% Conduction of heat in a sphere of radius R.
% The exterior face (r=R) has a sinusoidal variation in temperature
% T0*cos(omega*t).  The steady-state sinusoidal temperature interior
% to the sphere is illustrated.
% The diffusivity is kappa=k/(rho*c) and the input quantity is the
% dimensionless frequency (R^2)*omega/kappa.
% 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);
freqpar=-1;
while freqpar<=0
    freqpar=input('The dimensionless frequency is (R^2)*omega/kappa = ');
end
alpha=sqrt(freqpar/2)*(1+1i);
roverR=linspace(0,1,101);
roverR(1,1)=.0001;
omegat=linspace(0,6*pi,301);
T1=(ones(1,101));
F=(ones(1,101)./roverR)/(exp(alpha)-exp(-alpha));
G=F.*(exp(alpha*roverR)-exp(-alpha*roverR));
phase=angle(G);
mag=abs(G);
framedata=[];
for t=1:301
    framedata(t,:)=mag.*cos(omegat(1,t)*ones(1,101)+phase);
end

figure(1);clf;% Plot T(r,t) vs r for various t
axis([0 1 -1.2 1.2])
hold on
xlabel('Dimensionless Radius, r/R')
hold on
ylabel('Dimensionless Temperature, T(r,t)/T_0')
for t=1:10:101% time loop
    plot(roverR,framedata(t,:),'k')
    hold on
    if t<=21
        xline=[.9-.01*(t-1) .8-.01*(t-1)];
        yline=[framedata(t,91-t+1) framedata(t,91-t+1)+.1];
        plot(xline,yline,'k')
        text1=num2str(.02*(t-1));
        text(.65-.01*(t-1),framedata(t,91-t+1)+.1,['\omegat = ' text1 '\pi'])
    end
end
text(.2, -1.1,['Dimensionless Frequency = R^2\omega/\kappa = ' num2str(freqpar)]);
set(gca,'Box','on')
text(.2,1.1,'Press Enter to Continue')
hold off
pause

figure(2);clf;%Plot T(x,t) vs t for various x
axis([0 6*pi -1.2 1.2]);
hold on
for j=1:10:101% space loop
    plot(omegat,framedata(:,j))
    hold on
end
for j=101:-10:71
    plot([omegat(1,201) 15.8], [framedata(201,j) framedata(201,j)+.1],'k')
    text(15.9,framedata(201,j)+.1,['r/R = ' num2str(roverR(1,j))])
end
text(1,-1.1,['Dimensionless Frequency = R^2\omega/\kappa = ' num2str(freqpar)])
xlabel('Dimensionless Time, \omegat')
ylabel('Dimensionless Temperature, T(r,t)/T_0')
text(1,1.1,'Press Enter to Continue')
set(gca,'Box','on')
hold off
pause

figure(3);clf;%now do animation
axis([0 1 -1.2 1.2])
hold on
box on
xlabel('Dimensionless Radius, r/R')
hold on
ylabel('Dimensionless Temperature, T(r,t)/T_0')
hold on
L=plot(roverR, framedata(1,:),'k','EraseMode','xor');
text(.2,-1.1,['Dimensionless Frequency = R^2\omega/\kappa = ' num2str(freqpar)])
texthandl=text(.2,1.1,'Press Enter to Animate');
pause
set(texthandl,'String','                       ');
for k=2:301% time loop
    set(L,'Ydata',framedata(k,:))
    pause(.05)
end
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,omegat(1,1:101));
mesh(X,Y,framedata(1:101,:))
colormap(cool)
text(.2,0,-1.5,'Dimensionless Radius, r/R','Rotation',12)
text(0,6,-1.45,'Dimensionless Time, \omegat','Rotation',-33)
zlabel('Dimensionless Temperature, T(r,t)/T_0')
view(330,30)
axis([0 1 0 2*pi -1.1 1.1])