% cell polarization   on a periodic ring
clear

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

 

%parameters
L  = 1; %perimeter wlog = 1
Du = 0.001;  % diffusion coefficients
Dv = 1;
del = 1;  % wlog
K = 1; % wlog
gam = 0.5;
%parameters governing the signal
ampmx = 0.1;
amprate = 0.005;
ampzero = pi/amprate;

 
formatSpecF = '%6.2f\n';
% This file looks at solutions of a diffusion reaction system
% on a periodic ring.  

 % specify V_T
 Vt =  4.2;

N=100;  % number of grid points
M = 270; % number of time steps
kstep = 24;  %number of intermediate time steps
 
dxc = L/N;  % grid spacing
dt = 0.1;  % time step
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);
  
K0s = U;
Us = U;
t=0;
plotlist = [204,58,0];
pk = 1;

% 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);
 
 for n=1:M
 
for k = 1:kstep	
   signal = ((1-cos(2*pi*X))/2).^10;
    
 kamp = (t<ampzero)*ampmx*sin(amprate*t)^2; 
ko =  kamp*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
% save the amplitude of the solution
 svu(n) = max(U);
 svamp(n) = kamp;
 tm(n) = t;
 Us = [Us,U]; % save the solution
K0s = [K0s,ko];


  
%wwname1 = strcat('t=',sprintf(formatSpecF,t));

% title(strcat('Stimulus Amplitude = ',sprintf(formatSpecF,kamp)),'fontsize',18)
 


 end
hold off
% 
 
color = ['b','g'];
colorstar = ['b*','g*'];
% Plot the hysteretic response curve
figure(1)
plot(svamp,svu,'linewidth',2)
xlabel('Stimulus Amplitude','fontsize',20)
ylabel('Response Amplitude','fontsize',20)
axis([0 0.1 0 1.8])
set(gca,'Ytick',0:0.25:1.5)
hold on

target = 0.04;  % locate places where stimulus amplitude is target
bt = find(svamp>=target);
plotlist = [max(bt),min(bt)];
 
for j = 1:2
n = plotlist(j);
 
     plot(svamp(n),svu(n),'*','linewidth',4)
      
end

 hold off
 
 
 
 
 for j = 1:2
 figure(2)
 n  = plotlist(j);
 plot(X,Us(:,n),color(j), X,K0s(:,n),'r--','linewidth',2)
 hold on
 plot(X(N/2),Us(N/2,n),'*',X(N/2),K0s(N/2,n),'*','linewidth',4)
  
  hold on
  axis([0 L 0 1.6])
  set(gca,'Ytick',0:0.25:1.5)
  
  xlabel('\xi','fontsize',20)
  ylabel ('u, S','fontsize',20)
 end
 hold off
 
  
 
 
 figure(3)
 plot(tm,svu,tm,svamp,'linewidth',2)
 ylabel('S, u Amplitude','fontsize',20)
xlabel('t','fontsize',20)
set(gca,'Ytick',0:0.25:2)
hold on

for j = 1:2 
    n  = plotlist(j);
    plot(tm(n),svamp(n),'k*',tm(n),svu(n),'k*','linewidth',4)
end
legend('u','S','fontsize',20)
hold off
  
 %plots to make a  movie
for n=1:M
figure(4)
subplot(2,1,1)
plot(svamp,svu,svamp(n),svu(n),'*','linewidth',2)
xlabel('Stimulus Amplitude','fontsize',20)
ylabel('Response Amplitude','fontsize',20)

% 
subplot(2,1,2)
  plot(X,Us(:,n), X,K0s(:,n),'r--','linewidth',2)
  axis([0 L 0 1.75])
  xlabel('\xi','fontsize',20)
  ylabel ('u,S','fontsize',20)
  legend('U','S')
 pause(0.1)
 
%   
% %this command saves files for later making a movie with ffmpeg

%the simulation runs much faster if this  next 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
    
 
