// PROGRAM ising
// メトロポリス・アルゴリズムを用いた
// 正方格子上のイジング・モデルのモンテカルロ・シミュレーション

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

int main();
void initial(int *N, int *L, double *T, int spin[][33], int *mcs, int *nequil, 
             double w[], double *E, double *M);
void plotspins(int L, int spin[][33]);
void initialize_correl(double Ce[], double Cm[], double esave[], 
                       double msave[], int *nsave);
void initialize(double accum[], double *accept);
void Metropolis(int N, int L, int spin[][33], double *E, double *M, double w[], 
                double *accept);
int DeltaE(int x, int y, int L, int spin[][33]);
void show_change(int i, int j, int L, int spin[][33]);
void data(double E, double M, double accum[]);
void output(double T, int N, int mcs, double accum[], double accept);
void save_config(int L, double T, int spin[][33]);
void read_config(int L, double *T, int spin[][33]);
void correl(double Ce[], double Cm[], double E, double M, double esave[], 
            double msave[], int pass, int nsave);
void c_output(double Ce[], double Cm[], double accum[], int mcs, int nsave);

int main()
{
    int L, N, i, mcs, nsave, nequil, pass, spin[33][33];
    double Ce[21], Cm[21], E, M, T, accept, accum[11], esave[101], msave[101], 
           w[17];

    GWopen(0);

    initial(&N, &L, &T, spin, &mcs, &nequil, w, &E, &M);
    plotspins(L, spin);
    for(i = 1; i <= nequil; ++i)         // 系を平衡化する
        Metropolis(N, L, spin, &E, &M, w, &accept);
    initialize(accum, &accept);
    initialize_correl(Ce, Cm, esave, msave, &nsave);
    for(pass = 1; pass <= mcs; ++pass)      // スピンを更新してデータを積算する
    {
        Metropolis(N, L, spin, &E, &M, w, &accept);
        data(E, M, accum);
        correl(Ce, Cm, E, M, esave, msave, pass, nsave);
    }
    output(T, N, mcs, accum, accept);
    c_output(Ce, Cm, accum, mcs, nsave);

    GWquit();
    return 0;
}

void initial(int *N, int *L, double *T, int spin[][33], int *mcs, int *nequil, 
             double w[], double *E, double *M)
{
    int dE, right, up, x, y;
    double sum;
    char Stmp1_[_LBUFF_];

    rnd(-1);
    printf("linear dimension of lattice = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", L);
    (*N) = (*L)*(*L);                       // スピン数
    // 温度は J/k を単位に測られる
    printf("temperature = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", T);
    printf("# MC steps per spin for equilibrium = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", nequil);
    printf("# MC steps per spin for data = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", mcs);
    (*M) = 0;
    for(y = 1; y <= (*L); ++y)              // ランダムな初期配置
    {
        for(x = 1; x <= (*L); ++x)
        {
            if(rnd(0) < 0.5)
                spin[x][y] = 1;             // 上向きスピン
            else
                spin[x][y] = -1;
            (*M) += spin[x][y];             // 全磁化
        }
    }
    (*E) = 0;
    for(y = 1; y <= (*L); ++y)              // 初期エネルギーの計算
    {
        if(y == (*L))
            up = 1;                         // 周期的境界条件
        else
            up = y + 1;
        for(x = 1; x <= (*L); ++x)
        {
            if(x == (*L))
                right = 1;
            else
                right = x + 1;
            sum = spin[x][up] + spin[right][y];
            (*E) -= spin[x][y]*sum;         // 全エネルギー
        }
    }
    // ボルツマン確率の比の計算
    for(dE = -8; dE <= 8; dE += 4)
        w[dE + 8] = exp(-dE/(*T));
}

void plotspins(int L, int spin[][33])
{
    int i, j;
    float xwin, ywin;

    GWclear(-1);
    compute_aspect_ratio((float)(L+1), &xwin, &ywin);
    GWindow(0.0f, 0.0f, xwin, ywin);
    for(i = 1; i <= L; ++i)
    {
        for(j = 1; j <= L; ++j)
        {
            GWerase((i-1)*L+j, 0);
            GWsetogn(0, (i-1)*L+j);
            if(spin[i][j] == 1)
                GWsrect((float)i, (float)j, i+1.0f, j+1.0f, RED);
            else
                GWsrect((float)i, (float)j, i+1.0f, j+1.0f, GREEN);
        }
    }
}

void initialize_correl(double Ce[], double Cm[], double esave[], 
                       double msave[], int *nsave)
{
    int i;

    (*nsave) = 10;
    for(i = 1; i <= (*nsave); ++i)
    {
        esave[i] = 0;
        msave[i] = 0;
    }
    for(i = 1; i <= (*nsave); ++i)
    {
        Ce[i] = 0;
        Cm[i] = 0;
    }
}

void initialize(double accum[], double *accept)
// 積算された磁化とエネルギーの値を保存するために配列を使う
// 配列を使うことで容易に他の量を加えることができる
{
    int i;

    for(i = 1; i <= 5; ++i)
        accum[i] = 0;
    (*accept) = 0;
}

void Metropolis(int N, int L, int spin[][33], double *E, double *M, double w[], 
                double *accept)
// 単位スピン当りモンテカルロ・ステップの実行
{
    int dE, ispin, x, y;

    for(ispin = 1; ispin <= N; ++ispin)
    {
        x = (int)(L*rnd(0) + 1);            // ランダムな x 座標
        y = (int)(L*rnd(0) + 1);            // ランダムな y 座標
        dE = DeltaE(x, y, L, spin);         // エネルギー変化の計算
        if(rnd(0) <= w[dE + 8])
        {
            spin[x][y] *= - 1;       // スピンを反転する
            ++(*accept);
            (*M) += 2*spin[x][y];
            (*E) += dE;
            show_change(x,y,L,spin);
        }
    }
}

