Fix stale material-availability tables; add ion-aware lookup and Bethe fallback for DEDX_ICRU

include/dedx.h

@@ -44,7 +44,23 @@ enum {
     DEDX_ICRU73,        /**< ICRU Report 73 (2005) */
     DEDX_ICRU49,        /**< ICRU Report 49 (1993) — protons and alpha particles */
     _DEDX_0008,         /**< Reserved */
-    DEDX_ICRU,          /**< Auto-selects ICRU49 or ICRU73 based on ion type */
+    DEDX_ICRU,          /**< Auto-selects ICRU49 or ICRU73 based on ion type. A real, tabulated
+                             report like the others: it does not fall back to the Bethe-Bloch
+                             formula, so ion/target combinations no ICRU sub-report covers still
+                             fail with DEDX_ERR_COMBINATION_NOT_FOUND. Use DEDX_AUTO for a variant
+                             that also falls back to Bethe-Bloch. */
+    DEDX_AUTO,          /**< Like DEDX_ICRU (auto-selects the best tabulated report for the ion,
+                             preferring directly-tabulated data for a target whenever it exists),
+                             but additionally falls back to the Bethe-Bloch formula (see
+                             DEDX_BETHE_EXT00) for elemental targets no tabulated report covers
+                             (e.g. Boron), rather than failing outright. This mixes a real,
+                             literature-backed dataset with an analytical model as a deliberate,
+                             documented "best effort" hybrid -- comparable to the AUTO mode in
+                             dedx_web -- and is distinct from DEDX_ICRU, which never does this.
+                             Deliberately grouped with the other tabulated report identifiers
+                             (not >= DEDX_DEFAULT) so that a compound target resolves through the
+                             same tabulated-first path as DEDX_ICRU instead of always decomposing
+                             into elements. */
     DEDX_DEFAULT = 100, /**< Default program (Bethe formula) */
     DEDX_BETHE_EXT00    /**< Bethe formula with extensions */
 };
@@ -162,6 +178,44 @@ const int *dedx_get_program_list(void);
  */
 const int *dedx_get_material_list(int program);
 
+/** @brief Recommended capacity for materials/material_list buffers passed to
+ *  dedx_get_material_list_for_ion() or its dedx_fill_material_list_for_ion() wrapper.
+ *  Comfortably exceeds the total number of known materials, so passing this size as
+ *  @p max_materials guarantees the result is never truncated. This is a recommended
+ *  buffer size, not a hard cap enforced by the API: dedx_get_material_list_for_ion()
+ *  writes at most @p max_materials entries, whatever value the caller passes for it.
+ */
+#define DEDX_MAX_MATERIAL_LIST 1000
+
+/** @brief Fill the materials actually usable with a specific program/ion combination.
+ *
+ *  Unlike dedx_get_material_list(), which reports a coarse per-program hint,
+ *  this consults the same embedded-data resolution dedx_load_config() uses (including
+ *  DEDX_AUTO's Bethe-Bloch fallback and Bragg-additivity decomposition of compounds),
+ *  so a material appears here if and only if dedx_load_config() is expected to succeed
+ *  for it with this exact ion. Thread-safe: writes only into the caller-supplied
+ *  buffer, no shared or static state.
+ *
+ *  The output is length-delimited via @p materials_len, not -1-terminated (unlike
+ *  dedx_get_material_list()); use dedx_fill_material_list_for_ion() if you need a
+ *  -1-terminated array instead.
+ *
+ *  @param[in]  program        Program identifier.
+ *  @param[in]  ion            Ion identifier.
+ *  @param[out] materials      Caller-allocated array of at least @p max_materials entries;
+ *                             DEDX_MAX_MATERIAL_LIST is a safe capacity to pass here.
+ *  @param[in]  max_materials  Capacity of @p materials; at most this many entries are written.
+ *  @param[out] materials_len  Number of material identifiers written (always <= max_materials).
+ *  @param[out] err            Error code; 0 on success. Set to DEDX_ERR_ION_NOT_SUPPORTED or
+ *                             DEDX_ERR_ESTAR_NOT_IMPL if the program/ion combination itself is
+ *                             invalid (in which case @p materials_len is set to 0). For
+ *                             DEDX_DEFAULT/DEDX_BETHE_EXT00/DEDX_AUTO, only ions 1-112 are
+ *                             accepted, since that is the actual range covered by the periodic
+ *                             table backing the Bethe-Bloch evaluation.
+ */
+void dedx_get_material_list_for_ion(
+    int program, int ion, int *materials, unsigned int max_materials, unsigned int *materials_len, int *err);
+
 /** @brief Return a null-terminated list of ions supported by a program.
  *  @param[in] program  Program identifier.
  *  @return Pointer to a static array terminated by 0; do not free.

include/dedx_wrappers.h

@@ -34,6 +34,23 @@ void dedx_fill_material_list(int program, int *material_list);
  */
 void dedx_fill_ion_list(int program, int *ion_list);
 
