diff --git a/include/dedx_tools.h b/include/dedx_tools.h index 9f761f5..5706235 100644 --- a/include/dedx_tools.h +++ b/include/dedx_tools.h @@ -36,18 +36,65 @@ double dedx_get_csda(dedx_workspace *ws, dedx_config *config, float energy, int /** @brief Find the energy corresponding to a given stopping power value. * - * Inverts the stopping power curve. Since stopping power is non-monotonic, - * the @p side parameter selects which branch to use. + * Inverts the stopping power curve. Stopping power is not simply unimodal: + * real tables can rise to a maximum (the Bragg peak) at low energy, fall + * through the minimum-ionizing point, and rise again at relativistic + * energies (e.g. proton tables extending to several GeV), so a requested + * @p stp can be reachable at more than one energy. Among all reachable + * energies, @p side selects which one to return: 0 selects the lowest; + * any other value (the documented convention is 1) selects the highest. + * When only one energy reaches @p stp, @p side has no effect. * - * @param[in] ws Workspace with a loaded configuration. - * @param[in] config Loaded configuration. - * @param[in] stp Target stopping power in MeV cm²/g. - * @param[in] side 0 = low-energy branch, 1 = high-energy branch. + * @param[in] ws Workspace to hold the loaded configuration. + * @param[in] config Configuration to invert against. Loaded automatically + * if not already (config->loaded == 0): on return, + * config->loaded and config->cfg_id reflect that load, + * same as after an explicit dedx_load_config() call. + * @param[in] stp Target stopping power in MeV cm²/g. In practice must + * lie between dedx_get_min_stp() and dedx_get_max_stp() + * (inclusive) for this configuration -- both computed + * from the same tabulated data points this function + * searches -- or the value is unreachable and an error + * is returned. + * @param[in] side 0 = lowest-energy solution, any other value = + * highest-energy solution. * @param[out] err Error code; 0 on success. - * @return Energy in MeV/nucl (MeV per nucleon). + * @return Energy in MeV/nucl (MeV per nucleon), or -1 on error. */ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, int side, int *err); +/** @brief Find the maximum stopping power over a program's whole tabulated + * energy range. + * + * This is the peak of the stopping-power-vs-energy curve, not the Bragg + * peak of a depth-dose curve (which also depends on range straggling and is + * not computed by this library). It is computed exactly from the tabulated + * data points backing the configuration, not sampled or estimated. + * + * @param[in] ws Workspace to hold the loaded configuration. + * @param[in] config Configuration to inspect. Loaded automatically if not + * already (config->loaded == 0); see dedx_get_inverse_stp(). + * @param[out] err Error code; 0 on success. + * @return Maximum stopping power in MeV cm²/g, or -1 on error. + */ +double dedx_get_max_stp(dedx_workspace *ws, dedx_config *config, int *err); + +/** @brief Find the minimum stopping power over a program's whole tabulated + * energy range. + * + * For tables that reach relativistic energies this is typically the + * minimum-ionizing point, not the value at either tabulated endpoint (the + * curve can rise again past it — see dedx_get_inverse_stp()). Computed + * exactly from the tabulated data points, like dedx_get_max_stp(). + * + * @param[in] ws Workspace to hold the loaded configuration. + * @param[in] config Configuration to inspect. Loaded automatically if not + * already (config->loaded == 0); see dedx_get_inverse_stp(). + * @param[out] err Error code; 0 on success. + * @return Minimum stopping power in MeV cm²/g, or -1 on error. + */ +double dedx_get_min_stp(dedx_workspace *ws, dedx_config *config, int *err); + /** @brief Find the energy corresponding to a given CSDA range. * * @param[in] ws Workspace with a loaded configuration. diff --git a/src/dedx_tools.c b/src/dedx_tools.c index bb74e48..adb4d89 100644 --- a/src/dedx_tools.c +++ b/src/dedx_tools.c @@ -5,6 +5,7 @@ #include #include "dedx_data_access.h" +#include "dedx_lookup_data.h" typedef struct { dedx_workspace *ws; @@ -59,42 +60,66 @@ static double adapt_stp(double energy, dedx_tools_settings *set) { return 1.0 / stp; } -static double find_min_stp_func(double x, dedx_tools_settings *set) { - int err = 0; - double stp = dedx_get_stp(set->ws, set->cfg, x, &err); - if (err != 0 || stp == 0.0) - return INFINITY; - return 1.0 / stp; +/* Load (if needed) and return the internal dataset backing a config, so + * callers can walk its exact tabulated (energy, STP) knots directly instead + * of re-sampling the curve at an arbitrary density. dedx_get_stp() evaluates + * a spline built from precisely these knots (dedx_spline.c: coef[i].a = + * stopping[i]), so scanning them finds every feature present in the + * tabulated values -- though not necessarily every feature of the + * interpolated curve dedx_get_stp() actually evaluates, since a spline can + * in principle overshoot slightly between two knots. That's a property of + * the interpolation, not something this scan (or any sampling density) can + * see; it hasn't been observed in the tables this library ships. + * + * `static` (file-private, not declared in dedx_tools.h): only this + * translation unit needs it, unlike the dedx_internal_* functions declared + * in headers such as dedx_data_access.h for use across multiple .c files. */ +static dedx_internal_lookup_data *get_loaded_dataset(dedx_workspace *ws, dedx_config *config, int *err) { + /* Public entrypoints reset *err themselves (e.g. dedx_get_stp()) rather + * than treating an incoming nonzero value as a precondition -- do the + * same here so a stale error code left over in the caller's variable + * from an unrelated earlier call can't make a valid request fail. */ + *err = DEDX_OK; + + /* config->loaded == 1 only means *some* workspace holds this config's + * data at config->cfg_id -- not necessarily this ws (the config could + * have been loaded into a workspace that was since freed and replaced, + * or the caller could be reusing one dedx_config across workspaces). + * Reload whenever the config hasn't been loaded at all, or its cfg_id + * isn't a valid slot in *this* workspace; otherwise skip the reload on + * the common already-loaded-here path. */ + int needs_load = 0; + if (config->loaded == 0) + needs_load = 1; + else if (config->cfg_id < 0 || config->cfg_id >= ws->active_datasets) + needs_load = 1; + if (needs_load) + dedx_load_config(ws, config, err); + if (*err != 0) + return NULL; + + int id = config->cfg_id; + if (id < 0 || id >= ws->active_datasets) { + /* Defensive: dedx_load_config() succeeding is expected to always + * leave a valid cfg_id, so this should be unreachable in practice. */ + *err = DEDX_ERR_INVALID_DATASET_ID; + return NULL; + } + return ws->loaded_data[id]; } -static double find_min(double (*func)(double x, dedx_tools_settings *set), dedx_tools_settings *set, double acc) { - double x[] = {0.01, 10}; - double f[] = {func(x[0], set), func(x[1], set)}; - int i = 0; - double x_temp = 0; - double f_temp = 0; - double h = 0; - while (fabs(x[1] - x[0]) > acc) { - i = 0; - h = x[1] - x[0]; - if (f[1] > f[0]) - i = 1; - // try flip - - x_temp = x[i] + 2 * h * pow(-1, i); - f_temp = func(x_temp, set); - if (f_temp > f[i]) { - // try shrink - x_temp = x[i] + 0.5 * h * pow(-1, i); - f_temp = func(x_temp, set); - - if (f_temp > f[i]) - return -1; - } - f[i] = f_temp; - x[i] = x_temp; +static void min_max_stp_over_table(const dedx_internal_lookup_data *data, double *min_stp, double *max_stp) { + double lo = data->base[0].a; + double hi = data->base[0].a; + for (int i = 1; i < data->n; i++) { + double v = data->base[i].a; + if (v < lo) + lo = v; + if (v > hi) + hi = v; } - return (x[0] + x[1]) / 2; + *min_stp = lo; + *max_stp = hi; } double dedx_get_inverse_csda(dedx_workspace *ws, dedx_config *config, float range, int *err) { @@ -124,44 +149,222 @@ double dedx_get_inverse_csda(dedx_workspace *ws, dedx_config *config, float rang return (min + max) / 2; } +/* Bisection convergence, as an absolute tolerance in *log-energy* space (see + * bisect_monotonic_run() below for why bisecting there instead of in the raw + * energy value). For small differences, d(ln x) = dx/x, i.e. a difference in + * log-space is directly a *relative* difference in the original energy -- + * so this one constant fixes the achieved relative precision on the energy + * result regardless of scale. 1e-6 is about ten times coarser than a + * `float`'s own precision (~1.19e-7 relative, i.e. 24 bits of mantissa): + * dedx_get_stp() truncates the search value to a float before evaluating + * the spline, so refining much further than the float itself can represent + * would just spend extra iterations re-measuring rounding noise. */ +#define DEDX_INVERSE_STP_LOG_ACC 1e-6 + +/* Bisect the monotonic run of knots [lo_idx, hi_idx] for the energy where + * the spline equals stp, refining via dedx_get_stp() (not linear knot + * interpolation) so the result matches what callers would measure with the + * public curve. Returns 0 and sets *solution on success; returns -1 without + * touching *err if stp is not bracketed by this run. */ +static int bisect_monotonic_run(dedx_workspace *ws, + dedx_config *config, + const dedx_internal_lookup_data *data, + int lo_idx, + int hi_idx, + float stp, + int *err, + double *solution) { + double lo_x = data->base[lo_idx].x; + double hi_x = data->base[hi_idx].x; + double lo_v = data->base[lo_idx].a; + double hi_v = data->base[hi_idx].a; + + /* This run's STP values span [range_min, range_max] regardless of + * whether the run is rising or falling; stp must land in that span for + * a solution to exist anywhere on this run. */ + double range_min; + double range_max; + if (lo_v < hi_v) { + range_min = lo_v; + range_max = hi_v; + } else { + range_min = hi_v; + range_max = lo_v; + } + if (stp < range_min || stp > range_max) + return -1; + + /* Whether STP rises (ascending) or falls (descending) from lo_x to + * hi_x tells us which half to keep at each bisection step below. */ + int ascending; + if (hi_v >= lo_v) + ascending = 1; + else + ascending = 0; + + /* Bisect in log-energy space rather than raw energy. Two reasons: + * + * 1. The tabulated energy grid is itself predominantly log-spaced (e.g. + * NIST PSTAR/ASTAR-derived tables step by a roughly constant energy + * *ratio*, not a constant energy difference, so that a fixed number + * of points can cover ~0.001 to 10000 MeV/nucleon with even relative + * resolution). Bisecting in the same log-scale as the data matches + * the coordinate the curve actually varies smoothly in. + * + * 2. It makes the convergence tolerance exactly a relative tolerance on + * the energy result (see DEDX_INVERSE_STP_LOG_ACC above), rather than + * an approximation of one. A run spanning many decades (e.g. the + * post-peak run for proton/water PSTAR, ~0.08 to 10000 MeV/nucleon) + * then converges in a number of steps set by log2(ln(hi_x/lo_x) / + * DEDX_INVERSE_STP_LOG_ACC) -- a few dozen steps regardless of scale + * or of where in the run the root happens to fall, since bisection + * always halves the bracket every step no matter which half contains + * the root. + * + * dedx_get_stp() evaluates the actual interpolating spline at any energy + * (not just at knots), so this is correct even on the tail of a table + * where the raw knot spacing itself happens to be linear rather than + * log-spaced (observed for a few tables at their highest energies) -- + * the spline doesn't care how its own knots were spaced, and neither + * does bisecting the search variable in log space. */ + double log_lo = log(lo_x); + double log_hi = log(hi_x); + while (fabs(log_lo - log_hi) > DEDX_INVERSE_STP_LOG_ACC) { + double log_mid = (log_lo + log_hi) / 2; + double x_temp = exp(log_mid); + double f_temp = dedx_get_stp(ws, config, (float) x_temp, err); + if (*err != 0) + return -1; + if (ascending) { + /* STP too low at the midpoint -> the root is further up. */ + if (f_temp <= stp) + log_lo = log_mid; + else + log_hi = log_mid; + } else { + /* STP still at/above target at the midpoint -> the root is + * further up (we're descending, so STP keeps falling as energy + * rises). */ + if (f_temp >= stp) + log_lo = log_mid; + else + log_hi = log_mid; + } + } + *solution = exp((log_lo + log_hi) / 2); + return 0; +} + double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, int side, int *err) { if (config->ion_a <= 0) { *err = DEDX_ERR_ION_A_REQUIRED; return -1; } - double acc = 1e-5; - dedx_tools_settings set; - set.ws = ws; - if (*err != 0) + dedx_internal_lookup_data *data = get_loaded_dataset(ws, config, err); + if (data == NULL) return -1; - dedx_load_config(ws, config, err); - set.cfg = config; - if (*err != 0) - return -1; - double max = find_min(find_min_stp_func, &set, acc * 100); - double x1; - double x2; - double x_temp; - double f_temp; - if (side < 0) { - x1 = dedx_get_min_energy(config->program, config->ion); - x2 = max; - } else { - x2 = max; - x1 = dedx_get_max_energy(config->program, config->ion); - } - while (fabs(x1 - x2) > acc) { - x_temp = (x1 + x2) / 2; - f_temp = dedx_get_stp(set.ws, set.cfg, x_temp, err); + /* Stopping power vs. energy is not simply unimodal: real tables can rise + * to the Bragg peak, fall through the minimum-ionizing point, and rise + * again at relativistic energies (e.g. PSTAR/ICRU protons up to 10 GeV) + * -- an arbitrary requested STP can be reachable on more than one of + * these monotonic runs. Walk the exact tabulated knots once, bisecting + * every run that brackets stp, and keep the lowest- and highest-energy + * solutions found across ALL of them; side == 0 returns the lowest, any + * other side value returns the highest. For the common single-peak case + * this is exactly the old ascending/descending branch choice. See #121. */ + int found = 0; + double x_min_found = 0; + double x_max_found = 0; + int seg_start = 0; + int prev_dir = 0; - if (f_temp >= stp) { - x2 = x_temp; - } else { - x1 = x_temp; + for (int i = 1; i <= data->n; i++) { + /* dir is the sign of the step from knot i-1 to knot i: +1 rising, + * -1 falling, 0 flat (a flat step doesn't end the current run -- + * prev_dir is left untouched so a flat plateau stays part of + * whichever run it interrupts). Only computed for i < data->n, + * since there is no knot i to compare against once i == data->n. */ + int is_turning = 0; + int dir = 0; /* stays 0 (flat) unless one of the branches below fires */ + if (i < data->n) { + double delta = (double) data->base[i].a - (double) data->base[i - 1].a; + if (delta > 0) + dir = 1; + else if (delta < 0) + dir = -1; + if (dir != 0) { + if (prev_dir == 0) + prev_dir = dir; + else if (dir != prev_dir) + is_turning = 1; + } + } + /* A run ends either at a direction reversal (is_turning, and the + * turning knot i-1 is shared with the next run) or at the last + * knot (i == data->n, since there's nowhere further to walk). */ + if (i == data->n || is_turning) { + int seg_end; + if (i == data->n) + seg_end = data->n - 1; + else + seg_end = i - 1; + if (seg_end > seg_start) { + double solution; + if (bisect_monotonic_run(ws, config, data, seg_start, seg_end, stp, err, &solution) == 0) { + /* This run brackets stp, so it contributes one solution. + * Track the lowest- and highest-energy solutions seen so + * far across every run processed so far in this loop -- + * not just this one run -- so that once the whole table + * has been walked, x_min_found/x_max_found hold the + * overall lowest/highest reachable energy regardless of + * how many separate runs bracketed stp along the way. + * `!found` covers the very first solution: there is + * nothing yet to compare it against, so it always + * becomes both the running min and the running max. */ + if (!found || solution < x_min_found) + x_min_found = solution; + if (!found || solution > x_max_found) + x_max_found = solution; + found = 1; + } else if (*err != 0) { + return -1; + } + } + if (is_turning) { + seg_start = i - 1; + prev_dir = dir; + } } } - return (x1 + x2) / 2; + + if (!found) { + *err = DEDX_ERR_ENERGY_OUT_OF_RANGE; + return -1; + } + if (side == 0) + return x_min_found; + return x_max_found; +} + +double dedx_get_max_stp(dedx_workspace *ws, dedx_config *config, int *err) { + dedx_internal_lookup_data *data = get_loaded_dataset(ws, config, err); + if (data == NULL) + return -1; + double min_stp; + double max_stp; + min_max_stp_over_table(data, &min_stp, &max_stp); + return max_stp; +} + +double dedx_get_min_stp(dedx_workspace *ws, dedx_config *config, int *err) { + dedx_internal_lookup_data *data = get_loaded_dataset(ws, config, err); + if (data == NULL) + return -1; + double min_stp; + double max_stp; + min_max_stp_over_table(data, &min_stp, &max_stp); + return min_stp; } double dedx_get_csda(dedx_workspace *ws, dedx_config *config, float energy, int *err) { diff --git a/tests/test_inverse_stp.c b/tests/test_inverse_stp.c new file mode 100644 index 0000000..013f033 --- /dev/null +++ b/tests/test_inverse_stp.c @@ -0,0 +1,298 @@ +#include +#include +#include +#include +#include +#include + +/* Regression tests for dedx_get_inverse_stp(), dedx_get_max_stp(), and + * dedx_get_min_stp() (issue #121). + * + * The stopping-power-vs-energy curve is not simply unimodal: PSTAR/ICRU + * proton tables extend to 10 GeV/nucleon, far enough to show the curve rise + * to the Bragg peak (~0.08 MeV), fall through a minimum-ionizing point + * (~3000-4000 MeV), and rise again (relativistic rise / Fermi plateau) up to + * the table's max energy. A requested STP in that high-energy dip can be + * reachable at two distinct energies, on top of the usual two solutions + * around the Bragg peak. dedx_get_inverse_stp() finds every monotonic run in + * the exact tabulated knots (no arbitrary sampling density) and returns the + * lowest-energy (side=0) or highest-energy (side=1) reachable solution. + */ + +static int failures = 0; + +static void expect_int(const char *label, long got, long expected) { + if (got != expected) { + fprintf(stderr, "FAIL %s: got %ld expected %ld\n", label, got, expected); + failures++; + } +} + +static void expect_near(const char *label, double got, double expected, double rel) { + double scale = fabs(expected) > 1e-12 ? fabs(expected) : 1.0; + if (fabs(got - expected) > rel * scale) { + fprintf(stderr, "FAIL %s: got %.8g expected %.8g\n", label, got, expected); + failures++; + } +} + +static void expect_true(const char *label, int condition) { + if (!condition) { + fprintf(stderr, "FAIL %s\n", label); + failures++; + } +} + +static dedx_config *make_config(int program, int ion, int ion_a, int target) { + dedx_config *cfg = calloc(1, sizeof(dedx_config)); + cfg->program = program; + cfg->ion = ion; + cfg->target = target; + cfg->ion_a = ion_a; + return cfg; +} + +int main(void) { + int err = 0; + + /* --- STP below the Bragg peak's ascending floor: only the post-peak + * descending run reaches it, so side=0 and side=1 must agree on the + * same (only reachable) energy, with no error. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + expect_int("workspace alloc err", err, DEDX_OK); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + err = 0; + double e0 = dedx_get_inverse_stp(ws, cfg, 100.0f, 0, &err); + expect_int("low-floor stp side=0 err", err, DEDX_OK); + expect_true("low-floor stp side=0 positive", e0 > 0.0); + + err = 0; + double e1 = dedx_get_inverse_stp(ws, cfg, 100.0f, 1, &err); + expect_int("low-floor stp side=1 err", err, DEDX_OK); + expect_near("low-floor stp side=0/1 agree", e0, e1, 1e-3); + + int verify_err = 0; + double stp_at_e0 = dedx_get_stp(ws, cfg, (float) e0, &verify_err); + expect_int("low-floor roundtrip err", verify_err, DEDX_OK); + expect_near("low-floor roundtrip stp", stp_at_e0, 100.0, 1e-2); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- Around the Bragg peak, STP reachable on both the ascending + * (pre-peak) and descending (post-peak) runs: side=0 and side=1 must + * select distinct energies and each round-trip back to the requested + * STP. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + const float target_stp = 500.0f; + + err = 0; + double e_low = dedx_get_inverse_stp(ws, cfg, target_stp, 0, &err); + expect_int("branch stp side=0 err", err, DEDX_OK); + + err = 0; + double e_high = dedx_get_inverse_stp(ws, cfg, target_stp, 1, &err); + expect_int("branch stp side=1 err", err, DEDX_OK); + + expect_true("side=0 selects the lower energy", e_low < e_high); + + int verify_err = 0; + double stp_low = dedx_get_stp(ws, cfg, (float) e_low, &verify_err); + expect_int("side=0 roundtrip err", verify_err, DEDX_OK); + expect_near("side=0 roundtrip stp", stp_low, target_stp, 1e-2); + + verify_err = 0; + double stp_high = dedx_get_stp(ws, cfg, (float) e_high, &verify_err); + expect_int("side=1 roundtrip err", verify_err, DEDX_OK); + expect_near("side=1 roundtrip stp", stp_high, target_stp, 1e-2); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- Relativistic-rise dip around the minimum-ionizing point (~3000-4000 + * MeV for proton/water PSTAR): a target STP in this dip is reachable at + * two widely separated energies, one on each side of the minimum. This + * is exactly the region the old find_min()-based implementation got + * wrong -- it assumed the whole post-peak range was a single descending + * branch and rejected every STP here as out of range. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + const float target_stp = 2.05f; + + err = 0; + double e_low = dedx_get_inverse_stp(ws, cfg, target_stp, 0, &err); + expect_int("dip stp side=0 err", err, DEDX_OK); + + err = 0; + double e_high = dedx_get_inverse_stp(ws, cfg, target_stp, 1, &err); + expect_int("dip stp side=1 err", err, DEDX_OK); + + expect_true("dip side=0 selects the lower energy", e_low < e_high); + expect_true("dip side=0 energy is in the pre-minimum range", e_low > 100.0 && e_low < 3000.0); + expect_true("dip side=1 energy is in the post-minimum range", e_high > 3000.0 && e_high < 10000.0); + + int verify_err = 0; + double stp_low = dedx_get_stp(ws, cfg, (float) e_low, &verify_err); + expect_int("dip side=0 roundtrip err", verify_err, DEDX_OK); + expect_near("dip side=0 roundtrip stp", stp_low, target_stp, 1e-2); + + verify_err = 0; + double stp_high = dedx_get_stp(ws, cfg, (float) e_high, &verify_err); + expect_int("dip side=1 roundtrip err", verify_err, DEDX_OK); + expect_near("dip side=1 roundtrip stp", stp_high, target_stp, 1e-2); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- Out-of-range STP values (above the global maximum, or below the + * global minimum) must return an error rather than looping to a bogus or + * negative energy. The bounds are derived from dedx_get_max_stp()/ + * dedx_get_min_stp() themselves rather than hardcoded constants, so this + * test stays valid if the underlying tabulated data is ever updated. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + + err = 0; + double max_stp = dedx_get_max_stp(ws, cfg, &err); + expect_int("out-of-range setup: max stp err", err, DEDX_OK); + err = 0; + double min_stp = dedx_get_min_stp(ws, cfg, &err); + expect_int("out-of-range setup: min stp err", err, DEDX_OK); + + float stp_above_max = (float) (max_stp * 2.0); + float stp_below_min = (float) (min_stp * 0.5); + + err = 0; + double e_too_high = dedx_get_inverse_stp(ws, cfg, stp_above_max, 0, &err); + expect_int("stp above global max err", err, DEDX_ERR_ENERGY_OUT_OF_RANGE); + expect_near("stp above global max sentinel", e_too_high, -1.0, 1e-9); + + err = 0; + double e_too_low = dedx_get_inverse_stp(ws, cfg, stp_below_min, 0, &err); + expect_int("stp below global min err", err, DEDX_ERR_ENERGY_OUT_OF_RANGE); + expect_near("stp below global min sentinel", e_too_low, -1.0, 1e-9); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- A config that fails to load (here: a Bethe-type program with a + * custom target and no density set, DEDX_ERR_RHO_REQUIRED -- the same + * failure test_error_codes.c uses) must propagate that error cleanly out + * of all three functions that route through get_loaded_dataset(), not + * crash or silently return a bogus value. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_BETHE_EXT00, DEDX_PROTON, 1, 0); + + err = 0; + double e = dedx_get_inverse_stp(ws, cfg, 100.0f, 0, &err); + expect_true("unloadable config: inverse_stp err set", err != DEDX_OK); + expect_near("unloadable config: inverse_stp sentinel", e, -1.0, 1e-9); + + err = 0; + double max_stp = dedx_get_max_stp(ws, cfg, &err); + expect_true("unloadable config: max_stp err set", err != DEDX_OK); + expect_near("unloadable config: max_stp sentinel", max_stp, -1.0, 1e-9); + + err = 0; + double min_stp = dedx_get_min_stp(ws, cfg, &err); + expect_true("unloadable config: min_stp err set", err != DEDX_OK); + expect_near("unloadable config: min_stp sentinel", min_stp, -1.0, 1e-9); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- A config already loaded into a workspace that was since freed and + * replaced (config->loaded == 1, but config->cfg_id is not a valid slot + * in the new workspace) must be transparently reloaded into the new + * workspace rather than erroring out -- get_loaded_dataset() detects the + * stale cfg_id and reloads instead of trusting the loaded flag alone. */ + { + dedx_workspace *ws1 = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + + err = 0; + double max_stp_ws1 = dedx_get_max_stp(ws1, cfg, &err); + expect_int("stale workspace: first load err", err, DEDX_OK); + expect_true("stale workspace: cfg marked loaded", cfg->loaded != 0); + + dedx_free_workspace(ws1, &err); + + dedx_workspace *ws2 = dedx_allocate_workspace(1, &err); + err = 0; + double max_stp_ws2 = dedx_get_max_stp(ws2, cfg, &err); + expect_int("stale workspace: reload into ws2 err", err, DEDX_OK); + expect_near("stale workspace: same result after reload", max_stp_ws2, max_stp_ws1, 1e-9); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws2, &err); + } + + /* --- ion_a <= 0 is rejected up front, same contract as before. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 0, DEDX_WATER); + err = 0; + double e = dedx_get_inverse_stp(ws, cfg, 100.0f, 0, &err); + expect_int("ion_a required err", err, DEDX_ERR_ION_A_REQUIRED); + expect_near("ion_a required sentinel", e, -1.0, 1e-9); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- dedx_get_max_stp() / dedx_get_min_stp() report the true extrema of + * the curve: the Bragg peak and the minimum-ionizing point. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + + err = 0; + double max_stp = dedx_get_max_stp(ws, cfg, &err); + expect_int("max stp err", err, DEDX_OK); + + err = 0; + double min_stp = dedx_get_min_stp(ws, cfg, &err); + expect_int("min stp err", err, DEDX_OK); + + expect_true("min stp is less than max stp", min_stp < max_stp); + expect_true("min stp is positive", min_stp > 0.0); + + /* Brute-force reference: sample far more densely than the table + * itself and take the min/max, as an independent cross-check. */ + float emin = dedx_get_min_energy(cfg->program, cfg->ion); + float emax = dedx_get_max_energy(cfg->program, cfg->ion); + double log_emin = log((double) emin); + double log_emax = log((double) emax); + double reference_max = 0.0; + double reference_min = HUGE_VAL; + const int n = 5000; + for (int i = 0; i < n; i++) { + float e = (float) exp(log_emin + (log_emax - log_emin) * i / (n - 1)); + int sample_err = 0; + double s = dedx_get_stp(ws, cfg, e, &sample_err); + if (sample_err != 0) + continue; + if (s > reference_max) + reference_max = s; + if (s < reference_min) + reference_min = s; + } + expect_near("max stp matches brute-force reference", max_stp, reference_max, 5e-2); + expect_near("min stp matches brute-force reference", min_stp, reference_min, 5e-2); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + return failures; +}