%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function workshop07
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2d diffusion
% Eli Tziperman, APM115

close all; clear all
cases=[2,3,1];
for ii=1:length(cases)
  the_case=cases(ii);
  [X,Y,framedata,NT]=solve_diffusion_eqn(the_case);
  plot_results(X,Y,framedata,NT);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X,Y,framedata,NT]=solve_diffusion_eqn(the_case)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=39;% number of internal nodes
kappa=0.3; % diffusion coeff
dt=.01; % time step
dx=1/(N+1); % grid spacing
A=kappa/dx;
B=dt/dx;
NT=3000; % number of time steps
[X,Y]=meshgrid(0:1:N,0:1:N); % grid for plotting
X=X*dx; Y=Y*dx;
T=initial_conditions(the_case,N);
diffusion=zeros(N+1,N+1);
Flux_x_in=zeros(N+1,N+1); Flux_x_out=zeros(N+1,N+1);
Flux_y_in=zeros(N+1,N+1); Flux_y_out=zeros(N+1,N+1);
% array with plot data:
framedata=zeros(NT,N+1,N+1);

% integrate equation in time:
for n=1:NT
  % Interior points are 2:N, all points are 1:N+1.  
  
  % If boundary conditions are of fixed temperature, need to calculate
  % the diffusion term only in interior points.  If boundary
  % conditions are of no flux, the term needs to be evaluated at all
  % grid points, such that the flux_in and flux_out terms are zero at
  % the appropriate bundary grid points.
  
  % diffusion term written in flux for that is easy to change to deal
  % with flux boundary conditions: flux=-kappa*dT/dx
  Flux_x_in(2:N+1,1:N+1)= -A*(T(2:N+1,1:N+1) - T(1:N,1:N+1));
  % (in above line, note that Flux_x_in(i=1,j)=0 already from
  % initialization above, and similarly for other boundary terms that
  % should vanish for zero flux boundary conditions)
  Flux_x_out(1:N,1:N+1)=-A*(T(2:N+1,1:N+1)-T(1:N,1:N+1));
  Flux_y_in(1:N+1,2:N+1)= -A*(T(1:N+1,2:N+1) - T(1:N+1,1:N));
  Flux_y_out(1:N+1,1:N)=-A*(T(1:N+1,2:N+1)-T(1:N+1,1:N));

  diffusion(1:N+1,1:N+1)=B*(Flux_x_in(1:N+1,1:N+1)-Flux_x_out(1:N+1,1:N+1) ...
                       +Flux_y_in(1:N+1,1:N+1)-Flux_y_out(1:N+1,1:N+1));

  %% time step temperature in time:
  T=T+dt*(diffusion+source(the_case,n,N));
  framedata(n,:,:)=T(:,:);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function T=initial_conditions(the_case,N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initial conditions:
T=zeros(N+1,N+1); T(2:N,2:N)=ones(N-1,N-1); 
if the_case==1
  T=0*ones(N+1,N+1); 
  T(2:N,2:N)=zeros(N-1,N-1); 
elseif the_case==2
  N2=round(N/2);
  T=zeros(N+1,N+1); 
  T(N2-2:N2+2,2:N)=T(N2-2:N2+2,2:N)*0+1;
  T(2:N,N2-2:N2+2)=T(2:N,N2-2:N2+2)*0+1;
elseif the_case==3
  T=zeros(N+1,N+1); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function my_source=source(the_case,n,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the source used in the diffusion equation
%% moving heat source that is rotating whithin the domain:
my_source=zeros(N+1,N+1); 
N2=round(N/2);dk=3;dl=3;
if the_case==1
  %% some parameters controlling heat source structure:
  omega=8*2*pi/1000; % speed of rotation in space of heat source
  kc=N/2+N/4*cos(omega*n);lc=N/2+N/4*sin(omega*n);
  for k=1:N-2; 
    for l=1:N-2;
      my_source(k,l)=exp(-(k-kc)^2/dk^2-(l-lc)^2/dl^2);
    end
  end
elseif the_case==2
  my_source=0*my_source;
elseif the_case==3
  kc=N/2;lc=N/2;
  for k=1:N-2; 
    for l=1:N-2;
      my_source(k,l)=exp(-(k-kc)^2/dk^2-(l-lc)^2/dl^2);
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plot_results(X,Y,framedata,NT)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot 3d animation of solution
figure(1);clf
Q=surf(X,Y,squeeze(framedata(1,:,:)));
axis([0 1 0 1 0 1]); 
shading interp;
xlabel('x');
ylabel('y');
pause

for n=2:2:NT
  Q=surf(X,Y,squeeze(framedata(n,:,:)));
  axis([0 1 0 1 0 1]); 
  shading interp;
  pause(.005)
end
disp('enter to continue');
pause
