diff --git a/include/dedx_tools.h b/include/dedx_tools.h index a453d4e..fb58430 100644 --- a/include/dedx_tools.h +++ b/include/dedx_tools.h @@ -35,10 +35,17 @@ double dedx_get_csda(dedx_workspace *ws, dedx_config *config, float energy, int * Inverts the stopping power curve. Since stopping power is non-monotonic, * the @p side parameter selects which branch to use. * - * @param[in] ws Workspace with a loaded configuration. - * @param[in] config Loaded configuration. + * @note This function (re)loads @p config into @p ws internally, so the + * configuration need not be pre-loaded; config->ion_a must be set. + * Because it always loads, do not pre-load @p config into a + * single-dataset workspace, or the dataset will be loaded twice. + * + * @param[in] ws Workspace into which the configuration is loaded. + * @param[in] config Configuration to invert (loaded internally). * @param[in] stp Target stopping power in MeV cm²/g. - * @param[in] side 0 = low-energy branch, 1 = high-energy branch. + * @param[in] side Branch selector: negative selects the low-energy branch + * (below the Bragg peak); zero or positive selects the + * high-energy branch (above the Bragg peak). * @param[out] err Error code; 0 on success. * @return Energy in MeV/nucl (MeV per nucleon). */ @@ -54,6 +61,22 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, */ double dedx_get_inverse_csda(dedx_workspace *ws, dedx_config *config, float range, int *err); +/** @brief Return the maximum (Bragg-peak) stopping power. + * + * Locates the energy at which the mass stopping power reaches its maximum + * (the Bragg peak) and returns the stopping power value there. + * + * @note If @p config is not already loaded into @p ws it is loaded + * automatically; config->ion_a must be set beforehand. + * + * @param[in] ws Workspace used to evaluate the configuration. + * @param[in] config Configuration to evaluate; loaded automatically if not + * already loaded. config->ion_a must be set. + * @param[out] err Error code; 0 on success. + * @return Peak mass stopping power in MeV cm²/g. + */ +double dedx_get_bragg_peak_stp(dedx_workspace *ws, dedx_config *config, int *err); + /** @brief Convert an array of stopping power values between unit systems. * * @param[in] old_unit Source unit (dedx_stp_units). diff --git a/include/dedx_wrappers.h b/include/dedx_wrappers.h index 31b59c2..1cab8c7 100644 --- a/include/dedx_wrappers.h +++ b/include/dedx_wrappers.h @@ -97,4 +97,49 @@ int dedx_get_csda_range_table(const int program, const float *energies, double *csda_ranges); +/** @brief One-call inverse CSDA: energy for a given CSDA range. + * + * Allocates a fresh workspace and configuration internally, sets the nucleon + * number automatically, and is safe to call without managing object + * lifetimes (suitable for JavaScript/WASM and other flat-API consumers). + * + * @param[in] program Program identifier. + * @param[in] ion Ion identifier. + * @param[in] target Target material identifier. + * @param[in] range CSDA range in g/cm². + * @param[out] err Error code; 0 on success. + * @return Kinetic energy in MeV/nucl (MeV per nucleon), or -1 on failure. + */ +double dedx_get_inverse_csda_simple(int program, int ion, int target, double range, int *err); + +/** @brief One-call inverse stopping power: energy for a given stopping power. + * + * Allocates a fresh workspace and configuration internally and sets the + * nucleon number automatically. Because stopping power is non-monotonic the + * @p side parameter selects which branch to invert (see dedx_get_inverse_stp). + * + * @param[in] program Program identifier. + * @param[in] ion Ion identifier. + * @param[in] target Target material identifier. + * @param[in] stp Target stopping power in MeV cm²/g. + * @param[in] side Branch selector (see dedx_get_inverse_stp). + * @param[out] err Error code; 0 on success. + * @return Kinetic energy in MeV/nucl (MeV per nucleon), or -1 on failure. + */ +double dedx_get_inverse_stp_simple(int program, int ion, int target, double stp, int side, int *err); + +/** @brief One-call Bragg peak: maximum stopping power for an ion/target pair. + * + * Allocates a fresh workspace and configuration internally and sets the + * nucleon number automatically, then returns the peak (maximum) mass + * stopping power. + * + * @param[in] program Program identifier. + * @param[in] ion Ion identifier. + * @param[in] target Target material identifier. + * @param[out] err Error code; 0 on success. + * @return Peak mass stopping power in MeV cm²/g, or -1 on failure. + */ +double dedx_get_bragg_peak_stp_simple(int program, int ion, int target, int *err); + #endif // DEDX_WRAPPERS_H diff --git a/src/dedx_tools.c b/src/dedx_tools.c index bb74e48..3323acf 100644 --- a/src/dedx_tools.c +++ b/src/dedx_tools.c @@ -164,6 +164,31 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, return (x1 + x2) / 2; } +double dedx_get_bragg_peak_stp(dedx_workspace *ws, dedx_config *config, int *err) { + if (config->ion_a <= 0) { + *err = DEDX_ERR_ION_A_REQUIRED; + return -1; + } + double acc = 1e-5; + dedx_tools_settings set; + + if (*err != 0) + return -1; + if (config->loaded == 0) + dedx_load_config(ws, config, err); + if (*err != 0) + return -1; + set.ws = ws; + set.cfg = config; + + double peak_energy = find_min(find_min_stp_func, &set, acc * 100); + if (peak_energy < 0) { + *err = DEDX_ERR_ENERGY_OUT_OF_RANGE; + return -1; + } + return dedx_get_stp(ws, config, (float) peak_energy, err); +} + double dedx_get_csda(dedx_workspace *ws, dedx_config *config, float energy, int *err) { if (config->ion_a <= 0) { *err = DEDX_ERR_ION_A_REQUIRED; diff --git a/src/dedx_wrappers.c b/src/dedx_wrappers.c index c993a54..d674d30 100644 --- a/src/dedx_wrappers.c +++ b/src/dedx_wrappers.c @@ -253,3 +253,98 @@ int dedx_get_csda_range_table(const int program, return err; } + +double dedx_get_inverse_csda_simple(int program, int ion, int target, double range, int *err) { + int cleanup_err = DEDX_OK; + double energy; + dedx_config *config = allocate_wrapper_config(program, ion, target, err); + dedx_workspace *ws; + + if (*err != 0) + return -1.0; + config->ion_a = dedx_internal_get_nucleon(config->ion, err); + if (*err != 0) { + dedx_free_config(config, &cleanup_err); + return -1.0; + } + ws = allocate_wrapper_workspace(err); + if (*err != 0) { + dedx_free_config(config, &cleanup_err); + return -1.0; + } + dedx_load_config(ws, config, err); + if (*err != 0) { + dedx_free_config(config, &cleanup_err); + dedx_free_workspace(ws, &cleanup_err); + return -1.0; + } + + energy = dedx_get_inverse_csda(ws, config, (float) range, err); + dedx_free_config(config, &cleanup_err); + dedx_free_workspace(ws, &cleanup_err); + if (*err != 0) + return -1.0; + return energy; +} + +double dedx_get_inverse_stp_simple(int program, int ion, int target, double stp, int side, int *err) { + int cleanup_err = DEDX_OK; + double energy; + dedx_config *config = allocate_wrapper_config(program, ion, target, err); + dedx_workspace *ws; + + if (*err != 0) + return -1.0; + config->ion_a = dedx_internal_get_nucleon(config->ion, err); + if (*err != 0) { + dedx_free_config(config, &cleanup_err); + return -1.0; + } + ws = allocate_wrapper_workspace(err); + if (*err != 0) { + dedx_free_config(config, &cleanup_err); + return -1.0; + } + + /* dedx_get_inverse_stp() loads the configuration internally, so the + workspace (allocated for a single dataset) must not be pre-loaded here. */ + energy = dedx_get_inverse_stp(ws, config, (float) stp, side, err); + dedx_free_config(config, &cleanup_err); + dedx_free_workspace(ws, &cleanup_err); + if (*err != 0) + return -1.0; + return energy; +} + +double dedx_get_bragg_peak_stp_simple(int program, int ion, int target, int *err) { + int cleanup_err = DEDX_OK; + double peak_stp; + dedx_config *config = allocate_wrapper_config(program, ion, target, err); + dedx_workspace *ws; + + if (*err != 0) + return -1.0; + config->ion_a = dedx_internal_get_nucleon(config->ion, err); + if (*err != 0) { + dedx_free_config(config, &cleanup_err); + return -1.0; + } + ws = allocate_wrapper_workspace(err); + if (*err != 0) { + dedx_free_config(config, &cleanup_err); + return -1.0; + } + dedx_load_config(ws, config, err); + if (*err != 0) { + dedx_free_config(config, &cleanup_err); + dedx_free_workspace(ws, &cleanup_err); + return -1.0; + } + + peak_stp = dedx_get_bragg_peak_stp(ws, config, err); + dedx_free_config(config, &cleanup_err); + dedx_free_workspace(ws, &cleanup_err); + if (*err != 0) + return -1.0; + return peak_stp; +} diff --git a/tests/test_tools.c b/tests/test_tools.c index 9e8269d..c1d8c9c 100644 --- a/tests/test_tools.c +++ b/tests/test_tools.c @@ -30,6 +30,7 @@ int main(void) { double csda; double inverse_csda; double stp; + double bragg_peak; float old_values[2] = {1.0f, 2.0f}; float new_values[2] = {0.0f, 0.0f}; @@ -73,6 +74,22 @@ int main(void) { failures += failf("dedx_get_stp value", stp, 1.0); } + /* The Bragg-peak stopping power is the maximum of the curve, so it must + exceed the stopping power at 100 MeV (far up the high-energy branch). */ + bragg_peak = dedx_get_bragg_peak_stp(ws, cfg, &err); + if (err != DEDX_OK) { + failures += faili("dedx_get_bragg_peak_stp err", err, DEDX_OK); + } else if (!(bragg_peak > stp)) { + failures += failf("dedx_get_bragg_peak_stp value", bragg_peak, stp); + } + + /* Without a nucleon number the Bragg-peak helper must refuse to run. */ + cfg->ion_a = 0; + err = DEDX_OK; + bragg_peak = dedx_get_bragg_peak_stp(ws, cfg, &err); + failures += faili("dedx_get_bragg_peak_stp ion_a guard", err, DEDX_ERR_ION_A_REQUIRED); + cfg->ion_a = 1; + err = convert_units(DEDX_MEVCM2G, DEDX_KEVUM, DEDX_WATER, 2, old_values, new_values); if (err != DEDX_OK) { failures += faili("convert_units err", err, DEDX_OK); @@ -97,6 +114,30 @@ int main(void) { failures += failf("convert_units reverse second", new_values[1], 20.0); } + /* The Bragg-peak helper loads an unloaded configuration automatically, as + documented; verify that documented behavior on a fresh config. */ + { + int err2 = 0; + dedx_config *cfg2 = calloc(1, sizeof(dedx_config)); + dedx_workspace *ws2 = dedx_allocate_workspace(1, &err2); + if (cfg2 != NULL && ws2 != NULL && err2 == DEDX_OK) { + cfg2->program = DEDX_PSTAR; + cfg2->ion = DEDX_PROTON; + cfg2->target = DEDX_WATER; + cfg2->ion_a = 1; + /* deliberately not loaded; the helper must load it itself */ + bragg_peak = dedx_get_bragg_peak_stp(ws2, cfg2, &err2); + if (err2 != DEDX_OK) + failures += faili("bragg auto-load err", err2, DEDX_OK); + else if (!(bragg_peak > 0.0)) + failures += failf("bragg auto-load value", bragg_peak, 1.0); + } else { + failures += faili("bragg auto-load setup", err2, DEDX_OK); + } + dedx_free_config(cfg2, &err2); + dedx_free_workspace(ws2, &err2); + } + dedx_free_config(cfg, &err); dedx_free_workspace(ws, &err); return failures; diff --git a/tests/test_wrappers.c b/tests/test_wrappers.c new file mode 100644 index 0000000..7ede98f --- /dev/null +++ b/tests/test_wrappers.c @@ -0,0 +1,103 @@ +#include +#include +#include +#include +#include + +static int failf(const char *label, double got, double expected) { + fprintf(stderr, "FAIL %s: got %.8g expected %.8g\n", label, got, expected); + return 1; +} + +static int faili(const char *label, int got, int expected) { + if (got == expected) + return 0; + fprintf(stderr, "FAIL %s: got %d expected %d\n", label, got, expected); + return 1; +} + +static int approx(double a, double b, double rel) { + double scale = fabs(b) > 1e-12 ? fabs(b) : 1.0; + return fabs(a - b) <= rel * scale; +} + +int main(void) { + int failures = 0; + int err = 0; + const float energy = 100.0f; + double range = 0.0; + float stp_high = 0.0f; + double recovered; + double peak; + + /* Reference CSDA range at 100 MeV/nucl for protons in water. */ + err = dedx_get_csda_range_table(DEDX_PSTAR, DEDX_PROTON, DEDX_WATER, 1, &energy, &range); + failures += faili("csda_range_table err", err, DEDX_OK); + + /* Inverse CSDA flat wrapper should recover the original energy. */ + err = 0; + recovered = dedx_get_inverse_csda_simple(DEDX_PSTAR, DEDX_PROTON, DEDX_WATER, range, &err); + if (err != DEDX_OK) { + failures += faili("inverse_csda_simple err", err, DEDX_OK); + } else if (!approx(recovered, energy, 1e-3)) { + failures += failf("inverse_csda_simple energy", recovered, energy); + } + + /* Reference stopping power at 100 MeV/nucl (falling, high-energy branch). */ + err = 0; + stp_high = dedx_get_simple_stp_for_program(DEDX_PSTAR, DEDX_PROTON, DEDX_WATER, energy, &err); + failures += faili("simple_stp err", err, DEDX_OK); + + /* Inverse STP flat wrapper on the high-energy branch recovers the energy. */ + err = 0; + recovered = dedx_get_inverse_stp_simple(DEDX_PSTAR, DEDX_PROTON, DEDX_WATER, stp_high, 1, &err); + if (err != DEDX_OK) { + failures += faili("inverse_stp_simple err", err, DEDX_OK); + } else if (!approx(recovered, energy, 1e-2)) { + failures += failf("inverse_stp_simple energy", recovered, energy); + } + + /* The low-energy branch (side < 0) yields the rising-edge solution, which + lies below the high-energy solution for the same stopping power. */ + err = 0; + { + double low_branch = dedx_get_inverse_stp_simple(DEDX_PSTAR, DEDX_PROTON, DEDX_WATER, stp_high, -1, &err); + if (err != DEDX_OK) { + failures += faili("inverse_stp_simple low-branch err", err, DEDX_OK); + } else if (!(low_branch > 0.0 && low_branch < recovered)) { + failures += failf("inverse_stp_simple low-branch energy", low_branch, recovered); + } + } + + /* Bragg peak stopping power must exceed the stopping power at 100 MeV. */ + err = 0; + peak = dedx_get_bragg_peak_stp_simple(DEDX_PSTAR, DEDX_PROTON, DEDX_WATER, &err); + if (err != DEDX_OK) { + failures += faili("bragg_peak_stp_simple err", err, DEDX_OK); + } else if (!(peak > stp_high)) { + failures += failf("bragg_peak_stp_simple value", peak, stp_high); + } + + /* Error path: an unsupported program/ion combination (PSTAR is + proton-only) must report DEDX_ERR_ION_NOT_SUPPORTED and return -1 from + every flat wrapper, exercising their allocation-cleanup branches. */ + err = 0; + recovered = dedx_get_inverse_csda_simple(DEDX_PSTAR, DEDX_CARBON, DEDX_WATER, 1.0, &err); + failures += faili("inverse_csda_simple unsupported err", err, DEDX_ERR_ION_NOT_SUPPORTED); + if (!approx(recovered, -1.0, 1e-9)) + failures += failf("inverse_csda_simple unsupported return", recovered, -1.0); + + err = 0; + recovered = dedx_get_inverse_stp_simple(DEDX_PSTAR, DEDX_CARBON, DEDX_WATER, 10.0, 1, &err); + failures += faili("inverse_stp_simple unsupported err", err, DEDX_ERR_ION_NOT_SUPPORTED); + if (!approx(recovered, -1.0, 1e-9)) + failures += failf("inverse_stp_simple unsupported return", recovered, -1.0); + + err = 0; + recovered = dedx_get_bragg_peak_stp_simple(DEDX_PSTAR, DEDX_CARBON, DEDX_WATER, &err); + failures += faili("bragg_peak_stp_simple unsupported err", err, DEDX_ERR_ION_NOT_SUPPORTED); + if (!approx(recovered, -1.0, 1e-9)) + failures += failf("bragg_peak_stp_simple unsupported return", recovered, -1.0); + + return failures; +}