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

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

#Get the Heat Capacity/ spin from 5th column
C = data_16x16_provided[:,5]  # Get Values for C-  Heat Capacity / spin

#first we fit the polynomial to the data
fit = np.polyfit(T, C, 30) # fit a 30 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 = 'Provided 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()