%
% Ricardo Carretero 7/Dec/2005
% Integrates the following PDE using finite differences in time and space
%
%   u_t = f0(u,x) + f1(u,x) u_x + f2(u,x) u_xx + f3(u,x) u_xxx + f4(u,x) u_xxxx
%
% Ex:  PDE_integrator(N,dt,tf,R,disp,savesnaps)
%
% npts=N, dt:time step, tf:t_final, domain:[0,R]
%
%Ex: 
%global oodx oodxx oodxxx oodxxxx BCs;
%BCs=0		%% 0:fixed, 1:periodic. 2:no flux
%N=100;dt=0.005;tf=30;R=10;snaps=100;disp=100;
%[X,T,U]=PDE_integrator(N,dt,tf,R,disp,snaps);
%wave_plotsurf(T,X,U',R,1,snaps+1,'no')
%

function [X,T,U] = PDE_integrator(N,dt,tf,R,exposh,saveframes)

 global oodx oodxx oodxxx oodxxxx BCs;
 
 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

 L=0;
 dx=(R-L)/N;
 x=L+(0:N-1)*dx;
 oodx = 1/dx;
 oodxx = 1/dx^2;
 oodxxx = 1/dx^3;
 oodxxxx = 1/dx^4;
 oo3=1/3;
 oo4=1/4;
 oo5=1/5;
 oo6=1/6;
 
 rand('state',0)       %this fixes the seed for the random generation
 
 strlength=0;

 T = 0:stopsave*dt:tf;
 U = zeros(saveframes+1,N);
 X = 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=10*(0.1*rand(1,N)+cos(1*pi*(x-L)/(R-L)));
 uini=100+x*0; uini(1)=0; uini(N)=0;
 save=0;
 t=0;
 u=uini;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot IC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure(1); clf;
 set(gca,'FontSize',[16]);
 plot(x,u,'LineWidth',3); 
 title('u(x,0)');
 ax0=axis;
 if(ax0(3)>0)
  ax0(3)=0;
 end
 axis(ax0)
 fprintf('Press any key to continue\n');
 pause
 
 save=save+1;
 U(save,:) = u;
     
 fprintf('\nt=%5.2f ',t);
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for kk=1:maxstep
    
  t=t+dt;
  uN = u+dt*deriv_u(u,x,N);

switch BCs	% 0:fixed, 1:periodic. 2:no flux
  case 0			%% Fixed BCs
   uN(1)=uini(1);
   uN(N)=uini(N);
  case 2			%% no flux
   uN(1)=uN(2);
   uN(N)=uN(N-1);
  end

  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;
    U(save,:) = u;
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(fix(kk/stopdisp)==kk/stopdisp&kk>1)
     figure(1); clf;
     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);
     text(R-0.2*(R-L),ax0(4)*0.9,mytext); axis(ax0);
     pause(0.0001)
     %pause
   end

 end

 fprintf('\n');
 return;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FINITE DIFFERENCE SCHEME FROM 
% http://mathworld.wolfram.com/FiniteDifference.html
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function du = deriv_u(u,x,N)
 global oodx oodxx oodxxx oodxxxx BCs;


switch BCs	% 0:fixed, 1:periodic. 2:no flux
case 0			%% Fixed BCs
 up = [ 0 u(1:N-1) ];
 um = [ u(2:N) 0 ];
 upp = [ 0 0 u(1:N-2) ];
 umm = [ u(3:N) 0 0 ];
case 1			% Periodic BCs:
 up = [ u(N) u(1:N-1) ];
 um = [ u(2:N) u(1) ];
 upp = [ u(N-1:N) u(1:N-2) ];
 umm = [ u(3:N) u(1:2) ];
case 2			% no flux
 up = [ 0 u(1:N-2) 0 ];
 um = [0  u(3:N) 0 ];
 upp = [ 0 0 u(1:N-2) ];
 umm = [ u(3:N) 0 0 ];
 u(1)=0;
 u(N)=0;
end

ff0=f0(u,x);
ff1=f1(u,x);
ff2=f2(u,x);
ff3=f3(u,x);
ff4=f4(u,x);

 du = ff0 + ...
      ff1.*oodx.*(up-um)/2 + ...
      ff2.*oodxx.*(up-2*u+um) + ...
      ff3.*oodxxx.*(upp-2*up+2*um-umm)/2 + ...
      ff4.*oodxxxx.*(upp-4*up+6*u-4*um+umm);

if(BCs==0)
% Fixed BCs
 du(1)=0;
 du(N)=0;
end
 
 return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fu0 = f0(u,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global oodx oodxx oodxxx oodxxxx BCs;
 %fu0 = u.*(1-u);
 fu0 = 0;
 return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fu1 = f1(u,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global oodx oodxx oodxxx oodxxxx BCs;
 fu1 = 0;
 return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fu2 = f2(u,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global oodx oodxx oodxxx oodxxxx BCs;
 fu2 = 1;
 return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fu3 = f3(u,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global oodx oodxx oodxxx oodxxxx BCs;
 fu3 = 0;
 return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fu4 = f4(u,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global oodx oodxx oodxxx oodxxxx BCs;
 fu4 = 0;
 return