!====================================================================

 module global
   integer, parameter :: kdp=kind(0.d0)
   real (kind=kdp) dt, dt2, dt6
   real (kind=kdp) a, mu, nu                        ! some parameters
   real (kind=kdp), parameter :: pi = 4.d0*atan(1.d0)
 end module

!====================================================================
!
program main
!    
!====================================================================
!
!  This code integrates a 2D nonlinear system of the kind
!   
!  \dot(x1) = F1(x1,x2)
!  \dot(x2) = F2(x1,x2)
!
!  The code uses a 4th order Runge-Kutta algorithm.
!
!  This example illustrates the case of Kermack and McKendrick's 1927
!  epidemics spread model. This example is exercise 6.5.6 in S. Strogatz's 
!  book "Nonlinear dynamics and chaos"
!
!  Copyright 3/2017. F. Piazza, University of Orléans,
!                    Francesco.Piazza@gmail.com 
!
!====================================================================
!
use global
integer is, N
real (kind=kdp)  time, x0, y0, xm, xm1, beta, alpha
real (kind=kdp), dimension(2) :: x    ! the state vector
!
!====================================================================
!  Parameters

 dt   = 0.01d0            ! integration step 
 dt2  = dt/2.d0           !
 dt6  = dt/6.d0           !
 time = 0.d0              !
 is   = 10                ! sampling interval
 N    = 10000             ! Number of integration steps 
 mu   =  0.22d0           ! contagion rate
 nu   =  0.2              ! death rate
                     
! 
!====================================================================
! Initial condition 
!
 x0 = 10.9d0
 y0 = 0.3d0
 z0 = 0.d0
 beta = y0 / x0
 alpha = mu * x0 / nu
 xm = x0*(1. - beta*alpha - beta*alpha**2*(1.-beta/2.) ) !- &
           !beta*alpha**3*(1. - 1.5*beta + beta**2/6.) )

 xm1 = x0 * (alpha - alpha**2*beta/(1.-alpha) ) 

 if ( x0 > nu/mu )  then 
   print*, ' Epidemics spreads'
   print*, ' alpha ', alpha
   print*, ' x(t-->inf) = ', xm
   open (2,file='ep_spread.dat')   ! Trajectory file
 else
   print*, ' Epidemics does not spread'
   print*, ' alpha ', alpha
   print*, ' x(t-->inf) = ', xm
   open (2,file='ep_nospread.dat') ! Trajectory file
 end if

 x = (/ x0, y0 /)
!
!====================================================================
! Time integration
!
 do i=1, N
    call one_RK4_step(x)
    time = i * dt
    if ( mod(i,is) == 0) then 
        !print*, time, "   ", x(1), x(2)
        write (2,100)  time, x(1), x(2), x0+y0+z0-x(1)-x(2), xm
       end if
 end do
!
!====================================================================
 100 format(20(F16.8,3x)) 
!====================================================================
!
contains
!
!====================================================================

    subroutine one_RK4_step(x)   ! advances the time map one timestep
    use global 
    real (kind=kdp), dimension(2) :: x
    real (kind=kdp), dimension(2) :: k1, k2, k3, k4

    call Fvec(k1,x)
    call Fvec(k2,x + k1 * dt2)
    call Fvec(k3,x + k2 * dt2)
    call Fvec(k4,x + k3 * dt)
    x = x + dt6*(k1 + 2.d0*k2 + 2.d0*k3 + k4)

    end subroutine	
!
!--------------------------------------------------------------------
!    
    subroutine Fvec(F,x)   
    use global 
    real (kind=kdp), dimension(2) :: x, F

    F(1) = -mu * x(1) * x(2)                        
    F(2) =  mu * x(1) * x(2) - nu * x(2) 

    end subroutine
!
!====================================================================
!
end program  
