from numpy import *
from MatrixIO import *
from pylab import *

#reads in the data files with energy and magnetization and prints them as array
Energ=Readarr('data.txt')
Energ1=Readarr('data1.txt')
Energ2=Readarr('data2.txt')
Energ3=Readarr('data3.txt')
Energ4=Readarr('data4.txt')
Energ5=Readarr('data5.txt')
Energ6=Readarr('data6.txt')
Energ7=Readarr('data7.txt')
Energ8=Readarr('data8.txt')
Energ9=Readarr('data9.txt')
Energ10=Readarr('data10.txt')
Energ11=Readarr('data11.txt')

#reads in the final configuation of the system.
config=Readarr('stru.txt')
config1=Readarr('stru1.txt')
config2=Readarr('stru2.txt')
config3=Readarr('stru3.txt')
config4=Readarr('stru4.txt')
config5=Readarr('stru5.txt')
config6=Readarr('stru6.txt')
config7=Readarr('stru7.txt')
config8=Readarr('stru8.txt')
config9=Readarr('stru9.txt')
config10=Readarr('stru10.txt')
config11=Readarr('stru11.txt')

#Calculate the average magnetization after taking out the equilibration time

Av=[0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]
Av=[average(Energ[200:,1]),average(Energ1[200:,1]),average(Energ2[200:,1]),average(Energ3[200:,1]),average(Energ4[200:,1]),average(Energ5[200:,1]),average(Energ6[400:,1]),average(Energ7[200:,1]),average(Energ8[200:,1]),average(Energ9[200:,1]),average(Energ10[200:,1]),average(Energ11[200:,1])]

#Put temperature in a list
t=[0.5,1.0,1.5,2.0,2.1,2.2,2.269,2.4,2.5,2.6,2.8,3.0]
figure()
title('Magnetization per site vs Temperature')
xlabel('Temperature')
ylabel('Magnetisation')
plot(t,Av,'r')
show()

J=4.0
S=[abs(J*t[0]),abs(J*t[1]),abs(J*t[2]),abs(J*t[3]),abs(J*t[4]),abs(J*t[5]),abs(J*t[6]),abs(J*t[7]),abs(J*t[8]),abs(J*t[9]),abs(J*t[10]),abs(J*t[11])]

#number of lattice sites
N=100*100

#specific heat calculation
beta=[1./S[0],1./S[1],1./S[2],1./S[3],1./S[4],1./S[5],1./S[6],1./S[7],1./S[8],1./S[9],1./S[10],1./S[11]]

#calculate <E**2>
Mean_sqE=[mean(Energ[200:,0]**2),mean(Energ1[200:,0]**2),mean(Energ2[200:,0]**2),mean(Energ3[200:,0]**2),mean(Energ4[200:,0]**2),mean(Energ5[200:,0]**2),mean(Energ6[400:,0]**2),mean(Energ7[200:,0]**2),mean(Energ8[200:,0]**2),mean(Energ9[200:,0]**2),mean(Energ10[200:,0]**2),mean(Energ11[200:,0]**2)]

#calculate <E>
Mean_E=[mean(Energ[200:,0]),mean(Energ1[200:,0]),mean(Energ2[200:,0]),mean(Energ3[200:,0]),mean(Energ4[200:,0]),mean(Energ5[200:,0]),mean(Energ6[400:,0]),mean(Energ7[200:,0]),mean(Energ8[200:,0]),mean(Energ9[200:,0]),mean(Energ10[200:,0]),mean(Energ11[200:,0])]

#calculate beta**2*(<E2>-<E>**2
Spec_heat=[((beta[0]**2)*(Mean_sqE[0]-Mean_E[0]**2))/float(N),((beta[1]**2)*(Mean_sqE[1]-Mean_E[1]**2))/float(N),((beta[2]**2)*(Mean_sqE[2]-Mean_E[2]**2))/float(N),((beta[3]**2)*(Mean_sqE[3]-Mean_E[3]**2))/float(N),((beta[4]**2)*(Mean_sqE[4]-Mean_E[4]**2))/float(N),((beta[5]**2)*(Mean_sqE[5]-Mean_E[5]**2))/float(N),((beta[6]**2)*(Mean_sqE[6]-Mean_E[6]**2))/float(N),((beta[7]**2)*(Mean_sqE[7]-Mean_E[7]**2))/float(N),((beta[8]**2)*(Mean_sqE[8]-Mean_E[8]**2))/float(N),((beta[9]**2)*(Mean_sqE[9]-Mean_E[9]**2))/float(N),((beta[10]**2)*(Mean_sqE[10]-Mean_E[10]**2))/float(N),((beta[11]**2)*(Mean_sqE[11]-Mean_E[11]**2))/float(N)]


figure()
title('Specific heat per site vs Temperature')
xlabel('Temperature')
ylabel('Specific heat')
plot(t,Spec_heat,'g')
show()