%% Stochastic Hunter-Prey model in 2D
% This script is a model of two species, one hunter and one prey. The prey
% can reproduce and die naturally, and both the hunter and prey occupy a
% physical location in 2D. Both the hunter and prey move according to
% Browninan motion and if the hunter occupies the same space as the prey it
% is considered to have 'eaten' them.  The prey are an example of branching
% diffusion. The hunter is an example of a diffusion with lingering effects
% so that if the hunter eats, it slows down to "look around". There is also
% a background current to mimic a passing internal wave.
%%
clear all,close all
%%
% The maximum humber of particles.
numparts=1000;
%%
% Number of particles to start with.
numactive=500;
%% Parameters for domain size
L=20;
H=10;
%% Brownian motion parameters
dt=1e-3;
dtrt=sqrt(dt);
sigx=1e-2*H;
sigz=1e-2*H;
%% Birth process parameters
% Uniform random number larger than this gives birth.
birth_thresh=0.975;
%%
% Uniform random number larger than this gives death.
death_thresh=0.97575;
%% Hunter parameters
hunter_rad=0.05*min(L,H);
hunter_sigx0=2e-1*H;
hunter_sigz0=2e-1*H;
hunter_sigx=hunter_sigx0;
hunter_sigz=hunter_sigz0;
xhunter=0.25*L;
zhunter=0.5*H;
%%
% How many did I have to eat to slow down and eat more.
eattostop=10;
%%
% The number of time steps to slowdown for.
lengthstop=200;
%%
whatiate=0;
eating=0;
%% Deterministic motion parameters
% Mode number.
modenum=1;
%%
% Wave number.
k=2*pi/L;
%%
% Vertical wave number.
m=modenum*pi/H;
%%
% Buoyancy frequency.
N0=2e-2;
%%
% Buoyancy frequency squared.
N02=N0^2;
%%
% Frequency.
sigma=(N0*k*H)/sqrt(modenum^2*pi^2+k*H*k*H);
%%
% Period.
myper=2*pi/sigma;
%%
% Wave amplitude.
amp=0.025*100;
%%
% Flow solution to Internal Wave problem.
psi=@(x,z,t) amp*cos(k*x).*sin(m*z)*cos(sigma*t);
u=@(x,z,t) amp*cos(k*x).*cos(m*z)*m*cos(sigma*t);
w=@(x,z,t) amp*k*sin(k*x).*sin(m*z)*cos(sigma*t);

