from visual import *
from visual.graph import *   #for the plots
from random import *
from math import *



L=10
spin=[]

###
#
MAXITS=2000 #needed to limit the simulation time spent at one temperature
WARM=50     #needed to ignore transient effects at the start of the simulation
#
###


### VPython displays
#
print "L=",L
scene=display(x=0,y=0,width=400,height=400,background=color.white,
              autoscale=0,range=(.65*L,0.65*L,0.65*L),fov=1e-8,
              title="simulation") #finer control of the visualization window

pic2=gdisplay(x=400,y=0,width=400,height=400,xmin=0.0,xmax=4.0,ymin=0.0,ymax=1.5,
              xtitle="Temp", ytitle="Magnetization",title="Magnetization (L=%d)" % L,
              foreground=color.black, background=color.white)
magnetization=gdots(gdisplay=pic2,color=color.blue)
#
###



for i in range(0,L):
    spin.append([])
    for j in range(0,L):
        spin[i].append(sphere(pos=vector(i-L/2,j-L/2,0),radius=0.45,color=color.red,state=1))


### Replace eternally simulating at T=3.0
### with a sequence a temperatures with MAXITS iterations each
#T=3.0
#while(1):
#
        
for T in [0.5,1.0,1.5,2.0,2.125,2.25,2.5,2.75,3.0,3.5,4.0]:
    print "temperature is",T
    magnet=0.0
    sus=0.0
    for its in range (0,MAXITS):
        scene.title="simulation: T=%5.3f (%d)" % (T,its)  #show our progress
        
        ### This block is indented because of the new outer loop over temperatures
        ###
        ### GENERATE MONTE-CARLO CONFIGURATION
        for i in range (0,L):
            for j in range (0,L):
                #if T>1.: rate(100) #slow down so we can see what is happening

                deltaE=(2.0/T)*spin[i][j].state*(
                                      spin[i][(j+1)%L].state+spin[i][(j-1)%L].state+
                                      spin[(i+1)%L][j].state+spin[(i-1)%L][j].state
                                                )
                if (exp(-deltaE)>random()):
                    spin[i][j].state=-spin[i][j].state
                    
                if (spin[i][j].state==1):
                    spin[i][j].color=color.red
                else:
                    spin[i][j].color=color.blue

        #
        ###



        if (its>=WARM): #ignore transients in the beginning
            dum=0.0
            for i in range (0,L):
                for j in range (0,L):
                    dum+=spin[i][j].state
            magnet+=abs(dum)


    magnet=magnet/((MAXITS-WARM))

    print "magnetization", magnet/(L*L)

    print
    magnetization.plot(pos=(T,abs(magnet)/(L*L)))