%% Function name
%       tsaiwu

%% Revised:
%       23 January 2014

%% Author
%       Abiodun Yakub, 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-Wu failure theory

%% Usage
%       function [SR] = tsaiwu(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-Wu failure theory
%       strength ratio
%       angle ply

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

%% Code
function [SR] = tsaiwu(strain_glo,angle,moduli,strength)
% Tsai-Wu failure theory parameters
H1=(1/((strength(1))) - (1/(strength(2))));
H2=(1/((strength(3))) - (1/(strength(4))));
H11=(1/((strength(1)*strength(2))));
H22=(1/((strength(3)*strength(4))));
H66=(1/((strength(5))^2));
H6=0;
% Finding H12 using Mises-Hencky criterion
H12=(-1/2)*sqrt(1/((strength(1))*(strength(2))*(strength(3))*(strength(4))));
% 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-Wu criterion in quadratic form.
b=(H1*stress_loc(1))+(H2*stress_loc(2))+(H6*stress_loc(3));
a=(H11*(stress_loc(1)^2)+(H22*(stress_loc(2)^2)+(H66*(stress_loc(3)^2)+(2*H12*stress_loc(1)*stress_loc(2)))));
c=-1;
% Quadratic formula
R(1)=((-b)+sqrt((b^2)-(4*a*c)))/(2*a);
R(2)=((-b)-sqrt((b^2)-(4*a*c)))/(2*a);
% Strength ratio per Tsai-Wu criterion
if(R(1)<0)
    SR=R(2);
else
    SR=R(1);
end
end