! this Fortran 90 code is free for all to use and extend as they wish
! it requires a Fortran compiler (e.g. gfortran) to create an executable file
! for clarity it is best viewed with an editor that interprets Fortran (e.g. emacs)
! for those unfamiliar with Fortran 90, the ! symbol is used for comments

! *** numerical integration of PCF equations ***
! *** Thea Newman, Solaravus, June 18th 2020 ***

! supplementary file to Solaravus Technical Report 003
! "Epidemic suppression via an emergent pre-conditioning field"
! available on www.solaravus.com/news

! this code integrates the boxed differential equations in Appendix B
! using 2nd order Runge-Kutta for integration

! note: the choice of dt must be much smaller than the reciprocal of the smallest non-zero rate constant
! an automated definition of dt is not used since this code is designed for "hands-on" rather than "blind" use

! code still managable in size (ca 60 lines of instruction) as single file
! if significant additional code required, then modular form is recommended

! main variables:

! s - density of susceptibles
! d - density of diseased (i.e. infecteds - not using i due to its common use as an integer symbol)
! p - density of pre-conditioned
! r - density of recovereds (and/or dead)
! f - pre-conditioning field (phi)

! main parameters (Greek symbols used for rates):

! alpha  - rate of S -> P (pre-conditioning)
! beta   - rate of S -> D (infection)
! gamma  - rate of D -> R (recovery/death)
! eps    - rate of P -> S (reversion)
! kappa  - rate of production of phi (f) from D
! lambda - rate of decay of phi (f) field

program pcf

  ! all variables to be explicitly defined
  implicit none
  
  ! time increment and output
  integer::it,it_output 
  real*8::dt,time,time_max,time_output

  ! density variables w/ 'total' to check conservation of density (total density = 1)
  real*8::s,d,p,r,f,total
  ! "original values" for use in RK algorithm
  real*8::s_o,d_o,p_o,r_o,f_o
  ! explicit increments of densities for clarity in algorithm code
  real*8::delta_s,delta_d,delta_p,delta_r,delta_f
  ! change in recovered (or dead) on output time-scale (providing epidemic curve)
  real*8::r_prior,d_r

  ! initial size of infection
  real*8::d_init

  ! rate parameters
  real*8::alpha,beta,gamma,eps,kappa,lambda

  ! open data file for recording time evolution of densities
  open(unit=101,file='PCF_output_001',status='unknown')

  ! set values of rate parameters
  alpha=0.4d0
  beta=0.25d0
  gamma=0.1d0
  eps=0.01d0
  kappa=1.d0 ! kappa generally set to unity to fix time-scale
  lambda=0.05d0

  ! set time parameters (where one time unit = 1/kappa)
  dt=0.001d0 ! check this is sufficiently small for stable, reproducible results
  it_output=100
  time_output=dt*it_output
  time_max=1000.d0

  ! set initial size of infection density (relative to total population density = 1.0)
  d_init=1.d-6
  
  ! set initial values of densities - fixed by net population density constant at unity
  s=1.d0-d_init 
  d=d_init 
  p=0.d0
  r=0.d0
  f=0.d0

  ! initialise r_prior (for calculating epidemic curve)
  r_prior=0.d0

  ! initialise time
  it=0
  time=0.d0
  
  ! iterative loop - integrate equations until time_max reached
  do while (time.le.time_max)

     ! increment time
     it=it+1
     time=dt*it

     !! integrate forward using 2nd order RK !!

     ! record original values
     s_o=s
     d_o=d
     p_o=p
     r_o=r
     f_o=f

     ! increments based on 1st order differential equations
     delta_d=beta*d*s-gamma*d
     delta_p=alpha*f*s-eps*p
     delta_r=gamma*d
     delta_f=kappa*d-lambda*f

     ! from conservation of density:
     delta_s=-delta_d-delta_p-delta_r

     ! integrate forward 1/2 time step from original values
     s=s_o+0.5d0*delta_s*dt
     d=d_o+0.5d0*delta_d*dt
     p=p_o+0.5d0*delta_p*dt
     r=r_o+0.5d0*delta_r*dt
     f=f_o+0.5d0*delta_f*dt

     ! re-calculate increments from mid-point values (2nd order RK)
     delta_d=beta*d*s-gamma*d
     delta_p=alpha*f*s-eps*p
     delta_r=gamma*d
     delta_f=kappa*d-lambda*f

     ! from conservation of density:
     delta_s=-delta_d-delta_p-delta_r

     ! integrate forward whole time step from original values
     s=s_o+delta_s*dt
     d=d_o+delta_d*dt
     p=p_o+delta_p*dt
     r=r_o+delta_r*dt
     f=f_o+delta_f*dt

     ! record total density as check of conservation
     total=s+d+p+r
     
     ! output data periodically
     if (mod(it,it_output).eq.0) then
        d_r=(r-r_prior)/time_output ! epidemic curve scaled from output time step to unit time step
        r_prior=r
        write(101,*)time,s,d,p,r,f,d_r,total
     end if
     
  end do

  ! close data file
  close(unit=101)

end program pcf