// PROGRAM oscillators
// 1次元の線形連成振動子系のシミュレーション

#include "TrueBASIC.h"
#define loop 4
int main();
void initial(int *L, double u[], double v[], double *t, double *dt, double usave[]);
void update(int L, double u[], double v[], double *ke, double *t, double dt, double usave[]);
void animate(int L, double u[], double ke, double t, double usave[]);

int main()
{
    int L, i;
    double dt, ke, t, u[22], v[22], usave[22];

    GWopen(0);

    initial(&L, u, v, &t, &dt, usave);
    do
    {
        for(i = 0; i < loop; i++)
            update(L, u, v, &ke, &t, dt, usave);
        animate(L, u, ke, t, usave);
    } while(!TBkeyinput());

    GWquit();
    return 0;
}

void initial(int *L, double u[], double v[], double *t, double *dt, double usave[])
{
    int j;
    int DP_ = 0;
    double DATA_[] = {
        0,0.5,0,0,0,0,0,0,0,0
      , 0,0,0,0,0,0,0,0,0,0
    };

    (*t) = 0;
    (*dt) = 0.025;
    (*L) = 10;                              // 粒子の数
    GWindow(-1.0f, -4.0f, (float)((*L)+1), 4.0f);
    GWline(0.5f, 0.0f, (*L)+0.5f, 0.0f);
    GWanchor(1);
    GWsetogn(0, 1);
    for(j = 1; j <= (*L); ++j)
    {
        u[j] = DATA_[DP_++];                // 初期変位
        GWsrect((float)(j), (float)(u[j]), 
                (float)(j+0.2), (float)(u[j]+0.2), RED);
    }
    for(j = 1; j <= (*L); ++j)
        v[j] = DATA_[DP_++];                // 初期速度
    u[0] = 0;                               // 固定壁の条件
    u[(*L) + 1] = 0;
    for(j = 0; j <= 21; ++j)
        usave[j] = u[j];
}

void update(int L, double u[], double v[], double *ke, double *t, double dt, double usave[])
{
    int j;
    double a[21], amid[21], umid[22], vmid[21];

    // オイラー-リチャードソンのアルゴリズム
    (*ke) = 0;
    // K/M = 1 としてある
    for(j = 1; j <= L; ++j)
    {
        usave[j] = u[j];
        a[j] = -2*u[j] + u[j + 1] + u[j - 1];
        umid[j] = u[j] + 0.5*v[j]*dt;
        vmid[j] = v[j] + 0.5*a[j]*dt;
    }
    umid[0] = 0;
    umid[L + 1] = 0;
    for(j = 1; j <= L; ++j)
    {
        amid[j] = -2*umid[j] + umid[j + 1] + umid[j - 1];
        u[j] = u[j] + vmid[j]*dt;
        v[j] = v[j] + amid[j]*dt;
        (*ke) += v[j]*v[j];
    }
    (*t) += dt;
}

void animate(int L, double u[], double ke, double t, double usave[])
{
    int j;
    double E, pe;
    char Stmp1_[_LBUFF_];

    pe = pow(u[1] - u[0], 2.0);             // 左のばねの相互作用
    GWsetogn(0, -1);
    for(j = 1; j <= L; ++j)
    {
        // 位置エネルギーの計算
        pe += pow(u[j + 1] - u[j], 2.0);
        // 横振動
        GWsrect((float)(j), (float)(usave[j]), 
                (float)(j+0.2), (float)(usave[j]+0.2), RED);
    }
    GWflush(-1);
    GWsetogn(0, -2);
    sprintf(Stmp1_, "t = %6.2f", t);
    GWputtxt(0.0f, 3.7f, Stmp1_);
    E = 0.5*(ke + pe);
    sprintf(Stmp1_, "E = %6.4f", E);
    GWputtxt(0.0f, 3.4f, Stmp1_);
    GWflush(-2);
}