from numpy import *
from scipy import *
from visual import *
import time


class waveFunction:
    def __init__(self,n,l,m,elec):
        self.n = n
        self.l = l
        self.m = m
        self.elec = elec
        self.tonorm = 1.0
        self.vSign = 1.0

        self.c = 300000000
        
        self.tolerance = 1e-6
        self.maxIterations = 50
        self.maxSize = 10.0*self.n**2
        
        self.gamma0 = array([[1.0,0.0,0.0,0.0],[0.0,1.0,0.0,0.0],[0.0,0.0,-1.0,0.0],[0.0,0.0,0.0,-1.0]])
##        self.gamma0 = array([[0.0,0.0,1.0,0.0],[0.0,0.0,0.0,1.0],[1.0,0.0,0.0,0.0],[0.0,1.0,0.0,0.0]])
        self.gamma1 = array([[0.0,0.0,0.0,1.0],[0.0,0.0,1.0,0.0],[0.0,-1.0,0.0,0.0],[-1.0,0.0,0.0,0.0]])
        self.gamma2 = array([[0.0,0.0,0.0,-1.0j],[0.0,0.0,1.0j,0.0],[0.0,1.0j,0.0,0.0],[-1.0j,0.0,0.0,0.0]])
        self.gamma3 = array([[0.0,0.0,1.0,0.0],[0.0,0.0,0.0,-1.0],[-1.0,0.0,0.0,0.0],[0.0,1.0,0.0,0.0]])

    def radialWaveFunction(self,r):
        return special.genlaguerre(self.n-self.l-1,2*self.l+1)(2.0*r/(self.n))*e**(-r/(self.n))*(2.0*r/(self.n))**self.l

    def thetaphiWaveFunction(self,theta,phi):
        return special.sph_harm(self.m,self.l,theta,phi)

    def returnWaveFunction(self,r,theta,phi):
        return self.tonorm*self.radialWaveFunction(r)*self.thetaphiWaveFunction(theta,phi)

    def integrateThetaPhiSeparately(self):
        return integrate.dblquad(lambda phi, theta: (conjugate(self.thetaphiWaveFunction(theta,phi))*self.thetaphiWaveFunction(theta,phi)).real*sin(phi), 0, 2*pi, lambda phi: 0, lambda phi: pi, epsrel=self.tolerance)[0]

    def integrateRadialSeparately(self):
        return integrate.quad(lambda r: conjugate(self.radialWaveFunction(r))*self.radialWaveFunction(r)*r**2, 0, Inf, epsrel=self.tolerance, limit=self.maxIterations)[0]
   
    def normalize(self):
        norm = self.integrateRadialSeparately()*self.integrateThetaPhiSeparately()
        self.tonorm = 1/sqrt(norm)

    def returnSpinor(self,r,theta,phi):
        waveFunction = self.returnWaveFunction(r,theta,phi)

        if self.elec == -1:
            elec1coeff = 0
            elec2coeff = 1
        if self.elec == 1:
            elec1coeff = 1
            elec2coeff = 0
        if self.elec == 0:
            elec1coeff = 1
            elec2coeff = 1

        electrons = array([elec1coeff,elec2coeff])
        spins = array([[0.0,-sin(theta)-1j*cos(theta)],[-sin(theta)+1j*cos(theta),0.0]])

        coeffs1 = array(matrixmultiply(spins,electrons))
        coeffs2 = array(matrixmultiply(conjugate(spins),electrons))

        if elec1coeff == 1 and elec2coeff == 0:
            return 1/sqrt(2.0)*array([coeffs1[1]*waveFunction,coeffs2[1]*waveFunction,0.0,0.0])
        if elec2coeff == 1 and elec1coeff == 0:
            return 1/sqrt(2.0)*array([coeffs1[0]*waveFunction,coeffs2[0]*waveFunction,0.0,0.0])
        if elec1coeff == 1 and elec2coeff == 1:
            return 1/sqrt(4.0)*array([coeffs1[1]*waveFunction,coeffs2[1]*waveFunction,coeffs2[0]*waveFunction,coeffs1[0]*waveFunction])

    def getProbCurrent(self,r,phi):
        derivTheta = self.tonorm*derivative(self.thetaphiWaveFunction,0.0,dx=2*pi/self.maxIterations,n=1,order=5,args=(phi,))
        part1theta = conjugate(self.returnWaveFunction(r,0.0,phi))*derivTheta*self.radialWaveFunction(r)
        part2theta = self.returnWaveFunction(r,0.0,phi)*conjugate(derivTheta)*self.radialWaveFunction(r)
        thetaCurrent = (1/(2j)*(part1theta-part2theta)).real
        return vector(0.0,thetaCurrent,0.0)

    def getDensity(self,r,phi):
        return self.tonorm**2*(dot(conjugate(self.returnSpinor(r,0.0,phi)),self.returnSpinor(r,0.0,phi))).real

    def getProbVel(self,r,phi):
        return self.getProbCurrent(r,phi)/self.getDensity(r,phi)

    def getSpeed(self,r,phi):
        vel = mag(self.getProbVel(r,phi))
        return vel*sqrt(1-vel**2/self.c**2)

    def getvSign(self):
        self.vSign = sign(self.getProbVel(1.0,pi/2)[1])
        return self.vSign

    def integratePhiForMeanVel(self,r):
        return integrate.quad(lambda phi: self.getDensity(r,phi)*self.getSpeed(r,phi)*r**2*sin(phi), 0, pi, limit=self.maxIterations, epsrel=1)[0]

    def integrateRadialForMeanVel(self):
        return 2*pi*integrate.quad(lambda r: self.integratePhiForMeanVel(r), 0, Inf, limit=self.maxIterations, epsrel=1)[0]

    def getMeanVelocity(self):
        return self.integrateRadialForMeanVel()
   
    def getVelocityNorm(self,theta):
        return self.vSign*vector(-sin(theta),cos(theta),0.0)

    def getVelocity(self,r,theta,phi):
        return self.getSpeed(r,phi)*self.getVelocityNorm(theta)

    def getTransformation(self,r,theta,phi):
        speedOverc = self.getSpeed(r,phi)/self.c
        nVector = norm(self.getVelocity(r,theta,phi))
        
        phi = arctanh(speedOverc)
        coshExp = cosh(phi/2)
        sinhCon = sinh(phi/2)

        return array([[coshExp,0.0,-nVector[2]*sinhCon,-(nVector[0]-1j*nVector[1])*sinhCon],[0.0,coshExp,-(nVector[0]+1j*nVector[1])*sinhCon,nVector[2]*sinhCon],[-nVector[2]*sinhCon,-(nVector[0]-1j*nVector[1])*sinhCon,coshExp,0.0],[-(nVector[0]+1j*nVector[1])*sinhCon,nVector[2]*sinhCon,0.0,coshExp]])

    def getTransformedSpinor(self,r,theta,phi):
        return array(matrixmultiply(self.getTransformation(r,theta,phi),self.returnSpinor(r,theta,phi)))

    def getTransformedSpinorRadialFirst(self,r,theta,phi):
        return self.getTransformedSpinor(r,theta,phi)

    def getTransformedSpinorThetaFirst(self,theta,r,phi):
        return self.getTransformedSpinor(r,theta,phi)

    def getTransformedSpinorPhiFirst(self,phi,r,theta):
        return self.getTransformedSpinor(r,theta,phi)

    def returnDifferentialSpinor(self,r,theta,phi):
        spinorRadial = self.getTransformedSpinorRadialFirst
        spinorTheta = self.getTransformedSpinorThetaFirst
        spinorPhi = self.getTransformedSpinorPhiFirst

        radialDifferential = derivative(spinorRadial,r,dx=self.maxSize/self.maxIterations,n=1,order=5,args=(theta,phi))*r*sin(phi)
        thetaDifferential = derivative(spinorTheta,theta,dx=2*pi/self.maxIterations,n=1,order=5,args=(r,phi))
        phiDifferential = derivative(spinorPhi,phi,dx=pi/self.maxIterations,n=1,order=5,args=(r,theta))*sin(phi)

        self.gammaR = self.gamma1*sin(phi)*cos(theta) + self.gamma2*sin(phi)*sin(theta) + self.gamma3*cos(phi)
        self.gammaT = - self.gamma1*r*sin(phi)*sin(theta) + self.gamma2*r*sin(phi)*cos(theta)
        self.gammaP = self.gamma1*r*cos(phi)*cos(theta) + self.gamma2*r*cos(phi)*sin(theta) - self.gamma3*r*sin(phi)

        differentialSpinor = array(matrixmultiply(self.gammaR,radialDifferential)+matrixmultiply(self.gammaT,thetaDifferential)+matrixmultiply(self.gammaP,phiDifferential))

