From 5ca4e1063abbd282ff92da01fde92f2ffbe8493c Mon Sep 17 00:00:00 2001 From: Leszek Grzanka Date: Mon, 13 Jul 2026 17:08:42 +0200 Subject: [PATCH 1/4] Fix stale material-availability tables, add ion-aware lookup and Bethe fallback for DEDX_ICRU dedx_program_available_materials falsely claimed elements (Boron, Nickel, Zirconium, Gadolinium, Tantalum as applicable) as available for ASTAR, PSTAR, ICRU49, ICRU73 and ICRU73_OLD when the embedded data for those programs never tabulated them, e.g. HYDROGEN + PSTAR + BORON (#51). The rows are corrected against the actual embedded target-id tables; MSTAR, which genuinely tabulates all of these, is left untouched. Add dedx_get_material_list_for_ion() / dedx_fill_material_list_for_ion(), which answer "what materials work with this program and this ion" by consulting the same embedded-data resolution dedx_load_config() uses, including Bragg-additivity decomposition of compounds -- so the answer can't drift out of sync with the underlying data the way the hand-maintained table did. DEDX_ICRU (the auto-selecting report) now falls back to the Bethe-Bloch formula for elements no tabulated report covers, instead of failing with DEDX_ERR_COMBINATION_NOT_FOUND. Report-specific programs (PSTAR, ASTAR, ICRU49, ICRU73, ...) intentionally keep failing those combinations rather than silently blending in another program's data. Co-Authored-By: Claude Sonnet 5 --- include/dedx.h | 28 +++- include/dedx_wrappers.h | 15 ++ src/dedx.c | 134 +++++++++++++++- src/dedx_program_const.h | 31 ++-- src/dedx_wrappers.c | 9 ++ tests/test_material_availability.c | 238 +++++++++++++++++++++++++++++ 6 files changed, 438 insertions(+), 17 deletions(-) create mode 100644 tests/test_material_availability.c diff --git a/include/dedx.h b/include/dedx.h index 14bc5df..4bdc509 100644 --- a/include/dedx.h +++ b/include/dedx.h @@ -44,7 +44,9 @@ 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. Falls back to the + Bethe-Bloch formula (see DEDX_BETHE_EXT00) for elemental targets that + no tabulated report covers (e.g. Boron), rather than failing outright. */ DEDX_DEFAULT = 100, /**< Default program (Bethe formula) */ DEDX_BETHE_EXT00 /**< Bethe formula with extensions */ }; @@ -162,6 +164,30 @@ const int *dedx_get_program_list(void); */ const int *dedx_get_material_list(int program); +/** @brief Maximum number of entries dedx_get_material_list_for_ion() can write. */ +#define DEDX_MAX_MATERIAL_LIST 300 + +/** @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_ICRU'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. + * + * @param[in] program Program identifier. + * @param[in] ion Ion identifier. + * @param[out] materials Caller-allocated array of at least DEDX_MAX_MATERIAL_LIST entries. + * @param[in] max_materials Capacity of @p materials. + * @param[out] materials_len Number of material identifiers written. + * @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). + */ +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. diff --git a/include/dedx_wrappers.h b/include/dedx_wrappers.h index edb85ba..51b1d07 100644 --- a/include/dedx_wrappers.h +++ b/include/dedx_wrappers.h @@ -34,6 +34,21 @@ 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. + * @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. diff --git a/src/dedx.c b/src/dedx.c index 75b5981..5b8d6ca 100644 --- a/src/dedx.c +++ b/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_icru_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; @@ -234,6 +237,89 @@ const int *dedx_get_material_list(int program) { 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_ICRU'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. */ +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) { + return 1; /* Bethe formula: any known element, given caller-supplied rho */ + } + if (program == DEDX_ESTAR) { + return 0; + } + if (program == DEDX_MSTAR && ion > 1) { + ion_load = 2; + } + if (program == DEDX_ICRU) { + /* Falls back to Bethe-Bloch for any element find_data() can't resolve + * via a tabulated sub-program -- see load_icru_bethe_fallback(). */ + return 1; + } + if (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. */ +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) { + return 0; + } + 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) { + return 0; + } + 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)) { + *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 @@ -599,10 +685,22 @@ 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_ICRU is documented as auto-selecting the best tabulated report for the + * ion; when no tabulated ICRU/PSTAR/ASTAR table covers this element (e.g. + * Boron, which none of them tabulate), fall back to the Bethe-Bloch formula + * 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. Report-specific programs + * (PSTAR, ASTAR, ICRU49, ICRU73, ...) intentionally stay literal to their + * published data and do not get this fallback. */ + if (prog == DEDX_ICRU && target_load > 0 && target_load <= 99) { + return load_icru_bethe_fallback(data, config, energy, err); + } + if (prog == DEDX_ICRU || prog == DEDX_ICRU49 || prog == DEDX_ICRU73) { + *err = DEDX_ERR_COMBINATION_NOT_FOUND; + return -1; + } } dedx_internal_read_binary_data(data, prog_load, ion_load, target_load, err); if (*err != 0) @@ -729,3 +827,31 @@ static int load_bethe_2(stopping_data *data, dedx_config *config, float *energy, free(bethe); return 0; } + +/* Last-resort tier for DEDX_ICRU: 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() before find_data() runs; only the + * I-value needs fetching here since tabulated (program < 100) configs skip the + * dedx_internal_evaluate_i_pot() step that normally supplies it. */ +static int load_icru_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) + return -1; + } + + config->program = DEDX_BETHE_EXT00; + load_bethe_2(data, config, energy, err); + config->program = prog_temp; + + return (*err == DEDX_OK) ? 0 : -1; +} diff --git a/src/dedx_program_const.h b/src/dedx_program_const.h index f6d2227..5b76293 100644 --- a/src/dedx_program_const.h +++ b/src/dedx_program_const.h @@ -77,6 +77,13 @@ static const int dedx_full_ion_list[113] = { 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, -1 }; +/* Per-program materials lists below are a coarse "generally usable with this program" + * hint (elements the program's own tables lack are omitted; compounds are listed + * unconditionally since Bragg additivity can serve any compound whose elements the + * program covers). They are not ion-specific and don't account for per-compound + * decomposition feasibility. For an exact, ion-specific, decomposition-aware answer, + * use dedx_get_material_list_for_ion() (see dedx.h), which consults the embedded + * target-id tables directly instead of a hand-maintained list. */ static const int dedx_program_available_materials[110][290] = { // currently excluding BETHE, which is handled in dedx.c {1,2,3,4,5,6,7,8,9, @@ -107,12 +114,12 @@ static const int dedx_program_available_materials[110][290] = { 250,251,252,253,254,255,256,257,258,259, 260,261,262,263,264,265,266,267,268,269, 270,271,272,273,274,275,276,277,278,906,-1}, // ALL materials - {1,2,4,5,6,7,8,10,13,14,18,22,26,28,29,32,36,40,42,47,50,54,64, // ASTAR - 73,74,78,79,82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, + {1,2,4,6,7,8,10,13,14,18,22,26,29,32,36,42,47,50,54,64,74,78,79, // ASTAR (no Boron/Nickel/Zirconium/Tantalum: not in dedx_astar_target_ids) + 82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, 141,155,160,169,179,185,189,191,197,200,201,202,203,204,209,213,215,216,219, 221,222,223,225,226,227,232,238,245,252,255,263,264,266,276,277,906,-1}, - {1,2,4,5,6,7,8,10,13,14,18,22,26,28,29,32,36,40,42,47,50,54,64, // PSTAR - 73,74,78,79,82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, + {1,2,4,6,7,8,10,13,14,18,22,26,29,32,36,42,47,50,54,64,74,78,79, // PSTAR (no Boron/Nickel/Zirconium/Tantalum: not in dedx_pstar_target_ids) + 82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, 141,155,160,169,179,185,189,191,197,200,201,202,203,204,209,213,215,216,219, 221,222,223,225,226,227,232,238,245,252,255,263,264,266,276,277,906,-1}, { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, @@ -134,21 +141,21 @@ static const int dedx_program_available_materials[110][290] = { 73,74,78,79,82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, 141,155,160,169,179,185,189,191,197,200,201,202,203,204,209,213,215,216,219, 221,222,223,225,226,227,232,238,245,252,255,263,264,266,276,277,906,-1}, - {1,2,4,5,6,7,8,10,13,14,18,22,26,28,29,32,36,40,42,47,50,54,64, // ICRU73_OLD - 73,74,78,79,82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, + {1,2,4,6,7,8,10,13,14,18,22,26,28,29,32,36,42,47,50,54,74,78,79, // ICRU73_OLD (no Boron/Zirconium/Gadolinium/Tantalum: not in dedx_icru73_target_ids) + 82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, 141,155,160,169,179,185,189,191,197,200,201,202,203,204,209,213,215,216,219, 221,222,223,225,226,227,232,238,245,252,255,263,264,266,276,277,906,-1}, - {1,2,4,5,6,7,8,10,13,14,18,22,26,28,29,32,36,40,42,47,50,54,64, // ICRU73 - 73,74,78,79,82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, + {1,2,4,6,7,8,10,13,14,18,22,26,28,29,32,36,42,47,50,54,74,78,79, // ICRU73 (union of dedx_icru73new_target_ids and dedx_icru73_target_ids; same exclusions as ICRU73_OLD) + 82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, 141,155,160,169,179,185,189,191,197,200,201,202,203,204,209,213,215,216,219, 221,222,223,225,226,227,232,238,245,252,255,263,264,266,276,277,906,-1}, - {1,2,4,5,6,7,8,10,13,14,18,22,26,28,29,32,36,40,42,47,50,54,64, // ICRU49 - 73,74,78,79,82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, + {1,2,4,6,7,8,10,13,14,18,22,26,29,32,36,42,47,50,54,64,74,78,79, // ICRU49 (resolves to icru_pstar/icru_astar tables; no Boron/Nickel/Zirconium/Tantalum) + 82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, 141,155,160,169,179,185,189,191,197,200,201,202,203,204,209,213,215,216,219, 221,222,223,225,226,227,232,238,245,252,255,263,264,266,276,277,906,-1}, {-1}, // reserved internal slot - {1,2,4,5,6,7,8,10,13,14,18,22,26,28,29,32,36,40,42,47,50,54,64, // ICRU - 73,74,78,79,82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, + {1,2,4,6,7,8,10,13,14,18,22,26,28,29,32,36,42,47,50,54,64,74,78,79, // ICRU (union of icru90_p/a/C, icru_pstar/astar, icru73(new+old); no Boron/Zirconium/Tantalum) + 82,92,99,101,103,104,106,111,119,120,126,130,134,138,139, 141,155,160,169,179,185,189,191,197,200,201,202,203,204,209,213,215,216,219, 221,222,223,225,226,227,232,238,245,252,255,263,264,266,276,277,906,-1} }; diff --git a/src/dedx_wrappers.c b/src/dedx_wrappers.c index c993a54..173ac86 100644 --- a/src/dedx_wrappers.c +++ b/src/dedx_wrappers.c @@ -59,6 +59,15 @@ void dedx_fill_ion_list(int program, int *ion_list) { ion_list[i] = -1; } +int dedx_fill_material_list_for_ion(int program, int ion, int *material_list) { + unsigned int len = 0; + int err = 0; + + dedx_get_material_list_for_ion(program, ion, material_list, DEDX_MAX_MATERIAL_LIST, &len, &err); + material_list[len] = -1; + return err; +} + static dedx_config *allocate_wrapper_config(int program, int ion, int target, int *err) { dedx_config *config; diff --git a/tests/test_material_availability.c b/tests/test_material_availability.c new file mode 100644 index 0000000..df81743 --- /dev/null +++ b/tests/test_material_availability.c @@ -0,0 +1,238 @@ +#include + +#include "test_helpers.h" + +/* Regression coverage for issue #51: some ion+material+program combinations were + * advertised as available (e.g. HYDROGEN + PSTAR + BORON) when the embedded data + * doesn't actually cover them, and there was no per-ion way to ask which materials + * a program really supports. These tests lock in: + * - the corrected dedx_get_material_list() tables no longer claim elements the + * embedded data doesn't have, + * - the new dedx_get_material_list_for_ion() / dedx_fill_material_list_for_ion() + * give an accurate, ion-specific, decomposition-aware answer, + * - DEDX_ICRU falls back to the Bethe-Bloch formula for elements/compounds no + * tabulated report covers, while report-specific programs (PSTAR, ASTAR, ICRU49, + * ICRU73, ...) still fail those combinations as they should. */ + +static int contains(const int *list, int value) { + int i; + for (i = 0; list[i] != -1; i++) { + if (list[i] == value) + return 1; + } + return 0; +} + +static int check_absent(const int *list, int value, const char *label) { + if (contains(list, value)) { + fprintf(stderr, "FAIL %s: value %d unexpectedly present\n", label, value); + return 1; + } + return 0; +} + +static int check_present(const int *list, int value, const char *label) { + if (!contains(list, value)) { + fprintf(stderr, "FAIL %s: value %d unexpectedly absent\n", label, value); + return 1; + } + return 0; +} + +static int check_err(int got, int expected, const char *label) { + if (got != expected) { + fprintf(stderr, "FAIL %s: got err=%d, expected %d\n", label, got, expected); + return 1; + } + return 0; +} + +/* dedx_get_material_list() rows must no longer claim elements the embedded data + * doesn't back (issue #51's literal report: PSTAR + Boron). */ +static int test_corrected_static_tables(void) { + int failures = 0; + + failures += check_absent(dedx_get_material_list(DEDX_PSTAR), DEDX_BORON, "PSTAR material list"); + failures += check_absent(dedx_get_material_list(DEDX_ASTAR), DEDX_BORON, "ASTAR material list"); + failures += check_absent(dedx_get_material_list(DEDX_ICRU49), DEDX_BORON, "ICRU49 material list"); + failures += check_absent(dedx_get_material_list(DEDX_ICRU73), DEDX_BORON, "ICRU73 material list"); + failures += check_absent(dedx_get_material_list(DEDX_ICRU73_OLD), DEDX_BORON, "ICRU73_OLD material list"); + failures += check_absent(dedx_get_material_list(DEDX_ICRU), DEDX_BORON, "ICRU material list"); + + /* MSTAR genuinely tabulates Boron -- must not be over-corrected. */ + failures += check_present(dedx_get_material_list(DEDX_MSTAR), DEDX_BORON, "MSTAR material list"); + /* Nickel is genuinely tabulated by ICRU73/ICRU73_OLD -- must not be over-corrected. */ + failures += check_present(dedx_get_material_list(DEDX_ICRU73), DEDX_NICKEL, "ICRU73 material list"); + failures += check_present(dedx_get_material_list(DEDX_ICRU73_OLD), DEDX_NICKEL, "ICRU73_OLD material list"); + + return failures; +} + +/* dedx_load_config() must keep rejecting Boron for report-specific programs... */ +static int test_report_specific_programs_reject_boron(void) { + int failures = 0; + int err; + dedx_workspace *ws; + dedx_config *cfg; + + ws = dedx_allocate_workspace(1, &err); + cfg = calloc(1, sizeof(dedx_config)); + cfg->program = DEDX_PSTAR; + cfg->ion = DEDX_PROTON; + cfg->target = DEDX_BORON; + err = 0; + dedx_load_config(ws, cfg, &err); + failures += check_err(err, DEDX_ERR_COMBINATION_NOT_FOUND, "PSTAR+proton+Boron"); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + + ws = dedx_allocate_workspace(1, &err); + cfg = calloc(1, sizeof(dedx_config)); + cfg->program = DEDX_ICRU73; + cfg->ion = DEDX_LITHIUM; + cfg->target = DEDX_BORON; + err = 0; + dedx_load_config(ws, cfg, &err); + failures += check_err(err, DEDX_ERR_COMBINATION_NOT_FOUND, "ICRU73+lithium+Boron"); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + + return failures; +} + +/* ...while DEDX_ICRU falls back to the Bethe-Bloch formula for the same element, + * both directly and via Bragg-additivity decomposition of a Boron compound. */ +static int test_icru_bethe_fallback(void) { + int failures = 0; + int err; + dedx_workspace *ws; + dedx_config *cfg; + float stp; + + ws = dedx_allocate_workspace(1, &err); + cfg = calloc(1, sizeof(dedx_config)); + cfg->program = DEDX_ICRU; + cfg->ion = DEDX_PROTON; + cfg->target = DEDX_BORON; + err = 0; + dedx_load_config(ws, cfg, &err); + failures += check_err(err, DEDX_OK, "ICRU+proton+Boron load"); + if (err == DEDX_OK) { + stp = dedx_get_stp(ws, cfg, 100.0f, &err); + if (err != DEDX_OK || stp <= 0.0f) { + fprintf(stderr, "FAIL ICRU+proton+Boron stp: err=%d stp=%f\n", err, (double) stp); + failures++; + } + } + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + + ws = dedx_allocate_workspace(1, &err); + cfg = calloc(1, sizeof(dedx_config)); + cfg->program = DEDX_ICRU; + cfg->ion = DEDX_PROTON; + cfg->target = DEDX_BORON_CARBIDE; + err = 0; + dedx_load_config(ws, cfg, &err); + failures += check_err(err, DEDX_OK, "ICRU+proton+BoronCarbide load"); + failures += check_err(cfg->bragg_used, 1, "ICRU+proton+BoronCarbide should decompose via Bragg additivity"); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + + /* PSTAR still can't serve the same compound: Boron is not in its table and + * PSTAR does not get the Bethe fallback. */ + ws = dedx_allocate_workspace(1, &err); + cfg = calloc(1, sizeof(dedx_config)); + cfg->program = DEDX_PSTAR; + cfg->ion = DEDX_PROTON; + cfg->target = DEDX_BORON_CARBIDE; + err = 0; + dedx_load_config(ws, cfg, &err); + failures += check_err(err, DEDX_ERR_COMBINATION_NOT_FOUND, "PSTAR+proton+BoronCarbide"); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + + return failures; +} + +static int test_material_list_for_ion(void) { + int failures = 0; + int materials[DEDX_MAX_MATERIAL_LIST]; + unsigned int len; + int err; + + dedx_get_material_list_for_ion(DEDX_PSTAR, DEDX_PROTON, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); + failures += check_err(err, DEDX_OK, "PSTAR+proton for_ion err"); + materials[len] = -1; + failures += check_absent(materials, DEDX_BORON, "PSTAR+proton for_ion"); + failures += check_present(materials, DEDX_WATER, "PSTAR+proton for_ion"); + + /* ASTAR only supports helium; requesting it with a proton is an ion mismatch, + * not a material-availability question. */ + dedx_get_material_list_for_ion(DEDX_ASTAR, DEDX_PROTON, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); + failures += check_err(err, DEDX_ERR_ION_NOT_SUPPORTED, "ASTAR+proton for_ion"); + failures += check_err((int) len, 0, "ASTAR+proton for_ion count"); + + dedx_get_material_list_for_ion(DEDX_ICRU, DEDX_PROTON, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); + failures += check_err(err, DEDX_OK, "ICRU+proton for_ion err"); + materials[len] = -1; + failures += check_present(materials, DEDX_BORON, "ICRU+proton for_ion"); + failures += check_present(materials, DEDX_BORON_CARBIDE, "ICRU+proton for_ion"); + + dedx_get_material_list_for_ion(DEDX_ICRU73, DEDX_LITHIUM, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); + failures += check_err(err, DEDX_OK, "ICRU73+lithium for_ion err"); + materials[len] = -1; + failures += check_absent(materials, DEDX_BORON, "ICRU73+lithium for_ion"); + failures += check_absent(materials, DEDX_BORON_CARBIDE, "ICRU73+lithium for_ion"); + + dedx_get_material_list_for_ion(DEDX_ESTAR, 1001, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); + failures += check_err(err, DEDX_ERR_ESTAR_NOT_IMPL, "ESTAR+electron for_ion"); + failures += check_err((int) len, 0, "ESTAR+electron for_ion count"); + + return failures; +} + +/* The wrapper must agree with the core API and remain -1-terminated. */ +static int test_fill_material_list_for_ion(void) { + int failures = 0; + int wrapped[DEDX_MAX_MATERIAL_LIST + 1]; + int direct[DEDX_MAX_MATERIAL_LIST]; + unsigned int direct_len; + int err; + int i; + + int wrap_err = dedx_fill_material_list_for_ion(DEDX_ICRU, DEDX_PROTON, wrapped); + dedx_get_material_list_for_ion(DEDX_ICRU, DEDX_PROTON, direct, DEDX_MAX_MATERIAL_LIST, &direct_len, &err); + + failures += check_err(wrap_err, DEDX_OK, "fill_material_list_for_ion err"); + for (i = 0; i < (int) direct_len; i++) { + if (wrapped[i] != direct[i]) { + fprintf(stderr, "FAIL fill_material_list_for_ion: entry %d differs (%d vs %d)\n", i, wrapped[i], direct[i]); + failures++; + break; + } + } + failures += check_err(wrapped[direct_len], -1, "fill_material_list_for_ion terminator"); + + /* Invalid ion/program combination: wrapper reports the error and still + * produces an empty, -1-terminated list. */ + wrap_err = dedx_fill_material_list_for_ion(DEDX_ASTAR, DEDX_PROTON, wrapped); + failures += check_err(wrap_err, DEDX_ERR_ION_NOT_SUPPORTED, "fill_material_list_for_ion invalid combo err"); + failures += check_err(wrapped[0], -1, "fill_material_list_for_ion invalid combo terminator"); + + return failures; +} + +int main(void) { + int failures = 0; + + failures += test_corrected_static_tables(); + failures += test_report_specific_programs_reject_boron(); + failures += test_icru_bethe_fallback(); + failures += test_material_list_for_ion(); + failures += test_fill_material_list_for_ion(); + + if (failures == 0) + printf("test_material_availability: all checks passed\n"); + return failures; +} From 5749f9f246bd92b1f41e6a20b2deb7ef66d99ff1 Mon Sep 17 00:00:00 2001 From: Leszek Grzanka Date: Mon, 13 Jul 2026 17:43:10 +0200 Subject: [PATCH 2/4] Split Bethe-Bloch fallback into DEDX_AUTO, address review feedback DEDX_ICRU and DEDX_AUTO are different models and shouldn't be conflated: DEDX_ICRU stays a real, literal, tabulated auto-select between ICRU49/ICRU73 (no fallback), while the new DEDX_AUTO program layers a Bethe-Bloch tier on top for elements no tabulated report covers -- comparable to the AUTO mode in dedx_web. DEDX_AUTO is deliberately grouped with the other tabulated report identifiers (not >= DEDX_DEFAULT) so a compound target still resolves through the same tabulated-first path as DEDX_ICRU instead of always decomposing into elements, which was the case in an earlier draft of this change and gave DEDX_AUTO wrong stopping powers for directly-tabulated compounds like water. Also addresses Copilot's review comments on PR #144: reserves the terminator slot in dedx_fill_material_list_for_ion() and its test, clarifies the DEDX_MAX_MATERIAL_LIST/output-format docs, and rejects ions the periodic table doesn't actually cover (>112) for the Bethe-capable programs in dedx_get_material_list_for_ion(), consistent with what dedx_load_config() does. Bumps DEDX_MAX_MATERIAL_LIST to 1000 per request. Closes the codecov patch-coverage gaps: two of the added guards are truly unreachable given today's embedded data (verified by probing every candidate element/compound) and are marked LCOV_EXCL with a note why; the rest needed a DEDX_MSTAR ion-remap test case and fell out of removing now-dead i-value-fetch logic that direct testing showed was never reachable in the first draft. Co-Authored-By: Claude Sonnet 5 --- include/dedx.h | 48 ++++++++++++--- include/dedx_wrappers.h | 2 + src/dedx.c | 97 +++++++++++++++++++----------- src/dedx_program_const.h | 6 +- src/dedx_wrappers.c | 3 +- tests/test_material_availability.c | 94 +++++++++++++++++++++++------ 6 files changed, 181 insertions(+), 69 deletions(-) diff --git a/include/dedx.h b/include/dedx.h index 4bdc509..add7237 100644 --- a/include/dedx.h +++ b/include/dedx.h @@ -44,9 +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. Falls back to the - Bethe-Bloch formula (see DEDX_BETHE_EXT00) for elemental targets that - no tabulated report covers (e.g. Boron), rather than failing outright. */ + 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 */ }; @@ -164,26 +178,40 @@ const int *dedx_get_program_list(void); */ const int *dedx_get_material_list(int program); -/** @brief Maximum number of entries dedx_get_material_list_for_ion() can write. */ -#define DEDX_MAX_MATERIAL_LIST 300 +/** @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_ICRU's Bethe-Bloch fallback and Bragg-additivity decomposition of compounds), + * 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 DEDX_MAX_MATERIAL_LIST entries. - * @param[in] max_materials Capacity of @p materials. - * @param[out] materials_len Number of material identifiers written. + * @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). + * 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); diff --git a/include/dedx_wrappers.h b/include/dedx_wrappers.h index 51b1d07..513f319 100644 --- a/include/dedx_wrappers.h +++ b/include/dedx_wrappers.h @@ -44,6 +44,8 @@ void dedx_fill_ion_list(int program, int *ion_list); * @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). */ diff --git a/src/dedx.c b/src/dedx.c index 5b8d6ca..4604aaa 100644 --- a/src/dedx.c +++ b/src/dedx.c @@ -57,7 +57,7 @@ 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_icru_bethe_fallback(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); @@ -231,35 +231,31 @@ 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_ICRU's Bethe-Bloch + * 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. */ 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) { + if (program == DEDX_DEFAULT || program == DEDX_BETHE_EXT00 || program == DEDX_AUTO) { return 1; /* Bethe formula: any known element, given caller-supplied rho */ } - if (program == DEDX_ESTAR) { - return 0; - } + /* 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) { - /* Falls back to Bethe-Bloch for any element find_data() can't resolve - * via a tabulated sub-program -- see load_icru_bethe_fallback(). */ - return 1; - } - if (program == DEDX_ICRU49 || program == DEDX_ICRU73) { + 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); @@ -275,17 +271,19 @@ static int material_id_supported(int program, int ion, int material) { unsigned int i; int comp_err = DEDX_OK; - if (material <= 0) { + 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) { + 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; @@ -311,6 +309,14 @@ void dedx_get_material_list_for_ion( *err = DEDX_ERR_ESTAR_NOT_IMPL; return; } + if ((program == DEDX_DEFAULT || program == DEDX_BETHE_EXT00 || program == DEDX_AUTO) && (ion < 1 || ion > 112)) { + /* check_ion() accepts any ion in [1, 120] for these programs, but the Bethe + * evaluation path (dedx_internal_get_atom_mass/charge, via dedx_periodic_table) + * only has data for ions 1-112. Keep this function consistent with what + * dedx_load_config() will actually do for the same ion. */ + *err = DEDX_ERR_ION_NOT_SUPPORTED; + return; + } for (i = 0; candidates[i] != -1 && count < max_materials; i++) { if (material_id_supported(program, ion, candidates[i])) { @@ -325,7 +331,7 @@ const int *dedx_get_ion_list(int program) { * 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]; @@ -357,6 +363,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 @@ -398,6 +411,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 @@ -601,7 +615,7 @@ static int check_ion(int prog, int ion) { const int *ion_list; int i = 0; - if (prog >= DEDX_DEFAULT) { + if (prog >= DEDX_DEFAULT || prog == DEDX_AUTO) { if ((ion < 1) || (ion > 120)) return 0; return 1; @@ -669,6 +683,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; @@ -686,18 +704,19 @@ static int find_data(stopping_data *data, dedx_config *config, float *energy, in } if (dedx_embedded_resolve_program(prog_load, ion_load, target_load, &prog_load) != 0) { - /* DEDX_ICRU is documented as auto-selecting the best tabulated report for the - * ion; when no tabulated ICRU/PSTAR/ASTAR table covers this element (e.g. - * Boron, which none of them tabulate), fall back to the Bethe-Bloch formula - * 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. Report-specific programs - * (PSTAR, ASTAR, ICRU49, ICRU73, ...) intentionally stay literal to their - * published data and do not get this fallback. */ - if (prog == DEDX_ICRU && target_load > 0 && target_load <= 99) { - return load_icru_bethe_fallback(data, config, energy, err); + /* 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) { + if (prog == DEDX_ICRU || prog == DEDX_ICRU49 || prog == DEDX_ICRU73 || prog == DEDX_AUTO) { *err = DEDX_ERR_COMBINATION_NOT_FOUND; return -1; } @@ -828,12 +847,16 @@ static int load_bethe_2(stopping_data *data, dedx_config *config, float *energy, return 0; } -/* Last-resort tier for DEDX_ICRU: compute via the Bethe-Bloch formula using the +/* 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() before find_data() runs; only the - * I-value needs fetching here since tabulated (program < 100) configs skip the - * dedx_internal_evaluate_i_pot() step that normally supplies it. */ -static int load_icru_bethe_fallback(stopping_data *data, dedx_config *config, float *energy, int *err) { + * 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. */ +static int load_bethe_fallback(stopping_data *data, dedx_config *config, float *energy, int *err) { int prog_temp = config->program; int state; @@ -845,8 +868,10 @@ static int load_icru_bethe_fallback(stopping_data *data, dedx_config *config, fl 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) + 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; diff --git a/src/dedx_program_const.h b/src/dedx_program_const.h index 5b76293..ff4171e 100644 --- a/src/dedx_program_const.h +++ b/src/dedx_program_const.h @@ -6,7 +6,7 @@ static const char dedx_program_table[110][20] = { // clang-format off "(N/A)","ASTAR","PSTAR","ESTAR","MSTAR","ICRU73_OLD","ICRU73","ICRU49","","ICRU", // 0 - 9 - "","","","","","","","","","", // 10-19 + "AUTO","","","","","","","","","", // 10-19 "","","","","","","","","","", // 20-29 "","","","","","","","","","", // 30-39 "","","","","","","","","","", // 40-49 @@ -19,7 +19,7 @@ static const char dedx_program_table[110][20] = { }; static const char dedx_program_version_table[110][20] = { // clang-format off "(N/A)","1.2","1.2","ESTAR","3.12","2005/06","2005/06-REV","ICRU Report 49 D","","1.0", // 0 - 9 - "","","","","","","","","","", // 10-19 + "1.0","","","","","","","","","", // 10-19 "","","","","","","","","","", // 20-29 "","","","","","","","","","", // 30-39 "","","","","","","","","","", // 40-49 @@ -35,7 +35,7 @@ static const int dedx_available_programs[20] = { DEDX_ASTAR, DEDX_PSTAR, DEDX_ESTAR, DEDX_MSTAR, DEDX_ICRU73_OLD, DEDX_ICRU73, - DEDX_ICRU49, DEDX_ICRU, DEDX_DEFAULT, + DEDX_ICRU49, DEDX_ICRU, DEDX_AUTO, DEDX_DEFAULT, DEDX_BETHE_EXT00, -1}; static const int dedx_program_available_ions[110][20] = { diff --git a/src/dedx_wrappers.c b/src/dedx_wrappers.c index 173ac86..b526ab1 100644 --- a/src/dedx_wrappers.c +++ b/src/dedx_wrappers.c @@ -63,7 +63,8 @@ int dedx_fill_material_list_for_ion(int program, int ion, int *material_list) { unsigned int len = 0; int err = 0; - dedx_get_material_list_for_ion(program, ion, material_list, DEDX_MAX_MATERIAL_LIST, &len, &err); + /* Reserve the last slot for the -1 terminator. */ + dedx_get_material_list_for_ion(program, ion, material_list, DEDX_MAX_MATERIAL_LIST - 1, &len, &err); material_list[len] = -1; return err; } diff --git a/tests/test_material_availability.c b/tests/test_material_availability.c index df81743..31d6b79 100644 --- a/tests/test_material_availability.c +++ b/tests/test_material_availability.c @@ -10,9 +10,10 @@ * embedded data doesn't have, * - the new dedx_get_material_list_for_ion() / dedx_fill_material_list_for_ion() * give an accurate, ion-specific, decomposition-aware answer, - * - DEDX_ICRU falls back to the Bethe-Bloch formula for elements/compounds no - * tabulated report covers, while report-specific programs (PSTAR, ASTAR, ICRU49, - * ICRU73, ...) still fail those combinations as they should. */ + * - DEDX_AUTO falls back to the Bethe-Bloch formula for elements/compounds no + * tabulated report covers, while DEDX_ICRU and the other report-specific + * programs (PSTAR, ASTAR, ICRU49, ICRU73, ...) still fail those combinations, + * as they should -- DEDX_ICRU and DEDX_AUTO are deliberately different models. */ static int contains(const int *list, int value) { int i; @@ -65,10 +66,16 @@ static int test_corrected_static_tables(void) { failures += check_present(dedx_get_material_list(DEDX_ICRU73), DEDX_NICKEL, "ICRU73 material list"); failures += check_present(dedx_get_material_list(DEDX_ICRU73_OLD), DEDX_NICKEL, "ICRU73_OLD material list"); + /* DEDX_AUTO can reach any element via its Bethe-Bloch tier, same as DEDX_DEFAULT + * and DEDX_BETHE_EXT00. */ + failures += check_present(dedx_get_material_list(DEDX_AUTO), DEDX_BORON, "AUTO material list"); + return failures; } -/* dedx_load_config() must keep rejecting Boron for report-specific programs... */ +/* dedx_load_config() must keep rejecting Boron for report-specific programs, + * including DEDX_ICRU itself: it is a real, tabulated auto-selecting report, not a + * Bethe-Bloch hybrid, so a combination no ICRU sub-report covers still fails. */ static int test_report_specific_programs_reject_boron(void) { int failures = 0; int err; @@ -97,12 +104,23 @@ static int test_report_specific_programs_reject_boron(void) { dedx_free_config(cfg, &err); dedx_free_workspace(ws, &err); + ws = dedx_allocate_workspace(1, &err); + cfg = calloc(1, sizeof(dedx_config)); + cfg->program = DEDX_ICRU; + cfg->ion = DEDX_PROTON; + cfg->target = DEDX_BORON; + err = 0; + dedx_load_config(ws, cfg, &err); + failures += check_err(err, DEDX_ERR_COMBINATION_NOT_FOUND, "ICRU+proton+Boron (ICRU must not fall back)"); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + return failures; } -/* ...while DEDX_ICRU falls back to the Bethe-Bloch formula for the same element, +/* ...while DEDX_AUTO falls back to the Bethe-Bloch formula for the same element, * both directly and via Bragg-additivity decomposition of a Boron compound. */ -static int test_icru_bethe_fallback(void) { +static int test_auto_bethe_fallback(void) { int failures = 0; int err; dedx_workspace *ws; @@ -111,22 +129,36 @@ static int test_icru_bethe_fallback(void) { ws = dedx_allocate_workspace(1, &err); cfg = calloc(1, sizeof(dedx_config)); - cfg->program = DEDX_ICRU; + cfg->program = DEDX_AUTO; cfg->ion = DEDX_PROTON; cfg->target = DEDX_BORON; err = 0; dedx_load_config(ws, cfg, &err); - failures += check_err(err, DEDX_OK, "ICRU+proton+Boron load"); + failures += check_err(err, DEDX_OK, "AUTO+proton+Boron load"); if (err == DEDX_OK) { stp = dedx_get_stp(ws, cfg, 100.0f, &err); if (err != DEDX_OK || stp <= 0.0f) { - fprintf(stderr, "FAIL ICRU+proton+Boron stp: err=%d stp=%f\n", err, (double) stp); + fprintf(stderr, "FAIL AUTO+proton+Boron stp: err=%d stp=%f\n", err, (double) stp); failures++; } } dedx_free_config(cfg, &err); dedx_free_workspace(ws, &err); + ws = dedx_allocate_workspace(1, &err); + cfg = calloc(1, sizeof(dedx_config)); + cfg->program = DEDX_AUTO; + cfg->ion = DEDX_PROTON; + cfg->target = DEDX_BORON_CARBIDE; + err = 0; + dedx_load_config(ws, cfg, &err); + failures += check_err(err, DEDX_OK, "AUTO+proton+BoronCarbide load"); + failures += check_err(cfg->bragg_used, 1, "AUTO+proton+BoronCarbide should decompose via Bragg additivity"); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + + /* DEDX_ICRU still can't serve the same compound: it does not get the Bethe + * fallback, unlike DEDX_AUTO. */ ws = dedx_allocate_workspace(1, &err); cfg = calloc(1, sizeof(dedx_config)); cfg->program = DEDX_ICRU; @@ -134,13 +166,12 @@ static int test_icru_bethe_fallback(void) { cfg->target = DEDX_BORON_CARBIDE; err = 0; dedx_load_config(ws, cfg, &err); - failures += check_err(err, DEDX_OK, "ICRU+proton+BoronCarbide load"); - failures += check_err(cfg->bragg_used, 1, "ICRU+proton+BoronCarbide should decompose via Bragg additivity"); + failures += check_err(err, DEDX_ERR_COMBINATION_NOT_FOUND, "ICRU+proton+BoronCarbide"); dedx_free_config(cfg, &err); dedx_free_workspace(ws, &err); - /* PSTAR still can't serve the same compound: Boron is not in its table and - * PSTAR does not get the Bethe fallback. */ + /* PSTAR still can't serve the same compound either: Boron is not in its table + * and PSTAR does not get the Bethe fallback. */ ws = dedx_allocate_workspace(1, &err); cfg = calloc(1, sizeof(dedx_config)); cfg->program = DEDX_PSTAR; @@ -157,7 +188,7 @@ static int test_icru_bethe_fallback(void) { static int test_material_list_for_ion(void) { int failures = 0; - int materials[DEDX_MAX_MATERIAL_LIST]; + int materials[DEDX_MAX_MATERIAL_LIST + 1]; unsigned int len; int err; @@ -167,17 +198,34 @@ static int test_material_list_for_ion(void) { failures += check_absent(materials, DEDX_BORON, "PSTAR+proton for_ion"); failures += check_present(materials, DEDX_WATER, "PSTAR+proton for_ion"); + /* MSTAR remaps every ion above helium onto its single reference table; helium + * itself (ion == 2) already takes that remap path, since MSTAR has no ion == 1 + * entry. MSTAR genuinely tabulates Boron (unlike PSTAR/ASTAR/ICRU49/ICRU73). */ + dedx_get_material_list_for_ion(DEDX_MSTAR, DEDX_HELIUM, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); + failures += check_err(err, DEDX_OK, "MSTAR+helium for_ion err"); + materials[len] = -1; + failures += check_present(materials, DEDX_BORON, "MSTAR+helium for_ion"); + failures += check_present(materials, DEDX_WATER, "MSTAR+helium for_ion"); + /* ASTAR only supports helium; requesting it with a proton is an ion mismatch, * not a material-availability question. */ dedx_get_material_list_for_ion(DEDX_ASTAR, DEDX_PROTON, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); failures += check_err(err, DEDX_ERR_ION_NOT_SUPPORTED, "ASTAR+proton for_ion"); failures += check_err((int) len, 0, "ASTAR+proton for_ion count"); + /* DEDX_ICRU never falls back -- Boron and Boron Carbide stay absent. */ dedx_get_material_list_for_ion(DEDX_ICRU, DEDX_PROTON, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); failures += check_err(err, DEDX_OK, "ICRU+proton for_ion err"); materials[len] = -1; - failures += check_present(materials, DEDX_BORON, "ICRU+proton for_ion"); - failures += check_present(materials, DEDX_BORON_CARBIDE, "ICRU+proton for_ion"); + failures += check_absent(materials, DEDX_BORON, "ICRU+proton for_ion"); + failures += check_absent(materials, DEDX_BORON_CARBIDE, "ICRU+proton for_ion"); + + /* DEDX_AUTO does fall back -- both appear. */ + dedx_get_material_list_for_ion(DEDX_AUTO, DEDX_PROTON, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); + failures += check_err(err, DEDX_OK, "AUTO+proton for_ion err"); + materials[len] = -1; + failures += check_present(materials, DEDX_BORON, "AUTO+proton for_ion"); + failures += check_present(materials, DEDX_BORON_CARBIDE, "AUTO+proton for_ion"); dedx_get_material_list_for_ion(DEDX_ICRU73, DEDX_LITHIUM, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); failures += check_err(err, DEDX_OK, "ICRU73+lithium for_ion err"); @@ -189,6 +237,14 @@ static int test_material_list_for_ion(void) { failures += check_err(err, DEDX_ERR_ESTAR_NOT_IMPL, "ESTAR+electron for_ion"); failures += check_err((int) len, 0, "ESTAR+electron for_ion count"); + /* check_ion() alone accepts any ion in [1, 120] for DEDX_AUTO/DEFAULT/BETHE_EXT00, + * but the periodic-table data backing the Bethe evaluation only covers ions up to + * 112; an ion in (112, 120] must still be rejected, consistently with what + * dedx_load_config() would do. */ + dedx_get_material_list_for_ion(DEDX_AUTO, 115, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); + failures += check_err(err, DEDX_ERR_ION_NOT_SUPPORTED, "AUTO+ion115 for_ion"); + failures += check_err((int) len, 0, "AUTO+ion115 for_ion count"); + return failures; } @@ -201,8 +257,8 @@ static int test_fill_material_list_for_ion(void) { int err; int i; - int wrap_err = dedx_fill_material_list_for_ion(DEDX_ICRU, DEDX_PROTON, wrapped); - dedx_get_material_list_for_ion(DEDX_ICRU, DEDX_PROTON, direct, DEDX_MAX_MATERIAL_LIST, &direct_len, &err); + int wrap_err = dedx_fill_material_list_for_ion(DEDX_AUTO, DEDX_PROTON, wrapped); + dedx_get_material_list_for_ion(DEDX_AUTO, DEDX_PROTON, direct, DEDX_MAX_MATERIAL_LIST, &direct_len, &err); failures += check_err(wrap_err, DEDX_OK, "fill_material_list_for_ion err"); for (i = 0; i < (int) direct_len; i++) { @@ -228,7 +284,7 @@ int main(void) { failures += test_corrected_static_tables(); failures += test_report_specific_programs_reject_boron(); - failures += test_icru_bethe_fallback(); + failures += test_auto_bethe_fallback(); failures += test_material_list_for_ion(); failures += test_fill_material_list_for_ion(); From fea3ff8c0bb1b467c8536571f2be8bcd5cc01a13 Mon Sep 17 00:00:00 2001 From: Leszek Grzanka Date: Mon, 13 Jul 2026 17:59:44 +0200 Subject: [PATCH 3/4] Document return-value contracts for material-availability helpers Clarify what element_supported_for_ion(), material_id_supported(), and load_bethe_fallback() return, for readers who only skim the function signatures. Comment-only change, no behavior difference. Co-Authored-By: Claude Sonnet 5 --- src/dedx.c | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/src/dedx.c b/src/dedx.c index 4604aaa..d483266 100644 --- a/src/dedx.c +++ b/src/dedx.c @@ -240,7 +240,12 @@ const int *dedx_get_material_list(int 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. */ + * 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; @@ -264,7 +269,11 @@ static int element_supported_for_ion(int program, int ion, int 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. */ + * 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; @@ -855,7 +864,12 @@ static int load_bethe_2(stopping_data *data, dedx_config *config, float *energy, * 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. */ + * 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; From 3bbfe2cae43d1ac9680b320da8349d0c5c963947 Mon Sep 17 00:00:00 2001 From: Leszek Grzanka Date: Mon, 13 Jul 2026 18:37:28 +0200 Subject: [PATCH 4/4] Constrain ion range to 1-112 for Bethe-driven programs, fix silent NaN Copilot pointed out that check_ion() accepted ions up to 120 for DEDX_DEFAULT/DEDX_BETHE_EXT00/DEDX_AUTO, wider than the 1-112 range the periodic-table data backing the Bethe-Bloch formula actually supports. Tightened check_ion() itself to 1-112 (matching dedx_get_ion_list()) and dropped the now-redundant duplicate check in dedx_get_material_list_for_ion(). Verifying that surfaced a worse, pre-existing bug: for DEDX_DEFAULT/ DEDX_BETHE_EXT00 with a compound target, the compound decomposes into elements before check_ion() ever runs, so an invalid ion had to be caught inside load_bethe_2() instead -- and it wasn't. dedx_internal_get_atom_charge()/ get_atom_mass() each reset *err to DEDX_OK on entry, so the ion's DEDX_ERR_NOT_AN_ELEMENT was silently overwritten by the (valid) target element's successful lookup right after it, yielding a "successful" load with a NaN stopping power. Fixed by checking *err right after the ion's charge/mass pair, before the target's pair can clobber it. Co-Authored-By: Claude Sonnet 5 --- src/dedx.c | 28 ++++++++++++++-------- tests/test_error_codes.c | 37 ++++++++++++++++++++++++++++++ tests/test_material_availability.c | 7 +++--- 3 files changed, 59 insertions(+), 13 deletions(-) diff --git a/src/dedx.c b/src/dedx.c index d483266..81a5d45 100644 --- a/src/dedx.c +++ b/src/dedx.c @@ -311,6 +311,9 @@ void dedx_get_material_list_for_ion( *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; } @@ -318,14 +321,6 @@ void dedx_get_material_list_for_ion( *err = DEDX_ERR_ESTAR_NOT_IMPL; return; } - if ((program == DEDX_DEFAULT || program == DEDX_BETHE_EXT00 || program == DEDX_AUTO) && (ion < 1 || ion > 112)) { - /* check_ion() accepts any ion in [1, 120] for these programs, but the Bethe - * evaluation path (dedx_internal_get_atom_mass/charge, via dedx_periodic_table) - * only has data for ions 1-112. Keep this function consistent with what - * dedx_load_config() will actually do for the same ion. */ - *err = DEDX_ERR_ION_NOT_SUPPORTED; - return; - } for (i = 0; candidates[i] != -1 && count < max_materials; i++) { if (material_id_supported(program, ion, candidates[i])) { @@ -625,7 +620,11 @@ static int check_ion(int prog, int ion) { int i = 0; if (prog >= DEDX_DEFAULT || prog == DEDX_AUTO) { - if ((ion < 1) || (ion > 120)) + /* 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; } @@ -832,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; diff --git a/tests/test_error_codes.c b/tests/test_error_codes.c index a3b9d1c..065ca7a 100644 --- a/tests/test_error_codes.c +++ b/tests/test_error_codes.c @@ -130,5 +130,42 @@ int main(void) { dedx_free_config(cfg, &err); dedx_free_workspace(ws, &err); + /* DEDX_DEFAULT/DEDX_BETHE_EXT00 only support ions the periodic table covers + * (1-112); an out-of-range ion must fail cleanly, not silently. */ + ws = dedx_allocate_workspace(1, &err); + cfg = calloc(1, sizeof(dedx_config)); + cfg->program = DEDX_DEFAULT; + cfg->ion = 115; /* beyond dedx_periodic_table's 1-112 range */ + cfg->target = DEDX_WATER; + err = 0; + dedx_load_config(ws, cfg, &err); + /* Regression test: DEDX_DEFAULT/DEDX_BETHE_EXT00 with a compound target + * decomposes into elements (dedx_internal_evaluate_compound(), called eagerly for + * program >= 100) before check_ion() would ever run, so an invalid ion has to be + * caught inside the Bethe evaluation itself. It previously wasn't: load_bethe_2() + * called dedx_internal_get_atom_charge()/get_atom_mass() for the ion and then + * again for the (valid) target element, and the target's success silently + * overwrote the ion's DEDX_ERR_NOT_AN_ELEMENT, yielding a loaded config with a NaN + * stopping power instead of a load failure. */ + failures += check_err(err, DEDX_ERR_NOT_AN_ELEMENT, "DEFAULT with invalid ion and compound target"); + failures += check_err(cfg->loaded, 0, "failed load should not mark config loaded"); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + + /* Same invalid ion, but with a plain elemental target: this path goes through + * load_config_clean()'s check_ion() instead, which independently must also reject + * ions the periodic table doesn't cover for DEDX_DEFAULT/DEDX_BETHE_EXT00/DEDX_AUTO. */ + ws = dedx_allocate_workspace(1, &err); + cfg = calloc(1, sizeof(dedx_config)); + cfg->program = DEDX_DEFAULT; + cfg->ion = 115; + cfg->target = DEDX_CARBON; + cfg->rho = 1.0f; + err = 0; + dedx_load_config(ws, cfg, &err); + failures += check_err(err, DEDX_ERR_ION_NOT_SUPPORTED, "DEFAULT with invalid ion and elemental target"); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + return failures; } diff --git a/tests/test_material_availability.c b/tests/test_material_availability.c index 31d6b79..3c7e855 100644 --- a/tests/test_material_availability.c +++ b/tests/test_material_availability.c @@ -237,10 +237,9 @@ static int test_material_list_for_ion(void) { failures += check_err(err, DEDX_ERR_ESTAR_NOT_IMPL, "ESTAR+electron for_ion"); failures += check_err((int) len, 0, "ESTAR+electron for_ion count"); - /* check_ion() alone accepts any ion in [1, 120] for DEDX_AUTO/DEFAULT/BETHE_EXT00, - * but the periodic-table data backing the Bethe evaluation only covers ions up to - * 112; an ion in (112, 120] must still be rejected, consistently with what - * dedx_load_config() would do. */ + /* The periodic-table data backing the Bethe evaluation only covers ions up to 112 + * for DEDX_AUTO/DEFAULT/BETHE_EXT00 (see check_ion()); an ion beyond that must be + * rejected, consistently with what dedx_load_config() would do. */ dedx_get_material_list_for_ion(DEDX_AUTO, 115, materials, DEDX_MAX_MATERIAL_LIST, &len, &err); failures += check_err(err, DEDX_ERR_ION_NOT_SUPPORTED, "AUTO+ion115 for_ion"); failures += check_err((int) len, 0, "AUTO+ion115 for_ion count");