%
%
% Ricardo Carretero 20/Oct/2003
% Integrates a reaction diffusion eq in 2D using finite differences
%
%        u_t = f(u) + D (u_xx+u_yy)
%
% Ex:  PDE_2D(Nx,Ny,dt,tf,Lx,Rx,Ly,Ry,disp,savesnaps)
%
% npts=N, dt:time step, tf:t_final, domain:[L,R]
%
%Ex:
%global D allt allx ally allu ran;
%Nx=20;Ny=20;dt=0.0005;tf=5;Lx=0;Rx=pi+1;Ly=0;Ry=pi+1;snaps=100;disp=100;D=1;
%PDE_2D(Nx,Ny,dt,tf,Lx,Rx,Ly,Ry,disp,snaps)
%

function [] = PDE_2D(Nx,Ny,dt,tf,Lx,Rx,Ly,Ry,exposh,saveframes)

 global allt allx ally allu;

 allx=[]; ally=[]; allt=[]; allu=[];

 ran=1;
 noise=0;

 %This next line is here so I do not forget to globalize 
 if(isempty(ran))
   fprintf('\n');
   error('Error: first do: global global dx dxx dy dyy D oos allt allx ally allu ran;')
 end

 if(Nx~=Ny)
   fprintf('\n');
   error('You need Nx=Ny, Nx not= Ny has not been implemented yet')
 end

 maxstep=tf/dt;
 numdispscreen=100;
 dispscreen= fix(maxstep/numdispscreen);
 if(dispscreen==0) dispscreen=1;end;

 if(exposh>0) stopdisp= fix(maxstep/exposh); 
 else stopdisp=9e+99;
 end

 if(saveframes>0)stopsave= fix(maxstep/saveframes);
 else stopsave=9e+99;
 end
 if(stopsave==0) stopsave=9e+99;end;

 % Spatial computational domain defined

 x=linspace(-Lx,Rx,Nx);
 y=linspace(-Ly,Ry,Ny);
 X=kron(x,ones(Ny,1));
 Y=kron(y',ones(1,Nx));
 dx=x(2)-x(1);
 dy=y(2)-y(1);
 dxx=1/dx^2;
 dyy=1/dy^2;
 oos=1/6;

 rand('state',0)       %this fixes the seed for the random generation
 strlength=0;
 ifig=5;

 allx = x;
 ally = y;
 allt = 0:stopsave*dt:tf;
 allt = allt;

 allu = zeros(saveframes+1,Ny,Nx);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% uini=      0.01*exp(-(X-(Rx-Lx-2)/2).^2).*(0.1*rand(Ny,Nx)+1+sin(20*pi*(X-Lx-0.1)/(Rx-Lx)));
% uini=uini.*0.01*exp(-(Y-(Ry-Ly-2)/2).^2).*(0.1*rand(Ny,Nx)+1+sin(20*pi*(Y-Ly-0.1)/(Ry-Ly)));

 uini = 0.01*(0.1*rand(Ny,Nx)+sin(3*pi*(X-Lx)/(Rx-Lx)));
 uini = uini .*sin(2*pi*(Y-Ly)/(Ry-Ly));
 uini=uini+0*0.001;
 uini(1,1:Nx)=0;
 uini(Ny,1:Nx)=0;
 uini(1:Ny,1)=0;
 uini(1:Ny,Nx)=0;

 save=0;
 t=0;
 u=uini;
 allu(1,:,:) = u;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 save=1;
 kk=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 fprintf('Lx \t\t= %f\n',Lx);
 fprintf('Rx \t\t= %f\n',Rx);
 fprintf('Ly \t\t= %f\n',Ly);
 fprintf('Ry \t\t= %f\n',Ry);
 fprintf('dt \t\t= %f\n',dt);
 fprintf('dx \t\t= %f\n',dx);
 fprintf('dy \t\t= %f\n',dy);
 fprintf('Nx \t\t= %f\n',Nx);
 fprintf('Ny \t\t= %f\n',Ny);
 fprintf('tf \t\t= %f\n',tf);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot initial data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 high = max(max(u));
 low = min(min(u));
 viewUtop=[0 90];
 viewU=[-38 30];
 
  figure(1+ifig); clf;
  set(gca,'FontSize',[16]);
  zl='u(x,y,t)';title='';
  mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
  plotopU='mag'; axU=[Lx Rx Ly Ry low high];
  plotU(x,y,u,axU,zl,title,mytext,plotopU);
  view(viewU);
 
  figure(2+ifig); clf;
  set(gca,'FontSize',[16]);
  zl='u(x,y,t)';title='';
  mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
  plotopU='mag'; axU=[Lx Rx Ly Ry low high];
  plotU(x,y,u,axU,zl,title,mytext,plotopU);
  view(viewU);
 
  figure(3+ifig); clf;
  set(gca,'FontSize',[16]);
  zltop='u(x,y,t)';title='';
  plotopUtop='height'; axUtop=[Lx Rx Ly Ry low-5*(high-low) high+5*(high-low)];
  axis(axUtop);
  plotU(x,y,u,[],zltop,title,mytext,plotopUtop);
  view(viewUtop);
  colorbar

 pause(0.001);
% fprintf('Press any key...');pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for kk=1:maxstep
  t=kk*dt;

   u = rk4_2d(u,Nx,Ny,dt,oos,dxx,dyy);
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% include noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(noise>0)
    u  = u+noise*(rand(Ny,Nx)-0.5);
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display time on screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if (fix(kk/dispscreen)==kk/dispscreen)
    fprintf('%5.2f ',kk*dt);
    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+ifig); clf;
       set(gca,'FontSize',[16]);
       mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
       plotU(x,y,u,[],zl,title,mytext,plotopU);
       view(viewU);
       
       figure(2+ifig); clf;
       set(gca,'FontSize',[16]);
       mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
       plotU(x,y,u,axU,zl,title,mytext,plotopU);
       view(viewU);

       figure(3+ifig); clf;
       set(gca,'FontSize',[16]);
       plotU(x,y,u,[],zltop,title,mytext,plotopUtop);
       view(viewUtop);
%       colorbar
       
       pause(0.001);
    end

 end
figure(3+ifig);colorbar
 fprintf('\n');
 return;



