%
% Ricardo Carretero 20/Oct/2003
% Integrates a reaction diffusion eq in 1D using finite differences
%
%        u_t = f(u) + D u_xx
%
% Ex:  diffusion_PDE(N,dt,tf,L,R,disp,savesnaps)
%
% npts=N, dt:time step, tf:t_final, domain:[L,R]
%
%Ex: 
%global dx dxx D oos allt allx allu;
%N=100;dt=0.005;tf=5;L=0;R=10;snaps=100;disp=100;D=1;
%diffusion_PDE(N,dt,tf,L,R,disp,snaps)
%wave_plotsurf(allt,allx,allu',R,1,N+1,'no')
%

function [] = diffusion_PDE(N,dt,tf,L,R,exposh,saveframes)

 global allt allx allu;
 global dx dxx oos;
 
 ran=1;
 noise=0;

 maxstep=tf/dt;
 dispscreen= fix(maxstep/10); 

 if(exposh>0) stopdisp= fix(maxstep/exposh); 
 else stopdisp=9e+99;
 end
 if(saveframes>0)stopsave= fix(maxstep/saveframes);
 else stopsave=9e+99;
 end
 
 % Spatial computational domain defined

 dx=(R-L)/N;
 x=L+(0:N-1)*dx;
 dxx = 1/dx^2;
 oos=1/6;
 
 rand('state',0)       %this fixes the seed for the random generation
 
 strlength=0;

 allt = 0:stopsave*dt:tf;
 allu = zeros(saveframes+1,N);
 allx = x;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 fprintf('\nL \t\t= %f\n',L);
 fprintf('R \t\t= %f\n',R);
 fprintf('dt \t\t= %f\n',dt);
 fprintf('dx \t\t= %f\n',dx);
 fprintf('N \t\t= %f\n',N);
 fprintf('tf \t\t= %f\n',tf);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %uini=0.01*exp(-(x-(R-L)/2).^2).*(0.1*rand(1,N)+1+sin(20*pi*(x-L)/(R-L)));
 %uini=0.01*(1.0*rand(1,N)+sin(2*pi*(x-L)/(R-L)));
 uini=0.1*sin(1*pi*(x-L)/(R-L))+0.1*sin(2*pi*(x-L)/(R-L))+ ...
      0.3*sin(3*pi*(x-L)/(R-L))+0.3*sin(4*pi*(x-L)/(R-L))+ ...
      0.5*sin(5*pi*(x-L)/(R-L))+0.5*sin(6*pi*(x-L)/(R-L))+ ...
      0.9*sin(7*pi*(x-L)/(R-L))+0.9*sin(8*pi*(x-L)/(R-L));
 save=0;
 t=0;
 u=uini;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot IC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure(1); clf;
 set(gca,'FontSize',[16]);
 plot(x,u,'LineWidth',3); 
 title('u(x,0)');
 axis tight;
 ax=axis;
 kk=0;
 subplot(2,1,1)
 set(gca,'FontSize',[16]);
 mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
 plot(x,u,'LineWidth',3); 
 text(R-0.2*(R-L),max(u)*0.9,mytext);
%     axis(ax);
 subplot(2,1,2)
 set(gca,'FontSize',[16]);
 mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
 plot(x,u,'LineWidth',3); 
 text(R-0.2*(R-L),ax(4)*0.9,mytext);
 axis(ax);
 pause(0.0001)
 fprintf('Press any key\n');
 pause
 
 save=save+1;
 allu(save,:) = u;
     
 fprintf('\nt=%5.2f ',t);
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for kk=1:maxstep
    
  t=t+dt;
  uN = u+dt*deriv_u(u,N);

  %Fixed BCs
  uN(1)=0;
  uN(N)=0;

  u = uN;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% include noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(noise>0)
    u  = u+noise*(rand(N)-0.5);
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display time on screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if (fix(kk/dispscreen)==kk/dispscreen)
    fprintf('%5.2f ',t);
    strlength=strlength+1;
    if strlength>=10
     strlength=0;
     fprintf('\n');
    end
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% store solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if (fix(kk/stopsave)==kk/stopsave)
    save=save+1;
    allu(save,:) = u;
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(fix(kk/stopdisp)==kk/stopdisp&kk>1)
     figure(1); clf;
     subplot(2,1,1)
     set(gca,'FontSize',[16]);
     mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
     plot(x,u,'LineWidth',3); 
     text(R-0.2*(R-L),max(u)*0.9,mytext);
%     axis(ax);
     subplot(2,1,2)
     set(gca,'FontSize',[16]);
     mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
     plot(x,u,'LineWidth',3); 
     text(R-0.2*(R-L),ax(4)*0.9,mytext);
     axis(ax);
     drawnow
     pause(0.0001)
   end

 end

 fprintf('\n');
 return;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function rkZn = ode_rk(rkZ)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% global dx dxx oos;
%
% k1 = dt*deriv_u(rkZ       );
% k2 = dt*deriv_u(rkZ+0.5*k1);
% k3 = dt*deriv_u(rkZ+0.5*k2);
% k4 = dt*deriv_u(rkZ+    k3);
% 
% rkZn = rkZ + oos*(k1 + 2*k2 + 2*k3 + k4);
% 
% return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function du = deriv_u(u,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global dx dxx D oos;
 
%% Periodic BCs:
% up = [ u(1) u(1:N-1) ];
% um = [ u(2:N) u(N) ];

% Fixed BCs
 up = [ 0 u(1:N-1) ];
 um = [ u(2:N) 0 ];

 du = f(u) + D*dxx*(up-2*u+um);
 
 return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fu = f(u)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global dx dxx D oos;

 fu = 0*1*u;

