Rep:Mod:explosivesheep:Scripts

IsingLattice.py (unoptimized)

import numpy as np

class IsingLattice:

	E = 0.0
	E2 = 0.0
	M = 0.0
	M2 = 0.0

	n_cycles = 0

	def __init__(self, n_rows, n_cols):
		self.n_rows = n_rows
		self.n_cols = n_cols
		self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

	def energy(self):
		"Return the total energy of the current lattice configuration."
		energy = 0.0
		for i in len(self.lattice) - 1: #loop through indexes of all rows
			row = self.lattice[i]
			for j in len(row) - 1: #loop through indexes of all columns
				spin = self.lattice[i][j] #spin at the current position
				energy += spin * self.lattice[i][j-1] #multiply with its neighbour to the left
				# at the leftmost edge case, [i][-1] is the rightmost spin

				energy += spin * self.lattice[i-1][j] #multiply with its neighbour above
				# at the top edge case, [i-1][j] is the bottom spin

				# only multiply in one direction to negate the half factor
				# also preventing [j+1] or [i+1] which returns index errors
		return energy

	def magnetisation(self):
		"Return the total magnetisation of the current lattice configuration."
		magnetisation = 0
		for i in len(self.lattice) - 1: #loop through indexes of all rows
			for j in len(self.lattice[i]) - 1: #loop through indexes of all columns
				magnetisation += self.lattice[i][j]
				# loop through each spin and add to magnetisation
		return magnetisation

	def montecarlostep(self, T):
		# complete this function so that it performs a single Monte Carlo step
		# 1
		E_0 = self.energy()
	
		# 2 the following two lines will select the coordinates of the random spin for you
		random_i = np.random.choice(range(0, self.n_rows))
		random_j = np.random.choice(range(0, self.n_cols))
		self.lattice[random_i][random_j] = - self.lattice[random_i][random_j] #flips a spin at a random position
		
		# 3 calculate the energy difference
		delta_E = self.energy() - E_0
		print(delta_E)
		
		# 4 decide if the new config should be rejected
		if delta_E > 0.0: #if >0, coin flip; else do nothing and keep new configuration
			random_number = np.random.random() #choose a random %runnumber in the range[0,1) for you
			boltzmann_dist = np.exp(-delta_E/T)

			if random_number > boltzmann_dist: #if R>dist, reject new configuration; else do nothing and keep new configuration
				self.lattice[random_i][random_j] = - self.lattice[random_i][random_j] #reject new configuration by reversing the spin flip
		
		# 5 keep a sum of energies etc
		self.E += self.energy()
		self.E2 += self.energy()**2
		self.M += self.magnetisation()
		self.M2 += self.magnetisation()**2
		self.n_cycles += 1
		
		return self.energy(), self.magnetisation()

	def statistics(self):
		# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
		return self.E/self.n_cycles, self.E2/self.n_cycles, self.M/self.n_cycles, self.M2/self.n_cycles

Note: The ILanim.py provided only took 4 values from self.statistics(), omitting the number of cycles. But later on in other scripts given to us, a fifth value, n_cycles, had to be returned. The following optimized IsingLattice.py was changed based on this fact.

IsingLattice.py (optimized)

import numpy as np