##        print "differential spinor", matrixmultiply(self.gamma0,differentialSpinor)
##        print "conjugate spinor", conjugate(self.getTransformedSpinor(r,theta,phi))
##        print "product", (dot(conjugate(self.getTransformedSpinor(r,theta,phi)),matrixmultiply(self.gamma0,differentialSpinor))).imag

        return differentialSpinor






class integrations:
    def __init__(self,spinorOne,spinorTwo):
        self.spinorOne = spinorOne
        self.spinorTwo = spinorTwo

        self.maxIterations = 50
        self.tolerance = 5e-6
        self.maxSize = 5.0*max(self.spinorOne.n**2,self.spinorTwo.n**2)
        self.c = 300000000
        
        self.gamma0 = array([[1.0,0.0,0.0,0.0],[0.0,1.0,0.0,0.0],[0.0,0.0,-1.0,0.0],[0.0,0.0,0.0,-1.0]])
##        self.gamma0 = array([[0.0,0.0,1.0,0.0],[0.0,0.0,0.0,1.0],[1.0,0.0,0.0,0.0],[0.0,1.0,0.0,0.0]])

##    def integrateOverlapTheta(self,r):
##        return integrate.quad(lambda theta: (dot(conjugate(self.spinorOne.getTransformedSpinor(r,theta)),array(matrixmultiply(self.gamma0,self.spinorTwo.getTransformedSpinor(r,theta))))).real, 0, 2*pi, epsrel=self.tolerance, limit=self.maxIterations)[0]
##
##    def integrateOverlapRadial(self):
##        return integrate.quad(lambda r: self.integrateOverlapTheta(r)*r, 0, self.maxSize, epsrel=self.tolerance, limit=self.maxIterations)[0]
##
##    def integrateOverlap(self):
##        return self.integrateOverlapRadial()
##
##
##    def integrateNuclearTheta(self,r):
##        return integrate.quad(lambda theta: (dot(conjugate(self.spinorOne.getTransformedSpinor(r,theta)),array(matrixmultiply(self.gamma0,self.spinorTwo.getTransformedSpinor(r,theta))))).real, 0, 2*pi, epsrel=self.tolerance, limit=self.maxIterations)[0]
##
##    def integrateNuclearRadial(self):
##        return integrate.quad(lambda r: self.integrateNuclearTheta(r), 0, self.maxSize, epsrel=self.tolerance, limit=self.maxIterations)[0]
##
##    def integrateNuclear(self):
##        return self.integrateNuclearRadial()

    def integrateKineticPhi(self,r,theta):
        return integrate.quadrature(lambda phi: (1j*dot(conjugate(self.spinorOne.getTransformedSpinor(r,theta,phi)),matrixmultiply(self.gamma0,self.spinorTwo.returnDifferentialSpinor(r,theta,phi)))).real*r, 0, pi, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]
 
    def integrateKineticTheta(self,r):
        return integrate.quadrature(lambda theta: self.integrateKineticPhi(r,theta), 0, 2*pi, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]
       
    def integrateKineticRadial(self):
        return integrate.quadrature(lambda r: self.integrateKineticTheta(r), 0, self.maxSize, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]
        
    def integrateKinetic(self):
        return self.integrateKineticRadial()


