#!/usr/bin/python
# syntax ./Ising1D.py N T J h
import sys
import numpy as np
import math

def Ham(sig,J,h,N):
    H=0
    for j in range(N):
       H+=-J*sig[j]*sig[(j+1)%N]-h*sig[j] #%N for PBC
    return H


N=int(sys.argv[1]); T=float(sys.argv[2]); J=float(sys.argv[3]); h=float(sys.argv[4])

Z=0; M=0; M2=0; E=0; E2=0;

for k in range(2**N):
    s=[2*int(x)-1 for x in bin(k)[2:].zfill(N)]
    H=Ham(s,J,h,N)
    Boltz=math.exp(-H/T)
    Z+=Boltz
    M+=sum(s)*Boltz
    M2+=sum(s)**2*Boltz
    E+=H*Boltz
    E2+=H**2*Boltz
    
M=M/Z; M2=M2/Z; E=E/Z; E2=E2/Z;
c=(E2-E*E)/T**2
chi=(M2-M*M)/T
print "      E            C            M           chi"
print '{0:11.4f}  {1:11.4f}  {2:11.4f}  {3:11.4f}'.format(E,c,M,chi)
