/* Program that sums over all configurations of an Ising lattice of the 
   size  L0*L1 , at the moment the external field vanishes: h=0 ;
   The program computes the the energy density, the magnetisation
   per site defined as <|M|>/V and the probability distribution of the 
   magnetisation ; for 5^2 it takes less than a minute on a 1.8 Ghz PC
   it should not be used for larger lattices*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define D   2   /* should not be changed, program does not work correctly 
                   with other values */
#define L0  5   /* the lattice is of size L0*L1 */
#define L1  5
#define V  L0*L1

int spin[V];   /* contains the Ising spins */
double pmag[V+1];  /* the propability for a given magnetisation */
int nach[V][2*D]; /* encodes the neighbourhood relation on the lattice */
void main(void)
{
double boltz,z,xm,xe,beta;
int x,x0,x1,x0l,x1l,x0r,x1r,in;
int i,imax;
int i1,i2,i3;
int mag,ene;
printf("beta ? \n");                           scanf("%lf",&beta);
/* initialize the array for the propabilities */
for(i=0;i<=V;i++)
  {
  pmag[i]=0.;
  }

/*
the lattice points: example 4 x 4 lattice

   12  13  14  15

    8   9  10  11 

    4   5   6   7   

    0   1   2   3
*/
/* compute the neighbourhood relation */
for(x0=0;x0<L0;x0++)
  {
  x0l=x0-1; if(x0==0) x0l=L0-1;
  x0r=x0+1; if(x0r==L0) x0r=0;
  for(x1=0;x1<L1;x1++)
    {
    x1l=x1-1; if(x1==0) x1l=L1-1;
    x1r=x1+1; if(x1r==L1) x1r=0;
    x=x0*L1+x1;
    nach[x][0]=x0l*L1+x1;
    nach[x][1]=x0*L1+x1l;
    nach[x][2]=x0r*L1+x1;
    nach[x][3]=x0*L1+x1r;
    }
  } 

/* loop over the configurations */
z=0.;
xm=0.;
xe=0.;
imax=pow(2,V);
for(i=0;i<imax;i++)
  {
  i1=i;
  for(x=0;x<V;x++)
    {
    i2=i1/2;
    i3=i2*2;
    spin[x]=i1-i3; if(spin[x]==0) spin[x]=-1;
    i1=i2;
    }
/* compute the quantities of interest */
  mag=0;
  ene=0;
  for(x=0;x<V;x++)
    {
    mag+=spin[x];
    for(in=0;in<D;in++) ene-=spin[x]*spin[nach[x][in]];
    }
  boltz=exp(-beta*ene);
  z+=boltz;
  xm+=boltz*abs(mag);
  xe+=boltz*ene;
  pmag[(mag+V)/2]+=boltz;
  }
 printf("%f %f \n",xm/(z*L0*L1),-xe/(z*L0*L1));  
for(i=0;i<=V;i++)
  {
  pmag[i]=pmag[i]/z;
  }
for(i=0;i<=V;i++)
  {
  printf("%d  %f \n", 2*i-V, pmag[i]);
  }
}