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

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

  /* initial conditions */
  beta = 0.25,Gamma = 0.1;
  t0 = 0.0;
  x = 1.0;
  y = 1.0/N; 
  dx = 0;
  h = 1.0;  
  /********************/
  
  t = t0;
  fout1 = fopen("dat1.txt","w");
  fout2 = fopen("dat2.txt","w");
  fout3 = fopen("dat3.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,-dx);  
  for(i = 1;i <= IMAX*1.6;i++){
    x0 = x;
    y0 = y;
    t += h;
    /* begin runge-kutta */
    kx1 = fx(x0,y0);
    ky1 = fy(x0,y0);
    x = x0 + .5 * h * kx1;
    y = y0 + .5 * h * ky1;
    kx2 = fx(x,y);
    ky2 = fy(x,y);
    x = x0 + .5 * h * kx2;
    y = y0 + .5 * h * ky2;
    kx3 = fx(x,y);
    ky3 = fy(x,y);
    x = x0 + h * kx3;
    y = y0 + h * ky3;
    kx4 = fx(x,y);
    ky4 = fy(x,y);
    dx = (1.0/6.0) * h * (kx1 + 2 * kx2 + 2 * kx3 + kx4);
    dy = (1.0/6.0) * h * (ky1 + 2 * ky2 + 2 * ky3 + ky4);
    x = x0 + dx;
    y = y0 + dy;
    /* 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,-dx);
  }

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

  fin1 = fopen("dat1.txt","r");
  fin2 = fopen("dat2.txt","r");
  fin3 = fopen("dat3.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' with lines linetype 2, 'dat2.txt' with lines linetype 3, 'dat3.txt' with lines linetype 4\n");
  fprintf(gp, "plot 'dat1.txt' title 'S' with lines linetype 2, 'dat2.txt' title 'I' with lines linetype 3, 'dat3.txt' title 'S(t-1)-S(t)' with lines linetype 4\n"); 
  pclose(gp);
  fclose(fin1);
  fclose(fin2);
  fclose(fin3);  

  return 0;
}

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

double fy(double x, double y)/* Infectious */
{
  return beta * y * x- Gamma * y;
}