// PROGRAM demon
// 零磁場下の d = 1 イジング・モデルにおけるデーモン・アルゴリズム

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

int main();
void initial(int *N, int spin[], double *E, double *Ed, double *M, int *mcs, 
             double *Ecum, double *Edcum, double *Mcum, double *M2cum);
void setupscreen(int N, int spin[]);
void change(int N, int spin[], double *E, double *Ed, double *M, int *accept);
void data(double E, double Ed, double M, double *Ecum, double *Edcum, 
          double *Mcum, double *M2cum);
void averages(int N, double Ecum, double Edcum, double Mcum, double M2cum, 
              int mcs, int accept);
void showspin(int N, int dir, int isite);

int main()
{
    int N, accept = 0, imcs, mcs, spin[1001];
    double E, Ecum, Ed, Edcum, M, M2cum, Mcum;

    GWopen(0);

    initial(&N, spin, &E, &Ed, &M, &mcs, &Ecum, &Edcum, &Mcum, &M2cum);
    setupscreen(N, spin);
    for(imcs = 1; imcs <= mcs; ++imcs)
    {
        change(N, spin, &E, &Ed, &M, &accept);
        data(E, Ed, M, &Ecum, &Edcum, &Mcum, &M2cum);
    }
    averages(N, Ecum, Edcum, Mcum, M2cum, mcs, accept);

    GWquit();
    return 0;
}

void initial(int *N, int spin[], double *E, double *Ed, double *M, int *mcs, double *Ecum, double *Edcum, double *Mcum, double *M2cum)
{
    int isite;
    double Etot;
    char Stmp1_[_LBUFF_];

    rnd(-1);
    printf("number of spins = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", N);
    // 全エネルギーを 4J の整数倍になるように選ぶ
    // 結合定数 J は 1 にとる
    printf("desired total energy = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", &Etot);
    Etot = 4*(int)(Etot/4.0);
    printf("number of Monte Carlo steps per spin = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", mcs);
    // 最低エネルギー状態でのスピンの初期配置
    for(isite = 1; isite <= (*N); ++isite)
        spin[isite] = 1;
    (*M) = (*N);                            // 正味の磁化
    // 系の初期エネルギーの計算
    (*E) = -(*N);                           // 周期的境界条件
    (*Ed) = (Etot - (*E));
    printf("total energy =  %f\n", (*E) + (*Ed));
    // 和の変数の初期化
    (*Ecum) = 0;
    (*Edcum) = 0;
    (*Mcum) = 0;
    (*M2cum) = 0;
}

void setupscreen(int N, int spin[])
{
    int i;
    double dtheta;
    float r, xwin, ywin;

    r = N/(2*pi);
    compute_aspect_ratio(r+2, &xwin, &ywin);
    GWindow(-xwin, -ywin, xwin, ywin);
    dtheta = 2*pi/N;

    for(i = 1; i <= N; ++i)
        showspin(N, spin[i], i);
}

void change(int N, int spin[], double *E, double *Ed, double *M, int *accept)
{
    int i, isite;
    double de, left, right;

    for(i = 1; i <= N; ++i)
    {
        isite = (int)(rnd(0)*N + 1);        // ランダムにスピンを選ぶ
        // 隣接するスピンの値を決定する
        if(isite == 1)
            left = spin[N];
        else
            left = spin[isite - 1];
        if(isite == N)
            right = spin[1];
        else
            right = spin[isite + 1];
        // 試行エネルギー変化
        de =  2*spin[isite]*(left + right);
        if(de <= (*Ed))
        {
            // スピン反転動力学
            spin[isite] = -spin[isite];
            (*M) = (*M) + 2*spin[isite];
            ++(*accept);                    // 受け入れられた試行変化の回数
            (*Ed) -= de;
            (*E) += de;
        }
        showspin(N, spin[isite], isite);
    }
}

void data(double E, double Ed, double M, double *Ecum, double *Edcum, 
          double *Mcum, double *M2cum)
{
    // データを積算する
    (*Ecum) += E;
    (*Edcum) += Ed;
    (*Mcum) += M;
    (*M2cum) += M*M;
}

void averages(int N, double Ecum, double Edcum, double Mcum, double M2cum, 
              int mcs, int accept)
{
    double Eave, Edave, M2ave, Mave, T, accept_prob, norm;

    norm = 1.0/mcs;                       // データは各試行の後に集められている
    Edave = Edcum*norm;
    printf("mean demon energy = %f\n", Edave);
    T = 4/log(1 + 4/Edave);
    printf("T = %f\n", T);
    Eave = Ecum*norm;
    printf("<E> =  %f\n", Eave);
    Mave = Mcum*norm;
    printf("<M> = %f\n", Mave);
    M2ave = M2cum*norm;
    printf("<M*M> = %f\n", M2ave);
    accept_prob = accept*norm/N;
    printf("acceptance probability = %f\n", accept_prob);
}

void showspin(int N, int dir, int isite)
{
    float r, theta, x, y;

    r = N/(2*pi);
    theta = isite/r;
    x = r*cos(theta);
    y = r*sin(theta);
    GWerase(isite, 0);
    GWsetogn(0, isite);
    if(dir == 1)
        GWsrect(x-0.25f, y-0.25f, x+0.25f, y+0.25f, RED);
    else
        GWsrect(x-0.25f, y-0.25f, x+0.25f, y+0.25f, BLUE);
}