3-body XZ in p-space (not working)

.gitignore

@@ -1,6 +1,3 @@
-# probably not recommended
-Project.toml
-
 # Compiled FORTRAN libraries
 *.so
 
@@ -11,7 +8,7 @@ hpc/*
 
 # VS Code files for those working on multiple tools
 .vscode/*
-!.vscode/settings.json
+# .vscode/settings.json
 !.vscode/tasks.json
 !.vscode/launch.json
 !.vscode/extensions.json

B2R_comparison.jl

@@ -1,51 +0,0 @@
-using Plots, LinearAlgebra
-include("p_space.jl")
-
-berggren_mesh = get_mesh((0, 0.4 - 0.15im, 0.8, 6), 128)
-csm_mesh = get_mesh((0, 8 - 3im), 512)
-
-for mesh in (berggren_mesh, csm_mesh)
-    p, w = mesh
-    mesh_E = p.*p ./ (2*0.5)
-
-    # ResonanceEC: Eq. (20)
-    V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
-
-    training_points = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
-    training_E = Vector{ComplexF64}(undef, length(training_points))
-    EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
-
-    for (j, c) in enumerate(training_points)
-        evals, evecs = eigen(get_H_matrix(V_system(c), p, w))
-        i = identify_pole_i(p, evals)
-        training_E[j] = evals[i]
-        EC_basis[:, j] = evecs[:, i]
-    end
-
-    EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC
-    EC_basis_w = EC_basis .* w
-    N_EC = transpose(EC_basis_w) * EC_basis
-
-    extrapolate_points = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
-
-    exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
-    extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))
-
-    for (j, c) in enumerate(extrapolate_points)
-        exact_E[j] = quick_pole_E(V_system(c))
-
-        H = get_H_matrix(V_system(c), p, w)
-        H_EC = transpose(EC_basis_w) * H * EC_basis
-        evals = eigvals(H_EC, N_EC)
-        i = argmin(abs.(evals .- exact_E[j]))
-        extrapolate_E[j] = evals[i]
-    end
-
-    scatter(real.(training_E), imag.(training_E), label="training")
-    scatter!(real.(exact_E), imag.(exact_E), label="exact")
-    scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
-    plot!(real.(mesh_E), imag.(mesh_E), label="contour")
-    xlims!(-0.3,0.3)
-    ylims!(-0.120,0.020)
-    savefig("temp/" * string(rand(UInt16)) * ".pdf")
-end

EC.jl

@@ -0,0 +1,173 @@
+using Statistics, SparseArrays, LinearAlgebra, Arpack, Plots
+include("common.jl")
+
+"EC model for a Hamiltonian family H(c) = H0 + c * H1"
+mutable struct affine_EC
+    H0::AbstractMatrix{ComplexF64}
+    H1::AbstractMatrix{ComplexF64}
+    weights::Vector{ComplexF64}
+    trained::Bool
+
+    H0_EC
+    H1_EC
+    N_EC
+
+    ensemble_size::Int
+    H0_EC_ensemble
+    H1_EC_ensemble
+    N_EC_ensemble
+
+    training_E::Vector{ComplexF64}
+    exact_E::Vector{ComplexF64}
+    extrapolated_E::Vector{ComplexF64}
+    extrapolated_CI::Vector{ComplexF64}
+
+    affine_EC(H0::AbstractMatrix{ComplexF64}, H1::AbstractMatrix{ComplexF64}, weights::Vector{ComplexF64}=ones(ComplexF64, size(H0, 1)); ensemble_size=0) =
+        new(H0, H1, weights, false, nothing, nothing, nothing, ensemble_size, Matrix[], Matrix[], Matrix[], ComplexF64[], ComplexF64[], ComplexF64[], ComplexF64[])
+end
+
+"Train an EC model for a given range of c values.
+If a list is provided for ref_eval, they are used as reference values for picking the closest eigenvalues at each sampling point.
+If a single number is provided for ref_eval, it is used as a reference for the first point, and the previous eigenvalue is used as the reference for each successive point.
+If orthonormalize_threshold > 0, Gram-Schmidt orthonormalization is performed, using this value as the threshold for dropping redundant vectors."
+function train!(EC::affine_EC, c_vals; ref_eval=-10.0, CAEC=false, gram_schmidt_threshold=0, verbose=true, tol=1e-5)
+    training_vecs = Vector{ComplexF64}[]
+
+    for c in c_vals
+        verbose && println("Training for c = $c")
+
+        global current_E
+        if ref_eval isa Number
+            current_E = ref_eval
+            ref_eval = nothing
+        elseif !isnothing(ref_eval)
+            current_E = popfirst!(ref_eval)
+        end
+
+        H = EC.H0 + c .* EC.H1
+
+        evals, evecs = eigs(H, sigma=current_E, maxiter=5000, tol=tol, ritzvec=true, check=1)
+        current_E = nearest(evals, current_E)
+
+        push!(EC.training_E, current_E)
+        push!(training_vecs, evecs[:, nearestIndex(evals, current_E)])
+    end
+
+    CAEC && append!(training_vecs, conj.(training_vecs))
+
+    (EC.H0_EC, EC.H1_EC, EC.N_EC) = get_reduced_matrices(EC, training_vecs, gram_schmidt_threshold; verbose=true)
+
+    for _ in 1:EC.ensemble_size
+        subsample = resample(length(training_vecs))
+        if gram_schmidt_threshold > 0
+            (H0_EC, H1_EC, N_EC) = get_reduced_matrices(EC, training_vecs, gram_schmidt_threshold, subsample; verbose=false)
+            push!(EC.H0_EC_ensemble, H0_EC)
+            push!(EC.H1_EC_ensemble, H1_EC)
+            push!(EC.N_EC_ensemble, N_EC)
+        else
+            push!(EC.H0_EC_ensemble, EC.H0_EC[subsample, subsample])
+            push!(EC.H1_EC_ensemble, EC.H1_EC[subsample, subsample])
+            push!(EC.N_EC_ensemble, EC.N_EC[subsample, subsample])
+        end
+    end
+
+    EC.trained = true
+end
+
+function get_reduced_matrices(EC::affine_EC, training_vecs, gram_schmidt_threshold, subsample=1:length(training_vecs); verbose=false)
+    vecs = deepcopy(training_vecs[subsample])
+
+    if gram_schmidt_threshold > 0; vecs = gram_schmidt!(vecs, EC.weights, gram_schmidt_threshold; verbose=verbose); end
+
+    EC_basis = hcat(vecs...)
+    weights_mat = spdiagm(EC.weights)
+    H0_EC = transpose(EC_basis) * weights_mat * EC.H0 * EC_basis
+    H1_EC = transpose(EC_basis) * weights_mat * EC.H1 * EC_basis
+    N_EC = transpose(EC_basis) * weights_mat * EC_basis
+
+    return (H0_EC, H1_EC, N_EC)
+end
+
+resample(n::Int) = rand(1:n, n) |> unique |> sort
+
+"Extrapolate using a trained EC model for a given range of c values
+If a list is provided for ref_eval, they are used as reference values for picking the closest eigenvalues at each point.
+If a single number is provided for ref_eval, it is used as a reference for the first point, and the previous eigenvalue is used as the reference for each successive point.
+If precalculated_exact_E is provided, ref_eval is ignored.
+If pseudo_inv_tol > 0, the GEVP is avoided using Moore-Penrose psuedoinverse, using this value as the relative tolerance for dropping redundant vectors."
+function extrapolate!(EC::affine_EC, c_vals; ref_eval=EC.training_E[end], pseudo_inv_tol=0, verbose=true, tol=1e-5, precalculated_exact_E=nothing)
+    @assert EC.trained "EC model must be trained using train() before extrapolation"
+
+    for c in c_vals
+        global current_E
+
+        if isnothing(precalculated_exact_E)
+            if ref_eval isa Number
+                current_E = ref_eval
+                ref_eval = nothing
+            elseif !isnothing(ref_eval)
+                current_E = popfirst!(ref_eval)
+            end
+
+            verbose && println("Extact solution for c = $c")
+            
+            H = EC.H0 + c .* EC.H1
+            evals, _ = eigs(H, sigma=current_E, maxiter=5000, tol=tol, ritzvec=false, check=1)
+            current_E = nearest(evals, current_E)
+        else
+            current_E = popfirst!(precalculated_exact_E)
+        end
+
+        push!(EC.exact_E, current_E)
+
+        verbose && println("Extrapolating for c = $c")
+
+        evals = get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, pseudo_inv_tol)
+        push!(EC.extrapolated_E, nearest(evals, current_E))
+
+        if EC.ensemble_size > 0
+            E_ensemble = zeros(ComplexF64, EC.ensemble_size)
+            for i in 1:EC.ensemble_size
+                evals = get_extrapolated_evals(EC.H0_EC_ensemble[i], EC.H1_EC_ensemble[i], EC.N_EC_ensemble[i], c, pseudo_inv_tol)
+                E_ensemble[i] = nearest(evals, current_E)
+            end
+            re_CI = std(real.(E_ensemble))
+            im_CI = std(imag.(E_ensemble))
+            push!(EC.extrapolated_CI, complex(re_CI, im_CI))
+        end
+    end
+end
+
+"Solve the GEVP with or without Moore-Penrose psuedoinverse"
+function get_extrapolated_evals(H0_EC, H1_EC, N_EC, c, pseudo_inv_tol)
+    H_EC = H0_EC + c .* H1_EC
+    if pseudo_inv_tol > 0
+        inv_N_EC = pinv(N_EC; atol=pseudo_inv_tol)
+        H_EC = inv_N_EC * H_EC
+        return eigvals(H_EC)
+    else
+        return eigvals(H_EC, N_EC)
+    end
+end
+
+"Export EC data as CSV file"
+exportCSV(EC::affine_EC, filename) = exportCSV(filename, (EC.training_E, EC.exact_E, EC.extrapolated_E), ("training", "exact", "extrapolated"))
+
+"Plot EC data and optionally save figure to a file"
+function plot(EC::affine_EC, save_fig_filename=nothing; basis_points=nothing, basis_contour=nothing, xlims=nothing, ylims=nothing)
+    scatter(real.(EC.training_E), imag.(EC.training_E), label="training")
+    scatter!(real.(EC.exact_E), imag.(EC.exact_E), label="exact", markercolor=:white)
+    if EC.ensemble_size > 0
+        scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), xerror=real.(EC.extrapolated_CI), yerror=imag.(EC.extrapolated_CI), label="extrapolated", m=:x)
+    else
+        scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), label="extrapolated", m=:x)
+    end
+
+    isnothing(basis_points) || scatter!(real.(basis_points), imag.(basis_points), m=:x, label="basis")
+    isnothing(basis_contour) || plot!(real.(basis_contour), imag.(basis_contour), label="contour")
+
+    isnothing(xlims) || xlims!(xlims...)
+    isnothing(ylims) || ylims!(ylims...)
+
+    isnothing(save_fig_filename) || savefig(save_fig_filename)
+end

EC_test.ipynb

@@ -1,100 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "using Plots, LinearAlgebra\n",
-    "include(\"p_space.jl\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "vertices = (0, 0.4 - 0.15im, 0.8, 6)\n",
-    "subdivisions = 128\n",
-    "p, w = get_mesh(vertices, subdivisions)\n",
-    "mesh_E = p.*p ./ (2*0.5)\n",
-    "\n",
-    "# ResonanceEC: Eq. (20)\n",
-    "V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "training_points = range(0.75, 0.45, 5)\n",
-    "training_E = Vector{ComplexF64}(undef, length(training_points))\n",
-    "EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))\n",
-    "\n",
-    "for (j, c) in enumerate(training_points)\n",
-    "    evals, evecs = eigen(get_H_matrix(V_system(c), p, w))\n",
-    "    i = identify_pole_i(p, evals)\n",
-    "    training_E[j] = evals[i]\n",
-    "    EC_basis[:, j] = evecs[:, i]\n",
-    "end\n",
-    "\n",
-    "scatter(real.(training_E), imag.(training_E), label=\"training\")\n",
-    "plot!(real.(mesh_E), imag.(mesh_E), label=\"contour\")\n",
-    "xlims!(0,1)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "extrapolate_points = range(0.40, 0.25, 5)\n",
-    "ref_E = 0.2 - 0.1im\n",
-    "\n",
-    "exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))\n",
-    "extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))\n",
-    "\n",
-    "EC_basis_w = EC_basis .* w\n",
-    "N_EC = transpose(EC_basis_w) * EC_basis\n",
-    "\n",
-    "for (j, c) in enumerate(extrapolate_points)\n",
-    "    exact_E[j] = quick_pole_E(V_system(c))\n",
-    "\n",
-    "    EC_basis_w = EC_basis .* w\n",
-    "    H = get_H_matrix(V_system(c), p, w)\n",
-    "    H_EC = transpose(EC_basis_w) * H * EC_basis\n",
-    "    evals = eigvals(H_EC, N_EC)\n",
-    "    i = argmin(abs.(evals .- ref_E))\n",
-    "    ref_E = evals[i]\n",
-    "    extrapolate_E[j] = evals[i]\n",
-    "end\n",
-    "\n",
-    "scatter(real.(training_E), imag.(training_E), label=\"training\")\n",
-    "scatter!(real.(exact_E), imag.(exact_E), label=\"exact\")\n",
-    "scatter!(real.(extrapolate_E), imag.(extrapolate_E), label=\"extrapolated\")\n",
-    "plot!(real.(mesh_E), imag.(mesh_E), label=\"contour\")\n",
-    "xlims!(0,1)"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Julia 1.9.0",
-   "language": "julia",
-   "name": "julia-1.9"
-  },
-  "language_info": {
-   "file_extension": ".jl",
-   "mimetype": "application/julia",
-   "name": "julia",
-   "version": "1.9.0"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}

OSBRACKETS/README

@@ -1,71 +0,0 @@
-The programs are written in Fortran-90. The double precision is adopted.
-
-The following six files make up the package: 
-
-allosbrac.f90
-testallosbrac.f90
-allosoutput
-  
-osbrac.f90
-testosbrac.f90
-osoutput
-
-These are two independent groups of files that correspond to the two versions of the present
-program. 
-
-In one version, the subroutine producing the brackets is named ALLOSBRAC. It is contained in 
-the file allosbrac.f90 along with the routines this subroutine uses. 
-
-In the other version, the subroutine producing the brackets is named OSBRAC. It is contained
-in the file osbrac.f90 along with the routines this subroutine uses. 
-
-To test the set of programs contained in the file allosbrac.f90, one may say, e.g.,
-gfortran -O2 testallosbrac.f90 allosbrac.f90. To test the set of programs contained in the 
-file osbrac.f90, one may say, e.g., gfortran -O2 testosbrac.f90 osbrac.f90. The files
-allosoutput and osoutput are respective output files that contain results of the tests. 
- 
-The files testallosbrac.f90 and testosbrac.f90 contain, respectively, only the programs
-TESTALLOSBRAC and TESTOSBRAC. The tests performed are described in comment lines 
-in these programs and in the accompanying CPC paper. 
-
-The parameters that are set at the moment in the programs TESTALLOSBRAC and TESTOSBRAC  
-are just those for which the results in the files allosoutput and osoutput are listed. With these
-parameters, the programs run less than a second on a notebook. 
--------------------------------------------------------------------------------------------------------------------------------------------
-As said above, the brackets are produced either by the subroutine ALLOSBRAC or by the
-subroutine OSBRAC. The subroutines contain generous comments.
--------------------------------------------------------------------------------------------------------------------------------------------
-In the case when ALLOSBRAC is employed, all the routines used to calculate the brackets, 
-which are contained in the file allosbrac.f90, are as follows: the subroutine ALLOSBRAC (with
-comments), the subroutines ARR, FLPHI (with comments), COEFREL (with comments), and the
-function WIGMOD. 
-
-The subroutine ALLOSBRAC calls for the ARR, COE, FLPHI, and COEFREL subroutines. FLPHI calls 
-for the function WIGMOD. The meaning of the parameters of the ALLOSBRAC subroutine is
-explained in comments at its beginning.
-------------------------------------------------------------------------------------------------------------------------------------------
-In the case when OSBRAC is employed, all the routines used to calculate the brackets, 
-which are contained in the file osbrac.f90, are as follows: the subroutine OSBRAC (with
-comments), the subroutines ARR, FLPHI (with comments), COEFREL (with comments), and the
-function WIGMOD. 
-
-The subroutine OSBRAC calls for  the ARR, COE, FLPHI, and COEFREL subroutines. FLPHI calls 
-for the function WIGMOD. The meaning of the parameters of the OSBRAC subroutine is explained
-in comments at its beginning. 
--------------------------------------------------------------------------------------------------------------------------------------------
-While it is clear from the comments in the beginnings of ALLOSBRAC and OSBRAC how to
-implement these subroutines, this is also obvious from examples listed in the above mentioned
-TESTALLOSBRAC and TESTOSBRAC programs. 
-
-In general, to compile codes that include the present programs one may say e.g., 
-f95 main.f90 ... allosbrac.f90 ... or, alternatively, f95 main.f90 ... osbrac.f90 ....
-
-
-
-
-
-
-
-
-
-

OSBRACKETS/allosbrac.f90

@@ -1,451 +0,0 @@
-    SUBROUTINE ALLOSBRAC(NQMAX,LMIN,LMAX,CO,SI,BRAC)
-    
-! COMMENTS ON USE:
-!---------------------------------------------------------------------------------------------------------------------------------------
-! THE OSCILLATOR BRACKETS <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\VARPHI ARE 
-! CALCULATED  HERE. THIS IS DONE WITH THE HELP OF EQS. (18), (20)-TYPE, AND (22) IN 
-! THE ACCOMPANYING CPC PAPER. THE QUANTITIES L1, L2, L1P, AND L2P IN THE ABOVE
-! DEFINITION OF THE BRACKETS ARE THE PARTIAL ANGULAR MOMENTA, L IS THE TOTAL 
-! ANGULAR MOMENTUM, AND N1, N2, N1P, N2P ARE  THE RADIAL QUANTUM NUMBERS.
-! IN THE NOTATION LIKE L1P, ETC., "P" SYMBOLIZES "PRIMED" HERE AND BELOW.
-
-! THE SUBROUTINE RETURNS THE ARRAY OF ALL THE BRACKETS SUCH THAT
-! L1+L2+2*(N1+N2).LE.NQMAX, LMIN \LE L \LE LMAX, AND THE BRACKETS PERTAIN TO STATES
-! OF THE SAME PARITY (-1)^(L1+L2) WHICH IS THE PARITY OF NQMAX. 
- 
-! THE L1, L2, L1P, AND L2P ORBITAL MOMENTA ARE EXPRESSED IN THE SUBROUTINE IN 
-! TERMS OF THE M, N, MP, AND NP VARIABLES DEFINED AS FOLLOWS
-! M=(L1+L2-L-NN)/2, N=(L1-L2+L-NN)/2, MP=(L1P+L2P-L-NN)/2, NP=(L1P-L2P+L-NN)/2
-! WHERE NN EQUALS 0 OR 1 WHEN, RESPECTIVELY, NQMAX-L IS EVEN OR ODD. ONE THEN HAS
-! L1 = M+N+NN,  L2 = M-N+L, L1P = MP+NP+NN,  L2P = MP-NP+L.
-! WHEN L1, L2, L1P AND L2P TAKE ALL THE VALUES ALLOWED AT GIVEN L AND GIVEN PARITY 
-! THE N AND NP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO L-NN AND 
-! THE M AND MP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO (NQMAX-L-NN)/2.
-! ONE ALSO HAS NQ=2*(M+N1+N2)+L+NN=2*(MP+N1P+N2P)+L+NN.  
-
-! ALL THE PARAMETERS OF THE SUBROUTINE BUT BRAC ARE INPUT ONES. THE MEANING
-! OF THE NQMAX, LMIN, AND LMAX PARAMETERS IS SEEN FROM THE ABOVE DESCRIPTION.
-
-! THE BRAC PARAMETER DESIGNATES THE ARRAY OF OUTPUT BRACKETS. IT IS OF THE FORM
-! BRAC(NP,N1P,MP,N1,N2,N,M,L) . (AS SAID ABOVE, L1=L1(N,M), L2=L2(N,M), L1P=L1P(NP.MP),
-! AND L2P=L2P(NP,MP). THE QUANTITY N2P IS DETERMINED BY THE EQUALITY 
-! MP+N1P+N2P=M+N1+N2.) THE ORDER OF THE ARGUMENTS OF BRAC CORRESPONDS TO 
-! NESTING OF THE LOOPS AT ITS COMPUTATION.
-
-! THE ROUTINE PARAMETERS CO AND SI ARE AS FOLLOWS, CO = COS(PHI) AND SI = SIN(PHI). 
-! THESE QUANTITIES DEFINE THE PSEUDO(!!!)ORTHOGONAL TRANSFORMATION 
-! XI1P = CO*XI1+SI*XI2, XI2P = SI*XI1-CO*XI2. 
-
-! (BRACKETS PERTAINING TO THE CASE OF THE ORTHOGONAL TRANSFORMATION 
-! XI1P = CO*XI1+SI*XI2, XI2P = -SI*XI1+CO*XI2  ARE SIMPLY EXPRESSED IN TERMS OF 
-! THOSE CALCULATED IN THE PRESENT PROGRAM, SEE THE ACCOMPANYING CPC PAPER.
-
-! THE MEANING OF THE ARRAYS, OTHER THAN BRAC, ENTERING THE DIMENSION LIST IS AS
-! FOLLOWS. B IS A SUBSIDIARY ARRAY TO PERFORM THE RECURSION. 
-
-! FAC(I)=I!, DFAC(I)=(2I+1)!!, AND DEFAC(I)=(2I)!!. RFAC(I)=SQRT((2I)!)/I! 
-! THESE ARRAYS ARE PRODUCED BY THE "ARR" ROUTINE WHICH IS CALLED FROM THE PRESENT
-! ROUTINE AND WHICH IS CONTAINED IN THE PRESENT FILE. THE SET UPPER BOUNDS OF THESE
-! ARRAYS ARE SUFFICIENT FOR THE COMPUTATION. THESE ARRAYS ARE USED IN THE PRESENT
-! ROUTINE AND IN THE "FLPHI" ROUTINE.
-
-! THE A ARRAY IS A(N,L)=(-1)**N/SQRT(DEFAC(N)*DFAC(N+L)). IT APPEARS BOTH IN THE N1=N2=0
-! BRACKETS AND IN THE RELATION BETWEEN THE < | > AND [ | ] TYPE BRACKETS.
-
-! THE BI ARRAY IS BI(M,N)=FAC(N)/[FAC(M)*FAC(N-M)]. IT ENTERS THE 3J SYMBOLS. IT IS USED
-! IN THE "FLPHI" ROUTINE AND IN THE FUNCTION WIGMOD.
-
-! THE FL(NP,MP,N) ARRAY REPRESENTS THE QUANTITY 
-! [(2*L1P+1)*(2*L2P+1)]^{1/2}*F_L^VARPHI WHERE F_L^VARPHI IS GIVEN BY EQ. (23) IN THE
-! ACCOMPANYING CPC PAPER. THIS ARRAY IS PRODUCED IN ADVANCE BY THE
-! "FLPHI" ROUTINE WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS
-! CONTAINED IN THE PRESENT FILE. THE VARIABLES MP, NP, AND N ARE DEFINED ABOVE ALONG
-! WITH THEIR UPPER BOUNDS. 
-
-! THE PSIP(P,Q) AND PSIM(P,Q) ARRAYS REPRESENT THE QUANTITIES (4) IN THE ACCOMPANYING
-! CPC PAPER. THEY ARE PRODUCED IN ADVANCE BY THE "COEFREC" ROUTINE
-! WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS CONTAINED IN THE PRESENT
-! FILE. THEIR ARGUMENTS P AND Q TAKE THE VALUES L1P, L1P+/-1 AND L2P, L2P+/-1,
-! RESPECTIVELY. THE BOUNDS OF THE PSIP AND PSIM ARRAYS ARE SUCH THAT ALL THE P AND Q
-! VALUES REQUIRED TO PERFORM THE RECURSION ARE PROVIDED.
-
-! THE MENTIONED "FLPHI" ROUTINE USES THE FUNCTION WIGMOD WHICH IS ALSO CONTAINED
-! IN THE PRESENT FILE.
-!----------------------------------------------------------------------------------------------------------------------------------------
-    DOUBLE PRECISION BRAC,FL,PSIP,PSIM,A,FAC,DFAC,DEFAC,RFAC,BI,CO,SI,T,CO2,SI2,SC,FA,&
-    PLFACT,D2L12P,TN,PA,S,FA1,FA2 
-    DIMENSION& 
-    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),&
-! THE ARRAY IS OF THE FORM BRAC(NP,N1P,MP,N1,N2,N,M,L)      
-    FL(0:LMAX,0:(NQMAX-LMIN)/2,0:LMAX),&
-    PSIP((NQMAX+LMAX)/2+1,(NQMAX+LMAX)/2+1),&
-    PSIM((NQMAX+LMAX)/2+1,(NQMAX+LMAX)/2+1),A(0:NQMAX/2+1,0:NQMAX),&
-    FAC(0:2*NQMAX+1),DFAC(0:NQMAX+1),DEFAC(0:NQMAX/2+1),RFAC(0:NQMAX),&
-    BI(0:2*NQMAX+1,0:2*NQMAX+1) 
-    IF (NQMAX.GT.84) THEN
-        WRITE(6,*)'IN THE R(*8) COMPUTATION NQMAX SHOULD NOT EXCEED 84' 
-        STOP
-    ENDIF   
-    IF (CO.EQ.0.D0.OR.SI.EQ.0.D0) THEN
-        WRITE(6,*)'THE PROGRAM IS OF NO USE AT ZERO COS(PHI) OR SIN(PHI) VALUES'
-        WRITE(6,*)'COS(PHI)=',CO,' SIN(PHI)=',SI
-        STOP
-    ENDIF
-    IF (NQMAX.LT.LMAX) THEN
-        WRITE(6,*)'NQMAX=',NQMAX,' L=',LMAX
-        WRITE(6,*)'L SHOULD NOT EXCEED NQMAX'
-        STOP
-    ENDIF
-    T=SI/CO
-    CO2=CO**2 
-    SI2=SI**2
-    SC=SI*CO
-    CALL ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
-    CALL COE(2*NQMAX+1,FAC,BI)
-!   THE A_{NL} ARRAY:
-    DO NI=0,NQMAX/2+1
-        NA=NI-2*(NI/2)
-        K=1
-        IF(NA.EQ.1)K=-1
-        DO LI=0,NQMAX-NI
-            A(NI,LI)=K/SQRT(DEFAC(NI)*DFAC(NI+LI))
-        ENDDO
-    ENDDO             
-! BRACKETS [N_1'L_1'N_2'L_2'|0L_10L_2]_L^\PHI TO START THE RECURSION, EQ. (22) IN  
-! THE ACCOMPANYING CPC PAPER. 
-! THESE BRACKETS ARE REPRESENTED AS BRAC(NP,N1P,MP,0,0,N,M,L). 
-! SUBSIDIARY QUANTITIES: 
-    DO L=LMIN,LMAX        
-        CALL FLPHI(NQMAX,LMAX,LMIN,L,T,BI,RFAC,FL)
-! THIS SUBROUTINE PRODUCES THE ARRAY FL USED BELOW.    
-        NQML=NQMAX-L
-        NQMLD2=NQML/2
-        NN=NQML-2*(NQMLD2)
-        CALL COEFREC((NQMAX+LMAX)/2+1,L,NN,NQMAX,PSIP,PSIM) 
-        NMAX=L-NN 
-        MMAX=(NQML-NN)/2
-        DO M=0,MMAX
-            DO N=0,NMAX     
-                L1=M+N+NN 
-                L2=M-N+L
-                L1L2=L1+L2
-                L1L2ML=L1L2-L    
-                FA=SQRT((2.D0*L1+1)*(2*L2+1))*CO**L1L2
-                PLFACT=SQRT(FAC(L1L2+L+1)*FAC(L1L2ML))
-! COMPUTATION OF THE N1=N2=0 BRACKETS. THEY ARE REPRESENTED AS 
-! BRAC(NP,N1P,MP,0,0,N,M,L). 
-                MPMAX=(L1L2ML-NN)/2
-                DO MP=0,MPMAX
-                    L12P=2*MP+L+NN
-! THIS IS BY DEFINITION, L12P=L1P+L2P
-                    D2L12P=1.D0/2.D0**L12P
-                    N12P=MPMAX-MP
-! INDEED, 2*N12P=2*(N1P+N2P)=NQ-L12P. FOR THE BRACKETS WE COMPUTE NOW 
-! WE HAVE NQ=L1L2. THEN 2*N12P=L1L2-L12P. SINCE L1L2=2*MPMAX-L-NN AND 
-! L12P=2*MP+L+NN, ONE GETS N12P=MPMAX-MP.
-                    DO N1P=0,N12P
-                        N1PA=N1P-2*(N1P/2)
-                        IF (N1PA.EQ.0) THEN
-                            K=1
-                        ELSE
-                            K=-1
-                        ENDIF
-! K=(-1)**N1P
-                        N2P=N12P-N1P
-                        TN=T**N12P*K
-                        DO NP=0,NMAX
-                            L1P=MP+NP+NN
-                            L2P=L12P-L1P  
-! BY DEFINITION ONE ALSO HAS: L2P=L2P(MP,NP)=MP-NP+L
-                            PA=A(N1P,L1P)*A(N2P,L2P)
-                            BRAC(NP,N1P,MP,0,0,N,M,L)=FL(NP,MP,N)*FA*TN*D2L12P&
-                            *PLFACT*PA*PA
-                        ENDDO ! NP
-                    ENDDO ! N1P
-                ENDDO ! MP
-! RECURSION TO OBTAIN THE GENERAL FORM BRACKETS 
-! [N_1'L_1'N_2'L_2'|N_1L_1N_2L_2]_L^\PHI. 
-                N12MAX=(NQMAX-L1L2)/2
-                MPMAX0=MPMAX
-! THE N1=0, N2-1-->N2 RECURSION, EQ. (20)  IN THE ACCOMPANYING CPC PAPER WITH THE
-! MODIFICATION POINTED OUT THERE: 
-                DO N2=1,N12MAX
-                    MPMAX=MPMAX+1
-! THIS MPMAX VALUE IS (NQ-L-NN)/2 WHERE NQ=L1+L2+2*N2.
-                    DO MP=0,MPMAX
-                        N12P=MPMAX-MP
-                        DO N1P=0,N12P
-! MP=(L1P+L2P-L-NN)/2. THEREFORE, N12P=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.
-                            DO NP=0,NMAX
-                                L1P=MP+NP+NN
-                                L2P=MP-NP+L  
-                                S=0.D0 
-! BELOW THE RESTRICTIONS ARE IMPOSED ON THE BRAC ARRAY ENTERING THE RIGHT-HAND
-! SIDE OF THE RECURRENCE FORMULAE. THIS ARRAY IS OF THE FORM BRAC(K1,K2,K3,...) WHERE
-! K1=K1(NP), K2=N2(N1P), AND K3=K3(MP). NAMELY, K1=NP, OR NP+1, OR NP-1; 
-! K2=N1P OR N1P-1; K3=MP, OR MP-1, OR MP+1. (THE COMBINATION K2=N1P AND K3=MP+1
-! DOES NOT ARISE IN THE RECURRENCE FORMULAE.) THE QUANTITIES K1, K2, AND K3
-! REPRESENT, RESPECTIELY, NP, N1P, AND MP VALUES AT THE PRECEDING STAGE OF THE
-! RECURSION. THE RESTRICTIONS ENSURE THAT K1, K2, AND K3 RANGE WITHIN THE LIMITS
-! PERTAINING TO THE BRAC ARRAY OBTAINED AT THAT PRECEDING STAGE OF THE RECURSION. 
-
-! THUS K1, K2, AND K3 SHOULD BE NON-NEGATIVE. THEREFORE, WHEN K1=NP-1 THE VALUE
-! NP=0 IS TO BE EXCLUDED, WHEN K2=N1P-1 THE VALUE N1P=0 IS TO BE EXCLUDED, AND
-! WHEN K3=MP-1 THE VALUE MP=0 IS TO BE EXCLUDED. 
-
-! FURTHERMORE, IT SHOULD BE K1 \LE NMAX AND AT THE SAME TIME ONE HAS NP \LE NMAX.
-! THEREFORE, WHEN K1=NP+1 THE VALUE NP=NPMAX IS TO BE EXCLUDED.
-
-! IN ADDITION, IT SHOULD BE K2+K3 \LE MPMAX-1 WHILE ONE HAS N1P+MP \LE MPMAX. WHEN
-! K2=N1P AND K3=MP-1, OR K2=N1P-1 AND K3=MP-1, OR K2=N1P-1 AND K3=MP ONE
-! AUTOMATICALLY HAS N1P+MP \LE MPMAX. THEREFORE, IN THESE CASES THE CONDITION
-! K2+K3 \LE MPMAX-1 DOES NOT CREATE ANY RESTRICTIONS. 
-! BUT WHEN K2=N1P AND K3=MP, OR K2=N1P-1 AND K3=MP+1 THE CASE N1P+MP=MPMAX
-! IS TO BE FORBIDDEN. 
-
-! IN THE RECURRENCE FORMULAE BELOW THE DESCRIBED RESTRICTIONS ARE
-! IMPOSED.   
-                                IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,0,N2-1,N,M,L)*PSIP(L1P,L2P)
-                                IF (N1P.NE.0.AND.N1P.NE.N12P) S=S+&
-                                BRAC(NP,N1P-1,MP+1,0,N2-1,N,M,L)*PSIP(L1P+1,L2P+1)
-! RECALL THAT N12P=MPMAX-MP.                                
-                                IF (NP.NE.NMAX.AND.N1P.NE.0) S=S-&                               
-                                BRAC(NP+1,N1P-1,MP,0,N2-1,N,M,L)*PSIM(L1P+1,L2P)                                
-                                IF (NP.NE.0.AND.N1P.NE.N12P) S=S-&
-                                BRAC(NP-1,N1P,MP,0,N2-1,N,M,L)*PSIM(L1P,L2P+1)                                
-                                S=S*SC                                
-                                IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,0,N2-1,N,M,L)*SI2
-                                IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,0,N2-1,N,M,L)*CO2
-                                BRAC(NP,N1P,MP,0,N2,N,M,L)=S
-                            ENDDO ! NP
-                        ENDDO ! N1P
-                    ENDDO ! MP
-                ENDDO ! N2
-! N1-1-->N1 RECURSION, EQ. (20) IN THE ACCOMPANYING CPC PAPER:
-                DO N2=0,N12MAX
-                    MPMAX=MPMAX0+N2
-! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2=MPMAX0+N2 SINCE NQ=L1+L2+2*N2.
-! THE RECURSION:        
-                    DO N1=1,N12MAX-N2
-                        MPMAX=MPMAX+1
-! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2 AND NQ=L1+L2+2*N2+2*N1.            
-                        DO MP=0,MPMAX
-                            N12P=MPMAX-MP
-                            DO N1P=0,N12P
-! MP=(L1P+L2P-L-NN)/2. THEREFORE, N1PMAX=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.              
-                                DO NP=0,NMAX
-                                    L1P=MP+NP+NN
-                                    L2P=MP-NP+L  
-                                    S=0.D0                 
-                                    IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,N1-1,N2,N,M,L)*PSIP(L1P,L2P)
-                                    IF (N1P.NE.0.AND.N1P.NE.N12P) S=S+&
-                                    BRAC(NP,N1P-1,MP+1,N1-1,N2,N,M,L)*PSIP(L1P+1,L2P+1)
-                                    IF (NP.NE.NMAX.AND.N1P.NE.0) S=S-&
-                                    BRAC(NP+1,N1P-1,MP,N1-1,N2,N,M,L)*PSIM(L1P+1,L2P)
-                                    IF (NP.NE.0.AND.N1P.NE.N12P) S=S-&
-                                    BRAC(NP-1,N1P,MP,N1-1,N2,N,M,L)*PSIM(L1P,L2P+1)
-! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
-                                    S=-S*SC
-                                    IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,N1-1,N2,N,M,L)*CO2
-                                    IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,N1-1,N2,N,M,L)*SI2
-! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
-                                    BRAC(NP,N1P,MP,N1,N2,N,M,L)=S
-                                ENDDO ! NP
-                            ENDDO ! N1P
-                        ENDDO ! MP
-                    ENDDO ! N1
-                ENDDO ! N2
-! RENORMALIZATION OF THE BRACKETS, EQ. (18) IN THE ACCOMPANYING CPC PAPER:
-                DO N2=0,N12MAX
-                    MPMAX1=MPMAX0+N2
-                    DO N1=0,N12MAX-N2
-                        MPMAX=MPMAX1+N1
-                        FA1=A(N1,L1)*A(N2,L2)
-                        DO MP=0,MPMAX
-                            N12P=MPMAX-MP
-                            DO N1P=0,MPMAX-MP
-                                N2P=N12P-N1P
-                                DO NP=0,NMAX
-                                    L1P=MP+NP+NN
-                                    L2P=MP-NP+L                       
-                                    FA2=A(N1P,L1P)*A(N2P,L2P)
-                                    BRAC(NP,N1P,MP,N1,N2,N,M,L)=BRAC(NP,N1P,MP,N1,N2,N,M,L) &
-                                    *FA1/FA2
-                                ENDDO ! NP
-                            ENDDO ! N1P 
-                        ENDDO ! MP
-                    ENDDO ! N1
-                ENDDO ! N2                      
-            ENDDO ! N
-        ENDDO ! M
-    ENDDO ! L
-    RETURN
-    END 
- 
-    SUBROUTINE ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
-! FAC(I), DFAC(I),DEFAC(I), AND RFAC(I) ARE, RESPECTIVELY, THE QUANTITIES
-! I!, (2I+1)!!, (2I)!!, AND SQRT((2*I)!)/I!
-    DOUBLE PRECISION FAC,DFAC,DEFAC,RFAC
-    DIMENSION FAC(0:2*NQMAX+1),DFAC(0:NQMAX+1),DEFAC(0:NQMAX/2+1),RFAC(0:NQMAX)
-    FAC(0)=1.D0
-    DFAC(0)=1.D0
-    DEFAC(0)=1.D0
-    RFAC(0)=1.D0
-    DO I=1,2*NQMAX+1 
-        FAC(I)=FAC(I-1)*I 
-    ENDDO            
-    DO I=1,NQMAX+1           
-        DFAC(I)=DFAC(I-1)*(2*I+1)
-    ENDDO    
-    DO I=1,NQMAX/2+1
-        DEFAC(I)=DEFAC(I-1)*2*I
-    ENDDO
-    DO I=1,NQMAX    
-        RFAC(I)=RFAC(I-1)*2*SQRT(1-0.5D0/I)    
-    ENDDO
-    RETURN
-    END
-    
-
-    SUBROUTINE FLPHI(NQMAX,LMAX,LMIN,L,T,BI,RFAC,FL)
-! PROVIDES THE QUANTITY SQRT((2L1P+1)*(2L2P+1))*F_L^\VARPHI, SEE EQ. (23) IN THE 
-! ACCOMPANYING PAPER, IN THE FORM OF THE FL ARRAY.
-! USES THE FUNCTION WIGMOD. 
-! THIS FLPHI DIFFERS FROM THAT IN THE OTHER FILE osbrac.f90.
-    DOUBLE PRECISION BI,RFAC,FL,T,T2,SQP,F,WIGMOD,z
-    DIMENSION&
-    FL(0:LMAX,0:(NQMAX-LMIN)/2,0:LMAX),BI(0:2*NQMAX+1,0:2*NQMAX+1),RFAC(0:NQMAX)
-! THE OUTPUT IS FL(NP,MP,N) WHERE MP=(L1P+L2P-L-NN)/2, L1P+L2P=L1T+L2T,
-! NP=(L1P-L2P+L-NN)/2, AND N=(L1-L2+L-NN)/2, L1-L2=L1T-L2T. 
-    T2=T*T    
-    NQMAL=NQMAX-L
-    MPMAX=NQMAL/2
-    NN=NQMAL-2*MPMAX
-    LMNN=L-NN 
-    LL3=2*L
-    DO N=0,LMNN
-        L1ML2=2*N-LMNN 
-        DO MP=0,MPMAX
-            L1PL2P=2*MP+L+NN
-            DO  NP=0,LMNN
-                L1PML2P=2*NP-LMNN            
-                L1P=(L1PL2P+L1PML2P)/2
-                L2P=(L1PL2P-L1PML2P)/2
-                M3=L2P-L1P
-                SQP=SQRT((2*L1P+1.D0)*(2*L2P+1))
-                L1T=(L1PL2P+L1ML2)/2
-                L2T=(L1PL2P-L1ML2)/2
-                N3=L1T+L2T-L
-                IMIN=ABS(L1T-L1P)
-                IMAX=MIN(L1P+L1T,L2P+L2T)
-                NAL1=(L1T+L1P-IMIN)/2
-                NAL2=(L1T-L1P+IMIN)/2
-                NAL3=(L2T-L2P+IMIN)/2
-                NAL4=(L2T+L2P-IMIN)/2
-                NAL4P=NAL4-2*(NAL4/2)
-                F=0.D0
-                z=t**imin
-                if(nal4p.ne.0)z=-z
-                DO I=IMIN,IMAX,2                          
-                    F=F+RFAC(NAL1)*RFAC(NAL2)*RFAC(NAL3)*RFAC(NAL4)*&
-                    WIGMOD(L1T,L2T,L,L1P-I,I-L2P,BI,2*NQMAX+1)*z                   
-                    NAL1=NAL1-1
-                    NAL2=NAL2+1
-                    NAL3=NAL3+1
-                    NAL4=NAL4-1    
-                    z=-z*t2
-                ENDDO
-                FL(NP,MP,N)=F*SQRT((2*L1P+1)*(2*L2P+1)&
-                *BI(L1T+L-L2T,2*L1T)*BI(N3,2*L2T)/((LL3+1)*BI(N3,L1T+L2T+L+1)*BI(L+M3,LL3)))
-            ENDDO
-        ENDDO
-    ENDDO
-    RETURN
-    END
-          
-   FUNCTION WIGMOD(L1,L2,L3,M1,M2,BI,NMAX)
-! THE 3J SYMBOL IN TERMS OF BINOMIAL COEFFICIENTS WITHOUT THE FACTOR
-! SQRT(F) WHERE F IS AS FOLLOWS,
-! F=BI(L1+L3-L2,2*L1)*BI(N3,2*L2)/((LL3+1)*BI(N3,L1+L2+L3+1)*BI(L3+M3,LL3)) 
-! WITH N3=L1+L2-L3, LL3=2*L3, AND M3=-M1-M2.   
-    DOUBLE PRECISION WIGMOD,BI,S
-    DIMENSION BI(0:NMAX,0:NMAX)
-    M3=-M1-M2
-    N3=L1+L2-L3
-    LM1=L1-M1
-    LP2=L2+M2
-    KMIN=MAX(0,L1-L3+M2,L2-L3-M1)
-    KMAX=MIN(N3,LM1,LP2)
-    S=0.d0
-    NPH=1
-    DO K=KMIN,KMAX
-        S=S+NPH*BI(K,N3)*BI(LM1-K,L1+L3-L2)*BI(LP2-K,L2+L3-L1)
-        NPH=-NPH
-    ENDDO 
-    WIGMOD=S/SQRT(BI(LM1,2*L1)*BI(LP2,2*L2))
-    NY=KMIN+L1-L2-M3
-    NYP=NY-2*(NY/2)
-    IF(NYP.NE.0)WIGMOD=-WIGMOD
-    RETURN
-    END    
-
-    SUBROUTINE COE(NMAX,FAC,BI)
-    DOUBLE PRECISION FAC,BI
-    DIMENSION FAC(0:NMAX),BI(0:NMAX,0:NMAX)
-    DO N=0,NMAX
-        DO M=0,N/2
-            BI(M,N)=FAC(N)/(FAC(M)*FAC(N-M))
-            BI(N-M,N)=BI(M,N)
-        ENDDO
-    ENDDO  
-    RETURN
-    END
-        
-    SUBROUTINE COEFREC(MAXCOE,L,NN,NQMAX,PSIP,PSIM)
-! CALCULATES THE COEFFICIENTS OF THE RECURSION FORMULA, EQ. (21) IN THE ACCOMPANYING
-! CPC PAPER
-! MAXCOE EQUALS INT((NQMAX+LMAX)/2)+1 AT CALLS OF THIS ROUTINE.
-    DOUBLE PRECISION PSIP,PSIM
-    DIMENSION PSIP(MAXCOE,MAXCOE),PSIM(MAXCOE,MAXCOE)
-    LP=L+1
-    LM=L-1
-    DO M1=1,(NQMAX+L)/2+1
-        DO M2=1,M1
-            M1P2=M1+M2
-            NA=M1P2+L
-            NNA=NA-2*(NA/2)
-            M1M2=M1-M2
-! CALCULATION OF PSIP. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P AND L2P OR 
-! L1P+1 AND L2P+1.)
-            IF (NNA.EQ.NN.AND.M1P2.GE.L.AND.M1P2.LE.NQMAX.AND.&
-            ABS(M1M2).LE.L) THEN
-                IF (M1P2.GT.LP) THEN
-                    M1P2L=M1P2+L
-                    M1P2ML=M1P2-L
-                    PSIP(M1,M2)=SQRT(M1P2L*(M1P2L+1.D0)*M1P2ML*(M1P2ML-1)/((4*M1*M1-1)&
-                    *(4*M2*M2-1)))
-                ELSE
-                    PSIP(M1,M2)=0.D0
-                ENDIF
-                PSIP(M2,M1)=PSIP(M1,M2) 
-            ENDIF
-! CALCULATION OF PSIM. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P+1 AND L2P
-! OR L2P+1 AND L1P.)
-            IF (NNA.NE.NN.AND.M1P2.GE.LP.AND.M1P2.LE.NQMAX&
-            .AND.M1M2.LE.LP.AND.M1M2.GE.-LM) THEN
-                IF (M1M2.LT.L) THEN
-                    M1M2L=M1M2+L
-                    M1M2ML=M1M2-L 
-                    PSIM(M1,M2)=SQRT(M1M2L*(M1M2L+1.D0)*M1M2ML*(M1M2ML-1)/&
-                    ((4*M1*M1-1)*(4*M2*M2-1)))
-                ELSE
-                    PSIM(M1,M2)=0.D0
-                ENDIF
-                PSIM(M2,M1)=PSIM(M1,M2) 
-            ENDIF                                               
-        ENDDO
-    ENDDO 
-    RETURN
-    END             
-    

OSBRACKETS/osbrac.f90

@@ -1,467 +0,0 @@
-    SUBROUTINE OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC) 
-    
-! COMMENTS ON USE:
-!---------------------------------------------------------------------------------------------------------------------------------------
-! THE OSCILLATOR BRACKETS <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\VARPHI ARE 
-! CALCULATED  HERE. THIS IS DONE WITH THE HELP OF EQS. (18), (20)-TYPE, AND (22) IN 
-! THE ACCOMPANYING CPC PAPER. THE QUANTITIES L1, L2, L1P, AND L2P IN THE ABOVE
-! DEFINITION OF THE BRACKETS ARE THE PARTIAL ANGULAR MOMENTA, L IS THE TOTAL 
-! ANGULAR MOMENTUM, AND N1, N2, N1P, N2P ARE  THE RADIAL QUANTUM NUMBERS.
-! IN THE NOTATION LIKE L1P, ETC., "P" SYMBOLIZES "PRIMED" HERE AND BELOW. 
-
-! N,M,L, AND NQMAX ARE INPUT PARAMETERS OF THE SUBROUTINE, AND BRAC IS THE
-! OUTPUT PARAMETER REPRESENTING THE ARRAY OF CALCULATED BRACKETS. THE PARAMETERS
-! N AND M REPRESENT THE L1 AND L2 QUANTUM NUMBERS, SEE BELOW, SO THAT L1=L1(N,M)
-! AND L2=L2(N,M).
-
-! LET US USE THE NOTATION NQ=L1+L2+2*(N1+N2).  THE ARRAY BRAC CONTAINS ALL BRACKETS
-! EXISTING AT GIVEN INPUT L, L1, AND L2 VALUES AND HAVING N1 AND N2 SUCH THAT
-! NQ.LE.NQMAX. IT IS CLEAR THAT NQ AND NQMAX PERTAINING TO ALL THESE BRACKETS ARE OF
-! THE SAME PARITY EQUAL TO (-1)^(L1+L2).
-
-! IN OTHER WORDS, THE COMPUTED BRACKETS PERTAIN TO STATES WITH GIVEN N, M, AND L
-! VALUES SUCH THAT NQ.LE.NQMAX. AUTOMATICALLY, THESE STATES ARE OF THE SAME PARITY
-! AS THE PARITY OF NQMAX. (AND THIS PARITY EQUALS (-1)^(L1+L2).) 
- 
-! THE L1, L2, L1P, AND L2P ORBITAL MOMENTA ARE EXPRESSED IN TERMS OF THE N, M, NP, AND
-! MP VARIABLES AS FOLLOWS
-! N=(L1-L2+L-NN)/2, M=(L1+L2-L-NN)/2, NP=(L1P-L2P+L-NN)/2, MP=(L1P+L2P-L-NN)/2 
-! WHERE NN EQUALS 0 OR 1 WHEN, RESPECTIVELY, NQMAX-L IS EVEN OR ODD. ONE THEN HAS
-! L1 = M+N+NN,  L2 = M-N+L, L1P = MP+NP+NN,  L2P = MP-NP+L.
-! WHEN L1, L2, L1P AND L2P TAKE ALL THE VALUES ALLOWED AT GIVEN L AND GIVEN PARITY 
-! THE N AND NP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO L-NN AND 
-! THE M AND MP VARIABLES TAKE ALL THE INTEGER VALUES FROM ZERO UP TO (NQMAX-L-NN)/2.
-! ONE ALSO HAS NQ=2*(M+N1+N2)+L+NN=2*(MP+N1P+N2P)+L+NN. 
-
-! THE SAID ABOVE DETERMINES THE MEANING OF THE N,M, AND NQMAX PARAMETERS OF
-! THE SUBROUTINE.
-
-! THE BRAC ARRAY IS OF THE FORM BRAC(NP,N1P,MP,N1,N2). AS SAID ABOVE, NP AND MP
-! RERESENT L1P AND L2P. THE MISSING N2P VALUE IS DETERMNED BY THE EQUALITY
-! MP+N1P+N2P=M+N1+N2 THAT IS ANOTHER FORM OF THE RELATION 
-! 2(N1+N2)+L1+L2=2(N1P+N2P)+L1P+L2P.
-
-! THE ROUTINE PARAMETERS CO AND SI ARE AS FOLLOWS, CO = COS(PHI) AND SI = SIN(PHI). 
-! THESE QUANTITIES DEFINE THE PSEUDO(!)ORTHOGONAL TRANSFORMATION 
-! XI1P = CO*XI1+SI*XI2, XI2P = SI*XI1-CO*XI2. 
-
-! (BRACKETS PERTAINING TO THE CASE OF THE ORTHOGONAL TRANSFORMATION 
-! XI1P = CO*XI1+SI*XI2, XI2P = -SI*XI1+CO*XI2  ARE SIMPLY EXPRESSED IN TERMS OF 
-! THOSE CALCULATED IN THE PRESENT PROGRAM, SEE THE ACCOMPANYING CPC PAPER.
-
-! THE MEANING OF THE FIRSTCALL PARAMETER IS EXPLAINED BELOW.
-
-! THE MEANING OF THE ARRAYS, OTHER THAN BRAC, ENTERING THE DIMENSION LIST IS AS
-! FOLLOWS. B IS A SUBSIDIARY ARRAY TO PERFORM THE RECURSION. 
-
-! FAC(I)=I!, DFAC(I)=(2I+1)!!, AND DEFAC(I)=(2I)!!, AND RFAC(I)=SQRT((2I)!)/I! 
-! IS A  SUBSIDIARY ARRAY TO CALCULATE THE FL ARRAY, SEE BELOW. THESE ARRAYS ARE
-! PRODUCED BY THE "ARR" ROUTINE WHICH IS CALLED FROM THE PRESENT ROUTINE AND
-! WHICH IS CONTAINED IN THE PRESENT FILE. THE SET UPPER BOUNDS OF THESE
-! ARRAYS ARE SUFFICIENT FOR THE COMPUTATION. THESE ARRAYS ARE USED IN THE PRESENT
-! ROUTINE AND IN THE "FLPHI" ROUTINE.
-
-! THE A ARRAY IS A(N,L)=(-1)**N/SQRT(DEFAC(N)*DFAC(N+L)). IT APPEARS BOTH 
-! IN THE N1=N2=0 BRACKETS AND IN THE RELATION BETWEEN THE < | > AND [ | ] TYPE
-! BRACKETS.
-
-! THE BI ARRAY IS BI(M,N)=FAC(N)/[FAC(M)*FAC(N-M)]. IT ENTERS THE 3J SYMBOLS. IT IS USED
-! IN THE "FLPHI" ROUTINE AND IN THE FUNCTION WIGMOD.
-
-! THE FL(NP,MP,N) ARRAY REPRESENTS THE QUANTITY 
-! [(2*L1P+1)*(2*L2P+1)]^{1/2}*F_L^VARPHI WHERE F_L^VARPHI IS GIVEN BY EQ. (23) IN THE
-! ACCOMPANYING CPC PAPER. THIS ARRAY IS PRODUCED IN ADVANCE BY THE
-! "FLPHI" ROUTINE WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS
-! CONTAINED IN THE PRESENT FILE. THE VARIABLES MP, NP, AND N ARE DEFINED ABOVE ALONG
-! WITH THEIR UPPER BOUNDS.
-
-! THE PSIP(P,Q) AND PSIM(P,Q) ARRAYS REPRESENT THE QUANTITIES (4) IN THE ACCOMPANYING
-! CPC PAPER. THEY ARE PRODUCED IN ADVANCE BY THE "COEFREC" ROUTINE
-! WHICH IS CALLED FROM THE PRESENT ROUTINE AND WHICH IS CONTAINED IN THE PRESENT
-! FILE. THEIR ARGUMENTS P AND Q TAKE THE VALUES L1P, L1P+/-1 AND L2P, L2P+/-1,
-! RESPECTIVELY. THE BOUNDS OF THE PSIP AND PSIM ARRAYS ARE SUCH THAT ALL THE P AND Q
-! VALUES REQUIRED TO PERFORM THE RECURSION ARE PROVIDED.
-
-! THE MENTIONED "FLPHI" ROUTINE USES THE FUNCTION WIGMOD WHICH IS ALSO CONTAINED
-! IN THE PRESENT FILE.
-
-! THE PARAMETER FIRSTCALL IS A LOGICAL VARIABLE WHICH ENSURES THAT FLPHI
-! IS CALLED ONLY ONCE AT A GIVEN L VALUE PROVIDING THAT ALL THE COMPUTATIONS ARE
-! PERFORMED WITH THE SAME NQMAX. UNNECESSARY REPEATED CALLS OF FLPHI ARE
-! SUPPRESSED VIA SAYING FIRSTCALL=.FALSE. IN A PROPER PLACE. IF THE CALLING PROGRAM
-! INCLUDES NESTED LOOPS OVER N-TYPE VARIABLES FROM N=0 UP TO N=NMAX AND THESE
-! LOOPS CONTAIN CALLS OF OSBRAC THEN FIRSTCALL=.FALSE. IS TO BE SAID AFTER THE MOST
-! INNER OF THESE LOOPS WITH THIS PURPOSE. AND AFTER A CHANGE OF L ONE SHOULD
-! CALL FLPHI ANEW VIA SAYING FIRSTCALL=.TRUE. (SEE ALSO THE APPENDED PROGRAM
-! TESTOSBRAC AS AN EXAMPLE.)
-!----------------------------------------------------------------------------------------------------------------------------------------
-    DOUBLE PRECISION BRAC,FL,PSIP,PSIM,A,FAC,DFAC,DEFAC,RFAC,BI,CO,SI,T,CO2,SI2,SC,FA,&
-    PLFACT,D2L12P,TN,PA,S,FA1,FA2 
-    LOGICAL FIRSTCALL
-    DIMENSION BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2),&
-! THIS ARRAY IS OF THE FORM: BRAC(NP,N1P,MP,N1,N2)    
-    A(0:42,0:84),PSIP(85,85),PSIM(85,85),&    
-    FL(0:84,0:42,0:84),&
-! THIS ARRAY IS OF THE FORM: FL(NP,MP,N)
-    FAC(0:169),DFAC(0:84),DEFAC(0:42),RFAC(0:84),BI(0:169,0:169)
-    SAVE PSIP,PSIM,FL,FAC,DFAC,RFAC,A,BI
-    IF (NQMAX.GT.84) THEN
-        WRITE(6,*)'IN THE R(*8) COMPUTATION NQMAX SHOULD NOT EXCEED 84' 
-        STOP
-    ENDIF 
-    IF (CO.EQ.0.D0.OR.SI.EQ.0.D0) THEN
-        WRITE(6,*)'THE PROGRAM IS OF NO USE AT ZERO COS(PHI) OR SIN(PHI) VALUES'
-        WRITE(6,*)'COS(PHI)=',CO,' SIN(PHI)=',SI
-        STOP
-    ENDIF
-    IF (NQMAX.LT.L) THEN
-        WRITE(6,*)'NQMAX=',NQMAX,' L=',L
-        WRITE(6,*)'L SHOULD NOT EXCEED NQMAX'
-        STOP
-    ENDIF
-! BRACKETS [N_1'L_1'N_2'L_2'|0L_10L_2]_L^\PHI TO START THE RECURSION, EQ. (22) IN  
-! THE ACCOMPANYING CPC PAPER. 
-! SUBSIDIARY QUANTITIES:
-    NQML=NQMAX-L
-    NQMLD2=NQML/2
-    NN=NQML-2*(NQMLD2)
-    NPMAX=L-NN      
-    IF (M.GT.NQMLD2.OR.N.GT.NPMAX) THEN
-        WRITE(6,*) 'M SHOULD NOT EXCEED (NQ-L)2 AND N SHOULD NOT EXCEED L-NN. BUT'
-        WRITE(6,*) 'M=',M,' MMAX=',MPMAX,' N=',N,' NMAX=',NPMAX
-        STOP
-    ENDIF
-    L1=M+N+NN 
-    L2=M-N+L
-    L1L2=L1+L2
-    L1L2ML=L1L2-L    
-    MPMAX=(L1L2ML-NN)/2
-    FA=SQRT((2.D0*L1+1)*(2*L2+1))*CO**L1L2
-    T=SI/CO
-        IF (FIRSTCALL) THEN
-        CALL ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
-        CALL COE(2*NQMAX+1,FAC,BI)
-!   THE A_{NL} ARRAY:
-    DO NI=0,NQMAX/2
-        NA=NI-2*(NI/2)
-        K=1
-        IF(NA.EQ.1)K=-1
-        DO LI=0,NQMAX-2*NI
-            A(NI,LI)=K/SQRT(DEFAC(NI)*DFAC(NI+LI))
-        ENDDO
-    ENDDO         
-    ENDIF 
-     PLFACT=SQRT(FAC(L1L2+L+1)*FAC(L1L2ML))     
-     IF (FIRSTCALL) CALL FLPHI(NQMAX,L,N,T,BI,RFAC,FL)
-! THIS SUBROUTINE PRODUCES THE ARRAY FL USED BELOW.
-! COMPUTATION OF THE N1=N2=0 BRACKETS. THEY ARE THE QUANTITIES BRAC(NP,N1P,MP,0,0): 
-    DO MP=0,MPMAX
-        L12P=2*MP+L+NN
-! L12P=L1P+L2P
-        D2L12P=1.D0/2.D0**L12P
-        N12P=MPMAX-MP
-! N12P=N1P+N2P=(L1L2-L12P)/2
-        DO N1P=0,N12P
-            N1PA=N1P-2*(N1P/2)
-            IF (N1PA.EQ.0) THEN
-                K=1
-            ELSE
-                K=-1
-            ENDIF
-! K=(-1)**N1P
-            TN=T**N12P*K
-            N2P=N12P-N1P
-            DO NP=0,NPMAX
-                L1P=MP+NP+NN
-                L2P=L12P-L1P  
-! BY DEFINITION ONE ALSO HAS: L2P=L2P(MP,NP)=MP-NP+L
-                PA=A(N1P,L1P)*A(N2P,L2P)
-                BRAC(NP,N1P,MP,0,0)=FL(NP,MP,N)*FA*TN*D2L12P*PLFACT*PA*PA
-            ENDDO
-        ENDDO
-    ENDDO
-! RECURSION TO OBTAIN THE GENERAL FORM BRACKETS 
-! [N_1'L_1'N_2'L_2'|N_1L_1N_2L_2]_L^\PHI. 
-    IF (FIRSTCALL) CALL COEFREC((NQMAX+L)/2+1,L,NN,NQMAX,PSIP,PSIM) 
-    N12MAX=(NQMAX-L1L2)/2
-    MPMAX0=MPMAX
-    CO2=CO**2 
-    SI2=SI**2
-    SC=SI*CO
-! THE N1=0, N2-1-->N2 RECURSION, EQ. (20)  IN THE ACCOMPANYING CPC PAPER WITH THE
-! MODIFICATION POINTED OUT THERE: 
-    DO N2=1,N12MAX
-        MPMAX=MPMAX+1
-! THIS MPMAX VALUE IS (NQ-L-NN)/2 WHERE NQ=L1+L2+2*N2.
-        DO MP=0,MPMAX
-            N12P=MPMAX-MP
-            DO N1P=0,N12P
-! MP=(L1P+L2P-L-NN)/2. THEREFORE, N12P=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.
-                DO NP=0,NPMAX
-                    L1P=MP+NP+NN
-                    L2P=MP-NP+L  
-                    S=0.D0 
-! BELOW THE RESTRICTIONS ARE IMPOSED ON THE BRAC ARRAY ENTERING THE RIGHT-HAND
-! SIDE OF THE RECURRENCE FORMULAE. THIS ARRAY IS OF THE FORM BRAC(K1,K2,K3,...) WHERE
-! K1=K1(NP), K2=N2(N1P), AND K3=K3(MP). NAMELY, K1=NP, OR NP+1, OR NP-1; 
-! K2=N1P OR N1P-1; K3=MP, OR MP-1, OR MP+1. (THE COMBINATION K2=N1P AND K3=MP+1
-! DOES NOT ARISE IN THE RECURRENCE FORMULAE.) THE QUANTITIES K1, K2, AND K3
-! REPRESENT, RESPECTIELY, NP, N1P, AND MP VALUES AT THE PRECEDING STAGE OF THE
-! RECURSION. THE RESTRICTIONS ENSURE THAT K1, K2, AND K3 RANGE WITHIN THE LIMITS
-! PERTAINING TO THE BRAC ARRAY OBTAINED AT THAT PRECEDING STAGE OF THE RECURSION. 
-
-! THUS K1, K2, AND K3 SHOULD BE NON-NEGATIVE. THEREFORE, WHEN K1=NP-1 THE VALUE
-! NP=0 IS TO BE EXCLUDED, WHEN K2=N1P-1 THE VALUE N1P=0 IS TO BE EXCLUDED, AND
-! WHEN K3=MP-1 THE VALUE MP=0 IS TO BE EXCLUDED. 
-
-! FURTHERMORE, IT SHOULD BE K1 \LE NMAX AND AT THE SAME TIME ONE HAS NP \LE NMAX.
-! THEREFORE, WHEN K1=NP+1 THE VALUE NP=NPMAX IS TO BE EXCLUDED.
-
-! IN ADDITION, IT SHOULD BE K2+K3 \LE MPMAX-1 WHILE ONE HAS N1P+MP \LE MPMAX. WHEN
-! K2=N1P AND K3=MP-1, OR K2=N1P-1 AND K3=MP-1, OR K2=N1P-1 AND K3=MP ONE
-! AUTOMATICALLY HAS N1P+MP \LE MPMAX. THEREFORE, IN THESE CASES THE CONDITION
-! K2+K3 \LE MPMAX-1 DOES NOT CREATE ANY RESTRICTIONS. 
-! BUT WHEN K2=N1P AND K3=MP, OR K2=N1P-1 AND K3=MP+1 THE CASE N1P+MP=MPMAX
-! IS TO BE FORBIDDEN. 
-
-! IN THE RECURRENCE FORMULAE BELOW THE DESCRIBED RESTRICTIONS ARE
-! IMPOSED.                            
-                    IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,0,N2-1)*PSIP(L1P,L2P)
-                    IF (N1P.NE.0.AND.N1P.NE.N12P) &
-                    S=S+BRAC(NP,N1P-1,MP+1,0,N2-1)*PSIP(L1P+1,L2P+1)
-! RECALL THAT N12P=MPMAX-MP.                     
-                    IF (NP.NE.NPMAX.AND.N1P.NE.0) &
-                    S=S-BRAC(NP+1,N1P-1,MP,0,N2-1)*PSIM(L1P+1,L2P)
-                    IF (NP.NE.0.AND.N1P.NE.N12P) &
-                    S=S-BRAC(NP-1,N1P,MP,0,N2-1)*PSIM(L1P,L2P+1)
-                    S=S*SC
-                    IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,0,N2-1)*SI2
-                    IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,0,N2-1)*CO2
-                    BRAC(NP,N1P,MP,0,N2)=S
-                ENDDO
-            ENDDO
-        ENDDO
-    ENDDO ! N2
-! N1-1-->N1 RECURSION, EQ. (20) IN THE ACCOMPANYING CPC PAPER:
-    DO N2=0,N12MAX
-        MPMAX=MPMAX0+N2
-! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2=MPMAX0+N2 SINCE NQ=L1+L2+2*N2.
-! THE RECURSION:        
-        DO N1=1,N12MAX-N2
-            MPMAX=MPMAX+1
-! THE CURRENT MPMAX VALUE IS (NQ-L-NN)/2 AND NQ=L1+L2+2*N2+2*N1.            
-            DO MP=0,MPMAX
-                N12P=MPMAX-MP
-                DO N1P=0,N12P
-! MP=(L1P+L2P-L-NN)/2. THEREFORE, N1PMAX=MPMAX-MP=(NQ-L1P-L2P)/2=N1P+N2P.              
-                    DO NP=0,NPMAX
-                        L1P=MP+NP+NN
-                        L2P=MP-NP+L  
-                        S=0.D0                 
-                        IF (MP.NE.0) S=BRAC(NP,N1P,MP-1,N1-1,N2)*PSIP(L1P,L2P)
-                        IF (N1P.NE.0.AND.N1P.NE.N12P) & 
-                        S=S+BRAC(NP,N1P-1,MP+1,N1-1,N2)*PSIP(L1P+1,L2P+1)
-                        IF (NP.NE.NPMAX.AND.N1P.NE.0) &
-                        S=S-BRAC(NP+1,N1P-1,MP,N1-1,N2)*PSIM(L1P+1,L2P)
-                        IF (NP.NE.0.AND.N1P.NE.N12P) &
-                        S=S-BRAC(NP-1,N1P,MP,N1-1,N2)*PSIM(L1P,L2P+1)
-! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
-                        S=-S*SC
-                        IF (N1P.NE.0) S=S+BRAC(NP,N1P-1,MP,N1-1,N2)*CO2
-                        IF (N1P.NE.N12P) S=S+BRAC(NP,N1P,MP,N1-1,N2)*SI2
-! THESE RELATIONS ARE THE SAME AS THE CORRESPONDING ONES ABOVE.
-                        BRAC(NP,N1P,MP,N1,N2)=S
-                    ENDDO
-                ENDDO
-            ENDDO
-        ENDDO ! N1
-    ENDDO ! N2  
-! RENORMALIZATION OF THE BRACKETS, EQ. (18) IN THE ACCOMPANYING CPC PAPER:
-    DO N2=0,N12MAX
-        MPMAX1=MPMAX0+N2
-        DO N1=0,N12MAX-N2
-            MPMAX=MPMAX1+N1
-            FA1=A(N1,L1)*A(N2,L2)
-            DO MP=0,MPMAX
-                N12P=MPMAX-MP
-                DO N1P=0,N12P
-                    N2P=N12P-N1P
-                    DO NP=0,NPMAX
-                        L1P=MP+NP+NN
-                        L2P=MP-NP+L 
-                        FA2=A(N1P,L1P)*A(N2P,L2P)
-                        BRAC(NP,N1P,MP,N1,N2)=BRAC(NP,N1P,MP,N1,N2) &
-                        *FA1/FA2                           
-                    ENDDO ! NP
-                ENDDO ! N1P
-            ENDDO ! MP
-        ENDDO ! N1
-    ENDDO ! N2
-    RETURN
-    END 
- 
-    SUBROUTINE ARR(FAC,DFAC,DEFAC,RFAC,NQMAX)
-! FAC(I), DFAC(I),DEFAC(I), AND RFAC(I) ARE, RESPECTIVELY, THE QUANTITIES
-! I!, (2I+1)!!, (2I)!!, AND SQRT((2*I)!)/I!
-    DOUBLE PRECISION FAC,DFAC,DEFAC,RFAC
-    DIMENSION FAC(0:169),DFAC(0:84),DEFAC(0:42),RFAC(0:84)
-    FAC(0)=1.D0
-    DFAC(0)=1.D0
-    DEFAC(0)=1.D0
-    RFAC(0)=1.D0
-    DO I=1,2*NQMAX+1 
-        FAC(I)=FAC(I-1)*I 
-    ENDDO            
-    DO I=1,NQMAX+1           
-        DFAC(I)=DFAC(I-1)*(2*I+1)
-    ENDDO    
-    DO I=1,NQMAX/2+1
-        DEFAC(I)=DEFAC(I-1)*2*I
-    ENDDO
-    DO I=1,NQMAX    
-        RFAC(I)=RFAC(I-1)*2*SQRT(1-0.5D0/I)    
-    ENDDO
-    RETURN
-    END
-     
-    SUBROUTINE FLPHI(NQMAX,L,N,T,BI,RFAC,FL)
-! PROVIDES THE QUANTITY SQRT((2L1P+1)*(2L2P+1))*F_L^\VARPHI, SEE EQ. (23) IN THE CPC
-! PAPER, IN THE FORM OF THE FL ARRAY.
-! USES THE FUNCTION WIGMOD. 
-! THIS FLPHI ROUTINE IS NOT IDENTICAL TO THAT IN THE OTHER FILE allosbrac.f90.
-!    IMPLICIT DOUBLE PRECISION(A-H,O-Z)
-    DOUBLE PRECISION BI,RFAC,FL,T,T2,SQP,F,WIGMOD,z
-    DIMENSION BI(0:169,0:169),RFAC(0:84),FL(0:84,0:42,0:84)
-! THE OUTPUT IS FL(NP,MP,N) WHERE NP=(L1P-L2P+L-NN)/2,
-! MP=(L1P+L2P-(L+NN))/2, L1P+L2P=L1T+L2T, AND N=(L1-L2+L-NN)/2, L1-L2=L1T-L2T. 
-    T2=T*T    
-    NQMAL=NQMAX-L
-    MPMAX=NQMAL/2
-    NN=NQMAL-2*MPMAX
-    LMNN=L-NN 
-    L1ML2=2*N-LMNN
-    LL3=2*L 
-       DO MP=0,MPMAX
-        L1PL2P=2*MP+L+NN
-           DO NP=0,LMNN
-            L1PML2P=2*NP-LMNN   
-            L1P=(L1PL2P+L1PML2P)/2
-            L2P=(L1PL2P-L1PML2P)/2
-            M3=L2P-L1P
-            SQP=SQRT((2*L1P+1.D0)*(2*L2P+1))
-            L1T=(L1PL2P+L1ML2)/2
-            L2T=(L1PL2P-L1ML2)/2
-            N3=L1T+L2T-L
-            IMIN=ABS(L1T-L1P)
-            IMAX=MIN(L1P+L1T,L2P+L2T)
-            NAL1=(L1T+L1P-IMIN)/2
-            NAL2=(L1T-L1P+IMIN)/2
-            NAL3=(L2T-L2P+IMIN)/2
-            NAL4=(L2T+L2P-IMIN)/2
-            NAL4P=NAL4-2*(NAL4/2)
-            F=0.D0
-            z=t**imin
-            if(nal4p.ne.0)z=-z
-            DO I=IMIN,IMAX,2
-                F=F+RFAC(NAL1)*RFAC(NAL2)*RFAC(NAL3)*RFAC(NAL4)*&
-                WIGMOD(L1T,L2T,L,L1P-I,I-L2P,BI)*z
-                NAL1=NAL1-1
-                NAL2=NAL2+1
-                NAL3=NAL3+1
-                NAL4=NAL4-1                             
-                z=-z*t2
-            ENDDO
-                FL(NP,MP,N)=F*SQRT((2*L1P+1)*(2*L2P+1)&
-                *BI(L1T+L-L2T,2*L1T)*BI(N3,2*L2T)/((LL3+1)*BI(N3,L1T+L2T+L+1)*BI(L+M3,LL3)))
-        ENDDO
-    ENDDO
-    RETURN
-    END
-
-   FUNCTION WIGMOD(L1,L2,L3,M1,M2,BI)
-! THE 3J SYMBOL IN TERMS OF BINOMIAL COEFFICIENTS WITHOUT THE FACTOR
-! SQRT(F) WHERE F IS AS FOLLOWS,
-! F=BI(L1+L3-L2,2*L1)*BI(N3,2*L2)/((LL3+1)*BI(N3,L1+L2+L3+1)*BI(L3+M3,LL3)) 
-! WITH N3=L1+L2-L3, LL3=2*L3, AND M3=-M1-M2. 
-! THIS WIGMOD ROUTINE IS NOT IDENTICAL TO THAT IN THE OTHER FILE allosbrac.f90.  
-    DOUBLE PRECISION WIGMOD,BI,S
-    DIMENSION BI(0:169,0:169)
-    M3=-M1-M2
-    N3=L1+L2-L3
-    LM1=L1-M1
-    LP2=L2+M2
-    KMIN=MAX(0,L1-L3+M2,L2-L3-M1)
-    KMAX=MIN(N3,LM1,LP2)
-    S=0.d0
-    NPH=1
-    DO K=KMIN,KMAX
-        S=S+NPH*BI(K,N3)*BI(LM1-K,L1+L3-L2)*BI(LP2-K,L2+L3-L1)
-        NPH=-NPH
-    ENDDO 
-    WIGMOD=S/SQRT(BI(LM1,2*L1)*BI(LP2,2*L2))
-    NY=KMIN+L1-L2-M3
-    NYP=NY-2*(NY/2)
-    IF(NYP.NE.0)WIGMOD=-WIGMOD
-    RETURN
-    END    
-
-    SUBROUTINE COE(NMAX,FAC,BI)
-!    IMPLICIT DOUBLE PRECISION(A-H,O-Z)
-    DOUBLE PRECISION FAC,BI
-    DIMENSION FAC(0:169),BI(0:169,0:169)
-    DO N=0,NMAX
-        DO M=0,N/2
-            BI(M,N)=FAC(N)/(FAC(M)*FAC(N-M))
-            BI(N-M,N)=BI(M,N)
-        ENDDO
-    ENDDO  
-    RETURN
-    END
-    
-    SUBROUTINE COEFREC(MAXCOE,L,NN,NQMAX,PSIP,PSIM)
-! CALCULATES THE COEFFICIENTS OF THE RECURSION FORMULA, EQ. (21) IN THE ACCOMPANYING
-! CPC PAPER
-    DOUBLE PRECISION PSIP,PSIM
-    DIMENSION PSIP(85,85),PSIM(85,85)
-! MAXCOE EQUALS INT((NQMAX+LMAX)/2)+1 AT CALLS OF THIS ROUTINE.
-    LP=L+1
-    LM=L-1
-    DO M1=1,MAXCOE
-        DO M2=1,M1
-            M1P2=M1+M2
-            NA=M1P2+L
-            NNA=NA-2*(NA/2)
-            M1M2=M1-M2
-! CALCULATION OF PSIP. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P AND L2P OR 
-! L1P+1 AND L2P+1.)
-            IF (NNA.EQ.NN.AND.M1P2.GE.L.AND.M1P2.LE.NQMAX.AND.&
-            ABS(M1M2).LE.L) THEN
-                IF (M1P2.GT.LP) THEN
-                    M1P2L=M1P2+L
-                    M1P2ML=M1P2-L
-                    PSIP(M1,M2)=SQRT(M1P2L*(M1P2L+1.D0)*M1P2ML*(M1P2ML-1)/((4*M1*M1-1)&
-                    *(4*M2*M2-1)))
-                ELSE
-                    PSIP(M1,M2)=0.D0
-                ENDIF
-                PSIP(M2,M1)=PSIP(M1,M2) 
-            ENDIF
-! CALCULATION OF PSIM. (IN THIS CASE M1 AND M2 REPRESENT, RESPECTIVELY, L1P+1 AND L2P
-! OR L2P+1 AND L1P.)
-            IF (NNA.NE.NN.AND.M1P2.GE.LP.AND.M1P2.LE.NQMAX&
-            .AND.M1M2.LE.LP.AND.M1M2.GE.-LM) THEN
-                IF (M1M2.LT.L) THEN
-                    M1M2L=M1M2+L
-                    M1M2ML=M1M2-L 
-                    PSIM(M1,M2)=SQRT(M1M2L*(M1M2L+1.D0)*M1M2ML*(M1M2ML-1)/&
-                    ((4*M1*M1-1)*(4*M2*M2-1)))
-                ELSE
-                    PSIM(M1,M2)=0.D0
-                ENDIF
-                PSIM(M2,M1)=PSIM(M1,M2) 
-            ENDIF                                               
-        ENDDO
-    ENDDO 
-    RETURN
-    END             
-

OSBRACKETS/testallosbrac.f90

@@ -1,260 +0,0 @@
-    program testallosbrac
-! Tests of the program that calculates the <N_1'L_1'N_2'L_2'|N_1L_1N_2L_2>_L^\varphi brackets.
-    implicit double precision(a-h,o-z)
-    dimension fac(0:201),dfac(0:201),defac(0:201),rfac(0:100)
-! These arrays enter the routine ARR called below. The first three arrays are used merely to
-! compute the exact expression for the bracket below. The array rfac is of no use.    
-    allocatable :: brac(:,:,:,:,:,:,:,:),bi(:,:)
-! Below the allocation/deallocation actions are done independently for each of the tests
-! performed as it would be required in the case when only a single test is being performed. 
-    open(13,file='allosoutput')
-! The input angular parameters are CO=COS(VARPHI) and SI=SIN(VARPHI). Set, for example,
-    si=1/sqrt(2.d0)
-    co=-si
-! This corresponds to Eq. (13) in the accompanying CPC paper in the equal mass case. The
-! corresponding transformation (7) in that paper arises at calculating standard two-body 
-! shell-model matrix elements.
-! One might also test the program for COS(VARPHI) and SIN(VARPHI) values pertaining to
-! nonequal mass cases. For example, one might set
-!    si=1/sqrt(3.d0)
-!    co=-sqrt(2/3.d0)
-! 1. Comparison of the value of the <n1p,L00|n1,l1,n2,l2>_L^\varphi bracket with the result 
-! of its analytic calculation. Here l1+l2 should be of the same parity as L.
-    l=4
-    l1=2
-    l2=4
-    n1=0
-    n2=0
-    nq=2*n1+2*n2+l1+l2
-! also nq=2*n1p+L
-    NN=0
-! this is because nq-l=2*n1p is even
-    mp=0
-! here l1p=l, mp=(l1p+l2p-l-nn)/2=0    
-    np=l
-! np=(l1p-l2p+l-nn)/2=l    
-    n1p=(nq-l)/2
-! Computation of the exact bracket.
-    nfac=2*nq+1
-    CALL ARR(FAC,DFAC,DEFAC,RFAC,NQ)
-! The array of binomial coefficients to calculate the Wigner coefficients:                
-    allocate (BI(0:NFAC,0:NFAC),STAT=istatus)
-    IF (istatus /= 0) STOP "TESTOSBRAC: allocate BI failed"
-    call coe(NFAC,FAC,BI)
-    nm=n1p-n1-n2
-    nfa=nm-2*(nm/2)
-    LL3=2*L
-    N3=L1+L2-L
-! exact bracket:        
-    b=co**(2*n1+l1)*si**(2*n2+l2)*(-1)**l*sqrt((2.d0*l1+1)*(2*l2+1))*(-1)**nfa*&
-    WIGMOD(l1,l2,l,0,0,bi,nfac)*SQRT(BI(L1+L-L2,2*L1)*BI(N3,2*L2)/&
-    ((LL3+1)*BI(N3,L1+L2+L+1)*BI(L,LL3)))&
-    *sqrt(defac(n1p)*dfac(n1p+l)/(defac(n1)*defac(n2)*dfac(n1+l1)*dfac(n2+l2)))
-! bracket from the code:                
-    ndel=4
-    nqmax=nq+ndel        
-! The result should not depend on ndel, provided that ndel is even, i.e., nq and nqmax are of the
-! same parity.      
-    NQML=NQMAX-L
-    m=(l1+l2-l-nn)/2
-    n=(l1-l2+l-nn)/2
-    allocate (BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,&
-    0:(NQMAX-L)/2,0:L,0:(NQMAX-L)/2,L:L),STAT=istatus)
-    IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
-    CALL ALLOSBRAC(NQMAX,L,L,CO,SI,BRAC) 
-    a=brac(np,n1p,mp,n1,n2,n,m,l)
-    write(6,*)'1. bracket <n1p L 0 0|n1 l1 n2 l2>_L^varphi' 
-    write(6,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L   =',l
-    write(6,*)'exact bracket',b,' bracket from the code',a
-    write(13,*)'1. bracket <n1p L 0 0|n1 l1 n2 l2>_L^varphi'
-    write(13,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l 
-    write(13,*)'exact bracket',b,' bracket from the code',a
-!    write(6,*)co,si
-    deallocate (BRAC,STAT=istatus)
-! 2. Test of the symmetry of the brackets. The relation 
-! <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\varphi=<N1,L1,N2,L2|N1P,L1P,N2P,L2P>_L^\varphi
-! should hold true. This symmetry relation may be rewritten as
-! BRAC(NP,N1P,MP,N1,N2,N,M,L)=BRAC1(N,N1,M,N1P,N2P,NP,MP,L) where N and
-! M represent the angular momenta L1 and L2 and NP and MP represent L1P and L2P. 
-! An example:
-    nq=30
-    l=9
-! nn=1 in the present example.
-! calculate: nmax=l-nn.
-! nmax=8  in the present example. Any n and np that do not exceed nmax may be employed.
-    n=4 
-    np=3      
-! calculate: nqmld2=(nq-l-nn)/2.
-! nqmld2=10 in the present example. Any m, n1, n2, mp, n1p, and n2p may be employed
-! provided that m+n1+n2=mp+n1p+n2p=nqmld2.
-    m=5
-    n1=3
-    n2=2
-    mp=4
-    n1p=3
-    n2p=3
-    nqmax=32
-! Recall that nq=30. Any nqmax of the same parity as nq and such that nqmax.ge.nq  
-! may be employed. Then its choice does not influence the results.  
-    allocate (BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,&
-    0:(NQMAX-L)/2,0:L,0:(NQMAX-L)/2,L:L),STAT=istatus)
-    IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed" 
-    CALL ALLOSBRAC(NQMAX,L,L,CO,SI,BRAC)  
-    write(6,*)'2. Symmetry test, L=',l
-    write(6,*)'m =',m,' n =',n,' n1 =',n1,' n2 =',n2
-    write(6,*)'mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
-    write(6,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =',BRAC(np,n1p,mp,n1,n2,n,m,l)
-    write(6,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =',BRAC(n,n1,m,n1p,n2p,np,mp,l)
-    write(13,*)'2. Symmetry test, L=',l
-    write(13,*)'m =',m,' n =',n,' n1 =',n1,' n2 =',n2
-    write(13,*)'mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
-    write(13,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =',BRAC(np,n1p,mp,n1,n2,n,m,l)
-    write(13,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =',BRAC(n,n1,m,n1p,n2p,np,mp,l)  
-    deallocate (BRAC,STAT=istatus)
-! 3. Test of the relation 
-! \sum_{i,i',L}\sum_j<j|i><j|i'>=\sum_L N_0(L).
-! Here i=(n_1,l_1,n_2,l_2), i'=(n_3,l_3,n_4,l_4), and j=(n_1p,l_1p,n_2p,l_2p). The l.h.s. sum 
-! runs over all i and i' values such that nq(i)=nq(i')=nq. 
-! In the r.h.s. sum, N_0(L) is the number of the (n_1,l_1,n_2,l_2) states pertaining to
-! given nq and L values. The sum runs over all the L values compatible with a given nq value
-! and thus represents the total number of states with this nq.
-    nq=14
-    ndel=4
-    nqmax=nq+ndel    
-! The result should not depend on ndel, i.e., on nqmax, provided that ndel is even.         
-    numb=0
-    s=0.d0
-    LMIN=0
-    LMAX=NQ   
-    allocate (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),STAT=istatus)
-    IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
-    CALL ALLOSBRAC(NQMAX,LMIN,LMAX,CO,SI,BRAC)                   
-    do L=0,nq
-        nqml=nq-l
-        nqmld2=nqml/2  
-        nn=nqml-2*nqmld2          
-        mmax=(nq-(l+nn))/2
-        nmax=l-nn
-        nst=(nmax+1)*(mmax+1)*(mmax+2)/2
-        do na=0,nmax
-            do ma=0,mmax
-! IN THE ABOVE NOTATION NA=NA(L1,L2), MA=MA(L1,L2).
-                n12=mmax-ma               
-                do nb=0,nmax
-                    do mb=0,mmax
-! IN THE ABOVE NOTATION NB=NB(L3,L4), MB=MB(L4,L4).
-                        n34=mmax-mb               
-                        ds=0.d0
-                        do n2=0,n12
-                            n1=n12-n2
-                            do n4=0,n34
-                                n3=n34-n4
-                                DO MP=0,mmax
-                                    do n1p=0,mmax-mp
-                                        DO NP=0,nmax
-                                            DS=DS+BRAC(np,n1p,mp,n1,n2,na,ma,l)&
-                                            *BRAC(np,n1p,mp,n3,n4,nb,mb,l)
-                                        enddo ! NP
-                                    ENDDO ! n1p
-                                ENDDO ! MP
-                            enddo ! n4
-                        enddo ! n2
-                        s=s+ds
-                    enddo ! MB
-                enddo ! NB
-            enddo ! MA
-        enddo    ! NA
-        numb=numb+nst 
-    ENDDO !l 
-    deallocate (BRAC,STAT=istatus)   
-    write(6,*)'3. Test of the relation sum_{i,ip,L}sum_j<j|i><j|ip> = the number of states.'
-    write(6,*)'See the comments in the text of the program and in more detail in the' 
-    write(6,*)'accompanying CPC paper for further explanations.'  
-    write(6,*)'nqmax =',nqmax,' nq =',nq
-    write(6,*)'sum_{i1,i2,L}sum_j<j|i1><j|i2> =',s,' exact value (number of states) =', numb
-    write(13,*)'3. Test of the relation sum_{i,ip,L}sum_j<j|i><j|ip> = the number of states.'
-    write(13,*)'See the comments in the text of the program and in more detail in the' 
-    write(13,*)'accompanying CPC paper for further explanations.'  
-    write(13,*)'nqmax =',nqmax,' nq =',nq
-   write(13,*)'sum_{i1,i2,L}sum_j<j|i1><j|i2> =',s,' exact value (number of states)=',numb         
-! 4. Test of the relation \sum_{i,i',L}|\sum_j<j|i><j|i'>-\delta_{i,i'}|=0.
-! Here i=(n_1,l_1,n_2,l_2), i'=(n_1',l_1',n_2',l_2'), and j=(n_1'',l_1'',n_2'',l_2''). The outer
-! sum runs over all i and i' values such that nq(i)=nq(i') \le nqmax.  
-! Also one has nq(i)=nq(i')=nq(j). 
-    nqmax=12
-    lmin=0
-    lmax=nqmax
-    nqml=nqmax-lmin
-    nqmld2=nqml/2        
-    allocate (BRAC(0:LMAX,0:nqmld2,0:nqmld2,0:nqmld2,0:nqmld2,0:LMAX,& 
-    0:nqmld2,LMIN:LMAX),STAT=istatus)
-    IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
-    CALL ALLOSBRAC(NQMAX,LMIN,LMAX,CO,SI,BRAC) 
-! When nqmax is sufficiently high, this calculation is to be performed in somewhat another way
-! because of the memory restrictions. In such cases one may use ALLOSBRAC at lmin=lmax=l.
-! The BRAC array is to be allocated and deallocated inside the loop over l below, and
-! ALLOSBRAC is to be called inside this loop.   
-    s=0.d0
-! The quantity s represents the expression \sum_{i,i',L}|\sum_j<j|i><j|i'>-\delta_{i,i'}| we are
-! calculating.           
-    do l=0,nqmax
-    nqml=nqmax-l
-    nqmld2=nqml/2
-    nn=nqml-2*nqmld2
-    nmax=l-nn
-    lpnn=l+nn
-    do nq=lpnn,nqmax,2
-        mmax=(nq-lpnn)/2                       
-            do ma=0,mmax
-! The equality mmax=(nq-lpnn)/2 follows from the definition of m-type variables.             
-                n12=mmax-ma
-! The quantity n12 designates n1+n2 and here it is taken into account that
-! n1+n2+ma=mmax.                
-                do na=0,nmax
-! In the case of odd nqmax and l equal to zero, one has nmax=-1 and no contribution arises as it
-! should be. 
-                    do n1=0,n12                
-                        n2=n12-n1                           
-                        do mb=0,mmax
-                            n34=mmax-mb
-                            do nb=0,nmax
-                                do n3=0,n34
-                                    n4=n34-n3
-                                    ds=0.d0 
-! The quantity ds represents the contribution of \sum_j|<j|i>_L<j|i'>_L-\delta_{i,i'}| to the net
-! result.                                     
-                                    DO MP=0,MMAX
-                                        do n1p=0,mmax-mp
-                                            DO NP=0,NMAX
-                                            DS=DS+BRAC(np,n1p,mp,n1,n2,na,ma,l)&
-                                            *BRAC(np,n1p,mp,n3,n4,nb,mb,l)
-                                            enddo ! NP
-                                        ENDDO ! n1p
-                                    ENDDO ! MP
-                                    if(ma.eq.mb.and.na.eq.nb.and.n1.eq.n3.and.n2.eq.n4)then
-                                        ds=abs(ds-1.d0)
-                                    else
-                                        ds=abs(ds)
-                                    endif
-                                    s=s+ds
-                                enddo ! n3
-                            enddo ! nb
-                        enddo ! mb
-                    enddo ! n1
-                enddo ! na
-            enddo ! ma
-        enddo ! nq
-    ENDDO ! l
-    deallocate (BRAC,STAT=istatus)
-    write(6,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
-    write(6,*)'For further explanations see the comments in the text of the program'
-    write(6,*)'and, in more detail, in the accompanying CPC paper.'      
-    write(6,*)'nqmax =',nqmax
-    write(6,*)'sum_{i1,i2,L}|sum_j<j|i1><j|i2>-delta_{i1,i2}| =',s
-    write(13,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
-    write(13,*)'For further explanations see the comments in the text of the program'
-    write(13,*)'and, in more detail, in the accompanying CPC paper.'       
-    write(13,*)'nqmax =',nqmax
-    write(13,*)'sum_{i1,i2,L}|sum_j<j|i1><j|i2>-delta_{i1,i2}| =',s
-    end

OSBRACKETS/testosbrac.f90

@@ -1,266 +0,0 @@
-    program testosbrac
-! Tests of the program that calculates the <N_1'L_1'N_2'L_2'|N_1L_1N_2L_2>_L^\varphi brackets.
-    implicit double precision(a-h,o-z)
-    logical firstcall
-    dimension fac(0:169),dfac(0:84),defac(0:42),rfac(0:84),bi(0:169,0:169)
-! These arrays enter the routine ARR called below. The first three arrays are used merely to
-! compute the exact expression for the bracket below. The array rfac is of no use. 
-    allocatable :: brac(:,:,:,:,:),brac1(:,:,:,:,:)
-    open(13,file='osoutput')
-! The input angular parameters are CO=COS(VARPHI) and SI=SIN(VARPHI). Set, for example,
-    si=1/sqrt(2.d0)
-    co=-si
-! This corresponds to Eq. (13) in the accompanying CPC paper in the equal mass case. The
-! corresponding transformation (7) in that paper arises at calculating standard two-body 
-! shell-model matrix elements.
-! One might also test the program for COS(VARPHI) and SIN(VARPHI) values pertaining to
-! nonequal mass cases. For example, one might set
-!   si=1/sqrt(3.d0)
-!   co=-sqrt(2/3.d0)
-! 1. Comparison of the value of the <n1p,L00|n1,l1,n2,l2>_L^\varphi bracket with the result 
-! of its analytic calculation. Here l1+l2 should be of the same parity as L.
-    l=4
-    l1=2
-    l2=4
-    n1=3
-    n2=2
-    mp=0
-    np=l
-    nq=2*n1+2*n2+l1+l2 
-    n1p=(nq-l)/2
-! Computation of the exact bracket.
-    nfac=2*nq+1 
-    CALL ARR(FAC,DFAC,DEFAC,RFAC,NQ)
-! The array of binomial coefficients to calculate the Wigner coefficients:                
-    call coe(NFAC,FAC,BI)      
-    nm=n1p-n1-n2
-    nfa=nm-2*(nm/2)
-    LL3=2*L
-    N3=L1+L2-L
-! exact bracket:        
-    b=co**(2*n1+l1)*si**(2*n2+l2)*(-1)**l*sqrt((2.d0*l1+1)*(2*l2+1))*(-1)**nfa*&
-    WIGMOD(l1,l2,l,0,0,bi)*SQRT(BI(L1+L-L2,2*L1)*BI(N3,2*L2)/&
-    ((LL3+1)*BI(N3,L1+L2+L+1)*BI(L,LL3)))&
-    *sqrt(defac(n1p)*dfac(n1p+l)/(defac(n1)*defac(n2)*dfac(n1+l1)*dfac(n2+l2)))
-! bracket from the code:                
-    ndel=4
-    nqmax=nq+ndel        
-! The result should not depend on ndel, provided that ndel is even, i.e., nq and nqmax are of the
-! same parity.        
-    NQML=NQMAX-L
-    NN=0
-    m=(l1+l2-l-nn)/2
-    n=(l1-l2+l-nn)/2
-    allocate (BRAC(0:L,0:NQML/2,0:NQML/2,0:NQML/2,0:NQML/2),STAT=istatus)
-    IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
-    FIRSTCALL=.TRUE.
-    CALL OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
-    a=brac(np,n1p,mp,n1,n2)
-    write(6,*)'1. bracket <n1p,L00|n1l1n2l2>_L^phi'
-    write(6,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l 
-    write(13,*)'1. bracket <n1p,L00|n1l1n2l2>_L^phi'
-    write(13,*)'n1 =',n1,' l1 =',l1,' n2 =',n2,' l2 =',l2,' L =',l 
-    write(6,*)'exact bracket',b,' bracket from the code',a
-    write(13,*)'exact bracket',b,' bracket from the code',a
-    deallocate (BRAC,STAT=istatus)
-! 2. Test of the symmetry of the brackets. The relation 
-! <N1P,L1P,N2P,L2P|N1,L1,N2,L2>_L^\varphi=<N1,L1,N2,L2|N1P,L1P,N2P,L2P>_L^\varphi
-! should hold true. Let the bracket in the left-hand side of this relation is contained in the BRAC
-! array and that in its right-hand side in the BRAC1 array. Then the above symmetry relation may
-! be rewritten as BRAC(NP,N1P,MP,N1,N2)=BRAC1(N,N1,M,N1P,N2P) where N and M represent
-! the angular momenta L1 and L2 and NP and MP represent L1P and L2P.   
-! An example:
-    nq=30
-    l=9
-! nn=1 at this choice.
-! calculate: nmax=l-nn.
-! nmax=8  in the present example. Any n and np values 
-! that do not exceed nmax may be employed.
-    n=4 
-    np=3      
-! calculate: nqmld2=(nq-l-nn)/2.
-! nqmld2=10 in the present example. Any m, n1, n2, mp, n1p, and n2p may be employed
-! provided that m+n1+n2=mp+n1p+n2p=nqmld2.
-    m=5
-    n1=3
-    n2=2
-    mp=4
-    n1p=3
-    n2p=3
-    nqmax=36
-! Recall that nq=30. Any nqmax of the same parity as nq and such that nqmax.ge.nq  
-! may be employed. Then its choice does not influence the results.  
-    nqml=nqmax-l 
-    allocate (BRAC(0:L,0:NQML/2,0:NQML/2,0:NQML/2,0:NQML/2),STAT=istatus)
-    IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
-    CALL OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
-    allocate (BRAC1(0:L,0:NQML/2,0:NQML/2,0:NQML/2,0:NQML/2),STAT=istatus)
-    IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC1) failed"
-    CALL OSBRAC(NP,MP,L,NQMAX,CO,SI,FIRSTCALL,BRAC1)
-    write(6,*)'2. Symmetry test, L =',l
-    write(6,*)'m =',m,' n =',n,' n1=',n1,' n2=',n2,' mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
-    write(6,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =', BRAC(np,n1p,mp,n1,n2)
-    write(6,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =', BRAC1(n,n1,m,n1p,n2p)
-    write(13,*)'2. Symmetry test, L=',l
-    write(13,*)'m =',m,' n =',n,' n1=',n1,' n2=',n2,' mp =',mp,' np =',np,' n1p =',n1p,' n2p =',n2p
-    write(13,*)'bracket(mp,np,n1p,n2p;m,n,n1,n2) =', BRAC(np,n1p,mp,n1,n2)
-    write(13,*)'bracket(m,n,n1,n2;mp,np,n1p,n2p) =', BRAC1(n,n1,m,n1p,n2p)  
-    deallocate (BRAC1,STAT=istatus)
-    deallocate (BRAC,STAT=istatus)         
-! 3. Test of the relation 
-! \sum_{i,i',L}\sum_j<j|i>_L<j|i'>_L=\sum_L N_0(L).
-! Here i=(n_1,l_1,n_2,l_2), i'=(n_3,l_3,n_4,l_4), and j=(n_1p,l_1p,n_2p,l_2p). The l.h.s. sum 
-! runs over all i and i' values such that nq(i)=nq(i')=nq. 
-! In the r.h.s. sum, N_0(L) is the number of the (n_1,l_1,n_2,l_2) states pertaining to
-! given nq and L values. The sum runs over all the L values compatible with a given nq value
-! and thus represents the total number of states with this nq.
-    nq=14
-    ndel=6
-    nqmax=nq+ndel
-! The result should not depend on ndel, i.e., on nqmax, provided that ndel is even.         
-    numb=0
-    s=0.d0                   
-    do l=0,nq
-        FIRSTCALL=.true.
-        nqml=nqmax-l
-        nqmld2=nqml/2    
-        nn=nqml-2*nqmld2
-        mmax=(nq-(l+nn))/2
-        nmax=l-nn
-        nst=(nmax+1)*(mmax+1)*(mmax+2)/2
-        allocate (BRAC(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
-        IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
-        allocate (BRAC1(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
-        IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC1) failed"
-! The allocation/deallocation of BRAC and BRAC1, at nqmld2-->nqmax/2, might also be done
-! outside the loop.        
-        do ma=0,mmax
-            n12=mmax-ma 
-            do na=0,nmax 
-! IN THE ABOVE NOTATION NA=NA(L1,L2) AND MA=MA(L1,L2).                      
-                CALL OSBRAC(NA,MA,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
-! The order of summations below aims to avoid unnecessary calls of osbrac.                 
-                 do mb=0,mmax
-                    n34=mmax-mb 
-                    do nb=0,nmax 
-! IN THE ABOVE NOTATION NB=NB(L3,L4) AND MB=MB(L3,L4).                                  
-                        CALL OSBRAC(NB,MB,L,NQMAX,CO,SI,FIRSTCALL,BRAC1)
-                        ds=0.d0
-                        do n2=0,n12
-                            n1=n12-n2
-                            do n4=0,n34
-                                n3=n34-n4
-                                DO MP=0,mmax
-                                    do n1p=0,mmax-mp
-                                        DO NP=0,nmax
-                                            DS=DS+BRAC(NP,n1p,MP,n1,n2)*BRAC1(NP,n1p,MP,n3,n4)
-                                        enddo ! NP
-                                    ENDDO ! n1p
-                                ENDDO ! MP
-                            enddo ! n4
-                        enddo ! n2
-                        s=s+ds
-                    enddo ! nb
-                    FIRSTCALL=.false.
-! If the nb loop is nested inside the mb loop the FIRSTCALL=.false. command is more efficient.
-! In practice, usually the gain in the running time is relatively not large.                     
-                enddo ! mb
-            enddo ! ma
-        enddo    ! na
-        numb=numb+nst
-        deallocate (BRAC,STAT=istatus)
-        deallocate (BRAC1,STAT=istatus)
-    ENDDO ! l
-    write(6,*)'3. Test of the relation sum_{i,ip,L,j}<j|i>_L<j|ip>_L = the number of states.'
-    write(6,*)'See the comments in the text of the program and in more detail in the' 
-    write(6,*)'accompanying CPC paper for further explanations.'  
-    write(6,*)'nqmax =',nqmax,' nq =',nq
-    write(6,*)'sum_{i,ip,L,j}<j|i>_L<j|ip>_L  =',s,' exact value (number of states) =',numb
-    write(13,*)'3. Test of the relation sum_{i,ip,L}sum_j<j|i><j|ip> = the number of states.'
-    write(13,*)'See the comments in the text of the program and in more detail in the' 
-    write(13,*)'accompanying CPC paper for further explanations.'  
-    write(13,*)'nqmax =',nqmax,' nq =',nq         
-    write(13,*)'sum_{i,ip,L,j}<j|i>_L<j|ip>_L  =',s,' exact value (number of states) =',numb                      
-! 4. Test of the relation \sum_{i,i',L}\sum_j|<j|i>_L<j|i'>_L-\delta_{i,i'}|=0.
-! Here i=(n_1,l_1,n_2,l_2), i'=(n_1',l_1',n_2',l_2'), and j=(n_1'',l_1'',n_2'',l_2''). The outer
-! sum runs over all i and i' values such that nq(i)=nq(i') \le nqmax.  
-! Also one has nq(i)=nq(i')=nq(j).        
-    nqmax=12
-    s=0.d0 
-! The quantity s represents the expression \sum_{i,i',L}|\sum_j<j|i><j|i'>-\delta_{i,i'}| we are
-! calculating.                
-    do l=0,nqmax
-        FIRSTCALL=.true.   
-        nqml=nqmax-l
-        nqmld2=nqml/2    
-        nn=nqml-2*nqmld2
-        nmax=l-nn
-        lpnn=l+nn
-        allocate (BRAC(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
-        IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC) failed"
-        allocate (BRAC1(0:L,0:NQMLD2,0:NQMLD2,0:NQMLD2,0:NQMLD2),STAT=istatus)
-        IF (istatus /= 0) STOP "TESTOSBRAC: allocate(BRAC1) failed"
-! The allocation/deallocation of BRAC and BRAC1, at nqmld2-->nqmax/2, might also be done
-! outside the loop.          
-        do ma=0,nqmld2 
-! The above upper limit of summation corresponds to the definition of m-type
-! variables.                          
-            do na=0,nmax            
-! In the case of odd nqmax and l equal to zero, one has nmax=-1 and no contribution arises as it
-! should be.     
-                CALL OSBRAC(NA,MA,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
-! The order of summations below aims to avoid unnecessary calls of osbrac.                 
-                do mb=0,nqmld2
-                    do nb=0,nmax                    
-                        CALL OSBRAC(NB,MB,L,NQMAX,CO,SI,FIRSTCALL,BRAC1)                        
-                        do nq=l+nn,nqmax,2
-                            nqcmld2=(nq-l-nn)/2
-                            n12=nqcmld2-ma
-                            n34=nqcmld2-mb
-! We consider nq to be the number of quanta pertaining to states entering the brackets in our
-! sum. Then one has: ma+n1+n2=mb+n3+n4=mp+n1p+n2p=(nq-l-nn)/2 that is
-! nqcmld2.                           
-                            do n1=0,n12
-                                n2=n12-n1
-                                do n3=0,n34
-                                    n4=n34-n3                                     
-                                    ds=0.d0 
-! The quantity ds represents the contribution of \sum_j|<j|i>_L<j|i'>_L-\delta_{i,i'}| to the net
-! result.                                              
-                                    DO MP=0,nqcmld2
-                                        do n1p=0,nqcmld2-mp
-                                            DO NP=0,NMAX
-                                                DS=DS+BRAC(NP,n1p,MP,n1,n2)*BRAC1(NP,n1p,MP,n3,n4)
-                                            enddo ! NP
-                                        ENDDO ! n1p
-                                    ENDDO ! MP
-                                    if(ma.eq.mb.and.na.eq.nb.and.n1.eq.n3.and.n2.eq.n4)then
-                                        ds=abs(ds-1.d0)
-                                    else
-                                        ds=abs(ds)
-                                    endif
-                                    s=s+ds
-                                enddo ! n3
-                            enddo ! n1
-                        enddo ! nq
-                    enddo ! nb
-                    FIRSTCALL=.false.
-! If the nb loop is nested inside the mb loop the FIRSTCALL=.false. command is more efficient.
-! In practice, usually the gain in the running time is relatively not large. 
-                enddo ! mb
-            enddo ! ma
-        enddo ! na
-        deallocate (BRAC,STAT=istatus)
-        deallocate (BRAC1,STAT=istatus)
-    ENDDO ! l
-    write(6,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
-    write(6,*)'For further explanations see the comments in the text of the program'
-    write(6,*)'and, in more detail, in the accompanying CPC paper.'  
-    write(6,*)'nqmax =',nqmax
-    write(6,*)'sum_{i1,i2,L}|sum_j<j|i1>_L<j|i2>_L-delta_{i1,i2}| =',s
-    write(13,*)'4. Test of the relation sum_{i,ip,L}sum_j|<j|i>_L<j|ip>_L-delta_{i,ip}| = 0.'
-    write(13,*)'For further explanations see the comments in the text of the program'
-    write(13,*)'and, in more detail, in the accompanying CPC paper.'  
-    write(13,*)'nqmax =',nqmax
-    write(13,*)'sum_{i1,i2,L}|sum_j<j|i1>_L<j|i2>_L-delta_{i1,i2}| =',s
-    end

Project.toml

@@ -0,0 +1,15 @@
+[deps]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"

XZ_technique.jl

@@ -1,55 +0,0 @@
-using LinearAlgebra, Plots
-include("ho_basis.jl")
-include("p_space.jl")
-
-ω = 0.5 * exp(-1im * 0.47 * pi)
-μ = 0.5
-l = 0
-V1 = -5
-R1 = sqrt(3)
-V2 = 2
-R2 = sqrt(10)
-n_max = 15
-
-ns = collect(0:n_max)
-ls = fill(l, n_max + 1)
-
-T = sp_T_matrix(ns, ls; ω=ω, μ=μ)
-V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); ω=ω) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); ω=ω)
-
-n_EC = 8
-train_cs = (0.7 .+ 0.05 * randn(n_EC)) - 1im * (0.2 .+ 0.05 * randn(n_EC))
-target_cs = range(0.77, 0.22, 6)
-
-train_E = zeros(ComplexF64, n_EC)
-EC_basis = zeros(ComplexF64, (n_max + 1, length(train_cs)))
-exact_E = zeros(ComplexF64, length(target_cs))
-extrapolate_E = similar(exact_E)
-
-near_E = 0.2 + 0.2im
-
-for (j, c) in enumerate(train_cs)
-    H = T + c .* V
-    evals, evecs = eigen(H)
-    i = argmin(abs.(evals .- near_E))
-    train_E[j] = evals[i]
-    EC_basis[:, j] = evecs[:, i]
-end
-
-N_EC = transpose(EC_basis) * EC_basis
-
-for (j, c) in enumerate(target_cs)
-    exact_E[j] = quick_pole_E((p, q) -> c*(V1*g0(R1, p, q) + V2*g0(R2, p, q)), μ; cs_angle=0.5)
-
-    H = T + c .* V
-    H_EC = transpose(EC_basis) * H * EC_basis
-    evals = eigvals(H_EC, N_EC)
-    i = argmin(abs.(evals .- exact_E[j]))
-    extrapolate_E[j] = evals[i]
-end
-
-scatter(real.(train_E), imag.(train_E), label="training")
-scatter!(real.(exact_E), imag.(exact_E), label="exact")
-scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated")
-xlims!(-0.2,0.3)
-ylims!(-0.3,0.3)

berggren.jl

@@ -0,0 +1,22 @@
+using SparseArrays
+include("math.jl")
+
+"berg_bases1/2 are lists of (1+l_max) matrices containing the eigenbases corresponding to 1st and 2nd DOFs respectively,
+js are a list of tuples (j1, j2) corresponding to 1st and 2nd DOFs respectively,
+and ws are the weights needed to evaluate the inner products"
+function get_2p_p1p2_matrix(mesh_size, js, Λ, berg_bases1::Vector{Matrix{ComplexF64}}, berg_bases2::Vector{Matrix{ComplexF64}}, ks::Vector{ComplexF64}, ws::Vector{ComplexF64})
+    # TODO: Cache / precalculate
+    integral1(np, lp, n, l) = sum(ks .* berg_bases1[1 + lp][:, np] .* ws .* berg_bases1[1 + l][:, n])
+    integral2(np, lp, n, l) = sum(ks .* berg_bases2[1 + lp][:, np] .* ws .* berg_bases2[1 + l][:, n])
+
+    basis = iter_prod(js, 1:mesh_size, 1:mesh_size)
+    mat = zeros(ComplexF64, length(basis), length(basis))
+    Threads.@threads for idx in CartesianIndices(mat)
+        (ip, i) = Tuple(idx)
+        ((j1p, j2p), n1p, n2p) = basis[ip]
+        ((j1, j2), n1, n2) = basis[i]
+        val = racahs_reduction_formula(n1p, j1p, n2p, j2p, n1, j1, n2, j2, Λ, integral1, integral2)
+        if !(val ≈ 0); mat[idx] = val; end
+    end
+    return sparse(mat)
+end

berggren_3body_resonance.jl

@@ -0,0 +1,75 @@
+using LinearAlgebra, SparseArrays, Arpack
+include("common.jl")
+include("p_space.jl")
+include("ho_basis.jl")
+include("berggren.jl")
+
+println("No of threads = ", Threads.nthreads())
+atol = 10^-5
+maxevals = 10^5
+R_cutoff = 16
+
+Λ = 0
+m = 1.0
+μ = m/2 # due to simple relative coordinates
+
+target = 4.0766890719636875 - 0.012758927741074495im
+
+V_of_r(r) =  2 * exp(-(r-3)^2 / (1.5)^2)
+V_l(j, k, kp) = Vl_mat_elem(V_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
+
+vertices = [0, 2 - 0.2im, 3, 4]
+subdivisions = [15, 10, 10]
+ks, ws = get_mesh(vertices, subdivisions)
+
+jmax = 4
+tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
+js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
+
+basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
+basis_size = length(js) * length(ks)^2
+@assert length(basis) == basis_size "Something wrong with the basis"
+println("Basis size = $basis_size")
+
+# generate Berggren bases
+@time "berg_bases" begin
+    berg_bases = Vector{Matrix{ComplexF64}}(undef, jmax + 1)
+    for j in 0:jmax
+        _, berg_basis = eigen(get_H_matrix((k, kp) -> V_l(j, k, kp), ks, ws, μ); permute=false, scale=false)
+        N_berg = sum(berg_basis.^2 .* ws, dims=1)
+        berg_bases[1 + j] = berg_basis ./ transpose(sqrt.(N_berg))
+    end
+    to_berg_basis(mat, j) = transpose(berg_bases[1 + j] .* ws) * mat * berg_bases[1 + j]
+end
+
+@time "U_berggren" begin
+    U_blocks = [kron(berg_bases[1 + j1], berg_bases[1 + j2]) for (j1, j2) in js]
+    U = blockdiag(sparse.(U_blocks)...)
+end
+
+@time "Block diagonal part" begin
+    Hb_blocks = [kron_sum(to_berg_basis(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ), j1), to_berg_basis(get_H_matrix((k, kp) -> V_l(j2, k, kp), ks, ws, μ), j2)) for (j1, j2) in js]
+    Hb = blockdiag(sparse.(Hb_blocks)...)
+end
+
+@time "T_cross" T_cross = get_2p_p1p2_matrix(length(ks), js, Λ, berg_bases, berg_bases, ks, ws) ./ (2*μ)
+
+E_max = 30
+μω_global = 0.5
+# due to simple relative coordinates
+μω = μω_global * 2
+μ = m/2
+
+basis_ho = ho_basis_2B(E_max, Λ)
+
+@time "V12_HO" V12_HO =  get_src_V12_matrix(V_of_r, basis_ho, μω_global; atol=10^-6, maxevals=10^5)
+
+@time "W" W = get_W_matrix(basis, basis_ho, μω, μω; weights=true)
+
+@time "V12_p" V12_p = W * V12_HO * transpose(W)
+@time "V12" V12 = transpose(U) * V12_p * U
+
+@time "H" H = Hb + T_cross + V12
+
+@time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
+display(evals)

calculations/2body_data.jl

@@ -0,0 +1,33 @@
+using Roots, DelimitedFiles
+
+include("../EC.jl")
+include("../common.jl")
+include("../p_space.jl")
+
+μ = 0.5
+V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
+
+# determining c0 with EC
+temp_c = range(1.1, 0.9, 3)
+p, w = get_mesh([0, 8], [256])
+H0 = get_T_matrix(p, μ)
+V = get_V_matrix(V_system(1), p, w)
+EC = affine_EC(H0, V, w)
+train!(EC, temp_c; ref_eval=-0.2, CAEC=false, verbose=false)
+quick_extrapolate(c) = minimum(abs2, get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, 0))
+c0 = find_zero(quick_extrapolate, 0.85)
+
+training_c = range(1.2, 0.9, 9) # original: range(1.35, 0.9, 5)
+extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
+
+data_c = vcat(training_c, extrapolating_c)
+data_E = [quick_pole_E(V_system(c)) for c in data_c]
+
+# export to CSV
+file = "temp/2body_data.csv"
+delim = ','
+open(file, "w") do f
+    writedlm(f, ["c" "re_E" "im_E"], delim)
+    writedlm(f, [c0 0 0], delim) # first entry for the threshold
+    writedlm(f, hcat(data_c, real.(data_E), imag.(data_E)), delim)
+end

calculations/3body_Berggren_B2R_EC.jl

@@ -0,0 +1,40 @@
+include("../p_space.jl")
+include("../EC.jl")
+
+Λ = 0
+m = 1.0
+V_of_r(r) =  2 * exp(-(r-3)^2 / (1.5)^2)
+
+ϕ = 0.1
+vertices = [0, 4 * exp(-1im * ϕ), 5, 6]
+subdivisions = [40, 12, 15]
+jmax = 6
+
+E_max = 40
+μω_global = 0.5
+
+H0, weights = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m)
+
+# Vp = perturbation to make the state artificially bound
+Vp_of_r(r) = -exp(-(r/3)^2)
+Vp, _ = get_3b_H_matrix(jacobi, Vp_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m, false, true)
+
+training_c = [2.6, 2.4, 2.2, 2.0, 1.8]
+extrapolating_c = 0.0 : 0.2 : 1.2
+
+training_ref = -2.22 # complete list not needed because identification is simple
+
+extrapolating_ref = [4.076662025307587-0.012709842443350328im,
+    3.613318119833891-0.007335804709990623im,
+    3.1453431847006783-0.004030580410326795im,
+    2.672967129943755-0.00211498327461944im,
+    2.196542557810288-0.0010719835443437104im,
+    1.7164583929199813-0.0005455212208182736im,
+    1.233088227541505-0.0003070320106485624im]
+
+EC = affine_EC(H0, Vp, weights)
+train!(EC, training_c; ref_eval=training_ref, CAEC=true)
+extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
+
+exportCSV(EC, "temp/Berggren_B2R.csv")
+plot(EC, "temp/Berggren_B2R.pdf")

calculations/3body_Berggren_R2R_EC.jl

@@ -0,0 +1,42 @@
+include("../p_space.jl")
+include("../EC.jl")
+
+Λ = 0
+m = 1.0
+V_of_r(r) =  2 * exp(-(r-3)^2 / (1.5)^2)
+
+vertices = [0, 2 - 0.2im, 3, 4]
+subdivisions = [15, 10, 10]
+jmax = 4
+
+E_max = 40
+μω_global = 0.5
+
+H0, weights = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m)
+
+# Vp = perturbation to make the state artificially bound
+Vp_of_r(r) = -exp(-(r/3)^2)
+Vp, _ = get_3b_H_matrix(jacobi, Vp_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m, false, true)
+
+training_c = [1.1, 0.9, 0.7, 0.5]
+extrapolating_c = 0.0 : 0.2 : 1.2
+
+training_ref = [1.4750633616275919 - 0.0003021770706749637im
+    1.9567078295375822 - 0.0007646829108872369im
+    2.4351117758403076 - 0.001281037843108658im
+    2.9096543462392357 - 0.002962488527470604im]
+
+extrapolating_ref = [4.076662025307587-0.012709842443350328im,
+    3.613318119833891-0.007335804709990623im,
+    3.1453431847006783-0.004030580410326795im,
+    2.672967129943755-0.00211498327461944im,
+    2.196542557810288-0.0010719835443437104im,
+    1.7164583929199813-0.0005455212208182736im,
+    1.233088227541505-0.0003070320106485624im]
+
+EC = affine_EC(H0, Vp, weights)
+train!(EC, training_c; ref_eval=training_ref, CAEC=false)
+extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
+
+exportCSV(EC, "temp/Berggren_R2R.csv")
+plot(EC, "temp/Berggren_R2R.pdf")

calculations/3body_CSM_B2R_EC.jl

@@ -0,0 +1,40 @@
+include("../p_space.jl")
+include("../EC.jl")
+
+Λ = 0
+m = 1.0
+V_of_r(r) =  2 * exp(-(r-3)^2 / (1.5)^2)
+
+ϕ = 0.1
+vertices = [0, 6 * exp(-1im * ϕ)]
+subdivisions = [50]
+jmax = 4
+
+E_max = 40
+μω_global = 0.5
+
+H0, weights = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m)
+
+# Vp = perturbation to make the state artificially bound
+Vp_of_r(r) = -exp(-(r/3)^2)
+Vp, _ = get_3b_H_matrix(jacobi, Vp_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m, false, true)
+
+training_c = [2.6, 2.4, 2.2, 2.0, 1.8]
+extrapolating_c = 0.0 : 0.2 : 1.2
+
+training_ref = -2.22 # complete list not needed because identification is simple
+
+exact_E = [4.076662025307587-0.012709842443350328im,
+    3.613318119833891-0.007335804709990623im,
+    3.1453431847006783-0.004030580410326795im,
+    2.672967129943755-0.00211498327461944im,
+    2.196542557810288-0.0010719835443437104im,
+    1.7164583929199813-0.0005455212208182736im,
+    1.233088227541505-0.0003070320106485624im]
+
+EC = affine_EC(H0, Vp, weights)
+train!(EC, training_c; ref_eval=training_ref, CAEC=true)
+extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
+
+exportCSV(EC, "temp/CSM_B2R.csv")
+plot(EC, "temp/CSM_B2R.pdf")

calculations/3body_HO_B2R_ACCC.jl

@@ -0,0 +1,68 @@
+using Roots, LinearAlgebra, Plots
+
+include("../EC.jl")
+include("../common.jl")
+include("../ho_basis.jl")
+
+V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
+Λ = 0
+m = 1.0
+
+ϕ = 0.1
+μω_global = 0.5 * exp(-2im * ϕ)
+E_max = 40
+
+H0 = get_3b_H_matrix(jacobi, V_of_r, μω_global, E_max, Λ, m, true, true)
+
+# Vp = perturbation to make the state artificially bound
+Vp_of_r(r) = -exp(-(r/3)^2)
+@time "Vp" Vp = get_3b_H_matrix(jacobi, Vp_of_r, μω_global, E_max, Λ, m, false, true)
+
+training_ref = -2.22
+extrapolating_ref = [4.076662025307587-0.012709842443350328im,
+    3.613318119833891-0.007335804709990623im,
+    3.1453431847006783-0.004030580410326795im,
+    2.672967129943755-0.00211498327461944im,
+    2.196542557810288-0.0010719835443437104im,
+    1.7164583929199813-0.0005455212208182736im,
+    1.233088227541505-0.0003070320106485624im]
+
+training_c = range(2.8, 1.8, 5)
+extrapolating_c = 0.0 : 0.2 : 1.2
+
+EC = affine_EC(H0, Vp)
+train!(EC, training_c; ref_eval=training_ref, CAEC=true)
+extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
+
+# determining c0 with EC
+approx_c0 = 1.5
+quick_extrapolate(c) = minimum(abs2, get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, 1e-14))
+c0 = find_zero(quick_extrapolate, approx_c0)
+
+order::Int = ceil((length(training_c) - 1) / 2) # order of the Pade approximant
+
+# Solve coefficients as a linear system
+training_k = alt_sqrt.(EC.training_E)
+M_left_element(c, i) = alt_sqrt(c - c0)^i
+M_left = M_left_element.(training_c, (0:order)')
+M_right = -training_k .* M_left[:, 2:end] # remove the first column
+M = hcat(M_left, M_right) # M = [M_left | M_right]
+sol = M \ training_k
+a = sol[1:order+1]
+b = [1; sol[order+2:end]]
+
+# Pade approximant
+polynomial(a, c) = sum(i -> a[i+1] * alt_sqrt(c - c0)^i, 0:order)
+pade_approx(c) = polynomial(a, c) / polynomial(b, c)
+
+# Extrapolate
+extrapolated_k = pade_approx.([extrapolating_c; training_c])
+extrapolated_E = extrapolated_k .^ 2
+
+# Plotting
+scatter(real.(EC.training_E), imag.(EC.training_E), label="training")
+scatter!(real.(EC.exact_E), imag.(EC.exact_E), label="exact")
+scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), label="CAEC", m=:x)
+scatter!(real.(extrapolated_E), imag.(extrapolated_E), label="ACCC", m=:+)
+title!("3-body extrapolation with $(length(training_c)) training points")
+savefig("temp/3body_HO_B2R_ACCC-$(length(training_c)).pdf")

calculations/3body_HO_B2R_EC.jl

@@ -0,0 +1,36 @@
+include("../ho_basis.jl")
+include("../EC.jl")
+
+V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
+Λ = 0
+m = 1.0
+
+ϕ = 0.1
+μω_global = 0.5 * exp(-2im * ϕ)
+E_max = 40
+
+H0 = get_3b_H_matrix(jacobi, V_of_r, μω_global, E_max, Λ, m, true, true)
+
+# Vp = perturbation to make the state artificially bound
+Vp_of_r(r) = -exp(-(r/3)^2)
+@time "Vp" Vp = get_3b_H_matrix(jacobi, Vp_of_r, μω_global, E_max, Λ, m, false, true)
+
+training_ref = -2.22
+
+extrapolating_ref = [4.076662025307587-0.012709842443350328im,
+    3.613318119833891-0.007335804709990623im,
+    3.1453431847006783-0.004030580410326795im,
+    2.672967129943755-0.00211498327461944im,
+    2.196542557810288-0.0010719835443437104im,
+    1.7164583929199813-0.0005455212208182736im,
+    1.233088227541505-0.0003070320106485624im]
+
+training_c = [2.6, 2.4, 2.2, 2.0, 1.8]
+extrapolating_c = 0.0 : 0.2 : 1.2
+
+EC = affine_EC(H0, Vp)
+train!(EC, training_c; ref_eval=training_ref, CAEC=true)
+extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
+
+exportCSV(EC, "temp/HO_B2R.csv")
+plot(EC, "temp/HO_B2R.pdf")

calculations/3body_HO_R2R_EC.jl

@@ -0,0 +1,24 @@
+include("../ho_basis.jl")
+include("../EC.jl")
+
+V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
+Λ = 0
+m = 1.0
+
+μω_global = 0.5 * exp(-2im * pi / 9)
+E_max = 40
+
+T = get_3b_H_matrix(src, V_of_r, μω_global, E_max, Λ, m, true, false)
+V = get_3b_H_matrix(src, V_of_r, μω_global, E_max, Λ, m, false, true)
+
+ref_E = 5.9673 - 0.0006im
+
+training_c = 2.0 : -0.2 : 1.2
+extrapolating_c = 1.05 .- [0.0 : 0.1 : 0.4; 0.45 : 0.05 : 0.60]
+
+EC = affine_EC(T, V)
+train!(EC, training_c; ref_eval=ref_E, CAEC=false)
+extrapolate!(EC, extrapolating_c)
+
+exportCSV(EC, "temp/HO_R2R.csv")
+plot(EC, "temp/HO_R2R.pdf")

calculations/3body_dis_CSM_EC.jl

@@ -0,0 +1,77 @@
+include("../p_space.jl")
+include("../EC.jl")
+
+Λ = 0
+m = 1.0
+
+# Distinguishable particles: V12 = bound and V13 & V23 = resonant
+Vsubsystem_of_r(r) = -2 * exp(-r^2/4)
+Vdecay_of_r(r) = -exp(-r^2 / 3) + exp(-r^2 / 10)
+
+ϕ = 0.1
+vertices = [0, 6 * exp(-1im * ϕ)]
+subdivisions = [50]
+jmax = 4
+
+E_max = 40
+μω_global = 0.4
+
+# Jacobi coordinates
+μ1ω1 = μω_global * 1/2
+μ2ω2 = μω_global * 2
+μ1 = m * 1/2
+μ2 = m * 2/3
+
+atol=10^-5; maxevals=10^5; R_cutoff=16; verbose=true;
+
+###########
+
+verbose && println("No of threads = ", Threads.nthreads())
+
+V_l(j, k, kp) = Vl_mat_elem(Vsubsystem_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
+
+ks, ws = get_mesh(vertices, subdivisions)
+weights = repeat(kron(ws, ws), jmax + 1)
+block_size = length(ks)
+tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
+js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
+basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
+basis_size = length(js) * length(ks)^2
+@assert length(basis) == basis_size "Something wrong with the basis"
+verbose && println("Basis size = $basis_size")
+
+@time "Block diagonal part" begin
+    blocks = [kron_sum(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ1), get_T_matrix(ks, μ2)) for (j1, _) in js]
+    Ha = blockdiag(sparse.(blocks)...)
+end
+
+basis_ho = ho_basis_2B(E_max, Λ)
+verbose && println("HO basis size = ", basis_ho.dim)
+
+@time "V2_HO" V2_HO =  get_jacobi_V2_matrix(Vdecay_of_r, basis_ho, μω_global)
+
+@time "W_right" W_right = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=true)
+@time "W_left" W_left = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=false)
+
+@time "V2" Vb = W_left * V2_HO * transpose(W_right)
+
+###########
+
+training_c = [-0.55, -0.7, -0.85, -1, -1.2]
+extrapolating_c = [0.2, 0.1, 0.0, -0.1, -0.2, -0.3]
+
+ref_E = -0.5173809356244544
+
+exact_ref = [-0.31360746615280954-0.07689284936870341im
+    -0.3233372403877718-0.06011323914565665im
+    -0.339615582074795-0.04239442037174759im
+    -0.36333816534241997-0.02648721825958402im
+    -0.39376650561322885-0.014382935339817332im	
+    -0.4299479825535172-0.006510710745123606im]
+
+EC = affine_EC(Ha, Vb)
+train!(EC, training_c; ref_eval=ref_E, CAEC=true)
+extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_ref)
+
+exportCSV(EC, "temp/dis_CSM_B2R.csv")
+plot(EC, "temp/dis_CSM_B2R.pdf")

calculations/3body_dis_HO_EC.jl

@@ -0,0 +1,52 @@
+include("../ho_basis.jl")
+include("../EC.jl")
+
+Λ = 0
+m = 1.0
+
+# Distinguishable particles: V12 = bound and V13 & V23 = resonant
+Va_of_r(r) = -2 * exp(-r^2/4)
+Vb_of_r(r) = -exp(-r^2 / 3) + exp(-r^2 / 10)
+
+E_max = 40
+μω_global = 0.4 * exp(-2im * pi / 9)
+
+# due to Jacobi coordinates
+μ1ω1 = μω_global * 1/2
+μ2ω2 = μω_global * 2
+μ1 = m * 1/2
+μ2 = m * 2/3
+
+println("No of threads = ", Threads.nthreads())
+
+basis = ho_basis_2B(E_max, Λ)
+println("Basis size = ", basis.dim)
+
+@time "T1" T1 = get_sp_T_matrix(basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; μω_gen=μ1ω1, μ=μ1)
+@time "T2" T2 = get_sp_T_matrix(basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; μω_gen=μ2ω2, μ=μ2)
+
+@time "Va" Va = get_jacobi_V1_matrix(Va_of_r, basis, μ1ω1)
+@time "Vb" Vb = get_jacobi_V2_matrix(Vb_of_r, basis, μω_global)
+
+@time "Ha" Ha = T1 + T2 + Va
+@time "Eigenvalues" target_evals, _ = eigs(Ha, nev=5, ncv=50, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
+display(target_evals)
+
+training_c = [-0.55, -0.7, -0.85, -1, -1.2]
+extrapolating_c = [0.2, 0.1, 0.0, -0.1, -0.2, -0.3]
+
+ref_E = -0.5173809356244544
+
+exact_ref = [-0.3136074661528041-0.07689284936868852im
+    -0.323337240387771-0.06011323914564878im
+    -0.33961558207479553-0.04239442037174764im
+    -0.3633381653424224-0.026487218259589693im
+    -0.393766505613234-0.014382935339825854im
+    -0.42994798255352606-0.006510710745131777im]
+
+EC = affine_EC(Ha, Vb)
+train!(EC, training_c; ref_eval=ref_E, CAEC=true) # try CAEC=false !!!
+extrapolate!(EC, extrapolating_c, precalculated_exact_E = exact_ref)
+
+exportCSV(EC, "temp/dis_HO_B2R.csv")
+plot(EC, "temp/dis_HO_B2R.pdf")

calculations/ACCC.jl

@@ -0,0 +1,50 @@
+using Roots, LinearAlgebra, Plots
+
+include("../EC.jl")
+include("../common.jl")
+include("../p_space.jl")
+
+μ = 0.5
+V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
+
+# determining c0 with EC
+temp_c = range(1.1, 0.9, 3)
+p, w = get_mesh([0, 8], [256])
+H0 = get_T_matrix(p, μ)
+V = get_V_matrix(V_system(1), p, w)
+EC = affine_EC(H0, V, w)
+train!(EC, temp_c; ref_eval=-0.2, CAEC=false, verbose=false)
+quick_extrapolate(c) = minimum(abs2, get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, 0))
+c0 = find_zero(quick_extrapolate, 0.85)
+
+# Calculation of training and extrapolating E
+training_c = range(1.2, 0.9, 9) # original: range(1.35, 0.9, 5)
+training_E = [quick_pole_E(V_system(c)) for c in training_c]
+training_k = alt_sqrt.(2μ .* training_E)
+
+extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
+exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
+
+order::Int = ceil((length(training_c) - 1) / 2) # order of the Pade approximant
+
+# Solve coefficients as a linear system
+M_left_element(c, i) = alt_sqrt(c - c0)^i
+M_left = M_left_element.(training_c, (0:order)')
+M_right = -training_k .* M_left[:, 2:end] # remove the first column
+M = hcat(M_left, M_right) # M = [M_left | M_right]
+sol = M \ training_k
+a = sol[1:order+1]
+b = [1; sol[order+2:end]]
+
+# Pade approximant
+polynomial(a, c) = sum(i -> a[i+1] * alt_sqrt(c - c0)^i, 0:order)
+pade_approx(c) = polynomial(a, c) / polynomial(b, c)
+
+# Extrapolate
+extrapolated_k = pade_approx.([training_c; extrapolating_c])
+extrapolated_E = (extrapolated_k .^ 2) / (2μ)
+
+# Plotting
+scatter(real.(training_E), imag.(training_E), label="training")
+scatter!(real.(exact_E), imag.(exact_E), label="exact")
+scatter!(real.(extrapolated_E), imag.(extrapolated_E), label="extrapolated", m=:star5)

calculations/B2R_Berggren_poles.jl

@@ -0,0 +1,33 @@
+include("../EC.jl")
+include("../common.jl")
+include("../p_space.jl")
+
+# contour
+p, w = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128])
+
+μ = 0.5
+V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
+
+# generating a Berggren basis with a pole using the same system
+basis_c = 0.6
+basis_E, berg_basis = eigen(get_H_matrix(V_system(basis_c), p, w); permute=false, scale=false)
+basis_p = sqrt.(basis_E)
+N_berg = sqrt.(diag(transpose(berg_basis .* w) * berg_basis))
+berg_basis = berg_basis ./ transpose(N_berg)
+berg_basis_w = berg_basis .* w
+
+H0 = transpose(berg_basis_w) * get_T_matrix(p, μ) * berg_basis
+V = transpose(berg_basis_w) * get_V_matrix(V_system(1), p, w) * berg_basis
+
+training_c = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
+extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
+
+training_ref = -0.26
+extrapolating_ref = [quick_pole_E(V_system(c)) for c in extrapolating_c]
+
+EC = affine_EC(H0, V)
+train!(EC, training_c; ref_eval=training_ref, CAEC=true)
+extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
+
+exportCSV(EC, "temp/2b_GSM_B2R.csv")
+plot(EC, "temp/2b_GSM_B2R.pdf"; basis_points=basis_E, xlims=(0, 0.3), ylims=(-0.120, 0.020))

calculations/B2R_comparison.jl

@@ -0,0 +1,30 @@
+include("../EC.jl")
+include("../common.jl")
+include("../p_space.jl")
+
+berggren_mesh = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128])
+csm_mesh = get_mesh([0, 8 - 3im], [512])
+
+for (mesh, name) in zip((berggren_mesh, csm_mesh), ("beggren", "csm"))
+    p, w = mesh
+    mesh_E = p.*p ./ (2*0.5)
+
+    μ = 0.5
+    V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
+    
+    H0 = get_T_matrix(p, μ)
+    V = get_V_matrix(V_system(1), p, w)
+    
+    training_c = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
+    extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
+
+    training_ref = [quick_pole_E(V_system(c)) for c in training_c]
+    exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
+
+    EC = affine_EC(H0, V, w)
+    train!(EC, training_c; ref_eval=training_ref, CAEC=true)
+    extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
+    
+    #exportCSV(EC, "temp/2b_comparison_$name.csv")
+    plot(EC, "temp/2b_comparison_$name.pdf"; basis_contour=mesh_E, xlims=(-0.3,0.3), ylims=(-0.120,0.020))
+end

calculations/PMM.py

@@ -0,0 +1,81 @@
+#%%
+import pandas as pd
+import torch
+import numpy as np
+
+#%%
+df = pd.read_csv('../temp/2body_data.csv').sort_values(by='c')
+df['E'] = df['re_E'] + 1j * df['im_E']
+train_data = df[df['re_E'] < 0]
+target_data = df[df['re_E'] > 0]
+
+train_cs = train_data['c'].to_numpy()
+train_Es = torch.tensor(train_data['E'].to_numpy(), dtype=torch.complex128)
+
+#%%
+# hyperparameters
+N = 9
+
+# initialize random Hamiltonians
+H0 = torch.randn(N, N, dtype=torch.complex128)
+H0 = (H0 + torch.transpose(H0, 0, 1)).requires_grad_() # symmetric
+H1 = torch.randn(N, N, dtype=torch.complex128)
+H1 = (H1 + torch.transpose(H1, 0, 1)).requires_grad_() # symmetric
+
+#%%
+# training
+
+# generate a set of c values to follow by subdividing the training cs
+subdivisions = 3
+c_steps = np.concatenate([np.linspace(start, stop, subdivisions, endpoint=False) for (start, stop) in zip(train_cs, train_cs[1:])])
+c_steps = np.append(c_steps, train_cs[-1])
+
+lr = 0.05
+epochs = 100000
+for epoch in range(epochs):
+    Es = torch.empty(len(train_data), dtype=torch.complex128)
+    current_E = 0.0 # start at the threshold
+    for c in c_steps:
+        H = H0 + c * H1
+        evals = torch.linalg.eigvals(H)
+        current_E = evals[torch.argmin(torch.abs(evals - current_E))]
+        if np.any(c == train_cs):
+            index = np.where(c == train_cs)[0][0]
+            Es[index] = current_E
+
+    loss = ((Es - train_Es).abs() ** 2).sum()
+
+    if epoch % 1000 == 0:
+        print(f"Training {(epoch+1)/epochs:.1%} \t Loss: {loss}")
+
+    if H0.grad is not None:
+        H0.grad.zero_()
+    if H1.grad is not None:
+        H1.grad.zero_()
+    loss.backward()
+
+    with torch.no_grad():
+        H0 -= lr * H0.grad
+        H1 -= lr * H1.grad
+
+# %%
+# evaluate for all points
+all_c = torch.tensor(df['c'].values, dtype=torch.float64)
+exact_E = torch.tensor(df['E'].values, dtype=torch.complex128)
+pred_Es = torch.empty(len(df), dtype=torch.complex128)
+with torch.no_grad():
+    for (index, (c, E)) in enumerate(zip(all_c, exact_E)):
+        H = H0 + c * H1
+        evals = torch.linalg.eigvals(H)
+        i = torch.argmin(torch.abs(evals - E)) # TODO: more robust way to identify the eigenvector
+        pred_Es[index]= evals[i]
+
+# %%
+# plot the results
+import matplotlib.pyplot as plt
+plt.scatter(train_data['re_E'], train_data['im_E'], label='training')
+plt.scatter(target_data['re_E'], target_data['im_E'], label='target')
+plt.scatter(pred_Es.real, pred_Es.imag, marker='x', label='predicted')
+plt.legend()
+
+# %%

calculations/R2R_Berggren.jl

@@ -0,0 +1,28 @@
+include("../EC.jl")
+include("../common.jl")
+include("../p_space.jl")
+
+# contour
+vertices = [0, 0.4 - 0.15im, 0.8, 6]
+subdivisions = [128, 128, 128]
+p, w = get_mesh(vertices, subdivisions)
+mesh_E = p.*p ./ (2*0.5)
+
+μ = 0.5
+V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
+
+H0 = get_T_matrix(p, μ)
+V = get_V_matrix(V_system(1), p, w)
+
+training_c = range(0.75, 0.45, 5)
+extrapolating_c = range(0.40, 0.25, 5)
+
+training_ref = [quick_pole_E(V_system(c)) for c in training_c]
+exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
+
+EC = affine_EC(H0, V, w)
+train!(EC, training_c; ref_eval=training_ref, CAEC=false)
+extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
+
+#exportCSV(EC, "temp/2b_R2R.csv")
+plot(EC, "temp/2b_R2R.pdf"; basis_contour=mesh_E, xlims=(0, 1))

calculations/R2R_Berggren_poles.jl

@@ -0,0 +1,33 @@
+include("../EC.jl")
+include("../common.jl")
+include("../p_space.jl")
+
+# contour
+p, w = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128])
+
+μ = 0.5
+V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
+
+# generating a Berggren basis with a pole using the same system
+basis_c = 0.6
+basis_E, berg_basis = eigen(get_H_matrix(V_system(basis_c), p, w); permute=false, scale=false)
+basis_p = sqrt.(basis_E)
+N_berg = sqrt.(diag(transpose(berg_basis .* w) * berg_basis))
+berg_basis = berg_basis ./ transpose(N_berg)
+berg_basis_w = berg_basis .* w
+
+H0 = transpose(berg_basis_w) * get_T_matrix(p, μ) * berg_basis
+V = transpose(berg_basis_w) * get_V_matrix(V_system(1), p, w) * berg_basis
+
+training_c = range(0.79, 0.66, 4) # original: range(1.35, 0.9, 5)
+extrapolating_c = range(0.62, 0.40, 6) # original: range(0.75, 0.40, 8)
+
+training_ref = [quick_pole_E(V_system(c)) for c in training_c]
+extrapolating_ref = [quick_pole_E(V_system(c)) for c in extrapolating_c]
+
+EC = affine_EC(H0, V)
+train!(EC, training_c; ref_eval=training_ref, CAEC=false)
+extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
+
+exportCSV(EC, "temp/2b_GSM_R2R.csv")
+plot(EC, "temp/2b_GSM_R2R.pdf"; basis_points=basis_E, xlims=(0, 0.3), ylims=(-0.120, 0.020))

calculations/XZ/2body_HO.jl

@@ -0,0 +1,39 @@
+using Plots
+
+include("../../EC.jl")
+include("../../ho_basis.jl")
+include("../../p_space.jl")
+
+angle = 0.25 * pi # DOESN'T WORK WITHOUT ROTATION
+μω_gen = 0.5 * exp(-2im * angle)
+μ = 0.5
+l = 0
+V1 = -5
+R1 = sqrt(3)
+V2 = 2
+R2 = sqrt(10)
+n_max = 15
+
+ns = collect(0:n_max)
+ls = fill(l, n_max + 1)
+
+T = get_sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ)
+V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); μω_gen=μω_gen) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); μω_gen=μω_gen)
+
+n_EC = 8
+train_cs = (0.7 .+ 0.05 * randn(n_EC)) - 1im * (0.2 .+ 0.05 * randn(n_EC))
+target_cs = [0.5]
+
+near_E = 0.2 + 0.2im
+exact_E = [0.20845136860234303 - 0.07100640993695649im]
+
+EC = affine_EC(T, V; ensemble_size=32)
+train!(EC, train_cs; ref_eval=near_E, CAEC=false)
+extrapolate!(EC, target_cs; precalculated_exact_E=exact_E)
+
+plot(EC; xlims=(0,0.3), ylims=(-0.3,0.3))
+hline!([0], color=:red, label="continuum")
+xlabel!("Re(E)")
+ylabel!("Im(E)")
+plot!(legend=:bottomleft)
+savefig("temp/2b_HO_XZ.pdf")

calculations/XZ/2body_p_space.jl

@@ -0,0 +1,50 @@
+using Plots
+
+include("../../EC.jl")
+include("../../ho_basis.jl")
+include("../../p_space.jl")
+
+# paramters of the system
+
+angle = 0.0
+μ = 0.5
+l = 0
+V1 = -5
+R1 = sqrt(3)
+V2 = 2
+R2 = sqrt(10)
+
+n_EC = 8
+train_cs = (0.7 .+ 0.03 * randn(n_EC)) - 1im * (0.2 .+ 0.03 * randn(n_EC))
+near_E = 0.2 + 0.2im
+
+target_c = 0.5
+exact_E = 0.20845136860234303 - 0.07100640993695649im
+
+vertices = [0, 4 * exp(-1im * angle)]
+subdivisions = [256]
+ks, ws = get_mesh(vertices, subdivisions)
+
+V_of_r(r) = V1 * exp(-r^2 / R1^2) + V2 * exp(-r^2 / R2^2)
+V_mat_elem(k, kp) = Vl_mat_elem(V_of_r, l, k, kp; atol=10^-5, maxevals=10^5, R_cutoff=16)
+V = get_V_matrix(V_mat_elem, ks, ws)
+T = get_T_matrix(ks, μ)
+
+EC_p_space = affine_EC(T, V)
+train!(EC_p_space, train_cs; ref_eval=near_E, CAEC=false)
+extrapolate!(EC_p_space, [target_c]; precalculated_exact_E=[exact_E])
+
+# Plotting
+
+theme(:dark) # Set the global theme to dark
+
+scatter([real(exact_E)], [imag(exact_E)], label="exact", marker=:circle, markercolor=:white, bg = :black) # black background
+scatter!(real.(EC_p_space.training_E), imag.(EC_p_space.training_E), label="training", marker=:circle, color=:blue)
+scatter!(real.(EC_p_space.extrapolated_E), imag.(EC_p_space.extrapolated_E), label="extrapolated", marker=:x, color=:green)
+hline!([0], color=:red, label="continuum")
+plot!(legend=:bottomleft)
+xlabel!("Re(E)")
+ylabel!("Im(E)")
+xlims!(0, 0.3)
+ylims!(-0.3, 0.3)
+savefig("temp/2body_p_space.pdf")

calculations/XZ/3body_HO.jl

@@ -0,0 +1,36 @@
+include("../../ho_basis.jl")
+include("../../EC.jl")
+
+V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
+Λ = 0
+m = 1.0
+
+ϕ = 0.1 # DOESN'T WORK WITHOUT ROTATION
+μω_global = 0.5 * exp(-2im * ϕ)
+E_max = 40
+
+H0 = get_3b_H_matrix(jacobi, V_of_r, μω_global, E_max, Λ, m, true, true)
+
+# Vp = perturbation to make the state artificially bound
+Vp_of_r(r) = -exp(-(r/3)^2)
+@time "Vp" Vp = get_3b_H_matrix(jacobi, Vp_of_r, μω_global, E_max, Λ, m, false, true)
+
+training_ref = 2 + 0.5im
+
+exact_E = [4.076642792419057-0.012998408352259658im,
+    3.6129849325287-0.007397677539402868im,
+    3.145212908643357-0.0038660337822150753im,	
+    2.6729225739451596-0.0021090370393881063im,	
+    2.196385760253282-0.0010430088245526555im,	
+    1.7162659936896967-0.0004515351140200029,	
+    1.2329926791785895-0.00017698044022813525im]
+
+training_c = [0.6 - 0.16im] .+ 0.04 .* (randn(8) .+ 0.5im * randn(8))
+extrapolating_c = 0.0 : 0.2 : 1.2
+
+EC = affine_EC(H0, Vp)
+train!(EC, training_c; ref_eval=training_ref, CAEC=false)
+extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
+
+exportCSV(EC, "temp/3b_HO_XZ.csv")
+plot(EC, "temp/3b_HO_XZ.pdf")

calculations/XZ/3body_p_space.jl

@@ -0,0 +1,45 @@
+include("../../p_space.jl")
+include("../../EC.jl")
+
+using Arpack
+
+# target = 4.0766890719636875 - 0.012758927741074495im
+
+Λ = 0
+m = 1.0
+V_of_r(r) =  2 * exp(-(r-3)^2 / (1.5)^2)
+
+vertices = [0, 2 - 0.2im, 3, 4] # TODO: real contour instead of Berggren basis
+subdivisions = [15, 10, 10]
+jmax = 4
+
+E_max = 40
+μω_global = 0.5
+
+@time "H0" H0, _ = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m, true, true)
+
+# Vp = perturbation to make the state artificially bound
+Vp_of_r(r) = -exp(-(r/3)^2)
+@time "Vp" Vp, _ = get_3b_H_matrix(jacobi, Vp_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m, false, true)
+
+training_ref = 2 + 0.5im
+
+exact_E = [4.076642792419057-0.012998408352259658im,
+    3.6129849325287-0.007397677539402868im,
+    3.145212908643357-0.0038660337822150753im,	
+    2.6729225739451596-0.0021090370393881063im,	
+    2.196385760253282-0.0010430088245526555im,	
+    1.7162659936896967-0.0004515351140200029,	
+    1.2329926791785895-0.00017698044022813525im]
+
+training_c = [-1.5 - 0.5im] .+ (randn(8) .+ 0.05im * randn(8))
+extrapolating_c = 0.0 : 0.2 : 1.2
+
+EC = affine_EC(H0, Vp)
+train!(EC, training_c; ref_eval=training_ref, CAEC=false)
+extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
+
+exportCSV(EC, "temp/3b_p_space_XZ.csv")
+plot(EC, "temp/3b_p_space_XZ.pdf")
+
+# Results: training points are all over the place, and extrapolated values are garbage.

calculations/XZ/p_space_vs_HO.jl

@@ -0,0 +1,71 @@
+using Plots
+
+include("../../EC.jl")
+include("../../ho_basis.jl")
+include("../../p_space.jl")
+
+# paramters of the system
+
+angle = 0.0 * pi
+μ = 0.5
+l = 0
+V1 = -5
+R1 = sqrt(3)
+V2 = 2
+R2 = sqrt(10)
+
+n_EC = 8
+train_cs = (0.7 .+ 0.03 * randn(n_EC)) - 1im * (0.2 .+ 0.03 * randn(n_EC))
+near_E = 0.2 + 0.2im
+
+target_c = 0.5
+exact_E = 0.20845136860234303 - 0.07100640993695649im
+
+# HO basis
+
+global EC_HO
+begin
+    println("HO basis calculation")
+    μω_gen = 0.5 * exp(-1im * angle)
+    n_max = 40
+    ns = collect(0:n_max)
+    ls = fill(l, n_max + 1)
+
+    T = get_sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ)
+    V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); μω_gen=μω_gen) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); μω_gen=μω_gen)
+
+    global EC_HO = affine_EC(T, V)
+    train!(EC_HO, train_cs; ref_eval=near_E, CAEC=false)
+    extrapolate!(EC_HO, [target_c]; precalculated_exact_E=[exact_E])
+end
+
+# p-space
+
+global EC_p_space
+begin
+    println("p-space calculation")
+    vertices = [0, 4 * exp(-1im * angle)]
+    subdivisions = [256]
+    ks, ws = get_mesh(vertices, subdivisions)
+
+    V_of_r(r) = V1 * exp(-r^2 / R1^2) + V2 * exp(-r^2 / R2^2)
+    V_mat_elem(k, kp) = Vl_mat_elem(V_of_r, l, k, kp; atol=10^-5, maxevals=10^5, R_cutoff=16)
+    V = get_V_matrix(V_mat_elem, ks, ws)
+    T = get_T_matrix(ks, μ)
+
+    global EC_p_space = affine_EC(T, V)
+    train!(EC_p_space, train_cs; ref_eval=near_E, CAEC=false)
+    extrapolate!(EC_p_space, [target_c]; precalculated_exact_E=[exact_E])
+end
+
+# Plotting
+
+scatter([real(exact_E)], [imag(exact_E)], label="Exact", marker=:circle, markercolor=:white)
+scatter!(real.(EC_HO.training_E), imag.(EC_HO.training_E), label="HO basis training", marker=:circle, color=:blue)
+scatter!(real.(EC_HO.extrapolated_E), imag.(EC_HO.extrapolated_E), label="HO basis extrapolated", marker=:x, color=:blue)
+scatter!(real.(EC_p_space.training_E), imag.(EC_p_space.training_E), label="p-space training", marker=:circle, color=:red)
+scatter!(real.(EC_p_space.extrapolated_E), imag.(EC_p_space.extrapolated_E), label="p-space extrapolated", marker=:x, color=:red)
+plot!(legend=:bottomleft)
+xlabel!("Re(E)")
+ylabel!("Im(E)")
+savefig("temp/2b_p_space_vs_HO.pdf")

common.jl

@@ -0,0 +1,106 @@
+using LinearAlgebra, DelimitedFiles, SparseArrays
+
+@enum coordinate_system jacobi src
+
+"Square root function with the branch cut along the postive imaginary axis"
+alt_sqrt(x::Number)::ComplexF64 = sqrt(im * x) / sqrt(im)
+
+"Sum over array while minimizing catastrophic cancellation as much as possible"
+function better_sum(arr::Array{T}) where T<:Real
+    pos_arr = arr[arr .> 0]
+    neg_arr = arr[arr .< 0]
+
+    sort!(pos_arr)
+    sort!(neg_arr, rev=true)
+
+    return sum(pos_arr) + sum(neg_arr)
+end
+
+better_sum(arr::Array{ComplexF64}) = better_sum(real.(arr)) + 1im * better_sum(imag.(arr))
+
+"The triangle inequality for angular momenta"
+triangle_ineq(l1, l2, L) = abs(l1 - l2) ≤ L && L ≤ (l1 + l2)
+
+"Index of the nearest value in a list to a given reference point"
+nearestIndex(list::Array, ref) = argmin(norm.(list .- ref))
+
+"Nearest value in a list to a given reference point"
+nearest(list::Array, ref) = list[nearestIndex(list, ref)]
+
+"Simple implementation of the Kronecker sum"
+function kron_sum(A::AbstractMatrix, B::AbstractMatrix)
+    @assert size(A, 1) == size(A, 2) && size(B, 1) == size(B, 2) "Matrices should be square"
+    return kron(A, I(size(B, 1))) + kron(I(size(A, 1)), B)
+end
+
+"Flattened vector version of Iterators.product(...) with index hierachy reversed -- leftmost index has the highest hierachy"
+iter_prod(args...) = reverse.(collect(Iterators.product(reverse(args)...))[:])
+
+"Export CSV data for a scatter plot taking in data as a list of complex vectors (x=Re, y=Im), and a list of corresponding labels, typically used for EC results"
+function exportCSV(file::String, data, labels=nothing)
+    columns = ["x" "y" "label"]
+
+    if isnothing(labels)
+        columns[end] = ""
+        labels = fill("", length(data))
+    end
+    
+    open(file, "w") do f
+        writedlm(f, columns)
+        for (d, label) in zip(data, labels)
+            l = fill(label, length(d))
+            writedlm(f, hcat(reim(d)..., l))
+        end
+    end
+end
+
+"In-place c-orthogonalization via (modified) Gram-Schmidt. Only significant vectors are returned (c-normalized).
+The number of significant vectors to return are determined by the original singular values (compared to the threshold)."
+function gram_schmidt!(vecs::Vector{Vector{T}}, ws=ones(length(vecs[1])), threshold=1e-5; verbose=false) where T<: Number
+    c_product(i, j) = sum(vecs[i] .* ws .* vecs[j])
+    norm(i) = c_product(i, i)
+    proj(i, j) = (c_product(i, j) / norm(i)) .* vecs[i] # component of vec[i] in vec[j]
+
+    # initial normalization
+    for i in eachindex(vecs)
+        vecs[i] ./= sqrt(norm(i))
+    end
+
+    verbose && println("Absolute singular values = $(round.(c_singular_values(vecs, ws); sigdigits=1))")
+    target_dim = c_rank(vecs, ws, threshold)
+    verbose && println("Target dimensionality = $target_dim")
+
+    selected_vecs_i = Integer[]
+
+    while length(selected_vecs_i) < target_dim
+        i = argmax(i -> abs(norm(i)), setdiff(eachindex(vecs), selected_vecs_i)) # find the largest vector from the remaining
+        push!(selected_vecs_i, i)
+        for j in setdiff(eachindex(vecs), selected_vecs_i)
+            vecs[j] .-= proj(i, j)
+        end
+    end
+    
+    verbose && println("Absolute norms of selected vectors = $(round.(selected_vecs_i .|> norm .|> abs; sigdigits=1))")
+    verbose && println("Absolute norms of dropped vectors = $(round.(setdiff(eachindex(vecs), selected_vecs_i) .|> norm .|> abs; sigdigits=1))")
+
+    # final normalization
+    for i in selected_vecs_i
+        vecs[i] ./= sqrt(norm(i))
+    end
+
+    return vecs[selected_vecs_i]
+end
+
+"Same as above but the basis is provided as a matrix instead of an array of vectors"
+gram_schmidt!(vecs::Matrix{T}, ws=ones(size(vecs, 1)), threshold=1e-5; verbose=false) where T<: Number = hcat(gram_schmidt!([vecs[:, i] for i in axes(vecs, 2)], ws, threshold; verbose=verbose)...)
+
+"Rank of a basis set (pre-normalized) under c-product"
+c_rank(vecs, ws, threshold=1e-8) = count(c_singular_values(vecs, ws) .> threshold)
+
+"Singular values (magnitudes) of a basis set (pre-normalized) under c-product"
+function c_singular_values(vecs, ws)
+    basis = hcat(vecs...)
+    N = transpose(basis) * spdiagm(ws) * basis
+    singular_values = eigvals(N) .|> abs .|> sqrt
+    return singular_values
+end

helper.jl

@@ -1,12 +0,0 @@
-"Sum over array while minimizing catastrophic cancellation as much as possible"
-function better_sum(arr::Array{Float64})
-    pos_arr = arr[arr .> 0]
-    neg_arr = arr[arr .< 0]
-
-    sort!(pos_arr)
-    sort!(neg_arr, rev=true)
-
-    return sum(pos_arr) + sum(neg_arr)
-end
-
-better_sum(arr::Array{ComplexF64}) = better_sum(real.(arr)) + 1im * better_sum(imag.(arr))

ho_basis.jl

@@ -1,64 +1,74 @@
 using SparseArrays
-using SpecialFunctions
-include("helper.jl")
-include("osbrackets.jl")
+using QuadGK
+using LRUCache
+include("common.jl")
+include("math.jl")
 
-# Gaussian potentials in HO space
-inv_factorial(n) = Iterators.prod(inv.(1:n))
-sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
-N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + l + 3/2)
-Talmi(l, R, k1, k2; ω=1.0) = (-1)^(k1 + k2) * (1 + 1/(ω * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
-V_Gaussian(R, l, n1, n2; ω=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; ω=ω) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
+"2-body HO basis (used for 3-body systems without the CM DOF)"
+struct ho_basis_2B
+    E_max::Int
+    Λ::Int
+    dim::Int # dimensionality of the basis
+    Es::Vector{Int}
+    n1s::Vector{Int}
+    l1s::Vector{Int}
+    n2s::Vector{Int}
+    l2s::Vector{Int}
 
-function get_sp_basis(E_max)
-    Es = Int[]
-    ns = Int[]
-    ls = Int[]
+    function ho_basis_2B(E_max, Λ=-1)
+        Es = Int[]
+        n1s = Int[]
+        l1s = Int[]
+        n2s = Int[]
+        l2s = Int[]
 
-    # E = 2*n + l
-    for E in 0 : E_max
-        for n in 0 : E ÷ 2
-            l = E - 2*n
-            push!(Es, E)
-            push!(ns, n)
-            push!(ls, l)
-        end
-    end
-
-    return (Es, ns, ls)
-end
-
-function get_2p_basis(E_max)
-    Es = Int[]
-    n1s = Int[]
-    l1s = Int[]
-    n2s = Int[]
-    l2s = Int[]
-
-    # E = 2*n1 + l1 + 2*n2 + l2
-    for E in 2*E_max : -2 : 0
-        for n1 in 0 : E ÷ 2
-            for n2 in 0 : (E - 2*n1) ÷ 2
-                for l1 in 0 : (E - 2*n1 - 2*n2)
-                    l2 = E - 2*n1 - 2*n2 - l1
-                    push!(Es, E)
-                    push!(n1s, n1)
-                    push!(l1s, l1)
-                    push!(n2s, n2)
-                    push!(l2s, l2)
+        # E = 2*n1 + l1 + 2*n2 + l2
+        for E in E_max : -2 : 0 # same parity states only
+            for n1 in 0 : E ÷ 2
+                for n2 in 0 : (E - 2*n1) ÷ 2
+                    for l1 in 0 : (E - 2*n1 - 2*n2)
+                        l2 = E - 2*n1 - 2*n2 - l1
+                        if Λ≥0 && !triangle_ineq(l1, l2, Λ); continue; end
+                        push!(Es, E)
+                        push!(n1s, n1)
+                        push!(l1s, l1)
+                        push!(n2s, n2)
+                        push!(l2s, l2)
+                    end
                 end
             end
         end
-    end
 
-    return (Es, n1s, l1s, n2s, l2s)
+        return new(E_max, Λ, length(Es), Es, n1s, l1s, n2s, l2s)
+    end
 end
 
-function sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), ω=1.0, μ=1.0)
+"Number of possible distinct matrix elements for a given l."
+function cache_size_l(E_max::Int, l::Int)::Int
+    n_max = (E_max - l) ÷ 2
+    return (n_max * (n_max + 1)) ÷ 2
+end
+
+"Number of possible distinct matrix elements for all l."
+cache_size(E_max::Int)::Int = sum(l -> cache_size_l(E_max, l), 0 : E_max)
+
+"Preallocation of cache for PE matrix elements."
+prealloc_V_cache(E_max::Int, dtype::DataType=Float64) = LRU{Tuple{UInt8, UInt8, UInt8}, dtype}(maxsize=cache_size(E_max))
+
+function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
+    const_part = sqrt(μω_gen) * ho_basis_const(l, n1) * ho_basis_const(l, n2)
+    integrand(r) = ho_basis_func(l, n1, sqrt(μω_gen) * r) * ho_basis_func(l, n2, sqrt(μω_gen) * r) * V_of_r(r)
+    (integral, _) = quadgk(integrand, 0, Inf; atol=atol, maxevals=maxevals)
+    return const_part * integral
+end
+
+"KE matrix for a single DOF. Set kron_deltas=[other quantum numbers] for other DOFs which the operator does not act on.
+    E.g. get_sp_T_matrix(n1s, l1s, kron_deltas=[n2s l2s])"
+function get_sp_T_matrix(ns, ls, kron_deltas=[]; μω_gen=1.0, μ=1.0)
     mat = spzeros(length(ns), length(ns))
     for idx in CartesianIndices(mat)
-        if !mask[idx]; continue; end
         (i, j) = Tuple(idx)
+        all(arr -> arr[i]==arr[j], kron_deltas) || continue # check if all Kronecker deltas are non-zero
         if ls[i] == ls[j]
             if ns[i] == ns[j]
                 mat[idx] = ns[j] + ls[i]/2 + 3/4
@@ -68,44 +78,174 @@ function sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), ω=1.0, μ=1.0)
             end
         end
     end
-    return (ω / μ) .* mat
+    return (μω_gen / μ) .* mat
 end
 
-function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float64)
+"PE matrix for a single DOF. Set kron_deltas=[other quantum numbers] for other DOFs which the operator does not act on.
+    E.g. get_sp_V_matrix(n1s, l1s, kron_deltas=[n2s l2s])
+    Providing a preallocated cache is optional. Otherwise, provided E_max value is used to initialize one (defaults to 100)."
+function get_sp_V_matrix(V_l, ns, ls, kron_deltas=[]; dtype=Float64, E_max=100, cache=prealloc_V_cache(E_max, dtype))
     mat = zeros(dtype, length(ns), length(ns))
-    Threads.@threads for idx in CartesianIndices(mat)
-        if !mask[idx]; continue; end
+    Threads.@threads for idx in CartesianIndices(mat) 
         (i, j) = Tuple(idx)
+        all(arr -> arr[i]==arr[j], kron_deltas) || continue # check if all Kronecker deltas are non-zero
         if ls[i] == ls[j]
-            mat[idx] = V_l(ls[i], ns[i], ns[j])
+            l = UInt8(ls[i])
+            n1, n2 = UInt8.(minmax(ns[i], ns[j])) # assuming transpose symmetry
+            mat[idx] = (get!(cache, (l, n1, n2)) do; V_l(l, n1, n2); end)
         end
     end
     return sparse(mat)
 end
 
-function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
-    Emax = maximum(Es) ÷ 2      # TODO: Too many steps. Simplify.
+function Moshinsky_transform(basis::ho_basis_2B)
+    NQMAX = maximum(basis.Es)
+    @assert all(mod.(basis.Es, 2) .== mod(NQMAX, 2)) "Can only admit basis states with same parity"
 
-    ul1s = unique(l1s)
-    ul2s = unique(l2s)
-    BRACs = Matrix{Array}(undef, length(ul1s), length(ul2s))
-    for l1 in ul1s
-        for l2 in ul2s
-            BRACs[l1, l2] = cal_BRAC(Emax, Λ, l1, l2)
-        end
-    end
+    LMIN = basis.Λ
+    LMAX = basis.Λ
+    CO = 1/sqrt(2)
+    SI = 1/sqrt(2)
 
-    mat = spzeros(length(Es), length(Es))
-    s = hcat(Es, n1s, l1s, n2s, l2s)
-    for idx in CartesianIndices(mat)
+    # 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)
+    BRAC = 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},BRAC::Ptr{Array{Float64}})::Cvoid
+    
+    mat = zeros(basis.dim, basis.dim)
+
+    s = hcat(basis.Es, basis.n1s, basis.l1s, basis.n2s, basis.l2s)
+    Threads.@threads for idx in CartesianIndices(mat)
         (i, j) = Tuple(idx)
         (Elhs, N, L, n, l) = s[i, :]
         (Erhs, n1, l1, n2, l2) = s[j, :]
-        if Elhs == Erhs
-            phase = (-1)^(N + n + n1 + n2)
-            mat[i, j] = phase * get_bracket(BRACs[l1, l2], Emax, Λ, N, L, n, l, n1, n2)
+        if Elhs == Erhs && triangle_ineq(L, l, basis.Λ) && triangle_ineq(l1, l2, basis.Λ)
+            mat[i, j] = (-1)^(n1 + n2 + N + n) * pick_Moshinsky_bracket(BRAC, n1, l1, n2, l2, N, L, n, l, basis.Λ)
         end
     end
 
-    return mat
+    return sparse(mat)
 end
+
+function pick_Moshinsky_bracket(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
+
+function get_jacobi_V_matrix(V_of_r, basis::ho_basis_2B, μ1ω1, μω_global; atol=10^-6, maxevals=10^5)
+    V1 = get_jacobi_V1_matrix(V_of_r, basis, μ1ω1; atol=atol, maxevals=maxevals)
+    V2 = get_jacobi_V2_matrix(V_of_r, basis, μω_global; atol=atol, maxevals=maxevals)
+    return V1 + V2
+end
+
+function get_jacobi_V1_matrix(V_of_r, basis::ho_basis_2B, μ1ω1; atol=10^-6, maxevals=10^5)
+    V1_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μ1ω1, atol=atol, maxevals=maxevals)
+    V1 = get_sp_V_matrix(V1_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, E_max=basis.E_max)
+    return V1
+end
+
+function get_jacobi_V2_matrix(V_of_r, basis::ho_basis_2B, μω_global; atol=10^-6, maxevals=10^5)    
+    V_relative_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω_global, atol=atol, maxevals=maxevals)
+    V_relative_cache = prealloc_V_cache(basis.E_max, ComplexF64)
+
+    V_relative = get_sp_V_matrix(V_relative_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, cache=V_relative_cache) + get_sp_V_matrix(V_relative_elem, basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; dtype=ComplexF64, cache=V_relative_cache)
+    U = Moshinsky_transform(basis)
+    V2 =  transpose(U) * V_relative * U
+
+    return V2
+end
+
+function get_2p_p1p2_matrix(basis::ho_basis_2B, μ1ω1, μ2ω2; dtype=Float64)
+    # TODO: Cache for integrals
+    integral1(np, lp, n, l) = integral_HO(np, lp, n, l, μ1ω1)
+    integral2(np, lp, n, l) = integral_HO(np, lp, n, l, μ2ω2)
+
+    mat = zeros(dtype, basis.dim, basis.dim)
+    Threads.@threads for idx in CartesianIndices(mat)
+        (i, j) = Tuple(idx)
+        val = racahs_reduction_formula(basis.n1s[i], basis.l1s[i], basis.n2s[i], basis.l2s[i], basis.n1s[j], basis.l1s[j], basis.n2s[j], basis.l2s[j], basis.Λ, integral1, integral2)
+        if !(val ≈ 0); mat[idx] = val; end
+    end
+    return sparse(mat)
+end
+
+function get_src_V_matrix(V_of_r, basis::ho_basis_2B, μω, μω_global; atol=10^-6, maxevals=10^5)    
+    V_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω, atol=atol, maxevals=maxevals)
+    V_cache = prealloc_V_cache(basis.E_max, ComplexF64)
+
+    V1 = get_sp_V_matrix(V_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, cache=V_cache)
+    V2 = get_sp_V_matrix(V_elem, basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; dtype=ComplexF64, cache=V_cache)
+
+    V12 = get_src_V12_matrix(V_of_r, basis, μω_global; atol=atol, maxevals=maxevals)
+
+    return V1 + V2 + V12
+end
+
+function get_src_V12_matrix(V_of_r, basis::ho_basis_2B, μω_global; atol=10^-6, maxevals=10^5)
+    V_relative_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω_global, atol=atol, maxevals=maxevals)
+    V_relative = get_sp_V_matrix(V_relative_elem, basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; dtype=ComplexF64, E_max=basis.E_max)
+
+    U = Moshinsky_transform(basis)
+    V12 =  transpose(U) * V_relative * U
+
+    return V12
+end
+
+"Basis transformation from HO to momentum space"
+function get_W_matrix(basis_p, basis::ho_basis_2B, μ1ω1, μ2ω2=μ1ω1; weights=true)
+    W = zeros(ComplexF64, length(basis_p), basis.dim)
+    Threads.@threads for idx in CartesianIndices(W)
+        (i1, i2) = Tuple(idx)
+        ((j1, j2), (k1, w1), (k2, w2)) = basis_p[i1]
+        if j1 == basis.l1s[i2] && j2 == basis.l2s[i2]
+            elem1 = 1/sqrt(sqrt(μ1ω1)) * (-1)^basis.n1s[i2] * ho_basis(j1, basis.n1s[i2], 1/sqrt(μ1ω1) * k1)
+            elem2 = 1/sqrt(sqrt(μ2ω2)) * (-1)^basis.n2s[i2] * ho_basis(j2, basis.n2s[i2], 1/sqrt(μ2ω2) * k2)
+            W[idx] = elem1 * elem2 * (weights ? w1 * w2 : 1)
+        end
+    end
+    return sparse(W)
+end
+
+function get_3b_H_matrix(coord_system::coordinate_system, V_of_r, μω_global, E_max, Λ, m=1.0, kinetic_part=true, potential_part=true; atol=10^-5, maxevals=10^5, verbose=true)
+    if coord_system == jacobi
+        μ1ω1 = μω_global * 1/2
+        μ2ω2 = μω_global * 2
+        μ1 = m * 1/2
+        μ2 = m * 2/3        
+    elseif coord_system == src
+        μ1ω1 = μ2ω2 = μω = μω_global * 2
+        μ1 = μ2 = μ = m/2
+    end
+
+    verbose && println("No of threads = ", Threads.nthreads())
+
+    @time "Basis" basis = ho_basis_2B(E_max, Λ)
+    verbose && println("Basis size = ", basis.dim)
+
+    out = spzeros(basis.dim, basis.dim)
+
+    if kinetic_part
+        verbose && println("Constructing KE matrices")
+        @time "T1" out += get_sp_T_matrix(basis.n1s, basis.l1s, [basis.n2s, basis.l2s]; μω_gen=μ1ω1, μ=μ1)
+        @time "T2" out += get_sp_T_matrix(basis.n2s, basis.l2s, [basis.n1s, basis.l1s]; μω_gen=μ2ω2, μ=μ2)
+        if coord_system == src
+            @time "T_cross" out += get_2p_p1p2_matrix(basis, μ1ω1, μ2ω2; dtype=ComplexF64) ./ (2*μ)
+        end
+    end
+
+    if potential_part
+        verbose && println("Constructing PE matrices")
+        if coord_system == jacobi
+            @time "V" out += get_jacobi_V_matrix(V_of_r, basis, μ1ω1, μω_global; atol=atol, maxevals=maxevals)
+        elseif coord_system == src
+            @time "V" out += get_src_V_matrix(V_of_r, basis, μω, μω_global; atol=atol, maxevals=maxevals)
+        end
+    end
+
+    return out
+end

ho_basis_2body.jl

@@ -1,31 +0,0 @@
-using Arpack, SparseArrays
-include("ho_basis.jl")
-include("p_space.jl")
-
-E_max = 30
-ω = 0.2
-Λ = 0
-m = 1.0
-
-Va = -2
-Ra = 2
-
-μ1 = m * 1/2
-
-println("No of threads = ", Threads.nthreads())
-
-Es, n1s, l1s = get_sp_basis(E_max)
-println("Basis size = ", length(Es))
-
-println("Constructing KE matrices")
-@time "T1" T1 = sp_T_matrix(n1s, l1s; ω=ω, μ=μ1)
-
-println("Constructing PE matrices")
-V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; ω=ω)
-@time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s)
-
-println("Calculating spectrum")
-@time "H" H = T1 + V1
-@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=30, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
-
-display(evals)

ho_basis_3body.jl

@@ -1,43 +1,14 @@
-using Arpack, SparseArrays
+using Arpack
 include("ho_basis.jl")
-include("p_space.jl")
 
-E_max = 20
-ω = 0.2
 Λ = 0
 m = 1.0
+V_of_r(r) = -2 * exp(-r^2 / 4)
 
-Va = -2
-Ra = 2
+E_max = 40
+μω_global = 0.3
 
-μ1 = m * 1/2
-μ2 = m * 2/3 * 4
-c = sqrt(2)
+H = get_3b_H_matrix(src, V_of_r, μω_global, E_max, Λ, m)
 
-println("No of threads = ", Threads.nthreads())
-
-@time "Basis" Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max)
-@time "Masks" begin
-    mask1 = (n2s .== n2s') .&& (l2s .== l2s')
-    mask2 = (n1s .== n1s') .&& (l1s .== l1s')
-end
-
-println("Basis size = ", length(Es))
-
-println("Constructing KE matrices")
-@time "T1" T1 = sp_T_matrix(n1s, l1s; mask=mask1, ω=ω, μ=μ1)
-@time "T2" T2 = sp_T_matrix(n2s, l2s; mask=mask2, ω=ω, μ=μ2)
-
-println("Constructing PE matrices")
-V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; ω=ω)
-V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra / c, l, n1, n2; ω=ω)
-@time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s; mask=mask1)
-@time "V relative" V_relative = sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1) + sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2)
-@time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
-@time "V2" V2 =  U' * V_relative * U
-
-println("Calculating spectrum")
-@time "H" H = T1 + T2 + V1 + V2
 @time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=30, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
-
 display(evals)

ho_basis_3body_resonance.jl

@@ -0,0 +1,15 @@
+using Arpack
+include("ho_basis.jl")
+
+target_E = 4.07656088827514 - 0.012743522750966718im
+V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
+Λ = 0
+m = 1.0
+
+μω_global = 0.5 * exp(-2im * pi / 9)
+E_max = 30
+
+H = get_3b_H_matrix(src, V_of_r, μω_global, E_max, Λ, m)
+
+@time "Eigenvalues" evals, _ = eigs(H, nev=5, ncv=50, sigma=target_E, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
+display(evals)

ho_basis_benchmark.jl

@@ -2,7 +2,7 @@ using LinearAlgebra, Plots
 include("ho_basis.jl")
 include("p_space.jl")
 
-ω = 0.5 * exp(-1im * 0.45 * pi)
+μω_gen = 0.5 * exp(-1im * 0.45 * pi)
 μ = 0.5
 l = 0
 V1 = -5
@@ -14,10 +14,10 @@ n_max = 15
 ns = collect(0:n_max)
 ls = fill(l, n_max + 1)
 
-T = sp_T_matrix(ns, ls; ω=ω, μ=μ)
+T = get_sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ)
 
-V_l(l, n1, n2) = V1 * V_Gaussian(R1, l, n1, n2; ω=ω) + V2 * V_Gaussian(R2, l, n1, n2; ω=ω)
-V = sp_V_matrix(V_l, ns, ls; dtype=ComplexF64)
+V_l(l, n1, n2) = V1 * V_Gaussian(R1, l, n1, n2; μω_gen=μω_gen) + V2 * V_Gaussian(R2, l, n1, n2; μω_gen=μω_gen)
+V = get_sp_V_matrix(V_l, ns, ls; dtype=ComplexF64)
 
 cs = range(1.25, 0.25, 10)
 

math.jl

@@ -0,0 +1,47 @@
+using SpecialFunctions, WignerSymbols
+include("common.jl")
+
+# Gaussian potentials in HO space
+inv_factorial(n) = Iterators.prod(inv.(1:n))
+sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
+N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + l + 3/2)
+Talmi(l, R, k1, k2; μω_gen=1.0) = (-1)^(k1 + k2) * (1 + 1/(μω_gen * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
+V_Gaussian(R, l, n1, n2; μω_gen=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; μω_gen=μω_gen) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
+
+# for numerical evaluation of V matrix elements
+sqrt_double_factorial(n) = Iterators.prod(sqrt.(n:-2:1))
+sqrt_sqrt_pi = sqrt(sqrt(pi))
+laguerre(l, n, x) = gamma(n + l + 3/2) * better_sum([(-x * x)^k / gamma(k + l + 3/2) * inv_factorial(n - k) * inv_factorial(k)  for k in 0:n])
+ho_basis_const(l, n) = (-1)^n / sqrt_sqrt_pi * (2.0)^((n + l + 2) / 2) * sqrt_factorial(n) / sqrt_double_factorial(2*n + 2*l + 1)
+ho_basis_func(l, n, x) = x^(l + 1) * exp(-x^2 / 2) * laguerre(l, n, x)
+ho_basis(l, n, x) = ho_basis_const(l, n) * ho_basis_func(l, n, x)
+
+# for implementation of simple relative coordinates
+double_factorial(n::Int) = Iterators.prod(big, n:-2:1)
+
+"Gaussian integral for n ∈ Integers (Ref: Wolfram MathWorld + simplifications)"
+gauss_int(a, n) = double_factorial(n - 1) / (2.0 * a)^((n + 1)/2) * (iseven(n) ? sqrt(π / 2) : 1)
+
+"Gives ∫dp p u' u where u' and u are HO functions with different l (Ref: worked out in Mathematica)"
+function integral_HO(np, lp, n, l, μω)
+    s = [(-1)^(m + mp) * gauss_int(1, 2 * m + 2 * mp + l + lp + 3) * N_lnk(l, n, m) * N_lnk(lp, np, mp) for (m, mp) in Iterators.product(0:n, 0:np)]
+    return 2 * sqrt(μω) * better_sum(s)
+end
+
+"Gives <n' l'|| p ||n l> for the HO basis, where integral(np, lp, n, l) is a function that returns ∫dp p u' u"
+function reduced_matrix_element(np, lp, n, l, integral)::ComplexF64
+    wig::Float64 = wigner3j(Float64, lp, 1, l, 0, 0, 0)
+    if wig == 0
+        return 0
+    else
+        return (-1)^lp * sqrt(2*lp + 1) * sqrt(2*l + 1) * wig * integral(np, lp, n, l)
+    end
+end
+
+"Matrix element <n1p l1p n2p l2p| p1⋅p2 |n1 l1 n2 l2> (Ref: de-Shalit & Talmi, Eq 15.5), where integral(np, lp, n, l) is a function that returns ∫dp p u' u"
+function racahs_reduction_formula(n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, integral1, integral2)
+    val = wigner6j(Float64, l1p, l2p, Λ, l2, l1, 1)
+    if val != 0; val *= reduced_matrix_element(n1p, l1p, n1, l1, integral1); end
+    if val != 0; val *= reduced_matrix_element(n2p, l2p, n2, l2, integral2); end
+    return (-1)^(l1 + l2p + Λ) * val
+end

osbrackets.jl

@@ -1,42 +0,0 @@
-# compile using
-# gfortran -shared -fPIC osbrac.f90 -o osbrac.so
-# gfortran -shared -fPIC allosbrac.f90 -o allosbrac.so
-
-using Libdl
-
-function cal_BRAC(Emax, Λ, l1, l2)
-    ϵ = (2 * Emax - Λ) % 2
-    N = (l1 - l2 + Λ - ϵ) ÷ 2
-    M = (l1 + l2 - Λ - ϵ) ÷ 2
-    L = Λ
-    NQMAX = 2 * Emax
-    CO = 1/sqrt(2)
-    SI = 1/sqrt(2)
-    FIRSTCALL = true
-    # dimensions BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2)
-    BRAC = zeros(Float64, 1 + L, 1 + (NQMAX - L) ÷ 2, 1 + (NQMAX - L) ÷ 2, 1 + (NQMAX - L) ÷ 2, 1 + (NQMAX - L) ÷ 2)
-
-    lib = Libdl.dlopen("OSBRACKETS/osbrac.so") # Open the library explicitly.
-    sym = Libdl.dlsym(lib, :osbrac_)   # Get a symbol for the function to call.
-    # call signature OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
-    @ccall $sym(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
-    Libd.dlclose(lib)
-
-    return BRAC
-end
-
-function get_bracket(BRAC, Emax, Λ, n1′, l1′, n2′, l2′, n1, n2)
-    ϵ = (2 * Emax - Λ) % 2
-    NP = (l1′ - l2′ + Λ - ϵ) ÷ 2
-    N1P = n1′
-    MP = (l1′ + l2′ - Λ - ϵ) ÷ 2
-    N1 = n1
-    N2 = n2
-    N2P = n2′
-    if MP+N1P+N2P == M+N1+N2
-        # BRAC(NP,N1P,MP,N1,N2)
-        return BRAC[1 + NP, 1 + N1P, 1 + MP, 1 + N1, 1 + N2]
-    else
-        return 0
-    end
-end

p_space.jl

@@ -1,8 +1,6 @@
-using FastGaussQuadrature
-
-# Gaussian potentials in momentum space
-g0(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(-1 + exp(p*q*R^2))*R)/(2*sqrt(π))
-g1(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(2 + p*q*R^2 + exp(p*q*R^2)*(-2 + p*q*R^2)))/(2*p*sqrt(π)*q*R)
+using LinearAlgebra
+using SpecialFunctions, FastGaussQuadrature, QuadGK
+include("ho_basis.jl")
 
 function gausslegendre_shifted(a, b, n)
     scale = (b - a) / 2
@@ -13,11 +11,11 @@ function gausslegendre_shifted(a, b, n)
     return (p, w)
 end
 
-function get_mesh(vertices, subdivisions)
+function get_mesh(vertices::Vector, subdivs::Vector)
     p = Vector{ComplexF64}()
     w = Vector{ComplexF64}()
-    for (a, b) in zip(vertices, vertices[2:end])
-        p_new, w_new = gausslegendre_shifted(a, b, subdivisions)
+    for (a, b, subdiv) in zip(vertices, vertices[2:end], subdivs)
+        p_new, w_new = gausslegendre_shifted(a, b, subdiv)
         append!(p, p_new)
         append!(w, w_new)
     end
@@ -46,7 +44,72 @@ function identify_pole_i(p, evals, μ=0.5)
 end
 
 function quick_pole_E(V_pq, μ=0.5; cs_angle=0.4, cutoff=8.0, meshpoints=256)
-    p, w = get_mesh([0, cutoff * exp(-1im * cs_angle)], meshpoints)
+    p, w = get_mesh([0, cutoff * exp(-1im * cs_angle)], [meshpoints])
     evals = eigvals(get_H_matrix(V_pq, p, w, μ))
     return evals[identify_pole_i(p, evals, μ)]
+end
+
+# Gaussian potentials in momentum space
+g0(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(-1 + exp(p*q*R^2))*R)/(2*sqrt(π))
+g1(R, p, q) = (exp(-(1/4)*(p + q)^2*R^2)*(2 + p*q*R^2 + exp(p*q*R^2)*(-2 + p*q*R^2)))/(2*p*sqrt(π)*q*R)
+
+# general potential (numerical integration)
+jHat(l, z) = z * sphericalbesselj(l, z)
+function Vl_mat_elem(V_of_r, l, p, q; atol=0, maxevals=10^7, R_cutoff=Inf)
+    integrand(r) = jHat(l, p * r) * V_of_r(r) * jHat(l, q * r)
+    (integral, _) = quadgk(integrand, 0, R_cutoff; atol=atol, maxevals=maxevals)
+    return (2 / pi) * integral
+end
+
+"Return the Hamiltonian matrix (and the array of weights) for the given system."
+function get_3b_H_matrix(coord_system::coordinate_system, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m=1.0, kinetic_part=true, potential_part=true; atol=10^-5, maxevals=10^5, R_cutoff=16, verbose=true)
+    if coord_system == jacobi
+        μ1ω1 = μω_global * 1/2
+        μ2ω2 = μω_global * 2
+        μ1 = m * 1/2
+        μ2 = m * 2/3
+    else
+        error("Only Jacobi coordinates are implemented")
+    end
+
+    verbose && println("No of threads = ", Threads.nthreads())
+
+    V_l(j, k, kp) = Vl_mat_elem(V_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
+
+    ks, ws = get_mesh(vertices, subdivisions)
+    weights = repeat(kron(ws, ws), jmax + 1)
+    block_size = length(ks)
+    tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
+    js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
+    basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
+    basis_size = length(js) * length(ks)^2
+    @assert length(basis) == basis_size "Something wrong with the basis"
+    verbose && println("Basis size = $basis_size")
+
+    out = spzeros(basis_size, basis_size)
+
+    @time "Block diagonal part" begin
+        if kinetic_part & potential_part
+            Hb_blocks = [kron_sum(get_H_matrix((k, kp) -> V_l(j1, k, kp), ks, ws, μ1), get_T_matrix(ks, μ2)) for (j1, _) in js]
+        elseif kinetic_part
+            Hb_blocks = [kron_sum(get_T_matrix(ks, μ1), get_T_matrix(ks, μ2)) for _ in js]
+        elseif potential_part
+            Hb_blocks = [kron_sum(get_V_matrix((k, kp) -> V_l(j1, k, kp), ks, ws), spzeros(block_size, block_size)) for (j1, _) in js]
+        end
+        out += blockdiag(sparse.(Hb_blocks)...)
+    end
+
+    if potential_part
+        basis_ho = ho_basis_2B(E_max, Λ)
+        verbose && println("HO basis size = ", basis_ho.dim)
+    
+        @time "V2_HO" V2_HO =  get_jacobi_V2_matrix(V_of_r, basis_ho, μω_global)
+
+        @time "W_right" W_right = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=true)
+        @time "W_left" W_left = get_W_matrix(basis, basis_ho, μ1ω1, μ2ω2; weights=false)
+
+        @time "V2" out += W_left * V2_HO * transpose(W_right)
+    end
+
+    return (out, weights)
 end

p_space_2body.jl

@@ -0,0 +1,25 @@
+using SparseArrays, Arpack
+include("common.jl")
+include("p_space.jl")
+
+E_target = -0.3919
+
+μ = 0.5
+Va = -2
+Ra = 2
+V_of_r(r) = Va * exp(-r^2 / Ra^2)
+
+vertices = [0, 0.5 - 0.3im, 1, 4]
+subdivisions = [16, 16, 16]
+ks, ws = get_mesh(vertices, subdivisions)
+ls = collect(0:4)
+
+V_l(l, k, kp) = Vl_mat_elem(V_of_r, l, k, kp; atol=10^-5, maxevals=10^5, R_cutoff=16)
+
+H_l = [get_H_matrix((k, kp) -> V_l(l, k, kp), ks, ws, μ) for l in ls]
+
+H1 = blockdiag(sparse.(H_l)...)
+
+H = H1
+@time "Eigenvalues" evals, _ = eigs(H, nev=10, ncv=50, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
+E = nearest(evals, E_target)

p_space_3body.jl

@@ -0,0 +1,18 @@
+using Arpack
+include("p_space.jl")
+
+Λ = 0
+m = 1.0
+V_of_r(r) = -2 * exp(-r^2 / 4)
+
+vertices = [0, 0.5 - 0.3im, 1, 4]
+subdivisions = [10, 10, 10]
+jmax = 4
+
+E_max = 30
+μω_global = 0.5
+
+H, _ = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m)
+
+@time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
+display(evals)

p_space_3body_resonance.jl

@@ -0,0 +1,20 @@
+using Arpack
+include("p_space.jl")
+
+target = 4.0766890719636875 - 0.012758927741074495im
+
+Λ = 0
+m = 1.0
+V_of_r(r) =  2 * exp(-(r-3)^2 / (1.5)^2)
+
+vertices = [0, 2 - 0.2im, 3, 4]
+subdivisions = [15, 10, 10]
+jmax = 4
+
+E_max = 40
+μω_global = 0.5
+
+H, _ = get_3b_H_matrix(jacobi, V_of_r, vertices, subdivisions, jmax, μω_global, E_max, Λ, m)
+
+@time "Eigenvalues" evals, _ = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
+display(evals)

p_space_test.ipynb

@@ -1,61 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "using Plots, LinearAlgebra\n",
-    "include(\"p_space.jl\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "vertices = (0, 1 - 0.5im, 2, 6)\n",
-    "subdivisions = 64\n",
-    "p, w = get_mesh(vertices, subdivisions)\n",
-    "\n",
-    "scatter(real.(p), imag.(p))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# ComplexScaling-FV: Eq. (54)\n",
-    "V_pq(p, q) = -10 * g1(1, p, q)\n",
-    "\n",
-    "H = get_H_matrix(V_pq, p, w)\n",
-    "evals = eigen(H).values\n",
-    "\n",
-    "E_target = 0.258 - 0.164im\n",
-    "E = evals[argmin(norm.(evals .- E_target))]\n",
-    "\n",
-    "print(\"E = $E\")\n",
-    "scatter(real.(evals), imag.(evals), xlim = (0,1))"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Julia 1.9.0",
-   "language": "julia",
-   "name": "julia-1.9"
-  },
-  "language_info": {
-   "file_extension": ".jl",
-   "mimetype": "application/julia",
-   "name": "julia",
-   "version": "1.9.0"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}

simple_berggren.ipynb

@@ -1,69 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "using Plots, LinearAlgebra\n",
-    "include(\"p_space.jl\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "vertices = (0, 0.5 - 0.3im, 1, 6)\n",
-    "subdivisions = 64\n",
-    "p, w = get_mesh(vertices, subdivisions)\n",
-    "\n",
-    "# resonance example from my thesis\n",
-    "V_basis(p, q) = 2*g0(4, p, q) - 3*g0(2, p, q)\n",
-    "\n",
-    "basis_eig = eigen(get_H_matrix(V_basis, p, w))\n",
-    "basis = basis_eig.vectors .* w\n",
-    "\n",
-    "evals = basis_eig.values\n",
-    "E_target = 0.7\n",
-    "E = evals[argmin(norm.(evals .- E_target))]\n",
-    "\n",
-    "print(\"E = $E\")\n",
-    "scatter(real.(evals), imag.(evals), xlim = (0,2))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# ResonanceEC: Eq. (20)\n",
-    "V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))\n",
-    "\n",
-    "H = get_H_matrix(V_system(0.45), p, w)\n",
-    "H_berggren = transpose(basis) * H * basis\n",
-    "\n",
-    "evals = eigvals(H)\n",
-    "scatter(real.(evals), imag.(evals), xlim = (0, 0.5))"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "Julia 1.9.0",
-   "language": "julia",
-   "name": "julia-1.9"
-  },
-  "language_info": {
-   "file_extension": ".jl",
-   "mimetype": "application/julia",
-   "name": "julia",
-   "version": "1.9.0"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}

test/3body_orthogonality.jl

@@ -0,0 +1,8 @@
+include("../p_space_3body_resonance.jl")
+
+@time "Eigenvectors" evals, evecs = eigs(H, sigma=target, maxiter=5000, tol=1e-5, ritzvec=true, check=1)
+
+weights_mat = spdiagm(repeat(kron(ws, ws), jmax + 1))
+N = transpose(evecs) * weights_mat * evecs
+
+display(abs.(N))

test/OSBRACKETS.jl

@@ -2,27 +2,35 @@
 # gfortran -shared -fPIC osbrac.f90 -o osbrac.so
 # gfortran -shared -fPIC allosbrac.f90 -o allosbrac.so
 
-println("OSBRACKETS test against Brody et al.")
+println("OSBRACKETS test against Buck et al.")
 
-Emax = 3
-Λ = 1
-l1 = 1
-l2 = 2
+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];
 
-ϵ = (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)
+function calculate_single_bracket(n1′, l1′, n2′, l2′, n1, l1, n2, l2, Λ) # Efros notation -- DON'T CONFUSE
 
-@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
+    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
 
-function get_bracket(n1′, l1′, n2′, l2′, n1, n2)
     Nq = l1 + l2 + 2 * (n1 + n2)
     Nq′ = l1′ + l2′ + 2 * (n1′ + n2′)
     if Nq ≠ Nq′; return 0; end
@@ -31,13 +39,55 @@ function get_bracket(n1′, l1′, n2′, l2′, n1, n2)
     M′ = (l1′ + l2′ - Λ - ϵ) ÷ 2
     N1 = n1
     N2 = n2
+
     return BRAC[N′ + 1, N1′ + 1, M′ + 1, N1 + 1, N2 + 1]
 end
 
-test_n = [0 0 0 0 1 1]
-test_l = [0 1 1 2 0 1]
-test_N = [1 0 1 0 0 0]
-test_L = [1 2 0 1 1 0]
+println("OSBRAC results")
+osbracs = calculate_single_bracket.(n1, l1, n2, l2, N, L, n, l, Λ)
+display(osbracs)
 
-bracs = get_bracket.(test_n, test_l, test_N, test_L, 0, 0)
-display(bracs)
+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))

test/buck_et_al.jl

@@ -0,0 +1,25 @@
+using NuclearToolkit
+
+# Table 1 from 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];
+
+Elhs = l1 .+ l2 .+ 2 .* (n1 .+ n2)
+Erhs = l .+ L .+ 2 .* (n .+ N)
+E_max = max(maximum(Elhs), maximum(Erhs))
+j_max = 2 * E_max + 1
+l_max = j_max
+
+dtri = NuclearToolkit.prep_dtri(l_max + 1);
+dcgm0 = NuclearToolkit.prep_dcgm0(l_max);
+d6j = nothing # NuclearToolkit.prep_d6j_int(E_max, j_max, to);
+
+bracs = NuclearToolkit.HObracket_d6j.(N, L, n, l, n1, l1, n2, l2, Λ, 1.0, Ref(dtri), Ref(dcgm0), Ref(d6j))
+display(bracs)

test/cosmo.jl

@@ -0,0 +1,19 @@
+# depends on a fork of alvolya/cosmo
+
+println("COSMO test against Brody et al.")
+
+get_bracket(n, l, N, L, n1, l1, n2, l2, lambda, d) = @ccall "../cosmo/bin/shared.so"._Z24MoshinskyBracketsWrapperiiiiiiiiid(n::Cint,l::Cint,N::Cint,L::Cint,n1::Cint,l1::Cint,n2::Cint,l2::Cint,lambda::Cint,d::Cdouble)::Cdouble
+
+# Table 1 from 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];
+
+bracs = get_bracket.(n, l, N, L, n1, l1, n2, l2, Λ, 1.0)
+display(bracs)

test/ho_basis_test.jl

@@ -1,5 +1,5 @@
 using LinearAlgebra
-include("ho_basis.jl")
+include("../ho_basis.jl")
 
 l = 0
 V0 = -10
@@ -9,10 +9,10 @@ n_max = 20
 ns = collect(0:n_max)
 ls = fill(l, n_max + 1)
 
-T = sp_T_matrix(ns, ls)
+T = get_sp_T_matrix(ns, ls)
 
 V_l(l, n1, n2) = V0 * V_Gaussian(R, l, n1, n2)
-V = sp_V_matrix(V_l, ns, ls)
+V = get_sp_V_matrix(V_l, ns, ls)
 
 H = T + V
 

test/misc.jl

@@ -0,0 +1,65 @@
+using SparseArrays, LinearAlgebra
+include("../common.jl")
+
+#### gram_schmidt ####
+ 
+n = 64
+N = 8
+
+vecs = [rand(n) + 1im .* rand(n) for _ in 1:N]
+
+ws = rand(n)
+ws_mat = spdiagm(ws)
+
+basis = hcat(vecs...)
+
+H = rand(n, n)
+H += transpose(H) # complex symmetric
+
+H_EC = transpose(basis) * ws_mat * H * ws_mat * basis
+N_EC = transpose(basis) * ws_mat * basis
+
+evals = eigvals(H_EC, N_EC)
+println("Eigenvalues with GEVP:")
+display(evals)
+
+ortho_basis = hcat(gram_schmidt!(vecs, ws)...)
+N_ortho = transpose(ortho_basis) * ws_mat * ortho_basis
+println("Norm matrix after Gram-Schmidt:")
+display(round.(N_ortho; digits=2)) # should be ≈I
+@assert N_ortho ≈ I(N) "Gram-Schmidt did not yield an orthogonal basis"
+
+H_EC_ortho = transpose(ortho_basis) * ws_mat * H * ws_mat * ortho_basis
+
+evals_ortho = eigvals(H_EC_ortho)
+println("Eigenvalues after Gram-Schmidt:")
+display(evals_ortho)
+
+@assert evals ≈ evals_ortho "Gram-Schmidt did not preserve the eigenvalues"
+
+############
+println("\nRepeat with a redundant basis\n")
+
+println("Original dimensionality = $(length(vecs))")
+
+for pow in 6:9
+    noise = rand(n) ./ 10^pow
+    new_vec = vecs[1] + noise
+    push!(vecs, new_vec)
+end
+
+println("Dimensionality before Gram-Schmidt = $(length(vecs))")
+
+ortho_vecs = gram_schmidt!(vecs, ws; verbose=true)
+ortho_basis = hcat(ortho_vecs...)
+println("Dimensionality after Gram-Schmidt = $(length(ortho_vecs))")
+
+H_EC_ortho = transpose(ortho_basis) * ws_mat * H * ws_mat * ortho_basis
+
+evals_ortho = eigvals(H_EC_ortho)
+println("Eigenvalues after Gram-Schmidt:")
+display(evals_ortho)
+
+@assert isapprox(evals, evals_ortho; atol=1e-3) "Gram-Schmidt did not approximately preserve the eigenvalues"
+
+######################

test/moshinsky.jl

@@ -1,21 +1,21 @@
-println("### Test: U' * U == identity")
+println("### Test: transpose(U) * U == identity")
 
 using LinearAlgebra
 
 include("../ho_basis.jl")
 
-E_max = 15
+E_max = 30
 Λ = 0
 
 println("No of threads = ", Threads.nthreads())
 
-@time "Basis" Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max)
+@time "Basis" basis = ho_basis_2B(E_max, Λ)
 
-println("Basis size = ", length(Es))
+println("Basis size = ", basis.dim)
 
-@time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
+@time "Moshinsky brackets" U = Moshinsky_transform(basis)
 
-check = conj(U) * U - I
+check = transpose(U) * U - I
 
 println("Maximum = ", maximum(abs.(check)))
 println("Norm = ", sum(check .* conj(check)))

test/numerical_V.jl

@@ -0,0 +1,20 @@
+using Statistics
+
+include("../ho_basis.jl")
+
+ls = 0:4
+ns = 0:15
+
+R = 2
+μω_gen = 0.3
+V_of_r(r) = exp(-r^2 / R^2)
+
+for l in ls
+    println("Testing l=", l)
+
+    @time "Analytical" V_ana = [V_Gaussian(R, l, n1, n2; μω_gen=μω_gen) for (n1, n2) in Iterators.product(ns, ns)]
+    @time "Numerical" V_num = [V_numerical(V_of_r, l, n1, n2; μω_gen=μω_gen) for (n1, n2) in Iterators.product(ns, ns)]
+
+    mean_diff = Statistics.mean(abs.((V_num .- V_ana) ./ V_ana))
+    @assert mean_diff < 1e-5 "Mean absoute difference of $(mean_diff * 100)% between analytical and numerical calculations"
+end

test/p1p2_matrix.jl

@@ -0,0 +1,57 @@
+using LinearAlgebra, SparseArrays, Plots
+include("../p_space.jl")
+include("../ho_basis.jl")
+include("../berggren.jl")
+
+println("No of threads = ", Threads.nthreads())
+atol = 10^-5
+maxevals = 10^5
+R_cutoff = 16
+
+Λ = 0
+m = 1.0
+μ = m/2 # due to simple relative coordinates
+
+vertices = [0, 2 - 0.2im, 3, 4]
+subdivisions = [15, 10, 10]
+ks, ws = get_mesh(vertices, subdivisions)
+
+jmax = 4
+tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
+js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
+
+basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
+basis_size = length(js) * length(ks)^2
+@assert length(basis) == basis_size "Something wrong with the basis"
+println("Basis size = $basis_size")
+
+V_of_r(r) =  2 * exp(-(r-3)^2 / (1.5)^2)
+V_l(j, k, kp) = Vl_mat_elem(V_of_r, j, k, kp; atol=atol, maxevals=maxevals, R_cutoff=R_cutoff)
+
+# generate p-space bases (actually identity matrices)
+@time "p-space bases" ps_bases = [Matrix(spdiagm(1 ./ sqrt.(ws))) for _ in 0:jmax]
+
+# generate Berggren bases
+@time "Berggren bases" begin
+    berg_bases = Vector{Matrix{ComplexF64}}(undef, jmax + 1)
+    for j in 0:jmax
+        _, berg_basis = eigen(get_H_matrix((k, kp) -> V_l(j, k, kp), ks, ws, μ); permute=false, scale=false)
+        N_berg = sum(berg_basis.^2 .* ws, dims=1)
+        berg_bases[1 + j] = berg_basis ./ transpose(sqrt.(N_berg))
+    end
+end
+
+@time "BB tranform matrix" begin
+    U_blocks = [kron(berg_bases[1 + j1] .* ws, berg_bases[1 + j2] .* ws) for (j1, j2) in js]
+    U = blockdiag(sparse.(U_blocks)...)
+end
+
+@time "In p-space" T_cross_PS = get_2p_p1p2_matrix(length(ks), js, Λ, ps_bases, ps_bases, ws)
+@time "In BB" T_cross_BB = get_2p_p1p2_matrix(length(ks), js, Λ, berg_bases, berg_bases, ws)
+
+@time "Basis transform" T_cross_transformed = transpose(U) * T_cross_PS * U
+
+diff = abs.(T_cross_transformed .- T_cross_BB)
+println("Max error = $(maximum(diff))")
+
+#@time "In in HO" T_cross_HO = get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μω, μω; dtype=ComplexF64)

test/p_space_basic.jl

@@ -0,0 +1,21 @@
+using Plots
+include("../common.jl")
+include("../p_space.jl")
+
+vertices = [0, 1 - 0.5im, 2, 6]
+subdivisions = [64, 64, 64]
+p, w = get_mesh(vertices, subdivisions)
+
+scatter(real.(p), imag.(p))
+
+# ComplexScaling-FV: Eq. (54)
+V_pq(p, q) = -10 * g1(1, p, q)
+
+H = get_H_matrix(V_pq, p, w)
+evals = eigen(H).values
+
+E_target = 0.258 - 0.164im
+E = nearest(evals, E_target)
+
+print("E = $E")
+scatter(real.(evals), imag.(evals), xlim = (0,1))

test/racah_test_data.csv

@@ -0,0 +1,384 @@
+2,0,2,3,7,3,5,1,3,2.2274833794413214,0.45049378982090493,0.
+8,1,3,3,2,2,0,0,2,2.202162451867089,2.3467827533568713,0.
+4,1,7,1,0,2,0,3,2,1.1507774013957701,1.7553231930689552,0.
+3,3,7,1,7,0,3,3,3,1.7207018027964498,0.5156168080476098,0.
+10,3,2,3,2,2,8,1,2,1.8402526798628767,2.8981293881514993,0.
+4,3,8,1,0,3,0,3,3,2.471579266664868,0.9515272650259292,0.
+9,2,8,3,10,3,1,2,1,2.1675261170354094,0.38546777353112605,0.
+2,2,7,3,4,3,7,3,3,2.76860332907484,0.7916486882071214,0.
+2,3,3,3,5,1,1,0,1,1.6146747022554235,2.8377955838430013,0.
+4,3,5,0,7,3,2,3,3,0.9198951336314427,0.9950080709108984,0.
+7,3,4,2,2,2,8,0,2,0.2364795732399818,1.839193911487496,0.
+7,3,3,3,5,1,6,2,3,0.29514953289458434,0.3780481641738338,0.
+5,2,4,1,9,3,1,0,3,2.8086957499249543,2.1260743936378903,0.
+4,1,6,1,1,0,10,0,0,2.6628402589565523,2.425950546585856,0.
+8,3,1,0,6,3,0,3,3,2.1328726647005265,0.9431361454676961,0.
+4,3,5,0,5,1,0,3,3,2.576495013458395,0.37293785437117544,0.
+10,3,4,0,1,2,2,3,3,2.0030278179077214,1.1187270744696427,0.
+9,1,8,2,1,1,5,3,2,2.789813153068671,2.6763038014457665,0.
+5,0,6,1,5,3,9,3,1,1.9319197502299517,0.6216330166798851,0.
+8,0,6,2,7,3,8,2,2,1.8084358691372344,2.0181750785519634,0.
+3,1,5,1,9,0,7,1,1,1.599191638971793,1.552252944957802,0.
+1,0,7,2,6,3,4,2,2,2.4217191310288078,1.2381047980553532,0.
+9,3,2,0,1,2,4,2,3,1.261927284387386,0.39708064622835604,0.
+0,1,5,1,5,0,0,2,2,0.407989632316895,2.4021451804173477,0.
+3,1,6,3,9,2,2,1,2,1.2053920704239278,0.25375827636999704,0.
+9,2,7,0,5,2,5,3,2,2.539812963473447,0.7382167735661622,0.
+10,2,7,3,4,3,0,0,3,0.8427632747449358,0.7622120453149561,0.
+5,2,0,2,7,3,5,0,3,0.460843656858851,0.583650962079461,0.
+8,1,8,0,6,2,0,2,1,2.056601089239093,0.22262414697020771,0.
+4,2,4,2,6,1,0,3,3,1.4808637079563942,1.358717819227813,0.
+10,3,4,3,4,0,5,2,2,1.5284998408295385,2.925365977036159,0.
+10,1,1,1,5,1,10,1,0,1.1846342029623314,1.186468360104492,0.
+9,0,5,2,1,2,2,3,2,1.733114516642603,2.8283467464705456,0.
+8,3,8,1,10,3,4,3,3,0.8918790824078111,0.9129770894549143,0.
+10,0,9,3,5,3,8,3,3,0.6417078137732939,1.965954991753998,0.
+10,2,8,2,8,2,8,3,4,0.6600917142427227,0.987969186819146,0.
+8,2,4,1,3,2,9,1,3,0.8935573241720247,2.1162172596314806,0.
+4,2,10,3,7,2,7,2,2,1.5052737158592535,2.7171501665361504,0.
+0,1,3,1,5,2,5,0,2,1.577080591236359,0.34164973402053134,0.
+0,2,10,3,6,3,6,0,3,2.7925439641102736,1.831313118698632,0.
+1,0,9,3,9,0,0,3,3,2.3331462089933934,1.6945470525813482,0.
+8,2,7,3,7,2,3,1,1,1.4988059394628408,1.3530552540483156,0.
+9,2,8,1,8,3,6,1,2,0.6337531145346973,1.5268564417495032,0.
+10,1,0,0,7,2,7,2,1,1.3790583258395,1.5669277734464737,0.
+3,0,0,2,6,0,4,2,2,2.4730048216795737,2.70759381187844,0.
+8,1,3,1,5,3,8,2,2,0.8399302875607697,0.32436518041750384,0.
+4,1,1,2,4,3,10,1,3,0.8437127682376069,1.7940889253407324,0.
+9,0,3,1,4,3,2,2,1,0.8755496194878809,0.9675651711989977,0.
+9,0,2,2,3,1,8,2,2,1.4950146011492311,2.8152388609955157,0.
+7,0,8,3,2,3,1,3,3,2.2587603015484863,2.0666085413750848,0.
+10,3,9,1,4,3,0,2,4,1.5489892079797567,1.7151221413658782,0.
+3,0,1,2,2,3,4,1,2,0.25514014111567285,2.06281752109991,0.
+5,1,1,0,9,1,3,2,1,1.3329695370285548,0.5298097905264418,0.
+8,1,1,0,5,0,9,1,1,2.95254173854617,2.1687131876307753,0.
+6,1,8,0,5,2,1,2,1,1.5801597077900533,2.257296177275017,0.
+10,0,9,1,9,3,4,3,1,2.5352867964035397,1.7981404960338239,0.
+10,2,1,2,0,0,9,3,3,2.06463437201601,1.4295850747781111,0.
+1,0,1,3,6,3,6,3,3,0.8872777463742612,1.9002353924129451,0.
+6,3,10,0,2,3,2,1,3,1.2410196118280048,1.8568476669071874,0.
+8,2,1,2,2,3,6,3,2,0.5302021471052729,1.2216709678909177,0.
+1,2,8,1,0,3,3,0,3,2.8876975946711045,2.717322643102454,0.
+4,2,4,1,0,1,7,0,1,2.0026155982539864,0.7640140889966269,0.
+1,0,1,1,4,2,3,2,1,1.016582775534538,2.3841546013418045,0.
+2,0,3,1,10,1,1,1,1,1.6404532954231827,2.1972778830632356,0.
+7,3,1,3,2,2,4,3,4,2.648817905371118,2.798529790102262,0.
+9,0,4,1,6,0,2,1,1,2.3138662599243407,1.9450229178773055,0.
+0,1,8,1,1,1,1,0,1,0.6524221878214105,2.9034613327001972,0.
+6,1,6,3,0,1,6,2,3,1.831532710689828,2.890598375121809,0.
+4,1,0,0,0,2,7,3,1,0.626504505885948,2.6340691505629463,0.
+10,3,2,3,10,3,1,1,2,1.6591641837977482,2.452206213376977,0.
+3,1,6,2,2,2,7,1,3,1.1999874026225852,2.795586000355012,3.3573329331683617
+8,3,1,2,1,3,5,2,5,0.6830778416708307,1.1113967340568016,0.
+7,2,3,1,5,3,8,2,3,1.714181314975816,1.3518067184764888,0.
+0,2,4,2,0,1,4,1,2,1.060932248034836,0.8036155418186053,-1.1371377045290691
+4,3,9,3,10,3,10,3,5,0.42490400067847744,1.9191837942593741,0.
+0,0,6,1,2,2,7,2,1,1.4606282497051195,2.390102981429126,0.
+8,0,0,1,9,0,10,1,1,2.9779138627844803,1.482942802025069,0.
+1,2,2,3,3,3,8,3,4,1.540918746101882,0.9621051195430765,0.
+2,1,1,3,2,0,1,3,3,2.0584800719627525,1.7796724319652748,0.
+10,2,5,2,4,0,0,2,2,1.2759898818923578,0.6080747387544254,0.
+1,0,9,3,2,2,10,2,3,2.835152414163625,1.2470542223023808,0.
+1,2,2,0,6,0,7,2,2,1.8642054302207116,0.48490065380939873,0.
+5,1,10,2,5,0,10,3,3,0.6428303158043476,0.5501263433952537,2.105493662515726
+1,1,1,1,4,2,0,1,1,2.91586289712962,0.2030446107349677,0.
+4,2,0,2,1,3,1,0,3,0.31239590050870536,2.3530015462264577,0.
+10,0,0,3,0,1,4,3,3,1.4066735483371904,0.3492189542356794,0.
+3,1,9,3,6,2,4,2,4,2.039761320004641,2.0376049586912632,0.
+4,1,3,3,2,2,8,2,2,2.485347427066465,2.2425974670112305,0.
+8,0,10,3,7,3,3,0,3,1.145305638226675,0.7430070058758362,0.
+1,0,4,3,1,3,3,2,3,2.236729779510484,2.518871850291589,0.
+9,3,10,2,0,0,9,3,3,2.2231614378407896,2.934256387523143,0.
+1,3,7,1,6,2,7,1,3,1.739591417445606,2.852892536289553,0.
+9,0,0,2,7,1,0,2,2,2.628613403491599,1.948716181798507,0.
+2,0,5,3,3,2,4,3,3,0.7244831235767495,0.307815884940033,0.
+9,2,2,3,4,0,9,2,2,1.4712965684171477,0.7152116770614816,0.
+10,0,7,1,7,2,2,2,1,2.0944053463871253,0.7877074060643912,0.
+10,2,8,2,4,1,4,3,3,1.2400259243830605,1.533481685252715,0.
+7,2,2,0,6,3,2,1,2,0.7570369287424406,1.5727691945398687,2.4154000428935993
+9,0,9,2,4,2,5,0,2,2.9903195228568853,1.7342728120955062,0.
+0,2,2,1,2,2,2,2,3,2.871118456330721,1.8840932443901108,0.
+1,3,4,1,9,2,8,1,3,2.2820457592186356,2.8820550176739523,0.
+4,3,9,1,5,3,7,1,2,0.5157719636123437,2.2186070481213065,0.
+3,3,9,3,7,2,1,1,1,1.806277095079202,2.303389106635813,0.
+4,2,0,0,8,0,10,2,2,0.4100439145977339,1.0520725633546872,0.
+10,1,6,0,4,2,2,3,1,1.126129220884609,1.5550003535454513,0.
+3,2,10,1,6,1,10,1,1,0.2658540542443206,0.608227371386918,0.
+6,3,5,1,4,1,0,3,3,2.705056186004474,0.42463741187699755,0.
+10,0,5,2,5,1,1,2,2,0.8756297484714799,0.7584762754047034,0.
+7,1,7,3,5,1,0,3,4,2.914828464578812,2.261191790401713,0.
+10,0,3,1,10,2,4,3,1,2.50303873714987,1.5051451493521428,0.
+1,1,10,0,1,3,0,3,1,1.3072742851724795,1.9423953478998177,0.
+1,1,3,2,3,2,0,3,2,2.3584710757323544,1.1305990899916512,0.
+10,3,2,3,3,1,6,1,0,1.4654067038271226,1.985658804098854,0.
+7,2,0,2,9,2,0,1,1,1.3873373238347515,0.7703862068685314,0.
+8,3,7,0,8,2,5,1,3,0.8253295318767577,2.899984116328765,0.
+1,1,0,0,0,1,10,2,1,2.3641832402927276,2.1083504299969444,0.
+8,2,7,3,10,1,0,2,2,1.3652741505836983,2.0949097737438462,0.
+8,3,10,0,4,2,5,1,3,0.24020243315960244,1.9022194572611664,0.
+9,1,9,2,0,0,4,1,1,2.5376357730342507,1.8526159863451914,0.
+7,3,9,3,3,3,7,0,3,1.6807661985105629,2.687649747358564,0.
+7,2,2,0,9,2,8,1,2,2.4565884018951802,0.8982169296897196,0.
+9,2,0,3,10,3,5,1,2,0.6353524665165655,1.0371641317201128,0.
+3,0,7,3,2,2,1,1,3,2.7958304824728186,2.4767642915078376,0.
+9,0,7,1,10,1,4,0,1,0.25616825512317254,1.308344122398832,0.
+1,2,1,3,10,1,1,1,2,2.224763269329724,2.231499665082513,0.
+8,0,4,3,0,2,5,3,3,0.770816371156068,2.3322630361929546,0.
+5,1,1,2,2,2,3,1,2,0.8022984805467428,1.4426168425194534,0.
+7,0,5,0,0,2,10,2,0,1.9511722131189702,0.4950560622323814,0.
+9,3,7,1,7,0,10,2,2,2.775167770563364,0.45822632372936667,0.
+7,2,7,3,9,0,6,1,1,2.8779651825740737,2.4436285639972706,0.
+1,1,6,3,0,0,3,3,3,0.5442622594206528,2.9853682012412737,0.
+10,2,8,1,4,2,7,3,2,1.329367888136959,0.9413429526552286,0.
+7,3,0,1,4,1,9,1,2,1.3457614299418887,1.5531678736115984,0.
+4,2,8,2,10,1,3,3,3,0.8248650313839705,0.6616884801498313,0.
+0,0,8,3,3,1,10,2,3,2.0673228373578825,1.4463899838479248,0.
+0,3,0,2,9,3,0,3,5,0.906810420580336,2.811065359991593,0.
+2,1,2,0,4,3,8,2,1,1.4719829696395097,1.5045170591445123,0.
+8,3,4,2,7,0,6,2,2,2.934501293001735,2.475691074854602,0.
+4,1,7,3,9,1,2,2,3,1.5229908742702083,2.7910336359576577,0.
+3,0,4,3,9,2,8,2,3,0.7395750631044415,2.7437950070306316,0.
+0,3,10,0,8,3,5,0,3,2.6218495656463876,0.917324005589486,0.
+0,1,10,3,3,3,7,3,4,1.2601693603785153,2.8311848091672465,0.
+5,3,3,1,0,2,0,1,2,1.5322630072331016,2.4612811983034675,0.
+2,0,10,3,4,3,8,2,3,2.581504744416619,2.3352397148922988,0.
+8,2,10,1,8,2,10,1,1,1.1752269484482167,1.7646051998689414,0.
+0,2,4,1,10,2,7,0,2,1.1544888515337508,2.637958865288499,0.
+4,0,6,3,6,1,7,3,3,1.4304562327889254,1.498204529248202,0.
+8,2,6,1,3,3,6,0,3,2.3823736912605336,2.0477368833060705,0.
+4,2,2,3,10,1,3,1,1,2.5205545321391254,0.7220948301952252,0.
+10,0,7,1,8,3,9,2,1,1.2260538418022682,2.451257871532359,0.
+0,0,0,3,10,1,10,3,3,2.006553045177734,1.173388489802477,0.
+6,3,3,3,1,3,1,2,5,2.9612859268853233,1.4353422269551883,0.
+2,1,2,0,2,3,1,2,1,1.3130842056605396,1.1311416250752515,0.
+3,3,9,3,5,1,3,3,3,0.4024708536057622,2.0642879809097643,0.
+8,2,7,2,3,2,7,2,0,1.2953268563208105,0.5693758931801445,0.
+0,3,9,1,8,1,8,1,2,2.152069679767771,1.41039953394532,0.
+10,3,10,1,2,1,2,3,3,1.922953978768665,1.3543149100065444,0.
+4,2,5,1,10,3,8,0,3,0.6098612038176467,0.5440300294800884,0.
+8,2,7,1,3,1,8,2,2,0.608144934570964,2.548777086726803,0.
+9,2,3,0,6,0,10,2,2,1.7250046113327357,1.415947697499866,0.
+4,1,3,2,7,0,3,3,3,1.5376308517695847,2.8988547754211718,0.
+5,3,4,0,6,3,4,2,3,0.7636088061341804,0.7149031039865208,0.
+4,3,9,3,2,0,6,2,2,1.8376914709492826,1.479579605670744,0.
+2,3,10,2,5,2,10,0,2,0.8658425814507988,1.1218947490443987,0.
+6,0,9,1,2,3,10,2,1,2.7269511334702985,1.2873996903274714,0.
+5,3,1,1,4,1,2,2,3,2.0472755894627443,1.4404427131274424,0.
+7,3,5,1,5,3,1,0,3,2.1329904508839563,1.7676149035568383,0.
+0,1,6,1,3,2,1,3,1,1.1352856084425453,2.0153824999161527,0.
+9,1,4,2,9,2,1,0,2,0.5735267931557,0.883812410190215,0.
+10,2,8,0,2,2,0,0,2,2.286451077255511,2.3140961306516,0.
+3,3,3,3,2,2,0,2,2,0.7388292443690174,1.6936545691041545,0.
+10,2,2,1,10,0,10,2,2,1.4443453550869583,1.112093007841184,0.
+7,3,6,1,10,3,2,2,4,1.8267825436432625,1.0428273092046476,0.
+0,0,1,3,9,2,9,1,3,2.4360668218156167,1.0158379194857559,0.
+9,0,6,3,8,1,5,2,3,2.4459546670789223,2.5273734233477834,0.
+8,2,0,2,9,3,2,2,1,0.9849892765217048,0.31610510276469483,0.
+2,0,6,0,3,1,9,1,0,1.646375524767966,0.39139222602227974,0.
+2,3,0,2,6,0,10,1,1,0.3040006591059585,1.784871103068122,0.
+7,2,2,2,4,3,6,0,3,1.5667275423919191,1.7334318105948183,0.
+5,3,2,1,1,2,4,0,2,0.4724440845832545,2.0191544118568805,0.
+2,1,3,2,6,1,8,3,2,1.6569706450379993,1.22551530354268,0.
+4,0,8,2,10,3,6,2,2,1.5517249658736345,1.763047858243108,0.
+6,2,4,2,0,3,3,0,3,2.056076795383781,1.693356584470254,0.
+1,3,6,1,4,1,3,3,3,1.9168275535161614,2.2496624307522994,0.
+6,1,1,1,8,1,3,0,1,0.7505340273835843,1.1586206830861343,0.
+1,2,9,2,8,2,1,0,2,2.5465024812722747,0.40185310698150234,0.
+6,2,1,3,8,3,5,3,2,1.0806918086518902,2.785171921529357,0.
+2,3,9,1,0,1,7,3,2,2.649104146816554,0.22330923947941184,0.
+3,3,2,2,8,2,7,1,1,1.4153027357986057,0.7592362643035888,0.
+1,2,4,1,2,0,8,1,1,0.3716613874614527,2.741327801112848,0.
+8,2,4,0,7,3,8,2,2,2.772853450624477,0.7349370680695975,0.
+4,2,10,1,9,1,10,2,2,2.3205630100773416,0.7942049271730163,0.
+1,3,6,3,6,1,4,3,2,1.3281998811443967,1.0363866114967322,0.
+3,2,7,3,3,1,3,3,4,0.8164145918403629,1.064457528831932,0.
+5,1,6,0,6,1,5,0,1,2.7770950000390737,1.3710296853982933,0.
+4,3,3,3,3,1,2,2,3,2.714056705000763,2.4602516368227603,0.
+4,1,5,2,3,2,7,3,1,0.35146521178282075,0.3571498607132635,0.
+9,1,8,2,9,1,6,1,1,1.2274240317357688,1.6715623227186507,0.
+1,1,3,3,8,3,7,0,3,1.9317099576956442,1.552216114675943,0.
+9,0,3,1,8,0,2,1,1,2.9043448611385188,2.7338758740571345,0.
+6,0,8,2,2,2,4,1,2,1.675816848063465,0.8679109639840443,0.
+9,1,9,1,0,1,1,0,1,0.596382672070606,1.7119233696062492,0.
+7,2,3,1,5,1,0,3,3,1.374644910756329,2.8427817664580024,0.
+6,3,10,3,7,1,10,1,2,2.6167816000584327,0.8461412388948104,0.
+5,2,7,1,6,3,2,3,2,1.372612635683653,2.0961606258142194,0.
+8,1,7,2,8,3,1,2,3,0.9231592095990986,0.4730174634532238,0.
+3,1,1,1,10,1,10,3,2,0.3461817829038991,1.5964106156876117,0.
+8,0,7,3,0,2,4,1,3,2.096022436900264,1.9588093081620253,0.
+4,3,9,3,5,1,9,1,0,2.449961912097761,0.4920648568697943,0.
+4,0,8,1,6,2,1,3,1,2.7092656850999868,1.2113036872224647,0.
+2,2,4,3,2,2,1,0,2,0.7822138935161549,0.8243720955516407,0.
+5,3,9,3,8,0,7,3,3,1.245926765783301,2.995825714452799,0.
+9,1,9,0,5,3,1,3,1,0.9447956603373777,1.9688135363735153,0.
+5,2,2,2,3,3,7,1,2,2.2346399253675413,0.7417391778003699,0.
+0,1,8,3,0,2,2,3,2,1.0232322455743281,1.8925470263446593,0.
+5,3,2,3,8,2,9,0,2,2.9817835096507697,2.111249932049714,0.
+9,2,8,2,1,1,4,3,4,0.4083625562135289,0.40305961637541055,0.
+4,3,9,2,8,2,10,1,3,2.7964161333721353,2.78601941028056,0.
+1,1,10,0,10,0,10,1,1,2.8918303445485822,2.7282191833543576,0.
+3,1,1,3,0,1,9,2,3,0.5883986968362152,1.1131593470756491,0.
+3,1,5,3,6,3,10,0,3,1.29122181890919,2.1687719665184835,0.
+5,1,1,0,10,1,9,1,1,2.7757481174210463,2.3835512920725517,0.
+3,3,1,3,10,0,6,1,1,1.127749467977892,1.3592596107053452,0.
+0,1,10,2,10,2,4,0,2,1.403190963591391,1.7426120757596841,0.
+8,2,0,3,2,3,4,2,3,2.21376917027845,1.2373565371204496,0.
+6,3,2,2,0,2,6,2,4,2.0909192232076848,2.9595518129553477,0.
+2,0,4,3,0,3,2,0,3,2.9389053651691066,0.2571903449038695,0.
+3,0,8,2,2,2,8,3,2,1.1292930350844212,2.662958886705143,0.
+4,2,9,3,3,1,8,2,2,1.2069732981501304,0.6603448818737543,0.
+8,1,1,2,9,2,2,3,2,1.176331991745739,2.9348246139214273,0.
+8,2,9,1,3,2,6,0,2,1.2579465199621849,2.8708297284639936,0.
+2,1,0,0,4,2,0,2,1,0.5068840995027317,1.4926314796061328,0.
+8,0,5,2,0,1,6,1,2,2.083459207268281,2.044350098956075,0.
+8,1,3,1,8,3,3,1,2,0.27182349339562206,0.25278077977555524,0.
+9,0,9,3,8,2,7,1,3,2.0842851912303697,1.4617238592678117,0.
+6,0,7,3,4,2,5,3,3,0.5438422711246416,1.1960642527951837,0.
+3,0,0,3,1,3,5,2,3,2.142935683266975,0.3991912324294198,0.
+5,2,1,1,6,3,6,2,2,1.7422616583207757,1.1999207722763399,0.
+1,2,1,0,8,2,9,0,2,0.32666758600821844,2.6968272363747605,0.
+3,2,0,0,10,1,6,1,2,1.6221941314416704,1.1434410020543084,0.
+0,3,1,2,10,1,3,1,2,0.9726067645442806,0.5905895434956969,0.
+0,1,0,0,1,0,7,1,1,2.900837043628563,2.949017267115262,0.
+8,3,8,1,2,0,0,2,2,1.5234948494086824,0.3368477864264441,0.
+7,2,4,2,8,1,6,1,1,0.650243384315508,0.5631520583276046,0.
+3,3,9,3,2,0,9,0,0,2.3807738065036146,0.579866727058429,0.
+4,3,7,0,10,2,6,3,3,1.3173017138545768,0.5962764653299519,0.
+8,3,0,3,0,3,2,1,4,0.6053424778386112,0.8254011764061571,0.
+9,1,3,1,1,0,2,2,2,0.9352316673114225,2.648521223427407,0.
+4,3,6,2,1,3,8,3,3,1.1811263754353778,0.27949826874059935,0.
+6,2,9,1,0,3,8,3,2,1.7891686343587905,2.5475273097794897,0.
+3,3,6,3,1,1,4,2,1,0.8652651806579512,2.144133642324886,0.
+8,2,5,1,10,3,3,1,3,1.3153722390010687,2.370685736109543,0.
+9,1,5,2,7,3,2,1,2,1.4002015422875607,1.5704174280943155,0.
+4,3,8,3,9,1,4,2,1,1.1067444816569871,2.7633308770019767,0.
+8,1,5,1,7,0,7,1,1,1.3033564335609653,1.8205597872081913,0.
+6,2,10,2,0,2,9,0,2,1.417347861506408,0.8228017043867935,0.
+3,3,5,2,8,1,5,1,1,2.3178061623921504,2.446585257893128,0.
+0,2,5,2,5,1,0,0,1,0.417745684372842,0.962045539464802,0.
+8,1,8,1,8,1,5,0,1,2.241624712766236,2.117102285343411,0.
+5,2,1,1,2,0,2,3,3,1.452790633495452,2.3926515680383504,0.
+7,1,0,3,9,1,9,1,2,1.320602559814661,0.5022153836743799,0.
+10,2,8,3,1,1,4,3,4,0.8992253656072662,2.3674150876262807,0.
+2,3,9,0,10,1,9,2,3,1.0158467671907352,1.376506498457696,0.
+1,3,1,1,4,1,4,1,2,1.8583872135299577,2.4330946543470144,0.
+8,0,9,1,3,0,7,1,1,2.9549925009110574,0.7408650898965572,0.
+5,1,3,1,2,2,0,0,2,0.7685503537061233,1.8557635823419316,0.
+10,1,7,0,3,1,7,1,1,2.91588044339084,0.4955239475168276,0.
+6,2,10,0,10,3,1,1,2,1.2165093404056426,0.48447018353308113,0.
+8,3,0,2,9,0,9,3,3,1.2458573416728456,1.8984850604019927,0.
+5,1,6,1,5,2,2,3,2,1.6678698181298168,0.27271548312207194,0.
+1,0,5,1,6,3,7,3,1,2.4353366447686176,2.251622116348118,0.
+4,1,6,2,10,3,8,1,3,0.9785359746307756,2.2569876645542513,0.
+7,1,8,0,3,0,8,1,1,2.6254607678979696,1.4870987629995316,0.
+7,2,4,0,4,1,6,2,2,2.4860925929709445,2.078702116156186,0.
+4,3,10,0,2,3,4,1,3,2.7833993206816183,2.4091684481860574,0.
+10,2,3,2,10,0,3,2,2,2.600869455586518,2.7023313585291886,0.
+10,0,3,2,1,1,5,1,2,0.6308140012941346,2.192072324433969,0.
+0,2,1,3,8,3,9,2,5,0.47211349059914465,1.2835269519982067,0.
+1,3,5,3,9,0,8,1,1,0.39051019613873894,0.6191713045924181,0.
+8,2,1,2,5,2,10,1,1,2.8083327076207043,2.7889141797722985,0.
+0,1,6,2,4,3,1,3,3,2.576880910721891,2.8940089415951142,0.
+10,2,9,2,3,1,7,1,2,0.4344721828159619,2.6268647157878835,0.
+4,1,0,3,8,0,10,3,3,2.3791812481389023,0.41047900744269716,0.
+9,0,0,3,9,2,2,1,3,0.3591915419041567,2.665840985111095,0.
+9,3,6,2,5,1,5,2,3,2.1661396893720006,1.1610985541712369,0.
+7,3,2,3,2,1,5,2,2,1.2178868457103484,1.906937613481415,0.
+9,1,4,2,1,2,5,2,1,1.76351959987508,0.310172698567285,0.
+10,1,4,2,5,1,5,1,2,2.9057790676388597,2.246935216027948,0.
+9,3,6,1,7,0,5,2,2,1.8228283848861624,0.5136268936144317,0.
+4,1,2,1,10,2,1,2,2,2.4146521536083805,0.4173254666656194,0.
+9,1,2,1,3,2,8,0,2,2.091075002425989,1.6190364383142475,0.
+5,2,7,1,2,2,7,1,2,1.9881551763498084,0.9065836043584428,0.
+5,3,6,2,4,1,7,2,1,0.2437619648775846,0.4633502649952983,0.
+6,1,8,1,10,2,5,1,1,1.979325271197399,1.1385489265368367,0.
+3,2,6,2,4,0,8,3,3,2.761222760478871,2.281187219578926,0.
+7,3,2,2,3,0,4,1,1,1.387675975897055,2.978542623739007,0.
+0,2,8,0,0,1,2,2,2,0.7224763006498796,1.564980597384034,0.
+3,3,7,1,10,2,10,2,3,2.9525286612356716,0.8408929108171121,0.
+0,1,7,2,2,1,9,0,1,0.9126033819378172,2.8275866134399896,0.
+9,2,1,0,8,0,0,2,2,1.0314013900322792,0.6732978424388585,0.
+8,2,9,1,7,1,9,0,1,1.3865069980063334,2.540865404640523,0.
+6,2,0,2,2,1,9,2,3,0.4235665273207312,1.6846937203861359,0.
+5,2,5,2,5,0,3,1,1,2.7659439645673256,1.0115836172242267,0.
+0,1,9,3,6,3,2,1,3,1.2065541331332832,1.0405234817364937,0.
+8,1,6,1,4,0,1,1,1,0.7227473981225692,1.155456043230863,0.
+2,3,8,0,1,0,3,3,3,1.9299731710816523,0.6103164598521023,0.
+4,0,2,3,0,2,1,2,3,2.8621863382346,2.9310316681902524,0.
+8,3,5,1,8,2,6,3,4,0.45896263960829176,1.9911804294619255,0.
+1,0,3,2,6,2,4,0,2,2.335589089980126,1.0652452461608548,0.
+3,1,1,3,1,2,10,3,3,1.0154942249838212,1.8603629108124498,0.
+2,1,6,0,3,2,2,2,1,0.6444259802836747,1.4787144338375233,0.
+7,1,10,2,2,1,3,2,3,0.7147755411709733,0.741536139217259,0.
+0,2,1,3,2,1,8,1,1,1.8857744159033487,2.4763969783476716,0.
+2,3,8,0,10,3,4,0,3,1.9960210966818677,2.9863811890951855,0.
+2,3,2,3,2,1,5,2,2,1.4578550110452344,2.0655616115569337,0.
+8,3,7,2,0,0,0,2,2,2.1690558898141488,1.6547046613975827,0.
+7,2,3,1,6,3,3,0,3,2.553795269242924,2.98657012675977,-5.8584894789029045
+5,2,2,2,8,0,7,0,0,2.111885434964618,1.9937513659932842,0.
+8,3,8,0,8,2,0,3,3,2.8705609121506805,1.3331637714091094,0.
+3,2,6,3,1,2,0,1,1,2.34251648574082,0.5435767229545014,0.
+4,0,6,2,6,0,9,2,2,1.20338522636726,2.788327210207373,0.
+7,1,1,0,10,1,0,0,1,2.3923670262306302,2.4029578275837586,0.
+0,1,8,3,6,2,0,3,2,1.203543583831212,2.036143250730075,0.
+5,3,3,1,10,3,9,1,2,2.7314904241739955,0.21670230017016534,0.
+7,3,0,3,4,0,9,3,3,0.7221556673410454,2.504574860198268,0.
+0,1,9,0,7,2,2,1,1,0.6129040077484231,0.9377073179769542,0.
+8,3,8,3,8,0,6,1,1,0.2916837472356941,1.6706587000955677,0.
+5,0,8,1,2,1,9,1,1,1.1071388188914164,2.907305251123775,0.
+4,2,3,3,2,1,9,0,1,2.6049414020919137,1.1423144552963707,0.
+0,2,3,2,6,0,3,1,1,0.6589475462254843,0.3458167300233934,0.
+5,1,7,2,1,1,9,1,2,0.37152940436379245,2.257037185691721,0.
+5,2,5,2,9,3,3,2,4,1.1411180980151068,1.360756047176233,0.
+7,2,5,3,4,3,6,3,5,2.885010245830152,0.6362691751312788,0.
+5,2,9,3,10,1,0,1,2,0.37940369466801593,1.2334735865845734,0.
+0,0,7,2,1,2,7,1,2,1.5479759698453543,1.9801733480094859,0.
+2,2,4,0,10,0,10,2,2,1.3293709201965624,1.2689927819499371,0.
+9,0,3,1,0,2,8,2,1,2.164068098414817,1.2453603926890962,0.
+0,2,9,2,3,1,4,0,1,0.6696321137903438,2.6084879932864853,0.
+8,1,3,3,4,1,6,1,2,2.042466846990008,1.8892178256161216,0.
+2,3,2,1,10,3,10,2,3,1.4323483837153708,1.0105488574564596,0.
+5,0,10,3,7,3,10,2,3,2.726039671844071,2.8549137956975468,0.
+2,3,1,2,3,3,10,1,3,2.351151357868736,0.36675214929021926,0.
+2,0,7,0,9,1,3,1,0,0.7394282604433373,0.2998036281738994,0.
+4,3,0,3,6,2,9,2,4,1.4189068677052719,2.9473449090470467,0.
+9,2,1,3,9,2,7,0,2,2.2125653578966116,0.6449925605947713,0.
+3,3,9,3,7,0,2,0,0,2.6527765476432785,0.7017535912522956,0.
+3,2,1,3,10,3,10,2,3,1.3545764993952725,1.4463176767351857,0.
+6,0,6,3,7,3,10,1,3,0.5231215491635299,1.5951002259931988,0.
+10,1,0,3,0,0,3,2,2,0.932146690038874,0.6396108542280041,0.
+1,3,4,3,3,1,5,1,2,2.6568428894223652,1.0041895351349086,0.
+2,0,6,3,7,2,5,2,3,2.1388490384614993,1.084038962927921,0.
+1,0,4,3,3,3,6,1,3,1.1954487430499983,1.3480169360596266,0.
+4,2,10,3,7,0,8,3,3,0.9989698852671967,0.31492500891608577,0.
+0,1,5,1,4,1,7,1,2,0.2915533255281222,0.8488169939496006,0.
+5,3,8,2,6,3,10,0,3,1.3273506939181288,1.9576844972216119,0.
+4,1,1,1,8,1,3,1,0,2.5512768516580167,1.1776908814863933,0.
+0,1,7,2,4,3,5,2,3,2.55835426097922,0.7999378464474947,0.
+7,1,3,2,4,3,4,3,3,2.9469964384408645,1.0130439824706348,0.
+10,1,10,1,4,0,9,2,2,2.29282320065841,1.5895923699679324,0.
+3,3,8,2,9,2,8,3,5,1.4288368272508416,0.32914872565251896,0.
+5,2,10,0,7,1,4,2,2,2.0963238426388164,2.5184909838350205,0.
+7,3,4,0,1,2,10,3,3,0.7462285283562706,1.9678766341114837,0.
+5,1,6,0,6,1,1,0,1,2.742103144450633,1.0100914398280834,0.
+6,2,8,2,1,0,8,0,0,2.9686558639239102,2.571359181752765,0.
+6,2,1,0,8,3,0,2,2,2.8725089493368277,2.890444927773988,0.
+3,2,6,3,9,1,2,2,1,0.9140224473432359,1.7614216539656447,0.
+9,1,3,1,2,3,8,2,2,1.5224977087087135,1.7130595192975742,0.
+2,1,5,1,1,1,2,1,0,0.620215788594543,2.125090314126055,0.
+3,3,2,2,9,3,10,3,5,2.1571418399260214,0.2352162581649102,0.
+4,1,3,3,8,2,1,0,2,2.5710875717760064,1.8947511401060053,0.
+8,2,5,2,6,1,8,1,0,1.5775238404669114,2.9176906027237717,0.
+10,2,1,2,2,3,3,3,2,2.9869659703459304,0.6494176627258623,0.
+9,1,7,1,7,0,5,0,0,0.35843660067087235,1.1434109912448518,0.
+8,0,0,3,0,1,7,3,3,2.5846819654664843,0.30192303129821685,0.
+9,0,10,1,3,0,4,1,1,2.2400726780346583,2.76156913953512,0.
+0,2,6,1,3,1,3,0,1,1.2200817282296477,0.8137895784475973,0.
+4,3,0,1,10,1,3,2,3,1.1171419500565256,2.8469536184972846,0.
+6,1,3,2,2,3,3,0,3,0.8843975271072293,2.0058645383535723,0.
+0,2,4,3,3,1,1,2,1,2.612997115865329,0.34104461891724425,0.
+3,0,1,3,6,3,10,2,3,2.716076032829566,1.791883854821867,0.
+0,3,1,3,1,2,7,2,2,2.989190167535284,0.6804444348264811,0.
+7,0,7,0,6,1,8,1,0,1.4624896265248113,1.7583815597493881,0.
+6,3,6,3,1,1,9,2,3,2.4493465867634194,1.3172740005089132,0.
+3,2,1,2,5,3,8,1,2,2.9637052407119135,1.386678142676876,0.

test/racahs_reduction_formula.jl

@@ -0,0 +1,43 @@
+using QuadGK
+include("../math.jl")
+
+# Testing ∫dp p u' u where u' and u' may have different l
+
+function explicit_integral(np, lp, n, l, μω)
+    integrand(p) = p * (-1)^(n + np) * 1/sqrt(μω) * ho_basis(l, n, 1/sqrt(μω) * p) * ho_basis(lp, np, 1/sqrt(μω) * p)
+    (integral, _) = quadgk(integrand, 0, Inf; atol=1e-10)
+    return integral
+end
+
+n_list = 0 : 10
+l_list = 0 : 5
+μω_list = [1, 2.66]
+
+for np in n_list
+    for n in n_list
+        for lp in l_list
+            for l in l_list
+                for μω in μω_list
+                    println("Checking n=$n, np=$np, l=$l, lp=$lp, μω=$μω")
+                    expl = explicit_integral(np, lp, n, l, μω)
+                    anal = integral(np, lp, n, l, μω)
+                    @assert isapprox(expl, anal; atol=1e-8) "Explicit=$(expl) and analytical=$(anal)"
+                end
+            end
+        end
+    end
+end
+
+# Testing the full Racah's reduction formula
+
+test_data = DelimitedFiles.readdlm("test/racah_test_data.csv", ',') # format = n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, μ1ω1, μ2ω2, racah_matrix_element
+test_data_ints = Int.(test_data[:, 1:(end - 3)])
+test_data_floats = test_data[:, (end - 2):end]
+
+for i in 1 : size(test_data)[1]
+    n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ = test_data_ints[i, :]
+    μ1ω1, μ2ω2, racah_matrix_element = test_data_floats[i, :]
+    println("Checking n1'=$n1p, l1'=$l1p, n2'=$n2p, l2'=$l2p, n1=$n1, l1=$l1, n2=$n2, l2=$l2, Λ=$Λ, μ1ω1=$μ1ω1, μ2ω2=$μ2ω2")
+    val = racahs_reduction_formula(n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, μ1ω1, μ2ω2)
+    @assert isapprox(val, racah_matrix_element; atol=1e-8) "Formula=$(val) and test reference=$racah_matrix_element"
+end

test/simple_berggren.jl

@@ -0,0 +1,30 @@
+using Plots
+include("../common.jl")
+include("../p_space.jl")
+
+vertices = [0, 0.5 - 0.3im, 1, 6]
+subdivisions = [64, 64, 64]
+p, w = get_mesh(vertices, subdivisions)
+
+# resonance example from my thesis
+V_basis(p, q) = 2*g0(4, p, q) - 3*g0(2, p, q)
+
+basis_eig = eigen(get_H_matrix(V_basis, p, w))
+basis = basis_eig.vectors .* w
+
+basis_evals = basis_eig.values
+E_target = 0.7
+E = nearest(basis_evals, E_target)
+
+print("pole E = $E")
+
+# ResonanceEC: Eq. (20)
+V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
+
+H = get_H_matrix(V_system(0.45), p, w)
+H_berggren = transpose(basis) * H * basis
+
+evals = eigvals(H)
+scatter(real.(evals), imag.(evals), label="E")
+scatter!(real.(basis_evals), imag.(basis_evals), label="basis", markershape=:x)
+xlims!((0, 1))