// PROGRAM md

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

int main();
void initial(double *t, double *ke, double *kecum, double *pecum, double *vcum, 
             double *area);
void set_up_windows(int IC1_, int IC2_);
void headings();
void Verlet(double *t, double *ke, double *pe, double *virial);
double pbc(double pos, double L);
void accel(double *pe, double *virial);
void force(double dx, double dy, double *fx, double *fy, double *pot);
double separation(double ds, double L);
void show_positions(char* flag, int IC2_);
void show_output(double t, double ke, double pe, double virial, double *kecum, 
                 double *vcum, int *ncum, double area, int IC1_);
void save_config();

double Lx;
double Ly;
int N;
double R2cum[101];
double ax[37];
double ay[37];
double dr;
double dt;
double dt2;
double gcum[1001];
double vx[37];
double vy[37];
double x[37];
double xsave[101];
double y[37];
double ysave[101];

int main()
{
    int IC1_ = 1, IC2_ = 2, ncum;
    double E, area, ke, kecum, pe, pecum, t, vcum, virial;
    char flag[_LBUFF_];

    GWopen(0);

    initial(&t, &ke, &kecum, &pecum, &vcum, &area);
    set_up_windows(IC1_, IC2_);
    accel(&pe, &virial);
    E = ke + pe;                            // 全エネルギー
    ncum = 0;                               // データ積算回数
    strcpy(flag, "");
    GWselvp(IC2_);
    do
    {
        show_positions(flag, IC2_);
        Verlet(&t, &ke, &pe, &virial);
        show_output(t, ke, pe, virial, &kecum, &vcum, &ncum, area, IC1_);
    } while(!(strcmp(flag, "stop") == 0));
    save_config();

    GWquit();
    return 0;
}

void initial(double *t, double *ke, double *kecum, double *pecum, double *vcum, 
             double *area)
{
    int i;
    char file[_LBUFF_], heading[_LBUFF_], response[_LBUFF_], start[_LBUFF_], 
         Stmp1_[_LBUFF_];
    FILE *FP1_;
    int DP_ = 0;
    double DATA_[] = {
        1.09,0.98,-0.33,0.78,3.12,5.25,0.12,-1.19
      , 0.08,2.38,-0.08,-0.10,0.54,4.08,-1.94,-0.56
      , 2.52,4.39,0.75,0.34,3.03,2.94,1.70,-1.08
      , 4.25,3.01,0.84,0.47,0.89,3.11,-1.04,0.06
      , 2.76,0.31,1.64,1.36,3.14,1.91,0.38,-1.24
      , 0.23,5.71,-1.58,0.55,1.91,2.46,-1.55,-0.16
      , 4.77,0.96,-0.23,-0.83,5.10,4.63,-0.31,0.65
      , 4.97,5.88,1.18,1.48,3.90,0.20,0.46,-0.51
    };

    dt = 0.01;
    dt2 = dt*dt;
    strcpy(response, "");
    do
    {
        printf("read data statements (d) or file (f)? ");
        fgets(Stmp1_, _LBUFF_, stdin);
        sscanf(Stmp1_, "%[^\n]", start);
        if((strcmp(start, "d") == 0) || (strcmp(start, "D") == 0))
        {
            strcpy(response, "ok");
            N = 16;
            Lx = 6;
            Ly = 6;
            for(i = 1; i <= N; ++i)
            {
                x[i] = DATA_[DP_++];
                y[i] = DATA_[DP_++];
                vx[i] = DATA_[DP_++];
                vy[i] = DATA_[DP_++];
            }
        }
        else if((strcmp(start, "f") == 0) || (strcmp(start, "F") == 0))
        {
            strcpy(response, "ok");
            printf("file name = ");
            fgets(Stmp1_, _LBUFF_, stdin);
            sscanf(Stmp1_, "%[^\n]", file);
            FP1_ = fopen(file, "r");
            fgets(Stmp1_, _LBUFF_, FP1_);
            sscanf(Stmp1_, "%d", &N);
            fgets(Stmp1_, _LBUFF_, FP1_);
            sscanf(Stmp1_, "%lg", &Lx);
            fgets(Stmp1_, _LBUFF_, FP1_);
            sscanf(Stmp1_, "%lg", &Ly);
            fgets(Stmp1_, _LBUFF_, FP1_);
            sscanf(Stmp1_, "%[^\n]", heading);
            for(i = 1; i <= N; ++i)
            {
                fgets(Stmp1_, _LBUFF_, FP1_);
                sscanf(Stmp1_, " %lg, %lg", &x[i],&y[i]);
            }
            fgets(Stmp1_, _LBUFF_, FP1_);
            sscanf(Stmp1_, "%[^\n]", heading);
            for(i = 1; i <= N; ++i)
            {
                fgets(Stmp1_, _LBUFF_, FP1_);
                sscanf(Stmp1_, " %lg, %lg", &vx[i],&vy[i]);
            }
            fclose(FP1_);
        }
        else
        {
            printf("\n");
            printf("d or f are the only acceptable responses.\n");
        }
    } while(!(strcmp(response, "ok") == 0));
    (*ke) = 0;                              // 運動エネルギー
    for(i = 1; i <= N; ++i)
        (*ke) += vx[i]*vx[i] + vy[i]*vy[i];
    (*ke) *= 0.5;
    (*area) = Lx*Ly;
    (*t) = 0;                               // 時刻
    // 和の初期化
    (*kecum) = 0;
    (*pecum) = 0;
    (*vcum) = 0;
}

