% Conduction of heat in a sphere of radius R
% with exterior face (r=R) at temperature of T0
% and initial interior at zero temperature.
% The diffusivity is kappa=k/(rho*c)
% This is the well-known "Cooking a Turkey Problem" problem posed when
% investigating how long it will take to cook a turkey starting at
% temperature T0 and with a constant external temperature TH.
% 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);
roverR=linspace(0,1,101);roverR(1,1)=.0001;
nterm=20;% number of terms in the Fourier series solution.
for i=1:nterm
    phi(i,:)=sin(i*pi*roverR);
end
tau=linspace(0,.2,201);
T1=(ones(1,101));
F=ones(1,101)./roverR;
framedata=[];
for k=1:201%time loop
    sum=zeros(1,101);
    for i=1:nterm
        sum=sum+(((-1)^i)/i)*phi(i,:)*exp(-(i*i*pi*pi*tau(1,k)));
    end
    sum=sum*(2/pi).*F;
    T=T1+sum;
    framedata=[framedata;T];
end
framedata(1,:)=zeros(1,101);
framedata(1,101)=1;

figure(1);clf;% Plot T(x,t) vs x for various t
axis([-.2 1.2 -.2 1.2])
hold on
plot([0 0],[-.2,1.2]);
hold on
plot([1 1],[-.2,1.2]);
hold on
plot(roverR,framedata(1,:),'k')
hold on
xlabel('Dimensionless Radius, r/R')
hold on
ylabel('Dimensionless Temperature, [T(r,t)-T_0]/[T_H-T_0]')
for k=21:20:201% time loop
    plot(roverR,framedata(k,:),'k')
    hold on
end
text(.7,.05,['\kappat/R^2 = ' num2str(0)])
text(.65,.2,num2str(0.02))
text(.6,.3,num2str(0.04))
text(0.4,.85,num2str(0.2))
set(gca,'Box','on')
text(.1,1.13,'Press Enter to Continue')
pause
hold off

figure(2);clf;%Plot T(x,t) vs t for various x
axis([0 .2 -.2 1.2]);
hold on
plot(tau,ones(1,201))
hold on
for j=1:10:101% space loop
    plot(tau,framedata(:,j))
    hold on
end
text(.02,.95,['r/R = ' num2str(1)])
text(.04,.87,['r/R = ' num2str(0.9)])
text(.125,.35,['r/R = ' num2str(0)])
xlabel('Dimensionless Time, \kappat/R^2')
ylabel('Dimensionless Temperature, [T(r,t)-T_0]/[T_H-T_0]')
text(.02,1.1,'Press Enter to Continue')
set(gca,'Box','on')
pause
hold off

figure(3);clf;%now do animation
axis([-.2 1.2 -.2 1.2])
hold on
plot([0 0],[-.2,1.2]);
plot([1 1],[-.2,1.2]);
box on
xlabel('Dimensionless Radius, r/R')
hold on
ylabel('Dimensionless Temperature, [T(r,t)-T_0]/[T_H-T_0]')
hold on
L=plot(roverR, framedata(1,:),'k','EraseMode','xor');
hold on
texthandl=text(.4,1.1,'Press Enter to Animate');
pause
set(texthandl,'String','                       ');
for k=2:201% time loop
    set(L,'Ydata',framedata(k,:))
    pause(.05)
end
hold on
plot(roverR,framedata(1,:),'k')
hold on
set(texthandl,'String','Press Enter to Continue');
pause

figure(4);clf;%3-D Plot of T(r,t) vs x and t
%tau1=linspace(0,.2,101);
[X,Y]=meshgrid(roverR,tau);
mesh(X,Y,framedata)
colormap(cool)
text(.2,0,-.2,'Dimensionless Radius, r/R','Rotation',10)
text(0,.17,-.25,'Time, \kappat/R^2','Rotation',-34)
zlabel('[T(r,t) - T_0]/[T_H - T_0]')
view(330,30)
axis([0 1 0 .2 0 1])
hold off
text(.05,.2,.93,'Press Enter to Quit','Rotation',10)
pause