int DeltaE(int x, int y, int L, int spin[][33])
// 周期的境界条件
{
    int DeltaE_, down, left, right, up;

    if(x == 1)
        left = spin[L][y];
    else
        left = spin[x - 1][y];
    if(x == L)
        right = spin[1][y];
    else
        right = spin[x + 1][y];
    if(y == 1)
        down = spin[x][L];
    else
        down = spin[x][y - 1];
    if(y == L)
        up = spin[x][1];
    else
        up = spin[x][y + 1];
    DeltaE_ = 2*spin[x][y]*(left + right + up + down);
    return DeltaE_;
}

void show_change(int i, int j, int L, int spin[][33])
{
    static int show = 1;

    if(GWmouse(NULL, NULL, NULL)) {
       if(show)
          show = 0;
       else {
          show = 1;
          plotspins(L, spin);
       }
       while(GWmouse(NULL, NULL, NULL));
    }
    if(!show) return;
    GWerase((i-1)*L+j, 0);
    GWsetogn(0, (i-1)*L+j);
    if(spin[i][j] == 1)
        GWsrect((float)i, (float)j, i+1.0f, j+1.0f, RED);
     else
        GWsrect((float)i, (float)j, i+1.0f, j+1.0f, GREEN);
}

void data(double E, double M, double accum[])
// 各スピン当りモンテカルロ・ステップの後にデータを積算する
{
    accum[1] += E;
    accum[2] += (E*E);
    accum[3] += M;
    accum[4] += (M*M);
    accum[5] += fabs(M);
}

void output(double T, int N, int mcs, double accum[], double accept)
{
    double abs_mave, e2ave, eave, m2ave, mave, norm;

    norm = 1.0/(mcs*N);                       // スピン当りの平均
    accept *= norm;
    eave = accum[1]*norm;
    e2ave = accum[2]*norm;
    mave = accum[3]*norm;
    m2ave = accum[4]*norm;
    abs_mave = accum[5]*norm;
    printf("temperature =  %f\n", T);
    printf("acceptance probability =  %f\n", accept);
    printf("mean energy per spin =  %f\n", eave);
    printf("mean squared energy per spin =  %f\n", e2ave);
    printf("mean magnetization per spin =  %f\n", mave);
    printf("mean of absolute magnetization per spin =  %f\n", abs_mave);
    printf("mean squared magnetization per spin =  %f\n", m2ave);
}

void save_config(int L, double T, int spin[][33])
{
    int x, y;
    char file[_LBUFF_], Stmp1_[_LBUFF_];
    FILE *FP2_;

    printf("name of file for last configuration = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%[^\n]", file);
    FP2_ = fopen(file, "w");
    fprintf(FP2_, "%f\n", T);
    for(y = 1; y <= L; ++y)
        for(x = 1; x <= L; ++x)
            fprintf(FP2_, "%d\n", spin[x][y]);
    fclose(FP2_);
}

void read_config(int L, double *T, int spin[][33])
{
    int x, y;
    char file[_LBUFF_], Stmp1_[_LBUFF_];
    FILE *FP1_;

    printf("filename?");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%[^\n]", file);
    FP1_ = fopen(file, "r");
    fgets(Stmp1_, _LBUFF_, FP1_);
    sscanf(Stmp1_, "%lg", T);
    for(y = 1; y <= L; ++y)
    {
        for(x = 1; x <= L; ++x)
        {
            fgets(Stmp1_, _LBUFF_, FP1_);
            sscanf(Stmp1_, "%d", &spin[x][y]);
        }
    }
    fclose(FP1_);
}

void correl(double Ce[], double Cm[], double E, double M, double esave[], 
            double msave[], int pass, int nsave)
// 時間相関関数用にデータを積算する
// 最後の nsave 個の M と E の値を格納する
{
    int index, index0, tdiff;

    // index0 = 最も古い値の配列要素番号
    index0 = ((pass-1) % (nsave)) + 1;
    if(pass > nsave)
    {
        // nsave 個の値が格納された後の Ce と Cm を計算する
        index = index0;
        for(tdiff = nsave; tdiff >= 1; --tdiff)
        {
            Ce[tdiff] += E*esave[index];
            Cm[tdiff] += M*msave[index];
            if(++index > nsave)
                index = 1;
        }
    }
    // 最新の値を最も古い値の配列位置に格納する
    esave[index0] = E;
    msave[index0] = M;
}

void c_output(double Ce[], double Cm[], double accum[], int mcs, int nsave)
// 時間相関関数の計算
{
    int tdiff;
    double e2bar, ebar, m2bar, mbar, norm;

    ebar = accum[1]/mcs;
    e2bar = accum[2]/mcs;
    Ce[0] = e2bar - ebar*ebar;
    mbar = accum[3]/mcs;
    m2bar = accum[4]/mcs;
    Cm[0] = m2bar - mbar*mbar;
    norm = 1.0/(mcs - nsave);
    printf("\n");
    printf("           t        Ce(t)        Cm(t)\n");
    printf("\n");
    for(tdiff = 1; tdiff <= nsave; ++tdiff)
    {
        // 相関関数は C(t=0) = 1 及び C(infinity) = 0 
        // となるように定義される
        Ce[tdiff] = (Ce[tdiff]*norm - ebar*ebar)/Ce[0];
        Cm[tdiff] = (Cm[tdiff]*norm - mbar*mbar)/Cm[0];
        printf("%12d %12f %12f\n", tdiff, Ce[tdiff], Cm[tdiff]);
    }
}