from IsingLattice import *
from matplotlib import pylab as pl
import numpy as np

n_rows = 8
n_cols = 8
il = IsingLattice(n_rows, n_cols)
il.lattice = np.ones((n_rows, n_cols))
spins = n_rows*n_cols
runtime = 200000
times = range(runtime)
temps = np.arange(0.25, 5.0, 0.1)
#default was 1.5,3.0,0.1
energies = []
magnetisations = []
energysq = []
magnetisationsq = []
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 = il.statistics()
    energies.append(aveE)
    energysq.append(aveE2)
    magnetisations.append(aveM)
    magnetisationsq.append(aveM2)
    print(n_cycles)
 
    #Variance of energies is equal to the difference between avgE2 and (avgE)^2
    var_energies=np.subtract(energysq,np.multiply(energies,energies))
    
    #Standard Deviation of the energies is equal to the square root of the Variance
    SD_energies=np.sqrt(var_energies) 
    
    #Standard Error of the Mean (SEM) is equal to the Standard Deviation divided by the square root of the Number of Cycles
    SEM_energies=np.divide(SD_energies,np.sqrt(n_cycles))
    
    #Taking into account that we want the SEM of energy/spin
    SEM_energies_per_spin=SEM_energies/spins
    
    #Applying similar logic to calculate the SEM for the values of magnetisation:
    var_mag=np.subtract(magnetisationsq,np.multiply(magnetisations,magnetisations))
    SD_mag=np.sqrt(var_mag) 
    SEM_mag=np.divide(SD_mag,np.sqrt(n_cycles))
    SEM_mag_per_spin=SEM_mag/spins
    
    #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.errorbar(temps, np.array(energies)/spins,yerr=SEM_energies_per_spin) 
magax.errorbar(temps, np.array(magnetisations)/spins,yerr=SEM_mag_per_spin)
pl.show()

final_data = np.column_stack((temps, energies, energysq, magnetisations, magnetisationsq))
np.savetxt("8x8_mine.dat",final_data)