Use energy loss for small CSDA steps

niess · Feb 27, 2023 · 529976e · 529976e
1 parent 955149b
commit 529976e
Show file tree

Hide file tree

Showing 2 changed files with 42 additions and 6 deletions.
diff --git a/include/pumas.h b/include/pumas.h
@@ -37,7 +37,7 @@ extern "C" {
 /* PUMAS library version. */
 #define PUMAS_VERSION_MAJOR 1
 #define PUMAS_VERSION_MINOR 2
-#define PUMAS_VERSION_PATCH 2
+#define PUMAS_VERSION_PATCH 3
 
 /**
  * Projectiles supported by PUMAS.

diff --git a/src/pumas.c b/src/pumas.c
@@ -93,6 +93,13 @@
  * Maximum deflection angle for a soft scattering event, in degrees.
  */
 #define MAX_SOFT_ANGLE 1E+00
+/**
+ * Grammage ratio for small steps, in CSDA mode.
+ *
+ * Below this ratio, the step energy loss is computed from the stopping power,
+ * using finite differences. Above, the interpolated grammage integral is used.
+ */
+#define EPSILON_X 3E-03
 /**
  * Minimum step size.
  */
@@ -4842,8 +4849,14 @@ enum pumas_event transport_with_csda(const struct pumas_physics * physics,
         if (context->mode.direction == PUMAS_MODE_FORWARD) {
                 if (xi > xB) {
                         xf = xi - xB;
-                        kf = cel_kinetic_energy(
-                            physics, context, PUMAS_MODE_CSDA, material, xf);
+                        if (xB < EPSILON_X * xi) {
+                                kf = ki - cel_energy_loss(physics, context,
+                                    PUMAS_MODE_CSDA, material, ki) * xB;
+                                if (kf < 0.) kf = 0.;
+                        } else {
+                                kf = cel_kinetic_energy(physics, context,
+                                    PUMAS_MODE_CSDA, material, xf);
+                        }
                 } else {
                         xf = kf = 0.;
                         event = PUMAS_EVENT_LIMIT_ENERGY;
@@ -4865,8 +4878,13 @@ enum pumas_event transport_with_csda(const struct pumas_physics * physics,
                         kf = DBL_MAX;
                 } else {
                         xf = xB + xi;
-                        kf = cel_kinetic_energy(
-                            physics, context, PUMAS_MODE_CSDA, material, xf);
+                        if (xB < EPSILON_X * xi) {
+                                kf = ki + cel_energy_loss(physics, context,
+                                    PUMAS_MODE_CSDA, material, ki) * xB;
+                        } else {
+                                kf = cel_kinetic_energy(physics, context,
+                                    PUMAS_MODE_CSDA, material, xf);
+                        }
                 }
                 if (context->event & PUMAS_EVENT_LIMIT_ENERGY) {
                         if (ki >= context->limit.energy)
@@ -6521,9 +6539,27 @@ enum pumas_return step_transport(const struct pumas_physics * physics,
                                 grammage_max = grammage;
                                 event = PUMAS_EVENT_LIMIT_ENERGY;
                         }
-                } else
+                } else if (dX < EPSILON_X * Xtot) {
+                        k1 -= sgn * cel_energy_loss(physics, context, scheme,
+                            material, k1) * dX;
+                        if ((context->event & PUMAS_EVENT_LIMIT_ENERGY) &&
+                            (context->limit.energy > 0.)) {
+                                if (((context->mode.direction ==
+                                     PUMAS_MODE_FORWARD) &&
+                                     (k1 < context->limit.energy)) ||
+                                    ((context->mode.direction ==
+                                    PUMAS_MODE_BACKWARD) &&
+                                     (k1 > context->limit.energy))) {
+                                        k1 = context->limit.energy;
+                                        event = PUMAS_EVENT_LIMIT_ENERGY;
+                                }
+                        } else if (k1 < 0.) {
+                                k1 = 0.;
+                        }
+                } else {
                         k1 = cel_kinetic_energy(
                             physics, context, scheme, material, X);
+                }
         } else if (scheme == PUMAS_MODE_STRAGGLED) {
                 /* Fluctuate the CEL around its average value. */
                 double ratio;