% Conduction of heat in a slab with both surfaces x=0 and x-L
% suddenly elevated fron zero to T0 degrees and slab initially
% at temperature zero.
% The diffusivity is kappa=k/(rho*c).
% 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,101);
sn=[];
nterm=40;
for i=1:nterm
    sn(i,:)=sin((2*i-1)*pi*xoverL);
end
tau=linspace(0,.2,101);
T1=(ones(1,101));
framedata=[];
for k=1:101
    sum=zeros(1,101);
    for i=1:nterm
        sum=sum-((4/((2*i-1)*pi))*exp(-(((2*i-1)*pi)^2)*tau(1,k))*sn(i,:));
    end
    T=T1+sum;
    framedata=[framedata;T];
end
framedata(1,:)=zeros(1,101);
framedata(:,1)=ones(101,1);
framedata(:,101)=ones(101,1);
figure(1);clf;
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:10:101
    plot(xoverL,framedata(k,:),'k')
    hold on
    axis([-.2 1.2 -.2 1.2])
end
text(.05,.05,['\kappat/L^2 = ' num2str(0)])
text(.16,.2,num2str(0.02))
text(.4,.85,num2str(0.2))
text(1.04,.8,'T = T_0')
text(-.15,.8,'T = T_0')
set(gca,'Box','on')
text(.1,1.13,'Press Enter to Continue')
pause

figure(2);clf;
plot(tau,ones(1,101))
hold on
for j=11:10:51
    plot(tau,framedata(:,j))
    hold on
end
text(.02,.97,['x/L = ' num2str(0)])
text(.03,.8,['x/L = ' num2str(0.1)])
text(.05,.6, '0.2')
text(0.06,0.5,'0.3')
text(.068,.44,'0.4')
text(.075,.38,'0.5')
xlabel('Dimensionless Time, \kappat/L^2')
ylabel('Dimensionless Temperature, T(x,t)/T_0')
text(.1,.2,'Press Enter to Continue')
pause

figure(3);clf;%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 = T_0')
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:101
    set(L,'Ydata',framedata(i,:))
    pause(.05)
end
hold on
plot(xoverL,framedata(1,:),'k')
hold on
set(texthandl,'String','Press Enter to Continue');
pause
figure(4);clf;%Do dynamic color plot
axis([-.2 1.2 -.2 1.2 ]);
hold on
framedata(:,1)=ones(101,1);
framedata(:,100)=ones(101,1);
colormap(cool);
zdata=zeros(101,101);
for j=1:101
    zdata(j,:)=framedata(1,:);
end
Q=pcolor(xoverL,xoverL,zdata);
xlabel('Dimensionless Distance, x/L');
hold on
ylabel('Dimension y');
texthandl2=text(0,1.1,'Press Enter to Continue');
pause
set(texthandl2,'String','        ');
for i=2:60
    for j=1:101
        zdata(j,:)=framedata(i,:);
    end
    Q=pcolor(xoverL,xoverL,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)
xlabel('Distance, x/L','Rotation',-34)
ylabel('Time, \kappat/L^2','Rotation',12)
zlabel('Dimensionless Temperature, T(x,t)/T_0')
text(0.8,.05,.2,'Press Enter to Continue','Rotation',12)
view(60,30)
axis([0 1 0 .2 0 1])
pause