NPM/trapping
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

default branch written in kotlin

Browse Source
This commit is contained in:
Alexander Nozik 2018-02-17 13:38:33 +03:00
parent ca59a325f3
commit ef77e6a134
12 changed files with 1709 additions and 1813 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.hgignore
View File

@ -3,7 +3,7 @@
\.chg\..*$
\.rej$
\.conflict\~$
.gradle
.gradle/
build/
.idea/*
!gradle-wrapper.jar

47
build.gradle
View File

@ -1,26 +1,15 @@
apply plugin: 'java'
//apply plugin: 'c'
apply plugin: 'idea'
apply plugin: 'application'
mainClassName = "inr.numass.trapping.Trapping"
plugins {
id "org.jetbrains.kotlin.jvm" version "1.2.21"
id "java"
id "application"
}
group = 'inr.numass'
version = 'dev'
description = "Numass trapping simulation"
sourceCompatibility = 1.8
targetCompatibility = 1.8
tasks.withType(JavaCompile) {
options.encoding = 'UTF-8'
}
tasks.withType(JavaExec) {
enableAssertions = true
}
mainClassName = "inr.numass.trapping.TrappingKt"
repositories {
mavenCentral()
@ -28,9 +17,14 @@ repositories {
dependencies {
compile group: 'org.apache.commons', name: 'commons-math3', version: '3.6.1'
compile "org.jetbrains.kotlin:kotlin-stdlib-jdk8"
testCompile group: 'junit', name: 'junit', version: '4.12'
}
compileKotlin {
kotlinOptions.jvmTarget = "1.8"
}
task runTrap(type: JavaExec) {
classpath = sourceSets.main.runtimeClasspath
@ -40,22 +34,3 @@ task runTrap(type: JavaExec) {
args 'D:\\Work\\Numass\\trapping\\trap.out'
}
// native component
//model {
// components {
// scatter(NativeLibrarySpec)
// }
// toolChains {
// visualCpp(VisualCpp) {
// // Specify the installDir if Visual Studio cannot be located
// // installDir "C:/Apps/Microsoft Visual Studio 10.0"
// installDir "C:\\Program Files (x86)\\Microsoft Visual Studio 12.0"
// }
// gcc(Gcc) {
// // Uncomment to use a GCC install that is not in the PATH
// // path "/usr/bin/gcc"
// }
//
// }
//}

104
src/main/java/inr/numass/trapping/MultiCounter.java
View File

@ -1,104 +0,0 @@
/*
* 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 static java.lang.Integer.valueOf;
import java.util.HashMap;
import java.util.Map;
/**
* TODO есть объект MultiDimensionalCounter, исползовать его?
*
* @author Alexander Nozik
* @version $Id: $Id
*/
public class MultiCounter {
private HashMap<String, Integer> counts = new HashMap<>();
String name;
/**
* <p>Constructor for MultiCounter.</p>
*
* @param name a {@link java.lang.String} object.
*/
public MultiCounter(String name) {
this.name = name;
}
/**
* <p>getCount.</p>
*
* @param name a {@link java.lang.String} object.
* @return a int.
*/
public int getCount(String name) {
if (counts.containsKey(name)) {
return counts.get(name);
} else {
return -1;
}
}
/**
* <p>increase.</p>
*
* @param name a {@link java.lang.String} object.
*/
public synchronized void increase(String name) {
if (counts.containsKey(name)) {
Integer count = counts.get(name);
counts.remove(name);
counts.put(name, count + 1);
} else {
counts.put(name, valueOf(1));
}
}
/**
* <p>print.</p>
*
* @param out a {@link java.io.PrintWriter} object.
*/
public void print(PrintStream out) {
out.printf("%nValues for counter %s%n%n", this.name);
for (Map.Entry<String, Integer> entry : counts.entrySet()) {
String keyName = entry.getKey();
Integer value = entry.getValue();
out.printf("%s : %d%n", keyName, value);
}
}
/**
* <p>reset.</p>
*
* @param name a {@link java.lang.String} object.
*/
public synchronized void reset(String name) {
if (counts.containsKey(name)) {
counts.remove(name);
}
}
/**
* <p>resetAll.</p>
*/
public synchronized void resetAll() {
this.counts = new HashMap<>();
}
}

904
src/main/java/inr/numass/trapping/Scatter.java
View File

