/* ****************************************************************

This is the C version of program 2.6 from page 41 of
"Modeling Infectious Disease in humans and animals"
by Keeling & Rohani.

It is the SEIR epidemic with equal births and deaths.

This code is written to be simple, transparent and readily compiled.
Far more elegant and efficient code can be written.

This code can be compiled using gnu or intel C compilers:
icc -o Program_2_6  Program_2_6.c  -lm
g++ -o Program_2_6  Program_2_6.c  -lm

If you have any difficulties, we suggest you set DEBUG to 1 or -1

******************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

// Set DEBUG to 1 to output some debugging information
// Set DEBUG to -1 to output more debugging information
#define DEBUG 0

// Set up basic parameters
double beta=520.0/365.0;
double sigma=1.0/14.0;
double gamm=1.0/7.0;
double mu=1.0/(70.0*365.0);
double S0=0.1;
double E0=1e-4;
double I0=1e-4;

double MaxTime=60*365;

// Set up variables and rates of change
// Note we no-longer explicitly model R.
double t,S,E,I,R,Pop[3];
double dPop[3];

// Initialise the equations and Runge-Kutta integration
void Diff(double Pop[3])
{
  // Set up temporary variables to make the equations look neater

  double tmpS, tmpE, tmpI;

  tmpS=Pop[0]; tmpE=Pop[1]; tmpI=Pop[2];

  /* The differential equations */

  dPop[0] = mu - beta*tmpS*tmpI - mu*tmpS;              // dS/dt
  dPop[1] = beta*tmpS*tmpI - sigma*tmpE -mu*tmpE;       // dE/dt
  dPop[2] = sigma*tmpE - gamm*tmpI - mu*tmpI;           // dI/dt

  return;
}

void Runge_Kutta(double);
void Read_in_parameters(FILE *);
void Perform_Checks();
void Output_Data(FILE *);


// The main program

int main(int argc, char** argv)
{
  char FileName[1000];
  double step, Every;
  FILE *in, *out;
   
  /* Find the given parameter  file name and open the file 
     otherwise use parameters set at top of program */
  if(argc>1)
    {
      strcpy(FileName, *(argv+1));  
      
      if((in = fopen(FileName,"r"))==NULL)
	{
	  printf("Cannot read file %s\n",FileName);
	  exit(1);
	}

      /* Read in the parameters */
      
      Read_in_parameters(in);

      if(DEBUG) printf("\nReading parameters from file %s\n",FileName);
    }
  else
    {
      if(DEBUG) printf("\nUsing default parameters\n");
    }

  if(DEBUG) printf("beta=%g\nsigma=%g\ngamma=%g\nmu=%g\nInitial S=%g\nInitial E=%g\n\nInitial I=%g\nRuns until time %g\n\n",beta,sigma,gamm,mu,S0,E0,I0,MaxTime);


  /* Check all parameters are OK & set up intitial conditions */

  Perform_Checks();
 
  S=S0; E= E0; I=I0; R=1-S-E-I;

  /* Find a suitable time-scale for outputs */

  step=0.01/((beta+sigma+gamm+mu)*S0);
  
  Every=pow(10,floor(log10((1.0/((beta+sigma+gamm+mu)*S0)))));
  while(MaxTime/Every>10000)
    {
      Every*=10.0;
    }
  

  if(DEBUG) printf("Using a time step of %g and outputing data every %g\n\n",step,Every);
    


  if((out = fopen("Output","w"))==NULL)
    {
      printf("Cannot open file Output to write data\n");
      exit(1);
    }

  /* The main iteration routine */

  t=0;

  Output_Data(out);

  do
    {
      Runge_Kutta(step);  t+=step;
      
      /* If time has moved on sufficiently, output the current data */
      if( floor(t/Every) > floor((t-step)/Every) )
	{
	  Output_Data(out);
	}
      else
	{
	  if(DEBUG<0) Output_Data(stdout);
	}
    }
  while(t<MaxTime);

  Output_Data(out);

  fclose(out);
}






