%% 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

%% Testing code
clc
clear all
%% Inputs
% Material: Graphite/Epoxy
% Lamina: 30/45/60
strainglobalplies(:,:,1)=[-7.1666 7.1921 -0.0244;-4.7734 4.799 -0.0244;-2.3803 2.4058 -0.0244];
strainglobalplies(:,:,2)=[-2.3803 2.4058 -0.0244; 1.276E-02 1.276E-02 -0.0244; 2.4058 -2.3803 -0.0244];
strainglobalplies(:,:,3)=[2.4058 -2.3803 -0.0244; 4.799 -4.7734 -0.0244; 7.1921 -7.1666 -0.0244];
qbarplies(:,:,1)=[109.4 32.46 54.19; 32.46 23.65 20.05; 54.19 20.05 36.74]*10^9;
qbarplies(:,:,2)=[56.66 42.32 42.87; 42.32 56.66 42.87;42.87 42.87 46.59]*10^9;
qbarplies(:,:,3)=[23.65 32.46 20.05; 32.46 109.4 54.19; 20.05 54.1936 74]*10^9;
nplies=3;
tplies=0.125;
alphaplies=[5.64E-06 16.88E-06 -19.4676E-06; 11.26E-06 11.26E-06 -22.48E-06; 16.88E-06 5.64E-06 -19.4676E-06];
deltat=100;
betaplies=[0.15 0.45 -0.5196; 0.3 0.3 0.6; 0.45 0.15 -0.5196];
deltac=0.04;
fprintf('Number of plies: %g\n',nplies)
fprintf('Total Thicknes of Laminate: %g\n',tplies)
fprintf('Change in Temperature: %g\n',deltat)
fprintf('Change in Moisture: %g\n',deltac)
for k=1:1:nplies
    fprintf('\n\nFor Ply: %g\n\n',k)
    disp('   Transformed Reduced Stiffness Matrix')
    disp('______________________________________________')
    disp(qbarplies(:,:,k))
    disp('   Coefficient of Thermal Expansion')
    disp('______________________________________________')
    fprintf('    alphax  | %g\n    alphay  | %g\n    alphaxy | %g\n\n',alphaplies(:,k))
    disp('   Coefficient of Moisture Expansion')
    disp('______________________________________________')
    fprintf('    betax  | %g\n    betay  | %g\n    betaxy | %g\n\n',betaplies(:,k))
    disp('     Global Strains')
    disp('___________________________________________')
    for i=1:1:3
        x=i;
        if x==1
            disp('    Top')
        end
        if x==2
            disp('    Middle')
        end
        if x==3
            disp('    Bottom')
        end
        for j=1:1:3
           y=j;
        if y==1
            fprintf('       Epsx    | %G\n',strainglobalplies(i,j,k))
        end
        if y==2
            fprintf('       Epsx    | %G\n',strainglobalplies(i,j,k))
        end
        if y==3
            fprintf('       Gammaxy | %G\n',strainglobalplies(i,j,k))
        end 
        end
    end
end
%% Test
% Call function: stressglobalplies
stressglobalplies = stressgloballaminate(strainglobalplies,qbarplies,nplies,alphaplies,deltat,betaplies,deltac);
for k=1:1:nplies
    fprintf('\n\nFor Ply: %g\n\n',k)
    disp('     Global Stresses')
    disp('___________________________________________')
    for i=1:1:3
        x=i;
        if x==1
            disp('    Top')
        end
        if x==2
            disp('    Middle')
        end
        if x==3
            disp('    Bottom')
        end
        for j=1:1:3
           y=j;
        if y==1
            fprintf('       Sx    | %G\n',stressglobalplies(i,j,k))
        end
        if y==2
            fprintf('       Sx    | %G\n',stressglobalplies(i,j,k))
        end
        if y==3
            fprintf('       Txy   | %G\n',stressglobalplies(i,j,k))
        end 
        end
    end
end