import numpy as np 

def ising(n, m, beta, J, B, S0):
    # Allocate vectors derived from a system's state
    E = np.zeros(n+1)
    M = np.zeros(n+1)

    # Initialize system 
    S = S0
    # Magnetic moment associated with the initial condition
    M[0] = S.sum()
    # Energy associated with the initial condition
    E[0] = - J * ( (S[:,:m-1] * S[:,1:]).sum() + (S[:m-1,:]*S[1:,:]).sum() ) - B * M[0]

    for k in range(n):
        # PROPOSAL : generate candidate state 
        Sc = S.copy() 
        # select randomly an atom on the lattice
        i = int( np.floor(m * np.random.random()) )
        j = int( np.floor(m * np.random.random()) )
        Sc[i,j] = - S[i,j] # flip the spin

        # Change in magnetic moment due to this flip 
        dM = - 2 * S[i,j]

        # Change in the energy due to this flip
        dE = 0
        if i>0:
                dE = dE + S[i-1,j]
        if i < m-1:
                dE = dE + S[i+1,j]
        if j>0:
                dE = dE + S[i,j-1]
        if j<m-1:
                dE = dE + S[i,j+1]

        dE = 2*S[i,j]*(J*dE + B)

        # ACCEPT-REJECT STEP
        alpha = np.min([np.exp(-dE*beta),1])
        U = np.random.random()
        if U < alpha:
                S = Sc # Accept proposed candidate state 
        else:
                dM = 0
                dE = 0

        # Update energy and magnetic moment 
        E[k+1] = E[k] + dE
        M[k+1] = M[k] + dM
            
    return E, M, S