Skip to content

Commit

Permalink
Add spherical harmonics (qojulia#226)
Browse files Browse the repository at this point in the history 
* add spherical harmonics

* add spherical harmonics

* Some style and effiency changes

* add extension to negative m

* Some more style/efficiency

* Enhancement for small l

* Correct checking for large numbers

* fix spherical harmonics and wignersu2

* fix domain of ylm

* add description and remove old GSL ylm

* add tests

* Final style changes

* Update tests

* Update tests again

* tune tests

* Tune tests again
  • Loading branch information
@karolpezet @david-pl
karolpezet authored and david-pl committed Jul 24, 2018
1 parent 79b4821 commit e2a407f
Show file tree
Hide file tree
Showing 4 changed files with 94 additions and 11 deletions.
1 change: 0 additions & 1 deletion REQUIRE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,5 +4,4 @@ OrdinaryDiffEq 3.19.1
DiffEqCallbacks 1.1
StochasticDiffEq 4.4.5
RecursiveArrayTools
GSL 0.3.6
WignerSymbols 0.1.1
2 changes: 1 addition & 1 deletion src/QuantumOptics.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
80 changes: 71 additions & 9 deletions src/phasespace.jl
View file
Original file line number Diff line number Diff line change
@@ -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, α)
Expand Down Expand Up @@ -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]]
Expand All @@ -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)))
Expand Down Expand Up @@ -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]]
Expand All @@ -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
Expand All @@ -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
22 changes: 22 additions & 0 deletions test/test_phasespace.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,4 +124,26 @@ 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
int = 0
l = rand(4:33)
m = rand(0:l-1)
for i = 0:2pi/res:2pi, j = 0:pi/res:pi
int += sin(j)*abs2(ylm(l,m,j,i))
end
t1 = abs(int*2*(pi/res)^2)
@test isapprox(t1, 1.00, atol=1e-2)

l1 = rand(33:40)
m1 = rand(8:30)
l2 = rand(77:80)
m2 = rand(0:10)
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

0 comments on commit e2a407f

Please sign in to comment.