Fixed propagation procedure

build.gradle

@@ -17,6 +17,10 @@ tasks.withType(JavaCompile) {
     options.encoding = 'UTF-8'
 }
 
+tasks.withType(JavaExec){
+    enableAssertions = true
+}
+
 repositories {
     mavenCentral()
 }

gradle/wrapper/gradle-wrapper.properties

@@ -1,6 +1,6 @@
-#Sun May 08 12:57:40 MSK 2016
+#Tue Jun 14 13:15:55 MSK 2016
 distributionBase=GRADLE_USER_HOME
 distributionPath=wrapper/dists
 zipStoreBase=GRADLE_USER_HOME
 zipStorePath=wrapper/dists
-distributionUrl=https\://services.gradle.org/distributions/gradle-2.13-bin.zip
+distributionUrl=https\://services.gradle.org/distributions/gradle-2.14-bin.zip

src/main/java/inr/numass/trapping/SimulationManager.java

@@ -108,6 +108,7 @@ public class SimulationManager {
         Counter counter = new Counter();
         Predicate<Simulator.SimulationResult> reportIf = (res) -> res.state == Simulator.EndState.ACCEPTED;
         System.out.printf("%nStarting sumulation with initial energy %g and %d electrons.%n%n", initialE, num);
+        output.printf("%s\t%s\t%s\t%s\t%s\t%s%n", "E", "theta", "theta_start", "colNum", "L", "state");
         Stream.generate(() -> getRandomTheta()).limit(num).parallel()
                 .forEach((theta) -> {
                     double initZ = (generator.nextDouble() - 0.5) * Simulator.SOURCE_LENGTH;
@@ -144,7 +145,7 @@ public class SimulationManager {
     }
 
     private void printOne(PrintStream out, Simulator.SimulationResult res) {
-        out.printf("%g\t%g\t%g\t%d\t%s%n", res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.state.toString());
+        out.printf("%g\t%g\t%g\t%d\t%g\t%s%n", res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.l, res.state.toString());
     }
 
 

src/main/java/inr/numass/trapping/Simulator.java

@@ -45,7 +45,7 @@ public class Simulator {
         setELow(Elow);
     }
 
-    public static enum EndState {
+    public enum EndState {
 
         ACCEPTED,//трэппинговый электрон попал в аксептанс
         REJECTED,//трэппинговый электрон вылетел через заднюю пробку
@@ -64,14 +64,30 @@ public class Simulator {
         this.thetaPinch = Math.asin(Math.sqrt(Bsource / Bpinch));
     }
 
+    /**
+     * Set gas density in 1/m^3
+     *
+     * @param gasDensity
+     */
     public void setGasDensity(double gasDensity) {
         this.gasDensity = gasDensity;
     }
 
+    /**
+     * Longitudal magnetic field distribution
+     *
+     * @param magneticField
+     */
     public void setFieldFunc(UnivariateFunction magneticField) {
         this.magneticField = magneticField;
     }
 
+    /**
+     * Perform scattering in the given position
+     *
+     * @param pos
+     * @return
+     */
     private State scatter(State pos) {
         //Вычисляем сечения и нормируем их на полное сечение
         double sigmaIon = scatter.sigmaion(pos.e);
@@ -116,13 +132,13 @@ public class Simulator {
      * @return
      */
     private double freePath(double e) {
-        //FIXME double cross-section calculation
+        //FIXME redundant cross-section calculation
         //All cross-sections are in m^2
-        return new ExponentialDistribution(generator, 1 / scatter.sigmaTotal(e) / gasDensity).sample();
+        return new ExponentialDistribution(generator, 1d / scatter.sigmaTotal(e) / gasDensity).sample();
     }
 
     /**
-     * Calculate propagated position
+     * Calculate propagated position before scattering
      *
      * @param deltaL
      * @return z shift and reflection counter
@@ -131,23 +147,27 @@ public class Simulator {
         // if magnetic field not defined, consider it to be uniform and equal bSource
         if (magneticField == null) {
             double deltaZ = deltaL * cos(pos.theta); // direction already included in cos(theta)
+            double z0 = pos.z;
+            pos.addZ(deltaZ);
+            pos.l += abs(pos.z - z0) / cos(pos.theta);
 
-            //if we are crossing source boundary, check for end condition
+//            //if we are crossing source boundary, check for end condition
-            if(abs(deltaZ + pos.z)>SOURCE_LENGTH/2d) {
+//            while (abs(deltaZ + pos.z) > SOURCE_LENGTH / 2d && !pos.isFinished()) {
-                checkEndState(pos);
+//
-            }
+//                pos.checkEndState();
-
+//            }
-            // if track is finished apply boundary position
+//
-            if (pos.isFinished()) {
+//            // if track is finished apply boundary position
-                // remembering old z to correctly calculate total l
+//            if (pos.isFinished()) {
-                double oldz = pos.z;
+//                // remembering old z to correctly calculate total l
-                pos.z = pos.direction() * SOURCE_LENGTH / 2d;
+//                double oldz = pos.z;
-                pos.l += (pos.z - oldz)/cos(pos.theta);
+//                pos.z = pos.direction() * SOURCE_LENGTH / 2d;
-            } else {
+//                pos.l += (pos.z - oldz) / cos(pos.theta);
-                //else just add z
+//            } else {
-                pos.l += deltaL;
+//                //else just add z
-                pos.addZ(deltaZ);
+//                pos.l += deltaL;
-            }
+//                pos.addZ(deltaZ);
+//            }
 
             return pos;
         } else {
@@ -161,31 +181,45 @@ public class Simulator {
                 double b = field(pos.z);
 
                 double root = 1 - sin2 * b / bSource;
-                // change direction in case of reflection. Loss of precision here?
+                //preliminary reflection
                 if (root < 0) {
-                    //flip direction
                     pos.flip();
-                    // check if end state occured
-                    checkEndState(pos);
-                    // finish if it does
-                    if(pos.isFinished()){
-                        return pos;
-                    } else {
-                        // move in reversed direction
-                        pos.z += pos.direction() * delta * sqrt(-root);
-                    }
-                } else {
-                    // move
-                    pos.z += pos.direction() * delta * sqrt(root);
                 }
+                pos.addZ(pos.direction() * delta * sqrt(abs(root)));
+//                // change direction in case of reflection. Loss of precision here?
+//                if (root < 0) {
+//                    // check if end state occurred. seem to never happen since it is reflection case
+//                    pos.checkEndState();
+//                    // finish if it does
+//                    if (pos.isFinished()) {
+//                        //TODO check if it happens
+//                        return pos;
+//                    } else {
+//                        //flip direction
+//                        pos.flip();
+//                        // move in reversed direction
+//                        pos.z += pos.direction() * delta * sqrt(-root);
+//                    }
+//
+//                } else {
+//                    // move forward
+//                    pos.z += pos.direction() * delta * sqrt(root);
+//                    //check if it is exit case
+//                    if (abs(pos.z) > SOURCE_LENGTH / 2d) {
+//                        // check if electron is out
+//                        pos.checkEndState();
+//                        // finish if it is
+//                        if (pos.isFinished()) {
+//                            return pos;
+//                        }
+//                        // PENDING no need to apply reflection, it is applied automatically when root < 0
+//                        pos.z = signum(pos.z) * SOURCE_LENGTH / 2d;
+//                        if (signum(pos.z) == pos.direction()) {
+//                            pos.flip();
+//                        }
+//                    }
+//                }
 
-                //normalize in case reflection is beyound the source
-                if (abs(pos.z) > SOURCE_LENGTH / 2d) {
-                    pos.z = signum(pos.z) * SOURCE_LENGTH / 2d;
-                    if(signum(pos.z)== pos.direction()){
-                        pos.flip();
-                    }
-                }
 
                 curL += delta;
                 pos.l += delta;
@@ -194,30 +228,6 @@ public class Simulator {
         }
     }
 
-    /**
-     * Check if this position is an end state and apply it if necessary
-     *
-     * @param pos
-     * @return
-     */
-    private State checkEndState(State pos) {
-        //accepted by spectrometer
-        if (pos.theta < thetaPinch) {
-            if (pos.colNum == 0) {
-                //counting pass electrons
-                pos.setEndState(EndState.PASS);
-            } else {
-                pos.setEndState(EndState.ACCEPTED);
-            }
-        }
-
-        //through the rear magnetic pinch
-        if (pos.theta >= PI - thetaTransport) {
-            pos.setEndState(EndState.REJECTED);
-        }
-        return pos;
-    }
-
     /**
      * Magnetic field in the z point
      *
@@ -246,11 +256,12 @@ public class Simulator {
         while (!pos.isFinished()) {
             double dl = freePath(pos.e); // path to next scattering
             // propagate to next scattering position
-            propagate(pos,dl);
+            propagate(pos, dl);
 
             if (!pos.isFinished()) {
                 // perform scatter
                 scatter(pos);
+                // increase collision number
                 pos.colNum++;
                 if (pos.e < eLow) {
                     //Если энергия стала слишком маленькой
@@ -363,20 +374,53 @@ public class Simulator {
          */
         double addZ(double dZ) {
             this.z += dZ;
-            while (abs(this.z) > SOURCE_LENGTH / 2d) {
+            while (abs(this.z) > SOURCE_LENGTH / 2d && !isFinished()) {
+                checkEndState();
+                if (!isFinished()) {
+                    flip();
+                }
                 // reflecting from back wall
                 if (z < 0) {
-                    // reflecting from rear pinch
+                    if (isFinished()) {
-                    z = -SOURCE_LENGTH - z;
+                        z = -SOURCE_LENGTH / 2d;
+                    } else {
+                        // reflecting from rear pinch
+                        z = -SOURCE_LENGTH - z;
+                    }
                 } else {
-                    // reflecting from forward transport magnet
+                    if (isFinished()) {
-                    z = SOURCE_LENGTH - z;
+                        z = SOURCE_LENGTH / 2d;
+                    } else {
+                        // reflecting from forward transport magnet
+                        z = SOURCE_LENGTH - z;
+                    }
                 }
-                flip();
             }
             return z;
         }
 
+        /**
+         * Check if this position is an end state and apply it if necessary. Does not check z position.
+         *
+         * @return
+         */
+        private void checkEndState() {
+            //accepted by spectrometer
+            if (theta < thetaPinch) {
+                if (colNum == 0) {
+                    //counting pass electrons
+                    setEndState(EndState.PASS);
+                } else {
+                    setEndState(EndState.ACCEPTED);
+                }
+            }
+
+            //through the rear magnetic pinch
+            if (theta >= PI - thetaTransport) {
+                setEndState(EndState.REJECTED);
+            }
+        }
+
         /**
          * Reverse electron direction
          */
@@ -406,9 +450,8 @@ public class Simulator {
                 if (theta > PI / 2) {
                     newTheta = PI - newTheta;
                 }
-                if (Double.isNaN(newTheta)) {
+
-                    throw new Error();
+                assert !Double.isNaN(newTheta);
-                }
                 return newTheta;
             }
         }
@@ -437,13 +480,13 @@ public class Simulator {
 
             double newtheta = acos(result.getZ());
 
-            //следим чтобы угол был от 0 до Pi
+//            //следим чтобы угол был от 0 до Pi
-            if (newtheta < 0) {
+//            if (newtheta < 0) {
-                newtheta = -newtheta;
+//                newtheta = -newtheta;
-            }
+//            }
-            if (newtheta > Math.PI) {
+//            if (newtheta > Math.PI) {
-                newtheta = 2 * Math.PI - newtheta;
+//                newtheta = 2 * Math.PI - newtheta;
-            }
+//            }
 
             //change back to virtual angles
             if (magneticField == null) {
@@ -455,9 +498,7 @@ public class Simulator {
                 }
             }
 
-            if (Double.isNaN(theta)) {
+            assert !Double.isNaN(theta);
-                throw new Error();
-            }
 
             return theta;
         }

src/main/java/inr/numass/trapping/Trapping.java

@@ -16,10 +16,10 @@ public class Trapping {
         System.out.printf("Starting at %s%n%n", startTime.toString());
         new SimulationManager()
                 .withParameters(0.6, 3.7, 4.84, 18000d, 4000)
-//                .withOutputFile("D:\\Work\\Numass\\trapping\\trap 18, pinch 100A, fields.out")
+                .withOutputFile("D:\\Work\\Numass\\trapping\\trap 18, pinch 100A.out")
 //                .withFieldMap(z, b)
 //                .withDensity(5e20)
-                .simulateAll((int) 1e4);
+                .simulateAll((int) 1e6);
         Instant finishTime = Instant.now();
         System.out.printf("%nFinished at %s%n", finishTime.toString());
         System.out.printf("Calculation took %s%n", Duration.between(startTime,finishTime).toString());