from visual import *
from random import random

from Numeric import *
from mpmath import *

class NumericLattice:
    def __init__(self, nSpinsX, nSpinsY, stripeSpacingX, stripeSpacingY):
        self.nSpinsX = nSpinsX
        self.nSpinsY = nSpinsY
       
        self.stripeSpacingX = stripeSpacingX
        self.stripeSpacingY = stripeSpacingY

        self.Ja = -1.0
        self.Jbx = -1.0
        self.Jby = -1.0
        self.Jad = -1.0

        self.phiScale = 1.0
##        self.magS2 = 1.0
        self.temp = 0.001
        self.eta = .9  # damping

     
        #  Initializing Arrays                                              # to move to array language, make array, then operate with array
        self.zArray = zeros((stripeSpacingX,stripeSpacingY),Float32)        # scalars
        self.phaseArray = zeros((stripeSpacingX,stripeSpacingY),Float32)    # scalars
        self.phiArray = zeros((stripeSpacingX,stripeSpacingY),Float32)      # scalars

        #  Computation Arrays
        self.couplingXArray = zeros((nSpinsX,nSpinsX),Float32)              # coupling array for X direction (Up and Down), scalars
        self.couplingYArray = zeros((nSpinsY,nSpinsY),Float32)              # coupling array for Y direction (Left and Right), scalars

        self.spinArray = zeros((nSpinsX,nSpinsY,1,3),Float32)               # spin vectors in a matrix, vectors
        self.meanFieldArray = zeros((nSpinsX,nSpinsY,1,3),Float32)          # vectors
        self.torqueArray = zeros((nSpinsX,nSpinsY,1,3),Float32)             # torque vectors in a matrix, vectors
        
        self.tempArray = zeros((nSpinsX,nSpinsY,1,3),Float32)
        self.torqueTempArray = zeros((nSpinsX,nSpinsY,1,3),Float32)         # torque vectors in a matrix, vectors


    def setState(self, k, zArray, phaseArray, phiArray, phiScale, phiRatio, magS2):
        self.magS2 = magS2
        self.phiScale = phiScale
        self.phiRatio = phiRatio

        self.zArray = zArray
        self.phaseArray = phaseArray
        self.phiArray = phiArray

#        k2 = 2.0*pi*vector(0.0,1.0,0.0)

        for x in range(self.nSpinsX):
            for y in range(self.nSpinsY):
                xModXSpacing = x%self.stripeSpacingX
                yModYSpacing = y%self.stripeSpacingY
                
                xR = x - xModXSpacing
                yR = y - yModYSpacing               

                if(self.zArray[xModXSpacing][yModYSpacing] == 1.0):
                    tSigma = sin(self.phiScale*self.phiArray[xModXSpacing][yModYSpacing])                              
                    tzSigma = cos(self.phiScale*self.phiArray[xModXSpacing][yModYSpacing])
                else:
                    tSigma = sin(self.phiScale*self.phiArray[xModXSpacing][yModYSpacing])                              
                    tzSigma = cos(self.phiScale*self.phiArray[xModXSpacing][yModYSpacing])

                    
                tAngle = (k.x*x)/self.nSpinsX + (k.y*y)/self.nSpinsY #+ (k2.x*x)/self.nSpinsX + (k2.y*y)/self.nSpinsY
                self.phaseArray[xModXSpacing][yModYSpacing] = tAngle
                
                
                tSpinX = self.zArray[xModXSpacing][yModYSpacing]*tSigma*cos(tAngle)
                tSpinY = self.zArray[xModXSpacing][yModYSpacing]*tSigma*sin(tAngle)
                tSpinZ = self.zArray[xModXSpacing][yModYSpacing]*tzSigma

                self.spinArray[x][y] = (tSpinX,tSpinY,tSpinZ)

    def setSolitonState(self, k, zArray, phaseArray, phiArray, phiScale, phiRatio, magS2):
        self.magS2 = magS2
        self.phiScale = phiScale
        self.phiRatio = phiRatio

        self.zArray = zArray
        self.phaseArray = phaseArray
        self.phiArray = phiArray

        for x in range(self.nSpinsX):
            for y in range(self.nSpinsY):
                xModXSpacing = x%self.stripeSpacingX
                yModYSpacing = y%self.stripeSpacingY
                
                sign = -1+2*((x+y)%2)

                lf = 8.0
                lg = 8.0



                C1 = 1.6
                D1 = .9
                A1 = C1
                B1 = sqrt(C1*(2*D1-C1))

                C2 = 1.9
                D2 = .8
                A2 = C2
                B2 = -sqrt(C1)*(2*D2-C2)/sqrt(2*D1-C1)

                    
                if sign == 1:
                    tSpinX = A1*sech((x-self.nSpinsX/2.0)/lg)*tanh((x-self.nSpinsX/2.0)/lg)
                    tSpinY = B1*sech((x-self.nSpinsX/2.0)/lg)
                    tSpinZ = D1-C1*sech((x-self.nSpinsX/2.0)/lg)**2

                if sign == -1:               
                    tSpinX = A2*sech((x-self.nSpinsX/2.0)/lf)*tanh((x-self.nSpinsX/2.0)/lf)
                    tSpinY = B2*sech((x-self.nSpinsX/2.0)/lf)
                    tSpinZ = D2-C2*sech((x-self.nSpinsX/2.0)/lf)**2

