# 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

bracs = calculate_single_bracket.(n1, l1, n2, l2, N, L, n, l, Λ)
display(bracs)