% This file solves the resource consumption problem with chemotaxis 
clear all;

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

%moviename = ['MOVIE.qt']
%MakeQTMovie('start', moviename)

L = 200; % length of the domain
du = .01; % diffusion coefficient for u
dg = 0.01; % diffusion coefficient for g
chi = 0.2;
%chi = 0.;
% parameters
alpha = 0.01;
N=100; % number of spatial grid points
M = 20;  % number of steps
h = L/N;
dt = h;
T = 0;

scalu = du*dt/h^2;
scalg = dg*dt/h^2;

X = h*(1:N)';
  %set initial data
 U =0.1* exp(-(X-L/2).^2/10);
 G = ones(N,1)/2-U;
 Us = U; 
 figure(1)
  subplot(2,1,2)
  plot(X,G,'--','linewidth',2)
  xlabel('x','fontsize',20)
   ylabel('g(x,t)','fontsize',20)
  hold on
  subplot(2,1,1)
  plot(X, U,'--','linewidth',2)
   xlabel('x','fontsize',20)
    ylabel('u(x,t)','fontsize',20)
  axis([0 L  0 0.5])
 % axis([0 L -.5 1])
 hold on 
kstep = 200; % number of intermediate steps
t=0;

% uses Crank Nicolson to solve the diffusion equation

% set up matrix

Atmu=(1+scalu)*ones(N,1);
Atmu(1,1) = 1+scalu/2;
Atmu(N,1) = 1+scalu/2;

Amu=diag(Atmu)-diag(scalu/2*ones(N-1,1),1)-diag(scalu/2*ones(N-1,1),-1);

Atmg=(1+scalg)*ones(N,1);
Atmg(1,1) = 1+scalg/2;
Atmg(N,1) = 1+scalg/2;

Amg=diag(Atmg)-diag(scalg/2*ones(N-1,1),1)-diag(scalg/2*ones(N-1,1),-1);

% The algorithm
 
for n=2:M
for k = 1:kstep
F= alpha*U.*G ;
 
G_t =  scalg/2 * ([G(2:N );G(N )] - 2*G + [G(1);G(1:N-1 )]) -F ;  % no flux BC's

 % here is the code that does upwinding:
vel = chi*(G(2:N)-G(1:N-1))/h;  %calculate the velocity  % G is the chemoattractant
msk = (vel>0);
 Flx = [0;vel.*(msk.*U(1:N-1)+(1-msk).*U(2:N));0];% calculate the upwinded flux
Gdiff = (Flx(2:N+1)-Flx(1:N))/h;  % this is the gradient of the flux
 

U_t =  scalu/2 * ([U(2:N );U(N )] - 2*U + [U(1 );U(1:N-1 )]) -dt*Gdiff +F;
 
G = Amg\(G+G_t);
U = Amu\(U+U_t);

t=t + dt;
end
Us = [Us,U];
T = [T,t];
figure(1)
  subplot(2,1,2)
  plot(X,G)
  xlabel('x','fontsize',20)
 %  ylabel('g(x,t)','fontsize',20)
  hold on
  subplot(2,1,1)
  plot(X, U)
   xlabel('x','fontsize',20)
    ylabel('u(x,t)','fontsize',20)
  %axis([0 L  0 0.5])
 % axis([0 L -.5 1])

hold on
%MakeQTMovie('addframe','moviename')
%MakeQTMovie('quality', 1.0)

  pause(0.05)

end
hold off
%MakeQTMovie('finish','moviename')
figure(2)
 [xx,tt]=meshgrid(X,T); 
contour(xx,tt,Us',[0.25 0.25],'linewidth',2)
xlabel('X','fontsize',20)
ylabel('T','fontsize',20)

TL = t;
time=cputime-time ;       % outputs the cpu time taken to solve the equations

time_steps = M;
 