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

%   Define parameters

tfinal=20;
mu=1; nu=1; d=1; t=4;
nsteps=tfinal/t;
V=exp(-(X-5).^2./4-(Y-5).^2);      % these are the victims' odors

f=V./(d+V);					% pheromone response function

%   Set initial population of predators

Pt=((X.^2+Y.^2)<=4); 

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

%   Define dispersal kernel

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

for i=1:nsteps,

    P=Pt;
    numP=dx^2*trapz(trapz(P));

     theta=P.*exp(-nu/mu*f);     % transform population for jim-taxis

%   Disperse transformed population 
    ftheta=fft2(theta);
    thetat=real(fftshift(ifft2(ftheta.*fK)));
    
    %   transform back
    Pt=exp(nu/mu*f).*thetat;

%   Normalize the dispersed population

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

%   Plot the results in comparison with the victims

pcolor(X,Y,Pt), shading interp, colormap hot, axis square
hold on, contour(X,Y,V,'g'),  hold off
pause

end  % of time step iteration


