// PROGRAM hd
// 固い円盤の系の動力学
// 部分的に Allen と Tildesley の Fortran のプログラムを基礎にしている

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

int main();
void initial(double *vsum, double *rho, double *area);
void check_momentum();
void minimum_collision_time(int *i, int *j, double *tij);
void move(double tij);
void reset_list(int i, int j);
double separation(double ds, double L);
double pbc(double pos, double L);
void check_overlap();
void kinetic_energy(double *ke, int IC1_);
void uplist(int i);
void downlist(int j);
void check_collision(int i, int j);
void contact(int i, int j, double *virial);
void show_output(double t, int collisions, double temperature, double vsum, 
                 double rho, double area);
void save_config();
void set_up_windows(int IC1_, int IC2_);
void headings();
void show_positions(char* flag);
void show_disks(char* flag);
void getflag(char *flag);

double Lx;
double Ly;
int N;
double collision_time[101];
int partner[101];
double t;
double timebig;
double vx[101];
double vy[101];
double x[101];
double y[101];

int ishow = 1;

int main()
{
    int IC1_ = 1, IC2_ = 2, collisions, i, j;
    double area, ke, rho, temperature, tij, virial, vsum;
    char flag[_LBUFF_];

    GWopen(0);

    initial(&vsum, &rho, &area);
    set_up_windows(IC1_, IC2_);
    kinetic_energy(&ke, IC1_);
    GWanchor(1);
    GWselvp(IC2_);
    if(ishow)
        show_positions(flag);
    else
        show_disks(flag);
    temperature = ke/N;
    strcpy(flag, "show_output");
    collisions = 0;                         // 衝突回数
    do
    {
        minimum_collision_time(&i, &j, &tij);
        // 粒子を前方に動かし，衝突時間を tij だけ小さくする
        move(tij);
        t += tij;
        ++collisions;
        getflag(flag);
        if(ishow)
            show_positions(flag);
        else
            show_disks(flag);
        contact(i, j, &virial);             // 衝突の動力学を計算する
        vsum += virial;
        // 衝突相手のリストをリセットする
        reset_list(i, j);
        if(strlen(flag) && strcmp(flag, "no_show")) {
            GWselvp(IC1_);
            show_output(t, collisions, temperature, vsum, rho, area);
            GWselvp(IC2_);
            if((strcmp(flag, "show_output") == 0) || 
               (strcmp(flag, "change") == 0)) strcpy(flag, "");
        }
    } while(!(strcmp(flag, "stop") == 0));
    save_config();

    GWquit();
    return 0;
}