##                vec = vector(float(tSpinX),float(tSpinY),float(tSpinZ))
##                print mag(vec)

##                if sign == -1:               
##                    tSpinX = -A1*sech((x-self.nSpinsX/3.0)/lf)*tanh((x-self.nSpinsX/3.0)/lf)-A2*sech((x-2*self.nSpinsX/3.0)/lf)*tanh((x-2*self.nSpinsX/3.0)/lf)
##                    tSpinY = -B1*sech((x-self.nSpinsX/3.0)/lf)+B2*sech((x-2*self.nSpinsX/3.0)/lf)
##                    tSpinZ = -(1-A1*sech((x-self.nSpinsX/3.0)/lf)**2)+A2*sech((x-2*self.nSpinsX/3.0)/lf)**2
##                    
##                if sign == 1:
##                    tSpinX = A1*sech((x-self.nSpinsX/3.0)/lg)*tanh((x-self.nSpinsX/3.0)/lg)+A2*sech((x-2*self.nSpinsX/3.0)/lg)*tanh((x-2*self.nSpinsX/3.0)/lg)
##                    tSpinY = -B1*sech((x-self.nSpinsX/3.0)/lg)+B2*sech((x-2*self.nSpinsX/3.0)/lg)
##                    tSpinZ = (1-A1*sech((x-self.nSpinsX/3.0)/lg)**2)-A2*sech((x-2*self.nSpinsX/3.0)/lg)**2

##                A = 1.6
##                B = .1
##                C = 1.8

##                if sign == -1:               
##                    tSpinX = -A*sech((x-self.nSpinsX/2.0)/lf)
##                    tSpinY = B*sech((x-self.nSpinsX/2.0)/lf)
##                    tSpinZ = -C*tanh((x-self.nSpinsX/2.0)/lf)
##                    
##                if sign == 1:
##                    tSpinX = A*sech((x-self.nSpinsX/2.0)/lg)
##                    tSpinY = B*sech((x-self.nSpinsX/2.0)/lg)
##                    tSpinZ = C*tanh((x-self.nSpinsX/2.0)/lg)

##
##                if sign == -1:               
##                    tSpinX = -A*sech((x-self.nSpinsX/2.0)/lf)*tanh((x-self.nSpinsX/2.0)/lf)
##                    tSpinY = B*sech((x-self.nSpinsX/2.0)/lf)
##                    tSpinZ = -(1-C*sech((x-self.nSpinsX/2.0)/lf)**2)
##                    
##                if sign == 1:
##                    tSpinX = A*sech((x-self.nSpinsX/2.0)/lg)*tanh((x-self.nSpinsX/2.0)/lg)
##                    tSpinY = B*sech((x-self.nSpinsX/2.0)/lg)
##                    tSpinZ = (1-C*sech((x-self.nSpinsX/2.0)/lg)**2)
##
##                if x == self.nSpinsX/2.0+4:
##                    print tSpinX
##                if x == self.nSpinsX/2.0:
##                    print tSpinY
##                    print tSpinZ

                


