diff --git a/examples/generate_Yb174_database.jl b/examples/generate_Yb174_database.jl
index f4522bd..ce66bf4 100644
--- a/examples/generate_Yb174_database.jl
+++ b/examples/generate_Yb174_database.jl
@@ -1,52 +1,47 @@
 using MQDT
 
-# load Yb174 data
-parameters = MQDT.Yb174.PARA
-low_n_models = [
-    MQDT.Yb174.FMODEL_LOWN_S0,
-    MQDT.Yb174.FMODEL_LOWN_S1,
-    MQDT.Yb174.FMODEL_LOWN_P0,
-    MQDT.Yb174.FMODEL_LOWEST_P1,
-    MQDT.Yb174.FMODEL_LOWN_P1,
-    MQDT.Yb174.FMODEL_LOWN_P2,
-    MQDT.Yb174.FMODEL_LOWN_D1,
-    MQDT.Yb174.FMODEL_LOWN_D2,
-    MQDT.Yb174.FMODEL_LOWN_D3,
-]
-low_l_models = [
-    MQDT.Yb174.FMODEL_HIGHN_S0,
-    MQDT.Yb174.FMODEL_HIGHN_S1,
-    MQDT.Yb174.FMODEL_HIGHN_P0,
-    MQDT.Yb174.FMODEL_HIGHN_P1,
-    MQDT.Yb174.FMODEL_HIGHN_P2,
-    MQDT.Yb174.FMODEL_HIGHN_D1,
-    MQDT.Yb174.FMODEL_HIGHN_D2,
-    MQDT.Yb174.FMODEL_HIGHN_D3,
-    MQDT.Yb174.FMODEL_HIGHN_F2,
-    MQDT.Yb174.FMODEL_HIGHN_F3,
-    MQDT.Yb174.FMODEL_HIGHN_F4,
-    MQDT.Yb174.FMODEL_HIGHN_G3,
-    MQDT.Yb174.FMODEL_HIGHN_G4,
-    MQDT.Yb174.FMODEL_HIGHN_G5,
-]
-
+species = :Yb174
+parameters = MQDT.get_species_parameters(species)
 n_max = 30
+n_min_high_l = 25 # minimum n for high-l states
 
-# calculate low \nu MQDT states
-low_n_states = [eigenstates(NaN, n_max, model, parameters) for model in low_n_models]
+all_models = Vector{MQDT.fModel}()
 
-# calculate high \nu, low \ell MQDT states
-low_l_states = [eigenstates(NaN, n_max, model, parameters) for model in low_l_models]
+sr = 1 / 2
+Jc = 1 / 2
+ic = parameters.spin
+for lr in 0:(n_max - 1)
+    for Jr in abs(lr - sr):1:(lr + sr)
+        for Fc in abs(Jc - ic):1:(Jc + ic)
+            for F in abs(Fc - Jr):1:(Fc + Jr)
+                models = MQDT.get_fmodels(species, lr, Jr, Fc, F, parameters)
+                for model in models
+                    if !any(m -> m.name == model.name, all_models)
+                        push!(all_models, model)
+                    end
+                end
+            end
+        end
+    end
+end
 
-# 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(n_min_high_l, n_max, M, parameters) for M in high_l_models]
+@info "Calculating MQDT states..."
+states = Vector{MQDT.EigenStates}(undef, length(all_models))
+for (i, M) in enumerate(all_models)
+    n_min = NaN
+    if startswith(M.name, "SQDT")
+        n_min = n_min_high_l
+    end
+    @info "$(M.name)"
+    states[i] = MQDT.eigenstates(n_min, n_max, M, parameters)
+    @info "  found nu_min=$(minimum(states[i].n)), nu_max=$(maximum(states[i].n)), total states=$(length(states[i].n))"
+end
+
+# generate basis
+basis = basisarray(states, all_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))
+# calculate matrix elements
+MQDT.wigner_init_float(n_max, "Jmax", 9) # initialize Wigner symbol calculation
 @time me = matrix_elements(basis, parameters)
 
 # prepare tables
diff --git a/examples/generate_lu_fano.jl b/examples/generate_lu_fano.jl
index 3f31d89..9b41b8a 100644
--- a/examples/generate_lu_fano.jl
+++ b/examples/generate_lu_fano.jl
@@ -35,20 +35,20 @@ end
 
 # S series of Yb171
 f = 0.5
-s_S05 = eigenstates(24, 129, MQDT.Yb171.FMODEL_HIGHN_S05, MQDT.Yb171.PARA)
+s_S05 = eigenstates(24, 129, MQDT.Yb171.FMODEL_HIGHN_S05, MQDT.Yb171.PARAMETERS)
 b_S05 = basisarray([s_S05], [MQDT.Yb171.FMODEL_HIGHN_S05])
-m_S05 = [MQDT.multipole_moments(b, b, MQDT.Yb171.PARA)[4] for b in b_S05.states]
+m_S05 = [MQDT.multipole_moments(b, b, MQDT.Yb171.PARAMETERS)[4] for b in b_S05.states]
 g_S05 = -2m_S05 / f * wigner3j(1, 0, f, f, f, f)
 
 f = 1.5
-s_S15 = eigenstates(24, 129, MQDT.Yb171.FMODEL_HIGHN_S15, MQDT.Yb171.PARA)
+s_S15 = eigenstates(24, 129, MQDT.Yb171.FMODEL_HIGHN_S15, MQDT.Yb171.PARAMETERS)
 b_S15 = basisarray([s_S15], [MQDT.Yb171.FMODEL_HIGHN_S15])
-m_S15 = [MQDT.multipole_moments(b, b, MQDT.Yb171.PARA)[4] for b in b_S15.states]
+m_S15 = [MQDT.multipole_moments(b, b, MQDT.Yb171.PARAMETERS)[4] for b in b_S15.states]
 g_S15 = -2m_S15 / f * wigner3j(1, 0, f, f, f, f)
 
 scatter(; layout=(2, 1), size=(400, 500), legend=:topleft)
-scatter!(s_S05.n, mod.(-s_S05.nu[7, :], 1); title="¹⁷¹Yb S series", subplot=1, label="F=1/2")
-scatter!(s_S15.n, mod.(-s_S15.nu[1, :], 1); xlabel="ν", ylabel="μ", subplot=1, label="F=3/2")
+scatter!(s_S05.n, mod.(-s_S05.n, 1); title="¹⁷¹Yb S series", subplot=1, label="F=1/2")
+scatter!(s_S15.n, mod.(-s_S15.n, 1); xlabel="ν", ylabel="μ", subplot=1, label="F=3/2")
 scatter!(s_S05.n, g_S05; subplot=2)
 scatter!(s_S15.n, g_S15; xlabel="ν", ylabel="g factor", subplot=2)
 hline!([g_F(1/2, 1/2, 0, 0, 0)]; l=:dash, label="¹S₁", subplot=2)
@@ -59,21 +59,21 @@ savefig("Yb171_S_series.pdf")
 
 # P series of Yb171
 f = 0.5
-s_P05 = eigenstates(10, 70, MQDT.Yb171.FMODEL_HIGHN_P05, MQDT.Yb171.PARA)
+s_P05 = eigenstates(10, 70, MQDT.Yb171.FMODEL_HIGHN_P05, MQDT.Yb171.PARAMETERS)
 b_P05 = basisarray([s_P05], [MQDT.Yb171.FMODEL_HIGHN_P05])
-m_P05 = [MQDT.multipole_moments(b, b, MQDT.Yb171.PARA)[4] for b in b_P05.states]
+m_P05 = [MQDT.multipole_moments(b, b, MQDT.Yb171.PARAMETERS)[4] for b in b_P05.states]
 g_P05 = -2m_P05 / f * wigner3j(1, 0, f, f, f, f)
 
 f = 1.5
-s_P15 = eigenstates(10, 70, MQDT.Yb171.FMODEL_HIGHN_P15, MQDT.Yb171.PARA)
+s_P15 = eigenstates(10, 70, MQDT.Yb171.FMODEL_HIGHN_P15, MQDT.Yb171.PARAMETERS)
 b_P15 = basisarray([s_P15], [MQDT.Yb171.FMODEL_HIGHN_P15])
-m_P15 = [MQDT.multipole_moments(b, b, MQDT.Yb171.PARA)[4] for b in b_P15.states]
+m_P15 = [MQDT.multipole_moments(b, b, MQDT.Yb171.PARAMETERS)[4] for b in b_P15.states]
 g_P15 = -2m_P15 / f * wigner3j(1, 0, f, f, f, f)
 
 scatter(; layout=(2, 2), size=(700, 500))
-scatter!(s_P05.n, mod.(-s_P05.nu[1, :], 1); ylabel="μ", title="¹⁷¹Yb P F=1/2", subplot=1)
+scatter!(s_P05.n, mod.(-s_P05.n, 1); ylabel="μ", title="¹⁷¹Yb P F=1/2", subplot=1)
 scatter!(s_P05.n, g_P05; xlabel="ν", ylabel="g factor", subplot=3)
-scatter!(s_P15.n, mod.(-s_P15.nu[1, :], 1); ylabel="μ", title="¹⁷¹Yb P F=3/2", subplot=2)
+scatter!(s_P15.n, mod.(-s_P15.n, 1); ylabel="μ", title="¹⁷¹Yb P F=3/2", subplot=2)
 scatter!(s_P15.n, g_P15; xlabel="ν", ylabel="g factor", subplot=4)
 hline!([g_F(1/2, 1/2, 1, 1, 0)]; l=:dash, label="¹P₁", subplot=3)
 hline!([g_F(1/2, 1/2, 1, 1, 1)]; l=:dash, label="³P₁", subplot=3)
diff --git a/src/MQDT.jl b/src/MQDT.jl
index c984de8..c88b060 100644
--- a/src/MQDT.jl
+++ b/src/MQDT.jl
@@ -25,7 +25,9 @@ export lsQuantumNumbers,
     basisarray,
     matrix_elements,
     get_nu_limits_from_model,
-    single_channel_models
+    get_species_parameters,
+    get_fmodels,
+    single_channel_model
 
 include("general.jl")
 include("boundstates.jl")
