Merge cc584b8 into 79b4821

REQUIRE

@@ -4,5 +4,4 @@ OrdinaryDiffEq 3.19.1
 DiffEqCallbacks 1.1
 StochasticDiffEq 4.4.5
 RecursiveArrayTools
-GSL 0.3.6
 WignerSymbols 0.1.1

src/QuantumOptics.jl

@@ -24,7 +24,7 @@ export bases, Basis, GenericBasis, CompositeBasis, basis,
         nlevel, NLevelBasis, transition, nlevelstate,
         manybody, ManyBodyBasis, fermionstates, bosonstates,
                 manybodyoperator, onebodyexpect, occupation,
-        phasespace, qfunc, wigner, coherentspinstate, qfuncsu2, wignersu2,
+        phasespace, qfunc, wigner, coherentspinstate, qfuncsu2, wignersu2, ylm,
         metrics, tracenorm, tracenorm_h, tracenorm_nh,
                 tracedistance, tracedistance_h, tracedistance_nh,
                 entropy_vn, fidelity, ptranspose, PPT,

src/phasespace.jl

@@ -1,11 +1,10 @@
 module phasespace
 
-export qfunc, wigner, coherentspinstate, qfuncsu2, wignersu2
+export qfunc, wigner, coherentspinstate, qfuncsu2, wignersu2, ylm
 
 using ..bases, ..states, ..operators, ..operators_dense, ..fock, ..spin
 
 import WignerSymbols: clebschgordan
