#define ANISOTROPIC

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#define LATTICE_X   (30)
#define LATTICE_Y   (30)
#define LATTICE_Z   (30)

#define MIN(a,b)    (((a)>(b))?(b):(a))

struct configuration_t
{
    double phase[LATTICE_X][LATTICE_Y][LATTICE_Z];
};

void clear_configuration(struct configuration_t *cfgt)
{
    int i,j,k;

    for(i=0;i<LATTICE_X;i++)
        for(j=0;j<LATTICE_Y;j++)
            for(k=0;k<LATTICE_Z;k++)
                cfgt->phase[i][j][k]=0.0f;
}

double gen_random_number(void)
{
    return drand48();
}

int gen_random_int(void)
{
    return lrand48();
}

#define INIT_TO_ZERO    (0x1)

struct configuration_t *init_configuration(short flags)
{
    struct configuration_t *ret=malloc(sizeof(struct configuration_t));

    if(!ret)
        return NULL;

    if((flags&INIT_TO_ZERO)!=0)
        clear_configuration(ret);

    return ret;
}

void fini_configuration(struct configuration_t *cfgt)
{
    if(cfgt)
        free(cfgt);
}

double lsite_energy(struct configuration_t *cfgt,int x,int y,int z)
{
    double theta1,theta2,xy=0.0f;

    theta1=cfgt->phase[x][y][z];
    theta2=cfgt->phase[(x+1)%LATTICE_X][y][z];

#ifdef  ANISOTROPIC
    xy-=2.0f*cos(theta1-theta2);
#else
    xy-=cos(theta1-theta2);
#endif

    theta1=cfgt->phase[x][y][z];
    theta2=cfgt->phase[x][(y+1)%LATTICE_Y][z];

    xy-=cos(theta1-theta2);

    theta1=cfgt->phase[x][y][z];
    theta2=cfgt->phase[x][y][(z+1)%LATTICE_Z];

    xy-=cos(theta1-theta2);

    return xy;
}

double configuration_energy(struct configuration_t *cfgt)
{
    int i,j,k;
    double ret=0.0f;

    for(i=0;i<LATTICE_X;i++)
        for(j=0;j<LATTICE_Y;j++)
                for(k=0;k<LATTICE_Z;k++)
                    ret+=lsite_energy(cfgt,i,j,k);

    return ret;
}

void update_lsite(struct configuration_t *cfgt,int x,int y,int z)
{

#define AMPLITUDE   (2*M_PI)

    cfgt->phase[x][y][z]+=(gen_random_number()*AMPLITUDE);
}

void push_lsite(struct configuration_t *cfgt,int x,int y,int z,double *s)
{
    s[0]=cfgt->phase[x][y][z];
}

void pop_lsite(struct configuration_t *cfgt,int x,int y,int z,double *s)
{
    cfgt->phase[x][y][z]=s[0];
}

int update_configuration(struct configuration_t *cfgt,double beta)
{
    int i,j,k,x,y,z;
    double deltae,boltzmann_factor,p;
    double f[4];

    i=gen_random_int()%LATTICE_X;
    j=gen_random_int()%LATTICE_Y;
    k=gen_random_int()%LATTICE_Z;

    deltae=0.0f;

    for(x=i-1;x<=(i+1);x++)
        for(y=j-1;y<=(j+1);y++)
            for(z=k-1;z<=(k+1);z++)
                deltae-=lsite_energy(cfgt,(x+LATTICE_X)%LATTICE_X,(y+LATTICE_Y)%LATTICE_Y,(z+LATTICE_Z)%LATTICE_Z);

    push_lsite(cfgt,i,j,k,f);
    update_lsite(cfgt,i,j,k);

    for(x=i-1;x<=(i+1);x++)
        for(y=j-1;y<=(j+1);y++)
            for(z=k-1;z<=(k+1);z++)
                deltae+=lsite_energy(cfgt,(x+LATTICE_X)%LATTICE_X,(y+LATTICE_Y)%LATTICE_Y,(z+LATTICE_Z)%LATTICE_Z);

    boltzmann_factor=exp(-beta*deltae);
    p=MIN(1,boltzmann_factor);

    if(p<gen_random_number())
    {
        pop_lsite(cfgt,i,j,k,f);
        return 1;
    }

    return 2;
}

double magnetization(struct configuration_t *cfgt)
{
    int i,j,k;
    double a,b;

    a=b=0.0f;
    for(i=0;i<LATTICE_X;i++)
    {
        for(j=0;j<LATTICE_Y;j++)
        {
                for(k=0;k<LATTICE_Z;k++)
                {
                    a+=cos(cfgt->phase[i][j][k]);
                    b+=sin(cfgt->phase[i][j][k]);
                }
        }
    }

    return sqrt(pow(a,2.0f)+pow(b,2.0f));
}

double magnetization_per_spin(struct configuration_t *cfgt)
{
    double spins=LATTICE_X*LATTICE_Y*LATTICE_Z;

    return magnetization(cfgt)/spins;
}

double energy_per_spin(struct configuration_t *cfgt)
{
    double spins=LATTICE_X*LATTICE_Y*LATTICE_Z;

    return configuration_energy(cfgt)/spins;
}

double helicity_modulus(struct configuration_t *cfgt,double T)
{
    int i,j,k;
    double f,g;
    double spins=LATTICE_X*LATTICE_Y*LATTICE_Z;

    f=configuration_energy(cfgt);

    g=0.0f;
    for(i=0;i<LATTICE_X;i++)
        for(j=0;j<LATTICE_Y;j++)
            for(k=0;k<LATTICE_Z;k++)
                g+=sin(cfgt->phase[i][j][k]-cfgt->phase[(i+1)%LATTICE_X][j][k]);

    return (-0.5f*f-(1.0f/T)*pow(g,2.0f))/spins;
}

int main(void)
{
    int c,d;

    struct configuration_t *cfgt=init_configuration(INIT_TO_ZERO);
    double T,beta,mpsmean,epsmean,hmmean;

    int thermalization=100000;
    int iterations=150000;
    int valid;
    int runs;

    for(T=0.0f;T<=50.0f;T+=0.25f)
    {
        epsmean=mpsmean=hmmean=0.0f;
        beta=1.0/T;
        valid=0;
        runs=1;

        for(d=0;d<=runs;d++)
        {
            clear_configuration(cfgt);

            for(c=0;c<=iterations;c++)
            {
                double mps,eps,hm;
                int result=update_configuration(cfgt,beta);

                if((c>thermalization)&&(result==2))
                {
                    mps=eps=0.0f;
                    mps=magnetization_per_spin(cfgt);
                    eps=energy_per_spin(cfgt);
                    hm=helicity_modulus(cfgt,T);

                    mpsmean+=mps;
                    epsmean+=eps;
                    hmmean+=hm;
                    valid++;
                }
            }
        }

        epsmean/=((double)(valid));
        mpsmean/=((double)(valid));
        hmmean/=((double)(valid));

        printf("%f %d %f %f %f\n",T,valid,mpsmean,epsmean,hmmean);
        fflush(stdout);
    }

    fini_configuration(cfgt);

    return 0;
}