/* SEIAR model */
/* Naoki Osada 20200723 */
#include <stdio.h>
#include <math.h>
#define IMAX 300
#define N 1000000
double fx(double x, double y, double z, double u);
double fy(double x, double y, double z, double u);
double fz(double x, double y, double z, double u);
double fu(double x, double y, double z, double u);
double alpha,beta,Gamma,delta,epsilon,zeta;

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

  /* initial conditions */
  alpha = 2.0, beta = 0.25, Gamma = 0.1;
  delta = 0.1, epsilon = 0.2, zeta = 0.6;
  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");
  fout5 = fopen("dat5.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,u);
  fprintf(fout5, "%10.6f %10.6f\n",t,epsilon*y);  
  for(i = 1;i <= IMAX*1.4;i++){
    x0 = x;
    y0 = y;
    z0 = z;
    u0 = u;    
    t += h;
    /* begin runge-kutta */
    kx1 = fx(x0,y0,z0,u0);
    ky1 = fy(x0,y0,z0,u0);
    kz1 = fz(x0,y0,z0,u0);
    ku1 = fz(x0,y0,z0,u0);
    x = x0 + .5 * h * kx1;
    y = y0 + .5 * h * ky1;
    z = z0 + .5 * h * kz1;
    u = u0 + .5 * h * ku1;
    kx2 = fx(x,y,z,u);
    ky2 = fy(x,y,z,u);
    kz2 = fz(x,y,z,u);
    ku2 = fu(x,y,z,u);
    x = x0 + .5 * h * kx2;
    y = y0 + .5 * h * ky2;
    z = z0 + .5 * h * kz2;
    u = u0 + .5 * h * ku2;
    kx3 = fx(x,y,z,u);
    ky3 = fy(x,y,z,u);
    kz3 = fz(x,y,z,u);
    ku3 = fu(x,y,z,u);
    x = x0 + h * kx3;
    y = y0 + h * ky3;
    z = z0 + h * kz3;
    u = u0 + h * ku3;
    kx4 = fx(x,y,z,u);
    ky4 = fy(x,y,z,u);
    kz4 = fz(x,y,z,u);
    ku4 = fu(x,y,z,u);
    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);
    du = (1.0/6.0) * h * (ku1 + 2 * ku2 + 2 * ku3 + ku4);
    x = x0 + dx;
    y = y0 + dy;
    z = z0 + dz;
    u = u0 + du;
    /* 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,u);
    fprintf(fout5, "%10.6f %10.6f\n",t,epsilon*y0);
    printf("%4d %10.2f %10.2f %10.2f %10.2f %10.2f\n",i,N*x,N*y,N*z,N*u,N*epsilon*y0);
  }
  fclose(fout1);
  fclose(fout2);
  fclose(fout3);
  fclose(fout4);
  fclose(fout5);

  fin1 = fopen("dat1.txt","r");
  fin2 = fopen("dat2.txt","r");
  fin3 = fopen("dat3.txt","r");
  fin4 = fopen("dat4.txt","r");
  fin5 = fopen("dat5.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 'A' with lines linetype 5, 'dat5.txt' title 'newly infected' with lines linetype 6\n"); 
  pclose(gp);
  fclose(fin1);
  fclose(fin2);
  fclose(fin3);
  fclose(fin4);
  return 0;
}

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

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

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

double fu(double x, double y, double z, double u)/* Asymptomatic */
{
  return epsilon * (1.0 - zeta) * y - delta * u;
}