##                g = .40
##               
##                lg = 4.5*g
##                lf = 8.0*g
##                la = 10.0*g

##                A = .8
##                B = .2
##
##                lg = 4.0
##                lf = 4.0
##                la = 4.0
##
##                offsetY =  0.0
                
##                tSpinY = -B*sech((x-self.nSpinsX/2.0+offsetY)/lf)
##                tSpinZ = -sign*(1-2*sech((x-self.nSpinsX/2.0)/la)**2)
##                tSpinX = ((x-self.nSpinsX/2.0)/abs(x-self.nSpinsX/2.0+.001))*sign*sqrt(1-tSpinZ**2)

##                tSpinX = sign*A*sech((x-self.nSpinsX/2.0)/lf)
##                tSpinY = sqrt(1-A**2)*sech((x-self.nSpinsX/2.0)/lf)
##                tSpinZ = sign*tanh((x-self.nSpinsX/2.0)/la)

##                fx = exp(-abs(x-self.nSpinsX/2.0)/lf)
##                gx = exp(-abs(x-self.nSpinsX/2.0)/lg)
##                arg = (fx**2-3*gx**2)/(fx**2+gx**2)*1/gamma
##                theta = acos(arg)
##                theta = 2*3.1415926535/(self.nSpinsX-1.0)

##                Lambda = -1/(lf*lg)
##                b = -Lambda
##                gamma = 1+b

##                if sign == 1:
##                    tSpinX = exp(-abs(x-self.nSpinsX/2.0)/lf)*cos(theta)-exp(-abs(x-self.nSpinsX/2.0)/lg)*sin(theta)
##                    tSpinY = exp(-abs(x-self.nSpinsX/2.0)/lg)*cos(theta)+exp(-abs(x-self.nSpinsX/2.0)/lf)*sin(theta)
##                    tSpinZ = sqrt(1-exp(-2*abs(x-self.nSpinsX/2.0)/lg))
##                    
##                if sign == -1:
##                    tSpinX = sign*(exp(-abs(x-self.nSpinsX/2.0)/lf)*cos(theta)-exp(-abs(x-self.nSpinsX/2.0)/lg)*sin(theta))
##                    tSpinY = exp(-abs(x-self.nSpinsX/2.0)/lg)*cos(theta)+exp(-abs(x-self.nSpinsX/2.0)/lf)*sin(theta)
##                    tSpinZ = sign*sqrt(1-exp(-2*abs(x-self.nSpinsX/2.0)/lg))
                               
##                if sign == 1:
##                    tSpinX = A*exp(-abs(x-self.nSpinsX/2.0)/lf)*cos(theta)-B*exp(-abs(x-self.nSpinsX/2.0)/lg)*sin(theta)
##                    tSpinY = B*exp(-abs(x-self.nSpinsX/2.0)/lf)*sin(theta)+B*exp(-abs(x-self.nSpinsX/2.0)/lg)*cos(theta)
##                    tSpinZ = sqrt(1-B**2*exp(-2*abs(x-self.nSpinsX/2.0)/lg)-A**2*exp(-2*abs(x-self.nSpinsX/2.0)/lf))
##
##                if sign == -1:
##                    tSpinX = sign*(A*exp(-abs(x-self.nSpinsX/2.0)/lf)*cos(theta)-B*exp(-abs(x-self.nSpinsX/2.0)/lg)*sin(theta))
##                    tSpinY = sign*(B*exp(-abs(x-self.nSpinsX/2.0)/lf)*sin(theta)+B*exp(-abs(x-self.nSpinsX/2.0)/lg)*cos(theta))
##                    tSpinZ = sign*sqrt(1-B**2*exp(-2*abs(x-self.nSpinsX/2.0)/lg)-A**2*exp(-2*abs(x-self.nSpinsX/2.0)/lf))

##                tSpinX = sign*(A*sech((x-self.nSpinsX/2.0)/lf)*cos(theta)-B*sech((x-self.nSpinsX/2.0)/lf)*sin(theta))
##                tSpinY = (B*sech((x-self.nSpinsX/2.0)/lf)*cos(theta)+A*sech((x-self.nSpinsX/2.0)/lf)*sin(theta))
##                tSpinZ = sign*tanh((x-self.nSpinsX/2.0)/lf)

