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 4x4 lattice
data_4x4_provided= np.loadtxt("4x4.dat")

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

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

#first we fit the polynomial to the data


Tmin = 2.0 
Tmax = 3.0 #Setting Tmin and Tmax as desired
selection = np.logical_and(T > Tmin, T < Tmax) #choose only those rows where both conditions are true
peak_T_values = T[selection]
peak_C_values = C[selection]

fit = np.polyfit(peak_T_values, peak_C_values, 10) # fit a N order polynomial (this seems to give a good fit)

#now we generate interpolated values of the fitted polynomial over the chosen range
T_min = np.min(peak_T_values)
T_max = np.max(peak_T_values)
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

Cmax = np.max(fitted_C_values)
Tmax = T_range[fitted_C_values == Cmax]

print(Tmax)

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

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

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

pl.legend(loc='best')


pl.show()