void set_up_windows(int IC1_, int IC2_)
{
    float xwin, ywin;

    GWvport(0.0f, 0.90f, GWaspect(1), 1.0f);    // 数値の出力
    GWindow(0.0f, 3.0f, 80.0f, 0.0f);
    GWsavevp(IC1_);
    headings();
    GWvport(0.02f*GWaspect(1), 0.02f, GWaspect(1), 0.90f);
                                            // 粒子の軌跡
    compute_aspect_ratio((float)Lx, &xwin, &ywin);
    GWindow(0.0f, 0.0f, xwin, ywin);
    GWsavevp(IC2_);
    GWsetpen(BLACK, -1, -1, 4);
    GWrect(0.0f, 0.0f, (float)(Lx), (float)(Ly));
    GWanchor(1);
}

void headings()
{
    char Stmp1_[_LBUFF_];

    sprintf(Stmp1_, "time steps ");
    GWputtxt(1.0f, 1.0f, Stmp1_);
    sprintf(Stmp1_, "time ");
    GWputtxt(12.0f, 1.0f, Stmp1_);
    sprintf(Stmp1_, "energy ");
    GWputtxt(24.0f, 1.0f, Stmp1_);
    sprintf(Stmp1_, "<T> ");
    GWputtxt(36.0f, 1.0f, Stmp1_);
    sprintf(Stmp1_, "<P>");
    GWputtxt(48.0f, 1.0f, Stmp1_);
}

void Verlet(double *t, double *ke, double *pe, double *virial)
{
    int i;
    double xnew, ynew;

    for(i = 1; i <= N; ++i)
    {
        xnew = x[i] + vx[i]*dt + 0.5*ax[i]*dt2;
        ynew = y[i] + vy[i]*dt + 0.5*ay[i]*dt2;
        // 古い加速度を用いて速度の更新を部分的に行う
        vx[i] += 0.5*ax[i]*dt;
        vy[i] += 0.5*ay[i]*dt;
        x[i] = pbc(xnew, Lx);
        y[i] = pbc(ynew, Ly);
    }
    accel(pe, virial);                      // 新しい加速度
    (*ke) = 0;
    for(i = 1; i <= N; ++i)
    {
        // 新しい加速度を用いて速度の更新を行う
        vx[i] += 0.5*ax[i]*dt;
        vy[i] += 0.5*ay[i]*dt;
        (*ke) += vx[i]*vx[i] + vy[i]*vy[i];
    }
    (*ke) *= 0.5;
    (*t) += dt;
}

double pbc(double pos, double L)
{
    double pbc_;

    if(pos < 0)
        pbc_ = pos + L;
    else if(pos > L)
        pbc_ = pos - L;
    else
        pbc_ = pos;
    return pbc_;
}