class IsingLattice:

	E = 0.0
	E2 = 0.0
	M = 0.0
	M2 = 0.0

	n_cycles = 0
	equilibrium = 0

	def __init__(self, n_rows, n_cols):
		self.n_rows = n_rows
		self.n_cols = n_cols
		self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))
		self.equilibrium = n_rows * n_cols * 25
		#scale the equilibrium cutoff cycle according to the number of spins

	def energy(self):
		"Return the total energy of the current lattice configuration."
		energy = 0.0

		#shift the spins by moving down individual rows
		roll_row = np.roll(self.lattice, 1, axis=1) #axis 1 is a list of rows
		# the bottom row will be shifted to become the top row, satisfying periodic boundary conditions
		
		#shift the spins by moving individual columns to the right
		roll_col = np.roll(self.lattice, 1, axis=0) #axis 0 is a list of spins in a column
		# the rightmost spin will be shifted to become the leftmost spin, satisfying periodic boundary conditions
		
		#multiply the spins and sum the result
		#the 1/2 factor is negated as we avoided double counting by rolling in one direction only
		energy += -np.sum(np.multiply(self.lattice, roll_row)) #multiply neighbours' spins in horizontal direction and sum the energies
		energy += -np.sum(np.multiply(self.lattice, roll_col)) #multiply neighbours' spins in vertical direction and sum the energies
		return energy

	def magnetisation(self):
		"Return the total magnetisation of the current lattice configuration."
		return np.sum(self.lattice) #total magnetisation is sum of all spins in the lattice

	def montecarlostep(self, T):
		# complete this function so that it performs a single Monte Carlo step
		# Step 1: store the original energy
		E_0 = self.energy()
		
		# Step 2: the following two lines will select the coordinates of the random spin for you
		random_i = np.random.choice(range(0, self.n_rows))
		random_j = np.random.choice(range(0, self.n_cols))
		self.lattice[random_i][random_j] = - self.lattice[random_i][random_j] #apply the spin flip at random i,j
		
		# Step 3: calculate the energy difference
		delta_E = self.energy() - E_0 #energies are in units of Boltzmann constant
		# print(delta_E)
		
		# Step 4: decide if the new config should be rejected
		if delta_E > 0: #if spin is energetically unfavored, coin flip; else do nothing and keep new config
			random_number = np.random.random() #choose a random %runnumber in the range[0,1)
			exponential = np.exp(-delta_E/T) #calculates likelihood of spin flip by boltzmann function
			#if random number > boltmann distribution, reject new config; else do nothing and keep new config
			if random_number > exponential: 
				self.lattice[random_i][random_j] = -self.lattice[random_i][random_j] #reverse the spin flip
		
		# 5 keep a sum of energies etc. after equilibrium state is reached
		self.n_cycles += 1
		E = self.energy()
		M = self.magnetisation()

		if self.n_cycles > self.equilibrium: #check if it has reached a predetermined number of cycles
			self.E += E
			self.E2 += E**2
			self.M += M
			self.M2 += M**2
		
		return E, M

	def statistics(self):
		# complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them
		n = self.n_cycles - self.equilibrium
		return self.E/n, self.E2/n, self.M/n, self.M2/n, self.n_cycles

ILtemperaturerange.py

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

# Defining the Ising lattice
  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
  
  # Setting the number of Monte Carlo iterations
runtime = 150000
times = range(runtime)

# Setting the temperature range
  temps = np.arange(1.5, 5.0, 0.05) #start, end, increment
  
  # Lists to store statistics
  energies = []
  magnetisations = []
  energysq = []
  magnetisationsq = []
  
  for t in temps:
      for i in times:
          if i % 1000 == 0:
              print(str(i)+"th iteration at "+str(t)+"K")
        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)
      #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
  
  #Plot
  fig = pl.figure()
  enerax = fig.add_subplot(2,1,1) #y-axis on top figure
  enerax.set_ylabel("Energy per spin")
  enerax.set_xlabel("Temperature")
  enerax.set_ylim([-2.1, 2.1])
  magax = fig.add_subplot(2,1,2) #y-axis on bottom figure
  magax.set_ylabel("Magnetisation per spin")
  magax.set_xlabel("Temperature")
  magax.set_ylim([-1.1, 1.1])
  enerax.plot(temps, np.array(energies)/spins)
  magax.plot(temps, np.array(magnetisations)/spins)
  
  #Calculate errors: std error = square root of (variance/number of cycles) where variance = <E2> - <E>**2
energyerrors = np.sqrt((np.array(energysq)-np.array(energies)**2)/n_cycles)
magnetisationerrors = np.sqrt((np.array(magnetisationsq)-np.array(magnetisations)**2)/n_cycles)

#Plot error bars
  enerax.errorbar(temps, np.array(energies)/spins, yerr=energyerrors/2, linestyle='o') #energy error bars
 magax.errorbar(temps, np.array(magnetisations)/spins, yerr=magnetisationerrors/2, linestyle='o') #mag error bars
 
  pl.show()
  
  final_data = np.column_stack((temps, energies, energysq, magnetisations, magnetisationsq))
 file_name = str(n_rows) + "x" + str(n_cols) + ".dat" #dynamic file name
 np.savetxt(file_name, final_data)

IL.finalframe.py (modified)

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

