// PROGRAM ideal
// 1次元の古典理想気体にたいするデーモン・アルゴリズム

#include "TrueBASIC.h"
int main();
void initial(int *N, double v[], int *mcs, double *E, double *Ed, double *Ecum, double *Edcum, int *accept, double *dvmax);
void change(int N, double v[], double *E, double *Ed, double *Ecum, double *Edcum, int *accept, double dvmax);
void averages(int N, double Ecum, double Edcum, int mcs, int accept);

int main()
{
    int N, accept, imcs, mcs;
    double E, Ecum, Ed, Edcum, dvmax, v[101];

    initial(&N, v, &mcs, &E, &Ed, &Ecum, &Edcum, &accept, &dvmax);
    for(imcs = 1; imcs <= mcs; ++imcs)
    {
        change(N, v, &E, &Ed, &Ecum, &Edcum, &accept, dvmax);
    }
    averages(N, Ecum, Edcum, mcs, accept);

    return 0;
}

void initial(int *N, double v[], int *mcs, double *E, double *Ed, double *Ecum, 
             double *Edcum, int *accept, double *dvmax)
{
    int i;
    double vinitial;
    char Stmp1_[_LBUFF_];

    rnd(-1);
    printf("number of particles = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", N);
    printf("initial energy of system = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", E);
    (*Ed) = 0;                              // 初期デーモン・エネルギー
    printf("number of MC steps per particle = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", mcs);
    printf("maximum change in velocity = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", dvmax);
    // エネルギーを各粒子に等しく分配する
    vinitial = sqrt(2*(*E)/(*N));           // 質量は 1
    // 各粒子の初速度は全て同じ
    for(i = 1; i <= (*N); ++i)
    {
        v[i] = vinitial;
    }
    // 和の初期化
    (*Ecum) = 0;
    (*Edcum) = 0;
    (*accept) = 0;
}

void change(int N, double v[], double *E, double *Ed, double *Ecum, 
            double *Edcum, int *accept, double dvmax)
{
    int i, ipart;
    double de, dv, vtrial;

    for(i = 1; i <= N; ++i)
    {
        dv = (2*rnd(0) - 1)*dvmax;          // 速度の試行変化
        ipart = (int)(rnd(0)*N + 1);        // ランダムに粒子を選択する
        vtrial = v[ipart] + dv;             // 試行速度
        // 試行のエネルギー変化
        de = 0.5*(vtrial*vtrial - v[ipart]*v[ipart]);
        if(de <= (*Ed))
        {
            v[ipart] = vtrial;
            ++(*accept);
            (*Ed) -= de;
            (*E) += de;
        }
    }
    // 粒子当りのモンテカルロ・ステップごとにデータを積算する
    (*Ecum) += (*E);
    (*Edcum) += (*Ed);
}

void averages(int N, double Ecum, double Edcum, int mcs, int accept)
{
    double Edave, Esave, accept_prob, norm;

    norm = 1.0/mcs;
    Edave = Edcum*norm;                     // 平均のデーモン・エネルギー
    norm /= N;
    accept_prob = accept*norm;              // 受容確率
    // １粒子当りの平均値
    Esave = Ecum*norm;                      // １粒子当りの平均エネルギー
    printf("mean demon energy = %f\n", Edave);
    printf("mean system energy per particle = %f\n", Esave);
    printf("acceptance probability = %f\n", accept_prob);
}