##    def getRelativisticDistance(self,r1,r2,theta1,theta2):
##        velElectronOne = self.spinorOne.getVelocity(r1,theta1)
##        positionElectronTwoRelativeToElectronOneInLattice = vector(r1*cos(theta1),r1*sin(theta1),0.0)-vector(r2*cos(theta2),r2*sin(theta2),0.0)
##        gamma = 1/sqrt(1-dot(velElectronOne,velElectronOne)/self.c**2)
##        
##        if mag(velElectronOne) == 0.0:
##            positionElectronTwoRelativeToElectronOne = positionElectronTwoRelativeToElectronOneInLattice
##        else:
##            positionElectronTwoRelativeToElectronOne = positionElectronTwoRelativeToElectronOneInLattice+(gamma-1)/dot(velElectronOne,velElectronOne)*dot(velElectronOne,positionElectronTwoRelativeToElectronOneInLattice)*velElectronOne
##        dist = mag(positionElectronTwoRelativeToElectronOne)
##        if dist == 0.0:
##            return self.tolerance
##        else:
##            return dist
##
##
##    def integrandForExchange(self,r1,r2,theta1,theta2):
##        return (dot(matrixmultiply(conjugate(self.spinorOne.getTransformedSpinor(r1,theta1)),self.gamma0),self.spinorOne.getTransformedSpinor(r2,theta2))*r1*r2/self.getRelativisticDistance(r1,r2,theta1,theta2)*dot(matrixmultiply(conjugate(self.spinorTwo.getTransformedSpinor(r2,theta2)),self.gamma0),self.spinorTwo.getTransformedSpinor(r1,theta1))).real
##
##    def integrateExchangeTheta2(self,r1,r2,theta1):
##        return integrate.quadrature(lambda theta2: self.integrandForExchange(r1,r2,theta1,theta2), 0, 2*pi, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]   
##
##    def integrateExchangeTheta1(self,r1,r2):
##        return integrate.quadrature(lambda theta1: self.integrateExchangeTheta2(r1,r2,theta1), 0, 2*pi, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]
##
##    def integrateExchangeRadial2(self,r1):
##        return integrate.quadrature(lambda r2: self.integrateExchangeTheta1(r1,r2), 0, self.maxSize, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]
##    
##    def integrateExchangeRadial1(self):
##        return integrate.quadrature(lambda r1: self.integrateExchangeRadial2(r1), 0, self.maxSize, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]
##
##    def integrateExchange(self):
##        return self.integrateExchangeRadial1()
##
##
##    def integrandForCoulombic(self,r1,r2,theta1,theta2):
##        return (dot(matrixmultiply(conjugate(self.spinorOne.getTransformedSpinor(r1,theta1)),self.gamma0),self.spinorOne.getTransformedSpinor(r1,theta1))*r1*r2/self.getRelativisticDistance(r1,r2,theta1,theta2)*dot(matrixmultiply(conjugate(self.spinorTwo.getTransformedSpinor(r2,theta2)),self.gamma0),self.spinorTwo.getTransformedSpinor(r2,theta2))).real
##
##    def integrateCoulombicTheta2(self,r1,r2,theta1):
##        return integrate.quadrature(lambda theta2: self.integrandForCoulombic(r1,r2,theta1,theta2), 0, 2*pi, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]   
##
##    def integrateCoulombicTheta1(self,r1,r2):
##        return integrate.quadrature(lambda theta1: self.integrateCoulombicTheta2(r1,r2,theta1), 0, 2*pi, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]   
##
##    def integrateCoulombicRadial2(self,r1):
##        return integrate.quadrature(lambda r2: self.integrateCoulombicTheta1(r1,r2), 0, self.maxSize, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]
##    
##    def integrateCoulombicRadial1(self):
##        return integrate.quadrature(lambda r1: self.integrateCoulombicRadial2(r1), 0, self.maxSize, tol=self.tolerance, maxiter=self.maxIterations, vec_func=False)[0]
##
##    def integrateCoulombic(self):
##        return self.integrateCoulombicRadial1()