##      twist soliton
##                tSpinX = sign*(A*sech((x-self.nSpinsX/2.0)/lf)*cos(theta)-B*sech((x-self.nSpinsX/2.0)/lf)*sin(theta))
##                tSpinY = sign*(B*sech((x-self.nSpinsX/2.0)/lf)*cos(theta)+A*sech((x-self.nSpinsX/2.0)/lf)*sin(theta))
##                tSpinZ = sign*sqrt(1-tSpinX**2-tSpinY**2)

##      twist with pb  
##                tSpinX = Gamma*sign*sech((x-self.nSpinsX/2.0)/lf)*sin(theta)
##                tSpinY = Gamma*sign*sech((x-self.nSpinsX/2.0)/lf)*cos(theta)
##                tSpinZ = sign*tanh(abs(x-self.nSpinsX/2.0)/lf)

##      twist soliton
##                tSpinX = sign*sech((x-self.nSpinsX/2.0)/lg)*tanh((x-self.nSpinsX/2.0)/lg)
##                tSpinZ = sign*(1-2*sech((x-self.nSpinsX/2.0)/lf))
##                tSpinY = sign*sqrt(1-tSpinX**2-tSpinZ**2)

##      k =2.0
##                tSpinX = sign*.5*(cos((x-self.nSpinsX/2.0)*k)-.5*cos(2*(x-self.nSpinsX/2.0)*k))
##                tSpinY = sign*.5*(sin((x-self.nSpinsX/2.0)*k)-.5*sin(2*(x-self.nSpinsX/2.0)*k))
##                tSpinZ = sign*sqrt(1-tSpinX**2-tSpinY**2)

##      excellent twist (tacoma narrows w/o pb, complete pendular swinging)
##                tSpinX = sign*sech((x-self.nSpinsX/2.0)/lf)*sin(theta)
##                tSpinY = sign*sech((x-self.nSpinsX/2.0)/lf)*cos(theta)
##                tSpinZ = sign*tanh((x-self.nSpinsX/2.0)/lf)

##      new soliton w/o pb
##                tSpinX = sign*sech((x-self.nSpinsX/2.0)/lg)
##                tSpinY = sign*sech((x-self.nSpinsX/2.0)/lg)*tanh((x-self.nSpinsX/2.0)/lf)
##                tSpinZ = sign*tanh((x-self.nSpinsX/2.0)/lf)
               
##
##                tSpinX = sign*A*exp(-abs(x-self.nSpinsX/2.0)/lf)*cos(theta)-B*exp(-abs(x-self.nSpinsX/2.0)/lg)*sin(theta)
##                tSpinY = B*exp(-abs(x-self.nSpinsX/2.0)/lg)*cos(theta)+A*sign*exp(-abs(x-self.nSpinsX/2.0)/lf)*sin(theta)
##                tSpinZ = sign*sqrt(1-B**2*exp(-2*abs(x-self.nSpinsX/2.0)/lg)-A**2*exp(-2*abs(x-self.nSpinsX/2.0)/lf))

##                if sign == 1:
##                    tSpinX = A*sech((x-6.0)/la)*cos(phi)
##                    tSpinY = A*sech((x-6.0)/la)*sin(phi)
##                    tSpinZ = tanh((x-6.0)/la)
##                    
##                if sign == -1:
##                    tSpinX = sign*A*sech((x-6.0)/la)*cos(phi)
##                    tSpinY = sign*A*sech((x-6.0)/la)*sin(phi)
##                    tSpinZ = sign*tanh((x-6.0)/la)

##                if sign == 1:
##                    tSpinX = 1/sqrt(1+(x-self.nSpinsX/2.0)**2/la**2)
##                    tSpinY = 0.000
##                    tSpinZ = sqrt(1-1/(1+(x-self.nSpinsX/2.0)**2/la**2))
##                    
##                if sign == -1:
##                    tSpinX = sign*1/sqrt(1+(x-self.nSpinsX/2.0)**2/la**2)
##                    tSpinY = 0.000
##                    tSpinZ = sign*sqrt(1-1/(1+(x-self.nSpinsX/2.0)**2/la**2))

