%
%	LOGISTIC-GAUSSIAN including ALLEE Effect
%
%	matlab code to iterate a logistic growth single population model
%	with random (diffusive) movement.
%
%	The model itself looks like thi
%		h* = k h(t)^2/(A^2 + h(t)^2)
%	
%	This population is then allowed to disperse using a spatial
%	convolution, * below, (like a long step of  a heat equation):
%	
%		h(t+1) = p_h(x,y) * h
%	
%	p_h can be any probability function representing the
%	dispersal of the species; I have used Gaussians below.
%
%	Parameters appearing (and their interpretation):
%
%		r	reproduction rate
%		k	carrying capacity
%	
%	Parameters used for the code:
%
%		np	number of points in x and y directions (a power of 2)
%		xl	length of domain in x and y directions
%		dx	grid spacing
%		dt	a time-step length
%		muh -- diffusion parameter for the species
%		ngens --	number of steps to run the code (number of 
%			generations

%
%	set parameters and spatial grids
% 
%	use an imported map (already made available)
landmap=a;	%make sure a is created already!!!
%landmap will contain values from 0 to 10000 
%representing the ha forest per 10*10 km gridcell

ncols=length(landmap(1,:));	%isn't there a simpler way to get these?
nrows=length(landmap(:,1));
%concatenate to obtain a 512 * 512 size landscape
fillX=zeros(nrows,512-ncols);
land2=[landmap fillX];
fillY=zeros(512-nrows,512);
landscape=[land2;fillY];


np=512; muh=.2; dt=1; xl=10; dx=2*xl/np; allee=0.1; 
trunc=dx*5;%truncate dispersal kernel at distance of 5 cells = 50 km
x=linspace(-xl,xl-dx,np);
y=x;
[X,Y]=meshgrid(x,y);
ngens=15;	%not too high for the movie!
allZero=zeros(np);
%
%	set up spatial parameters
%

r=2.75*ones(np);
%first experiment; let k depend on the landscape
%to avoid Divide by Zero, add 0.001 as a minimum value
k=0.001+landscape/2500;
	%now k should be between 0.001 and 4.0
%k=4*ones(np);

%	Set up stuff for movie
M=moviein(ngens); %	If you want movies

%
%	Set up initial conditions
% 
%h0=.5*rand(np);
h0=(rand(np)<=0.025);	%only in 2.5 percent of the cells
%h0=(-1<= X & X <= 1 & -1 <= Y & Y <= 1);

%
%	Define movement kernel for the species (Gaussian)	
% 

hhker=exp(-(X.^2+Y.^2)/(4*muh*dt))/(4*pi*dt*muh);
%truncate the kernel at distance 1 (NB 
hhker=hhker.*(sqrt(X.^2+Y.^2)<=trunc);
hhker=hhker/(dx^2*trapz(trapz(hhker)));
Fhker=fft2(hhker)*dx^2;	%	FFT is taken because we will use it often,
% there are two factors of dx; one for each dimension

%	Now we are in a position to iterate for
% 	a number of generations equal to ngens
%	but first plot a cross-section of carrying capacity
%plot(k(225,:),'r'), axis([200 400 0 4]);

for j=1:ngens

%
%	Advance the life cycles using logistic model
%

%hn=h0+h0.*r.*(1-h0./k).*(h0-allee)./k;
%this (N-A)/K allee effect was removed, as it caused NUM-instability
hn=k.*h0.^2./(allee.^2+h0.^2);	
%this works ok, what bothers me is the absence of intrinsic growth rate r in the eq.
hn=max(hn,allZero);

%
% 	Do the dispersal st
  fhn=fft2(hn);		%	First take fft of the species
%	Now do the convolutions.
%	The fftshift serves to center the probability functions
  h0=real(fftshift(ifft2(Fhker.*fhn)));
%	Now plot the solutio
pcolor(X,Y,h0), shading flat, colormap(hot), axis square, colorbar;
%hold on,plot(h0(225,:)), axis([200 400 0 4]), hold off;
pause(0.5);
M(:,j) = getframe; %	If you want map+movies

end	
%	Here ends the iteration
%  NB To play the movie,  use this command:
%  movie(M,2,8)