#Initiate plot
 figure = pl.figure()
 subplotnumber = 0
 
 sizes = [7, 16, 32]
 #original setting is 7
 #16 and 32 used to find out equilibrium steps for larger system sizes
temps = [1, 2.69, 5]#temperatures set to below, around and above curie temperature
runtime = 50000 #should be enough to find equilibrium step, not interested in data atm

for temperature in temps:
	for n in sizes:
		il = IsingLattice(n,n)
		spins = n**2

		#subplot settings
		subplotnumber += 1
		# print(subplotnumber)
		cols = len(sizes) #
		rows = len(temps) * 3

		#Monte Carlo Simulation
		times = range(runtime)
		E=[]
		M=[]
		for i in times:
			if i %1000 == 0:
				print("Step: ", i, " Size:", n, " Temp: ", temperature)
			energy, magnetisation = il.montecarlostep(temperature)
			E.append(energy)
			M.append(magnetisation)

		#diagram subplot
		diag = figure.add_subplot(rows,cols,subplotnumber)
		diag.matshow(il.lattice)
		diag.tick_params(axis='both', which='major', labelsize=8)
		subplotnumber += cols
		# print(subplotnumber)

		#energy subplot
		enerax = figure.add_subplot(rows,cols,subplotnumber)
		enerax.set_ylabel(r"Energy per spin ($k_B$ Joules)", fontsize=8)
		enerax.set_xlabel("Monte Carlo Steps", fontsize=8)
		enerax.tick_params(axis='both', which='major', labelsize=8)
		enerax.plot(times, np.array(E)/spins)
		enerax.set_title(str(n) + " x " + str(n) + " at " + str(temperature) + " K")
		subplotnumber += cols
		# print(subplotnumber)

		#magentisation subplot
		magax = figure.add_subplot(rows,cols,subplotnumber)
		magax.set_ylabel("Magnetisation per spin", fontsize=8)
		magax.set_xlabel("Monte Carlo Steps", fontsize=8)
		magax.plot(times, np.array(M)/spins)
		magax.tick_params(axis='both', which='major', labelsize=8)
		subplotnumber -= 2*cols
		# print(subplotnumber)
	subplotnumber += 2*cols
	# print(subplotnumber)


pl.show()

ILsavedata.py

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

# Defining the Ising lattice
  for n in [2,4,8,16,32]: #looping the predetermined sizes
      il = IsingLattice(n, n)
      il.lattice = np.ones((n, n))
  
      # Setting the number of Monte Carlo iterations
    runtime = 150000
     times = range(runtime)
 
     # Setting the temperature range
     temps = np.arange(1.2, 5.0000001, 0.2) #start, end, increment

    energies = []
      magnetisations = []
      energysq = []
      magnetisationsq = []
  
      for t in temps:
          for i in times:
              if i % 1000 == 0:
                  print(str(i)+"th Monte Carlo cycle at "+str(t)+"K on a "+str(n)+"x"+str(n)+" lattice")
            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)
          #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
  
      #Calculate errors: std error = square root of variance/number of cycles
    energyerrors = np.sqrt((np.array(energysq)-np.array(energies)**2)/n_cycles)
    magnetisationerrors = np.sqrt((np.array(magnetisationsq)-np.array(magnetisations)**2)/n_cycles)

    final_data = np.column_stack((temps, energies, energysq, magnetisations, magnetisationsq))
    file_name = str(n) + "x" + str(n) + ".dat"
    np.savetxt(file_name, final_data)

ILplotsystemsize.py

from IsingLattice import *

from matplotlib import pylab as pl

import numpy as np



enfig = pl.subplot(6,2,1)

enfig.set_ylabel("Energy \n per spin ($k_B$)")

enfig.set_xlabel("Temperature")

enfig.set_title("Comparing Energies in Different System Sizes")

magfig = pl.subplot(6,2,2)

magfig.set_ylabel("Magnetisation \n per spin")

magfig.set_xlabel("Temperature")

magfig.set_title("Comparing Magnetisations in Different System Sizes")

subplotnumber = 3



# Defining the Ising lattice

