/* 
************************************************************************
*  ルンゲ・クッタ法による1次元運動方程式の解法                         *
*  投げ上げ                                                            *
************************************************************************
*/
       
#include <stdio.h>
#include <math.h>

#define N 2                                   /* number of equations */
 
main()
{
   void rungen(double t, double x[], double dt);
   double f(double t, double x[], int i);
   
   double t, x[N], v0, dt, h, xreal, xerror;
   int j;
   
   FILE *output;
   output = fopen("fall.dat","w");

   printf("Input: v0\n");
   scanf("%lf", &v0);
   printf("v0=%g\n", v0);
   printf("Input: dt\n");
   scanf("%lf", &dt);
   printf("dt=%g\n", dt);
   
   t = 0;
   x[0] = 0.0;                                /* initial position  */
   x[1] = v0;                                /* initial velocity  */
   
   fprintf(output, "%f\t%f\n", t, x[0]);

   while(x[0] >= 0)
   {
      rungen(t, x, dt);
      t += dt;    
      xreal=v0*t-0.5*9.8*t*t;
      xerror=x[0]-xreal;      
      fprintf(output, "%f\t%f\t%f\t%e\n", t, x[0], xreal, xerror);
   }
   printf("data stored in fall.dat\n");
   fclose(output);
}
/*-----------------------end of main program--------------------------*/

/* Runge-Kutta subroutine */
void rungen(double x, double y[], double step)
{
   double f(double x, double y[], int i);
   double h=step/2.0,                         /* the midpoint */
          t1[N], t2[N], t3[N],                /* temporary storage */
          k1[N], k2[N], k3[N],k4[N];          /* for Runge-Kutta  */
   int i;
 
   for (i=0; i<N; i++) t1[i] = y[i]+0.5*(k1[i]=step*f(x, y, i));
   for (i=0; i<N; i++) t2[i] = y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
   for (i=0; i<N; i++) t3[i] = y[i]+    (k3[i]=step*f(x+h, t2, i));
   for (i=0; i<N; i++) k4[i] =                 step*f(x + step, t3, i);

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

/* definition of equations - this is the harmonic oscillator */
double  f(double x, double y[], int i)
{
   if (i == 0) return(y[1]);               /* RHS of first equation */
   if (i == 1) return(-9.8);              /* RHS of second equation */
}