//This program is an implementation of the Metropolis Algorithm for a 2D Ising lattice
//50x50 with J=1, k=1, and periodic boundary condition
//physical grid is [0,Nmax-1]x[0,Nmax-1]
//spins are either -1 or 1

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

#define Nmax 5 //max size of lattice

/*constants for ran1(long *)*/
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

int lattice[Nmax][Nmax];//lattice[0..Nmax-1][0..Nmax-1]
int endi=Nmax-1;
double eng[5];//lookup table for energy changes
long seed1,seed2,seed3;//seed for ran1()

void metropolis(double ttemp,double *m, double *c);
int moveforward(int x_i,int x_j,double *delE);
void getboltzfactor(double bbeta);
int totalmag(void);
double hamiltonian(void);
void initialize_lattice(int choice);
double ran1(long *);

//
//stuffs for dianogistic
// 
void printflattice(int);
void printfm(double, double);

FILE *digfile2;

void printlattice(int index)
{
	int i,j;
	char filename[25];

	sprintf(filename,"lattice%d.dat",index);
	digfile2=fopen(filename,"w");
	for(i=0;i<Nmax;i++)
	{
		for(j=0;j<Nmax;j++)
			fprintf(digfile2,"%d\t",lattice[i][j]);
		fprintf(digfile2,"\n");
	}
	fclose(digfile2);
}

void printfm(double time,double m)
{
	fprintf(digfile2,"%lg\t%lg\n",time,m);
}
//
// end of stuff for dianogistic
//

main()
{
	double temp,tstep;
	double mag;//avg magnetization
	double spec_heat; //specific heat
	int i,Ntemp;

	FILE *outfile;

	//priming output data file
	outfile=fopen("ising.dat","w");
	fprintf(outfile,"temp\tmag\tspecheat\n");

	//variables for temperature increments
	Ntemp=25;
	tstep=5.0/Ntemp;

	//calculation
	initialize_lattice(1);

//printlattice(0);
//hamiltonian();

	temp=5.0;

	for(i=1;i<=Ntemp;i++)
	{
		metropolis(temp,&mag,&spec_heat);
		fprintf(outfile,"%lg\t%lg\t%lg\n",temp,mag,spec_heat);
		fflush(outfile);
		printf("Temperature done: %lg\n",temp);
		temp-=tstep;
	}

	return 0;
}

void metropolis(double ttemp,double *m, double *c)
{
	long N,n;
	double beta,dE;
	int M;
	double Msum;
	double E,Esum,E2sum;
	int x,y;

	//reseeding random number generators for each run
	//this is not absolutely necessary
	seed1=-(int)(ran1(&seed1)*777);

	//
	//take out transient; wait for lattice to reach equilibrium at ttemp
	//
	N=(long int)10000*(Nmax*Nmax);

	beta=1.0/ttemp;
	getboltzfactor(beta);//this will generate the boltzman factor used for the Metro step

	for(n=1;n<=N;n++)
	{
		//pick a random spin to flip in the 2D grid
		x=(int)(Nmax*ran1(&seed1));
		y=(int)(Nmax*ran1(&seed1));
		moveforward(x,y,&dE);

/*
if(n==1000*Nmax*Nmax)
printlattice(1000);
if(n==2000*Nmax*Nmax)
printlattice(2000);
if(n==3000*Nmax*Nmax)
printlattice(3000);
if(n==5000*Nmax*Nmax)
printlattice(5000);
*/

	}

	//
	//take data
	//
	//initial values of M and E fro the lattice
	M=totalmag();
	E=hamiltonian();
	Msum=0.0;
	Esum=E2sum=0.0;

	N=(long int)450000*(Nmax*Nmax);

	//N is also used for the time to take temporal average
	for(n=1;n<=N;n++)
	{
		//pick a random spin to flip in the 2D grid
		x=(int)(Nmax*ran1(&seed1));
		y=(int)(Nmax*ran1(&seed1));
		if(moveforward(x,y,&dE))
		{
			//if step taken
			//update E and M
			E+=dE;
			if(lattice[x][y]>0)
				M+=2;//if new spin is up
			else
				M-=2;//if new spin is down
			
		}

		//collecting statistics
		//avg magnetization
		Msum+=fabs((double)M);
		//for specific heat, we need to calculate <E> and <E^2>
		Esum+=E;
		E2sum+=E*E;

	}

	Esum/=N;
	E2sum/=N;
	*m=Msum/N/(Nmax*Nmax);//average magnetization per spin
	*c=beta*beta*(E2sum-Esum*Esum)/(Nmax*Nmax);//specific heat of the lattice at temperature ttemp

}