##                if sign == 1:
##                    tSpinX = (x-5.99)/abs(x-5.99)*sech((x-6.0)/la)*cos(phi)
##                    tSpinY = (x-5.99)/abs(x-5.99)*sech((x-6.0)/la)*sin(phi)
##                    tSpinZ = tanh((x-6.0)/la)
##                    
##                if sign == -1:
##                    tSpinX = (x-5.99)/abs(x-5.99)*sech((x-6.0)/la)*cos(phi)
##                    tSpinY = sign*(x-5.99)/abs(x-5.99)*sech((x-6.0)/la)*sin(phi)
##                    tSpinZ = sign*tanh((x-6.0)/la)

##                tSpinX = sech((x-4.0)/la)
##                tSpinY = tanh((x-4.0)/la)
##                tSpinZ = sech((x-4.0)/la)


                self.spinArray[x][y] = (tSpinX,tSpinY,tSpinZ)

##                print self.spinArray[x][y]

    def setCouplings(self, Ja, Jbx, Jby, Jad):          # sets the couplings of neighboring spins and creates two coupling arrays.
        self.Ja = Ja                                    # Unit Cell interior region coupling. Unit Cell Object specific. Currently a constant.
        self.Jbx = Jbx                                  # Unit Cell boundary coupling., between unit cells, in the x direction.
        self.Jby = Jby                                  # in the y direction...
        self.Jad = Jad                                  # separation of couplings between middle spins.

        self.couplingXArray = zeros((self.nSpinsX,self.nSpinsX),Float32)        # resets to all zeros
        self.couplingYArray = zeros((self.nSpinsY,self.nSpinsY),Float32)

        for i in range(self.nSpinsX):
            for j in range(self.nSpinsX):
                if      i - j == 1:
                    if      i%self.stripeSpacingX == 0 and j%self.stripeSpacingX == (self.stripeSpacingX-1):    self.couplingXArray[i][j] += self.Jbx
                    else:                                                                                       self.couplingXArray[i][j] += self.Ja 
                if      j - i == 1:
                    if      j%self.stripeSpacingX == 0 and i%self.stripeSpacingX == (self.stripeSpacingX-1):    self.couplingXArray[i][j] += self.Jbx
                    else:                                                                                       self.couplingXArray[i][j] += self.Ja
                if      i - j == self.nSpinsX - 1 or j - i == self.nSpinsX - 1:                                 self.couplingXArray[i][j] += self.Jbx

        i=0
        j=0

        for i in range(self.nSpinsY):
            for j in range(self.nSpinsY):
                if      i - j == 1:
                    if      i%self.stripeSpacingY == 0 and j%self.stripeSpacingY == (self.stripeSpacingY-1):    self.couplingYArray[i][j] += self.Jby
                    else:                                                                                       self.couplingYArray[i][j] += self.Ja
                if      j - i == 1:
                    if      j%self.stripeSpacingY == 0 and i%self.stripeSpacingY == (self.stripeSpacingY-1):    self.couplingYArray[i][j] += self.Jby
                    else:                                                                                       self.couplingYArray[i][j] += self.Ja
                if      i - j == self.nSpinsY - 1 or j - i == self.nSpinsY - 1:                                 self.couplingYArray[i][j] += self.Jby

##        self.couplingXArray[0][self.nSpinsX-1] = 0
##        self.couplingXArray[self.nSpinsX-1][0] = 0
        # periodic boundary conditions in the x direction


    def getTorques(self, spinArray, t):
        transSpinArray = transpose(spinArray)                                             # array preparation for coupling computation
        twiceTransSpinArray = transpose(transpose(spinArray,(1,0,2,3)))                   #

        meanFieldXArray = matrixmultiply(transSpinArray,self.couplingXArray)                   # computation of mean field
        meanFieldYArray = matrixmultiply(twiceTransSpinArray,self.couplingYArray)              #
      
        self.meanFieldArray = transpose(meanFieldXArray) + transpose(meanFieldYArray, (2,3,1,0))    # sum of x and y directions in correct shape
        
        for x in range(self.nSpinsX):
            for y in range(self.nSpinsY):
                self.torqueArray[x][y] = -cross(spinArray[x][y][0],self.meanFieldArray[x][y][0])

