// PROGRAM ideal
// demon algorithm for the one-dimensional, ideal classical gas 

#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;                              // initial demon energy
    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);
    // divide energy equally among particles
    vinitial = sqrt(2*(*E)/(*N));           // mass unity
    // all particles have same initial velocities
    for(i = 1; i <= (*N); ++i)
    {
        v[i] = vinitial;
    }
    // initialize sums
    (*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;          // trial change in velocity
        ipart = (int)(rnd(0)*N + 1);        // select random particle
        vtrial = v[ipart] + dv;             // trial velocity
        // trial energy change
        de = 0.5*(vtrial*vtrial - v[ipart]*v[ipart]);
        if(de <= (*Ed))
        {
            v[ipart] = vtrial;
            ++(*accept);
            (*Ed) -= de;
            (*E) += de;
        }
    }
    // accumulate data after each Monte Carlo step per particle
    (*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;                     // mean demon energy
    norm /= N;
    accept_prob = accept*norm;              // acceptance probability
    // system averages per particle
    Esave = Ecum*norm;                      // mean energy per system particle
    printf("mean demon energy = %f\n", Edave);
    printf("mean system energy per particle = %f\n", Esave);
    printf("acceptance probability = %f\n", accept_prob);
}