from visual import *
from random import random
from time import clock

scene.width = 900
scene.height = 500
scene.autoscale = 0
scene.up = (0,0,1)
scene.forward = (0,-1,0)

Ntotx = 200
Ntoty = 1
Ntotz = 1
Ntotal = Ntotx*Ntoty*Ntotz

spins = []

def lattice(N=3, delta=1.0):
    xmin = -Lx/2.
    nspin = 0
    ny = 0
    nz = 0
    
    for nx in range(Ntotx):
          x = xmin+nx*delta
          spins.append(frame())
          spins[-1].pos = vector(x,0,0)
          
          spins[-1].spin = vector()
          spins[-1].spinvec = arrow(pos=(x,0,0), axis=(0,0,1), color=(.9,.9,.9), length=2, shaftwidth=0.2)
          
          spins[-1].torque = vector()
          spins[-1].torquevec = arrow(pos=(x,0,0), axis=(0,0,1), color=(.3,.3,1), length=1, shaftwidth=0.05)
          
          spins[-1].near = range(2)
          
          spins[-1].nspin = nspin
          
          spins[-1].indices = (nx,ny,nz)
          
          nspin = nspin+1
          
    for s in spins:
          nx, ny, nz = s.indices
          if nx == 0: # left 'end'
              nspinl = nspin-1
              s.near[0] = nspinl
          else:
              nspinl = nx-1
              s.near[0] = nspinl
          if nx == Ntotx-1: # right 'end'
              nspinr = 0
              s.near[1] = nspinr
          else:
              nspinr = nx+1
              s.near[1] = nspinr
    return spins

dx = 1
Lx = dx*Ntotx

spins = lattice(N=Ntotx, delta=dx)

dt = .1
J = 1
k = 2.*(10)/Lx
spinvec_length = 2
#cons_torque_mag = 0
sigma = .4

# J = -1, k = 2.*(3)/Lx, sigma = .4, Nx = 80, results in clear illustration of spontaneous realignment


for s in spins:
    pos = s.pos
    
    spinx = sigma*cos(pi*k*pos.x)
    spiny = sigma*sin(pi*k*pos.x)
    spinz = 1*(J/abs(J))**s.nspin
    
    s.spin = vector(spinx,spiny,spinz)

    s.pos.z = s.pos.z + 0*(-1)**s.nspin
    
    s.spinvec.axis = s.spin*spinvec_length
    s.spinvec.pos = pos - s.spinvec.axis/2.


tt = clock()
Nsteps = 0
while 1:
    for s in spins:
        cons_spinz = s.spin.z
        s.torque = vector(0,0,0)

        for nn in range(2):
            nspin = s.near[nn]
            if nspin != None:
                s.torque = s.torque+cross(s.spin,spins[nspin].spin)

        s.torque.z = 0

#        if Nsteps == 0:
#            cons_torque_mag = mag(s.torque)

#        s.torque = cons_torque_mag*s.torque/mag(s.torque)
        
        s.spin = s.spin+J*s.torque*dt
        s.spin.z = 0
        s.spin = sigma*s.spin/mag(s.spin)
        s.spin.z = cons_spinz

        s.spin = s.spin/mag(s.spin)
        
        s.spinvec.axis = s.spin*spinvec_length
        s.spinvec.pos = s.pos - s.spinvec.axis/2.

        s.torquevec.axis = J*s.torque/mag(s.torque)
        s.torquevec.pos = s.pos + s.spinvec.axis/2.


    if Nsteps == 1000:
        tt = clock()-tt
        print '%0.1f' % tt, 'seconds for', Nsteps, 'steps with', Ntotal, 'spins'
    Nsteps = Nsteps+1


#Times:

#dt = .05,
#w/ torque vectors:   9.4 seconds for 1000 steps with 40 spins
#w/o torque vectors:  7.7 seconds for 1000 steps with 40 spins
