""" From "COMPUTATIONAL PHYSICS" & "COMPUTER PROBLEMS in PHYSICS"
    by RH Landau, MJ Paez, and CC Bordeianu (deceased)
    Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, 
    C Bordeianu, Univ Bucharest, 2018. 
    Please respect copyright & acknowledge our work."""

# IsingVis.py: Ising model with Visual

"""Dirichlet boundary conditions surrounding four walls 
 Domain dimensions: WxH, with 2 triangles per square  
 Based on FEM2DL_Box Matlab program in Polycarpou, Intro Finite 
 Element Method in Electromagnetics, Morgan & Claypool (2006) """

from vpython import *
import numpy as np
import random

# Display for the arrows
scene = canvas(x=0, y=0, width=700, height=200,
               range=40, title='Spins')
engraph = graph(y=200, width=700, height=300,
                title='E of Spin System', xtitle='iteration',
                ytitle='E', xmax=500, xmin=0, ymax=5, ymin=-5)
enplot = gcurve(color=color.yellow)
N = 30
B = 1.
mu = .33                                      # g mu
J = .20
k = 1.                                   # Boltmann
T = 100.
state = np.zeros((N))               # spins up(1), down (0)
S = np.zeros((N), float)
test = state
random.seed()                           # Seed generator


def energy(S):
    FirstTerm = 0.
    SecondTerm = 0.
    for i in range(0, N - 2):
        FirstTerm += S[i] * S[i + 1]
    FirstTerm *= -J
    for i in range(0, N - 1):
        SecondTerm += S[i]
    SecondTerm *= -B * mu
    return (FirstTerm + SecondTerm)


ES = energy(state)


def spstate(state):                         # Plots spins
    for obj in scene.objects:
        obj.visible = 0  # Erase arrows
    j = 0
    for i in range(-N, N, 2):
        if state[j] == -1:
            ypos = 5             # Spin down
        else:
            ypos = 0
        if 5 * state[j] < 0:
            arrowcol = (1, 1, 1)  # White = down
        else:
            arrowcol = (0.7, 0.8, 0)
        arrow(pos=vector(i, ypos, 0),
              axis=vector(0, 5 * state[j], 0))
        j += 1


for i in range(0, N):
    state[i] = -1  # Initial spins all down

for obj in scene.objects:
    obj.visible = 0
spstate(state)
ES = energy(state)

for j in range(1, 500):
    rate(3)
    test = state
    r = int(N * random.random())   # Flip spin randomly
    test[r] *= -1
    ET = energy(test)
    p = np.exp((ES - ET) / (k * T))  # Boltzmann test
    enplot.plot(pos=vector(j, ES, 0))       # Adds segment to curve
    if p >= random.random():
        state = test
        spstate(state)
        ES = ET
