/* 
************************************************************************
*  ルンゲ・クッタ法による2次元運動方程式の解法                         *
*  速度の１乗および２乗に比例する抵抗がある斜方投射                                    *
************************************************************************
  t：時間、x[0]=x、x[1]=y：位置、x[2]=vx、x[3]=vy：速度、dt：時間の刻み幅
  t0、x0、y0、vx0、vy0：初期値
  v0：初速度、theta：投げ上げ角度
  alpha：速度に比例する抵抗の係数
  beta：速度の２乗に比例する抵抗の係数

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

#define N 4                                   /* 方程式の数 */

double alpha, beta;
 
main()
{
   void rungen(double t, double x[], double dt);
   double f(double t, double x[], int i);
   
   double t, x[N], v0, vx0, vy0, theta, dt, h;
   double rho, rball, mass, eta;
   double pi;
   int j;
   
/* 出力ファイルの指定 */
   FILE *output;
   output = fopen("ball.dat","w");          /* ball.datにデータをセーブ */

/*  初速度 */
   printf("Input: v0 (km/h)\n");
   scanf("%lf", &v0);
   printf("v0=%g\n", v0);
/* 投げ上げ角度 */
   printf("Input: theta\n");
   scanf("%lf", &theta);
   printf("theta=%g\n", theta);
   pi=3.1415926535897932385;

/* 時間の刻み幅 */
   printf("Input: dt\n");
   scanf("%lf", &dt);
   printf("dt=%g\n", dt);

   rho=1.2; /* 空気の密度  */
   rball=0.0715/2.0; /* ボールの半径  */
   mass=0.1445;  /* ボールの質量  */
   eta=1.8e-5;  /* 空気の粘性係数 */
   v0=v0*1000.0/3600.0;
   vx0=v0*cos(pi*theta/180);
   vy0=v0*sin(pi*theta/180);

   alpha=6.0*pi*rball*eta/mass;
   beta=0.25*pi*rho*rball*rball/mass;

   
/* 初期値 */
   t = 0.0;
   x[0] = 0.0;                                /* 初期位置 */
   x[1] = 0.0;                                /* 初期位置 */
   x[2] = vx0;                                /* 初期速度 */
   x[3] = vy0;                                /* 初期速度 */
   
   fprintf(output, "%f\t%f\t%f\n", t, x[0], x[1]);

   while(x[1] >= 0.0)
   {
/* ルンゲ・クッタ法ルーチンを呼び出す */
      rungen(t, x, dt);
      t += dt;    
      fprintf(output, "%f\t%f\t%f\n", t, x[0], x[1]);
   }
   printf("data stored in ball.dat\n");
   fclose(output);
}
/*-----------------------end of main program--------------------------*/

void rungen(double t, double x[], double dt)
{
   double f(double t, double x[], int i);
   double h=dt/2.0,
          t1[N], t2[N], t3[N],                /* Runge-Kutta法のための */
          k1[N], k2[N], k3[N],k4[N];          /* 一時変数              */
   int i;
 
   for (i=0; i<N; i++) 
     {
       k1[i] = dt*f(t, x, i);
       t1[i] = x[i]+0.5*k1[i];
     }
   for (i=0; i<N; i++) 
     {
       k2[i] = dt*f(t+h, t1, i);
       t2[i] = x[i]+0.5*k2[i];
     }
   for (i=0; i<N; i++) 
     {
       k3[i] = dt*f(t+h, t2, i);
       t3[i] = x[i]+    k3[i];
     }
   for (i=0; i<N; i++) 
     {
       k4[i] = dt*f(t + dt, t3, i);
     }

   for (i=0; i<N; i++) x[i] += (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
/*--------------------------------------------------------------------*/

/* 微分方程式の右辺 */
double  f(double t, double x[], int i)
{
  double vabs;                         /* 速度の絶対値 */
  vabs=sqrt(x[2]*x[2]+x[3]*x[3]);
   if (i == 0) return(x[2]);                            /* dx[0]/dx = dx/dt */
   if (i == 1) return(x[3]);                            /* dx[1]/dx = dy/dt */
   if (i == 2) return(-alpha*x[2] -beta*vabs*x[2]);     /* dx[2]/dx = dvx/dt */ 
   if (i == 3) return(-9.8-alpha*x[3] -beta*vabs*x[3]); /* dx[3]/dx = dvy/dt */
}