void accel(double *pe, double *virial)
{
    int i, j;
    double dx, dy, fxij, fyij, pot;

    for(i = 1; i <= N; ++i)
    {
        ax[i] = 0;
        ay[i] = 0;
    }
    (*pe) = 0;
    (*virial) = 0;
    for(i = 1; i <= N - 1; ++i)             // 粒子 i に働く合力の計算
    {
        for(j = i + 1; j <= N; ++j)         // 粒子 j > i による寄与
        {
            dx = separation(x[i] - x[j], Lx);
            dy = separation(y[i] - y[j], Ly);
            // 換算単位系では 質量 = 1 なので 加速度 = 力 である
            force(dx, dy, &fxij, &fyij, &pot);
            ax[i] += fxij;
            ay[i] += fyij;
            ax[j] -= fxij;                  // ニュートンの第3法則
            ay[j] -= fyij;
            (*pe) += pot;
            (*virial) += dx*fxij + dy*fyij;
        }
    }
}

void force(double dx, double dy, double *fx, double *fy, double *pot)
{
    double f_over_r, r2, rm2, rm6;

    r2 = dx*dx + dy*dy;
    rm2 = 1/r2;
    rm6 = rm2*rm2*rm2;
    f_over_r = 24*rm6*(2*rm6 - 1)*rm2;
    (*fx) = f_over_r*dx;
    (*fy) = f_over_r*dy;
    (*pot) = 4*(rm6*rm6 - rm6);
}

double separation(double ds, double L)
{
    double separation_;

    if(ds > 0.5*L)
        separation_ = ds - L;
    else if(ds < -0.5*L)
        separation_ = ds + L;
    else
        separation_ = ds;
    return separation_;
}

void show_positions(char* flag, int IC2_)
{
    int i, k;

    GWselvp(IC2_);
    if(TBkeyinput())
    {
        k = TBgetkey();
        if(k == 'r')
        {
            GWerase(1, 0);
            GWrefresh();
            strcpy(flag, "");
        }
        else if(k == 's')
            strcpy(flag, "stop");
        else if(k == 'n')
            strcpy(flag, "no_show");
    }
    if(strcmp(flag, "no_show") != 0)
    {
        GWerase(1, 0);
        for(i = 1; i <= N; ++i)
            GWsetpxl((float)(x[i]), (float)(y[i]), RED);
    }
}

void show_output(double t, double ke, double pe, double virial, double *kecum, 
                 double *vcum, int *ncum, double area, int IC1_)
{
    double E, mean_ke, p;
    char Stmp1_[_LBUFF_];

    GWselvp(IC1_);
    GWsetogn(0, 2);
    GWerase(2, -1);
    ++(*ncum);
    sprintf(Stmp1_, "%6d ", (*ncum));
    GWputtxt(1.0f, 2.0f, Stmp1_);
    sprintf(Stmp1_, "%.2f ", t);            // 時間
    GWputtxt(12.0f, 2.0f, Stmp1_);
    E = ke + pe;                            // 全エネルギー
    sprintf(Stmp1_, "%.3f ", E);
    GWputtxt(24.0f, 2.0f, Stmp1_);
    (*kecum) += ke;
    (*vcum) += virial;
    mean_ke = (*kecum)/(*ncum);             // さらに N で割る必要がある
    p = mean_ke + (0.5*(*vcum))/(*ncum);    // 平均圧力 * 面積
    p /= area;
    sprintf(Stmp1_, "%.2f ", mean_ke/N);    // 平均の運動学的温度
    GWputtxt(36.0f, 2.0f, Stmp1_);
    sprintf(Stmp1_, "%.2f ", p);
    GWputtxt(48.0f, 2.0f, Stmp1_);
    GWsetogn(0, N + 1);
}

void save_config()
{
    int i;
    char file[_LBUFF_];
    FILE *FP1_;

    if(!GWinput("name of saved configuration = ", file, _LBUFF_)) return;
    FP1_ = fopen(file, "w");
    fprintf(FP1_, "%d\n", N);
    fprintf(FP1_, "%f\n", Lx);
    fprintf(FP1_, "%f\n", Ly);
    fprintf(FP1_, "        x(i)        y(i)\n");
    // 出力と出力の間のコンマは
    // True BASIC での読み込みを可能にするため
    for(i = 1; i <= N; ++i)
        fprintf(FP1_, "   %9.4f, %9.5f\n", x[i], y[i]);
    fprintf(FP1_, "       vx(i)       vy(i)\n");
    for(i = 1; i <= N; ++i)
        fprintf(FP1_, "   %9.4f, %9.5f\n", vx[i], vy[i]);
    fclose(FP1_);
}