%
%
% 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_conic(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 gg alpha rho K
%gg=1;alpha=1.5;K=0.1;rho=18.5;aa=92;bb=64;D1=1;D2=10;
%Nx=50;Ny=20;dt=0.002;tf=150;Lx=0;Rx=75;Ly=0;Ry=4;snaps=100;disp=100;
%PDE_conic(Nx,Ny,dt,tf,Lx,Rx,Ly,Ry,disp,snaps)
%

function [] = PDE_conic(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 gg alpha rho K

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

 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 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

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


 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 u0 = 9.933832844;
 v0 = 9.289221896;

 circumference=Ry-Ly;
 radius=circumference/(2*pi);
 theta=linspace(0,2*pi,Ny-1)'+pi/2;

% exite a parituclar mode (n,m):
 n=0;m=9*2;
 randnoise=1*0.25;

 uini = u0+1*0.01*(randnoise*(rand(Ny,Nx)-0.5) + ...
        0.1*(cos(m*pi*(X-Lx)/(Rx-Lx)) .* ...
           cos(n*(Y-Ly)/radius))) ;

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

 save=0;
 t=0;

% fix BCs
 uini(1,1:Nx)=u0;
 uini(Ny,1:Nx)=u0;
 uini(1:Ny,1)=u0;
 uini(1:Ny,Nx)=u0;
 vini(1,1:Nx)=v0;
 vini(Ny,1:Nx)=v0;
 vini(1:Ny,1)=v0;
 vini(1:Ny,Nx)=v0;

% no flux BCs
 uini(1,1:Nx) = uini(2,1:Nx);
 uini(Ny,1:Nx)= uini(Ny-1,1:Nx);
 uini(1:Ny,1) = uini(1:Ny,2);
 uini(1:Ny,Nx)= uini(1:Ny,Nx-1);
 vini(1,1:Nx) = vini(2,1:Nx);
 vini(Ny,1:Nx)= vini(Ny-1,1:Nx);
 vini(1:Ny,1) = vini(1:Ny,2);
 vini(1:Ny,Nx)= vini(1:Ny,Nx-1);

% no flux + predioic BCs
 uini(1,1:Nx) = uini(Ny,1:Nx);
 uini(Ny,1:Nx)= uini(Ny-1,1:Nx);
 uini(1:Ny,1) = uini(1:Ny,2);
 uini(1:Ny,Nx)= uini(1:Ny,Nx-1);
 vini(1,1:Nx) = vini(Ny,1:Nx);
 vini(Ny,1:Nx)= vini(Ny-1,1:Nx);
 vini(1:Ny,1) = vini(1:Ny,2);
 vini(1:Ny,Nx)= vini(1:Ny,Nx-1);

 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('Rymin \t\t= %f\n',Ry);
 fprintf('Rymax \t\t= %f\n',Ry*(1+slope));
 fprintf('dt \t\t= %f\n',dt);
 fprintf('dx \t\t= %f\n',dx);
 fprintf('dy~ \t\t= %f\n',sum(dy)/length(dy));
 fprintf('Nx \t\t= %f\n',Nx);
 fprintf('Ny \t\t= %f\n',Ny);
 fprintf('tf \t\t= %f\n',tf);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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];
 tailview=[56 26];
 mycolor=1-[zeros(1.5*32,3);gray(32)];
 
  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*(1+slope) 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*(1+slope) lowv highv];
  plotU(X,Y,v,axV,zl,title,mytext,plotopU);
  view(viewU);
 
  figure(2+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*(1+slope) 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*(1+slope) lowv highv];
  plotU(X,Y,v,axV,zl,title,mytext,plotopU);
  view(viewU);
 
  figure(3+ifig); clf;
  subplot(2,1,1)
  set(gca,'FontSize',[16]);
  zltop='u(x,y,t)';title='';
  plotopUtop='height'; axUtop=[Lx Rx Ly Ry*(1+slope) lowu-5*(highu-lowu) highu+5*(highu-lowu)];
  axis(axUtop);
  plotU(X,Y,u,[],zltop,title,mytext,plotopUtop);
  view(viewUtop);
  colormap(mycolor);
%  colorbar
  subplot(2,1,2)
  set(gca,'FontSize',[16]);
  zltop='v(x,y,t)';title='';
  plotopUtop='height'; axVtop=[Lx Rx Ly Ry*(1+slope) lowv-5*(highv-lowv) highv+5*(highv-lowv)];
  axis(axVtop);
  plotU(X,Y,v,[],zltop,title,mytext,plotopUtop);
  view(viewUtop);
  colormap(mycolor);
%  colorbar

  figure(4+ifig);clf
  ZZ = radius*cos(theta)*linspace(1,1+slope,Nx);
  XX = radius*sin(theta)*linspace(1,1+slope,Nx);
  YY = ones(size(theta))*x;
  surf (XX,YY,ZZ,u(1:Ny-1,:))
  shading interp
  colormap(mycolor)
  axis equal
  view(tailview)
  hold on
  plot3(radius*sin(theta),Lx*ones(length(y)-1,1),radius*cos(theta),'k','linewidth',1)
  plot3(radius*(1+slope*1.01)*sin(theta),Rx*ones(length(y)-1,1),radius*(1+slope*1.01)*cos(theta),'k','linewidth',1)
  hold off

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

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

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

% no flux + predioic BCs
 u(1,1:Nx) = u(Ny,1:Nx);
 u(Ny,1:Nx)= u(Ny-1,1:Nx);
 u(1:Ny,1) = u(1:Ny,2);
 u(1:Ny,Nx)= u(1:Ny,Nx-1);
 v(1,1:Nx) = v(Ny,1:Nx);
 v(Ny,1:Nx)= v(Ny-1,1:Nx);
 v(1:Ny,1) = v(1:Ny,2);
 v(1:Ny,Nx)= v(1:Ny,Nx-1);
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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(2+ifig); clf;
       subplot(1,2,1)
       set(gca,'FontSize',[16]);
       mytext=['T=(' num2str(kk*dt) ')/(' num2str(tf) ')'];
       plotU(X,Y,u,axU,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,axV,zl,title,mytext,plotopU);
       view(viewU);

       figure(3+ifig); clf;
       subplot(2,1,1)
       set(gca,'FontSize',[16]);
       plotU(X,Y,u,[],zltop,title,mytext,plotopUtop);
       view(viewUtop);
       colormap(mycolor);
%       colorbar
       axis tight;
       subplot(2,1,2)
       set(gca,'FontSize',[16]);
       plotU(X,Y,v,[],zltop,title,mytext,plotopUtop);
       view(viewUtop);
%       colorbar
       colormap(mycolor);
       axis tight;

       figure(4+ifig);clf
       ZZ = radius*cos(theta)*linspace(1,1+slope,Nx);
       XX = radius*sin(theta)*linspace(1,1+slope,Nx);
       YY = ones(size(theta))*x;
       surf (XX,YY,ZZ,u(1:Ny-1,:))
       shading interp
       colormap(mycolor)
       axis equal
       view(tailview)
       hold on
       plot3(radius*sin(theta),Lx*ones(length(y)-1,1),radius*cos(theta),'k','linewidth',1)
       plot3(radius*(1+slope*1.01)*sin(theta),Rx*ones(length(y)-1,1),radius*(1+slope*1.01)*cos(theta),'k','linewidth',1)
       hold off


       pause(0.001);
    end

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



