%Conduction of heat in a slab with convection at both surfaces
% solved by finite differeences in space and state variables in time,
% using 21 nodes counting the faces of the slab.
% This m-file was written at the University of Wyoming in the Electrical
% and Computer Engineering Department and is to be distributed without
% cost.
clear all
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
figure(1);clf
figure(2);clf
figure(3);clf
figure(4);clf
figure(5);clf
xoverL=linspace(0,1,21);
xoverL=[-.1 xoverL 1.1];
N=20;% 19 internal nodes
h1Loverk=1;
h2Loverk=1;
A=zeros(21,21);
B=zeros(21,1);
B(1,1)=2*h1Loverk*N;
for i=2:20
    A(i,i)=-2*N*N;
end
for i=2:20
    A(i,i+1)=N*N;
    A(i+1,i)=N*N;
end
A(1,1)=-2*N*N-2*N*h1Loverk;
A(1,2)=2*N*N;
A(2,1)=N*N;


A(21,21)=-2*N*N-2*N*h2Loverk;
A(21,20)=A(1,2);
%A(20,21)=A(2,1);
T=.001;
[Ad,Bd]=c2d(A,B,T);
x=zeros(21,1);
framedata=x';
for k=1:2000
    x=Ad*x+Bd;
    if fix(k/10)==k/10;
        framedata=[framedata;x'];
    end
end
framedata=[ones(201,1) framedata zeros(201,1)];
framedata(1,1)=0;


figure(1);% Plot T(x,t) vs x for various t
f1=[ 0 0.08 .2 .31 .43 .51 .58 .64 .76 .88 1];
xpatch1=[f1 ones(1,11)-f1 f1(1,1)];
y1=[-.1 -.04 -.16 -.12 -.05 -.14 -.07 -.11 -.16 -.12 -.08];
ypatch1=[y1 1.1*ones(1,11)+(y1) y1(1,1)];
patch(xpatch1,ypatch1,'c')
hold on
plot(xoverL,framedata(1,:),'k')
hold on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Temperature, T(x,t)/T_0')
for k=11:20:201
    plot(xoverL,framedata(k,:),'k')
    hold on
    axis([-.2 1.2 -.2 1.2])
end
text(.05,.05,['\kappat/L^2 = ' num2str(0)])
text(.18,.2,num2str(0.2))
text(.22,.6,num2str(2))
text(1.04,.8,'T = 0')
text(-.15,.8,'T = T_0')
set(gca,'Box','on')
text(.5,.7,['h_1L/k = ' num2str(h1Loverk)]);
text(.5,.6,['h_2L/k = ' num2str(h2Loverk)]);
text(-.05,.4,'h_1');
text(1.02,.4,'h_2');
text(.1,1.13,'Press Enter to Continue')
pause

figure(2)%Plot T(x,t) vs t for various x
tau=linspace(0,2,201);
plot(tau,ones(1,201))
hold on
for j=2:2:22
    plot(tau,framedata(:,j)')
    hold on
end
text(1.2,.25,['x/L = ' num2str(1)])
text(1.2,.66,['x/L = ' num2str(0)])
%text(.082,.42,['x/L = ' num2str(0.2)])
%text(0.09,0.4,['x/L = ' num2str(.2)])
xlabel('Dimensionless Time, \kappat/L^2')
ylabel('Dimensionless Temperature, T(x,t)/T_0')
text(.4,.9,'Press Enter to Continue')
pause

figure(3)%now do animation
xpatch=[0 0 .2 .3 .5 .6 .85 1 1 .82 .66 .43 .2 0];
ypatch=[-.1 1.1 1.15 1.05 1.13 1.1 1.14 1.1 -.1 -.03 -.14 -.1 -.16 -.1];
patch(xpatch,ypatch,'c')
axis([-.2 1.2 -.2 1.2])
hold on
text(-.15,.8,'T = T_0')
text(1.05,.8,'T = 0')
text(-.05,.4,'h_1');
text(1.02,.4,'h_2');
box on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Dimensionless Temperature, T(x,t)/T_0')
hold on
L=plot(xoverL, framedata(1,:),'k','EraseMode','xor');
hold on
texthandl=text(.4,1,'Press Enter to Animate');
pause
set(texthandl,'String','                       ');
for i=2:201
    set(L,'Ydata',framedata(i,:))
    pause(.02)
end
hold on
plot(xoverL,framedata(1,:),'k')
hold on

text(.5,.7,['h_1L/k = ' num2str(h1Loverk)]);
text(.5,.6,['h_2L/k = ' num2str(h2Loverk)]);
set(texthandl,'String','Press Enter to Continue');
pause

figure(4)%Color plot
axis([-.2 1.2 -.2 1.2 ]);
hold on
framedata(:,1)=ones(201,1);
colormap(cool);
zdata=zeros(21,23);
for j=1:21
    zdata(j,:)=framedata(1,:);
end
Q=pcolor(xoverL,linspace(0,1,21),zdata);
xlabel('Dimensionless Distance, x/L');
hold on
ylabel('Dimension y');
texthandl2=text(0,1.1,'Press Enter to Continue');
pause

set(texthandl2,'String','        ');
y=linspace(0,1,21);
for i=2:2:200
    for j=1:21
        zdata(j,:)=framedata(i,:);
    end
    Q=pcolor(xoverL,y,zdata);
    pause(.005)
end
set(texthandl2,'String','Press Enter to Continue');
pause
    
figure(5);clf;%3-D Plot of T(x,t) vs x and t
[X,Y]=meshgrid(xoverL,tau);
mesh(X,Y,framedata)
colormap(cool)
axis([-.1 1.1 0 2 0 1.1]);
hold on
text(0,.2,.6,'Press Enter to Continue');
text(0,0,-.2,'Location, x/L','Rotation',-34)
text(1.1,.4,-.2,'Dimensionless Time, \kappat/L^2','Rotation',12)
zlabel('Dimensionless Temperature, T(x,t)/T_0')
view(60,30)