from IsingLattice import *
from matplotlib import pylab as pl
import numpy as np
n_rows = 32
n_cols = 32
il = IsingLattice(n_rows, n_cols)
il.lattice = np.ones((n_rows, n_cols))
spins = n_rows*n_cols
runtime = 1000
times = range(runtime)
temps = np.arange(0.25, 5.25, 0.25)
energies = []
magnetisations = []
energysq = []
magnetisationsq = []
E_errors = []
M_errors = []
for t in temps:
for i in times:
if i % 100 == 0:
print(t, i)
energy, magnetisation = il.montecarlostep(t)
aveE, aveE2, aveM, aveM2, n_cycles, std_errorE, std_errorM = il.statistics()
energies.append(aveE)
energysq.append(aveE2)
magnetisations.append(aveM)
magnetisationsq.append(aveM2)
E_errors.append(std_errorE)
M_errors.append(std_errorM)
#reset the IL object for the next cycle
il.E = 0.0
il.E2 = 0.0
il.M = 0.0
il.M2 = 0.0
il.n_cycles = 0
fig = pl.figure()
enerax = fig.add_subplot(2,1,1)
enerax.set_ylabel("Energy per spin")
enerax.set_xlabel("Temperature")
enerax.set_ylim([-2.1, 2.1])
magax = fig.add_subplot(2,1,2)
magax.set_ylabel("Magnetisation per spin")
magax.set_xlabel("Temperature")
magax.set_ylim([-1.1, 1.1])
enerax.plot(temps, np.array(energies)/spins)
enerax.errorbar(temps, np.array(energies)/spins, xerr = 0, yerr = np.array(E_errors)/spins)
magax.plot(temps, np.array(magnetisations)/spins)
magax.errorbar(temps, np.array(magnetisations)/spins, xerr = 0, yerr = np.array(M_errors)/spins)
pl.show()
final_data = np.column_stack((temps, energies, energysq, magnetisations, magnetisationsq))
np.savetxt("32x32.dat", final_data)