diff --git a/src/boundstates.jl b/src/boundstates.jl
index c1d825d..c9bda6a 100644
--- a/src/boundstates.jl
+++ b/src/boundstates.jl
@@ -112,7 +112,7 @@ function nu(N::Number, M::Model, P::Parameters)
 end
 
 function nu(N::Vector, M::Model, P::Parameters)
-    n = Matrix{Float64}(undef, M.size, length(N))
+    n = Matrix{Float64}(undef, length(M), length(N))
     for i in eachindex(N)
         n[:, i] = nu(N[i], M, P)
     end
@@ -232,7 +232,7 @@ end
 
 function rot(N::Number, M::fModel, P::Parameters)
     if iszero(M.angles) # check presence of channel mixing
-        return diagm(ones(M.size))
+        return diagm(ones(length(M)))
     else
         t = M.angles[:, 1]
         i = findall(!iszero, M.angles[:, 2]) # check energy dependence
@@ -241,7 +241,7 @@ function rot(N::Number, M::fModel, P::Parameters)
             n = nu(N, M, P)[j] # find nu relative to threshold
             t = theta(n, M.angles)
         end
-        return rot(t, M.mixing, M.size)
+        return rot(t, M.mixing, length(M))
     end
 end
 
@@ -380,7 +380,7 @@ If overwrite_model_limits is True, we will ignore the model nu limits and always
 """
 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)
+        nu_min_model, nu_max_model = M.range
         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)
@@ -434,8 +434,8 @@ Returns a `BasisArray`, which is a list of `BasisState` instances.
 """
 function basisarray(T::EigenStates, M::fModel)
     c = M.outer_channels
-    p = unique_parity(c)
-    f = good_quantum_number(c)
+    p = parity(M)
+    f = M.F
     lr_list = Vector{Union{Int,Nothing}}(nothing, length(M.core))
     for (idx, lr) in zip(findall(M.core), get_lr(c))
         lr_list[idx] = lr
@@ -445,7 +445,7 @@ function basisarray(T::EigenStates, M::fModel)
 
     B = Vector{BasisState}()
     for i in eachindex(T.n)
-        if !isone(M.size) || lr_list[first_relevant_channel] < T.nu[first_relevant_channel, i]
+        if !isone(length(M)) || lr_list[first_relevant_channel] < T.nu[first_relevant_channel, i]
             ei = T.n[i]
             if abs(ei - round(Int, ei)) < 1e2eps()
                 ei = round(ei)
diff --git a/src/general.jl b/src/general.jl
index b64336e..643acdd 100644
--- a/src/general.jl
+++ b/src/general.jl
@@ -323,42 +323,7 @@ function unique_parity(T::Channels)
     if length(f) == 1
         return f[1]
     else
-        error("unique_parity: channels have different parity")
-    end
-end
-
-"""
-    good_quantum_number(T::lsChannels)
-    good_quantum_number(T::jjChannels)
-    good_quantum_number(T::fjChannels)
-
-Given a ls or jj (fj) channel representation, this function checks whether each channel
-has the same total (hyperfine) angular momentum J (F) and, if true, returns it.
-"""
-function good_quantum_number(T::lsChannels)
-    f = unique(get_J(T))
-    if length(f) == 1
-        return f[1]
-    else
-        error("good_quantum_number: channels have different J")
-    end
-end
-
-function good_quantum_number(T::jjChannels)
-    f = unique(get_J(T))
-    if length(f) == 1
-        return f[1]
-    else
-        error("good_quantum_number: channels have different J")
-    end
-end
-
-function good_quantum_number(T::fjChannels)
-    f = unique(get_F(T))
-    if length(f) == 1
-        return f[1]
-    else
-        error("good_quantum_number: channels have different F")
+        error("Channels have different parity")
     end
 end
 
