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 16x16 lattice
data_16x16_mine= np.loadtxt("16x16_mine.dat")

#Get the temperatures from the first column
T = data_16x16_mine[:,0] 

#Calculate the Heat Capacity from my Data

#Calculating Heat Capacity:

#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_16x16=np.subtract(data_16x16_mine[:,2],np.multiply(data_16x16_mine[:,1],data_16x16_mine[:,1]))

#Heat Capacity is equal to the Variance / T^2  
Heat_Capacity_16x16 = np.divide(var_energies_16x16, np.multiply(data_16x16_mine[:,0], data_16x16_mine[:,0]))

C = Heat_Capacity_16x16 /spins  # Get Values for C-  Heat Capacity / spin

#first we fit the polynomial to the data
fit = np.polyfit(T, C, 20) # fit a 20 order polynomial (this seems to give a good fit)

#now we generate interpolated values of the fitted polynomial over the range of our function
T_min = np.min(T)
T_max = np.max(T)
T_range = np.linspace(T_min, T_max, 1000) #generate 1000 evenly spaced points between T_min and T_max
fitted_C_values = np.polyval(fit, T_range) # use the fit object to generate the corresponding values of C



#Plotting data recorded through simulation
pl.plot(T, C,'bo-', markersize = 3, label = 'My Data: 16x16 Lattice')

#Plotting 'fitted' data
pl.plot(T_range,fitted_C_values,'ro-', markersize=3, label = 'Fitted curve: 16x16 Lattice') 

pl.title('Plot of Heat Capacity vs Temperature for Lattice Size: 16x16')
pl.ylabel("Heat Capacity (in units of kB)")
pl.xlabel("Temperature")

pl.legend(loc='best')


pl.show()