// PROGRAM sho
// 調和振動子

#include "TrueBASIC.h"
int main();
void initial(double *x0, double *v, double *w2, double *t, double *dt, 
             double *dt2, int *nshow);
void update(double *x, double *v, double w2, double *t, double dt, double dt2);
void output(double x, double v, double t);

int main()
{
    int counter, nshow;
    double dt, dt2, t, v, w2, x;

    initial(&x, &v, &w2, &t, &dt, &dt2, &nshow);
    counter = 0;
    do
    {
        update(&x, &v, w2, &t, dt, dt2);
        ++counter;
        if(counter == nshow)
        {
            output(x, v, t);
            counter = 0;
        }
    } while(!(kbhit()));

    return 0;
}

void initial(double *x0, double *v, double *w2, double *t, double *dt, 
             double *dt2, int *nshow)
{
    double show_time;
    char Stmp1_[_LBUFF_];

    printf("initial position = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", x0);              // メートル
    (*t) = 0;
    (*v) = 0;                               // 初速度
    // もとの (角) 振動数
    printf("ratio of k/m = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", w2);
    printf("time step = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", dt);
    (*dt2) = 0.5*(*dt);
    show_time = 0.1;                        // 出力の時間間隔
    (*nshow) = (int)(show_time/(*dt));
    printf("         time     position     velocity\n");
    printf("\n");                           // 空白行
}

void update(double *x, double *v, double w2, double *t, double dt, double dt2)
{
    double a, amid, vmid, xmid;

    // オイラー-リチャードソンのアルゴリズム
    a = -w2*(*x);
    vmid = (*v) + a*dt2;                    // 中点での速度
    xmid = (*x) + (*v)*dt2;                 // 中点の位置
    amid = -w2*xmid;                        // 中点での加速度
    (*v) += amid*dt;
    (*x) += vmid*dt;
    (*t) += dt;
}

void output(double x, double v, double t)
{
    printf("%13.4f %12.4f %12.4f\n", t, x, v);
}