(* :Title: PredatorPrey *)

(* :Author: H. Edward Donley
  <hedonley@iup.edu> or
  <http://www.math.iup.edu/people/hedonley.html>
*)

(* :Summary:
  Functional response simulations for predator-prey interactions.
*)

(* :Package Version: 1.0 *)

(* :Keywords: 
  differential equations, population modeling 
*)

(* :Mathematica Version: 2.2 or higher *)

BeginPackage["PredatorPrey`"]

Unprotect[FunctionalResponse, PreyRadius, PredatorDelay, PreyDisplayFlag, PredatorIterations, PreyCaptureProbability];

(* Function usage statements *)

FunctionalResponse::usage = "Simulates functional response for predator-prey interactions.
Can model probabilistic capture of prey and handling time for prey.";

(* Option usage statements *)

PreyRadius::usage = "PreyRadius is an option for FunctionalResponse.  It is a positive number specifying
the maximum distance from which predators can capture prey.  The radius should be a fraction of the length and width
of the 1 by 1 square environment.  The default value is 0.1.";

PredatorDelay::usage = "PredatorDelay is an option for FunctionalResponse.  It is a non-negative number specifying
the time that is required for each predator to process one captured prey.  The default value is 0.";

PreyDisplayFlag::usage = "PreyDisplayFlag is an option for FunctionalResponse.  If the value is True, then the
simulation will be displayed as a graphical animation.  Otherwise, only the data points for the functional response
will be returned from FunctionalResponse.  The default value is True.";

PredatorIterations::usage = "PredatorIterations is an option for FunctionalResponse.  It is a positive integer specifying
the number of times that the predator tries to capture prey.  The default value is 30.";

PreyCaptureProbability::usage = "PreyCaptureProbability is an option for FunctionalResponse.  It is the probability of
the predator capturing a prey when it is within PreyRadius of that individual. Its value should be a number between 0 and 1,
inclusive.  The default value is 1.";

(* Error messages *)

FunctionalResponse::prey =  "Value of argument nPrey, ``, is not a positive integer.";

PreyRadius::prad =  "Value of argument PreyRadius, ``, is not a positive number.";

PredatorDelay::dlay =  "Value of argument delay, ``, is not a non-negative number.";

PreyDisplayFlag::disp = "Value is PreyDisplayFlag, ``, is neither True nor False.";

PredatorIterations::iter = "Value of PredatorIterations, ``, is not a positive integer.";

PreyCaptureProbability::nano = "Value of PreyCaptureProbability, ``, is not a number.";

PreyCaptureProbability::ubnd = "Value of PreyCaptureProbability, ``, is not between 0 and 1, inclusive.";


(* Option default values *)

Options[FunctionalResponse] = {
  PreyRadius -> 0.1,
  PredatorDelay -> 0,
  PreyDisplayFlag -> True,
  PredatorIterations -> 30,
  PreyCaptureProbability -> 1
};

Begin["`Private`"]

(* Function definitions *)

FunctionalResponse[nPrey_Integer, opts___?OptionQ] :=
  Module[
    {radius, delay, displayFlag, iterations, captureProbabilityprey, preyplot, predatorplot, predator, i, j, hits, previousHit, response, eatenpreyplot, numHits},

    (* Process the arguments and options *)

    If[nPrey <= 0,
      Message[FunctionalResponse::prey, nPrey];
      Return[$Failed]
    ];

    radius = PreyRadius /. {opts} /. Options[FunctionalResponse];
    If[radius <= 0,
      Message[PreyRadius::prad, radius];
      Return[$Failed]
    ];

    delay = PredatorDelay /. {opts} /. Options[FunctionalResponse];
    If[delay < 0,
      Message[PredatorDelay::dlay, delay];
      Return[$Failed]
    ];

    displayFlag = PreyDisplayFlag /. {opts} /. Options[FunctionalResponse];
    If[displayFlag =!= False && displayFlag =!= True,
      Message[PreyDisplayFlag::disp, displayFlag];
      Return[$Failed]
    ];

    iterations = PredatorIterations /. {opts} /. Options[FunctionalResponse];
    If[iterations <= 0,
      Message[PredatorIterations::iter, iterations];
      Return[$Failed]
    ];

    captureProbability = PreyCaptureProbability /. {opts} /. Options[FunctionalResponse];
    If[NumberQ[captureProbability] === False,
      Message[PreyCaptureProbability::nano, captureProbability];
      Return[$Failed]
    ];
    If[captureProbability < 0 || captureProbability > 1,
      Message[PreyCaptureProbability::ubnd, captureProbability];
      Return[$Failed]
    ];


    prey = Table[{Random[], Random[]}, {nPrey}];
    If[displayFlag,
      preyplot = Graphics[Map[(({RGBColor[0, 0, 1], Disk[#, radius]})) &, prey]];
      Show[preyplot, PlotRange -> {{0, 1}, {0, 1}}, AspectRatio -> Automatic,  Frame -> True];
    ];
    previousHit = 0;
    response = {};
    numHits = 0;
    For[i = 1, i <= iterations, i++,
      predator = {Random[], Random[]};
      If[displayFlag,
        predatorplot = Graphics[{RGBColor[1, 0, 0], Thickness[.02], 
          {Line[{predator - {radius, 0}, predator + {radius, 0}}], 
           Line[{predator - {0, radius}, predator + {0, radius}}]}}];
      ];
      hits = Position[Map[( ((# - predator) . (# - predator) <= radius^2) && (Random[] <= captureProbability) ) &, prey], True, 1];
      If[hits =!= {},  (* predator is near at least one prey *)
        AppendTo[response, {Length[prey], 1.0/(i - previousHit + delay)}];
        numHits++;
        previousHit = i;
        If[displayFlag,
          eatenpreyplot = Graphics[{RGBColor[0, 1, 1], Disk[prey[[hits[[1]]]][[1]], radius]}];
        ];
        prey = Delete[prey, hits[[1]]];
        If[displayFlag,
          preyplot = Graphics[Map[(({RGBColor[0, 0, 1], Disk[#, radius]})) &, prey]];
          Show[preyplot, eatenpreyplot, predatorplot, PlotRange -> {{0, 1}, {0, 1}}, AspectRatio -> Automatic, Frame -> True];
          Table[Show[preyplot, eatenpreyplot, predatorplot, PlotRange -> {{0, 1}, {0, 1}}, AspectRatio -> Automatic, Frame -> True], {j, 1, Floor[delay]}];
        ],
      (* Else *)
        If[displayFlag,
            Show[preyplot, predatorplot, PlotRange -> {{0, 1}, {0, 1}}, AspectRatio -> Automatic, Frame -> True];
        ];
      ];
    ];
    Return[response];
  ];

End[];

Protect[FunctionalResponse, PreyRadius, PredatorDelay, PreyDisplayFlag, PredatorIterations, PreyCaptureProbability];

EndPackage[];