##        self.torqueArray[0][0] = vector(0.0,0.0,0.0)
##        self.torqueArray[-1][0] = vector(0.0,0.0,0.0)
##        self.torqueArray[0][1] = vector(0.0,0.0,0.0)
##        self.torqueArray[-1][1] = vector(0.0,0.0,0.0)
        ## + vector(0.0,0.0,2.0)


    def timeEvolve(self, t, dt, temp, EBu, EBd):                      # Takes us one step forward in time by adding the torque to our spins and renormalizing.
        self.temp = temp                                    # adds a random torque <--> kick to each spin, with a scale set by temp.
                           
##        for x in range(self.nSpinsX):
##            for y in range(self.nSpinsY):
##                randomU = random()                          # perturbed by unit ball
##                randomV = random()                          # we produce a uniform distribution on the unit ball
##                theta = 2.0*3.1416*randomU
##                phi = arccos(2.0*randomV-1.0)
##                kick = self.temp*norm(vector(sin(theta)*cos(phi),cos(theta)*sin(phi),cos(phi)))
##                            
##                if x+y == 0:
##                    self.tempArray[x][y][0] = kick        # this kick is constant in magnitude but random in direction
##                if x+y == 1:
##                    self.tempArray[x][y][0] = kick        # this kick is constant in magnitude but random in direction

##        self.Jbx = self.Jbx - .0001
##        self.setCouplings(self.Ja, self.Jbx, self.Jby, self.Jad)
        

        for x in range(self.nSpinsX):
            for y in range(self.nSpinsY):
                self.torqueTempArray[x][y] = -cross(self.spinArray[x][y][0],self.tempArray[x][y][0])

        k0 = dt*self.torqueTempArray

        self.getTorques(self.spinArray, t)              # for time dep solns: set t = tn
        k1 = dt*self.torqueArray
        self.getTorques(self.spinArray + k1/2.0, t + dt/2.0)        # for time dep solns: set t = tn + dt/2.0
        k2 = dt*self.torqueArray
        self.getTorques(self.spinArray + k2/2.0, t + dt/2.0)        # for time dep solns: set t = tn + dt/2.0
        k3 = dt*self.torqueArray
        self.getTorques(self.spinArray + k3, t + dt)                # for time dep solns: set t = tn + dt
        k4 = dt*self.torqueArray


        self.spinArray = self.spinArray + k0 + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0
        
##        print (self.spinArray[0][0][0][0],self.spinArray[0][0][0][1],self.spinArray[0][0][0][2])
##        print (self.spinArray[1][0][0][0],self.spinArray[1][0][0][1],self.spinArray[1][0][0][2])
##        print (self.spinArray[2][0][0][0],self.spinArray[2][0][0][1],self.spinArray[2][0][0][2])
##        print (self.spinArray[0][1][0][0],self.spinArray[0][1][0][1],self.spinArray[0][1][0][2])
##        print (self.spinArray[1][1][0][0],self.spinArray[1][1][0][1],self.spinArray[1][1][0][2])
##        print (self.spinArray[2][1][0][0],self.spinArray[2][1][0][1],self.spinArray[2][1][0][2])
##        print ""
##        print self.Jbx
##        print ""

##        print dot(self.spinArray[1][0][0],self.meanFieldArray[1][0][0])
##        print self.spinArray[0][0][0]-self.spinArray[1][0][0]+self.spinArray[2][0][0]

##        print self.spinArray[0][0][0][1]+self.spinArray[1][0][0][1]+self.spinArray[2][0][0][1]
##        print dot(self.spinArray[0][0][0],self.meanFieldArray[0][0][0])+dot(self.spinArray[1][0][0],self.meanFieldArray[1][0][0])+dot(self.spinArray[2][0][0],self.meanFieldArray[2][0][0])

##        print self.returnTotalEnergy()
        
##        print self.returnSpinYSum()

        for x in range(self.nSpinsX):           # normalizes the state
            for y in range(self.nSpinsY):
