#include <stdio.h>
#include <math.h>

// 離心率
#define e 0.5
// 刻み幅
#define h 0.0001
//define h 0.00125663706143592
// 総ステップ数
#define N 10000 //20000 // 10000
// 変数ｘの要素数 (x,y,u,vで４つ)
#define VAR_N 4

void dfunc(double*,double*);
void rk4(int,double*,void (*)(double*,double*));
double E(double*);

int main(void){
  double x[2*VAR_N];
  x[0] = 1.0-e             ;//x
  x[1] = 0.0               ;//y
  x[2] = 0.0               ;//u
  x[3] = sqrt((1.0+e)/x[0]);//v
  //時間t
  double t=0.0;
  //カウンタ
  int i;
  for(i=0;i<N;i++){
    printf("%e %e %e %e %e %e\n",t,x[0],x[1],x[2],x[3],E(x));
    rk4(VAR_N,x,dfunc);
    t += h;
  }
  printf("%e %e %e %e %e %e\n",t,x[0],x[1],x[2],x[3],E(x));
  return 0;
}

// 微係数
void dfunc(double* x,double* dx){
  double a = - pow(x[0]*x[0]+x[1]*x[1],-1.5);
  dx[0]=h* x[2];   // dx = u dt
  dx[1]=h* x[3];   // dy = v dt
  dx[2]=h* a*x[0]; // du = -x/r^3 dt
  dx[3]=h* a*x[1]; // dv = -y/r^3 dt
}

/* ４次ルンゲ＝クッタ法 */
/* n: 与える変数の個数（配列ｘの要素数）               */
/* x: 解きたい微分方程式の変数の配列                   */
/* f: （微分方程式から）微係数を求める関数へのポインタ */
/* return: なし。ただし、xの各値を書き換える。         */
/* なお、刻み幅ｈは固定（定数）                        */
void rk4(int n,double *x,void (*f)(double*,double*)){
  double buf[n];
  double g[4][n];
  int i;
  (*f)(x,g[0]);
  for(i=0;i<n;i++)
    buf[i] = x[i] + 0.5*g[0][i] ;
  (*f)(buf,g[1]);
  for(i=0;i<n;i++)
    buf[i] = x[i] + 0.5*g[1][i] ;
  (*f)(buf,g[2]);
  for(i=0;i<n;i++)
    buf[i] = x[i] + g[2][i] ;
  (*f)(buf,g[3]);
  for(i=0;i<n;i++)
    x[i] += (g[0][i]+2.0*(g[1][i]+g[2][i])+g[3][i]) / 6.0 ;
  return;
}

double E(double* arr){
  double x=arr[0];
  double y=arr[1];
  double u=arr[2];
  double v=arr[3];
  return 0.5*(u*u+v*v)-pow(x*x+y*y,-0.5);
}