conic.h: avoid spurious re-hit (and absorption) on the surface just reflected from

mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr

@@ -12,18 +12,18 @@
 * %D
  * Example instrument that shows some ways of using the Conics pair components
 *
- * %Example: Test_Conics_pairs OPTIC=1 Detector: psd_i_I=6.71864e-13
- * %Example: Test_Conics_pairs OPTIC=1 Detector: psd_i_reflections_2_I=6.71864e-13
- * %Example: Test_Conics_pairs OPTIC=2 fs=100000 Detector: psd_i_I=2.75574e-21
- * %Example: Test_Conics_pairs OPTIC=2 fs=100000 Detector: psd_i_reflections_2_I=2.75574e-21
- * %Example: Test_Conics_pairs OPTIC=3 Detector: psd_i_I=6.71864e-13
- * %Example: Test_Conics_pairs OPTIC=4 ssize=1e-4 Detector: psd_i_I=8.42691e-09
- * %Example: Test_Conics_pairs quadratic=1 Detector: psd_i_I=4.93027e-13
+ * %Example: Test_Conics_pairs OPTIC=1 Detector: psd_i_I=7.13655e-13
+ * %Example: Test_Conics_pairs OPTIC=1 Detector: psd_i_reflections_2_I=7.13655e-13
+ * %Example: Test_Conics_pairs OPTIC=2 fs=100000 Detector: psd_i_I=3.65133e-21
+ * %Example: Test_Conics_pairs OPTIC=2 fs=100000 Detector: psd_i_reflections_2_I=3.65072e-21
+ * %Example: Test_Conics_pairs OPTIC=3 Detector: psd_i_I=7.13655e-13
+ * %Example: Test_Conics_pairs OPTIC=4 ssize=1e-4 Detector: psd_i_I=8.44676e-09
+ * %Example: Test_Conics_pairs quadratic=1 Detector: psd_i_I=5.23944e-13
  * %Example: Test_Conics_pairs quadratic=1 fs=8 fi=6 Detector: psd_i_I=1.2236e-12
- * %Example: Test_Conics_pairs OPTIC=2 quadratic=1 fs=100000 Detector: psd_i_I=2.03618e-21
- * %Example: Test_Conics_pairs OPTIC=3 quadratic=1 Detector: psd_i_I=4.93027e-13
+ * %Example: Test_Conics_pairs OPTIC=2 quadratic=1 fs=100000 Detector: psd_i_I=2.71421e-21
+ * %Example: Test_Conics_pairs OPTIC=3 quadratic=1 Detector: psd_i_I=5.23944e-13
  * %Example: Test_Conics_pairs OPTIC=3 quadratic=1 fs=8 fi=6 Detector: psd_i_I=8.80918e-13
- * %Example: Test_Conics_pairs OPTIC=4 quadratic=1 ssize=1e-4 Detector: psd_i_I=6.41135e-09
+ * %Example: Test_Conics_pairs OPTIC=4 quadratic=1 ssize=1e-4 Detector: psd_i_I=6.42631e-09
 *
 * %P
 * OPTIC: [ ]     Flag to choose between 1: EH pair, 2: PH pair, 3: HE pair, 4: PP pair

mcstas-comps/share/conic.h

@@ -1449,6 +1449,47 @@ double getTimeOfFirstCollisionFlat(_class_particle p, FlatSurf s) {
         return -1;
     return t+tz;
 }
