from visual import *
from visual.graph import *
from random import *
from math import *

##lattice length
L=4
# create a list 
spin=[]

for i in range(0,L):
# add a new list for each element of original list
# allows us to create a 2D lattice or array of spins 
    spin.append([])
    for j in range(0,L):
# i and j now give lattice coordinates of spin object
        spin[i].append(sphere(pos=vector(i-L/2,j-L/2,0),
                              radius=0.2,color=color.red,state=1))

        
## set ratio of energy to thermal energy kT
T=3.0
MAXITS=1000
WARM=10

## code for loop over T goes in here
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]:
    magnetization=0.0
    energy=0.0

## remember to indent all following code after adding the T loop.
    for its in range(0,MAXITS):

    # simulation code to update configuration of spins using Monte Carlo
    # algorithm

        for i in range (0,L):
            for j in range (0,L):
                #rate(200)
    #using periodic boundary conditions in 2D
                deltaE=(2.0/T)*spin[i][j].state*(
                                      spin[i][(j+1)%L].state+
                                      spin[i][(j-1+L)%L].state+
                                      spin[(i+1)%L][j].state+
                                      spin[(i+L-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
    #  end simulation part of code
        m=0
        e=0
        for i in range(0,L):
            for j in range(0,L):
                m+=spin[i][j].state
                e+=spin[i][j].state*(spin[(i+1)%L][j].state+
                spin[(i-1+L)%L][j].state+
                spin[i][(j+1)%L].state+
                spin[i][(j-1+L)%L].state)
        if(its>WARM):
            magnetization+=abs(m)
            energy+=e
            
    magnetization=magnetization/(MAXITS-WARM)
    energy=energy/(MAXITS-WARM)
    print magnetization, energy