#!/usr/bin/python
# syntax ./MC1D.py N T J h S
# N size, S steps

import sys
import numpy as np
import math

# Total energy
def Ham(sig,J,h,N):
    H=0
    for j in range(N):
       H+=-J*sig[j]*sig[(j+1)%N]-h*sig[j] # modulo "%" for periodic boundaries
    return H

# Local potential
def dE(sig,J,h,N,k):
    Vk=-J*(sig[(k-1)%N]+sig[(k+1)%N])*sig[k]-h*sig[k]
    return -2*Vk # flip dE=-2*V

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

# Note initial sigma is out of equilibrium
s0=np.random.choice((1,-1),N)
E0=Ham(s0,1,h,N)
sig=s0; E=E0;

for step in range(S):
    k=np.random.randint(N)
    dEk=dE(sig,J,h,N,k)
    Boltz=math.exp(-dEk/T)
    r=np.random.random()
    if(r < Boltz):
        sig[k]=-sig[k]
        E+=dEk
        M=sum(sig)

    print '{0:9d}  {1:11.0f}  {2:11d}'.format(step,E,M)