/*--------------------------------------------*/
/*　4次ルンゲクッタ法の復習(2変数・調和振動子)*/
/*--------------------------------------------*/
#include<stdio.h>
#include <stdlib.h>
#include<math.h>

#define max_it     1000000     /*　最大時間ステップ数　*/
#define t_end      30e0         /*  計算終了の t */
void func(double t, double x[], double f[]);

main(void)
{
  int i,num;
  double t,x[2],f[2],k1[2],k2[2],k3[2],k4[2],f_fin[2],dt;
  double x_tmp[2];
  FILE *output_file;
  
  output_file=fopen("output.dat","w"); /* 出力ファイルオープン*/ 

  /* ----------------  初期条件 ----------------- */

  x[0]=0e0;
  x[1]=1e0;
  dt=1e-3;
  t=0e0;

  /* -------------- 初期条件終わり ------------- */


  
  /* -------- 4次ルンゲクッタの時間積分 --------- */    

  for(num=0; num < max_it; num++){

    func(t,x,f);
    for(i=0; i < 2 ; i++){
      k1[i]=f[i]*dt; 
      x_tmp[i]=x[i]+k1[i]*5e-1;
    }

    func(t+dt*5e-1,x_tmp,f);
    for(i=0; i < 2 ; i++){
      k2[i]=f[i]*dt; 
      x_tmp[i]=x[i]+k2[i]*5e-1;
    }

    func(t+dt*5e-1,x_tmp,f);
    for(i=0; i < 2 ; i++){
      k3[i]=f[i]*dt; 
      x_tmp[i]=x[i]+k3[i];
    }

    func(t+dt,x_tmp,f);
    for(i=0; i < 2 ; i++){
      k4[i]=f[i]*dt; 
    }

    for(i=0; i < 2 ; i++){
    x[i] = x[i] + 1e0 / 6e0 * ( k1[i] + 2e0*k2[i] + 2e0*k3[i] + k4[i] );
    f_fin[i] = 1e0 / 6e0 * ( k1[i] + 2e0*k2[i] + 2e0*k3[i] + k4[i] ) / dt;
    }

    t=t+dt;
    
    /* --------            出力ルーチン              -----------*/

    if(num%2==0){   /* 2 step ごとに出力 */
      fprintf(output_file,"%f %f %f \n",(float)t,(float)x[0],(float)x[1]);
    }

    /* --------        出力ルーチン終わり              -----------*/

    /* --------        time step control ---------- */

    dt =1e-1;

    for(i=0; i < 2 ; i++){
      if( abs(f_fin[i]) > 1e-4 && fabs(x[i]) > 1e-4){
	if( dt  > 1e-1 * fabs(x[i] / f_fin[i])){
	  dt =  1e-1 * fabs(x[i] / f_fin[i]) ;
	}
      }
    }

    /* --------        time step control 終わり ---------- */

    printf("t = %f\n",t);

    if(t > t_end) break; /* 計算終了判定 */
  }

}

void func(double t, double x[], double f[])
{

  f[0] = x[1];
  f[1] = -x[0];
}