% 
%    Testing out Turchin-taxis
% 
%   Set up space:
%
np=64;
xl=10; dx=2*xl/np;
x=linspace(-xl,xl-dx,np);
[X,Y]=meshgrid(x,x);

%   Define parameters

tfinal=20;
D2=4; D1=1.5; d=1; t=1; %t=tfinal;
nsteps=tfinal/t;
%V=100*exp(-(X-5).^2/4-(Y-5).^2);      % these are the victims
H=1.48; rockies;
V=10*B.^2;

V= 10*exp(-(X-4).^2/9-(Y-4).^2/9).*(1-cos(20*pi/xl*X));
figure(2), surf(X,Y,V),shading interp
%numV=dx^2*trapz(trapz(V));             % this is the number  of victims
%muhat=4*del2((D2-D1)*(D2*V)./(d*D1+ D2*V) ,dx,dx);
%muhat=(D2-D1)*(D2*V)./(d*D1+ D2*V);
muhat=(D2-D1)*(V)./(d*D1+ D2*V);

%   Set initial population of predators

Pt=((X.^2+Y.^2)<=1); 
%Pt=100*exp(-(X-5).^2-(Y-5).^2);

%   Define dispersal kernel -- remember we will use 1/D2 for diffusion

K=.25/(pi*(1/D2)*t)*exp(-(X.^2+Y.^2)/(4*(1/D2)*t));
K=.25/(pi*D2*t)*exp(-(X.^2+Y.^2)/(4*D2*t));
fK=dx.^2*fft2(K);

for i=1:nsteps,

    P=Pt;
    numP=dx^2*trapz(trapz(P));
    C=P.*(D2-muhat);
 C=P.*(1-muhat);

%   Disperse population according to turchin-taxis

    fC=fft2(C);
    Ct=real(fftshift(ifft2(fC.*fK)));
    Pt=Ct./(1-muhat);

%   Normalize the dispersed population

    numPt=dx^2*trapz(trapz(Pt));
    Pt=numP/numPt*Pt;

%   Plot the results in comparison with the victims

figure(1),pcolor(X,Y,Pt), shading interp, caxis('manual'), colormap hot, axis square
hold on, contour(X,Y,V,'b'), caxis('manual'),  hold off
pause(.1)

end  % of time step iteration


