// PROGRAM lattice_gas
// 格子気体中の粒子の拡散のシミュレーション

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

int main();
void initial(int *L, double *L_2, int *N, double *tmax, int *im2, int *IC2_);
void lattice(double site[][41], int x[], int y[], double x0[], double y0[], 
             int L, int N, int *im1, int *IC1_);
void move(double site[][41], int x[], int y[], int L, int N, int im1, int IC1_);
void direction(int *xtemp, int *ytemp, int L);
void data(int x[], int y[], double x0[], double y0[], int L, double L_2, int N, 
          double t, int im2, int IC2_);
void separation(double *dx, double *dy, int L, double L_2);

int main()
{
    int im1, im2;
    int IC1_ = 1, IC2_ = 2, L, N, x[1201], y[1201];
    double L_2, site[41][41], t, tmax, x0[1201], y0[1201];

    GWopen(0);

    initial(&L, &L_2, &N, &tmax, &im2, &IC2_);
    lattice(site, x, y, x0, y0, L, N, &im1, &IC1_);
    for(t = 1; t <= tmax; ++t)
    {
        move(site, x, y, L, N, im1, IC1_);
        data(x, y, x0, y0, L, L_2, N, t, im2, IC2_);
    }

    GWquit();
    return 0;
}

void initial(int *L, double *L_2, int *N, double *tmax, int *im2, int *IC2_)
{
    char Stmp1_[_LBUFF_];
    rnd(-1);
    (*L) = 40;                              // 格子の1辺の長さ
    (*L_2) = 0.5*(*L);
    (*N) = 500;                             // 粒子数
    (*tmax) = 50;                   // 1粒子あたりのモンテカルロ・ステップ数
    GWcolor(BLACK, CLR_BACKGRND);
    GWcolor(WHITE, CLR_TEXT);
    GWvport((float)(0.52*GWaspect(1)), 0.0f, GWaspect(1), 1.0f);
    GWindow(-1.0f, -1.0f, (float)((*tmax) + 1), (float)((*L_2) + 1));
    GWsetogn(0, (*N)+1);
    GWsetpen(WHITE, -1, -1, MIX_COPYPEN);
    *im2 = GWcmbmrk(0, 1, (float)((*L_2)/(*tmax)/4), WHITE, -1, MIX_COPYPEN);
    sprintf(Stmp1_, "Plot of <R^2(t)> versus t");
    GWputtxt(1.0f, (float)((*L_2)), Stmp1_);
    GWrect(0.0f, 0.0f, (float)((*tmax)), (float)((*L_2)));
    GWsavevp(*IC2_);
}

void lattice(double site[][41], int x[], int y[], double x0[], double y0[], 
             int L, int N, int *im1, int *IC1_)
{
    int i, ix, iy, xadd, yadd;
    float r, xwin, ywin;
    char Stmp1_[_LBUFF_];

    GWvport(0.0f, 0.2f, 0.5f*GWaspect(1), 1.0f);
    r = 0.5f;                                // 箱のサイズ
    compute_aspect_ratio(L+r, &xwin, &ywin);
    GWindow(-0.1f*xwin, -0.1f*ywin, xwin, ywin);
    GWsetogn(0, N+1);
    GWcolor(WHITE, CLR_TEXT);
    sprintf(Stmp1_, "density =  %f", (double)N/L/L);
    GWputtxt(0.0f, ywin*0.93f, Stmp1_);
    GWsetpen(RED, -1, -1, MIX_COPYPEN);
    GWrect(1-r, 1-r, L+r, L+r);
    for(iy = 1; iy <= L; ++iy)              // 格子点の描画
    {
        for(ix = 1; ix <= L; ++ix)
        {
            site[ix][iy] = 0;
            GWsetpxl((float)(ix), (float)(iy), RED);
        }
    }
    i = 0;
    GWanchor(1);

    GWsetpen(GWC_BLUE, -1, -1, 4);
    *im1 = GWcmbmrk(0, 6, 1.0f, GWC_BLUE, -1, MIX_COPYPEN);

    do
    {
        xadd = (int)(L*rnd(0)) + 1;
        yadd = (int)(L*rnd(0)) + 1;
        if(site[xadd][yadd] == 0)
        {
            ++i;                            // 粒子数への加算
            site[xadd][yadd] = -1;          // 占有された格子点
            GWsetogn(0, i);
            GWputcmb(*im1, (float)(xadd), (float)(yadd), 0.0f, 0.0f, 0);
            x[i] = xadd;
            y[i] = yadd;
            x0[i] = x[i];                   // t = 0 における x 座標
            y0[i] = y[i];
        }
    } while(i != N);
    GWsavevp(*IC1_);
}

void move(double site[][41], int x[], int y[], int L, int N, int im1, int IC1_)
{
    int i, particle, xtrial, ytrial;

    GWselvp(IC1_);
    GWsetogn(0, N+1);
    for(particle = 1; particle <= N; ++particle)
    {
        i = (int)(N*rnd(0)) + 1;
        xtrial = x[i];
        ytrial = y[i];
        direction(&xtrial, &ytrial, L);
        if(site[xtrial][ytrial] >= 0)       // 空の格子点
        {
            site[x[i]][y[i]] = 0;
            GWsrect((float)(x[i]-0.5), (float)(y[i]-0.5), 
                    (float)(x[i]+0.5), (float)(y[i]+0.5), -1);
            GWerase(i, 0);
            GWsetogn(0, i);
            x[i] = xtrial;
            y[i] = ytrial;
            site[x[i]][y[i]] = -1;          // 新しく占有された格子点
            GWputcmb(im1, (float)(x[i]), (float)(y[i]), 0.0f, 0.0f, 0);
        }
    }
}

void direction(int *xtemp, int *ytemp, int L)
// 方向をランダムに選び，周期境界条件を使う
{
    int dir;

    dir = (int)(4*rnd(0)) + 1;
    switch(dir)
    {
        case 1:
            if(++(*xtemp) > L) (*xtemp) = 1;
            break;
        case 2:
            if(--(*xtemp) < 1) (*xtemp) = L;
            break;
        case 3:
            if(++(*ytemp) > L) (*ytemp) = 1;
            break;
        case 4:
            if(--(*ytemp) < 1) (*ytemp) = L;
    }
}

void data(int x[], int y[], double x0[], double y0[], int L, double L_2, int N, 
          double t, int im2, int IC2_)
{
    int i;
    double R2, R2bar, dx, dy;

    R2 = 0;
    for(i = 1; i <= N; ++i)
    {
        dx = x[i] - x0[i];
        dy = y[i] - y0[i];
        separation(&dx, &dy, L, L_2);
        R2 += dx*dx + dy*dy;
    }
    R2bar = R2/N;
    GWselvp(IC2_);
    GWsetogn(0, N+1);
    GWputcmb(im2, (float)t, (float)R2bar, 0.0f, 0.0f, 0);
}

void separation(double *dx, double *dy, int L, double L_2)
// 距離の計算に周期境界条件を使う
{
    if(fabs((*dx)) > L_2)
        (*dx) -= sgn((*dx))*L;
    if(fabs((*dy)) > L_2)
        (*dy) -= sgn((*dy))*L;
}