%% Function name
%       anglemoduli

%% Revised:
%       13 January 2014

%% Author
%       Peter Vinje, Trey Moore, and Autar Kaw
%       Section: All
%       Semester: Fall 2012

%% Purpose
%       Given elastic modulii of a unidirectional lamina and the angle of the 
%       ply in degrees, output elastic modulii in the x and y directions, 
%       Poisson's ratio in the x-y plane, the shear modulus in the x-y
%       plane, as well as the shear coupling terms with axial load in both
%       the x and y directions

%% Usage
%       function [moduli_ang] = anglemoduli(moduli,angle)
% Input variables
%       moduli=vector with four elastic moduli of unidirectional lamina
%       [moduli]=[E1 E2 nu12 G12]
%       angle=angle of ply in degrees
%       E1=longitudinal elastic modulus
%       E2=transverse elastic modulus
%       nu12=major Poisson's ratio
%       G12=in-plane shear modulus
% Output variables
%       Ex=elastic modulus in the x-direction
%       Ey=elastic modulus in the y-direction
%       nuxy=Poisson's ratio in the x-y plane
%       Gxy=shear modulus in x-y plane
%       Mx=shear coupling with axial load in x-direction
%       My=shear coupling with axial load in y-direction
%       moduli=vector with four elastic moduli of unidirectional lamina in
%              x,y direction
%       [moduli_ang]=[nuxy Gxy Mx My]
% Keywords
%       elastic modulus
%       Poisson's ratio
%       shear modulus
%       shear coupling
%       angle ply

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

%% Code
function [moduli_ang] = anglemoduli(moduli,angle)
% Calculating the sine of the desired angle
s=sind(angle);
% Calculating the cosine of the desired angle
c=cosd(angle);
% Creating the S matrix from the moduli vector
S=[1/moduli(1) -moduli(3)/moduli(1) 0; -moduli(3)/moduli(1) 1/moduli(2) 0;
    0 0 1/moduli(4)];
% Defining elements of the transformed reduced compliance matrix, Sbar
sbar11 = S(1,1)*c^4 + (2*S(1,2)+S(3,3))*s^2*c^2 + S(2,2)*s^4;
sbar12 = S(1,2)*(s^4+c^4) + (S(1,1)+S(2,2)-S(3,3))*s^2*c^2;
sbar22 = S(1,1)*s^4 + (2*S(1,2)+S(3,3))*s^2*c^2 + S(2,2)*c^4;
sbar16 = (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;
sbar26 = (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;
sbar66 = 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);
% Constructing the Sbar matrix from the elements
Sbar =[sbar11 sbar12 sbar16; sbar12 sbar22 sbar26; sbar16 sbar26 sbar66];
% Engineering constants of the angle ply
Ex = 1/Sbar(1,1);
Ey = 1/Sbar(2,2);
nuxy = -Sbar(1,2)/Sbar(1,1);
Gxy = 1/Sbar(3,3);
Mx = -Sbar(1,3)*moduli(1);
My = -Sbar(2,3)*moduli(1);
% Elements of engineering constants create the vector moduli_ang matrix
[moduli_ang] = [Ex Ey nuxy Gxy Mx My];
end