%% Function name
%       qqbar

%% Revised
%       30 January 2014

%% Author
%       Shyam Ramanath Thillainathan, Trey Moore, & Autar Kaw
%       Section: All
%       Semester: Fall 2013

%% Purpose
%       Given the four elastic moduli of plies, E1, E2, 
%       nu12,G12, angle of plies and number of plies this 
%       function calculates the Reduced Stiffness Matrix, Q 
%       and the Transformed Reduced Stiffness Matrix Qbar of the plies.

%% Usage
%       function [Qplies,Qbarplies] = qqbar(nplies,moduliplies,angleplies)
% Input variables
%       nplies=number of plies
%       moduliplies=vector of elastic modulii and Poisson's ratio for each lamina
%       E1=longitudinal elastic modulus
%       E2=transverse elastic modulus
%       nu12=major Poisson's ratio
%       G12=in-plane shear modulus
%       angleplies=angle of ply in degrees for each lamina   
% Output variables
%       Qplies=reduced stiffness matrix for each lamina
%       Qbarplies=transformed reduced stiffness matrix for each lamina
% Keyword
%       reduced stiffness matrix
%       transformed reduced stiffness matrix

%% License Agreement
%       http://www.eng.usf.edu/~kaw/OCW/composites/license/limiteduse.pdf

%% CODE
function [Qplies,Qbarplies] = qqbar(nplies,moduliplies,angleplies)
for k=1:1:nplies
    % Angle of each ply
    angle=angleplies(k);
    % Longitudinal elastic modulus
    E1=moduliplies(k,1);
    % Transverse elastic modulus
    E2=moduliplies(k,2);
    % Major Poisson's ratio
    nu12=moduliplies(k,3);
    % In-plane shear modulus
    G12=moduliplies(k,4);
    % Calculating the sine of the desired angle
    s=sind(angle);
    % Calculating the cosine of the desired angle
    c=cosd(angle);
    % Compliance matrix
    [S]=[1/E1 -nu12/E1 0;-nu12/E1 1/E2 0;0 0 1/G12];
    % Transformed compliance matrix values
    Sb11=(S(1,1)*c^4)+((2*S(1,2)+S(3,3))*s^2*c^2)+(S(2,2)*s^4);
    Sb12=(S(1,2)*(s^4+c^4))+((S(1,1)+S(2,2)-S(3,3))*s^2*c^2);
    Sb16=((2*S(1,1)-2*S(1,2)-S(3,3))*s*c^3)-((2*S(2,2)-2*S(1,2)-S(3,3))*s^3*c);
    Sb22=(S(1,1)*s^4)+((2*S(1,2)+S(3,3))*s^2*c^2)+(S(2,2)*c^4);
    Sb26=((2*S(1,1)-2*S(1,2)-S(3,3))*s^3*c)-((2*S(2,2)-2*S(1,2)-S(3,3))*s*c^3);
    Sb66=(2*(2*S(1,1)+2*S(2,2)-4*S(1,2)-S(3,3))*s^2*c^2)+(S(3,3)*(s^4+c^4));
    % Transformed compliance matrix
    Sbar=[Sb11 Sb12 Sb16;Sb12 Sb22 Sb26;Sb16 Sb26 Sb66];
    % Reduced Stiffness Matrix
    Q=inv(S);
    % Transformed Reduced Stiffness Matrix
    Qbar = inv(Sbar);
    % Reduced stiffness and transformed reduced stiffness matrices for each lamina
    for i=1:1:3
        for j=1:1:3
            Qplies(i,j,k) = Q(i,j);
            Qbarplies(i,j,k) = Qbar(i,j);
        end
    end
end