clc;
length=500e-6;
width=100e-6;
tal=10e-6;tsio2=10e-6;%thickness
m=tsio2/tal; % ratio of thickness
L=linspace(0,length,11);
cpal=900;cpsio2=703; %Specific Heat capacity
gal=160;gsio2=1.38;% thermal conductivity
pal=2700;psio2=2200; % density (Kg/m^3)
aal=23e-6;asio2=0.55e-6;% THermal expansion coefficient
Eal=70e9;Esio2=73.1e9; %Young's Modulus
n=Esio2/Eal; % ratio of young's modulus
time=linspace(0,28e-3,29);
dz_sim=[8.18E-08	1.25E-06	1.53E-06	1.71E-06	1.83E-06	1.92E-06	1.97E-06	2.00E-06	2.03E-06	2.04E-06	2.06E-06	2.06E-06	2.07E-06	2.07E-06	2.07E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06	2.08E-06
];
T_tip=[3.30000000000000e-18,3.26804100000000,8.53223300000000,12.2809725000000,14.8053250000000,16.4773940000000,17.5580220000000,18.2597890000000,18.7835160000000,19.1567360000000,19.4120770000000,19.5605800000000,19.6689430000000,19.7284450000000,19.7655980000000,19.7884480000000,19.8045770000000,19.8096940000000,19.8065850000000,19.8032460000000,19.8000870000000,19.7985080000000,19.7969280000000,19.7956000000000,19.7948550000000,19.7941070000000,19.7935490000000,19.7931750000000,19.7927990000000];
%plotyy(time,T_tip, time, dz_tip);
for j=1:11;
    massal(j)=L(j)*width*tal*pal; %mass
masssio2(j)=L(j)*width*tsio2*psio2; %mass
Gcal(j)=gal*width*tal/L(j); Gcsio2(j)=gsio2*width*tsio2/L(j);%thermal conductance
;
end
Cal=massal*cpal;
Csio2=masssio2*cpsio2; %heat capacitance
C=(Cal+Csio2)/2; % total equivalent capacitance
%Gcal=gal*width*tal/length; Gcsio2=gsio2*width*tsio2/length;%thermal conductance
G=Gcal+Gcsio2; % total equivalent thermal conductance.

Thot=313; Tcold=293;
h=tal+tsio2; %height
tl2h=3*length^2/(2*h);
opm2=(1+m)^2;
opmn=1+m*n;
msp1mn=m^2+1/(m*n);
da=aal-asio2;
delT=zeros(2,29); Tana=zeros(2,29); Tsim=zeros(2,29); dzana=zeros(2,29);
dz(1,1)=0;
for t=1:29
    for k=1:11
        tau(k)=C(k)/G(k);
        dummy2=(t-1)*1e-3/tau(k);
        dummy=exp(-dummy2);
        delT(k,t)=(Thot-Tcold)*dummy;
        Tana(k,t)=Thot-delT(k,t);
        Tsim(k,t)=Thot-delT(k,t);
        dzana(k,t)=tl2h*opm2*(aal-asio2)*(Tana(k,t)-Tcold)/(3*opm2+opmn*msp1mn);
        dztip(k,t)=0.5e-3*dzana(k,t)/L(k);
        if (t*1e-3)<tau(k); 
           dz(k,t)=dzana(k,t);
        elseif t==1 
            dz(k,t)=dz(k,t);
        else t>1 && (t*1e-3)>tau(k);
            dz(k,t)=dz(k,t-1);
        end
    end
end
plot(time,Tana,time,Tsim,'rd')
figure;
plot(time,dz(1:29),time,dz_sim,'rd')