%% Function name
%       tsaihill

%% Revised:
%       23 January 2014

%% Author
%       Srikanth Gunti, Trey Moore, & Autar Kaw
%       Section: All
%       Semester: Fall 2012

%% Purpose
%       Given elastic modulii and ultimate strengths of a unidirectional 
%       lamina, angle of the ply in degrees, and the global strains acting
%       on the lamina, output the strength ratio of the lamina based on
%       Tsai-Hill failure theory

%% Usage
%       funciton [SR] = tsaihill(strain_glo,angle,moduli,strength)
% 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
%       strain_glo=vector of global strain applied to unidirectional lamina
%       [strain_glo]=[epsx;epsy;epsxy]
%       epsx=longitudinal global strain
%       epsy=transverse global strain
%       epsxy=in-plane global strain
%       strength=vector of the five ultimate strain strengths of the 
%                unidirectional lamina
%       [strength]=[s1tu s1cu s2tu s2cu s12u] 
%       s1tu=ultimite longitudinal tensile strength
%       s1cu=ultimite longitudinal compressive strength
%       s2tu=ultimite transverse tensile strength
%       s2cu=ultimite transverse compressive strength
%       s12u=ultimite in-plane shear strength
% Output variable
%       SR=Strength Ratio
% Keyword
%       Tsai-Hill failure theory
%       strength ratio
%       angle ply

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

%% Code
function [SR] = tsaihill(strain_glo,angle,moduli,strength)
% Compliance matrix values 
S11=1/moduli(1);
S12=-1*moduli(3)/moduli(1);
S22=1/moduli(2);
S66=1/moduli(4);
% Sine of the angle of the lamina
s=sind(angle);
% Cosine of the angle of the lamina
c=cosd(angle);
% Transformed compliance matrix values 
S11bar=S11*(c^4)+(2*S12+S66)*(s^2)*(c^2)+S22*(s^4);
S12bar=S12*((s^4)+(c^4))+(S11+S22-S66)*(s^2)*(c^2);
S22bar=S11*(s^4)+(2*S12+S66)*(s^2)*(c^2)+S22*(c^4);
S16bar=(2*S11-2*S12-S66)*s*(c^3)-(2*S22-2*S12-S66)*(s^3)*c;
S26bar=(2*S11-2*S12-S66)*(s^3)*c-(2*S22-2*S12-S66)*s*(c^3);
S66bar=2*(2*S11+2*S22-4*S12-S66)*(s^2)*(c^2)+S66*((s^4)+(c^4));
% Transformed compliance matrix
Sbar=[S11bar S12bar S16bar; S12bar S22bar S26bar; S16bar S26bar S66bar];
% Transformed reduced stiffness matrix
Qbar=inv(Sbar);
% Global Stress
stress_glo=Qbar*strain_glo;
% Transformation matrix
T=[c^2 s^2 2*s*c; s^2 c^2 -2*s*c; -s*c s*c c^2-s^2];
% Local Stress
stress_loc=T*stress_glo;
% Tsai-Hill criterion
TSH=((stress_loc(1))/(strength(1)))^2 - (stress_loc(1))*(stress_loc(2))...
    /(strength(1))^2 + ((stress_loc(2))/(strength(3)))^2 + ((stress_loc(3))...
    /(strength(5)))^2;
%Strength ratio per Tsai-Hill criterion
SR=sqrt(1/TSH);
end