(* :Title: SurplusProduction              *)
(* :Author: Arne Eide            *)
(* :Summary:
This package contains functions to the ECONMULT model.
*)
(* :Context: EconMult`SurplusProduction`  *)
(* :Package Version: 1.0         *)
(* :Copyright: Arne Eide         *)
(* :History:
	Version 1.0 by Arne Eide 2004
*)
(* :Keywords:
	Surplus Production Models,
	Bioeconomics
*)
(* :Source:                      *)
(* :Warning: None.               *)
(* :Mathematica Version: 5.0     *)
(* :Limitation: None.            *)
(* :Discussion:                  *)


BeginPackage["EconMult`SurplusProduction`"]

$GrowthModels::usage                   = "$GrowthModels is a package variable containing the currently implemented growth models in the SurplusProduction package"

BevertonHolt::usage                    = "BevertonHolt is the valid value of the SurplusProduction option GrowthModel."
BiomassMaximum::usage                  = "BiomassMaximum is an option of SurplusProduction with default value 'k'"
CatchabilityCoefficient                = "CatchabilityCoefficient is an option of SurplusProduction with default value 'q'"
GompertzFox::usage                     = "GompertzFox is the valid value of the SurplusProduction option GrowthModel."
GrowthModel::usage                     = "GrowthModel is an option of SurplusProduction with default value 'VerhulstSchaefer'"
IntrinsicGrowthRate::usage             = "IntrinsicGrowthRate is an option of SurplusProduction with default value 'r'"
LengthWeightRelation::usage            = "LengthWeightRelation is an option of SurplusProduction with default value 'b'"
MaximumSustainableYield::usage         = "MaximumSustainableYield is an option of SurplusProduction with default value 'msy'"
RichardsPellaTomlinson::usage          = "RichardsPellaTomlinson is the valid value of the SurplusProduction option GrowthModel."
RichardsPellaTomlinsonParameter::usage = "RichardsPellaTomlinsonParameter is an option of SurplusProduction with default value 'm'"
SurplusProduction::usage               = "SurplusProduction[x, opts] gives the surplus growth related to stock biomass x."
UseMSY::usage                          = "UseMSY is an option of SurplusProduction with default value 'True'"
VerhulstSchaefer::usage                = "VerhulstSchaefer is the default value of the SurplusProduction option GrowthModel."
b::usage                               = "b is the default value of the SurplusProduction option LengthWeightRelation."
k::usage                               = "k is the default value of the SurplusProduction option BiomassMaximum."
m::usage                               = "m is the default value of the SurplusProduction option RichardsPellaTomlinsonParameter."
msy::usage                             = "msy is the default value of the SurplusProduction option MaximumSustainableYield."
q::usage                               = "q is the default value of the SurplusProduction option CatchabilityCoefficient."
r::usage                               = "r is the default value of the SurplusProduction option IntrinsicGrowthRate."
BevertonHoltCatch::usage               = "BevertonHoltCatch[f, opts] gives the equilibrium BevertonHoltCatch of fishing effort 'f'. Default options is the SurplusProduction options."
BevertonHoltGrowth::usage              = "BevertonHoltGrowth[f, opts] gives the equilibrium BevertonHoltGrowth of fishing effort 'f'. Default options is the SurplusProduction options."
GompertzFoxCatch::usage                = "GompertzFoxCatch[f, opts] gives the equilibrium GompertzFoxCatch of fishing effort 'f'. Default options is the SurplusProduction options."
GompertzFoxGrowth::usage               = "GompertzFoxGrowth[f, opts] gives the equilibrium BevertonHoltCatch of fishing effort 'f'. Default options is the SurplusProduction options."
RichardsPellaTomlinsonCatch::usage     = "BevertonHoltCatch[f, opts] gives the equilibrium BevertonHoltCatch of fishing effort 'f'. Default options is the SurplusProduction options."
RichardsPellaTomlinsonGrowth::usage    = "BevertonHoltCatch[f, opts] gives the equilibrium BevertonHoltCatch of fishing effort 'f'. Default options is the SurplusProduction options."
VerhulstSchaeferCatch::usage           = "BevertonHoltCatch[f, opts] gives the equilibrium BevertonHoltCatch of fishing effort 'f'. Default options is the SurplusProduction options."
VerhulstSchaeferGrowth::usage          = "BevertonHoltCatch[f, opts] gives the equilibrium BevertonHoltCatch of fishing effort 'f'. Default options is the SurplusProduction options."


