Direction, density effect and pinch scattering

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

@@ -8,15 +8,7 @@
  */
 package inr.numass.trapping;
 
-import static java.lang.Math.abs;
+import static  org.apache.commons.math3.util.FastMath.*;
-import static java.lang.Math.acos;
-import static java.lang.Math.atan;
-import static java.lang.Math.cos;
-import static java.lang.Math.exp;
-import static java.lang.Math.log;
-import static java.lang.Math.sin;
-import static java.lang.Math.sqrt;
-import static java.lang.Math.tan;
 import org.apache.commons.math3.random.RandomGenerator;
 import org.apache.commons.math3.util.Pair;
 
@@ -31,6 +23,9 @@ public class Scatter {
     static final double emass = 18780; // Electron mass in atomic units
     static final double R = 13.6; // Ryberg energy in eV
     private final RandomGenerator generator;
+
+    boolean debug = false;
+
     MultiCounter counter = new MultiCounter("Accept-reject calls");
 
     private double fmax = 0;
@@ -39,6 +34,12 @@ public class Scatter {
         this.generator = generator;
     }
 
+    private void count(String counterName){
+        if(debug){
+            counter.increase(counterName);
+        }
+    }
+
     public Pair<Double, Double> randomel(double E) {
         // This subroutine generates  energy loss and polar scatt. angle according to
         // electron elastic scattering in molecular hydrogen.
@@ -52,7 +53,7 @@ public class Scatter {
         double[] u = new double[3];
         int i;
 
-        counter.increase("randomel-calls");
+        count("randomel-calls");
         if (E >= 250.) {
             Gmax = 1.e-19;
         } else if (E < 250. && E >= 150.) {
@@ -64,7 +65,7 @@ public class Scatter {
         gam = 1. + T / (clight * clight); // relativistic correction factor
         b = 2. / (1. + gam) / T;
         for (i = 1; i < 5000; i++) {
-            counter.increase("randomel");
+            count("randomel");
             c = 1. + b - b * (2. + b) / (b + 2. * generator.nextDouble());
             K2 = 2. * T * (1. + gam) * abs(1d - c); // momentum transfer squared
             a = (4. + K2) * (4. + K2) / (gam * gam);
@@ -117,7 +118,7 @@ public class Scatter {
         //   (from: Phys. Rev. A51 (1995) 3745 )
         double[] pC = {1.2e-1, 1.9e-1, 1.9e-1, 1.5e-1, 1.1e-1, 7.5e-2, 5.0e-2,
             3.3e-2, 2.2e-2, 1.4e-2, 9.3e-3, 6.0e-3, 3.7e-3, 1.8e-3};
-        counter.increase("randomexc-calls");
+        count("randomexc-calls");
         T = 20000. / 27.2;
         //
         xmin = Ecen * Ecen / (2. * T);
@@ -153,7 +154,7 @@ public class Scatter {
             // Generation of y values with the Neumann acceptance-rejection method:
             y = ymin;
             for (j = 1; j < 5000; j++) {
-                counter.increase("randomexc1");
+                count("randomexc1");
                 y = ymin + (ymax - ymin) * generator.nextDouble();
                 K = exp(y / 2.);
                 fy = sumexc(K);
@@ -176,7 +177,7 @@ public class Scatter {
                 Dmax = 400.;
             }
             for (j = 1; j < 5000; j++) {
-                counter.increase("randomexc2");
+                count("randomexc2");
                 c = -1. + 2. * generator.nextDouble();
                 D = Dexc(E, c) * 1.e22;
                 if (Dmax * generator.nextDouble() < D) {
@@ -192,7 +193,7 @@ public class Scatter {
         N = 7; // the number of electronic states in our calculation
         pmax = p[1]; // the maximum of the p[] values
         for (j = 1; j < 5000; j++) {
-            counter.increase("randomexc3");
+            count("randomexc3");
             n = (int) (N * generator.nextDouble());
             if (generator.nextDouble() * pmax < p[n]) {
                 break;
@@ -215,7 +216,7 @@ public class Scatter {
                 N = 28; // the number of B vibrational states in our calculation
                 pmax = pB[7]; // maximum of the pB[] values
                 for (j = 1; j < 5000; j++) {
-                    counter.increase("randomexc4");
+                    count("randomexc4");
                     v = (int) (N * generator.nextDouble());
                     if (generator.nextDouble() * pmax < pB[v]) {
                         break;
@@ -235,7 +236,7 @@ public class Scatter {
                 N = 14; // the number of C vibrational states in our calculation
                 pmax = pC[1]; // maximum of the pC[] values
                 for (j = 1; j < 5000; j++) {
-                    counter.increase("randomexc4");
+                    count("randomexc4");
                     v = (int) (N * generator.nextDouble());
                     if (generator.nextDouble() * pmax < pC[v]) {
                         break;
@@ -275,7 +276,7 @@ public class Scatter {
         double K2, KK, fE, kej, ki, kf, Rex, arg, arctg;
         int i;
         double st1, st2;
-        counter.increase("randomion-calls");
+        count("randomion-calls");
         //
         // I. Generation of theta
         // -----------------------
@@ -291,7 +292,7 @@ public class Scatter {
         // Generation of y values with the Neumann acceptance-rejection method:
         y = ymin;
         for (j = 1; j < 5000; j++) {
-            counter.increase("randomion1");
+            count("randomion1");
             y = ymin + (ymax - ymin) * generator.nextDouble();
             K = exp(y / 2.);
             c = 1. + b - K * K / (4. * T);
@@ -350,7 +351,7 @@ public class Scatter {
             // Generation of q with inverse transform method
             // (we use the Newton-Raphson method in order to solve the nonlinear eq.
             // for the inversion) :
-            counter.increase("randomion2");
+            count("randomion2");
             F = Fmin + h * generator.nextDouble();
             y = 0.;
             for (i = 1; i <= 30; i++) {

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

@@ -1,9 +1,15 @@
 package inr.numass.trapping;
 
+import org.apache.commons.math3.analysis.UnivariateFunction;
+import org.apache.commons.math3.analysis.interpolation.LinearInterpolator;
+import org.apache.commons.math3.random.ISAACRandom;
 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.util.Map;
+import java.util.function.Function;
 import java.util.function.Predicate;
 import java.util.stream.Stream;
 
@@ -35,6 +41,21 @@ public class SimulationManager {
         return this;
     }
 
+    public SimulationManager withFieldMap(UnivariateFunction fieldMap) {
+        this.simulator.setFieldFunc(fieldMap);
+        return this;
+    }
+
+    public SimulationManager withFieldMap(double[] z, double[] b) {
+        this.simulator.setFieldFunc(new LinearInterpolator().interpolate(z, b));
+        return this;
+    }
+
+    public SimulationManager withDensity(double density){
+        this.simulator.setGasDensity(density);
+        return this;
+    }
+
     /**
      * Симулируем пролет num электронов.
      *
@@ -47,7 +68,8 @@ 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) -> {
-                    Simulator.SimulationResult res = simulator.simulate(initialE, theta);
+                    double initZ = generator.nextDouble() * Simulator.SOURCE_LENGTH;
+                    Simulator.SimulationResult res = simulator.simulate(initialE, theta, initZ);
                     if (reportIf.test(res)) {
                         if (output != null) {
                             printOne(output, res);
@@ -61,7 +83,8 @@ public class SimulationManager {
 
     private double getRandomTheta() {
         double x = generator.nextDouble();
-        return Math.acos(x);
+        // from 0 to 2 Pi
+        return Math.acos(1 - 2 * x);
     }
 
     private void printStatistics(Counter counter) {
@@ -70,7 +93,7 @@ public class SimulationManager {
         output.printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI);
         output.printf("The transport reflection angle is %g.%n", simulator.thetaTransport * 180 / Math.PI);
         output.printf("The starting energy is %g.%n", initialE);
-        output.printf("The lower energy boundary is %g.%n%n", simulator.Elow);
+        output.printf("The lower energy boundary is %g.%n%n", simulator.eLow);
 
         output.printf("The total number of ACCEPTED events is %d.%n", counter.accepted);
         output.printf("The total number of PASS events is %d.%n", counter.pass);

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

@@ -4,40 +4,47 @@
  */
 package inr.numass.trapping;
 
-import java.io.PrintStream;
+import org.apache.commons.math3.analysis.UnivariateFunction;
-import java.util.function.Predicate;
+import org.apache.commons.math3.distribution.ExponentialDistribution;
-import java.util.stream.Stream;
-
 import org.apache.commons.math3.geometry.euclidean.threed.Rotation;
 import org.apache.commons.math3.geometry.euclidean.threed.SphericalCoordinates;
 import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
-import org.apache.commons.math3.random.MersenneTwister;
+import org.apache.commons.math3.random.JDKRandomGenerator;
 import org.apache.commons.math3.random.RandomGenerator;
-import org.apache.commons.math3.random.SynchronizedRandomGenerator;
 import org.apache.commons.math3.util.Pair;
 import org.apache.commons.math3.util.Precision;
 
+import java.util.function.Function;
+
+import static org.apache.commons.math3.util.FastMath.*;
+
 /**
  * @author Darksnake
  */
 public class Simulator {
 
+    public static final double SOURCE_LENGTH = 3d;
+    private static final double DELTA_L = 0.1; //step for dZ calculation
+
     private RandomGenerator generator;
     private Scatter scatter;
-    double Elow = 14000d;//default value
+    double eLow = 14000d;//default value
     double thetaTransport = 24.107064 / 180 * Math.PI;// default value
     double thetaPinch = 19.481097 / 180 * Math.PI;// default value
+    double gasDensity = 5e18;// m^-3
+    double bSource = 0.6;
+    private UnivariateFunction magneticField;
 
 
     public Simulator() {
-        generator = new SynchronizedRandomGenerator(new MersenneTwister());
+        generator = new JDKRandomGenerator();// Be careful here since it is not synchronized
         scatter = new Scatter(generator);
     }
 
     public Simulator(double Bsource, double Btransport, double Bpinch, double Elow) {
         this();
         setFields(Bsource, Btransport, Bpinch);
-        setElow(Elow);
+        setELow(Elow);
     }
 
     public static enum EndState {
@@ -49,163 +56,223 @@ public class Simulator {
         NONE
     }
 
-    public final void setElow(double elow) {
+    public final void setELow(double eLow) {
-        Elow = elow;
+        this.eLow = eLow;
     }
 
     public final void setFields(double Bsource, double Btransport, double Bpinch) {
+        this.bSource = Bsource;
         this.thetaTransport = Math.asin(Math.sqrt(Bsource / Btransport));
         this.thetaPinch = Math.asin(Math.sqrt(Bsource / Bpinch));
     }
 
+    public void setGasDensity(double gasDensity) {
+        this.gasDensity = gasDensity;
+    }
+
+    public void setFieldFunc(UnivariateFunction magneticField) {
+        this.magneticField = magneticField;
+    }
+
+    private State scatter(State pos) {
+        //Вычисляем сечения и нормируем их на полное сечение
+        double sigmaIon = scatter.sigmaion(pos.e);
+        double sigmaEl = scatter.sigmael(pos.e);
+        double sigmaExc = scatter.sigmaexc(pos.e);
+        double sigmaTotal = sigmaEl + sigmaIon + sigmaExc;
+        sigmaIon /= sigmaTotal;
+        sigmaEl /= sigmaTotal;
+        sigmaExc /= sigmaTotal;
+
+        //проверяем нормировку
+        assert Precision.equals(sigmaEl + sigmaExc + sigmaIon, 1, 1e-2);
+
+        double alpha = generator.nextDouble();
+
+        Pair<Double, Double> delta;
+        if (alpha > sigmaEl) {
+            if (alpha > sigmaEl + sigmaExc) {
+                //ionization case
+                delta = scatter.randomion(pos.e);
+            } else {
+                //excitation case
+                delta = scatter.randomexc(pos.e);
+            }
+        } else {
+            // elastic
+            delta = scatter.randomel(pos.e);
+        }
+
+        //Обновляем значени угла и энергии независимо ни от чего
+        pos.substractE(delta.getFirst());
+        //Изменение угла
+        pos.addTheta(delta.getSecond() / 180 * Math.PI);
+
+        return pos;
+    }
+
+    /**
+     * Calculate distance to the next scattering
+     *
+     * @param e
+     * @return
+     */
+    private double freePath(double e) {
+        //FIXME double cross-section calculation
+        //All cross-sections are in m^2
+        return new ExponentialDistribution(generator, 1 / scatter.sigmaTotal(e) / gasDensity).sample();
+    }
+
+    /**
+     * @param 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 dl
+     * @return
+     */
+    private double dZ(State position, double dl) {
+        // if magnetic field not defined, consider it to be uniform and equal bSource
+        if (magneticField == null) {
+            return dl / cos(position.theta);
+        } else {
+            double dz = 0;
+            double curL = 0;
+            while (curL <= dl) {
+                double delta = min(dl - curL, DELTA_L);
+                double b = field(position.z + dz);
+
+                dz += delta / sqrt(1 - pow(sin(position.theta), 2) * b / bSource);
+                curL += DELTA_L;
+            }
+            return dz;
+        }
+    }
+
+    /**
+     * Magnetic field with reflection taken into account
+     *
+     * @param z
+     * @return
+     */
+    private double field(double z) {
+        if (magneticField == null) {
+            return bSource;
+        } else {
+            while (!(z > 0 && z < SOURCE_LENGTH)) {
+                // reflecting from back wall
+                if (z < 0) {
+                    z = -z;
+                } else {
+                    z -= SOURCE_LENGTH;
+                }
+            }
+            return magneticField.value(z);
+        }
+    }
+
     /**
      * Симулируем один пробег электрона от начального значения и до вылетания из
-     * иточника или до того момента, как энергия становится меньше Elow.
+     * иточника или до того момента, как энергия становится меньше eLow.
-     * Считаем, что за один проход источника происходит в среднем существенно
-     * меньше одного столкновения, поэтому выход вперед или назад совершенно
-     * симментричны.
      */
-    public SimulationResult simulate(double initEnergy, double initTheta) {
+    public SimulationResult simulate(double initEnergy, double initTheta, double initZ) {
         assert initEnergy > 0;
-        assert initTheta > 0 && initTheta < Math.PI / 2;
+        assert initTheta > 0 && initTheta < Math.PI;
+        assert initZ >= 0 && initZ < SOURCE_LENGTH;
 
-        if (initTheta < this.thetaPinch) {
+//        if (initTheta < this.thetaPinch) {
-            if (generator.nextBoolean()) {
+//            if (generator.nextBoolean()) {
-                return new SimulationResult(EndState.PASS, initEnergy, initTheta, initTheta, 0);
+//                return new SimulationResult(EndState.PASS, initEnergy, initTheta, initTheta, 0);
-            } else {
+//            } else {
-                return new SimulationResult(EndState.REJECTED, initEnergy, initTheta, initTheta, 0);
+//                return new SimulationResult(EndState.REJECTED, initEnergy, initTheta, initTheta, 0);
-            }
+//            }
-        } else if (initTheta < this.thetaTransport) {
+//        } else if (initTheta < this.thetaTransport) {
-            return new SimulationResult(EndState.REJECTED, initEnergy, initTheta, initTheta, 0);
+//            return new SimulationResult(EndState.REJECTED, initEnergy, initTheta, initTheta, 0);
-        }
+//        }
 
-        double E = initEnergy;
+        State pos = new State(initEnergy, initTheta, initZ);
-        double theta = initTheta;
+        EndState endState = EndState.NONE;
-        int colNum = 0;
-        EndState state = EndState.NONE;
 
         boolean stopflag = false;
 
         while (!stopflag) {
-            colNum++;
-            Pair<Double, Double> delta;
 
-            //Вычисляем сечения и нормируем их на полное сечение
+            double dl = freePath(pos.e); // path to next scattering
-            double sigmaIon = scatter.sigmaion(E);
+            double dz = dZ(pos, dl); // z coordinate to next scattering
-            double sigmaEl = scatter.sigmael(E);
+            double expectedZ = pos.z + dz; // expected z position of next scattering
-            double sigmaexc = scatter.sigmaexc(E);
-            double sigmaTotal = sigmaEl + sigmaIon + sigmaexc;
-            sigmaIon /= sigmaTotal;
-            sigmaEl /= sigmaTotal;
-            sigmaexc /= sigmaTotal;
-            //проверяем нормировку
-            assert Precision.equals(sigmaEl + sigmaexc + sigmaIon, 1, 1e-2);
 
-            double alpha = generator.nextDouble();
+            //if no scattering on current source pass
-
+            if (expectedZ < 0 || expectedZ > SOURCE_LENGTH) {
-            if (alpha > sigmaEl) {
+                //accepted by spectrometer
-                if (alpha > sigmaEl + sigmaexc) {
+                if (pos.theta < thetaPinch) {
-                    //ionization case
-                    delta = scatter.randomion(E);
-                } else {
-                    //excitation case
-                    delta = scatter.randomexc(E);
-                }
-            } else {
-                // elastic
-                delta = scatter.randomel(E);
-            }
-
-            //Обновляем значени угла и энергии независимо ни от чего
-            E -= delta.getFirst();
-            //Изменение угла
-            theta = addTheta(theta, delta.getSecond() / 180 * Math.PI);
-            //следим чтобы угол был от 0 до 90, если он перекинется через границу, считаем что электрон остается в потоке
-            theta = normalizeTheta(theta);
-
-            if (theta < thetaPinch) {
                     stopflag = true;
-                if (generator.nextBoolean()) {
                     //Учитываем тот факт, что электрон мог вылететь в правильный угол, но назад
-                    state = EndState.ACCEPTED;
+                    if (pos.colNum == 0) {
+                        //counting pass electrons
+                        endState = EndState.PASS;
                     } else {
-                    state = EndState.REJECTED;
+                        endState = EndState.ACCEPTED;
                     }
                 }
 
-            if (theta >= thetaPinch && theta <= thetaTransport) {
+
-                //Вылет через заднюю пробку эквивалентен отражению от передней.
+                //through the rear magnetic trap
-                //это верно при низких концентрациях
+                if (pos.theta >= PI - thetaTransport) {
                     stopflag = true;
-                state = EndState.REJECTED;
+                    endState = EndState.REJECTED;
                 }
 
-            if (E < Elow) {
+                if (pos.e < eLow) {
                     //Если энергия стала слишком маленькой
                     stopflag = true;
-                state = EndState.LOWENERGY;
+                    endState = EndState.LOWENERGY;
                 }
             }
+            if (!stopflag) {
+                // perform scatter
+                propagate(pos, dl, dz);
+                pos.colNum++;
+                scatter(pos);
+            }
 
-        SimulationResult res = new SimulationResult(state, E, theta, initTheta, colNum);
+        }
+
+        SimulationResult res = new SimulationResult(endState, pos.e, pos.theta, initTheta, pos.colNum, pos.l);
         return res;
     }
 
-    /**
-     * Сложение вектора с учетом случайно распределения по фи
-     *
-     * @param theta
-     * @param dTheta
-     * @return
-     */
-    private double addTheta(double theta, double dTheta) {
-        //Генерируем случайный фи
-        double phi = generator.nextDouble() * 2 * Math.PI;
-        //Создаем начальный вектор в сферических координатах
-        SphericalCoordinates init = new SphericalCoordinates(1, 0, theta + dTheta);
-        // Задаем вращение относительно оси, перпендикулярной исходному вектору
-        SphericalCoordinates rotate = new SphericalCoordinates(1, 0, theta);
-        // поворачиваем исходный вектор на dTheta
-        Rotation rot = new Rotation(rotate.getCartesian(), phi, null);
-
-        Vector3D result = rot.applyTo(init.getCartesian());
-
-        //      assert Vector3D.angle(result, rotate.getCartesian()) == dTheta;
-        return Math.acos(result.getZ());
-    }
-
-    private double normalizeTheta(double theta) {
-        double res = theta;
-        if (theta < 0) {
-            res = -theta;
-        }
-        if (res >= 0 && res <= Math.PI / 2) {
-            return res;
-        } else if (res > Math.PI / 2 && res <= Math.PI) {
-            return Math.PI - res;
-        } else if (res > Math.PI && res <= Math.PI * 3 / 2) {
-            return res - Math.PI;
-        } else {
-            throw new IllegalStateException();
-        }
-    }
-
     public void resetDebugCounters() {
+        scatter.debug = true;
         scatter.counter.resetAll();
     }
 
     public void printDebugCounters() {
+        if (scatter.debug) {
             scatter.counter.print(System.out);
+        } else {
+            throw new RuntimeException("Debug not initiated");
+        }
     }
 
     public static class SimulationResult {
 
-        public SimulationResult(EndState state, double E, double theta, double initTheta, int collisionNumber) {
+        public SimulationResult(EndState state, double E, double theta, double initTheta, int collisionNumber, double l) {
             this.state = state;
             this.E = E;
             this.theta = theta;
             this.initTheta = initTheta;
             this.collisionNumber = collisionNumber;
+            this.l = l;
         }
 
         public EndState state;
@@ -213,6 +280,145 @@ public class Simulator {
         public double initTheta;
         public double theta;
         public int collisionNumber;
+        public double l;
+    }
+
+    /**
+     * Current electron position in simulation. Not thread safe!
+     */
+    private class State {
+        public State(double e, double theta, double z) {
+            this.e = e;
+            this.theta = theta;
+            this.z = z;
+        }
+
+        /**
+         * Current energy
+         */
+        double e;
+        /**
+         * Current theta recalculated to the field in the center of the source
+         */
+        double theta;
+        /**
+         * Current total path
+         */
+        double l = 0;
+
+        /**
+         * Number of scatterings
+         */
+        int colNum = 0;
+
+        /**
+         * current z;
+         */
+        double z;
+
+        /**
+         * @param dE
+         * @return resulting E
+         */
+        double substractE(double dE) {
+            this.e -= dE;
+            return e;
+        }
+
+        boolean isForward() {
+            return theta <= Math.PI / 2;
+        }
+
+        /**
+         * add Z and calculate direction change
+         *
+         * @param dZ
+         * @return
+         */
+        double addZ(double dZ) {
+            this.z += dZ;
+            while (!(this.z > 0 && this.z < SOURCE_LENGTH)) {
+                // reflecting from back wall
+                if (z < 0) {
+                    z = -z;
+                    flip();
+                } else {
+                    z -= SOURCE_LENGTH;
+                    flip();
+                }
+            }
+            return z;
+        }
+
+        /**
+         * Reverse electron direction
+         */
+        void flip() {
+            theta = PI - theta;
+        }
+
+        /**
+         * Magnetic field in the current point
+         *
+         * @return
+         */
+        double field() {
+            return Simulator.this.field(z);
+        }
+
+        /**
+         * Real theta angle in current point
+         *
+         * @return
+         */
+        double realTheta() {
+            if (magneticField == null) {
+                return theta;
+            } else {
+                return asin(sin(theta) * sqrt(field() / bSource));
+            }
+        }
+
+        /**
+         * Сложение вектора с учетом случайно распределения по фи
+         *
+         * @param dTheta
+         * @return resulting angle
+         */
+        double addTheta(double dTheta) {
+            //Генерируем случайный фи
+            double phi = generator.nextDouble() * 2 * Math.PI;
+
+            //change to real angles
+            double realTheta = realTheta();
+            //Создаем начальный вектор в сферических координатах
+            SphericalCoordinates init = new SphericalCoordinates(1, 0, realTheta + dTheta);
+            // Задаем вращение относительно оси, перпендикулярной исходному вектору
+            SphericalCoordinates rotate = new SphericalCoordinates(1, 0, realTheta);
+            // поворачиваем исходный вектор на dTheta
+            Rotation rot = new Rotation(rotate.getCartesian(), phi, null);
+
+            Vector3D result = rot.applyTo(init.getCartesian());
+
+            double newtheta = acos(result.getZ());
+
+            //следим чтобы угол был от 0 до Pi
+            if (newtheta < 0) {
+                newtheta = -newtheta;
+            }
+            if (newtheta > Math.PI && newtheta <= Math.PI * 3 / 2) {
+                newtheta = 2 * Math.PI - newtheta;
+            }
+
+            //change back to virtual angles
+            if (magneticField == null) {
+                theta = newtheta;
+            } else {
+                theta = asin(sin(newtheta) * sqrt(bSource / field()));
+            }
+            return theta;
+        }
+
     }
 
 

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

@@ -1,13 +1,34 @@
 package inr.numass.trapping;
 
 import java.io.FileNotFoundException;
+import java.io.IOException;
+import java.time.Duration;
+import java.time.Instant;
 
-/**
- * Hello world!
- */
 public class Trapping {
-    public static void main(String[] args) throws FileNotFoundException {
+    public static void main(String[] args) throws IOException {
 //        new SimulationManager().withParameters(0.6, 3.7, 7.2, 18000d, 4000).simulateAll((int) 1e6);
-        new SimulationManager().withParameters(0.6, 3.7, 4.84, 14000d, 4000).simulateAll((int) 1e6);
+
+//        -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();
+        Instant startTime = Instant.now();
+        System.out.printf("Starting at %s%n%n", startTime.toString());
+        new SimulationManager()
+                .withParameters(0.6, 3.7, 4.84, 18000d, 4000)
+//                .withFieldMap(z, b)
+                .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());
     }
 }