%% Time stepping parameters
numsteps=20;
numouts=20;
%% Allocate and Initialize
% Allocate arrays particles.
x=-1e4*ones(numparts,1);
z=-1e4*ones(numparts,1);
alive=zeros(numparts,1);
%%
% Initialize particles.
x(1:numactive,1)=L*rand(numactive,1);
z(1:numactive,1)=H*rand(numactive,1);
alive(1:numactive)=1;
%%
% Grid for velocity field pictures.
myxs=linspace(0,1*L,1000);
myzs=linspace(0,H,100);
[xx ,zz]=meshgrid(myxs,myzs);
%%
t=0;
totdead=zeros(numouts);
totborn=zeros(numouts);
toteaten=zeros(numouts);
%% Primary Loop
for ii=1:numouts
    totdead(ii)=0;
    totborn(ii)=0;
    toteaten(ii)=0;
    for jj=1:numsteps
        %%
        % Advect and diffuse the active particles.
        pertx=sigx*randn(numactive,1);
        pertz=sigz*randn(numactive,1);
        x(1:numactive,1)=mod(x(1:numactive,1)+alive(1:numactive).*(dt*u(x(1:numactive,1),z(1:numactive,1),t)+dtrt*pertx),L);
        z(1:numactive,1)=z(1:numactive,1)+alive(1:numactive).*(dt*w(x(1:numactive,1),z(1:numactive,1),t)+dtrt*pertz);
        z(1:numactive,1)=H-abs(z(1:numactive,1)-H);
        perthunterx=hunter_sigx*randn(1,1);
        perthunterz=hunter_sigz*randn(1,1);
        xhunter=mod(xhunter+dt*u(xhunter,zhunter,t)+dtrt*perthunterx,L);
        zhunter=zhunter+dt*w(xhunter,zhunter,t)+dtrt*perthunterz;
        zhunter=H-abs(zhunter-H);
        t=t+dt;
        %%
        % Natural death, do this first so you're not always killing the
        % newborns.
        deathpert=rand(numactive,1);
        deathcheck=(deathpert>death_thresh);
        numdead=sum(deathcheck);
        numactive=numactive-numdead;
        totdead(ii)=totdead(ii)+numdead;
        %%
        % Birth new particles.
        birthpert=rand(numactive,1);
        birthcheck=(birthpert>birth_thresh);
        numborn=sum(birthcheck);
        totborn(ii)=totborn(ii)+numborn;
        %%
        % Get born next to parent where parents are automatically at the end
        % Put the newborns at the front of the queue so they are not the
        % first to die.
        dummyx=x(1:numactive);
        dummyz=z(1:numactive);
        x(1:numborn)=x(numactive-numborn+1:numactive);
        z(1:numborn)=z(numactive-numborn+1:numactive);
        x(numborn+1:numactive+numborn)=dummyx;
        z(numborn+1:numactive+numborn)=dummyz;
        numactive=numactive+numborn;
        alive=zeros(numparts,1);
        alive(1:numactive)=1;
        %%
        % Can you eat it?
        huntdist=sqrt((x(1:numactive,1)-xhunter*ones(numactive,1)).^2+(z(1:numactive,1)-zhunter*ones(numactive,1)).^2);
        eaten=find(huntdist<hunter_rad);
        toteaten(ii)=toteaten(ii)+length(eaten);
        alive(eaten)=0;
        x(eaten,1)=-1e10;
        z(eaten,1)=-1e10;
        deadind=find(alive==0);
        didIeat=length(eaten);
        whatiate=whatiate+didIeat;
        eating=eating*(didIeat<eattostop+1)+lengthstop*(didIeat>eattostop);
        sigfactor=0.01*(eating>0)+1*(eating==0);
        %%
        % If the hunter has just eaten it will stop to look around. This is
        % achieved by lowering its standard deviation.
        hunter_sigx=hunter_sigx0*sigfactor;
        hunter_sigz=hunter_sigz0*sigfactor;
        eating=max(eating-1,0);
    end
    %% Plot
    figure(1)
    set(gcf,'DefaultLineLineWidth',1,'DefaultTextFontSize',14,...
        'DefaultAxesFontSize',14);
    clf
    diagfact=1/sqrt(2);
    pcolor(xx,zz,psi(xx,zz,t)/amp),shading flat
    hold on
    plot(x(1:numactive,1),z(1:numactive,1),'k.','Markersize',4)
    plot(xhunter,zhunter,'ko','Markersize',8)
    plot(mod(xhunter+hunter_rad,L),(zhunter),'ko','Markersize',4)
    plot(mod(xhunter-hunter_rad,L),(zhunter),'ko','Markersize',4)
    plot(mod(xhunter,L),(zhunter+hunter_rad),'ko','Markersize',4)
    plot(mod(xhunter,L),(zhunter-hunter_rad),'ko','Markersize',4)
    plot(mod(xhunter+diagfact*hunter_rad,L),(zhunter-diagfact*hunter_rad),'ko','Markersize',4)
    plot(mod(xhunter-diagfact*hunter_rad,L),(zhunter-diagfact*hunter_rad),'ko','Markersize',4)
    plot(mod(xhunter-diagfact*hunter_rad,L),(zhunter+diagfact*hunter_rad),'ko','Markersize',4)
    plot(mod(xhunter+diagfact*hunter_rad,L),(zhunter+diagfact*hunter_rad),'ko','Markersize',4)
    myindalive=(alive(1:numactive)==1);
    title(['Number of active prey ' int2str(sum(myindalive))])
    drawnow
    %% Update
    % Now remove the eaten particles.
    livingind=find(alive==1);
    dummyx=x(livingind);
    dummyz=z(livingind);
    numactive=length(dummyx);
    %%
    % Allocate arrays particles.
    x=-1e4*ones(numparts,1);
    z=-1e4*ones(numparts,1);
    alive=zeros(numparts,1);
    alive(1:numactive)=1;
    x(1:numactive)=dummyx;
    z(1:numactive)=dummyz;
end

