%One dimensional conduction of heat in a slab with left face at T0 and 
% right face at zero and slab initially at zero temperature.
% 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(i*pi*xoverL);
end
tau=linspace(0,.2,101);
T1=(ones(1,101)-xoverL);
framedata=[];
for k=1:101
    sum=zeros(1,101);
    for i=1:nterm
        sum=sum-((2/((i)*pi))*exp(-(((i)*pi)^2)*tau(1,k))*sn(i,:));
    end
    T=T1+sum;
    framedata=[framedata;T];
end
framedata(1,:)=zeros(1,101);

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: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,.54,num2str(0.2))
text(1.04,.8,'T = 0')
text(-.15,.8,'T = T_0')
set(gca,'Box','on')
text(.1,1.13,'Press Enter to Continue')
pause

figure(2)%Plot T(x,t) vs t for various x
plot(tau,ones(1,101))
hold on
for j=11:10:101
    plot(tau,framedata(:,j))
    hold on
end
text(.02,.97,['x/L = ' num2str(0)])
text(.03,.8,['x/L = ' num2str(0.1)])
text(.1,-.08,['x/L = ' num2str(1)])
text(0.04,0.6,['x/L = ' num2str(.2)])
xlabel('Dimensionless Time, \kappat/L^2')
ylabel('Dimensionless Temperature, T(x,t)/T_0')
text(.02,1.1,'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')
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)%Color plot
axis([-.2 1.2 -.2 1.2 ]);
hold on
framedata(:,1)=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=3:2:61
    for j=1:101
        zdata(j,:)=framedata(i,:);
    end
    Q=pcolor(xoverL,xoverL,zdata);
    pause(.01)
end
set(texthandl2,'String','Press Enter to Continue');
pause
    
figure(5)%3-D Plot of T(x,t) vs x and t
tau1=linspace(0,.1,51);
[X,Y]=meshgrid(xoverL,tau1);
mesh(X,Y,framedata(1:51,:))
colormap(cool)
text(.2,0,-.2,'Distance, x/L','Rotation',-34)
text(1,.05,-.2,'Time, \kappat/L^2','Rotation',12)
zlabel('Dimensionless Temperature, T(x,t)/T_0')
view(60,30)
axis([0 1 0 .1 0 1])