AS4/PEEK-PEEK liner: 
(1)=PEEK; (2)=AS4/PEEK: 
E1(1)=3.64e9;E2(1)=3.64e9;v12(1)=0.4;G12(1)=1.3e9;v21(1)=v12(1)*E2(1)/E1(1);
E1(2)=131e9;E2(2)=8.7e9;v12(2)=0.28;G12(2)=5e9;v21(2)=v12(2)*E2(2)/E1(2);

AS4/Epoxy-steel liner:
(1)=steel; (2)=AS4/epoxy:
E1(1)=2.07e11;E2(1)=2.07e11;v12(1)=0.3;G12(1)=7.96154e10;v21(1)=v12(1)*E2(1)/E1(1);
E1(2)=1.354e11;E2(2)=9.37e9;v12(2)=0.32;G12(2)=4.96e9;v21(2)=v12(2)*E2(2)/E1(2);


AS4/Epoxy-Ti liner:
(1)=Ti; (2)=AS4/epoxy:
E1(1)=1.138e11;E2(1)=1.138e11;v12(1)=0.342;G12(1)=4.24e10;v21(1)=v12(1)*E2(1)/E1(1);
E1(2)=1.354e11;E2(2)=9.37e9;v12(2)=0.32;G12(2)=4.96e9;v21(2)=v12(2)*E2(2)/E1(2);


AS4/Epoxy-Al liner:
(1)=Al; (2)=AS4/epoxy:
E1(1)=7.1e10;E2(1)=7.1e10;v12(1)=0.3;G12(1)=2.731e10;v21(1)=v12(1)*E2(1)/E1(1);
E1(2)=1.354e11;E2(2)=9.37e9;v12(2)=0.32;G12(2)=4.96e9;v21(2)=v12(2)*E2(2)/E1(2);



E1(1)=1.138e11;E2(1)=1.138e11;v12(1)=0.342;G12(1)=4.24e10;v21(1)=v12(1)*E2(1)/E1(1);
E1(2)=1.354e11;E2(2)=9.37e9;v12(2)=0.32;G12(2)=4.96e9;v21(2)=v12(2)*E2(2)/E1(2);
format long
ABD_global_total=zeros(6,6);Q_original=zeros(4,1);Q=zeros(6,1);
N=input('number of layers=');M=input('number of material=');t_total=0;z=0;
for m=1:1:N;
    t(m)=input('t=');
   angle(m)=input('angle=');
  material(m)=input('material=');
   t_total=t_total+t(m);
end
z(1)=-(t_total/2-t(1)/2);
for k=2:1:N;
  z(k)=z(k-1)+t(k-1)/2+t(k)/2;
end

for m=1:1:N;
    s(m)=sin(angle(m)*pi/180);
    c(m)=cos(angle(m)*pi/180);
end

for i=1:1:N;
T(:,:,i)=[(c(i))^4           2*(c(i))^2*(s(i))^2             (s(i))^4          4*(c(i))^2*(s(i))^2
       (c(i))^2*(s(i))^2     (c(i))^4+(s(i))^4             (c(i))^2*(s(i))^2    -4*(c(i))^2*(s(i))^2
          (s(i))^4           2*(c(i))^2*(s(i))^2             (c(i))^4          4*(c(i))^2*(s(i))^2
        (c(i))^3*(s(i))    -c(i)*s(i)*((c(i))^2-(s(i))^2)   -(c(i))*(s(i))^3   -2*c(i)*s(i)*((c(i))^2-(s(i))^2)
        (c(i))*(s(i))^3     c(i)*s(i)*((c(i))^2-(s(i))^2)   -(c(i))^3*(s(i))   2*c(i)*s(i)*((c(i))^2-(s(i))^2)
       (c(i))^2*(s(i))^2    -2*(c(i))^2*(s(i))^2           (c(i))^2*(s(i))^2   ((c(i))^2-(s(i))^2)^2];
end

for i=1:1:N;                         %%%%%%%%%%%%%%%%%using if commond to add different materials,in this example 8 materials can be defined.
        if  material(i)==1
       Q11_original(i)=E1(1)/(1-v12(1)*v21(1));
       Q22_original(i)=E2(1)/(1-v12(1)*v21(1));
       Q12_original(i)=v21(1)*E1(1)/(1-v12(1)*v21(1));
       Q66_original(i)=G12(1);
        end
                
        if material(i)==2
       Q11_original(i)=E1(2)/(1-v12(2)*v21(2));
       Q22_original(i)=E2(2)/(1-v12(2)*v21(2));
       Q12_original(i)=v21(2)*E1(2)/(1-v12(2)*v21(2));
       Q66_original(i)=G12(2);
        end
        
end       
for i=1:1:N
          Q_original(:,:,i)=[Q11_original(i)
                  Q12_original(i)
                  Q22_original(i)
                  Q66_original(i)];
       Q(:,:,i)=T(:,:,i)*Q_original(:,:,i);
end


A11=0; A12=0;  A22=0; A16=0; A26=0; A66=0;
B11=0; B12=0;  B22=0; B16=0; B26=0; B66=0;
D11=0; D12=0;  D22=0; D16=0; D26=0; D66=0;
for i=1:N
    A11=A11+Q(1,1,i)*t(i);
    A12=A12+Q(2,1,i)*t(i);
    A22=A22+Q(3,1,i)*t(i);
    A16=A16+Q(4,1,i)*t(i);
    A26=A26+Q(5,1,i)*t(i);
    A66=A66+Q(6,1,i)*t(i);
    
    
    B11=B11+Q(1,1,i)*t(i)*z(i);
    B12=B12+Q(2,1,i)*t(i)*z(i);
    B22=B22+Q(3,1,i)*t(i)*z(i);
    B16=B16+Q(4,1,i)*t(i)*z(i);
    B26=B26+Q(5,1,i)*t(i)*z(i);
    B66=B66+Q(6,1,i)*t(i)*z(i);
    
    D11=D11+Q(1,1,i)*(t(i)*(z(i))^2+(t(i))^3/12);
    D12=D12+Q(2,1,i)*(t(i)*(z(i))^2+(t(i))^3/12);
    D22=D22+Q(3,1,i)*(t(i)*(z(i))^2+(t(i))^3/12);
    D16=D16+Q(4,1,i)*(t(i)*(z(i))^2+(t(i))^3/12);
    D26=D26+Q(5,1,i)*(t(i)*(z(i))^2+(t(i))^3/12);
    D66=D66+Q(6,1,i)*(t(i)*(z(i))^2+(t(i))^3/12);
end    
    
ABD_global_total=[A11  A12  A16  B11  B12  B16
                  A12  A22  A26  B12  B22  B26      
                  A16  A26  A66  B16  B26  B66 
                  B11  B12  B16  D11  D12  D16
                  B12  B22  B26  D12  D22  D26     
                  B16  B26  B66  D16  D26  D66 ];
              
              


Ex=(A11*A22-A12^2)/A22/t_total;
Ey=(A11*A22-A12^2)/A11/t_total;
Gxy=A66/t_total;
v_xy=A12/A22;
Exb=12*(D11*D22-D12^2)/t_total^3/D22;
Gxyb=12*D66/t_total^3;
Ex
Exb
Ey
Gxy
Gxyb
v_xy