mirror of
https://github.com/tgorordo/WignerSymbols.jl.git
synced 2026-06-05 15:42:15 -07:00
3j and 6j working and tested
This commit is contained in:
Jutho
parent
b9eb46d5cd
commit
56b02ebcf2
3 changed files with 137 additions and 20 deletions
|
|
@ -10,10 +10,10 @@ clebschgordan(j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) = clebschgordan(Fl
|
|||
function clebschgordan(T::Type{<:AbstractFloat}, j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂)
|
||||
s = wigner3j(T, j₁, j₂, j₃, m₁, m₂, -m₃)
|
||||
iszero(s) && return s
|
||||
if isodd(j₁ - j₂ + m₃)
|
||||
if isodd(convert(Int,j₁ - j₂ + m₃))
|
||||
s *= -one(s)
|
||||
end
|
||||
s /= sqrt(convert(T, j₃+j₃+1))
|
||||
s *= sqrt(convert(T, j₃+j₃+one(j₃)))
|
||||
return s
|
||||
end
|
||||
|
||||
|
|
@ -34,8 +34,8 @@ function wigner3j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, m₃ =
|
|||
β₂ = convert(Int, j₁ - m₁ )
|
||||
β₃ = convert(Int, j₂ + m₂ )
|
||||
|
||||
# extra sign in definition
|
||||
sgn = isodd(j₁ - j₃ - m₃) ? -sgn : sgn
|
||||
# extra sign in definition: α₁ - α₂ = j₁ + m₂ - j₂ + m₁ = j₁ - j₂ + m₃
|
||||
sgn = isodd(α₁ - α₂) ? -sgn : sgn
|
||||
|
||||
# dictionary lookup or compute
|
||||
if haskey(Wigner3j, (β₁, β₂, β₃, α₁, α₂))
|
||||
|
|
@ -44,7 +44,8 @@ function wigner3j(T::Type{<:AbstractFloat}, j₁, j₂, j₃, m₁, m₂, m₃ =
|
|||
s1 = Δ²(j₁, j₂, j₃)
|
||||
s2 = prod(map(primefactorial, (β₂, β₁ - α₁, β₁ - α₂, β₃, β₃ - α₁, β₂ - α₂)))
|
||||
|
||||
snum, rnum = splitsquare(s1.num * s2)
|
||||
|
||||
snum, rnum = splitsquare(s1.num*s2)
|
||||
sden, rden = splitsquare(s1.den)
|
||||
s = convert(BigInt, snum) // convert(BigInt, sden)
|
||||
r = convert(BigInt, rnum) // convert(BigInt, rden)
|
||||
|
|
@ -118,6 +119,7 @@ function δ(j₁, j₂, j₃)
|
|||
end
|
||||
|
||||
# triangle coefficient
|
||||
Δ(j₁, j₂, j₃) = Δ(Float64, j₁, j₂, j₃)
|
||||
function Δ(T::Type{<:AbstractFloat}, j₁, j₂, j₃)
|
||||
if !δ(j₁, j₂, j₃)
|
||||
throw(DomainError())
|
||||
|
|
@ -139,7 +141,7 @@ end
|
|||
|
||||
# reorder parameters determining the 3j symbol to canonical order:
|
||||
# j₁ >= j₂ >= j₃ and m₁ >= 0 or m₁ == 0 && m₂ >= 0
|
||||
function reorder3j(j₁, j₂, j₃, m₁, m₂, m₃, sign = one(UInt8))
|
||||
function reorder3j(j₁, j₂, j₃, m₁, m₂, m₃, sign = one(Int8))
|
||||
if j₁ < j₂
|
||||
return reorder3j(j₂, j₁, j₃, m₂, m₁, m₃, -sign)
|
||||
elseif j₂ < j₃
|
||||
|
|
@ -150,7 +152,7 @@ function reorder3j(j₁, j₂, j₃, m₁, m₂, m₃, sign = one(UInt8))
|
|||
return reorder3j(j₁, j₂, j₃, -m₁, -m₂, -m₃, -sign)
|
||||
else
|
||||
# sign doesn't matter if total J=j₁ + j₂ + j₃ is even
|
||||
if iseven(j₁ + j₂ + j₃)
|
||||
if iseven(convert(UInt,j₁ + j₂ + j₃))
|
||||
sign = one(sign)
|
||||
end
|
||||
return (j₁, j₂, j₃, m₁, m₂, m₃, sign)
|
||||
|
|
@ -183,14 +185,14 @@ function compute3jseries(β₁, β₂, β₃, α₁, α₂)
|
|||
nums = Vector{T}(length(krange))
|
||||
dens = Vector{T}(length(krange))
|
||||
for (i, k) in enumerate(krange)
|
||||
num = iseven(k) ? primefactorial(1) : -primefactorial(1)
|
||||
num = iseven(k) ? one(T) : -one(T)
|
||||
den = primefactorial(k)*primefactorial(k-α₁)*primefactorial(k-α₂)*
|
||||
primefactorial(β₁-k)*primefactorial(β₂-k)*primefactorial(β₃-k)
|
||||
nums[i], dens[i] = divgcd!(num, den)
|
||||
end
|
||||
den = commondenominator!(nums, dens)
|
||||
totalnum = sumlist!(nums)
|
||||
totalden = convert(BigInt, PrimeFactorization(den))
|
||||
totalden = convert(BigInt, den)
|
||||
return totalnum//totalden
|
||||
end
|
||||
|
||||
|
|
@ -209,7 +211,7 @@ function compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄)
|
|||
end
|
||||
den = commondenominator!(nums, dens)
|
||||
totalnum = sumlist!(nums)
|
||||
totalden = convert(BigInt, PrimeFactorization(den))
|
||||
totalden = convert(BigInt, den)
|
||||
return totalnum//totalden
|
||||
end
|
||||
|
||||
|
|
|
|||
|
|
@ -56,12 +56,9 @@ struct PrimeFactorization{T<:Unsigned} <: Integer
|
|||
end
|
||||
PrimeFactorization(powers::Vector{T}) where {T<:Unsigned} = PrimeFactorization{T}(powers, one(Int8))
|
||||
|
||||
Base.copy(a::PrimeFactorization) = PrimeFactorization(copy(a.powers), a.sign)
|
||||
Base.:(==)(a::P, b::P) where {P<:PrimeFactorization} = a.powers == b.powers && a.sign == b.sign
|
||||
|
||||
# define our own factor function, returning an instance of PrimeFactorization
|
||||
function primefactor(n::Integer)
|
||||
iszero(n) && return PrimeFactorization(UInt8[], zero(UInt8))
|
||||
iszero(n) && return PrimeFactorization(UInt8[], zero(Int8))
|
||||
sn = n < 0 ? -one(Int8) : one(Int8)
|
||||
n = abs(n)
|
||||
m = length(factortable)
|
||||
|
|
@ -104,7 +101,16 @@ function primefactorial(n::Integer)
|
|||
@inbounds return PrimeFactorization(copy(factorialtable[n+1]))
|
||||
end
|
||||
|
||||
# Conversion from PrimeFactorization to `BigInt` using `bigprime`
|
||||
# Methods for PrimeFactorization:
|
||||
Base.copy(a::PrimeFactorization) = PrimeFactorization(copy(a.powers), a.sign)
|
||||
|
||||
Base.one(::Type{PrimeFactorization{T}}) where {T<:Unsigned} = PrimeFactorization(Vector{T}())
|
||||
Base.zero(::Type{PrimeFactorization{T}}) where {T<:Unsigned} = PrimeFactorization(Vector{T}(), zero(Int8))
|
||||
|
||||
Base.promote_rule(P::Type{<:PrimeFactorization},::Type{<:Integer}) = P
|
||||
Base.promote_rule(P::Type{<:PrimeFactorization},::Type{BigInt}) = BigInt
|
||||
|
||||
Base.convert(P::Type{<:PrimeFactorization}, n::Integer) = convert(P, primefactor(n))
|
||||
function Base.convert(::Type{BigInt}, a::PrimeFactorization)
|
||||
A = big(a.sign)
|
||||
for (n, e) in enumerate(a.powers)
|
||||
|
|
@ -114,6 +120,24 @@ function Base.convert(::Type{BigInt}, a::PrimeFactorization)
|
|||
end
|
||||
return A
|
||||
end
|
||||
Base.convert(::Type{PrimeFactorization{T}}, a::PrimeFactorization{T}) where {T<:Unsigned} = a
|
||||
Base.convert(::Type{PrimeFactorization{T1}}, a::PrimeFactorization{T2}) where {T1<:Unsigned, T2<:Unsigned} = PrimeFactorization(map(T1, a.powers), a.sign)
|
||||
|
||||
Base.:(==)(a::PrimeFactorization, b::PrimeFactorization) = a.powers == b.powers && a.sign == b.sign
|
||||
function Base.:<(a::PrimeFactorization, b::PrimeFactorization)
|
||||
if a.sign != b.sign
|
||||
return a.sign < b.sign
|
||||
elseif a.sign < 0
|
||||
return <(-b, -a)
|
||||
else
|
||||
ag, bg = divgcd(a, b)
|
||||
if length(ag.powers) <= length(bg.powers) && all(k->ag.powers[k]<bg.powers[k], 1:length(ag.powers))
|
||||
return true
|
||||
else
|
||||
return convert(BigInt, ag) < convert(BigInt, bg)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# Methods for PrimeFactorization: only fast multiplication, and lcm and gcd.
|
||||
# Addition and subtraction will require conversion to BigInt
|
||||
|
|
@ -181,15 +205,15 @@ end
|
|||
# given a list of numerators and denominators, compute the common denominator and
|
||||
# the rescaled numerator after putting all fractions at the same common denominator
|
||||
function commondenominator!(nums::Vector{P}, dens::Vector{P}) where {P<:PrimeFactorization}
|
||||
isempty(nums) && return one(P)
|
||||
# accumulate lcm of denominator
|
||||
den = copy(dens[1])
|
||||
den = PrimeFactorization(copy(dens[1].powers))
|
||||
for i = 2:length(dens)
|
||||
_vmax!(den.powers, dens[i].powers)
|
||||
end
|
||||
# rescale numerators
|
||||
for i = 1:length(nums)
|
||||
powers = _vsub!(_vadd!(nums[i].powers, den.powers), dens[i].powers)
|
||||
nums[i] = PrimeFactorization(powers, nums[i].sign)
|
||||
_vsub!(_vadd!(nums[i].powers, den.powers), dens[i].powers)
|
||||
end
|
||||
return den
|
||||
end
|
||||
|
|
|
|||
|
|
@ -1,5 +1,96 @@
|
|||
using WignerSymbols
|
||||
using Base.Test
|
||||
|
||||
# write your own tests here
|
||||
@test 1 == 2
|
||||
smalljlist = 0:1//2:10
|
||||
largejlist = 0:1//2:1000
|
||||
# test triangle coefficient
|
||||
for j1 in smalljlist, j2 in smalljlist
|
||||
for j3 = abs(j1-j2):(j1+j2)
|
||||
@test Δ(j1,j2,j3) ≈ sqrt(factorial(float(j1+j2-j3))*factorial(float(j1-j2+j3))*
|
||||
factorial(float(j2+j3-j1))/factorial(float(j1+j2+j3+1)))
|
||||
end
|
||||
end
|
||||
|
||||
# test 3j:
|
||||
#--------
|
||||
# test cg orthogonality
|
||||
for j1 in smalljlist, j2 in smalljlist
|
||||
d1 = 2*j1+1
|
||||
d2 = 2*j2+1
|
||||
M = zeros(d1*d2, d1*d2)
|
||||
ind1 = 1
|
||||
for m1 in -j1:j1, m2 in -j2:j2
|
||||
ind2 = 1
|
||||
for j3 in abs(j1-j2):(j1+j2)
|
||||
for m3 in -j3:j3
|
||||
M[ind1,ind2] = clebschgordan(j1,m1,j2,m2,j3,m3)
|
||||
ind2 += 1
|
||||
end
|
||||
end
|
||||
ind1 += 1
|
||||
end
|
||||
@test M'*M ≈ one(M)
|
||||
end
|
||||
|
||||
# test recurrence relations: Phys Rev E 57, 7274 (1998)
|
||||
for k = 1:10
|
||||
j2 = convert(BigFloat,rand(0:1//2:1000))
|
||||
j3 = convert(BigFloat,rand(0:1//2:1000))
|
||||
m2 = convert(BigFloat,rand(-j2:j2))
|
||||
m3 = convert(BigFloat,rand(-j3:j3))
|
||||
|
||||
for j in max(abs(j2-j3),abs(m2+m3))+1:(j2+j3)-1
|
||||
X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(m2+m3)^2))
|
||||
Y = (2*j+1)*((m2+m3)*(j2*(j2+1)-j3*(j3+1)) - (m2-m3)*j*(j+1))
|
||||
Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(m2+m3)^2))
|
||||
tol = 10*max(abs(X),abs(Y),abs(Z))*eps(BigFloat)
|
||||
@test (X*wigner3j(BigFloat,j+1,j2,j3,-m2-m3,m2,m3) + Z*wigner3j(BigFloat,j-1,j2,j3,-m2-m3,m2,m3))≈(-Y*wigner3j(BigFloat,j,j2,j3,-m2-m3,m2,m3)) atol=tol
|
||||
end
|
||||
end
|
||||
|
||||
# test 6j
|
||||
#----------
|
||||
# test orthogonality
|
||||
for j1 in smalljlist, j2 in smalljlist, j4 in smalljlist
|
||||
for j5 in max(abs(j1-j2-j4),abs(j1-j2+j4),abs(j1+j2-j4)):(j1+j2+j4)
|
||||
j6range = max(abs(j2-j4),abs(j1-j5)):min((j2+j4),(j1+j5))
|
||||
j3range = max(abs(j1-j2),abs(j4-j5)):min((j1+j2),(j4+j5))
|
||||
@test length(j6range) == length(j3range)
|
||||
M = zeros(length(j3range),length(j6range))
|
||||
for (k2,j6) in enumerate(j6range)
|
||||
for (k1,j3) in enumerate(j3range)
|
||||
M[k1,k2] = sqrt(2*j3+1)*sqrt(2*j6+1)*wigner6j(j1,j2,j3,j4,j5,j6)
|
||||
end
|
||||
end
|
||||
@test M'*M ≈ one(M)
|
||||
end
|
||||
end
|
||||
|
||||
# test special case
|
||||
for j1 in smalljlist, j2 in smalljlist
|
||||
j6 = 0
|
||||
j4 = j2
|
||||
j5 = j1
|
||||
for j3 in abs(j1-j2):(j1+j2)
|
||||
@test wigner6j(j1,j2,j3,j4,j5,j6) ≈ (-1)^(j1+j2+j3)/sqrt((2*j1+1)*(2*j2+1))
|
||||
end
|
||||
end
|
||||
|
||||
# test recurrence relations: Phys Rev E 57, 7274 (1998)
|
||||
for k = 1:10
|
||||
j2 = convert(BigFloat,rand(largejlist))
|
||||
j3 = convert(BigFloat,rand(largejlist))
|
||||
l1 = convert(BigFloat,rand(largejlist))
|
||||
l2 = convert(BigFloat,rand(abs(l1-j3):(l1+j3)))
|
||||
l3 = convert(BigFloat,rand(abs(l1-j2):min(l1+j2)))
|
||||
|
||||
for j in intersect(abs(j2-j3):(j2+j3), abs(l2-l3):(l2+l3))
|
||||
X = j*sqrt(((j+1)^2-(j2-j3)^2)*((j2+j3+1)^2-(j+1)^2)*((j+1)^2-(l2-l3)^2)*((l2+l3+1)^2 - (j+1)^2))
|
||||
Y = (2*j+1)*( j*(j+1)*( -j*(j+1) + j2*(j2+1) + j3*(j3+1) - 2*l1*(l1+1)) +
|
||||
l2*(l2+1)*( j*(j+1) + j2*(j2+1) - j3*(j3+1) ) +
|
||||
l3*(l3+1)*( j*(j+1) - j2*(j2+1) + j3*(j3+1) ) )
|
||||
Z = (j+1)*sqrt((j^2-(j2-j3)^2)*((j2+j3+1)^2-j^2)*(j^2-(l2-l3)^2)*((l2+l3+1)^2 - j^2))
|
||||
tol = 10*max(abs(X),abs(Y),abs(Z))*eps(BigFloat)
|
||||
@test (X*wigner6j(BigFloat,j+1,j2,j3,l1,l2,l3) + Z*wigner6j(BigFloat,j-1,j2,j3,l1,l2,l3))≈(-Y*wigner6j(BigFloat,j,j2,j3,l1,l2,l3)) atol=tol
|
||||
end
|
||||
end
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue