Improve eigenstates

examples/generate_Yb174_database.jl

@@ -30,25 +30,20 @@ low_l_models = [
     MQDT.Yb174.FMODEL_HIGHN_G5,
 ]
 
-# bounds (as found in the parameters file)
-low_n_min = [1, 2, 1.5, 1.7, 2.7, 1.5, 2, 2, 2] .+ 0.5
-low_n_max = [2, 26, 5.5, 2.7, 5.7, 4.5, 26, 5, 18] .- 0.5
+n_max = 30
 
 # calculate low \nu MQDT states
-low_n_states = [eigenstates(low_n_min[i], low_n_max[i], low_n_models[i], parameters) for i in eachindex(low_n_min)]
-
-# bounds (as found in the parameters file)
-n_min = [2, 26, 6, 6, 5, 26, 5, 18, 25, 7, 25, 25, 25, 25] .+ 0.5
-n_max = 30
+low_n_states = [eigenstates(NaN, n_max, model, parameters) for model in low_n_models]
 
 # calculate high \nu, low \ell MQDT states
-low_l_states = [eigenstates(n_min[i], n_max, low_l_models[i], parameters) for i in eachindex(n_min)]
+low_l_states = [eigenstates(NaN, n_max, model, parameters) for model in low_l_models]
 
 # calculate high \ell SQDT states
+n_min_high_l = 25 # minimum n for high-l states
 l_max = n_max - 1
 MQDT.wigner_init_float(n_max, "Jmax", 9) # initialize Wigner symbol calculation
 high_l_models = single_channel_models(:Yb174, 5:l_max)
-high_l_states = [eigenstates(25, n_max, M, parameters) for M in high_l_models]
+high_l_states = [eigenstates(n_min_high_l, n_max, M, parameters) for M in high_l_models]
 
 # generate basis and calculate matrix elements
 basis = basisarray(vcat(low_n_states, low_l_states, high_l_states), vcat(low_n_models, low_l_models, high_l_models))

src/boundstates.jl

@@ -294,30 +294,38 @@ end
     mroots(N::Number, M::Model, P::Parameters)
     mroots(N1::Number, N2::Number, M::Model, P::Parameters)
 
-Function that returns the roots of the M matrix in the range from `N` to `N+1` or, if provided, from `N1` to `N2`.
+Function that returns the roots of the M matrix in the range from `N-0.5` to `N+0.5` or, if provided, from `N1` to `N2`.
 Root finding algorithm used from the package `Roots`.
 """
 function mroots(N::Number, M::Model, P::Parameters)
-    m(n) = det(mmat(n, M, P))
-    z = find_zeros(m, (N-0.5, N+0.5))
-    c = Int[]
-    for i in eachindex(z)
-        mz = m(z[i])
-        if abs(mz) > 1e-10
-            println("Warning: skipped a root that seems inaccurate. Value was $(mz) at n=$(z[i]) for $(M.name).")
-        else
-            push!(c, i)
-        end
-    end
-    return z[c]
+    return mroots(N-0.5, N+0.5, M, P)
 end
 
 function mroots(N1::Number, N2::Number, M::Model, P::Parameters)
-    z = Float64[]
-    for i in N1:N2
-        append!(z, mroots(i, M, P))
+    nus = Float64[]
+    m(n) = det(mmat(n, M, P))
+    N_list = [N1; (N1 + 0.5):N2; N2]
+    for (i, _N1) in enumerate(N_list[1:(end - 1)])
+        _N2 = N_list[i + 1]
+        if abs(_N1 - _N2) < 1e-10
+            _N2 = _N2 + 1e-10
+        end
+        z = find_zeros(m, (_N1, _N2))
+        new_nus = Float64[]
+        for nu in z
+            if abs(m(nu)) > 1e-10
+                println("Warning: skipped a root that seems inaccurate. Value was $(m(nu)) at n=$(nu) for $(M.name).")
+            else
+                push!(new_nus, nu)
+            end
+        end
+        if i != 1 && length(new_nus) > 0 && abs(new_nus[1] - _N1) < eps()
+            # if the first root is at the lower boundary, skip it unless it's the very first interval
+            popfirst!(new_nus)
+        end
+        append!(nus, new_nus)
     end
-    return z
+    return nus
 end
 
 """
@@ -363,8 +371,25 @@ See also [`eigenenergies`](@ref)
 
 Function that returns the bound states corresponding to an MQDT model in the form of an instance of `EigenStates`,
 which contains the global reference principal quantum number, the channel-dependent principal quantum numbers, as well as the channel coefficients.
-"""
-function eigenstates(N1::Number, N2::Number, M::Model, P::Parameters)
+N1 can be NaN to use the lowest possible effective principal quantum number for the model M.
+N2 can only be NaN if the model M has a defined maximum effective principal quantum number.
+We will always use the higher of N1 and the M's minimum effective principal quantum number
+and the lower of N2 and the M's maximum effective principal quantum number.
+
+If overwrite_model_limits is True, we will ignore the model nu limits and always use the given bounds.
+"""
+function eigenstates(N1::Number, N2::Number, M::Model, P::Parameters; overwrite_model_limits::Bool=false)
+    if !overwrite_model_limits
+        nu_min_model, nu_max_model = get_nu_limits_from_model(M)
+        N1 = isnan(N1) ? nu_min_model : max(N1, nu_min_model)
+        N2 = isnan(N2) ? nu_max_model : min(N2, nu_max_model)
+    elseif isnan(N1) || isnan(N2)
+        error("When overwrite_model_limits is true, N1 and N2 cannot be NaN.")
+    end
+    if isinf(N2)
+        error("Cannot use infinite N2 if the model does not define a maximum effective principal quantum number.")
+    end
+
     z, n = eigenenergies(N1, N2, M, P)
     mfunc(e) = mmat(e, M, P)
     a = similar(n)

src/general.jl

@@ -570,7 +570,7 @@ end
 Get the minimum and maximum effective principal quantum number (nu_min, nu_max) from a given Model.
 
 The name of a Model usually is something like "S J=0, ν > 2" or "S F=1/2, ν > 26",
-from this we can extract nu_min = 2 and nu_max NaN.
+from this we can extract nu_min = 2 and nu_max Inf.
 
 Some Model names also contain a upper limit, e.g. "P J=1, 1.7 < ν < 2.7",
 from which we extract nu_min = 1.7 and nu_max = 2.7.
@@ -586,7 +586,7 @@ function get_nu_limits_from_model(model::Model)
     m = match(r"ν\s*>\s*([0-9/\.]+)", model.name)
     if m !== nothing
         nu_min = parse(Float64, m.captures[1])
-        return nu_min, NaN
+        return nu_min, Inf
     end
 
     throw(ArgumentError("No match found for 'ν > ...' or '... < ν < ...' in model: $(model.name)"))
@@ -678,7 +678,7 @@ function single_channel_models(species::Symbol, l::Integer)
     for i in eachindex(jt)
         m[i] = fModel(
             species,
-            "L=$l, J=$(jt[i]), Jr=$(jr[i])",
+            "L=$l, J=$(jt[i]), Jr=$(jr[i]), ν > $(l+1)",
             1,
             [""],
             Bool[1],