Begin["`Private`"]


(*****************************************)
(*       Growth Model Options            *)
(*****************************************)


Options[SurplusProduction] =
     {
      GrowthModel -> VerhulstSchaefer,
      RichardsPellaTomlinsonParameter -> m,
      UseMSY -> True,
      MaximumSustainableYield -> msy,
      CatchabilityCoefficient -> q,
      BiomassMaximum -> k,
      IntrinsicGrowthRate -> r,
      LengthWeightRelation -> b,
      UnitPrice -> p,
      EffortUnitCost -> c,
      CatchUnitCost -> 0
      };


$GrowthModels=
    {
      VerhulstSchaefer,
      GompertzFox,
      RichardsPellaTomlinson,
      BevertonHolt
      };


(*****************************************)
(*  Growth  and Catch Model Algorithms   *)
(*****************************************)


(*           GompertzFox                 *)

GompertzFoxGrowth[x_, opts___] :=
	Module[{usemsy, msyopt, maxbio, ropts},

		{
		usemsy, msyopt, maxbio, ropts
		} =
		{
		UseMSY,
		MaximumSustainableYield,
		BiomassMaximum,
		IntrinsicGrowthRate
		}  /. {opts} /. Options[SurplusProduction];
		If[usemsy,
			-E msyopt (x/maxbio) Log[x/maxbio],
			-ropts x Log[x/maxbio]
		]
	];


GompertzFoxCatch[ff_,opts___]:=
	Module[{usemsy,msyopt,ccopt,maxbio,ropts},

		{
		usemsy,msyopt,ccopt,maxbio,ropts
		} =
		{
		UseMSY,
		MaximumSustainableYield,
		CatchabilityCoefficient,
		BiomassMaximum,
		IntrinsicGrowthRate
		}  /.{opts}/. Options[SurplusProduction];
		If[usemsy,
			ccopt*ff*maxbio * E^(-ff*maxbio*ccopt/(E*msyopt)),
			ccopt*ff*maxbio * E^(-ff*ccopt/ropts)
		]
	];



(*            VerhulstSchaefer            *)

VerhulstSchaeferGrowth[x_, opts___] :=
	Module[{usemsy, msyopt, maxbio, ropts},

		{
		usemsy, msyopt, maxbio, ropts
		} =
		{
		UseMSY,
		MaximumSustainableYield,
		BiomassMaximum,
		IntrinsicGrowthRate
		}  /. {opts} /. Options[SurplusProduction];

		If[usemsy,
			4 msyopt (x/maxbio)(1-x/maxbio),
			ropts x (1-x/maxbio)
		]
	];


VerhulstSchaeferCatch[ff_,opts___]:=
	Module[{usemsy,msyopt,ccopt,maxbio,ropts},

		{
		usemsy,msyopt,ccopt,maxbio,ropts
		} =
		{
		UseMSY,
		MaximumSustainableYield,
		CatchabilityCoefficient,
		BiomassMaximum,
		IntrinsicGrowthRate
		}  /.{opts}/. Options[SurplusProduction];
		If[usemsy,
			ccopt*ff*maxbio - ff maxbio^2 ccopt/(4 msyopt),
			ccopt*ff*maxbio - ff maxbio ccopt/ropts
		]
	];




(*           RichardsPellaTomlinson        *)

RichardsPellaTomlinsonGrowth[x_, opts___] :=
	Module[
		{rptpar, usemsy, msyopt, maxbio, ropts},

		{
		rptpar, usemsy, msyopt, maxbio, ropts
		} =
		{
		RichardsPellaTomlinsonParameter,
		UseMSY,
		MaximumSustainableYield,
		BiomassMaximum,
		IntrinsicGrowthRate
		}  /. {opts} /. Options[SurplusProduction];

		If[rptpar === 0,
			If[usemsy,

				msyopt -msyopt x/maxbio,

				ropts (-maxbio + x)
			],

			If[rptpar === 1 ,

				GompertzFoxGrowth[x, opts],

				If[usemsy,
					msyopt/(rptpar^(1/(1-rptpar))-rptpar^(rptpar/(1-rptpar))) * (x/maxbio) * (1-(x/maxbio)^(rptpar-1)),
					ropts*x*(1-(x/maxbio)^(rptpar-1))

				]
			]
		]
	];