void initial(double *vsum, double *rho, double *area)
{
    int i;
    double ax, ay, col, nx, ny, row, vmax;
    char file[_LBUFF_], heading[_LBUFF_], start[_LBUFF_], Stmp1_[_LBUFF_];
    FILE *FP1_;

    t = 0;
    printf("read file (f) or lattice start (l) = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%[^\n]", start);
    if(strcmp(start, "f") == 0 || strcmp(start, "F") == 0)
    {
        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 if(strcmp(start, "l") == 0 || strcmp(start, "L") == 0)
    {
        rnd(-1);
        printf("N = ");
        fgets(Stmp1_, _LBUFF_, stdin);
        sscanf(Stmp1_, "%d", &N);       // sqr(N) が整数になるように N を選ぶ
        printf("Lx = ");
        fgets(Stmp1_, _LBUFF_, stdin);
        sscanf(Stmp1_, "%lg", &Lx);
        printf("Ly = ");
        fgets(Stmp1_, _LBUFF_, stdin);
        sscanf(Stmp1_, "%lg", &Ly);
        printf("vmax = ");
        fgets(Stmp1_, _LBUFF_, stdin);
        sscanf(Stmp1_, "%lg", &vmax);
        nx = sqrt(N);
        ny = nx;
        if(nx >= Lx || nx >= Ly)
        {
            printf("box too small\n");
            GWquit();
            exit(0);
        }
        ax = Lx/nx;                         // "格子" 間隔
        ay = Ly/ny;
        i = 0;
        for(col = 1; col <= nx; ++col)
        {
            for(row = 1; row <= ny; ++row)
            {
                ++i;
                x[i] = (col - 0.5)*ax;
                y[i] = (row - 0.5)*ay;
                // 位置と速度をランダムに選ぶ
                vx[i] = (2*rnd(0) - 1)*vmax;
                vy[i] = (2*rnd(0) - 1)*vmax;
            }
        }
    }
    check_overlap();                        // 2つの円盤が重なっているか調べる
    check_momentum();
    (*area) = Lx*Ly;
    (*rho) = N/(*area);
    timebig = 1.0e10;
    (*vsum) = 0;                            // ビリアル和
    for(i = 1; i <= N; ++i)
        partner[i] = N;
    collision_time[N] = timebig;
    // 初期状態の衝突のリストを構成する
    for(i = 1; i <= N; ++i)
        uplist(i);
    GWsetogn(0, N + 1);
}

void check_momentum()
{
    int i;
    double vxsum, vysum;

    vxsum = 0;
    vysum = 0;
    for(i = 1; i <= N; ++i)
    {
        vxsum += vx[i];      // 重心の全速度 (全運動量)
        vysum += vy[i];
    }
    vxsum = vxsum/N;
    vysum = vysum/N;
    for(i = 1; i <= N; ++i)
    {
        vx[i] -= vxsum;
        vy[i] -= vysum;
    }
}

void minimum_collision_time(int *i, int *j, double *tij)
// 最短衝突時間を見つける
{
    int k;

    (*tij) = timebig;
    for(k = 1; k <= N; ++k)
    {
        if(collision_time[k] < (*tij))
        {
            (*tij) = collision_time[k];
            (*i) = k;
        }
    }
    (*j) = partner[(*i)];
}

void move(double tij)
{
    int k;

    for(k = 1; k <= N; ++k)
    {
        collision_time[k] -= tij;
        x[k] += vx[k]*tij;
        y[k] += vy[k]*tij;
        x[k] = pbc(x[k], Lx);
        y[k] = pbc(y[k], Ly);
    }
}

void reset_list(int i, int j)
// 関係する粒子についての衝突リストをリセットする
{
    int k, test;

    for(k = 1; k <= N; ++k)
    {
        test = partner[k];
        if(k == i || test == i || k == j || test == j)
            uplist(k);
    }
    downlist(i);
    downlist(j);
}

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_;
}

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

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

void check_overlap()
{
    int i, j;
    double dx, dy, r, r2, tol;

    tol = 1.0e-4;
    for(i = 1; i <= N - 1; ++i)
    {
        for(j = i + 1; j <= N; ++j)
        {
            dx = separation(x[i] - x[j], Lx);
            dy = separation(y[i] - y[j], Ly);
            r2 = dx*dx + dy*dy;
            if(r2 < 1)
            {
                r = sqrt(r2);
                if((1 - r) > tol)
                {
                    printf("particles  %d and  %doverlap\n", i, j);
                    GWquit();
                    exit(0);
                }
            }
        }
    }
}

void kinetic_energy(double *ke, int IC1_)
{
    int i;
    char Stmp1_[_LBUFF_];

    GWselvp(IC1_);
    (*ke) = 0;
    for(i = 1; i <= N; ++i)
        (*ke) += vx[i]*vx[i] + vy[i]*vy[i];
    (*ke) *= 0.5;
    sprintf(Stmp1_, "%.3f", (*ke)/N);
    GWputtxt(36.0f, 2.0f, Stmp1_);
}

void uplist(int i)
// j > i の粒子との衝突を探す
{
    int j;

    if(i == N) return;
    collision_time[i] = timebig;
    for(j = i + 1; j <= N; ++j)
        check_collision(i, j);
}

void downlist(int j)
// 粒子 i < j の衝突を調べる
{
    int i;

    if(j == 1) return;
    for(i = 1; i <= j - 1; ++i)
        check_collision(i, j);
}

void check_collision(int i, int j)
{
    double bij, discr, dvx, dvy, dx, dy, r2, tij, v2, xcell, ycell;

    // i と，j の周期的な像との衝突を考慮する
    for(xcell = -1; xcell <= 1; ++xcell)
    {
        for(ycell = -1; ycell <= 1; ++ycell)
        {
            dx = x[i] - x[j] + xcell*Lx;
            dy = y[i] - y[j] + ycell*Ly;
            dvx = vx[i] - vx[j];
            dvy = vy[i] - vy[j];
            bij = dx*dvx + dy*dvy;
            if(bij < 0)
            {
                r2 = dx*dx + dy*dy;
                v2 = dvx*dvx + dvy*dvy;
                discr = bij*bij - v2*(r2 - 1);
                if(discr > 0)
                {
                    tij = (-bij - sqrt(discr))/v2;
                    if(tij < collision_time[i])
                    {
                        collision_time[i] = tij;
                        partner[i] = j;
                    }
                }
            }
        }
    }
}