for n in [2,4,8,16,32]:

	il = IsingLattice(n, n)

	il.lattice = np.ones((n, n))

	spins = n*n

	string = str(n) + "x" + str(n)

	


	array = np.loadtxt(string + ".dat") #load and store data

	temps = array[:,0]

	energies = array[:,1]

	energysq = array[:,2]

	magnetisations = array[:,3]

	magnetisationssq = array[:,4]



	#Plots

	#energy

	pl.subplot(6,2,subplotnumber) #y-axis on top figure

	pl.title(string)

	pl.ylabel("Energy per spin")

	pl.xlabel("Temperature")

	pl.plot(temps, np.array(energies)/spins)

	


	energyerrors = np.sqrt((np.array(energysq)-np.array(energies)**2)/n_cycles) * 0.5

	pl.errorbar(temps, np.array(energies)/spins, yerr=energyerrors, linestyle='o') #energy error bars

	subplotnumber += 1

	


	#plot on top graph for comparisons

	enfig.plot(temps, np.array(energies)/spins, label = string)

	enfig.errorbar(temps, np.array(energies)/spins, yerr=energyerrors, linestyle='o')



	#magnetisation

	pl.subplot(6,2,subplotnumber) #y-axis on bottom figure

	pl.title(string)

	pl.ylabel("Magnetisation \n per spin")

	pl.xlabel("Temperature")

	pl.plot(temps, np.array(magnetisations)/spins)

	


	magnetisationerrors = np.sqrt((np.array(magnetisationsq)-np.array(magnetisations)**2)/n_cycles) * 0.5

	pl.errorbar(temps, np.array(magnetisations)/spins, yerr=magnetisationerrors, linestyle='o') #mag error bars

	subplotnumber += 1

	#Calculate errors: std error = square root of (variance/number of cycles) where variance = <E2> - <E>**2



	#plot on top graph for comparisons

	magfig.plot(temps, np.array(magnetisations)/spins, label = string)

	magfig.errorbar(temps, np.array(magnetisations)/spins, yerr=magnetisationerrors, linestyle='o')

	


	final_data = np.column_stack((temps, energies, energysq, magnetisations, magnetisationsq))

    file_name = str(n) + "x" + str(n) + ".dat"

    np.savetxt(file_name, final_data)



magfig.legend(loc='upper right')

enfig.legend(loc='upper right')



pl.show()

ILcompare.py

from IsingLattice import *

from matplotlib import pylab as pl

from matplotlib import gridspec as gridspec

import numpy as np

	


#Define subplots

comparison_p = pl.subplot(4,2,1)

pl.title("A Comparison Of System Sizes (Python)")

pl.ylabel("Heat Capacity \n per spin ($k_B$)")

pl.xlabel("Temperature (K)")



comparison_c = pl.subplot(4,2,2)

pl.title("A Comparison Of System Sizes (C++)")

pl.ylabel("Heat Capacity \n per spin ($k_B$)")

pl.xlabel("Temperature (K)")



pl.rc('legend', **{'fontsize':9})



#Importing Data from C++ Files

subplotnumber = 3

for n in [2, 4, 8, 16, 32]:

	spins = n**2



	#Load Python data

	file_name = str(n) + "x" + str(n) + ".dat"

	array_p = np.loadtxt(file_name)

	temps_p = array_p[:,0]

	energies = array_p[:,1]

	energysq = array_p[:,2]

	#Calculate heat capacity per spins for array_p

	heatcap_p = (energysq - energies ** 2)/(temps_p**2 * spins)



	#Load C++ data

	array_c = np.loadtxt("C++_data/" + file_name)

	temps_c = array_c[:,0]

	heatcap_c = array_c[:,5]

	# array has these columns: T, E, E2, M, M2, C per spins



	#Plot temp v. heat capacity

	#on top comparison subplot

	pl.subplot(4,2,1)

	comparison_p.plot(temps_p, heatcap_p, label = file_name[:-4]) #add to comparison subplot

	pl.subplot(4,2,2)

	comparison_c.plot(temps_c, heatcap_c, label = file_name[:-4])



	#new subplot showing this system size only

	pl.subplot(4,2,subplotnumber) 

	subplotnumber +=1

	pl.plot(temps_c, heatcap_c, label = "C++ Simulation Data")

	pl.plot(temps_p, heatcap_p, label = "Python Simulation Data")

	pl.title(file_name[:-4])

	pl.grid(True)



	#Construct plot axes and legend

	pl.ylabel("Heat Capacity \n per spin ($k_B$)")

	pl.xlabel("Temperature (K)")

	pl.legend(loc="center left", bbox_to_anchor=(0.6, 0.5))



	#Set a reasonable scale

	pl.axis([min(temps_c), max(temps_c), 0, 2])



