yapanuwan
/
BergEC-jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

OSBRACKETS implemented and NuclearToolkit.jl removed

This commit is contained in:
Nuwan Yapa 2024-04-19 15:18:36 -04:00
parent 278974f5b6
commit fce2f997c8
1 changed files with 29 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
ho_basis.jl
View File

@ -1,5 +1,4 @@
using SparseArrays using SparseArrays
using NuclearToolkit
using SpecialFunctions using SpecialFunctions
include("helper.jl") include("helper.jl")
@ -84,25 +83,45 @@ function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float
end end
function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ) function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
E_max = maximum(Es) Emax = maximum(Es)
j_max = 2 * E_max + 1
l_max = j_max
to = 0 # unused
dtri = NuclearToolkit.prep_dtri(l_max + 1); LMIN = Λ
dcgm0 = NuclearToolkit.prep_dcgm0(l_max); LMAX = Λ
d6j = nothing # NuclearToolkit.prep_d6j_int(E_max, j_max, to); CO = 1/sqrt(2)
SI = 1/sqrt(2)
BRACs = Vector{Array{Float64}}(undef, 2) # two arrays for both parities
for parity in 1:2
NQMAX = Emax - mod1(Emax, 2) + parity
# dimensions BRAC(0:LMAX,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:LMAX,0:(NQMAX-LMIN)/2,LMIN:LMAX)
BRACs[parity] = zeros(Float64, 1 + LMAX, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + LMAX, 1 + (NQMAX-LMIN) ÷ 2, 1 + LMAX-LMIN)
@ccall "../OSBRACKETS/allosbrac.so".allosbrac_(NQMAX::Ref{Int32},LMIN::Ref{Int32},LMAX::Ref{Int32},CO::Ref{Float64},SI::Ref{Float64},BRACs[parity]::Ptr{Array{Float64}})::Cvoid
end
mat = spzeros(length(Es), length(Es)) mat = spzeros(length(Es), length(Es))
s = hcat(Es, n1s, l1s, n2s, l2s) s = hcat(Es, n1s, l1s, n2s, l2s)
for idx in CartesianIndices(mat) for idx in CartesianIndices(mat)
(i, j) = Tuple(idx) (i, j) = Tuple(idx)
(Elhs, N, L, n, l) = s[i, :] (Elhs, N, L, n, l) = s[i, :]
(Erhs, n1, l1, n2, l2) = s[j, :] (Erhs, n1, l1, n2, l2) = s[j, :]
if Elhs == Erhs && triangle_ineq(L, l, Λ) && triangle_ineq(l1, l2, Λ) if Elhs == Erhs && triangle_ineq(L, l, Λ) && triangle_ineq(l1, l2, Λ)
mat[i, j] = NuclearToolkit.HObracket_d6j(N, L, n, l, n1, l1, n2, l2, Λ, 1.0, dtri, dcgm0, d6j) parity = mod1(Elhs, 2)
mat[i, j] = (-1)^(n1 + n2 + N + n) * pick_Moshinsky_brakcet(BRACs[parity], n1, l1, n2, l2, N, L, n, l, Λ)
end end
end end
return mat return mat
end end
function pick_Moshinsky_brakcet(BRAC, n1′, l1′, n2′, l2′, n1, l1, n2, l2, Λ) # Efros notation -- don't confuse
ϵ = (l1 + l2 - Λ) % 2
NP = (l1′ - l2′ + Λ - ϵ) ÷ 2
MP = (l1′ + l2′ - Λ - ϵ) ÷ 2
N = (l1 - l2 + Λ - ϵ) ÷ 2
M = (l1 + l2 - Λ - ϵ) ÷ 2
# BRAC(NP,N1P,MP,N1,N2,N,M,L)
return BRAC[1 + NP, 1 + n1′, 1 + MP, 1 + n1, 1 + n2, 1 + N, 1 + M, 1]
end