// PROGRAM drag

#include "TrueBASIC.h"
int main();
void initial(double *y, double *y0, double *v, double *t0, double *g, 
             double *vt2, double *dt, double *dt_2);
void Euler_Richardson(double *y, double *v, double *t, double g, double vt2, 
                      double dt, double dt_2);
void print_table(double y, double y0, double t, double dt, int *nshow);

int main()
{
    int counter, nshow;
    double dt, dt_2, g, t, v, vt2, y, y0;

    // 抵抗力は v^2 に比例すると仮定
    initial(&y, &y0, &v, &t, &g, &vt2, &dt, &dt_2);
    print_table(y, y0, t, dt, &nshow);
    while(!(t >= 0))
        Euler_Richardson(&y, &v, &t, g, vt2, dt, dt_2);
    print_table(y, y0, t, dt, &nshow);      // t = 0 での値
    counter = 0;
    while(t <= 0.80)
    {
        Euler_Richardson(&y, &v, &t, g, vt2, dt, dt_2);
        if((++counter % nshow) == 0)
            print_table(y, y0, t, dt, &nshow);
    }

    return 0;
}

void initial(double *y, double *y0, double *v, double *t0, double *g, 
             double *vt2, double *dt, double *dt_2)
{
    double vt;
    char Stmp1_[_LBUFF_];

    (*t0) = -0.132;                         // 時刻の初期値 (表3.1を見よ)
    (*y0) = 0;
    (*y) = (*y0);
    (*v) = 0;                               // 初速度
    printf("terminal velocity (m/s) = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", &vt);
    (*vt2) = vt*vt;
    (*dt) = 0.001;
    (*dt_2) = 0.5*(*dt);
    (*g) = 9.8;
}

void Euler_Richardson(double *y, double *v, double *t, double g, double vt2, 
                      double dt, double dt_2)
// オイラー-リチャードソン法
{
    double a, amid, v2, vmid, vmid2, ymid;

    v2 = (*v)*(*v);
    a = -g*(1 - v2/vt2);
    vmid = (*v) + a*dt_2;                   // 中点での速度
    ymid = (*y) + (*v)*dt_2;                // 中点の位置
    vmid2 = vmid*vmid;
    amid = -g*(1 - vmid2/vt2);              // 中点での加速度
    (*v) += amid*dt;
    (*y) += vmid*dt;
    (*t) += dt;
}

void print_table(double y, double y0, double t, double dt, int *nshow)
{
    double distance_fallen, show_time;

    if(t < 0)
    {
        show_time = 0.1;                    // 出力の時間間隔
        // dt = 0.001 としてみるとよい
        (*nshow) = (int)(show_time/dt);
        printf("\n        time  displacement\n\n");
    }
    distance_fallen = y0 - y;
    printf("%12f %12f\n", t, distance_fallen);
}