class basis:
    def __init__(self,maxn):
        self.maxn = maxn
        self.tonorm = 1

        count = 0
        for n in xrange(maxn):
            for l in xrange(n+1):
                for m in xrange(2*l+1):
                    count += 1
        self.maxNbasis = count
        
        print "Total modes:", self.maxNbasis
        
        i = 0
        self.basis = [0]*self.maxNbasis
        for n in xrange(maxn):
            for l in xrange(n+1):
                for m in xrange(2*l+1):
                    if m-l<0:
                        elec = -1
                    if m-l>0:
                        elec = 1
                    if m-l==0:
                        elec = 0
                        
                    print "Radial QN", n+1
                    print "Angular QN", l
                    print "Magnetic QN", m-l
                    print "Elec index", elec
                    
                    self.basis[i] = waveFunction(n+1,l,m-l,elec)
                    self.basis[i].normalize()
                    
                    print "Normalization was:", self.basis[i].tonorm
##                    print "Mean velocity was:", self.basis[i].getvSign()*self.basis[i].getMeanVelocity()
                    i += 1

        self.integrations = [0]*self.maxNbasis
        for i in xrange(self.maxNbasis):
            self.integrations[i] = [0]*self.maxNbasis
            for j in xrange(self.maxNbasis):       
                self.integrations[i][j] = integrations(self.basis[i],self.basis[j])

        self.overlaps = ones((self.maxNbasis,self.maxNbasis),Float)
        self.Noverlaps = ones((self.maxNbasis,self.maxNbasis),Float)
        self.KEoverlaps = ones((self.maxNbasis,self.maxNbasis),Float)
        self.XCoverlaps = ones((self.maxNbasis,self.maxNbasis),Float)
        self.COoverlaps = ones((self.maxNbasis,self.maxNbasis),Float)

        print "Basis and Integrations Set Up and Normalized at:", time.ctime()

    def fillOverlaps(self):
        for i in xrange(self.maxNbasis):
            for j in xrange(i+1):
                overlap = self.integrations[i][j].integrateOverlap()
                if abs(overlap) < 1e-15:
                    overlap = 0.0
                self.overlaps[i][j] = overlap
                self.overlaps[j][i] = overlap
                print "Overlap",i,j,"is:",overlap,"at:",time.ctime()

    def fillNoverlaps(self):
        for i in xrange(self.maxNbasis):
            for j in xrange(i+1):
                Noverlap = self.integrations[i][j].integrateNuclear()
                if abs(Noverlap) < 1e-15:
                    Noverlap = 0.0
                self.Noverlaps[i][j] = Noverlap
                self.Noverlaps[j][i] = Noverlap
                print "Nuclear",i,j,"is:",Noverlap,"at:",time.ctime()

    def fillKEoverlaps(self):
        for i in xrange(self.maxNbasis):
            for j in xrange(i+1):
                KEoverlap = self.integrations[i][j].integrateKinetic()
                if abs(KEoverlap) < 1e-15:
                    KEoverlap = 0.0
                self.KEoverlaps[i][j] = KEoverlap
                self.KEoverlaps[j][i] = KEoverlap
                print "Kinetic",i,j,"is:",KEoverlap,"at:",time.ctime()

    def fillXCEnergies(self):
        for i in xrange(self.maxNbasis):
            for j in xrange(i+1):
                XCoverlap = self.integrations[i][j].integrateExchange()
                if abs(XCoverlap) < 1e-15:
                    XCoverlap = 0.0
                self.XCoverlaps[i][j] = XCoverlap
                self.XCoverlaps[j][i] = XCoverlap
                print "Exchange Overlap",i,j,"is:",XCoverlap,"at:",time.ctime()

    def fillCOEnergies(self):
        for i in xrange(self.maxNbasis):
            for j in xrange(i+1):
                COoverlap = self.integrations[i][j].integrateCoulombic()
                if abs(COoverlap) < 1e-15:
                    COoverlap = 0.0
                self.COoverlaps[i][j] = COoverlap
                self.COoverlaps[j][i] = COoverlap
                print "Coulombic Overlap",i,j,"is:",COoverlap,"at:",time.ctime()

    def short(self):
