Improve fmodel

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

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)

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")

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)