WIP : Spin operators.

Author: amitjamadagni

src/arrays/arraymath.jl

@@ -61,15 +61,15 @@ end
 
 *(dm1::DualMatrix, dm2::DualMatrix) = (dm1.qarr*dm2.qarr)'
 
-function +{B<:OrthonormalBasis,N}(qarr1::Union(QuArray{B,N},CTranspose{B,N}), qarr2::Union(QuArray{B,N},CTranspose{B,N}))
+function +{B<:OrthonormalBasis,N}(qarr1::AbstractQuArray{B,N}, qarr2::AbstractQuArray{B,N})
     if bases(qarr1) == bases(qarr2)
         return QuArray(coeffs(qarr1)+coeffs(qarr2), bases(qarr1))
     else
         error("Bases not compatible")
   end
 end
 
-function -{B<:OrthonormalBasis,N}(qarr1::Union(QuArray{B,N},CTranspose{B,N}), qarr2::Union(QuArray{B,N},CTranspose{B,N}))
+function -{B<:OrthonormalBasis,N}(qarr1::AbstractQuArray{B,N}, qarr2::AbstractQuArray{B,N})
     if bases(qarr1) == bases(qarr2)
         return QuArray(coeffs(qarr1)-coeffs(qarr2), bases(qarr1))
     else

src/arrays/ladderops.jl

@@ -34,10 +34,48 @@ end
 
 function momentumop(n::Int)
     cop = raiseop(n)
-    return scale!(1/sqrt(2.), im*(cop-cop'))
+    return scale(im/sqrt(2.), cop-cop')
 end
 
+function spin_h1(j::Float64)
+    if (isinteger(2*j) && j>0)
+        m = reverse([-j:j-1])
+        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)
+    else
+        return error("j must be non-negative integer or half integer")
+  end
+end
+
+function spin_h2(j::Float64)
+    if (isinteger(2*j) && j>0)
+        m = reverse([-j:j])
+        N = length(m)
+        return spdiagm(m,0,N,N)
+    else
+        return error("j must be non-negative integer or half integer")
+  end
+end
+
+spin_jx(j::Float64) = QuArray(0.5*(spin_h1(j)+spin_h1(j)'))
+spin_jy(j::Float64) = QuArray(-0.5*im*(spin_h1(j)-spin_h1(j)'))
+spin_jz(j::Float64) = QuArray(spin_h2(j))
+spin_jp(j::Float64) = QuArray(spin_h1(j))
+spin_jm(j::Float64) = QuArray(spin_h1(j)')
+sigmax = 2.0 * spin_jx(0.5)
+sigmay = 2.0 * spin_jy(0.5)
+sigmaz = 2.0 * spin_jz(0.5)
+
 export raiseop,
     lowerop,
     positionop,
-    momentumop
+    momentumop,
+    spin_jx,
+    spin_jy,
+    spin_jz,
+    spin_jp,
+    spin_jm,
+    sigmax,
+    sigmay,
+    sigmaz

src/arrays/quarray.jl

@@ -35,7 +35,7 @@
 
     rawbases(qarr::QuArray, i) = qarr.bases[i]
     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 +105,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 #
 ######################