/* 
************************************************************************
*  ルンゲ・クッタ法による1次元運動方程式の解法                         *
*  速度に比例する抵抗がある調和振動子                                  *
************************************************************************
  t：時間、x：位置、v：速度、dt：時間の刻み幅
  t0、x0、v0：初期値
  tlim：時間の最大値
  alpha：速度に比例する抵抗の係数

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

double alpha;
 
main()
{
   double f(double t, double x, double v, int i);
   
   double t, x, v, x0, v0, dt, h;
   double pi,tlim;
   int j;
   double t10, t11, t20, t21, t30, t31,       /* Runge-Kutta法のための */
     k10, k11, k20, k21, k30, k31 ,k40, k41;  /* 一時変数     */
   
   FILE *output;
   output = fopen("harm.dat","w");            /* harm.datにデータをセーブ */

/* 速度に比例する抵抗の係数 */
   printf("Input: alpha\n");
   scanf("%lf", &alpha);
   printf("alpha=%g\n", alpha);
/* 時間の刻み幅 */
   printf("Input: dt\n");
   scanf("%lf", &dt);
   printf("dt=%g\n", dt);

   
   t = 0;
   x0 = 1.0;                                /* 初期位置  */
   v0 = 0.0;                                /* 初期速度  */
   x = x0;
   v = v0; 
   
   fprintf(output, "%f\t%f%f\n", t, x, v);

      pi=3.1415926535897932385;
      tlim=6.0*pi;

   while(t <= tlim)
   {

      k10 = dt*f(t, x, v, 0);
      t10 = x+0.5*k10;
      k11 = dt*f(t, x, v, 1);
      t11 = v+0.5*k11;
      
      k20 = dt*f(t+h, t10, t11, 0);
      t20 = x+0.5*k20;
      k21 = dt*f(t+h, t10, t11, 1);
      t21 = v+0.5*k21;
      
      k30 = dt*f(t+h, t20, t21, 0);
      t30 = x+    k30;
      k31 = dt*f(t+h, t20, t21, 1);
      t31 = v+    k31;
      
      k40 = dt*f(t + dt, t30, t31, 0);
      k41 = dt*f(t + dt, t30, t31, 1);
      
      x += (k10+2*k20+2*k30+k40)/6.0;
      v += (k11+2*k21+2*k31+k41)/6.0;

      t += dt;    
      fprintf(output, "%f\t%f\t%f\t\n", t, x, v);
   }
   printf("data stored in harm.dat\n");
   fclose(output);
}
/*--------------------------------------------------------------------*/

/* 微分方程式の右辺 */
double  f(double t, double x, double v, int i)
{
   if (i == 0) return(v);                  /* dx/dt*/
   if (i == 1) return(-x-alpha*v);         /* dv/dt*/}