%% Function Name
%       NTMT

%% Revised
%       30 January 2014

%% Author
%       Srujana Anne, Trey Moore, & Autar Kaw
%       Semester: Fall 2013

%% Purpose
%       Given the transformed reduced stiffness matrix for each ply, the
%       temperature change , the number of plies, the thickness of plies,
%       and the coefficients of thermal expansion for each lamina. output
%       the fictitious thermal loads

%% Usage
%       function [NT,MT] = NTMT(qbarplies,deltat,alphaplies,nplies,tplies)
% Input Variables
%       nplies=number of plies
%       qbarplies=transformed reduced stiffness matrix for each ply
%       deltat=temperature change
%       alphaplies=coefficients of thermal expansion for each ply
%       tplies=thickness of each ply
% Output Variables
%       NT=vector of fictitious thermal forces
%       MT=vector of fictitious thermal moments
% Keywords
%       transformed reduced stiffness matrix
%       temperature change
%       coefficients of thermal expansion
%       thickness of plies
%       fictitious thermal forces
%       fictitious thermal moments

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

%% Code
function [NT,MT]=NTMT(qbarplies,deltat,alphaplies,nplies,tplies)
% Thickness of lamina
t=sum(tplies);
% Position of top and bottom of each lamina
h(1)=-t/2;
for i=2:1:nplies+1
    h(i)=h(i-1)+tplies(i-1);
end
% Fictitious thermal loads 
NT=0;
MT=0;
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
        alpha(i)=alphaplies(i,k);
    end
    NT=NT+(deltat*qbar*alpha'*(h(k+1)-h(k)));
    MT=MT+(0.5*deltat*qbar*alpha'*(h(k+1)^2-h(k)^2));
end
end