-import GSL: sf_legendre_sphPlm
 
 """
     qfunc(a, α)
@@ -323,7 +322,7 @@ function wignersu2(rho::DenseOperator, theta::Real, phi::Real)
         BandT[S+1,S+1] = clebschgordan(1,0,S,S,S+1,S)BandT[1,1][1:N+1-S].*BandT[S,S+1]+
             clebschgordan(1,1,S,S-1,S+1,S)*BandT[1,2][1:N+1-S].*BandT[S,S][2:end]
         BandT[S+1,S+2] = BandT[1,2][1:N+1-(S+1)].*BandT[S,S+1][2:end]
-        for  M=1:S-1
+        @inbounds for M=1:S-1
             BandT[S+1,M+1] = clebschgordan(1, 0, S, M, S+1,M)*BandT[1,1][1:N+1-M].*BandT[S,M+1] +
                 clebschgordan(1,1,S,M-1,S+1,M)*BandT[1,2][1:N+1-M].*BandT[S,M][2:end] -
                 clebschgordan(1,-1,S,M+1,S+1,M)*[zeros(1); BandT[1,2][1:N-M].*BandT[S,M+2][1:N-M]]
@@ -341,7 +340,7 @@ function wignersu2(rho::DenseOperator, theta::Real, phi::Real)
     end
 
     ### State decomposition ###
-    c = rho.data;
+    c = rho.data
     EVT = Array{Complex128}(N,N+1)
     @inbounds for S = 1:N, M = 0:S
         EVT[S,M+1] = conj(sum(BandT[S,M+1].*diag(c,M)))
@@ -373,7 +372,7 @@ function wignersu2(rho::DenseOperator, Ntheta::Int; Nphi::Int=2Ntheta)
         BandT[S+1,S+1] = clebschgordan(1,0,S,S,S+1,S)BandT[1,1][1:N+1-S].*BandT[S,S+1]+
             clebschgordan(1,1,S,S-1,S+1,S)*BandT[1,2][1:N+1-S].*BandT[S,S][2:end]
         BandT[S+1,S+2] = BandT[1,2][1:N+1-(S+1)].*BandT[S,S+1][2:end]
-        for  M=1:S-1
+        @inbounds for M=1:S-1
             BandT[S+1,M+1] = clebschgordan(1, 0, S, M, S+1,M)*BandT[1,1][1:N+1-M].*BandT[S,M+1] +
                 clebschgordan(1,1,S,M-1,S+1,M)*BandT[1,2][1:N+1-M].*BandT[S,M][2:end] -
                 clebschgordan(1,-1,S,M+1,S+1,M)*[0.0im; BandT[1,2][1:N-M].*BandT[S,M+2][1:N-M]]
@@ -389,15 +388,15 @@ function wignersu2(rho::DenseOperator, Ntheta::Int; Nphi::Int=2Ntheta)
     end
 
     ### State decomposition ###
-    c = rho.data;
+    c = rho.data
     EVT = Array{Complex128}(N,N+1)
     @inbounds for S = 1:N, M = 0:S
         EVT[S,M+1] = conj(sum(BandT[S,M+1].*diag(c,M)))
     end
 
     wignermap = Array{Float64}(Ntheta,Nphi)
     @inbounds for i = 1:Ntheta, j = 1:Nphi
-        wignermap[i,j] = _wignersu2int(N,i*1pi/(Ntheta-1),j*2pi/(Nphi-1)-pi, EVT)
+        wignermap[i,j] = _wignersu2int(N,i*1pi/(Ntheta-1)-1pi/(Ntheta-1),j*2pi/(Nphi-1)-2pi/(Nphi-1)-pi, EVT)
     end
     return wignermap*sqrt((N+1)/(4pi))
 end
@@ -414,8 +413,71 @@ function _wignersu2int(N::Integer, theta::Real, phi::Real, EVT::Array{Complex128
 end
 wignersu2(psi::Ket, args...) = wignersu2(dm(psi), args...)
 
-function ylm(l::Integer, m::Integer, theta::Real, phi::Real)
-    sf_legendre_sphPlm(l,m,cos(theta))*e^(1im*m*phi)
+"""
+    ylm(l::Integer,m::Integer,theta::Real,phi::Real)
+
+Spherical harmonics Y(l,m)(θ,ϕ) where l ∈ N,  m = -l,-l+1,...,l-1,l, θ ∈ [0,π],
+and ϕ ∈ [0,2π).
+
+This function calculates the value of Y(l,m) spherical harmonic at position θ and ϕ.
+"""
+function ylm(l::Integer,m::Integer,theta::Real,phi::Real)
+    phi_ = mod(phi,2pi)
+    theta_ = mod(theta,2pi)
+    phase = e^(1.0im*m*phi_)
+    if theta_ ≈ 0
+        if m == 0
+            return @. phase*sqrt((2*l+1)/pi)/2
+        else
+            return 0
+        end
+    elseif theta_ ≈ pi
+        if m == 0
+            return @. phase*(-1)^l*sqrt((2*l+1)/pi)/2
+        else
+            return 0
+        end
+    else
+        if l == 0
+            return 1.0/sqrt(4pi)
+        else
+            m_ = abs(m)
+            norm = _calc_ylm_norm(l, m_)
+            sign = m > 0 ? (-1)^m_ : 1
+            arg = cos(theta_)
+            p_ll = 1.0
+            @inbounds for fact = 1.0:l
+                p_ll *= @. 1.0/((2*fact))*sqrt(1-arg^2)
+            end
+
+            if m_ == l
+                return @. p_ll/norm*phase*sign
+            elseif l-m_ == 1
+                p_llp1 = @. 2*l*arg/sqrt(1-arg^2)*p_ll
+                return @. p_llp1/norm*phase*sign
+            else
+                p_llp1 = @. 2*l*arg/sqrt(1-arg^2)*p_ll
+                @inbounds for mr = -l:-m_-2
+                    p_llp2 = @. -2*(mr+1)*arg/sqrt(1-arg^2)*p_llp1-(l-mr)*(l+mr+1)*p_ll
+                    p_ll = p_llp1
+                    p_llp1 = p_llp2
+                end
+                return @. p_llp1/norm*phase*sign
+            end
+        end
+    end
+end
+
+function _calc_ylm_norm(l::Int, m_::Int)
+    # TODO: Clean up Int types
+    if 0 < l+m_ <= 33
+        norm = @. Float64(sqrt(4pi)/sqrt(2*l+1)*sqrt(factorial(Int128(l-m_)))/sqrt(factorial(Int128(l+m_))))
+    elseif 0 < l-m_ <= 33 && l+m_ > 33
+        norm = @. Float64(sqrt(4pi)/sqrt(2*l+1)*sqrt(factorial(Int128(l-m_)))/sqrt(factorial(BigInt(l+m_))))
+    else
+        norm = @. Float64(sqrt(4pi)/sqrt(2*l+1)*sqrt(factorial(BigInt(l-m_)))/sqrt(factorial(BigInt(l+m_))))
+    end
+    norm
 end
 
 end #module

test/test_phasespace.jl

@@ -124,4 +124,24 @@ qsu2dm = sum(qfuncsu2(dmrs,res).*costhetam)*(π/res)^2
 @test isapprox(wsu2, 1.0, atol=1e-2)
 @test isapprox(wsu2dm, 1.0, atol=1e-2)
 
+########### YLM test #############
+res = 1000
+l1 = rand(0:50)
+m1 = rand(0:l1)
+int = 0
+for i = 0:2pi/res:2pi, j = 0:pi/res:pi
+    int += sin(j)*abs2(ylm(l1,m1,j,i))
+end
+t1 = abs(int*2*(pi/res)^2)
+@test isapprox(t1, 1.00, atol=1e-2)
+
+l2 = rand(67:80)
+m2 = rand(0:30)
+int = 0
+for i = 0:2pi/res:2pi, j = 0:pi/res:pi
+    int += sin(j)*ylm(l1,m1,j,i)*conj(ylm(l2,m2,j,i))
+end
+t2 = abs(int*2*(pi/res)^2)
+@test isapprox(t2, 0, atol=1e-2)
+
 end # testset