! 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