C----------------------------------------------------------------------
C----------------------------------------------------------------------
C   PP.f :     Predator prey reaction-diffusion model
C
C
C Part of the Electronic Supplementary Material for 
C Sherratt, J.A. & Smith M.J. (2007) 
C Periodic Travelling Waves in Cyclic Populations: 
C Field studies and reaction-diffusion models
C
C
C This is the AUTO code used to calculate the wave families in
C Figure 2(a) of the paper. Instructions for running this code
C are given in the online appendix. 
C
C SUBROUTINE FUNC defines the four ODEs that have travelling wave 
C solutions; these are equations A3 in the Electronic Supplementary
C Material.
C
C SUBROUTINE STPNT defines the parameter values
C
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C
      SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP)
C     ---------- ----
C
C Evaluates the algebraic equations or ODE right hand side
C
C Input arguments :
C      NDIM   :   Dimension of the ODE system 
C      U      :   State variables
C      ICP    :   Array indicating the free parameter(s)
C      PAR    :   Equation parameters
C
C Values to be returned :
C      F      :   ODE right hand side values
C
C Normally unused Jacobian arguments : IJAC, DFDU, DFDP (see manual)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(NDIM), PAR(*), F(NDIM), ICP(*)
C
C Define the parameters other than the dispesal rates
       psigma=0.15d0 !prey to predator conversion rate
       pmu=0.05d0 !predator death rate
       pkappa=0.2d0 !half-saturation constant in prey consumption by predators

C Do some pre-calculations and work out the dispersal rates from PAR(2)
       U1=U(1)
       U2=U(2)
       U3=U(3)
       U4=U(4)
       fkinet=(U1*(1-U1))-((U1*U3)/(U1+pkappa)) !the prey kinetic equation
       gkinet=((psigma*U3*U1)/(U1+pkappa))-(pmu*U3) !the predator kinetic equation
       Du=(PAR(2))**0.5 !The prey dispersal rate
       Dv=(1/PAR(2))**0.5 !The predator dispersal rate
C
C These four equations calculate the ODEs
       F(1)=U2
       F(2)=((-PAR(1)**2)/Du)*(U2+fkinet)
       F(3)=U4
       F(4)=((-PAR(1)**2)/Dv)*(U4+gkinet)
C
      RETURN
      END
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C
      SUBROUTINE STPNT(NDIM,U,PAR)
C     ---------- -----
C
C Input arguments :
C      NDIM   :   Dimension of the ODE system 
C
C Values to be returned :
C      U      :   A starting solution vector
C      PAR    :   The corresponding equation-parameter values
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(NDIM), PAR(*)
C
C Initialise wave speed and the ratio of the dispersal rates
C These are continuation parameters, i.e. they can be varied
C by auto to study their effects on the travelling wave family
C
       PAR(1)=1.e-1 !wave speed
       PAR(2)=1.0d0 !alpha; the ratio of the dispersal rates
C
C The solution gets initialised with the steady state values of the 
C kinetic equations. We need to provide the parameter values again.

       psigma=0.15d0 !prey to predator conversion rate
       pmu=0.05d0 !predator death rate
       pkappa=0.2d0 !half-saturation constant in prey consumption by predators

       U(1)=(pmu*pkappa)/(psigma-pmu)
       U(2)=0.
       U(3)=(1-U(1))*(pkappa+U(1))
       U(4)=0.
C
      RETURN
      END
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C The following subroutines are not used here,
C but they must be supplied as dummy routines
C
      SUBROUTINE BCND 
      RETURN 
      END 
C 
      SUBROUTINE ICND 
      RETURN 
      END 
C 
      SUBROUTINE FOPT 
      RETURN 
      END 
C 
      SUBROUTINE PVLS
      RETURN
      END
C----------------------------------------------------------------------
C----------------------------------------------------------------------