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

#First, we use the NumPy loadtxt function to read our previously saved file for the 2x2 lattice
data_2x2= np.loadtxt("2x2_mine.dat")


#Define number of spins in lattice (so that we can use this to convert values of energy to energy/spin)
spins=4

#Variance of energy/spin is equal to the difference between avg(E^2) and (avgE)^2 where E is energy/spin
var_energies_2x2=np.subtract(data_2x2[:,2],np.multiply(data_2x2[:,1],data_2x2[:,1]))

#Heat Capacity is equal to the Variance / T^2  

Heat_Capacity_2x2 = np.divide(var_energies_2x2, np.multiply(data_2x2[:,0], data_2x2[:,0]))

#Plotting heat capacity/spin vs temperature
pl.plot(data_2x2[:,0],(Heat_Capacity_2x2),'bo-', markersize = 3, label = 'Lattice Size = 2x2')

#We can use a similar method for the rest of the lattice sizes

#4x4 lattice:

data_4x4= np.loadtxt("4x4_mine.dat")
spins=16
var_energies_4x4=np.subtract(data_4x4[:,2],np.multiply(data_4x4[:,1],data_4x4[:,1]))
Heat_Capacity_4x4 = np.divide(var_energies_4x4, np.multiply(data_4x4[:,0], data_4x4[:,0]))
pl.plot(data_4x4[:,0],(Heat_Capacity_4x4),'go-', markersize = 3, label = 'Lattice Size = 4x4')


#8x8 lattice:

data_8x8= np.loadtxt("8x8_mine.dat")
spins=64
var_energies_8x8=np.subtract(data_8x8[:,2],np.multiply(data_8x8[:,1],data_8x8[:,1]))
Heat_Capacity_8x8 = np.divide(var_energies_8x8, np.multiply(data_8x8[:,0], data_8x8[:,0]))
pl.plot(data_8x8[:,0],(Heat_Capacity_8x8),'ro-', markersize = 3, label = 'Lattice Size = 8x8')

#16x16 lattice:

data_16x16= np.loadtxt("16x16_mine.dat")
spins=256
var_energies_16x16=np.subtract(data_16x16[:,2],np.multiply(data_16x16[:,1],data_16x16[:,1]))
Heat_Capacity_16x16 = np.divide(var_energies_16x16, np.multiply(data_16x16[:,0], data_16x16[:,0]))
pl.plot(data_16x16[:,0],(Heat_Capacity_16x16),'co-', markersize = 3, label = 'Lattice Size = 16x16')

#32x32 lattice:

data_32x32= np.loadtxt("32x32_mine.dat")
spins=1024
var_energies_32x32=np.subtract(data_32x32[:,2],np.multiply(data_32x32[:,1],data_32x32[:,1]))
Heat_Capacity_32x32 = np.divide(var_energies_32x32, np.multiply(data_32x32[:,0], data_32x32[:,0]))
pl.plot(data_32x32[:,0],(Heat_Capacity_32x32),'mo-', markersize = 3, label = 'Lattice Size = 32x32')


pl.title('Plots of Heat Capacity vs Temperature for Lattice Sizes Tested')
pl.ylabel("Heat Capacity (in units of kB)")
pl.xlabel("Temperature")

pl.legend(loc='upper right')

pl.show()