#!/usr/bin/python3
#
# Plots 2d Ising model data

import numpy as np
import matplotlib.pyplot as plt

filename = "ising.dat"

data = np.transpose(np.loadtxt(filename))
Epoints = data[0]
E2points = data[1]
Mpoints = data[2]
M2points = data[3]
Tpoints = data[4]

chipoints =  256 * (M2points - Mpoints**2) / Tpoints
CVpoints = 256 * (E2points - Epoints**2) / Tpoints**2

plt.subplot(221)
plt.ylim(-2.2,0)
plt.xlabel('T/T0')
plt.ylabel('Energy density')
plt.plot(Tpoints, Epoints, 'bo')

plt.subplot(222)
plt.ylim(-1.2,1.2)
plt.xlabel('T/T0')
plt.ylabel('Magnetisation')
plt.plot(Tpoints, Mpoints, 'g^')

plt.subplot(223)
plt.ylim(-1,100)
plt.xlabel('T/T0')
plt.ylabel('chi')
plt.plot(Tpoints, chipoints, 'rv')

plt.subplot(224)
plt.ylim(-0.1,2)
plt.plot(Tpoints, CVpoints, 'm*')
plt.xlabel('T/T0')
plt.ylabel('Specific heat')

plt.show()