int moveforward(int x_i,int x_j,double *delE)
//if step taken, spin will be flipped here
//other output is delE and the count for the Metro step
{
	int nn_i,nn_j;
	int spinsum;

	//
	//calculate E_x-E_x'
	//
	//finding energy change from interaction with neighbors of (x_i,x_j)
	//nn of (i,j) are (i+-1,j+-1) with periodic boundary
	//left
	if(!x_i)
		nn_i=Nmax-1;
	else
		nn_i=x_i-1;
	spinsum=lattice[nn_i][x_j];
	//right
	if(x_i==Nmax-1)
		nn_i=0;
	else
		nn_i=x_i+1;
	spinsum+=lattice[nn_i][x_j];
	//top
	if(!x_j)
		nn_j=Nmax-1;
	else
		nn_j=x_j-1;
	spinsum+=lattice[x_i][nn_j];
	//bottom
	if(x_j==Nmax-1)
		nn_j=0;
	else
		nn_j=x_j+1;
	spinsum+=lattice[x_i][nn_j];


	//Now, we take the Metropolis step
	//delE=2J*s_x*SUM(nn)
	*delE=2*lattice[x_i][x_j]*spinsum;

	if(*delE>0)
	{
		//delE>0 take the Metropolis step only with prob e^(-beta*delE)
		if(ran1(&seed1)<eng[abs(spinsum)])
		{
			//if random <= prob take the Metropolis step
			lattice[x_i][x_j]=-lattice[x_i][x_j];//flip the spin
			return 1;//count this step
		}
	}
	else
	{
		//delE<=0 take the Metropolis step
		lattice[x_i][x_j]=-lattice[x_i][x_j];//flip the spin
		return 1;//count this step
	}

	return 0;//not taking step
}

int totalmag(void)
{
	int i,j;
	int spinsum;
	
	spinsum=0;
	for(i=0;i<Nmax;i++)
		for(j=0;j<Nmax;j++)
		{
			spinsum+=lattice[i][j];
		}
	return spinsum;
}

double hamiltonian(void)
{
	int i,j;
	double energy;

	energy=0.0;
	for(i=0;i<endi;i++)
	{	
		for(j=0;j<endi;j++)
		{
			energy-=(lattice[i][j]*(lattice[i+1][j]+lattice[i][j+1]));			
		}
		energy-=(lattice[i][endi]*(lattice[i+1][endi]+lattice[i][0]));//wrap around pts
	}
	for(j=0;j<endi;j++)
		energy-=(lattice[endi][j]*(lattice[0][j]+lattice[endi][j+1]));//wrap around pts
	energy-=(lattice[endi][endi]*(lattice[0][endi]+lattice[endi][0]));//wrap around pts
	return energy;
}

void initialize_lattice(int choice)
{
	int i,j;

	seed1=-1;//initialize seed for random generator ran1()

	//to initialize lattice either at T=0 or T=infy
	if(choice)
	{
		//lattice at T=infy will have random uncorrelated spins
		for(i=0;i<Nmax;i++)
			for(j=0;j<Nmax;j++)
				{
				if(ran1(&seed1)>0.5)
					lattice[i][j]=1;
				else
					lattice[i][j]=-1;
				}
	}
	else
	{
		//lattice at T=0 will have all spins in the same direction
		for(i=0;i<Nmax;i++)
			for(j=0;j<Nmax;j++)
				lattice[i][j]=1;
	}
}

void getboltzfactor(double bbeta)
{

	//we will also initialize a look up table for the energy changes due to spin flips
	//E_x-E_x'=2Js_xSum(nn) [see note in text for this reduction]
	//Sum(nn) can only take on z+1 values: -z,-z+2,-z+4,,,,+z
	//for our regular 2D lattice, z (num of nn neighbors) = 4
	//only positive values needed to be evaluated

	eng[2]=exp(-bbeta*4.0);//e^(-beta*delE) -> delE=2J*2
	eng[4]=exp(-bbeta*8.0);//e^(-beta*delE) -> delE=2J*4

}

//random number generators

double ran1(long *idum)
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	double temp;
	
	if(*idum <= 0 || !iy)
		{
		if(-(*idum) < 1)
			*idum = 1;
		else
			*idum = -(*idum);
		for(j=NTAB+7;j>=0;j--)
			{
			k=(*idum)/IQ;
			*idum=IA*(*idum-k*IQ)-IR*k;
			if(*idum<0)
				*idum += IM;
			if(j<NTAB)
				iv[j]= *idum;
			}
		iy=iv[0];
		}
	k=(*idum)/IQ;
	*idum=IA*(*idum-k*IQ)-IR*k;
	if(*idum < 0)
		*idum += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j]= *idum;
	if((temp=AM*iy) >RNMX)
		return RNMX;
	else
		return temp;
}