@ -1,904 +0,0 @@
/*
* written by Sebastian Voecking <seb.voeck@uni-muenster.de>
*
* See scatter.h for details
*
* Included in this file are function from Ferenc Glueck for calculation of
* cross sections.
*/
package inr.numass.trapping;
import static org.apache.commons.math3.util.FastMath.*;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.util.Pair;
/**
*
* @author Darksnake
*/
public class Scatter {
static final double a02 = 28e-22; // Bohr radius squared
static final double clight = 137; // velocity of light in atomic units
static final double emass = 18780; // Electron mass in atomic units
static final double R = 13.6; // Ryberg energy in eV
private final RandomGenerator generator;
boolean debug = false;
MultiCounter counter = new MultiCounter("Accept-reject calls");
private double fmax = 0;
public Scatter(RandomGenerator generator) {
this.generator = generator;
}
private void count(String counterName){
if(debug){
counter.increase(counterName);
}
}
public Pair<Double, Double> randomel(double E) {
// This subroutine generates energy loss and polar scatt. angle according to
// electron elastic scattering in molecular hydrogen.
// Input:
// E: incident electron energy in eV.
// Output:
// *Eloss: energy loss in eV
// *theta: change of polar angle in degrees
double H2molmass = 69.e6;
double T, c = 1, b, G, a, gam, K2, Gmax;
double[] u = new double[3];
int i;
count("randomel-calls");
if (E >= 250.) {
Gmax = 1.e-19;
} else if (E < 250. && E >= 150.) {
Gmax = 2.5e-19;
} else {
Gmax = 1.e-18;
}
T = E / 27.2;
gam = 1. + T / (clight * clight); // relativistic correction factor
b = 2. / (1. + gam) / T;
for (i = 1; i < 5000; i++) {
count("randomel");
c = 1. + b - b * (2. + b) / (b + 2. * generator.nextDouble());
K2 = 2. * T * (1. + gam) * abs(1d - c); // momentum transfer squared
a = (4. + K2) * (4. + K2) / (gam * gam);
G = a * Del(E, c);
if (G > Gmax * generator.nextDouble()) {
break;
}
}
return new Pair<>(2d * emass / H2molmass * (1d - c) * E, acos(c) * 180d / Math.PI);
}
Pair<Double, Double> randomexc(double E) {
// 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
double Ecen = 12.6 / 27.21;
double[] sum = new double[1001];
double T, c = 0., K, xmin, ymin, ymax, x, y, fy, dy, pmax;
double D, Dmax;
int i, j, n = 0, N, v = 0;
// Energy values of the excited electronic states:
// (from Mol. Phys. 41 (1980) 1501, in Hartree atomic units)
double[] En = {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 )
double[] p = {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)
double[] EB = {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)
double[] EC = {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 )
double[] pB = {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 )
double[] pC = {1.2e-1, 1.9e-1, 1.9e-1, 1.5e-1, 1.1e-1, 7.5e-2, 5.0e-2,
3.3e-2, 2.2e-2, 1.4e-2, 9.3e-3, 6.0e-3, 3.7e-3, 1.8e-3};
count("randomexc-calls");
T = 20000. / 27.2;
//
xmin = Ecen * Ecen / (2. * T);
ymin = log(xmin);
ymax = log(8. * T + xmin);
dy = (ymax - ymin) / 1000.;
// Initialization of the sum[] vector, and fmax calculation:
if (fmax == 0) {
synchronized (this) {
for (i = 0; i <= 1000; i++) {
y = ymin + dy * i;
K = exp(y / 2.);
sum[i] = sumexc(K);
if (sum[i] > fmax) {
fmax = sum[i];
}
}
fmax = 1.05 * fmax;
}
}
//
// Scattering angle *theta generation:
//
T = E / 27.2;
double theta;
if (E >= 100.) {
xmin = Ecen * Ecen / (2. * T);
ymin = log(xmin);
ymax = log(8. * T + xmin);
// dy = (ymax - ymin) / 1000.;
// Generation of y values with the Neumann acceptance-rejection method:
y = ymin;
for (j = 1; j < 5000; j++) {
count("randomexc1");
y = ymin + (ymax - ymin) * generator.nextDouble();
K = exp(y / 2.);
fy = sumexc(K);
if (fmax * generator.nextDouble() < fy) {
break;
}
}
// Calculation of c=cos(theta) and theta:
x = exp(y);
c = 1. - (x - xmin) / (4. * T);
theta = acos(c) * 180. / Math.PI;
} else {
if (E <= 25.) {
Dmax = 60.;
} else if (E > 25. && E <= 35.) {
Dmax = 95.;
} else if (E > 35. && E <= 50.) {
Dmax = 150.;
} else {
Dmax = 400.;
}
for (j = 1; j < 5000; j++) {
count("randomexc2");
c = -1. + 2. * generator.nextDouble();
D = Dexc(E, c) * 1.e22;
if (Dmax * generator.nextDouble() < D) {
break;
}
}
theta = acos(c) * 180. / 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
for (j = 1; j < 5000; j++) {
count("randomexc3");
n = (int) (N * generator.nextDouble());
if (generator.nextDouble() * pmax < p[n]) {
break;
}
}
if (n < 0) {
n = 0;
}
if (n > 6) {
n = 6;
}
if (n > 1) {
}
double Eloss;
switch (n) {
case 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
for (j = 1; j < 5000; j++) {
count("randomexc4");
v = (int) (N * generator.nextDouble());
if (generator.nextDouble() * pmax < pB[v]) {
break;
}
}
if (v < 0) {
v = 0;
}
if (v > 27) {
v = 27;
}
Eloss = EB[v] * 27.2;
break;
case 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
for (j = 1; j < 5000; j++) {
count("randomexc4");
v = (int) (N * generator.nextDouble());
if (generator.nextDouble() * pmax < pC[v]) {
break;
}
}
if (v < 0) {
v = 0;
}
if (v > 13) {
v = 13;
}
Eloss = EC[v] * 27.2;
break;
default:
// Bp, Bpp, D, Dp, EF states
Eloss = En[n] * 27.2;
break;
}
return new Pair<>(Eloss, theta);
}
Pair<Double, Double> randomion(double E) {
// 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
//
double Ei = 15.45 / 27.21;
double c, b, K, xmin, ymin, ymax, x, y, T, G, W, Gmax;
double q, h, F, Fmin, Fmax, Gp, Elmin, Elmax, qmin, qmax, El, wmax;
double WcE, Jstarq, WcstarE, w, D2ion;
int j;
double K2, KK, fE, kej, ki, kf, Rex, arg, arctg;
int i;
double st1, st2;
count("randomion-calls");
//
// I. Generation of theta
// -----------------------
Gmax = 1.e-20;
if (E < 200.) {
Gmax = 2.e-20;
}
T = E / 27.2;
xmin = Ei * Ei / (2. * T);
b = xmin / (4. * T);
ymin = log(xmin);
ymax = log(8. * T + xmin);
// Generation of y values with the Neumann acceptance-rejection method:
y = ymin;
for (j = 1; j < 5000; j++) {
count("randomion1");
y = ymin + (ymax - ymin) * generator.nextDouble();
K = exp(y / 2.);
c = 1. + b - K * K / (4. * T);
G = K * K * (Dinel(E, c) - Dexc(E, c));
if (Gmax * generator.nextDouble() < G) {
break;
}
}
// y --> x --> c --> theta
x = exp(y);
c = 1. - (x - xmin) / (4. * T);
double theta = acos(c) * 180. / 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.) {
return new Pair<>(15.45 + gensecelen(E), theta);
}
// For theta>=20 the free electron model is used
// (with full correlation between theta and Eloss)
if (theta >= 20.) {
return new Pair<>(E * (1. - 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.;
} else if (theta <= 0.2 && theta > 0.05) {
wmax = 4.;
} else {
wmax = 8.;
}
// 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. * T * (1. - Ei / (2. * T) - sqrt(1. - Ei / T) * c));
Elmin = Ei;
Elmax = (E + 15.45) / 2. / 27.2;
qmin = Elmin / K - K / 2.;
qmax = Elmax / K - K / 2.;
//
q = qmax;
Fmax = 1. / 2. + 1. / Math.PI * (q / (1. + q * q) + atan(q));
q = qmin;
Fmin = 1. / 2. + 1. / Math.PI * (q / (1. + q * q) + atan(q));
h = Fmax - Fmin;
// Generation of Eloss with the Neumann acceptance-rejection method:
El = 0;
for (j = 1; j < 5000; j++) {
// 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.;
for (i = 1; i <= 30; i++) {
G = 1. / 2. + (y + sin(2. * y) / 2.) / Math.PI;
Gp = (1. + cos(2. * y)) / Math.PI;
y = y - (G - F) / Gp;
if (abs(G - F) < 1.e-8) {
break;
}
}
q = tan(y);
// We have the q value, so we can define El, and calculate the weight:
El = q * K + K * K / 2.;
// First Born approximation formula of Bethe for e-H ionization:
KK = K;
ki = sqrt(2. * T);
kf = sqrt(2. * (T - El));
K2 = 4. * T * (1. - El / (2. * T) - sqrt(1. - El / T) * c);
if (K2 < 1.e-9) {
K2 = 1.e-9;
}
K = sqrt(K2); // momentum transfer
Rex = 1. - K * K / (kf * kf) + K2 * K2 / (kf * kf * kf * kf);
kej = sqrt(2. * abs(El - Ei) + 1.e-8);
st1 = K2 - 2. * El + 2.;
if (abs(st1) < 1.e-9) {
st1 = 1.e-9;
}
arg = 2. * kej / st1;
if (arg >= 0.) {
arctg = atan(arg);
} else {
arctg = atan(arg) + Math.PI;
}
st1 = (K + kej) * (K + kej) + 1.;
st2 = (K - kej) * (K - kej) + 1.;
fE = 1024. * El * (K2 + 2. / 3. * El) / (st1 * st1 * st1 * st2 * st2 * st2)
* exp(-2. / kej * arctg) / (1. - exp(-2. * Math.PI / kej));
D2ion = 2. * kf / ki * Rex / (El * K2) * fE;
K = KK;
//
WcE = D2ion;
Jstarq = 16. / (3. * Math.PI * (1. + q * q) * (1. + q * q));
WcstarE = 4. / (K * K * K * K * K) * Jstarq;
w = WcE / WcstarE;
if (wmax * generator.nextDouble() < w) {
break;
}
}
return new Pair<>(El * 27.2, theta);
}
double gensecelen(double E) {
// 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.
double Ei = 15.45, eps2 = 14.3, b = 6.25;
double B;
double D, A, eps, a, u, epsmax;
int iff = 0;
B = 0;
if (iff == 0) {
B = atan((Ei - eps2) / b);
iff = 1;
}
epsmax = (E + Ei) / 2.;
A = atan((epsmax - eps2) / b);
D = b / (A - B);
u = generator.nextDouble();
a = b / D * (u + D / b * B);
eps = eps2 + b * tan(a);
return eps - Ei;
}
double Del(double E, double c) {
// This subroutine computes the differential cross section
// Del= 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
// Del: in m^2/steradian
double[] Cel = {
-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
};
double[] e = {0., 3., 6., 12., 20., 32., 55., 85., 150., 250.};
double[] t = {0., 10., 20., 30., 40., 60., 80., 100., 140., 180.};
double[][] D = {
{2.9, 2.7, 2.5, 2.1, 1.8, 1.2, 0.9, 1., 1.6, 1.9},
{4.2, 3.6, 3.1, 2.5, 1.9, 1.1, 0.8, 0.9, 1.3, 1.4},
{6., 4.4, 3.2, 2.3, 1.8, 1.1, 0.7, 0.54, 0.5, 0.6},
{6., 4.1, 2.8, 1.9, 1.3, 0.6, 0.3, 0.17, 0.16, 0.23},
{4.9, 3.2, 2., 1.2, 0.8, 0.3, 0.15, 0.09, 0.05, 0.05},
{5.2, 2.5, 1.2, 0.64, 0.36, 0.13, 0.05, 0.03, 0.016, 0.02},
{4., 1.7, 0.7, 0.3, 0.16, 0.05, 0.02, 0.013, 0.01, 0.01},
{2.8, 1.1, 0.4, 0.15, 0.07, 0.02, 0.01, 0.007, 0.004, 0.003},
{1.2, 0.53, 0.2, 0.08, 0.03, 0.0074, 0.003, 0.0016, 0.001, 0.0008}
};
double T, K2, K, d, st1, st2, DH, gam, Delreturn = 0., CelK, Ki, theta;
int i, j;
T = E / 27.2;
if (E >= 250.) {
gam = 1. + T / (clight * clight); // relativistic correction factor
K2 = 2. * T * (1. + gam) * (1. - c);
if (K2 < 0.) {
K2 = 1.e-30;
}
K = sqrt(K2);
if (K < 1.e-9) {
K = 1.e-9; // momentum transfer
}
d = 1.4009; // distance of protons in H2
st1 = 8. + K2;
st2 = 4. + K2;
// DH is the diff. cross section for elastic electron scatt.
// on atomic hydrogen within the first Born approximation :
DH = 4. * 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.
if (K < 3.) {
i = (int) (K / 0.1);
Ki = i * 0.1;
CelK = Cel[i] + (K - Ki) / 0.1 * (Cel[i + 1] - Cel[i]);
} else if (K >= 3. && K < 5.) {
i = (int) (30 + (K - 3.) / 0.2);
Ki = 3. + (i - 30) * 0.2;
CelK = Cel[i] + (K - Ki) / 0.2 * (Cel[i + 1] - Cel[i]);
} else if (K >= 5. && K < 9.49) {
i = (int) (40 + (K - 5.) / 0.5);
Ki = 5. + (i - 40) * 0.5;
CelK = Cel[i] + (K - Ki) / 0.5 * (Cel[i + 1] - Cel[i]);
} else {
CelK = 0.;
}
Delreturn = 2. * gam * gam * DH * (1. + sin(K * d) / (K * d)) * (1. + CelK);
} else {
theta = acos(c) * 180. / Math.PI;
for (i = 0; i <= 8; i++) {
if (E >= e[i] && E < e[i + 1]) {
for (j = 0; j <= 8; j++) {
if (theta >= t[j] && theta < t[j + 1]) {
Delreturn = 1.e-20 * (D[i][j] + (D[i][j + 1] - D[i][j])
* (theta - t[j]) / (t[j + 1] - t[j]));
}
}
}
}
}
return Delreturn;
}
double Dexc(double E, double c) {
// This subroutine computes the differential cross section
// Del= 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
// Dexc: in m^2/steradian
double K2, K, sigma = 0., T, theta;
double EE = 12.6 / 27.2;
double[] e = {0., 25., 35., 50., 100.};
double[] t = {0., 10., 20., 30., 40., 60., 80., 100., 180.};
double[][] D = {
{60., 43., 27., 18., 13., 8., 6., 6., 6.},
{95., 70., 21., 9., 6., 3., 2., 2., 2.,},
{150., 120., 32., 8., 3.7, 1.9, 1.2, 0.8, 0.8},
{400., 200., 12., 2., 1.4, 0.7, 0.3, 0.2, 0.2}
};
int i, j;
//
T = E / 27.2;
if (E >= 100.) {
K2 = 4. * T * (1. - EE / (2. * T) - sqrt(1. - EE / T) * c);
if (K2 < 1.e-9) {
K2 = 1.e-9;
}
K = sqrt(K2); // momentum transfer
sigma = 2. / K2 * sumexc(K) * a02;
} else if (E <= 10.) {
sigma = 0.;
} else {
theta = acos(c) * 180. / Math.PI;
for (i = 0; i <= 3; i++) {
if (E >= e[i] && E < e[i + 1]) {
for (j = 0; j <= 7; j++) {
if (theta >= t[j] && theta < t[j + 1]) {
sigma = 1.e-22 * (D[i][j] + (D[i][j + 1] - D[i][j])
* (theta - t[j]) / (t[j + 1] - t[j]));
}
}
}
}
}
return sigma;
}
double Dinel(double E, double c) {
// This subroutine computes the differential cross section
// Dinel= 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
// Dinel: in m2/steradian
double[] Cinel = {
-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
};
double Ei = 0.568;
double T, K2, K, st1, F, DH, Dinelreturn, CinelK, Ki;
int i;
if (E < 16.) {
return Dexc(E, c);
}
T = E / 27.2;
K2 = 4. * T * (1. - Ei / (2. * T) - sqrt(1. - Ei / T) * c);
if (K2 < 1.e-9) {
K2 = 1.e-9;
}
K = sqrt(K2); // momentum transfer
st1 = 1. + K2 / 4.;
F = 1. / (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. / (K2 * K2) * (1. - 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.) {
i = (int) (K / 0.1);
Ki = i * 0.1;
CinelK = Cinel[i] + (K - Ki) / 0.1 * (Cinel[i + 1] - Cinel[i]);
} else if (K >= 3. && K < 5.) {
i = (int) (30 + (K - 3.) / 0.2);
Ki = 3. + (i - 30) * 0.2;
CinelK = Cinel[i] + (K - Ki) / 0.2 * (Cinel[i + 1] - Cinel[i]);
} else if (K >= 5. && K < 9.49) {
i = (int) (40 + (K - 5.) / 0.5);
Ki = 5. + (i - 40) * 0.5;
CinelK = Cinel[i] + (K - Ki) / 0.5 * (Cinel[i + 1] - Cinel[i]);
} else {
CinelK = 0.;
}
Dinelreturn = 2. * DH * (1. + CinelK);
return Dinelreturn;
}
double sumexc(double K) {
double[] Kvec = {0., 0.1, 0.2, 0.4, 0.6, 0.8, 1., 1.2, 1.5, 1.8, 2., 2.5, 3., 4., 5.};
double[][] fvec = {
{2.907e-1, 2.845e-1, 2.665e-1, 2.072e-1, 1.389e-1, // B
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},
{3.492e-1, 3.367e-1, 3.124e-1, 2.351e-1, 1.507e-1, // C
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},
{6.112e-2, 5.945e-2, 5.830e-2, 5.072e-2, 3.821e-2, // Bp
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},
{2.066e-2, 2.127e-2, 2.137e-2, 1.928e-2, 1.552e-2, // Bpp
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},
{9.405e-2, 9.049e-2, 8.613e-2, 7.301e-2, 5.144e-2, // D
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},
{4.273e-2, 3.862e-2, 3.985e-2, 3.362e-2, 2.486e-2, // Dp
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},
{0.000e-3, 2.042e-3, 7.439e-3, 2.200e-2, 3.164e-2, // EF
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}};
double[] EeV = {12.73, 13.20, 14.77, 15.3, 14.93, 15.4, 13.06};
int n, j, jmin = 0, nmax;
double En, sum;
double[] f = new double[7];
double[] x4 = new double[4];
double[] f4 = new double[4];
//
sum = 0.;
nmax = 6;
for (n = 0; n <= nmax; n++) {
En = EeV[n] / 27.21; // En is the excitation energy in Hartree atomic units
if (K >= 5.) {
f[n] = 0.;
} else if (K >= 3. && K <= 4.) {
f[n] = fvec[n][12] + (K - 3.) * (fvec[n][13] - fvec[n][12]);
} else if (K >= 4. && K <= 5.) {
f[n] = fvec[n][13] + (K - 4.) * (fvec[n][14] - fvec[n][13]);
} else {
for (j = 0; j < 14; j++) {
if (K >= Kvec[j] && K <= Kvec[j + 1]) {
jmin = j - 1;
}
}
if (jmin < 0) {
jmin = 0;
}
if (jmin > 11) {
jmin = 11;
}
for (j = 0; j <= 3; j++) {
x4[j] = Kvec[jmin + j];
f4[j] = fvec[n][jmin + j];
}
f[n] = lagrange(4, x4, f4, K);
}
sum += f[n] / En;
}
return sum;
}
double lagrange(int n, double[] xn, double[] fn, double x) {
int i, j;
double f, aa, bb;
double[] a = new double[100];
double[] b = new double[100];
f = 0.;
for (j = 0; j < n; j++) {
for (i = 0; i < n; i++) {
a[i] = x - xn[i];
b[i] = xn[j] - xn[i];
}
a[j] = b[j] = aa = bb = 1.;
for (i = 0; i < n; i++) {
aa = aa * a[i];
bb = bb * b[i];
}
f += fn[j] * aa / bb;
}
return f;
}
double sigmael(double E) {
// 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
double[] e = {0., 1.5, 5., 7., 10., 15., 20., 30., 60., 100., 150., 200., 300., 400.};
double[] s = {9.6, 13., 15., 12., 10., 7., 5.6, 3.3, 1.1, 0.9, 0.5, 0.36, 0.23, 0.15};
double gam, sigma = 0., T;
int i;
T = E / 27.2;
if (E >= 400.) {
gam = (emass + T) / emass;
sigma = gam * gam * Math.PI / (2. * T) * (4.2106 - 1. / T) * a02;
} else {
for (i = 0; i <= 12; i++) {
if (E >= e[i] && E < e[i + 1]) {
sigma = 1.e-20 * (s[i] + (s[i + 1] - s[i]) * (E - e[i]) / (e[i + 1] - e[i]));
}
}
}
return sigma;
}
double sigmaexc(double E) {
// 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
double sigma;
if (E < 9.8) {
sigma = 1.e-40;
} else if (E >= 9.8 && E <= 250.) {
sigma = sigmaBC(E) + sigmadiss10(E) + sigmadiss15(E);
} else {
sigma = 4. * Math.PI * a02 * R / E * (0.80 * log(E / R) + 0.28);
}
// sigma=sigmainel(E)-sigmaion(E);
return sigma;
}
double sigmaion(double E) {
// 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.
//
double B = 15.43, U = 15.98;
double sigma, t, u, S, r, lnt;
if (E < 16.) {
sigma = 1.e-40;
} else if (E >= 16. && E
<= 250.) {
t = E / B;
u = U / B;
r = R / B;
S = 4. * Math.PI * a02 * 2. * r * r;
lnt = log(t);
sigma = S / (t + u + 1.) * (lnt / 2. * (1. - 1. / (t * t)) + 1. - 1. / t - lnt / (t + 1.));
} else {
sigma = 4. * Math.PI * a02 * R / E * (0.82 * log(E / R) + 1.3);
}
return sigma;
}
double sigmaBC(double E) {
// 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
double[] aB = {-4.2935194e2, 5.1122109e2, -2.8481279e2,
8.8310338e1, -1.6659591e1, 1.9579609,
-1.4012824e-1, 5.5911348e-3, -9.5370103e-5};
double[] aC = {-8.1942684e2, 9.8705099e2, -5.3095543e2,
1.5917023e2, -2.9121036e1, 3.3321027,
-2.3305961e-1, 9.1191781e-3, -1.5298950e-4};
double lnsigma, lnE, lnEn, sigmaB, Emin, sigma, sigmaC;
int n;
sigma = 0.;
Emin = 12.5;
lnE = log(E);
lnEn = 1.;
lnsigma = 0.;
if (E < Emin) {
sigmaB = 0.;
} else {
for (n = 0; n <= 8; n++) {
lnsigma += aB[n] * lnEn;
lnEn = lnEn * lnE;
}
sigmaB = exp(lnsigma);
}
sigma += sigmaB;
// sigma=0.;
// C state:
Emin = 15.8;
lnE = log(E);
lnEn = 1.;
lnsigma = 0.;
if (E < Emin) {
sigmaC = 0.;
} else {
for (n = 0; n <= 8; n++) {
lnsigma += aC[n] * lnEn;
lnEn = lnEn * lnE;
}
sigmaC = exp(lnsigma);
}
sigma += sigmaC;
return sigma * 1.e-4;
}
//////////////////////////////////////////////////////////////////
double sigmadiss10(double E) {
// 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
double[] a = {-2.297914361e5, 5.303988579e5, -5.316636672e5,
3.022690779e5, -1.066224144e5, 2.389841369e4,
-3.324526406e3, 2.624761592e2, -9.006246604};
double lnsigma, lnE, lnEn, Emin, sigma;
int n;
// E is in eV
sigma = 0.;
Emin = 9.8;
lnE = log(E);
lnEn = 1.;
lnsigma = 0.;
if (E < Emin) {
sigma = 0.;
} else {
for (n = 0; n <= 8; n++) {
lnsigma += a[n] * lnEn;
lnEn = lnEn * lnE;
}
sigma = exp(lnsigma);
}
return sigma * 1.e-4;
}
//////////////////////////////////////////////////////////////////
double sigmadiss15(double E) {
// 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
double[] a = {-1.157041752e3, 1.501936271e3, -8.6119387e2,
2.754926257e2, -5.380465012e1, 6.573972423,
-4.912318139e-1, 2.054926773e-2, -3.689035889e-4};
double lnsigma, lnE, lnEn, Emin, sigma;
int n;
// E is in eV
sigma = 0.;
Emin = 16.5;
lnE = log(E);
lnEn = 1.;
lnsigma = 0.;
if (E < Emin) {
sigma = 0.;
} else {
for (n = 0; n <= 8; n++) {
lnsigma += a[n] * lnEn;
lnEn = lnEn * lnE;
}
sigma = exp(lnsigma);
}
return sigma * 1.e-4;
}
/**
* Полное сечение с учетом квазиупругих столкновений
*
* @param E
* @return
*/
public double sigmaTotal(double E) {
return sigmael(E) + sigmaexc(E) + sigmaion(E);
}
}

219
src/main/java/inr/numass/trapping/SimulationManager.java
View File

@ -1,219 +0,0 @@
package inr.numass.trapping;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.interpolation.LinearInterpolator;
import org.apache.commons.math3.random.ISAACRandom;
import org.apache.commons.math3.random.JDKRandomGenerator;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.random.SynchronizedRandomGenerator;
import java.io.*;
import java.util.Map;
import java.util.function.Function;
import java.util.function.Predicate;
import java.util.stream.Stream;
/**
* Created by darksnake on 04-Jun-16.
*/
public class SimulationManager {
RandomGenerator generator = new JDKRandomGenerator();
Simulator simulator = new Simulator();
private double initialE = 18000;
private double range = 4000;
private PrintStream output = System.out;
private PrintStream statisticOutput = System.out;
private Predicate<Simulator.SimulationResult> reportFileter = (res) -> res.state == Simulator.EndState.ACCEPTED;
// public SimulationManager withParameters(double bSource, double bTransport, double bPinch, double initialE, double energyRange) {
// this.simulator = new Simulator(bSource, bTransport, bPinch, initialE - energyRange);
// this.initialE = initialE;
// return this;
// }
public SimulationManager withInitialE(double initialE){
this.initialE = initialE;
simulator.setELow(initialE-range);
return this;
}
public SimulationManager withRange(double range){
this.range = range;
simulator.setELow(initialE-range);
return this;
}
public SimulationManager withFields(double bSource, double bTransport, double bPinch){
this.simulator.setFields(bSource,bTransport,bPinch);
return this;
}
public SimulationManager withOutputFile(String fileName) throws IOException {
File outputFile = new File(fileName);
if (!outputFile.exists()) {
outputFile.getParentFile().mkdirs();
outputFile.createNewFile();
}
this.output = new PrintStream(new FileOutputStream(outputFile));
return this;
}
/**
* Select output for accepted events
*
* @param output
* @return
*/
public SimulationManager withOutput(PrintStream output) {
this.output = output;
return this;
}
/**
* Select output for statistics
*
* @param statisticOutput
* @return
*/
public SimulationManager withStatisticOutput(PrintStream statisticOutput) {
this.statisticOutput = statisticOutput;
return this;
}
public SimulationManager withReportFilter(Predicate<Simulator.SimulationResult> filter) {
this.reportFileter = filter;
return this;
}
/**
* Set field map as function
*
* @param fieldMap
* @return
*/
public SimulationManager withFieldMap(UnivariateFunction fieldMap) {
this.simulator.setFieldFunc(fieldMap);
return this;
}
/**
* Set field map from values
*
* @param z
* @param b
* @return
*/
public SimulationManager withFieldMap(double[] z, double[] b) {
this.simulator.setFieldFunc(new LinearInterpolator().interpolate(z, b));
return this;
}
/**
* Set source density
* PENDING replace by source thickness?
*
* @param density
* @return
*/
public SimulationManager withGasDensity(double density) {
this.simulator.setGasDensity(density);
return this;
}
/**
* Симулируем пролет num электронов.
*
* @param num
* @return
*/
public synchronized Counter simulateAll(int num) {
Counter counter = new Counter();
System.out.printf("%nStarting sumulation with initial energy %g and %d electrons.%n%n", initialE, num);
output.printf("%s\t%s\t%s\t%s\t%s\t%s%n", "E", "theta", "theta_start", "colNum", "L", "state");
Stream.generate(() -> getRandomTheta()).limit(num).parallel()
.forEach((theta) -> {
double initZ = (generator.nextDouble() - 0.5) * Simulator.SOURCE_LENGTH;
Simulator.SimulationResult res = simulator.simulate(initialE, theta, initZ);
if (reportFileter.test(res)) {
if (output != null) {
printOne(output, res);
}
}
counter.count(res);
});
printStatistics(counter);
return counter;
}
private double getRandomTheta() {
double x = generator.nextDouble();
// from 0 to Pi
return Math.acos(1 - 2 * x);
}
private void printStatistics(Counter counter) {
output.println();
output.println("***RESULT***");
output.printf("The total number of events is %d.%n%n", counter.count);
output.printf("The spectrometer acceptance angle is %g.%n", simulator.thetaPinch * 180 / Math.PI);
output.printf("The transport reflection angle is %g.%n", simulator.thetaTransport * 180 / Math.PI);
output.printf("The starting energy is %g.%n", initialE);
output.printf("The lower energy boundary is %g.%n", simulator.eLow);
output.printf("The source density is %g.%n", simulator.gasDensity);
output.println();
output.printf("The total number of ACCEPTED events is %d.%n", counter.accepted);
output.printf("The total number of PASS events is %d.%n", counter.pass);
output.printf("The total number of REJECTED events is %d.%n", counter.rejected);
output.printf("The total number of LOWENERGY events is %d.%n%n", counter.lowE);
}
private void printOne(PrintStream out, Simulator.SimulationResult res) {
out.printf("%g\t%g\t%g\t%d\t%g\t%s%n", res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.l, res.state.toString());
}
/**
* Statistic counter
*/
public static class Counter {
int accepted = 0;
int pass = 0;
int rejected = 0;
int lowE = 0;
int count = 0;
public synchronized Simulator.SimulationResult count(Simulator.SimulationResult res) {
count++;
switch (res.state) {
case ACCEPTED:
accepted++;
break;
case LOWENERGY:
lowE++;
break;
case PASS:
pass++;
break;
case REJECTED:
rejected++;
break;
case NONE:
throw new IllegalStateException();
}
return res;
}
public synchronized void reset() {
count = 0;
accepted = 0;
pass = 0;
rejected = 0;
lowE = 0;
}
}
}

522
src/main/java/inr/numass/trapping/Simulator.java
View File

@ -1,522 +0,0 @@
/*
* To change this template, choose Tools | Templates
* and open the template in the editor.
*/
package inr.numass.trapping;
import org.apache.commons.math3.analysis.UnivariateFunction;
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.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math3.random.JDKRandomGenerator;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.util.Pair;
import org.apache.commons.math3.util.Precision;
import static org.apache.commons.math3.util.FastMath.*;
/**
* @author Darksnake
*/
public class Simulator {
public static final double SOURCE_LENGTH = 3d;
private static final double DELTA_L = 0.1; //step for propagate calculation
private RandomGenerator generator;
private Scatter scatter;
double eLow = 14000d;//default value
double thetaTransport = 24.107064 / 180 * Math.PI;// default value
double thetaPinch = 19.481097 / 180 * Math.PI;// default value
double gasDensity = 5e20;// m^-3
double bSource = 0.6;
private UnivariateFunction magneticField;
public Simulator() {
generator = new JDKRandomGenerator();// Be careful here since it is not synchronized
scatter = new Scatter(generator);
}
public Simulator(double Bsource, double Btransport, double Bpinch, double Elow) {
this();
setFields(Bsource, Btransport, Bpinch);
setELow(Elow);
}
public enum EndState {
ACCEPTED,//трэппинговый электрон попал в аксептанс
REJECTED,//трэппинговый электрон вылетел через заднюю пробку
LOWENERGY,//потерял слишком много энергии
PASS,//электрон никогда не запирался и прошел напрямую, нужно для нормировки
NONE
}
public final void setELow(double eLow) {
this.eLow = eLow;
}
public final void setFields(double Bsource, double Btransport, double Bpinch) {
this.bSource = Bsource;
this.thetaTransport = Math.asin(Math.sqrt(Bsource / Btransport));
this.thetaPinch = Math.asin(Math.sqrt(Bsource / Bpinch));
}
/**
* Set gas density in 1/m^3
*
* @param gasDensity
*/
public void setGasDensity(double gasDensity) {
this.gasDensity = gasDensity;
}
/**
* Longitudal magnetic field distribution
*
* @param magneticField
*/
public void setFieldFunc(UnivariateFunction magneticField) {
this.magneticField = magneticField;
}
/**
* Perform scattering in the given position
*
* @param pos
* @return
*/
private State scatter(State pos) {
//Вычисляем сечения и нормируем их на полное сечение
double sigmaIon = scatter.sigmaion(pos.e);
double sigmaEl = scatter.sigmael(pos.e);
double sigmaExc = scatter.sigmaexc(pos.e);
double sigmaTotal = sigmaEl + sigmaIon + sigmaExc;
sigmaIon /= sigmaTotal;
sigmaEl /= sigmaTotal;
sigmaExc /= sigmaTotal;
//проверяем нормировку
assert Precision.equals(sigmaEl + sigmaExc + sigmaIon, 1, 1e-2);
double alpha = generator.nextDouble();
Pair<Double, Double> delta;
if (alpha > sigmaEl) {
if (alpha > sigmaEl + sigmaExc) {
//ionization case
delta = scatter.randomion(pos.e);
} else {
//excitation case
delta = scatter.randomexc(pos.e);
}
} else {
// elastic
delta = scatter.randomel(pos.e);
}
//Обновляем значени угла и энергии независимо ни от чего
pos.substractE(delta.getFirst());
//Изменение угла
pos.addTheta(delta.getSecond() / 180 * Math.PI);
return pos;
}
/**
* Calculate distance to the next scattering
*
* @param e
* @return
*/
private double freePath(double e) {
//FIXME redundant cross-section calculation
//All cross-sections are in m^2
return new ExponentialDistribution(generator, 1d / scatter.sigmaTotal(e) / gasDensity).sample();
}
/**
* Calculate propagated position before scattering
*
* @param deltaL
* @return z shift and reflection counter
*/
private State propagate(State pos, double deltaL) {
// if magnetic field not defined, consider it to be uniform and equal bSource
if (magneticField == null) {
double 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 {
double curL = 0;
double sin2 = sin(pos.theta) * sin(pos.theta);
while (curL <= deltaL - 0.01 && !pos.isFinished()) {
//an l step
double delta = min(deltaL - curL, DELTA_L);
double b = field(pos.z);
double 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 double field(double z) {
if (magneticField == null) {
return bSource;
} else {
return magneticField.value(z);
}
}
/**
* Симулируем один пробег электрона от начального значения и до вылетания из
* иточника или до того момента, как энергия становится меньше eLow.
*/
public SimulationResult simulate(double initEnergy, double initTheta, double initZ) {
assert initEnergy > 0;
assert initTheta > 0 && initTheta < Math.PI;
assert abs(initZ) <= SOURCE_LENGTH / 2d;
State pos = new State(initEnergy, initTheta, initZ);
while (!pos.isFinished()) {
double dl = freePath(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);
}
}
}
SimulationResult res = new SimulationResult(pos.endState, pos.e, pos.theta, initTheta, pos.colNum, pos.l);
return res;
}
public void resetDebugCounters() {
scatter.debug = true;
scatter.counter.resetAll();
}
public void printDebugCounters() {
if (scatter.debug) {
scatter.counter.print(System.out);
} else {
throw new RuntimeException("Debug not initiated");
}
}
public static class SimulationResult {
public SimulationResult(EndState state, double E, double theta, double initTheta, int collisionNumber, double l) {
this.state = state;
this.E = E;
this.theta = theta;
this.initTheta = initTheta;
this.collisionNumber = collisionNumber;
this.l = l;
}
public EndState state;
public double E;
public double initTheta;
public double theta;
public int collisionNumber;
public double l;
}
/**
* Current electron position in simulation. Not thread safe!
*/
private class State {
public State(double e, double theta, double z) {
this.e = e;
this.theta = theta;
this.z = z;
}
/**
* Current energy
*/
double e;
/**
* Current theta recalculated to the field in the center of the source
*/
double theta;
/**
* Current total path
*/
double l = 0;
/**
* Number of scatterings
*/
int colNum = 0;
/**
* current z. Zero is the center of the source
*/
double z;
EndState endState = EndState.NONE;
/**
* @param dE
* @return resulting E
*/
double substractE(double dE) {
this.e -= dE;
return e;
}
boolean isForward() {
return theta <= Math.PI / 2;
}
double direction() {
return isForward() ? 1 : -1;
}
boolean isFinished() {
return this.endState != EndState.NONE;
}
void setEndState(EndState state) {
this.endState = state;
}
/**
* add Z and calculate direction change
*
* @param dZ
* @return
*/
double addZ(double dZ) {
this.z += dZ;
while (abs(this.z) > SOURCE_LENGTH / 2d && !isFinished()) {
// reflecting from back wall
if (z < 0) {
if (theta >= PI - thetaTransport) {
setEndState(EndState.REJECTED);
}
if (isFinished()) {
z = -SOURCE_LENGTH / 2d;
} else {
// reflecting from rear pinch
z = -SOURCE_LENGTH - z;
}
} else {
if (theta < thetaPinch) {
if (colNum == 0) {
//counting pass electrons
setEndState(EndState.PASS);
} else {
setEndState(EndState.ACCEPTED);
}
}
if (isFinished()) {
z = SOURCE_LENGTH / 2d;
} else {
// reflecting from forward transport magnet
z = 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
*/
void flip() {
if(theta<0||theta>PI){
throw new Error();
}
theta = PI - theta;
}
/**
* Magnetic field in the current point
*
* @return
*/
double field() {
return Simulator.this.field(z);
}
/**
* Real theta angle in current point
*
* @return
*/
double realTheta() {
if (magneticField == null) {
return theta;
} else {
double newTheta = asin(min(abs(sin(theta)) * sqrt(field() / bSource), 1));
if (theta > PI / 2) {
newTheta = PI - newTheta;
}
assert !Double.isNaN(newTheta);
return newTheta;
}
}
/**
* Сложение вектора с учетом случайно распределения по фи
*
* @param dTheta
* @return resulting angle
*/
double addTheta(double dTheta) {
//Генерируем случайный фи
double phi = generator.nextDouble() * 2 * Math.PI;
//change to real angles
double realTheta = realTheta();
//Создаем начальный вектор в сферических координатах
SphericalCoordinates init = new SphericalCoordinates(1, 0, realTheta + dTheta);
// Задаем вращение относительно оси, перпендикулярной исходному вектору
SphericalCoordinates rotate = new SphericalCoordinates(1, 0, realTheta);
// поворачиваем исходный вектор на dTheta
Rotation rot = new Rotation(rotate.getCartesian(), phi, null);
Vector3D result = rot.applyTo(init.getCartesian());
double newtheta = acos(result.getZ());
// //следим чтобы угол был от 0 до Pi
// if (newtheta < 0) {
// newtheta = -newtheta;
// }
// if (newtheta > Math.PI) {
// newtheta = 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 !Double.isNaN(theta);
return theta;
}
}
}

27
src/main/java/inr/numass/trapping/Trapping.java
View File

@ -1,27 +0,0 @@
package inr.numass.trapping;
import java.io.IOException;
import java.time.Duration;
import java.time.Instant;
import java.util.function.Predicate;
public class Trapping {
public static void main(String[] args) throws IOException {
double[] z = {-1.736, -1.27, -0.754, -0.238, 0.278, 0.794, 1.31, 1.776};
double[] b = {3.70754, 0.62786, 0.60474, 0.60325, 0.60333, 0.60503, 0.6285, 3.70478};
// System.out.println("Press any key to start...");
// System.in.read();
Instant startTime = Instant.now();
System.out.printf("Starting at %s%n%n", startTime.toString());
new SimulationManager()
.withOutputFile("D:\\Work\\Numass\\trapping\\test1.out")
.withFields(0.6, 3.7, 7.2)
// .withFieldMap(z, b)
.withGasDensity(1e19) // per m^3
.withReportFilter(res -> true)
.simulateAll((int) 1e6);
Instant finishTime = Instant.now();
System.out.printf("%nFinished at %s%n", finishTime.toString());
System.out.printf("Calculation took %s%n", Duration.between(startTime, finishTime).toString());
}
}

96
src/main/kotlin/inr/numass/trapping/MultiCounter.kt Normal file
View File

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

963
src/main/kotlin/inr/numass/trapping/Scatter.kt Normal file
View File

@ -0,0 +1,963 @@
/*
* written by Sebastian Voecking <seb.voeck@uni-muenster.de>
*
* See scatter.h for details
*
* Included in this file are function from Ferenc Glueck for calculation of
* cross sections.
*/
package inr.numass.trapping
import org.apache.commons.math3.random.JDKRandomGenerator
import org.apache.commons.math3.util.FastMath.*
import org.apache.commons.math3.random.RandomGenerator
import org.apache.commons.math3.util.Pair
//var generator: RandomGenerator = JDKRandomGenerator()
var debug = false
var counter = MultiCounter("Accept-reject calls")
/**
* Using top level object instead of class
* @author Darksnake
*/
object Scatter {
private var fmax = 0.0
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<Double, Double> {
// 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 * Del(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)
}
internal fun randomexc(E: Double, generator: RandomGenerator): Pair<Double, Double> {
// 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
val Ecen = 12.6 / 27.21
val sum = DoubleArray(1001)
var T: Double = 20000.0 / 27.2
var c = 0.0
var K: Double
var xmin: Double
var ymin: Double
var ymax: Double
val x: Double
var y: Double
var fy: Double
val dy: Double
var pmax: Double
var D: Double
val Dmax: Double
var i: Int
var j: Int
var n = 0
var N: Int
var v = 0
// Energy values of the excited electronic states:
// (from Mol. Phys. 41 (1980) 1501, in Hartree atomic units)
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 )
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)
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)
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 )
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 )
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)
count("randomexc-calls")
//
xmin = Ecen * Ecen / (2.0 * T)
ymin = log(xmin)
ymax = log(8.0 * T + xmin)
dy = (ymax - ymin) / 1000.0
// Initialization of the sum[] vector, and fmax calculation:
synchronized(this) {
if (fmax == 0.0) {
i = 0
while (i <= 1000) {
y = ymin + dy * i
K = exp(y / 2.0)
sum[i] = sumexc(K)
if (sum[i] > fmax) {
fmax = sum[i]
}
i++
}
fmax *= 1.05
}
}
//
// Scattering angle *theta generation:
//
T = E / 27.2
val theta: Double
if (E >= 100.0) {
xmin = Ecen * Ecen / (2.0 * T)
ymin = log(xmin)
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:
x = exp(y)
c = 1.0 - (x - xmin) / (4.0 * T)
theta = acos(c) * 180.0 / Math.PI
} else {
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 = Dexc(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++
}
if (n < 0) {
n = 0
}
if (n > 6) {
n = 6
}
if (n > 1) {
}
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)
}
internal fun randomion(E: Double, generator: RandomGenerator): Pair<Double, Double> {
// 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 * (Dinel(E, c) - Dexc(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 = 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)
}
internal 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
}
internal fun Del(E: Double, c: Double): Double {
// This subroutine computes the differential cross section
// Del= 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
// Del: in m^2/steradian
val Cel = 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)
val e = doubleArrayOf(0.0, 3.0, 6.0, 12.0, 20.0, 32.0, 55.0, 85.0, 150.0, 250.0)
val t = doubleArrayOf(0.0, 10.0, 20.0, 30.0, 40.0, 60.0, 80.0, 100.0, 140.0, 180.0)
val D = 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))
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 = Cel[i] + (K - Ki) / 0.1 * (Cel[i + 1] - Cel[i])
}
K >= 3.0 && K < 5.0 -> {
i = (30 + (K - 3.0) / 0.2).toInt()
Ki = 3.0 + (i - 30) * 0.2
CelK = Cel[i] + (K - Ki) / 0.2 * (Cel[i + 1] - Cel[i])
}
K >= 5.0 && K < 9.49 -> {
i = (40 + (K - 5.0) / 0.5).toInt()
Ki = 5.0 + (i - 40) * 0.5
CelK = Cel[i] + (K - Ki) / 0.5 * (Cel[i + 1] - Cel[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 >= e[i] && E < e[i + 1]) {
j = 0
while (j <= 8) {
if (theta >= t[j] && theta < t[j + 1]) {
Delreturn = 1e-20 * (D[i][j] + (D[i][j + 1] - D[i][j]) * (theta - t[j]) / (t[j + 1] - t[j]))
}
j++
}
}
i++
}
}
return Delreturn
}
internal fun Dexc(E: Double, c: Double): Double {
// This subroutine computes the differential cross section
// Del= 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
// Dexc: in m^2/steradian
var K2: Double
val K: Double
var sigma = 0.0
val T: Double = E / 27.2
val theta: Double
val EE = 12.6 / 27.2
val e = doubleArrayOf(0.0, 25.0, 35.0, 50.0, 100.0)
val t = doubleArrayOf(0.0, 10.0, 20.0, 30.0, 40.0, 60.0, 80.0, 100.0, 180.0)
val D = 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))
var i: Int
var j: Int
//
if (E >= 100.0) {
K2 = 4.0 * T * (1.0 - EE / (2.0 * T) - sqrt(1.0 - EE / 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 >= e[i] && E < e[i + 1]) {
j = 0
while (j <= 7) {
if (theta >= t[j] && theta < t[j + 1]) {
sigma = 1e-22 * (D[i][j] + (D[i][j + 1] - D[i][j]) * (theta - t[j]) / (t[j + 1] - t[j]))
}
j++
}
}
i++
}
}
return sigma
}
internal fun Dinel(E: Double, c: Double): Double {
// This subroutine computes the differential cross section
// Dinel= 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
// Dinel: in m2/steradian
val Cinel = 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)
val Ei = 0.568
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 Dexc(E, c)
}
K2 = 4.0 * T * (1.0 - Ei / (2.0 * T) - sqrt(1.0 - Ei / 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 = Cinel[i] + (K - Ki) / 0.1 * (Cinel[i + 1] - Cinel[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 = Cinel[i] + (K - Ki) / 0.2 * (Cinel[i + 1] - Cinel[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 = Cinel[i] + (K - Ki) / 0.5 * (Cinel[i + 1] - Cinel[i])
} else {
CinelK = 0.0
}
Dinelreturn = 2.0 * DH * (1.0 + CinelK)
return Dinelreturn
}
private fun sumexc(K: Double): Double {
val Kvec = 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)
val fvec = arrayOf(doubleArrayOf(2.907e-1, 2.845e-1, 2.665e-1, 2.072e-1, 1.389e-1, // B
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), doubleArrayOf(3.492e-1, 3.367e-1, 3.124e-1, 2.351e-1, 1.507e-1, // C
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), doubleArrayOf(6.112e-2, 5.945e-2, 5.830e-2, 5.072e-2, 3.821e-2, // Bp
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), doubleArrayOf(2.066e-2, 2.127e-2, 2.137e-2, 1.928e-2, 1.552e-2, // Bpp
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), doubleArrayOf(9.405e-2, 9.049e-2, 8.613e-2, 7.301e-2, 5.144e-2, // D
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), doubleArrayOf(4.273e-2, 3.862e-2, 3.985e-2, 3.362e-2, 2.486e-2, // Dp
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), doubleArrayOf(0.000e-3, 2.042e-3, 7.439e-3, 2.200e-2, 3.164e-2, // EF
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))
val EeV = doubleArrayOf(12.73, 13.20, 14.77, 15.3, 14.93, 15.4, 13.06)
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 = EeV[n] / 27.21 // En is the excitation energy in Hartree atomic units
if (K >= 5.0) {
f[n] = 0.0
} else if (K >= 3.0 && K <= 4.0) {
f[n] = fvec[n][12] + (K - 3.0) * (fvec[n][13] - fvec[n][12])
} else if (K >= 4.0 && K <= 5.0) {
f[n] = fvec[n][13] + (K - 4.0) * (fvec[n][14] - fvec[n][13])
} else {
j = 0
while (j < 14) {
if (K >= Kvec[j] && K <= Kvec[j + 1]) {
jmin = j - 1
}
j++
}
if (jmin < 0) {
jmin = 0
}
if (jmin > 11) {
jmin = 11
}
j = 0
while (j <= 3) {
x4[j] = Kvec[jmin + j]
f4[j] = fvec[n][jmin + j]
j++
}
f[n] = lagrange(4, x4, f4, K)
}
sum += f[n] / En
n++
}
return sum
}
internal fun lagrange(n: Int, xn: DoubleArray, fn: DoubleArray, x: Double): Double {
var i: Int
var j: Int
var f: Double = 0.0
var aa: Double
var bb: Double
val a = DoubleArray(100)
val b = DoubleArray(100)
j = 0
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
}
internal fun sigmael(E: Double): Double {
// 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
val e = 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)
val s = 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)
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 >= e[i] && E < e[i + 1]) {
sigma = 1e-20 * (s[i] + (s[i + 1] - s[i]) * (E - e[i]) / (e[i + 1] - e[i]))
}
i++
}
}
return sigma
}
internal 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
val sigma: Double
if (E < 9.8) {
sigma = 1e-40
} else if (E in 9.8..250.0) {
sigma = sigmaBC(E) + sigmadiss10(E) + sigmadiss15(E)
} else {
sigma = 4.0 * Math.PI * a02 * R / E * (0.80 * log(E / R) + 0.28)
}
// sigma=sigmainel(E)-sigmaion(E);
return sigma
}
internal 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
}
internal fun sigmaBC(E: Double): Double {
// 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
val aB = doubleArrayOf(-4.2935194e2, 5.1122109e2, -2.8481279e2, 8.8310338e1, -1.6659591e1, 1.9579609, -1.4012824e-1, 5.5911348e-3, -9.5370103e-5)
val aC = doubleArrayOf(-8.1942684e2, 9.8705099e2, -5.3095543e2, 1.5917023e2, -2.9121036e1, 3.3321027, -2.3305961e-1, 9.1191781e-3, -1.5298950e-4)
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
}
//////////////////////////////////////////////////////////////////
internal fun sigmadiss10(E: Double): Double {
// 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
val a = doubleArrayOf(-2.297914361e5, 5.303988579e5, -5.316636672e5, 3.022690779e5, -1.066224144e5, 2.389841369e4, -3.324526406e3, 2.624761592e2, -9.006246604)
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 += a[n] * lnEn
lnEn *= lnE
n++
}
exp(lnsigma)
}
return sigma * 1e-4
}
//////////////////////////////////////////////////////////////////
internal fun sigmadiss15(E: Double): Double {
// 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
val a = doubleArrayOf(-1.157041752e3, 1.501936271e3, -8.6119387e2, 2.754926257e2, -5.380465012e1, 6.573972423, -4.912318139e-1, 2.054926773e-2, -3.689035889e-4)
var lnsigma: Double = 0.0
val lnE: Double = log(E)
var lnEn: Double
val Emin: Double = 16.5
var sigma: Double
var n: Int
// E is in eV
sigma = 0.0
lnEn = 1.0
if (E < Emin) {
sigma = 0.0
} else {
n = 0
while (n <= 8) {
lnsigma += a[n] * lnEn
lnEn *= lnE
n++
}
sigma = exp(lnsigma)
}
return sigma * 1e-4
}
/**
* Полное сечение с учетом квазиупругих столкновений
*
* @param E
* @return
*/
fun sigmaTotal(E: Double): Double {
return sigmael(E) + sigmaexc(E) + sigmaion(E)
}
}