comparison_p.legend(loc="center left", bbox_to_anchor=(0.8, 0.5))

comparison_c.legend(loc="center left", bbox_to_anchor=(0.8, 0.5))



pl.show() #show the comparisons graph

ILplotheatcap.py

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

#Looping all the lattice sizes
 for n in [2,4,8,16,32]:
 	file_name = str(n) + "x" + str(n) + ".dat"
 	array = np.loadtxt(file_name) #load and store data
 	#array has these columns: temps, energies, energysq, magnetisations, magnetisationsq
 	spins = n**2
 
 	#Calculate heat capacity
 	pl.ylabel("Heat Capacity per spin")
 	pl.xlabel("Temperature (K)")
 	temps = array[:,0]
 	energies = array[:,1]
 	energysq = array[:,2]
 	heatcap = (energysq - energies**2)/(temps**2) #Var(E)/kT^2, where k is included in the units of heat cap
 
 	#plot heat capacity
 	pl.plot(temps, heatcap/spins, label=file_name[:-4])
 
 pl.legend(loc="center left", bbox_to_anchor=(1,0.5))
 
 pl.show()

ILtestfits.py

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

#a function to plot a given polynomial
   def PlotPolynomial(polynomial, x, description): #input actual x data and label of the curve
   	lin_x = np.linspace(np.min(x), np.max(x), 1000) #generate 1000 evenly spaced points between min and max of temps data given
   	lin_y = np.polyval(polynomial, lin_x) #generate a smooth fitted line using the best fit
   	pl.plot(lin_x, lin_y, label = description) #plot the fit on graph
   
   #a function to return statistical data and maxima of a fit
def AnalysePolynomial(polynomial, x, y): #input actual x and y data
	y_hat = np.polyval(polynomial, x)
	y_bar = np.sum(y)/len(y) #mean of observed data
	ssres = np.sum((y-y_hat)**2) #sum of sqaures of residuals, i.e. value between model and observed
	sstot = np.sum((y-y_bar)**2)  #sum of total variance in the observed data
	ssreg = np.sum((y_hat-y_bar)**2)

	if sstot != 0:
		R2 = 1 - ssres/sstot
	else:
		R2 = 0 #ssres/sstot must not be infinity since the minimum value of R2 is zero

	Extrema = np.roots(np.polyder(fit)) #returns a list of all extrema x values
	#Determine a suitable extrema around the curie temperature
	Curie = 0.0
	for point in Extrema:
		if np.imag(point) == 0 and np.real(point) > 1 and np.real(point) < 5.: #if root is real and in the range 1 < Tc < 5
			Curie = np.real(point) #sometimes numpy returns a real root with a very small imaginary component because how floats are handled in python
	return R2, ssreg, ssres, sstot, Curie	

#define the data set in the test fit
   n=64
   file_name = str(n) + "x" + str(n) + ".dat"
 array = np.loadtxt("C++_data/" + file_name)
 temps = array[:,0]
   heatcap = array[:,5]
   
   #Plot the simulation data
   pl.plot(temps, heatcap, label = "C++ Simulation Data")
 
   #loop through various degrees of polynomials
for deg in [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]:
	fit = np.polyfit(temps, heatcap, deg) #fit
	R2, ssreg, ssres, sstot, Curie = AnalysePolynomial(fit, temps, heatcap) #get coefficient of determiantion
	PlotPolynomial(fit, temps, "Test Fit, " + str(deg) + "th Polynomial \n R2:" + str(round(R2, 6))) #plot

#plot settings
   pl.rc('legend', **{'fontsize':10})
   pl.legend(loc="upper right")
   pl.ylabel("Heat Capacity per spin")
   pl.xlabel("Temperature")
   pl.show()

ILheatcapfit.py

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

