// PROGRAM planet
// 惑星の運動

#include "TrueBASIC.h"
#include "csgraphics.h"

int main();
void initial(double pos[], double vel[], double *t, double *GM, double *dt, 
             int *nshow);
void Euler(double pos[], double vel[], double *t, double GM, double dt);
void output(double pos[], double t);

int main()
{
    int counter, nshow;
    double GM, dt, t;
    double pos[3], vel[3];                  // ２つの配列の次元を定義する

    GWopen(0);

    initial(pos, vel, &t, &GM, &dt, &nshow);
    output(pos, t);
    counter = 0;
    do
    {
        Euler(pos, vel, &t, GM, dt);
        if((++counter % nshow) == 0)
            output(pos, t);                 // "地球"の位置を描く
    } while(!TBkeyinput());

    GWquit();
    return 0;
}

void initial(double pos[], double vel[], double *t, double *GM, double *dt, 
             int *nshow)
{
    float xwin, ywin;
    char Stmp1_[_LBUFF_];

    (*GM) = 4.0*(pi*pi);                    // 定数 GM (天文単位)
    (*t) = 0;
    printf("time step = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", dt);
    printf("number of time steps between plots = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", nshow);
    printf("initial x position = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", &pos[1]);
    pos[2] = 0;                             // 初期位置の y 座標
    vel[1] = 0;                             // 初速度の x 成分
    printf("initial y-velocity = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", &vel[2]);
    // 予想される半長軸の最大値
    compute_aspect_ratio((float)(2*pos[1]), &xwin, &ywin);  // ライブラリ csgraphics の中にある
    GWindow(-xwin, -ywin, xwin, ywin);
}

void Euler(double pos[], double vel[], double *t, double GM, double dt)
// オイラー-クロマー・アルゴリズム
{
    int i;
    double accel[3], r2, r3;

    r2 = pos[1]*pos[1] + pos[2]*pos[2];
    r3 = r2*sqrt(r2);
    for(i = 1; i <= 2; ++i)
    {
        accel[i] = -GM*pos[i]/r3;
        vel[i] += accel[i]*dt;
        pos[i] += vel[i]*dt;
    }
    (*t) += dt;
}

void output(double pos[], double t)
{
    double radius;

    if(t == 0)
    {
        GWclear(-1);
        radius = 0.05;                      // "太陽"の半径(比例してない)
        GWsetpen(RED, -1, -1, 4);
        GWsetmrk(6, (float)(2*radius), RED, -1, 4);
        GWputmrk(0.0f, 0.0f);   // 原点にある"太陽"
    }
    GWsetpxl((float)(pos[1]), (float)(pos[2]), BLUE);
}