% Conduction of heat in a slab, infinie in extent in y and z directions
% with left face at x=0 and right face x=L insulated. Initially left half 
% (0<x,L/2) of slab is initially at T0 and right half (L/2<x<L) 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);
xoverL=linspace(0,1,101);
cs=[];
a=[];
nterm=40;
for i=1:nterm;
    cs(i,:)=cos((2*i-1)*pi*xoverL);
    a(1,i)=(2*((-1)^(i+1)))/((2*i-1)*pi);
end
tau=linspace(0,.2,201);
framedata=[];
for k=1:201;%time loop
    sum=zeros(1,101);
    for i=1:nterm;%series loop
        sum=sum+a(1,i)*exp(-(((2*i-1)*pi)^2)*tau(1,k))*cs(i,:);
    end
    T=.5*ones(1,101)+sum;
    framedata=[framedata;T];
end
framedata(1,:)=[ones(1,50) .5 zeros(1,50)];
    
figure(1);clf;% Plot T(x,t) vs x for various t
axis([-.2 1.2 -.3 1.3]);
hold on
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.3*ones(1,11)+(y1) y1(1,1)];
patch(xpatch1,ypatch1,'c')
hold on
xinsul=[-.001 -.001 -.15 -.1 -.17 -.08 -.13 -.001];
yinsul=[-.10 1.22 1.0 .8 .68 .52 .18 -.1];
patch(xinsul, yinsul, 'r');
hold on
yinsul(1,2)=1.2;yinsul(1,1)=-.08;yinsul(1,8)=-.08;
patch(1.002*ones(1,8)-xinsul,yinsul,'r')
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:201
    plot(xoverL,framedata(k,:),'k')
    hold on
end
text(.52,.94,['\kappat/L^2 = ' num2str(0)])
text(.37,.92,num2str(0.01))
text(.15,.5,num2str(0.2))
text(1.04,.8,'k = 0')
text(-.15,.8,'k = 0')
set(gca,'Box','on')
text(.1,1.13,'Press Enter to Continue')
pause

figure(2);clf;%Plot T(x,t) vs t for various x
axis([0 .2 -.1 1.1])
hold on
plot(tau,framedata(:,1))
hold on
for j=11:10:101
    plot(tau,framedata(:,j))
    hold on
end
text(.04,.97,['x/L = ' num2str(0)])
text(.045,.8, num2str(0.2))
text(.05,.65, num2str(.4))
text(0.055,.4, num2str(.6))
text(.1,.25,['x/L = ',num2str(1)])
xlabel('Dimensionless Time, \kappat/L^2')
ylabel('Dimensionless Temperature, T(x,t)/T_0')
text(.1,1,'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(xpatch1,ypatch1,'c')
axis([-.2 1.2 -.3 1.3])
hold on
yinsul(1,2)=1.22;yinsul(1,1)=-.1;yinsul(1,8)=-.1;
patch(xinsul, yinsul, 'r');
hold on
yinsul(1,2)=1.2;yinsul(1,1)=-.08;yinsul(1,8)=-.08;
patch(1.002*ones(1,8)-xinsul,yinsul,'r')
hold on
text(-.15,.8,'k = 0')
text(1.05,.8,'k = 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.05,'Press Enter to Animate');
pause
set(texthandl,'String','                       ');
for i=2:201
    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;%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
yaxis=linspace(0,1,101);
Q=pcolor(xoverL,yaxis,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=5:4:201
    for j=1:101
        zdata(j,:)=framedata(i,:);
    end
    Q=pcolor(xoverL,yaxis,zdata);
    pause(.01)
end
set(texthandl2,'String','Press Enter to Continue');
pause
    
figure(5);clf;%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])