% Conduction of heat in an infinite cylinder of radius R
% with exterior face (r=R) at T=0 and and initial interior at zero
% temperature. Heat is generated at a uniform rate q per unit volume as 
% would be the case with ohmic heating with a uniform current density.
% 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);
% eigenvalues--zeros of the J0 Bessel function
lambdaR=[2.404826 5.520078 8.653728 11.791534 14.9309177 18.0761064];
lambdaR=[lambdaR 21.211637 24.352472 27.493479 30.634606];
roverR=linspace(0,1,101);
nterm=10;
for i=1:nterm
    J0(i,:)=besselj(0,lambdaR(1,i)*roverR);
    J1(1,i)=besselj(1,lambdaR(1,i));
end
tau=linspace(0,.8,201);
%T1=(ones(1,101));
framedata=[];
for k=1:201
    sum=zeros(1,101);
    for i=1:nterm
        sum=sum+(2/(((lambdaR(1,i))^3)*J1(1,i)))*(1-exp(-((lambdaR(1,i)^2)*tau(1,k))))*J0(i,:);
    end
    framedata=[framedata;sum];
end
framedata(1,:)=zeros(1,101);
%framedata(1,101)=1;

figure(1);clf;% Plot T(x,t) vs x for various t
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)k/qR^2')
for k=21:20:201
    plot(roverR,framedata(k,:),'k')
    hold on
    axis([-.2 1.2 -.1 .3])
end
text(.2,-.02,['\kappat/R^2 = ' num2str(0)])
text(.2,.06,num2str(0.08))
text(.2,.13,num2str(0.16))
text(0.2,.25,num2str(0.8))
set(gca,'Box','on')
text(.6,.26,'Press Enter to Continue')
pause

figure(2);clf;%Plot T(x,t) vs t for various x
axis([0 .3 -.1 .3]);
hold on
plot(tau,ones(1,201))
hold on
for j=11:10:101
    plot(tau,framedata(:,j))
    hold on
end
text(.23,-.02,['r/R = ' num2str(1)])
text(.2,.028,['r/R = ' num2str(0.9)])
text(.16,.17,['r/R = ' num2str(0)])
xlabel('Dimensionless Time, \kappat/R^2')
ylabel('Dimensionless Temperature, T(r,t)k/qR^2')
text(.02,.25,'Press Enter to Continue')
pause

figure(3);clf;%now do animation
axis([-.2 1.2 -.1 .3])
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)k/qR^2')
hold on
L=plot(roverR, framedata(1,:),'k','EraseMode','xor');
hold on
texthandl=text(.4,.25,'Press Enter to Animate');
pause
set(texthandl,'String','                       ');
for i=2:201
    set(L,'Ydata',framedata(i,:))
    pause(.03)
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,tau1);
mesh(X,Y,framedata(1:101,:))
colormap(cool)
text(.2,0,-.06,'Dimensionless Radius, r/R','Rotation',12)
text(0,.17,-.06,'Time, \kappat/R^2','Rotation',-34)
zlabel('Dimensionless Temperature, T(r,t)k/qR^2')
view(330,30)
axis([0 1 0 .2 0 .3])