% -------------------------------------------------------------------------
% 
% The Mathematics behind biological invasions:
%  
%  - University of Utah,
%    Department of Mathematics
%    June 2 - 13, 2003
%  
%  - Fred Adler
%    Mark Lewis
%    Mike Neubert
%    Nancy Sundell-Turner
% 
% Numerical solution of a 1-D IDE using FFTs.
% 
% -------------------------------------------------------------------------

%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Pre - Initialization:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
clf;
F_P = [370 240 900 700];
Fig = gcf;
set(Fig, 'position', F_P);


%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Initialization:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% System Constants:
Generation_Count = 20;


%%%%% Spatial Discretization:
% The physical domain is defined for -1 <= x <= 1.
% However, for the FFT, the computational domain is buffered with 
%  zeros for -2 <= x < -1 and 1 <= x <= 2.
% (The number of grid points should be a power of 2 for the FFT solve.)
P      = 10;                % 2^P data points for the Physical domain
Px     = 2^(P+1);           % Number of data points in the computational domain
x_step = 4 / (2^(P+1));     % 2^(P+1) data points for the computational domain
x      = [-2.0 : x_step : 2.0 - 0.001*x_step];
x_plot = [-1.0 : x_step : 1.0 - 0.001*x_step];


%%%%% Wave Front:
% Population threshold for the wave speed calculation.
wave_front = 0.01;

% To suppress numerical error, a week Allee effect is introduced.
Allee_tol = 10^(-15) * ones( size( x ) );




%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Dispersal Kernel:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Normal Distribution:
t_s = .01;
D   = 0.05;
K   = 1/2.0 * 1/sqrt(pi*D*t_s) * exp(-(x.^2)/(4*D*t_s));

%%%%% Laplace (or double-exponential) distribution:
% a = 600;
% D = 0.05;
% K = 1/2.0 * sqrt(a/D) * exp(-sqrt(a/D)*abs(x));

%%%%% Cauchy Distribution:
% beta = .001;
% K    = 1/pi * beta ./ (beta^2 + x.^2);

%%%%% Finite Moment Kernel:
% alpha = 150;
% K     = alpha^2/4 * exp( -alpha * sqrt(abs(x)) );


%%%%% Plot the dispersal Kernel:
subplot(3,2,1);
plot( x_plot, K(2^(P-1)+1:2^(P+1)-2^(P-1)), 'k-');
xlabel( 'x');
ylabel( 'y');
title(  'Dispersal Kernel');


%%%%% Compute the FFT of the Dispersal Kernel:
% Swap the right and left half buffers:
K_rr = fftshift( K );
K_w  = fft( K_rr );




%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Growth Functions:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Beverton-Holt:
% Growth parameters:
R_0   = 2.0;
C_Cap = 1.0;
%% Plot variables:
N_plot  = [0.0 : 0.05 : 1.5 * C_Cap];
FN_plot = (R_0 * N_plot) ./ (1 + ((R_0 - 1)/C_Cap)*N_plot);

%%%%% Ricker:
% Growth parameters:
 %R     = 1.4;
% C_Cap = 1.0;
% Plot variables:
%N_plot  = [0.0 : 0.05 : 1.5 * C_Cap];
 %FN_plot = N_plot .* exp(R * (1 - (N_plot/C_Cap)));

%%%%% Linear;
% Growth parameters:
% M     = 1.0;
% Plot variables:
% N_plot  = [0.0 : 0.05 : 1.5 * C_Cap];
% FN_plot = M * N_plot;


%%%%% Plot the dispersal Kernel:
subplot(3,2,2);
plot( N_plot, FN_plot, 'k-');
hold on
plot( N_plot, N_plot, 'k--');
hold off
xlabel( 'N');
ylabel( 'F(N)');
title(  'Growth Function');




%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Define the Initial Population:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Define the initial population:
% This is a delta function IC:
N         = zeros( 1, Px);
N(Px/2+1) = 1;


%%%%% FFT the initial population:
% Swap the right and left half buffers:
N_rr      = fftshift( N );
N_w       = fft( N_rr );


%%%%% Plot the initial population:
subplot(3,1,2);
plot( x_plot, N(2^(P-1)+1:2^(P+1)-2^(P-1)), 'k-')
xlabel( 'x');
ylabel( 'N_{t}');
title(  'Population Density');



%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Compute the solution to the IDE:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Compute N_{t}(x) for t = 1 to Generation_Count:
for i = 1:1:Generation_Count

%%%%% Dispersal Stage:
% Convolve the dispersal kernel K with N_k-1:
  Np1_w     = ( K_w .* N_w );

% IFFT the population:
% (Note, the factor 2^(-(P-1)) is needed for the convolution theorem.)
  Np1       = 2^(-(P-1)) * ifft( Np1_w );
  Np1       = real( Np1 );      % Zero any complex component of the FFT.
  
  
% Implement absorbing boundary conditions:
  Np1(2^(P-1)+1:Px - (2^(P-1))) = zeros(1,2^(P));
  
  
%%%%% Growth Stage:
% Week Allee effect:
% (Needed for any numerical noise at the leading edge of the population.)
  N_noise = le( Allee_tol, Np1);
  Np1     = Np1 .* N_noise;
  
% Beverton-Holt:
  Np1 = (R_0 * Np1) ./ (1 + ((R_0 - 1)/C_Cap)*Np1);
  
% Ricker:
 % Np1 = Np1 .* exp(R * (1 - (Np1/C_Cap)));
  
% Linear;
%  Np1 = M * Np1;
  
  
%%%%% FFT the new population:
  N_w = fft( Np1 );
  
  
%%%%% Plot the new generation:
  subplot(3,1,2);
  Np1_plot  = fftshift( Np1 );
  hold on
  plot( x_plot, Np1_plot(2^(P-1)+1:2^(P+1)-2^(P-1)), 'k-')
  hold off
  
  
%%%%% Compute the location of the wave front:
  T_1 = Np1_plot(2^(P):2^(P+1)-2^(P-1)) - wave_front;
  T_2 = sign( T_1);
  [Y,ndx(i)]  = min(T_2);
  
%%%%% Plot the location of the wave front:
  subplot(3,1,3);
  plot( [1:1:i], ndx*x_step, 'ko-');
  axis( [0, Generation_Count, 0, 1.1*max(ndx*x_step)] );
  xlabel( ' Generation');
  ylabel( 'x');
  title(  'Wave front location');
  
  
%%%%% Update the generation count:
  Gen_Num = i
  pause( 0.5 );    % Pause to refresh the screen.
  
end