New iteration

build.gradle

@@ -18,7 +18,7 @@ tasks.withType(JavaCompile) {
 }
 
 repositories {
-    mavencentral()
+    mavenCentral()
 }
 
 dependencies {

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

@@ -7,7 +7,7 @@ import org.apache.commons.math3.random.JDKRandomGenerator;
 import org.apache.commons.math3.random.RandomGenerator;
 import org.apache.commons.math3.random.SynchronizedRandomGenerator;
 
-import java.io.PrintStream;
+import java.io.*;
 import java.util.Map;
 import java.util.function.Function;
 import java.util.function.Predicate;
@@ -31,6 +31,16 @@ public class SimulationManager {
         return this;
     }
 
+    public SimulationManager withOutputFile(String fileName) throws IOException {
+        File outputFile = new File(fileName);
+        if (!outputFile.exists()) {
+            outputFile.getParentFile().mkdirs();
+            outputFile.createNewFile();
+        }
+        this.output = new PrintStream(new FileOutputStream(outputFile));
+        return this;
+    }
+
     /**
      * Select output for accepted events
      *
@@ -100,7 +110,7 @@ public class SimulationManager {
         System.out.printf("%nStarting sumulation with initial energy %g and %d electrons.%n%n", initialE, num);
         Stream.generate(() -> getRandomTheta()).limit(num).parallel()
                 .forEach((theta) -> {
-                    double initZ = (generator.nextDouble()-1) * Simulator.SOURCE_LENGTH;
+                    double initZ = (generator.nextDouble() - 0.5) * Simulator.SOURCE_LENGTH;
                     Simulator.SimulationResult res = simulator.simulate(initialE, theta, initZ);
                     if (reportIf.test(res)) {
                         if (output != null) {

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

@@ -11,7 +11,6 @@ import org.apache.commons.math3.geometry.euclidean.threed.SphericalCoordinates;
 import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
 import org.apache.commons.math3.random.JDKRandomGenerator;
 import org.apache.commons.math3.random.RandomGenerator;
-import org.apache.commons.math3.util.DoubleArray;
 import org.apache.commons.math3.util.Pair;
 import org.apache.commons.math3.util.Precision;
 
@@ -23,7 +22,7 @@ import static org.apache.commons.math3.util.FastMath.*;
 public class Simulator {
 
     public static final double SOURCE_LENGTH = 3d;
-    private static final double DELTA_L = 0.1; //step for deltaZ calculation
+    private static final double DELTA_L = 0.1; //step for propagate calculation
 
     private RandomGenerator generator;
     private Scatter scatter;
@@ -123,64 +122,104 @@ public class Simulator {
     }
 
     /**
-     * @param position
+     * Calculate propagated position
-     * @param dl
-     * @return
-     */
-    private State propagate(State position, double dl, double dz) {
-        //increase l
-        position.l += dl;
-        position.addZ(dz);
-        return position;
-    }
-
-    /**
-     * calculate z coordinate change with known path length. Does not change position.
      *
      * @param deltaL
      * @return z shift and reflection counter
      */
-    private Pair<Double, Integer> deltaZ(State position, double deltaL) {
+    private State propagate(State pos, double deltaL) {
         // if magnetic field not defined, consider it to be uniform and equal bSource
-        int direction = position.isForward() ? 1 : -1;
         if (magneticField == null) {
-            double deltaZ = direction * deltaL * cos(position.theta);
+            double deltaZ = deltaL * cos(pos.theta); // direction already included in cos(theta)
-            int reflectionCounter = (int) (deltaZ / SOURCE_LENGTH);
+
-            double curZ = normalizeZ(position.z + deltaZ);
+            //if we are crossing source boundary, check for end condition
-            return new Pair<>(curZ - position.z, reflectionCounter);
+            if(abs(deltaZ + pos.z)>SOURCE_LENGTH/2d) {
+                checkEndState(pos);
+            }
+
+            // if track is finished apply boundary position
+            if (pos.isFinished()) {
+                // remembering old z to correctly calculate total l
+                double oldz = pos.z;
+                pos.z = pos.direction() * SOURCE_LENGTH / 2d;
+                pos.l += (pos.z - oldz)/cos(pos.theta);
+            } else {
+                //else just add z
+                pos.l += deltaL;
+                pos.addZ(deltaZ);
+            }
+
+            return pos;
         } else {
-            int reflect = 0;
-            double curZ = position.z;
             double curL = 0;
+            double sin2 = sin(pos.theta) * sin(pos.theta);
 
+            while (curL <= deltaL - 0.01 && !pos.isFinished()) {
+                //an l step
+                double delta = min(deltaL - curL, DELTA_L);
 
-            double sin2 = sin(position.theta) * sin(position.theta);
+                double b = field(pos.z);
 
-            while (curL <= deltaL) {
-                double delta = DELTA_L; //min(deltaL - curL, DELTA_L);
-                double b = field(curZ);
-
-                //TODO check this!
                 double root = 1 - sin2 * b / bSource;
                 // change direction in case of reflection. Loss of precision here?
                 if (root < 0) {
-                    direction = -direction;
+                    //flip direction
-                    root = -root;
+                    pos.flip();
-                    reflect++;
+                    // 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);
                 }
 
-                curZ += direction * delta * sqrt(root);
+                //normalize in case reflection is beyound the source
-
+                if (abs(pos.z) > SOURCE_LENGTH / 2d) {
-                curZ = normalizeZ(curZ);
+                    pos.z = signum(pos.z) * SOURCE_LENGTH / 2d;
+                    if(signum(pos.z)== pos.direction()){
+                        pos.flip();
+                    }
+                }
 
                 curL += delta;
+                pos.l += delta;
             }
-            return new Pair<>(curZ - position.z, reflect);
+            return pos;
         }
     }
 
     /**
-     * Magnetic field with reflection taken into account
+     * 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
      *
      * @param z
      * @return
@@ -193,19 +232,6 @@ public class Simulator {
         }
     }
 
-    private double normalizeZ(double z) {
-        while ((abs(z) > SOURCE_LENGTH / 2)) {
-            if (z < 0) {
-                // reflecting from rear pinch
-                z += SOURCE_LENGTH / 2d;
-            } else {
-                // reflecting from forward transport magnet
-                z -= SOURCE_LENGTH / 2d;
-            }
-        }
-        return z;
-    }
-
     /**
      * Симулируем один пробег электрона от начального значения и до вылетания из
      * иточника или до того момента, как энергия становится меньше eLow.
@@ -216,52 +242,24 @@ public class Simulator {
         assert abs(initZ) <= SOURCE_LENGTH / 2d;
 
         State pos = new State(initEnergy, initTheta, initZ);
-        EndState endState = EndState.NONE;
-
-        boolean stopflag = false;
-
-        while (!stopflag) {
 
+        while (!pos.isFinished()) {
             double dl = freePath(pos.e); // path to next scattering
-            Pair<Double, Integer> dzCalc = deltaZ(pos, dl);
+            // propagate to next scattering position
-            double dz = dzCalc.getFirst();// z coordinate to next scattering
+            propagate(pos,dl);
-            int reflections = dzCalc.getSecond();
 
-            //if no scattering on current source pass
+            if (!pos.isFinished()) {
-            if (reflections > 0) {
-                //accepted by spectrometer
-                if (pos.theta < thetaPinch) {
-                    stopflag = true;
-                    if (pos.colNum == 0) {
-                        //counting pass electrons
-                        endState = EndState.PASS;
-                    } else {
-                        endState = EndState.ACCEPTED;
-                    }
-                }
-
-
-                //through the rear magnetic pinch
-                if (pos.theta >= PI - thetaTransport) {
-                    stopflag = true;
-                    endState = EndState.REJECTED;
-                }
-            }
-            if (!stopflag) {
                 // perform scatter
-                propagate(pos, dl, dz);
-                pos.colNum++;
                 scatter(pos);
+                pos.colNum++;
                 if (pos.e < eLow) {
                     //Если энергия стала слишком маленькой
-                    stopflag = true;
+                    pos.setEndState(EndState.LOWENERGY);
-                    endState = EndState.LOWENERGY;
                 }
             }
-
         }
 
-        SimulationResult res = new SimulationResult(endState, pos.e, pos.theta, initTheta, pos.colNum, pos.l);
+        SimulationResult res = new SimulationResult(pos.endState, pos.e, pos.theta, initTheta, pos.colNum, pos.l);
         return res;
     }
 
@@ -330,6 +328,8 @@ public class Simulator {
          */
         double z;
 
+        EndState endState = EndState.NONE;
+
         /**
          * @param dE
          * @return resulting E
@@ -343,6 +343,18 @@ public class Simulator {
             return theta <= Math.PI / 2;
         }
 
+        double direction() {
+            return isForward() ? 1 : -1;
+        }
+
+        boolean isFinished() {
+            return this.endState != EndState.NONE;
+        }
+
+        void setEndState(EndState state) {
+            this.endState = state;
+        }
+
         /**
          * add Z and calculate direction change
          *
@@ -354,12 +366,13 @@ public class Simulator {
             while (abs(this.z) > SOURCE_LENGTH / 2d) {
                 // reflecting from back wall
                 if (z < 0) {
-                    z += SOURCE_LENGTH / 2d;
+                    // reflecting from rear pinch
-                    flip();
+                    z = -SOURCE_LENGTH - z;
                 } else {
-                    z -= SOURCE_LENGTH / 2d;
+                    // reflecting from forward transport magnet
-                    flip();
+                    z = SOURCE_LENGTH - z;
                 }
+                flip();
             }
             return z;
         }
@@ -400,6 +413,7 @@ public class Simulator {
             }
         }
 
+
         /**
          * Сложение вектора с учетом случайно распределения по фи
          *

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

@@ -8,16 +8,7 @@ import java.time.Instant;
 public class Trapping {
     public static void main(String[] args) throws IOException {
 //        new SimulationManager().withParameters(0.6, 3.7, 7.2, 18000d, 4000).simulateAll((int) 1e6);
-
+        double[] z = {-1.736, -1.27, -0.754, -0.238, 0.278, 0.794, 1.31, 1.776};
-//        -0.236	3.70754
-//        0.23	0.62786
-//        0.746	0.60474
-//        1.262	0.60325
-//        1.778	0.60333
-//        2.294	0.60503
-//        2.81	0.6285
-//        3.276	3.70478
-        double[] z = {-0.236, 0.23, 0.746, 1.262, 1.778, 2.294, 2.81, 3.276};
         double[] b = {3.70754, 0.62786, 0.60474, 0.60325, 0.60333, 0.60503, 0.6285, 3.70478};
         System.out.println("Press any key to start...");
         System.in.read();
@@ -25,9 +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")
 //                .withFieldMap(z, b)
 //                .withDensity(5e20)
-                .simulateAll((int) 1e5);
+                .simulateAll((int) 1e4);
         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());