clear
echo?on
d1=43;d2=7;dx=0.15;dy=0.1;xy=dx/dy;yx=dy/dx;
t=zeros（d1d2）;t1=ones（d1d2）;t0=zeros（d1d2）;
x=zeros（d1）;y=zeros（d2）;
x（1）=0;x（2）=dx/2;
for?i=3:d1-1
????x（i）=x（i-1）+dx;
end
x（d1）=（d1-2）*dx;
y（1）=0;y（2）=dy/2;
for?i=3:d2-1
????y（i）=y（i-1）+dy;
end
y（d2）=（d2-2）*dy;
t1=20*ones（d1d2）;
t=zeros（d1）;
dt=0.1;ttt=30;
%nnn=tt/dt;
echo?off
%for?iii=1:nnn
%????ttt=iii*dt;
????tf=30;af=6.6;af=1/af;bta=-6.6;v=0.0625;
????tin=100;tout=100;d=0.05;l=6;
????for?i=1:42
????????if?x（i）????????????t1（i1）=tin+（300-tout）*x（i）/（v*ttt）;
????????elseif?x（i）>v*ttt+28*d
????????????t1（i1）=300-（300-tout）*（x（i）-v*ttt-28*d）/（l-28*d-v*ttt）;
????????else?
????????????zz=-0.123*（x（i）-v*ttt）/d-3.52*exp（-0.123*（x（i）-v*ttt）/d）;
????????????t1（i1）=10060*exp（zz）;
????????end
????end
????for?i=1:41
????????for?j=2:6
????????????t1（ij）=t1（i1）-10*（j-1）;
????????end
????end
for?iii=1:500
????t0=t1;
????cd=1300;a=0.0003;
????an=[1.69-0.5940.401-0.1680.027-0.0370.046-0.05?0.039-0.012];
????bn=[0???0.333?0.017-0.1310.0540.00030.007-0.0120.0260];
????fncp=zeros（d1d2）;
????for?i=1:10
????????fncp=an（i）*cos（i*pi/9*t1）+bn（i）*sin（i*pi/9*t1）+fncp;
????end
????fncp=fncp/4.1868;
????b=1.3e-6;c=1.5e-9;
????fnk=（a+b*t1+c*t1.^2）*360;
????for?i=2:42
????????for?j=2:6
????????????fnae（ij）=2*yx*fnk（ij）*fnk（i+1j）/（fnk（ij）+fnk（i+1j））;
????????????fnaw（ij）=2*yx*fnk（ij）*fnk（i-1j）/（fnk（ij）+fnk（i-1j））;
????????????fnan（ij）=2*xy*fnk（ij）*fnk（ij+1）/（fnk（ij）+fnk（ij+1））;
????????????fnas（ij）=2*xy*fnk（ij）*fnk（ij-1）/（fnk（ij）+fnk（ij-1））;
????????end
????end
????fnap0=cd*fncp*dx*dy/dt;
????fnbb=fnap0.*t1;
????%for?i=2:41
????%????for?j=1:5
????%????????t1（ij）=t（i1）-10*（j-1）;
????%????end
????%?end
???%?t0=t1;
????kk=af+0.5*dx/fnk（12）;
????bb=fnbb（22）+tf*dy/kk;
????ap=fnae（22）+fnan（22）+fnap0（22）+dy/kk+2*xy*fnk（21）-bta*dx*dy;
????fff=fnae（22）*t1（32）+fnan（22）*t1（23）+2*xy*fnk（21）*t1（21）+bb;
????t1（22）=fff/ap;
????for?j=3;5
????????kk=af+0.5*dx/fnk（1j）;
????????bb=fnbb（2j）+tf*dy/kk;
????????ap=fnae（2j）+fnas（2j）+fnan（2j）+fnap0（2j）+dy/kk-bta*dx*dy;
????????fff=fnae（2j）*t1（3j）+fnan（2j）*t1（2j+1）+fnas（2j）*t1（2j-1）+bb;
????????t1（2j）=fff/ap;
????