##                if x == 1:
##                    self.spinArray[x][y][0] = self.magS2*norm(self.spinArray[x][y][0])
##                else:
                self.spinArray[x][y][0] = norm(self.spinArray[x][y][0])




    # Data for Export Methods
   
    def returnState(self):
        return self.spinArray

    def returnPhaseArray(self):
        return self.phaseArray

    def returnPhiArray(self):
        return self.phiArray

    def returnzArray(self):
        return self.zArray

    def returnMeanFieldArray(self):
        return self.meanFieldArray

    def returnTorqueArray(self):
        return self.torqueArray

    def returnPhiScale(self):
        return self.phiScale

    def returnSpinYSum(self):
        spinYSum = 0.0
        for x in range(self.nSpinsX):
            for y in range(self.nSpinsY):         
                spinYSum += self.spinArray[x][y][0][1]
        return spinYSum
    
    def returnSpinSum(self):
        spinSum = vector(0.0,0.0,0.0)
        for x in range(self.nSpinsX):
            for y in range(self.nSpinsY):
                spinSum += self.spinArray[x][y][0]
        return spinSum

    def returnEnergies(self):
        energyArray = zeros((self.nSpinsX,self.nSpinsY),Float32)                 # scalars
        for x in range(self.nSpinsX):
            for y in range(self.nSpinsY):
                energyArray[x][y] = -dot(self.spinArray[x][y][0],self.meanFieldArray[x][y][0])
        return energyArray

    def returnBondXEnergies(self):
        spinArray = self.returnState()
        
        energyArrayX = zeros((self.nSpinsX+1,self.nSpinsY),Float32)
        for x in range(self.nSpinsX+1):
            for y in range(self.nSpinsY):
                energyArrayX[x][y] = -dot(spinArray[x%self.nSpinsX][y%self.nSpinsY][0],spinArray[(x-1)%self.nSpinsX][y%self.nSpinsY][0])

        return energyArrayX

    def returnBondYEnergies(self):
        spinArray = self.returnState()
        
        energyArrayY = zeros((self.nSpinsX,self.nSpinsY+1),Float32)
        for x in range(self.nSpinsX):
            for y in range(self.nSpinsY+1):              
                energyArrayY[x][y] = -dot(spinArray[x%self.nSpinsX][y%self.nSpinsY][0],spinArray[x%self.nSpinsX][(y-1)%self.nSpinsY][0])
        return energyArrayY

    def returnTotalEnergy(self):
        energyArray = self.returnEnergies()
        totEnergy = 0.0
        for x in range(self.nSpinsX):
            for y in range(self.nSpinsY):
                totEnergy = totEnergy + energyArray[x][y]
        return totEnergy/2.0

    def returnCouplingX(self):
        return self.couplingXArray

    def returnCouplingY(self):
        return self.couplingYArray

    def returnStripeSpacingX(self):
        return self.stripeSpacingX

    def returnStripeSpacingY(self):
        return self.stripeSpacingY









    # Perturbation Methods


    def randomizeState(self):
        for x in range(self.nSpinsX):
            for y in range(self.nSpinsY):
                randomU = random()
                randomV = random()
                
                theta = 2.0*3.1416*randomU
                phi = arccos(2.0*randomV-1.0)
                
                self.spinArray[x][y][0] = norm(vector(sin(theta)*cos(phi),cos(theta)*sin(phi),cos(phi)))





























##        self.temp = 0.00
##        self.EBu = 0.00
##        self.EBd = 0.00
##        self.EBaxis = vector(0.0,0.0,1.0)


##        w = (30.0/10.0)*2.0*pi
##        A = .7
##        arg = A*sin(w*t)
##        EField = self.EBu*vector(0.0,0.0,cos(w*t))
##        appliedFieldArray = [[EField],[EField],[EField]]

                
##                if self.temp != 0.0:                            # adds a random torque <--> kick to each spin, with a scale set by temp.
##                    randomU = random()                          # perturbed by unit ball
##                    randomV = random()                          # we produce a uniform distribution on the unit ball
##                    theta = 2.0*3.1416*randomU
##                    phi = arccos(2.0*randomV-1.0)
##                    kick = self.temp*norm(vector(sin(theta)*cos(phi),cos(theta)*sin(phi),cos(phi)))
##                    
##                    self.meanFieldArray[x][y] += kick        # this kick is constant in magnitude but random in direction

