diff --git a/build.gradle b/build.gradle index 129bee4..b223762 100644 --- a/build.gradle +++ b/build.gradle @@ -8,7 +8,7 @@ mainClassName = "inr.numass.trapping.Trapping" group = 'inr.numass' version = 'dev' -description = """trapping""" +description = "Numass trapping simulation" sourceCompatibility = 1.8 targetCompatibility = 1.8 @@ -18,9 +18,9 @@ tasks.withType(JavaCompile) { } repositories { - - maven { url "http://repo.maven.apache.org/maven2" } + mavencentral() } + dependencies { compile group: 'org.apache.commons', name: 'commons-math3', version:'3.6.1' testCompile group: 'junit', name: 'junit', version:'4.12' diff --git a/src/main/resources/csource/constants.h b/docs/csource/constants.h similarity index 100% rename from src/main/resources/csource/constants.h rename to docs/csource/constants.h diff --git a/src/main/resources/csource/random.c b/docs/csource/random.c similarity index 100% rename from src/main/resources/csource/random.c rename to docs/csource/random.c diff --git a/src/main/resources/csource/random.h b/docs/csource/random.h similarity index 100% rename from src/main/resources/csource/random.h rename to docs/csource/random.h diff --git a/docs/csource/readme.txt b/docs/csource/readme.txt new file mode 100644 index 0000000..88229bc --- /dev/null +++ b/docs/csource/readme.txt @@ -0,0 +1 @@ +Original C sources provided by Sebastian Voecking and Ferenc Glueck \ No newline at end of file diff --git a/src/main/resources/csource/scatter.c b/docs/csource/scatter.c similarity index 100% rename from src/main/resources/csource/scatter.c rename to docs/csource/scatter.c diff --git a/src/main/java/inr/numass/trapping/SimulationManager.java b/src/main/java/inr/numass/trapping/SimulationManager.java index 9eb1484..279434f 100644 --- a/src/main/java/inr/numass/trapping/SimulationManager.java +++ b/src/main/java/inr/numass/trapping/SimulationManager.java @@ -31,27 +31,59 @@ public class SimulationManager { return this; } + /** + * Select output for accepted events + * + * @param output + * @return + */ public SimulationManager withOutput(PrintStream output) { this.output = output; return this; } + /** + * Select output for statistics + * + * @param statisticOutput + * @return + */ public SimulationManager withStatisticOutput(PrintStream statisticOutput) { this.statisticOutput = statisticOutput; return this; } + /** + * Set field map as function + * + * @param fieldMap + * @return + */ public SimulationManager withFieldMap(UnivariateFunction fieldMap) { this.simulator.setFieldFunc(fieldMap); return this; } + /** + * Set field map from values + * + * @param z + * @param b + * @return + */ public SimulationManager withFieldMap(double[] z, double[] b) { this.simulator.setFieldFunc(new LinearInterpolator().interpolate(z, b)); return this; } - public SimulationManager withDensity(double density){ + /** + * Set source density + * PENDING replace by source thickness? + * + * @param density + * @return + */ + public SimulationManager withDensity(double density) { this.simulator.setGasDensity(density); return this; } @@ -68,7 +100,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() * Simulator.SOURCE_LENGTH; + double initZ = (generator.nextDouble()-1) * Simulator.SOURCE_LENGTH; Simulator.SimulationResult res = simulator.simulate(initialE, theta, initZ); if (reportIf.test(res)) { if (output != null) { @@ -83,7 +115,7 @@ public class SimulationManager { private double getRandomTheta() { double x = generator.nextDouble(); - // from 0 to 2 Pi + // from 0 to Pi return Math.acos(1 - 2 * x); } @@ -106,6 +138,9 @@ public class SimulationManager { } + /** + * Statistic counter + */ public static class Counter { int accepted = 0; int pass = 0; diff --git a/src/main/java/inr/numass/trapping/Simulator.java b/src/main/java/inr/numass/trapping/Simulator.java index dd4d24c..e6ad24a 100644 --- a/src/main/java/inr/numass/trapping/Simulator.java +++ b/src/main/java/inr/numass/trapping/Simulator.java @@ -11,11 +11,10 @@ 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; -import java.util.function.Function; - import static org.apache.commons.math3.util.FastMath.*; /** @@ -24,7 +23,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 dZ calculation + private static final double DELTA_L = 0.1; //step for deltaZ calculation private RandomGenerator generator; private Scatter scatter; @@ -138,24 +137,45 @@ public class Simulator { /** * calculate z coordinate change with known path length. Does not change position. * - * @param dl - * @return + * @param deltaL + * @return z shift and reflection counter */ - private double dZ(State position, double dl) { + private Pair deltaZ(State position, double deltaL) { // if magnetic field not defined, consider it to be uniform and equal bSource + int direction = position.isForward() ? 1 : -1; if (magneticField == null) { - return dl / cos(position.theta); + double deltaZ = direction * deltaL * cos(position.theta); + int reflectionCounter = (int) (deltaZ / SOURCE_LENGTH); + double curZ = normalizeZ(position.z + deltaZ); + return new Pair<>(curZ - position.z, reflectionCounter); } else { - double dz = 0; + int reflect = 0; + double curZ = position.z; 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; + + double sin2 = sin(position.theta) * sin(position.theta); + + 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; + root = -root; + reflect++; + } + + curZ += direction * delta * sqrt(root); + + curZ = normalizeZ(curZ); + + curL += delta; } - return dz; + return new Pair<>(curZ - position.z, reflect); } } @@ -169,18 +189,23 @@ public class Simulator { 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); } } + 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. @@ -188,17 +213,7 @@ public class Simulator { public SimulationResult simulate(double initEnergy, double initTheta, double initZ) { assert initEnergy > 0; assert initTheta > 0 && initTheta < Math.PI; - assert initZ >= 0 && initZ < SOURCE_LENGTH; - -// if (initTheta < this.thetaPinch) { -// if (generator.nextBoolean()) { -// return new SimulationResult(EndState.PASS, initEnergy, initTheta, initTheta, 0); -// } else { -// return new SimulationResult(EndState.REJECTED, initEnergy, initTheta, initTheta, 0); -// } -// } else if (initTheta < this.thetaTransport) { -// return new SimulationResult(EndState.REJECTED, initEnergy, initTheta, initTheta, 0); -// } + assert abs(initZ) <= SOURCE_LENGTH / 2d; State pos = new State(initEnergy, initTheta, initZ); EndState endState = EndState.NONE; @@ -208,15 +223,15 @@ public class Simulator { while (!stopflag) { double dl = freePath(pos.e); // path to next scattering - double dz = dZ(pos, dl); // z coordinate to next scattering - double expectedZ = pos.z + dz; // expected z position of next scattering + Pair dzCalc = deltaZ(pos, dl); + double dz = dzCalc.getFirst();// z coordinate to next scattering + int reflections = dzCalc.getSecond(); //if no scattering on current source pass - if (expectedZ < 0 || expectedZ > SOURCE_LENGTH) { + if (reflections > 0) { //accepted by spectrometer if (pos.theta < thetaPinch) { stopflag = true; - //Учитываем тот факт, что электрон мог вылететь в правильный угол, но назад if (pos.colNum == 0) { //counting pass electrons endState = EndState.PASS; @@ -226,23 +241,22 @@ public class Simulator { } - //through the rear magnetic trap + //through the rear magnetic pinch if (pos.theta >= PI - thetaTransport) { stopflag = true; endState = EndState.REJECTED; } - - if (pos.e < eLow) { - //Если энергия стала слишком маленькой - stopflag = true; - endState = EndState.LOWENERGY; - } } if (!stopflag) { // perform scatter propagate(pos, dl, dz); pos.colNum++; scatter(pos); + if (pos.e < eLow) { + //Если энергия стала слишком маленькой + stopflag = true; + endState = EndState.LOWENERGY; + } } } @@ -312,7 +326,7 @@ public class Simulator { int colNum = 0; /** - * current z; + * current z. Zero is the center of the source */ double z; @@ -337,13 +351,13 @@ public class Simulator { */ double addZ(double dZ) { this.z += dZ; - while (!(this.z > 0 && this.z < SOURCE_LENGTH)) { + while (abs(this.z) > SOURCE_LENGTH / 2d) { // reflecting from back wall if (z < 0) { - z = -z; + z += SOURCE_LENGTH / 2d; flip(); } else { - z -= SOURCE_LENGTH; + z -= SOURCE_LENGTH / 2d; flip(); } } @@ -375,7 +389,14 @@ public class Simulator { if (magneticField == null) { return theta; } else { - return asin(sin(theta) * sqrt(field() / bSource)); + double newTheta = asin(min(abs(sin(theta)) * sqrt(field() / bSource), 1)); + if (theta > PI / 2) { + newTheta = PI - newTheta; + } + if (Double.isNaN(newTheta)) { + throw new Error(); + } + return newTheta; } } @@ -406,7 +427,7 @@ public class Simulator { if (newtheta < 0) { newtheta = -newtheta; } - if (newtheta > Math.PI && newtheta <= Math.PI * 3 / 2) { + if (newtheta > Math.PI) { newtheta = 2 * Math.PI - newtheta; } @@ -415,7 +436,15 @@ public class Simulator { theta = newtheta; } else { theta = asin(sin(newtheta) * sqrt(bSource / field())); + if (newtheta > PI / 2) { + theta = PI - theta; + } } + + if (Double.isNaN(theta)) { + throw new Error(); + } + return theta; } diff --git a/src/main/java/inr/numass/trapping/Trapping.java b/src/main/java/inr/numass/trapping/Trapping.java index bf40655..a80b4b5 100644 --- a/src/main/java/inr/numass/trapping/Trapping.java +++ b/src/main/java/inr/numass/trapping/Trapping.java @@ -26,7 +26,8 @@ public class Trapping { new SimulationManager() .withParameters(0.6, 3.7, 4.84, 18000d, 4000) // .withFieldMap(z, b) - .simulateAll((int) 1e4); +// .withDensity(5e20) + .simulateAll((int) 1e5); 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()); diff --git a/src/main/resources/win32-amd64/libScatter.dll b/src/main/resources/win32-amd64/libScatter.dll deleted file mode 100644 index 6d41d27..0000000 Binary files a/src/main/resources/win32-amd64/libScatter.dll and /dev/null differ diff --git a/src/main/resources/win32-x86/libScatter.dll b/src/main/resources/win32-x86/libScatter.dll deleted file mode 100644 index be996e0..0000000 Binary files a/src/main/resources/win32-x86/libScatter.dll and /dev/null differ diff --git a/src/test/java/hep/dataforge/trapping/AppTest.java b/src/test/java/hep/dataforge/trapping/AppTest.java deleted file mode 100644 index 60af9b7..0000000 --- a/src/test/java/hep/dataforge/trapping/AppTest.java +++ /dev/null @@ -1,38 +0,0 @@ -package hep.dataforge.trapping; - -import junit.framework.Test; -import junit.framework.TestCase; -import junit.framework.TestSuite; - -/** - * Unit test for simple App. - */ -public class AppTest - extends TestCase -{ - /** - * Create the test case - * - * @param testName name of the test case - */ - public AppTest( String testName ) - { - super( testName ); - } - - /** - * @return the suite of tests being tested - */ - public static Test suite() - { - return new TestSuite( AppTest.class ); - } - - /** - * Rigourous Test :-) - */ - public void testApp() - { - assertTrue( true ); - } -}