+/** @brief Fill an array with the materials actually usable with a program/ion combination.
+ *
+ *  Unlike dedx_fill_material_list(), which reports a coarse per-program hint,
+ *  this reflects exactly what dedx_load_config() can resolve for this specific ion
+ *  (see dedx_get_material_list_for_ion() in dedx.h for details).
+ *
+ *  @param[in]  program        Program identifier.
+ *  @param[in]  ion            Ion identifier.
+ *  @param[out] material_list  Caller-allocated array of at least DEDX_MAX_MATERIAL_LIST
+ *                             entries; filled with material identifiers terminated by -1.
+ *                             One slot is reserved for the terminator, so at most
+ *                             DEDX_MAX_MATERIAL_LIST - 1 material IDs are ever written.
+ *  @return 0 on success, non-zero error code if the program/ion combination itself is invalid
+ *          (e.g. DEDX_ERR_ION_NOT_SUPPORTED); @p material_list is still -1-terminated (empty).
+ */
+int dedx_fill_material_list_for_ion(int program, int ion, int *material_list);
+
 /** @brief Evaluate stopping power at arbitrary energy points.
  *
  *  @param[in]  program      Program identifier.

src/dedx.c

@@ -57,6 +57,9 @@ static int load_config_clean(dedx_workspace *ws, dedx_config *config, int *err);
 static int find_data(stopping_data *data, dedx_config *config, float *energy, int *err);
 static int load_compound(dedx_workspace *ws, dedx_config *config, int *err);
 static int load_bethe_2(stopping_data *data, dedx_config *config, float *energy, int *err);
+static int load_bethe_fallback(stopping_data *data, dedx_config *config, float *energy, int *err);
+static int element_supported_for_ion(int program, int ion, int element);
+static int material_id_supported(int program, int ion, int material);
 
 dedx_workspace *dedx_allocate_workspace(unsigned int count, int *err) {
     int i = 0;
@@ -228,18 +231,111 @@ const int *dedx_get_program_list(void) {
 
 const int *dedx_get_material_list(int program) {
     /* returns a list of available materials, terminated with -1 */
-    if (program == DEDX_BETHE_EXT00 || program == DEDX_DEFAULT)
+    if (program == DEDX_BETHE_EXT00 || program == DEDX_DEFAULT || program == DEDX_AUTO)
         return dedx_program_available_materials[0];
     else
         return dedx_program_available_materials[program];
 }
 
