diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4bebd19 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +/.gradle/ +/.idea/ +/build/ +/output/ +!gradle-wrapper.jar diff --git a/.hgignore b/.hgignore deleted file mode 100644 index 2300bba..0000000 --- a/.hgignore +++ /dev/null @@ -1,13 +0,0 @@ - -\.orig$ -\.orig\..*$ -\.chg\..*$ -\.rej$ -\.conflict\~$ -.gradle/ -output/ -.ipynb_checkpoints/ -notebooks/images/ -build/ -.idea/* -!gradle-wrapper.jar diff --git a/build.gradle b/build.gradle index d464d5c..9b36dbc 100644 --- a/build.gradle +++ b/build.gradle @@ -1,15 +1,14 @@ plugins { - id "org.jetbrains.kotlin.jvm" version "1.3.10" - id "java" + id "org.jetbrains.kotlin.jvm" version "1.4.10" id "application" } -group = 'inr.numass' -version = 'dev' +group = 'ru.inr.mass' +version = '1.0.0' description = "Numass trapping simulation" -mainClassName = "inr.numass.trapping.TrappingKt" +mainClassName = "ru.inr.mass.trapping.TrappingKt" repositories { mavenCentral() diff --git a/gradle/wrapper/gradle-wrapper.jar b/gradle/wrapper/gradle-wrapper.jar index a5fe1cb..e708b1c 100644 Binary files a/gradle/wrapper/gradle-wrapper.jar and b/gradle/wrapper/gradle-wrapper.jar differ diff --git a/gradle/wrapper/gradle-wrapper.properties b/gradle/wrapper/gradle-wrapper.properties index 3ba40c1..be52383 100644 --- a/gradle/wrapper/gradle-wrapper.properties +++ b/gradle/wrapper/gradle-wrapper.properties @@ -1,5 +1,5 @@ -distributionBase=GRADLE_USER_HOME -distributionPath=wrapper/dists -zipStoreBase=GRADLE_USER_HOME -zipStorePath=wrapper/dists -distributionUrl=https\://services.gradle.org/distributions/gradle-4.5-bin.zip +distributionBase=GRADLE_USER_HOME +distributionPath=wrapper/dists +distributionUrl=https\://services.gradle.org/distributions/gradle-6.7-bin.zip +zipStoreBase=GRADLE_USER_HOME +zipStorePath=wrapper/dists diff --git a/gradlew b/gradlew index cccdd3d..4f906e0 100644 --- a/gradlew +++ b/gradlew @@ -1,5 +1,21 @@ #!/usr/bin/env sh +# +# Copyright 2015 the original author or authors. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + ############################################################################## ## ## Gradle start up script for UN*X @@ -28,7 +44,7 @@ APP_NAME="Gradle" APP_BASE_NAME=`basename "$0"` # Add default JVM options here. You can also use JAVA_OPTS and GRADLE_OPTS to pass JVM options to this script. -DEFAULT_JVM_OPTS="" +DEFAULT_JVM_OPTS='"-Xmx64m" "-Xms64m"' # Use the maximum available, or set MAX_FD != -1 to use that value. MAX_FD="maximum" @@ -66,6 +82,7 @@ esac CLASSPATH=$APP_HOME/gradle/wrapper/gradle-wrapper.jar + # Determine the Java command to use to start the JVM. if [ -n "$JAVA_HOME" ] ; then if [ -x "$JAVA_HOME/jre/sh/java" ] ; then @@ -109,10 +126,11 @@ if $darwin; then GRADLE_OPTS="$GRADLE_OPTS \"-Xdock:name=$APP_NAME\" \"-Xdock:icon=$APP_HOME/media/gradle.icns\"" fi -# For Cygwin, switch paths to Windows format before running java -if $cygwin ; then +# For Cygwin or MSYS, switch paths to Windows format before running java +if [ "$cygwin" = "true" -o "$msys" = "true" ] ; then APP_HOME=`cygpath --path --mixed "$APP_HOME"` CLASSPATH=`cygpath --path --mixed "$CLASSPATH"` + JAVACMD=`cygpath --unix "$JAVACMD"` # We build the pattern for arguments to be converted via cygpath @@ -138,19 +156,19 @@ if $cygwin ; then else eval `echo args$i`="\"$arg\"" fi - i=$((i+1)) + i=`expr $i + 1` done case $i in - (0) set -- ;; - (1) set -- "$args0" ;; - (2) set -- "$args0" "$args1" ;; - (3) set -- "$args0" "$args1" "$args2" ;; - (4) set -- "$args0" "$args1" "$args2" "$args3" ;; - (5) set -- "$args0" "$args1" "$args2" "$args3" "$args4" ;; - (6) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" ;; - (7) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" "$args6" ;; - (8) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" "$args6" "$args7" ;; - (9) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" "$args6" "$args7" "$args8" ;; + 0) set -- ;; + 1) set -- "$args0" ;; + 2) set -- "$args0" "$args1" ;; + 3) set -- "$args0" "$args1" "$args2" ;; + 4) set -- "$args0" "$args1" "$args2" "$args3" ;; + 5) set -- "$args0" "$args1" "$args2" "$args3" "$args4" ;; + 6) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" ;; + 7) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" "$args6" ;; + 8) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" "$args6" "$args7" ;; + 9) set -- "$args0" "$args1" "$args2" "$args3" "$args4" "$args5" "$args6" "$args7" "$args8" ;; esac fi @@ -159,14 +177,9 @@ save () { for i do printf %s\\n "$i" | sed "s/'/'\\\\''/g;1s/^/'/;\$s/\$/' \\\\/" ; done echo " " } -APP_ARGS=$(save "$@") +APP_ARGS=`save "$@"` # Collect all arguments for the java command, following the shell quoting and substitution rules eval set -- $DEFAULT_JVM_OPTS $JAVA_OPTS $GRADLE_OPTS "\"-Dorg.gradle.appname=$APP_BASE_NAME\"" -classpath "\"$CLASSPATH\"" org.gradle.wrapper.GradleWrapperMain "$APP_ARGS" -# by default we should be in the correct project dir, but when run from Finder on Mac, the cwd is wrong -if [ "$(uname)" = "Darwin" ] && [ "$HOME" = "$PWD" ]; then - cd "$(dirname "$0")" -fi - exec "$JAVACMD" "$@" diff --git a/gradlew.bat b/gradlew.bat index e95643d..ac1b06f 100644 --- a/gradlew.bat +++ b/gradlew.bat @@ -1,3 +1,19 @@ +@rem +@rem Copyright 2015 the original author or authors. +@rem +@rem Licensed under the Apache License, Version 2.0 (the "License"); +@rem you may not use this file except in compliance with the License. +@rem You may obtain a copy of the License at +@rem +@rem https://www.apache.org/licenses/LICENSE-2.0 +@rem +@rem Unless required by applicable law or agreed to in writing, software +@rem distributed under the License is distributed on an "AS IS" BASIS, +@rem WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +@rem See the License for the specific language governing permissions and +@rem limitations under the License. +@rem + @if "%DEBUG%" == "" @echo off @rem ########################################################################## @rem @@ -13,15 +29,18 @@ if "%DIRNAME%" == "" set DIRNAME=. set APP_BASE_NAME=%~n0 set APP_HOME=%DIRNAME% +@rem Resolve any "." and ".." in APP_HOME to make it shorter. +for %%i in ("%APP_HOME%") do set APP_HOME=%%~fi + @rem Add default JVM options here. You can also use JAVA_OPTS and GRADLE_OPTS to pass JVM options to this script. -set DEFAULT_JVM_OPTS= +set DEFAULT_JVM_OPTS="-Xmx64m" "-Xms64m" @rem Find java.exe if defined JAVA_HOME goto findJavaFromJavaHome set JAVA_EXE=java.exe %JAVA_EXE% -version >NUL 2>&1 -if "%ERRORLEVEL%" == "0" goto init +if "%ERRORLEVEL%" == "0" goto execute echo. echo ERROR: JAVA_HOME is not set and no 'java' command could be found in your PATH. @@ -35,7 +54,7 @@ goto fail set JAVA_HOME=%JAVA_HOME:"=% set JAVA_EXE=%JAVA_HOME%/bin/java.exe -if exist "%JAVA_EXE%" goto init +if exist "%JAVA_EXE%" goto execute echo. echo ERROR: JAVA_HOME is set to an invalid directory: %JAVA_HOME% @@ -45,28 +64,14 @@ echo location of your Java installation. goto fail -:init -@rem Get command-line arguments, handling Windows variants - -if not "%OS%" == "Windows_NT" goto win9xME_args - -:win9xME_args -@rem Slurp the command line arguments. -set CMD_LINE_ARGS= -set _SKIP=2 - -:win9xME_args_slurp -if "x%~1" == "x" goto execute - -set CMD_LINE_ARGS=%* - :execute @rem Setup the command line set CLASSPATH=%APP_HOME%\gradle\wrapper\gradle-wrapper.jar + @rem Execute Gradle -"%JAVA_EXE%" %DEFAULT_JVM_OPTS% %JAVA_OPTS% %GRADLE_OPTS% "-Dorg.gradle.appname=%APP_BASE_NAME%" -classpath "%CLASSPATH%" org.gradle.wrapper.GradleWrapperMain %CMD_LINE_ARGS% +"%JAVA_EXE%" %DEFAULT_JVM_OPTS% %JAVA_OPTS% %GRADLE_OPTS% "-Dorg.gradle.appname=%APP_BASE_NAME%" -classpath "%CLASSPATH%" org.gradle.wrapper.GradleWrapperMain %* :end @rem End local scope for the variables with windows NT shell diff --git a/src/main/kotlin/inr/numass/trapping/FunctionCache.kt b/src/main/kotlin/ru/inr/mass/trapping/FunctionCache.kt similarity index 63% rename from src/main/kotlin/inr/numass/trapping/FunctionCache.kt rename to src/main/kotlin/ru/inr/mass/trapping/FunctionCache.kt index 88a235e..75460cb 100644 --- a/src/main/kotlin/inr/numass/trapping/FunctionCache.kt +++ b/src/main/kotlin/ru/inr/mass/trapping/FunctionCache.kt @@ -1,23 +1,23 @@ -package inr.numass.trapping - -import java.util.* - -class FunctionCache(val xPrecision: Double, val function: (Double) -> Double) { - private val values = TreeMap() - - operator fun invoke(x: Double): Double { - val floor: MutableMap.MutableEntry? = values.floorEntry(x) - val ceil: MutableMap.MutableEntry? = values.ceilingEntry(x) - if (floor == null || ceil == null || ceil.key - floor.key <= xPrecision) { - val newValue = function(x) - values[x] = newValue - return newValue - } else { - val x0 = floor.key - val x1 = ceil.key - val y0 = floor.value - val y1 = ceil.value - return y0 + (x - x0) * (y1 - y0) / (x1 - x0) - } - } +package ru.inr.mass.trapping + +import java.util.* + +class FunctionCache(private val xPrecision: Double, val function: (Double) -> Double) { + private val values = TreeMap() + + operator fun invoke(x: Double): Double { + val floor: MutableMap.MutableEntry? = values.floorEntry(x) + val ceil: MutableMap.MutableEntry? = values.ceilingEntry(x) + return if (floor == null || ceil == null || ceil.key - floor.key <= xPrecision) { + val newValue = function(x) + values[x] = newValue + newValue + } else { + val x0 = floor.key + val x1 = ceil.key + val y0 = floor.value + val y1 = ceil.value + y0 + (x - x0) * (y1 - y0) / (x1 - x0) + } + } } \ No newline at end of file diff --git a/src/main/kotlin/inr/numass/trapping/MultiCounter.kt b/src/main/kotlin/ru/inr/mass/trapping/MultiCounter.kt similarity index 93% rename from src/main/kotlin/inr/numass/trapping/MultiCounter.kt rename to src/main/kotlin/ru/inr/mass/trapping/MultiCounter.kt index b6c3f2c..9194b97 100644 --- a/src/main/kotlin/inr/numass/trapping/MultiCounter.kt +++ b/src/main/kotlin/ru/inr/mass/trapping/MultiCounter.kt @@ -1,96 +1,96 @@ -/* - * Copyright 2015 Alexander Nozik. - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package inr.numass.trapping - -import java.io.PrintStream -import java.lang.Integer.valueOf -import java.util.HashMap - -/** - * TODO есть объект MultiDimensionalCounter, исползовать его? - * - * @author Alexander Nozik - * @version $Id: $Id - */ -class MultiCounter(private var name: String) { - - private var counts = HashMap() - - /** - * - * getCount. - * - * @param name a [java.lang.String] object. - * @return a int. - */ - fun getCount(name: String): Int? { - return counts[name] - } - - /** - * - * increase. - * - * @param name a [java.lang.String] object. - */ - @Synchronized - fun increase(name: String) { - if (counts.containsKey(name)) { - val count = counts[name] - counts.remove(name) - counts[name] = count!! + 1 - } else { - counts[name] = valueOf(1) - } - - } - - /** - * - * print. - * - * @param out a [java.io.PrintWriter] object. - */ - fun print(out: PrintStream) { - out.printf("%nValues for counter %s%n%n", this.name) - for ((keyName, value) in counts) { - - out.printf("%s : %d%n", keyName, value) - } - } - - /** - * - * reset. - * - * @param name a [java.lang.String] object. - */ - @Synchronized - fun reset(name: String) { - if (counts.containsKey(name)) { - counts.remove(name) - } - } - - /** - * - * resetAll. - */ - @Synchronized - fun resetAll() { - this.counts = HashMap() - } -} +/* + * Copyright 2015 Alexander Nozik. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package ru.inr.mass.trapping + +import java.io.PrintStream +import java.lang.Integer.valueOf +import java.util.* + +/** + * TODO есть объект MultiDimensionalCounter, исползовать его? + * + * @author Alexander Nozik + * @version $Id: $Id + */ +class MultiCounter(private var name: String) { + + private var counts = HashMap() + + /** + * + * getCount. + * + * @param name a [java.lang.String] object. + * @return a int. + */ + fun getCount(name: String): Int? { + return counts[name] + } + + /** + * + * increase. + * + * @param name a [java.lang.String] object. + */ + @Synchronized + fun increase(name: String) { + if (counts.containsKey(name)) { + val count = counts[name] + counts.remove(name) + counts[name] = count!! + 1 + } else { + counts[name] = valueOf(1) + } + + } + + /** + * + * print. + * + * @param out a [java.io.PrintWriter] object. + */ + fun print(out: PrintStream) { + out.printf("%nValues for counter %s%n%n", this.name) + for ((keyName, value) in counts) { + + out.printf("%s : %d%n", keyName, value) + } + } + + /** + * + * reset. + * + * @param name a [java.lang.String] object. + */ + @Synchronized + fun reset(name: String) { + if (counts.containsKey(name)) { + counts.remove(name) + } + } + + /** + * + * resetAll. + */ + @Synchronized + fun resetAll() { + this.counts = HashMap() + } +} diff --git a/src/main/kotlin/inr/numass/trapping/RandomGeneratorBridge.kt b/src/main/kotlin/ru/inr/mass/trapping/RandomGeneratorBridge.kt similarity index 85% rename from src/main/kotlin/inr/numass/trapping/RandomGeneratorBridge.kt rename to src/main/kotlin/ru/inr/mass/trapping/RandomGeneratorBridge.kt index 56bf65e..939c963 100644 --- a/src/main/kotlin/inr/numass/trapping/RandomGeneratorBridge.kt +++ b/src/main/kotlin/ru/inr/mass/trapping/RandomGeneratorBridge.kt @@ -1,28 +1,28 @@ -package inr.numass.trapping - -import org.apache.commons.math3.random.RandomGenerator -import org.apache.commons.rng.UniformRandomProvider - -class RandomGeneratorBridge(val provider: UniformRandomProvider) : RandomGenerator { - override fun nextBoolean(): Boolean = provider.nextBoolean() - - override fun nextFloat(): Float = provider.nextFloat() - - override fun setSeed(seed: Int) = error("Not supported") - - override fun setSeed(seed: IntArray?) = error("Not supported") - - override fun setSeed(seed: Long) = error("Not supported") - - override fun nextBytes(bytes: ByteArray) = provider.nextBytes(bytes) - - override fun nextInt(): Int = provider.nextInt() - - override fun nextInt(n: Int): Int = provider.nextInt(n) - - override fun nextGaussian(): Double = error("Not supported") - - override fun nextDouble(): Double = provider.nextDouble() - - override fun nextLong(): Long = provider.nextLong() +package ru.inr.mass.trapping + +import org.apache.commons.math3.random.RandomGenerator +import org.apache.commons.rng.UniformRandomProvider + +class RandomGeneratorBridge(private val provider: UniformRandomProvider) : RandomGenerator { + override fun nextBoolean(): Boolean = provider.nextBoolean() + + override fun nextFloat(): Float = provider.nextFloat() + + override fun setSeed(seed: Int) = error("Not supported") + + override fun setSeed(seed: IntArray?) = error("Not supported") + + override fun setSeed(seed: Long) = error("Not supported") + + override fun nextBytes(bytes: ByteArray) = provider.nextBytes(bytes) + + override fun nextInt(): Int = provider.nextInt() + + override fun nextInt(n: Int): Int = provider.nextInt(n) + + override fun nextGaussian(): Double = error("Not supported") + + override fun nextDouble(): Double = provider.nextDouble() + + override fun nextLong(): Long = provider.nextLong() } \ No newline at end of file diff --git a/src/main/kotlin/inr/numass/trapping/Scatter.kt b/src/main/kotlin/ru/inr/mass/trapping/Scatter.kt similarity index 97% rename from src/main/kotlin/inr/numass/trapping/Scatter.kt rename to src/main/kotlin/ru/inr/mass/trapping/Scatter.kt index cbe576f..bddd60f 100644 --- a/src/main/kotlin/inr/numass/trapping/Scatter.kt +++ b/src/main/kotlin/ru/inr/mass/trapping/Scatter.kt @@ -1,983 +1,982 @@ -package inr.numass.trapping - -import org.apache.commons.math3.random.RandomGenerator -import org.apache.commons.math3.util.FastMath.* -import org.apache.commons.math3.util.Pair - -/** - * - * The initial code was written by Ferenc Glueck and then rewritten by [Sebastian Voecking](mailto:seb.voeck@uni-muenster.de) - * into C++. To reference the code, don't forget to also include [kasiopea](http://stacks.iop.org/1367-2630/19/i=5/a=053012). - * @author Darksnake - */ -object Scatter { - - var debug = false - - val counter = MultiCounter("Accept-reject calls") - - - private val a02 = 28e-22 // Bohr radius squared - private val clight = 137.0 // velocity of light in atomic units - private val emass = 18780.0 // Electron mass in atomic units - private val R = 13.6 // Ryberg energy in eV - - private fun count(counterName: String) { - if (debug) { - counter.increase(counterName) - } - } - - fun randomel(E: Double, generator: RandomGenerator): Pair { - // This subroutine generates energy loss and polar scatt. angle according to - // electron elastic scattering in molecular hydrogen. - // Input: - // E: incident electron energy in eV. - // Output: - // *Eloss: energy loss in eV - // *theta: change of polar angle in degrees - val H2molmass = 69e6 - val T: Double = E / 27.2 - var c = 1.0 - val b: Double - var G: Double - var a: Double - val gam: Double - var K2: Double - val Gmax: Double = when { - E >= 250.0 -> 1e-19 - E < 250.0 && E >= 150.0 -> 2.5e-19 - else -> 1e-18 - } - var i: Int = 1 - - count("randomel-calls") - gam = 1.0 + T / (clight * clight) // relativistic correction factor - b = 2.0 / (1.0 + gam) / T - while (i < 5000) { - count("randomel") - c = 1.0 + b - b * (2.0 + b) / (b + 2.0 * generator.nextDouble()) - K2 = 2.0 * T * (1.0 + gam) * abs(1.0 - c) // momentum transfer squared - a = (4.0 + K2) * (4.0 + K2) / (gam * gam) - G = a * dElastic(E, c) - if (G > Gmax * generator.nextDouble()) { - break - } - i++ - } - return Pair(2.0 * emass / H2molmass * (1.0 - c) * E, acos(c) * 180.0 / Math.PI) - } - - - // Energy values of the excited electronic states: - // (from Mol. Phys. 41 (1980) 1501, in Hartree atomic units) - private val En = doubleArrayOf(12.73 / 27.2, 13.2 / 27.2, 14.77 / 27.2, 15.3 / 27.2, 14.93 / 27.2, 15.4 / 27.2, 13.06 / 27.2) - // Probability numbers of the electronic states: - // (from testelectron7.c calculation ) - private val p = doubleArrayOf(35.86, 40.05, 6.58, 2.26, 9.61, 4.08, 1.54) - // Energy values of the B vibrational states: - // (from: Phys. Rev. A51 (1995) 3745 , in Hartree atomic units) - private val EB = doubleArrayOf(0.411, 0.417, 0.423, 0.428, 0.434, 0.439, 0.444, 0.449, 0.454, 0.459, 0.464, 0.468, 0.473, 0.477, 0.481, 0.485, 0.489, 0.493, 0.496, 0.500, 0.503, 0.507, 0.510, 0.513, 0.516, 0.519, 0.521, 0.524) - // Energy values of the C vibrational states: - // (from: Phys. Rev. A51 (1995) 3745 , in Hartree atomic units) - private val EC = doubleArrayOf(0.452, 0.462, 0.472, 0.481, 0.490, 0.498, 0.506, 0.513, 0.519, 0.525, 0.530, 0.534, 0.537, 0.539) - // Franck-Condon factors of the B vibrational states: - // (from: Phys. Rev. A51 (1995) 3745 ) - private val pB = doubleArrayOf(4.2e-3, 1.5e-2, 3.0e-2, 4.7e-2, 6.3e-2, 7.3e-2, 7.9e-2, 8.0e-2, 7.8e-2, 7.3e-2, 6.6e-2, 5.8e-2, 5.1e-2, 4.4e-2, 3.7e-2, 3.1e-2, 2.6e-2, 2.2e-2, 1.8e-2, 1.5e-2, 1.3e-2, 1.1e-2, 8.9e-3, 7.4e-3, 6.2e-3, 5.2e-3, 4.3e-3, 3.6e-3) - // Franck-Condon factors of the C vibrational states: - // (from: Phys. Rev. A51 (1995) 3745 ) - private val pC = doubleArrayOf(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) - - private val Ecen = 12.6 / 27.21 - - private val fmax by lazy { - val sum = DoubleArray(1001) - - val T: Double = 20000.0 / 27.2 - val xmin = Ecen * Ecen / (2.0 * T) - val ymin = log(xmin) - val ymax = log(8.0 * T + xmin) - val dy = (ymax - ymin) / 1000.0 - - - var res = 0.0 - var i = 0 - while (i <= 1000) { - val y = ymin + dy * i - val K = exp(y / 2.0) - sum[i] = sumexc(K) - if (sum[i] > res) { - res = sum[i] - } - i++ - } - res *= 1.05 - res - } - - - fun randomexc(E: Double, generator: RandomGenerator): Pair { - // This subroutine generates energy loss and polar scatt. angle according to - // electron excitation scattering in molecular hydrogen. - // Input: - // E: incident electron energy in eV. - // Output: - // *Eloss: energy loss in eV - // *theta: change of polar angle in degrees - - - var c = 0.0 - var K: Double - var y: Double - var fy: Double - var pmax: Double - var D: Double - var j: Int - var n = 0 - var N: Int - var v = 0 - - count("randomexc-calls") - - // - // Scattering angle *theta generation: - // - val T = E / 27.2 - val theta: Double - if (E >= 100.0) { - val xmin = Ecen * Ecen / (2.0 * T) - val ymin = log(xmin) - val ymax = log(8.0 * T + xmin) - // dy = (ymax - ymin) / 1000.; - // Generation of y values with the Neumann acceptance-rejection method: - y = ymin - j = 1 - while (j < 5000) { - count("randomexc1") - y = ymin + (ymax - ymin) * generator.nextDouble() - K = exp(y / 2.0) - fy = sumexc(K) - if (fmax * generator.nextDouble() < fy) { - break - } - j++ - } - // Calculation of c=cos(theta) and theta: - val x = exp(y) - c = 1.0 - (x - xmin) / (4.0 * T) - theta = acos(c) * 180.0 / Math.PI - } else { - val Dmax = when { - E <= 25.0 -> 60.0 - E > 25.0 && E <= 35.0 -> 95.0 - E > 35.0 && E <= 50.0 -> 150.0 - else -> 400.0 - } - j = 1 - - while (j < 5000) { - count("randomexc2") - c = -1.0 + 2.0 * generator.nextDouble() - D = dExcitation(E, c) * 1e22 - if (Dmax * generator.nextDouble() < D) { - break - } - j++ - } - theta = acos(c) * 180.0 / Math.PI - } - // Energy loss *Eloss generation: - - // First we generate the electronic state, using the Neumann - // acceptance-rejection method for discrete distribution: - N = 7 // the number of electronic states in our calculation - pmax = p[1] // the maximum of the p[] values - j = 1 - while (j < 5000) { - count("randomexc3") - n = (N * generator.nextDouble()).toInt() - if (generator.nextDouble() * pmax < p[n]) { - break - } - j++ - }// Bp, Bpp, D, Dp, EF states - // C state; we generate now a vibrational state, - // using the Franck-Condon factors - // the number of C vibrational states in our calculation - // maximum of the pC[] values - // B state; we generate now a vibrational state, - // using the Frank-Condon factors - // the number of B vibrational states in our calculation - // maximum of the pB[] values - when { - n < 0 -> n = 0 - n > 6 -> n = 6 - } - val Eloss: Double - when (n) { - 0 -> { - // B state; we generate now a vibrational state, - // using the Frank-Condon factors - N = 28 // the number of B vibrational states in our calculation - pmax = pB[7] // maximum of the pB[] values - j = 1 - while (j < 5000) { - count("randomexc4") - v = (N * generator.nextDouble()).toInt() - if (generator.nextDouble() * pmax < pB[v]) { - break - } - j++ - } - if (v < 0) { - v = 0 - } - if (v > 27) { - v = 27 - } - Eloss = EB[v] * 27.2 - } - 1 -> { - // C state; we generate now a vibrational state, - // using the Franck-Condon factors - N = 14 // the number of C vibrational states in our calculation - pmax = pC[1] // maximum of the pC[] values - j = 1 - while (j < 5000) { - count("randomexc4") - v = (N * generator.nextDouble()).toInt() - if (generator.nextDouble() * pmax < pC[v]) { - break - } - j++ - } - if (v < 0) { - v = 0 - } - if (v > 13) { - v = 13 - } - Eloss = EC[v] * 27.2 - } - else -> - // Bp, Bpp, D, Dp, EF states - Eloss = En[n] * 27.2 - } - return Pair(Eloss, theta) - } - - fun randomion(E: Double, generator: RandomGenerator): Pair { - // This subroutine generates energy loss and polar scatt. angle according to - // electron ionization scattering in molecular hydrogen. - // Input: - // E: incident electron energy in eV. - // Output: - // *Eloss: energy loss in eV - // *theta: change of polar angle in degrees - // The kinetic energy of the secondary electron is: Eloss-15.4 eV - // - val Ei = 15.45 / 27.21 - var c: Double - val b: Double - var K: Double - val xmin: Double - val ymin: Double - val ymax: Double - val x: Double - var y: Double - val T: Double = E / 27.2 - var G: Double - var Gmax: Double - var q: Double - val h: Double - var F: Double - val Fmin: Double - val Fmax: Double - var Gp: Double - val Elmin: Double - val Elmax: Double = (E + 15.45) / 2.0 / 27.2 - val qmin: Double - val qmax: Double - var El: Double - val wmax: Double - var WcE: Double - var Jstarq: Double - var WcstarE: Double - var w: Double - var D2ion: Double - var j: Int = 1 - var K2: Double - var KK: Double - var fE: Double - var kej: Double - var ki: Double - var kf: Double - var Rex: Double - var arg: Double - var arctg: Double - var i: Int - var st1: Double - var st2: Double - count("randomion-calls") - // - // I. Generation of theta - // ----------------------- - Gmax = 1e-20 - if (E < 200.0) { - Gmax = 2e-20 - } - xmin = Ei * Ei / (2.0 * T) - b = xmin / (4.0 * T) - ymin = log(xmin) - ymax = log(8.0 * T + xmin) - // Generation of y values with the Neumann acceptance-rejection method: - y = ymin - while (j < 5000) { - count("randomion1") - y = ymin + (ymax - ymin) * generator.nextDouble() - K = exp(y / 2.0) - c = 1.0 + b - K * K / (4.0 * T) - G = K * K * (dInelastic(E, c) - dExcitation(E, c)) - if (Gmax * generator.nextDouble() < G) { - break - } - j++ - } - // y --> x --> c --> theta - x = exp(y) - c = 1.0 - (x - xmin) / (4.0 * T) - val theta = acos(c) * 180.0 / Math.PI - // - // II. Generation of Eloss, for fixed theta - // ---------------------------------------- - // - // For E<=100 eV we use subr. gensecelen - // (in this case no correlation between theta and Eloss) - if (E <= 100.0) { - return Pair(15.45 + gensecelen(E, generator), theta) - } - // For theta>=20 the free electron model is used - // (with full correlation between theta and Eloss) - if (theta >= 20.0) { - return Pair(E * (1.0 - c * c), theta) - } - // For E>100 eV and theta<20: analytical first Born approximation - // formula of Bethe for H atom (with modification for H2) - // - // Calc. of wmax: - if (theta >= 0.7) { - wmax = 1.1 - } else if (theta <= 0.7 && theta > 0.2) { - wmax = 2.0 - } else if (theta <= 0.2 && theta > 0.05) { - wmax = 4.0 - } else { - wmax = 8.0 - } - // We generate the q value according to the Jstarq pdf. We have to - // define the qmin and qmax limits for this generation: - K = sqrt(4.0 * T * (1.0 - Ei / (2.0 * T) - sqrt(1.0 - Ei / T) * c)) - Elmin = Ei - qmin = Elmin / K - K / 2.0 - qmax = Elmax / K - K / 2.0 - // - q = qmax - Fmax = 1.0 / 2.0 + 1.0 / Math.PI * (q / (1.0 + q * q) + atan(q)) - q = qmin - Fmin = 1.0 / 2.0 + 1.0 / Math.PI * (q / (1.0 + q * q) + atan(q)) - h = Fmax - Fmin - // Generation of Eloss with the Neumann acceptance-rejection method: - El = 0.0 - j = 1 - while (j < 5000) { - // Generation of q with inverse transform method - // (we use the Newton-Raphson method in order to solve the nonlinear eq. - // for the inversion) : - count("randomion2") - F = Fmin + h * generator.nextDouble() - y = 0.0 - i = 1 - while (i <= 30) { - G = 1.0 / 2.0 + (y + sin(2.0 * y) / 2.0) / Math.PI - Gp = (1.0 + cos(2.0 * y)) / Math.PI - y -= (G - F) / Gp - if (abs(G - F) < 1e-8) { - break - } - i++ - } - q = tan(y) - // We have the q value, so we can define El, and calculate the weight: - El = q * K + K * K / 2.0 - // First Born approximation formula of Bethe for e-H ionization: - KK = K - ki = sqrt(2.0 * T) - kf = sqrt(2.0 * (T - El)) - K2 = 4.0 * T * (1.0 - El / (2.0 * T) - sqrt(1.0 - El / T) * c) - if (K2 < 1e-9) { - K2 = 1e-9 - } - K = sqrt(K2) // momentum transfer - Rex = 1.0 - K * K / (kf * kf) + K2 * K2 / (kf * kf * kf * kf) - kej = sqrt(2.0 * abs(El - Ei) + 1e-8) - st1 = K2 - 2.0 * El + 2.0 - if (abs(st1) < 1e-9) { - st1 = 1e-9 - } - arg = 2.0 * kej / st1 - arctg = if (arg >= 0.0) { - atan(arg) - } else { - atan(arg) + Math.PI - } - st1 = (K + kej) * (K + kej) + 1.0 - st2 = (K - kej) * (K - kej) + 1.0 - fE = 1024.0 * El * (K2 + 2.0 / 3.0 * El) / (st1 * st1 * st1 * st2 * st2 * st2) * exp(-2.0 / kej * arctg) / (1.0 - exp(-2.0 * Math.PI / kej)) - D2ion = 2.0 * kf / ki * Rex / (El * K2) * fE - K = KK - // - WcE = D2ion - Jstarq = 16.0 / (3.0 * Math.PI * (1.0 + q * q) * (1.0 + q * q)) - WcstarE = 4.0 / (K * K * K * K * K) * Jstarq - w = WcE / WcstarE - if (wmax * generator.nextDouble() < w) { - break - } - j++ - } - - return Pair(El * 27.2, theta) - } - - private fun gensecelen(E: Double, generator: RandomGenerator): Double { - // This subroutine generates secondary electron energy W - // from ionization of incident electron energy E, by using - // the Lorentzian of Aseev et al. (Eq. 8). - // E and W in eV. - val Ei = 15.45 - val eps2 = 14.3 - val b = 6.25 - val B: Double = atan((Ei - eps2) / b) - val D: Double - val A: Double - val eps: Double - val a: Double - val u: Double = generator.nextDouble() - val epsmax: Double - - epsmax = (E + Ei) / 2.0 - A = atan((epsmax - eps2) / b) - D = b / (A - B) - a = b / D * (u + D / b * B) - eps = eps2 + b * tan(a) - return eps - Ei - } - - - // This subroutine computes the differential cross section - // dElastic= d sigma/d Omega of elastic electron scattering - // on molecular hydrogen. - // See: Nishimura et al., J. Phys. Soc. Jpn. 54 (1985) 1757. - // Input: E= electron kinetic energy in eV - // c= cos(theta), where theta is the polar scatt. angle - // dElastic: in m^2/steradian - private val elasticCel = doubleArrayOf(-0.512, -0.512, -0.509, -0.505, -0.499, -0.491, -0.476, -0.473, -0.462, -0.452, -0.438, -0.422, -0.406, -0.388, -0.370, -0.352, -0.333, -0.314, -0.296, -0.277, -0.258, -0.239, -0.221, -0.202, -0.185, -0.167, -0.151, -0.135, -0.120, -0.105, -0.092, -0.070, -0.053, -0.039, -0.030, -0.024, -0.019, -0.016, -0.014, -0.013, -0.012, -0.009, -0.008, -0.006, -0.005, -0.004, -0.003, -0.002, -0.002, -0.001) - private val elasticE = doubleArrayOf(0.0, 3.0, 6.0, 12.0, 20.0, 32.0, 55.0, 85.0, 150.0, 250.0) - private val elasticT = doubleArrayOf(0.0, 10.0, 20.0, 30.0, 40.0, 60.0, 80.0, 100.0, 140.0, 180.0) - private val elasticD = arrayOf( - doubleArrayOf(2.9, 2.7, 2.5, 2.1, 1.8, 1.2, 0.9, 1.0, 1.6, 1.9), - doubleArrayOf(4.2, 3.6, 3.1, 2.5, 1.9, 1.1, 0.8, 0.9, 1.3, 1.4), - doubleArrayOf(6.0, 4.4, 3.2, 2.3, 1.8, 1.1, 0.7, 0.54, 0.5, 0.6), - doubleArrayOf(6.0, 4.1, 2.8, 1.9, 1.3, 0.6, 0.3, 0.17, 0.16, 0.23), - doubleArrayOf(4.9, 3.2, 2.0, 1.2, 0.8, 0.3, 0.15, 0.09, 0.05, 0.05), - doubleArrayOf(5.2, 2.5, 1.2, 0.64, 0.36, 0.13, 0.05, 0.03, 0.016, 0.02), - doubleArrayOf(4.0, 1.7, 0.7, 0.3, 0.16, 0.05, 0.02, 0.013, 0.01, 0.01), - doubleArrayOf(2.8, 1.1, 0.4, 0.15, 0.07, 0.02, 0.01, 0.007, 0.004, 0.003), - doubleArrayOf(1.2, 0.53, 0.2, 0.08, 0.03, 0.0074, 0.003, 0.0016, 0.001, 0.0008) - ) - - - private fun dElastic(E: Double, c: Double): Double { - val T: Double = E / 27.2 - var K2: Double - var K: Double - val d: Double - val st1: Double - val st2: Double - val DH: Double - val gam: Double - var Delreturn = 0.0 - val CelK: Double - val Ki: Double - val theta: Double - var i: Int - var j: Int - if (E >= 250.0) { - gam = 1.0 + T / (clight * clight) // relativistic correction factor - K2 = 2.0 * T * (1.0 + gam) * (1.0 - c) - if (K2 < 0.0) { - K2 = 1e-30 - } - K = sqrt(K2) - if (K < 1e-9) { - K = 1e-9 // momentum transfer - } - d = 1.4009 // distance of protons in H2 - st1 = 8.0 + K2 - st2 = 4.0 + K2 - // DH is the diff. cross section for elastic electron scatt. - // on atomic hydrogen within the first Born approximation : - DH = 4.0 * st1 * st1 / (st2 * st2 * st2 * st2) * a02 - // CelK calculation with linear interpolation. - // CelK is the correction of the elastic electron - // scatt. on molecular hydrogen compared to the independent atom - // model. - when { - K < 3.0 -> { - i = (K / 0.1).toInt() - Ki = i * 0.1 - CelK = elasticCel[i] + (K - Ki) / 0.1 * (elasticCel[i + 1] - elasticCel[i]) - } - K >= 3.0 && K < 5.0 -> { - i = (30 + (K - 3.0) / 0.2).toInt() - Ki = 3.0 + (i - 30) * 0.2 - CelK = elasticCel[i] + (K - Ki) / 0.2 * (elasticCel[i + 1] - elasticCel[i]) - } - K >= 5.0 && K < 9.49 -> { - i = (40 + (K - 5.0) / 0.5).toInt() - Ki = 5.0 + (i - 40) * 0.5 - CelK = elasticCel[i] + (K - Ki) / 0.5 * (elasticCel[i + 1] - elasticCel[i]) - } - else -> CelK = 0.0 - } - Delreturn = 2.0 * gam * gam * DH * (1.0 + sin(K * d) / (K * d)) * (1.0 + CelK) - } else { - theta = acos(c) * 180.0 / Math.PI - i = 0 - while (i <= 8) { - if (E >= elasticE[i] && E < elasticE[i + 1]) { - j = 0 - while (j <= 8) { - if (theta >= elasticT[j] && theta < elasticT[j + 1]) { - Delreturn = 1e-20 * (elasticD[i][j] + (elasticD[i][j + 1] - elasticD[i][j]) * (theta - elasticT[j]) / (elasticT[j + 1] - elasticT[j])) - } - j++ - } - } - i++ - } - } - return Delreturn - } - - private val excEE = 12.6 / 27.2 - private val excitationE = doubleArrayOf(0.0, 25.0, 35.0, 50.0, 100.0) - private val excitationT = doubleArrayOf(0.0, 10.0, 20.0, 30.0, 40.0, 60.0, 80.0, 100.0, 180.0) - private val excitationD = arrayOf(doubleArrayOf(60.0, 43.0, 27.0, 18.0, 13.0, 8.0, 6.0, 6.0, 6.0), doubleArrayOf(95.0, 70.0, 21.0, 9.0, 6.0, 3.0, 2.0, 2.0, 2.0), doubleArrayOf(150.0, 120.0, 32.0, 8.0, 3.7, 1.9, 1.2, 0.8, 0.8), doubleArrayOf(400.0, 200.0, 12.0, 2.0, 1.4, 0.7, 0.3, 0.2, 0.2)) - - - private fun dExcitation(E: Double, c: Double): Double { - // This subroutine computes the differential cross section - // dElastic= d sigma/d Omega of excitation electron scattering - // on molecular hydrogen. - // Input: E= electron kinetic energy in eV - // c= cos(theta), where theta is the polar scatt. angle - // dExcitation: in m^2/steradian - var K2: Double - val K: Double - var sigma = 0.0 - val T: Double = E / 27.2 - val theta: Double - var i: Int - var j: Int - // - if (E >= 100.0) { - K2 = 4.0 * T * (1.0 - excEE / (2.0 * T) - sqrt(1.0 - excEE / T) * c) - if (K2 < 1e-9) { - K2 = 1e-9 - } - K = sqrt(K2) // momentum transfer - sigma = 2.0 / K2 * sumexc(K) * a02 - } else if (E <= 10.0) { - sigma = 0.0 - } else { - theta = acos(c) * 180.0 / Math.PI - i = 0 - while (i <= 3) { - if (E >= excitationE[i] && E < excitationE[i + 1]) { - j = 0 - while (j <= 7) { - if (theta >= excitationT[j] && theta < excitationT[j + 1]) { - sigma = 1e-22 * (excitationD[i][j] + (excitationD[i][j + 1] - excitationD[i][j]) * (theta - excitationT[j]) / (excitationT[j + 1] - excitationT[j])) - } - j++ - } - } - i++ - } - } - return sigma - } - - // This subroutine computes the differential cross section - // dInelastic= d sigma/d Omega of inelastic electron scattering - // on molecular hydrogen, within the first Born approximation. - // Input: E= electron kinetic energy in eV - // c= cos(theta), where theta is the polar scatt. angle - // dInelastic: in m2/steradian - private val cInelastic = doubleArrayOf(-0.246, -0.244, -0.239, -0.234, -0.227, -0.219, -0.211, -0.201, -0.190, -0.179, -0.167, -0.155, -0.142, -0.130, -0.118, -0.107, -0.096, -0.085, -0.076, -0.067, -0.059, -0.051, -0.045, -0.039, -0.034, -0.029, -0.025, -0.022, -0.019, -0.016, -0.014, -0.010, -0.008, -0.006, -0.004, -0.003, -0.003, -0.002, -0.002, -0.001, -0.001, -0.001, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000) - private val inelasticEi = 0.568 - - - private fun dInelastic(E: Double, c: Double): Double { - val T: Double = E / 27.2 - var K2: Double - val K: Double - val st1: Double - val F: Double - val DH: Double - val Dinelreturn: Double - val CinelK: Double - val Ki: Double - val i: Int - if (E < 16.0) { - return dExcitation(E, c) - } - K2 = 4.0 * T * (1.0 - inelasticEi / (2.0 * T) - sqrt(1.0 - inelasticEi / T) * c) - if (K2 < 1e-9) { - K2 = 1e-9 - } - K = sqrt(K2) // momentum transfer - st1 = 1.0 + K2 / 4.0 - F = 1.0 / (st1 * st1) // scatt. formfactor of hydrogen atom - // DH is the diff. cross section for inelastic electron scatt. - // on atomic hydrogen within the first Born approximation : - DH = 4.0 / (K2 * K2) * (1.0 - F * F) * a02 - // CinelK calculation with linear interpolation. - // CinelK is the correction of the inelastic electron - // scatt. on molecular hydrogen compared to the independent atom - // model. - if (K < 3.0) { - i = (K / 0.1).toInt() - Ki = i * 0.1 - CinelK = cInelastic[i] + (K - Ki) / 0.1 * (cInelastic[i + 1] - cInelastic[i]) - } else if (K >= 3.0 && K < 5.0) { - i = (30 + (K - 3.0) / 0.2).toInt() - Ki = 3.0 + (i - 30) * 0.2 - CinelK = cInelastic[i] + (K - Ki) / 0.2 * (cInelastic[i + 1] - cInelastic[i]) - } else if (K >= 5.0 && K < 9.49) { - i = (40 + (K - 5.0) / 0.5).toInt() - Ki = 5.0 + (i - 40) * 0.5 - CinelK = cInelastic[i] + (K - Ki) / 0.5 * (cInelastic[i + 1] - cInelastic[i]) - } else { - CinelK = 0.0 - } - Dinelreturn = 2.0 * DH * (1.0 + CinelK) - return Dinelreturn - } - - - private val sumexcKvec = doubleArrayOf(0.0, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.5, 1.8, 2.0, 2.5, 3.0, 4.0, 5.0) - private val sumexcfvec = arrayOf( - doubleArrayOf(2.907e-1, 2.845e-1, 2.665e-1, 2.072e-1, 1.389e-1, - 8.238e-2, 4.454e-2, 2.269e-2, 7.789e-3, 2.619e-3, 1.273e-3, 2.218e-4, 4.372e-5, 2.889e-6, 4.247e-7), // B - doubleArrayOf(3.492e-1, 3.367e-1, 3.124e-1, 2.351e-1, 1.507e-1, - 8.406e-2, 4.214e-2, 1.966e-2, 5.799e-3, 1.632e-3, 6.929e-4, 8.082e-5, 9.574e-6, 1.526e-7, 7.058e-9), // C - doubleArrayOf(6.112e-2, 5.945e-2, 5.830e-2, 5.072e-2, 3.821e-2, - 2.579e-2, 1.567e-2, 8.737e-3, 3.305e-3, 1.191e-3, 6.011e-4, 1.132e-4, 2.362e-5, 1.603e-6, 2.215e-7), // Bp - doubleArrayOf(2.066e-2, 2.127e-2, 2.137e-2, 1.928e-2, 1.552e-2, - 1.108e-2, 7.058e-3, 4.069e-3, 1.590e-3, 5.900e-4, 3.046e-4, 6.142e-5, 1.369e-5, 9.650e-7, 1.244e-7), // Bpp - doubleArrayOf(9.405e-2, 9.049e-2, 8.613e-2, 7.301e-2, 5.144e-2, - 3.201e-2, 1.775e-2, 8.952e-3, 2.855e-3, 8.429e-4, 3.655e-4, 4.389e-5, 5.252e-6, 9.010e-8, 7.130e-9), // D - doubleArrayOf(4.273e-2, 3.862e-2, 3.985e-2, 3.362e-2, 2.486e-2, - 1.612e-2, 9.309e-3, 4.856e-3, 1.602e-3, 4.811e-4, 2.096e-4, 2.498e-5, 2.905e-6, 5.077e-8, 6.583e-9), // Dp - doubleArrayOf(0.000e-3, 2.042e-3, 7.439e-3, 2.200e-2, 3.164e-2, - 3.161e-2, 2.486e-2, 1.664e-2, 7.562e-3, 3.044e-3, 1.608e-3, 3.225e-4, 7.120e-5, 6.290e-6, 1.066e-6))// EF - private val sumexcEeV = doubleArrayOf(12.73, 13.20, 14.77, 15.3, 14.93, 15.4, 13.06) - - private fun _sumexc(K: Double): Double { - var n: Int = 0 - var j: Int - var jmin = 0 - val nmax: Int = 6 - var En: Double - var sum: Double = 0.0 - val f = DoubleArray(7) - val x4 = DoubleArray(4) - val f4 = DoubleArray(4) - - // - while (n <= nmax) { - En = sumexcEeV[n] / 27.21 // En is the excitation energy in Hartree atomic units - when { - K >= 5.0 -> f[n] = 0.0 - K in 3.0..4.0 -> f[n] = sumexcfvec[n][12] + (K - 3.0) * (sumexcfvec[n][13] - sumexcfvec[n][12]) - K in 4.0..5.0 -> f[n] = sumexcfvec[n][13] + (K - 4.0) * (sumexcfvec[n][14] - sumexcfvec[n][13]) - else -> { - j = 0 - while (j < 14) { - if (K >= sumexcKvec[j] && K <= sumexcKvec[j + 1]) { - jmin = j - 1 - } - j++ - } - if (jmin < 0) { - jmin = 0 - } - if (jmin > 11) { - jmin = 11 - } - j = 0 - while (j <= 3) { - x4[j] = sumexcKvec[jmin + j] - f4[j] = sumexcfvec[n][jmin + j] - j++ - } - f[n] = lagrange(4, x4, f4, K) - } - } - sum += f[n] / En - n++ - } - return sum - } - - //Caching function for precision - private val sumexcCache = FunctionCache(0.01, this::_sumexc) - - private fun sumexc(K: Double): Double = sumexcCache(K) - - private fun lagrange(n: Int, xn: DoubleArray, fn: DoubleArray, x: Double): Double { - var i: Int - var j: Int = 0 - var f: Double = 0.0 - var aa: Double - var bb: Double - val a = DoubleArray(100) - val b = DoubleArray(100) - while (j < n) { - i = 0 - while (i < n) { - a[i] = x - xn[i] - b[i] = xn[j] - xn[i] - i++ - } - bb = 1.0 - aa = bb - b[j] = aa - a[j] = b[j] - i = 0 - while (i < n) { - aa *= a[i] - bb *= b[i] - i++ - } - f += fn[j] * aa / bb - j++ - } - return f - } - - // This function computes the total elastic cross section of - // electron scatt. on molecular hydrogen. - // See: Liu, Phys. Rev. A35 (1987) 591, - // Trajmar, Phys Reports 97 (1983) 221. - // E: incident electron energy in eV - // sigmael: cross section in m^2 - private val sigmaelE = doubleArrayOf(0.0, 1.5, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 60.0, 100.0, 150.0, 200.0, 300.0, 400.0) - private val sigmaelS = doubleArrayOf(9.6, 13.0, 15.0, 12.0, 10.0, 7.0, 5.6, 3.3, 1.1, 0.9, 0.5, 0.36, 0.23, 0.15) - - fun sigmael(E: Double): Double { - val gam: Double - var sigma = 0.0 - val T: Double = E / 27.2 - var i: Int - if (E >= 400.0) { - gam = (emass + T) / emass - sigma = gam * gam * Math.PI / (2.0 * T) * (4.2106 - 1.0 / T) * a02 - } else { - i = 0 - while (i <= 12) { - if (E >= sigmaelE[i] && E < sigmaelE[i + 1]) { - sigma = 1e-20 * (sigmaelS[i] + (sigmaelS[i + 1] - sigmaelS[i]) * (E - sigmaelE[i]) / (sigmaelE[i + 1] - sigmaelE[i])) - } - i++ - } - } - return sigma - } - - fun sigmaexc(E: Double): Double { - // This function computes the electronic excitation cross section of - // electron scatt. on molecular hydrogen. - // E: incident electron energy in eV, - // sigmaexc: cross section in m^2 - return when { - E < 9.8 -> 1e-40 - E in 9.8..250.0 -> sigmaBC(E) + sigmadiss10(E) + sigmadiss15(E) - else -> 4.0 * Math.PI * a02 * R / E * (0.80 * log(E / R) + 0.28) - } - // sigma=sigmainel(E)-sigmaion(E); - } - - fun sigmaion(E: Double): Double { - // This function computes the total ionization cross section of - // electron scatt. on molecular hydrogen. - // E: incident electron energy in eV, - // sigmaion: total ionization cross section of - // e+H2 --> e+e+H2^+ or e+e+H^+ +H - // process in m^2. - // - // E<250 eV: Eq. 5 of J. Chem. Phys. 104 (1996) 2956 - // E>250: sigma_i formula on page 107 in - // Phys. Rev. A7 (1973) 103. - // Good agreement with measured results of - // PR A 54 (1996) 2146, and - // Physica 31 (1965) 94. - // - val B = 15.43 - val U = 15.98 - val sigma: Double - val t: Double - val u: Double - val S: Double - val r: Double - val lnt: Double - if (E < 16.0) { - sigma = 1e-40 - } else if (E >= 16.0 && E <= 250.0) { - t = E / B - u = U / B - r = R / B - S = 4.0 * Math.PI * a02 * 2.0 * r * r - lnt = log(t) - sigma = S / (t + u + 1.0) * (lnt / 2.0 * (1.0 - 1.0 / (t * t)) + 1.0 - 1.0 / t - lnt / (t + 1.0)) - } else { - sigma = 4.0 * Math.PI * a02 * R / E * (0.82 * log(E / R) + 1.3) - } - return sigma - } - - // This function computes the sigmaexc electronic excitation - // cross section to the B and C states, with energy loss - // about 12.5 eV. - // E is incident electron energy in eV, - // sigmaexc in m^2 - private val aB = doubleArrayOf(-4.2935194e2, 5.1122109e2, -2.8481279e2, 8.8310338e1, -1.6659591e1, 1.9579609, -1.4012824e-1, 5.5911348e-3, -9.5370103e-5) - private val aC = doubleArrayOf(-8.1942684e2, 9.8705099e2, -5.3095543e2, 1.5917023e2, -2.9121036e1, 3.3321027, -2.3305961e-1, 9.1191781e-3, -1.5298950e-4) - - private fun sigmaBC(E: Double): Double { - var lnsigma: Double = 0.0 - var lnE: Double = log(E) - var lnEn: Double = 1.0 - val sigmaB: Double - var Emin: Double = 12.5 - var sigma: Double = 0.0 - val sigmaC: Double - var n: Int - if (E < Emin) { - sigmaB = 0.0 - } else { - n = 0 - while (n <= 8) { - lnsigma += aB[n] * lnEn - lnEn *= lnE - n++ - } - sigmaB = exp(lnsigma) - } - sigma += sigmaB - // sigma=0.; - // C state: - Emin = 15.8 - lnE = log(E) - lnEn = 1.0 - lnsigma = 0.0 - if (E < Emin) { - sigmaC = 0.0 - } else { - n = 0 - while (n <= 8) { - lnsigma += aC[n] * lnEn - lnEn *= lnE - n++ - } - sigmaC = exp(lnsigma) - } - sigma += sigmaC - return sigma * 1e-4 - } - - ////////////////////////////////////////////////////////////////// - - // This function computes the sigmadiss10 electronic - // dissociative excitation - // cross section, with energy loss - // about 10 eV. - // E is incident electron energy in eV, - // sigmadiss10 in m^2 - private val sigmadiss10a = doubleArrayOf(-2.297914361e5, 5.303988579e5, -5.316636672e5, 3.022690779e5, -1.066224144e5, 2.389841369e4, -3.324526406e3, 2.624761592e2, -9.006246604) - - private fun sigmadiss10(E: Double): Double { - var lnsigma: Double = 0.0 - val lnE: Double = log(E) - var lnEn: Double = 1.0 - val Emin = 9.8 - var n: Int - // E is in eV - val sigma = if (E < Emin) { - 0.0 - } else { - n = 0 - while (n <= 8) { - lnsigma += sigmadiss10a[n] * lnEn - lnEn *= lnE - n++ - } - exp(lnsigma) - } - return sigma * 1e-4 - } - - ////////////////////////////////////////////////////////////////// - // This function computes the sigmadiss15 electronic - // dissociative excitation - // cross section, with energy loss - // about 15 eV. - // E is incident electron energy in eV, - // sigmadiss15 in m^2 - private val sigmadiss15a = doubleArrayOf(-1.157041752e3, 1.501936271e3, -8.6119387e2, 2.754926257e2, -5.380465012e1, 6.573972423, -4.912318139e-1, 2.054926773e-2, -3.689035889e-4) - - private fun sigmadiss15(E: Double): Double { - var lnsigma: Double = 0.0 - val lnE: Double = log(E) - var lnEn: Double - val Emin: Double = 16.5 - var n: Int - // E is in eV - lnEn = 1.0 - val sigma = if (E < Emin) { - 0.0 - } else { - n = 0 - while (n <= 8) { - lnsigma += sigmadiss15a[n] * lnEn - lnEn *= lnE - n++ - } - exp(lnsigma) - } - return sigma * 1e-4 - } - - /** - * Полное сечение с учетом квазиупругих столкновений - * - * @param E - * @return - */ - fun sigmaTotal(E: Double): Double { - return sigmael(E) + sigmaexc(E) + sigmaion(E) - - } -} +package ru.inr.mass.trapping + +import org.apache.commons.math3.random.RandomGenerator +import org.apache.commons.math3.util.FastMath.* +import org.apache.commons.math3.util.Pair + +/** + * + * The initial code was written by Ferenc Glueck and then rewritten by [Sebastian Voecking](mailto:seb.voeck@uni-muenster.de) + * into C++. To reference the code, don't forget to also include [kasiopea](http://stacks.iop.org/1367-2630/19/i=5/a=053012). + * @author Darksnake + */ +object Scatter { + + var debug = false + + val counter = MultiCounter("Accept-reject calls") + + + private val a02 = 28e-22 // Bohr radius squared + private val clight = 137.0 // velocity of light in atomic units + private val emass = 18780.0 // Electron mass in atomic units + private val R = 13.6 // Ryberg energy in eV + + private fun count(counterName: String) { + if (debug) { + counter.increase(counterName) + } + } + + fun randomel(E: Double, generator: RandomGenerator): Pair { + // This subroutine generates energy loss and polar scatt. angle according to + // electron elastic scattering in molecular hydrogen. + // Input: + // E: incident electron energy in eV. + // Output: + // *Eloss: energy loss in eV + // *theta: change of polar angle in degrees + val H2molmass = 69e6 + val T: Double = E / 27.2 + var c = 1.0 + val b: Double + var G: Double + var a: Double + val gam: Double + var K2: Double + val Gmax: Double = when { + E >= 250.0 -> 1e-19 + E < 250.0 && E >= 150.0 -> 2.5e-19 + else -> 1e-18 + } + var i: Int = 1 + + count("randomel-calls") + gam = 1.0 + T / (clight * clight) // relativistic correction factor + b = 2.0 / (1.0 + gam) / T + while (i < 5000) { + count("randomel") + c = 1.0 + b - b * (2.0 + b) / (b + 2.0 * generator.nextDouble()) + K2 = 2.0 * T * (1.0 + gam) * abs(1.0 - c) // momentum transfer squared + a = (4.0 + K2) * (4.0 + K2) / (gam * gam) + G = a * dElastic(E, c) + if (G > Gmax * generator.nextDouble()) { + break + } + i++ + } + return Pair(2.0 * emass / H2molmass * (1.0 - c) * E, acos(c) * 180.0 / Math.PI) + } + + + // Energy values of the excited electronic states: + // (from Mol. Phys. 41 (1980) 1501, in Hartree atomic units) + private val En = doubleArrayOf(12.73 / 27.2, 13.2 / 27.2, 14.77 / 27.2, 15.3 / 27.2, 14.93 / 27.2, 15.4 / 27.2, 13.06 / 27.2) + // Probability numbers of the electronic states: + // (from testelectron7.c calculation ) + private val p = doubleArrayOf(35.86, 40.05, 6.58, 2.26, 9.61, 4.08, 1.54) + // Energy values of the B vibrational states: + // (from: Phys. Rev. A51 (1995) 3745 , in Hartree atomic units) + private val EB = doubleArrayOf(0.411, 0.417, 0.423, 0.428, 0.434, 0.439, 0.444, 0.449, 0.454, 0.459, 0.464, 0.468, 0.473, 0.477, 0.481, 0.485, 0.489, 0.493, 0.496, 0.500, 0.503, 0.507, 0.510, 0.513, 0.516, 0.519, 0.521, 0.524) + // Energy values of the C vibrational states: + // (from: Phys. Rev. A51 (1995) 3745 , in Hartree atomic units) + private val EC = doubleArrayOf(0.452, 0.462, 0.472, 0.481, 0.490, 0.498, 0.506, 0.513, 0.519, 0.525, 0.530, 0.534, 0.537, 0.539) + // Franck-Condon factors of the B vibrational states: + // (from: Phys. Rev. A51 (1995) 3745 ) + private val pB = doubleArrayOf(4.2e-3, 1.5e-2, 3.0e-2, 4.7e-2, 6.3e-2, 7.3e-2, 7.9e-2, 8.0e-2, 7.8e-2, 7.3e-2, 6.6e-2, 5.8e-2, 5.1e-2, 4.4e-2, 3.7e-2, 3.1e-2, 2.6e-2, 2.2e-2, 1.8e-2, 1.5e-2, 1.3e-2, 1.1e-2, 8.9e-3, 7.4e-3, 6.2e-3, 5.2e-3, 4.3e-3, 3.6e-3) + // Franck-Condon factors of the C vibrational states: + // (from: Phys. Rev. A51 (1995) 3745 ) + private val pC = doubleArrayOf(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) + + private val Ecen = 12.6 / 27.21 + + private val fmax by lazy { + val sum = DoubleArray(1001) + + val T: Double = 20000.0 / 27.2 + val xmin = Ecen * Ecen / (2.0 * T) + val ymin = log(xmin) + val ymax = log(8.0 * T + xmin) + val dy = (ymax - ymin) / 1000.0 + + + var res = 0.0 + var i = 0 + while (i <= 1000) { + val y = ymin + dy * i + val K = exp(y / 2.0) + sum[i] = sumexc(K) + if (sum[i] > res) { + res = sum[i] + } + i++ + } + res *= 1.05 + res + } + + + fun randomexc(E: Double, generator: RandomGenerator): Pair { + // This subroutine generates energy loss and polar scatt. angle according to + // electron excitation scattering in molecular hydrogen. + // Input: + // E: incident electron energy in eV. + // Output: + // *Eloss: energy loss in eV + // *theta: change of polar angle in degrees + + + var c = 0.0 + var K: Double + var y: Double + var fy: Double + var pmax: Double + var D: Double + var j: Int + var n = 0 + var N: Int + var v = 0 + + count("randomexc-calls") + + // + // Scattering angle *theta generation: + // + val T = E / 27.2 + val theta: Double + if (E >= 100.0) { + val xmin = Ecen * Ecen / (2.0 * T) + val ymin = log(xmin) + val ymax = log(8.0 * T + xmin) + // dy = (ymax - ymin) / 1000.; + // Generation of y values with the Neumann acceptance-rejection method: + y = ymin + j = 1 + while (j < 5000) { + count("randomexc1") + y = ymin + (ymax - ymin) * generator.nextDouble() + K = exp(y / 2.0) + fy = sumexc(K) + if (fmax * generator.nextDouble() < fy) { + break + } + j++ + } + // Calculation of c=cos(theta) and theta: + val x = exp(y) + c = 1.0 - (x - xmin) / (4.0 * T) + theta = acos(c) * 180.0 / Math.PI + } else { + val Dmax = when { + E <= 25.0 -> 60.0 + E > 25.0 && E <= 35.0 -> 95.0 + E > 35.0 && E <= 50.0 -> 150.0 + else -> 400.0 + } + j = 1 + + while (j < 5000) { + count("randomexc2") + c = -1.0 + 2.0 * generator.nextDouble() + D = dExcitation(E, c) * 1e22 + if (Dmax * generator.nextDouble() < D) { + break + } + j++ + } + theta = acos(c) * 180.0 / Math.PI + } + // Energy loss *Eloss generation: + + // First we generate the electronic state, using the Neumann + // acceptance-rejection method for discrete distribution: + N = 7 // the number of electronic states in our calculation + pmax = p[1] // the maximum of the p[] values + j = 1 + while (j < 5000) { + count("randomexc3") + n = (N * generator.nextDouble()).toInt() + if (generator.nextDouble() * pmax < p[n]) { + break + } + j++ + }// Bp, Bpp, D, Dp, EF states + // C state; we generate now a vibrational state, + // using the Franck-Condon factors + // the number of C vibrational states in our calculation + // maximum of the pC[] values + // B state; we generate now a vibrational state, + // using the Frank-Condon factors + // the number of B vibrational states in our calculation + // maximum of the pB[] values + when { + n < 0 -> n = 0 + n > 6 -> n = 6 + } + val Eloss: Double + when (n) { + 0 -> { + // B state; we generate now a vibrational state, + // using the Frank-Condon factors + N = 28 // the number of B vibrational states in our calculation + pmax = pB[7] // maximum of the pB[] values + j = 1 + while (j < 5000) { + count("randomexc4") + v = (N * generator.nextDouble()).toInt() + if (generator.nextDouble() * pmax < pB[v]) { + break + } + j++ + } + if (v < 0) { + v = 0 + } + if (v > 27) { + v = 27 + } + Eloss = EB[v] * 27.2 + } + 1 -> { + // C state; we generate now a vibrational state, + // using the Franck-Condon factors + N = 14 // the number of C vibrational states in our calculation + pmax = pC[1] // maximum of the pC[] values + j = 1 + while (j < 5000) { + count("randomexc4") + v = (N * generator.nextDouble()).toInt() + if (generator.nextDouble() * pmax < pC[v]) { + break + } + j++ + } + if (v < 0) { + v = 0 + } + if (v > 13) { + v = 13 + } + Eloss = EC[v] * 27.2 + } + else -> + // Bp, Bpp, D, Dp, EF states + Eloss = En[n] * 27.2 + } + return Pair(Eloss, theta) + } + + fun randomion(E: Double, generator: RandomGenerator): Pair { + // This subroutine generates energy loss and polar scatt. angle according to + // electron ionization scattering in molecular hydrogen. + // Input: + // E: incident electron energy in eV. + // Output: + // *Eloss: energy loss in eV + // *theta: change of polar angle in degrees + // The kinetic energy of the secondary electron is: Eloss-15.4 eV + // + val Ei = 15.45 / 27.21 + var c: Double + val b: Double + var K: Double + val xmin: Double + val ymin: Double + val ymax: Double + val x: Double + var y: Double + val T: Double = E / 27.2 + var G: Double + var Gmax: Double + var q: Double + val h: Double + var F: Double + val Fmin: Double + val Fmax: Double + var Gp: Double + val Elmin: Double + val Elmax: Double = (E + 15.45) / 2.0 / 27.2 + val qmin: Double + val qmax: Double + var El: Double + val wmax: Double + var WcE: Double + var Jstarq: Double + var WcstarE: Double + var w: Double + var D2ion: Double + var j: Int = 1 + var K2: Double + var KK: Double + var fE: Double + var kej: Double + var ki: Double + var kf: Double + var Rex: Double + var arg: Double + var arctg: Double + var i: Int + var st1: Double + var st2: Double + count("randomion-calls") + // + // I. Generation of theta + // ----------------------- + Gmax = 1e-20 + if (E < 200.0) { + Gmax = 2e-20 + } + xmin = Ei * Ei / (2.0 * T) + b = xmin / (4.0 * T) + ymin = log(xmin) + ymax = log(8.0 * T + xmin) + // Generation of y values with the Neumann acceptance-rejection method: + y = ymin + while (j < 5000) { + count("randomion1") + y = ymin + (ymax - ymin) * generator.nextDouble() + K = exp(y / 2.0) + c = 1.0 + b - K * K / (4.0 * T) + G = K * K * (dInelastic(E, c) - dExcitation(E, c)) + if (Gmax * generator.nextDouble() < G) { + break + } + j++ + } + // y --> x --> c --> theta + x = exp(y) + c = 1.0 - (x - xmin) / (4.0 * T) + val theta = acos(c) * 180.0 / Math.PI + // + // II. Generation of Eloss, for fixed theta + // ---------------------------------------- + // + // For E<=100 eV we use subr. gensecelen + // (in this case no correlation between theta and Eloss) + if (E <= 100.0) { + return Pair(15.45 + gensecelen(E, generator), theta) + } + // For theta>=20 the free electron model is used + // (with full correlation between theta and Eloss) + if (theta >= 20.0) { + return Pair(E * (1.0 - c * c), theta) + } + // For E>100 eV and theta<20: analytical first Born approximation + // formula of Bethe for H atom (with modification for H2) + // + // Calc. of wmax: + if (theta >= 0.7) { + wmax = 1.1 + } else if (theta <= 0.7 && theta > 0.2) { + wmax = 2.0 + } else if (theta <= 0.2 && theta > 0.05) { + wmax = 4.0 + } else { + wmax = 8.0 + } + // We generate the q value according to the Jstarq pdf. We have to + // define the qmin and qmax limits for this generation: + K = sqrt(4.0 * T * (1.0 - Ei / (2.0 * T) - sqrt(1.0 - Ei / T) * c)) + Elmin = Ei + qmin = Elmin / K - K / 2.0 + qmax = Elmax / K - K / 2.0 + // + q = qmax + Fmax = 1.0 / 2.0 + 1.0 / Math.PI * (q / (1.0 + q * q) + atan(q)) + q = qmin + Fmin = 1.0 / 2.0 + 1.0 / Math.PI * (q / (1.0 + q * q) + atan(q)) + h = Fmax - Fmin + // Generation of Eloss with the Neumann acceptance-rejection method: + El = 0.0 + j = 1 + while (j < 5000) { + // Generation of q with inverse transform method + // (we use the Newton-Raphson method in order to solve the nonlinear eq. + // for the inversion) : + count("randomion2") + F = Fmin + h * generator.nextDouble() + y = 0.0 + i = 1 + while (i <= 30) { + G = 1.0 / 2.0 + (y + sin(2.0 * y) / 2.0) / Math.PI + Gp = (1.0 + cos(2.0 * y)) / Math.PI + y -= (G - F) / Gp + if (abs(G - F) < 1e-8) { + break + } + i++ + } + q = tan(y) + // We have the q value, so we can define El, and calculate the weight: + El = q * K + K * K / 2.0 + // First Born approximation formula of Bethe for e-H ionization: + KK = K + ki = sqrt(2.0 * T) + kf = sqrt(2.0 * (T - El)) + K2 = 4.0 * T * (1.0 - El / (2.0 * T) - sqrt(1.0 - El / T) * c) + if (K2 < 1e-9) { + K2 = 1e-9 + } + K = sqrt(K2) // momentum transfer + Rex = 1.0 - K * K / (kf * kf) + K2 * K2 / (kf * kf * kf * kf) + kej = sqrt(2.0 * abs(El - Ei) + 1e-8) + st1 = K2 - 2.0 * El + 2.0 + if (abs(st1) < 1e-9) { + st1 = 1e-9 + } + arg = 2.0 * kej / st1 + arctg = if (arg >= 0.0) { + atan(arg) + } else { + atan(arg) + Math.PI + } + st1 = (K + kej) * (K + kej) + 1.0 + st2 = (K - kej) * (K - kej) + 1.0 + fE = 1024.0 * El * (K2 + 2.0 / 3.0 * El) / (st1 * st1 * st1 * st2 * st2 * st2) * exp(-2.0 / kej * arctg) / (1.0 - exp(-2.0 * Math.PI / kej)) + D2ion = 2.0 * kf / ki * Rex / (El * K2) * fE + K = KK + // + WcE = D2ion + Jstarq = 16.0 / (3.0 * Math.PI * (1.0 + q * q) * (1.0 + q * q)) + WcstarE = 4.0 / (K * K * K * K * K) * Jstarq + w = WcE / WcstarE + if (wmax * generator.nextDouble() < w) { + break + } + j++ + } + + return Pair(El * 27.2, theta) + } + + private fun gensecelen(E: Double, generator: RandomGenerator): Double { + // This subroutine generates secondary electron energy W + // from ionization of incident electron energy E, by using + // the Lorentzian of Aseev et al. (Eq. 8). + // E and W in eV. + val Ei = 15.45 + val eps2 = 14.3 + val b = 6.25 + val B: Double = atan((Ei - eps2) / b) + val D: Double + val A: Double + val eps: Double + val a: Double + val u: Double = generator.nextDouble() + val epsmax: Double + + epsmax = (E + Ei) / 2.0 + A = atan((epsmax - eps2) / b) + D = b / (A - B) + a = b / D * (u + D / b * B) + eps = eps2 + b * tan(a) + return eps - Ei + } + + + // This subroutine computes the differential cross section + // dElastic= d sigma/d Omega of elastic electron scattering + // on molecular hydrogen. + // See: Nishimura et al., J. Phys. Soc. Jpn. 54 (1985) 1757. + // Input: E= electron kinetic energy in eV + // c= cos(theta), where theta is the polar scatt. angle + // dElastic: in m^2/steradian + private val elasticCel = doubleArrayOf(-0.512, -0.512, -0.509, -0.505, -0.499, -0.491, -0.476, -0.473, -0.462, -0.452, -0.438, -0.422, -0.406, -0.388, -0.370, -0.352, -0.333, -0.314, -0.296, -0.277, -0.258, -0.239, -0.221, -0.202, -0.185, -0.167, -0.151, -0.135, -0.120, -0.105, -0.092, -0.070, -0.053, -0.039, -0.030, -0.024, -0.019, -0.016, -0.014, -0.013, -0.012, -0.009, -0.008, -0.006, -0.005, -0.004, -0.003, -0.002, -0.002, -0.001) + private val elasticE = doubleArrayOf(0.0, 3.0, 6.0, 12.0, 20.0, 32.0, 55.0, 85.0, 150.0, 250.0) + private val elasticT = doubleArrayOf(0.0, 10.0, 20.0, 30.0, 40.0, 60.0, 80.0, 100.0, 140.0, 180.0) + private val elasticD = arrayOf( + doubleArrayOf(2.9, 2.7, 2.5, 2.1, 1.8, 1.2, 0.9, 1.0, 1.6, 1.9), + doubleArrayOf(4.2, 3.6, 3.1, 2.5, 1.9, 1.1, 0.8, 0.9, 1.3, 1.4), + doubleArrayOf(6.0, 4.4, 3.2, 2.3, 1.8, 1.1, 0.7, 0.54, 0.5, 0.6), + doubleArrayOf(6.0, 4.1, 2.8, 1.9, 1.3, 0.6, 0.3, 0.17, 0.16, 0.23), + doubleArrayOf(4.9, 3.2, 2.0, 1.2, 0.8, 0.3, 0.15, 0.09, 0.05, 0.05), + doubleArrayOf(5.2, 2.5, 1.2, 0.64, 0.36, 0.13, 0.05, 0.03, 0.016, 0.02), + doubleArrayOf(4.0, 1.7, 0.7, 0.3, 0.16, 0.05, 0.02, 0.013, 0.01, 0.01), + doubleArrayOf(2.8, 1.1, 0.4, 0.15, 0.07, 0.02, 0.01, 0.007, 0.004, 0.003), + doubleArrayOf(1.2, 0.53, 0.2, 0.08, 0.03, 0.0074, 0.003, 0.0016, 0.001, 0.0008) + ) + + + private fun dElastic(E: Double, c: Double): Double { + val T: Double = E / 27.2 + var K2: Double + var K: Double + val d: Double + val st1: Double + val st2: Double + val DH: Double + val gam: Double + var Delreturn = 0.0 + val CelK: Double + val Ki: Double + val theta: Double + var i: Int + var j: Int + if (E >= 250.0) { + gam = 1.0 + T / (clight * clight) // relativistic correction factor + K2 = 2.0 * T * (1.0 + gam) * (1.0 - c) + if (K2 < 0.0) { + K2 = 1e-30 + } + K = sqrt(K2) + if (K < 1e-9) { + K = 1e-9 // momentum transfer + } + d = 1.4009 // distance of protons in H2 + st1 = 8.0 + K2 + st2 = 4.0 + K2 + // DH is the diff. cross section for elastic electron scatt. + // on atomic hydrogen within the first Born approximation : + DH = 4.0 * st1 * st1 / (st2 * st2 * st2 * st2) * a02 + // CelK calculation with linear interpolation. + // CelK is the correction of the elastic electron + // scatt. on molecular hydrogen compared to the independent atom + // model. + when { + K < 3.0 -> { + i = (K / 0.1).toInt() + Ki = i * 0.1 + CelK = elasticCel[i] + (K - Ki) / 0.1 * (elasticCel[i + 1] - elasticCel[i]) + } + K >= 3.0 && K < 5.0 -> { + i = (30 + (K - 3.0) / 0.2).toInt() + Ki = 3.0 + (i - 30) * 0.2 + CelK = elasticCel[i] + (K - Ki) / 0.2 * (elasticCel[i + 1] - elasticCel[i]) + } + K >= 5.0 && K < 9.49 -> { + i = (40 + (K - 5.0) / 0.5).toInt() + Ki = 5.0 + (i - 40) * 0.5 + CelK = elasticCel[i] + (K - Ki) / 0.5 * (elasticCel[i + 1] - elasticCel[i]) + } + else -> CelK = 0.0 + } + Delreturn = 2.0 * gam * gam * DH * (1.0 + sin(K * d) / (K * d)) * (1.0 + CelK) + } else { + theta = acos(c) * 180.0 / Math.PI + i = 0 + while (i <= 8) { + if (E >= elasticE[i] && E < elasticE[i + 1]) { + j = 0 + while (j <= 8) { + if (theta >= elasticT[j] && theta < elasticT[j + 1]) { + Delreturn = 1e-20 * (elasticD[i][j] + (elasticD[i][j + 1] - elasticD[i][j]) * (theta - elasticT[j]) / (elasticT[j + 1] - elasticT[j])) + } + j++ + } + } + i++ + } + } + return Delreturn + } + + private val excEE = 12.6 / 27.2 + private val excitationE = doubleArrayOf(0.0, 25.0, 35.0, 50.0, 100.0) + private val excitationT = doubleArrayOf(0.0, 10.0, 20.0, 30.0, 40.0, 60.0, 80.0, 100.0, 180.0) + private val excitationD = arrayOf(doubleArrayOf(60.0, 43.0, 27.0, 18.0, 13.0, 8.0, 6.0, 6.0, 6.0), doubleArrayOf(95.0, 70.0, 21.0, 9.0, 6.0, 3.0, 2.0, 2.0, 2.0), doubleArrayOf(150.0, 120.0, 32.0, 8.0, 3.7, 1.9, 1.2, 0.8, 0.8), doubleArrayOf(400.0, 200.0, 12.0, 2.0, 1.4, 0.7, 0.3, 0.2, 0.2)) + + + private fun dExcitation(E: Double, c: Double): Double { + // This subroutine computes the differential cross section + // dElastic= d sigma/d Omega of excitation electron scattering + // on molecular hydrogen. + // Input: E= electron kinetic energy in eV + // c= cos(theta), where theta is the polar scatt. angle + // dExcitation: in m^2/steradian + var K2: Double + val K: Double + var sigma = 0.0 + val T: Double = E / 27.2 + val theta: Double + var i: Int + var j: Int + // + if (E >= 100.0) { + K2 = 4.0 * T * (1.0 - excEE / (2.0 * T) - sqrt(1.0 - excEE / T) * c) + if (K2 < 1e-9) { + K2 = 1e-9 + } + K = sqrt(K2) // momentum transfer + sigma = 2.0 / K2 * sumexc(K) * a02 + } else if (E <= 10.0) { + sigma = 0.0 + } else { + theta = acos(c) * 180.0 / Math.PI + i = 0 + while (i <= 3) { + if (E >= excitationE[i] && E < excitationE[i + 1]) { + j = 0 + while (j <= 7) { + if (theta >= excitationT[j] && theta < excitationT[j + 1]) { + sigma = 1e-22 * (excitationD[i][j] + (excitationD[i][j + 1] - excitationD[i][j]) * (theta - excitationT[j]) / (excitationT[j + 1] - excitationT[j])) + } + j++ + } + } + i++ + } + } + return sigma + } + + // This subroutine computes the differential cross section + // dInelastic= d sigma/d Omega of inelastic electron scattering + // on molecular hydrogen, within the first Born approximation. + // Input: E= electron kinetic energy in eV + // c= cos(theta), where theta is the polar scatt. angle + // dInelastic: in m2/steradian + private val cInelastic = doubleArrayOf(-0.246, -0.244, -0.239, -0.234, -0.227, -0.219, -0.211, -0.201, -0.190, -0.179, -0.167, -0.155, -0.142, -0.130, -0.118, -0.107, -0.096, -0.085, -0.076, -0.067, -0.059, -0.051, -0.045, -0.039, -0.034, -0.029, -0.025, -0.022, -0.019, -0.016, -0.014, -0.010, -0.008, -0.006, -0.004, -0.003, -0.003, -0.002, -0.002, -0.001, -0.001, -0.001, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000) + private val inelasticEi = 0.568 + + + private fun dInelastic(E: Double, c: Double): Double { + val T: Double = E / 27.2 + var K2: Double + val K: Double + val st1: Double + val F: Double + val DH: Double + val Dinelreturn: Double + val CinelK: Double + val Ki: Double + val i: Int + if (E < 16.0) { + return dExcitation(E, c) + } + K2 = 4.0 * T * (1.0 - inelasticEi / (2.0 * T) - sqrt(1.0 - inelasticEi / T) * c) + if (K2 < 1e-9) { + K2 = 1e-9 + } + K = sqrt(K2) // momentum transfer + st1 = 1.0 + K2 / 4.0 + F = 1.0 / (st1 * st1) // scatt. formfactor of hydrogen atom + // DH is the diff. cross section for inelastic electron scatt. + // on atomic hydrogen within the first Born approximation : + DH = 4.0 / (K2 * K2) * (1.0 - F * F) * a02 + // CinelK calculation with linear interpolation. + // CinelK is the correction of the inelastic electron + // scatt. on molecular hydrogen compared to the independent atom + // model. + if (K < 3.0) { + i = (K / 0.1).toInt() + Ki = i * 0.1 + CinelK = cInelastic[i] + (K - Ki) / 0.1 * (cInelastic[i + 1] - cInelastic[i]) + } else if (K >= 3.0 && K < 5.0) { + i = (30 + (K - 3.0) / 0.2).toInt() + Ki = 3.0 + (i - 30) * 0.2 + CinelK = cInelastic[i] + (K - Ki) / 0.2 * (cInelastic[i + 1] - cInelastic[i]) + } else if (K >= 5.0 && K < 9.49) { + i = (40 + (K - 5.0) / 0.5).toInt() + Ki = 5.0 + (i - 40) * 0.5 + CinelK = cInelastic[i] + (K - Ki) / 0.5 * (cInelastic[i + 1] - cInelastic[i]) + } else { + CinelK = 0.0 + } + Dinelreturn = 2.0 * DH * (1.0 + CinelK) + return Dinelreturn + } + + + private val sumexcKvec = doubleArrayOf(0.0, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.5, 1.8, 2.0, 2.5, 3.0, 4.0, 5.0) + private val sumexcfvec = arrayOf( + doubleArrayOf(2.907e-1, 2.845e-1, 2.665e-1, 2.072e-1, 1.389e-1, + 8.238e-2, 4.454e-2, 2.269e-2, 7.789e-3, 2.619e-3, 1.273e-3, 2.218e-4, 4.372e-5, 2.889e-6, 4.247e-7), // B + doubleArrayOf(3.492e-1, 3.367e-1, 3.124e-1, 2.351e-1, 1.507e-1, + 8.406e-2, 4.214e-2, 1.966e-2, 5.799e-3, 1.632e-3, 6.929e-4, 8.082e-5, 9.574e-6, 1.526e-7, 7.058e-9), // C + doubleArrayOf(6.112e-2, 5.945e-2, 5.830e-2, 5.072e-2, 3.821e-2, + 2.579e-2, 1.567e-2, 8.737e-3, 3.305e-3, 1.191e-3, 6.011e-4, 1.132e-4, 2.362e-5, 1.603e-6, 2.215e-7), // Bp + doubleArrayOf(2.066e-2, 2.127e-2, 2.137e-2, 1.928e-2, 1.552e-2, + 1.108e-2, 7.058e-3, 4.069e-3, 1.590e-3, 5.900e-4, 3.046e-4, 6.142e-5, 1.369e-5, 9.650e-7, 1.244e-7), // Bpp + doubleArrayOf(9.405e-2, 9.049e-2, 8.613e-2, 7.301e-2, 5.144e-2, + 3.201e-2, 1.775e-2, 8.952e-3, 2.855e-3, 8.429e-4, 3.655e-4, 4.389e-5, 5.252e-6, 9.010e-8, 7.130e-9), // D + doubleArrayOf(4.273e-2, 3.862e-2, 3.985e-2, 3.362e-2, 2.486e-2, + 1.612e-2, 9.309e-3, 4.856e-3, 1.602e-3, 4.811e-4, 2.096e-4, 2.498e-5, 2.905e-6, 5.077e-8, 6.583e-9), // Dp + doubleArrayOf(0.000e-3, 2.042e-3, 7.439e-3, 2.200e-2, 3.164e-2, + 3.161e-2, 2.486e-2, 1.664e-2, 7.562e-3, 3.044e-3, 1.608e-3, 3.225e-4, 7.120e-5, 6.290e-6, 1.066e-6))// EF + private val sumexcEeV = doubleArrayOf(12.73, 13.20, 14.77, 15.3, 14.93, 15.4, 13.06) + + private fun _sumexc(K: Double): Double { + var n: Int = 0 + var j: Int + var jmin = 0 + val nmax: Int = 6 + var En: Double + var sum: Double = 0.0 + val f = DoubleArray(7) + val x4 = DoubleArray(4) + val f4 = DoubleArray(4) + + // + while (n <= nmax) { + En = sumexcEeV[n] / 27.21 // En is the excitation energy in Hartree atomic units + when { + K >= 5.0 -> f[n] = 0.0 + K in 3.0..4.0 -> f[n] = sumexcfvec[n][12] + (K - 3.0) * (sumexcfvec[n][13] - sumexcfvec[n][12]) + K in 4.0..5.0 -> f[n] = sumexcfvec[n][13] + (K - 4.0) * (sumexcfvec[n][14] - sumexcfvec[n][13]) + else -> { + j = 0 + while (j < 14) { + if (K >= sumexcKvec[j] && K <= sumexcKvec[j + 1]) { + jmin = j - 1 + } + j++ + } + if (jmin < 0) { + jmin = 0 + } + if (jmin > 11) { + jmin = 11 + } + j = 0 + while (j <= 3) { + x4[j] = sumexcKvec[jmin + j] + f4[j] = sumexcfvec[n][jmin + j] + j++ + } + f[n] = lagrange(4, x4, f4, K) + } + } + sum += f[n] / En + n++ + } + return sum + } + + //Caching function for precision + private val sumexcCache = FunctionCache(0.01, this::_sumexc) + + private fun sumexc(K: Double): Double = sumexcCache(K) + + private fun lagrange(n: Int, xn: DoubleArray, fn: DoubleArray, x: Double): Double { + var i: Int + var j: Int = 0 + var f: Double = 0.0 + var aa: Double + var bb: Double + val a = DoubleArray(100) + val b = DoubleArray(100) + while (j < n) { + i = 0 + while (i < n) { + a[i] = x - xn[i] + b[i] = xn[j] - xn[i] + i++ + } + bb = 1.0 + aa = bb + b[j] = aa + a[j] = b[j] + i = 0 + while (i < n) { + aa *= a[i] + bb *= b[i] + i++ + } + f += fn[j] * aa / bb + j++ + } + return f + } + + // This function computes the total elastic cross section of + // electron scatt. on molecular hydrogen. + // See: Liu, Phys. Rev. A35 (1987) 591, + // Trajmar, Phys Reports 97 (1983) 221. + // E: incident electron energy in eV + // sigmael: cross section in m^2 + private val sigmaelE = doubleArrayOf(0.0, 1.5, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 60.0, 100.0, 150.0, 200.0, 300.0, 400.0) + private val sigmaelS = doubleArrayOf(9.6, 13.0, 15.0, 12.0, 10.0, 7.0, 5.6, 3.3, 1.1, 0.9, 0.5, 0.36, 0.23, 0.15) + + fun sigmael(E: Double): Double { + val gam: Double + var sigma = 0.0 + val T: Double = E / 27.2 + var i: Int + if (E >= 400.0) { + gam = (emass + T) / emass + sigma = gam * gam * Math.PI / (2.0 * T) * (4.2106 - 1.0 / T) * a02 + } else { + i = 0 + while (i <= 12) { + if (E >= sigmaelE[i] && E < sigmaelE[i + 1]) { + sigma = 1e-20 * (sigmaelS[i] + (sigmaelS[i + 1] - sigmaelS[i]) * (E - sigmaelE[i]) / (sigmaelE[i + 1] - sigmaelE[i])) + } + i++ + } + } + return sigma + } + + fun sigmaexc(E: Double): Double { + // This function computes the electronic excitation cross section of + // electron scatt. on molecular hydrogen. + // E: incident electron energy in eV, + // sigmaexc: cross section in m^2 + return when { + E < 9.8 -> 1e-40 + E in 9.8..250.0 -> sigmaBC(E) + sigmadiss10(E) + sigmadiss15(E) + else -> 4.0 * Math.PI * a02 * R / E * (0.80 * log(E / R) + 0.28) + } + // sigma=sigmainel(E)-sigmaion(E); + } + + fun sigmaion(E: Double): Double { + // This function computes the total ionization cross section of + // electron scatt. on molecular hydrogen. + // E: incident electron energy in eV, + // sigmaion: total ionization cross section of + // e+H2 --> e+e+H2^+ or e+e+H^+ +H + // process in m^2. + // + // E<250 eV: Eq. 5 of J. Chem. Phys. 104 (1996) 2956 + // E>250: sigma_i formula on page 107 in + // Phys. Rev. A7 (1973) 103. + // Good agreement with measured results of + // PR A 54 (1996) 2146, and + // Physica 31 (1965) 94. + // + val B = 15.43 + val U = 15.98 + val sigma: Double + val t: Double + val u: Double + val S: Double + val r: Double + val lnt: Double + if (E < 16.0) { + sigma = 1e-40 + } else if (E in 16.0..250.0) { + t = E / B + u = U / B + r = R / B + S = 4.0 * Math.PI * a02 * 2.0 * r * r + lnt = log(t) + sigma = S / (t + u + 1.0) * (lnt / 2.0 * (1.0 - 1.0 / (t * t)) + 1.0 - 1.0 / t - lnt / (t + 1.0)) + } else { + sigma = 4.0 * Math.PI * a02 * R / E * (0.82 * log(E / R) + 1.3) + } + return sigma + } + + // This function computes the sigmaexc electronic excitation + // cross section to the B and C states, with energy loss + // about 12.5 eV. + // E is incident electron energy in eV, + // sigmaexc in m^2 + private val aB = doubleArrayOf(-4.2935194e2, 5.1122109e2, -2.8481279e2, 8.8310338e1, -1.6659591e1, 1.9579609, -1.4012824e-1, 5.5911348e-3, -9.5370103e-5) + private val aC = doubleArrayOf(-8.1942684e2, 9.8705099e2, -5.3095543e2, 1.5917023e2, -2.9121036e1, 3.3321027, -2.3305961e-1, 9.1191781e-3, -1.5298950e-4) + + private fun sigmaBC(E: Double): Double { + var lnsigma: Double = 0.0 + var lnE: Double = log(E) + var lnEn: Double = 1.0 + val sigmaB: Double + var Emin: Double = 12.5 + var sigma: Double = 0.0 + val sigmaC: Double + var n: Int + if (E < Emin) { + sigmaB = 0.0 + } else { + n = 0 + while (n <= 8) { + lnsigma += aB[n] * lnEn + lnEn *= lnE + n++ + } + sigmaB = exp(lnsigma) + } + sigma += sigmaB + // sigma=0.; + // C state: + Emin = 15.8 + lnE = log(E) + lnEn = 1.0 + lnsigma = 0.0 + if (E < Emin) { + sigmaC = 0.0 + } else { + n = 0 + while (n <= 8) { + lnsigma += aC[n] * lnEn + lnEn *= lnE + n++ + } + sigmaC = exp(lnsigma) + } + sigma += sigmaC + return sigma * 1e-4 + } + + ////////////////////////////////////////////////////////////////// + + // This function computes the sigmadiss10 electronic + // dissociative excitation + // cross section, with energy loss + // about 10 eV. + // E is incident electron energy in eV, + // sigmadiss10 in m^2 + private val sigmadiss10a = doubleArrayOf(-2.297914361e5, 5.303988579e5, -5.316636672e5, 3.022690779e5, -1.066224144e5, 2.389841369e4, -3.324526406e3, 2.624761592e2, -9.006246604) + + private fun sigmadiss10(E: Double): Double { + var lnsigma: Double = 0.0 + val lnE: Double = log(E) + var lnEn: Double = 1.0 + val Emin = 9.8 + var n: Int + // E is in eV + val sigma = if (E < Emin) { + 0.0 + } else { + n = 0 + while (n <= 8) { + lnsigma += sigmadiss10a[n] * lnEn + lnEn *= lnE + n++ + } + exp(lnsigma) + } + return sigma * 1e-4 + } + + ////////////////////////////////////////////////////////////////// + // This function computes the sigmadiss15 electronic + // dissociative excitation + // cross section, with energy loss + // about 15 eV. + // E is incident electron energy in eV, + // sigmadiss15 in m^2 + private val sigmadiss15a = doubleArrayOf(-1.157041752e3, 1.501936271e3, -8.6119387e2, 2.754926257e2, -5.380465012e1, 6.573972423, -4.912318139e-1, 2.054926773e-2, -3.689035889e-4) + + private fun sigmadiss15(E: Double): Double { + var lnsigma: Double = 0.0 + val lnE: Double = log(E) + val Emin: Double = 16.5 + var n: Int + // E is in eV + var lnEn: Double = 1.0 + val sigma = if (E < Emin) { + 0.0 + } else { + n = 0 + while (n <= 8) { + lnsigma += sigmadiss15a[n] * lnEn + lnEn *= lnE + n++ + } + exp(lnsigma) + } + return sigma * 1e-4 + } + + /** + * Полное сечение с учетом квазиупругих столкновений + * + * @param E + * @return + */ + fun sigmaTotal(E: Double): Double { + return sigmael(E) + sigmaexc(E) + sigmaion(E) + + } +} diff --git a/src/main/kotlin/inr/numass/trapping/SimulationManager.kt b/src/main/kotlin/ru/inr/mass/trapping/SimulationManager.kt similarity index 96% rename from src/main/kotlin/inr/numass/trapping/SimulationManager.kt rename to src/main/kotlin/ru/inr/mass/trapping/SimulationManager.kt index 8a62bf1..395e53c 100644 --- a/src/main/kotlin/inr/numass/trapping/SimulationManager.kt +++ b/src/main/kotlin/ru/inr/mass/trapping/SimulationManager.kt @@ -1,197 +1,197 @@ -package inr.numass.trapping - -import org.apache.commons.math3.analysis.interpolation.LinearInterpolator -import org.apache.commons.math3.random.JDKRandomGenerator -import org.apache.commons.rng.UniformRandomProvider -import org.apache.commons.rng.simple.RandomSource -import java.io.File -import java.io.PrintStream -import java.nio.file.Files -import java.nio.file.StandardOpenOption -import java.util.stream.LongStream - -private val seedGenerator = JDKRandomGenerator() - - -/** - * Created by darksnake on 04-Jun-16. - */ -class SimulationManager() { - - var outputDirectory: File = File("./output") - var fileName = "trap[test]" - - var comment = "" - - /** - * A supplier for random generator. Each track has its own generator - */ - var generatorFactory: (Long) -> UniformRandomProvider = { RandomSource.create(RandomSource.MT_64, seedGenerator.nextInt()) } - var reportFilter: (Simulator.SimulationResult) -> Boolean = { it.state == Simulator.EndState.ACCEPTED } - - var initialE = 18000.0 - var eLow: Double = 14000.0 - var thetaTransport: Double = 24.107064 / 180 * Math.PI - var thetaPinch: Double = 19.481097 / 180 * Math.PI - var gasDensity: Double = 5e20// m^-3 - var bSource: Double = 0.6 - var magneticField: ((Double) -> Double)? = null - - /** - * A syntentic property which allows to set lower energy by specifying range instead of energy itself - */ - var range: Double - get() = initialE - eLow - set(value) { - eLow = initialE - value - } - - -// fun setOutputFile(fileName: String) { -// val outputFile = File(fileName) -// if (!outputFile.exists()) { -// outputFile.parentFile.mkdirs() -// outputFile.createNewFile() -// } -// this.output = PrintStream(FileOutputStream(outputFile)) -// } - - - fun withReportFilter(filter: (Simulator.SimulationResult) -> Boolean): SimulationManager { - this.reportFilter = filter - return this - } - - fun setFields(Bsource: Double, Btransport: Double, Bpinch: Double) { - this.bSource = Bsource - this.thetaTransport = Math.asin(Math.sqrt(Bsource / Btransport)) - this.thetaPinch = Math.asin(Math.sqrt(Bsource / Bpinch)) - } - - /** - * Set field map from values - * - * @param z - * @param b - * @return - */ - fun setFieldMap(z: DoubleArray, b: DoubleArray) { - this.magneticField = { LinearInterpolator().interpolate(z, b).value(it) } - } - - /** - * Симулируем пролет num электронов. - * - * @param num - * @return - */ - @Synchronized - fun simulateAll(num: Number): Counter { - if (!outputDirectory.exists()) { - outputDirectory.mkdirs() - } - - - val simulator = Simulator(eLow, thetaTransport, thetaPinch, gasDensity, bSource, magneticField) - - val counter = Counter() - - val header = """ - E_init = $initialE; - E_low = $eLow; - theta_pinch = $thetaPinch; - theta_transport = $thetaTransport; - density = $gasDensity; - """.trimIndent() + "\n\n" + comment - - - val outputPath = outputDirectory.toPath().resolve("$fileName.out") - - PrintStream(Files.newOutputStream(outputPath, StandardOpenOption.CREATE, StandardOpenOption.TRUNCATE_EXISTING, StandardOpenOption.WRITE)).use { output -> - output.println("# " + header.replace("\n", "\n# "))//adding comment symbols - - System.out.printf("%nStarting simulation with initial energy %g and %d electrons.%n%n", initialE, num.toLong()) - output.printf("%s\t%s\t%s\t%s\t%s\t%s\t%s%n", "id", "E", "theta", "theta_start", "colNum", "L", "state") - LongStream.rangeClosed(1, num.toLong()).parallel() - .mapToObj { id -> - val generator = RandomGeneratorBridge(generatorFactory(id)) - val theta = Math.acos(1 - 2 * generator.nextDouble())// from 0 to Pi - val z = (generator.nextDouble() - 0.5) * Simulator.SOURCE_LENGTH - simulator.simulate(id, generator, initialE, theta, z).also { counter.count(it) } - } - .filter(reportFilter) - .forEach { res -> - printOne(output, res) - } - } - - val statisticsPath = outputDirectory.toPath().resolve("$fileName.stat") - PrintStream(Files.newOutputStream(statisticsPath, StandardOpenOption.CREATE, StandardOpenOption.TRUNCATE_EXISTING, StandardOpenOption.WRITE)).use { - it.println(header + "\n") - printStatistics(it, simulator, counter) - - } - return counter - } - - private fun printStatistics(output: PrintStream, simulator: Simulator, counter: Counter) { - output.apply { - println() - println("***RESULT***") - printf("The total number of events is %d.%n%n", counter.count) - printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI) - printf("The transport reflection angle is %g.%n", simulator.thetaTransport * 180 / Math.PI) - printf("The starting energy is %g.%n", initialE) - printf("The lower energy boundary is %g.%n", simulator.eLow) - printf("The source density is %g.%n", simulator.gasDensity) - - println() - printf("The total number of ACCEPTED events is %d.%n", counter.accepted) - printf("The total number of PASS events is %d.%n", counter.pass) - printf("The total number of REJECTED events is %d.%n", counter.rejected) - printf("The total number of LOWENERGY events is %d.%n%n", counter.lowE) - } - - } - - private fun printOne(out: PrintStream, res: Simulator.SimulationResult) { - out.printf("%d\t%g\t%g\t%g\t%d\t%g\t%s%n", res.id, res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.l, res.state.toString()) - out.flush() - } - - - /** - * Statistic counter - */ - class Counter { - internal var accepted = 0 - internal var pass = 0 - internal var rejected = 0 - internal var lowE = 0 - internal var count = 0 - - fun count(res: Simulator.SimulationResult): Simulator.SimulationResult { - synchronized(this) { - count++ - when (res.state) { - Simulator.EndState.ACCEPTED -> accepted++ - Simulator.EndState.LOWENERGY -> lowE++ - Simulator.EndState.PASS -> pass++ - Simulator.EndState.REJECTED -> rejected++ - Simulator.EndState.NONE -> throw IllegalStateException() - } - } - return res - } - - @Synchronized - fun reset() { - count = 0 - accepted = 0 - pass = 0 - rejected = 0 - lowE = 0 - } - - } -} +package ru.inr.mass.trapping + +import org.apache.commons.math3.analysis.interpolation.LinearInterpolator +import org.apache.commons.math3.random.JDKRandomGenerator +import org.apache.commons.rng.UniformRandomProvider +import org.apache.commons.rng.simple.RandomSource +import java.io.File +import java.io.PrintStream +import java.nio.file.Files +import java.nio.file.StandardOpenOption +import java.util.stream.LongStream + +private val seedGenerator = JDKRandomGenerator() + + +/** + * Created by darksnake on 04-Jun-16. + */ +class SimulationManager() { + + var outputDirectory: File = File("./output") + var fileName = "trap[test]" + + var comment = "" + + /** + * A supplier for random generator. Each track has its own generator + */ + var generatorFactory: (Long) -> UniformRandomProvider = { RandomSource.create(RandomSource.MT_64, seedGenerator.nextInt()) } + var reportFilter: (Simulator.SimulationResult) -> Boolean = { it.state == Simulator.EndState.ACCEPTED } + + var initialE = 18000.0 + var eLow: Double = 14000.0 + var thetaTransport: Double = 24.107064 / 180 * Math.PI + var thetaPinch: Double = 19.481097 / 180 * Math.PI + var gasDensity: Double = 5e20// m^-3 + var bSource: Double = 0.6 + var magneticField: ((Double) -> Double)? = null + + /** + * A syntentic property which allows to set lower energy by specifying range instead of energy itself + */ + var range: Double + get() = initialE - eLow + set(value) { + eLow = initialE - value + } + + +// fun setOutputFile(fileName: String) { +// val outputFile = File(fileName) +// if (!outputFile.exists()) { +// outputFile.parentFile.mkdirs() +// outputFile.createNewFile() +// } +// this.output = PrintStream(FileOutputStream(outputFile)) +// } + + + fun withReportFilter(filter: (Simulator.SimulationResult) -> Boolean): SimulationManager { + this.reportFilter = filter + return this + } + + fun setFields(Bsource: Double, Btransport: Double, Bpinch: Double) { + this.bSource = Bsource + this.thetaTransport = Math.asin(Math.sqrt(Bsource / Btransport)) + this.thetaPinch = Math.asin(Math.sqrt(Bsource / Bpinch)) + } + + /** + * Set field map from values + * + * @param z + * @param b + * @return + */ + fun setFieldMap(z: DoubleArray, b: DoubleArray) { + this.magneticField = { LinearInterpolator().interpolate(z, b).value(it) } + } + + /** + * Симулируем пролет num электронов. + * + * @param num + * @return + */ + @Synchronized + fun simulateAll(num: Number): Counter { + if (!outputDirectory.exists()) { + outputDirectory.mkdirs() + } + + + val simulator = Simulator(eLow, thetaTransport, thetaPinch, gasDensity, bSource, magneticField) + + val counter = Counter() + + val header = """ + E_init = $initialE; + E_low = $eLow; + theta_pinch = $thetaPinch; + theta_transport = $thetaTransport; + density = $gasDensity; + """.trimIndent() + "\n\n" + comment + + + val outputPath = outputDirectory.toPath().resolve("$fileName.out") + + PrintStream(Files.newOutputStream(outputPath, StandardOpenOption.CREATE, StandardOpenOption.TRUNCATE_EXISTING, StandardOpenOption.WRITE)).use { output -> + output.println("# " + header.replace("\n", "\n# "))//adding comment symbols + + System.out.printf("%nStarting simulation with initial energy %g and %d electrons.%n%n", initialE, num.toLong()) + output.printf("%s\t%s\t%s\t%s\t%s\t%s\t%s%n", "id", "E", "theta", "theta_start", "colNum", "L", "state") + LongStream.rangeClosed(1, num.toLong()).parallel() + .mapToObj { id -> + val generator = RandomGeneratorBridge(generatorFactory(id)) + val theta = Math.acos(1 - 2 * generator.nextDouble())// from 0 to Pi + val z = (generator.nextDouble() - 0.5) * Simulator.SOURCE_LENGTH + simulator.simulate(id, generator, initialE, theta, z).also { counter.count(it) } + } + .filter(reportFilter) + .forEach { res -> + printOne(output, res) + } + } + + val statisticsPath = outputDirectory.toPath().resolve("$fileName.stat") + PrintStream(Files.newOutputStream(statisticsPath, StandardOpenOption.CREATE, StandardOpenOption.TRUNCATE_EXISTING, StandardOpenOption.WRITE)).use { + it.println(header + "\n") + printStatistics(it, simulator, counter) + + } + return counter + } + + private fun printStatistics(output: PrintStream, simulator: Simulator, counter: Counter) { + output.apply { + println() + println("***RESULT***") + printf("The total number of events is %d.%n%n", counter.count) + printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI) + printf("The transport reflection angle is %g.%n", simulator.thetaTransport * 180 / Math.PI) + printf("The starting energy is %g.%n", initialE) + printf("The lower energy boundary is %g.%n", simulator.eLow) + printf("The source density is %g.%n", simulator.gasDensity) + + println() + printf("The total number of ACCEPTED events is %d.%n", counter.accepted) + printf("The total number of PASS events is %d.%n", counter.pass) + printf("The total number of REJECTED events is %d.%n", counter.rejected) + printf("The total number of LOWENERGY events is %d.%n%n", counter.lowE) + } + + } + + private fun printOne(out: PrintStream, res: Simulator.SimulationResult) { + out.printf("%d\t%g\t%g\t%g\t%d\t%g\t%s%n", res.id, res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.l, res.state.toString()) + out.flush() + } + + + /** + * Statistic counter + */ + class Counter { + internal var accepted = 0 + internal var pass = 0 + internal var rejected = 0 + internal var lowE = 0 + internal var count = 0 + + fun count(res: Simulator.SimulationResult): Simulator.SimulationResult { + synchronized(this) { + count++ + when (res.state) { + Simulator.EndState.ACCEPTED -> accepted++ + Simulator.EndState.LOWENERGY -> lowE++ + Simulator.EndState.PASS -> pass++ + Simulator.EndState.REJECTED -> rejected++ + Simulator.EndState.NONE -> throw IllegalStateException() + } + } + return res + } + + @Synchronized + fun reset() { + count = 0 + accepted = 0 + pass = 0 + rejected = 0 + lowE = 0 + } + + } +} diff --git a/src/main/kotlin/inr/numass/trapping/Simulator.kt b/src/main/kotlin/ru/inr/mass/trapping/Simulator.kt similarity index 97% rename from src/main/kotlin/inr/numass/trapping/Simulator.kt rename to src/main/kotlin/ru/inr/mass/trapping/Simulator.kt index 6ce096b..f248255 100644 --- a/src/main/kotlin/inr/numass/trapping/Simulator.kt +++ b/src/main/kotlin/ru/inr/mass/trapping/Simulator.kt @@ -1,435 +1,437 @@ -/* - * To change this template, choose Tools | Templates - * and open the template in the editor. - */ -package inr.numass.trapping - -import org.apache.commons.math3.distribution.ExponentialDistribution -import org.apache.commons.math3.geometry.euclidean.threed.Rotation -import org.apache.commons.math3.geometry.euclidean.threed.SphericalCoordinates -import org.apache.commons.math3.random.RandomGenerator -import org.apache.commons.math3.util.FastMath.* -import org.apache.commons.math3.util.Pair -import org.apache.commons.math3.util.Precision - -/** - * @author Darksnake - * @property magneticField Longitudal magnetic field distribution - * @property gasDensity gas density in 1/m^3 - */ -class Simulator( - val eLow: Double, - val thetaTransport: Double, - val thetaPinch: Double, - val gasDensity: Double,// m^-3 - val bSource: Double, - val magneticField: ((Double) -> Double)? -) { - - enum class EndState { - - ACCEPTED, //трэппинговый электрон попал в аксептанс - REJECTED, //трэппинговый электрон вылетел через заднюю пробку - LOWENERGY, //потерял слишком много энергии - PASS, //электрон никогда не запирался и прошел напрямую, нужно для нормировки - NONE - } - - - /** - * Perform scattering in the given position - * - * @param pos - * @return - */ - private fun scatter(pos: State): State { - //Вычисляем сечения и нормируем их на полное сечение - var sigmaIon = Scatter.sigmaion(pos.e) - var sigmaEl = Scatter.sigmael(pos.e) - var sigmaExc = Scatter.sigmaexc(pos.e) - val sigmaTotal = sigmaEl + sigmaIon + sigmaExc - sigmaIon /= sigmaTotal - sigmaEl /= sigmaTotal - sigmaExc /= sigmaTotal - - //проверяем нормировку - assert(Precision.equals(sigmaEl + sigmaExc + sigmaIon, 1.0, 1e-2)) - - val alpha = pos.generator.nextDouble() - - val delta: Pair = when { - alpha < sigmaEl -> Scatter.randomel(pos.e, pos.generator)// elastic - alpha > sigmaEl + sigmaExc -> Scatter.randomion(pos.e, pos.generator)//ionization case - else -> Scatter.randomexc(pos.e, pos.generator)//excitation case - } - - //Обновляем значени угла и энергии независимо ни от чего - pos.substractE(delta.first) - //Изменение угла - pos.addTheta(delta.second / 180 * Math.PI) - - return pos - } - - /** - * Calculate distance to the next scattering - * - * @param e - * @return - */ - private fun freePath(generator: RandomGenerator, e: Double): Double { - //FIXME redundant cross-section calculation - //All cross-sections are in m^2 - return ExponentialDistribution(generator, 1.0 / Scatter.sigmaTotal(e) / gasDensity).sample() - } - - /** - * Calculate propagated position before scattering - * - * @param deltaL - * @return z shift and reflection counter - */ - private fun propagate(pos: State, deltaL: Double): State { - // if magnetic field not defined, consider it to be uniform and equal bSource - if (magneticField == null) { - val deltaZ = deltaL * cos(pos.theta) // direction already included in cos(theta) - // double z0 = pos.z; - pos.addZ(deltaZ) - pos.l += deltaL - - // //if we are crossing source boundary, check for end condition - // while (abs(deltaZ + pos.z) > SOURCE_LENGTH / 2d && !pos.isFinished()) { - // - // pos.checkEndState(); - // } - // - // // 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 { - var curL = 0.0 - val sin2 = sin(pos.theta) * sin(pos.theta) - - while (curL <= deltaL - 0.01 && !pos.isFinished) { - //an l step - val delta = min(deltaL - curL, DELTA_L) - - val b = field(pos.z) - - val root = 1 - sin2 * b / bSource - //preliminary reflection - if (root < 0) { - pos.flip() - } - 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(); - // } - // } - // } - - - curL += delta - pos.l += delta - } - return pos - } - } - - /** - * Magnetic field in the z point - * - * @param z - * @return - */ - private fun field(z: Double): Double { - return magneticField?.invoke(z) ?: bSource - } - - /** - * Симулируем один пробег электрона от начального значения и до вылетания из - * иточника или до того момента, как энергия становится меньше eLow. - */ - fun simulate(id: Long, generator: RandomGenerator, initEnergy: Double, initTheta: Double, initZ: Double): SimulationResult { - assert(initEnergy > 0) - assert(initTheta > 0 && initTheta < Math.PI) - assert(abs(initZ) <= SOURCE_LENGTH / 2.0) - - val pos = State(generator, initEnergy, initTheta, initZ) - - while (!pos.isFinished) { - val dl = freePath(generator, pos.e) // path to next scattering - // propagate to next scattering position - propagate(pos, dl) - - if (!pos.isFinished) { - // perform scatter - scatter(pos) - // increase collision number - pos.colNum++ - if (pos.e < eLow) { - //Если энергия стала слишком маленькой - pos.setEndState(EndState.LOWENERGY) - } - } - } - - return SimulationResult(id, pos.endState, pos.e, pos.theta, initTheta, pos.colNum, pos.l) - } - - fun resetDebugCounters() { - debug = true - counter.resetAll() - } - - fun printDebugCounters() { - if (debug) { - counter.print(System.out) - } else { - throw RuntimeException("Debug not initiated") - } - } - - data class SimulationResult(val id: Long, val state: EndState, val E: Double, val theta: Double, val initTheta: Double, val collisionNumber: Int, val l: Double) - - /** - * Current electron position in simulation. Not thread safe! - * - * @property e Current energy - * @property theta Current theta recalculated to the field in the center of the source - * @property z current z. Zero is the center of the source - */ - private inner class State(val generator: RandomGenerator, internal var e: Double, internal var theta: Double, internal var z: Double) { - /** - * Current total path - */ - var l = 0.0 - - /** - * Number of scatterings - */ - var colNum = 0 - - var endState = EndState.NONE - - internal val isForward: Boolean - get() = theta <= Math.PI / 2 - - internal val isFinished: Boolean - get() = this.endState != EndState.NONE - - /** - * @param dE - * @return resulting E - */ - internal fun substractE(dE: Double): Double { - this.e -= dE - return e - } - - internal fun direction(): Double { - return (if (isForward) 1 else -1).toDouble() - } - - internal fun setEndState(state: EndState) { - this.endState = state - } - - /** - * add Z and calculate direction change - * - * @param dZ - * @return - */ - internal fun addZ(dZ: Double): Double { - this.z += dZ - while (abs(this.z) > SOURCE_LENGTH / 2.0 && !isFinished) { - // reflecting from back wall - if (z < 0) { - if (theta >= PI - thetaTransport) { - setEndState(EndState.REJECTED) - } - z = if (isFinished) { - -SOURCE_LENGTH / 2.0 - } else { - // reflecting from rear pinch - -SOURCE_LENGTH - z - } - } else { - if (theta < thetaPinch) { - if (colNum == 0) { - //counting pass electrons - setEndState(EndState.PASS) - } else { - setEndState(EndState.ACCEPTED) - } - } - z = if (isFinished) { - SOURCE_LENGTH / 2.0 - } else { - // reflecting from forward transport magnet - SOURCE_LENGTH - z - } - } - if (!isFinished) { - 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 - */ - fun flip() { - if (theta < 0 || theta > PI) { - throw Error() - } - theta = PI - theta - } - - /** - * Magnetic field in the current point - * - * @return - */ - fun field(): Double { - return this@Simulator.field(z) - } - - /** - * Real theta angle in current point - * - * @return - */ - fun realTheta(): Double { - if (magneticField == null) { - return theta - } else { - var newTheta = asin(min(abs(sin(theta)) * sqrt(field() / bSource), 1.0)) - if (theta > PI / 2) { - newTheta = PI - newTheta - } - - assert(!java.lang.Double.isNaN(newTheta)) - return newTheta - } - } - - - /** - * Сложение вектора с учетом случайно распределения по фи - * - * @param dTheta - * @return resulting angle - */ - fun addTheta(dTheta: Double): Double { - //Генерируем случайный фи - val phi = generator.nextDouble() * 2.0 * Math.PI - - //change to real angles - val realTheta = realTheta() - //Создаем начальный вектор в сферических координатах - val init = SphericalCoordinates(1.0, 0.0, realTheta + dTheta) - // Задаем вращение относительно оси, перпендикулярной исходному вектору - val rotate = SphericalCoordinates(1.0, 0.0, realTheta) - // поворачиваем исходный вектор на dTheta - val rot = Rotation(rotate.cartesian, phi, null) - - val result = rot.applyTo(init.cartesian) - - val newTheta = acos(result.z) - - // //следим чтобы угол был от 0 до Pi - // if (newtheta < 0) { - // newtheta = -newtheta; - // } - // if (newtheta > Math.PI) { - // newtheta = 2 * Math.PI - newtheta; - // } - - //change back to virtual angles - if (magneticField == null) { - theta = newTheta - } else { - theta = asin(sin(newTheta) * sqrt(bSource / field())) - if (newTheta > PI / 2) { - theta = PI - theta - } - } - - assert(!java.lang.Double.isNaN(theta)) - - return theta - } - - } - - companion object { - - const val SOURCE_LENGTH = 3.0 - private const val DELTA_L = 0.1 //step for propagate calculation - } - - -} +/* + * To change this template, choose Tools | Templates + * and open the template in the editor. + */ +package ru.inr.mass.trapping + +import org.apache.commons.math3.distribution.ExponentialDistribution +import org.apache.commons.math3.geometry.euclidean.threed.Rotation +import org.apache.commons.math3.geometry.euclidean.threed.SphericalCoordinates +import org.apache.commons.math3.random.RandomGenerator +import org.apache.commons.math3.util.FastMath.* +import org.apache.commons.math3.util.Pair +import org.apache.commons.math3.util.Precision +import ru.inr.mass.trapping.Scatter.counter +import ru.inr.mass.trapping.Scatter.debug + +/** + * @author Darksnake + * @property magneticField Longitudal magnetic field distribution + * @property gasDensity gas density in 1/m^3 + */ +class Simulator( + val eLow: Double, + val thetaTransport: Double, + val thetaPinch: Double, + val gasDensity: Double,// m^-3 + val bSource: Double, + val magneticField: ((Double) -> Double)? +) { + + enum class EndState { + + ACCEPTED, //трэппинговый электрон попал в аксептанс + REJECTED, //трэппинговый электрон вылетел через заднюю пробку + LOWENERGY, //потерял слишком много энергии + PASS, //электрон никогда не запирался и прошел напрямую, нужно для нормировки + NONE + } + + + /** + * Perform scattering in the given position + * + * @param pos + * @return + */ + private fun scatter(pos: State): State { + //Вычисляем сечения и нормируем их на полное сечение + var sigmaIon = Scatter.sigmaion(pos.e) + var sigmaEl = Scatter.sigmael(pos.e) + var sigmaExc = Scatter.sigmaexc(pos.e) + val sigmaTotal = sigmaEl + sigmaIon + sigmaExc + sigmaIon /= sigmaTotal + sigmaEl /= sigmaTotal + sigmaExc /= sigmaTotal + + //проверяем нормировку + assert(Precision.equals(sigmaEl + sigmaExc + sigmaIon, 1.0, 1e-2)) + + val alpha = pos.generator.nextDouble() + + val delta: Pair = when { + alpha < sigmaEl -> Scatter.randomel(pos.e, pos.generator)// elastic + alpha > sigmaEl + sigmaExc -> Scatter.randomion(pos.e, pos.generator)//ionization case + else -> Scatter.randomexc(pos.e, pos.generator)//excitation case + } + + //Обновляем значени угла и энергии независимо ни от чего + pos.substractE(delta.first) + //Изменение угла + pos.addTheta(delta.second / 180 * Math.PI) + + return pos + } + + /** + * Calculate distance to the next scattering + * + * @param e + * @return + */ + private fun freePath(generator: RandomGenerator, e: Double): Double { + //FIXME redundant cross-section calculation + //All cross-sections are in m^2 + return ExponentialDistribution(generator, 1.0 / Scatter.sigmaTotal(e) / gasDensity).sample() + } + + /** + * Calculate propagated position before scattering + * + * @param deltaL + * @return z shift and reflection counter + */ + private fun propagate(pos: State, deltaL: Double): State { + // if magnetic field not defined, consider it to be uniform and equal bSource + if (magneticField == null) { + val deltaZ = deltaL * cos(pos.theta) // direction already included in cos(theta) + // double z0 = pos.z; + pos.addZ(deltaZ) + pos.l += deltaL + + // //if we are crossing source boundary, check for end condition + // while (abs(deltaZ + pos.z) > SOURCE_LENGTH / 2d && !pos.isFinished()) { + // + // pos.checkEndState(); + // } + // + // // 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 { + var curL = 0.0 + val sin2 = sin(pos.theta) * sin(pos.theta) + + while (curL <= deltaL - 0.01 && !pos.isFinished) { + //an l step + val delta = min(deltaL - curL, DELTA_L) + + val b = field(pos.z) + + val root = 1 - sin2 * b / bSource + //preliminary reflection + if (root < 0) { + pos.flip() + } + 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(); + // } + // } + // } + + + curL += delta + pos.l += delta + } + return pos + } + } + + /** + * Magnetic field in the z point + * + * @param z + * @return + */ + private fun field(z: Double): Double { + return magneticField?.invoke(z) ?: bSource + } + + /** + * Симулируем один пробег электрона от начального значения и до вылетания из + * иточника или до того момента, как энергия становится меньше eLow. + */ + fun simulate(id: Long, generator: RandomGenerator, initEnergy: Double, initTheta: Double, initZ: Double): SimulationResult { + assert(initEnergy > 0) + assert(initTheta > 0 && initTheta < Math.PI) + assert(abs(initZ) <= SOURCE_LENGTH / 2.0) + + val pos = State(generator, initEnergy, initTheta, initZ) + + while (!pos.isFinished) { + val dl = freePath(generator, pos.e) // path to next scattering + // propagate to next scattering position + propagate(pos, dl) + + if (!pos.isFinished) { + // perform scatter + scatter(pos) + // increase collision number + pos.colNum++ + if (pos.e < eLow) { + //Если энергия стала слишком маленькой + pos.setEndState(EndState.LOWENERGY) + } + } + } + + return SimulationResult(id, pos.endState, pos.e, pos.theta, initTheta, pos.colNum, pos.l) + } + + fun resetDebugCounters() { + debug = true + counter.resetAll() + } + + fun printDebugCounters() { + if (debug) { + counter.print(System.out) + } else { + throw RuntimeException("Debug not initiated") + } + } + + data class SimulationResult(val id: Long, val state: EndState, val E: Double, val theta: Double, val initTheta: Double, val collisionNumber: Int, val l: Double) + + /** + * Current electron position in simulation. Not thread safe! + * + * @property e Current energy + * @property theta Current theta recalculated to the field in the center of the source + * @property z current z. Zero is the center of the source + */ + private inner class State(val generator: RandomGenerator, internal var e: Double, internal var theta: Double, internal var z: Double) { + /** + * Current total path + */ + var l = 0.0 + + /** + * Number of scatterings + */ + var colNum = 0 + + var endState = EndState.NONE + + internal val isForward: Boolean + get() = theta <= Math.PI / 2 + + internal val isFinished: Boolean + get() = this.endState != EndState.NONE + + /** + * @param dE + * @return resulting E + */ + internal fun substractE(dE: Double): Double { + this.e -= dE + return e + } + + internal fun direction(): Double { + return (if (isForward) 1 else -1).toDouble() + } + + internal fun setEndState(state: EndState) { + this.endState = state + } + + /** + * add Z and calculate direction change + * + * @param dZ + * @return + */ + internal fun addZ(dZ: Double): Double { + this.z += dZ + while (abs(this.z) > SOURCE_LENGTH / 2.0 && !isFinished) { + // reflecting from back wall + if (z < 0) { + if (theta >= PI - thetaTransport) { + setEndState(EndState.REJECTED) + } + z = if (isFinished) { + -SOURCE_LENGTH / 2.0 + } else { + // reflecting from rear pinch + -SOURCE_LENGTH - z + } + } else { + if (theta < thetaPinch) { + if (colNum == 0) { + //counting pass electrons + setEndState(EndState.PASS) + } else { + setEndState(EndState.ACCEPTED) + } + } + z = if (isFinished) { + SOURCE_LENGTH / 2.0 + } else { + // reflecting from forward transport magnet + SOURCE_LENGTH - z + } + } + if (!isFinished) { + 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 + */ + fun flip() { + if (theta < 0 || theta > PI) { + throw Error() + } + theta = PI - theta + } + + /** + * Magnetic field in the current point + * + * @return + */ + fun field(): Double { + return this@Simulator.field(z) + } + + /** + * Real theta angle in current point + * + * @return + */ + fun realTheta(): Double { + if (magneticField == null) { + return theta + } else { + var newTheta = asin(min(abs(sin(theta)) * sqrt(field() / bSource), 1.0)) + if (theta > PI / 2) { + newTheta = PI - newTheta + } + + assert(!java.lang.Double.isNaN(newTheta)) + return newTheta + } + } + + + /** + * Сложение вектора с учетом случайно распределения по фи + * + * @param dTheta + * @return resulting angle + */ + fun addTheta(dTheta: Double): Double { + //Генерируем случайный фи + val phi = generator.nextDouble() * 2.0 * Math.PI + + //change to real angles + val realTheta = realTheta() + //Создаем начальный вектор в сферических координатах + val init = SphericalCoordinates(1.0, 0.0, realTheta + dTheta) + // Задаем вращение относительно оси, перпендикулярной исходному вектору + val rotate = SphericalCoordinates(1.0, 0.0, realTheta) + // поворачиваем исходный вектор на dTheta + val rot = Rotation(rotate.cartesian, phi, null) + + val result = rot.applyTo(init.cartesian) + + val newTheta = acos(result.z) + + // //следим чтобы угол был от 0 до Pi + // if (newtheta < 0) { + // newtheta = -newtheta; + // } + // if (newtheta > Math.PI) { + // newtheta = 2 * Math.PI - newtheta; + // } + + //change back to virtual angles + if (magneticField == null) { + theta = newTheta + } else { + theta = asin(sin(newTheta) * sqrt(bSource / field())) + if (newTheta > PI / 2) { + theta = PI - theta + } + } + + assert(!java.lang.Double.isNaN(theta)) + + return theta + } + + } + + companion object { + + const val SOURCE_LENGTH = 3.0 + private const val DELTA_L = 0.1 //step for propagate calculation + } + + +} diff --git a/src/main/kotlin/inr/numass/trapping/Trapping.kt b/src/main/kotlin/ru/inr/mass/trapping/Trapping.kt similarity index 91% rename from src/main/kotlin/inr/numass/trapping/Trapping.kt rename to src/main/kotlin/ru/inr/mass/trapping/Trapping.kt index 9ed354d..327617a 100644 --- a/src/main/kotlin/inr/numass/trapping/Trapping.kt +++ b/src/main/kotlin/ru/inr/mass/trapping/Trapping.kt @@ -1,33 +1,33 @@ -package inr.numass.trapping - -import java.time.Duration -import java.time.Instant - - -fun main(args: Array) { - //val z = doubleArrayOf(-1.736, -1.27, -0.754, -0.238, 0.278, 0.794, 1.31, 1.776) - //val b = doubleArrayOf(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(); - - val es = listOf(12000.0, 14000.0, 16000.0, 18000.0) - - val startTime = Instant.now() - System.out.printf("Starting at %s%n%n", startTime.toString()) - - for (e in es) { - SimulationManager().apply { - comment = "Out of the box cross-sections" - fileName = "trap[$e]" - setFields(0.6, 3.7, 7.2) - gasDensity = 1e19 - initialE = e - range = 4000.0 - }.simulateAll(1_000_000) - } - - val finishTime = Instant.now() - System.out.printf("%nFinished at %s%n", finishTime.toString()) - System.out.printf("Calculation took %s%n", Duration.between(startTime, finishTime).toString()) -} - +package ru.inr.mass.trapping + +import java.time.Duration +import java.time.Instant + + +fun main() { + //val z = doubleArrayOf(-1.736, -1.27, -0.754, -0.238, 0.278, 0.794, 1.31, 1.776) + //val b = doubleArrayOf(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(); + + val es = listOf(12000.0, 14000.0, 16000.0, 18000.0) + + val startTime = Instant.now() + System.out.printf("Starting at %s%n%n", startTime.toString()) + + for (e in es) { + SimulationManager().apply { + comment = "Out of the box cross-sections" + fileName = "trap[$e]" + setFields(0.6, 3.7, 7.2) + gasDensity = 1e19 + initialE = e + range = 4000.0 + }.simulateAll(1_000_000) + } + + val finishTime = Instant.now() + System.out.printf("%nFinished at %s%n", finishTime.toString()) + System.out.printf("Calculation took %s%n", Duration.between(startTime, finishTime).toString()) +} +