clc;
clear all;
%responsivity calculation of multilayer cantilever beam
%material =[SU8 Al SiO2]; % from bottom to top.
E=[3.2e9 70e9 73.1e9];% Young's modulus [Pa]
rho=[1.5e6 2.7e6 2.22e6];%density [g/m^3]
alfa=[8.71e-5 23e-6 0.55e-6];%thermal expansion coefficient [m/K]
Gc=[0.2 160 1.38];% Thermal conductivity [W/m-k]
Cp=[2 0.9 0.703];%Specific heat capacity [J/g-K]
%Thickness matrix
t=[7.5e-6 50e-9 1e-6; 15e-6 100e-9 2e-6; 22.5e-6 150e-9 3e-6; 30e-6 200e-9 4e-6; 37.5e-6 250e-9 5e-6; 45e-6 300e-9 6e-6];
tal=[50e-9 100e-9 150e-9 200e-9 250e-9 300e-9];%Thickness of aluminum [m]
tsio2=[1e-6 2e-6 3e-6 4e-6 5e-6 6e-6];%Thickness of SiO2 [m]
tsu8=[7.5e-6 15e-6 22.5e-6 30e-6 37.5e-6 45e-6];%Thickness of SU8[m]
L=[0.5e-3 1e-3 2e-3 3e-3 4e-3 5e-3 6e-3 7e-3];%[m]
w=L/4;%[m]
% calculating radius of curvature.
%find Ai
 sumaEtw=0;sumEtw=0;sumCEtw=0;
 C(1)=t(6,1);
  for i=2:3
        C(i)=C(i-1)+t(6,i);
        CEtw=C(i)*E(i)*t(6,i)*w(1);
        sumCEtw=sumCEtw+CEtw;
  end
for i=1:3
        aEtw=alfa(i)*E(i)*t(6,i)*w(1);
        Etw=E(i)*t(6,i)*w(1);
        sumaEtw=sumaEtw+aEtw;
        sumEtw=sumEtw+Etw;
end
    dT=20; %[k] Temperature difference
for i=1:3
    A(i)=(-alfa(i)+sumaEtw/sumEtw)*dT*E(i)*t(6,i)*w(1);
    B(i)=(C(i)-sumCEtw/sumEtw)*E(i)*t(6,i)*w(1);
end
 sumEI=0;
for i=1:3
    I(i)=(t(6,i)^3)/12;%calculating Moment I for each layer
    EI=E(i)*I(i);% and multiplying with Young's modulus
    sumEI=sumEI+EI; %summing
end 
sumAC=0;sumBC=0;
for i=2:3
    AC=A(i)*C(i);
    sumAC=sumAC+AC;
    BC=B(i)*C(i);
    sumBC=sumBC+BC;
end
r=(2*sumAC)/(0+sumBC);% radius of curvature
dy=((L(1)^2)/2)*r % tip deflection


%finding thermal response
% for i=1:length(tal)
%     for j=1:length(L)
%         %Aluminum
%         Aal(i,j)=tal(i)*w(j);
%         Val=Aal*L(j);
%         Mal=Val*rho(2);
%         Thcapal=Mal*Cp(2);
%         thGal=Gc(2)*Aal/L(j);
%         %Su8
%         Asu8(i,j)=tsu8(i)*w(j);
%         Vsu8=Asu8*L(j);
%         Msu8=Vsu8*rho(1);
%         Thcapsu8=Msu8*Cp(1);
%         thGsu8=Gc(1)*Asu8/L(j);
%         %SiO2
%         Asio2(i,j)=tsio2(i)*w(j);
%         Vsio2=Asio2*L(j);
%         Msio2=Vsio2*rho(3);
%         Thcapsio2=Msio2*Cp(3);
%         thGsio2=Gc(3)*Asio2/L(j);
%         
%     end
% end
%time constant calculation
% Ceq=Thcapal+Thcapsu8+Thcapsio2;%Equivalent thermal capacitance
% Geq=thGal+thGsu8+thGsio2; %Equivalent thermal cconductance
% for i=1:6
%     for j=1:8
%         Tau_al(i,j)=Ceq(i,j)/Geq(i,j);
%     end
% end