##        print "0 0 overlap is:", self.integrations[0][0].integrateOverlap()
##        print "0 0 nuclear overlap is:",self.integrations[0][0].integrateNuclear()
        print "0 0 kinetic overlap is:",self.integrations[0][0].integrateKinetic()
##        print "0 0 exchange overlap is:",self.integrations[0][0].integrateExchange()
##        print "0 0 coulombic overlap is:",self.integrations[0][0].integrateCoulombic()

        print "1 1 kinetic overlap is:",self.integrations[1][1].integrateKinetic()
        print "2 2 kinetic overlap is:",self.integrations[2][2].integrateKinetic()

##        print "3 3 overlap is:", self.integrations[3][3].integrateOverlap()
##        print "3 3 nuclear overlap is:",self.integrations[3][3].integrateNuclear()
        print "3 3 kinetic overlap is:",self.integrations[3][3].integrateKinetic()
        print "4 4 kinetic overlap is:",self.integrations[4][4].integrateKinetic()
##        print "3 3 exchange overlap is:",self.integrations[3][3].integrateExchange()
##        print "3 3 coulombic overlap is:",self.integrations[3][3].integrateCoulombic()

        

maxn = 3

Basis = basis(maxn)

##Basis.fillOverlaps()
##dump = open("dumpOverlaps", mode='w')
##dump.write(str(Basis.overlaps))
##dump.close()
##
##Basis.fillNoverlaps()
##dump = open("dumpNoverlaps", mode='w')
##dump.write(str(Basis.Noverlaps))
##dump.close()