+/* Checks whether a single element (target <= 99) is reachable for program/ion,
+ * mirroring find_data()'s own dispatch (MSTAR's ion remap, DEDX_AUTO's Bethe-Bloch
+ * fallback, the resolve-across-sub-tables behavior of the ICRU-family programs) so
+ * the answer stays consistent with what dedx_load_config() will actually do.
+ *
+ * Returns 1 if dedx_load_config() is expected to succeed for this exact
+ * (program, ion, element) combination, 0 otherwise. `element` is assumed to already
+ * be a plain element id (1-99); callers are responsible for that precondition --
+ * see material_id_supported()'s `material <= 99` check below. */
+static int element_supported_for_ion(int program, int ion, int element) {
+    int ion_load = ion;
+    int resolved;
+
+    if (program == DEDX_DEFAULT || program == DEDX_BETHE_EXT00 || program == DEDX_AUTO) {
+        return 1; /* Bethe formula: any known element, given caller-supplied rho */
+    }
+    /* DEDX_ESTAR is not handled specially here: dedx_get_material_list_for_ion() -- the
+     * only caller -- already rejects it with DEDX_ERR_ESTAR_NOT_IMPL before reaching this
+     * function, and dedx_embedded_has_table() below returns 0 for it regardless (no
+     * embedded program data is registered for DEDX_ESTAR). */
+    if (program == DEDX_MSTAR && ion > 1) {
+        ion_load = 2;
+    }
+    if (program == DEDX_ICRU || program == DEDX_ICRU49 || program == DEDX_ICRU73) {
+        return dedx_embedded_resolve_program(program, ion_load, element, &resolved) == 0;
+    }
+    return dedx_embedded_has_table(program, ion_load, element);
+}
+
+/* Checks whether a material (element or compound) is reachable for program/ion.
+ * Compounds are available iff their composition is known and every constituent
+ * element is itself reachable, matching the Bragg-additivity decomposition
+ * load_compound() performs at calculation time.
+ *
+ * Returns 1 if dedx_load_config() is expected to succeed for this exact
+ * (program, ion, material) combination, 0 otherwise -- for a compound, this recurses
+ * once per constituent element and is 1 only if every one of them is reachable. */
+static int material_id_supported(int program, int ion, int material) {
+    float composition[20][2];
+    unsigned int comp_len;
+    unsigned int i;
+    int comp_err = DEDX_OK;
+
+    if (material <= 0) { /* LCOV_EXCL_START -- master material list never contains this, and
+                             composition rows are always positive element IDs in practice */
+        return 0;
+    } /* LCOV_EXCL_STOP */
+    if (material <= 99) {
+        return element_supported_for_ion(program, ion, material);
+    }
+
+    dedx_internal_get_composition(material, composition, &comp_len, &comp_err);
+    if (comp_err != DEDX_OK || comp_len == 0) { /* LCOV_EXCL_START -- every compound in the master
+                                                    material list has known embedded composition */
+        return 0;
+    } /* LCOV_EXCL_STOP */
+    for (i = 0; i < comp_len; i++) {
+        if (!material_id_supported(program, ion, (int) composition[i][0])) {
+            return 0;
+        }
+    }
+    return 1;
+}
+
+void dedx_get_material_list_for_ion(
+    int program, int ion, int *materials, unsigned int max_materials, unsigned int *materials_len, int *err) {
+    const int *candidates = dedx_program_available_materials[0];
+    unsigned int count = 0;
+    int i;
+
+    *err = DEDX_OK;
+    *materials_len = 0;
+
+    if (!check_ion(program, ion)) {
+        /* For DEDX_DEFAULT/DEDX_BETHE_EXT00/DEDX_AUTO, check_ion() itself already
+         * enforces the 1-112 range the Bethe evaluation path actually supports (see
+         * its comment), so no separate check is needed here. */
+        *err = DEDX_ERR_ION_NOT_SUPPORTED;
+        return;
+    }
+    if (program == DEDX_ESTAR) {
+        *err = DEDX_ERR_ESTAR_NOT_IMPL;
+        return;
+    }
+
+    for (i = 0; candidates[i] != -1 && count < max_materials; i++) {
+        if (material_id_supported(program, ion, candidates[i])) {
+            materials[count++] = candidates[i];
+        }
+    }
+    *materials_len = count;
+}
+
 const int *dedx_get_ion_list(int program) {
     /* Returns a -1-terminated list of ions available for the program. Both branches
      * return immutable const data -- the shared full-ion table for unrestricted
      * programs, or the program's row in dedx_program_available_ions -- so the function
      * performs no writes and is safe to call concurrently from multiple threads. */
-    if (program == DEDX_BETHE_EXT00 || program == DEDX_DEFAULT) /* any ion, no restrictions */
+    if (program == DEDX_BETHE_EXT00 || program == DEDX_DEFAULT || program == DEDX_AUTO) /* any ion, no restrictions */
         return dedx_full_ion_list;
     else
         return dedx_program_available_ions[program];
@@ -271,6 +367,13 @@ float dedx_get_min_energy(int program, int ion) {
         energy_min = get_min_energy_icru(ion);
         break;
     case DEDX_ICRU:
+    case DEDX_AUTO:
+        /* DEDX_AUTO resolves through the same tabulated sub-programs as DEDX_ICRU
+         * before falling back to Bethe-Bloch, so these advisory bounds mirror
+         * DEDX_ICRU's. These are best-effort hints, not authoritative: the actual
+         * bounds of whichever tier ends up serving a given target (a tabulated
+         * report, or the DEDX_BETHE_EXT00 fallback) may differ; dedx_get_stp()
+         * still enforces the real bounds of the dataset it loaded. */
         if (ion == DEDX_CARBON)
             energy_min = 0.001f / 12.0f;
         else
@@ -312,6 +415,7 @@ float dedx_get_max_energy(int program, int ion) {
         energy_max = get_max_energy_icru(ion);
         break;
     case DEDX_ICRU:
+    case DEDX_AUTO: /* see the comment in dedx_get_min_energy() */
         if (ion == DEDX_CARBON)
             energy_max = 10000.0f / 12.0f;
         else
@@ -515,8 +619,12 @@ static int check_ion(int prog, int ion) {
     const int *ion_list;
     int i = 0;
 
-    if (prog >= DEDX_DEFAULT) {
-        if ((ion < 1) || (ion > 120))
+    if (prog >= DEDX_DEFAULT || prog == DEDX_AUTO) {
+        /* The Bethe evaluation path (dedx_internal_get_atom_mass/charge, via
+         * dedx_periodic_table) only has data for ions 1-112, matching
+         * dedx_full_ion_list -- the same list dedx_get_ion_list() returns for these
+         * programs -- so accept nothing wider than that here. */
+        if ((ion < 1) || (ion > 112))
             return 0;
         return 1;
     }
@@ -583,6 +691,10 @@ static int find_data(stopping_data *data, dedx_config *config, float *energy, in
 
     if (prog == DEDX_ICRU) {
     } else if (prog == DEDX_ICRU49) {
+    } else if (prog == DEDX_AUTO) {
+        /* Resolve exactly like DEDX_ICRU (its tabulated sub-programs by ion type);
+         * the only difference is what happens below if that resolution fails. */
+        prog_load = DEDX_ICRU;
     } else if (prog == DEDX_ESTAR) {
         *err = DEDX_ERR_ESTAR_NOT_IMPL;
         return -1;
@@ -599,10 +711,23 @@ static int find_data(stopping_data *data, dedx_config *config, float *energy, in
         return 0;
     }
 
-    if (dedx_embedded_resolve_program(prog_load, ion_load, target_load, &prog_load) != 0
-        && (prog == DEDX_ICRU || prog == DEDX_ICRU49 || prog == DEDX_ICRU73)) {
-        *err = DEDX_ERR_COMBINATION_NOT_FOUND;
-        return -1;
+    if (dedx_embedded_resolve_program(prog_load, ion_load, target_load, &prog_load) != 0) {
+        /* DEDX_AUTO layers a Bethe-Bloch fallback on top of DEDX_ICRU's tabulated
+         * auto-selection: when no tabulated report covers this element (e.g. Boron,
+         * which none of them tabulate), compute it analytically instead of failing
+         * outright. Compounds (target > 99) are excluded here: they are handled below
+         * via Bragg-additivity decomposition into elements, each of which goes through
+         * this same fallback individually. DEDX_ICRU itself, and the report-specific
+         * programs (PSTAR, ASTAR, ICRU49, ICRU73, ...), intentionally stay literal to
+         * their published data and never get this fallback -- see DEDX_AUTO's docs
+         * in dedx.h for why. */
+        if (prog == DEDX_AUTO && target_load > 0 && target_load <= 99) {
+            return load_bethe_fallback(data, config, energy, err);
+        }
+        if (prog == DEDX_ICRU || prog == DEDX_ICRU49 || prog == DEDX_ICRU73 || prog == DEDX_AUTO) {
+            *err = DEDX_ERR_COMBINATION_NOT_FOUND;
+            return -1;
+        }
     }
     dedx_internal_read_binary_data(data, prog_load, ion_load, target_load, err);
     if (*err != 0)
@@ -706,8 +831,19 @@ static int load_bethe_2(stopping_data *data, dedx_config *config, float *energy,
         return -1;
     }
 
+    /* dedx_internal_get_atom_charge()/get_atom_mass() both reset *err to DEDX_OK and
+     * share the same id < 113 validity check, so this pair always succeeds or fails
+     * together -- one check right after both is enough. It must happen right away,
+     * though: config->target is already known <= 99 (hence < 113) at this point, so
+     * its own charge/mass lookup just below always succeeds and would otherwise
+     * silently overwrite an ion failure here, leaving *err == DEDX_OK with PZ/PA at
+     * their failure sentinel (this is reachable in practice: a compound target
+     * decomposes into elements -- see load_compound() -- before check_ion() ever runs
+     * for its ion). */
     PZ = dedx_internal_get_atom_charge(config->ion, err);
     PA = dedx_internal_get_atom_mass(config->ion, err);
+    if (*err != 0)
+        return -1;
     TZ = dedx_internal_get_atom_charge(config->target, err);
     TA = dedx_internal_get_atom_mass(config->target, err);
     rho = config->rho;
@@ -729,3 +865,42 @@ static int load_bethe_2(stopping_data *data, dedx_config *config, float *energy,
     free(bethe);
     return 0;
 }
+
+/* Last-resort tier for DEDX_AUTO: compute via the Bethe-Bloch formula using the
+ * element's embedded density and mean excitation potential. config->rho is already
+ * populated by dedx_internal_validate_rho(), which runs unconditionally before
+ * find_data(). config->_temp_i_value, however, is not: DEDX_AUTO is deliberately
+ * grouped with the tabulated report identifiers (< DEDX_DEFAULT), so
+ * dedx_internal_validate_config() skips dedx_internal_evaluate_i_pot() for it (that
+ * step only runs for program >= 100). For a compound element reached via
+ * load_compound(), _temp_i_value is already set per-element there; for a direct
+ * elemental target, it still defaults to 0 and must be fetched here.
+ *
+ * Returns 0 on success, with `data` filled in with the computed Bethe-Bloch
+ * stopping-power values and *err set to DEDX_OK, or -1 on failure, with *err set to
+ * the error code -- the same success/failure convention as find_data()/load_bethe_2(),
+ * which this function is a thin wrapper around. */
+static int load_bethe_fallback(stopping_data *data, dedx_config *config, float *energy, int *err) {
+    int prog_temp = config->program;
+    int state;
+
+    *err = DEDX_OK;
+
+    if (config->_temp_i_value == 0.0f) {
+        state = config->compound_state;
+        if (state == DEDX_DEFAULT_STATE) {
+            state = dedx_internal_target_is_gas(config->target, err) ? DEDX_GAS : DEDX_CONDENSED;
+        }
+        config->_temp_i_value = dedx_internal_get_i_value(config->target, state, err);
+        if (*err != DEDX_OK) { /* LCOV_EXCL_START -- every element in the valid 1-99
+                                   range that reaches this function has embedded I-value data */
+            return -1;
+        } /* LCOV_EXCL_STOP */
+    }
+
+    config->program = DEDX_BETHE_EXT00;
+    load_bethe_2(data, config, energy, err);
+    config->program = prog_temp;
+
+    return (*err == DEDX_OK) ? 0 : -1;
+}