%% Function name
%       QbarandSbar

%% Revised
%       September 29, 2013

%% Author
%       Britton Steele, Christian Muneton, Trey Moore, & Autar Kaw
%       Section: All
%       Semester: Fall 2012

%% Purpose
%       Given the angle of the ply, the the elastic modulii, 
%       Poisson's ratio, and shear modulus of a unidirectional lamina,
%       output the Qbar and Sbar matricies

%% Usage
%       function [Qbar,Sbar] = QbarandSbar(angle,moduli)
% Input variables
%       moduli=vector with four elastic modulii of unidirectional lamina
%       [moduli]=[E1 E2 nu12 G12]
%       E1=longitudinal elastic modulus
%       E2=transverse elastic modulus
%       nu12=major Poisson's ratio
%       G12=in-plane shear modulus
%       angle=angle of ply given in degrees.
% Output variables
%       [Qbar]=vector of transformed reduced stiffness matrix
%       [Sbar]=transformed compliance matrix
% Keyword
%       transformed compliance matrix 
%       transformed reduced stiffness matrix

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

%% Code
function [Qbar,Sbar] = QbarandSbar(angle,moduli)
% Sine of the angle of the lamina
s=sind(angle);
% Cosine of the angle of the lamina
c=cosd(angle);
% transformed compliance matrix values 
Sb11=((1/moduli(1))*c^4)+((2*(-moduli(3)/moduli(1))+(1/moduli(4)))*s^2*c^2)+((1/moduli(2))*s^4);
Sb12=((-moduli(3)/moduli(1))*(s^4+c^4))+(((1/moduli(1))+(1/moduli(2))-(1/moduli(4)))*s^2*c^2);
Sb16=((2*(1/moduli(1))-2*(-moduli(3)/moduli(1))-(1/moduli(4)))*s*c^3)-((2*(1/moduli(2))-2*(-moduli(3)/moduli(1))-(1/moduli(4)))*s^3*c);
Sb22=((1/moduli(1))*s^4)+((2*(-moduli(3)/moduli(1))+(1/moduli(4)))*s^2*c^2)+((1/moduli(2))*c^4);
Sb26=((2*(1/moduli(1))-2*(-moduli(3)/moduli(1))-(1/moduli(4)))*s^3*c)-((2*(1/moduli(2))-2*(-moduli(3)/moduli(1))-(1/moduli(4)))*s*c^3);
Sb66=(2*(2*(1/moduli(1))+2*(1/moduli(2))-4*(-moduli(3)/moduli(1))-(1/moduli(4)))*s^2*c^2)+((1/moduli(4))*(s^4+c^4));
% transformed compliance matrix and transformed reduced stiffness matrix
Sbar=[Sb11 Sb12 Sb16;Sb12 Sb22 Sb26;Sb16 Sb26 Sb66];
Qbar=inv(Sbar);
end