#a function to return stats and a positive extrema of a numpy polynomial
def AnalysePolynomial(polynomial, x, y): #input actual x and y data
	y_hat = np.polyval(polynomial, x) # y values from model
	y_bar = np.sum(y)/len(y) #mean of observed data
	ssres = np.sum((y-y_hat)**2) #sum of sqaures of residuals, i.e. value between model and observed
	sstot = np.sum((y-y_bar)**2)  #sum of total variance in the observed data
	ssreg = np.sum((y_hat-y_bar)**2)  #sum of squares of residuals in a regression model

	if sstot != 0:
		R2 = 1 - ssres/sstot
	else:
		R2 = 0 #if ssres/sstot approaches infinity, then R2 is at its minimum value: zero

	Extrema = np.roots(np.polyder(fit)) #returns a list of all extrema x values
	#Determine a suitable extrema around the curie temperature
	for point in Extrema:
		if np.imag(point) == 0 and np.real(point) > 1 and np.real(point) < 5.: #if root is real and in the range 1 < Tc < 5
			Curie = np.real(point) #sometimes numpy returns a real root with a very small imaginary component because how floats are handled in python
	return R2, ssreg, ssres, sstot, Curie

#a function to plot a given polynomial
def PlotPolynomial(polynomial, x, description): #input actual x data and label of the curve
	lin_x = np.linspace(np.min(x), np.max(x), 1000) #generate 1000 evenly spaced points between min and max of temps data given
	lin_y = np.polyval(polynomial, lin_x) #generate a smooth fitted line using the best fit
	pl.plot(lin_x, lin_y, label = description) #plot the fit on graph

#a function to return a text string from the polynomial (as given by np.polyfit())
def ToText(polynomial, limit=9999999999999, dp=1): #unlimited length of text string and 1dp by default
	equation = ""
	for j in range(len(polynomial)):
		i = len(polynomial) - 1 - j #i is the power of x because coefficients are given in decreasing powers of x
		coefficient = polynomial[j]
		if coefficient >= 0.005: #displays non-zero coefficients only
			equation += " + " + str(round(coefficient, dp)) #round coefficient 
			equation += r"x^{" + str(i) + r"}" #index of the coeffficient is the power to which T is raised
		elif coefficient <= 0.005: #displays non-zero coefficients only
			equation += str(round(coefficient, dp)) #round coefficient
			equation += r"x^{" + str(i) + r"}" 
		if len(equation) >= limit: #use this to avoid equation taking up too much space
			equation += " + ..."
			break
	return equation
	
#Define subplots
comparison = pl.subplot(4,2,(1,2))
pl.title("A Comparison Of System Sizes Without Polynomial Fit")
pl.ylabel("Heat Capacity per spin")
pl.xlabel("Temperature")

