Added : Operators

Position and momentum, Spin operators, Related Tests.

Author: amitjamadagni

src/arrays/arraymath.jl

@@ -107,6 +107,15 @@ end
 
 tensor{B<:OrthonormalBasis}(bra::DualVector{B}, ket::QuVector{B}) = tensor(ket, bra)
 
+##############
+# Commutator #
+##############
+
+function commutator{B<:OrthonormalBasis}(qm1::AbstractQuMatrix{B}, qm2::AbstractQuMatrix{B})
+    return qm1*qm2-qm2*qm1
+end
+
 export normalize,
     normalize!,
-    tensor
+    tensor,
+    commutator

src/arrays/arrays.jl

@@ -5,4 +5,4 @@
     include("arraymath.jl")
     include("constructors.jl")
     include("ladderops.jl")
-    
+    include("spinops.jl")

src/arrays/ladderops.jl

@@ -1,11 +1,11 @@
 ###########################
 # Ladder Operator Methods #
-########################### 
-    # n specifies which particle (a.k.a tensor product 
+###########################
+    # n specifies which particle (a.k.a tensor product
     # factor) the operator acts on
     raiseop(b::AbstractBasis, n=1) = QuArray(raisematrix(size(b), n), b, b)
     raiseop(lens, n=1) = raiseop(FiniteBasis(lens), n)
-    
+
     lowerop(b::AbstractBasis, n=1) = QuArray(lowermatrix(size(b), n), b, b)
     lowerop(lens, n=1) = lowerop(FiniteBasis(lens), n)
 
@@ -22,10 +22,23 @@
     before_eye(lens, pivot) = speye(prod(lens[[1:pivot-1]]))
     after_eye(lens, pivot) = speye(prod(lens[[pivot+1:length(lens)]]))
     function eye_sandwich(lens, pivot, op)
-        return kron(before_eye(lens, pivot), 
-                    op, 
+        return kron(before_eye(lens, pivot),
+                    op,
                     after_eye(lens, pivot))
     end
 
+function positionop(n::Int)
+    cop = raiseop(n)
+    return scale!(1/sqrt(2.),cop+cop')
+end
+
+function momentumop(n::Int)
+    cop = raiseop(n)
+    return scale(im/sqrt(2.), cop-cop')
+end
+
 export raiseop,
-    lowerop
+    lowerop,
+    positionop,
+    momentumop
+

src/arrays/quarray.jl

@@ -35,7 +35,10 @@
 
     rawbases(qarr::QuArray, i) = qarr.bases[i]
     bases(qarr::QuArray, i) = rawbases(qarr, i)
-    
+
+    rawbases(qarr::AbstractQuArray) =  qarr.bases
+
+    bases(qarr::QuArray, i) = rawbases(qarr, i)
     # works generally as long as the single index form is defined
     rawbases(qarr::AbstractQuArray) =  ntuple(ndims(qarr), i->rawbases(qarr, i))
     bases(qarr::AbstractQuArray) = ntuple(ndims(qarr), i->bases(qarr, i))
@@ -105,11 +108,11 @@
 # LabelQuArray #
 ################
     typealias LabelQuArray{B<:LabelBasis,T,N,A} QuArray{B,T,N,A}
-    typealias TupleArray{T<:Tuple,N} Array{T,N} 
+    typealias TupleArray{T<:Tuple,N} Array{T,N}
 
     Base.getindex(larr::LabelQuArray, tups::Union(Tuple,TupleArray)...) = getindex(larr, map(getindex, bases(larr), tups)...)
     Base.setindex!(larr::LabelQuArray, x, tups::Union(Tuple,TupleArray)...) = setindex!(larr, x, map(getindex, bases(larr), tups)...)
-    
+
 ######################
 # Printing Functions #
 ######################

src/arrays/spinops.jl

@@ -0,0 +1,84 @@
+#########################
+# Spin Operator Methods #
+#########################
+
+immutable HalfSpin{S}
+    function HalfSpin()
+        S > 0 ? new() : error("Spin must be positive")
+    end
+end
+
+# The below function `spin` takes in a value and converts its to a type parameter
+# i.e., from value domain to type domain making it type unstable.
+
+spin(j::Integer) = j > 0 ? HalfSpin{j}() : error("Spin must be positive")
+
+function spin(j::Float64)
+    if (j > 0) && isinteger(j)
+        spin(int(j))
+    elseif (j > 0) && isinteger(2*j)
+        spin(int(2*j))
+    else
+        error("Spin must be positive and a multiple of 1/2.")
+    end
+end
+
+spin_value{S}(::HalfSpin{S}) = S/2
+
+#####################
+# Auxiliary methods #
+#####################
+
+function mat_coeffs_1(j::Float64)
+    m = [j-1:-1:-j]
+    N = length(m)+1
+    m_coeffs = [sqrt(j*(j+1.0) - (x+1.0)*x) for x in m]
+    return spdiagm(m_coeffs,1,N,N)
+end
+
+function mat_coeffs_2(j::Float64)
+    m = [j:-1:-j]
+    N = length(m)
+    return spdiagm(m,0,N,N)
+end
+
+function spin_h1(k::HalfSpin)
+    j = spin_value(k)
+    return mat_coeffs_1(j)
+end
+
+function spin_h2(k::HalfSpin)
+    j = spin_value(k)
+    return mat_coeffs_2(j)
+end
+
+################################
+# Jx, Jy, Jz, Jp, Jm operators #
+################################
+
+spin_Jx(j) = QuArray(0.5*(spin_h1(spin(j))+spin_h1(spin(j))'))
+spin_Jy(j) = QuArray(-0.5*im*(spin_h1(spin(j))-spin_h1(spin(j))'))
+spin_Jz(j) = QuArray(spin_h2(spin(j)))
+spin_Jp(j) = QuArray(spin_h1(spin(j)))
+spin_Jm(j) = QuArray(spin_h1(spin(j))')
+
+############################
+#  Pauli spin 1/2 matrices #
+############################
+
+const sigma_x = 2.0 * spin_Jx(1/2)
+const sigma_y = 2.0 * spin_Jy(1/2)
+const sigma_z = 2.0 * spin_Jz(1/2)
+
+# TODO : unicode sigma_y, sigma_z
+const σₓ = sigma_x
+
+export spin_Jx,
+    spin_Jy,
+    spin_Jz,
+    spin_Jp,
+    spin_Jm,
+    sigma_x,
+    sigma_y,
+    sigma_z,
+    σₓ

test/operatortest.jl

@@ -0,0 +1,17 @@
+#######################
+# Spin Operators Test #
+#######################
+
+@assert commutator(sigma_x, sigma_y) == 2*im*sigma_z
+@assert commutator(sigma_y, sigma_z) == 2*im*sigma_x
+@assert commutator(sigma_z, sigma_x) == 2*im*sigma_y
+@assert coeffs(commutator(sigma_x, sigma_x)) == spzeros(2,2)
+
+######################################
+# Position & Momentum Operators Test #
+######################################
+
+p = positionop(2)
+m = momentumop(2)
+@assert coeffs(commutator(sigma_x, p)) == spzeros(2,2)
+@assert coeffs(commutator(sigma_y, m)) == spzeros(2,2)

test/runtests.jl

@@ -4,3 +4,4 @@ using Base.Test
 include("multest.jl")
 include("tensortest.jl")
 include("labelbasistest.jl")
+include("operatortest.jl")