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

#the following 5 lines read the data from the files written in the previous section, assigning each column a variable
temps2, energies2, energysq2, magnetisations2, magnetisationsq2 = np.loadtxt('2x2.dat',unpack=True)
temps4, energies4, energysq4, magnetisations4, magnetisationsq4 = np.loadtxt('4x4.dat',unpack=True)
temps8, energies8, energysq8, magnetisations8, magnetisationsq8 = np.loadtxt('8x8.dat',unpack=True)
temps16, energies16, energysq16, magnetisations16, magnetisationsq16 = np.loadtxt('16x16.dat',unpack=True)
temps32, energies32, energysq32, magnetisations32, magnetisationsq32 = np.loadtxt('32x32.dat',unpack=True)
#the following lines create two subplots for the magnetisation and energy, defining different parts of the graph
fig = pl.figure()
enerax = fig.add_subplot(2,1,1)
enerax.set_ylabel("Energy per spin")
enerax.set_xlabel("Temperature")
enerax.set_ylim([-2.1, 0.1])
magax = fig.add_subplot(2,1,2)
magax.set_ylabel("Magnetisation per spin")
magax.set_xlabel("Temperature")
magax.set_ylim([-1.1, 1.1])
#the following section of code plots the data and adds the error bars as in the previous code
#all five data sets are plotted on the same graph
enerax.plot(temps2, np.array(energies2)/4,c='red',label='2x2')
en2_stand_dev = np.sqrt(energysq2 - np.multiply(energies2,energies2))
en2_error = np.divide(en2_stand_dev,np.sqrt(100000))
enerax.errorbar(temps2, np.array(energies2)/4,yerr=en2_error,ecolor='red',linestyle="None")
enerax.plot(temps4, np.array(energies4)/16,c='orange',label='4x4')
en4_stand_dev = np.sqrt(energysq4 - np.multiply(energies4,energies4))
en4_error = np.divide(en4_stand_dev,np.sqrt(100000))
enerax.errorbar(temps4, np.array(energies4)/16,yerr=en4_error,ecolor='orange',linestyle="None")
enerax.plot(temps8, np.array(energies8)/64,c='yellow',label='8x8')
en8_stand_dev = np.sqrt(energysq8 - np.multiply(energies8,energies8))
en8_error = np.divide(en8_stand_dev,np.sqrt(100000))
enerax.errorbar(temps8, np.array(energies8)/64,yerr=en8_error,ecolor='yellow',linestyle="None")
enerax.plot(temps16, np.array(energies16)/256,c='green',label='16x16')
en16_stand_dev = np.sqrt(energysq16 - np.multiply(energies16,energies16))
en16_error = np.divide(en16_stand_dev,np.sqrt(100000))
enerax.errorbar(temps16, np.array(energies16)/256,yerr=en16_error,ecolor='green',linestyle="None")
enerax.plot(temps32, np.array(energies32)/1024,c='blue',label='32x32')
en32_stand_dev = np.sqrt(energysq32 - np.multiply(energies32,energies32))
en32_error = np.divide(en32_stand_dev,np.sqrt(100000))
enerax.errorbar(temps32, np.array(energies32)/1024,yerr=en2_error,ecolor='blue',linestyle="None")
magax.plot(temps2, np.array(magnetisations2)/4,c='red',label='2x2')
mag2_stand_dev = np.sqrt(magnetisationsq2 - np.multiply(magnetisations2,magnetisations2))
mag2_error = np.divide(mag2_stand_dev,np.sqrt(100000))
magax.errorbar(temps2, np.array(magnetisations2)/4,yerr=mag2_error,ecolor='red',linestyle="None")
magax.plot(temps4, np.array(magnetisations4)/16,c='orange',label='4x4')
mag4_stand_dev = np.sqrt(magnetisationsq4 - np.multiply(magnetisations4,magnetisations4))
mag4_error = np.divide(mag4_stand_dev,np.sqrt(100000))
magax.errorbar(temps4, np.array(magnetisations4)/16,yerr=mag4_error,ecolor='orange',linestyle="None")
magax.plot(temps8, np.array(magnetisations8)/64,c='yellow',label='8x8')
mag8_stand_dev = np.sqrt(magnetisationsq8 - np.multiply(magnetisations8,magnetisations8))
mag8_error = np.divide(mag8_stand_dev,np.sqrt(100000))
magax.errorbar(temps8, np.array(magnetisations8)/64,yerr=mag8_error,ecolor='yellow',linestyle="None")
magax.plot(temps16, np.array(magnetisations16)/256,c='green',label='16x16')
mag16_stand_dev = np.sqrt(magnetisationsq16 - np.multiply(magnetisations16,magnetisations16))
mag16_error = np.divide(mag16_stand_dev,np.sqrt(100000))
magax.errorbar(temps16, np.array(magnetisations16)/256,yerr=mag16_error,ecolor='green',linestyle="None")
magax.plot(temps32, np.array(magnetisations32)/1024,c='blue',label='32x32')
mag32_stand_dev = np.sqrt(magnetisationsq32 - np.multiply(magnetisations32,magnetisations32))
mag32_error = np.divide(mag32_stand_dev,np.sqrt(100000))
magax.errorbar(temps32, np.array(magnetisations32)/1024,yerr=mag32_error,ecolor='blue',linestyle="None")
#this adds the legend outside the graph for each plot
enerax.legend(loc='center left',prop={'size':15},bbox_to_anchor=(1, 0.5))
magax.legend(loc='center left',prop={'size':15},bbox_to_anchor=(1, 0.5))
pl.show()