From 84ef96cf8a40f6c14c36df6d1e909ffeacc5f381 Mon Sep 17 00:00:00 2001 From: Isaac Meyer Date: Fri, 6 Jul 2018 09:57:17 -0500 Subject: [PATCH] added reconstruction ability for resonance samples --- .../nuclear-data-resonance-covariance.ipynb | 106 +++++++++++++++--- openmc/data/resonance.py | 20 ++-- openmc/data/resonance_covariance.py | 66 ++++++++--- 3 files changed, 154 insertions(+), 38 deletions(-) diff --git a/examples/jupyter/nuclear-data-resonance-covariance.ipynb b/examples/jupyter/nuclear-data-resonance-covariance.ipynb index d1def693c6..88c421ce5f 100644 --- a/examples/jupyter/nuclear-data-resonance-covariance.ipynb +++ b/examples/jupyter/nuclear-data-resonance-covariance.ipynb @@ -118,7 +118,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 5, @@ -129,7 +129,7 @@ "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUkAAAD8CAYAAAD6+lbaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3X2wXHWd5/H3JzeEJ+UxghCySxyi\nu4GaVcxEVmstxwgEnTVulc6EXceU4mTXAccZZ0th3RUXlyqYcWW1RLYiZATHIaQYXVJb0RhAtKZK\nIJFhxIAZIihciUAMMC5IHu797h/n18m5ndN9Tz/d233O55U6le7T53Sfcx++9/f4/SkiMDOzYnNm\n+wLMzIaZg6SZWRsOkmZmbThImpm14SBpZtaGg6SZWRsDC5KSVkjaIWmnpMsH9TlmZoOkQYyTlDQG\n/CNwPjAObAUujoiH+/5hZmYDNKiS5DJgZ0Q8FhH7gPXAygF9lpnZwMwd0PsuAJ7MPR8H3tTq4LFj\njo0jTjhpQJdiZgB7d43vjohXdXv+hb97bPxqz0SpY3/4o72bI2JFt581TAYVJFWwb0q9XtIaYA3A\n3ONP5Mw/+viALsXMAHZc9fGf93L+r/ZMcP/mf1bq2LHTHp3fy2cNk0FVt8eBhbnnZwBP5Q+IiLUR\nsTQilo4dc2y2T9k2k/r5me3ep/m1ft9n0fu1+ozGPQ/i693J+83099p6E8BkyX9VMqgguRVYLGmR\npHnAKmDjgD7LzGZAEOyPiVJbGdONgJF0pKTb0uv3SToz99oVaf8OSRfm9v+ZpO2SfizpVklH9Xrf\nAwmSEXEAuAzYDDwCbIiI7dOdp8i2mSxhND6zX+9V9rV+fWa792v1GY177ue9T/eZvR5rw6FfJck0\nAuZ64CJgCXCxpCVNh10CPBcRZwHXAdemc5eQFbzOBlYAX5Y0JmkB8CfA0og4BxhLx/VkUG2SRMQm\nYFM35/qXx2z4BMFE/4YMHhwBAyCpMQImP0xwJfCZ9Ph24EuSlPavj4i9wOOSdqb3e4Isph0taT9w\nDE3NfN0Y6hk3brMyGy6TRKmthKIRMAtaHZNqpy8AJ7c6NyJ+AXyOLFjuAl6IiO90cZtTDHWQVMDE\nkXGwg2HiqGjbGWFmgxPABFFqA+ZL2pbb1jS93bQjYNocU7hf0olkpcxFwOnAsZLe39FNFhhYdbtf\nxvYe+nrM2au27WxmNlglS4kAuyNiaZvXpx0BkztmXNJc4HhgT5tz3wE8HhHPAkj6BvBm4K/LXnSR\noS5JmtnwCGB/RKmthDIjYDYCq9Pj9wJ3RzaPeiOwKvV+LwIWA/eTVbPPk3RMartcTtZx3JOhL0nm\nvbxwP0c9eQQAanSgpYJmo7pdVKIMTd3f/Nz8NbHpxaGqdO/vFXFAUmMEzBiwLiK2S7oK2BYRG4Gb\ngK+ljpk9pJ7qdNwGsk6eA8ClETEB3CfpduCBtP/vgbW9XutIBcmjnzji0C9zU6vEbA6/qQJ/TWxa\nARN9/DkpGgETEZ/OPX4ZeF+Lc68Gri7YfyVwZf+ucsSCJPiX2Wy2ZDNu6mfkgmRDu+p14/XGa+2q\nktO9T1W5em2dExOFHcvVNrJB0sxmVtZx4yA5MppLic0lwnwpqZP2yrqo631b97Jxkg6SI6coKJrZ\nYEy6JDm66tq22C23SVqnXJIccQqIORw+sckKOUBapwIxUcP5J5UJkpAbYG5mA+HqtplZC4HYF2Oz\nfRkzruuys6SFkr4r6ZGUCfhjaf9JkrZIejT9f2L/Lre8dhmDzFmTrHPZYPI5pbYq6eVuDgB/HhH/\nEjgPuDRlDL4cuCsiFgN3peczbuxlR4F23CZp3ZhIA8qn26qk6+p2ROwiS2xJRPxa0iNkyTBXAm9L\nh90M3AN8sqer7JIDgVn/RIiJqFYpsYy+tEmmBXreANwHnJoCKBGxS9Ip/fiMXsQYqNzaRJUy3XRM\n/xGxTk1WrJRYRs9BUtIrgL8F/jQi/ilL41bqvCnrbg+SJuo5jtIzjayfso6b+vX19lR2lnQEWYD8\nekR8I+1+WtJp6fXTgGeKzi1ad9vMhpc7bjqUMv/eBDwSEZ/PvZTPJrwauKP7y+ufQSyfalY3E6FS\nW5X0UnZ+C/CHwEOSHkz7/gtwDbBB0iVk6dQLk2bOJrfH+WtgnfOMmw5FxN9RvGoZZGtLDK1G1qCi\n7EF1Ubf7tf6YdO92fTh7kFlnsgQXDpK10/jDOOrzvouqz2UyskN//1B0Uo13lX+0BGJ/Dacl1j5I\nmlk5EdRyMHn97riJJrOtUbIa1Y65ohLZdOMkB9Hj38n7uRQ5asRkya3Uu0krJO2QtFPSYdOX07ra\nt6XX70uTVhqvXZH275B0YW7/CZJul/STlFfiX/d61y5JJm6jNGsv6F9JUtIYcD1wPjAObJW0MSIe\nzh12CfBcRJwlaRVwLfAHKUfEKuBs4HTgTkmvTWtvfwH4dkS8V9I84Jher7X2Jclmo1qSLNLuXqp0\nnzZzJphTaithGbAzIh6LiH3AerK8D3kryfI/ANwOLE/js1cC6yNib0Q8DuwElkk6Dngr2fhtImJf\nRDzf6z07SDZpDAsa9eo3eFqi9VcgJqPcBsyXtC23rWl6uwXAk7nn42lf4TERcQB4ATi5zbmvAZ4F\n/krS30u6UVLP0/lc3TazUrIlZUuHjN0RsbTN60XFj+Y/3a2OabV/LnAu8NGIuE/SF8hSNf63Etfb\nkkuSBRodGh6iYpZXLpdkyXyS48DC3PMzgKdaHSNpLnA8sKfNuePAeETcl/bfThY0e+Ig2UZ+Rk7V\nDOq+qvr1spTgIuaU2krYCiyWtCh1sKwiy/uQl88D8V7g7oiItH9V6v1eBCwG7o+IXwJPSnpdOmc5\n8DA9cnV7GlWdwjioe6jC18Za61fW8Yg4IOkyYDMwBqyLiO2SrgK2RcRGsg6Yr0naSVaCXJXO3S5p\nA1kAPABcmnq2AT4KfD0F3seAD/Z6rQ6SJbjqbZZlJu/n3O2I2ARsatr36dzjl2mRICcirgauLtj/\nINCuLbRjDpIlVa0kOaig72mJ1ZV13HhaoplZC17jxqbRKPXk2yhHtSQ0DG2So/q1q6us46Z+PXMO\nkl3wFEarqzqmSuv5jiWNpdHt/zc9X5Qmoz+aJqfP6/0yh9Mo/1Ed5Wu32dHhjJvK6MefhY8Bj+Se\nXwtcFxGLgefIJqlXkiLLR9mYxpifzjjsXAq2bnghsA5JOgN4F3Bjei7g7WQj3SGbnP6eXj5j2OWT\n9TrwWJVFwP7JOaW2Kum1TfJ/AZ8AXpmenww8nyajQ/GkdWBm1902s95l1e1qBcAyellS9veAZyLi\nh/ndBYcWlq+qtO72wbnes/jzU1TNn41UaaPS3GDd6ePc7ZHR65Ky75b0TuAo4DiykuUJkuam0mTR\npPXK0mQKlMFhuUoGPVyom8zkM3UdVg11HQLUddknIq6IiDMi4kyyOZV3R8R/AL5LNhkdssnpd/R8\nlSPkYBtl08+Sg4eNPvUzwcXIGMTdfBL4eJqUfjIpS3DdTM7wCFRXt20m9HONm1HRl1/liLgHuCc9\nfowsNXutzTkwda63q9s26rLebc/dNjMr1BhMXjcOkgOUL1W5hGVVULWqdBkOkjMk5kwdeN739y+o\nzrer4g9DqjQbLXXt3XaQnCGanBpA+p1JyG2SNhOq1nNdhoPkDCpK3OugYqMiQhxwkDQza83VbRu4\ng9XtxswcypcmO21jnI32QbdJVldd2yTrV3YeEo1OnOZla9v9DHbaxjgbwcoBstr6mU9S0gpJOyTt\nlHR5wetHppy0O1OO2jNzr12R9u+QdGHTeVNy3PbKQXIWFXXcOMjYsOpn0l1JY8D1wEXAEuBiSUua\nDrsEeC4izgKuI8tVSzpuFXA2sAL4cnq/huYctz1xkJxljUDZKmnvdIl8G691My2xH0mCm9+jzGeW\nOdaGUx+nJS4DdkbEYxGxD1gPrGw6ZiVZTlrIctQuTzlrVwLrI2JvRDwO7Ezvd1iO235wm+QQaFXl\nbqRgm+7c/P9Fr7U7r1fN79PJZ7rUPFoi4ED/EuouAJ7MPR8H3tTqmIg4IOkFsnwQC4B7m85t5K1t\nznHbM5ckzay0Dqrb8yVty21rmt6qTO7ZVscU7m+R47ZnLkkOieZB5mbDpsO527sjYmmb18eBhbnn\nRblnG8eMS5oLHA/saXPuu2nKcSvpryPi/WUvuohLkkMm30bZrN2+2WqT7OQzeznWhkOESm0lbAUW\np9VV55F1xGxsOmYjWU5ayHLU3h0RkfavSr3fi4DFwP0tctz2FCDBJcmhdHApiBJteLPdJtnL+7rE\nPHr6leAitTFeBmwGxoB1EbFd0lXAtojYSJaL9mspN+0essBHOm4D8DBwALg0Iib6cmEFegqSkk4g\n60U6h+xX+kPADuA24EzgZ8DvR8RzPV1lDQ0yGYZZNyL6O5g8IjYBm5r2fTr3+GXgfS3OvRq4us17\n30PKcdurXqvbXwC+HRH/AvhXZGOTLgfuSutu35WeW5cm52ZbJ1XwMq+ZdU5MTM4ptVVJL6slHge8\nlbQ8Q0Tsi4jnmTq2qfLrbpvVSR/bJEdGLyH/NcCzwF+lKUA3SjoWODUidgGk/08pOlnSmsbwgImX\nXuzhMqptzoFsg8PHUzb2QXGp0W1+1k+Nudv9mpY4KnoJknOBc4EbIuINwIt0ULWu0rrbM2G63JMO\niDZwkbVLltmqpJcgOQ6MR8R96fntZEHzaUmnAaT/n+ntEq2heQpjw0y3V7b7bA8BqrY6rpbYy7rb\nvwSelPS6tGs5WZd8fmxT7dbdHrRO0qENaghQc0KObhN0uPQ7WqKmHTe9jpP8KPD1NBj0MeCDZIF3\ng6RLgCdo0YVvZqOnalXpMnoKkhHxIFA09Wh5L+9r7RVNYXSpzGZC1Xquy/CMmxHWPDNnGDKTW3Vl\nnTIOkjZiNHkoUOZTrHXKAdXKqNrwnjIcJCugebnag/unCXqd5K00A7dJ2ggrWq52utJh/rV2pdDm\n98k/n+4zOznWhlsgJivWc12Gg6SZlVbHv2kOkhXSKJVNHhHM2a9sMNbE9G2VIUBTMw+16zlvNy6y\nl2NtyLnjxqpizv7sB1kpw16ZdXIOy5vvAGZFavhzUb8Ghhp5+VVZlJwydbBROjwAc/ZlWwgmjp76\n099uqmE32dBrWACppDpmAXJJssKOenaMGANSNboxVKgxvjJyKxXP2Vf+B7ubbOgumY6+ACYnqxUA\ny3CQrDjl2yQncyW6OYcf1/I9HOAMsihZsVJiGQ6SZlZaHcdJuk2yBhqDxWMsV1WeyNollRL6Ts6L\n0u2JnaRK8/ISFRMltwpxkKwRTcCBo7NgOHFU8JvTJ/jN6RMo4Jhfakq+yrJDd6ZLlTYbKzbaoJTr\ntCnbcSNphaQdknZKOixhd1oy9rb0+n2Szsy9dkXav0PShWnfQknflfSIpO2SPtaPu3Z1u2bm/kZZ\nEJwQ85479Ddy3yuz/x24rK0+/XxIGgOuB84nS+C9VdLGiHg4d9glwHMRcZakVcC1wB9IWkK2vOzZ\nwOnAnZJeS7a87J9HxAOSXgn8UNKWpvfsmEuSNaSA/cdNsvdVE+x91QQEHPfzrAu8k0zjZYYAdTNc\nyIZUQEyq1FbCMmBnRDwWEfuA9WSLCOblFxW8HVguSWn/+ojYGxGPAzuBZRGxKyIeAIiIX5Ot3rqg\n19vuKUhK+rNUrP2xpFslHSVpUSoaP5qKyvN6vUjrv3nPz+F77/o833vX59EEHPc39xYe12nJ0tXt\nqlPJbVoLgCdzz8c5PKAdPCYiDgAvACeXOTdVzd8A3EePellSdgHwJ8DSiDgHGCMrAl8LXJfW3X6O\nrMhsZlVQvuNmfmM11LStaXqnokh62MSvFse0PVfSK4C/Bf40Iv5puluaTq9tknOBoyXtB44BdgFv\nB/59ev1m4DPADT1+jg3Aiq98AoA4An7+39/MnAPFGXpaZevxejY1VP77uDsiilYtaBgHFuaenwE8\n1eKYcUlzgeOBPe3OlXQEWYD8ekR8o/TVttHLQmC/AD5Hto7NLrKi8A+B51PRGIqL0DZkFPCTP/ry\ntO2PRe2LnQwB6mVlRRsCjcHkZbbpbQUWp+a5eWS10I1Nx+QXFXwvcHdERNq/KvV+LwIWA/en9sqb\ngEci4vO933Cml+r2iWQNqIvIepiOBS4qOLTwb4+kNY2i+MRLL3Z7GdYnS274Y6D4Z7wxzKdoOmIn\nQ4B6WVnRhkO/1t1OBanLgM1kHSwbImK7pKskvTsddhNwsqSdwMeBy9O524ENZKuzfhu4NCImgLcA\nfwi8XdKDaXtnr/fcS3X7HcDjEfEsgKRvAG8GTpA0N30RiorQAETEWmAtwFGnL/SvyxAoStxrNkUf\n525HxCZgU9O+T+cev0yL1VYj4mrg6qZ9f0fJXqNO9NK7/QRwnqRjUjG3se72d8mKxuB1t80qJV+r\naLdVSS9tkveRjV16AHgovdda4JPAx1MR+WSyIrONiINTGEuOj/S0xBop27NdsSDZ67rbVwJXNu1+\njGygqI2wRjq1okXGus1MXvQZNkpKd8pUiqclWkutVmG0Gqvhz4KDpLWVL1GW0UlQdQAeQSV/DqrE\nQdKmlU/W22iv7EcGcgfIEeOku2Zm7dXxD5uDpJUy3RRFq4kafu+dKs060mijbMVDgKxqXJK0jmly\n6vCgg/s7HEjsEunoqeP3zEHSulJ2eFAdf6kqK+jrtMRR4SBpXZtS9W5Rimw3F9ztmyOoht8vB0kz\nK62Of9QcJK0njUHmk3Oz5WmbB55rkoN5WZpLlXX8hRt5NfyeOUhaX8xJaZYPm5nT1LFjI66G30MP\nAbK+2n9sTO3xTsHTw31GX9k0aVX7Y+iSpPXVES/qYNUbYHIeLTt1bAS5d9usd42qN5RPjGGjoY5/\n7KatbktaJ+kZST/O7TtJ0pa0tvaWtN4NynxR0k5JP5J07iAv3sxmWA2T7pZpk/wqsKJp3+XAXWlt\n7bvSc8gWAluctjV4Kdna2/+KOJS9vIK/QLVS0zbJaYNkRHyfbK3bvJVka2qT/n9Pbv8tkbmXbFGw\n0/p1sTZ6jvh/gsagcTGAZZpsRrkkWdqpEbELIP1/Stq/AHgyd5zX3bbD5njb6NJkua3Ue0krJO1I\nzXOXF7x+pKTb0uv3SToz99oVaf8OSReWfc9u9HsIUNGvQuHfFa+7XS+NKYztMghZfUgaA64na6Jb\nAlwsaUnTYZcAz0XEWcB1wLXp3CXAKuBssqbAL0saK/meHev2R/bpRjU6/f9M2j8OLMwd13bd7YhY\nGhFLx445tsvLMLMZ1b/q9jJgZ0Q8FhH7gPVkzXV5+Wa924HlafnqlcD6iNgbEY8DO9P7lXnPjnUb\nJDeSrakNU9fW3gh8IPVynwe80KiWmzWqYo3SpKvgI6azjpv5jZpi2tY0vVuZprmDx0TEAeAFsmWq\nW507kOa+acdJSroVeBvZTY+TLSF7DbBB0iXAE8D70uGbgHeSRfaXgA/2eoFWPY02q6r1gtZC+e/Z\n7ohY2ub1Mk1zrY5ptb+o0NfzT9m0QTIiLm7x0vKCYwO4tNeLsnqIMdDEbF+FdaR/f9jKNM01jhmX\nNBc4nmykTbtzSzX3dcLN6DZrHCBHi+hr7/ZWYLGkRZLmkXXEbGw6Jt+s917g7lQQ2wisSr3fi8jG\nZd9f8j075mmJZlZOHweKR8QBSZcBm4ExYF1EbJd0FbAtIjYCNwFfk7STrAS5Kp27XdIG4GHgAHBp\nREwAFL1nr9fqIGlm5fWxHTkiNpH1Y+T3fTr3+GUO9Xc0n3s1cHWZ9+yVg6SZlVfDzjYHSTMrrY4j\nEhwkzaw8B0kzsxainvlBHSTNrDyXJM3MWnObpJlZOw6SZmYtVDChbhkOkmZWinB128ysLQdJM7N2\nHCTNzNqoYZDsdt3tv5T0k7S29jclnZB7rXCBHrN+a2Q2d4bzGeIlZVv6Koevu70FOCcifhv4R+AK\naL1AT9+u1iyn8ctYtV/KoeYlZQ9XtO52RHwnrTkBcC9ZBmBovUCP2cC4JDlz+rmk7KjoR2byDwHf\nSo+97rbNOEUWKF39Hrw6Vrd76riR9CmyzMBfb+wqOKzlutvAGoC5x5/Yy2WY2UyoYFW6jK6DpKTV\nwO8By9O6E9DhutvAWoCjTl9Ywy+99VO+9FK1ksxQqeHXtqvqtqQVwCeBd0fES7mXWi3QYzZjXN0e\njMaMG1e3m7RYd/sK4EhgiySAeyPiP7VboMdsplTtl3SYaLJ+X9xu192+qc3xhQv0mNmIm6E2SUkn\nAbcBZwI/A34/Ip4rOG418F/T0/8RETen/W8kG7p4NNmiYB+LiJD0l8C/BfYBPwU+GBHPT3c9Xnfb\nzEqboer25cBdEbEYuCs9n3odWSC9EngT2TDDKyU1eoBvIOsUXpy2xjjvwvHd03GQNLPyZmYw+Urg\n5vT4ZuA9BcdcCGyJiD2plLkFWCHpNOC4iPhB6lC+pXF+m/HdbXnutpmVNkPtvadGxC6AiNgl6ZSC\nY1qNyV6QHjfvb/Yhsir9tBwkzay88kFyvqRtuedr07A/ACTdCby64LxPlXz/VmOypx2rXTC+uy0H\nSTMrp7PVEndHxNKWbxXxjlavSXpa0mmpFHka8EzBYeNko24azgDuSfvPaNp/cKx2i/HdbblN0sxK\nmcFxkhuB1enxauCOgmM2AxdIOjF12FwAbE7V9F9LOk/Z+MQPNM5vM767LQdJMysvotzWm2uA8yU9\nCpyfniNpqaQbs8uIPcBnga1puyrtA/gIcCNZgp2fcii3xJeAV5KN735Q0v8uczGubptZaTPRcRMR\nvwKWF+zfBnw493wdsK7FcecU7D+rm+txkDSzcpzgwsysvarliizDQdLMSnOQNDNrJehHp8zIcZA0\ns9LqmGHJQdLMynOQNDMr1hhMXjddrbude+0/SwpJ89NzSfpiWnf7R5LOHcRFm9ksiECT5bYq6Xbd\nbSQtJBsN/0Ru90UcyuG2hiyvm9nI8lIQTbzu9uGK1t1OrgM+wdQvyUrglsjcC5yQJqibjaQ6Vi/b\nqeMaN90uBPZu4BcR8Q9NL3ndbbOqCmAyym0V0nHHjaRjyHK+XVD0csE+r7ttlRCqXimpYzW8/25K\nkr8FLAL+QdLPyPK1PSDp1XS47nZELI2IpWPHHNvFZZjNrNoHSOpZ3e64JBkRDwEH06mnQLk0InZL\n2ghcJmk92QI9LzTSsJvZ6Ktaz3UZZYYA3Qr8AHidpHFJl7Q5fBPwGFket68Af9yXqzSz2Ve2Z7ti\ncbTbdbfzr5+ZexzApb1flpkNm2wwecUiYAmecWNm5TkLkJlZay5Jmpm1UsH2xjIcJM36pGgKY7WG\nw1RvXnYZDpJmfdIcECs577uG1W0vKWvWZ6GqBshs+YYyWy8knSRpi6RH0/+FU/IkrU7HPCppdW7/\nGyU9lLKRfTGtv50/b0r2suk4SJpZeTOz7vblwF0RsRi4Kz2fQtJJwJVkk1aWAVfmgukNZFOeGxnJ\nVuTOK8pe1paDpFmfVXFq3kEzM5h8JXBzenwz8J6CYy4EtkTEnoh4DtgCrEhZx46LiB+kcdu3NJ1f\nlL2sLbdJmg1IFQOlJkvXpedL2pZ7vjYi1pY899TGdOaI2CXplIJjWmUcW5AeN++fkr2sqQbeloOk\n2YBVJntQ0Mlg8t0RsbTVi5LuBF5d8NKnSr5/q4xjhfunyV7WloOk2YBVIkACIvo2mDwi3tHyc6Sn\nJZ2WSpGnAc8UHDYOvC33/AzgnrT/jKb9TzE1e1lj/wOSlkXEL9tdq9skzay8mem42Qg0eqtXA3cU\nHLMZuEDSianD5gJgc6qm/1rSealX+wPAHRHxUEScEhFnpnwT48C50wVIcJA0s07MTJC8Bjhf0qNk\nPdHXAEhaKunG7DJiD/BZYGvarkr7AD4C3EiWjeynwLd6uRhXt82snM7aJLv/mIhfAcsL9m8DPpx7\nvg5Y1+K4c6b5jDPLXo+DpJmV1kHvdmU4SJpZSX2pSo+cMpnJ10l6RtKPm/Z/VNIOSdsl/UVu/xVp\nOtAOSRcO4qLNbBYEM9UmOVTKlCS/CnyJbOQ6AJJ+l2xU/G9HxN7GYE9JS4BVwNnA6cCdkl4bERP9\nvnAzmwX1q21PX5KMiO8De5p2fwS4JiL2pmMa45hWAusjYm9EPE7Wu7Ssj9drZrNIEaW2Kul2CNBr\ngX8j6T5J35P0O2l/q6lCh5G0RtI2SdsmXnqxy8swsxnl6nZH550InAf8DrBB0mtoPVXo8J3ZPM61\nAEedvrBaX1WzKoqAifrVt7sNkuPAN1KWjfslTQLz0/6FueMaU4LMrAoqVkoso9vq9v8B3g4g6bXA\nPGA32XSiVZKOlLSILJfb/f24UDMbAq5uH07SrWQTyedLGidLdLkOWJeGBe0DVqdS5XZJG4CHgQPA\npe7ZNquIALzGzeEi4uIWL72/xfFXA1f3clFmNowCwm2SZmbFAnfcmJm1VbH2xjIcJM2sPAdJM7NW\nqtdzXYaDpJmVE4BTpZmZteGSpJlZK56WaGbWWkB4nKSZWRs1nHHj1RLNrLwZmLst6SRJWyQ9mv4/\nscVxq9Mxj0pandv/RkkPpRUSvpiWlm28VriiQjsOkmZWTkTWu11m683lwF0RsRi4Kz2fQtJJZHkk\n3kSW2PvKXDC9AVhDlmBnMbAinZNfUeFs4HNlLsZB0szKm5ksQCuBm9Pjm4H3FBxzIbAlIvZExHPA\nFmCFpNOA4yLiBynpzi2581utqNCWg6SZlRTExESprUenRsQugPT/KQXHtFoFYUF63LwfWq+o0JY7\nbsysnM5Spc2XtC33fG1ajQAASXcCry4471Ml37/VKgjtVkcoXFEhlThbcpA0q6hI4UJx6HHvb1q6\nvXF3RCxt+TYR72j1mqSnJZ0WEbtS9bmoWjxOlue24QzgnrT/jKb9T+XOKVpR4dl2N+LqtpmVEkBM\nRqmtRxuBRm/1auCOgmM2AxdIOjF12FwAbE7V819LOi/1an8gd36rFRXacpA0qyhFtjUe9yxS0t0y\nW2+uAc6X9ChwfnqOpKWSbswuJfYAnwW2pu2qtA+yDpobyZa0/inwrbR/HfCatKLCeg6tqNCWq9tm\nNdCv6nYfOmWm/4yIXwHLC/ZMZZ8RAAAC5ElEQVRvAz6ce76OLPAVHXdOwf59tFhRoR2VCKQDJ+lZ\n4EVKFH0rZD71ul+o3z0P2/3+84h4VbcnS/o22T2VsTsiVnT7WcNkKIIkgKRt7Rp6q6Zu9wv1u+e6\n3W9VuU3SzKwNB0kzszaGKUiunf6QSqnb/UL97rlu91tJQ9MmaWY2jIapJGlmNnRmPUhKWpHyu+2U\ndFhKpKqQ9LOU4+7BxpzWsnnzRoGkdZKeSQN1G/sK70+ZL6bv+Y8knTt7V969Fvf8GUm/SN/nByW9\nM/faFemed0i6cHau2jo1q0FS0hhwPXARsAS4WNKS2bymAfvdiHh9bljItHnzRshXSXn7clrd30Uc\nyvW3hiz/3yj6KoffM8B16fv8+ojYBJB+rlcBZ6dzvpx+/m3IzXZJchmwMyIeS6Ph15PlkquLMnnz\nRkJEfB/Y07S71f2tBG6JzL3ACSmRwUhpcc+trATWR8TeiHicbMrcsoFdnPXNbAfJVjnhqiiA70j6\noaQ1aV+ZvHmjrNX9Vf37fllqRliXa0Kp+j1X1mwHyXa536rmLRFxLllV81JJb53tC5pFVf6+3wD8\nFvB6YBfwP9P+Kt9zpc12kBwHFuae53O/VUpEPJX+fwb4JllV6+lGNbNN3rxR1ur+Kvt9j4inI2Ii\nsrVXv8KhKnVl77nqZjtIbgUWS1okaR5Zw/bGWb6mvpN0rKRXNh6T5b77MeXy5o2yVve3EfhA6uU+\nD3ihUS0fdU1tq/+O7PsM2T2vknSkpEVknVb3z/T1WedmNVVaRByQdBlZAs0xYF1EbJ/NaxqQU4Fv\nppUt5wJ/ExHflrSVLIX8JcATwPtm8Rp7IulWskzR8yWNk61kdw3F97cJeCdZ58VLwAdn/IL7oMU9\nv03S68mq0j8D/iNARGyXtAF4GDgAXBoRg887Zj3zjBszszZmu7ptZjbUHCTNzNpwkDQza8NB0sys\nDQdJM7M2HCTNzNpwkDQza8NB0sysjf8PkgcWMPD86KAAAAAASUVORK5CYII=\n", "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -212,7 +212,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The final function of the covariance module is the ability to sample a new set of parameters using the covariance matrix. Currently the sampling uses np.multivariate_normal(). Because parameters are assumed to have a multivariate normal distribution this method doesn't not currently guarantee that sampled parameters will be positive. " + "The covariance module also has the ability to sample a new set of parameters using the covariance matrix. Currently the sampling uses np.multivariate_normal(). Because parameters are assumed to have a multivariate normal distribution this method doesn't not currently guarantee that sampled parameters will be positive. " ] }, { @@ -225,19 +225,19 @@ "output_type": "stream", "text": [ "Sample 1\n", - " energy J neutronWidth captureWidth fissionWidthA fissionWidthB\n", - "0 0.031860 2.0 0.000475 0.106353 0.0 0.0\n", - "1 2.823805 2.0 0.000326 0.103970 0.0 0.0\n", - "2 16.242208 1.0 0.000430 0.110773 0.0 0.0\n", - "3 16.770072 2.0 0.012445 0.083466 0.0 0.0\n", - "4 20.559242 2.0 0.011852 0.082108 0.0 0.0\n", + " energy L J neutronWidth captureWidth fissionWidthA fissionWidthB\n", + "0 0.032744 0 2.0 0.000476 0.104726 0.0 0.0\n", + "1 2.824910 0 2.0 0.000367 0.093434 0.0 0.0\n", + "2 16.246703 0 1.0 0.000428 0.128670 0.0 0.0\n", + "3 16.772205 0 2.0 0.012712 0.087642 0.0 0.0\n", + "4 20.558978 0 2.0 0.011754 0.088571 0.0 0.0\n", "Sample 2\n", - " energy J neutronWidth captureWidth fissionWidthA fissionWidthB\n", - "0 0.032548 2.0 0.000476 0.105254 0.0 0.0\n", - "1 2.829121 2.0 0.000375 0.086731 0.0 0.0\n", - "2 16.234402 1.0 0.000438 0.144065 0.0 0.0\n", - "3 16.771254 2.0 0.012171 0.088258 0.0 0.0\n", - "4 20.571634 2.0 0.011214 0.094332 0.0 0.0\n" + " energy L J neutronWidth captureWidth fissionWidthA fissionWidthB\n", + "0 0.030677 0 2.0 0.000474 0.108886 0.0 0.0\n", + "1 2.822705 0 2.0 0.000357 0.099604 0.0 0.0\n", + "2 16.250280 0 1.0 0.000527 0.128495 0.0 0.0\n", + "3 16.770029 0 2.0 0.013519 0.076974 0.0 0.0\n", + "4 20.555458 0 2.0 0.011091 0.094415 0.0 0.0\n" ] } ], @@ -252,6 +252,82 @@ "print('Sample 2')\n", "print(first_five_sample_2)" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can reconstruct the cross section from the sampled parameters using the reconstruct method. This method also required the equivalent openmc.data.IncidentNeutron.resonance.ResonanceRange object. " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ]" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "gd157_endf.resonances.ranges" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Reconstructing in the ReichMoore region using our previously generated samples. " + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0,0.5,'Cross section (b)')" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XmcW3W5+PHPc7JMZutM95YutLSU\ntQhYlrJWdkEW2VQUERFERa8/dy9eFZeL+wUFZVFEvIJXEGWRVXYFpJTN0kJpS+lK95a2syQ55/n9\nkZOZJJNMlkkmyfR5v15hkpNzvnk6LefJdxdVxRhjjMnkVDsAY4wxtckShDHGmKwsQRhjjMnKEoQx\nxpisLEEYY4zJyhKEMcaYrCxBGGOMycoShDHGmKwsQRhjjMnKEoQxxpisgtUOYCBGjRqlU6ZMqXYY\nxhhTV+bNm7dBVUfnO6+uE8SUKVN4/vnnqx2GMcbUFRF5q5DzrInJGGNMVpYgjDHGZGUJwhhjTFZ1\nmSBE5FQRuWHr1q3VDsUYY4asukwQqnqPql7S1tZW7VCMMWbIqssEYYwxpvIsQRhjzCDavCPKmq2d\n1Q6jIJYgjDFmEJ39/d9zyQ9+Xe0wClLXE+WMMabePBL+gv/ss1WNoxCWIIwxZhBFtwfw4kKk2oEU\nwBKEMcYMoiX3jgVgryrHUQjrgzDGGJOV1SCMMWYQzTvgC2wdNtVqEMYYY9JtbZsGkrj1/n3BWq59\nbHGVI8rNahDGGFMln7glsV3BZ94zvcqRZGc1CGOMqZJA4zKCbbW7p40lCGOMqZKmKdfRuMsd1Q4j\nJ0sQxhhjsrIEYYwxJquaSRAispeIXCcid4jIp6odjzHG7OwqmiBE5CYRWSci8zOOnyQir4vIYhH5\nGoCqLlTVS4FzgVmVjMuYUkU746xfvo3Vb2zh7aVb2bG1G1WtdljGVESlh7neDFwD3JI8ICIB4Frg\neGAlMFdE7lbVBSJyGvA1/xpjqm7H1i6euf1frHpxGV3dw4gH+25SFZAuJs0cyVEffBetI+phhR1j\nClPRBKGqT4rIlIzDBwOLVXUpgIj8ETgdWKCqdwN3i8jfgFsrGZsxuSxdtZQnbr6L6JuNxAPTQAKE\nou0073iDbllOV2ADMaeTlliAkR2jiIb35K2X9+GWl59kryNHM+e8/XEcqfYfw9SBgKsE3WpHkVs1\nJspNAFakvF4JHCIic4AzgQbgvlwXi8glwCUAkydPrlyUZqegqqx4Zzn/+ve9rPjba0TWTofgXqjz\nLiKx9TTu+DsdY94gfMxY2qbsQVPzoXRrnE071vHahvncuul54lv/wTnPtzJm21ks/MeBrFrwMB/4\n5jGEIzYP1fTve7d4TF0rcHG1I8muGv+Cs321UlV9HHg838WqegNwA8CsWbOs8dcUxVOPJVuW8MLy\nx3lh4YME58LEDQcioXfRENyLkG6mqeOftO3ayWGf+Cgjpn6g3/Jcz+WpFY9xzej/Zuvqm/nM44tY\npefyh8sf4vwfnEgwFBikP5mpR1vHfpbH99yDfaodSA7VSBArgUkprycCq6sQh9kJdLvdvL7pdV5c\n8RTPr3icpSsWs9+iCUzb+C6mh88j2tCOhDppiS9i3F4jmX3xWTQPO6vg8gNOgDm7HscRk+Zw49yf\ncWXTLXzxgS7W81H+74q/c953T0DEmptMdluG71HtEPpVjQQxF9hdRKYCq4APAucVU4CInAqcOn16\nba5fYqqjK97F4i2LWbBhPq+ueoYFG17lzR1rmbYqwP7Lp3PwO/sxs/nDxMLDiDbHadI1TN6/laM+\nciyNLacM6LODTpBPHfIVZoyeyTfky3z+gXY2yBk8cvNzHHfhIWX6ExozuCqaIETkNmAOMEpEVgLf\nUtXfiMhlwINAALhJVV8tplxVvQe4Z9asWTXacmcqJe7FWd+xntU7VvPm1jd5c9PrLN2wgDe3LWd1\ndAvhqDL1beHAFbtw5ua98Jz3886waXiBMNuHdROJbGPG7BYOPX1/Is3hssd37G7vJXRCiG92f55L\nnprKomf2Yb9jNjNm1+Fl/yxjKq3So5g+lOP4ffTTEW12PnEvzpbuLWzu2symrk1s7trMxq6NbOzc\nwOoty1izbTlrOtayLvoOLomup+ZOZep6jwPebuLojbsS7j6ErtAUtrVOIR5qYtNICHrv0D7OZZ9j\np7DPYVMIhCo/N/SoKcdxyRGX8eiamzlo7eXcc+WjfPyXZyI2ssnUmbocZmFNTINPVYlrHNdz6Xa7\nibrR9J9etOd56vGueBfbY9vZEdvB9uh2dkS3saN7Kzui7yRex7azJbadrfEO+ow4UKWtS5mxyWP3\nzS7HbAnTvm0cwa4JxLwJdDWMZ3vzeLojI9jYDqhH0OmgdXSA6YeMZ58jdqO5raEKvy0474BLmf/2\nc8T/cCfx0Pk8cetzzPmINTWZ+lKXCWKgTUx3vnEnT69+umcGbO9M2ORrz3+V/r6iWc7VnOeml+3l\nOd//r6a92/tTs5WZ4zOznOd6Lp66xNXF81xcdXHVw1MvcUw9XPVw/edx/z0XxVW37827SEFVmj1l\nWNxjRKfHiG6PCR3KiB0wvENp2xGiqauVQLQdLz6cmDeSqIygOzyCroYRdEVGsD3cyvb25B/JJRTo\npnF4kF1ntLHXoVMZO6WNUEPtjBr6+rFXc/6yObz/H0t5/fHRHHZ23Ia+mrqyU/5rXfvSLSza8hrQ\nO+ZW/FtgZiNAz/va+1oy3us9N72MbNdmLTtfDF7KZ2t6DE7GndvR9DIcL/EIKIRdJaAQ9CDoJZ6H\n3EQZPce8xPPksYCnOJ4k3nMhGBca4krIhVAcAvEA4kYQrwG8BtQLg4ZxvQiu24irjYnnRHClAc9p\nIB5sJB5sJBZqIRoeRkeohW2RRsichKwuwUCUcLND+8hGRk4eweQ9xjJ6YivDRkVwAjWzlFhWreFW\nvnryj7nn1f9hQvz/8befPcT7//PkaodlTMF2ygQx9aEZfGr5BLJPyUgc67nvJocoavK5gqbf2pMv\nRUH7DGns7zMk7QzteVd6PyJHOkr/nOzxgKASQCXjJw4qDkjiZ58H6a9dcYg5QTwnhNfzM4QXCaJS\nxDd2VUTiOI7iBJRQRIi0hBk5oon2scMYPm4Yre0RmlrDtI6M0NgaqvshorN3PYbb51xL659fZu2b\nM+ja3k2kpTrNXsYUqy4TxED7IDaMOoKtgZbyBlU22Rpz+mvgSX+v7+1U/RyXqHqIX50R8fOdJLbH\nFRH/ufS8FofEkhGOEAoIwXCQYEPiEYoECUfCNDSHaYiEaYgECIUDBIIOwZBDOBIk1BggHAnS0Bgk\nFAkQagjU/Q2/FF858Rf89J8fZ0rnu3jw6gc4/fLTqx2SMQWpywQx0D6Ic791Ml68b7tPz71L0r/Z\np58jaccyzxGkz5f+3nNSyk2+txPeMHc241p3YexJu9N0y0LeXjqBWDROKFyX/+uZnUxtN+JWSDgS\nJNISSjyaex8NTf6jMfGtN5x8RHofoYbEN+FQOPEIJh+hxCMQcggE/Ucg8XCSD0dwHEl8Sxex5LAT\n+cRR32HxqL8TDw3j4WseqnY4xhRkp0wQxgy29sbhjD9xMk3bl7F6/jbbQ8LUBUsQxgySC479Luua\nnqI7PJp5d82tdjjG5FWXCUJEThWRG7Zu3VrtUIwpWGvDMJgDwdh2Ftz9YrXDMSavukwQqnqPql7S\n1tZ3dy9jatlH33cFcfcZtoV2Y+3StdUOx5h+1WWCMKZejR82mc17vAEIT1xzT7XDMaZfliCMGWTv\nP/+LNO5YyJatI3Bdr9rhmEGgqlz72GLWbeuqdihFsQRhzCDbZ+JsNrW/QCzUznO3PlrtcMwgmL/q\nHX784Ot8/o8vVTuUotRlgrBOalPvJp26J4F4B288WtRWKKZOxb1ETXFHd7zKkRSnLhOEdVKbenf6\nMZ/Di89jR3AP3llnX3SGunqdFFuXCcKYehcJNRKdshwvEOaxn99W7XCMycoShDFV8t6LLiXctZaN\nK2tnDwtjUlmCMKZK9ph0ELHg83RGprH8hdeqHY4ZBPW2wIolCGOqqOXQRgD+edNfqxyJKUVsyUI2\nX/O9vOfVZw9EnSYIG8Vkhoqzzv8a4c4ldGyfUO1QTAne+vC5vH3NH3DXr+73PE9dGsbdSUzWZX3/\n6P87mqdXP12JEAekLhOEjWIyQ0VzQwvR5gV0RSbwyt0PVzscUyS3I5Z4Eo/2e97Sd14jPPw51jf+\nNuv7m7o2cdW8q8od3oDVZYIwZiiZfMreoB6v3P1MtUMxRVo+/j08OudaXLew3oVcZ50x/z9o3Tym\nfIGViSUIY6rslJMvItT1Bl3x3fE8W3qjniyddAoAHR39L6HhdSmXPnM1e6zbN+v747btxu7zjyx7\nfANlCcKYKnPEwWtfQndkLE/fcmu1wzFFUL9OsLFje7/nOcvfBuDdK/aoeEzlZAnCmBqw/wePRzyX\nJU8uqXYoQ9PcX8O326BjU0WKV6//JqaGNYlO7Lbt9VVDtARhTA04dPaJhLpfI6r74MXra72euvD8\nzYmfW1dWpPh6m99QKEsQxtSKUauINozgkeuuq3YkQ053V5ylK5qhUmsiDdE9xi1BGFMjDvv4BxAv\nxqp5m6sdypDz9MPb6P5nGytfn1/egv3EkD8/SMbP+lCXCcImypmhaJ99DiDcvYCo7EO8o6Pa4Qwp\n3rZE2//idesr9AlWg6gZNlHODFUybgOxcDsPXffLaocypEjy/i3lXRgxWR/w8iSI3pYtq0EYY0p0\n5MXnI26Uta90VjuUIUmcytzyNE8bkzj5E4NTgwOcLEEYU0NmTN+TcGwBMWcm3dvfqXY4Q4eC64QR\nKfctz08MBU9wzJ1I2jf3v1xHNViCMKbGOOO3EAsP48Frf1HtUIaMjaOO54mj/gc3FqxI+YUPYspd\nk3Dc2qtCWIIwpsYcc8nHcdwoGxfU3g2jXm1tnwVAPBquSPmeluPvqvb6JyxBGFNjpkydQii2gGhg\nJp1bKzPzd+eTuPlWqg+iPKOYam8klCUIY2pQaMI24uFh3H+NNTOVk1NAZ3FxEjd1L89SG4UVVYcJ\nQkRmi8i1IvKKiKwXkeUicp+IfEZEbJypMRVw7Kc/ieN2s2WR7VddHn4NolI34bzl1l7zUSH6TRAi\ncj/wCeBB4CRgPLA38A0gAtwlIqdVOkhjdjYTJ4wjFHuVWHAm2za9Xe1whoDEDVqlWgki2yW1V2PI\nlK8Gcb6qXqSqd6vqalWNq+p2VX1BVX+qqnOA2tsnz5ghoGFSlHiolfuvubbaodS/5Ey1snQmlyDL\nnTbfCrC1oN8Eoaobks9FZJyInOYvczEu2zmDxZbaMDuD4z5zCY7bxfYlkWqHMnSoW+by/J955kFI\nliametgcqqBOahH5BPAccCZwNvCsiHy8koH1x5baMDuD8eNGJSbNBWeyed1b1Q6nrmmyiancCSJZ\nMSm4ZtKbKLQGRy1lKnQU05eBA1T1Y6p6AfBu4KuVC8sYAxCZ4hEPtfDAtb+qdih1LnknL3OxydVc\n89UGsiwz7sZjZQ6m/ApNECuBbSmvtwEryh+OMSbV8Z+6GMftpGPZsGqHUuf8GoRX5hqEz80zCzrb\nGKZKxVJO/c47F5Ev+E9XAf8SkbtI5ODTSTQ5GWMqaMzY4T3NTGtXvsbYiXtWO6T65N+hPbe8VQhJ\nVknyNjENzT6IVv+xBPgrvRW0u4A1FYzLGONrmhbEDTXz8PW/qXYodayyNQivhHWUyjK5rsL6rUGo\n6hWDFYgxJrvjPnURf/7ak3S/NbzaodSxRIIoz5pJfSn5+iD6PMneb9G5BRrbyxbXQOWbKHeDiOyb\n471mEfm4iHy4MqEZYwBGjxpGOL6AWHg/Vr71UrXDqVPJGkRlSs/bdJWlEyLryKeHv1megMokXxPT\nL4FvishCEbldRH4pIjeJyFMkJsi1AndUPEpjdnLNM5pwg008ev0t1Q6lLvXcvsve7p/ck7r4pTY8\nN0tzl1tbe0Lka2J6CThXRFqAWSSW2ugEFqrq64MQnzEGOPaTF3DHV58gunJMtUOpT/4w04p1DJe0\n1Ebtd1IXtHuGqm4HHq9sKMaYXEaNbCXsLiAansniRU8xfcaR1Q6pziQnylWm9Hyd3z0zqVMqEn1r\nELW3oJ8t921MnWjZsx032MiTN91e7VDqV6VGMeXJPNlmTQ+FxfqMMTXiuE9+hEB8B+7qXeri5lJb\nkk1Mlfm95V94TzJ+gufFM0spZ0hlYQnCmDoxor2JcHwhbng/Frz6cLXDqSvau2hShT6g0GGuvR5a\n/khlYimjQhfrmyEiN4rIQyLyaPJR6eCMMemG7TsGNxjhmd/dXe1Q6oskh7mWfTGmxH8LTDzBlEpD\n16vr094ThLc27ihbZOVQUCc1cDtwHXAjUPsLiBgzRB1zyQe4/f89jK6djKoiWRaBM7lVqmku72J9\nWUTmh0gd1BqMw6YX17Jr+cIasEITRFxVK7qcpIicAZwCjAGuVdWHKvl5xtSjEcMaCXmv090wkxde\n/AvvPvDMaodUJ5KLMZV5aGlyKaYCE8+Olgkp1/a9JrBiezmiKptC+yDuEZFPi8h4ERmRfOS7yJ9U\nt05E5mccP0lEXheRxSLyNQBV/auqXgx8DPhAsX8QY3YWbfvtghdoYN6tf692KHVDe5baqE4NorBP\nrd9O6gtI7AnxNDDPfzxfwHU3k9jLuoeIBIBrgfeS2N/6QyKyd8op3/DfN8Zk8Z6LziIY24asm4pb\nB0tG14RkH0TZE0RxfRD9E+L51nQaZAUlCFWdmuWxWwHXPQlsyjh8MLBYVZeqahT4I3C6JPwQuF9V\nXyj2D2LMzmLEsAhh73XiDfvyr2dvrXY49aVaw1yzvR3vTnu5vXUSS8It5QuqDAodxRQSkc+JyB3+\n4zIRCZX4mRNI32xopX/ss8BxwNkicmk/sVwiIs+LyPPr16/PdZoxQ9qwA6biBRp4+Y6nqx1KnUiv\nQTy4cDEPLnyjbKXnbbrK8rZkWXfpnebjyxRReRTaxPQrEtuM/tJ/vNs/Voqsmyup6s9V9d2qeqmq\nXpfrYlW9QVVnqeqs0aNHlxiCMfXt2ItOIxh7h+CGacTitbXAWy1SSd9y9EvPvZ8vPVfGDv48CcIr\n0+BPz9NB3Uei0ARxkKpeoKqP+o8LgYNK/MyVwKSU1xOB1SWWZcxOqb25gQZ9g3jDvjz55K+rHU4d\nqNQ8iIS8C+9p3+/FpUTy249/g99eOHhLgheaIFwRmZZ8ISK7Ufp8iLnA7iIyVUTCwAeBomb9iMip\nInLD1q1bSwzBmPo3+pA98QJhFt5le0QUKnkjn7Zhf6ZtOKAcJSb+20/iWbVoM1vW9N3sKVbC9qdd\nkWPpanxP0deVqtAE8WXgMRF5XESeAB4FvpjvIhG5DXgG2ENEVorIRaoaBy4DHgQWAn9S1VeLCVpV\n71HVS9ra2oq5zJgh5ZiPvZdgbCuRTXuwrTNzLIjJyr8nH//GhRz/xsfKV2w/TUx//dmLbFldn83h\nhS73/YiI7A7sQaKu9pqqdue5DFX9UI7j9wH3FROoMSZdY0OQBuctOiP78OC9P+bsc66sdkg1T1Xp\njpdvaHBPw1Ed7C9dinxbjh7j/zyTxCzn6cA04BT/mDGmiiYffwheIMyaR9ZWO5Qal5yvANF4+eca\nlDYPovaTSr4mpqP9n6dmebyvgnH1y/ogjEk4+sxDCMa20PrOTFasX1DtcGqfaplnUycTTylJJ+uA\nzgFFU279JghV/Zb/9DuqemHqA/hu5cPLGZf1QRgDBIIODY0b6GjZmwdv/1G1w6l9qnR1vVOBcstf\nZC0otJP6z1mO3VHOQIwxpdn7nPegTojuuWHbSCgPVSUe7ahAweUvshbk64PYU0TOAtpE5MyUx8eA\nyKBEaIzp10FHTScUXUtT7GBeXnx/tcOpaeop8XjmTm4DN5iT1wZTvhrEHiT6GtpJ7384ELi4sqHl\nZn0QxvQSERrHxOho2Y1/3Pq7aodTkyRZs1LBjZW/DyLvjnL9XFrL+h3mqqp3AXeJyGxVfWaQYspL\nVe8B7pk1a1bVkpQxtWT2RSfx4E9eIbh4V7rjXTQErYKfnaJu+WsQJbXs1cFeT4X2QVwqIu3JFyIy\nXERuqlBMxpgiTd99DJHomwQCB/P4szdWO5yapQpeV2f5y62xZbrLpdAEsZ+qbkm+UNXNQDnmqRtj\nyqR93zFEIyNY8Gdb4TUnVbwK9EGUVIWogyamQhOEIyI9i4n4u8kVul1p2VkfhDF9nfDJk3DcToat\nm8nbW5dVO5zapOAN8M68cnMH976SXF80UVa5OqlrrdWp0ATxU+BpEfmuiHyHxM5yVRt0bfMgjOmr\ntaWBiLOCWNMB3Hvn96sdTk1SBS1hkbxUZ1z7NJfd+mJmycUXJIWlg79e8xNuuPh7uPFY76eVe2/t\nHArdUe4W4CxgLbAeOFNVf1/JwIwxxZt2+mG4gQY6norjlTS7d4hTxdOBrcW0NfwQLXt+Pe2YlPKr\nLjCnrH5lX2KBw9i4pndXhJu+fHkJH1i8QmsQACOAHar6C2C9iEytUEzGmBIdcdK+hGIbGNY9m2cX\n2FzWPnTge0JExtyPQ/r6SwNttuqPJEfSpib8DbtW7PNSFbrl6LeArwLJtBkC/rdSQRljSuM4Do3j\n4uxonc6zt95W7XBqknoDq0EctfRsPvnsVbie9nZOl2sGez+tTqWt9zQwhdYg3g+cBuwAUNXVQGul\ngsrHOqmNye3oT74X8Vxal+3N5h22b3sqVYUBtt/vve5IAOIp5VR2hZNE4bGYm3mo4gpNEFFN1KcU\nQESaKxdSftZJbUxukycPp8F9Ey9yCPc+YHtEJPTOpC71i7jG43S88ELP69S9wCu7BlaibDd17/EC\nO7gHqtAE8ScRuR5oF5GLgb8DNhvHmBo19ui9iYda2PDQBlvAL5VScif1mh9dwVvnfbjndSyWMp+i\nAr/ju/50HS+/+I+e1/FYrJ+zK6PQUUw/IbF6659JrM/0Tb+z2hhTg0780GxC0Y20b5/NS28+XO1w\naodqyQnz1WfvTnsdj3XTuyf1QAPzpcS28tEZzL1qUc86UvFoSg2ilpqY/CalR1X1yyRqDo0iEqpo\nZMaYkoVCARpHbKejdQ8e/f311Q6ndgxgFFNM0q+Lx1Ju2BW6Y3c3Tun9vGhKDaLGmpieBBpEZAKJ\n5qULgZsrFZQxZuCO+PSpoC6ti3Znw/a3qx1OTdABzIPIzCtu2jf6yvdBeF6NNjEBoqodwJnAL1T1\n/cDelQvLGDNQU3cbRcRdBg2z+evd3652OLVBQd3S2oO6aWXhjPN6XieamPxiB6Gfx03p8+hqnF7x\nz4MiEoSIzAY+DPzNP2ZrMRlT4yYevz+xcCvbH+kiVoVvoDVHvZKXqdgy/GzW7HJ4z+t4NLUmUnyC\nKDanxOO1W4P4DxKT5P6iqq+KyG7AY5ULq382zNWYwhx/zrsJRdfR1n00j86zFfrLOdfMjUd75rWV\nVoEo9CJ/mGusAqvQ5lHoKKYnVfU0Vf2h/3qpqn6usqEZYwbKcRyaJjt0tOzKy398qNrhVJGfGbR8\nK68mvtGXPpNatMCOZr9sjaf3nUS7uor+zGIVsxaTMaYOnfDZU3DcTka9fTCL3n4h/wVDUM+Wo56U\nrUPZTdlXorQiixuJlLmXdrclCGPMQI0Z2UwkuJKu5gP5221VW6W/ynpnUpMyiskrqj8iPQt48WhK\nuSVEVGANInmW56b3QUS7dxT/oUWyBGHMTuCAj5+AitDw4kQ2dex86zOJ3/mgKngp+0G4A6hNpHZS\nlzKKSQquQSSX2khvYop110gNQkR+JCLDRCQkIo+IyAYR+UilgzPGlMf+B00lEltGMHgEd9z1jWqH\nUwWJm6xo+qqosQGMDPJSry0p0fS9/fZXirqZfRDRHGeWT6E1iBNU9R3gfcBKYAbw5YpFZYwpu/En\nzCQWHkbHwy5dsc5qhzOoemsQTlondWa7fjE8d4CjigrtpE5+Xkas8Wh3jjPLp9AEkVxW42TgNlXd\nVKF4CmLzIIwp3knnHkQ4uppW91jufuIn1Q5nkKWs5qqpCaKYb+Hp3+8TTT7+CKMy1SD6+1w3Y4Jf\nLFo7NYh7ROQ1YBbwiIiMBirfAJaDzYMwpniO49C673C6Gsfy1p2v4Q5w45z6Iik/e2/msaISRDo3\ntQZRyVFMyQFYbmYfRI0kCFX9GjAbmKWqMRIbB51eycCMMeV32qdPJBTdwMhtx/DYK7dWO5xB5N+M\nVdJ2lIt2ln6TTWtiKqUGUXATk19LyahBuANIboUqtJP6HCCuqq6IfIPEdqO7VDQyY0zZNUWCNO4S\np7N5KnN/f/dOs1dEz59SnbTVXOPx0tvxvfjAamCa9fbb9+9Dkov1ubXbB/FfqrpNRI4ATgR+B/yq\ncmEZYyrlvV86g2BsG+PWHMGzi+7Of8GQkPy27qStytrVXcxNNnMeRHwgE6mLGObqf17GnI3B2ECo\n0ASRTJWnAL9S1buAcGVCMsZU0qj2JhqGb6Fj2Eweu+mmnaMWkdw/QZ20BZmiXR0lF5lWgyjpd1hs\nE1PGfhTR2kkQq/wtR88F7hORhiKuNcbUmBO+9H4C8Q7Gr5jD80seqHY4g0jSNnbo7iwmQWSMYnK9\nPseKU9wtNHMV2prpgyCRGB4ETlLVLcAIbB6EMXVrl3HDCLeup2PYu/j7b3aG7eV7m5g05aY+oAXv\n3NSZ1AOJqTBen07qMi5Nm0Oho5g6gCXAiSJyGTBGVXfmpSGNqXvHfukMAvEOxq44mnlvPlLtcCqs\nt4kp9Zv4QBJEYtjpADohCGQ5logzvX8hue91eqe4Vyt9ECLyH8AfgDH+439F5LOVDMwYU1m7ThhO\nuHkNnS0zefDGoT3mRCWlBpFyM491FjOjPKOJqcSNh3pLy337TZvzkMxBGX0Q7gBmgReq0Cami4BD\nVPWbqvpN4FDg4sqFZYwZDMd/9RwC8R2Me+sonl50V7XDGQTpw1zdAUw20wFOlJN+br+u23cIrWbs\ndjQYGwgVvOUovSOZ8J8X14BmjKk5k3ZpJ9K+kc7WfXn8ulvwyrnlWgZVLdtmPcXrnUmdtlhfV+lz\nCVITTSkjwTRLE1P2OQ/ZRzG289sfAAAfZUlEQVRlzqyuhEITxG+Bf4nIt0Xk28CzwG8qFlUethaT\nMeVz8tfPIRjbwrj1J/PAvMr9b/3+/7qDD36/yrUUddK2HS2qHT/jK3Ha6qolNTf108SUknwkx3pP\nmYv3VUKhndQ/Ay4ENgGbgQtV9apKBpYnHluLyZgyGTOymcjEKJ3NU5l/01PE3Mp0fh6/YSRHrxpW\nkbLzE/+/AVLbg+LdRfxZM2/QnpLcylRKqhn11wfR9+afOcw1c1RTJeRNECLiiMh8VX1BVX+uqler\n6osVj8wYM2jO+OqZNHSvYUTHKdzx6A+rHU7ZpXZSp1Yh4gNox0/csBO1iFLWPVTJNoop4a23lmb5\nvBpsYtJEg93LIjK54tEYY6qirSlM8wFtdEdGs+62NbzTPdSab5MJIpjWxKSxOHE3xk9/dAavL3m8\nqBITi+f5+0yUdK/Ofvud+8hfeeVHL6V+UuK/GRWGmqhB+MYDr/q7yd2dfFQyMGPM4DrnMycT6VpC\nRN7Lzbf8v2qHUxFKKK2T2ou5vPr4nUSWfo6nrri//4szh+V4SrIGkdmBXBAJZj289Jn5bG2fnnIk\nsR1PZo0hc4e5SsgeYV9XVDQKY0zVBQMOu37oEBb9+W0an5jGotNeZsbYd1U7rDJJ3N1VgiR2LEjw\n4i6b1vjLbTgz85SRngTU85Dk4M6SVvsOZT0ukp6Jog3tAHixjM8vJSkVqd8ahIhMF5HDVfWJ1AeJ\nX8fKikdnjBlUx524Hw2yiGjLbO64+sqKLORXjcUBk30Q6oTSOpu9uIvbcxsscukLTW1iyriVrn4R\nrj8KojtyX+/kWO9UssehGUtrZM6sroR8TUxXAduyHO/w3zPGDDFHXn4eoegmxqw6mQefv6ns5b/T\nUfoKqiVRpbcGEUpLUP/uWkowlGhI6W9mc285vTy3t4kpc/MfffByWPMyrHohd3E5mpgkRxia0Z8e\namjqN9xyyJcgpqjqK5kHVfV5YEpFIjLGVNWMySNomNJNd+NEXr/+hbJ3WK9bubqs5eWVcmNXCab1\n9jbFdsEJ+LfBHN/cc3FTO6m9INde+iiPffNPACx7diML/7gLq5aty3l90TUIN/346V+ofD9RvgQR\n6ee9xnIGYoypHed+7QM0dr1BWN7Hr6/7XFnL3rhqRVnLy097brqeEwKv96v4nmt2QQLJb/J5bocZ\n923XdXvnZ7uJMl5f0wLAltcTTUvrly7PHZWTfZir5KpCpDRjjRz5z97EVkH5PmGuiPRZc0lELgLm\nVSYkY0y1NYaDTLvoaEQ9hs2bzb8W5RnhU4Tt694uW1mFUM/t6UP2nBCasoucSiPqj1HVIvsgNNa7\ncLiqk1bGkrb9eXTOtSzveqf4gHPVZLze2/WYfXctvtwS5EsQnwcuFJHHReSn/uMJ4BPAf1Q+PGNM\ntRx95J4Ehy2nu2kG//zZPXTGi1n5NLeODZvLUk6hPM8leavznCBud+LPIV6caDiSMly0yAQRV8Tv\ne9CeAaGJ11vbjgQgFm0pOt6cLV1eb41DnMFZCq/fBKGqa1X1MBLDXJf5jytUdbaqDu7XAGPMoHv/\nFRfS1LGIRvdUbry+PE1N3e9kG/dSOa4bTVmrzyHur+AajHcSCzb2bryTq2knh7Rhppq4eSdHS/Xc\nvovs11Bga0eO5T+0N0E4gdyzsMup0LWYHlPVX/iPRysdlDGmNoxsjTDxgsMQdWn616E89codJZWT\nuvicu20Au7iV8tmuiyKIl0gMbneixuB4najTiOuvPZVvFJNkTnZIG1WUnNOQXobTz3Ia2XQ1Hkrn\n6tnZ39TeeRM1UYMwxpjjj92PwOjVdDdN5eWr5rGpY0PRZbgpi891d1R+FdK0z/b3bu5JENHeBOEG\nIuzwm5zyfdvXzASRNi3BHyrbU0byZ/nmfKhWflhrJksQxpi8zvv2x2nqeAEJH8+vf/DVoie7xWO9\nG/PEuwf3tuO6MRAH8WdQqz8jWbQTxGH7tuS8jHzfyjPmOqR0GisN/ilO8oD/o3zrJam0pkSyk9Ug\nRGQ3EfmNiJRWhzXGVExzQ5CZX/0Qkc7VtK06lf+78/tFXR+PdffMP/BiDZUIMfdnR/0ahJ8gkqNc\nhcRopq2vvRsAzdsHkVmDkJ75cerk+HZfxtnO3Y2p66UOzq27op8iIjeJyDoRmZ9x/CQReV1EFovI\n1wBUdamqXlTJeIwxpZu1zwQaj2jHcyJsv2sML75ReHdkLNrV8+1aveZKhZhVcu9mwa9B+BPORDL7\nQor8Vu45PaOY3ECOP9MgrLhaSZVOQzcDJ6UeEJEAcC3wXmBv4EMisneF4zDGlMGHPnEKgeZX6W6a\nztNXPsXGHesLuq6zM3V5jeKHfg5Esg8i2auscT8ROOn7UeevQWTQ3vPdQPYaRCErropXQp9MkaOj\nSlXRBKGqT5LYhS7VwcBiv8YQBf4InF7JOIwx5SEifODKz9Hc8QyEj+Tm//ou8QJucLHulG/rMri7\nysV6NgXyE4I/n8AJZsbd/01XJLOJKdBzTeas6J6SCqlBaDT/OX2uGZwFD6vRBzEBSJ1rvxKYICIj\nReQ64AAR+Xqui0XkEhF5XkSeX7++sG8vxpjyaWsMMeuKS2jZ/hqN29/Hr6/5z7zXRLsSCUK8OPHg\ncN7ZtqXSYfZw4/68AvETgp8gAuGM5bPzfSvXzJfZl+v2PyzxUQXsVe021kxXcB/ViCzb34Kq6kZV\nvVRVp6nqlbkuVtUbVHWWqs4aPXp0BcM0xuSy79TRjDn/IMLRzcjLh3HfI/2v+hrzl7eIdC5HnQD/\nuP++wQgTANevQaifIMRf0yjckr6aqjr93fCzDFjVhryDWLWAGsQhJ4/Me05mM5QEhkATUw4rgUkp\nrycCg7y8ozFmoN574kE07NuJSogV/+vw2vKXcp7r+gkiqIn/1de8+MagxAgQj/k1CCdZg0gkgkBL\n/wkhH+l3LVP/owpIECMnhDnnC3v1e446he7tVl7VSBBzgd1FZKqIhIEPAkVtXyoip4rIDVu3DrV9\nc42pLx/+wgVEGv5FPDyRR77zd7Z2ZW86isYSCUKHdROM7cBdl//mWi7JUUwEkh3G/jDbhiJvuhlf\n2pVI7l6L5DwIL0cdI2UTaxGHMTPGZz0t0vli4hzPZdeZS3OFUjGVHuZ6G/AMsIeIrBSRi1Q1DlwG\nPAgsBP6kqq8WU66q3qOql7S1tZU/aGNMwRxHOOdn36R9xyMQPJDf/ueP/Z3W0iXXPyLk0Ni1AE/2\noitlVdVKSs7idgJ+U1MyQTgw9c2/pZ27YfOqfkrKuNlLJHcTU3IHuxzbgjpuNOXU/m73yVnf3ex/\n0tGpJfRzTflUehTTh1R1vKqGVHWiqv7GP36fqs7w+xuKm3FjjKkpzQ1Bjvjxl2nfOpdQx7H8+qrv\n9jnHjSaSgTgQHr4WN9jCnddfPyjxeX4ndSDZxERiox5BaPrAnmnnPnvvXwov12nMO6PZ7c6+cq2T\nsi+2OLlvwyKJZOt43UyctnvBsZVL7Xaf98OamIypLbuOH8GUTx5L8/ZluAsO4v4Hbk17v7cfQDj8\nsouIdK7nnZdaCxrlM1BuLPEZ4nigLip+85YDE2ekJ4h1rxS+mZEXSNQgpJ/Z0u2vtWY97ni9taf+\nFt5Tf7E/8Xd+Dncn1sEapGkQ9ZkgrInJmNpz+JxZtB7sEYx3sez2IG8sW9jznusnCHGESXvuQ6Pz\nD9zwrvzvT6+peFyuPwJIRHG8btRJ9n8Ik/fdL+3c+ObC+0bUCaESJBjPvSnQ8sknZD0uKXMfJMfO\ncgAS3Mr4SSs4/b+OSz9e7KS+EtVlgjDG1KYzP3sRTc3PQKCNR77/UM+ifsmhpslvy4d/5dO0bFvE\njkXTmfvM0xWNSWM9iy8hXjeek1xYTwi3937JDMQ7wN2DaDx730if5b4BlQiOtyPtWNe27eSddJea\nIPLEf+blFzBm6kQAgs2JLUynZCS2SrEEYYwpGxHhfT/+b5o7HsYNzOQvV18NpNQg/PH7u+61F8P3\nW0so3snzv1nLawsW5ixzoFx/uQsRQTSKGwj3vE7V1P0KsYZx3HvrdXnLTDYrqTT2mdW85JUX8weV\nkiD631s6PcYLfv5FLvzREYyfOi3/Z5RBXSYI64Mwpna1NYaYfNlHiHSuYsPL4+iOdvcumJfS3n7q\nV75Ja/sDBDzhsf95g7lPPlOReLxkH4EAdOMGIv7L9Jtv83TFcbtZ92QH2WjK6QE3cY46jWSOblry\n3Es597fuHb3U20nt+J3U5337EBpi/feBOI7QNCzc7znlVJcJwvogjKltcw7dH6/5eWINY7jvqmvx\nYsmhpr23HBHhzJ9cz/DWvxKOdvDcH7bzp6tuzT13oEQ98yAcgGjPpLOeCoQ/LLdl0ihaOufiOe9m\n7r8e7luQ9HaoB+OJZiU3yzLfm5fn3nNb/PkPkrIdneP3QQwf18wZX5+Tdv7grLiUW10mCGNM7Tvw\nU5+moWsdGxc04yV3lMsYsRMKOJz1P7+lbfpzDN/8MutfG8cNn/0db6/aWLY4krOZEzWGlIXx/FiS\nN20nEGDScbsBytzfvNCnHEmZ3+F42wFQv7kqEO/th4hvC/fTsZAsI/vIp1HTd03/zFzFDBJLEMaY\nijhw3xk4vER3wzS6NiduzE6w74gdxxHO/q+fscsF0xmz4Va0ezR3XjGXu355F2584MNgexKEA0is\nz/uOP8rJCTocfcGHGbb9H6i8m8cfvDPjzN7btaPpHdOO2+WXFcNzR5Hru78bTOwb0RXZo7fUQHWW\n0SiEJQhjTEWICExzQBw63vK/rfczKWzOKWdw6C/+mxGhG2jfPJ+Vr7Ry42V38Nq8pTmvKYQmk4xI\nWoJIdlL31iCCiAjTPngwwXgnr92+pp9S0/spHE0kiEjHCuLB8ZBnq1F1evfEcPqZB1FtdZkgrJPa\nmPrw7rM+QiDeRdxNrM8ZyFKDSDVp/BjO+eVfGfWh0YzZ9EtCncojNy7j99/4P7ZtztwBrjC98yBA\nndQEkbj9iSb7RxKxzT7tVFq6H8MN7sUD//e/KSX13si9oEvATZns5k98c3QdXqCRWHBEwfH1mx60\nusmjLhOEdVIbUx/22WNPQtG3iEYmACB5EgQkvtmfcOYFHH39r2ne5Q5Gr72X7eva+f1XH+OJPz6D\nW+Q2ntrTxCQQ6L1Wevog/OUsUlZM3esTpxLu3syyh7p65nKklSkQjCWbmRT8/a1xEl9a4+Hh/f8Z\nU2ZfZza7Td+/o2f/7GqrywRhjKkPwYCDOr2r+QeChd9yxgxv54NX3s6kL76P4dErad/yBvMf7+Sm\nz/2F5a+uLbicnmGujkAwNbn43841OaehN7YD57yHxvhjuKHdeOzPf068n1qoaNoEOZVlALhtG3tG\nRfWvtzSR9ARx4qXvI+I9VUAZlWcJwhhTURLZ1vPcCRbfITv7iGM4/fqHCR62mFHrr8PZ4XHPL17l\nz1c+QNeO/N+0ezqpRSCUcpv380OyiSkWT9/6c/J5JxOMbWfxgyv909Obexwv2Q+hHP29i+lq+h2n\nfOt7RLoKGYGVkiD6WWqj2gNdLUEYYyoq1N57AywlQQA0NoQ49/M/58AffIeGYVczZs39rH3T4Xdf\neYC3l/a/fWlqE5OEU+ZhJJuY/CGn8c70Po4jT3ovTR3PEZN92LRuY8a9WpGUjurdJ83giz/7HWNH\njyMQz78VcmqqcbJ13Eu2MwdfXSYI66Q2pn40TxzV8zwQGtiQzt1335PzrnqK1rNHMGr9Twnu6OLP\nP3yORf9amfMa9XprEIFI7+cnb70j1z8IQDC0Pe06EYFJ28EJ8OtfXk1XWlOQgqTXOJIcLaAGkdKv\n4fS3NGuVZ8rVZYKwTmpj6sfIfffpeR4IDmybT0gMCz3p/K8z84pvE3Z+xLB3lvPwTa/x6uPZh8P2\nLCnuCE5TOLUgAPYIb+KYxz/Dofse2efafc46m2Csg+a356Dhw3rfEHImCA1mX6ojTUpSyL4ya7Xn\nUCfUZYIwxtSPyXse3PPcGWANItVe+83mfT+8G4b/kvYti3j8tiW89XLfzuvk0h3iCKGW5p7jyXkQ\n06+/mVGf+QxNMw/oc+27Djqchq43s36+4+To/wjnH47rpdRGJJClDyLZf25NTMaYoWyXkaN7ngfD\n5V1obvjo8Zz2/Yfwhv+Glh1ruO9X89ixNX257t4EESAyrHcDn+Q399DYsYz+7GVZt/4MOII6WfoU\nRBF/C1PPSa8VBRoLGMWUUmsIZOmDkIyf1WIJwhhTUZFQ7zfkYGjgTUyZ2oaP4oTLbyUY/zW4Dvf8\n7NG0uQvqL/ftOEJkeMr8hALvvtqQvUYg/pBZlYa048HWxpxlOe6Ovsf6qUFUmyUIY0zFBWKJDuBQ\n+/iKlD920h5MO/9cRq29h41rG1j+yrqe97zk6rCO0DqqtzbT31afaUJZzhNF/Fynkl4rCrc29z3f\nN/nEzbROdDj5UzNTyupbvvZ5Uh11mSBsFJMx9WVTS+KGPXKX6RX7jENPuYyu3V8i0rWRv9/0dE8t\nIvnTcRzaU5q7Cv2SLjm6TSTsf/PPmOjW0E8NYkRzOx/9xhymviul2a3Eob+DoS4ThI1iMqa+RE48\nnEXTGpg6tb2in3PspVfRvPVhurrbWPdWYoJezzBXRxg5vHeNpEL3dU7ugpdOcfymM83IIE1tw3KX\nla2oLAcly7NqqMsEYYypL184eS+u/vLhWTuCy2mX3Q8iNvE1HDfKP297EujtpHYch+Gpo5j6WVk2\nTdabOgTCicSgGTWI5lEjcxfV76zp1BOTbUs2k9oYY8pm3zMuYPiml1j7ZqL/oWcUU8ChMaXDPHPz\nolxyJbXkpD910msQ7WNGZzs9UVaeZcCTtEbGMVmCMMYMKfvP+Sjx0L/xnCbeXrypZ2sGcQIEnNQJ\naqXffAUlGMk+ZHfEqMmMWft89usKrBBYE5MxxlSABELohHWgHi/e/1xKJ3XGzbbAPoisNQ2FUCSS\n9fTW0WPZd+FvezYiSrvMy77VaKZpp88BYOxx+xYWY4XUbve5McaUaMrhx7PkzjWsft2leYQ/DyJj\n34VCR7lmr2g4hBqzJ4iGxibGHLA10X2QmZMK3MriqJOO5+AjYkRayj9vpBhWgzDGDDmHHnM+oe43\nicVH4/lbjjrh9OGnoWwT1LLQkZv7HvOEcHNT9gtEaPnF02TbdjQ5oqoQ1U4OYAnCGDMEhdrH44be\nQp1G3GgiMQQa0m/ooVBhCWLLzEk8PfazacdEhUhza44roGHUlKzHi0kQtaAuE4RNlDPG5NXib//Z\nkZiX0NCYniCaQ9mbiDKdc9SFvLJb+q1S1SHSmjtBAIjf9zFt8dW917m1sUproeoyQdhEOWNMPo0T\nEqOM3GhiP4rGjBt6W3vu4aipdmkdz78v+Hf6QQ0Qac09IS5VSDwaOxIzyV0vXtA1taIuE4QxxuQz\n/l0HEu7eQiw8AYDmjBnO7aPGlF64Ck3t+WaFJ2oLbaNGIP5zq0EYY0wN2PegU4h09S7V3dqWPsPZ\niXdnXlKEAK0jRuU/DXA+elbPDnLJlWXrhSUIY8yQ1D5+Bngbel63Dk/UILYHVwEQ3nXXAZTupO0t\nkY1ookM60NpKPJBIEF6BM6lrhSUIY8zQ5ATA6d0furkl0Sl943CXn49ZRKClpfSyNUC4wFFQTiBI\n1B+xGhvWW+vY/3ilqa2A/auryCbKGWOGLA33jnQM+jf0n5x5PMs3FbBvdL8K+W6dqDUEQkFemrCE\n/d8ej06e2PPu4Wcdy+FnDTCMCrMEYYwZsgLDopCxdfQZB0woR8mFnxkMMnfik7w08T5+2HxjGT57\n8FgTkzFmyGoZ338/Qam0iAThSGIviq7QDjy1PghjjKkJE/aZmf+kEsQK6n/wRy55yjh/G9JxbYVN\nzqsV1sRkjBmydp91DDt+cgXB6NvAMQMq6/Dju3j9xYWs3zyCdVPy92E4XgwXcF2PG078Bbe9fht7\nj9p9QDEMtrpMECJyKnDq9OmV29/WGFP/2sZPZ8y7HmNUUyz/yXnsf9bJjD3uPcy+8lG+/578yWZM\n9HpCq/djxLjLGDtsEl856CsDjmGw1WWCUNV7gHtmzZp1cbVjMcbUMCfAtIuuwhm7d1mKG9/WyKLv\nvZdQ1n2q0z1y3OWsen0eR43apSyfXQ11mSCMMaZQ4QM+UN7ygoV13f7ovMNZtnF/WiPVX7a7VJYg\njDGmRFtPmU84FCZb/0ZjOMBe4wtb0K9WWYIwxpgS/eepn6t2CBVlw1yNMcZkZQnCGGNMVpYgjDHG\nZGUJwhhjTFaWIIwxxmRlCcIYY0xWliCMMcZkZQnCGGNMVqL+Ztr1SETWA2/5L9uArf08z/w5Cujd\nsDa/1DILfT/zWDVjLEd8mbGGioxvMGKs1b/n/o7Z33Nt/T3nem8o/T23q+rovBGo6pB4ADf09zzL\nz+dLLb/Q9zOPVTPGcsSXGWOx8Q1GjLX695wnVvt7rqG/51zvDcW/53yPodTEdE+e55k/B1J+oe9n\nHqtmjOWIL/V5rcZYq3/P9fQ7TH1eqzEOdnzZjg+Fv+d+1XUT00CIyPOqOqvacfSn1mOs9fjAYiyH\nWo8PLMZKGUo1iGLdUO0AClDrMdZ6fGAxlkOtxwcWY0XstDUIY4wx/duZaxDGGGP6YQnCGGNMVpYg\njDHGZGUJIgsRmSMiT4nIdSIyp9rxZCMizSIyT0TeV+1YshGRvfzf3x0i8qlqx5ONiJwhIjeKyF0i\nckK148lGRHYTkd+IyB3VjiXJ/7f3O/939+Fqx5NNLf7eUtXDvz0YgglCRG4SkXUiMj/j+Eki8rqI\nLBaRr+UpRoHtQARYWYPxAXwV+FM5YytnjKq6UFUvBc4Fyj60r0wx/lVVLwY+BpR3Z/vyxbhUVS8q\nd2yZioz1TOAO/3d3WqVjKyXGwfq9DSC+iv7bK5tiZ0jW+gM4CjgQmJ9yLAAsAXYDwsDLwN7ATODe\njMcYwPGvGwv8oQbjOw74IIl/XO+rxd+hf81pwNPAebUao3/dT4EDazzGO2ro/5uvA/v759xaybhK\njXGwfm9liK8i//bK9QgyxKjqkyIyJePwwcBiVV0KICJ/BE5X1SuB/ppoNgMNtRafiLwHaCbxP2un\niNynql4txeiXczdwt4j8Dbi1XPGVK0YREeAHwP2q+kI54ytXjIOlmFhJ1KonAi8xiK0QRca4YLDi\nSiomPhFZSAX/7ZXLkGtiymECsCLl9Ur/WFYicqaIXA/8HrimwrFBkfGp6uWq+nkSN90by5kc+lHs\n73COiPzc/z3eV+ngfEXFCHyWRG3sbBG5tJKBpSj29zhSRK4DDhCRr1c6uAy5Yr0TOEtEfkXpy0iU\nS9YYq/x7S5Xrd1iNf3tFG3I1iBwky7GcMwRV9U4S/xMMlqLi6zlB9ebyh5JTsb/Dx4HHKxVMDsXG\n+HPg55ULJ6tiY9wIVOsGkjVWVd0BXDjYweSQK8Zq/t5S5YqvGv/2iraz1CBWApNSXk8EVlcplmxq\nPT6wGMulHmJMqodYaz3GWo+vXztLgpgL7C4iU0UkTKKD9+4qx5Sq1uMDi7Fc6iHGpHqItdZjrPX4\n+lftXvJyP4DbgDVAjET2vsg/fjKwiMSIgsstPovRYqyvWGs9xlqPr5SHLdZnjDEmq52lickYY0yR\nLEEYY4zJyhKEMcaYrCxBGGOMycoShDHGmKwsQRhjjMnKEoTZKYiIKyIvpTwKWVJ9UEhiz4zd+nn/\n2yJyZcax/f0F3xCRv4vI8ErHaXY+liDMzqJTVfdPefxgoAWKyIDXMhORfYCA+qt95nAbffcM+CC9\nK+T+Hvj0QGMxJpMlCLNTE5FlInKFiLwgIv8WkT39483+BjBzReRFETndP/4xEbldRO4BHhIRR0R+\nKSKvisi9InKfiJwtIseKyF9SPud4Ecm2AOSHgbtSzjtBRJ7x47ldRFpU9XVgi4gcknLducAf/ed3\nAx8q72/GGEsQZufRmNHElPqNfIOqHgj8CviSf+xy4FFVPQh4D/BjEWn235sNXKCqx5DYXW0KiQ1/\nPuG/B/AosJeIjPZfXwj8NktchwPzAERkFPAN4Dg/nueBL/jn3Uai1oCIHApsVNU3AFR1M9AgIiNL\n+L0Yk9POsty3MZ2qun+O95Lf7OeRuOEDnACcJiLJhBEBJvvPH1bVTf7zI4DbNbEnx9si8hgk1nMW\nkd8DHxGR35JIHB/N8tnjgfX+80NJbAL1z8ReRoSBZ/z3/gg8LSJfJJEobssoZx2wC7Axx5/RmKJZ\ngjAGuv2fLr3/Twhwlt+808Nv5tmReqifcn9LYkOdLhJJJJ7lnE4SySdZ1sOq2qe5SFVXiMgy4Gjg\nLHprKkkRvyxjysaamIzJ7kHgs/62pIjIATnO+weJ3dUcERkLzEm+oaqrSaz9/w3g5hzXLwSm+8+f\nBQ4Xken+ZzaJyIyUc28D/gdYoqorkwf9GMcBy4r48xmTlyUIs7PI7IPIN4rpu0AIeEVE5vuvs/kz\niaWd5wPXA/8Ctqa8/wdgharm2iP5b/hJRVXXAx8DbhORV0gkjD1Tzr0d2IfezumkdwPP5qihGFMy\nW+7bmAHyRxpt9zuJnwMOV9W3/feuAV5U1d/kuLYReMy/xi3x868G7lbVR0r7ExiTnfVBGDNw94pI\nO4lO5e+mJId5JPorvpjrQlXtFJFvkdjIfnmJnz/fkoOpBKtBGGOMycr6IIwxxmRlCcIYY0xWliCM\nMcZkZQnCGGNMVpYgjDHGZGUJwhhjTFb/HwGMuUlD/oGOAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "rm_resonance = gd157_endf.resonances.ranges[0]\n", + "energy_range = [rm_resonance.energy_min, rm_resonance.energy_max]\n", + "energies = np.logspace(np.log10(energy_range[0]),\n", + " np.log10(energy_range[1]), 10000)\n", + "for sample in range(n_samples):\n", + " xs = gd157_endf.res_covariance.ranges[0].reconstruct(energies, rm_resonance, sample)\n", + " elastic_xs = xs[2]\n", + " plt.loglog(energies, elastic_xs)\n", + "plt.xlabel('Energy (eV)')\n", + "plt.ylabel('Cross section (b)')\n", + "\n", + " " + ] } ], "metadata": { diff --git a/openmc/data/resonance.py b/openmc/data/resonance.py index d58f706eb9..95a145da56 100644 --- a/openmc/data/resonance.py +++ b/openmc/data/resonance.py @@ -202,7 +202,7 @@ class ResonanceRange(object): return cls(target_spin, energy_min, energy_max, {0: a}, {0: ap}) - def reconstruct(self, energies): + def reconstruct(self, energies, use_sample = False, sample_parameters = None): """Evaluate cross section at specified energies. Parameters @@ -221,8 +221,8 @@ class ResonanceRange(object): raise RuntimeError("Resonance reconstruction not available.") # Pre-calculate penetrations and shifts for resonances - if not self._prepared: - self._prepare_resonances() + if not self._prepared or use_sample: + self._prepare_resonances(use_sample, sample_parameters) if isinstance(energies, Iterable): elastic = np.zeros_like(energies) @@ -394,8 +394,11 @@ class MultiLevelBreitWigner(ResonanceRange): return mlbw - def _prepare_resonances(self): - df = self.parameters.copy() + def _prepare_resonances(self, use_sample = False, sample_parameters = None): + if use_sample == False: + df = self.parameters.copy() + else: + df = sample_parameters.copy() # Penetration and shift factors p = np.zeros(len(df)) @@ -653,8 +656,11 @@ class ReichMoore(ResonanceRange): return rm - def _prepare_resonances(self): - df = self.parameters.copy() + def _prepare_resonances(self, use_sample = False, sample_parameters = None): + if use_sample == False: + df = self.parameters.copy() + else: + df = sample_parameters.copy() # Penetration and shift factors p = np.zeros(len(df)) diff --git a/openmc/data/resonance_covariance.py b/openmc/data/resonance_covariance.py index e761a3756e..6ed3db424e 100644 --- a/openmc/data/resonance_covariance.py +++ b/openmc/data/resonance_covariance.py @@ -113,7 +113,7 @@ class ResonanceCovariances(Resonances): # Throw error for unsupported formalisms if formalism in [0,7]: - raise TypeError('LRF= ', formalism, + raise NotImplementedError('LRF= ', formalism, 'covariance not supported for this formalism') if unresolved_flag in (0,1): @@ -122,14 +122,14 @@ class ResonanceCovariances(Resonances): erange = _FORMALISMS[formalism].from_endf(ev, file_obj, items, file2params) elif unresolved_flag == 2: - warnings.warn('Unresolved resonance not supported.' - 'Covariance values for the' - 'unresolved region not imported.') + warn_str = 'Unresolved resonance not supported.'\ + 'Covariance values for the unresolved region not imported.' + warnings.warn(warn_str) + ranges.append(erange) return cls(ranges) - class ResonanceCovarianceRange(object): """Resonace covariance range. Base class for different formalisms. Parameters @@ -243,6 +243,7 @@ class ResonanceCovarianceRange(object): param_list = ['energy','neutronWidth','captureWidth'] mean_array = pd.DataFrame.as_matrix(parameters[param_list]) spin = pd.DataFrame.as_matrix(parameters['J']) + l_value = pd.DataFrame.as_matrix(parameters['L']) gf = pd.DataFrame.as_matrix(parameters['fissionWidth']) mean = mean_array.flatten() for i in range(n_samples): @@ -253,9 +254,9 @@ class ResonanceCovarianceRange(object): gt = gn + gg + gf records = [] for j, E in enumerate(energy): - records.append([energy[j], spin[j], gt[j], gn[j], + records.append([energy[j], l_value[j], spin[j], gt[j], gn[j], gg[j], gf[j]]) - columns = ['energy', 'J', 'totalWidth', 'neutronWidth', + columns = ['energy', 'L', 'J', 'totalWidth', 'neutronWidth', 'captureWidth', 'fissionWidth'] sample_params = pd.DataFrame.from_records(records, columns=columns) samples.append(sample_params) @@ -264,6 +265,7 @@ class ResonanceCovarianceRange(object): param_list = ['energy','neutronWidth','captureWidth','fissionWidth'] mean_array = pd.DataFrame.as_matrix(parameters[param_list]) spin = pd.DataFrame.as_matrix(parameters['J']) + l_value = pd.DataFrame.as_matrix(parameters['L']) mean = mean_array.flatten() for i in range(n_samples): sample = np.random.multivariate_normal(mean,cov) @@ -274,9 +276,9 @@ class ResonanceCovarianceRange(object): gt = gn + gg + gf records = [] for j, E in enumerate(energy): - records.append([energy[j], spin[j], gt[j], gn[j], + records.append([energy[j], l_value[j], spin[j], gt[j], gn[j], gg[j], gf[j]]) - columns = ['energy', 'J', 'totalWidth', 'neutronWidth', + columns = ['energy', 'L', 'J', 'totalWidth', 'neutronWidth', 'captureWidth', 'fissionWidth'] sample_params = pd.DataFrame.from_records(records, columns=columns) samples.append(sample_params) @@ -286,6 +288,7 @@ class ResonanceCovarianceRange(object): 'fissionWidth', 'competitiveWidth'] mean_array = pd.DataFrame.as_matrix(parameters[param_list]) spin = pd.DataFrame.as_matrix(parameters['J']) + l_value = pd.DataFrame.as_matrix(parameters['L']) mean = mean_array.flatten() for i in range(n_samples): sample = np.random.multivariate_normal(mean,cov) @@ -297,9 +300,9 @@ class ResonanceCovarianceRange(object): gt = gn + gg + gf records = [] for j, E in enumerate(energy): - records.append([energy[j], spin[j], gt[j], gn[j], + records.append([energy[j], l_value[j], spin[j], gt[j], gn[j], gg[j], gf[j], gx[j]]) - columns = ['energy', 'J', 'totalWidth', 'neutronWidth', + columns = ['energy', 'L', 'J', 'totalWidth', 'neutronWidth', 'captureWidth', 'fissionWidth', 'competitveWidth'] sample_params = pd.DataFrame.from_records(records, columns=columns) samples.append(sample_params) @@ -310,6 +313,7 @@ class ResonanceCovarianceRange(object): param_list = ['energy','neutronWidth','captureWidth'] mean_array = pd.DataFrame.as_matrix(parameters[param_list]) spin = pd.DataFrame.as_matrix(parameters['J']) + l_value = pd.DataFrame.as_matrix(parameters['L']) gfa = pd.DataFrame.as_matrix(parameters['fissionWidthA']) gfb = pd.DataFrame.as_matrix(parameters['fissionWidthB']) mean = mean_array.flatten() @@ -320,17 +324,19 @@ class ResonanceCovarianceRange(object): gg = sample[2::3] records = [] for j, E in enumerate(energy): - records.append([energy[j], spin[j], gn[j], + records.append([energy[j], l_value[j], spin[j], gn[j], gg[j], gfa[j], gfb[j]]) - columns = ['energy', 'J', 'neutronWidth', + columns = ['energy', 'L', 'J', 'neutronWidth', 'captureWidth', 'fissionWidthA','fissionWidthB'] sample_params = pd.DataFrame.from_records(records, columns=columns) samples.append(sample_params) elif mpar == 5: - param_list = ['energy','neutronWidth','captureWidth','fissionWidthA','fissionWidthB'] + param_list = ['energy','neutronWidth','captureWidth', + 'fissionWidthA','fissionWidthB'] mean_array = pd.DataFrame.as_matrix(parameters[param_list]) spin = pd.DataFrame.as_matrix(parameters['J']) + l_value = pd.DataFrame.as_matrix(parameters['L']) mean = mean_array.flatten() for i in range(n_samples): sample = np.random.multivariate_normal(mean,cov) @@ -341,15 +347,43 @@ class ResonanceCovarianceRange(object): gfb = sample[4::5] records = [] for j, E in enumerate(energy): - records.append([energy[j], spin[j], gn[j], + records.append([energy[j], l_value[j], spin[j], gn[j], gg[j], gfa[j], gfb[j]]) - columns = ['energy', 'J', 'neutronWidth', + columns = ['energy', 'L', 'J', 'neutronWidth', 'captureWidth', 'fissionWidthA','fissionWidthB'] sample_params = pd.DataFrame.from_records(records, columns=columns) samples.append(sample_params) self.samples = samples + def reconstruct(self, energies, resonances, sampleN): + """Evaluate the cross section at specified energies for an already + sampled set of resonance parameters. + + Parameters + ---------- + energies : float or Iterable of float + Energies at which the cross section should be evaluated + resonances : openmc.data.Resonance object + Corresponding resonance range with File 2 data. Used for + reconstruction method + sampleN : int + Index of sample of resonance parameters to be used + + Returns + ------- + 3-tuple of float or numpy.ndarray + Elastic, capture, and fission cross sections at the specified + energies + + """ + if self.samples[sampleN] is None: + raise ValueError("Sample of resonance parameters has not been set.") + sample_parameters = self.samples[sampleN] + xs_array = resonances.reconstruct(energies, use_sample = True, + sample_parameters = sample_parameters) + return xs_array + class MultiLevelBreitWignerCovariance(ResonanceCovarianceRange): """Multi-level Breit-Wigner resolved resonance formalism covariance data. Parameters