%1.Parameters
promptEf='Value of Ef=';
promptvf='Value of vf=';
promptEm='Value of Em=';
promptvm='Value of vm=';
promptVf='Value of Vf=';
promptVm='Value of Vm=';
promptGf='Value of Gf=';
promptGm='Value of Gm=';
Ef=input(promptEf);
vf=input(promptvf);
Em=input(promptEm);
vm=input(promptvm);
Vf=input(promptVf);
Vm=input(promptVm);
Gf=input(promptGf);
Gm=input(promptGm);
E11=Ef*Vf+Em*Vm;
v12=vf*Vf+vm*Vm;
E22=Ef*Em/(Ef*Vm+Em*Vf);
v21=E22*v12/E11;
G12=Gf*Gm/(Gf*Vm+Gm*Vf);
%2&3-Number of layers,Thickness
prompt_thickness='layer thickness=';
t=input(prompt_thickness);
promptN = 'Number of layers=';
N=input(promptN);
orientation=zeros(1,N);
fprintf('Enter values of layer orientation from bottom to top\n')
%Orientation of fibers
 for i=1:N
 fprintf('Angle of fibers in layers %i\n',i);
 prompt_theta='Theta=';
 theta=input(prompt_theta);
 orientation(1,i)=theta;
 end
disp(orientation);
%4.Qij matrix
Q11=E11/(1-v12*v21);
Q22=E22/(1-v12*v21);
Q12=v12*E22/(1-v12*v21);
Q21=v21*E11/(1-v12*v21);
Q66=G12;
Q=[Q11 Q12 0;Q21 Q22 0;0 0 Q66];
fprintf('Matrix Q=\n');
disp(Q);
R=[1 0 0;0 1 0;0 0 2];
Rinv=inv(R);
%loop starting for different layers
m=cell(N,1);
for n=1:N
 c=cosd(orientation(1,n));%degree
 s=sind(orientation(1,n));%degree
 %Transformation matrix
 a11=c.^2; a12=s.^2; a13=2*c*s;
 a21=s.^2; a22=c.^2; a23=(-2*c*s);
 a31=(-c*s); a32=c*s; a33=c.^2-s.^2;
 T=[a11 a12 a13;a21 a22 a23;a31 a32 a33];
 Tinv=inv(T);
 %5.Qij-bar matrix
 Q_bar=Tinv*Q*R*T*Rinv;
 %Q_bar=R*inv(T)*Rinv*Q*T;
 fprintf('Q_bar Matrix=\n');
 disp(Q_bar);
 m{n}=Q_bar;
 end
%6.z values
z=zeros(1,(N+1));
h=N*t;
for j=1:N+1
 z(1,j)=h/2-t*(j-1);
end
%disp(z);
aa=0;
bb=0;
dd=0;
A=zeros(3,3);
B=zeros(3,3);
D=zeros(3,3);
for yy=1:3
 for zz=1:3
 for xx=1:n
 aa=aa+(m{xx}(yy,zz))*(z(1,(xx))-z(1,xx+1));
 bb=bb+(m{xx}(yy,zz))*(z(1,(xx))^2-(z(1,xx+1))^2);
 dd=dd+(m{xx}(yy,zz))*(z(1,(xx))^3-(z(1,xx+1))^3);
 end
 A(zz,yy)=aa;
 B(zz,yy)=0.5*bb;
 D(zz,yy)=(1/3)*dd;
 aa=0;
 bb=0;
 dd=0;
 end
end
%7. ABD Matrix
ABBD=[A(1,1) A(1,2) A(1,3) B(1,1) B(1,2) B(1,3);
 A(2,1) A(2,2) A(2,3) B(2,1) B(2,2) B(2,3);
 A(3,1) A(3,2) A(3,3) B(3,1) B(3,2) B(3,3);
 B(1,1) B(1,2) B(1,3) D(1,1) D(1,2) D(1,3);
 B(2,1) B(2,2) B(2,3) D(2,1) D(2,2) D(2,3);
 B(3,1) B(3,2) B(3,3) D(3,1) D(3,2) D(3,3)];
fprintf('Matrix ABBD= \n');
disp(ABBD);
A = [A(1,1) A(1,2) A(1,3);A(2,1) A(2,2) A(2,3);A(3,1) A(3,2) A(3,3)];
B = [B(1,1) B(1,2) B(1,3);B(2,1) B(2,2) B(2,3);B(3,1) B(3,2) B(3,3)];
D = [D(1,1) D(1,2) D(1,3); D(2,1) D(2,2) D(2,3);D(3,1) D(3,2) D(3,3)];
Ainv = inv(A);
Bx = (B)*(Ainv)*(B);
Dstar = (D)- Bx;
D1 = inv(Dstar);
B1 = -(Ainv)*(B)*(D1);
C1 = -(D1)*(B)*(Ainv);
A1 = (Ainv)+((Ainv)*(B)*(Dstar)*(B)*(Ainv));
fprintf ('matrix A1= \n');
disp (A1);
fprintf ('matrix B1= \n');
disp (B1);
fprintf ('matrix C1= \n');
disp (C1);
fprintf ('matrix D1= \n');
disp (D1);