counts = new HashMap<>();
+ String name;
+
+ /**
+ * Constructor for MultiCounter.
+ *
+ * @param name a {@link java.lang.String} object.
+ */
+ public MultiCounter(String name) {
+ this.name = name;
+ }
+
+ /**
+ * getCount.
+ *
+ * @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;
+ }
+ }
+
+ /**
+ * increase.
+ *
+ * @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));
+ }
+
+ }
+
+ /**
+ * print.
+ *
+ * @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 entry : counts.entrySet()) {
+
+ String keyName = entry.getKey();
+ Integer value = entry.getValue();
+ out.printf("%s : %d%n", keyName, value);
+ }
+ }
+
+ /**
+ * reset.
+ *
+ * @param name a {@link java.lang.String} object.
+ */
+ public synchronized void reset(String name) {
+ if (counts.containsKey(name)) {
+ counts.remove(name);
+ }
+ }
+ /**
+ * resetAll.
+ */
+ public synchronized void resetAll() {
+ this.counts = new HashMap<>();
+ }
+}
diff --git a/src/main/java/inr/numass/trapping/Scatter.java b/src/main/java/inr/numass/trapping/Scatter.java
new file mode 100644
index 0000000..4efe192
--- /dev/null
+++ b/src/main/java/inr/numass/trapping/Scatter.java
@@ -0,0 +1,903 @@
+/*
+ * written by Sebastian Voecking
+ *
+ * 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 java.lang.Math.abs;
+import static java.lang.Math.acos;
+import static java.lang.Math.atan;
+import static java.lang.Math.cos;
+import static java.lang.Math.exp;
+import static java.lang.Math.log;
+import static java.lang.Math.sin;
+import static java.lang.Math.sqrt;
+import static java.lang.Math.tan;
+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;
+ MultiCounter counter = new MultiCounter("Accept-reject calls");
+
+ private double fmax = 0;
+
+ public Scatter(RandomGenerator generator) {
+ this.generator = generator;
+ }
+
+ public Pair 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;
+
+ counter.increase("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++) {
+ counter.increase("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 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};
+ counter.increase("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++) {
+ counter.increase("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++) {
+ counter.increase("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++) {
+ counter.increase("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++) {
+ counter.increase("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++) {
+ counter.increase("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 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;
+ counter.increase("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++) {
+ counter.increase("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) :
+ counter.increase("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);
+
+ }
+}
diff --git a/src/main/java/hep/dataforge/trapping/Trapping.java b/src/main/java/inr/numass/trapping/Trapping.java
similarity index 91%
rename from src/main/java/hep/dataforge/trapping/Trapping.java
rename to src/main/java/inr/numass/trapping/Trapping.java
index 70cdc9b..58d1b24 100644
--- a/src/main/java/hep/dataforge/trapping/Trapping.java
+++ b/src/main/java/inr/numass/trapping/Trapping.java
@@ -1,4 +1,4 @@
-package hep.dataforge.trapping;
+package inr.numass.trapping;
import java.io.File;
import java.io.FileNotFoundException;
@@ -32,7 +32,11 @@ public class Trapping {
int rejected = 0;
int lowE = 0;
- List results = simulator.simulateAll(E, 500000);
+ simulator.scatter.counter.resetAll();
+ List results = simulator.simulateAll(E, 10000);
+ simulator.scatter.counter.print(System.out);
+
+ System.out.printf("%nSimulation complete.%n%n");
for (ElectronTrappingSimulator.SimulaionResult res : results) {
if (out != null) {
diff --git a/src/main/resources/csource/constants.h b/src/main/resources/csource/constants.h
new file mode 100644
index 0000000..346eb29
--- /dev/null
+++ b/src/main/resources/csource/constants.h
@@ -0,0 +1,15 @@
+/* $Id$ */
+/*
+ * constants.h
+ *
+ * written by Sebastian Voecking
+ *
+ * Constants which are used through out the scatter package
+ */
+
+#define C 299792458. /* velocity of light in SI units */
+#define CHARGEUNIT 1.602177e-19 /* electron charge (in SI, without sign) */
+#define BOLTZMANN 1.380658e-23 /* Boltzmann constant */
+#define ELECTRON_MASS 9.109389e-31 /* Electron mass */
+
+
diff --git a/src/main/resources/csource/random.c b/src/main/resources/csource/random.c
new file mode 100644
index 0000000..76d77ba
--- /dev/null
+++ b/src/main/resources/csource/random.c
@@ -0,0 +1,184 @@
+/* $Id$ */
+/*
+ * random.c
+ *
+ * written by Sebastian Voecking
+ *
+ * For details see random.h
+ */
+
+#include "random.h"
+
+#include
+
+static int method = RANDOM_STDLIB;
+
+static int rand_seed=-1; // for random calculation
+static double get_random_1_cw();
+static void subrn(double* u, int len);
+static double random_james();
+
+void random_set_method(RandomMethode meth)
+{
+ method = meth;
+}
+
+double random_get()
+{
+ switch(method) {
+ case RANDOM_STDLIB:
+ return ((double)rand())/RAND_MAX;
+ case RANDOM_CW:
+ return get_random_1_cw();
+ case RANDOM_JAMES:
+ return random_james();
+ default:
+ return 0;
+ }
+}
+
+void random_seed(int seed)
+{
+ switch(method) {
+ case RANDOM_STDLIB:
+ srand(seed);
+ break;
+ case RANDOM_CW:
+ rand_seed = seed;
+ break;
+ case RANDOM_JAMES:
+ break;
+ }
+}
+
+double get_random_1_cw()
+ // function gives a random value between 0 and 1.
+ // NO imput value needed.
+ // numerical rec. version provided by Ch. Weinheimer
+{
+#define IA 16807
+#define IM 2147483647
+#define AM (1./IM)
+#define IQ 127773
+#define IR 2836
+#define NTAB 32
+#define NDIV (1+(IM-1)/NTAB)
+#define RNMX (1.-1.2e-7)
+
+ int j, k;
+ static int iy=0, iv[NTAB];
+ double temp;
+
+ if (rand_seed<=0 || !iy){
+ if (-rand_seed<1) rand_seed=1;
+ else rand_seed=-rand_seed;
+ for (j=NTAB+7; j>=0; j--){
+ k=rand_seed/IQ;
+ rand_seed=IA*(rand_seed-k*IQ)-IR*k;
+ if (rand_seed<0) rand_seed+=IM;
+ if (jRNMX) return RNMX;
+ else return temp;
+
+#undef IA
+#undef IM
+#undef AM
+#undef IQ
+#undef IR
+#undef NTAB
+#undef NDIV
+#undef RNMX
+}
+
+void subrn(double *u,int len)
+{
+// This subroutine computes random numbers u[1],...,u[len]
+// in the (0,1) interval. It uses the 0=900000000) ijkl=1;
+ ij=ijkl/30082;
+ kl=ijkl-30082*ij;
+ i=((ij/177)%177)+2;
+ j=(ij%177)+2;
+ k=((kl/169)%178)+1;
+ l=kl%169;
+ for(ii=1;ii<=97;ii++)
+ { s=0; t=0.5;
+ for(jj=1;jj<=24;jj++)
+ { m=(((i*j)%179)*k)%179;
+ i=j; j=k; k=m;
+ l=(53*l+1)%169;
+ if((l*m)%64 >= 32) s=s+t;
+ t=0.5*t;
+ }
+ uu[ii]=s;
+ }
+ c=362436./16777216.;
+ cd=7654321./16777216.;
+ cm=16777213./16777216.;
+ i97=97;
+ j97=33;
+ iff=1;
+ }
+ for(ivec=1;ivec<=len;ivec++)
+ { uni=uu[i97]-uu[j97];
+ if(uni<0.) uni=uni+1.;
+ uu[i97]=uni;
+ i97=i97-1;
+ if(i97==0) i97=97;
+ j97=j97-1;
+ if(j97==0) j97=97;
+ c=c-cd;
+ if(c<0.) c=c+cm;
+ uni=uni-c;
+ if(uni<0.) uni=uni+1.;
+ if(uni==0.)
+ { uni=uu[j97]*0.59604644775391e-07;
+ if(uni==0.) uni=0.35527136788005e-14;
+ }
+ u[ivec]=uni;
+ }
+
+ // cout << endl<< "random: " << u[1] << endl << flush;
+
+ return;
+}
+
+double random_james()
+{
+// This function computes 1 random number in the (0,1) interval,
+// using the subrn subroutine.
+ double u[2];
+ subrn(u,1);
+ return u[1];
+}
+
+///////////////////////////////////////////////////////////
+
+
+
diff --git a/src/main/resources/csource/random.h b/src/main/resources/csource/random.h
new file mode 100644
index 0000000..3d1931e
--- /dev/null
+++ b/src/main/resources/csource/random.h
@@ -0,0 +1,33 @@
+/* $Id$ */
+/*
+ * random.h
+ *
+ * written by Sebastian Vц╤cking
+ *
+ * The random module is a frontend for different random generators.
+ */
+
+/*
+ * The different generators
+ */
+typedef enum
+{
+ RANDOM_STDLIB, /* standard c random generator (the default one)*/
+ RANDOM_CW,
+ RANDOM_JAMES
+} RandomMethode;
+
+/*
+ * Selects a generator
+ */
+void random_set_method(RandomMethode method);
+
+/*
+ * Returns a number between 0 and 1 generated with the selected generator.
+ */
+double random_get();
+
+/*
+ * Seeds the random generator if it supports it.
+ */
+void random_seed(int seed);
diff --git a/src/main/resources/csource/scatter.c b/src/main/resources/csource/scatter.c
new file mode 100644
index 0000000..edc65fd
--- /dev/null
+++ b/src/main/resources/csource/scatter.c
@@ -0,0 +1,1094 @@
+/* $Id$ */
+/*
+ * scatter.c
+ *
+ * written by Sebastian Voecking
+ *
+ * See scatter.h for details
+ *
+ * Included in this file are function from Ferenc Glueck for calculation of
+ * cross sections.
+ */
+
+//#include "scatter.h"
+
+#include
+#include
+#include
+#include "constants.h"
+#include "random.h"
+
+double sigmael(double E);
+double sigmaexc(double E);
+double sigmaion(double E);
+void randomel(double E, double *Eloss, double *theta);
+void randomexc(double E, double *Eloss, double *theta);
+void randomion(double E, double *Eloss, double *theta);
+
+static double sigmaBC(double E);
+static double sigmadiss10(double E);
+static double sigmadiss15(double E);
+
+static double Del(double E, double c);
+static double Dexc(double E, double c);
+static double sumexc(double K);
+static double Dinel(double E, double c);
+static void gensecelen(double E, double *W);
+
+static void subrn(double* u, int len);
+static double lagrange(int n, double *xn, double *fn, double x);
+
+static const double emass = 18780; // Electron mass in atomic units
+static const double a02 = 28e-22; // Bohr radius squared
+static const double clight = 137; // velocity of light in atomic units
+static const double R = 13.6; // Ryberg energy in eV
+
+static double density; /* Density of the residual gas */
+static double event_x; /* Distance to the interaction point */
+static double x; /* Already moved distance */
+static double totalx;
+static double energy_loss;
+static double scatter_angle;
+
+static FILE* output = NULL;
+
+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
+ static double e[14] = {0., 1.5, 5., 7., 10., 15., 20., 30., 60., 100., 150., 200., 300., 400.};
+ static double s[14] = {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 * M_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. * M_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.
+ //
+ static 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. * M_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. * M_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
+ static double aB[9] = {-4.2935194e2, 5.1122109e2, -2.8481279e2,
+ 8.8310338e1, -1.6659591e1, 1.9579609,
+ -1.4012824e-1, 5.5911348e-3, -9.5370103e-5};
+ static double aC[9] = {-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
+ static double a[9] = {-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
+ static double a[9] = {-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;
+}
+
+void randomel(double E, double *Eloss, double *theta) {
+ // 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
+ static double H2molmass = 69.e6;
+ double T, c, b, u[3], G, a, gam, K2, Gmax;
+ int i;
+ 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++) {
+ subrn(u, 2);
+ c = 1. + b - b * (2. + b) / (b + 2. * u[1]);
+ K2 = 2. * T * (1. + gam) * fabs(1. - c); // momentum transfer squared
+ a = (4. + K2)*(4. + K2) / (gam * gam);
+ G = a * Del(E, c);
+ if (G > Gmax * u[2]) break;
+ }
+ *theta = acos(c)*180. / M_PI;
+ *Eloss = 2. * emass / H2molmass * (1. - c) * E;
+ return;
+}
+
+void randomexc(double E, double *Eloss, double *theta) {
+ // 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
+ static int iff = 0;
+ static double sum[1001], fmax, Ecen = 12.6 / 27.21;
+ double T, c = 0., u[3], 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)
+ static double En[7] = {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 )
+ static double p[7] = {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)
+ static double EB[28] = {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)
+ static double EC[14] = {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 )
+ static double pB[28] = {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 )
+ static double pC[14] = {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};
+
+ 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 (iff == 0) {
+ fmax = 0;
+ 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;
+ iff = 1;
+ }
+ //
+ // Scattering angle *theta generation:
+ //
+ T = E / 27.2;
+ 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:
+ for (j = 1; j < 5000; j++) {
+ subrn(u, 2);
+ y = ymin + (ymax - ymin) * u[1];
+ K = exp(y / 2.);
+ fy = sumexc(K);
+ if (fmax * u[2] < fy) break;
+ }
+ // Calculation of c=cos(theta) and theta:
+ x = exp(y);
+ c = 1. - (x - xmin) / (4. * T);
+ *theta = acos(c)*180. / M_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++) {
+ subrn(u, 2);
+ c = -1. + 2. * u[1];
+ D = Dexc(E, c)*1.e22;
+ if (Dmax * u[2] < D) break;
+ }
+ *theta = acos(c)*180. / M_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++) {
+ subrn(u, 2);
+ n = (int) (N * u[1]);
+ if (u[2] * pmax < p[n]) break;
+ }
+ if (n < 0) n = 0;
+ if (n > 6) n = 6;
+ if (n > 1) // Bp, Bpp, D, Dp, EF states
+ {
+ *Eloss = En[n]*27.2;
+ return;
+ }
+ if (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
+ for (j = 1; j < 5000; j++) {
+ subrn(u, 2);
+ v = (int) (N * u[1]);
+ if (u[2] * pmax < pB[v]) break;
+ }
+ if (v < 0) v = 0;
+ if (v > 27) v = 27;
+ *Eloss = EB[v]*27.2;
+ }
+ if (n == 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++) {
+ subrn(u, 2);
+ v = (int) (N * u[1]);
+ if (u[2] * pmax < pC[v]) break;
+ }
+ if (v < 0) v = 0;
+ if (v > 13) v = 13;
+ *Eloss = EC[v]*27.2;
+ }
+ return;
+}
+
+void randomion(double E, double *Eloss, double *theta) {
+ // 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
+ //
+ static double Ei = 15.45 / 27.21;
+ double c, b, u[3], 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;
+ //
+ // 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:
+ for (j = 1; j < 5000; j++) {
+ subrn(u, 2);
+ y = ymin + (ymax - ymin) * u[1];
+ K = exp(y / 2.);
+ c = 1. + b - K * K / (4. * T);
+ G = K * K * (Dinel(E, c) - Dexc(E, c));
+ if (Gmax * u[2] < G) break;
+ }
+ // y --> x --> c --> theta
+ x = exp(y);
+ c = 1. - (x - xmin) / (4. * T);
+ *theta = acos(c)*180. / M_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.) {
+ gensecelen(E, &W);
+ *Eloss = 15.45 + W;
+ return;
+ }
+ // For theta>=20 the free electron model is used
+ // (with full correlation between theta and Eloss)
+ if (*theta >= 20.) {
+ *Eloss = E * (1. - c * c);
+ return;
+ }
+ // 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. / M_PI * (q / (1. + q * q) + atan(q));
+ q = qmin;
+ Fmin = 1. / 2. + 1. / M_PI * (q / (1. + q * q) + atan(q));
+ h = Fmax - Fmin;
+ // Generation of Eloss with the Neumann acceptance-rejection method:
+ 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) :
+ subrn(u, 2);
+ F = Fmin + h * u[1];
+ y = 0.;
+ for (i = 1; i <= 30; i++) {
+ G = 1. / 2. + (y + sin(2. * y) / 2.) / M_PI;
+ Gp = (1. + cos(2. * y)) / M_PI;
+ y = y - (G - F) / Gp;
+ if (fabs(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. * fabs(El - Ei) + 1.e-8);
+ st1 = K2 - 2. * El + 2.;
+ if (fabs(st1) < 1.e-9) st1 = 1.e-9;
+ arg = 2. * kej / st1;
+ if (arg >= 0.)
+ arctg = atan(arg);
+ else
+ arctg = atan(arg) + M_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. * M_PI / kej));
+ D2ion = 2. * kf / ki * Rex / (El * K2) * fE;
+ K = KK;
+ //
+ WcE = D2ion;
+ Jstarq = 16. / (3. * M_PI * (1. + q * q)*(1. + q * q));
+ WcstarE = 4. / (K * K * K * K * K) * Jstarq;
+ w = WcE / WcstarE;
+ if (wmax * u[2] < w) break;
+ }
+ //
+ *Eloss = El * 27.2;
+ //
+ return;
+}
+
+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
+ static double Cel[50] = {
+ -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
+ };
+ static double e[10] = {0., 3., 6., 12., 20., 32., 55., 85., 150., 250.};
+ static double t[10] = {0., 10., 20., 30., 40., 60., 80., 100., 140., 180.};
+ static double D[9][10] = {
+ {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. / M_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;
+ static double EE = 12.6 / 27.2;
+ static double e[5] = {0., 25., 35., 50., 100.};
+ static double t[9] = {0., 10., 20., 30., 40., 60., 80., 100., 180.};
+ static double D[4][9] = {
+ {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. / M_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 sumexc(double K) {
+ static double Kvec[15] = {0., 0.1, 0.2, 0.4, 0.6, 0.8, 1., 1.2, 1.5, 1.8, 2., 2.5, 3., 4., 5.};
+ static double fvec[7][15] ={
+ {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}};
+ static double EeV[7] = {12.73, 13.20, 14.77, 15.3, 14.93, 15.4, 13.06};
+ int n, j, jmin = 0, nmax;
+ double En, f[7], x4[4], f4[4], sum;
+ //
+ 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 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
+ static double Cinel[50] = {
+ -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
+ };
+ static 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;
+}
+
+
+////////////////////////////////////////////////////////////////////
+
+void gensecelen(double E, double *W) {
+ // 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.
+ static double Ei = 15.45, eps2 = 14.3, b = 6.25;
+ static double B;
+ double D, A, eps, a, u, epsmax;
+ static int iff = 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 = random_get();
+ a = b / D * (u + D / b * B);
+ eps = eps2 + b * tan(a);
+ *W = eps - Ei;
+ return;
+}
+
+void subrn(double *u, int len) {
+ // This subroutine computes random numbers u[1],...,u[len]
+ // in the (0,1) interval. It uses the 0= 900000000) ijkl = 1;
+ ij = ijkl / 30082;
+ kl = ijkl - 30082 * ij;
+ i = ((ij / 177) % 177) + 2;
+ j = (ij % 177) + 2;
+ k = ((kl / 169) % 178) + 1;
+ l = kl % 169;
+ for (ii = 1; ii <= 97; ii++) {
+ s = 0;
+ t = 0.5;
+ for (jj = 1; jj <= 24; jj++) {
+ m = (((i * j) % 179) * k) % 179;
+ i = j;
+ j = k;
+ k = m;
+ l = (53 * l + 1) % 169;
+ if ((l * m) % 64 >= 32) s = s + t;
+ t = 0.5 * t;
+ }
+ uu[ii] = s;
+ }
+ c = 362436. / 16777216.;
+ cd = 7654321. / 16777216.;
+ cm = 16777213. / 16777216.;
+ i97 = 97;
+ j97 = 33;
+ iff = 1;
+ }
+ for (ivec = 1; ivec <= len; ivec++) {
+ uni = uu[i97] - uu[j97];
+ if (uni < 0.) uni = uni + 1.;
+ uu[i97] = uni;
+ i97 = i97 - 1;
+ if (i97 == 0) i97 = 97;
+ j97 = j97 - 1;
+ if (j97 == 0) j97 = 97;
+ c = c - cd;
+ if (c < 0.) c = c + cm;
+ uni = uni - c;
+ if (uni < 0.) uni = uni + 1.;
+ if (uni == 0.) {
+ uni = uu[j97]*0.59604644775391e-07;
+ if (uni == 0.) uni = 0.35527136788005e-14;
+ }
+ u[ivec] = uni;
+ }
+
+ // cout << endl<< "random: " << u[1] << endl << flush;
+
+ return;
+}
+
+double lagrange(int n, double *xn, double *fn, double x) {
+ int i, j;
+ double f, a[100], b[100], aa, bb;
+ 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;
+}
+
+/* The following functions are used to calculate the new impulse after an
+ * interaction.
+ * In these functions a vector is a double[3] array and a matrix a double[9]
+ * array.*/
+
+/* Sets the matrix to a rotation matrix which performs a rotation around the
+ * y axis of the angle theta and then around the z axis of phi. theta and
+ * phi are given in radians */
+static void matrix_rot(double theta, double phi, double* matrix) {
+ double cos_theta = cos(theta);
+ double sin_theta = sin(theta);
+ double cos_phi = cos(phi);
+ double sin_phi = sin(phi);
+
+ matrix[0] = cos_theta*cos_phi;
+ matrix[1] = -sin_phi;
+ matrix[2] = sin_theta*cos_phi;
+ matrix[3] = cos_theta*sin_phi;
+ matrix[4] = cos_phi;
+ matrix[5] = sin_theta*sin_phi;
+ matrix[6] = -sin_theta;
+ matrix[7] = 0;
+ matrix[8] = cos_theta;
+}
+
+/* Multiplicates the matrix a with the vector x and returns the result as
+ * vector y*/
+static void matrix_vec(double* a, double* x, double* y) {
+ int i, j;
+
+ for (i = 0; i < 3; i++) {
+ y[i] = 0;
+ for (j = 0; j < 3; j++) {
+ y[i] += a[3 * i + j] * x[j];
+ }
+ }
+}
+
+/* Returns the length of the vector x */
+static double vector_abs(double* x) {
+ double x2 = 0;
+ int i;
+
+ for (i = 0; i < 3; i++)
+ x2 += x[i] * x[i];
+
+ return sqrt(x2);
+}
+
+/* Returns the polar angle phi of the vector x in radians */
+static double vector_phi(double* x) {
+ return atan2(x[1], x[0]);
+}
+
+/* Returns the azimutal angle theta of the vector x in radians */
+static double vector_theta(double* x) {
+ return acos(x[2] / vector_abs(x));
+}
+
+/* Creates a vector from spheric coordinates and returns it as x */
+static void vector_spheric(double r, double theta, double phi, double* x) {
+ double sin_theta = sin(theta);
+ double cos_theta = cos(theta);
+ double sin_phi = sin(phi);
+ double cos_phi = cos(phi);
+
+ x[0] = r * sin_theta * cos_phi;
+ x[1] = r * sin_theta * sin_phi;
+ x[2] = r * cos_theta;
+}
+
+static double sigma_total(double E) {
+ return 1e4 * (sigmaion(E) // sigma for ionisation of H2
+ + sigmaexc(E) // sigma for excitation of H2
+ + sigmael(E)); // sigma for elastic scattering
+}
+
+/*
+void scatter_init(double pressure, double temperature)
+{
+ random_set_method(RANDOM_CW);
+ //random_seed(time(0));
+ //output = fopen("scatter.dat", "a");
+
+ density = pressure / BOLTZMANN / temperature;
+ energy_loss = 0;
+ scatter_angle = 0;
+
+ scatter_set_output_file(NULL);
+ scatter_reset();
+}
+ */
+
+void scatter_reset() {
+ x = 0;
+ totalx = 0;
+ event_x = -log(random_get());
+}
+
+void scatter_move(double dx, double E) {
+ double sigma = sigma_total(E);
+ x += dx * density * sigma;
+ totalx += dx;
+ //fprintf(stderr, "%e\n", x);
+}
+
+int scatter_has_event() {
+ return (x >= event_x);
+}
+
+/*
+Interaction scatter_event(double *P, double E)
+{
+ double el = sigmael(E);
+ double exc = sigmaexc(E);
+ double ion = sigmaion(E);
+ double random = (el + exc + ion) * random_get();
+
+ double phi_old = vector_phi(P);
+ double theta_old = vector_theta(P);
+
+ double eloss, angle, theta, phi;
+ double matrix[9];
+ double y[3];
+ double p2;
+ double Enew;
+
+ Interaction interaction;
+
+ if(random < el) {
+ randomel(E, &eloss, &angle);
+ interaction = ELASTIC;
+ }
+ else if(random < exc) {
+ randomexc(E, &eloss, &angle);
+ interaction = EXCITATION;
+ }
+ else {
+ randomion(E, &eloss, &angle);
+ interaction = IONIZATION;
+ }
+
+ Enew = (E - eloss) * CHARGEUNIT;
+ p2 = Enew*Enew/C/C + 2*Enew*ELECTRON_MASS;
+ theta = M_PI / 180 * angle;
+ phi = 2 * M_PI * random_get();
+ vector_spheric(sqrt(p2), theta, phi, y);
+
+ matrix_rot(theta_old, phi_old, matrix);
+ matrix_vec(matrix, y, P);
+
+ energy_loss = eloss;
+ scatter_angle = angle;
+ scatter_reset(E);
+
+ if(output)
+ fprintf(output, "%f\t%d\t%f\t%f\n", E, interaction, eloss, angle);
+
+ return interaction;
+}
+ */
+
+double scatter_get_energy_loss() {
+ return energy_loss;
+}
+
+double scatter_get_scatter_angle() {
+ return scatter_angle;
+}
+
+double scatter_free_path(double E) {
+ return 1 / density / sigma_total(E);
+}
+
+void scatter_set_output_file(const char* filename) {
+ if (output) {
+ fclose(output);
+ output = NULL;
+ }
+
+ if (filename) {
+ output = fopen(filename, "w");
+ if (output)
+ fprintf(output, "#Energy\tEvent\tEloss\tAngle\n");
+ else
+ fprintf(stderr,
+ "Warning: Cannot write to file '%s'. Debugging output will be"
+ "omitted.\n", filename);
+ }
+}