%% Function name
%       stressgloballaminate

%% Revised
%       3 February 2014

%% Authors
%       Dustin Smith, Eshan Patil, Trey Moore, & Autar Kaw 
%       Section: All
%       Semester: Fall 2013

%% Purpose 
%       Given the change in temperature, the change in moisture content, the number of angle plies, the qbar matrix for
%       each ply, the coefficients of moisture and thermal expansion, and the
%       global strains of all of the plies, output the global stresses for each ply 

%% Usage
%       function [stressglobalplies] = stressgloballaminate(strainglobalplies,
%                                      qbarplies,nplies,alphaplies,deltat,betaplies,deltac)
% Input variables
%       qbarplies=transformed reduced stiffness matrix for each ply
%       deltat=temperature change
%       alphaplies=coefficient of thermal expansion for each ply
%       betaplies=coefficient of moisture expansion for each ply
%       deltac=moisture change
%       nplies=number of plies
%       angleplies=angle of each ply
%       strainglobalplies=global strains at top, middle, and bottom of each laminate
% Output variables
%       stressglobalplies=global stress at top, middle, and bottom of each laminate
% Keyword
%       temperature change 
%       moisture content
%       mechancial strain
%       mechancial stress
%       hygrothermal strain
%       mositure strain
%       coefficents of thermal expansion
%       coefficents of moisture expansion
%       global strain
%       global stress

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

%% Code
function [stressglobalplies] = stressgloballaminate(strainglobalplies,qbarplies,nplies,alphaplies,deltat,betaplies,deltac)
for k=1:1:nplies
% Strain from differeint components
    for i=1:1:3
        for j=1:1:3
            strainply(j)=strainglobalplies(i,j,k);
            alpha(j)=alphaplies(j,k);
            beta(j)=betaplies(j,k);
        end
% Mechanical strain
       strainmechply=strainply-(alpha*deltat)-(beta*deltac);
       for j=1:1:3
            strainmechglobalplies(i,j,k)=strainmechply(j);
        end  
    end
end
 % Transformed reduced stiffness matrix for each ply
   for k=1:1:nplies
       for i=1:1:3
          for j=1:1:3
             qbar(i,j)=qbarplies(i,j,k);
          end
       end
       for i=1:1:3
    for j=1:1:3
            strainmechglobal(j)=strainmechglobalplies(i,j,k);
    end 
            stressglobal=qbar*strainmechglobal';
            stressglobalplies(i,1,k)=stressglobal(1);
            stressglobalplies(i,2,k)=stressglobal(2);
            stressglobalplies(i,3,k)=stressglobal(3);
       end
   end 
end