94 lines
3.3 KiB
Julia
94 lines
3.3 KiB
Julia
# compile using
|
||
# gfortran -shared -fPIC osbrac.f90 -o osbrac.so
|
||
# gfortran -shared -fPIC allosbrac.f90 -o allosbrac.so
|
||
|
||
println("OSBRACKETS test against Buck et al.")
|
||
|
||
n1 = [0,0,0,0,0,0,2,2,2,2];
|
||
l1 = [0,1,1,2,2,2,2,2,2,2];
|
||
n2 = [0,0,0,0,0,0,1,1,1,1];
|
||
l2 = [0,3,5,2,4,5,3,3,4,4];
|
||
N = [0,0,0,0,1,0,0,1,0,3];
|
||
L = [0,2,1,1,3,5,3,0,2,2];
|
||
n = [0,1,0,0,0,0,1,2,4,0];
|
||
l = [0,0,5,3,1,2,6,5,2,4];
|
||
Λ = [0,2,6,4,3,4,4,5,2,4];
|
||
|
||
function calculate_single_bracket(n1′, l1′, n2′, l2′, n1, l1, n2, l2, Λ) # Efros notation -- DON'T CONFUSE
|
||
|
||
Emax = max(l1′ + l2′ + 2 * (n1′ + n2′), l1 + l2 + 2 * (n1 + n2))
|
||
|
||
ϵ = (Emax - Λ) % 2
|
||
N = (l1 - l2 + Λ - ϵ) ÷ 2
|
||
M = (l1 + l2 - Λ - ϵ) ÷ 2
|
||
L = Λ
|
||
NQMAX = Emax
|
||
CO = 1/sqrt(2)
|
||
SI = 1/sqrt(2)
|
||
FIRSTCALL = true
|
||
# from source: BRAC(NP,N1P,MP,N1,N2) with dimensions BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2)
|
||
BRAC = zeros(Float64, L + 1, (NQMAX - L) ÷ 2 + 1,(NQMAX - L) ÷ 2 + 1,(NQMAX - L) ÷ 2 + 1,(NQMAX - L) ÷ 2 + 1)
|
||
|
||
@ccall "../OSBRACKETS/osbrac.so".osbrac_(N::Ref{Int32},M::Ref{Int32},L::Ref{Int32},NQMAX::Ref{Int32},CO::Ref{Float64},SI::Ref{Float64},FIRSTCALL::Ref{UInt8},BRAC::Ptr{Array{Float64}})::Cvoid
|
||
|
||
Nq = l1 + l2 + 2 * (n1 + n2)
|
||
Nq′ = l1′ + l2′ + 2 * (n1′ + n2′)
|
||
if Nq ≠ Nq′; return 0; end
|
||
N′ = (l1′ - l2′ + Λ - ϵ) ÷ 2
|
||
N1′ = n1′
|
||
M′ = (l1′ + l2′ - Λ - ϵ) ÷ 2
|
||
N1 = n1
|
||
N2 = n2
|
||
|
||
return BRAC[N′ + 1, N1′ + 1, M′ + 1, N1 + 1, N2 + 1]
|
||
end
|
||
|
||
println("OSBRAC results")
|
||
osbracs = calculate_single_bracket.(n1, l1, n2, l2, N, L, n, l, Λ)
|
||
display(osbracs)
|
||
|
||
function calculate_all_brackets(n1′, l1′, n2′, l2′, n1, l1, n2, l2, Λ) # Efros notation -- DON'T CONFUSE
|
||
Emax = max(maximum(l1′ + l2′ + 2 * (n1′ + n2′)), maximum(l1 + l2 + 2 * (n1 + n2)))
|
||
|
||
LMIN = minimum(Λ)
|
||
LMAX = maximum(Λ)
|
||
CO = 1/sqrt(2)
|
||
SI = 1/sqrt(2)
|
||
|
||
BRACs = Vector{Any}(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
|
||
|
||
out = Float64[]
|
||
|
||
for (Λ, n1′, l1′, n2′, l2′, n1, l1, n2, l2) in zip(Λ, n1′, l1′, n2′, l2′, n1, l1, n2, l2)
|
||
# BRAC(NP,N1P,MP,N1,N2,N,M,L)
|
||
|
||
Nq = l1 + l2 + 2 * (n1 + n2)
|
||
ϵ = (Nq - Λ) % 2
|
||
NP = (l1′ - l2′ + Λ - ϵ) ÷ 2
|
||
N1P = n1′
|
||
MP = (l1′ + l2′ - Λ - ϵ) ÷ 2
|
||
N1 = n1
|
||
N2 = n2
|
||
N = (l1 - l2 + Λ - ϵ) ÷ 2
|
||
M = (l1 + l2 - Λ - ϵ) ÷ 2
|
||
|
||
parity = mod1(l1 + l2, 2)
|
||
push!(out, BRACs[parity][1 + NP, 1 + N1P, 1 + MP, 1 + N1, 1 + N2, 1 + N, 1 + M, 1 + LMIN + Λ])
|
||
end
|
||
|
||
return out
|
||
end
|
||
|
||
println("ALLOSBRAC results")
|
||
allosbracs = calculate_all_brackets(n1, l1, n2, l2, N, L, n, l, Λ)
|
||
display(allosbracs)
|
||
|
||
println("Difference")
|
||
display(abs.(allosbracs - osbracs))
|