diff --git a/build.gradle b/build.gradle index 1233762..d464d5c 100644 --- a/build.gradle +++ b/build.gradle @@ -1,5 +1,5 @@ plugins { - id "org.jetbrains.kotlin.jvm" version "1.2.30" + id "org.jetbrains.kotlin.jvm" version "1.3.10" id "java" id "application" } @@ -18,6 +18,10 @@ repositories { dependencies { compile group: 'org.apache.commons', name: 'commons-math3', version: '3.6.1' compile "org.jetbrains.kotlin:kotlin-stdlib-jdk8" + + // https://mvnrepository.com/artifact/org.apache.commons/commons-rng-simple + compile group: 'org.apache.commons', name: 'commons-rng-simple', version: '1.1' + testCompile group: 'junit', name: 'junit', version: '4.12' } @@ -25,12 +29,12 @@ compileKotlin { kotlinOptions.jvmTarget = "1.8" } -task runTrap(type: JavaExec) { - classpath = sourceSets.main.runtimeClasspath - - main = mainClassName - - // arguments to pass to the application - args 'D:\\Work\\Numass\\trapping\\trap.out' -} +//task runTrap(type: JavaExec) { +// classpath = sourceSets.main.runtimeClasspath +// +// main = mainClassName +// +// // arguments to pass to the application +// args 'D:\\Work\\Numass\\trapping\\trap.out' +//} diff --git a/notebooks/cross-sections.ipynb b/notebooks/cross-sections.ipynb index 4bf4d34..38b4963 100644 --- a/notebooks/cross-sections.ipynb +++ b/notebooks/cross-sections.ipynb @@ -11,33 +11,21 @@ }, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Added jar: [commons-math3-3.6.1.jar]\n", - "Added jar: [trapping-dev.jar]\n" - ] - } - ], - "source": [ - "%classpath add jar ../build/install/trapping/lib/commons-math3-3.6.1.jar\n", - "%classpath add jar ../build/install/trapping/lib/trapping-dev.jar" - ] - }, - { - "cell_type": "code", - "execution_count": 17, - "metadata": { - "ExecuteTime": { - "end_time": "2018-04-08T18:58+0000", - "start_time": "2018-04-08T18:58+0000" - } - }, - "outputs": [ + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "3161978b-3c71-4d18-a97e-f972beece3ff", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + }, { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "30d035de-65b2-4c03-96c6-25b897f10de4", + "model_id": "6884d267-e5ee-4841-8c60-9c66eea437c5", "version_major": 2, "version_minor": 0 }, @@ -48,11 +36,39 @@ } ], "source": [ - "def scatter = inr.numass.trapping.Scatter.INSTANCE\n", + "%classpath add jar ../build/install/trapping/lib/commons-math3-3.6.1.jar\n", + "%classpath add jar ../build/install/trapping/lib/trapping-dev.jar" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "ExecuteTime": { + "end_time": "2018-04-08T18:58+0000", + "start_time": "2018-04-08T18:58+0000" + } + }, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "582da53b-a518-4a8a-b7ce-84d38c580cf3", + "version_major": 2, + "version_minor": 0 + }, + "method": "display_data" + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "scatter = inr.numass.trapping.Scatter.INSTANCE\n", "\n", "def e = (1..20) // energy in keV\n", "\n", - "def factor = 1e-22\n", + "factor = 1e-22\n", "def sigmaEl = e.collect{scatter.sigmael(it*1000)/factor}\n", "def sigmaIon = e.collect{scatter.sigmaion(it*1000)/factor}\n", "def sigmaExc = e.collect{scatter.sigmaexc(it*1000)/factor}\n", @@ -63,6 +79,112 @@ "plot << new Line(x: e, y: sigmaExc, displayName: \"Excitation\")\n" ] }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "E(keV)\telastic\tion\texc\n", + "12.0\t0.43948540252301355\t2.7362696996880764\t2.2754271467658227\n", + "12.1\t0.4360219815420227\t2.716347129996631\t2.259247553470829\n", + "12.2\t0.4326153394877464\t2.6967291910800784\t2.2433117665045716\n", + "12.3\t0.42926409223681455\t2.6774088086825643\t2.227614186457129\n", + "12.4\t0.42596690029168804\t2.6583791281775273\t2.212149386181278\n", + "12.5\t0.4227224669965498\t2.639633506055748\t2.1969121041685504\n", + "12.6\t0.41952953683810995\t2.621165501807934\t2.1818972382301722\n", + "12.7\t0.4163868938266485\t2.602968870180603\t2.1670998394665926\n", + "12.8\t0.4132933599529065\t2.585037553785318\t2.1525151065102723\n", + "12.9\t0.4102477937167139\t2.567365676042546\t2.138138380027344\n", + "13.0\t0.4072490887234942\t2.5499475344425586\t2.123965137464613\n", + "13.1\t0.4042961723450181\t2.532777594106803\t2.1099909880291814\n", + "13.2\t0.40138800444100226\t2.515850481634202\t2.096211667888705\n", + "13.3\t0.3985235761383543\t2.4991609792177183\t2.082623035581028\n", + "13.4\t0.3957019086650494\t2.4827040190174023\t2.069221067622563\n", + "13.5\t0.39292205223581\t2.466474677776929\t2.0560018543054204\n", + "13.6\t0.39018308498692494\t2.4504681716713828\t2.042961595673853\n", + "13.7\t0.3874841119576961\t2.434679851374739\t2.0300965976711156\n", + "13.8\t0.3848242641161482\t2.419105197336177\t2.0174032684483527\n", + "13.9\t0.3822026974267757\t2.403739815254919\t2.00487811482759\n", + "14.0\t0.3796185919582249\t2.3885794317439313\t1.9925177389113418\n", + "14.1\t0.3770711510289289\t2.373619890173301\t1.9803188348317768\n", + "14.2\t0.3745596003888268\t2.358857146684646\t1.9682781856327505\n", + "14.3\t0.3720831874354002\t2.3442872663683785\t1.9563926602783916\n", + "14.4\t0.36964118046236066\t2.32990641959609\t1.9446592107822769\n", + "14.5\t0.36723286793941123\t2.315710878500726\t1.9330748694515214\n", + "14.6\t0.3648575578215951\t2.3016970135976464\t1.9216367462404509\n", + "14.7\t0.36251457688682004\t2.287861290540006\t1.910342026208781\n", + "14.8\t0.3602032701002296\t2.274200267002247\t1.8991879670794969\n", + "14.9\t0.357923000004158\t2.26071058968582\t1.8881718968918926\n", + "15.0\t0.3556731461324769\t2.2473889914415754\t1.8772912117454528\n", + "15.1\t0.353453104448203\t2.234232288503513\t1.8665433736304802\n", + "15.2\t0.35126228680329497\t2.221237377828907\t1.8559259083415982\n", + "15.3\t0.3491001204196282\t2.2084012345400286\t1.845436403470438\n", + "15.4\t0.3469660473901797\t2.1957209094629615\t1.835072506474014\n", + "15.5\t0.3448595241995173\t2.1831935267592297\t1.824831922815469\n", + "15.6\t0.3427800212627204\t2.170816281646157\t1.8147124141740314\n", + "15.7\t0.34072702248191666\t2.158586438202101\t1.8047117967211865\n", + "15.8\t0.33870002481965195\t2.1465013272528792\t1.7948279394602165\n", + "15.9\t0.3366985378883526\t2.134558344335901\t1.7850587626263867\n", + "16.0\t0.33472208355517563\t2.122754947738674\t1.7754022361452126\n", + "16.1\t0.33277019556158083\t2.111088656608532\t1.7658563781463463\n", + "16.2\t0.33084241915698676\t2.0995570491305764\t1.7564192535307461\n", + "16.3\t0.32893831074590485\t2.0881577607709536\t1.7470889725889123\n", + "16.4\t0.3270574375479805\t2.0768884825827687\t1.7378636896680617\n", + "16.5\t0.3251993772703896\t2.0657469595720115\t1.7287416018862318\n", + "16.6\t0.32336371779207285\t2.054730989121034\t1.7197209478913829\n", + "16.7\t0.3215500568593091\t2.043838419467228\t1.7108000066636726\n", + "16.8\t0.31975800179215863\t2.0330671482346325\t1.7019770963591483\n", + "16.9\t0.3179871692013229\t2.022415121016355\t1.693250573193194\n", + "17.0\t0.3162371847149949\t2.011880330005736\t1.684618830362136\n", + "17.1\t0.3145076827152896\t2.0014608126743285\t1.6760802970014965\n", + "17.2\t0.31279830608386505\t1.991154650494809\t1.6676334371794341\n", + "17.3\t0.31110870595636253\t1.9809599677070577\t1.6592767489239986\n", + "17.4\t0.3094385414853119\t1.9708749301257074\t1.6510087632828718\n", + "17.5\t0.3077874796111622\t1.9608977439875288\t1.6428280434143254\n", + "17.6\t0.30615519484111486\t1.9510266548371218\t1.6347331837082066\n", + "17.7\t0.3045413690354524\t1.9412599464494087\t1.6267228089357717\n", + "17.8\t0.3029456912010658\t1.9315959397875393\t1.6187955734272899\n", + "17.9\t0.3013678572919011\t1.9220329919948316\t1.6109501602763423\n", + "18.0\t0.29980757001605496\t1.9125694954194674\t1.603185280569817\n", + "18.1\t0.29826453864926233\t1.9032038766707016\t1.5954996726426307\n", + "18.2\t0.29673847885453164\t1.8939345957053948\t1.5878921013562524\n", + "18.3\t0.29522911250768985\t1.8847601449437406\t1.5803613574001418\n", + "18.4\t0.2937361675286164\t1.8756790484131034\t1.5729062566152603\n", + "18.5\t0.2922593777179467\t1.866689860918917\t1.565525639338836\n" + ] + }, + { + "data": { + "text/plain": [ + "null" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "println \"E(keV)\\telastic\\tion\\texc\"\n", + "for(d = 12.0; d < 18.6; d += 0.1){\n", + " def sigmaEl = scatter.sigmael(d*1000)/factor\n", + " def sigmaIon = scatter.sigmaion(d*1000)/factor\n", + " def sigmaExc = scatter.sigmaexc(d*1000)/factor\n", + " println \"$d\\t$sigmaEl\\t$sigmaIon\\t$sigmaExc\"\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": null, diff --git a/src/main/kotlin/inr/numass/trapping/FunctionCache.kt b/src/main/kotlin/inr/numass/trapping/FunctionCache.kt new file mode 100644 index 0000000..88a235e --- /dev/null +++ b/src/main/kotlin/inr/numass/trapping/FunctionCache.kt @@ -0,0 +1,23 @@ +package inr.numass.trapping + +import java.util.* + +class FunctionCache(val xPrecision: Double, val function: (Double) -> Double) { + private val values = TreeMap() + + operator fun invoke(x: Double): Double { + val floor: MutableMap.MutableEntry? = values.floorEntry(x) + val ceil: MutableMap.MutableEntry? = values.ceilingEntry(x) + if (floor == null || ceil == null || ceil.key - floor.key <= xPrecision) { + val newValue = function(x) + values[x] = newValue + return newValue + } else { + val x0 = floor.key + val x1 = ceil.key + val y0 = floor.value + val y1 = ceil.value + return y0 + (x - x0) * (y1 - y0) / (x1 - x0) + } + } +} \ No newline at end of file diff --git a/src/main/kotlin/inr/numass/trapping/RandomGeneratorBridge.kt b/src/main/kotlin/inr/numass/trapping/RandomGeneratorBridge.kt new file mode 100644 index 0000000..56bf65e --- /dev/null +++ b/src/main/kotlin/inr/numass/trapping/RandomGeneratorBridge.kt @@ -0,0 +1,28 @@ +package inr.numass.trapping + +import org.apache.commons.math3.random.RandomGenerator +import org.apache.commons.rng.UniformRandomProvider + +class RandomGeneratorBridge(val provider: UniformRandomProvider) : RandomGenerator { + override fun nextBoolean(): Boolean = provider.nextBoolean() + + override fun nextFloat(): Float = provider.nextFloat() + + override fun setSeed(seed: Int) = error("Not supported") + + override fun setSeed(seed: IntArray?) = error("Not supported") + + override fun setSeed(seed: Long) = error("Not supported") + + override fun nextBytes(bytes: ByteArray) = provider.nextBytes(bytes) + + override fun nextInt(): Int = provider.nextInt() + + override fun nextInt(n: Int): Int = provider.nextInt(n) + + override fun nextGaussian(): Double = error("Not supported") + + override fun nextDouble(): Double = provider.nextDouble() + + override fun nextLong(): Long = provider.nextLong() +} \ No newline at end of file diff --git a/src/main/kotlin/inr/numass/trapping/Scatter.kt b/src/main/kotlin/inr/numass/trapping/Scatter.kt index 78657fa..c49a803 100644 --- a/src/main/kotlin/inr/numass/trapping/Scatter.kt +++ b/src/main/kotlin/inr/numass/trapping/Scatter.kt @@ -407,7 +407,7 @@ object Scatter { 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 + y -= (G - F) / Gp if (abs(G - F) < 1e-8) { break } @@ -666,16 +666,23 @@ object Scatter { return Dinelreturn } - private fun sumexc(K: Double): Double { + 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 fvec = arrayOf( + doubleArrayOf(2.907e-1, 2.845e-1, 2.665e-1, 2.072e-1, 1.389e-1, + 8.238e-2, 4.454e-2, 2.269e-2, 7.789e-3, 2.619e-3, 1.273e-3, 2.218e-4, 4.372e-5, 2.889e-6, 4.247e-7), // B + doubleArrayOf(3.492e-1, 3.367e-1, 3.124e-1, 2.351e-1, 1.507e-1, + 8.406e-2, 4.214e-2, 1.966e-2, 5.799e-3, 1.632e-3, 6.929e-4, 8.082e-5, 9.574e-6, 1.526e-7, 7.058e-9), // C + doubleArrayOf(6.112e-2, 5.945e-2, 5.830e-2, 5.072e-2, 3.821e-2, + 2.579e-2, 1.567e-2, 8.737e-3, 3.305e-3, 1.191e-3, 6.011e-4, 1.132e-4, 2.362e-5, 1.603e-6, 2.215e-7), // Bp + doubleArrayOf(2.066e-2, 2.127e-2, 2.137e-2, 1.928e-2, 1.552e-2, + 1.108e-2, 7.058e-3, 4.069e-3, 1.590e-3, 5.900e-4, 3.046e-4, 6.142e-5, 1.369e-5, 9.650e-7, 1.244e-7), // Bpp + doubleArrayOf(9.405e-2, 9.049e-2, 8.613e-2, 7.301e-2, 5.144e-2, + 3.201e-2, 1.775e-2, 8.952e-3, 2.855e-3, 8.429e-4, 3.655e-4, 4.389e-5, 5.252e-6, 9.010e-8, 7.130e-9), // D + doubleArrayOf(4.273e-2, 3.862e-2, 3.985e-2, 3.362e-2, 2.486e-2, + 1.612e-2, 9.309e-3, 4.856e-3, 1.602e-3, 4.811e-4, 2.096e-4, 2.498e-5, 2.905e-6, 5.077e-8, 6.583e-9), // Dp + doubleArrayOf(0.000e-3, 2.042e-3, 7.439e-3, 2.200e-2, 3.164e-2, + 3.161e-2, 2.486e-2, 1.664e-2, 7.562e-3, 3.044e-3, 1.608e-3, 3.225e-4, 7.120e-5, 6.290e-6, 1.066e-6))// EF val EeV = doubleArrayOf(12.73, 13.20, 14.77, 15.3, 14.93, 15.4, 13.06) var n: Int = 0 var j: Int @@ -690,33 +697,32 @@ object Scatter { // 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 + when { + K >= 5.0 -> f[n] = 0.0 + K in 3.0..4.0 -> f[n] = fvec[n][12] + (K - 3.0) * (fvec[n][13] - fvec[n][12]) + K in 4.0..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++ } - 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) } - 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++ @@ -724,15 +730,18 @@ object Scatter { return sum } + //Caching function for precision + private val sumexcCache = FunctionCache(0.01, this::_sumexc) + private fun sumexc(K: Double): Double = sumexcCache(K) + fun lagrange(n: Int, xn: DoubleArray, fn: DoubleArray, x: Double): Double { var i: Int - var j: Int + var j: Int = 0 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) { diff --git a/src/main/kotlin/inr/numass/trapping/SimulationManager.kt b/src/main/kotlin/inr/numass/trapping/SimulationManager.kt index c4e7276..6b937aa 100644 --- a/src/main/kotlin/inr/numass/trapping/SimulationManager.kt +++ b/src/main/kotlin/inr/numass/trapping/SimulationManager.kt @@ -4,6 +4,8 @@ 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 org.apache.commons.rng.UniformRandomProvider +import org.apache.commons.rng.simple.RandomSource import java.io.* import java.util.stream.Stream @@ -13,7 +15,7 @@ import java.util.stream.Stream */ class SimulationManager { - var generator: RandomGenerator = JDKRandomGenerator() + var generator: UniformRandomProvider = RandomSource.create(RandomSource.SPLIT_MIX_64) /** * output for accepted events @@ -88,11 +90,11 @@ class SimulationManager { * @return */ @Synchronized - fun simulateAll(num: Int): Counter { + fun simulateAll(num: Number): Counter { val counter = Counter() - val simulator = Simulator(eLow, thetaTransport, thetaPinch, gasDensity, bSource, magneticField, generator) + val simulator = Simulator(eLow, thetaTransport, thetaPinch, gasDensity, bSource, magneticField, RandomGeneratorBridge(generator)) - System.out.printf("%nStarting sumulation with initial energy %g and %d electrons.%n%n", initialE, num) + System.out.printf("%nStarting sumulation with initial energy %g and %d electrons.%n%n", initialE, num.toLong()) output.printf("%s\t%s\t%s\t%s\t%s\t%s%n", "E", "theta", "theta_start", "colNum", "L", "state") Stream.generate { getRandomTheta() }.limit(num.toLong()).parallel() .forEach { theta -> diff --git a/src/main/kotlin/inr/numass/trapping/Trapping.kt b/src/main/kotlin/inr/numass/trapping/Trapping.kt index 9191061..b718209 100644 --- a/src/main/kotlin/inr/numass/trapping/Trapping.kt +++ b/src/main/kotlin/inr/numass/trapping/Trapping.kt @@ -1,5 +1,6 @@ package inr.numass.trapping +import org.apache.commons.rng.simple.RandomSource import java.time.Duration import java.time.Instant @@ -13,11 +14,13 @@ fun main(args: Array) { System.out.printf("Starting at %s%n%n", startTime.toString()) SimulationManager().apply { + generator = RandomSource.create(RandomSource.SPLIT_MIX_64) setOutputFile("D:\\Work\\Numass\\trapping\\test1.out") setFields(0.6, 3.7, 7.2) gasDensity = 1e19 - reportFilter = {true} - }.simulateAll(1e5.toInt()) + initialE = 14000.0 + range = 4000.0 + }.simulateAll(1e6) val finishTime = Instant.now() System.out.printf("%nFinished at %s%n", finishTime.toString())