/* Lotka-Volterra equation */
/* Runge-Kutta法 */
/* lotka_voltera_euler.c */
#include <stdio.h>
#include <math.h>
#define N 800
double fx(double x, double y);
double fy(double x, double y);
double alpha,beta,Gamma,delta;

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

 /* ここから初期条件 */
  alpha = 4.0,beta = 2.0,Gamma = 1.0,delta = 3.0;
  t0 = 0.0;
  x = 2.0;//2.0;
  y = 1.0;
 /* ここまで初期条件 */
  h = 20.0 / N;

  t = t0;
  fout1 = fopen("dat1.txt","w");
  fout2 = fopen("dat2.txt","w");
  fprintf(fout1, "%10.6f %10.6f\n",t,x);
  fprintf(fout2, "%10.6f %10.6f\n",t,y);
  kx2=ky2=kx3=ky3=kx4=ky4=0; //dummy
  for(i = 1;i < N;i++){
    x0 = x;
    y0 = y;
    t += h;
    /*  begin Euler */ 
    //     x = x0; //dummy
    //     y = y0; //dummy
    //     kx1 = fx(x,y);
    //     ky1 = fy(x,y);
    //     x = x0 + h * kx1;
    //     y = y0 + h * ky1;
     /*   end Euler */
    /* begin Huen */
    kx1 = fx(x0,y0); 
    ky1 = fy(x0,y0); 
    x = x0 + h * kx1;
    y = y0 + h * ky1;
    kx2 = fx(x,y);
    ky2 = fy(x,y);
    x = x0 + .5 * h * (kx1 + kx2);
    y = y0 + .5 * h * (ky1 + ky2);
    /* end Heun */
    /* 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);
        x = x0 + (1.0/6.0) * h * (kx1 + 2 * kx2 + 2 * kx3 + kx4);
        y = y0 + (1.0/6.0) * h * (ky1 + 2 * ky2 + 2 * ky3 + ky4);
	end runge-kutta */
    fprintf(fout1, "%10.6f %10.6f\n",t,x);
    fprintf(fout2, "%10.6f %10.6f\n",t,y);
  }

  fclose(fout1);
  fclose(fout2);

  fin1 = fopen("dat1.txt","r");
  fin2 = fopen("dat2.txt","r");

  gp = popen("gnuplot -persist","w");
  fprintf(gp, "set zeroaxis\n");
  fprintf(gp, "set xrange [%f:%f] \n",t0,t0+N*h);
  fprintf(gp, "set yrange [%f:%f] \n",0.0,10.0);
  //  fprintf(gp, "plot 'dat1.txt' with lines linetype 2, 'dat2.txt' with lines linetype 3\n");
    fprintf(gp, "plot 'dat1.txt' with lines linetype -1, 'dat2.txt' with lines linetype 1\n");
  pclose(gp);
  fclose(fin1);
  fclose(fin2);

  return 0;
}

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

double fy(double x, double y)
{
  return (- Gamma + delta * x) * y;
}