void Read_in_parameters(FILE *in)
{
  char str[200];
  
  fscanf(in,"%s",str);
  fscanf(in,"%lf",&beta);
  
  fscanf(in,"%s",str);
  fscanf(in,"%lf",&sigma);
  
  fscanf(in,"%s",str);
  fscanf(in,"%lf",&gamm);
  
  fscanf(in,"%s",str);
  fscanf(in,"%lf",&mu);
  
  fscanf(in,"%s",str);
  fscanf(in,"%lf",&S0);
  
  fscanf(in,"%s",str);
  fscanf(in,"%lf",&I0);
  
  fscanf(in,"%s",str);
  fscanf(in,"%lf",&E0);
  
  fscanf(in,"%s",str);
  fscanf(in,"%lf",&MaxTime);
  
  fclose(in);
}

void Perform_Checks()
{
   if(S0<=0) {printf("ERROR: Initial level of susceptibles (%g) is less than or equal to zero\n",S0); exit(1);}

  if(I0<=0) {printf("ERROR: Initial level of exposed (%g) is less than or equal to zero\n",E0); exit(1);}

  if(I0<=0) {printf("ERROR: Initial level of infectious (%g) is less than or equal to zero\n",I0); exit(1);}

  if(beta<=0) {printf("ERROR: Transmission rate beta (%g) is less than or equal to zero\n",beta); exit(1);}

  if(sigma<=0) {printf("ERROR: Exposed to infectious rate sigma (%g) is less than or equal to zero\n",sigma); exit(1);}

  if(gamm<=0) {printf("ERROR: Recovery rate gamma (%g) is less than or equal to zero\n",gamm); exit(1);}

  if(mu<0) {printf("ERROR: Birth / Death rate mu (%g) is less than zero\n",mu); exit(1);}

  if(MaxTime<=0) {printf("ERROR: Maximum run time (%g) is less than or equal to zero\n",MaxTime); exit(1);}

  if(S0+E0+I0>1) {printf("WARNING: Initial level of susceptibles+infecteds (%g+%g+%g=%g) is greater than one\n",S0,E0,I0,S0+E0+I0);}

  if(beta*sigma<(gamm+mu)*(sigma+mu)) {printf("WARNING: Basic reproductive ratio (R_0=%g) is less than one\n",beta*sigma/((gamm+mu)*(sigma+mu)));}

}



void Output_Data(FILE *out)
{
  if(out!=stdout) fprintf(out,"%g \t %g %g %g\n",t,S,E,I);

  if(DEBUG) printf("%g \t %g %g %g\n",t,S,E,I);
    
}




void Runge_Kutta(double step)
{
  int i;
  double dPop1[3], dPop2[3], dPop3[3], dPop4[3];
  double tmpPop[3], initialPop[3];

  /* Integrates the equations one step, using Runge-Kutta 4
     Note: we work with arrays rather than variables to make the
     coding easier */

  initialPop[0]=S; initialPop[1]=E; initialPop[2]=I;

  Diff(initialPop);
  for(i=0;i<3;i++)
    {
      dPop1[i]=dPop[i];
      tmpPop[i]=initialPop[i]+step*dPop1[i]/2;
    }

  Diff(tmpPop);
  for(i=0;i<3;i++)
    {
      dPop2[i]=dPop[i];
      tmpPop[i]=initialPop[i]+step*dPop2[i]/2;  
    }

  Diff(tmpPop);
  for(i=0;i<3;i++)
    {
      dPop3[i]=dPop[i];
      tmpPop[i]=initialPop[i]+step*dPop3[i]; 
    }

  Diff(tmpPop);

  for(i=0;i<3;i++)
    {
      dPop4[i]=dPop[i];

      tmpPop[i]=initialPop[i]+(dPop1[i]/6 + dPop2[i]/3 + dPop3[i]/3 + dPop4[i]/6)*step;
    }


  S=tmpPop[0]; E=tmpPop[1]; I=tmpPop[2]; R=1-S-E-I;

  return;
}