Migrate to git. Update versions

This commit is contained in:
Alexander Nozik 2020-10-18 19:33:32 +03:00
parent d1b66778b4
commit 981cb5f821
14 changed files with 1865 additions and 1855 deletions

5
.gitignore vendored Normal file
View File

@ -0,0 +1,5 @@
/.gradle/
/.idea/
/build/
/output/
!gradle-wrapper.jar

View File

@ -1,13 +0,0 @@
\.orig$
\.orig\..*$
\.chg\..*$
\.rej$
\.conflict\~$
.gradle/
output/
.ipynb_checkpoints/
notebooks/images/
build/
.idea/*
!gradle-wrapper.jar

View File

@ -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()

Binary file not shown.

View File

@ -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

53
gradlew vendored
View File

@ -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" "$@"

43
gradlew.bat vendored
View File

@ -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

View File

@ -1,23 +1,23 @@
package inr.numass.trapping
import java.util.*
class FunctionCache(val xPrecision: Double, val function: (Double) -> Double) {
private val values = TreeMap<Double, Double>()
operator fun invoke(x: Double): Double {
val floor: MutableMap.MutableEntry<Double, Double>? = values.floorEntry(x)
val ceil: MutableMap.MutableEntry<Double, Double>? = 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<Double, Double>()
operator fun invoke(x: Double): Double {
val floor: MutableMap.MutableEntry<Double, Double>? = values.floorEntry(x)
val ceil: MutableMap.MutableEntry<Double, Double>? = 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)
}
}
}

View File

@ -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<String, Int>()
/**
*
* 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<String, Int>()
/**
*
* 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()
}
}

View File

@ -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()
}

View File

@ -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
}
}
}

View File

@ -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<Double, Double> = 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<Double, Double> = 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
}
}

View File

@ -1,33 +1,33 @@
package inr.numass.trapping
import java.time.Duration
import java.time.Instant
fun main(args: Array<String>) {
//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())
}