169
src/main/kotlin/inr/numass/trapping/SimulationManager.kt Normal file
View File

@ -0,0 +1,169 @@
package inr.numass.trapping
import org.apache.commons.math3.analysis.UnivariateFunction
import org.apache.commons.math3.analysis.interpolation.LinearInterpolator
import org.apache.commons.math3.random.JDKRandomGenerator
import org.apache.commons.math3.random.RandomGenerator
import java.io.*
import java.util.stream.Stream
/**
* Created by darksnake on 04-Jun-16.
*/
class SimulationManager {
var generator: RandomGenerator = JDKRandomGenerator()
/**
* output for accepted events
*/
var output: PrintStream = System.out
/**
* output for statistics
*/
var statisticOutput = System.out
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
}
// from 0 to Pi
private fun getRandomTheta(): Double {
val x = generator.nextDouble()
return Math.acos(1 - 2 * x)
}
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: Int): Counter {
val counter = Counter()
val simulator = Simulator(eLow, thetaTransport, thetaPinch, gasDensity, bSource, magneticField, generator)
System.out.printf("%nStarting sumulation with initial energy %g and %d electrons.%n%n", initialE, num)
output.printf("%s\t%s\t%s\t%s\t%s\t%s%n", "E", "theta", "theta_start", "colNum", "L", "state")
Stream.generate { getRandomTheta() }.limit(num.toLong()).parallel()
.forEach { theta ->
val initZ = (generator.nextDouble() - 0.5) * Simulator.SOURCE_LENGTH
val res = simulator.simulate(initialE, theta, initZ)
if (reportFilter(res)) {
printOne(output, res)
}
counter.count(res)
}
printStatistics(statisticOutput, 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("%g\t%g\t%g\t%d\t%g\t%s%n", res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.l, res.state.toString())
}
/**
* 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
}
}
}

443
src/main/kotlin/inr/numass/trapping/Simulator.kt Normal file
View File

@ -0,0 +1,443 @@
/*
* 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)?,
val generator: RandomGenerator) {
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 = generator.nextDouble()
val delta: Pair<Double, Double>
if (alpha > sigmaEl) {
if (alpha > sigmaEl + sigmaExc) {
//ionization case
delta = Scatter.randomion(pos.e, generator)
} else {
//excitation case
delta = Scatter.randomexc(pos.e, generator)
}
} else {
// elastic
delta = Scatter.randomel(pos.e, generator)
}
//Обновляем значени угла и энергии независимо ни от чего
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(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(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(initEnergy, initTheta, initZ)
while (!pos.isFinished) {
val dl = freePath(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(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")
}
}
class SimulationResult(var state: EndState, var E: Double, var theta: Double, var initTheta: Double, var collisionNumber: Int, var 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(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
*/
internal fun flip() {
if (theta < 0 || theta > PI) {
throw Error()
}
theta = PI - theta
}
/**
* Magnetic field in the current point
*
* @return
*/
internal fun field(): Double {
return this@Simulator.field(z)
}
/**
* Real theta angle in current point
*
* @return
*/
internal 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
*/
internal 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 {
val SOURCE_LENGTH = 3.0
private val DELTA_L = 0.1 //step for propagate calculation
}
}

26
src/main/kotlin/inr/numass/trapping/Trapping.kt Normal file
View File

@ -0,0 +1,26 @@
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 startTime = Instant.now()
System.out.printf("Starting at %s%n%n", startTime.toString())
SimulationManager().apply {
setOutputFile("D:\\Work\\Numass\\trapping\\test1.out")
setFields(0.6, 3.7, 7.2)
gasDensity = 1e19
reportFilter = {true}
}.simulateAll(1e5.toInt())
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())
}