%
%
% Ricardo Carretero 20/Oct/2003
% Integrates a reaction diffusion eq in 2D with 2 populations using finite differences
%
%        u_t = f(u,v) + D1 (u_xx+u_yy)
%        v_t = g(u,v) + D2 (v_xx+v_yy)
%
% Ex:  PDE_2D_2pop(Nx,Ny,dt,tf,Lx,Rx,Ly,Ry,disp,savesnaps)
%
% npts=N, dt:time step, tf:t_final, domain:[L,R]
%
%Ex:
%global D1 D2 allt allx ally allu allv ran;
%global aa bb cc dd hh kk
%aa=1;bb=1;cc=0.3;dd=1;hh=0.4;kk=2;D1=1;D2=1;
%Nx=40;Ny=40;dt=0.1;tf=1000;Lx=0;Rx=500;Ly=0;Ry=500;snaps=100;disp=100;
%PDE_2D_2pop(Nx,Ny,dt,tf,Lx,Rx,Ly,Ry,disp,snaps)
%film2d
%

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

 global allt allx ally allu allv;
 global D1 D2 allt allx ally allu allv ran;
 global aa bb cc dd hh kk

 format compact
 allx=[]; ally=[]; allt=[]; allu=[]; allv=[];

 ran=1;
 noise=0*1e-5;

 %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 D1 D2 oos allt allx ally allu allv 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/dy^2;
 dyy=1/dx^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);
 allv = zeros(saveframes+1,Ny,Nx);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %fixed point

% u0 = dd/cc;
% v0 = aa/bb;
% u0 = cc/(1-cc);
% v0 = (aa-aa*cc-cc*bb)/((cc-1)^2);
 u0 = cc*hh/(kk-cc);
 v0 = hh*kk*(aa*kk-aa*cc-cc*bb*hh)/((cc-kk)^2);

%check if fixed point is fixed:
reactionLV_u(u0,v0)
reactionLV_v(u0,v0)

% exite a parituclar mode (n,m):
 n=0;m=0;
 randu=1*1*u0;
 randv=1*1*v0;

 uini = u0;
 uini = uini+randu*(rand(Ny,Nx)-0.5);
 uini = uini+0*(cos(n*pi*(X-Lx)/(Rx-Lx)).*cos(m*pi*(Y-Ly)/(Ry-Ly)));

 vini = v0;
 vini = vini+randv*(rand(Ny,Nx)-0.5);
 vini = vini+1*(cos(n*pi*(X-Lx)/(Rx-Lx)).*cos(m*pi*(Y-Ly)/(Ry-Ly)));

 save=0;
 t=0;

 u=uini;
 v=vini;

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);
 fprintf('aa \t\t= %f\n',aa);
 fprintf('bb \t\t= %f\n',bb);
 fprintf('cc \t\t= %f\n',cc);
 fprintf('dd \t\t= %f\n',dd);
 fprintf('hh \t\t= %f\n',hh);
 fprintf('kk \t\t= %f\n',kk);
 fprintf('D1 \t\t= %f\n',D1);
 fprintf('D2 \t\t= %f\n',D2);

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

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

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

   [u,v] = rk4_2d_periodic(u,v,Nx,Ny,dt,oos,dxx,dyy);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% include noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(noise>0)
    u  = u+noise*(rand(Ny,Nx)-0.5);
    v  = v+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;
    allv(save,:,:) = v;
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if(fix(kk/stopdisp)==kk/stopdisp&kk>1)

       figure(1+ifig); clf;
       subplot(1,2,1)
       set(gca,'FontSize',[16]);
       mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
       plotU(x,y,u,[],zl,title,mytext,plotopU);
       view(viewU);
       subplot(1,2,2)
       set(gca,'FontSize',[16]);
       mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
       plotU(x,y,v,[],zl,title,mytext,plotopU);
       view(viewU);
       
       figure(3+ifig); clf;
       subplot(1,2,1)
       set(gca,'FontSize',[16]);
       plotU(x,y,u,[],zltop,title,mytext,plotopUtop);
       view(viewUtop);
%       colorbar
       axis tight;
       subplot(1,2,2)
       set(gca,'FontSize',[16]);
       plotU(x,y,v,[],zltop,title,mytext,plotopUtop);
       view(viewUtop);
       axis tight;
%       colorbar

       pause(0.001);
    end

 end
figure(3+ifig);subplot(1,2,1);colorbar
figure(3+ifig);subplot(1,2,2);colorbar
 fprintf('\n');
 return;



