diff --git a/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr b/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr index fe92c5028..d14a24a44 100644 --- a/mcstas-comps/examples/Tests_optics/Test_Conics_pairs/Test_Conics_pairs.instr +++ b/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 diff --git a/mcstas-comps/share/conic.h b/mcstas-comps/share/conic.h index ce0543031..e84b00a9d 100755 --- a/mcstas-comps/share/conic.h +++ b/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,6 +1813,8 @@ 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; @@ -1779,7 +1822,9 @@ void traceSingleNeutron(_class_particle* _particle, Scene* s) { 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); }