%% Function name
%      ABD

%% Revised
%      28 January 2014

%% Author
%      Brent Savage, Willie Valentin, Trey Moore, & Autar Kaw
%      Section: All
%      Semester: Fall 2013

%% Purpose
%       Given the number of plies, the reduced stiffness matrix for each 
%       ply, the angle of the ply in degrees, and the thickness of each ply
%       output the extensional, coupling, bending, normalized extensional,
%       normalized coupling, and normalized bending matrices

%% Usage
%      function [A,B,D,An,Bn,Dn] = ABD(nplies,Qplies,angleplies,tplies)
% Input variables
%       nplies=number of plies
%       Qplies=reduced stiffness matrix for each ply
%       angleplies= angle of each ply in degrees
%       tplies=thickness of each ply
% Output variables
%       A=extensional stiffness matrix
%       B=coupling stiffness matrix
%       D=bending stiffness matrix
%       An=normalized extensional stiffness matriix
%       Bn=normalized coupling stiffness matrix
%       Dn=normalized bending stiffness matrix
% Keywords
%       extensional stiffness matrix 
%       coupling stiffness matrix
%       bending stiffness matrix
%       normalized extensional stiffness matrix
%       normalized coupling stiffness matrix
%       normalized bending stiffness matrix

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

%% Code
function[A,B,D,An,Bn,Dn]=ABD(nplies,Qplies,angleplies,tplies)
% Three dimensional reduced stiffness matrix consisting of a reduced
% stiffness matrix for each ply
for i=1:1:nplies
    Qbarplies(1,1,i)=(Qplies(1,1,i)*((cosd(angleplies(1,i)))^4))...
        +(Qplies(2,2,i)*((sind(angleplies(1,i)))^4))+(2*(Qplies(1,2,i)...
        +(2*Qplies(3,3,i)))*((sind(angleplies(1,i)))^2)...
        *((cosd(angleplies(1,i)))^2));
    Qbarplies(1,2,i)=((Qplies(1,1,i)+Qplies(2,2,i)-4*Qplies(3,3,i))...
        *((sind(angleplies(1,i)))^2)*((cosd(angleplies(1,i)))^2))...
        +(Qplies(1,2,i)*(((cosd(angleplies(1,i)))^4)...
        +((sind(angleplies(1,i)))^4)));
    Qbarplies(2,1,i)=((Qplies(1,1,i)+Qplies(2,2,i)-4*Qplies(3,3,i))...
        *((sind(angleplies(1,i)))^2)*((cosd(angleplies(1,i)))^2))...
        +(Qplies(1,2,i)*(((cosd(angleplies(1,i)))^4)...
        +((sind(angleplies(1,i)))^4)));
    Qbarplies(2,2,i)=((Qplies(1,1,i)*((sind(angleplies(1,i)))^4))...
        +(Qplies(2,2,i)*((cosd(angleplies(1,i)))^4))+(2*(Qplies(1,2,i)...
        +2*Qplies(3,3,i))*((sind(angleplies(1,i)))^2)...
        *((cosd(angleplies(1,i)))^2)));
    Qbarplies(1,3,i)=(Qplies(1,1,i)-Qplies(1,2,i)-2*Qplies(3,3,i))...
        *cosd(angleplies(1,i))^3*sind(angleplies(1,i))...
        -(Qplies(2,2,i)-Qplies(1,2,i)-2*Qplies(3,3,i))...
        *sind(angleplies(1,i))^3*cosd(angleplies(1,i));
    Qbarplies(3,1,i)=(Qplies(1,1,i)-Qplies(1,2,i)-2*Qplies(3,3,i))...
        *cosd(angleplies(1,i))^3*sind(angleplies(1,i))-(Qplies(2,2,i)...
        -Qplies(1,2,i)-2*Qplies(3,3,i))*(sind(angleplies(1,i))^3)...
        *cosd(angleplies(1,i));
    Qbarplies(2,3,i)=((Qplies(1,1,i)-Qplies(1,2,i)-2*Qplies(3,3,i))...
        *cosd(angleplies(1,i))*((sind(angleplies(1,i)))^3))...
        -((Qplies(2,2,i)-Qplies(1,2,i)-2*Qplies(3,3,i))...
        *((cosd(angleplies(1,i)))^3)*sind(angleplies(1,i)));
    Qbarplies(3,2,i)=((Qplies(1,1,i)-Qplies(1,2,i)-2*Qplies(3,3,i))...
        *cosd(angleplies(1,i))*((sind(angleplies(1,i)))^3))...
        -((Qplies(2,2,i)-Qplies(1,2,i)-2*Qplies(3,3,i))...
        *((cosd(angleplies(1,i)))^3)*sind(angleplies(1,i)));
    Qbarplies(3,3,i)=((Qplies(1,1,i)+Qplies(2,2,i)-2*Qplies(1,2,i)...
        -2*Qplies(3,3,i))*((sind(angleplies(1,i)))^2)...
        *((cosd(angleplies(1,i)))^2))+(Qplies(3,3,i)...
        *(((sind(angleplies(1,i)))^4)+((cosd(angleplies(1,i)))^4)));
end
% Total thickness of the laminate
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
% Extensional stiffness matrix, coupling stiffness matrix and bending stiffness matrix
A=0;
d=0;
b=0;
for i=1:1:nplies
    h1(1,i)=(h(i+1)-h(i));
    h2(1,i)=(h(i+1)^2-h(i)^2);
    h3(1,i)=(h(i+1)^3-h(i)^3);
    A=A+(Qbarplies(:,:,i)*h1(1,i));
    b=b+(Qbarplies(:,:,i)*h2(1,i));
    d=d+(Qbarplies(:,:,i)*h3(1,i));
end
% Multiplication due to integration
B=0.5*b;
D=(1/3)*d;    
% Normalized extensional stiffness matrix
An=A*(1/t);
% Normalized coupling stiffness matrix
Bn=B*(2/t^2);
% Normalized bending stiffness matrix
Dn=D*(12/t^3);
end