#pragma rtGlobals=1		// Use modern global access method.

Menu "Ising"
	"Metropolis Initialize"
	"MetroploisChange"
End

Macro MetroploisChange(GridSize,GridName,NoChange,Temp,makeplot,updatestep)
	variable/g gGridsize,gNoChange,gTemp,gmakeplot,gupdatestep
	string/g gGridName
	variable gridsize=gGridSize,NoChange=gNoChange,Temp=gTemp,makeplot=gmakeplot,updatestep=gupdatestep
	String GridName=gGridName
	Prompt GridSize,"Enter the size of the simulation grid: (100)"
	Prompt GridName,"Enter the wave name for the simulation:"
	Prompt NoChange,"Enter the number of times to randomly Change:"
	Prompt Temp,"Enter the temperature: (0 to 10 say)"
	Prompt makeplot,"Make a new image plot?",popup,"Yes;No"
	
	silent 1
	gGridsize=GridSize;gGridName=GridName;gNoChange=NoChange;gTemp=Temp;gmakeplot=makeplot
	gupdatestep=updatestep
	Make/n=(Gridsize,Gridsize)/o $GridName
	Initialize($GridName,GridSize)
	variable energy=getEnergy($gridname,gridsize)
	print "Energy = "+num2str(energy)+"; Energy per spin = "+num2str(energy/gridsize^2)
	
	if(Makeplot==1)
		//display the ising thing
		Display /W=(5.25,41.75,357.75,371)
		AppendImage $gridname
		ModifyImage $gridname ctab= {*,*,PlanetEarth256,0}
		ModifyGraph mirror=2
	endif
	
	ChangeIsing($GridName,GridSize,NoChange,Temp,updatestep)
EndMacro

Macro MetropolisInitialize(GridSize,GridName)
	variable/g gGridsize
	string/g gGridName
	variable gridsize=gGridSize
	String GridName=gGridName
	Prompt GridSize,"Enter the size of the simulation grid: (100)"
	Prompt GridName,"Enter the wave name for the simulation:"
	
	silent 1
	gGridsize=GridSize;gGridName=GridName
	Make/n=(Gridsize,Gridsize)/o $GridName
	Initialize($GridName,GridSize)
	variable energy=getEnergy($gridname,gridsize)
	print "Energy = "+num2str(energy)+"; Energy per spin = "+num2str(energy/gridsize^2)
	
	//display the ising thing
	Display /W=(5.25,41.75,357.75,371)
	AppendImage $gridname
	ModifyImage $gridname ctab= {*,*,PlanetEarth256,0}
	ModifyGraph mirror=2
EndMacro

Function ChangeIsing(GridName,GridSize,NoChange,Temp,updatestep)
	wave gridname
	variable gridsize,NoChange,Temp,updatestep
	
	variable counter=0,Energy0=0,Energy,DeltaEnergy,value=1,xval,yval,test,intg=1
	//First Get the energy of the matrix
	Energy0=getEnergy(Gridname,gridsize)
	do
		//Then change a random position's spin
		xval=round(abs(enoise(gridsize-1)))
		yval=round(abs(enoise(gridsize-1)))
		gridname[xval][yval]*=-1
		//Then recalculate the energy and energy difference
		Energy=getEnergy(Gridname,gridsize)
		deltaEnergy=Energy-energy0
		//Then if delta energy is negative or 0 keep and repeat
//print energy,energy0,deltaenergy,xval,yval
		if(deltaEnergy>0)//so if you increase the energy
			value=abs(enoise(1)) //random number between 0 and 1
//print value,exp(-energy/temp)
			test=exp(-energy/temp)
			if(numtype(test)>0)//NAN
				counter-=1
				Energy=Energy0
				gridname[xval][yval]*=-1
				//beep
			else
				if(value>test)//then don't allow should be E/kT not just T
					counter-=1
					Energy=Energy0
					gridname[xval][yval]*=-1
					//beep
					beep
				endif
			endif
			
		endif
		Energy0=Energy
		counter+=1
		if(counter/updatestep==intg)
			doUpdate
			print counter
			intg+=1
		endif
	while(counter<noChange)
	//Return(-energy)
End

Function getEnergy(Gridname,gridsize)
	wave gridname
	variable gridsize
	
	variable counterx=0,countery=0,Energy=0
	//assign +-1 to each position from 0 to gridsize-1
	do
		do
			//gridname(counterx,countery)=value/abs(value)
			if(counterx<(gridsize-1)&&countery<(gridsize-1))
				//then this is a normal position and add d (2) associations
				Energy+=GridName[counterx][countery]*Gridname[counterx+1][countery]
				Energy+=GridName[counterx][countery]*Gridname[counterx][countery+1]
			else
				if(counterx==(gridsize-1)&&countery<(gridsize-1))
					//add just one or d-1
					Energy+=GridName[counterx][countery]*Gridname[counterx][countery+1]
				else//then countery==gridsize-1
					if(counterx<(gridsize-1)&&countery==(Gridsize-1))
						Energy+=GridName[counterx][countery]*Gridname[counterx+1][countery]
					endif
				endif
			endif
			countery+=1
		while(countery<gridsize)
		countery=0
		counterx+=1
	while(counterx<gridsize)
	Return(-energy)
End

Function Initialize(Gridname,gridsize)
	wave gridname
	variable gridsize
	
	variable counterx=0,countery=0,value
	//assign +-1 to each position from 0 to gridsize-1
	do
		do
			value=enoise(100)
			gridname[counterx][countery]=value/abs(value)
			countery+=1
		while(countery<gridsize)
		countery=0
		counterx+=1
	while(counterx<gridsize)
	
End