% cell polarization   on a periodic ring with a moving stimulus
clear

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

 

%parameters
L  = 1; %perimeter wlog = 1
Du = 0.0009;
Dv = 0.06;
del = 1;  % wlog
K = 1; % wlog
gam = 0.5;
ampmx = 0.1;
tgrate = 0.05;
tgmx = 0.7;
 
% This file looks at solutions of a diffusion reaction system
% on a periodic ring.  

u = [0 :0.01:2];

f= gam*u.^2./(1+u.^2);
 Vt =  4.2;
 

figure(1)
plot(u,f*Vt-u,'linewidth',2)
hold on
plot(u,zeros(length(u),1),'--')
hold off

xlabel('u','fontsize',20)
ylabel('f(u,V_T)','fontsize',20)
 
N=100; % number of grid points
M = 150;  % number of time steps
kstep = 3;
 
dxc = L/N;
dt = 0.2;
scalu = Du*dt/(2*dxc^2);

scalv = Dv*dt/(2*dxc^2);
 
%specify initial data
    X = dxc*(1:N)';
    V0 = Vt;
   V = V0*ones(N,1);
   U = zeros(N,1);
  
 
t=0;
T = 0;

% uses Crank Nicolson to solve the diffusion equation

% set up matrices
Agu = scalu*ones(N-1,1);
Agv = scalv*ones(N-1,1);
 

Aexpu = -diag([Agu;scalu]+[scalu;Agu]) +diag(Agu,-1)+diag(Agu,1) ;
Aexpu(N,1)=scalu;
Aexpu(1,N)=scalu;

Amu = sparse(diag(ones(N,1))-Aexpu);

Aexpv = -diag([Agv;scalv]+[scalv;Agv]) +diag(Agv,-1)+diag(Agv,1) ;
Aexpv(N,1)=scalv;
Aexpv(1,N)=scalv;

Amv = sparse(diag(ones(N,1))-Aexpv);

plotlist=[69,246,0];  %these were selected by examination of the answer
pk = 1;

Cx = [1;cos(2*pi*X)];
Sx = [0; sin(2*pi*X)];
 for n=1:M
 
for k = 1:kstep	
    midx =  tgmx*sin(tgrate*t)^2-0.25; %location of the moving stimulus center
    signal = ((1-cos(2*pi*(X-midx)))/2).^10;
    
 
ko =  ampmx*signal;
f =  ko+gam*U.^2./(K^2+U.^2);
Fv = V.*f-del*U;

U_t =  scalu*([U(2:N);U(1)] - 2*U+[U(N);U(1:N-1)]) +dt*(Fv);
V_t =  scalv*([V(2:N);V(1)] - 2*V+[V(N);V(1:N-1)]) -dt*(Fv);

U = Amu\(U+U_t);
V = Amv\(V+V_t);
t = t + dt;

end
 
%track the location of the maximum
mxstim(n) = midx+0.5;
[mxu,mj] = max(U);
mxU(n) = X(mj);
T =[T;t];

 
[Umx,jumx] = max(U);
[Vmn,jvmx] = min(V);
[Smx,jsmx] = max(ko);
 % uncomment the next 9 lines to see the time dependent simulation
 
  figure(2)
%    
  
subplot(2,1,2)
plot(X,U,X,V, X,ko,'r--','linewidth',2)
hold on
plot(X(jumx),Umx,'*',X(jsmx),Smx,'*',X(jvmx),Vmn,'*','linewidth',4)

hold off
  legend('U','V','S')

  axis([0 L 0 4])
  xlabel('x','fontsize',20)
  ylabel ('U,V,S','fontsize',20)
  
  
 R = 1+0.2*[U(end);U];
 R0 = 1.5;

subplot(2,1,1)

 plot3(R.*Cx,R.*Sx,[U(end);U],...
     R0*Cx,R0*Sx,[ko(end);ko],...
      R0*Cx,R0*Sx,zeros(N+1,1),'--',...
    Cx,Sx,zeros(N+1,1),'--',...
    R.*Cx,R.*Sx,zeros(N+1,1),'linewidth',2)
hold on
plot3(R(jumx)*cos(2*pi*X(jumx)),R(jumx)*sin(2*pi*X(jumx)),Umx,'*',...
    R(jumx)*cos(2*pi*X(jumx)),R(jumx)*sin(2*pi*X(jumx)),0,'*',...
    R0*cos(2*pi*X(jsmx)),R0*sin(2*pi*X(jsmx)),Smx,'*','linewidth',4)
hold off
legend('U','S')
axis([-1.5 1.5 -1.5 1.5 0 1.7])
zlabel('U','fontsize',20)
pause(0.1)
% %this command saves files for later making a movie with qt-7

% the simulation runs much faster if this line is commented out
   %saveas(gcf, sprintf('./movie_files/frame%03g',n),'jpeg')
 
% the command to make a movie from these frames is

% ffmpeg -f image2 -pattern_type glob -framerate 12 -i 'frame*.jpg' foo.avi
 end
    
  
hold off
        
 figure(4)
 plot(mxstim,mxU,'linewidth',2)
 xlabel('Stimulus location','fontsize',20)
 ylabel('Pole location','fontsize',20)
 annotation('textarrow',[0.5 .6],[0.42 .52],'linewidth',2)
 annotation('textarrow',[0.5 .4],[0.58 .48],'linewidth',2)

  
 
  figure(5)
 plot(T(2:end),mxU,T(2:end),mxstim,'linewidth',2)
 xlabel('t','fontsize',20)
 ylabel('Stimulus-Pole locations','fontsize',20)
 legend('max(U)','Max(S)','fontsize',16)
 
  