Added : Operators

Position and momentum, Spin operators, Related Tests.

Conflicts:

	src/arrays/quarray.jl

src/arrays/ladderops.jl

@@ -37,45 +37,7 @@ function momentumop(n::Int)
  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,
- spin_jx,
- spin_jy,
- spin_jz,
- spin_jp,
- spin_jm,
- sigmax,
- sigmay,
- sigmaz
+ momentumop

src/arrays/quarray.jl

@@ -34,8 +34,6 @@
  coeffs(qarr::QuArray) = rawcoeffs(qarr)
 
  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)

src/arrays/spinops.jl

@@ -0,0 +1,77 @@
+#########################
+# Spin Operator Methods #
+#########################
+
+immutable HalfSpin{S}
+ function HalfSpin()
+ S > 0 ? new() : error("Spin must be positive")
+ end
+end
+
+spin(j::Integer) = j > 0 ? HalfSpin{2*j}() : error("Spin must be positive")
+
+function spin(j::Float64)
+ if (j > 0) && isinteger(j)
+ spin(int(j))
+ elseif (j > 0) && isinteger(2*j)
+ HalfSpin{int(2*j)}()
+ else
+ error("Spin must be positive and a multiple of 1/2.")
+ end
+end
+
+spin_value{S}(::HalfSpin{S}) = S/2
+
+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)