forwiki = ""
inversesizes = []
curies = []
#Importing Data from C++ Files
subplotnumber = 3
for n in [2, 4, 8, 16, 32]:
	file_name = str(n) + "x" + str(n) + ".dat"
	array = np.loadtxt("C++_data/" + file_name)
	# array has these columns: T, E, E2, M, M2, C per spins

	#Plot temp v. heat capacity
	#on top comparison subplot
	pl.subplot(4,2,(1,2))
	temps = array[:,0]
	heatcap = array[:,5]
	comparison.plot(temps, heatcap, label = file_name[:-4]) #add to comparison subplot

	#new subplot showing this system size only
	pl.subplot(4,2,subplotnumber) 
	subplotnumber +=1
	pl.plot(temps, heatcap, label = "C++ Simulation Data")
	pl.title(file_name[:-4])
	pl.grid(True)

	#First General Fit Over All Region
	fit = np.polyfit(temps, heatcap, 270) #fit the data up to the 270th power
	R2, ssreg, ssres, sstot, Curie = AnalysePolynomial(fit, temps, heatcap) #analyse the fits's error and maxima
	PlotPolynomial(fit, temps, "Fit Over All Region")

	#Generate text on graphs
	pl.text(0.55, 0.5, r"Over All Region: $" + ToText(fit, limit=80) + r"$") #add equations text
	pl.text(0.55, 0.35, "Coefficient of Determiantion: " + str(R2)) #add the coefficient of determiantion
	pl.text(0.55, 0.2, "Maxima: " + "(" + str(Curie) + ", " + str(np.polyval(fit, Curie)) + ")") #add extrema points

	#Now we do another fit on the particular range around the estimated Curie temperature
	#We use the coefficient of determination to get the range of particular data for the more precise fit
	Tmin = Curie * (1 - (1-R2))
	Tmax = Curie * (1 + (1-R2))
	selection = np.logical_and(temps>Tmin, temps<Tmax)
	select_temps, select_heatcap = temps[selection], heatcap[selection]

	#The particular temperature range should include at least 10 values around the Curie temperature for a reasonable fit no matter how accurate the original fit was
	if len(select_temps) < 10:
		index = np.argmin(abs(temps - Curie)) #index of the temperature with lowest difference with Curie temp
		select_temps = temps[index-5:index+4] #select 10 closest temps around Curie
		select_heatcap = heatcap[index-5:index+4] #select 10 closest heat capacities around peak

	fit = np.polyfit(select_temps, select_heatcap, 2) #fit the data in the particular range as defined above
	R2, ssreg, ssres, sstot, Curie = AnalysePolynomial(fit, select_temps, select_heatcap)
	PlotPolynomial(fit, temps, "Fit Around Peak") #Plot the new particular fit

	#Generate text shown on graphs
	pl.text(0.55, 1.8, r"Around Peak: $" + ToText(fit, limit=80) + r"$") #add equations text
	pl.text(0.55, 1.65, "Coefficient of Determiantion: " + str(R2)) #add the coefficient of determiantion
	pl.text(0.55, 1.5, "Maxima: " + "(" + str(Curie) + ", " + str(np.polyval(fit, Curie)) + ")") #add extrema points
	pl.legend(loc="upper right")

	#Display texts for wiki writing
	forwiki += file_name[:-4] + ": " + ToText(fit, dp=6) + "\n R2: " + str(R2) + "\n Maxima: " + "(" + str(Curie) + ", " + str(np.polyval(fit, Curie)) + ") \n\n"

	#Append data for later graph plotting on the finite system size effect
	inversesizes.append(1/n)
	curies.append(Curie)

	#Construct plot axes
	pl.ylabel("Heat Capacity per spin")
	pl.xlabel("Temperature /K")

	#Set a reasonable scale
	pl.axis([min(temps), max(temps), 0, 2]) #same scale for all sizes for comparison purposes

#Add legend to top comparisons graph (legend must be generated after all data has been plotted to prevent matplotlib error)
pl.subplot(4,2,(1,2))
pl.legend(loc="center left", bbox_to_anchor=(0.85,0.5))
pl.grid(True)
pl.rc('legend', **{'fontsize':9})

#Plot the linear graph for the finite system size effect
pl.subplot(4,2,8)
pl.scatter(inversesizes, curies)
pl.title("Finite Size Effect")
pl.ylabel("Curie Temperature from Simulation (K)")
pl.xlabel(r"Lattice Size$^{-1}$")
linearfit = np.polyfit(inversesizes, curies, 1) #get a linear fit on the graph
PlotPolynomial(linearfit, inversesizes, None) #plot the fit
R2, ssreg, ssres, sstot, Curie = AnalysePolynomial(linearfit, inversesizes, curies) #get coefficient of determiantion

pl.text(-0.075, 2.50, r"$" + ToText(linearfit, limit=80, dp=6) + r"$") #add equations text
pl.text(-0.075, 2.485, "Coefficient of Determiantion: " + str(R2)) #add the coefficient of determiantion
forwiki += "Finite Size Effect: " + ToText(linearfit, dp=6) + "\n R2: " + str(R2) #display texts for wiki writing

print(forwiki)
pl.show()

To produce the graph in the final Conclusion, I sliced the Finite System Size Effect data:

#Plot the linear graph for the finite system size effect
......
curies = curies[2:]
inversesizes = inversesizes[2:]
linearfit = np.polyfit(inversesizes, curies, 1) #get a linear fit on the graph
PlotPolynomial(linearfit, inversesizes, None) #plot the fit
R2, ssreg, ssres, sstot, Curie = AnalysePolynomial(linearfit, inversesizes, curies) #get coefficient of determiantion

pl.text(-0.075, 2.50, r"$" + ToText(linearfit, limit=80, dp=6) + r"$") #add equations text
pl.text(-0.075, 2.485, "Coefficient of Determiantion: " + str(R2)) #add the coefficient of determiantion
forwiki += "Finite Size Effect: " + ToText(linearfit, dp=6) + "\n R2: " + str(R2) #display texts for wiki writing

print(forwiki)
pl.show()