// PROGRAM conduct
// many demon algorithm for Ising chain

#include "TrueBASIC.h"
int main();
void initial(int *N, int spin[], int *nmcs, int *nvisit);
void change(int N, int spin[], double Ed[], double *accept);
void heat(int N, int spin[], double Edsum[]);
void data(int N, int spin[], double Ed[], double Edsum[], double Msum[]);
void output(int N, double Edsum[], double Msum[], int nmcs, double accept);

int main()
{
    int N, imcs, nmcs, nvisit, spin[1001], Itmp1_;
    double Ed[1001], Edsum[1001], Msum[1001], accept = 0;
    for(Itmp1_ = 0; Itmp1_ <= 1000; ++Itmp1_) {
       Ed[Itmp1_] = 0;
       Edsum[Itmp1_] = 0;
       Msum[Itmp1_] = 0;
    }

    // heat added at spin 1 and subtracted at spin N
    initial(&N, spin, &nmcs, &nvisit);
    for(imcs = 1; imcs <= nmcs; ++imcs)
    {
        change(N, spin, Ed, &accept);
        if((imcs & nvisit) == 0)
            heat(N, spin, Edsum);
        data(N, spin, Ed, Edsum, Msum);
    }
    output(N, Edsum, Msum, nmcs, accept);

    return 0;
}

void initial(int *N, int spin[], int *nmcs, int *nvisit)
{
    int i;
    char Stmp1_[_LBUFF_];

    rnd(-1);
    printf("number of spins = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", N);
    printf("number of MC steps per spin = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", nmcs);
    printf("MC steps between updates of end spins = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", nvisit);
    // initial random configuration
    for(i = 1; i <= (*N); ++i)
    {
        if(rnd(0) > 0.5)
            spin[i] = 1;
        else
            spin[i] = -1;
    }
}

void change(int N, int spin[], double Ed[], double *accept)
// spin flip dynamics
{
    int isite, i;
    double de;

    // do one Monte Carlo step
    for(i = 2; i <= N - 1; ++i)
    {
        // pick spin at random from spins 2 to N - 1
        isite = (int)(rnd(0)*(N - 2) + 2);
        // trial energy change
        de = 2*spin[isite]*(spin[isite - 1] + spin[isite + 1]);
        if(de <= Ed[isite])
        {
            spin[isite] *= - 1;
            ++(*accept);
            Ed[isite] -= de;
        }
    }
}

void heat(int N, int spin[], double Edsum[])
{
    // attempt to add energy at spin 1 and remove it at spin N
    // possible only if spins 1 and 2 are aligned
    // and spins N and N - 1 are not aligned
    if((spin[1]*spin[2] == 1) && (spin[N]*spin[N - 1] == -1))
    {
        Edsum[1] +=  2;
        Edsum[N] -= 2;
        spin[1] *= -1;
        spin[N] *= -1;
    }
}

void data(int N, int spin[], double Ed[], double Edsum[], double Msum[])
{
    int i;

    for(i = 2; i <= N - 1; ++i)
    {
        Edsum[i] += Ed[i];
        Msum[i] += spin[i];
    }
}

void output(int N, double Edsum[], double Msum[], int nmcs, double accept)
{
    int i;
    double Edave, M, accept_prob, norm, temperature;

    norm = 1.0/nmcs;
    accept_prob = accept*norm/(N - 2);
    printf("acceptance probability =  %f\n", accept_prob);
    printf("\n");
    printf("           i        Ed(i)            T         M(i)\n");
    printf("\n");
    for(i = 2; i <= N-1; ++i)
    {
        Edave = Edsum[i]*norm;
        temperature = 0;
        if(Edave != 0)
            if((1 + 4/Edave) > 0)
                temperature = 4/log(1 + 4/Edave);
        M = Msum[i]*norm;
        printf("%12d %12f %12f %12f\n", i, Edave, temperature, M);
    }
}