void contact(int i, int j, double *virial)
// 粒子 i と j の接触時の衝突の動力学の計算
{
    double delvx, delvy, dvx, dvy, dx, dy, factor;

    dx = separation(x[i] - x[j], Lx);
    dy = separation(y[i] - y[j], Ly);
    dvx = vx[i] - vx[j];
    dvy = vy[i] - vy[j];
    factor = dx*dvx + dy*dvy;
    delvx = - factor*dx;
    delvy = - factor*dy;
    vx[i] += delvx;
    vx[j] -= delvx;
    vy[i] += delvy;
    vy[j] -= delvy;
    (*virial) = delvx*dx + delvy*dy;
}

void show_output(double t, int collisions, double temperature, double vsum, 
                 double rho, double area)
{
    double mean_pressure, mean_virial;
    char Stmp1_[_LBUFF_];

    GWsetogn(0, N + 2);
    GWerase(N + 2, -1);
    mean_virial = vsum/(2*t);
    mean_pressure = rho*temperature + mean_virial/area;
    sprintf(Stmp1_, "%6d", collisions);
    GWputtxt(1.0f, 2.0f, Stmp1_);
    sprintf(Stmp1_, "%.2f", t);
    GWputtxt(12.0f, 2.0f, Stmp1_);
    sprintf(Stmp1_, "%.2f", mean_pressure);
    GWputtxt(24.0f, 2.0f, Stmp1_);
    GWsetogn(0, N + 1);
}

void save_config()
{
    int i, j;
    double tij;
    char file[_LBUFF_], Stmp1_[_LBUFF_];
    FILE *FP1_;

    // 最終配置で粒子が接触していないように粒子を動かす
    minimum_collision_time(&i, &j, &tij);
    move(0.5*tij);
    if(!GWinput("name of saved configuration = ", Stmp1_, _LBUFF_)) return;
    sscanf(Stmp1_, "%[^\n]", file);
    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");
    for(i = 1; i <= N; ++i)
        fprintf(FP1_, "  %10.5f, %10.5f\n", x[i], y[i]);
    fprintf(FP1_, "       vx(i)       vy(i)\n");
    for(i = 1; i <= N; ++i)
        fprintf(FP1_, "  %10.5f, %10.5f\n", vx[i], vy[i]);
    fclose(FP1_);
}

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((float)(0.02*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));
    GWsetmrk(6, 1.0f, RED, -1, 4);
}

void headings()
{
    char Stmp1_[_LBUFF_];

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

void show_positions(char* flag)
{
    int i;

    if(strcmp(flag, "no_show") != 0)
    {
        GWerase(N + 1, 0);
        for(i = 1; i <= N; ++i)
            GWsetpxl((float)(x[i]), (float)(y[i]), RED);
    }
}

void show_disks(char* flag)
{
    static int ncolor = 0;
    int i;

    if(strcmp(flag, "no_show") != 0)
    {
        ncolor = (ncolor % 6) + 1;
        if(ncolor == 1) GWrefresh();
        ncolor = GWkrgb((int)(rnd(0)*256), (int)(rnd(0)*256), (int)(rnd(0)*256));
        for(i = 1; i <= N; ++i) {
            GWsetogn(0, -i);
            GWsetpen(ncolor, -1, -1, 4);
            GWputmrk((float)(x[i]), (float)(y[i]));
            GWflush(-i);
        }
    }
}

void getflag(char *flag)
{
    int k;

    if(TBkeyinput())
    {
        k = TBgetkey();
        strcpy(flag, "show_output");
        if(k == 'r' || k == 'R') {
            GWerase(0, -1);
            GWrefresh();
        } else if(k == 's' || k == 'S')
            strcpy(flag, "stop");
        else if(k == 'n' || k == 'V')
            strcpy(flag, "no_show");
        else if(k == 'c' || k == 'C' || k == '\t') {
            GWerase(0, -1);
            if(ishow)
                ishow = 0;
            else
                ishow = 1;
        }
    }
}