##                if self.EBu != 0.0 or self.EBd != 0.0:
##                    self.meanFieldArray[x][y] += appliedFieldArray[x][y]



    

##    def setStateManually(self, k, baseSigma, zArray, sigmaArray, phaseArray):
##        for x in range(self.nSpinsX):
##            for y in range(self.nSpinsY):
##                xModXSpacing = x%self.stripeSpacingX
##                yModYSpacing = y%self.stripeSpacingY
##                xR = x - xModXSpacing
##                yR = y - yModYSpacing
##
##                z = zArray[xModXSpacing][yModYSpacing]
##
##                tempSigma = baseSigma*sigmaArray[xModXSpacing][yModYSpacing]
##                tempAngle = phaseArray[xModXSpacing][yModYSpacing] + k.x*xR + k.y*yR
##                
##                tempSpinX = tempSigma*cos(tempAngle)
##                tempSpinY = tempSigma*sin(tempAngle)
##                tempSpinZ = z*(1.0-tempSigma**2)**.5
##               
##                self.spinArray[x][y] = (tempSpinX,tempSpinY,tempSpinZ)
##
##
##

##
##        # Output Arrays
##        self.energyArray = zeros((nSpinsX,nSpinsY),Float32)                 # scalars
##        self.maxEnergyArray = -4.0*ones((nSpinsX,nSpinsY),Float32)
##        self.minEnergyArray = zeros((nSpinsX,nSpinsY),Float32)
##
##        self.TrackSpins = TrackSpins




##    def setStateFromH(self, k, baseSigma, n):
##        self.Ham.makeHamiltonian(k)
##
##        if n == -1:     state = self.Ham.getGroundState()
##        else:           state = self.Ham.getState(n)
##
##        xComponentUCArray = state.real
##        yComponentUCArray = state.imaginary
##
##        for x in range(self.stripeSpacingX):
##            for y in range(self.stripeSpacingY):
##                z = (-1.0)**(x%2 + y%2)
##                self.phaseArray[x][y] = k.x*x + k.y*y + (1-z)*pi/2
##
##        for x in range(self.nSpinsX):
##            for y in range(self.nSpinsY):
##                xModXSpacing = x%self.stripeSpacingX
##                yModYSpacing = y%self.stripeSpacingY
##                xR = x - xModXSpacing
##                yR = y - yModYSpacing
##
##                z = (-1.0)**(xModXSpacing%2 + yModYSpacing%2)
##
##                self.sigmaArray[xModXSpacing][yModYSpacing] = baseSigma*(xComponentUCArray[xModXSpacing][yModYSpacing]**2 + yComponentUCArray[xModXSpacing][yModYSpacing]**2)**0.5
##                self.phaseArray[xModXSpacing][yModYSpacing] = k.x*xModXSpacing + k.y*yModYSpacing + (1-z)*pi/2
##
##                tempSigma = self.sigmaArray[xModXSpacing][yModYSpacing]
##                tempAngle = self.phaseArray[xModXSpacing][yModYSpacing] + k.x*xR + k.y*yR
##
##                tempSpinX = baseSigma*(xComponentUCArray[xModXSpacing][yModYSpacing]*cos(tempAngle) - yComponentUCArray[xModXSpacing][yModYSpacing]*sin(tempAngle))
####                tempSpinX = baseSigma*xComponentUCArray[xModXSpacing][yModYSpacing]
##                tempSpinY = baseSigma*(xComponentUCArray[xModXSpacing][yModYSpacing]*sin(tempAngle) + yComponentUCArray[xModXSpacing][yModYSpacing]*cos(tempAngle))
####                tempSpinY = baseSigma*yComponentUCArray[xModXSpacing][yModYSpacing]
##                tempSpinZ = z*(1.0-tempSigma**2)**.5
##               
##                self.spinArray[x][y] = (tempSpinX,tempSpinY,tempSpinZ)
##
##