##Basis.fillKEoverlaps()
##dump = open("dumpKEoverlaps", mode='w')
##dump.write(str(Basis.KEoverlaps))
##dump.close()

##Basis.fillXCEnergies()
##dump = open("dumpXCoverlaps", mode='w')
##dump.write(str(Basis.XCoverlaps))
##dump.close()
##
##Basis.fillCOEnergies()
##dump = open("dumpCOoverlaps", mode='w')
##dump.write(str(Basis.COoverlaps))
##dump.close()

Basis.short()









##class densityPoint:
##    def __init__(self):
##        self.visibility = 1
##        self.point = sphere(visible=self.visibility, radius=.001)
##
##    def setAttributes(self,position,size,color):
##        self.point.pos = position
##        self.point.radius = size
##        self.point.color = color
##
##    def toggleVisibility(self):
##        self.visibility = (1+self.visibility)%2
##
##class currentVector:
##    def __init__(self):
##        self.visibility = 1
##        self.currentArrow = arrow(visible=self.visibility)
##
##    def setAttributes(self,position,direction):
##        self.currentArrow.pos = position
##        self.currentArrow.axis = direction
##
##    def toggleVisibility(self):
##        self.visibility = (1+self.visibility)%2
##
##class visualLattice:
##    def __init__(self,electronOne,sizeOfLattice,center):
##        self.electronOne = electronOne
##        self.sizeOfLattice = sizeOfLattice
##        self.center = center
##        self.scale = sizeOfLattice
##
##        self.pointLattice=[None]*sizeOfLattice
##        self.currentLattice=[None]*sizeOfLattice
##        for x in xrange(sizeOfLattice):
##            self.pointLattice[x]=[None]*sizeOfLattice
##            self.currentLattice[x]=[None]*sizeOfLattice
##            for y in xrange(sizeOfLattice):
##                position = vector(x-sizeOfLattice/2.0,y-sizeOfLattice/2.0,0.0)
##                r = mag(position)
##                theta = atan2(position[1],position[0])
##                currentDirection = self.electronOne.getVelocityRelLattice(r,theta)
##                transformationMatrix = self.electronOne.getTransformation(r,theta)
##                
##                density = self.electronOne.getDensity(r,theta)
##                size = self.scale*density
##                phase = self.electronOne.getPhase(r,theta)
##                colorArg = phase/pi
##                color = ((1+abs(colorArg))/2.0,0.0,(1-abs(colorArg))/2.0)
##                
##                self.pointLattice[x][y] = densityPoint()
##                self.currentLattice[x][y] = currentVector()
##                self.pointLattice[x][y].setAttributes(position,size,color)
##                self.currentLattice[x][y].setAttributes(position,currentDirection)
##
##    def updateLattice(self,electronOne):
##        for x in xrange(self.sizeOfLattice):
##            for y in xrange(self.sizeOfLattice):
##                position = vector(x-sizeOfLattice/2.0,y-sizeOfLattice/2.0,0.0)
##                r = mag(position)
##                theta = atan2(position[1],position[0])
##                currentDirection = self.electronOne.getVelocityRelLattice(r,theta)
##                transformationMatrix = self.electronOne.getTransformation(r,theta)
##                
##                density = self.electronOne.getDensity(r,theta)
##                size = self.scale*density
##                phase = self.electronOne.getPhase(r,theta)
##                colorArg = phase/pi
##                color = ((1+abs(colorArg))/2.0,0.0,(1-abs(colorArg))/2.0)
##                
##                self.pointLattice[x][y].setAttributes(position,size,color)
##                self.currentLattice[x][y].setAttributes(position,currentDirection)


