// PROGRAM qmwalk

#include "TrueBASIC.h"
int main();
void initial(double x[], int *N, int *N0, double *Vref, double *binsize, 
             double *ds, double *dt, int *mcs, int *nequil, double *Esum);
void walk(double x[], int *N, int N0, double *Vref, double *Vave, double ds, 
          double dt);
void data(double x[], double psi[], int N, double Vref, double Vave, 
          double *Esum, int imcs, double binsize);
void plot_psi(double psi[], double binsize);
double V(double x);

int main()
{
    int N, N0, i, imcs, mcs, nequil;
    double Esum, Vave, Vref, binsize, ds, dt, psi[1001], x[2001];

    for(i = 0; i <= 1000; ++i) {
       psi[i] = 0;
    }
    for(i = 0; i <= 2000; ++i) {
       x[i] = 0;
    }

    GWopen(0);

    initial(x, &N, &N0, &Vref, &binsize, &ds, &dt, &mcs, &nequil, &Esum);
    for(i = 1; i <= nequil; ++i)            // 平衡化
        walk(x, &N, N0, &Vref, &Vave, ds, dt);
    for(imcs = 1; imcs <= mcs; ++imcs)
    {
        walk(x, &N, N0, &Vref, &Vave, ds, dt);
        data(x, psi, N, Vref, Vave, &Esum, imcs, binsize);
    }
    plot_psi(psi, binsize);

    GWquit();
    return 0;
}

void initial(double x[], int *N, int *N0, double *Vref, double *binsize, 
             double *ds, double *dt, int *mcs, int *nequil, double *Esum)
{
    int i;
    double width;
    char Stmp1_[_LBUFF_];

    rnd(-1);
    printf("desired number of walkers = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", N0);
    printf("step length = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", ds);
    (*dt) = (*ds)*(*ds);                    // 時間刻み
    printf("number of Monte Carlo steps per walker = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", mcs);
    (*nequil) = (int)(0.4*(*mcs));
    // 初期の粒子数は任意に選ぶ
    (*N) = (*N0);
    // 粒子の初期位置はランダムに選ぶ
    width = 1;                              // 粒子の領域の初期幅
    (*Vref) = 0;
    for(i = 1; i <= (*N); ++i)
    {
        x[i] = (2*rnd(0) - 1)*width;
        (*Vref) += V(x[i]);
    }
    (*Vref) /= (*N);
    (*binsize) = 2*(*ds);                   // 配列 psi に対応する区分け用の箱の幅
    (*Esum) = 0;
    GWclear(-1);
    printf("      MC steps         N            E         Vref\n");
    printf("\n");
}

void walk(double x[], int *N, int N0, double *Vref, double *Vave, double ds, 
          double dt)
{
    int Nin, i;
    double Vsum, dV, potential;

    Nin = (*N);                             // 試行の初めにおける粒子数
    Vsum = 0;
    for(i = Nin; i >= 1; --i)
    {
        if(rnd(0) < 0.5)
            x[i] += ds;
        else
            x[i] -= ds;
        potential = V(x[i]);                // x でのポテンシャル
        dV = potential - (*Vref);
        if(dV < 0)                          // 粒子を加えるかどうかの検査
        {
            if(rnd(0) < -dV*dt)
            {
                ++(*N);
                x[(*N)] = x[i];             // 新しい粒子
                // 係数 2 は x に 2 個の粒子があるため
                Vsum += 2*potential;
            }
            else
                Vsum += potential;          // 古い粒子のみ
        }
        else
        {
            if(rnd(0) < dV*dt)              // 粒子を取り除くかどうかの検査
            {
                x[i] = x[(*N)];
                --(*N);
            }
            else
                Vsum += potential;
        }
    }
    (*Vave) = Vsum/(*N);                    // 平均ポテンシャル
    (*Vref) = (*Vave) - ((*N) - N0)/(N0*dt);         // 新しい参照エネルギー
}

void data(double x[], double psi[], int N, double Vref, double Vave, 
          double *Esum, int imcs, double binsize)
{
    int bin, i;

    // データを積算する
    (*Esum) += Vave;
    // 粒子を'箱'に入れる
    for(i = 1; i <= N; ++i)
    {
        bin = (int)(x[i]/binsize);
        ++(psi[bin + 500]);
    }
    if(((imcs) % (10)) == 0)
    {
        printf("%12d ", imcs);
        printf("%12d ", N);
        printf("%12.3f ", (*Esum)/imcs);
        printf("%12.3f\n", Vref);
    }
}

void plot_psi(double psi[], double binsize)
{
    int i, imax;
    double P, norm, pmax, sum = 0, ymax;

    i = 0;
    pmax = psi[i + 500]*psi[i + 500];
    sum = pmax;
    do
    {
        ++i;
        P = psi[i + 500]*psi[i + 500] + psi[-i + 500]*psi[-i + 500];
        if(P > pmax)
            pmax = P;
        sum += P;
    } while(!(P == 0));
    imax = i;
    norm = sqrt(sum*binsize);
    ymax = sqrt(pmax)/norm;
    GWindow((float)(-imax), (float)(-0.1*ymax), (float)(imax), 
            (float)(1.1*ymax));
    GWline((float)(-imax), 0.0f, (float)(imax), 0.0f);
    for(i = -imax; i <= imax; ++i)
        GWline((float)(i-1), (float)(psi[i + 499]/norm), (float)(i), 
               (float)(psi[i + 500]/norm));
}

double V(double x)
{
    double V_;

    V_ = 0.5*x*x;
    return V_;
}