RichardsPellaTomlinsonCatch[ff_, opts___] :=
	Module[
		{rptpar, usemsy, msyopt, ccopt, maxbio, ropts},

		{
		rptpar, usemsy, msyopt, ccopt, maxbio, ropts
		} =
		{
		RichardsPellaTomlinsonParameter,
		UseMSY,
		MaximumSustainableYield,
		CatchabilityCoefficient,
		BiomassMaximum,
		IntrinsicGrowthRate
		}  /. {opts} /. Options[SurplusProduction];

		If[rptpar === 0,

			If[usemsy,
				(ff maxbio msyopt ccopt)/(msyopt + ff maxbio ccopt),
				(ff maxbio ccopt ropts)/(-ff ccopt + ropts)
			],


			If[rptpar === 1 && usemsy, Exp[-ff maxbio ccopt/(E msyopt)]ff maxbio ccopt,

				If[rptpar === 1 ,

					GompertzFoxCatch[x, opts],

					If[usemsy,
						ccopt*f*(maxbio*(1-ff maxbio ccopt (rptpar-1) rptpar^(-rptpar/(rptpar-1))/msyopt)^(1/(rptpar-1))),
						ccopt*f*(maxbio*(1-ff ccopt/ropts)^(1/(rptpar-1)))
					]
				]
			]
		]
	];





(*     BevertonHolt (when ind.gr.rate=mort.rate)   *)


BevertonHoltGrowth[x_, opts___] :=

	Module[{lwrel, usemsy, msyopt, maxbio, ropts},

		{
		lwrel, usemsy, msyopt, maxbio, ropts
		} =
		{
		LengthWeightRelation,
		UseMSY,
		MaximumSustainableYield,
		BiomassMaximum,
		IntrinsicGrowthRate
		}  /. {opts} /. Options[SurplusProduction];


		RichardsPellaTomlinsonGrowth[
			x,
			RichardsPellaTomlinsonParameter	-> lwrel/(1+lwrel),
			UseMSY				-> usemsy,
			MaximumSustainableYield		-> msyopt,
			BiomassMaximum			-> maxbio,
			IntrinsicGrowthRate		-> ropts
		]


	];


BevertonHoltCatch[ff_, opts___] :=

	Module[{lwrel, usemsy, msyopt, maxbio, ropts},

		{
		lwrel, usemsy, msyopt, maxbio, ropts
		} =
		{
		LengthWeightRelation,
		UseMSY,
		MaximumSustainableYield,
		BiomassMaximum,
		IntrinsicGrowthRate
		}  /. {opts} /. Options[SurplusProduction];


		RichardsPellaTomlinsonCatch[
			ff,
			RichardsPellaTomlinsonParameter	-> lwrel/(1+lwrel),
			UseMSY				-> usemsy,
			MaximumSustainableYield		-> msyopt,
			BiomassMaximum			-> maxbio,
			IntrinsicGrowthRate		-> ropts
		]


	];



(*****************************************)
(*          SurplusProduction            *)
(*****************************************)

SurplusProduction[x_,opts___]:=
	Module[{model,price,rptpar,usemsy,msyopt,ccopt,maxbio,ropts,lwrel},

		{
		model,price,rptpar,usemsy,msyopt,ccopt,maxbio,ropts,lwrel
		} =
		{
		GrowthModel,
		UnitPrice,
		RichardsPellaTomlinsonParameter,
		UseMSY,
		MaximumSustainableYield,
		CatchabilityCoefficient,
		BiomassMaximum,
		IntrinsicGrowthRate,
		LengthWeightRelation
		}  /.{opts}/. Options[SurplusProduction];
		ToExpression[ToString[model]<>"Growth"][x,opts]
	];



End[]


EndPackage[]