/* 惑星運動のRunge-Kutta法による近似 */
#include <stdio.h>
#include <math.h>
#define N 800
void f(double * u);

int main(void){
  FILE *gp,*fin,*fout;
  double u[4],k1[4],k2[4],k3[4],k4[4],h;
  int i,j;

  h = 0.02;
  u[0] =3;
  u[1] = 0;
  u[2] = 0.3;
  u[3] = 0.2;

  fout = fopen("dat.txt","w");
  fprintf(fout, "%10.5f %10.5f \n",u[0],u[1]);
  for(i = 0;i <= N;i++){
    for(j = 0;j <= 3;j++)
      k1[j] = u[j];
    f(k1);
    for(j = 0;j <= 3;j++){
       k1[j] = h * k1[j];
       k2[j] = u[j] + 0.5 * k1[j];
    }
    f(k2);
    for(j = 0;j <= 3;j++){
      k2[j] = h * k2[j];
      k3[j] = u[j] + 0.5 * k2[j];
    }
    f(k3);
    for(j = 0;j <= 3;j++){
      k3[j] = h * k3[j];
      k4[j] = u[j] + k3[j];
    }
    f(k4);
    for(j = 0;j <= 3;j++)
      k4[j] = h * k4[j];
    for(j = 0;j <= 3;j++)
      u[j] = u[j] + (1.0/6.0) * (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]);
    fprintf(fout, "%10.5f %10.5f \n",u[0],u[1]);
  }
  fclose(fout);

  fin = fopen("dat.txt","r");
  gp = popen("gnuplot -persist","w");
  fprintf(gp, "set zeroaxis\n");
  fprintf(gp, "set xrange [-1:4] \n");
  fprintf(gp, "set yrange [-1:2] \n");
  fprintf(gp, "plot 'dat.txt' with lines linetype -1\n");
  pclose(gp);
  fclose(fin);
  
  return 0;
}

void f(double * u)
{
  double r2,u0,u1;

  u0 = u[0];
  u1 = u[1];
  r2 = pow(u0 * u0 + u1 * u1,1.5);
  u[0] = u[2];
  u[1] = u[3];
  u[2] = -u0 / r2;
  u[3] = -u1 / r2;

  return;
}