/* SEIR model */
/* Naoki Osada 20200713 */
#include <stdio.h>
#include <math.h>
#define IMAX 300
#define N 1000000
double fx(double x, double y, double z);
double fy(double x, double y, double z);
double fz(double x, double y, double z);
double beta,Gamma,epsilon;

int main(void)
{
  FILE *gp,*fin1,*fin2,*fout1,*fout2,*fin3,*fout3,*fin4,*fout4;
  int i;
  double h,t,t0;
  double x,y,z,x0,y0,z0,dx,dy,dz;
  double kx1,kx2,kx3,kx4,ky1,ky2,ky3,ky4,kz1,kz2,kz3,kz4;

  /* initial conditions */
  beta = 0.25,Gamma = 0.1; /* R_0=2.5 */
  epsilon = 0.2;
  t0 = 0.0;
  x = 1.0;
  y = 0;
  z = (double) 1.0/N;
  /**********************/
  h = 1.0;

  t = t0;
  fout1 = fopen("dat1.txt","w");
  fout2 = fopen("dat2.txt","w");
  fout3 = fopen("dat3.txt","w");
  fout4 = fopen("dat4.txt","w");  
  fprintf(fout1, "%10.6f %10.6f\n",t,x);
  fprintf(fout2, "%10.6f %10.6f\n",t,y);
  fprintf(fout3, "%10.6f %10.6f\n",t,z);
  fprintf(fout4, "%10.6f %10.6f\n",t,epsilon*y);  
  for(i = 1;i <= IMAX;i++){
    x0 = x;
    y0 = y;
    z0 = z;
    t += h;
    /* begin runge-kutta */
    kx1 = fx(x0,y0,z0);
    ky1 = fy(x0,y0,z0);
    kz1 = fz(x0,y0,z0);
    x = x0 + .5 * h * kx1;
    y = y0 + .5 * h * ky1;
    z = z0 + .5 * h * kz1;
    kx2 = fx(x,y,z);
    ky2 = fy(x,y,z);
    kz2 = fz(x,y,z);
    x = x0 + .5 * h * kx2;
    y = y0 + .5 * h * ky2;
    z = z0 + .5 * h * kz2;
    kx3 = fx(x,y,z);
    ky3 = fy(x,y,z);
    kz3 = fz(x,y,z);
    x = x0 + h * kx3;
    y = y0 + h * ky3;
    z = z0 + h * kz3;
    kx4 = fx(x,y,z);
    ky4 = fy(x,y,z);
    kz4 = fz(x,y,z);
    dx = (1.0/6.0) * h * (kx1 + 2 * kx2 + 2 * kx3 + kx4);
    dy = (1.0/6.0) * h * (ky1 + 2 * ky2 + 2 * ky3 + ky4);
    dz = (1.0/6.0) * h * (kz1 + 2 * kz2 + 2 * kz3 + kz4);
    x = x0 + dx;
    y = y0 + dy;
    z = z0 + dz;
    /* end runge-kutta */
    fprintf(fout1, "%10.6f %10.6f\n",t,x);
    fprintf(fout2, "%10.6f %10.6f\n",t,y);
    fprintf(fout3, "%10.6f %10.6f\n",t,z);
    fprintf(fout4, "%10.6f %10.6f\n",t,epsilon*y);
  }

  fclose(fout1);
  fclose(fout2);
  fclose(fout3);
  fclose(fout4);

  fin1 = fopen("dat1.txt","r");
  fin2 = fopen("dat2.txt","r");
  fin3 = fopen("dat3.txt","r");
  fin4 = fopen("dat4.txt","r");  

  /*  gp = popen("../gnuplot.exe -persist","w"); */
  gp = popen("gnuplot -persist","w");
  fprintf(gp, "set zeroaxis\n");
  fprintf(gp, "set xrange [%f:%f] \n",0.0,(float)IMAX);
  fprintf(gp, "set yrange [%f:%f] \n",0.0,1.0);
  fprintf(gp, "plot 'dat1.txt' title 'S' with lines linetype 2, 'dat2.txt' title 'E' with lines linetype 3, 'dat3.txt' title 'I' with lines linetype 4, 'dat4.txt' title 'newly infected' with lines linetype 5\n"); 
  pclose(gp);
  fclose(fin1);
  fclose(fin2);
  fclose(fin3);
  fclose(fin4);
  return 0;
}

double fx(double x, double y, double z)/* Susceptible */
{
  return - beta * x * z;
}

double fy(double x, double y, double z)/* Exposed */
{
  return beta * x * z - epsilon * y;
}

double fz(double x, double y, double z)/* Infectious */
{
  return epsilon * y - Gamma * z;
}