// PROGRAM boltzmann
// 1次元の粒子に対するメトロポリス・アルゴリズム

#include "TrueBASIC.h"
int main();
void initial(double *v0, double *E0, double *beta, int *mcs, int *nequil, 
             double *delta, int *nbin, double *delv);
void initialize_sums(double P[], double accum[], double *accept, int nbin);
void metropolis(double *v, double *E, double beta, double delta, 
                double *accept);
void data(double P[], double accum[], double v, double E, double delv);
void output(double P[], double accum[], int mcs, double *accept, int nbin, 
            double delv);

int main()
{
    int imcs, mcs, nequil, nbin;
    double E, P[201], accept, accum[4], beta, delta, delv, v;

    initial(&v, &E, &beta, &mcs, &nequil, &delta, &nbin, &delv);
    for(imcs = 1; imcs <= nequil; ++imcs)    // 系を平衡化する
        metropolis(&v, &E, beta, delta, &accept);
    initialize_sums(P, accum, &accept, nbin);
    for(imcs = 1; imcs <= mcs; ++imcs)
    {
        metropolis(&v, &E, beta, delta, &accept);
        // 各試行変化の後にデータを積算する
        data(P, accum, v, E, delv);
    }
    output(P, accum, mcs, &accept, nbin, delv);

    return 0;
}

void initial(double *v0, double *E0, double *beta, int *mcs, int *nequil, double *delta, int *nbin, double *delv)
{
    double T, vmax;
    char Stmp1_[_LBUFF_];

    rnd(-1);
    printf("temperature = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", &T);
    (*beta) = 1/T;
    printf("number of Monte Carlo steps = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", mcs);
    (*nequil) = (int)(0.1*(*mcs));
    printf("initial speed = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", v0);
    (*E0) = 0.5*(*v0)*(*v0);                // 初期運動エネルギー
    printf("maximum change in velocity = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", delta);
    vmax = 10*sqrt(T);
    (*nbin) = 20;                           // 速度の区分数
    (*delv) = vmax/(*nbin);                 // 速度の間隔
}

void initialize_sums(double P[], double accum[], double *accept, int nbin)
{
    int i, ibin;

    for(ibin = -nbin; ibin <= nbin; ++ibin)
        P[ibin + 100] = 0;
    for(i = 1; i <= 3; ++i)
        accum[i] = 0;
    (*accept) = 0;
}

void metropolis(double *v, double *E, double beta, double delta, double *accept)
{
    double dE, dv, vtrial;

    dv = (2*rnd(0) - 1)*delta;              // 試行速度変化
    vtrial = (*v) + dv;                     // 試行速度
    dE = 0.5*(vtrial*vtrial - (*v)*(*v));   // 試行エネルギー変化
    if(dE > 0)
        if(exp(-beta*dE) < rnd(0))
            return;                         // 試行を受け入れない
    (*v) = vtrial;
    ++(*accept);
    (*E) += dE;
}

void data(double P[], double accum[], double v, double E, double delv)
{
    int ibin;

    accum[1] += E;
    accum[2] += E*E;
    accum[3] += v;
    ibin = (int)(v/delv);
    ++(P[ibin + 100]);
}

void output(double P[], double accum[], int mcs, double *accept, int nbin, 
            double delv)
{
    int ibin;
    double E2ave, Eave, prob, sigma2, v, vave;

    (*accept) /= mcs;
    printf("acceptance probability = %f\n", (*accept));
    vave = accum[3]/mcs;
    printf("mean velocity = %f\n", vave);
    Eave = accum[1]/mcs;
    printf("mean energy = %f\n", Eave);
    E2ave = accum[2]/mcs;
    sigma2 = E2ave - Eave*Eave;
    printf("sigma_E =  %f\n", sqrt(sigma2));
    printf("\n");
    printf("           v         P(v)\n");
    printf("\n");
    v = -nbin*delv;
    for(ibin = -nbin; ibin <= nbin; ++ibin)
    {
        if(P[ibin+100] > 0)
        {
            prob = P[ibin+100]/mcs;
            printf("%12f       %6.3f\n", v, prob);
        }
        v += delv;
    }
}