commit c57cbc39c4702bc849562804d6d1e7e1f3a58950 Author: Darksnake Date: Wed Jul 10 15:07:58 2013 +0400 [no commit message] diff --git a/.hgignore b/.hgignore new file mode 100644 index 0000000..2f945c3 --- /dev/null +++ b/.hgignore @@ -0,0 +1,5 @@ +\.orig$ +\.orig\..*$ +\.chg\..*$ +\.rej$ +\.conflict\~$ diff --git a/nbactions.xml b/nbactions.xml new file mode 100644 index 0000000..d44d98c --- /dev/null +++ b/nbactions.xml @@ -0,0 +1,37 @@ + + + + run + + process-classes + org.codehaus.mojo:exec-maven-plugin:1.2.1:exec + + + -ea -classpath %classpath hep.dataforge.trapping.Trapping D:\\PlayGround\\Trapping.res + java + + + + debug + + process-classes + org.codehaus.mojo:exec-maven-plugin:1.2.1:exec + + + -Xdebug -Xrunjdwp:transport=dt_socket,server=n,address=${jpda.address} -ea -classpath %classpath ${packageClassName} D:\\PlayGround\\Trapping.res + java + true + + + + profile + + process-classes + org.codehaus.mojo:exec-maven-plugin:1.2.1:exec + + + ${profiler.args} -ea -classpath %classpath ${packageClassName} D:\\PlayGround\\Trapping.res + ${profiler.java} + + + diff --git a/pom.xml b/pom.xml new file mode 100644 index 0000000..8ce882d --- /dev/null +++ b/pom.xml @@ -0,0 +1,57 @@ + + 4.0.0 + + hep.DataForge + Trapping + 1.0-SNAPSHOT + jar + + Trapping + http://maven.apache.org + + + UTF-8 + + + + + + + org.apache.maven.plugins + maven-resources-plugin + 2.6 + + + + org.apache.maven.plugins + maven-compiler-plugin + 2.3.2 + + 1.7 + 1.7 + true + + + + + + + + junit + junit + 3.8.1 + test + + + org.apache.commons + commons-math3 + 3.2 + + + net.java.dev.jna + jna + 3.5.2 + + + diff --git a/src/main/java/hep/dataforge/trapping/DoubleValue.java b/src/main/java/hep/dataforge/trapping/DoubleValue.java new file mode 100644 index 0000000..a0938ba --- /dev/null +++ b/src/main/java/hep/dataforge/trapping/DoubleValue.java @@ -0,0 +1,33 @@ +/* + * To change this template, choose Tools | Templates + * and open the template in the editor. + */ +package hep.dataforge.trapping; + +/** + * Класс нужен исключительно чтобы сделать простой доступ к Сишным экспортам + * @author Darksnake + */ +public class DoubleValue { + private double value; + + public DoubleValue(double value) { + this.value = value; + } + + + + /** + * @return the value + */ + public double getValue() { + return value; + } + + /** + * @param value the value to set + */ + public void setValue(double value) { + this.value = value; + } +} diff --git a/src/main/java/hep/dataforge/trapping/ElectronTrappingSimulator.java b/src/main/java/hep/dataforge/trapping/ElectronTrappingSimulator.java new file mode 100644 index 0000000..c4645f1 --- /dev/null +++ b/src/main/java/hep/dataforge/trapping/ElectronTrappingSimulator.java @@ -0,0 +1,244 @@ +/* + * To change this template, choose Tools | Templates + * and open the template in the editor. + */ +package hep.dataforge.trapping; + +import java.util.ArrayList; +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.util.Precision; + +/** + * + * @author Darksnake + */ +public class ElectronTrappingSimulator { + + private TrappingRandomGenerator generator = new TrappingRandomGenerator(); + double Elow = 14000d; + double thetaTransport = 24.107064 / 180 * Math.PI; + double thetaPinch = 19.481097 / 180 * Math.PI; + + public ElectronTrappingSimulator() { + } + + public static enum EndState { + + ACCEPTED,//трэппинговый электрон попал в аксептанс + REJECTED,//трэппинговый электрон вылетел через заднюю пробку + LOWENERGY,//потерял слишком много энергии + PASS,//электрон никогда не запирался и прошел напрямую, нужно для нормировки + NONE + } + + public void setFields(double Bsource, double Btransport, double Bpinch) { + this.thetaTransport = Math.asin(Math.sqrt(Bsource / Btransport)); + this.thetaPinch = Math.asin(Math.sqrt(Bsource / Bpinch)); + } + + /** + * Симулируем один пробег электрона от начального значения и до вылетания из + * иточника или до того момента, как энергия становится меньше Elow. + * Считаем, что за один проход источника происходит в среднем существенно + * меньше одного столкновения, поэтому выход вперед или назад совершенно + * симментричны. + * + */ + public SimulaionResult simulateOne(double initEnergy, double initTheta) { + assert initEnergy > 0; + assert initTheta > 0 && initTheta < Math.PI / 2; + + if (initTheta < this.thetaPinch) { + if (generator.heads()) { + return new SimulaionResult(EndState.PASS, initEnergy, initTheta, 0); + } else { + return new SimulaionResult(EndState.REJECTED, initEnergy, initTheta, 0); + } + } else if(initTheta sigmaEl) { + if (alpha > sigmaEl + sigmaexc) { + //ionization case + Scatter.randomion(E, dE, dTheta); + } else { + //excitation case + Scatter.randomexc(E, dE, dTheta); + } + } else { + // elastic + Scatter.randomel(E, dE, dTheta); + } + +// if (alpha < sigmaEl) { +// Scatter.randomel(E, dE, dTheta); +// } else if (alpha < sigmaexc) { +// Scatter.randomexc(E, dE, dTheta); +// } else { +// Scatter.randomion(E, dE, dTheta); +// } + + + + + //Обновляем значени угла и энергии независимо ни от чего + E -= dE.getValue(); + //Изменение угла + theta = addTheta(theta, dTheta.getValue() / 180 * Math.PI); + //следим чтобы угол был от 0 до 90, если он перекинется через границу, считаем что электрон остается в потоке + theta = normalizeTheta(theta); + + if (theta < thetaPinch) { + stopflag = true; + if (generator.heads()) { + //Учитываем тот факт, что электрон мог вылететь в правильный угол, но назад + state = EndState.ACCEPTED; + } else { + state = EndState.REJECTED; + } + } + + if (theta >= thetaPinch && theta <= thetaTransport) { + //Вылет через заднюю пробку эквивалентен отражению от передней. + //это верно при низких концентрациях + stopflag = true; + state = EndState.REJECTED; + } + + if (E < Elow) { + //Если энергия стала слишком маленькой + stopflag = true; + state = EndState.LOWENERGY; + } + } + return new SimulaionResult(state, E, theta, colNum); + + } + + /** + * Сложение вектора с учетом случайно распределения по фи + * + * @param theta + * @param dTheta + * @return + */ + double addTheta(double theta, double dTheta) { + //Генерируем случайный фи + double phi = generator.next() * 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); + + Vector3D result = rot.applyTo(init.getCartesian()); + + // assert Vector3D.angle(result, rotate.getCartesian()) == dTheta; + return Math.acos(result.getZ()); + + } + + /** + * Симулируем пролет num электронов. + * + * @param E + * @param num + * @return + */ + public ArrayList simulateAll(double E, int num) { + ArrayList res = new ArrayList(); + double theta; + + System.out.printf("%nStarting sumulation with initial energy %g and %d electrons.%n%n", E, num); + for (int i = 0; i < num; i++) { +// System.out.printf("Running electron number %d", i); + theta = this.getRandomTheta(); + res.add(this.simulateOne(E, theta)); + } + return res; + } + + 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 double getRandomTheta() { + double x = generator.next(); + return Math.acos(x); + } + +// /** +// * Генерируем случайный угол таким образом, чтобы электрон заведомо был +// * траппинговый +// * +// * @return +// */ +// public double getRandomTrappedTheta() { +// double res = 0; +// while (res < this.thetaTransport) { +// res = this.getRandomTheta(); +// +// } +// return res; +// } + public class SimulaionResult { + + public SimulaionResult(EndState state, double E, double theta, int collisionNumber) { + this.state = state; + this.E = E; + this.theta = theta; + this.collisionNumber = collisionNumber; + } + public EndState state; + public double E; + public double theta; + public int collisionNumber; + } +} diff --git a/src/main/java/hep/dataforge/trapping/Scatter.java b/src/main/java/hep/dataforge/trapping/Scatter.java new file mode 100644 index 0000000..1f77bb1 --- /dev/null +++ b/src/main/java/hep/dataforge/trapping/Scatter.java @@ -0,0 +1,100 @@ +/* + * To change this template, choose Tools | Templates + * and open the template in the editor. + */ +package hep.dataforge.trapping; + +import com.sun.jna.Library; +import com.sun.jna.Native; +import com.sun.jna.ptr.DoubleByReference; + +/** + * + * @author Darksnake + */ +public class Scatter { + + public interface LibScatter extends Library { + + public void randomel(double E, DoubleByReference Eloss, DoubleByReference theta); + + public void randomexc(double E, DoubleByReference Eloss, DoubleByReference theta); + + public void randomion(double E, DoubleByReference Eloss, DoubleByReference theta); + + public double sigmael(double E); + + public double sigmaexc(double E); + + public double sigmaion(double E); + } + + /** + * PENDING переделать, чтобы возвращались нормальные значения + * @param E + * @param Eloss + * @param theta + */ + static void randomel(double E, DoubleValue Eloss, DoubleValue theta) { + LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class); + + DoubleByReference ElossPointer = new DoubleByReference(Eloss.getValue()); + DoubleByReference thetaPointer = new DoubleByReference(theta.getValue()); + lib.randomel(E, ElossPointer, thetaPointer); + Eloss.setValue(ElossPointer.getValue()); + theta.setValue(thetaPointer.getValue()); + + } + + static void randomexc(double E, DoubleValue Eloss, DoubleValue theta) { + LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class); + + DoubleByReference ElossPointer = new DoubleByReference(Eloss.getValue()); + DoubleByReference thetaPointer = new DoubleByReference(theta.getValue()); + lib.randomexc(E, ElossPointer, thetaPointer); + Eloss.setValue(ElossPointer.getValue()); + theta.setValue(thetaPointer.getValue()); + + } + + static void randomion(double E, DoubleValue Eloss, DoubleValue theta) { + LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class); + + DoubleByReference ElossPointer = new DoubleByReference(Eloss.getValue()); + DoubleByReference thetaPointer = new DoubleByReference(theta.getValue()); + lib.randomion(E, ElossPointer, thetaPointer); + Eloss.setValue(ElossPointer.getValue()); + theta.setValue(thetaPointer.getValue()); + + } + + /** + * Все сечения в м^2 + * @param E + * @return + */ + public static double sigmael(double E) { + LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class); + return lib.sigmael(E); + } + + public static double sigmaexc(double E) { + LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class); + return lib.sigmaexc(E); + } + + public static double sigmaion(double E) { + LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class); + return lib.sigmaion(E); + } + + /** + * Полное сечение с учетом квазиупругих столкновений + * @param E + * @return + */ + public static double sigmaTotal(double E){ + LibScatter lib = (LibScatter) Native.loadLibrary("libScatter", LibScatter.class); + return lib.sigmael(E)+lib.sigmaexc(E)+lib.sigmaion(E); + } +} diff --git a/src/main/java/hep/dataforge/trapping/Trapping.java b/src/main/java/hep/dataforge/trapping/Trapping.java new file mode 100644 index 0000000..56ce7eb --- /dev/null +++ b/src/main/java/hep/dataforge/trapping/Trapping.java @@ -0,0 +1,93 @@ +package hep.dataforge.trapping; + +import java.io.File; +import java.io.FileNotFoundException; +import java.io.PrintWriter; +import java.util.ArrayList; +import java.util.Iterator; + +/** + * Hello world! + * + */ +public class Trapping { + + public static void main(String[] args) throws FileNotFoundException { + PrintWriter out = null; + if ((args != null) && (args[0] != null)) { + File file = new File(args[0]); + out = new PrintWriter(file); + } + + double E = 18000d; + System.out.println(); +// System.setProperty("jna.library.path", "d:\\projects\\Trapping\\target\\classes\\win32-amd64\\libScatter.dll"); + ElectronTrappingSimulator simulator = new ElectronTrappingSimulator(); + simulator.setFields(0.6, 3.6, 7.2); +// ElectronTrappingSimulator.SimulaionResult result; + + int accepted = 0; + int pass = 0; + int rejected = 0; + int lowE = 0; + + ArrayList results = simulator.simulateAll(E, 100000); + + for (Iterator it = results.iterator(); it.hasNext();) { + ElectronTrappingSimulator.SimulaionResult res = it.next(); + + if (out != null) { + out.printf("%g\t%g\t%d\t%s%n", res.E, res.theta * 180 / Math.PI, res.collisionNumber, res.state.toString()); + } + switch (res.state) { + + case ACCEPTED: + System.out.printf("%g\t%g\t%d\t%s%n", res.E, res.theta * 180 / Math.PI, res.collisionNumber, res.state.toString()); + accepted++; + break; + case LOWENERGY: + lowE++; + break; + case PASS: + pass++; + break; + case REJECTED: + rejected++; + break; + case NONE: + throw new IllegalStateException(); + } + } + + System.out.printf("%nThe total number of events is %d.%n%n", results.size()); + + System.out.printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI); + System.out.printf("The transport mirroring angle is %g.%n", simulator.thetaTransport * 180 / Math.PI); + System.out.printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI); + System.out.printf("The starting energy is %g.%n", E); + System.out.printf("The lower energy boundary is %g.%n%n", simulator.Elow); + + System.out.printf("The total number of ACCEPTED events is %d.%n", accepted); + System.out.printf("The total number of PASS events is %d.%n", pass); + System.out.printf("The total number of REJECTED events is %d.%n", rejected); + System.out.printf("The total number of LOWENERGY events is %d.%n%n", lowE); + + if (out != null) { + out.printf("The total number of events is %d.%n%n", results.size()); + + out.printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI); + out.printf("The transport mirroring angle is %g.%n", simulator.thetaTransport * 180 / Math.PI); + out.printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI); + out.printf("The starting energy is %g.%n", E); + out.printf("The lower energy boundary is %g.%n%n", simulator.Elow); + + out.printf("The total number of ACCEPTED events is %d.%n", accepted); + out.printf("The total number of PASS events is %d.%n", pass); + out.printf("The total number of REJECTED events is %d.%n", rejected); + out.printf("The total number of LOWENERGY events is %d.%n%n", lowE); + out.close(); + } + + + } +} diff --git a/src/main/java/hep/dataforge/trapping/TrappingRandomGenerator.java b/src/main/java/hep/dataforge/trapping/TrappingRandomGenerator.java new file mode 100644 index 0000000..761a1dc --- /dev/null +++ b/src/main/java/hep/dataforge/trapping/TrappingRandomGenerator.java @@ -0,0 +1,43 @@ +/* + * To change this template, choose Tools | Templates + * and open the template in the editor. + */ +package hep.dataforge.trapping; + +import org.apache.commons.math3.random.JDKRandomGenerator; +import org.apache.commons.math3.random.RandomGenerator; + +/** + * + * @author Darksnake + */ +public class TrappingRandomGenerator { + RandomGenerator generator; + + public TrappingRandomGenerator() { + this.generator = new JDKRandomGenerator(); + } + + public TrappingRandomGenerator(RandomGenerator generator) { + this.generator = generator; + } + + /** + * heads-tails random. + * @return + */ + public boolean heads(){ + return generator.nextBoolean(); + } + + /** + * next uniform in [0;1] + * @return + */ + public double next(){ + return generator.nextDouble(); + } + + + +} diff --git a/src/main/resources/win32-amd64/libScatter.dll b/src/main/resources/win32-amd64/libScatter.dll new file mode 100644 index 0000000..6d41d27 Binary files /dev/null and b/src/main/resources/win32-amd64/libScatter.dll differ diff --git a/src/main/resources/win32-x86/libScatter.dll b/src/main/resources/win32-x86/libScatter.dll new file mode 100644 index 0000000..be996e0 Binary files /dev/null and b/src/main/resources/win32-x86/libScatter.dll differ diff --git a/src/test/java/hep/dataforge/trapping/AppTest.java b/src/test/java/hep/dataforge/trapping/AppTest.java new file mode 100644 index 0000000..60af9b7 --- /dev/null +++ b/src/test/java/hep/dataforge/trapping/AppTest.java @@ -0,0 +1,38 @@ +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 ); + } +}