+
+/*! \brief Function to compute time of next collision for a ConicSurf the particle is on
+
+Used right after an interaction with s: the particle is then on the surface, i.e.
+C = 0 in the quadratic equation and the roots are t = 0 and t = -B/A. Using -B/A
+directly avoids spurious near-zero roots caused by round-off.
+
+@param p _class_particle to consider, located on s
+@param s ConicSurf to consider
+@return Time until the next collision or -1 if particle will not hit s again
+*/
+double getTimeOfNextCollisionConic(_class_particle p, ConicSurf s) {
+    double A = p.vx*p.vx+p.vy*p.vy-s.k3*p.vz*p.vz;
+    double B = 2*(p.vx*p.x+p.vy*p.y-s.k3*p.vz*p.z)-s.k2*p.vz;
+    if (A == 0)
+        return -1;
+    double t = -B/A;
+    if (t <= 0 || p.vz*t+p.z > s.ze || p.vz*t+p.z < s.zs)
+        return -1;
+    return t;
+}
+
+/*! \brief Function to compute time of next collision for a FlatSurf the particle is on
+
+@see getTimeOfNextCollisionConic()
+
+@param p _class_particle to consider, located on s
+@param s FlatSurf to consider
+@return Time until the next collision or -1 if particle will not hit s again
+*/
+double getTimeOfNextCollisionFlat(_class_particle p, FlatSurf s) {
+    double A = p.vx*p.vx-s.k3*p.vz*p.vz;
+    double B = 2*(p.vx*p.x-s.k3*p.vz*p.z)-s.k2*p.vz;
+    if (A == 0)
+        return -1;
+    double t = -B/A;
+    if (t <= 0 || p.vz*t+p.z > s.ze || p.vz*t+p.z < s.zs || p.vy*t+p.y < s.ll || p.vy*t+p.y > s.rl)
+        return -1;
+    return t;
+}
+
 /*! \brief Function to rotate a vector (by D. Liu)
 
 @param v vector to rotate
@@ -1684,10 +1725,10 @@ double reflectNeutronFlat(_class_particle* _particle, FlatSurf s) {
 
 @param p Pointer of particle to reflect
 @param c ConicSurf to use
+@param t Time until the collision, from traceSingleNeutron()
 
 */
-void traceNeutronConic(_class_particle* _particle, ConicSurf c) {
-    double t = getTimeOfFirstCollisionConic(*_particle, c);
+void traceNeutronConic(_class_particle* _particle, ConicSurf c, double t) {
     if (t < 0)
         return;
     else {
@@ -1707,10 +1748,10 @@ void traceNeutronConic(_class_particle* _particle, ConicSurf c) {
 
 @param p Pointer of particle to reflect
 @param f FlatSurf to use
+@param t Time until the collision, from traceSingleNeutron()
 
 */
-void traceNeutronFlat(_class_particle* _particle, FlatSurf f) {
-    double t = getTimeOfFirstCollisionFlat(*_particle, f);
+void traceNeutronFlat(_class_particle* _particle, FlatSurf f, double t) {
     if (t < 0)
         return;
     else {
@@ -1772,14 +1813,18 @@ void initSimulation(Scene* s) {
 void traceSingleNeutron(_class_particle* _particle, Scene* s) {
 
     int contact = 1;
+    enum GEO last_type = NONE; // surface the particle is on after the previous step
+    int last_index = -1;
     do {
         double t;
         enum  GEO type = NONE;
         int index = -1;
         int i;
 
         for (i = 0; i < s->num_c; i++) {
-            double t2 = getTimeOfFirstCollisionConic(*_particle,s->c[i]);
+            double t2 = (last_type == CONIC && last_index == i) ?
+                getTimeOfNextCollisionConic(*_particle,s->c[i]) :
+                getTimeOfFirstCollisionConic(*_particle,s->c[i]);
 
             if (t2 <= 0)
                 continue;
@@ -1791,7 +1836,9 @@ void traceSingleNeutron(_class_particle* _particle, Scene* s) {
         }
 
         for (i = 0; i < s->num_f; i++) {
-            double t2 = getTimeOfFirstCollisionFlat(*_particle,s->f[i]);
+            double t2 = (last_type == FLAT && last_index == i) ?
+                getTimeOfNextCollisionFlat(*_particle,s->f[i]) :
+                getTimeOfFirstCollisionFlat(*_particle,s->f[i]);
 
             if (t2 <= 0)
                 continue;
@@ -1831,18 +1878,20 @@ void traceSingleNeutron(_class_particle* _particle, Scene* s) {
                 traceNeutronDetector(_particle, s->d[index]);
                 break;
             case FLAT:
-	        traceNeutronFlat(_particle, s->f[index]);
+	        traceNeutronFlat(_particle, s->f[index], t);
                 break;
             case DISK:
                 traceNeutronDisk(_particle, s->di[index]);
                 break;
             case CONIC:
-                traceNeutronConic(_particle, s->c[index]);
+                traceNeutronConic(_particle, s->c[index], t);
                 break;
             default:
                 contact = 0;
                 break;
         }
+        last_type = type;
+        last_index = index;
     } while (contact && !_particle->_absorbed);
 
 }