0382 / CGcoefficient.jl

Compute CG coefficient, Racah coefficient, and Wigner 3j, 6j, 9j Symbols, and give the exact results.
https://0382.github.io/CGcoefficient.jl/
MIT License
1 stars 1 forks source link

f9j and gsl9j not correct? #3

Closed fgerick closed 2 years ago

fgerick commented 2 years ago
julia> using CGcoefficient

julia> using GSL_jll

julia> function gsl9j(dj1::Int, dj2::Int, dj3::Int,
                       dj4::Int, dj5::Int, dj6::Int,
                       dj7::Int, dj8::Int, dj9::Int)
           ccall(
               (:gsl_sf_coupling_9j, GSL_jll.libgsl),
               Cdouble,
               (Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint),
               dj1, dj2, dj3, dj4, dj5, dj6, dj7, dj8, dj9
           )
       end
gsl9j (generic function with 1 method)

julia> nineJ(1,2,3,4,5,6,3,6,9)
(1//1274)√(3//5)

julia> f9j(1,2,3,4,5,6,3,6,9)
0.0

julia> gsl9j(1,2,3,4,5,6,3,6,9)
0.0
0382 commented 2 years ago

In this package, the nineJ function's parameters are the exact angular momentum number. Because I think this function is always used for demonstration, I define its parameters can be Rational. As for the f9j and d9j functions, as well as the GSL library, the parameters are twice of the exact angular momentum number. Because the angular momentum number may be half integers, which may cause trouble in some calculation codes. People often use twice of the exact angular momentum number, in this way, we only need to deal with integers.

The right benchmark is like this:

julia> f9j(2,4,6,8,10,12,6,12,18)
0.0006080036650247122

julia> nineJ(1,2,3,4,5,6,3,6,9)
(1//1274)√(3//5)

julia> float(ans)
0.0006080036650247122

julia> f9j(1,2,3,4,5,6,3,6,9)
0.0

julia> nineJ(1//2,1,3//2,2,5//2,3,3//2,3,9//2)
ERROR: 2, 5/2 cannnot couple to 3
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:33
 [2] _d9j(dj1::BigInt, dj2::BigInt, dj3::BigInt, dj4::BigInt, dj5::BigInt, dj6::BigInt, dj7::BigInt, dj8::BigInt, dj9::B
igInt)
   @ CGcoefficient C:\Users\18322\.julia\packages\CGcoefficient\sUxev\src\WignerSymbols.jl:85
 [3] nineJ(j1::Rational{Int64}, j2::Int64, j3::Rational{Int64}, j4::Int64, j5::Rational{Int64}, j6::Int64, j7::Rational{
Int64}, j8::Int64, j9::Rational{Int64})
   @ CGcoefficient C:\Users\18322\.julia\packages\CGcoefficient\sUxev\src\WignerSymbols.jl:239
 [4] top-level scope
   @ REPL[8]:1
fgerick commented 2 years ago

Ah yes, I could've gotten this from the HalfInt notation... Thanks for clarifying.