C
C This is the FORTRAN version of program 2.6 from page 41 of
C "Modeling Infectious Disease in humans and animals"
C by Keeling & Rohani.
C 
C It is the SEIR epidemic including births and deaths at an equal rate.
C
C This code is written to be simple, transparent and readily compiled.
C Far more elegant and efficient code can be written.
C
C This code can be compiled using the intel fortran compiler:
C ifort -Vaxlib  -o Program_2_6  Program_2_6.f  
C 


C Main program starts here.

      program main

      REAL beta
      REAL gamma
      REAL mu
      REAL sigma

      REAL S,E,I,R

      REAL S0
      REAL E0
      REAL I0
      REAL MaxTime
      REAL EVERY, step, t

      INTEGER GivesName

      CHARACTER*2000 str, FileName

      COMMON /parameters/ beta, sigma, gamma, mu
      COMMON /variables/ S, E, I, R

      GivesName=iargc()

      if (GivesName.eq.0) then
         beta=1.4247
         sigma=0.07143
         gamma=0.14286
         mu=0.0000391
         S0=0.1
         E0=1.0d-4
         I0=1.0d-4
         MaxTime=21900
      else
 
c     
c     READ IN ALL THE VARIABLES
c

         call getarg(1,FileName)

         open(1,file=FileName,STATUS='OLD',ACCESS='SEQUENTIAL')
  
         
         read(1,*) str
         read(1,*) beta
         
         read(1,*) str
         read(1,*) gamma

         read(1,*) str
         read(1,*) sigma

         read(1,*) str
         read(1,*) mu
         
         read(1,*) str
         read(1,*) S0
         
         read(1,*) str
         read(1,*) E0

         read(1,*) str
         read(1,*) I0
         
         read(1,*) str
         read(1,*) MaxTime 
         
         close(1)
      endif
C
C     Check all variables are OK & set up intitial conditions */
C     
      
      if ( S0.le.0) then
         write(*,*) "ERROR: Initial level of susceptibles (",S0,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( E0.le.0) then
         write(*,*) "ERROR: Initial level of exposed (",E0,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( I0.le.0) then
         write(*,*) "ERROR: Initial level of infectious (",I0,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( beta.le.0) then
         write(*,*) "ERROR: Transmission rate beta (",beta,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( sigma.le.0) then
         write(*,*) "ERROR: Exposed to Infectious rate sigma (",sigma,") 
     .        is less than or equal to zero."
         STOP
      endif

      if ( gamma.le.0) then
         write(*,*) "ERROR: Recovery rate gamma (",gamma,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( mu.le.0) then
         write(*,*) "ERROR: Recovery rate mu (",mu,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( MaxTime.le.0) then
         write(*,*) "ERROR: Maximum run time (",MaxTime,") is 
     .        less than or equal to zero."
         STOP
      endif

      if (S0+E0+I0.ge.1) then
        write(*,*) "WARNING: Initial level of susceptibles+infecteds
     .        (",S0,"+",E0,"+",I0,"=",S0+E0+I0,") is greater than one."
      endif

      if (beta*sigma.lt.(gamma+mu)*(sigma+mu)) then
        write(*,*) "WARNING: Basic reproductive ratio (R_0=",
     .        beta*sigma/((gamma+mu)*(sigma+mu)),") is less than one."
      endif


      S=S0
      E=E0
      I=I0
      R=1-S0-E0-I0
     
C     
C     Find a suitable time-scale for outputs 
C     
      step=0.01/((beta+gamma+mu+sigma))

      Every=1.0/((beta+gamma+mu+sigma))
      if (Every.gt.1) then
         Every=10.0**INT(log10(Every))
      else
         Every=10.0**INT(log10(Every)-1)
      endif
      DO WHILE (MaxTime/Every.gt.10000)
         Every=Every*10.0  
      ENDDO

      open(1,recl=3000,file='Output',ACCESS='SEQUENTIAL')
C     for F77 use
C      open(1,file='Output_Risk',ACCESS='SEQUENTIAL')
C

C
C     The main iteration routine
C     

      t=0

      DO WHILE (t.lt.MaxTime)
         

         CALL Runge_Kutta(step)
         
         t=t+step
         
C     If time has moved on sufficiently, output the current data
         
         if( INT(t/Every).gt.INT((t-step)/Every) ) then
            write(1,*) t,S,E,I,R
         endif
      
      ENDDO
      
      END
  


      SUBROUTINE Runge_Kutta(step)

      REAL InitialPop(4), tmpPop(4)
      REAL dPop1(4), dPop2(4), dPop3(4), dPop4(4)

      REAL S,E,I,R, step
      COMMON /variables/ S, E, I, R

C
C     Integrates the equations one step, using Runge-Kutta 4 
C     Note: we work with arrays rather than variables to make the
C     coding easier
C
      InitialPop(1)=S
      InitialPop(2)=E
      InitialPop(3)=I
      InitialPop(4)=R

      CALL Diff(InitialPop,dPop1)
      do k=1,4
         tmpPop(k)=InitialPop(k)+step*dPop1(k)/2.0
      ENDDO 

      CALL Diff(tmpPop,dPop2)
      do k=1,4
         tmpPop(k)=InitialPop(k)+step*dPop2(k)/2.0
      ENDDO

      CALL Diff(tmpPop,dPop3)
      do k=1,4
         tmpPop(k)=InitialPop(k)+step*dPop3(k)
      ENDDO

      CALL Diff(tmpPop,dPop4)

      do k=1,4
         tmpPop(k)=InitialPop(k)+step*(dPop1(k)/6 + dPop2(k)/3 +
     .        dPop3(k)/3 + dPop4(k)/6)
      ENDDO

      S=tmpPop(1)
      E=tmpPop(2)
      I=tmpPop(3)
      R=tmpPop(4)

      RETURN
      END


 
C The Main Differential Equations.

      SUBROUTINE Diff(Pop, dPop)

      REAL Pop(4), dPop(4)

      REAL beta
      REAL sigma
      REAL gamma
      REAL mu

      COMMON /parameters/ beta, sigma, gamma, mu


C     Set up temporary variables to make the equations look neater

      REAL tmpS, tmpE, tmpI, tmpR

      tmpS=Pop(1)
      tmpE=Pop(2)
      tmpI=Pop(3)
      tmpR=Pop(4)

C
C     The differential equations
C

C     dS/dt =
      dPop(1) = mu - beta*tmpS*tmpI - mu*tmpS

C     dE/dt =
      dPop(2) = beta*tmpS*tmpI - sigma*tmpE - mu*tmpE

C     dI/dt =
      dPop(3) = sigma*tmpE - gamma*tmpI - mu*tmpI

C     dR/dt =
      dPop(4) = gamma*tmpI - mu*tmpR

      RETURN
      END