feat: add flat inverse / Bragg-peak convenience wrappers (#122)

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).

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

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;

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;
+}

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;

tests/test_wrappers.c

@@ -0,0 +1,103 @@
+#include <dedx.h>
+#include <dedx_wrappers.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+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;
+}