diff --git a/include/dedx.h b/include/dedx.h index 14bc5df..add7237 100644 --- a/include/dedx.h +++ b/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. diff --git a/include/dedx_wrappers.h b/include/dedx_wrappers.h index edb85ba..513f319 100644 --- a/include/dedx_wrappers.h +++ b/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. diff --git a/src/dedx.c b/src/dedx.c index 75b5981..81a5d45 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_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; +} diff --git a/src/dedx_program_const.h b/src/dedx_program_const.h index f6d2227..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] = { @@ -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..b526ab1 100644 --- a/src/dedx_wrappers.c +++ b/src/dedx_wrappers.c @@ -59,6 +59,16 @@ 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; + + /* 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; +} + static dedx_config *allocate_wrapper_config(int program, int ion, int target, int *err) { dedx_config *config; 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 new file mode 100644 index 0000000..3c7e855 --- /dev/null +++ b/tests/test_material_availability.c @@ -0,0 +1,293 @@ +#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_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; + 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"); + + /* 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, + * 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; + 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); + + 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_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_auto_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_AUTO; + cfg->ion = DEDX_PROTON; + cfg->target = DEDX_BORON; + err = 0; + dedx_load_config(ws, cfg, &err); + 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 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; + 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, "ICRU+proton+BoronCarbide"); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + + /* 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; + 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 + 1]; + 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"); + + /* 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_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"); + 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"); + + /* 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"); + + 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_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++) { + 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_auto_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; +}