@@ -379,7 +344,8 @@ See also [`kModel`](@ref)
     fModel(
         species::Symbol,
         name::String,
-        size::Int,
+        F::Float64,
+        range::Tuple{Float64,Float64},
         terms::Vector{String},
         core::Vector{Bool},
         defects::Matrix{Float64},
@@ -400,7 +366,8 @@ Contains all relevant parameters to form the K matrix and calculate the spectrum
 julia> RYDBERG_S0 = fModel(
            :Yb174,
            "S J=0, ν > 2", # fit for states 6s7s upward [Phys. Rev. X 15, 011009 (2025)]
-           6,
+           0,
+           (2, Inf),
            ["6sns 1S0", "4f13 5d 6snl a", "6pnp 1S0", "4f13 5d 6snl b", "6pnp 3P0", "4f13 5d 6snl c"],
            Bool[1, 0, 1, 0, 1, 0],
            [
@@ -438,7 +405,8 @@ julia> RYDBERG_S0 = fModel(
 struct fModel <: Model
     species::Symbol
     name::String
-    size::Int
+    F::Float64
+    range::Tuple{Float64,Float64}
     terms::Vector{String}
     core::Vector{Bool}
     defects::Matrix{Float64}
@@ -450,11 +418,12 @@ struct fModel <: Model
     thresholds_dict::Union{Dict{Union{coreQuantumNumbers,String},Float64},Nothing}
 end
 
-function fModel(species, name, size, terms, core, defects, mixing, angles, inner_channels, outer_channels, unitary)
+function fModel(species, name, F, range, terms, core, defects, mixing, angles, inner_channels, outer_channels, unitary)
     return fModel(
         species,
         name,
-        size,
+        F,
+        range,
         terms,
         core,
         defects,
@@ -537,15 +506,11 @@ end
 # --------------------------------------------------------
 
 function Base.length(T::Model)
-    return T.size
-end
-
-function Base.size(T::Model)
-    return (T.size,)
+    return length(T.terms)
 end
 
 function get_lr(T::fModel)
-    l = zeros(Int, T.size)
+    l = zeros(Int, length(T.terms))
     l[findall(T.core)] = get_lr(T.outer_channels)
     return l
 end
@@ -566,6 +531,16 @@ function get_J(T::fModel)
     return get_J(T.inner_channels)
 end
 
+function parity(T::Model)
+    p_inn = unique_parity(T.inner_channels)
+    p_out = unique_parity(T.outer_channels)
+    if p_inn == p_out
+        return p_out
+    else
+        error("Inner and outer channels have different parity")
+    end
+end
+
 """
 Get the minimum and maximum effective principal quantum number (nu_min, nu_max) from a given Model.
 
@@ -670,35 +645,77 @@ function get_thresholds(M::kModel, P::Parameters)
     return thresholds
 end
 
-function single_channel_models(species::Symbol, l::Integer)
-    @assert l > 0 "l must be positive and nonzero for this function"
-    jr = [l-1/2, l-1/2, l+1/2, l+1/2]
-    jt = [l-1, l, l, l+1]
-    m = Vector{fModel}(undef, 4)
-    for i in eachindex(jt)
-        m[i] = fModel(
-            species,
-            "L=$l, J=$(jt[i]), Jr=$(jr[i]), ν > $(l+1)",
-            1,
-            [""],
-            Bool[1],
-            [0;;],
-            [""],
-            [0;;],
-            jjChannels([jjQuantumNumbers(0.5, 0, 0.5, l, jr[i], jt[i])]),
-            jjChannels([jjQuantumNumbers(0.5, 0, 0.5, l, jr[i], jt[i])]),
-            [1;;],
-        )
+function get_species_parameters(species::Symbol)
+    species_module = getfield(MQDT, species)
+    for nm in names(species_module; all=true)
+        obj = getfield(species_module, nm)
+        if isa(obj, Parameters) && obj.species == species
+            return obj
+        end
+    end
+    return error("No Parameters found for species: $species")
+end
+
+"""
+    get_fmodels(species::Symbol, lr::Int, Jr::Float64, Fc::Float64, F::Float64, param::Parameters)
+
+Given atomic `species` and quantum numbers of a Rydberg channel (`lr`, `Jr`, `Fc`, `F`),
+this function searches for all existing `fModel` in the species module, which includes the specified channel.
+If no matching `fModel` is found, a single-channel `fModel` is generated and returned.
+"""
+function get_fmodels(species::Symbol, lr::Int, Jr::Float64, Fc::Float64, F::Float64, param::Parameters)
+    ic = param.spin
+
+    if ic == 0
+        qn = jjQuantumNumbers(0.5, 0, 0.5, lr, Jr, F)
+    else
+        qn = fjQuantumNumbers(0.5, 0, 0.5, Fc, lr, Jr, F)
+    end
+
+    species_module = getfield(MQDT, species)
+    models = Vector{fModel}()
+    for nm in names(species_module; all=true)
+        obj = getfield(species_module, nm)
+        if !isa(obj, fModel) || (obj.species != species) || (obj.F != F)
+            continue
+        end
+        for channel in obj.outer_channels.i
+            if channel == qn
+                push!(models, obj)
+            end
+        end
+    end
+    if !isempty(models)
+        return models
     end
-    return m
+
+    # else generate single-channel model
+    return [single_channel_model(species, lr, Jr, Fc, F, param)]
 end
 
-function single_channel_models(species::Symbol, l_list::UnitRange{Int64})
-    m = Vector{fModel}()
-    for l in l_list
-        append!(m, single_channel_models(species, l))
+function single_channel_model(species::Symbol, lr::Int, Jr::Float64, Fc::Float64, F::Float64, param::Parameters)
+    sr = 1 / 2
+    Jc = 1 / 2
+    ic = param.spin
+    if !(abs(lr - sr) <= Jr <= (lr + sr)) || !(abs(Jc - ic) <= Fc <= (Jc + ic)) || !(abs(Jr - Fc) <= F <= (Jr + Fc))
+        error("Quantum numbers not allowed: lr=$lr, Jr=$Jr, Fc=$Fc, F=$F")
     end
-    return m
+
+    if lr <= 0
+        error("expected only lr > 0: lr=$lr, Jr=$Jr, Fc=$Fc, F=$F")
+    end
+
+    if param.spin == 0
+        channel = jjChannels([jjQuantumNumbers(0.5, 0, 0.5, lr, Jr, F, F)])
+        name = "SQDT L=$lr, Jr=$Jr, J=$F, ν > $(lr+1)"
+        term = "JJ: L=$lr, Jr=$Jr, J=$F"
+    else
+        channel = fjChannels([fjQuantumNumbers(0.5, 0, 0.5, Fc, lr, Jr, F)])
+        name = "SQDT L=$lr, Jr=$Jr, F=$F, ν > $(lr+1)"
+        term = "FJ: I=$(param.spin), L=$lr, Jr=$Jr, F=$F"
+    end
+
+    return fModel(species, name, F, (lr+1, Inf), [term], Bool[1], [0;;], [""], [0;;], channel, channel, [1;;])
 end
 
 # --------------------------------------------------------
diff --git a/src/parameters/Sr87.jl b/src/parameters/Sr87.jl
index b6cc2da..6eeec9a 100644
--- a/src/parameters/Sr87.jl
+++ b/src/parameters/Sr87.jl
@@ -2,7 +2,7 @@ module Sr87
 
 using ..MQDT: Parameters, fModel, lsChannels, fjChannels, lsQuantumNumbers, fjQuantumNumbers, coreQuantumNumbers
 
-export PARA,
+export PARAMETERS,
     FMODEL_HIGHN_S35,
     FMODEL_HIGHN_S45,
     FMODEL_HIGHN_S55,
@@ -24,7 +24,7 @@ export PARA,
 # Isotope data
 THRESHOLDS = Dict(coreQuantumNumbers(0, 0.5, 4) => 45932.287373577, coreQuantumNumbers(0, 0.5, 5) => 45932.120512528)
 
-PARA = Parameters(
+PARAMETERS = Parameters(
     :Sr87,
     1822.888486192*86.9088774970, # nuclear mass
     4.5, # nuclear spin
@@ -40,14 +40,15 @@ PARA = Parameters(
 FMODEL_HIGHN_S35 = fModel(
     :Sr87,
     "S F=7/2, ν > 11",
-    1,
+    3.5,
+    (11.0, Inf),
     ["5sns 3S1"],
     # ["(6s1/2)(ns1/2)"],
     Bool[1],
     [3.370778 0.418 -0.3;],
     [""],
     [0;;],
-    lsChannels([lsQuantumNumbers(0.5, 1, 0, 0, 0, 1)]),
+    lsChannels([lsQuantumNumbers(0.5, 1, 0, 0, 0, 1, 3.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 4, 0, 0.5, 3.5)]),
     [-1;;],
 )
@@ -55,14 +56,15 @@ FMODEL_HIGHN_S35 = fModel(
 FMODEL_HIGHN_S45 = fModel(
     :Sr87,
     "S F=9/2, ν > 11",
-    2,
+    4.5,
+    (11.0, Inf),
     ["5sns 1S0", "5sns 3S1"],
     # ["(6s1/2)(ns1/2)", "(6s1/2)(ns1/2)"],
     Bool[1, 1],
     [3.26896 -0.138 0.9; 3.370778 0.418 -0.3],
     [""],
     [0;;],
-    lsChannels([lsQuantumNumbers(0.5, 0, 0, 0, 0, 0), lsQuantumNumbers(0.5, 1, 0, 0, 0, 1)]),
+    lsChannels([lsQuantumNumbers(0.5, 0, 0, 0, 0, 0, 4.5), lsQuantumNumbers(0.5, 1, 0, 0, 0, 1, 4.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 4, 0, 0.5, 4.5), fjQuantumNumbers(0.5, 0, 0.5, 5, 0, 0.5, 4.5)]),
     [3/2/sqrt(5) sqrt(11)/2/sqrt(5); sqrt(11)/2/sqrt(5) -3/2/sqrt(5)],
 )
@@ -70,22 +72,24 @@ FMODEL_HIGHN_S45 = fModel(
 FMODEL_HIGHN_S55 = fModel(
     :Sr87,
     "S F=11/2, ν > 11",
-    1,
+    5.5,
+    (11.0, Inf),
     ["5sns 3S1"],
     # ["(6s1/2)(ns1/2)"],
     Bool[1],
     [3.370778 0.418 -0.3;],
     [""],
     [0;;],
-    lsChannels([lsQuantumNumbers(0.5, 1, 0, 0, 0, 1)]),
+    lsChannels([lsQuantumNumbers(0.5, 1, 0, 0, 0, 1, 5.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 5, 0, 0.5, 5.5)]),
     [1;;],
 )
 
 FMODEL_LOWN_P45 = fModel(
     :Sr87,
-    "clock, 1.8 < ν < 2.2",
-    4,
+    "P F=9/2 (clock), 1.8 < ν < 2.2",
+    4.5,
+    (1.8, 2.2),
     ["5snp 1P1", "5snp 3P0", "5snp 3P1", "5snp 3P2"],
     # ["(6s1/2)(np1/2)", "(6s1/2)(np1/2)", "(6s1/2)(np3/2)", "(6s1/2)(np3/2)"],
     Bool[1, 1, 1, 1],
@@ -93,10 +97,10 @@ FMODEL_LOWN_P45 = fModel(
     [""],
     [0;;],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 0),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2),
+        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 0, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 4.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 4, 1, 0.5, 4.5),
@@ -115,14 +119,15 @@ FMODEL_LOWN_P45 = fModel(
 FMODEL_HIGHN_P25 = fModel(
     :Sr87,
     "P F=5/2, ν > 5",
-    1,
+    2.5,
+    (5.0, Inf),
     ["5snp 3P2"],
     # ["(6s1/2)(np3/2)"],
     Bool[1],
     [2.882 0.446 -1.9;],
     [""],
     [0;;],
-    lsChannels([lsQuantumNumbers(0.5, 1, 0, 1, 1, 2)]),
+    lsChannels([lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 2.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 4, 1, 1.5, 2.5)]),
     [-1;;],
 )
@@ -130,7 +135,8 @@ FMODEL_HIGHN_P25 = fModel(
 FMODEL_HIGHN_P35 = fModel(
     :Sr87,
     "P F=7/2, ν > 5",
-    3,
+    3.5,
+    (5.0, Inf),
     ["5snp 1P1", "5snp 3P1", "5snp 3P2"],
     # ["(6s1/2)(np1/2)", "(6s1/2)(np3/2)", "(6s1/2)(np3/2)"],
     Bool[1, 1, 1],
@@ -138,9 +144,9 @@ FMODEL_HIGHN_P35 = fModel(
     [""],
     [0;;],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2),
+        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1, 3.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1, 3.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 3.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 4, 1, 0.5, 3.5),
@@ -153,7 +159,8 @@ FMODEL_HIGHN_P35 = fModel(
 FMODEL_HIGHN_P45 = fModel(
     :Sr87,
     "P F=9/2, ν > 7",
-    4,
+    4.5,
+    (7.0, Inf),
     ["5snp 1P1", "5snp 3P0", "5snp 3P1", "5snp 3P2"],
     # ["(6s1/2)(np1/2)", "(6s1/2)(np1/2)", "(6s1/2)(np3/2)", "(6s1/2)(np3/2)"],
     Bool[1, 1, 1, 1],
@@ -161,10 +168,10 @@ FMODEL_HIGHN_P45 = fModel(
     [""],
     [0;;],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 0),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2),
+        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 0, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 4.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 4, 1, 0.5, 4.5),
@@ -183,7 +190,8 @@ FMODEL_HIGHN_P45 = fModel(
 FMODEL_HIGHN_P55 = fModel(
     :Sr87,
     "P F=11/2, ν > 5",
-    3,
+    5.5,
+    (5.0, Inf),
     ["5snp 1P1", "5snp 3P1", "5snp 3P2"],
     # ["(6s1/2)(np3/2)", "(6s1/2)(np1/2)", "(6s1/2)(np3/2)"],
     Bool[1, 1, 1],
@@ -191,9 +199,9 @@ FMODEL_HIGHN_P55 = fModel(
     [""],
     [0;;],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2),
+        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1, 5.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1, 5.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 5.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 4, 1, 1.5, 5.5),
@@ -206,14 +214,15 @@ FMODEL_HIGHN_P55 = fModel(
 FMODEL_HIGHN_P65 = fModel(
     :Sr87,
     "P F=13/2, ν > 5",
-    1,
+    6.5,
+    (5.0, Inf),
     ["5snp 3P2"],
     # ["(6s1/2)(np3/2)"],
     Bool[1],
     [2.882 0.446 -1.9;],
     [""],
     [0;;],
-    lsChannels([lsQuantumNumbers(0.5, 1, 0, 1, 1, 2)]),
+    lsChannels([lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 6.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 5, 1, 1.5, 6.5)]),
     [1;;],
 )
@@ -221,14 +230,15 @@ FMODEL_HIGHN_P65 = fModel(
 FMODEL_HIGHN_D15 = fModel(
     :Sr87,
     "D F=3/2, ν > 47",
-    1,
+    1.5,
+    (47.0, Inf),
     ["5snd 3D3"],
     # ["(6s1/2)(nd5/2)"],
     Bool[1],
     [2.655 -41.4 -15363;],
     [""],
     [0;;],
-    lsChannels([lsQuantumNumbers(0.5, 1, 0, 2, 2, 3)]),
+    lsChannels([lsQuantumNumbers(0.5, 1, 0, 2, 2, 3, 1.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 4, 2, 2.5, 1.5)]),
     [-1;;],
 )
@@ -236,7 +246,8 @@ FMODEL_HIGHN_D15 = fModel(
 FMODEL_HIGHN_D25 = fModel(
     :Sr87,
     "D F=5/2, ν > 47",
-    3,
+    2.5,
+    (47.0, Inf),
     ["5snd 1D2", "5snd 3D2", "5snd 3D3"],
     # ["(6s1/2)(nd3/2)", "(6s1/2)(nd5/2)", "(6s1/2)(nd5/2)"],
     Bool[1, 1, 1],
@@ -244,9 +255,9 @@ FMODEL_HIGHN_D25 = fModel(
     ["1.2"],
     [-0.14 0;],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 2, 2, 2),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 2),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 3),
+        lsQuantumNumbers(0.5, 0, 0, 2, 2, 2, 2.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 2, 2.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 3, 2.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 4, 2, 1.5, 2.5),
@@ -259,7 +270,8 @@ FMODEL_HIGHN_D25 = fModel(
 FMODEL_HIGHN_D35 = fModel(
     :Sr87,
     "D F=7/2, ν > 47",
-    4,
+    3.5,
+    (47.0, Inf),
     ["5snd 1D2", "5snd 3D1", "5snd 3D2", "5snd 3D3"],
     # ["(6s1/2)(nd3/2)", "(6s1/2)(nd5/2)", "(6s1/2)(nd3/2)", "(6s1/2)(nd5/2)"],
     Bool[1, 1, 1, 1],
@@ -267,10 +279,10 @@ FMODEL_HIGHN_D35 = fModel(
     ["1.3"],
     [-0.14 0;],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 2, 2, 2),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 1),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 2),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 3),
+        lsQuantumNumbers(0.5, 0, 0, 2, 2, 2, 3.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 1, 3.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 2, 3.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 3, 3.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 4, 2, 1.5, 3.5),
@@ -289,7 +301,8 @@ FMODEL_HIGHN_D35 = fModel(
 FMODEL_HIGHN_D45 = fModel(
     :Sr87,
     "D F=9/2, ν > 47",
-    4,
+    4.5,
+    (47.0, Inf),
     ["5snd 1D2", "5snd 3D1", "5snd 3D2", "5snd 3D3"],
     # ["(6s1/2)(nd3/2)", "(6s1/2)(nd5/2)", "(6s1/2)(nd3/2)", "(6s1/2)(nd5/2)"],
     Bool[1, 1, 1, 1],
@@ -297,10 +310,10 @@ FMODEL_HIGHN_D45 = fModel(
     ["1.3"],
     [-0.14 0;],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 2, 2, 2),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 1),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 2),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 3),
+        lsQuantumNumbers(0.5, 0, 0, 2, 2, 2, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 1, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 2, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 3, 4.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 4, 2, 1.5, 4.5),
@@ -319,7 +332,8 @@ FMODEL_HIGHN_D45 = fModel(
 FMODEL_HIGHN_D55 = fModel(
     :Sr87,
     "D F=11/2, ν > 47",
-    4,
+    5.5,
+    (47.0, Inf),
     ["5snd 1D2", "5snd 3D1", "5snd 3D2", "5snd 3D3"],
     # ["(6s1/2)(nd3/2)", "(6s1/2)(nd5/2)", "(6s1/2)(nd3/2)", "(6s1/2)(nd5/2)"],
     Bool[1, 1, 1, 1],
@@ -327,10 +341,10 @@ FMODEL_HIGHN_D55 = fModel(
     ["1.3"],
     [-0.14 0;],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 2, 2, 2),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 1),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 2),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 3),
+        lsQuantumNumbers(0.5, 0, 0, 2, 2, 2, 5.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 1, 5.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 2, 5.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 3, 5.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 4, 2, 1.5, 5.5),
@@ -349,7 +363,8 @@ FMODEL_HIGHN_D55 = fModel(
 FMODEL_HIGHN_D65 = fModel(
     :Sr87,
     "D F=13/2, ν > 47",
-    3,
+    6.5,
+    (47.0, Inf),
     ["5snd 1D2", "5snd 3D2", "5snd 3D3"],
     # ["(6s1/2)(nd3/2)", "(6s1/2)(nd5/2)", "(6s1/2)(nd5/2)"],
     Bool[1, 1, 1],
@@ -357,9 +372,9 @@ FMODEL_HIGHN_D65 = fModel(
     ["1.2"],
     [-0.14 0;],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 2, 2, 2),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 2),
-        lsQuantumNumbers(0.5, 1, 0, 2, 2, 3),
+        lsQuantumNumbers(0.5, 0, 0, 2, 2, 2, 6.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 2, 6.5),
+        lsQuantumNumbers(0.5, 1, 0, 2, 2, 3, 6.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 4, 2, 2.5, 6.5),
@@ -372,14 +387,15 @@ FMODEL_HIGHN_D65 = fModel(
 FMODEL_HIGHN_D75 = fModel(
     :Sr87,
     "D F=15/2, ν > 47",
-    1,
+    7.5,
+    (47.0, Inf),
     ["5snd 3D3"],
     # ["(6s1/2)(nd5/2)"],
     Bool[1],
     [2.655 -41.4 -15363;],
     [""],
     [0;;],
-    lsChannels([lsQuantumNumbers(0.5, 1, 0, 2, 2, 3)]),
+    lsChannels([lsQuantumNumbers(0.5, 1, 0, 2, 2, 3, 7.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 5, 2, 2.5, 7.5)]),
     [1;;],
 )
@@ -387,7 +403,8 @@ FMODEL_HIGHN_D75 = fModel(
 FMODEL_HIGHN_F45 = fModel(
     :Sr87,
     "F F=9/2, ν > 9",
-    4,
+    4.5,
+    (9.0, Inf),
     ["5snf 1F3", "5snf 3F2", "5snf 3F3", "5snf 3F4"],
     # ["(6s1/2)(nf5/2)", "(6s1/2)(nf5/2)", "(6s1/2)(nf7/2)", "(6s1/2)(nf7/2)"],
     Bool[1, 1, 1, 1],
@@ -395,10 +412,10 @@ FMODEL_HIGHN_F45 = fModel(
     [""],
     [0;;],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 3, 3, 3),
-        lsQuantumNumbers(0.5, 1, 0, 3, 3, 2),
-        lsQuantumNumbers(0.5, 1, 0, 3, 3, 3),
-        lsQuantumNumbers(0.5, 1, 0, 3, 3, 4),
+        lsQuantumNumbers(0.5, 0, 0, 3, 3, 3, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 3, 3, 2, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 3, 3, 3, 4.5),
+        lsQuantumNumbers(0.5, 1, 0, 3, 3, 4, 4.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 4, 3, 2.5, 4.5),
diff --git a/src/parameters/Sr88.jl b/src/parameters/Sr88.jl
index 838e223..586e6b6 100644
--- a/src/parameters/Sr88.jl
+++ b/src/parameters/Sr88.jl
@@ -2,7 +2,7 @@ module Sr88
 
 using ..MQDT: Parameters, fModel, kModel, lsChannels, jjChannels, lsQuantumNumbers, jjQuantumNumbers, coreQuantumNumbers
 
-export PARA,
+export PARAMETERS,
     FMODEL_LOWN_P1,
     FMODEL_HIGHN_S0,
     FMODEL_HIGHN_S1,
@@ -19,7 +19,7 @@ export PARA,
 # Isotope data
 THRESHOLDS = Dict(coreQuantumNumbers(0, 0.5) => 45932.1956,)
 
-PARA = Parameters(
+PARAMETERS = Parameters(
     :Sr88,
     1822.888486192*87.9056122571, # nuclear mass
     0, # nuclear spin
@@ -35,7 +35,8 @@ PARA = Parameters(
 FMODEL_HIGHN_S0 = fModel(
     :Sr88,
     "S J=0, ν > 10",
-    1,
+    0,
+    (10.0, Inf),
     ["5sns 1S0"],
     Bool[1],
     [3.26896 -0.138 0.9;],
@@ -50,6 +51,7 @@ FMODEL_HIGHN_S1 = fModel(
     :Sr88,
     "S J=1, ν > 11",
     1,
+    (11.0, Inf),
     ["5sns 3S1"],
     Bool[1],
     [3.370778 0.418 -0.3;],
@@ -62,8 +64,9 @@ FMODEL_HIGHN_S1 = fModel(
 
 FMODEL_LOWN_P1 = fModel(
     :Sr88,
-    "recombination, 1.8 < ν < 2.2",
-    2,
+    "P J=1 (recombination), 1.8 < ν < 2.2",
+    1,
+    (1.8, 2.2),
     ["5snp 1P1", "5snp 3P1"],
     Bool[1, 1],
     [0.87199081 0; 0.13140955 0],
@@ -77,7 +80,8 @@ FMODEL_LOWN_P1 = fModel(
 FMODEL_HIGHN_P0 = fModel(
     :Sr88,
     "P J=0, ν > 7",
-    1,
+    0,
+    (7.0, Inf),
     ["5snp 3P0"],
     Bool[1],
     [2.8867 0.44 -1.9],
@@ -91,7 +95,8 @@ FMODEL_HIGHN_P0 = fModel(
 FMODEL_HIGHN_P1 = fModel(
     :Sr88,
     "P J=1, ν > 5",
-    2,
+    1,
+    (5.0, Inf),
     ["5snp 1P1", "5snp 3P1"],
     Bool[1, 1],
     [2.724 -4.67 -157; 2.8826 0.407 -1.3],
@@ -105,7 +110,8 @@ FMODEL_HIGHN_P1 = fModel(
 FMODEL_HIGHN_P2 = fModel(
     :Sr88,
     "P J=2, ν > 5",
-    1,
+    2,
+    (5.0, Inf),
     ["5snp 3P2"],
     Bool[1],
     [2.882 0.446 -1.9],
@@ -120,6 +126,7 @@ FMODEL_HIGHN_D1 = fModel(
     :Sr88,
     "D J=1, ν > 17",
     1,
+    (17.0, Inf),
     ["5snd 3D1"],
     Bool[1],
     [2.67524 -13.15 -4444],
@@ -134,6 +141,7 @@ FMODEL_HIGHN_D2 = fModel(
     :Sr88,
     "D J=2, ν > 25",
     2,
+    (25.0, Inf),
     ["5snd 1D2", "5snd 3D2"],
     Bool[1, 1],
     [2.3847 -39.41 -1090; 2.66149 -16.77 -6656],
@@ -146,8 +154,9 @@ FMODEL_HIGHN_D2 = fModel(
 
 FMODEL_HIGHN_D3 = fModel(
     :Sr88,
-    "D J=3, ν > 47",
-    1,
+    "D J=3, ν > 47", # This is a conservative estimate, based on the considered experimental data. In practise the fit may be reasonable good for lower \nu.
+    3,
+    (47.0, Inf),
     ["5snd 3D3"],
     Bool[1],
     [2.655 -41.4 -15363],
@@ -161,7 +170,8 @@ FMODEL_HIGHN_D3 = fModel(
 FMODEL_HIGHN_F2 = fModel(
     :Sr88,
     "F J=2, ν > 9",
-    1,
+    2,
+    (9.0, Inf),
     ["5snf 3F2"],
     Bool[1],
     [0.12 -2.2 120],
@@ -175,7 +185,8 @@ FMODEL_HIGHN_F2 = fModel(
 FMODEL_HIGHN_F3 = fModel(
     :Sr88,
     "F J=3, ν > 9",
-    2,
+    3,
+    (9.0, Inf),
     ["5snf 1F3", "5snf 3F3"],
     Bool[1, 1],
     [0.089 -2 30; 0.12 -2.2 120],
@@ -189,7 +200,8 @@ FMODEL_HIGHN_F3 = fModel(
 FMODEL_HIGHN_F4 = fModel(
     :Sr88,
     "F J=4, ν > 9",
-    1,
+    4,
+    (9.0, Inf),
     ["5snf 3F4"],
     Bool[1],
     [0.12 -2.2 120],
diff --git a/src/parameters/Yb171.jl b/src/parameters/Yb171.jl
index 5126f1e..c5abb50 100644
--- a/src/parameters/Yb171.jl
+++ b/src/parameters/Yb171.jl
@@ -11,7 +11,7 @@ using ..MQDT:
     fjQuantumNumbers,
     coreQuantumNumbers
 
-export PARA,
+export PARAMETERS,
     FMODEL_HIGHN_S05,
     FMODEL_HIGHN_S15,
     FMODEL_HIGHN_P05,
@@ -50,7 +50,7 @@ THRESHOLDS = Dict(
     "4f13 5d 6s" => 83967.7,
 )
 
-PARA = Parameters(
+PARAMETERS = Parameters(
     :Yb171,
     1822.88848628*170.9363258, # nuclear mass
     0.5, # nuclear spin
@@ -69,7 +69,8 @@ PARA = Parameters(
 FMODEL_HIGHN_S05 = fModel(
     :Yb171,
     "S F=1/2, ν > 26", # [Phys. Rev. A 112, 042817 (2025), Phys. Rev. X 15, 011009 (2025)]
-    7,
+    0.5,
+    (26.0, Inf),
     ["6sns 1S0", "4f13 5d 6snl a", "6pnp 1S0", "4f13 5d 6snl b", "6pnp 3P0", "4f13 5d 6snl c", "6sns 3S1"],
     Bool[1, 0, 1, 0, 1, 0, 1],
     [
@@ -116,7 +117,8 @@ FMODEL_HIGHN_S05 = fModel(
 FMODEL_HIGHN_S15 = fModel(
     :Yb171,
     "S F=3/2, ν > 26", # [Phys. Rev. X 15, 011009 (2025)]
-    1,
+    1.5,
+    (26.0, Inf),
     ["6sns 3S1"],
     Bool[1],
     [0.438426851 3.91762642 -10612.6828 8017432.38 -2582622910.0;],
@@ -130,7 +132,8 @@ FMODEL_HIGHN_S15 = fModel(
 FMODEL_HIGHN_P05 = fModel(
     :Yb171,
     "P F=1/2, ν > 5.7", # [Phys. Rev. A 112, 042817 (2025), Phys. Rev. X 15, 011009 (2025)] # fit for ν > 28, but extrapolates nicely down for ν > 5.7
-    8,
+    0.5,
+    (5.7, Inf),
     [
         "6snp 1P1",
         "6snp 3P1",
@@ -191,7 +194,8 @@ FMODEL_HIGHN_P05 = fModel(
 FMODEL_HIGHN_P15 = fModel(
     :Yb171,
     "P F=3/2, ν > 10", # [Phys. Rev. A 112, 042817 (2025)] # fit for ν > 10
-    11,
+    1.5,
+    (10.0, Inf),
     [
         "6snp 1P1",
         "6snp 3P1",
@@ -265,7 +269,8 @@ FMODEL_HIGHN_P15 = fModel(
 FMODEL_HIGHN_D05 = fModel(
     :Yb171,
     "D F=1/2, ν > 30", # [Phys. Rev. X 15, 011009 (2025)]
-    1,
+    0.5,
+    (30.0, Inf),
     ["6snd 3D1"],
     Bool[1],
     [0.75258093 0.382628525 -483.120633;],
@@ -279,7 +284,8 @@ FMODEL_HIGHN_D05 = fModel(
 FMODEL_HIGHN_D15 = fModel(
     :Yb171,
     "D F=3/2, ν > 30", # [Phys. Rev. A 112, 042817 (2025), Phys. Rev. X 15, 011009 (2025)]
-    6,
+    1.5,
+    (30.0, Inf),
     ["6snd 1D2", "6snd 3D2", "4f13 5d 6snl a", "4f13 5d 6snl b", "6pnp 1D2", "6snd 3D1"],
     Bool[1, 1, 0, 0, 1, 1],
     [
@@ -324,7 +330,8 @@ FMODEL_HIGHN_D15 = fModel(
 FMODEL_HIGHN_D25 = fModel(
     :Yb171,
     "D F=5/2, ν > 30", # [Phys. Rev. A 112, 042817 (2025), Phys. Rev. X 15, 011009 (2025)]
-    6,
+    2.5,
+    (30.0, Inf),
     ["6snd 1D2", "6snd 3D2", "4f13 5d 6snl a", "4f13 5d 6snl b", "6pnp 1D2", "6snd 3D3"],
     Bool[1, 1, 0, 0, 1, 1],
     [
@@ -369,7 +376,8 @@ FMODEL_HIGHN_D25 = fModel(
 FMODEL_HIGHN_D35 = fModel(
     :Yb171,
     "D F=7/2, ν > 14", # fit only valid from 30d upward [Phys. Rev. X 15, 011009 (2025)], provides good match around 18d
-    1,
+    3.5,
+    (14.0, Inf),
     ["6snd 3D3"],
     Bool[1],
     [0.72895315 -0.20653489 220.484722;],
@@ -383,7 +391,8 @@ FMODEL_HIGHN_D35 = fModel(
 FMODEL_HIGHN_F25 = fModel(
     :Yb171,
     "F F=5/2, ν > 20", # [Phys. Rev. A 112, 042817 (2025)]
-    12,
+    2.5,
+    (20.0, Inf),
     [
         "6snf 1F3",
         "6snf 3F3",
@@ -462,7 +471,8 @@ FMODEL_HIGHN_F25 = fModel(
 FMODEL_HIGHN_F35 = fModel(
     :Yb171,
     "F F=7/2, ν > 20", # [Phys. Rev. A 112, 042817 (2025)]
-    8,
+    3.5,
+    (20.0, Inf),
     [
         "6snf 1F3",
         "6snf 3F3",
@@ -523,7 +533,8 @@ FMODEL_HIGHN_F35 = fModel(
 FMODEL_HIGHN_F45 = fModel(
     :Yb171,
     "F F=9/2, ν > 20", # extracted from FMODEL_HIGHN_F35
-    1,
+    4.5,
+    (20.0, Inf),
     ["6snf 3F4"],
     Bool[1],
     [0.0834193873 -1.11453386 -1545.71844;],
@@ -537,7 +548,8 @@ FMODEL_HIGHN_F45 = fModel(
 FMODEL_HIGHN_G25 = fModel(
     :Yb171,
     "G F=5/2, ν > 25", # [Phys. Rev. A 112, 042817 (2025)]
-    1,
+    2.5,
+    (25.0, Inf),
     ["6sng 3G3"],
     Bool[1],
     [0.02613255 -0.14203905;],
@@ -551,7 +563,8 @@ FMODEL_HIGHN_G25 = fModel(
 FMODEL_HIGHN_G35 = fModel(
     :Yb171,
     "G F=7/2, ν > 25", # [Phys. Rev. A 112, 042817 (2025)]
-    3,
+    3.5,
+    (25.0, Inf),
     ["6sng +G4", "6sng -G4", "6sng 3G3"],
     Bool[1, 1, 1],
     [
@@ -581,7 +594,8 @@ FMODEL_HIGHN_G35 = fModel(
 FMODEL_HIGHN_G45 = fModel(
     :Yb171,
     "G F=9/2, ν > 25", # [Phys. Rev. A 112, 042817 (2025)]
-    3,
+    4.5,
+    (25.0, Inf),
     ["6sng +G4", "6sng -G4", "6sng 3G5"],
     Bool[1, 1, 1],
     [
@@ -611,7 +625,8 @@ FMODEL_HIGHN_G45 = fModel(
 FMODEL_HIGHN_G55 = fModel(
     :Yb171,
     "G F=11/2, ν > 25", # extracted from FMODEL_HIGHN_G45
-    1,
+    5.5,
+    (25.0, Inf),
     ["6sng 3G5"],
     Bool[1],
     [0.02536571 -0.18507079;],
@@ -629,7 +644,8 @@ FMODEL_HIGHN_G55 = fModel(
 FMODEL_LOWN_S05 = fModel(
     :Yb171,
     "S F=1/2, 2 < ν < 26", # model from [Phys. Rev. X 15, 011009 (2025)] with modified 3S1 data taken from Yb174 fit to NIST
-    7,
+    0.5,
+    (2.0, 26.0),
     ["6sns 1S0", "4f13 5d 6snl a", "6pnp 1S0", "4f13 5d 6snl b", "6pnp 3P0", "4f13 5d 6snl c", "6sns 3S1"],
     Bool[1, 0, 1, 0, 1, 0, 1],
     [
@@ -676,7 +692,8 @@ FMODEL_LOWN_S05 = fModel(
 FMODEL_LOWN_S15 = fModel(
     :Yb171,
     "S F=3/2, 2 < ν < 26", # taken from Yb174 fit to NIST
-    1,
+    1.5,
+    (2.0, 26.0),
     ["6sns 3S1"],
     Bool[1],
     [0.432841 0.724559 -1.95424;],
@@ -690,7 +707,8 @@ FMODEL_LOWN_S15 = fModel(
 FMODEL_LOWEST_P05 = fModel(
     :Yb171,
     "P F=1/2, 1.5 < ν < 2.5", # fit to Yb174 NIST data
-    3,
+    0.5,
+    (1.5, 2.5),
     ["6snp 1P1", "6snp 3P1", "6snp 3P0"],
     Bool[1, 1, 1],
     [0.161083 0; 0.920424 0; 0.180701 0],
@@ -716,7 +734,8 @@ FMODEL_LOWEST_P05 = fModel(
 FMODEL_LOWN_P05 = fModel(
     :Yb171,
     "P F=1/2, 2.9 < ν < 5.9", # data taken from Yb174 NIST fits
-    3,
+    0.5,
+    (2.9, 5.9),
     ["6snp 1P1", "6snp 3P1", "6snp 3P0"],
     Bool[1, 1, 1],
     [
@@ -746,7 +765,8 @@ FMODEL_LOWN_P05 = fModel(
 FMODEL_LOWEST_P15 = fModel(
     :Yb171,
     "P F=3/2, 1.5 < ν < 2.5", # fit to Yb174 NIST data
-    3,
+    1.5,
+    (1.5, 2.5),
     ["6snp 1P1", "6snp 3P1", "6snp 3P2"],
     Bool[1, 1, 1],
     [0.161083 0; 0.920424 0; 0.110501 0],
@@ -772,7 +792,8 @@ FMODEL_LOWEST_P15 = fModel(
 FMODEL_LOWN_P15 = fModel(
     :Yb171,
     "P F=3/2, 3 < ν < 10", # data taken from Yb174 NIST fits
-    10,
+    1.5,
+    (3.0, 10.0),
     [
         "6snp 1P1",
         "6snp 3P1",
@@ -840,7 +861,8 @@ FMODEL_LOWN_P15 = fModel(
 FMODEL_LOWEST_P25 = fModel(
     :Yb171,
     "P F=5/2, 1.5 < ν < 4.5", # data taken from Yb174 NIST fit
-    1,
+    2.5,
+    (1.5, 4.5),
     ["6snp 3P2"],
     Bool[1],
     [0.906105 0.383471 1.23512;],
@@ -854,7 +876,8 @@ FMODEL_LOWEST_P25 = fModel(
 FMODEL_LOWN_P25 = fModel(
     :Yb171,
     "P F=5/2, 5 < ν < 20", # data taken from Yb174 [Phys. Rev. X 15, 011009 (2025)] fits
-    4,
+    2.5,
+    (5.0, 20.0),
     ["6snp 3P2", "4f13 5d 6snl a", "4f13 5d 6snl b", "4f13 5d 6snl c"],
     Bool[1, 0, 0, 0],
     [
@@ -882,7 +905,8 @@ FMODEL_LOWN_P25 = fModel(
 FMODEL_LOWN_D05 = fModel(
     :Yb171,
     "D F=1/2, 2 < ν < 30", # data taken from Yb174 NIST fit for 2 < ν < 5
-    1,
+    0.5,
+    (2.0, 30.0),
     ["6snd 3D1"],
     Bool[1],
     [0.758222 -0.017906 3.392161;],
@@ -896,7 +920,8 @@ FMODEL_LOWN_D05 = fModel(
 FMODEL_LOWN_D15 = fModel(
     :Yb171,
     "D F=3/2, 2 < ν < 30", # model from [arXiv:2507.11487v1] with modified 3D1 data taken from Yb174 fit to NIST
-    6,
+    1.5,
+    (2.0, 30.0),
     ["6snd 1D2", "6snd 3D2", "4f13 5d 6snl a", "4f13 5d 6snl b", "6pnp 1D2", "6snd 3D1"],
     Bool[1, 1, 0, 0, 1, 1],
     [
@@ -941,7 +966,8 @@ FMODEL_LOWN_D15 = fModel(
 FMODEL_LOWN_D25 = fModel(
     :Yb171,
     "D F=5/2, 2 < ν < 30", # model from [arXiv:2507.11487v1] with modified 3D3 data taken from Yb174 fit to NIST
-    6,
+    2.5,
+    (2.0, 30.0),
     ["6snd 1D2", "6snd 3D2", "4f13 5d 6snl a", "4f13 5d 6snl b", "6pnp 1D2", "6snd 3D3"],
     Bool[1, 1, 0, 0, 1, 1],
     [
@@ -986,7 +1012,8 @@ FMODEL_LOWN_D25 = fModel(
 FMODEL_LOWN_D35 = fModel(
     :Yb171,
     "D F=7/2, 2 < ν < 14", # data taken from Yb174 NIST fit for 2 < ν < 5, provides good match around 18d
-    1,
+    3.5,
+    (2.0, 14.0),
     ["6snd 3D3"],
     Bool[1],
     [0.734512 -0.019501 3.459114;],
diff --git a/src/parameters/Yb173.jl b/src/parameters/Yb173.jl
index 3a24c9e..df8f86c 100644
--- a/src/parameters/Yb173.jl
+++ b/src/parameters/Yb173.jl
@@ -2,7 +2,7 @@ module Yb173
 
 using ..MQDT: Parameters, fModel, lsChannels, fjChannels, lsQuantumNumbers, fjQuantumNumbers, coreQuantumNumbers
 
-export PARA,
+export PARAMETERS,
     FMODEL_HIGHN_S15,
     FMODEL_HIGHN_S25,
     FMODEL_HIGHN_S35,
@@ -21,7 +21,7 @@ THRESHOLDS = Dict(
     "4f13 5d 6s" => 83967.7,
 )
 
-PARA = Parameters(
+PARAMETERS = Parameters(
     :Yb173,
     1822.88848628*172.938216212, # nuclear mass
     2.5, # nuclear spin
@@ -43,14 +43,15 @@ PARA = Parameters(
 FMODEL_HIGHN_S15 = fModel(
     :Yb173,
     "S F=3/2, ν > 26",
-    1,
+    1.5,
+    (26.0, Inf),
     ["6sns 3S1"],
     # ["(6s1/2)(ns1/2)"],
     Bool[1],
     [0.438426851 3.91762642 -10612.6828 8017432.38 -2582622910.0],
     [""],
     [0;;],
-    lsChannels([lsQuantumNumbers(0.5, 1, 0, 0, 0, 1)]),
+    lsChannels([lsQuantumNumbers(0.5, 1, 0, 0, 0, 1, 1.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 2, 0, 0.5, 1.5)]),
     [-1;;],
 )
@@ -58,7 +59,8 @@ FMODEL_HIGHN_S15 = fModel(
 FMODEL_HIGHN_S25 = fModel(
     :Yb173,
     "S F=5/2, ν > 26",
-    7,
+    2.5,
+    (26.0, Inf),
     ["6sns 1S0", "4f13 5d 6snl a", "6pnp 1S0", "4f13 5d 6snl b", "6pnp 3P0", "4f13 5d 6snl c", "6sns 3S1"],
     # ["(6s1/2)(ns1/2)", "5d a", "(6p3/2)(np3/2)", "5d b", "(6p1/2)(np1/2)", "5d c", "(6s1/2)(ns1/2)"],
     Bool[1, 0, 1, 0, 1, 0, 1],
@@ -81,10 +83,10 @@ FMODEL_HIGHN_S25 = fModel(
         0.137723736 0
     ],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 0, 0, 0),
-        lsQuantumNumbers(0.5, 0, 1, 1, 0, 0),
-        lsQuantumNumbers(0.5, 1, 1, 1, 1, 0),
-        lsQuantumNumbers(0.5, 1, 0, 0, 0, 1),
+        lsQuantumNumbers(0.5, 0, 0, 0, 0, 0, 2.5),
+        lsQuantumNumbers(0.5, 0, 1, 1, 0, 0, 2.5),
+        lsQuantumNumbers(0.5, 1, 1, 1, 1, 0, 2.5),
+        lsQuantumNumbers(0.5, 1, 0, 0, 0, 1, 2.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 2, 0, 0.5, 2.5),
@@ -106,14 +108,15 @@ FMODEL_HIGHN_S25 = fModel(
 FMODEL_HIGHN_S35 = fModel(
     :Yb173,
     "S F=7/2, ν > 26",
-    1,
+    3.5,
+    (26.0, Inf),
     ["6sns 3S1"],
     # ["(6s1/2)(ns1/2)"],
     Bool[1],
     [0.438426851 3.91762642 -10612.6828 8017432.38 -2582622910.0],
     [""],
     [0;;],
-    lsChannels([lsQuantumNumbers(0.5, 1, 0, 0, 0, 1)]),
+    lsChannels([lsQuantumNumbers(0.5, 1, 0, 0, 0, 1, 3.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 3, 0, 0.5, 3.5)]),
     [1;;],
 )
@@ -121,7 +124,8 @@ FMODEL_HIGHN_S35 = fModel(
 FMODEL_HIGHN_P05 = fModel(
     :Yb173,
     "P F=1/2, ν > 10",
-    4,
+    0.5,
+    (10.0, Inf),
     ["6snp 3P2", "4f13 5d 6snl a", "4f13 5d 6snl b", "4f13 5d 6snl c"],
     # ["(6s1/2)(np3/2)", "5d a", "5d b", "5d c"],
     Bool[1, 0, 0, 0],
@@ -133,7 +137,7 @@ FMODEL_HIGHN_P05 = fModel(
     ],
     ["1.2", "1.3", "1.4"],
     [0.071426685 0; 0.027464110 0; -0.029741862 0],
-    lsChannels([lsQuantumNumbers(0.5, 1, 0, 1, 1, 2)]),
+    lsChannels([lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 0.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 2, 1, 1.5, 0.5)]),
     [-1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1],
 )
@@ -141,7 +145,8 @@ FMODEL_HIGHN_P05 = fModel(
 FMODEL_HIGHN_P15 = fModel(
     :Yb173,
     "P F=3/2, ν > 10",
-    10,
+    1.5,
+    (10.0, Inf),
     [
         "6snp 1P1",
         "6snp 3P1",
@@ -184,9 +189,9 @@ FMODEL_HIGHN_P15 = fModel(
         -0.029741862 0 0
     ],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2),
+        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1, 1.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1, 1.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 1.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 2, 1, 0.5, 1.5),
@@ -210,7 +215,8 @@ FMODEL_HIGHN_P15 = fModel(
 FMODEL_HIGHN_P25 = fModel(
     :Yb173,
     "P F=5/2, ν > 10",
-    12,
+    2.5,
+    (10.0, Inf),
     [
         "6snp 1P1",
         "6snp 3P1",
@@ -258,10 +264,10 @@ FMODEL_HIGHN_P25 = fModel(
         0.029741862 0 0
     ],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 0),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2),
+        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1, 2.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1, 2.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 0, 2.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 2.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 2, 1, 0.5, 2.5),
@@ -288,7 +294,8 @@ FMODEL_HIGHN_P25 = fModel(
 FMODEL_HIGHN_P35 = fModel(
     :Yb173,
     "P F=7/2, ν > 10",
-    10,
+    3.5,
+    (10.0, Inf),
     [
         "6snp 1P1",
         "6snp 3P1",
@@ -331,9 +338,9 @@ FMODEL_HIGHN_P35 = fModel(
         -0.029741862 0 0
     ],
     lsChannels([
-        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1),
-        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2),
+        lsQuantumNumbers(0.5, 0, 0, 1, 1, 1, 3.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 1, 3.5),
+        lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 3.5),
     ]),
     fjChannels([
         fjQuantumNumbers(0.5, 0, 0.5, 2, 1, 1.5, 3.5),
@@ -357,7 +364,8 @@ FMODEL_HIGHN_P35 = fModel(
 FMODEL_HIGHN_P45 = fModel(
     :Yb173,
     "P F=9/2, ν > 10",
-    4,
+    4.5,
+    (10.0, Inf),
     ["6snp 3P2", "4f13 5d 6snl a", "4f13 5d 6snl b", "4f13 5d 6snl c"],
     # ["(6s1/2)(np3/2)", "5d a", "5d b", "5d c"],
     Bool[1, 0, 0, 0],
@@ -369,7 +377,7 @@ FMODEL_HIGHN_P45 = fModel(
     ],
     ["1.2", "1.3", "1.4"],
     [0.071426685 0; 0.027464110 0; -0.029741862 0],
-    lsChannels([lsQuantumNumbers(0.5, 1, 0, 1, 1, 2)]),
+    lsChannels([lsQuantumNumbers(0.5, 1, 0, 1, 1, 2, 4.5)]),
     fjChannels([fjQuantumNumbers(0.5, 0, 0.5, 3, 1, 1.5, 4.5)]),
     [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1],
 )
diff --git a/src/parameters/Yb174.jl b/src/parameters/Yb174.jl
index f79b57f..b560d31 100644
--- a/src/parameters/Yb174.jl
+++ b/src/parameters/Yb174.jl
@@ -2,7 +2,7 @@ module Yb174
 
 using ..MQDT: Parameters, fModel, lsChannels, jjChannels, lsQuantumNumbers, jjQuantumNumbers, coreQuantumNumbers
 
-export PARA,
+export PARAMETERS,
     FMODEL_HIGHN_S0,
     FMODEL_HIGHN_S1,
     FMODEL_HIGHN_P0,
@@ -36,7 +36,7 @@ THRESHOLDS = Dict(
     "4f13 5d 6s" => 83967.7,
 )
 
-PARA = Parameters(
+PARAMETERS = Parameters(
     :Yb174,
     1822.88848628*173.9388621, # nuclear mass
     0, # nuclear spin
@@ -54,7 +54,8 @@ PARA = Parameters(
 FMODEL_HIGHN_S0 = fModel(
     :Yb174,
     "S J=0, ν > 2", # fit for states 6s7s upward [Phys. Rev. X 15, 011009 (2025)]
-    6,
+    0,
+    (2.0, Inf),
     ["6sns 1S0", "4f13 5d 6snl a", "6pnp 1S0", "4f13 5d 6snl b", "6pnp 3P0", "4f13 5d 6snl c"],
     Bool[1, 0, 1, 0, 1, 0],
     [
@@ -98,6 +99,7 @@ FMODEL_HIGHN_S1 = fModel(
     :Yb174,
     "S J=1, ν > 26", # fit only valid from 28s upward [Phys. Rev. Lett. 128, 033201 (2022)]
     1,
+    (26.0, Inf),
     ["6sns 3S1"],
     Bool[1],
     [0.4382 4 -1e4 8e6 -3e9;], # adjusted for non-rydberg ritz
@@ -111,7 +113,8 @@ FMODEL_HIGHN_S1 = fModel(
 FMODEL_HIGHN_P0 = fModel(
     :Yb174,
     "P J=0, ν > 6", # [Phys. Rev. X 15, 011009 (2025)]
-    2,
+    0,
+    (6.0, Inf),
     ["6snp 3P0", "4f13 5d 6snd"],
     Bool[1, 0],
     [0.953661478 -0.287531374; 0.198460766 0],
@@ -125,7 +128,8 @@ FMODEL_HIGHN_P0 = fModel(
 FMODEL_HIGHN_P1 = fModel(
     :Yb174,
     "P J=1, ν > 6", # [Phys. Rev. X 15, 011009 (2025)]
-    6,
+    1,
+    (6.0, Inf),
     ["6snp 1P1", "6snp 3P1", "4f13 5d 6snl a", "4f13 5d 6snl b", "4f13 5d 6snl c", "4f13 5d 6snl d"],
     Bool[1, 1, 0, 0, 0, 0],
     [
@@ -163,7 +167,8 @@ FMODEL_HIGHN_P1 = fModel(
 FMODEL_HIGHN_P2 = fModel(
     :Yb174,
     "P J=2, ν > 5", # [Phys. Rev. X 15, 011009 (2025)]
-    4,
+    2,
+    (5.0, Inf),
     ["6snp 3P2", "4f13 5d 6snl a", "4f13 5d 6snl b", "4f13 5d 6snl c"],
     Bool[1, 0, 0, 0],
     [
@@ -187,6 +192,7 @@ FMODEL_HIGHN_D1 = fModel(
     :Yb174,
     "D J=1, ν > 26", # fit only valid from 30d upward [Phys. Rev. X 15, 011009 (2025)]
     1,
+    (26.0, Inf),
     ["6snd 3D1"],
     Bool[1],
     [0.75258093 0.3826 -483.1;],
@@ -200,7 +206,8 @@ FMODEL_HIGHN_D1 = fModel(
 FMODEL_HIGHN_D2 = fModel(
     :Yb174,
     "D J=2, ν > 5", # [Phys. Rev. X 15, 011009 (2025)]
-    5,
+    2,
+    (5.0, Inf),
     ["6snd 1D2", "6snd 3D2", "4f13 5d 6snl a", "4f13 5d 6snl b", "6pnp 1D2"],
     Bool[1, 1, 0, 0, 1],
     [
@@ -241,7 +248,8 @@ FMODEL_HIGHN_D2 = fModel(
 FMODEL_HIGHN_D3 = fModel(
     :Yb174,
     "D J=3, ν > 18", # fit only valid from 30d upward [Phys. Rev. X 15, 011009 (2025)], provides good match around 21d
-    1,
+    3,
+    (18.0, Inf),
     ["6snd 3D3"],
     Bool[1],
     [0.72902016 -0.705328923 829.238844;],
@@ -255,7 +263,8 @@ FMODEL_HIGHN_D3 = fModel(
 FMODEL_HIGHN_F2 = fModel(
     :Yb174,
     "F J=2, ν > 25", # [arXiv:2507.11487v1]
-    1,
+    2,
+    (25.0, Inf),
     ["6snf 3F2"],
     Bool[1],
     [0.0718252326 -1.00091963 -106.291066;],
@@ -269,7 +278,8 @@ FMODEL_HIGHN_F2 = fModel(
 FMODEL_HIGHN_F3 = fModel(
     :Yb174,
     "F J=3, ν > 7", # [Phys. Rev. A 112, 042817 (2025)]
-    7,
+    3,
+    (7.0, Inf),
     ["6snf 1F3", "6snf 3F3", "4f13 5d 6snl a", "4f13 5d 6snl b", "4f13 5d 6snl c", "4f13 5d 6snl d", "4f13 5d 6snl e"],
     Bool[1, 1, 0, 0, 0, 0, 0],
     [
@@ -311,7 +321,8 @@ FMODEL_HIGHN_F3 = fModel(
 FMODEL_HIGHN_F4 = fModel(
     :Yb174,
     "F J=4, ν > 25", # [arXiv:2507.11487v1]
-    1,
+    4,
+    (25.0, Inf),
     ["6snf 3F4"],
     Bool[1],
     [0.0839027969 -2.91009023;],
@@ -325,7 +336,8 @@ FMODEL_HIGHN_F4 = fModel(
 FMODEL_HIGHN_G3 = fModel(
     :Yb174,
     "G J=3, ν > 25", # [arXiv:2507.11487v1]
-    1,
+    3,
+    (25.0, Inf),
     ["6sng 3G3"],
     Bool[1],
     [0.0260964574 -0.14139526;],
@@ -339,7 +351,8 @@ FMODEL_HIGHN_G3 = fModel(
 FMODEL_HIGHN_G4 = fModel(
     :Yb174,
     "G J=4, ν > 25", # [Phys. Rev. A 112, 042817 (2025)]
-    2,
+    4,
+    (25.0, Inf),
     ["6sng +G4", "6sng -G4"],
     Bool[1, 1],
     [0.0262659964 -0.148808463; 0.0254568575 -0.134219071;],
@@ -353,7 +366,8 @@ FMODEL_HIGHN_G4 = fModel(
 FMODEL_HIGHN_G5 = fModel(
     :Yb174,
     "G J=5, ν > 25", # [Phys. Rev. A 112, 042817 (2025)]
-    1,
+    5,
+    (25.0, Inf),
     ["6sng 3G5"],
     Bool[1],
     [0.02536571 -0.18507079;], # updated since [arXiv:2507.11487v1]
@@ -371,7 +385,8 @@ FMODEL_HIGHN_G5 = fModel(
 FMODEL_LOWN_S0 = fModel(
     :Yb174,
     "S J=0, 1 < ν < 2", # fit to the 6s^2 ground state
-    1,
+    0,
+    (1.0, 2.0),
     ["6sns 1S0"],
     Bool[1],
     [0.525055 0;],
@@ -386,6 +401,7 @@ FMODEL_LOWN_S1 = fModel(
     :Yb174,
     "S J=1, 2 < ν < 26", # fit to NIST data between 7s and 13s, extrapolation seems good up to 30s
     1,
+    (2.0, 26.0),
     ["6sns 3S1"],
     Bool[1],
     [0.432841 0.724559 -1.95424;],
@@ -399,7 +415,8 @@ FMODEL_LOWN_S1 = fModel(
 FMODEL_LOWN_P0 = fModel(
     :Yb174,
     "P J=0, 1.5 < ν < 5.5", # fit to NIST data between 6p and 9p
-    1,
+    0,
+    (1.5, 5.5),
     ["6snp 3P0"],
     Bool[1],
     [0.969279 0.288219 1.36228;],
@@ -413,7 +430,8 @@ FMODEL_LOWN_P0 = fModel(
 FMODEL_LOWEST_P1 = fModel(
     :Yb174,
     "P J=1, 1.7 < ν < 2.7", # fit to NIST data for the 6p states
-    2,
+    1,
+    (1.7, 2.7),
     ["6snp 1P1", "6snp 3P1"],
     Bool[1, 1],
     [0.161083 0; 0.920424 0],
@@ -427,7 +445,8 @@ FMODEL_LOWEST_P1 = fModel(
 FMODEL_LOWN_P1 = fModel(
     :Yb174,
     "P J=1, 2.7 < ν < 5.7", # fit to NIST data between 7p and 9p
-    2,
+    1,
+    (2.7, 5.7),
     ["6snp 1P1", "6snp 3P1"],
     Bool[1, 1],
     [0.967223 -3.03997 0.569205; 0.967918 0.25116 0.868505;],
@@ -441,7 +460,8 @@ FMODEL_LOWN_P1 = fModel(
 FMODEL_LOWN_P2 = fModel(
     :Yb174,
     "P J=2, 1.5 < ν < 4.5", # fit to NIST data between 6p and 8p
-    1,
+    2,
+    (1.5, 4.5),
     ["6snp 3P2"],
     Bool[1],
     [0.906105 0.383471 1.23512;],
@@ -456,6 +476,7 @@ FMODEL_LOWN_D1 = fModel(
     :Yb174,
     "D J=1, 2 < ν < 26", # fit to NIST data between 5d and 8d, causes a µ=0.005 difference at the 30d state
     1,
+    (2.0, 26.0),
     ["6snd 3D1"],
     Bool[1],
     [0.758222 -0.017906 3.392161;],
@@ -470,6 +491,7 @@ FMODEL_LOWN_D2 = fModel(
     :Yb174,
     "D J=2, 2 < ν < 5", # fit to NIST data between 5d and 7d
     2,
+    (2.0, 5.0),
     ["6snd 1D2", "6snd 3D2"],
     Bool[1, 1],
     [0.703156 0.973192; 0.724546 0.372621],
@@ -483,7 +505,8 @@ FMODEL_LOWN_D2 = fModel(
 FMODEL_LOWN_D3 = fModel(
     :Yb174,
     "D J=3, 2 < ν < 18", # fit to NIST data between 5d and 8d, provides good match around 21d
-    1,
+    3,
+    (2.0, 18.0),
     ["6snd 3D3"],
     Bool[1],
     [0.734512 -0.019501 3.459114;],
diff --git a/test/small_test.jl b/test/small_test.jl
index 32d4de0..97a88d3 100644
--- a/test/small_test.jl
+++ b/test/small_test.jl
@@ -1,7 +1,8 @@
 using MQDT
 
 @testset "small test" begin
-    s = MQDT.eigenstates(24, 30, MQDT.Yb171.FMODEL_HIGHN_S15, MQDT.Yb171.PARA)
-    b = MQDT.basisarray([s], [MQDT.Yb171.FMODEL_HIGHN_S15])
-    m = [MQDT.multipole_moments(i, i, MQDT.Yb171.PARA)[4] for i in b.states]
+    s = eigenstates(24, 30, MQDT.Yb171.FMODEL_HIGHN_S15, MQDT.Yb171.PARAMETERS)
+    b = basisarray([s], [MQDT.Yb171.FMODEL_HIGHN_S15])
+    m = [MQDT.multipole_moments(i, i, MQDT.Yb171.PARAMETERS)[3] for i in b.states]
+    @test isapprox(m, 1e6*[2.491, 2.887, 3.329, 3.82], atol=1e3)
 end
diff --git a/test/test_model.jl b/test/test_model.jl
index 594de4e..d742474 100644
--- a/test/test_model.jl
+++ b/test/test_model.jl
@@ -2,10 +2,10 @@ using MQDT
 using PythonCall
 
 function test_model_struct(T::fModel)
-    if (T.size != length(T.terms)) || (T.size != length(T.core)) || (T.size != size(T.defects, 1))
+    if (length(T.terms) != length(T.core)) || (length(T.terms) != size(T.defects, 1))
         println("Model size does not correspond to provided parameters.")
         @test false
-    elseif size(T.unitary) != (T.size, T.size)
+    elseif size(T.unitary) != (length(T.terms), length(T.terms))
         println("Frame transformation matrix does not have the correct dimensions.")
         @test false
     elseif length(findall(T.core)) != size(T.outer_channels) != size(T.inner_channels)
@@ -14,6 +14,12 @@ function test_model_struct(T::fModel)
     elseif length(T.mixing) != size(T.angles, 1)
         println("Wrong number of close-coupling rotations.")
         @test false
+    elseif MQDT.parity(T) != 1 && MQDT.parity(T) != -1
+        println("Channels have inconsistent parity.")
+        @test false
+    elseif [T.F] != unique(MQDT.get_F(T.inner_channels)) || [T.F] != unique(MQDT.get_F(T.outer_channels))
+        println("Channels don't have the same good quantum number.")
+        @test false
     else
         println("Model $(T.name) passed.")
         @test true
@@ -21,14 +27,13 @@ function test_model_struct(T::fModel)
 end
 
 function test_model_struct(T::kModel)
-    if (T.size != length(T.terms)) ||
-       (T.size != length(T.jjscheme)) ||
-       (T.size != length(T.lschannels.i)) ||
-       (T.size != length(T.jjchannels.i)) ||
-       (T.size != length(T.K1))
+    if (length(T.terms) != length(T.jjscheme)) ||
+       (length(T.terms) != length(T.lschannels.i)) ||
+       (length(T.terms) != length(T.jjchannels.i)) ||
+       (length(T.terms) != length(T.K1))
         println("Model size does not correspond to provided parameters.")
         @test false
-    elseif size(T.K0) != (T.size, T.size)
+    elseif size(T.K0) != (length(T.terms), length(T.terms))
         println("Frame transformation matrix does not have the correct dimensions.")
         @test false
     else
@@ -67,6 +72,7 @@ function test_model_unitary(T::fModel, P::Parameters)
             lc = unique(MQDT.get_lc(T.inner_channels))
             lr = unique(MQDT.get_lr(T.inner_channels))
             ft = unique(MQDT.get_F(T.outer_channels))[1]
+            @test ft == T.F
             qn = MQDT.AngularMomenta(lc, lr, P.spin)
             jj = MQDT.jj_channels(qn)
             fj = [unique.(MQDT.get_F.(jj))[i][1] for i in 1:length(jj)]
@@ -89,9 +95,29 @@ function test_model_unitary(T::kModel)
     @test isapprox(t1, t2; rtol=0.001)
 end
 
-function test_model_states_and_matrix_elements(model::fModel, param::Parameters)
+function test_model_name(model::fModel)
+    f_tot_str = match(r"F=([\d]+/?[\d]*)", model.name)
+    if f_tot_str == nothing
+        f_tot_str = match(r"J=([\d]+)", model.name)
+    end
+
+    @test f_tot_str != nothing
+    s = f_tot_str.captures[1]
+    if occursin("/", s)
+        parts = split(s, "/")
+        f_tot = parse(Float64, parts[1]) / parse(Float64, parts[2])
+    else
+        f_tot = parse(Float64, s)
+    end
+    @test f_tot == model.F
+
     nu_min, nu_max = get_nu_limits_from_model(model)
-    nu_max = isnan(nu_max) ? nu_min + 5 : min(nu_max, nu_min + 5)
+    @test (nu_min, nu_max) == model.range
+end
+
+function test_model_states_and_matrix_elements(model::fModel, param::Parameters)
+    nu_min, nu_max = model.range
+    nu_max = min(nu_max, nu_min + 5)
     states = eigenstates(nu_min, nu_max, model, param)
     basis = basisarray(states, model)
     @test size(basis) == length(states)
@@ -113,9 +139,6 @@ end
 @testset "Models" begin
     for species in (MQDT.Sr87, MQDT.Sr88, MQDT.Yb171, MQDT.Yb173, MQDT.Yb174)
         for nm in names(species, all=true)
-            if !isdefined(species, nm)
-                continue
-            end
             obj = getfield(species, nm)
             if !isa(obj, MQDT.Model)
                 continue
@@ -124,9 +147,10 @@ end
             println("Testing $(species) $(nm)")
             @testset "$(species) $(nm)" begin
                 if isa(obj, fModel)
+                    test_model_name(obj)
                     test_model_struct(obj)
-                    test_model_unitary(obj, species.PARA)
-                    test_model_states_and_matrix_elements(obj, species.PARA)
+                    test_model_unitary(obj, species.PARAMETERS)
+                    test_model_states_and_matrix_elements(obj, species.PARAMETERS)
                 elseif isa(obj, kModel)
                     test_model_struct(obj)
                     test_model_unitary(obj)