##count = 0
##for n in xrange(maxn):
##    for l in xrange(n+1):
##        for m in xrange(2*l+1):
##            if l-1 == abs(int(m-l)):
##                pass
##            else:
##                count += 1
##
##coeffsOne = [None]*count
##coeffsTwo = [None]*count
##i = 0
##for n in xrange(maxn):
##    for l in xrange(n+1):
##        for m in xrange(2*l+1):
##            if l-1 == abs(int(m-l)):
##                pass
##            else:
##                if m-l==-2:
##                    coeffsOne[i] = 1.0
##                else:
##                    coeffsOne[i] = 0.0
##                if m-l==2:
##                    coeffsTwo[i] = 1.0
##                else:
##                    coeffsTwo[i] = 0.0
##                i += 1
##
##print "Electron One Coefficients", coeffsOne
##print "Electron Two Coefficients", coeffsTwo


##sizeOfLattice = 24
##center = (0.0,0.0,0.0)
##
##scene = display(title='Density Plot', x=0, y=0, width=600, height=600, center=center, background=(0,0,0))
##visualLattice = visualLattice(electronOne, sizeOfLattice, center)
##visualLattice.updateLattice(electronOne)

##results = list()
##
##for v in xrange(0,21,1):
##    velScale = float(v/3000.0)
##
##    integrations.normalize()
##    
##    meanVel = integrations.getMeanVelocity()
##    xcInt = integrations.integrateExchange()
##    coInt = integrations.integrateCoulombic()
##
##    print "Mean Velocity is:", meanVel
##    print "Exchange Integral is:", xcInt
##    print "Coulombic Integral is:", coInt
##
##    results.append([meanVel,xcInt,coInt])
##
##    print "Time was:", time.ctime()
##
##dump = open("dumpAnalysis", mode='w')
##dump.write(str(results))
##dump.close()


##    def getPhase(self,r,theta):
##        return atan2(self.returnSpinor(r,theta)[0].imag,self.returnSpinor(r,theta)[0].real)


