Patch a bug in the change of medium locating.

niess · Mar 21, 2017 · 34a498b · 34a498b
1 parent d3773b5
commit 34a498b
Show file tree

Hide file tree

Showing 2 changed files with 21 additions and 12 deletions.
diff --git a/include/ent.h b/include/ent.h
@@ -216,7 +216,7 @@ struct ent_context;
  * signs that the corresponding medium has no boundaries.
  *
  * **Warning** : it is an error to return zero or less for any state if the
- * extension is finite.
+ * initial medium is finite, i.e. returns a stricly positive step value.
  */
 typedef double ent_medium_cb(struct ent_context * context,
     struct ent_state * state, struct ent_medium ** medium);
@@ -438,13 +438,14 @@ enum ent_return ent_vertex(struct ent_physics * physics,
  * code is returned as detailed below.
  *
  * Transport the given particle in a medium. The Monte-Carlo transport ends
- * whenever an interaction occurs, or a medium boundary is crossed, or a user
+ * whenever an interaction occurs, or the neutrino escapes, or a user
  * specified limit is reached. Note that if *physics* is `NULL` neutrino
  * interactions are disabled.
  *
  * __Error codes__
  *
- *     ENT_RETURN_DOMAIN_ERROR     Some input parameter is inconsistent.
+ *     ENT_RETURN_DOMAIN_ERROR     Some input parameter is inconsistent or an
+ * inconsistent step value was returned.
  */
 enum ent_return ent_transport(struct ent_physics * physics,
     struct ent_context * context, struct ent_state * state,

diff --git a/src/ent.c b/src/ent.c
@@ -1299,9 +1299,9 @@ static enum ent_return transport_step(struct ent_context * context,
                         const double tmp_step =
                             context->medium(context, state, &tmp_medium);
                         if (tmp_medium == *medium) {
-                                s1 = s3;
-                        } else {
                                 s2 = s3;
+                        } else {
+                                s1 = s3;
                                 step1 = tmp_step;
                                 if (tmp_medium != medium1) medium1 = tmp_medium;
                         }
@@ -1314,21 +1314,17 @@ static enum ent_return transport_step(struct ent_context * context,
         }
 
         /* Get the end step local properties. */
-        double density1;
-        const double step2 = (*medium)->density(*medium, state, &density1);
-        if ((density1 < 0.) || ((*medium)->A <= 0.))
+        double density1, step2;
+        step2 = (*medium)->density(*medium, state, &density1);
+        if ((density1 <= 0.) || ((*medium)->A <= 0.))
                 return ENT_RETURN_DOMAIN_ERROR;
-        if (step2 > 0.) {
-                if ((step1 <= 0.) || (step2 < step1)) step1 = step2;
-        }
 
         /* Offset the end step position if a boundary crossing occured. */
         if (ds_boundary != 0.) {
                 state->position[0] += ds_boundary * sgn * state->direction[0];
                 state->position[1] += ds_boundary * sgn * state->direction[1];
                 state->position[2] += ds_boundary * sgn * state->direction[2];
                 *medium = medium1;
-                *step += ds_boundary;
 
                 /* Reset any distance limit if no more valid. */
                 if ((context->distance_max > 0.) &&
@@ -1366,7 +1362,15 @@ static enum ent_return transport_step(struct ent_context * context,
                 state->distance += *step;
         }
 
+        /* Update the end step density if a valid change of medium occured. */
+        if ((ds_boundary != 0.) && (*medium != NULL)) {
+                step2 = (*medium)->density(*medium, state, &density1);
+                if ((density1 <= 0.) || (medium1->A <= 0.))
+                        return ENT_RETURN_DOMAIN_ERROR;
+        }
+
         /* Update and return. */
+        if ((step2 > 0.) && ((step1 <= 0.) || (step2 < step1))) step1 = step2;
         *step = (step1 < 0) ? 0. : step1;
         *density = density1;
         return ENT_RETURN_SUCCESS;
@@ -2228,6 +2232,10 @@ enum ent_return ent_transport(struct ent_physics * physics,
                             ENT_RETURN_SUCCESS)
                                 goto exit;
                         if (event_ != ENT_EVENT_NONE) break;
+                        if ((medium != NULL) && (step == 0.)) {
+                                rc = ENT_RETURN_DOMAIN_ERROR;
+                                goto exit;
+                        }
                 }
         else {
                 /* This is a uniform medium of infinite extension.