diff --git a/docs/source/examples/index.rst b/docs/source/examples/index.rst index 89a1f1fe0b..bcb1b1ad8f 100644 --- a/docs/source/examples/index.rst +++ b/docs/source/examples/index.rst @@ -24,6 +24,7 @@ Basic Usage triso candu nuclear-data + nuclear-data-resonance-covariance ------------------------------------ Multi-Group Cross Section Generation diff --git a/docs/source/examples/nuclear-data-resonance-covariance.rst b/docs/source/examples/nuclear-data-resonance-covariance.rst new file mode 100644 index 0000000000..4b505c9a58 --- /dev/null +++ b/docs/source/examples/nuclear-data-resonance-covariance.rst @@ -0,0 +1,13 @@ +.. _notebook_nuclear_data_resonance_covariance: + +================================== +Nuclear Data: Resonance Covariance +================================== + +.. only:: html + + .. notebook:: ../../../examples/jupyter/nuclear-data-resonance-covariance.ipynb + +.. only:: latex + + IPython notebooks must be viewed in the online HTML documentation. diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 7feaa8d608..e7af5e273e 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -79,6 +79,11 @@ Resonance Data openmc.data.MultiLevelBreitWigner openmc.data.ReichMoore openmc.data.RMatrixLimited + openmc.data.ResonanceCovariances + openmc.data.ResonanceCovarianceRange + openmc.data.SingleLevelBreitWignerCovariance + openmc.data.MultiLevelBreitWignerCovariance + openmc.data.ReichMooreCovariance openmc.data.ParticlePair openmc.data.SpinGroup openmc.data.Unresolved diff --git a/examples/jupyter/nuclear-data-resonance-covariance.ipynb b/examples/jupyter/nuclear-data-resonance-covariance.ipynb new file mode 100644 index 0000000000..2c44b6f7fa --- /dev/null +++ b/examples/jupyter/nuclear-data-resonance-covariance.ipynb @@ -0,0 +1,964 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this notebook we will explore features of the Python API that allow us to import and manipulate resonance covariance data. A full description of the ENDF-VI and ENDF-VII formats can be found in the [ENDF102 manual](https://www.oecd-nea.org/dbdata/data/manual-endf/endf102.pdf)." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import os\n", + "from pprint import pprint\n", + "import shutil\n", + "import subprocess\n", + "import urllib.request\n", + "\n", + "import h5py\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import openmc.data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### ENDF: Resonance Covariance Data\n", + "\n", + "Let's download the ENDF/B-VII.1 evaluation for $^{157}$Gd and load it in:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Download ENDF file\n", + "url = 'https://t2.lanl.gov/nis/data/data/ENDFB-VII.1-neutron/Gd/157'\n", + "filename, headers = urllib.request.urlretrieve(url, 'gd157.endf')\n", + "\n", + "# Load into memory\n", + "gd157_endf = openmc.data.IncidentNeutron.from_endf(filename, covariance=True)\n", + "gd157_endf" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can access the parameters contained within File 32 in a similar manner to the File 2 parameters from before. " + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
energyJneutronWidthcaptureWidthfissionWidthAfissionWidthBL
00.03142.00.0004740.10720.00.00
12.82502.00.0003450.09700.00.00
216.24001.00.0004000.09100.00.00
316.77002.00.0128000.08050.00.00
420.56002.00.0113600.08800.00.00
\n", + "
" + ], + "text/plain": [ + " energy J neutronWidth captureWidth fissionWidthA fissionWidthB L\n", + "0 0.0314 2.0 0.000474 0.1072 0.0 0.0 0\n", + "1 2.8250 2.0 0.000345 0.0970 0.0 0.0 0\n", + "2 16.2400 1.0 0.000400 0.0910 0.0 0.0 0\n", + "3 16.7700 2.0 0.012800 0.0805 0.0 0.0 0\n", + "4 20.5600 2.0 0.011360 0.0880 0.0 0.0 0" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "gd157_endf.resonance_covariance.ranges[0].parameters[:5]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The newly created object will contain multiple resonance regions within `gd157_endf.resonance_covariance.ranges`. We can access the full covariance matrix from File 32 for a given range by:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "covariance = gd157_endf.resonance_covariance.ranges[0].covariance" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This covariance matrix currently only stores the upper triangular portion as covariance matrices are symmetric. Plotting the covariance matrix:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUIAAAD8CAYAAAACGq0tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXvwZVV15z9fugWCTuTRQBBwIKGdDDqOml/QTDIJEXk5iW0mWkEzCTFYJCnIJDEphbFGDWoNZpKQODFW9WhHzEMkThw7Y4+kRYkzKUUaNGCrpDtotLWHh01ICIXQsuaPcy4cNnvvs8/j3t99rE/Vrd895+yz99rnd++667EfMjMcx3FWmUPWWwDHcZz1xhWh4zgrjytCx3FWHleEjuOsPK4IHcdZeVwROo6z8rgidBxnXZB0nqTbJe2VdFnk+mGS3l9fv1HSKfX5J0m6WtJtkr4g6fKhskxNEbZ10nGc1UXSBuAdwPnA6cArJJ0eFLsIuNfMTgOuAt5Wn385cJiZ/Svge4CfmyjJvkxFERZ20nGc1eUMYK+Z3WFmDwHXAFuCMluAq+v3HwDOkiTAgCdL2gh8G/AQ8A9DhNk45OYMj3YSQNKkk5+PFd60aZOdcsopUxLFcRyAm2+++R4zO7bv/adJ9kBh2f2wG3iwcWqrmW1tHJ8IfLVxvA94flDNo2XM7KCk+4BjqJTilqoZjgB+xcwOFHckwrQUYWsnJV0MXAzw9Kc/nV2f/jSPBAbqITwC8Oj5Q3jkce8nPMIhTyi7XjRlnDcZcuch/+zayoT/m2k8g7De2Odh8lno0n6ufNu1SZu5suFndWwZSp+JNmz4u2glhTwA/Fxh2TfBg2a2limiyLlwvm+qzBnAt4CnAUcB/0fSRyeGVx+mpQhbO1n/OmwFWFtbM3j8BwZI/kOb1ybXYwpyPWjKOG8y5M5D/tm1lQn/N9N4BmG9sc/D5LPQpf1c+bZrkzZzZcPP6tgylD6ToYhRY2n7gJMbxycBX0+U2Ve7wU8FDgCvBD5iZg8Dd0n6K2AN6K0Ip2W6lHTScZwFQlSWU8mrgJuAzZJOlXQocAGwPSizHbiwfv8y4GNWrRLzFeCFqngy8ALgi707xvQUYUkns8QskLZftpQL8giHPPoKr6XazdXjrC6r/v8/pPDVhpkdBC4FrgO+AFxrZrslXSHpJXWxdwPHSNoLvAaYjD55B/AU4HNUuuYPzOzWIf2aimtcBzYnndwAbDOz3V3qmCi9lMsVugqhS9Q8bisbXksdr7fL66w/q/4ZGPNnwMx2ADuCc29ovH+QaqhMeN/9sfNDmFaMMNrJruSUUixGmIqVDCEVc1kU5iFm6SwHI8cI54qpKcJpklOCoSs7KdM81yXDnMrClZRfL0IZSrOaJXKvZ9Y4bH9RssYlZfvK0KwzPJ5G8soV4RzR9QMXux67lmuv5FzXepuMaW2mwglds8YxpdIna5yrp4slH7tWmjXO9b1NUZRkjZv1xBRRKuYdk2/IM1mgrPFcsZCK0HGc9WHDegswJRZLET74IBx+eLZIzE0Y2zXO3Tv013eMX+4xXeNcFj3VZqydtmx8SRvhta6ucV/3t6+7m5OhzTUudZtj7XQJ5XRBuCKcD2olGH5gYu5Q25c+dTxEcebcv1zZmKwlrmTJl6XE/SqVeREo+QzMI4si52JI2Z3FUoQ1KcXR/NumMFIKMFRQzfJ9Y4qp+9rilW3ncvGyNlm7Xm+LYZbEImNlh15LxYTDvznGiMHl6knFWIfERrvK5zHCPAupCKH9i5pTYCkLMlV3qQKcR5r9HKIcx3wGY7WT+kGL1ZNTJl3a7PMcSu7p+0z6/OgOwRWh4zgrzWSK3TKyFAo+dGvDoQx9g/HhcayurnJ2aW9sFiUONQ9M+1kt6v+ibWrd5LVoLKLMT6A0s5eK1cTqSo37GuJilMQBx3bDu7qauWtjfXm7/BDlSD3P0nhrrs1cTK4rzR/nVP19n0npD+lYMcLS16KxFIoQyuJfYfwwJJYkCa/FjkuVROwLkaq3y/CPVD3hK2UhT8hdb3smbc+/TQnkAvttX+7U/yIlayrRVtJm3xhhqKDbZEj9GIf3x364uyj2rmwofC0ay+ryA+lMcCyz1rwWO05da37A2zKnsexert6mfM2yKYsn98UJ64opppKMY9jvlOUcttv2bEvajPUj1d+wzdT/JfW5iD3r1OcmPE71K1Zvl2fSpMvzi8nbB88aLxC5L0WzTKjEUvfnvuCpD2eJC11abxfXOaeEwzIpZZCSIWXRlfwIpBRW137nlFpOobW12aZsYnWVKPYUqefXJl+sntIf7rFwReg4zkqzzFnjZe3Xo3SNtcXKtv1qzzM567ALi9TneSPlquYs4C51546H1v/EepaT3v2SdLKkj9cbLO+W9Ev1+aMl7ZS0p/571HjijkdTuYVJjFSCIXR5cx/usM6QWIA/F/xPJQRix7kvQ6qPubhSiXxtMcJUPal+9GkzjO+l2ojJlrqW+n+31ZOSPay35Jnk2sw9v5R8fZnECEtei8YQmQ8Cv2pm/5Jqz4BL6r2LLwOuN7PNwPU8trz2upBSCk3rrvkKzzWP2+oOYzS5WGHsS5yKRYX1tB2nlGIY72oqwJjlWBJniynR8EtX+lzb/gdtbTaVQK6NXF/Ca6lnWfJMuvzvc8clbbaVHQsfPhNgZvvN7Jb6/T9S7TtwIo/flPlq4KVDhRxKiRJLkbO8VoncF9vJE3tuY32O2pRj2/mujDl8RtJ5km6XtFfSEwwmSYdJen99/UZJpzSuPVvSJ2tv9DZJ+WWpWhjlv1EL+FzgRuB4M9sPlbIEjkvcc7GkXZJ23X333WOIUUzuQ9Hm0sbKz5Jce2Mq6nlU+POgiPs+l7YfknnoWxtj7mInaQPVJkznA6cDr6g9yiYXAfea2WnAVcDb6ns3An8E/LyZPRM4E3h4SN8Gf9olPQX4H8Avm9k/lN5nZlvNbM3M1o499tihYjiOM2VGjhGeAew1szvM7CHgGipvsknTu/wAcJYkAecAt5rZXwOY2TfM7Fu9O1YucxxJT6JSgn9sZn9Wn75T0gn19ROAu4a0MQ1yv+olsamw/CxpsyrGkieW8FhvK3G924d+/++YlxF7vn0oTYqM54qPpghPBL7aON5Xn4uWqbf/vA84BngGYJKuk3SLpNf2681jDMkai2rf0S+Y2W83LjU3Zb4Q+FB/8aZDGARvUpINXEU8RtifZYoRdlCEmyahr/p1cVBVLKdihWU2Aj8A/GT998ckndWnPxOGjCP8fuCngNskfbY+95+AK4FrJV1EtSP9qPuPjsVEGTaHMTTfN2lmgsPzsaxs7EOeuze0FNrqSX2JUpnj2L2hpRsOuwjliD2nVPmwnbA/ubpTfUz1LexPrB9NuWPPNpQjrDtWNlYu15+wfEqGsL+pz0Lsf5mTYwxF2HGK3T1mtpa5vg84uXF8EvD1RJl9dVzwqcCB+vxfmtk9AJJ2AM+jGqXSi96K0Mz+L+lM+SDtPCtyimFC7kOUGm5R8qErGRYx5Nc9ZYW0WScxuVLPqXmc++Klwgypcyn529psytr2TEv+x211pZRtqlxKlj4yNK/FjlPvhzLi0JibgM2STgW+BlwAvDIoM/EuPwm8DPiYmZmk64DXSjoCeAj4IapkSm+WfmZJGznFkPpVdhfRGcpYVtosEfCkkeoys4OSLgWuoxpxs83Mdku6AthlZtupQm9/KGkvlSV4QX3vvZJ+m0qZGrDDzD48RJ6VV4SO45QzZoTczHYAO4Jzb2i8f5BEaM3M/ohqCM0orGbkP0HXX+n1SJxMcxxhLPY2DaaVKGg7PwarmiyD5Z5i5xZhg9KgdPParCmNnQ2te5p9G6vuXMxz3t3OeZcvxSIquRJcEQaUZHdLlNGysaz96kqf5zBmsiJW96zwhVlXlNgQj9J7pi1T6pqzGEzrczLtcMGyfsJcETqOU8QyL8y6rAp+FNpcmjYLbBozUdrc8iG//NNMxHRpp5RFSZaUDLDvywJPsZsrllXBj0bXeOCsEg6lNF2wNnmGxD679HWsGGtuNkgu1jukzb4xwiFlulyb5o/DMscIl7VfoxObPhazkqYdp+tqtYWzF9ZDrlnSRxHMQvbwf9PlczMvzxaW1yJcRJnXhYkCTE2jijEN1ziHu8b582PgrvFyKkJ3jTuQmlMaKxO+H1uGUtw1Tpft0+aqu8bLqjCWtV+O40yBRdyPpIRFtGLXnZSF1WWKWmm8aGz3elpu2rzEsfr0ZZlihNPsixh3z5J5wi3CAeTihV1cyZx7M1RJhcmS3Go7Oaa18k4oU5eBxtN2aWP0GQg97RhhSV1jxgiXkcH9krRB0mck/a/6+NR6x6k99Q5Uhw4Xcz4JxxnOKqHQt54hX76hiZhcvbnj3LX1sEL7xgjD55fL7g+RZxbDZ5YxWTKGzL9EtZXnhLcBV9X7Gt9LtRPV0hJab83scpOYaxQOw4kprS4uVdfhM7EseIm8MfoOTxniInZNrJS0MQ1y/9+hMnnWeBwGySzpJODfAe+qjwW8kGrHKZiTfY1nRd+hH0NcvJgi7lJXaNX2rauvkhxC3y/3NIfXLDNjbuc5bwxV3r8DvBYe/WQdA/x9veMUxHemAtZ3X2PHcbrjrnEEST8C3GVmNzdPR4qGO1NVJ5dwX+OuY+1KrLm22FKTaQ7gnscY4ZB2ZklJjHBI3bnjtvPd21tORTh0F7uXSHoxcDjw7VQW4pGSNtZWYWxnqpVjEocLY0TNTGwqPhieS9XdZFJfW9wuN+g49UVNyd32vq3/zTZzZVP9DH9UYm3miMUf2+QO7+3SXq6eyfXmuZR8bW3E7h3CIiq5Enr3y8wuN7OTzOwUqk1VPmZmPwl8nGrHKZjTfY2nTZjoaI47bH5pY+9jVmIunphK1nS1TlN1xSzBmNxNhRQqp5L+tz2P2PlQcTeVRqx86pXqb9iXnOKNlS19dm2fgZx8uXrCz99QxnaNJZ0n6XZJeyVdFrl+WD3yZG89EuWU4PrTJd0v6dd6d6pmGgr+dcBr6p2njqHaiWrlCD+cuV//Lr/Y03J924h9qVeRPs8gpuBiFt4Y8qTqGet/N5YilLQBeAdwPnA68ApJpwfFLgLuNbPTqLbrfFtw/Srgf/fryeMZJcFjZjcAN9Tv7wDOGKPeZSL2QYy5ol0/sH3rKHGNY1/cMb5Qfa2U0i992FZp2RL6yt4mx1iW25gWYMjIc43PAPbW+gJJ1wBbgM83ymwB3lS//wDwe5JU7238UuAO4J/GEGZZXX7HcaaApKIXsGkyKqR+XRxUdSLw1cZxbITJo2XqnMN9wDGSnkzlef76WP1axCE/C0n4S52L++VoixHG2iqpZ4hMXelbd8xK7dPPWZNK9IxVdxfLeBASbCxUGQ8/fI+ZreVqi5wLR5ikyvw61aSN+2ulOxhXhDOiqxuacndzH/xpxg/nQaHMA31d+mn9b2b+fylXhG0l9gEnN45jI0wmZfZJ2gg8FTgAPB94maTfAI4EHpH0oJn9XplwT8QV4YyJDSsJz+fIKdRYYmZIzCg1fCQs0yZ3WE+p0m+jS7+aMpTeN1YSq6sS7CrjzJRhF4uwnZuAzZJOBb5GNfLklUGZ7VQjTz5JNRLlY2ZmwL99TCS9Cbh/iBIEV4QzJ6as2pRNSX3TJGfRlFiqoaLu006MPkqg67MtlWdMq29ure9DDoHDDy8r+4//mL1sZgclXQpcR7Vy1zYz2y3pCmCXmW2nGnHyh/UIlANUynIquCJ0HKeMcS1CzGwHsCM494bG+weBl7fU8aYxZHFFuE6UBv1zrnDI0KB8Lm45LboMiRkiz3oOn4lZi2M93z5DigYxoiKcJ5azVwtELH7WJW7W9QsVKx+LnZUq6j5xtJJ4YtdrJW1PQ/Gsdz19Mum9GdkinCeWs1cLRMxaGDrLIKechliIYRuls03asuJjzaooZRoKNTZTqOT/MbZ8U8UVoTNNcpnYLlnkti9cV6sul+To4rJ3qaeLa9xFYYw1RjLXx5TFm1KcKZmGJMtKQyy9cEXozJJUZjlG1+EzKStvveYwd2EeXOO+jD3VL9fOVF3j0qzxguGK0HGcMtwidGZJlyxyGznXuc3VzA2ozg2SDimpJyVDSvZcm22yhO3niD2/VBa4Sxu5Zz+NAdUzn2K3YCxnrxacpqvbNuxizBhhKt6Ui+2VfHnHiBGWxNly51MytIUFYqGH2H1DY4Rt5VN4jHAclrNXS0LqCxqzTEqC7W1WR0n8MGyzTVGX1tOnzVyyIWWxpuRLXRtiEeb6krIWc9ZkTp4u8vVmiRXhoKcj6UhJH5D0RUlfkPR9ko6WtLPe13inpKPGEnbVyH3om2WaFktp4qPNOgvrbWuzeS384odlc/fmyoblQiUxOQ6tuNyPRNu1VMKqrZ6w3bayueOSNtvKjsbGjWWvBWPoz8TvAh8xs+8G/jXV/saXAdfX+xpfXx87PQm/VDEl11QAOSWYUh4lX5rcF3tWrHf7fWWIKca2eG1feabqGk/mGpe8FozeqlvStwM/CPwMgJk9BDwkaQtwZl3saqqVq183REjHceYAd42jfCdwN/AHkj4j6V31yrHHm9l+gPrvcbGbfV/jfpRYiF1ihE0rMRXneqTRaqrNFDkLtgt9s8ZtVvIYcpTc0yZHXxlLY4GjxgjdNX4cG4HnAe80s+dS7R1Q7AYv477GsyD3pWpzX7skV9rqCeUpKduVZt2lsc9Y+/PgUi8FS6wIh0i8D9hnZjfWxx+gUoR3SjrBzPZLOgG4a6iQzmPklFYqsD8hN9wkdU/ufNu1ofe0lUtlvmP3Tmu4St97SmN7Xa9NNUYIC6nkSuhtEZrZ/wO+Kulf1KfOotqBarKqLKzovsazIhwiMqHkQ9/mspa6xm1uc9+yzfKT9yVf+pi73zZcJefyp6zfrv2MHZc+ky4u9VRdY0+WJPlF4I8lHUq1td6rqJTrtZIuAr5Cy8KKjuMsCEucLBnUKzP7LBDbqeqsIfU6ZaTGuaXcwEn5krhZGENsG0+XiuW1lU21mZIhtNDaQgAlbeRCA2FduX62hRhS/Zq00ybfurPEinBOnrDTlzCD3OZGdUmChPcNjTN1Gbc4r0wzJroQjJgskXSepNsl7ZX0hESrpMMkvb++fqOkU+rzZ0u6WdJt9d8XDu7W0Aqc+aHLF66PJdl8P0QhDL0/dTwLxlLkY9UzVuKpiBEtQkkbgHcAZ1MlXm+StN3MPt8odhFwr5mdJukC4G3ATwD3AD9qZl+X9CyqDaDCzeE74RbhEhFahLEkSlOhNRMLuYRLm5WZazOWRCgpGytfSs6dLE0uNOuKlR3SzzEprXsUGcYdPnMGsNfM7qgnY1wDbAnKbKGalAHVqJSzJMnMPmNmkz2QdwOHSzpsSNfcIlwiQosrNbQkfJ+yFlL1xGJczeO2GGHf45IYYdvwma5WZSre2CVGWJLZ7zJEpk89o1iEXbbzhE2SdjWOt5rZ1sbxicBXG8f7qDZuJ1am3v7zPuAYKotwwo8DnzGzb5YKFsMVoeM45ZS7xveYWSyROkGRc9aljKRnUrnL55QKlcIV4RJSajXFCF3nyf0p6yfV/jTcwRIrLCZvm0xDrk3abCs7YZYZ4FC+wYybNd4HnNw4Pgn4eqLMPkkbgadSbfSOpJOADwI/bWZ/O1QYV4TOo8SUZuo457a2JQL6Jgpy8pUkc9pCACXXwn6GMqSG4IwpQx/5Ro0RjsNNwGZJpwJfAy4AXhmUmUzO+CTwMuBjZmaSjgQ+DFxuZn81hjCeLFlmDh6cehNdEholw3sm1x9L5aRnlqRofvGHJC5iSiZVT9d6x0hyDJGhFyMmS8zsIHApVcb3C8C1ZrZb0hWSXlIXezdwjKS9wGt4bC2DS4HTgP8s6bP1K7q4SyluES4zGzc+IcAP5a7x5LjN4gjvjV1LWUp9LcNJnakERe6+mHy5etrK9rW8YmVL5eiS5BnVNR5x+pyZ7QB2BOfe0Hj/IJGZaWb2FuAtowmCK8KVoo/iacsYx76cpW5yzn0tzT6n3NQ2edvuC2WZ0IxBhgowjE+mFFHMfS/JOJdk6FM/EnPoGs8Vy9krx3HGxxWhswx0zfbGxtCVWHu5bG5Yd86iSWWCcy5fKmmSOp/LNseOY5n00ErMZe0nxzE5UvLG2gyvxcqmZOiNK0JnWXjcF+LBBwF45PAjOHiw+oynFF7z/hKXqy1+1uYe5sqWxAVzyjPXZpsSy8nXdq20bPgjUPK8SssOwhWhs5Q0At9dP99tlmBb2bZ6S5Rsk9KkSe56KgYXk6lN4eXa6hKrHCtGOAquCJ2l4+DBx4bXHH4E998PT3nK44t0cam6Jjm6Hndx+WKub07uLtZSyfNoU1ipMn0orWc013gBF10tYZC9LOlXJO2W9DlJ75N0uKRT6yVz9tRL6Bw6lrDOiGzcCF/8YvUC9u17YpE21zN2HFNWzWux49KyTespp6TD+F3YRqofJQoxVSbVj9Q9OTm6UFrPqK7xEu5Z0vvpSDoR+I/Ampk9C9hANTr8bcBV9b7G91ItpeM4zqKzxIpwqMQbgW+T9DBwBLAfeCGPTZW5GngT8M6B7TjT4DnPASoL6vTvrk61xZjassYp1zjn1pYkNHJxtSalWeOUfLmsbFhvKrGUi/uVZI3bZEhdix3HnmdvljhG2NsiNLOvAb9JtS/JfuA+4Gbg7+vpM1BNmh60YKIzI3btai/TQp/gfKkrl3M9J23GsqQlbmpJm0PKxu7t47qPKUNv3CJ8PJKOolo48VTg74E/Bc6PFA2X1pncfzFwMcDTn/70vmI4Y7FWrZgUS260WUSx8qWJkFT2uaTenMUX1tFmFaXK5/rSZnnl+lbyjPpea8tE92aJLcIhvXoR8CUzuxtA0p8B/wY4UtLG2iqMLa0DVBu8A1sB1tbWosrSWT9K3alUhjhWT1vZVN0ppdmWLGlTSjEZ28r3cTNL+p2rt8u1XCZ9MN0WZl0ohjydrwAvkHSEJPHYvsYfp1oyB3xfY8dZHjxZ8kTM7EZJHwBuAQ4Cn6Gy8D4MXCPpLfW5d48hqDNbYsH85rXY+1Q9pWVz90KZ9ZNyB3NJni5y9XEzu7j/XdvsYiWPwgIquRKG7mv8RuCNwek7qDZmcRacWCYzlY2dXIuRyobGysXaaKs/JnPsfKqeqbqTQV0l8c7YvUNd41HwGKGzqrQlH4ZaU7EhJCVlU2VStA37Sck7hjJsU8yzYNQB1UvIcvbKGZWYguqSlR3DtSttM/WFzynQEqU7hKFKfW5c4yVOlrgidIoIramYuxuWjZVvkhso3WU4S47YeMPSAdW5OpvWcBeZcpZpbjB2n7JTsTqX1CKcwQhMx3GWgpGzxpLOk3S7pL2SLotcP6xer2BvvX7BKY1rl9fnb5d07tCuLad6d6ZGiVXWlugILZYS9650bGBJUibnYrdZUTlLLHdP2Habm9+WAOlqKc5bjFDSBuAdwNlUM9BukrTdzD7fKHYRcK+ZnSZpso7BT0g6nWpdg2cCTwM+KukZZvatvvK4Reh0pu2LlftCjxl76zqtbQzGGj7TpczcxAjHtQjPAPaa2R1m9hBwDdVMtSZbqNYrAPgAcFY9ZnkLcI2ZfdPMvgTsZeBIFbcInV7ELJuUJdY2PCVlwcw6sxoyzbm7KzB8ZpOk5gT2rfVssgknAl9tHO8Dnh/U8WgZMzso6T7gmPr8p4J7B61p4IrQGURTeQ39AuaGz4w5BMWHz/TDDB46WFzPPWa2lrmuWBOFZUru7YQrQmcQJRZhlzGIXVy7LtZP3/GHYzD2mMhc3dN0jc0eW9R8BPYBJzeOY+sSTMrsk7QReCpwoPDeTniM0HGcIiaKsORVwE3A5npF+0Opkh/bgzLbqdYrgGr9go+ZmdXnL6izyqcCm4FPD+mbW4TOaDRjfbHxe2EsMDfGr23MYbN8aF22jXOMyZ0biN02SLskOxuWDd/3lS9XNtbmEMa0COuY36XAdVSr228zs92SrgB2mdl2qnUK/lDSXipL8IL63t2SrqVa5OUgcMmQjDG4InRGJvclbJudkjo3Od9lqEjpcJiwnq7Xmm3kkkJtQ22mkSwZffgMo7rGmNkOYEdw7g2N9w8CL0/c+1bgrWPJ4orQGZVDDj4EGzcWffFiymFaWc/1zELHFNFY8vSJkfZl5BjhXOGK0BmXenhF6cDkmLWUU6Ix1zi8N2UB9XV/x3CNQxlSyrGrDG1t5OrtyiOPwIMPjlLV3OGK0JkqMcsn57rmFGgf17ivjF3bzNVTkk1fBNfYLULHcRyWVxG2/kxI2ibpLkmfa5w7WtLOehP3nfVGTqji7fVk6FslPW+awjvzzyEHH3riucCCaZKyXLoMSg7rD7PYQ6yj9Z7tsp6MPHxmrij5RLwHOC84dxlwfb2J+/X1MVS72G2uXxfj+xk7dcwwpohySjA3/CZVJnb+kIYqbJ4rYZrzk2P9G6vuIeVyrLQiNLNPUI3hadKcDH018NLG+fdaxaeodrQ7YSxhncWlRBG1zUYZojz6jt8bgxJrtq/FWlp2lKmCdbKk5LVo9I0RHm9m+wHMbL+k4+rzsYnUJ1JtAP84fF/j1SOWVBgyRWyaU+xK6h5ijeWm2HVRwLOcYgeLae2VMPYUu+LJ0Ga21czWzGzt2GOPHVkMx3HGZpld474W4Z2STqitwROAu+rzo0+GdpaPlDWYGwITs9BKkyVdyreRGz4zJs26+1rDY8u1zMNn+lqEzcnQzU3ctwM/XWePXwDcN3GhHSckFzuLDbbuEw+LJUu6yDEGscx1W+KnKXNOvraEUyox1YeVtgglvQ84k2qhxX1U+xhfCVwr6SLgKzw2H3AH8GKqFWMfAF41BZmdJaHN8ssNs2mrN3bfLKe0hfe0yT9W1niaMcJltghbFaGZvSJx6axIWQMuGSqUszqE0+RiFkyoGNuUW3Oa3RAFMNbwGUivptPn3hzTnFNttpgZ4RJ8Zomz7nRVcH0USBjbK6ljrKlpOaUeky92f59rKRn6stIWoeM4DrgidJyZ0BbbK7WiYvX0Hac3Bl3HEQ6p22OE/XBF6MwduaW0Jtfb7p8wJDbX1Z2cloLtsvJMs7ymEGShAAARBUlEQVQPnynHFaEzd6QSHV1ihLl7pxUj7LIM1xgxwrAvobxjL8MFy6sIffMmZy4JxxFOzpXeO2FSx1iucR+FEhvLOK0FHabpGs9qrnFqdatIuQvrMnskXVifO0LShyV9UdJuSVeWtOmK0JlbmsojptxSZSdlmow5FGZVmeGA6tTqVo8i6WiqMc3PB84A3thQmL9pZt8NPBf4fknntzXonw7HcYqYoSJMrW7V5Fxgp5kdMLN7gZ3AeWb2gJl9vJLXHgJuoZrqm8VjhM5C0BZDm1DqKs6SsWe3hHW3zckOZRhCByW3SdKuxvFWM9taeG9qdasmqZWuHkXSkcCPAr/b1qArQmchyCU5cudiSmhIzHCsKXZdFHaXa3M0fOYeM1tLXZT0UeA7IpdeX1h/dqUrSRuB9wFvN7M72ipzRegsDKm5yZNrTdoyw+F9qSx1829ThvDemHw52Ur62pQrVlcXOcZMloyBmb0odU1SanWrJvuAMxvHJwE3NI63AnvM7HdK5PEYobNQhNnk5pCYpuJq++I3kyttlmaY9Q3bTA3LaZ6PtZ2qN6wjNWwodxzKt2Crz6RWt2pyHXCOpKPqJMk59TkkvQV4KvDLpQ26InQWjliGOFQo08oSr/LwGZiZIrwSOFvSHuDs+hhJa5LeBWBmB4A3AzfVryvM7ICkk6jc69OBWyR9VtKr2xp019hxnCJmNbPEzL5BfHWrXcCrG8fbgG1BmX3E44dZXBE6C0+pBZiL38VicSm3si02F5MvvDd8Hx6Hdcfc666xyqEs8xS7vvsa/9d65Patkj5Yp6kn1y6v9zW+XdK50xLccboSi9+F12IudtcYYSz+F5Zt1tMWIwyvNc/l+hebVTOEZV6huu++xjuBZ5nZs4G/AS4HkHQ6cAHwzPqe35e0YTRpHWeGdLWk+mSHF4nJwqwruZ2nmX1C0inBub9oHH4KeFn9fgtwjZl9E/iSpL1U018+OYq0jpMgZvXEXM+YG5lzg/u4xrlETso1zrn3ba5/iXxjWIXL7BqPESP8WeD99fsTqRTjhCeM9p7g+xo7zmLhijCBpNcDB4E/npyKFEvua0w16JG1tbVoGccZQmhJxWJvkN9wvnkuZwnG6m27NyVDk7Y4X1sCaHI85jjCZaS3IqyXvfkR4Kx60ybwfY2dOabENY6VH8M1Du8N34fHbfHGEnfZs8bl9FKEks4DXgf8kJk90Li0HfgTSb8NPA3YDHx6sJSO05OUcmseT8rFZobkFFublddmXZZYhCklHas/JntKjj6s9C52iX2NLwcOA3ZKAviUmf28me2WdC3weSqX+RIz+9a0hHecNnKKrnl+UeniYg9lpS3CxL7G786Ufyvw1iFCOc5YlFh4fetLkXKdc/Xk6u1ybSzrL8ZKK0LHcRxwReg4S0Gb9ZdKUJSM94tlnsNZJ7GyqbpLxgB2Sdr4OMI8rgidlSE27axkIHPbkJhU2ZRrnCuTimHmEimx49T7obgidJwloMs0uLGURyp5MYbFVhojHMMiHHNh1nnDFaGzkrS5t33uzw3VabPOwrJju8Zj4K6x4ywZKfc2LNNl/m8f1zgnW7OdeXCNXRE6juPgitBxlpKmtVS6GkyurtLzixgjXGaLcHGH1DvOSLQNtJ5kmmPXYmVT96WuNxVx+GrKl7qeOo4N3RnCrBZmlXS0pJ2S9tR/j0qUu7Aus6de+yC8vr25oHQOtwgdh3QsrmQKW9tc41QcsmRqXNdrXRI0XZlh1vgy4Hozu1LSZfXx65oFJB1NNd13jWqFq5slbTeze+vr/x64v7RBtwgdp6ZpeUFZZrjkWolbOu0pdmNlkGe0VP8W4Or6/dXASyNlzgV2mtmBWvntpF5JX9JTgNcAbylt0C1Cx2lQmvVtW8Umd39uZknqeirOGGsnV3YIHWOEmyTtahxvrdcgLeF4M9tftWn7JR0XKXMi8NXGcXMR6DcDvwU8EN6UwhWh4zjFmBVblveY2VrqoqSPAt8RufT6wvqji0BLeg5wmpn9SrjFSA5XhI6ToM01blpgKUuya8IiZX2m2hya6e6GAeOsqmdmL0pdk3SnpBNqa/AE4K5IsX3AmY3jk4AbgO8DvkfSl6n023GSbjCzM8ngMULHSRDmaOGJrmzzfPiKKcTJfbE6YscpmZrH4fXpYcBDha9BbAcmWeALgQ9FylwHnCPpqDqrfA5wnZm908yeZmanAD8A/E2bEoQCRRjb17hx7dckmaRN9bEkvb3e1/hWSc9rq99x5p1wOEp4vnmt+QoVYvO+8P7YcaxsTK7Y8TSGz0xaKXsN4krgbEl7gLPrYyStSXoXgJkdoIoF3lS/rqjP9aLENX4P8HvAe5snJZ1cC/mVxunzqZbn3ww8H3hn/ddxFpbU8Ji26zk3umRYzqT8vEyxG9M1zrZi9g3grMj5XcCrG8fbgG2Zer4MPKukzdanY2afAGKa9irgtTx+l7otwHut4lPAkbWP7zjOwjNRhCWvxaLXz4SklwBfM7O/Di7lUtphHRdL2iVp1913391HDMeZKUPibyWxv771ltQzXuxwORVh56yxpCOoUtznxC5Hzvm+xs5SkXIzU7HA2NjC0rpziZrYuVQmexxm4xqvB32Gz3wXcCrw1/UOdicBt0g6A9/X2FkBUjG3XIywjVgypuRarJ1pTbGrFOHDI9Qzf3RWhGZ2G/DoSO96vM6amd0jaTtwqaRrqJIk901GiDvOMpHLBIdlSqbqLcYUuxW2CGP7GptZajvPHcCLgb1U01teNZKcjjO3dLW2UooxN92uj0XXlu3ux4oqwsS+xs3rpzTeG3DJcLEcx5k/VtgidBwnT2qsX2pQcy4OWFp36r5pLsM1aWEZcUXoOCMQUza5GOGQxEqsjbZ6PEaYxxWh44xE16WyShdM6HItNlRnPGtwMtd4+XBF6DgzIOYal0yxa7uWq2dRp9itB64IHWeGpIbdhMepjG/bwOzpzyzxGKHjOCuNW4SO44xAW5Kk63Fb/ePjitBxnI6UJFAm59qG4cQoVXzjxQg9WeI4zog0lVgzoREbhpOajRJbaKF0rcPuGB4jdBxnMKnkR5PU9dzqNuFxLHs8Du4aO46z0ixvssQ3b3KcGRIbAB2+YgOkc1sExI4n9w1ZtOGJzGaFaklHS9opaU/996hEuQvrMnskXdg4f6ikrZL+RtIXJf14W5uuCB1nxoSub2zTpT51xY7bzndnJps3XQZcb2abgevr48ch6WjgjVTL/Z0BvLGhMF8P3GVmzwBOB/6yrUFXhI6zTjQXZSiZOZKas9x8H9vFLnVvd2a2necW4Or6/dXASyNlzgV2mtkBM7sX2AmcV1/7WeC/AJjZI2Z2T1uDrggdZ51ouq85RVW6zH/KNR6PmW3edPxkQef673GRMtH9kSQdWR+/WdItkv5U0vFtDfbe11jSL0q6XdJuSb/ROH95va/x7ZLObavfcVaZNouwTZGVurzjLsxapAg3TTZnq18XN2uR9FFJn4u8thQKktofaSPVFiF/ZWbPAz4J/GZbZb32NZb0w1Tm67PN7JuSjqvPnw5cADwTeBrwUUnPMLPlTDU5zkrRaRzhPWa2lqzJ7EWpa5LulHSCme2vtwO+K1JsH3Bm4/gk4AbgG1Sr43+wPv+nwEVtwvbd1/gXgCvN7Jt1mYmgW4BrzOybZvYlqiX7z2hrw3Gcx+970nRrm3G+8NU8PykbG2A9DjNzjbcDkyzwhcCHImWuA86RdFSdJDkHuK5eJf/PeUxJngV8vq3Bvk/oGcC/lXSjpL+U9L31ed/X2HF60pxlEhtOkxpqEyZFQsU5rjKciSK8Ejhb0h7g7PoYSWuS3gVgZgeANwM31a8r6nMArwPeJOlW4KeAX21rsO+A6o3AUcALgO8FrpX0nfi+xo4ziGnNChlvit305xqb2TeoLLnw/C7g1Y3jbcC2SLm/A36wS5t9FeE+4M9qM/TTkh4BNuH7GjvOYEIrMDXGMLWWYezaOCzvXOO+T+t/Ai8EkPQM4FDgHirf/gJJh0k6FdgMfHoMQR1nlcgtoJAaThMOn2kbn9iPmbjGM6fXvsZU5ui2ekjNQ8CFtXW4W9K1VMHJg8AlnjF2nGVheecaD9nX+D8kyr8VeOsQoRzHqQitwrbltVLL+vueJXl89RnHmXNyCy/ElunKZZKH4QuzOo6zjqSUYW5VmtKped1YzmSJK0LHWRBSFmGJ4nPXOI8rQsdZQNrc3bZNovrjitBxnJXGLULHceaIcNB17Pp08Bih4zhzRJsCHH/4zCN41thxnLkktnvd9KbbuWvsOM4cklqgYXz32GOEjuM4eIzQcZy5JrcBvI8jzOOK0HGWjDBO6DHCdlwROs6SkZubPAzPGjuOs0CMuzx/E7cIHcdZIKaTNV7OZIlv8O44Tgemv0K1pKMl7ZS0p/57VKLchXWZPZIubJx/haTbJN0q6SOSNrW16YrQcZxCZrad52XA9Wa2Gbi+Pn4cko6mWi3/+VRbBr+x3tpzI/C7wA+b2bOBW4FL2xp0Reg4TiEGPFz4GsQW4Or6/dXASyNlzgV2mtkBM7sX2AmcR7WTpoAnSxLw7RRsIDcXMcKbb775Hm3Y8E9UG0CtCptYrf7C6vV53vr7z4fdft918OetbmbN4ZJ2NY631lv4lnC8me0HMLP9ko6LlInuoW5mD0v6BeA24J+APcAlbQ3OhSI0s2Ml7TKztfWWZVasWn9h9fq8bP01s/PGqkvSR4HviFx6fWkVkXMm6UnALwDPBe4A/htwOfCWXGVzoQgdx1ktzOxFqWuS7pR0Qm0NngDcFSm2j2p3zQknATcAz6nr/9u6rmuJxBhDPEboOM68sR2YZIEvBD4UKXMdcE6dIDkKOKc+9zXgdEnH1uXOBr7Q1uA8WYSl8YNlYdX6C6vX51Xr71hcCVwr6SLgK8DLASStAT9vZq82swOS3gzcVN9zhZkdqMv9OvAJSQ8Dfwf8TFuDqvZldxzHWV3cNXYcZ+VxReg4zsqz7opQ0nmSbpe0V1JrdmdRkfTletrPZyfjq0qnEi0CkrZJukvS5xrnov1Txdvr//mtkp63fpL3J9HnN0n6Wv1//qykFzeuXV73+XZJ566P1E6MdVWEkjYA7wDOB04HXiHp9PWUacr8sJk9pzG2rHUq0QLxHqqR/U1S/Tsf2Fy/LgbeOSMZx+Y9PLHPAFfV/+fnmNkOgPpzfQHwzPqe368//84csN4W4RnAXjO7w8weAq6hml6zKpRMJVoIzOwTwIHgdKp/W4D3WsWngCPr8WILRaLPKbYA15jZN83sS8Beqs+/MwestyKMTpNZJ1mmjQF/IelmSRfX5x43lQiITSVaZFL9W/b/+6W1y7+tEe5Y9j4vNOutCKPTZGYuxWz4fjN7HpVbeImkH1xvgdaRZf6/vxP4LqoZDvuB36rPL3OfF571VoT7gJMbxydRsFLEImJmX6//3gV8kMotunPiEmamEi0yqf4t7f/dzO40s2+Z2SPAf+cx93dp+7wMrLcivAnYLOlUSYdSBZO3r7NMoyPpyZL+2eQ91XSgz1E2lWiRSfVvO/DTdfb4BcB9Exd60QlinT9G9X+Gqs8XSDpM0qlUiaJPz1o+J866TrEzs4OSLqWaI7gB2GZmu9dTpilxPPDBank0NgJ/YmYfkXQTkalEi4ik91FNgt8kaR/VopnRqVLADuDFVAmDB4BXzVzgEUj0+UxJz6Fye78M/ByAme2uFwD4PHAQuMTMlnMDkAXEp9g5jrPyrLdr7DiOs+64InQcZ+VxReg4zsrjitBxnJXHFaHjOCuPK0LHcVYeV4SO46w8/x+h4IgRZE+IUQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.imshow(covariance,cmap='seismic',vmin=-0.08, vmax=0.08)\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The correlation matrix can be constructed using the covariance matrix and also give some insight into the relations among the parameters." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAD8CAYAAADt2MYTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvX+cZUV55/9+bl/uNE07NEMDA46TAQckBFmio7AxX2P8FUxYcHdJvujXhPjjReJX47p+80OTXXNDQl4km9VdY6JBQ9D8kLjsJsGVXYO/1nUVBZTAgAiTccQWmqEdmknT9LT33vr+UVXnPqdunXPP7Xt7+sfU5/W6rz516jlVT9U9t7qep54fYowhISEh4VhHba0ZSEhISFgPSIthQkJCAmkxTEhISADSYpiQkJAApMUwISEhAUiLYUJCQgKQFsOEhIQ1gIjcICIHRWRvQb2IyPtEZJ+I3CMiz1N1V4nIQ+5z1ah4WrXFUEQuEZFvusG8c7X6SUhI2JC4EbikpP5VwNnuczXwAQAR2Qb8JnAR8ELgN0XkpFEwtCqLoYiMAX+EHdB5wGtE5LzV6CshIWHjwRjzBeBQCcnlwEeNxe3AlIicDvwEcJsx5pAx5gngNsoX1cqoj6KRCF4I7DPG7AcQkZuwg7s/RjwhYk4GZjgRgOc9bzfLy7Cl8zQcfzydjqWruaX7qafghBO6z+vy8rL922jYv08/DePjIGLLi4tw/PG2vLjYpavXLe2WLd1+9LOLi/ba8+FpReDIEXv/uONgbAyWlrp1vs8tW7p1Ia2/Bluu1+3Ht+Nx5Ei+Xc+fv9ZjWVqy5Vqtd06+/3173/cJgDEgwuKiLY6PW5pwLL6fsbHuHPg50bR+nL6sefXz6efk6ae7/Pk58XPg58/zo8c5NmbnpNEAwXpSdYxQq9nhaP7C7yXGX6PR+50dd1yXVs+1fqfCOSqD/l78tZ8//f75PvW7oN+TI0fs/aJ3ql6H+pidk1ZbqNfhrrvumjPGnNKfyzh2i5jFirSPwn3Akrp1vTHm+gG6eybwHVWecfeK7g+N1VoMYwxfpAlE5Grs9pedO3fy+ocf5rd4HQBf/vL7mZ2FnQv3w7nnsrhkV6eJ8Q4dauzdCxec38naumdvLSs/PGNpd+ywdfv2wa5d3R/W3r1w7rm2vHdvl25qCh54wNL6l/6BB2D37i7t7t32/vh4l7ZehwMHyPqcGO+w/0CNHTvyfe7aBZOTlh9POz5uy/4abFtTU/bjefc4cCA/Fs+fv/a0WyctD9u323ZnZ+39M7bb+Xv8cXjGM9TiTgdaLTr1BnudBmf37i5/u3ZBo27n9/4H7NgmJ/N9Tox3eHBfLaPdf8B+D75e86rn08+lp/V9bt9u+9BjmxjvZP37cT48Y8fZwK74i61G7j3x/fu59/PnvwfP7/0P1HJz+/Wv27rTTikei5+TB/fV2L3bzmMnImzV6L6rD+6r5b5736bmwdcBTE/Dtik7t3pOZmbsO+K/a0/ry1NTsG3SzsmhhQbbpjrI2Ni3e5gbAIvAL1SkbcKSMWbPEN3F/rWYkvtDY7UWw74Mu/8S1wPs2bPHNL/1rezfXoP3sLM1wz2t8zgfmMD+P+owQY0OFyzdSYcXZm1dsPTVrLyztd/d3QXAOUt7gXPxQ72gdTdwPlB317sc/RTntfYCuwG7Spy3FNKe62jHFW2dc1ruzW3toMNWzmo96NrVfe4GJjmn9UBGS2vcPtvaAdi3/KzWPmhNAdOOdnc2Tku7C+qOv5YfG1zQ2mvrgA7bHA87gHHOaM24+zup0eG09iy0JrM+O9Ro0aBBx/GKa3c8m78ODdfnPW5sk2oO7HdzTut+YDcdGq5/HL/1HK/ZfLZ2q3YcrZ+j1nZgSn0P59JhQvUPHba673sHy46/CRbte9JatnOSzd+kms9693tw/J7Xusf1b7+zi9p3g+ymw6ndsbi2/Fj8nJzTuh9au+nUG7mFz0MvkOcs3YN/x85Z2uvmADp1PzZXl/G3nQ7b3Nx235OdrQPAtHrfPK0vT7PMNgC2cYiOux4GwlE9cZ0BnqXKO4BH3P2XBPc/P4oOV2tsRQNJSEjYoBDsv4oqnxHgFuDn3KnyxcCTxphHgU8BrxSRk9zBySvdvaEhqxG1RkTqwIPAy4DvAncArzXG3Bej37Nnj7nzq1/Nys2xMX79iGFpyW77NTrUesQRX9b/lXW5H22s/bJ2ip6NPb8WKOKh7D70H1cZTThPqzEHVb+bQfsf1fvQb3417TA8rIS/DjXGxuSuYUTXHSLmlyrSvhNK+xKRj2F3eNPAY9gT4uMAjDEfFBEB3o89HFkEXm+MudM9+wbg111T1xpj/mwFw+nBquwMjTEt4K3YFfsbwMeLFsIYmu02v7tFOPHER3MvUvZFLyxQWzhsP3Rgbo5aa9nSzs/D/Hx3wZydpdZapkbH0rpyh5pVvCwtdTXQc3PUlpSK2NHqa12m1cqeY24u105W52l1XRltWKeh+wzLkbFkdW5OMiws5NvRmJmxH18f69P34etitH6cIa+aX/2cflb3Ebbr63z9/HyeP30dzolXnuo+Y/y1WlZBu7DQMxb/DuUQlssQjqXf/Gne9ZyE70lY1rT6ekjUKn76wRjzGmPM6caY44wxO4wxf2qM+aAx5oOu3hhj3mKMebYx5rl+IXR1NxhjdrvPSBZCWKWd4aAId4YezbExmu127wMLCyzWtwLulHFpkeX6RKbQhpL/pq0W1Ou2vrXc1ZirOl3O9EDuGroHDjnasn4LaIeGblddF+3SKu0M/Y8x4DejCedM1Q+0MxxkTvR3NsjOsGIfZe9Cv12Z7zPkJ0anaaM8rPOd4bNEzNsr0v5yn53hesRqHaCMBM12m+bYGFe68i5g/E/+hM6bruaD/8nee8fbOywywcTSYZic7NlJQv4F1Epuv7hliJSzHaZ+zpVjiJ4oFtAODd2uug55CH8UIaqMK2urwriL+inkvR8UbUzl4a99vf7OqqBGp/BdKBtL0aKzEtp+z5bNZ1hXdZFeCTazy9q6XgyhuyACXAmce/fd1OYP8aEP2dOxd7zdbWYWFmByMvrSaehy1RckRlf27EpevCp6u0HaWvHOMNLOoDubot1KrJ2ynU1ZH7F2fDlGVzb2fotFlZ2hbie2GMXmuGxHt9I5KdsZDoujfJp81LHuF8OEhIT1g7H+JBsWG2Ix9HrD5tgYzV27oNXiGc8IiJwFcUxkqHLKXITYf/qyZ4f9LzyK/+BDicl9+Kky7lg/ZeWqfYR1fcXkPmOJtTNonW63aN5jPPQTk6uK0LF+YnMyCghpMVw38CLzr771V/ngB+29DjW2TnbosK1HBCn74ReVh1k8y0TBMtoYr1XEyio/mCqiWFWeNwKqvAPrERuFz43B5cqw4cbWbLf5/ROE5z//Zp7//JutvnBpidrCYYDM/GVxyf7YawuHWW7VunqbuTmWW65utmsH7uu0/ofZ2e4Pyl1nn9lHutdzB6nNHczaqs0dzP8oZx/pmvd42lYrKwP5tpz5jm8na8vxkGvX18/MZNcdatTmD3XNgCKmNVldCGdak+OnaCyq/w61jDa7VvzneI214+qzcStzqGidnyM3zuw7c6fh4ZzosejvIePX8+farh3Yn7Vbmz9kP5oH9Z7k2tH96e8u8m6E32+UP89DOCfhexKW1ZzkzMWGgNcZjsK0Zj1iXZvWlMEfqjTbbdi3j3uWzuH888kWRep1OuMT9kXQXv9LSzA+bl8SvSDU61HTmpjpSlbWz/k+fbtlJjshivpdyXMEJ+aUHw54Ggh2yUU8hKY3uk837tzcls1Z0EZPuyV9hGOL7azLno2OI5zrpSX7DoVjUWY3QH7cfea6cA6KeCh7F/q8rzlTMEDGxoYyd9klYv59Rdo3JdOahISEzQrvjrdZsWHHpg9Vrv6u4YLZL8DcuTw4fypg/0med26HzvgErRY0/H9stysEev/Thv/ddVnvBFxZPxe2EyvHdg4datQUfa7cz2C4X78BKu1cwrYqjq30erVpy1C1nYiRec2H9RlkDiryE37vhe2UvQtVv4sRYqOKwFWw4cfWbLe5/pmC/Ng3YWqK6WkbysiHOqotLbKwoGy9ZmdptZRbVeAOldMZOp0TkNM/+XKG0J0sdKNzrmdZn4pW6y2zsnbP0igr+zhiHs4dr0anV2cYurA5+DH3jI38/OVc6sI+9bUuK9qsH92Oq+/pI+xHfWfZ2BRtbWkx/x1pWj1fir9cn55+34NR18KMVqOPO57WIHrakJ+edtV7Ep0T9/16XXdPWdMuLfWoEVYCGeCzEbFhdYYhmmNjNL/9be5f2AnAP/0TXPQC+6NYbtVyrnqgdkmhvk/tnIpOX7Oy0hVBscGvfsajX7uxcFBlbS+3allMv6peIiFfRfVFPBTxU7gD7jMvpXMSmVt9v2isXt+3Ej1q0fdQle++7ZfQhvrIsueL+AjvD6szPEvEXFuR9rVJZ7h28GY3nX9vF/fpaXjBC+xL0aj3mrFkL06JG1v4QvXUuR9K7McQe9ljP+Sydj2dv1+22Nox5hXmIcrGEuO/bE5C3nWfRQth2dxW6TM2jqLxZn0GhwixBTUcV1YOFsKyA5sy/ormr6gu92w9/p0OMn8xfleCze6BsqnG1my3qf22UPtt4W17vmRfkoWFrikNnZz5DPPzVqzyotX8fD7CjYt+E15rWsC6AjpzlQ41e61fXi3KKVpfzvW5sNA1hXDtZG05HnLtqmtvFuJ5yKAit9ToWFMk1acflRapw3GH/Wh+svl0fWbXrlw2f7k+fTmck7DPyDgzMdSJirk+An6zuQvKOesCF8FIj6Xsewj/GRTu+grmr6cceU9yom9YF3mnMhRFKFoBNrNpzabZGSYkJKwuNvtp8qbRGYZojo3RfOopOuMTzM7a/B9A1P6qDGW6mX70a4Uicb1MxwfV5mPQPkeJfjpd6K+3rdLuIHW+z1HSDsJDv+d0n8PqDHeLmP9YkfbVG1BnuOIdrYg8S0Q+JyLfEJH7ROTfuPvbROQ2l+D5tlHlNB0UzXab5gkn9B6c1uvd02TyJ33+5dHlUIyM6Zt0OdZuDKGOKey/qE7fKyoX/XC0LikcY5meqQp//XSGRe0UjWMlfeb0fSV9xHgrqiv6vvu1U8R72G6VOSnrs2z+ivhbKTa7B8owfLeA/88Y84PAxcBbXG7kdwKfMcacDXzGldcEzXab950sXHGFuhmYcGidoTdryV5PZ/ai63JmOF4v40wbdLTtnH4sEsFYP9ujH/M6Q6ffzOqUuUzWrubP6QY9D5lJidMpZT/CQCeXW1yUSYfvM7e4FUT/7lCLm3vQNSPJXN10WfWZlZ2+To+zZ750ZHMKdIZhdHI9Fv89+LLvw3/HMzPUFg53+9B6wXBOgnLuOwsR++4j70LGn34XtFudfk+UzrWnDkYa6TqZ1lRpSOTvsDkL3g+8xBjzqEv6/HljzHPKnl0NMVkjFzF7YYGH57eyc0envzvUWke6rsBfv0jXfXkjEKcKIl1X6rPMLa0f7xX61H2E49FjyO5X7DNst6fPgrFkfVZwx4vShmMpM6MpqCtD7ntleDH5bBHzvoq0P1lBTBaRS4D/jA2G82FjzHVB/XuBH3fFCeBUY8yUq2sD97q6h40xl1VkrRAj0YeKyC7gh4GvAKe5LFa4BfHUgmdyeZNXEzpAbPOpp9i5fRmo975c6gUH8h4CA3iV9GCl3gAVPGTCuqHRb1yDeN6U1Q3iLeHqa3RW5O1Rqc+KY+lQ0XMkQNF3lqsfZP7WAKM8QBGRMeCPgFdgs2neISK3GGPu9zTGmH+r6H8Ju8Z4PG2MuXBE7AAjEO9FZBL4r8DbjTGHqz5njLneGLPHGLPnlFNOGZaNhISEVcaIdYYvBPYZY/YbY5aBm4DLS+hfA3xshaxXwlCLoYgch10I/9IY89/c7ceceIz7e7Do+aOJZrudHarIFpe0PJI9rSeTnv+PPjMTzY6n9Ts5vWDoThZmMiuiDTOiaWjdmeIh48/pP3vaXUF2vNw4UTq5NcyOl+t/A2XHyx2S6PdG9anrcn32yY6X08/qcpAdbyRSAyNdDJ8JfEeVZ9y9HojIDwBnAp9Vt8dF5E4RuV1EXj3QIAqw4l2vy2v6p8A3jDHvUVW3AFcB17m/fzcUhyNGs92GsTE6GGrbt+fFj+npvNfB1FRXXHR1QFbOnp2asn91Wber++lHOz0dv/blkNaXw7Hodn2fHt5xO4bt2/P8OR6yOdF9lvET9hnSan79OKu0G9aF/RT1ofsJr2Pt7thR3Efs2RLkTupdPx1q1KrMX9mcuLY1bVbWtOPj/VU5FTHAkjotIneq8vXGmOtVOXbOUnSAcSVwszFGp8rcaYx5RETOAj4rIvcaY/6xOnu9GEYF8CLgZ4F7ReRud+/XsYvgx0XkjcDDwE8Pw+BqwOoQhV8/Ymi4+Iedya3UnNlNve6+8vEJW0cNxifyL1Q9H9HG09WIZ1oLdUb62dAspl87GX2QHS9nBhLxT9bPhqefmUlGkBWuFugpO8FYcgcXyj2xqJ2QVuvgsvpQ/1nQZ+6Z4CClQ40wo6EeS48bZuiDrPqJ0Wr0jDVyqBO2W8RDVu7zTuW+y4L5y419BIvhgO54c30OUGaAZ6nyDuCRAtorgbfoG8aYR9zf/SLyeaw+cajFcMV7Z2PMF40xYoy5wBhzofvcaoz5njHmZcaYs93fQ8MwuFrwierlxOORE4+3L8uBA9TryjTiwH7AlWcezjdQYlrTE5V45uG8uBMzFfGiuTKhyEXTJm76kz3r+/Bud5oHF6nZ/yhqC4dz4mRuMZ552LalRLmc+Y034VFzkI3b8ZfNgYriHZq5hH2EpjbatMbTZn0qdUAY6TprP5iTUO2ho2tr3rWYnEVPV1Gx9feSzZluJ/gnk/vnFbw32fen5i9mWpPVaXOjwOwmZoaTzYl3KxwBRmhacwdwtoicKSIN7IJ3S09/Is8BTgK+rO6dJCJb3PU0dmN2f/jsoFj7I6o1hBeZLewOfN8+2L3b/0eHxx+HU06pUZucZLll7zfqHStqelHEi526HIpyPjZeKP6EtJOTvbS6nTJR2D8XthuKxePjxaeTMTE+rNf8FfET9hlTBwyiOqgyTt1GjFaPpd/cTk93n42NpaxcBtdPhxq1KvNXNu6QP10e5LS+IgQ4biQtgTGmJSJvBT6FNa25wRhzn4hcA9xpjPEL42uAm0zeBvAHgT8RkQ72Z3qdPoVeKY7pxTAhIWEwjNK7xBhzK3BrcO/dQbkZee5LwHNHyAqQFsNcxOx3tw3nHNhPh7Ns5dRU9z/h5GT+H+zkZFcs0v+xobesdxb9aHU5UpcTxYKdQk4U0nWD2NA5XiuNbYV1HVQU6X7tlLRbZJuXzU9sR1y1z4LvrIf32LNlGMH8DUw7wp3hRnW1q4LNPLaB0Gy3uWZMkGd/qXuzXmfbeNflKeY65a97TFk0YiYTRaY18/N5MxyNiBuY7iNnWhMJG5ahwLQm1N8V8VAY0bss0rW6zpmVaNqiOQlptTlKONf+vn4uxm+Z6Y9zx8uZ7+h+ysbdDzoCeb+51SYysfckNJ3S5cC0ZlTYzL7Jx/zOUMPrEGu81t6YmeGe1nlccL7dIS274KkNAp1haNpQpi/rpx+LtavbKSqX6dJiO4cBdYba3KgSf/14H0ZnGM6XRtFz/fgrM3NRfeRMWYraCpA7ROk3fwO+J7n3rUgPPcjOtQSbfWeYFsMA2nXvzbOGC+77LB1eSo3uP95GnUIxNOquNYjI2sdlrLI73iCuZ0W8roC/yn2uFu2o2lnpd9YP2mRmFeakr+pgSKTFMCEh4ZjHZg/uupkX+hXDu+59YLsgL7vP7ggXFhgfVxuGubnuTqxAP5ZzwdJ1oa5Ko0wnp90Dw2ddXVYf0yl5FOgMO9R6s+PFxlKUva+fTk73E7aj3P56MsEVuM3V6ER1htqFciD+Qne8yPxl7m5F7UaQ27Frt74qOtcyt82jnB3P8r15dYabNtL1qNAcG6N55AgsLLA4vg2AifEg1JV+8WIiUyxEVUEorFw5Fvqr33OwokxwEIhSmgfNe1k5DLUVG7OmHXBOOvVG3xBeHtlY9JzoMFj9xgV2ERmf6B3LCkJ4Fc5B0fjDe/3ehbDc6o3oPmwIr+eKZAEI+uGcDRjpejPvekcCr0P81acME5/7JACdV/0UNWCZBnXoTQIf/iBiC2QV/VDkudAdLoN6+TVyLnpVfqz9eCgo9+jABh1zv/qSOdHo8fhw9D26tH56XuiKASPQGcZc53rem4I5KasrLA+ixxwAG3XXVwWbeWwjQ7Pd5vdPEOTSGeTSmZx4pl2nMpSJybFoKRplyejLoqWESe5XELUmKibHotYU8bNaYvIASeQz17vwuSr8VRWTy8YdQeiOBwVickES+ayPJCavKpKYPAD8KfO728aGhJ/cal/qUEwuEF8LRerVEJMHEN1yYnJsLCGvRWMpE5PL2u03J+G4PN/BODU2pZgc8neUxeQLRMwnK9LuTGJyQkLCZsZGzW9SBWkxHABd1z2xCev33gPnn5/VZ/q8cJcQ0eUUhsWKtAOUthUtF+iMYruXnJ5thfZtfXWGg/JaVN+n3SKd4UC8e2xAnWH2/a6CzlCwERU2KzaqeL+m8Icq8s9c+LQw453OpAaVsuNlCOu0HrCPzjDneqaz4zmdV/ZDCXSGucWjis5wSNMaPSdZm7occ8cLniszrYnNSY+Orkhn6HHgQNQdbySmNWXueCF/A+oMo7RJZ1gJQ+sMXWKXO4HvGmMuFZEzsfkMtgFfA37W5TgoxEbRGYbwwR384pXpaJYW8x4MS0vdckxXFdET9ei4YrQF+sNMZzeMzrCKbk0vqKGesorOsIo+MdZnwTj9GLIxV+mzjN9Wa9V0hpW+X3cvpyst0i1H5m/UOsMLRcxtFWlP3YA6w1Es4v8G+IYq/x7wXpc3+QngjSPoY13CB3fgxhvhxhupveTF1F59GX984wSvfV1eVFtuuZ1VvZ79cHzZX2tRVYvJ+n4oAupySAu9pjUaWfBP/QOO8Bcra1rNb9hO1q+i1fyFdT19BDzp+33FwVi7YZ8RfrOyFpMHUQEEyH2/BfMVjq2n7ZAHd6/S/NH73a8Um3lnOBTfIrID+Cngw64swEuBmx3JR4CRJGtZr2i22zR/4Rdo/sIvsPXuL8AVV3DxxXDllYpofr77bsYSfmsxWUccCUXAUPzWUWzCqDphlJpI1JpoJGYvxnvxvGLUmp4IMkq0rByhpUwkDJ/zIqIeF/SK16HqwPMTi0QTithVo9aE31MZqkb9iUWi0fz1i1oTiMmjgHfHq/LZiBh2Ef9PwK9C9gs9GZg3xvhvYobijFdXu+xWdz7++ONDspGQkLDaGHGq0HWHFS/iInIpcNAYc5eIvMTfjpBGlZIuU9b1YHWGK+VjPcCfMjMmMPcfmZqC666Dyy51BOPjtFoqXQB0RSUdigl6gormdoJhcNeANueB0i/4Z5FoVxDivlIIL/1sv1D5etz9srtpHiKZ33IIn9Pqgn786Wd9lsCqYylB7jsMw/4XtRu2GUsPUe8N4RUNIbc22fE2HIbZ0b4IuExEfhIYB7Zid4pTIlJ3u8MdFGe82nTwp8yHDryDT3yie78zuZXG0iKd+gQEmdRiGe8ykxuXnc/Tl0a+DmLuFbrtuXJM2d+hG8VZ86ORy3IXHlzosZVkxwOoOd59O7rPHlpdLshqp8eS9an66BlLQYa7sE89lnBs4bM9KocYInMU4yGckyL+YmMhODgZlb7Qtrl5seKxGWPeZYzZYYzZhc1s9VljzP8DfA64wpGtu7zJq41mu822PxSuvVYdUMwdzP47h1njctd0qM0fyvQ9tflD+R96mc7QZUTLsqnpDHihznBpKa+ndNCuhZqfXJ9eN+lotd7R02bX84dy7XaPETqZCVGuT19WfWhaoJv1z2eU86leHbXXJeb60FkC9bjU2MKMhj5yuB5LOM5s3nRd7rgkv0hntH4swfebKwe6Zf9cVr+wkK8L29UZA3308yExajFZRC4RkW+KyD4ReWek/udF5HERudt93qTqrhKRh9znqiGHZtschTueE5N/2ZnWnEXXtObrwOuMMUfKnt+opjVlaI6NdcXnhQUee3orp53SyZtW0OmaTvgfTR93uyq0RSYeZfehwLQmRMz8g14zl0omJlUNg0M3v8g4C01rykyR3DjDdnV9Tz7lgB8Y0LTGobDPorldgZib+14Z3rTm+SLm9oq0jT6mNc4k70HgFdizhTuA1+gsdyLy88AeY8xbg2e3Yc359mDVcHcBzzfGPDHAcHowkoMfY8zngc+76/3AC0fR7kaGjpjdPHKE02YfhFN295hTxETEDP0Wi+AHGL78ZSgUkwsW1Ch/EX6qQPdTCVrvF5rW9HkuNyeDmMcU8BflvcJYot9NwbswDFa6cFbBiIO7vhDY59YLROQm4HKq5T/+CeA2n5NdRG4DLgE+NgxDm1kFkJCQMGKISKUPMO2tRdzn6qCpZwLfUeUiy5N/LSL3iMjNIvKsAZ8dCGkxXEX4iNnNLVuQ5zxmb/rscxHdWaY/q5odT0e+VtnxOgQubZEQXoU6Q+eOl9WXhRiLZZ8LbQf1c7qfWAgvXda2caHrnrYzDG3qAlvBnO1lv4jjfiwtF8LL61mruBaWIHeA4frJzUE4n7pP/S7oOXHl7LsO5z7IjjeSQxSRvPF82QfmjDF71Of6sLVID6HO7hPALmPMBcCnsXbLVZ8dGBvVPnJDwWfdg3Z26rvcqllTm3qdZRo24x7Y/MfjE/gT4px4pXM1Q48phj59zpmKxPIHF4lk09P2b2A+U2haU2DKUimzX0jry7GsdkV9hiYn+jqW8a4KP2BNayJjqdHJtxO2G0HulFnz1G+cZfy5cmiaFH12RNnxcn30w/e/349iBniWKvdYnhhjvqeKH8J6t/lnXxI8+/lqjBUj7QyPErwOcf/cVvbPbaWx7/5sZzMzk3elyzY3kcTmMbtDAKam8qeX/VzPiqASvXSoZX1k/AX2i9mJZUwnp3lQtHpsWZu6PuRvwCTyMVrfbtccJZiTmNmSN1lxy3TbAAAgAElEQVSp4I5XtvPKm/eUtKPnJL/LKtRpFn6/I9A/9mCwnWE/3AGcLSJnikgDa5FyS747OV0VL6Pr9vsp4JUicpKInAS80t0bCmkxPIpottt89NnCR58tyA/ts1LY9u2cNX24Ky0tLdGg66qXQ1gO3PFYWuqKR2XueEtLeXEyaDN0fwvNebQ4pk1r+rnj5RbriEidS1wfidCSu16JO144X33c8TKTlApJ5CsfWlRwx9OmNYUuimvgjket+4+r76cPnB3yW7GL2DeAjxtj7hORa0TkMkf2NhG5T0T+AXgb8PPu2UPAb2MX1DuAa/xhyjBIYnJCQkI1+J3hiGCMuRW4Nbj3bnX9LuBdBc/eANwwMmZIi+FRR9d1b4yt9afggQf4+7nn8cqXu52F1ueFblStFotLNZudD/K6wFaLRazOcNy3Q4HnSoEo06FGzbWZib4FbnTA0XPHC2nLxNsCHVzPWPq542l9nqLt0MeNLoKcznBQd7w+ulKvS/a0WblsjobBaojf6wSbd2TrHF6H+I4nDa9c+AIdXgxAbXaWh+tnsXNHJ6cT9Ir73Gut3fHCugLdVMyeLWeorBdR92zu0Eb/sEJdWUW9Ws8C3c/eTunSBrLxC3SGA0fijujronrKEdhM9rRT0Q4yHNeo7BWjGPHOcL0h6QzXEM12m/ecKMiP7ePxx+Hxx4HJSXZu74bhyqn7Qv3Y3FxXNaRcz3poVSM1OuU6pDBTXZjpT5V7XOwKdIb6uVi70eyCZSG8HLxbWg5FEbI9T77PfnpA/ayavx7dqGtX/1PIHRQRHKCEYbo0yuYkFsJLX4flInOjYTDaA5R1h43J9SaCN7s57Xjnzr0ED8822LHD/qfaWvd+pePQarnoN+7W5KQ7bLEvYL+oNeFJs0fux6rMRLxomXMbm5qydnL1enYNzm2tQEzOteP7D9uB3rJHWZ9lYvL0dNfPu17PeMqJlkEf2fPbt7sI5RO9J9gFUWGKPIByZddPOF9ajM+J9K0W1BuF/Opx9ZSLeFgpNvnOcPOObAMh57r37W+z89M3wBVX0Jncyr59luac3fYHkvvCtJjsTGsgHzW5TIwK/XmByuLZUHUrpB1YTC4TH1coJkfLFXRyldwtB+FP3+snJo9qARMZrf5xnSEthgkJCdWwyXeGSWe4TpC57v3ADyBvtLamrVaviWBO3Jmf76r/fMY7rx9yFZlebWGhu/tTdT3BGJSuyuvHckbXgR4wZ3StGY3oDHPtFOjHojo592zWrp6LiM4wZ8uoocaiMwaG7QJZhsBw/vSchO3qeSzUGYapBzRNmU6zwOa0Eq23yRwWSWeYcDSRue6NH6Ex9wi7d5/Rrdy3j0PT57Btqqv7m6g7neHCAp3JrfnGvA4spjMMQn5l9aH+aXw8rzMcH8/ppnRmwJx4FtGz5XSGk5PFOi/XhzZtyXjwdbE+XdvZ4jw52V34nTFwrh2UiiDsc3o6m1Pfhz5Rz6kWBtEZun4yfWcWSqzRM7dAXm8ZzEmpztA/F+NhpdjkO8PNO7INDK9D/Jn7DLh/8uedC0xP8/3vq8VL+y6HBx8R0xXtNhcL9dShN9J1GC06F6E6jCStdZgDtrOSPqOmQMQjXXvaIv56UioMMBaCco53OtGxaH5i7WR6XX8oFvJXr/fSqnLuMC00gVopNvliONQMiciUC63zgIh8Q0T+uYhsE5HbXATa25zvYMKAaLbbfPyHhEsvhUsvtYEdaLVsgFgP7eI2P5+PWqPNK1YS6TqIZp3rU9dpcx5Hq0W3XHQeJSiXRrqORPLJ9aEiXWuTlw61wkjXmj9/3RMpvLXcVTU4d7xwLD3tRMq5cYYoGEvP3Po+vUugfg56Miz2ZFzU0a1H5Y4Hm1pMHirStYh8BPjfxpgPO2frCeDXgUPGmOtcKO+TjDG/VtbOZox0PSroALHs28ct+87jUpdoSp8e9+xACsQivbsIT5OHEaXCdkcilg3YJ6xsLGX0/ep8n/1oB2l3GNoy/oaNdL1nasrc+eIXV6KVT3xiwyWRX/ESLiJbgRfTdZ5eBpZF5HLIwut8BBtap3QxTEhI2ABIYnIhzgIeB/5MRL4uIh8WkROA04wxjwK4v6fGHk55k6tBB4hl1y4uO39/VwRbWODwgvsKZ2ezQ8WeXYSK8hLuYvSO0ouZ4elj7hSYvE6sH2K0K9FfxZ4p2i1l9904RqIvK+GjyjP9+Fgpj+FzRe2k0+T+GGaG6sDzgA8YY34YeAroyXBVBGPM9T4K7imnnDIEG8cGmu02zRNOQJ79vzm8ULOL4NJS98C0Xmdb/TDb6jZbXMwdL1v4lpZ6dFwZQv1YWdgwZS7TQxvRs4Xt5K4jLndZ2y6itzbJyY1NI3Djy+nuylzYwnZnZoojhUdMdiqjwLUwWhcLXab7LOJvtUJ4bfLFcBiuZ4AZY8xXXPlm7GL4mIicbox51AVnPDgskwkW3uxmK//S3hgf5/HH4bRnLEK9zmGsac0k5N3UnGtXeDrqkdMtRaLLxNzxvEtdj7mMei4H5UaXc6ur17sugLrszUTAusZ5t7QyF0DyOYyz52K0MVc4bbqyY4ddRAJXuKI5KkPMHS90O+xpV81JjU7OtTDXjn4u4rIIjDSJ/EZd6KpgxSMzxsyKyHdE5DnGmG8CL8Nmtrofmy/5Oo7BvMmrDe26d/V3DWfM3UPnlAuozc4yxzYAtk5GIqIEKDStKTDD8dDJ1YvqgEJznmg7vv8g8X3uUMTXqXHFfuA9h0p6cQx493RF/GSLSAFthqBcaloTmMRUmRPNfw4hfWR8MZ5WDB/cdZNi2GX+l4C/dCfJ+4HXY0Xvj4vIG4GHgZ8eso+EhIT1gHSAUgxjzN1O73eBMebVxpgnjDHfM8a8zBhztvs7dDjuhDz8ocr1zxTkn/2TlYimpjhr8iBnTR4ks43ziOkMI3aGQDd7X8w2z5eDuqj9XWjHpzMC+rLOWhf06UXAGp2sLtMZhtnxNPRBkXuukp2h71/ZGWZZ5VbBzjC0vYzNLXNzOTvDQp1hqF/0KSA8ks6wEjYm1wlAV4fY4Ag88AC3zP0IYI20c3osHRzVi1CRFzbTTaFE1qmpvChWNQNeJINcTjcZZomDfLuaP00byxKnUfRcSBs+F9KWZccbQaTr6HP9MgYW9VlWB+s20rWIXAL8Z2AM+LAx5rqg/h3Am4AW1nLlDcaYb7u6NnCvI33YGHMZQyIthhscXofYfPJJLpt90N3dnSfS/q8OMQPtDjVqIW29nvsx1wgWVY2yH4prJ+un7LmCdqroDEuh2+3XZwEPWq/aj7aMh77tVGyz6hyM5ABlhGKyiIwBfwS8AnsYe4eI3GKMuV+RfR3YY4xZFJE3A78P/N+u7mljzIUjYcYhRa3ZBGi22zRPPBF5zgHkOQfszaUl68IHXTHZi5pKBPULXYealaaUSA3A/DzLrVo3IE6QTH25Vev2E0aB0eXZ2a6pDOQj0TgxOVcmb1qj+8whFJOd2VDPc2GfEf6ArpjsAklq17go75FyaFOYow0T2WuE4yyj7RcpXGNtksj3wwuBfcaY/c5h4ybgck1gjPmcMcb7Fd6OzY+8akg7w02CbqJ66GCoQTcKthd9/Q9ifLy7gGV0MEGrK1L5yCr1Oo3WYvf+1JTVY42Pw9SUrQOo20RWtdayPQENAyj451z/YTtAtzw9bRc1bVqztERnfMImUXI6sM74RK9pjYvUUmstd6NVe1odtca1nZV9u36cu3bZBXBqW2amlI3FmRvpIAq6HO7CesTkpSXbpos8AyXmRgsLdoxZJBrHn58jNX+58sIC+ChGalc+FAY7TZ4WkTtV+XpjzPWq/EzgO6o8A1xU0t4bgf+hyuOu/RZwnTHmb6syVoS0GCYkJFRHdTF5ro9vskTuRQMliMjrgD3Aj6nbO40xj4jIWcBnReReY8w/VmUuhrQYbiL4NKTNMaH5kY/Ad77D4V/6jSyPyszcBFu2wPHH11hY6P6T39Zyoub0dLZjbGBPPDuTW/M7m/GJblnFywOynUiNTrZDy1BWDuvqjZx+098L2+2xI6TXvhBVn4v3qK6j/NYbMLWtdCxZn3rcfcTRTtm4ddnzP9mI0+p7wVitoXjXZrMTzuVKMVrTmhngWaq8A3ikt0t5OfAbwI8ZY474+8aYR9zf/SLyeeCHgaEWw6Qz3IRotts0r7oKTj/dLnhf/CJ88YssLMB3vwtPPw2f/jQcOGA/7N0LDzxAhxqNhUM0Fg51s+nNHcyZdNTmD3V1jnMHc6YhtflD2YJQmz+UlaO0cwe7YbrmDubqa/OHciG8fF3WrjKXqS0czuk9vV40e07Tan4Uf7l2ff3sI922Fg5n/YTjzJ6lGy0n/Gj08BCZvxx/mnfNX1gXoc340+G8hsFodYZ3AGeLyJnOTvlK4JZ8d/LDwJ8AlxljDqr7J4nIFnc9DbwI6+wxFNLOcJPCnzL/6pVvoPH97wP2HGBhwb6rDzxgVWoAPP4oHHecvfb6J38wUZC83NfldkJl7ngqEG1Pu8pdEFzQWrWbqWnznjBBetiPOv3W7fhnc+0E/IVmOJl5kaPN+A/HEpbLoAPyxngIrrN5iIw7V1dGOyrTmhHuDI0xLRF5K/AprGnNDcaY+0TkGuBOY8wtwH/Aepf+FxGBrgnNDwJ/IiId7IbuuuAUekVIi+Emhl0QheZznwvAj8y+Hq67jj/e9wa+/nX4xV90hFNnw0MPAbC4ZH+oE9qnV+1udDkUv3Q5rIul9KxsohOYnOTc9NR1D0If7BL+eg4/ggT0ZbS63HdRVGPpEV/rvdGrY3U99X3c9EZykgwjz45njLkVuDW49251/fKC574EPHdkjDgkMXmTo9lu07z3Xpr33sv1v3MQLryQuTn4F/9CEU1OwimnsLSkckr53CGBp0pt4XD0GuxpcCbyuajb4BaNMLqMNt8Jk9zrOsh7f7i6XD9KZA29NHK0mnfFn68r8irxtGE7mRiqyjFPlKJx+z6jYwtNnFRdzxwV0GrVwUiQPFASEhIS2PS+yZt3ZAkZ9Cnzl/6P4YEH4GMf28fsrPVUuebOX4ZWi4mLL+axp+3J6Na5OTo7dlpxzNuvYU9OvT2gP2nOoDPORU5cc+K2PqUOT2NdXfasP9VVdeDEv+CUVYusy5PbaGha1WdoO9hzaj59aresxhI7YdflcHfYU1Zj03MSlvUJdfRUP2gnRqvnJKpGGBRpMUzYLPA6xA8Cf3Xhhdz/23cDIPwVcAVP1Y/jtOOtyPfg5PM4x/2AFpnAa4pqdLJymY4s1FeFtPr01aMqbameULXVqPfqFPv1WbqQ9eGvX11V2tw/gj78hWPrRzsUNvlimHSGxxia7TbvAVp3380CYDV1TwB380//ROZ6lqn4lpZyarhQt6fz1vdkZAsT2dPbTnbdR2eoaXN6Nu9hQlxnmIOO8K340+3qZ7PFytGGOrlouwrRxbpADxiri9EWtZMrK5fEnjkZBptcZ5gWw2MQzXab38EmsTkLsIEdvslJJ5Ethg891N1N5H5LwSmmf/f9qXD2g3QV+nQ2dwAQ1OWiMSsRtoY1bemKrMEpb9lJqnsux5MeR8iDhnKVy7nv0emlDXgID06KVAm6nZ523ZxkPIR9+oTz/jldrtfzX9qo3PH8aXKVzwbEUEu4iPxbbIgdgw2n83rgdKzT9Tbga8DPOkfshHUEHTHbZmx4vr0891wAXnJa9weae7XHx6m38gnnoRuFJTSX0aYiOZEwljC9iLYej66t+9HIxMQIbRF/oTlKUQQZ324O9Xik61B81W2F7fS0W9HcqG87AU9DIYnJcYjIM4G3YUPsnI81nLwS+D3gvcaYs7Hy1xtHwWhCQsIaY5OLycNyXQeOF5HvYxPIPwq8FHitq/8I0AQ+MGQ/CasAf8rM2BgX8yEarZ/lcMuenJ72vfvpnHIeABOtw1hHAHcw0VqEut3Z5RJCET8U8M8VHU5UOeQIDx70cxpFhwhlhxaxPnxd1UOOomvdh64P+dPlfjwU1cXKsflcMdLOMA5jzHeBP8DmOXkUeBK4C5g3xniFxQw2VE/COkaz3eZ2gHq9q/JxIaSiuveI7g2qnfKGKBLfwvuhGBrWh3q6zOwloC0TF8v6HIY29mxV/laLhxVjE+8MhxGTT8IGYzwTOAM4AXhVhLQoLE9KIr+O4BPVN+YP0pg/CB/6EDU6NFjmsae3Zj/aWmuZZYrj9kF+9xPWheXcoYqqKyqH9Ho3FO60imjL+gzp9VjK+gx3erFxVB1LVd4HnZOhscnF5GH+lbwc+JYx5nFjzPeB/wb8CDAlIn42omF5ICWRX49otts0TzuN5mmnIf/uldm28PjjlZteEI1Z/zR7Eq0H0axztDMzuXZi0axD2hodW+eYqdEpj3R94EB2u0anMHF9h1pPpOvavgejSdpztLqtEuTGXRC1O1rnxprV6W26SwDVoZZLCJWVgyTyI9k1+uCum/Q0eZgZehi4WEQmxIaU8HmTPwdc4WhS3uSEhM2CTb4zXDHXxpiviMjNWPOZFjZ5y/XAJ4GbROR33L0/HQWjCUcH+lClUzfUZh7m6S072TruDkq2b88OUJiezh+geLu4eqMbyt+LaCHt9u29Ie592dVlJiI7dnTTEKhMcB1qxWH/6djQ/Tr6jgurn1071OhYHnz/ALt3W2PreiMXWiuj1aiQHS+D66emxx3jwcdX8+kPItkGs92xmr9cOx7ahnNYbNCFrgqGGpkx5jeB3wxu78cme0nYwPCue+940nAah+nUlf+w93WN2A5miNi+5WgjNnW5E9/Qn9YbDvexv+s5edXG0pFrzV+ILCp1RZu/MmRjG9A+EGePGa2L2EZm5dVYtNJpcsKxima7zXtOFOTEb3cV9E5/5/VYOeV9gU4O6KVVesFo4noddkrTqiTyUZ2h3hEpvWT2rL4O+dXPHjjQdS9UOsOMVqOCzjDsp2e+wnbCxPBlSeQ9f7q8GjrDJCYnHMvwWfc6ziigNj3N0hJMjHdUqGy3sylIbF4jT5vVe3HR1xWJeVpsDsXkQETNBUx14jXQzT7n+Yklhtc7tV1ndRfDMCL1CpLIh+MOeQ95ioq+YZ++zj+ny6shJg+WHW/DIe0ME/qi2W5zzZhwzZhwuDXBxN1fynYaOTs5LYY6aLoc6vVs4cra0OJysLvQtIOgU2/YD7V8f+5+rP2sD//DD3Y7OkyYF+kr86PHWRYpPDYnRe0UzF/MtGdobOKdYVoMExISqmHEYrKIXCIi3xSRfSLyzkj9FhH5a1f/FRHZpere5e5/U0R+YhTDS4thQiU02+2uDvFFU3Gdl7bFC/RstbksuVlO3wjYTHSzj+RptZjnaDvUunUFOsOcrnFmJsuOV6OT8eCvfblDzfbfanXpD+zvhvgK7Az1c7pcBG1n6PnX/ITt+jnwtLlx+3qXpdDX5WjnD2V6wg61dakzFJEx4I+wjhrnAa8RkfMCsjcCTxhjdgPvxcY9wNFdCfwQcAnwx669obAx97MJawatQ6zRDUs4OUk3jJQ3s0GJzDpKTagfUzovlG6vpz7ryCJmWpP7IU5P50THWln2Pt9HYM4T0tboxDPrlSB3wl0wlp52y/SJ6tmo/rMgY+DQGO1p8guBfcaY/bZpuQnr0aaz3F2OjW0ANrTS+51N8+XATS6P8rdEZJ9r78vDMJR2hgkDw+sQmZpi6wNfZesDX7U/+PEJHpmz+jntieAXglw8Vbc4Zroup0P05Zw+LYzzp3cfZT/OWAxAfR1rVyFnWqPrytrtBz9W8ocRYVnvsvQcRdupQjsKDLYznPbutu5zddDaM4HvqPIMvXEMMhoX7+BJ4OSKzw6MtBgmrAjel1kuWkQuclnwFg5zxnZlZuPE6KUlYGEhv/5odzxHm+1etJmNo9VufZVNa5Q7XtYPBaY1obvbgQNdHkZtWkPEPdDzoOnKTGvCujLaEYnJxsByq1bpA8x5d1v3uT5oTmJdVKSp8uzASGJyworhRWYAWkfoTG7l9tvh4ovpiqStFhN1YHIy2xWOj9eobd/eXfy0GYn2TvHYvt0uTDHTmjIxuaJpTY4HjwIPFGA405p+XjAhP34eYvweZQ8UY0aXQQC7m3uWKsfiGHiaGRfv4ETgUMVnB0baGSYkJFSCXwyrfCrgDuBsETlTRBrYA5FbAppbsPENwMY7+Kwxxrj7V7rT5jOBs4GvDju+tDNMGArdNKRjNJ96ih+Z+zQdLmMRq2+7+Sa44gqYmKwz4YyYO0zkXfmU3V7otmbru25mWpeo3fbC0FkxF0C0HjLSTg/GJ+J19XwIs06E/xAhbXjdw5+6V8hf5PmMB+XOOJKTZEa7MzTGtETkrcCnsFHybzDG3Cci1wB3GmNuwcY1+HN3QHIIu2Di6D6OPWxpAW8xxrSH5Skthgkjgc+p8jP3GXa3YGLpEAC7dm2zkusDD7C4y0XOxrrj1bznideVTU/bRXBuLu8R4mgzExMnMnbqDWtGMjnZXZBay12/5tnZbmAH364ve33d9HRPuxlP09N0xidsH9Ctd+1kC/fcwW47EXFUG3Jn/ehx+rHoOq+jHB+3POg6sPVu3Bl/uuyeA6gtLXYPg4bECMVkjDG3ArcG996trpeAny549lrg2tFxkxbDhBHCB3doHjnCIbYB8P73w4UXQmNqKn+oqRcTtWBlbmqhiYxHWBfq0spMa2LmPEXtbt/eLffTEVbQGeZMa4rGqcv9dJpl+s+wvD51husOaTFMGCn8DrH5la8A8PFdN8Psm2DHDv77f7c0/+rVZLs37SYXEyWBnDjYIyZr20WCXCYl7m59xeQy0XcFYnLIQ6zPIre+Ku3m2onMySjQ6eRTTW82pMUwYeTQaUiv/q7hjPpB2LePev0CQHl8aLEYYGrKlp0IqMVkisTkhcNWPPQLQCAmo8Xk+fnursmLoVNT+XZ1n1NTPWJyRqt3tvOHunUFYrKH7se3k82JKveIyU4dsJZictoZJiQkJDgc04uhiNwAXAocdPmREZFtwF8Du4ADwM8YY55wrjL/GfhJYBH4eWPM11aH9YT1jO4ps9C86y6YnOTTn7Z1l13ayfR5HWq90asj+rJsdxWGBitzPZuezovVfkdZ0EcOOuRYmT1gUN8XMVvCWHkQPWUFneEosNl3hlWUCTdinaE13gl8xiWK/4wrg3W6Ptt9riblSz7m0Wy3aT7/+SzuOIc/+AP4gz/ohpjKFjhnnKbFybxJTDwcVZEuTIcUywWUdah6mDAqXVs/0XkQnqq0PQxdGUZsZ7ju0PfbNsZ8AWvjo3E5NkE87u+r1f2PGovbsZnyTh8VswkbE812m98/Qdiy5Wts2fI1fMTs7IcTc8fTP17njpdFxdbueD5ShEeBO17fSNehO96+favjjufMiEJ3vGxsIX993PEq0Y7IHc8foFT5bESsVGd4mjHmUQBjzKMicqq7X+RA/WjYgHPcvhpg586dK2QjYaNAu+51MFnE7K2TnV7RVyeLAtixo9gdLxQBtamKSgg1qDteZ/c5NtJ1iTtetmAP4o7nxlbmjhe62GX86QVujRJCbdRdXxWM2h2vsgN1ypuckLCxcMyLyQV4zIu/7q+PSrkqDtQJmwM+QOw1YwLj42zd+yVb4cI+LS51PTVyNnduJ6Vps4OSSHiqXNoBFfY/FgIro48FJY2E/S8Lt18ZReHHVFm3nY0nqC97fjXC/qfFMA7tQK0Txd8C/JxYXAw86cXphAQPb4coLzrZ3nDZ8cbH6brRtZZzesLM8FjpDIFoCK9MPxZEui7K3tfTbquVj3Qd0y9qhOUAOtK17yezg9SHR0pvmUWv9vq+QKeZq1M6w6y8CpGuN/tiWMW05mPAS7DBGmeweZKvAz4uIm8EHqbrP3gr1qxmH9a05vWrwHPCJkAYMRvseceuXV03Na2Ti0aKhnLTmqmpqGlN9LmK7nhR05qwHCDnQaJcBPV11m6ZO17EDTEa6Tq5460IfRdDY8xrCqpeFqE1wFuGZSrh2ID3ZX7zrFUrn8V+OpwF2CCi/rdf0yHuVfoAIJ783bvkTW7NLwKDRKvW/ZSJo7FyBKGo66PqFLajrqPqgALaaLsjgjEb96S4ClI8w4Q1RbPd5gPbhQ9sF+TZ83bnMT9Po67sA52JjBcBc6JviWlN7cD+7HqtTWsykbUg0rU2rcnaDEXfNTatOebF5ISEhARIYnJCwqrDu+4xNkaj3obJSQ4v1DK1Xm16mlYLGvW8PqxI1MxEUhfcAZQuzacW6OeO520Uj3bYf29fOUjY/346wxFlx9vsi2ESkxPWDfwp8yMLW9n6vz/ZrajXeeKJ7rU3lwGifreZWOoiVWt9XY+JjuojLBdmxyvQ31VacHRbsXaKFvdIXTiWnAmRNq0ZkQ5xs4vJaTFMWFdotttc/0xBLt3d/THPz/OMZyidofY3DjX63jSFTpYYXuvWspPdmB4wcOUbVGdY5itdSWcYM/0p0hnOz2fmR56/HL/uuVGa1sDmXgyTmJyw7uDNbpZb9pS5ceAAX1w4g1e+3EWp8b+2er376wuCtfr6Zez9Bp184vVwRxmKwjpqTZDfuCcjX4Uk8jlaLyaXtdPPhEiPpV+S+xGJyUcruGtRVKyA5kJsIJitQBu41hjz167uRuDHsHmWwUbPurtfv2lnmLAu0Wy3+d0twu9uETj3XF5Z/6ytGB/vFZMjImEHm5Tdr5W+nKFfEnldjtFqDBIiK2y3qJ18QvZe/kJxO1bWtCPAURSTi6JiaSwCP2eM+SFsVK3/JCL6P8ivGGMudJ++CyGknWFCQkJFHMUDlMuxjh5go2J9Hvi1PC/mQXX9iIgcBE4BApek6kg7w4R1C+/L3Dz5ZORl7jBjfj6/+1C6Ma8fc7dhfp6J8Q4T413dY4ZQD+ht/HzjBw50ZcJR2hmGtoQaBXaGWbmCnWFUv7g2OsNpEblTfa4eoJtcVEEjvlkAABjWSURBVCzg1DJiEXkh0AD+Ud2+VkTuEZH3isiWKp2mnWHCukfOda9ez0t9rVY3vBdkfxssw9QUi0t2ERgfx4b78vo6JZJqs5usHe2OF+rk+unzAvS440XMe3pMa/qZywxiWrM27nhzxpg9RZUi8mlge6TqNwbhyQWK+XPgKmOMH+S7gFnsAnk9dld5Tb+20mKYsCHgXffe9j3Dtn1OQtq9m+Xxrd24iB7ZggfoH68+CPHJktwi4c1wNG0sgo7OaqefzeroJoUKs9QB+QjfJVkAa0r3GetT86VpfTnkY70doBhjXl5UJyKPicjpLlaqjooV0m0FPgn8OxdM2rftg8McEZE/A365Ck9JTE7YMGi227zvZOGsS87hrEvOAaCxdJitk52uGOoXQjqwd2++ASU+ZlFhUNn6QtOapUVLO3/IftyCUps7mC18PqtduBACmXCcE1G9OyGd7Lms3dlH8hF25uZsZjvfhzedgbwpjadVZf9ch1p2PSyO4gFKUVSsDCLSAP4GG1n/vwR1PrygYKPw7w2fjyHtDBM2FKzI7GIILz0F4+PsP1DjrO124eLAAWpTU3S2n0Ftfp6//Vt7+8orUaJwoytK+nIsIVRETI4mkxrUA0Unkdcoi1oTRtVxaUNzYnFY9hiRaQ0ctQOUaFQsEdkD/KIx5k3AzwAvBk4WkZ93z3kTmr8UkVOwwabvBn6xSqdpMUxISKiEo3WabIz5HvGoWHcCb3LXfwH8RcHzL11Jv2kxTNhw6KYhHaP5F3/BWUeO8MglbwBg39x5fPkT8OY3w9Y9e3i1esM79UZXvxbq/ZxonendnM5Q6+t0O+GzvhwTR2N6vvA6LId6wJC2pnSPfXWGoW5yhTjmfZNF5AYROSgie9W9/yAiD7ij67/Rxo4i8i4R2Sci3xSRn1gtxhMSmu02zde9Dm66KbOLfs5z4Ljj3I/2z/+cP/xD+MM/tPS1hcOZiUxt4bAtuxOB2vyh7HS2RieLkt2hltFmuj2nP8x0cq7OL0DhB7rmNTndo7uOlgMXu9r8oa6uERu6LNMRLizYsl90XV3Wp9Y1DoHkmxzPm3wbcL4x5gLgQexRNiJyHnAl4K3C/1hExkbGbUJCgGa7TfO225ictGq0f/xHeOghV/msZ3H88XD88a7s9GzZtS6HOjmdgN7TehSY1oR6ueihRVk7YSTuMv76eciskgfKMZ0q1BjzBRHZFdz7e1W8HbjCXV8O3GSMOQJ8S0T2AS8EvjwSbhMSIvBmNwDN7dv5kaeegks/xmMXXcaUspW2O7fuNdjdgN7VZQtYYC6jEZrM6HJ4oqx3ixCIyardnn5Cc55wYa0gbo/atGazi8mj+JfxBqxTNdgcyberOp83uQcpb3JCwsZCWgxLICK/gTVr/Ut/K0JWmDcZax3Onj17ojQJCVWhD1XeBmwDTnvyQZpNa4945ZUuOOz8POjcKAsLdKa2dXVzXqRcWspEz54doPN4Ccua1qPsWf1cWParjo9RmKvz9bod1I4w5K/gYGdQpMWwACJyFXAp8DKXCApS3uSENYYPENs8/niYmuLii+39ej0ffkufGvsT45gHSlRM7nOaHBOTw2fD67AcRvAuO02OBYSNecIMi7QYRiAil2D9/X7MGLOoqm4B/kpE3gOcAZwNfHVoLhMSBkC2ID71FH/1/kPurj2cWBzfxoRenApc6YrulbnbxZ4PUURbpNeL8VDUfoz3Ij5WgmM+O57Lm/xl4DkiMuOswt8PPAO4TUTuFpEPAhhj7gM+DtwP/E/gLcaY9qpxn5BQgGa7TfOEE5CTH0dOfhyw6UcnWodzdN7lLosg421DZmerR60JE9mXoWrUmtnZfMa+MFH9GkSt2eymNSvNm/ynJfTXAtcOw1RCwijgo90A0DpCY3aWL83s5Ecu7obc379wKrt20TVN0aY3/l4Y4SaMUhMEaS06bQbySalU1O4sio2PWqMSQtV8WUfVmZ6GpaV49BvtdqjHNCSSmJyQkJBAWgwTEjY09Cnzm2cN557rKtzOa4ffVE1O5nVrk5NdVzifJQ+ni1NlIHfY4mkgb2PYtV9UtCpMGHRdAIGc26BOZ6CR6Ttj7nj6QCilCq2EFMIr4ZhAs93mA9uFk09u5/RsjZbTGWr9HFQK4ZVBZ91TiIrJYait1nJvmd4QXkXZ8fx1WNYhvFJ2vGpIO8OEYwZeh1hbeirTs33l6w1e8ALy4bMAdu3qLiChjjAM9zVICC/3bIcaNdVOFuk67MOvLAV9Fka61ivSiCJdH63seGuFtBgmHFPwZjfveNKaxp59tqsIbPNyiImoOde5uMF1kRlObJdWFJm6x/4xwkO0j5LnVookJickbDI0223ec6LwnhOFCy90N2dn89FdlJgcM60JxWToNW4uE5MzUVjLla5cWUwO+FttMfmYN61JSEhI8NioC10VpMUw4ZiEP2W2KQQMTE9zeKmRqQdrO3bYMsQz00XKRWLpUDrDonI/naHG2mTH23BIYnLCMY1mu801Y8Jiq8HWA/dk9zv1RmZLnYmbdMsaRXXhc0X1mQmPKuv7OpBsrL6o3GPWMySOlpgsIttE5DYRecj9PamAru084O4WkVvU/TNF5Cvu+b92yaP6Ii2GCcc8mu02v3+CsPVHL8hMXWoH9tNYsq573rQGVCY9hVyWvcA/OFzIdBa+LBuej6itsuNl2fucvq82d5Da0mK3pflD+cx5rl7TZn2OKDueP00+CsFd3wl8xhhzNvAZV47haWPMhe5zmbr/e8B73fNPAG+s0mkSkxMS6Gbd69TtKXNtfJx7Dmzl/PPpmtZo1zgN7WIXQW53ForC+iS6IDteJvquNDveiMRkOGpi8uXAS9z1R4DPYwPD9IVLD/pS4LXq+SbwgX7Ppp1hQoKDF5mvGROYm+OC7da4epkGyzQ4tGBjC2ozHO8dskyjR0z1iIUA88+FIm5I7+MZxhB7riis1xqcJk+LyJ3qc/UAXZ3mE8G7v6cW0I27tm8XkVe7eycD88YYv2wXBpgOkXaGCQkJlWFM5R3mnDFmT1GliHwa2B6p+o0B2NlpjHlERM4CPisi9wKHI3SVgkenxTAhQUH7Ml8jhlbLJp8DK31mtnxqx1ZbWqRRr+N/TjHf5OWWpW3MzsIOl+Zidha2n9HtfHaWzo6dtFrQmJmxbXlPmNlZ2LGj2+fsLDVfduG9ajt2sNyq0Zibg+1qnVlYgMmtI5gdA4wmIp8x5uVFdSLymIicbox5VEROBw4WtPGI+7tfRD4P/DDwX4EpEam73WHlANNJTE5IiKDZbvNuI9QWDrONQ2zjEHfcAcs0YG6OhYXuAQkLCzwy12C5Ze8sLtmPXziXWzUa9Y5NOzA11V1QJyfzvsmursEy7NhhP34l9iG7fJ87dmRltm+3n4UF28f0dPc56AkxtnIYYLniZyjcAlzlrq8C/i4kEJGTRGSLu54GXgTc76Luf45ukrro8zGsKG+yqvtlETGOGcTifS5v8j0i8rwqTCQkrEc0222aJ55IZ2obnaltXDT3SRp3fgmmprj5ZmXSMn0qZ0wt0qjbOxPj9uP1fY16EA3HxRnsTG61NLFyvW51im5H1xmf6ImW0xmfyHSYnjbTabo8L941b3ToVPwMheuAV4jIQ8ArXBkR2SMiH3Y0PwjcKSL/gF38rjPG3O/qfg14h8vOeTIl8Vc1qojJN2IjW39U3xSRZzlGH1a3X4UN9X82cBH2BOeiKowkJKxH6DSkz/s7w2W7HoHZWc4996yMpjZ3MAv5VRT2PxNv5w5Sm5qyhyjzh7KT6k69Ydvxp8Y+2b1bPDNafzAyP5/VhbTMz9twYD7E19Jib9ixFWF0YnJpL8Z8D3hZ5P6dwJvc9ZeA5xY8vx+bongg9N0ZGmO+AByKVL0X+FXyysnLgY8ai9uxsvvpgzKVkJCwHuEXwyqfjYcV6QxF5DLgu8aYfwiqngl8R5VL8yb7Y/fHH398JWwkJBwVNNttmu02X7tc4H/+T9ixg/e/XxH4kPwlyHR9oS2hE3Gjddq2MLQzdOXMrjAsaz3hCO0M02KoICIT2OPvd8eqI/cK8yYbY/YYY/accsopg7KRkHDU0Wy3ab7xjTz2RIPLLw8qIwtOzEWuH4rc8cK6GH3YT7YAjwxpZxji2cCZwD+IyAHs0fXXRGQ7KW9ywiaHj5j9K7+ibs7PdzPTKUQXo4LseN5EpscdT2W8y+ro5EJ6aVpf9m58o410bYDvV/xsPAxsZ2iMuRdlEe4WxD3GmDnnLP1WEbkJe3DypLckT0jYLPCue7WFJ+2N8XFu/fwEl1yS9/SIxjN09n9hlJrMRMajLKKNqo+2E9KOTEw+Ogcoa4W+i6HLm/wSrHvNDPCbxpiio+pbgZ8E9gGLwOtHxGdCwrqCj5gN8O624eWFJsS9iO3SwugyoeH2SnZ2scjZw+MYXgwL8ibr+l3q2gBvGZ6thISE9YfNvTNMHigJCSuEP2W+ZkzYvburIwzjCOZ2ZaFesCCEVxj2P6P17c8fyjxXshBevqzCe40yhJfFUTG6XhMk3+SEhCHhdYiH5g2Tk90zkp3b7eKU8wBxesFopOsCnWG0rCNdh+2sWgivtDNMSEjog2a7zftOFt63Rdh5/lZ2nr+Vv92yBdnyC1xxRT7Q63KrqwfsF+la6wyHiXQ9utPko+KbvCZIi2FCwojQbLdZANi1C3btYj8Ac8zM5O3/fOL6UKSOL3fldWXt6EVxNNjcdoZJTE5IGCH0KXPzFa/gxtlPcMcd1+Ej19fo5LxDegK/Rgytoz7OBQvcoPcHx8bUB1ZBWgwTEhIqYnPrDNNimJAwYugAsffyc5xySjfgkz8ZLkovOmg5xGjF4hjSYpiQkDAg7CnzGM1nPwT8n+z+4uSpjNMb4kvr/DTKFriqi99oD1A2J9IBSkLCKqLZbtO8/fbujVaLAwfsZehxUnSqXHSarJ8Lyxqj2y0akp1hQkLCiqEPVd7xpOG88f3ALoAsN4qOzhU7MAnv67rQBzrcYY7O4BqSmJyQkJCwyQ9QkpickHAU4F333nOiIM92Hinz8y4js3LP82Kud8fzrnoqRBeQ1fnrXHl+PntutO54R8fOUES2ichtIvKQ+3tShObHReRu9VnyuZNF5EYR+Zaqu7BKv2lnmJBwFOEPVfyC8eABuzDu3u1ymHjx1rvU6UjXGmFka4canSyvSoaRRro+KvrAdwKfMcZcJyLvdOVf0wTGmM8BF4JdPLGRsv5ekfyKMebmQTpNi2FCwlGG1yG+6nbDRXOftDenXgD1Ol+8c4IX/yg5JWKNjs14pw9CXH2HGjV1DVALadlwp8mXAy9x1x8BPk+wGAa4AvgfxpjFYTpNYnJCwhqg2W7zPy4Wdr75p9j55p+iM30qLCzwoz/qCJzoG0a6zqAiXUfFZBW1hqWlEXF91NzxTvNBod3fU/vQXwl8LLh3rUtX/F6fX7kfVpw3WUR+SUS+KSL3icjvq/vvcnmTvykiP1GFiYSEYxHNdps3fEd4w3eE2t1fA+Dtb3eVLslTFvGmYkIowIrJa58QatonfHOfq3UrIvJpEdkb+YTZZUrhsm8+F/iUuv0u4FzgBcA2yneVGVaUN1lEfhy7lb3AGHNERE5198/DrtI/BJwBfFpEzjHGbN4jqISEYwbezrAS5owxewpbMqYwNriIPCYipxtjHnWL3cGSfn4G+BtjTJZ4RaUaOSIifwb8chWGV5o3+c3YDPZHHI1n9nLgJmPMEWPMt7BKzYGTOSckHCvwp8zN5z+fw1M7efvbndF1vZHtCrXOMPvUG10j7YBWX680ZUAcR01MvgW4yl1fBfxdCe1rCERkn6tdRAR4NbA38lwPVjpL5wD/l4h8RUT+l4i8wN1PeZMTElYAb3bz7Gd/v2tmM3ewaxYTZMfzka6zKNgqk16YHW+0ka6PymJ4HfAKEXkIeIUrIyJ7ROTDnkhEdmGzcf6v4Pm/FJF7gXuBaeB3qnS60tPkOnAScDFWLv+4iJzFgHmTgesB9uzZE6VJSDiW4M1uOhhqzjzmcGuCSejVGU5O9prdFCWcH2mk69U/TTbGfA94WeT+ncCbVPkAkc2WMealK+l3pf8uZoD/Ziy+ilUkTJPyJickDAWfU+Vwa4LDrQm2HrjHVtTrmeseAK1Wt9xq5cuM2gXPY3P7Jq90xv4WeCmAiJwDNIA5rKx/pYhsEZEzgbOBr46C0YSEYwVeZH7PiYL8M0OWIArry7zcqsHCAo16x1rXLC3B0lK3vLCQj4Y9MtMaOKYjXcfyJgM3ADc4c5tl4CqXJvQ+Efk4cD/QAt6STpITEjYLNrdv8jB5k19XQH8tcO0wTCUkHOvwAWIZG2O5ZWjU6zzxBJx2ktPZOde9Rj2fRqBRd+54rVb3VNmdRA+Pzb0YJg+UhIR1jGa7ze9uEQ7Xt3HaX7/PHoy4zwMPWN3gcn2C5fpE5pHSqTcyg2xvajPaA5SUHS8hIWEN4HWIp/7O27o3l5ayxPWNpcM0lg5ni2Ft/lCPac3okA5QEhIS1hDNdpv/93HJTo7Zu5cPftAFh3X3OuMTVhweH8+7441cTD5GD1ASEhLWB3TE7F8/Ynhb/atQ3wOzswDc/sA2Lr7Y0i63nD4Ruj7MI8HGXOiqIC2GCQkJFbG5D1DSYpiQsIHQTUMq/BaP8L3v1di2axcAe+pWYm60WjTqy2Q/71Yr75EyFDamPrAK0mKYkLAB4V33ti19F+YWALj5znN47ZWd7LTZe6Q0Qu+VFaPDRj0proK0GCYkbFB4HeLbvmdd+1/7ow/TYSc1YHGplp2hdKiNbmOYxOSEhIT1CLsg2vgov8WTPPUUTCwtMT410T00WVrKny6vGElnmJCQkOCQdIYJCQnrFNp1b6J+BLDWNtPTTmc4Pj4ineHm3hkmo+uEhE2CZrtNc8sWHpzbxhkLD2b22cut2ggD12xeo+u0GCYkbCI0223+6jmCPOc4JlhkgkUarUW21kfhkudPkzenb3ISkxMSNhm82c1X7rWnzKefDl/84qha35i7vipIi2FCwiaEPmV+NfDikbQ6UHa8DYckJickJAyA1dcZishPu3zsHREpTDcqIpe4/Oz7ROSd6v6ZLlndQyLy1yLSqNJvWgwTEjYpfBrSvwV2/NZvjaDFoxa1Zi/wr4AvFBGIyBjwR8CrgPOA17i87QC/B7zXGHM28ATwxiqdJjE5IWGTQ0e7GQ4G+H5fqqF7MeYbADbtcSFeCOwzxux3tDcBl4vIN7D5mV7r6D4CNIEP9Ot3XSyGd91115yMjT2FTSp1rGCaY2u8cOyNeb2N9weGe/zJT8EnpisSj4vInap8vUsPPCrEcrRfBJwMzBtjWup+NHd7iHWxGBpjThGRO40xhfqBzYZjbbxw7I15s43XGHPJqNoSkU8D2yNVv2GM+bsqTUTumZL7fbEuFsOEhIRjC8aYlw/ZRFGO9jlgSkTqbndYOXd7OkBJSEjYiLgDONudHDeAK4FbXMrizwFXOLqrgCo7zXW1GI5Sn7ARcKyNF469MR9r4x0JRORfuhzt/xz4pIh8yt0/Q0RuBXC7vrcCnwK+AXzcGHOfa+LXgHeIyD6sDvFPK/VrF9KEhISEYxvraWeYkJCQsGZIi2FCQkIC62AxLHKp2WwQkQMicq+I3O3tr0Rkm4jc5tyGbhORk9aaz5VCRG4QkYMislfdi45PLN7nvvN7ROR5a8f5ylEw5qaIfNd9z3eLyE+qune5MX9TRH5ibbhOKMKaLoZ9XGo2I37cGHOhsj17J/AZ5zb0GVfeqLgRCO3Qisb3KuBs97maCt4B6xQ30jtmsK5gF7rPrQDuvb4S+CH3zB+79z9hnWCtd4aZS40xZhm4Cbh8jXk6mrgc6y6E+/vqNeRlKBhjvgAcCm4Xje9y4KPG4nasXdjpR4fT0aFgzEW4HLjJGHPEGPMtYB/2/U9YJ1jrxTDmUlPJdWYDwgB/LyJ3icjV7t5pxphHAdzfU9eMu9VB0fg2+/f+Vif+36BUH5t9zBsea70Yrth1ZgPiRcaY52FFxLeIyGhCzG1MbObv/QPAs4ELgUeB/+jub+Yxbwqs9WJY5FKz6WCMecT9PQj8DVZEesyLh+7vwbXjcFVQNL5N+70bYx4zxrSNMR3gQ3RF4U075s2CtV4Moy41a8zTyCEiJ4jIM/w18EpszLZbsO5CMIDb0AZC0fhuAX7OnSpfDDzpxemNjkD3+S+x3zPYMV8pIltE5Ezs4dFXjzZ/CcVY00ANxpiWiHiXmjHgBuVSs5lwGvA3Lj5bHfgrY/7/9u7YJAIwBsPwkx0cxBGutncAucLCIWzdRLgBxB2sz1ocwkr4Lc5SK8VTeZ8JkuaDEELW/cw8YDczWzzj/Ig1fsnM3GKDk/dTqmvc+Li/O5w5LBFecPHjBX+DT3rezMypwwj8hEtYa+1nZodHvOJqrfV/H4r8QZ3jJYnjj8lJ8isUhkmiMEwSFIZJgsIwSVAYJgkKwyQBb1OGeeSNOEtDAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "corr = np.zeros([len(covariance),len(covariance)])\n", + "for i in range(len(covariance)):\n", + " for j in range(len(covariance)):\n", + " corr[i, j]=covariance[i, j]/covariance[i, i]**(0.5)/covariance[j, j]**(0.5)\n", + "plt.imshow(corr, cmap='seismic',vmin=-1.0, vmax=1.0)\n", + "plt.colorbar()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sampling and Reconstruction" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The covariance module also has the ability to sample a new set of parameters using the covariance matrix. Currently the sampling uses numpy.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." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/icmeyer/openmc/openmc/data/resonance_covariance.py:233: UserWarning: Sampling routine does not guarantee positive values for parameters. This can lead to undefined behavior in the reconstruction routine.\n", + " warnings.warn(warn_str)\n" + ] + }, + { + "data": { + "text/plain": [ + "openmc.data.resonance.ReichMoore" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rm_resonance = gd157_endf.resonances.ranges[0]\n", + "n_samples = 5\n", + "samples = gd157_endf.resonance_covariance.ranges[0].sample(n_samples)\n", + "type(samples[0])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The sampling routine requires the incorporation of the `openmc.data.ResonanceRange` for the same resonance range object. This allows each sample itself to be its own `openmc.data.ResonanceRange` with a new set of parameters. Looking at some of the sampled parameters below:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Sample 1\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
energyLJneutronWidthcaptureWidthfissionWidthAfissionWidthB
00.03115102.00.0004710.1070450.00.0
12.82092102.00.0003340.0988850.00.0
216.21740801.00.0004530.0683850.00.0
316.77102102.00.0136700.0712780.00.0
420.55968502.00.0106090.0975460.00.0
\n", + "
" + ], + "text/plain": [ + " energy L J neutronWidth captureWidth fissionWidthA fissionWidthB\n", + "0 0.031151 0 2.0 0.000471 0.107045 0.0 0.0\n", + "1 2.820921 0 2.0 0.000334 0.098885 0.0 0.0\n", + "2 16.217408 0 1.0 0.000453 0.068385 0.0 0.0\n", + "3 16.771021 0 2.0 0.013670 0.071278 0.0 0.0\n", + "4 20.559685 0 2.0 0.010609 0.097546 0.0 0.0" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "print('Sample 1')\n", + "samples[0].parameters[:5]" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Sample 2\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
energyLJneutronWidthcaptureWidthfissionWidthAfissionWidthB
00.03383802.00.0004800.1033250.00.0
12.82237002.00.0003670.0915990.00.0
216.24396801.00.0003110.0896550.00.0
316.77599302.00.0130500.0844760.00.0
420.56169002.00.0111630.0868020.00.0
\n", + "
" + ], + "text/plain": [ + " energy L J neutronWidth captureWidth fissionWidthA fissionWidthB\n", + "0 0.033838 0 2.0 0.000480 0.103325 0.0 0.0\n", + "1 2.822370 0 2.0 0.000367 0.091599 0.0 0.0\n", + "2 16.243968 0 1.0 0.000311 0.089655 0.0 0.0\n", + "3 16.775993 0 2.0 0.013050 0.084476 0.0 0.0\n", + "4 20.561690 0 2.0 0.011163 0.086802 0.0 0.0" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "print('Sample 2')\n", + "samples[1].parameters[:5]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can reconstruct the cross section from the sampled parameters using the reconstruct method of `openmc.data.ResonanceRange`. For more on reconstruction see the Nuclear Data example notebook. " + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ]" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "gd157_endf.resonances.ranges" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0,0.5,'Cross section (b)')" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XecXHX1+P/XmZntm+ym95CEEJLQAoTQkaaCEKqi2EBBREXxJxb8flDko4iNjxgLCAJiAykqhCIgxUiLSYCEQAjpJKT3bJ+Ze35/zMzu7OzMTr1zZ3bP8/GYx+7cenazuWfeXVQVY4wxJpHP6wCMMcaUJksQxhhjkrIEYYwxJilLEMYYY5KyBGGMMSYpSxDGGGOSsgRhjDEmKUsQxhhjkrIEYYwxJilLEMYYY5IKeB1APoYOHaoTJkzwOgxjjCkrixYt2q6qw9IdV5YJQkRmA7MnT57MwoULvQ7HGGPKioisy+S4sqxiUtW5qnpFQ0OD16EYY0yfVZYJwhhjjPssQRhjjEnKEoQxxpikyjJBiMhsEbl9z549XodijDF9VlkmCGukNsYY95VlgjDGmHIVbmqiY8MGr8PISFmOgzDGmHL11sVXsGd7Kye8/HevQ0nLEoQxxhTRvDGXwxg4wetAMlCWVUzWSG2MMe4rywRhjdTGGOO+skwQxhhj3GcJwhhjTFKWIIwxxiRlCcIYY0xSliCMMcYkVZYJwrq5GmOM+8oyQVg3V2OMcV9ZJghjjDHuswRhjDEmKUsQxhhjkrIEYYzp81Zva+KHjy9DVb0OpaxYgjDG9HmX3bOQ2+et5t2dLV6HUlYsQRhj+ryQ4wAgiMeRlBdbD8IY0+eVQs1Sx4b3CNSUQCBZKMsShA2UM8bkQjwqQDgdQf5w3Su8cvk3vQkgR2WZIGygnDEmG16XIMJtrbRXD+L1YZd4G0iWyjJBGGNMeSmvqqUYSxDGmD7vyJ2v8btlP0ajjdUmM5YgjDF93hefvIMxy7fh3/CG16Fkbdnra/jH7U97cm/rxWSM6fuiNTwi3nwmdjprmLq3kjttbYR376Zi5MiU5z7/m7dxfFXuBdcLK0EYY/oPj1qrQymqttZ/+qOsPPmUXs/1KjmAJQhjTD8Q9lfQXDvSu2Fy0cSkvu6VNi1L3vEimoxZgjDG9HlvTr6E+bO+QyhY2Ebq0I4d7H7wwYJes5RYgjDG9Hm7B04BwAnnf61dm5tZs3gbABuuvppN132HjvXr879wCbIEYYzpB7Tbl3z85XvzefzWSG+o8PYdkcsGQ72ek6pxXBEcKd2+QqUbmTHGFFo0QbS+sRT/wAFU7rdfXpfbWzGMVdPexwQnTeZJMcfHW1M/zZaRszgoryjcYwnCGNNvxNaDWPuRjwAw7e1leV3v9eHn0FIxmL27gtT0fuOkm7eMnJXX/d1WMlVMIjJNRG4TkQdF5Atex2OM6UNivYjSfdLPUjh6vb1twd5vn6RuSzs6ChqLG1xNECJyl4hsFZGlCdvPEJHlIrJSRK4FUNVlqnolcBEw0824jMlZOISzYz2hd+YTfncJ2rTN64hMBiT6gC70KIhYiaS5vfc2iGSaX3mlwNEUnttVTL8HfgX8IbZBRPzAr4H3AxuABSLyiKq+JSLnANdGzzHGc1sWvMLqf/yL7e8Je53htPkaaa9sQMWHT3dT3fYmDcH1jG3YyMEfOYraUz/u3ZzSJq1CLzkaIEwHEGjaBMxIeZxXI7jz5WqCUNV5IjIhYfMsYKWqrgYQkfuAc4G3VPUR4BEReQz4i5uxGZNIVVm1fC1r7v0Hzav3sk/HsK9uP1SOw1cZpL51PVXhDQQ61tAuftQRlBq2105hEzN57S9NTLnrGo772plUHfl+r38c000kMbg1V1+yKqS+wItG6jFAfKfhDcDRInIycAFQBTye6mQRuQK4AmD8+PHuRWn6vFDY4c0Nu1n2yL9h8XJCbY3sq5mA4z8MqQ0zoH09o/xLGDBjIvt96FQmjHo/FX5ft/NXbG1i3jtbWfz0m8zcsI23BpzNu7ds4P0zr2P0V75vpYkS0fmvEC1BhPxViBZgUETnddPtT3JAGfxteJEgkv1WVFWfB55Pd7Kq3g7cDjBz5sy+mbaNK1o6Qry2bheLX1yBb+FiavZBa+V4wv4GCMxigH89Y/R1RkwbxqGfOI+aob2XAgJ+H9NGDWTaqIE4J+7PE0s389jdz3FYaBBzl8zi1Ouu5oAf/KIsHgR9X/c2iHkn/h/Vrds5uEDX7av/xF4kiA3AuLj3Y4GN2VxARGYDsydPnlzIuEwf0tIR4q2Ne1mybhfvvr6OxreXM7DFR0fFGDRQS5gDCbOV0a2vMXq8w5QLTmXgkZ/O+X+6zyecdegojr/xw3znjpc47PWtPLPlg/i+93X2v+HmAv90JmuxcXJxdUxtNUMLcOHI30u6po1yTSBeJIgFwAEiMhF4D/gY8PFsLqCqc4G5M2fO/JwL8ZkyoarsaO5g7fZmVm5tYuV7e9m+ehsNa9YwoqmNKmmgrXI4w8QPHIjoFkY3v8rIIXuZeMx+DJ39UaThYwWNqbG2kp9fdRLX3/NfJr6wnmfWn8SA3/2Y4Zd/q6D3MeUleffa0s8ariYIEbkXOBkYKiIbgOtV9U4RuQp4EvADd6nqm27GYcpPWzDMjuYOdjZ1sL25ne372tm8p42NO1rYubkJZ9MOGnftYHR7G/VahfoHMaiigUEATMDvb6KuZR3j2t9gxNB2xhw9gSFnno8Mvtj12AN+H/976dFc19zKfkta+OfzI7joqHlUH3aS6/c2vdOCN1IXrpa7fc0adt//AMO/+Q2kRIocbvdiSvq/UVUfp5eG6HSsiqm0qCrtISf6CtMedOgIO7QHo+9DDm3BME3tIZraQjS1h9gX/drUFqKpLUhrc4j2fR2wu4mKvU0MbG1hiNPBQMehRv34pYZK/0DGBuoYC0A9UI/f30Z162YGtLxFQ8VOBg/qYPjUoQw96Xgqpl4GlbWe/E78PuF7XziJH3z3r4zcOZ5nb3qaM+85Eqmq8yQeowlfCyWzB3kmvWs3fPFLdKxZQ+NFH6Fq4kRefeApRh82Jc/48lOWU23kW8X0nxXbeGHl9h7zd8X6SGuP7bH33feTeF6a4xP3k7g/w/M0IfAe8fc4vmu/o0rYUcIOhMNhHAccx8EJK044sk9j2xzQsIPjKGh01KgDTjhMIBTC1xHEFwzjD4aochwq1aHKcajREFXqUKUOFapUqVKlEMCHHz8iFdRJBbW+Sob6qnF8lXG/zarIyz8EfA4VwWYqQ/uoCW6n1nmH+sBeBtS2MXCwn4bxjQw6/FAqpp+GNIwDX2n1Na8K+Lnym+fzj6/9njUN7+Pt713PtJt+5nVY/ZpT4JHUmUvfi0mjU83GSg8vPxOg5rH5UD3M9ehSKcsEka/37nyEIVsTP8kl+yQgqfdLiu1p30dmcBTi/mR6FCcze5/8/F6OFR+KRL6Kv0dcmfNHX0lWukp4RvvCHficIIFQK4FwG/5wGxXOHgLaSoW2UiHtVEo7lRUhaqqD1NZBXWMFtcPqqRszlIpR4wiMmYYMmQg1g8qutW9MYw2Tv3Iu63/xAi9vPIqJb75E9UHHeR1WvyMkfAorGLeuG9HqYXKAMk0Q+VYxTQmtZVVTY/RdtJtajwwf2971ffdtdP1RiMZ/SXvN7vtSXDvh3B5/4JK4PbJN0BTHCIKDoAhhfKKR96L4ol9FFJ84iBDZ5tPI9xL96ose44NABQQqBH+lj0CFEKgKEKjyE6gJEKiupKKmCn9tNb6aWnx1A/ENGoavcSJS2wBVAyKvyrqye+Dn4oMzRvPT/X3UbhjMCzf9g9P/eGy/+LlLkaPe/N4fuvktT+6br7JMEPlWMR1x4zc5oqM5+i72FJUk7xP30cux2V4rj/f2cCk7n/3qbB6+6m5W1JzCYQ/cxbCLLvM6pP6pkIPj4qX5P7lvZ5KJ+XatAUa7E0+BlGWCyFvsE6wxRTKkvoqqC4+mfe4WXvj7Ns67IIgEKrwOq/9xaaqNnKqYWneRLEEUer6ofJRWq54xfdhHP3QYleHlbBw4k3d/+3Ovw+mXtPD9XHOWbVfWcCj7GWPzVZYJQkRmi8jte/bs8ToUYzIW8PuY9IXz8DvtLHjBj4Z6X0PAFFJssj63ZuvL/lN/qlNKZQwElGmCUNW5qnpFQ0OD16EYk5VTZu5HVWA1WxoOZ/VvbvE6nP6ndGpvaGkv/Rr+skwQxpSz6VddSCDUzKL5FWi4+NUG/VIsMYRdKkHk8Kl/277uXc3bgpHYOsKlk8UsQRhTZEcfNJoa/0q2DTiUdXf+xutw+hXX2iByqGJqau9exbgxWMfKSeezeXdL8lu4ldx6UZYJwtogTLmb+oUL8IU7eO35JtcGWZku7g2Ui142h3NC4e5VTOv2/zjvjj+dlr0udcXNQVkmCGuDMOVu1oz9qNflbKqfyaaH7/U6nH7DtUZq4KWV27Nam/q97d2nG++ojr53ek8Qyx66kXULcp7KLitlmSCM6QvGf+o0EGHRQ+94HUo/kDAxWYHtaQ3y8d/N55r7F+d0/r6dbZ3f+1u3Jj0mltzevl/YfMufcrpPtixBGOORk046iIa2ZWyomsW+xQu8DqdfSL4uQ/46QpGH94qt+3I6v3VfkpHWKWwcdAwLB1ye032yZQnCGI+ICINOm0I4UM38Xz/sdTj9g0tVTD3nXUtP4yfWjOvN5lYSy0VZJghrpDZ9xQc+8QEGtK7l3Y6DCe/a4nU4fZ7j2kjqHAa3hePaGnau6rpSii6zXiSOskwQ1kht+oqA30f9ZKG1ejj//cmtXofT97n0kM3pqlu2pbhYiqt5ME1IWSYIY/qS91/zcSo7drPu3aGoi71sjIsT4akyqy1ATTDz62uoq92h1yVdPJQ2QYjIsSLyaxFZIiLbRORdEXlcRL4kIvYR3pg8DaitorF2HTsGTGfFH+/xOpy+LY8Eoaq0r1iRdF97q5/3tVVwzKZsrt+VCcKtbb0c551eE4SIPAFcDjwJnAGMAqYD1wHVwMMico7bQRrT18248jzECfLWvzZ6HUoflf9AuV1/+AOrZyd/3MUu68+qANiVINZd88Oua6UoQmia8RFuSDdb1KdUdXvCtibg1ejrZhEZ2vM0Y0w2Djh4IovaH2Zz1eHsfnspjVMP9jqkPimfKrzWN5am3JdLrVC3VBXfYF1CI+t7LUHEJwcRGSki50R7EI1MdowxJncjT51EOFDNgjkPeh1K3+XWVBudJYien/KzTUqb3mknHCqNtqiMGqlF5HLgv8AFwIeBV0Tks24GliYe6+Zq+pwTP3k29a3rea9pf5yOdq/D6ZPc6irqiz6LGpt399i36X+uS3t+fAnknZUDuO2q53sc45RwN9dvAIer6qWqeglwJPAt98LqnXVzNX2R3+9j0IjdNNeOYdGc27wOp2/K4xm7tyV1Q7I/GNnnp+dcTHv+/vfcb+qxTBPEBiB+DPk+YH3hwzGmfzvh65/GH2pl3WtWgnBDPtN9t61/PdkFAfA5rZGvSTJQ6pwkGRzjrV4bqUXka9Fv3wPmi8jDRH6Wc4lUORljCmjw8EEMcZazre4w1r84j3HHn+R1SH1KPuMgAppkvqToAIbeGqnfPvATyWMpoaVFU0lXghgQfa0C/kFXonsY2ORiXMb0W1M+fBSOr4I37v6X16H0GZ2P4nwaqXs7tZfrbhp1XIo92SWIkuvmqqo3FCsQY0zEYR86kSV/vZMtTKNj724qBzZ6HVKfUfBOTNEL5nvZ9qpBWZ/TvGcvdQ0D87xz79INlLtdRJJ2yBaROhH5rIgkLz8ZY3I2YlKIlpoRvHyzNVYXVD7zGfX2gT+n6qKuc9445PNZn922q2ePqUJLN1DuN8B3ReQQYCmwjcgI6gOAgcBdwJ9djdCYfuikr1/CmqufYsvKCq9D6Vvc6ipahMFtYQ/WpE5XxfQ6cJGI1AMziUy10QosU9XlRYjPmH6puq6aYf6VbK47lGVz5zJt9myvQ+oTNJ/KoCxOfeffq3j70cXM/vF5ud+vBGTUzVVVm1T1eVW9V1X/YcnBGPfNuORU1Odn5YPWYbBQ8qlhCkkl68adnmJv9yqmp+9dx/p9jbTMn586ljynbXVzfe2Yspzu20ZSm/5g0nEzaGxdyVbfoezeYp0G8xKtAhJyf6iuaTyTVfufn3xnimd9Qadv92Aq+LJMEDaS2vQX4w6tpq1qCAt+ervXofQJ+ZUgqnu5cC5VV+U/DsIY46HjvnwxlcE97Ng4yBYTKoRir8pWqkOkM5TpZH1TROQOEXlKRJ6NvdwOzpj+LlBZwYiatewYMJ0Ff7rP63DKnnudjYpfGijGB4Z03VxjHgBuA+4Aij+cz5h+7Ogvzmb9L9ay8el34NNeR1Pe8prNNUkOyCctZNtIreHiF0cyTRAhVbUV1Y3xwIjpkxnS/iRbKw7n3XdWMn7KZK9DKmNuTbWRw/WyzC6OFv+zeaZtEHNF5IsiMkpEBsderkZmjOk05YThBCsH8MacP3odSnkrYB1TZOK/6FQbqdo2er1ftiWIEhsoF+eS6NdvxG1TYFJhwzHGJDPjMxfy+ov3s6tpLG0dQaorbYR1LvKZzbXHtRwH6Sw6ZF/ZlHUVU6l2c1XViUlelhyMKRKf38fYIZvZU78/L/zmLq/DKV/55IeE57mGwukvWMApvR0PShCZ9mKqEJGviMiD0ddVImIfYYwpouO++lF84Q52L9zhdSjlK59G6oRTu41kzumyWZYgit1Fl8zbIG4lsszob6KvI6PbjDFFUj92FMNDy9heM4MlLy/yOpyyVNAqprhP9CmvW8j7hboniGJUOWWaII5S1UtU9dno6zPAUW4GZozp6dBzphL2V7P6rrleh1JmCt9FVDWuiil2eRd7opZyCSIsIvvH3ojIJFwYDyEi50UH5D0sIh8o9PWNKXcHnH8mA1vXsavjAHbsafY6nLJTyGdsuNsn+mhmSGynKGAvJjxYUS7TBPEN4DkReV5E/g08C1yTyYkicpeIbBWRpQnbzxCR5SKyUkSuBYjOFPs54FLgoxn/FMb0IxPG76OlZhQv3vI7r0MpH53P6cJlCMcJIbEEoMkf9lu2pk4QKU5Jfb9QiZYgVPUZIosEfSX6OlBVn8vwHr8HzojfICJ+4NfAmcB04GIRmR53yHXR/caYBEdfcwmBYDMtK8KE3VoAp68qYJuA060NIsUxvf77ZJchnMTiTxH+7dMtOXpq9OsFwFnAZGB/4KzotrRUdR6wM2HzLGClqq5W1Q7gPuBcifgx8ISqvprdj2JM/1DZ2MBoWcbO2kOY92imn9MMgGb7sb034a6Bcj3aIjKSYSzRxKCOw8a312Zzg7ylK0G8L/p1dpLX2XncdwywPu79hui2LwOnAx8WkSuTnSgiV4jIQhFZuG3btjxCMKZ8zfzUcajPz46/v+B1KOUlnxJEwvO8xyf6JMcU8jO+Osq2N5YU8IrppVty9Prot/+rqmvi94nIxDzumyx1qqrOAeakiel24HaAmTNnWvna9Euj3ncCQ+78DbsCB7Ni7WYOmDDS65DKQ0G7nYY6R1JLDtfN9AyJLZTqOEWfsC/TRuqHkmx7MI/7bgDGxb0fC2zM9GRbUc4YmHaEn47KRhbf8gevQykjhZxqo6uKKZfm44wru2K1V6072fD8WzncKXfp2iCmisiFQIOIXBD3uhToZXmltBYAB4jIRBGpBD4GPJLpybainDFwyFWfpa51E807h9PU1uF1OCUt9kk/n26uiaWEbrOrxto2EvNPKNjLFTNNVpHjtv35btZxdIbnFEa6EsSBRNoaGune/nAE8LlMbiAi9wIvAweKyAYRuUxVQ8BVwJPAMuB+VX0ztx/BmP7JF6hgv6Hraaodz/O//JPX4ZSFfEZSa8K8SprBOAh2v5vz/RLN911WsGtlKl0bxMPAwyJyrKq+nMsNVPXiFNsfBx7P5ZoiMhuYPXmyzYtv+rdjv3kJK/7fAvYtaUJVkQJODtc3FbCbq+OQtvdSb/8cGYeS/EAtwsC5TNsgrhSRxtgbERkkIp5NKWlVTMZEVI8YxWjeYGftdF6e+7zX4ZS8/EZSJzyo1el6/muKRNHrMIjMkrmXKT/TBHGoqu6OvVHVXcDh7oRkjMnGrM+diqjy3t/mex1Kyctvsr6EKqZw1yd4TX5I77FkfKB3nTUzTRA+ERkUexNdTS7TxYYKznoxGdNl+DHHMrxtCTsCh/DW0lVeh1PSCvms1bCmvWDTwtRNq+VQGZhpgrgZeElEvi8i/wu8BPzEvbB6Z1VMxnR3yGmjCAdqeOtXf/U6lBIVfZAXcHqKTNogWnsb+ZxntiqZ6b5V9Q/AhcAWYBtwgara4rjGlIgpn/oYDc2r2NM+mfe27vU6nJIjndNVaMHWhNC4NojOSyZcevGhX+olqEzLEKVfxQQwGGhW1V8C2/IcSW2MKSDx+Zg6tZm2qqG8/GOb5bWn2DgIdybrS9nNNW1E6UmpJwgRuR74FvDt6KYKwLOO19YGYUxPR1zzeepaN9K0fTi7W9q9DqekSGyss6PdVoLL8iLdRLqZ5jJJX4oLppLi2uLL5vN9bjK9w/nAOUAzgKpuBAa4FVQ61gZhTE++6homj1xPc81onvuJZ73QS1JsFLRqwlrSeeg+L1IuGSK/KqaSaYMAOjRScacAIlLnXkjGmFwd/f+upLp9B02rK2gPFX8FstIVq2JyCrasXGQJ0K7EE3ebzJRBN6ZME8T9IvJboFFEPgf8C7jDvbCMMbmoaBjEhPq32Vs7iSd/+RevwykZnfX4TuIUGblTR+MGynXeqE/JtBfTz4jM3voQkfmZvhttrPaEtUEYk9px115CRXAf+17dSzDX+va+JlZq0Ngn/0JcMq4NIpfzs+ojlOT8Ikz9nWkjdR3wrKp+g0jJoUZEKlyNrBfWBmFMajVjxrOfbwm766bx5N0ZT5Lcp8V3R43v5qqqNLWHmHDtY/xz6aasrtmtu2wOz2oV9xuZ85VphPOAKhEZQ6R66TNE1po2xpSg4685H3+olb3/XkfIShF0rtig4MRNkYHjsHZ7MwBznlmZ1RXjP8Frqum+ezs/wwRR8t1cAVHVFuAC4Jeqej4w3b2wjDH5qJ96MPvpq+yqO4R/3vOY1+F4Lr4XU7cRzMnmU8pQ9yqmHMZBiD/LOxZfxglCRI4FPgHE/to8m4vJGJPeidfMxh9qZc9zawgXcIBYWZKuBYPiu7k6wVDmA5oTJ+tzNK9P9/m2QSRdE7vAMo3waiKD5P6uqm+KyCTgOffCMsbkq376DPZjEbvrDubJe6wtIiZ+NHU42LUSX9opOHoMlHN6FCCyiqOvtEGo6jxVPUdVfxx9v1pVv+JuaKlZLyZjMnPCNecTCLaw67l3+3ePpvgHeVwiCLe2I8EQP3zxt4zbujbLa3b9PrXHNxmc3lcSRKmxXkzGZGbAtEOYIK+yu/YgHv/d370Ox0ORj/+q0m0ltlB7G7p6FaFx53HOgmezumKkkTo21YaLbRApSja+IiSYskwQxpjMnfCtDxMINrPvxc20dvTP0dWd60krOHFVTMG2Dpp3d7Bv4AS2jn5/r9dIfPY7xE33ncMIub7UBmGMKVN1B0xnUuBV9tROY+6c+7wOxyPREgQ+iKtqC7a1EWjeDECAdMkz4ZN8spHUWegzVUwi8hMRGSgiFSLyjIhsF5FPuh2cMaYwTvrOJVS176T1jSC79vXHmV5jVUwQ/zQPdbST63M6frI8zSFDZJ4gSn8cxAdUdS9wNrABmAJ8w7WojDEFVTV2AlMHvUVTzXieuulOr8Mpus7P+irdezG1ddD1GMyumqjbVBeaQxVTpgkiRVtFoRY+6k2mCSI2rcaHgHtVdadL8WTEejEZk71j//dq6lvW07RpKGs37fI6nOKSWBVT97r7UEcHkuO029qtDSKHh3WGCSIUqEkRTukkiLki8jYwE3hGRIYBbe6F1TvrxWRM9vz1DRw6dSttVUN56Qf9qxTR1Vog3SbrC3cEOwsO2Y+kzm8upnKQ6TiIa4FjgZmqGiSycNC5bgZmjCm8Gd/6GkOalrGvbQovLVjhdThF1NXNFSc+QXTEfZDvvSSRuDd+JHUfzQ8ZN1J/BAipalhEriOy3OhoVyMzxhSc+P0cc0Y94UANq389tx9NwdHViym+7SAcDJF520PC70qdzpKJdO4q3oIQhVxfO5VMq5i+o6r7ROQE4IPAPcCt7oVljHHLhE9ewriWV9hTdSiP3v2o1+EURec4CAQnfiR1MAhx+7K6pgOJbRDax1YMyjRBxDoInwXcqqoPA5XuhGSMcdup155FZXAfu+dtY1dT/+n2qirdpsgIdwQzfqQnTsynTrjbOhORg4qXIMTn/r0yTRDvRZccvQh4XESqsjjXGFNi6g46gmkDX6epZgKPf+92r8NxXeyTvUOgWxWT0xFC/L7oMemukfA+HO6aJTbhPn1Fpg/5i4AngTNUdTcwGBsHYUxZO+6HX6exaRXNu8Yzf1Ffb7COJQFft/ED4Y6u2VzTfvpP2O10BLsaqZ3cqqnyoSH3p03JtBdTC7AK+KCIXAUMV9WnXI3MGOMqX+0Ajj3VR8hfy4o5j9MR6vuzvaoEundz3b2PTB/qiUdFqqci14qtKFfMEkSoI+j6PTLtxXQ18GdgePT1JxH5spuBGWPcN+mzlzGp4z/sqTmEv/30z+7e7IlvwTPfd/ceKcQaqR38ENf7x3n0+SwaqbtXMoWDYUS6J4hitkF0K/24JNMqpsuAo1X1u6r6XeAY4HPuhdU7G0ltTOGceuNl1LVspGlFLW+t2OTejebfBv/5mXvX71W0ikkC3bqHtlc14M+wDSKREwrFlSBym64jH+FSKUEQ+anjK7zCFPM3kcBGUhtTOJVj9mPWETtor2jg1R/e1ycXFuosQUigc31qgGXTLulspM56oFww3DUK24MqplIqQdwNzBeR74nI94BXgP41Vt+YPmz6V69mQsuL7Kk6jL/+5C9eh+MaB3+3NakBfIHIVHPpJ89LrGKKK0FQ/EbqkilBqOr/AZ8BdgK7gM+o6i1uBmaMKa7Tb/wE9S3v0bKynv+6MA2iPtZDAAAblElEQVTHyrZjeaP5zKJMMpco9gBXqehx/65urtn13HfCisS6uWqsCquYCSLk+j3S/kZExCciS1X1VVWdo6q/UNXXXI/MGFNUVeMnc+L7Wgj7qln+q2doaivsJ9Qnd3+TefuugI6mgl43M9EqJl+g2zoO8fuyXhgi1FXTrmS4fGgBOcESqGLSSJ+wxSIy3vVojDGemnTZ55nKPPbWTOFv37rVnTUH2vYW/pppxLdBOAlzGMXeZtt+4ISduAFyxW+k3jjvHdfvkWnKHAW8GV1N7pHYy83AjDHeeN8t1zJs7xs0tR7Ig797ouDXDzfvLvg10+tKEKFwwgCzaBJM1wbRY6qNkENXCSJWxVS8CSaC4Yr0B+UpkOFxN7gahTGmZPjqBvLBr8zgodu2sPtlP4uOXMWRR+xfsOvv2bKNwUWeC7pzqg2fn46OhKVsNNbQnH0bRNckfbF2jCKOpFb3k1GvdxCRySJyvKr+O/5F5LeywfXojDGeaDjmfZwwczMhfx1L5zzP1t0tBbv2rnc3F+xaGYt+sldfBe0trUkPSfvpP+HZHxlP4V0bRDGa+tOloFuAfUm2t0T3GWP6qClXfYVDK5+nqXoiT3z9DtoLNPfPXo+XO23fl6KRPMseSE7cpH+dyaWIvZhyWQc7W+kSxARVXZK4UVUXAhNcicgYUzKO+8X32W/vPJoCh3DvN+8oSKN1y/bi92KKr/pJTBAabZNIV8WU2AYRP2VHWCozukZheZ8gqnvZl2IlbWNMXyGBCs6Y8zmG7V1Mc/Nk/vqT+3K+li8cWXeifZ/73TN7kK51IIJ7Igmiqj1SklmycDGQ/RgGdbSznifsq+q6T5FoCZQgFohIjzmXROQyYFEhAxGRSSJyp4g8WMjrGmPyExg6htnfOYmBzevYtXIQD9+T20TOgVCk7r+9ubhTeagqihAIRxqnnaZIHJUaqT3fvumIyIFZ9kBy4oaJxBJEMYcA+vzeD5T7KvAZEXleRG6Ovv4NXA5cne7iInKXiGwVkaUJ288QkeUislJErgVQ1dWqelmuP4gxxj01047kjEuHU9Wxhy3/aePJf7yS9TX8TqTkEGwv7jRukVoxIeBEEkSoOVKSCWhz9+MkXUNzwvgJRzpLDGF/rLKleFVMJ9/4adfv0etPo6pbVPU4It1c10ZfN6jqsaqaSVeE3wNnxG8QET/wa+BMYDpwsYhMzzpyY0xRDTt9Nqe/vxlfOMS6uZv599PZTagQa8jt6KhyI7zU9yVSfeR3IonBaY988g74slxqNaH6yAn7O1OG44vN51S85DdkzCjX75HpXEzPqeovo69nM724qs4jMn9TvFnAymiJoQO4Dzg344iNMZ4Z/8nPcuoxGxBV3rlvDS/PezOzEx2nM0G0McjFCHuKNKwLfo0miI5IFZe/Irsqmti8SzGO+pHEOqXoz1jduj2nWEuNF+tKjwHWx73fAIwRkSEichtwuIh8O9XJInKFiCwUkYXbtm1zO1ZjTIJJn/8KJx22CiXAm79fxosvvJ3+JA13VuG0BIZDqHgN1ZEaJiFA5J6xtgN/dXaPv8TSgUMgZZtDtRZ/OhE3eJEgkpXBVFV3qOqVqrq/qt6U6mRVvV1VZ6rqzGHDhrkYpjEmlQOv/jonHfgmYali2V1v8J8Xl/d6vBMOdSaIYEU96+e/UIwwgUgbhCL4JRiNJRJHZXV21UGSkCDC9DbVRd9YU8OLBLEBGBf3fiyw0YM4jDF5mPrNb3PSAYsJSw3L71zMvJdSTx4XCgVR8TGgeQ0Aq578b7HCxFEF8eHzRaqUHCeSIHw96ofS8CWWICpJNRahmFNuuMmLBLEAOEBEJopIJfAxIKuJ/2zJUWNKw7Rrr+OkSYsI+2pZ8bvXmffyyqTHhYJBHPFTV7mbQKiFneuL9+iJdHMFfyxBaOoG5e3vJp9BqOWFZ/EljCQPS0XKKibxYM0LN7j6ryQi9wIvAweKyAYRuUxVQ8BVwJPAMuB+Vc2wpSvClhw1pnRM+5/rOXHcAoK+OlbcsYgX5q/qcUyoow0VH74KGNz6FlsCh7JjRY9JGlwRmzNJxMHnBHFiVUNJPuSvee7FpNd49Ru/wtncfQ6nkK+ml3KCJYi0VPViVR2lqhWqOlZV74xuf1xVp0TbG250MwZjjPumX38DJ455haCvnuW/XchLi9Z02++0t6PiRwT2O3E/HF+Al24uzooBqg4qgojgD7cRJtLNVpI83rcv6llNpqosPuwqNow9pdv2oL+O1NNdWILwjFUxGVN6Drrh+5ww8iWCvgEs+/V8Fi7u6qzY0dYG4kN8wlGfv4hhe15jvR7NgjvucD0uDcfWbVD84XZCEk0QIsz64rRuxzbt6DmDUM8V6CJCFXWoowRC3We6DYbC+KyR2jtWxWRMaTr4Bzdy/PD/EPQ3sPgXL7J+fWS+o3BbZAyC+CIP5pO+9SGqW7fz6vyRvHbfX12NSZ0wiCACAaeNkC+aBASOOrRrsJk/1Mr22oNZ9UpCA3qKBAEQ1ip8TvelWTe8vRYfhV2u1StlmSCMMaXrkB/exKwBT9Hhb+SZ6/+G4zgEowki1hFo5OEHcfSZ9QSCrcx/ZiDP3DSnc2W3QnPCoc5eRX7tIOTvXkqoCEbmZBrpfwvHX8mbtz3abb8mrkAXJyi1JFYnbXltCT5xf56kYijLBGFVTMaUtiN++lMObH2C5sqJPDbnMcLtkU/U4u+qsz/o4tmcdEEj9c3v8va6g3ngijm8tyyDQXdZ0nA4WnQBP+2dCSI2N19lKDKorW6IMKRtGRsDx7DsyX91nu8kyVuxUkPQX48AgVDXvE67lq1DClzFNGLLgoJeL1NlmSCsismYEucPcPwNn6C+eQNbFzfT0RoZxZw4YeoB536A8358NuN3Pct2DmTu/63i0WtvoXlX4RYVckLRT/MCfoI4/sjaDbFerpVObH0I4ejLZ+H4Arz5h1c72x6StUFUBCPndPjrI9cIda2r1rq5veBLNYw7yJsSSVkmCGNM6avafxYTaxbSVjWc5f9ZDXQvQcTUjxvD2fd9n2MPe4/GPStYt/tQ/vKNf/PEDbfSkWJ50Gw44cjDVegaC0HnFghIdI3qUJiJJxzFGN8StgyYyXM3/jyyPUnVV0U4kiBCFXWAUuF0lSBag415x5zIGTWm2/uRrcUpUZRlgrAqJmPKwyGXnoY/1MbeNZF6fJ8v+UdrEWHGVZfxkbsu59ABL1HXtJHVmw7kD195gn/+8HedbRi5iPRiis7o6o8rDUjs3pHYnGDk6wd+9Hmq27ezbuVIdm7aFGnkTlChcT2XVAkQeR8INtNUNcb1xXzqhrp6+U5lmSCsismY8jDo6LMY2rKMvdX7AyD+3h85/vp6TvzpdXzkNxdzUOBZqpu3s+rdSdxz1WM89ZN7CLZn3zso9oAXlECgqzQgvkgsPokkjXAosq+mcQBTD22lpWYUL153W2SqjgQ+DeEPdZVu/L5IFVpNcDuhQC3t2jeeTWWZIIwxZcIfYEDVFtQXnf8okNkn64rBgzj5Vz/gw784n+k8RXXzdlasHsc9X3qEp2/+C6GOzOvku9oghEBV1/0lWpqRaIJwwl2J4Pj/7xKGtixlg+9Ylj32zx7XVKAy1LW2ta8qmlyItJ3sqx7X45x4oql7RqVyxlcn0dC2qiuAIrAEYYxx1cDRXSu1+dOUIBJVDx/GKbf9iPP/7xymhx+junkH76wYye+/+AjPzHmAcEf6B21nggACVYG4PdEE4Ys8bTUh5xz3lZMRdXj7n4lL2hAdU9HV7rDf8fWgDgecGMYX7uhsCE8phwf8/lMnUO3szv7EPFiCMMa4aujBEzu/l0Buj5y6USM55Y6bOfenH2JaxyNUN2/n7beGcM+XHuKdZ1/u9VyNNVILVNbFTdEdLUHEZnV1EtoNxs08lHEVS9hVe2Cyq1JBpIpJUGZ+6gtceuUeZlz2Nerb38vpZ8yKlSBSs0ZqY8rH2GNndn6faRVTKgPGjuXUu27h7Js+wIH7HkI74Om/NvPYd+/GCScfe+BEq6PEp1TVdw2Si/WoipUgkg14OParH04ZS0DiGs5FqDv8QvD5qZFMuuhm94RP7GqrRcoQZZkgrJHamPJRM3Jy5/f+gL+XIzPXOGECp//518w6Qxi17SXWbt2P+6/+PcHWnivVdSUIoWZQfef2zgWAol+T9TwafOBkGvclX+ciUJm8equyKoOG9BzXrpZoYijWahNlmSCMMWWkeiC+cOTTdqA2Td18lg65+COcOudyJm36OzuC43nomj/1aJcItkXGOYgPaofEjVGQhBJEisduXWBH0u0V0QbvkK+62/bqQYX9GaGrx1Us0mItN2EJwhjjuqpgZDqLyrqes6Xmq3HsaE7584+ZtPlv7HAmMPe6e7rtD3fERnH7GDi8awBBZy+m6PtUD92qgckfk4HomtZhf/cEUTcidc2GRLvcZrviXGcVU5EXqrMEYYxx3bhwZOTvgFGDXbl+dX0tJ951I+M2/Yv39k5g8d+e69wXao/OJOv3UTMsSYKIPgVVkz8Oqwb3TGqC4q+Orkzn615tVjN8UMo4x+uzsZv3+vOM9i3sdX+xlGWCsEZqY8rLifs/wFnV17P/4RPTH5yj+iGD2P+Lp1PftIGFj2/vrGpyoiUIn1+oGjyw64RotU3sU3mqWpuK+urk2+uqkm4fMHpEyhh9NQNT7os3ZP8hSbf7JdYjqzjrTZRlgrBGamPKS/WMs5jQuASpKfw8RfEO+uCpjGABbYEhvHz7QwCEO6Izyfp8VAzs2Ujd2VidYnqMQE2SajFVKuprkx4/aL9eBslFn7gNTat7+zG6Yuq8X+TLidfOZlLrU5x03WW9nl8oZZkgjDFl5qyb4coXodadKqZ4h1x/DQP2rWPl62FUlXAwmiACPqSyqwG5RxVTigr+QG2SEoQKlQ11SY8fNHZU0u2Rmwkz5DamntHW+w+R4sk8dOohnHnPj6gdkrqUUkiWIIwx7quogZEHF+VWYw6YwJDQEpoDI9iwYGlnCcLnD3Q7rrNnUCCyPVUbRKCmZ6+ksFRQ25i8BsNf1UsvJoXjb72fmZ+8qtefocekhkVunO6Mw5vbGmOMexpOm4FomDce/A8ajNbbVySMwYgmiOCI0QCEKnqOoQCoSPLAVwLUNA7IOi7NtH9qiiqmYrMEYYzpcw7/xEU07F3J1h0NnSWIxEF6sSqm9pNP5O2dz/L6aacnvdaEcT0TgSMBKuuTVzH1RnLMD16xBGGM6XPq6mqo1/U0V4wi2BJtg6is6HZMLEEcd8Aw5k46lguP3T/ptapGjOXU57/UbZsjAaoGZJ8gMilAnPLJqZ2xea0sE4R1czXGpFM5LFJi2PdedBR3tIE6tp60+CP7xw6qZe2PzmLWxOQN6IFxExn+ja932+ZIgOoB9UmP700mVUzTTxhN/ejk3VyLrSwThHVzNcakM/iEGaAOe5oiz4mKmu4JAn9m80KJCEMu696t1JFAbqPC4/ODph7LcMh5ZzNj2lsMbFmX/T0KqCwThDHGpHPgGWdQ07qdvRVjga5pPsSJNFr7AoGU56bjSIBAdWZzLo1uf6nrTVyCkF4SBMDxV1+FaPYr6BWSJQhjTJ/U2FBHbccWHF/kQV4VrRLyRVdz81fkkSB8gYzPl7jG8fhpugc0b8j4fhn3fiowSxDGmD6rQrpWYKscEJnmIjZhHgmN1tlwfNmcG9fgHDdae9yoNzq//8KvT845FjdZgjDG9Fn++q5P3rWDI9N8+KJrizp5DC5QX++lh8NfvyXFiV33PPlnP2fimkcZtGs5viyXYi2W0ozKGGMKoGJi17QXtYOijdUaWUs6UJv9QLcYld4buAftXkHj7hWRY2u6ShuJKWnypsc4fPEcAI45bxJjpnSfq6rOH1mdrnFscabWSGQJwhjTZw0/4ZjO72uGRqqYTvvyUUysf4sZpx+f3cXSNCrHm/DA/Z2lhcYhcQ3NCZfY/+wtHHDeZgCOPGMC533tiG77PzTni5wwaxsHXXh2drEWiCUIY0yfNe2oaZ3f+6PVOKNmzeBDP7sKny+7x9+xZ1czZsR7GR1bc8ghEO3kFKgbzIS255Ie5//yvwlc/kDK61QNHMBhn/1oVnEWUlkmCBsoZ4zJRH1NJYN2PsewfS+lPziNI2afwHk3fArUYXDFirTHd64frQ6hQHTMRGId06jDYHLyKT5KQe79vDykqnOBuTNnzvyc17EYY0rbx+//fkG7iX7pt6cD6R/qB9ct59Wmeg6eOYV/L38F8GzOvZyVZQnCGGOy0WMBniI4+OYbOe/CATScfErX5HtlliHKsgRhjDGlzldXR+P553XbVmb5wRKEMcbkavpZgwk5GTz2rQRhjDH9yymzZ3gdgqusDcIYY1zWub5DmZUgLEEYY4zrtNuXcmEJwhhjXBbrRaWUxkpxmbIEYYwxbst0MeoSYwnCGGNcFy1BqJUgjDHGxBErQRhjjElKyvNRWzLjIESkDvgN0AE8r6p/9jgkY4wpCCm37ktRrqY1EblLRLaKyNKE7WeIyHIRWSki10Y3XwA8qKqfA85xMy5jjCmmilGRz+JVjc0eR5Idt8s9vwfOiN8gIn7g18CZwHTgYhGZDowF1kcPC7sclzHGFM3AKT4+MuTrDJ241utQsuJqglDVecDOhM2zgJWqulpVO4D7gHOBDUSShOtxGWNMMa0dfRbzfcN4esRlXoeSFS8exGPoKilAJDGMAf4GXCgitwJzU50sIleIyEIRWbht2zZ3IzXGmAJ4/4xJ/Lzx21x82lFeh5IVLxqpk3UEVlVtBj6T7mRVvR24HWDmzJnl2fJjjOlXBtdV8uw1J3sdRta8KEFsAMbFvR8LbMzmArbkqDHGuM+LBLEAOEBEJopIJfAx4JFsLqCqc1X1ioaGBlcCNMYY434313uBl4EDRWSDiFymqiHgKuBJYBlwv6q+6WYcxhhjsudqG4SqXpxi++PA47leV0RmA7MnT56c6yWMMcakUZbdSa2KyRhj3FeWCcIYY4z7yjJBWC8mY4xxX1kmCKtiMsYY94lq+Y41E5FtwLro2wZgTy/fJ34dCmzP4nbx18x0f+I2L2PMNr5kcSXb5mWM9u+cf3zJ4kq2zf6dSyvGfONrVNVhaSNQ1T7xAm7v7fskXxfmev1M9ydu8zLGbONLFk+pxWj/zvbvbP/OuceXyassq5hSmJvm+8Sv+Vw/0/2J27yMMdv4UsVTSjHav3Nm++zfObMY0u0vpRgLEV9aZV3FlA8RWaiqM72OozcWY/5KPT6wGAuh1OOD8ogxUV8qQWTrdq8DyIDFmL9Sjw8sxkIo9figPGLspt+WIIwxxvSuP5cgjDHG9MIShDHGmKQsQRhjjEnKEkQSInKyiPxHRG4TkZO9jicVEakTkUUicrbXsSQSkWnR39+DIvIFr+NJRkTOE5E7RORhEfmA1/EkIyKTROROEXnQ61hion9390R/d5/wOp5kSvH3lqgc/v76XIIQkbtEZKuILE3YfoaILBeRlSJybZrLKNAEVBNZAa8UYwT4FnB/KcanqstU9UrgIqDgXfsKFOM/VPVzwKXAR0s0xtWq6vpK91nGegHwYPR3d47bseUSY7F+b3nG6OrfX0FkM7KvHF7AScARwNK4bX5gFTAJqAQWA9OBQ4BHE17DAV/0vBHAn0s0xtOJrMZ3KXB2qcUXPecc4CXg46X4O4w772bgiBKP8cES+n/zbWBG9Ji/uBlXrjEW6/dWoBhd+fsrxMvVBYO8oKrzRGRCwuZZwEpVXQ0gIvcB56rqTUBv1TO7gKpSjFFETgHqiPyHbRWRx1XVKZX4otd5BHhERB4D/lKI2AoZo4gI8CPgCVV9tZDxFSrGYskmViKl6rHA6xSxFiLLGN8qVlzxsolRRJbh4t9fIfS5KqYUxgDr495viG5LSkQuEJHfAn8EfuVybDFZxaiq/6OqXyXy4L2jUMmhUPFF23HmRH+POa8emKWsYgS+TKQk9mERudLNwOJk+3scIiK3AYeLyLfdDi5Bqlj/BlwoIreS+zQShZI0Ro9/b4lS/R69+PvLSp8rQaQgSbalHCGoqn8j8p+gmLKKsfMA1d8XPpSksv0dPg8871YwKWQb4xxgjnvhJJVtjDsArx4eSWNV1WbgM8UOJoVUMXr5e0uUKkYv/v6y0l9KEBuAcXHvxwIbPYollVKPsdTjA4ux0MohVovRRf0lQSwADhCRiSJSSaRx9xGPY0pU6jGWenxgMRZaOcRqMbrJ61byQr+Ae4FNQJBI5r4suv1DwDtEehP8j8VYvvFZjP0zVoux+C+brM8YY0xS/aWKyRhjTJYsQRhjjEnKEoQxxpikLEEYY4xJyhKEMcaYpCxBGGOMScoShOkXRCQsIq/HvTKZTr0oJLJmxqRe9n9PRG5K2DYjOtkbIvIvERnkdpym/7EEYfqLVlWdEff6Ub4XFJG85zITkYMAv0Zn+kzhXnquF/AxumbI/SPwxXxjMSaRJQjTr4nIWhG5QUReFZE3RGRqdHtddPGXBSLymoicG91+qYg8ICJzgadExCcivxGRN0XkURF5XEQ+LCKnicjf4+7zfhFJNgHkJ4CH4477gIi8HI3nARGpV9XlwG4ROTruvIuA+6LfPwJcXNjfjDGWIEz/UZNQxRT/iXy7qh4B3Ap8Pbrtf4BnVfUo4BTgpyJSF913LHCJqp5KZHW1CUQW/Lk8ug/gWWCaiAyLvv8McHeSuI4HFgGIyFDgOuD0aDwLga9Fj7uXSKkBETkG2KGqKwBUdRdQJSJDcvi9GJNSf5nu25hWVZ2RYl/sk/0iIg98gA8A54hILGFUA+Oj3z+tqjuj358APKCR9Tg2i8hzEJnLWUT+CHxSRO4mkjg+neTeo4Bt0e+PIbIA1IuRtYyoBF6O7rsPeElEriGSKO5NuM5WYDSwI8XPaEzWLEEYA+3Rr2G6/k8IcGG0eqdTtJqnOX5TL9e9m8iCOm1EkkgoyTGtRJJP7FpPq2qP6iJVXS8ia4H3ARfSVVKJqY5ey5iCsSomY5J7EvhydFlSROTwFMe9QGR1NZ+IjABOju1Q1Y1E5v2/Dvh9ivOXAZOj378CHC8ik6P3rBWRKXHH3gv8HFilqhtiG6MxjgTWZvHzGZOWJQjTXyS2QaTrxfR9oAJYIiJLo++TeYjItM5Lgd8C84E9cfv/DKxX1VRrJD9GNKmo6jbgUuBeEVlCJGFMjTv2AeAguhqnY44EXklRQjEmZzbdtzF5ivY0aoo2Ev8XOF5VN0f3/Qp4TVXvTHFuDfBc9Jxwjvf/BfCIqj6T209gTHLWBmFM/h4VkUYijcrfj0sOi4i0V1yT6kRVbRWR64ksYv9ujvdfasnBuMFKEMYYY5KyNghjjDFJWYIwxhiTlCUIY4wxSVmCMMYYk5QlCGOMMUlZgjDGGJPU/w9oMXzoHSpRCgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "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 samples:\n", + " xs = sample.reconstruct(energies)\n", + " elastic_xs = xs[2]\n", + " plt.loglog(energies, elastic_xs)\n", + "plt.xlabel('Energy (eV)')\n", + "plt.ylabel('Cross section (b)')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Subset Selection" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Another capability of the covariance module is selecting a subset of the resonance parameters and the corresponding subset of the covariance matrix. We can do this by specifying the value we want to discriminate and the bounds within one energy region. Selecting only resonances with J=2:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
energyLJneutronWidthcaptureWidthfissionWidthAfissionWidthB
00.031402.00.0004740.10720.00.0
12.825002.00.0003450.09700.00.0
316.770002.00.0128000.08050.00.0
420.560002.00.0113600.08800.00.0
521.650002.00.0003760.11400.00.0
\n", + "
" + ], + "text/plain": [ + " energy L J neutronWidth captureWidth fissionWidthA fissionWidthB\n", + "0 0.0314 0 2.0 0.000474 0.1072 0.0 0.0\n", + "1 2.8250 0 2.0 0.000345 0.0970 0.0 0.0\n", + "3 16.7700 0 2.0 0.012800 0.0805 0.0 0.0\n", + "4 20.5600 0 2.0 0.011360 0.0880 0.0 0.0\n", + "5 21.6500 0 2.0 0.000376 0.1140 0.0 0.0" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "lower_bound = 2; # inclusive\n", + "upper_bound = 2; # inclusive\n", + "rm_res_cov_sub = gd157_endf.resonance_covariance.ranges[0].subset('J',[lower_bound,upper_bound])\n", + "rm_res_cov_sub.file2res.parameters[:5]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The subset method will also store the corresponding subset of the covariance matrix" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(180, 180)" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rm_res_cov_sub.covariance\n", + "gd157_endf.resonance_covariance.ranges[0].covariance.shape\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Checking the size of the new covariance matrix to be sure it was sampled properly: " + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of parameters\n", + "Original: 60\n", + "Subet: 36\n", + "Covariance Size\n", + "Original: (180, 180)\n", + "Subset: (108, 108)\n" + ] + } + ], + "source": [ + "old_n_parameters = gd157_endf.resonance_covariance.ranges[0].parameters.shape[0]\n", + "old_shape = gd157_endf.resonance_covariance.ranges[0].covariance.shape\n", + "new_n_parameters = rm_res_cov_sub.file2res.parameters.shape[0]\n", + "new_shape = rm_res_cov_sub.covariance.shape\n", + "print('Number of parameters\\nOriginal: '+str(old_n_parameters)+'\\nSubet: '+str(new_n_parameters)+'\\nCovariance Size\\nOriginal: '+str(old_shape)+'\\nSubset: '+str(new_shape))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And finally, we can sample from the subset as well" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/icmeyer/openmc/openmc/data/resonance_covariance.py:233: UserWarning: Sampling routine does not guarantee positive values for parameters. This can lead to undefined behavior in the reconstruction routine.\n", + " warnings.warn(warn_str)\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
energyLJneutronWidthcaptureWidthfissionWidthAfissionWidthB
00.03106502.00.0004740.1079540.00.0
12.82380902.00.0003340.0982180.00.0
216.76918602.00.0139870.0729100.00.0
320.55664902.00.0108140.0987800.00.0
421.65459102.00.0003660.1176790.00.0
\n", + "
" + ], + "text/plain": [ + " energy L J neutronWidth captureWidth fissionWidthA fissionWidthB\n", + "0 0.031065 0 2.0 0.000474 0.107954 0.0 0.0\n", + "1 2.823809 0 2.0 0.000334 0.098218 0.0 0.0\n", + "2 16.769186 0 2.0 0.013987 0.072910 0.0 0.0\n", + "3 20.556649 0 2.0 0.010814 0.098780 0.0 0.0\n", + "4 21.654591 0 2.0 0.000366 0.117679 0.0 0.0" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "samples_sub = rm_res_cov_sub.sample(n_samples)\n", + "samples_sub[0].parameters[:5]" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index 7158e9fe3f..44c59628c9 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -27,5 +27,6 @@ from .urr import * from .library import * from .fission_energy import * from .resonance import * +from .resonance_covariance import * from .multipole import * from .grid import * diff --git a/openmc/data/endf.py b/openmc/data/endf.py index 0d1f402c2f..8a118f1695 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -71,6 +71,24 @@ def float_endf(s): return float(_ENDF_FLOAT_RE.sub(r'\1e\2', s)) +def _int_endf(s): + """Convert string to int. Used for INTG records where blank entries + indicate a 0. + + Parameters + ---------- + s : str + Integer or spaces + + Returns + ------- + integer + The number or 0 + """ + s = s.strip() + return int(s) if s else 0 + + def get_text_record(file_obj): """Return data from a TEXT record in an ENDF-6 file. @@ -250,6 +268,50 @@ def get_tab2_record(file_obj): return params, Tabulated2D(breakpoints, interpolation) +def get_intg_record(file_obj): + """ + Return data from an INTG record in an ENDF-6 file. Used to store the + covariance matrix in a compact format. + + Parameters + ---------- + file_obj : file-like object + ENDF-6 file to read from + + Returns + ------- + numpy.ndarray + The correlation matrix described in the INTG record + """ + # determine how many items are in list and NDIGIT + items = get_cont_record(file_obj) + ndigit = int(items[2]) + npar = int(items[3]) # Number of parameters + nlines = int(items[4]) # Lines to read + NROW_RULES = {2: 18, 3: 12, 4: 11, 5: 9, 6: 8} + nrow = NROW_RULES[ndigit] + + # read lines and build correlation matrix + corr = np.identity(npar) + for i in range(nlines): + line = file_obj.readline() + ii = _int_endf(line[:5]) - 1 # -1 to account for 0 indexing + jj = _int_endf(line[5:10]) - 1 + factor = 10**ndigit + for j in range(nrow): + if jj+j >= ii: + break + element = _int_endf(line[11+(ndigit+1)*j:11+(ndigit+1)*(j+1)]) + if element > 0: + corr[ii, jj] = (element+0.5)/factor + elif element < 0: + corr[ii, jj] = (element-0.5)/factor + + # Symmetrize the correlation matrix + corr = corr + corr.T - np.diag(corr.diagonal()) + return corr + + def get_evaluations(filename): """Return a list of all evaluations within an ENDF file. @@ -288,7 +350,7 @@ class Evaluation(object): Attributes ---------- info : dict - Miscallaneous information about the evaluation. + Miscellaneous information about the evaluation. target : dict Information about the target material, such as its mass, isomeric state, whether it's stable, and whether it's fissionable. diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 1cc0e38879..3917496b8e 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -24,6 +24,7 @@ from .njoy import make_ace from .product import Product from .reaction import Reaction, _get_photon_products_ace from . import resonance as res +from . import resonance_covariance as res_cov from .urr import ProbabilityTables import openmc.checkvalue as cv from openmc.mixin import EqualityMixin @@ -148,6 +149,8 @@ class IncidentNeutron(EqualityMixin): and the values are Reaction objects. resonances : openmc.data.Resonances or None Resonance parameters + resonance_covariance : openmc.data.ResonanceCovariance or None + Covariance for resonance parameters summed_reactions : collections.OrderedDict Contains summed cross sections, e.g., the total cross section. The keys are the MT values and the values are Reaction objects. @@ -228,6 +231,10 @@ class IncidentNeutron(EqualityMixin): def resonances(self): return self._resonances + @property + def resonance_covariance(self): + return self._resonance_covariance + @property def summed_reactions(self): return self._summed_reactions @@ -289,6 +296,12 @@ class IncidentNeutron(EqualityMixin): cv.check_type('resonances', resonances, res.Resonances) self._resonances = resonances + @resonance_covariance.setter + def resonance_covariance(self, resonance_covariance): + cv.check_type('resonance covariance', resonance_covariance, + res_cov.ResonanceCovariances) + self._resonance_covariance = resonance_covariance + @summed_reactions.setter def summed_reactions(self, summed_reactions): cv.check_type('summed reactions', summed_reactions, Mapping) @@ -744,7 +757,7 @@ class IncidentNeutron(EqualityMixin): return data @classmethod - def from_endf(cls, ev_or_filename): + def from_endf(cls, ev_or_filename, covariance=False): """Generate incident neutron continuous-energy data from an ENDF evaluation Parameters @@ -753,6 +766,10 @@ class IncidentNeutron(EqualityMixin): ENDF evaluation to read from. If given as a string, it is assumed to be the filename for the ENDF file. + covariance : bool + Flag to indicate whether or not covariance data from File 32 should be + retrieved + Returns ------- openmc.data.IncidentNeutron @@ -784,6 +801,11 @@ class IncidentNeutron(EqualityMixin): if (2, 151) in ev.section: data.resonances = res.Resonances.from_endf(ev) + if (32, 151) in ev.section and covariance: + data.resonance_covariance = ( + res_cov.ResonanceCovariances.from_endf(ev, data.resonances) + ) + # Read each reaction for mf, mt, nc, mod in ev.reaction_list: if mf == 3: diff --git a/openmc/data/resonance.py b/openmc/data/resonance.py index d58f706eb9..5e4bd7129e 100644 --- a/openmc/data/resonance.py +++ b/openmc/data/resonance.py @@ -16,6 +16,7 @@ except ImportError: _reconstruct = False import openmc.checkvalue as cv + class Resonances(object): """Resolved and unresolved resonance data @@ -90,14 +91,14 @@ class Resonances(object): # Determine whether discrete or continuous representation items = get_head_record(file_obj) - n_isotope = items[4] # Number of isotopes + n_isotope = items[4] # Number of isotopes ranges = [] for iso in range(n_isotope): items = get_cont_record(file_obj) abundance = items[1] - fission_widths = (items[3] == 1) # fission widths are given? - n_ranges = items[4] # number of resonance energy ranges + fission_widths = (items[3] == 1) # fission widths are given? + n_ranges = items[4] # number of resonance energy ranges for j in range(n_ranges): items = get_cont_record(file_obj) @@ -112,7 +113,7 @@ class Resonances(object): # unresolved resonance region erange = Unresolved.from_endf(file_obj, items, fission_widths) - #erange.material = self + # erange.material = self ranges.append(erange) return cls(ranges) @@ -162,6 +163,13 @@ class ResonanceRange(object): self._prepared = False self._parameter_matrix = {} + def __copy__(self): + cls = type(self) + new_copy = cls.__new__(cls) + new_copy.__dict__.update(self.__dict__) + new_copy._prepared = False + return new_copy + @classmethod def from_endf(cls, ev, file_obj, items): """Create resonance range from an ENDF evaluation. @@ -437,7 +445,7 @@ class MultiLevelBreitWigner(ResonanceRange): self._l_values = np.array(l_values) self._competitive = np.array(competitive) for l in l_values: - self._parameter_matrix[l] = df[df.L == l].as_matrix() + self._parameter_matrix[l] = df[df.L == l].values self._prepared = True @@ -682,7 +690,7 @@ class ReichMoore(ResonanceRange): self._l_values = np.array(l_values) for (l, J) in lj_values: self._parameter_matrix[l, J] = df[(df.L == l) & - (abs(df.J) == J)].as_matrix() + (abs(df.J) == J)].values self._prepared = True diff --git a/openmc/data/resonance_covariance.py b/openmc/data/resonance_covariance.py new file mode 100644 index 0000000000..300e6dbf60 --- /dev/null +++ b/openmc/data/resonance_covariance.py @@ -0,0 +1,708 @@ +from collections import MutableSequence +import warnings +import io +import copy + +import numpy as np +import pandas as pd + +from . import endf +import openmc.checkvalue as cv +from .resonance import Resonances + + +def _add_file2_contributions(file32params, file2params): + """Function for aiding in adding resonance parameters from File 2 that are + not always present in File 32. Uses already imported resonance data. + + Paramaters + ---------- + file32params : pandas.Dataframe + Incomplete set of resonance parameters contained in File 32. + file2params : pandas.Dataframe + Resonance parameters from File 2. Ordered by energy. + + Returns + ------- + parameters : pandas.Dataframe + Complete set of parameters ordered by L-values and then energy + + """ + # Use l-values and competitiveWidth from File 2 data + # Re-sort File 2 by energy to match File 32 + file2params = file2params.sort_values(by=['energy']) + file2params.reset_index(drop=True, inplace=True) + # Sort File 32 parameters by energy as well (maintaining index) + file32params.sort_values(by=['energy'], inplace=True) + # Add in values (.values converts to array first to ignore index) + file32params['L'] = file2params['L'].values + if 'competitiveWidth' in file2params.columns: + file32params['competitiveWidth'] = file2params['competitiveWidth'].values + # Resort to File 32 order (by L then by E) for use with covariance + file32params.sort_index(inplace=True) + return file32params + + +class ResonanceCovariances(Resonances): + """Resolved resonance covariance data + + Parameters + ---------- + ranges : list of openmc.data.ResonanceCovarianceRange + Distinct energy ranges for resonance data + + Attributes + ---------- + ranges : list of openmc.data.ResonanceCovarianceRange + Distinct energy ranges for resonance data + + """ + + @property + def ranges(self): + return self._ranges + + @ranges.setter + def ranges(self, ranges): + cv.check_type('resonance ranges', ranges, MutableSequence) + self._ranges = cv.CheckedList(ResonanceCovarianceRange, + 'resonance range', ranges) + + @classmethod + def from_endf(cls, ev, resonances): + """Generate resonance covariance data from an ENDF evaluation. + + Parameters + ---------- + ev : openmc.data.endf.Evaluation + ENDF evaluation + resonances : openmc.data.Resonance object + openmc.data.Resonanance object generated from the same evaluation + used to import values not contained in File 32 + + Returns + ------- + openmc.data.ResonanceCovariances + Resonance covariance data + + """ + file_obj = io.StringIO(ev.section[32, 151]) + + # Determine whether discrete or continuous representation + items = endf.get_head_record(file_obj) + n_isotope = items[4] # Number of isotopes + + ranges = [] + for iso in range(n_isotope): + items = endf.get_cont_record(file_obj) + abundance = items[1] + fission_widths = (items[3] == 1) # Flag for fission widths + n_ranges = items[4] # Number of resonance energy ranges + + for j in range(n_ranges): + items = endf.get_cont_record(file_obj) + # Unresolved flags - 0: only scattering radius given + # 1: resolved parameters given + # 2: unresolved parameters given + unresolved_flag = items[2] + formalism = items[3] # resonance formalism + + # Throw error for unsupported formalisms + if formalism in [0, 7]: + error = 'LRF='+str(formalism)+' covariance not supported '\ + 'for this formalism' + raise NotImplementedError(error) + + if unresolved_flag in (0, 1): + # Resolved resonance region + resonance = resonances.ranges[j] + erange = _FORMALISMS[formalism].from_endf(ev, file_obj, + items, resonance) + ranges.append(erange) + + elif unresolved_flag == 2: + warn = 'Unresolved resonance not supported. Covariance '\ + 'values for the unresolved region not imported.' + warnings.warn(warn) + + return cls(ranges) + + +class ResonanceCovarianceRange: + """Resonace covariance range. Base class for different formalisms. + + Parameters + ---------- + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + + Attributes + ---------- + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + parameters : pandas.DataFrame + Resonance parameters + covariance : numpy.array + The covariance matrix contained within the ENDF evaluation + lcomp : int + Flag indicating format of the covariance matrix within the ENDF file + file2res : openmc.data.ResonanceRange object + Corresponding resonance range with File 2 data. + mpar : int + Number of parameters in covariance matrix for each individual resonance + formalism : str + String descriptor of formalism + """ + def __init__(self, energy_min, energy_max): + self.energy_min = energy_min + self.energy_max = energy_max + + def subset(self, parameter_str, bounds): + """Produce a subset of resonance parameters and the corresponding + covariance matrix to an IncidentNeutron object. + + Parameters + ---------- + parameter_str : str + parameter to be discriminated + (i.e. 'energy', 'captureWidth', 'fissionWidthA'...) + bounds : np.array + [low numerical bound, high numerical bound] + + Returns + ------- + res_cov_range : openmc.data.ResonanceCovarianceRange + ResonanceCovarianceRange object that contains a subset of the + covariance matrix (upper triangular) as well as a subset parameters + within self.file2params + + """ + # Copy range and prevent change of original + res_cov_range = copy.deepcopy(self) + + parameters = self.file2res.parameters + cov = res_cov_range.covariance + mpar = res_cov_range.mpar + # Create mask + mask1 = parameters[parameter_str] >= bounds[0] + mask2 = parameters[parameter_str] <= bounds[1] + mask = mask1 & mask2 + res_cov_range.parameters = parameters[mask] + indices = res_cov_range.parameters.index.values + # Build subset of covariance + sub_cov_dim = len(indices)*mpar + cov_subset_vals = [] + for index1 in indices: + for i in range(mpar): + for index2 in indices: + for j in range(mpar): + if index2*mpar+j >= index1*mpar+i: + cov_subset_vals.append(cov[index1*mpar+i, + index2*mpar+j]) + + cov_subset = np.zeros([sub_cov_dim, sub_cov_dim]) + tri_indices = np.triu_indices(sub_cov_dim) + cov_subset[tri_indices] = cov_subset_vals + + res_cov_range.file2res.parameters = parameters[mask] + res_cov_range.covariance = cov_subset + return res_cov_range + + def sample(self, n_samples): + """Sample resonance parameters based on the covariances provided + within an ENDF evaluation. + + Parameters + ---------- + n_samples : int + The number of samples to produce + + Returns + ------- + samples : list of openmc.data.ResonanceCovarianceRange objects + List of samples size `n_samples` + + """ + warn_str = 'Sampling routine does not guarantee positive values for '\ + 'parameters. This can lead to undefined behavior in the '\ + 'reconstruction routine.' + warnings.warn(warn_str) + parameters = self.parameters + cov = self.covariance + + # Symmetrizing covariance matrix + cov = cov + cov.T - np.diag(cov.diagonal()) + formalism = self.formalism + mpar = self.mpar + samples = [] + + # Handling MLBW/SLBW sampling + if formalism == 'mlbw' or formalism == 'slbw': + params = ['energy', 'neutronWidth', 'captureWidth', 'fissionWidth', + 'competitiveWidth'] + param_list = params[:mpar] + mean_array = parameters[param_list].values + mean = mean_array.flatten() + par_samples = np.random.multivariate_normal(mean, cov, + size=n_samples) + spin = parameters['J'].values + l_value = parameters['L'].values + for sample in par_samples: + energy = sample[0::mpar] + gn = sample[1::mpar] + gg = sample[2::mpar] + gf = sample[3::mpar] if mpar > 3 else parameters['fissionWidth'].values + gx = sample[4::mpar] if mpar > 4 else parameters['competitiveWidth'].values + gt = gn + gg + gf + gx + + records = [] + for j, E in enumerate(energy): + records.append([energy[j], l_value[j], spin[j], gt[j], + gn[j], gg[j], gf[j], gx[j]]) + columns = ['energy', 'L', 'J', 'totalWidth', 'neutronWidth', + 'captureWidth', 'fissionWidth', 'competitiveWidth'] + sample_params = pd.DataFrame.from_records(records, + columns=columns) + # Copy ResonanceRange object + res_range = copy.copy(self.file2res) + res_range.parameters = sample_params + samples.append(res_range) + + # Handling RM sampling + elif formalism == 'rm': + params = ['energy', 'neutronWidth', 'captureWidth', + 'fissionWidthA', 'fissionWidthB'] + param_list = params[:mpar] + mean_array = parameters[param_list].values + mean = mean_array.flatten() + par_samples = np.random.multivariate_normal(mean, cov, + size=n_samples) + spin = parameters['J'].values + l_value = parameters['L'].values + for sample in par_samples: + energy = sample[0::mpar] + gn = sample[1::mpar] + gg = sample[2::mpar] + gfa = sample[3::mpar] if mpar > 3 else parameters['fissionWidthA'].values + gfb = sample[4::mpar] if mpar > 3 else parameters['fissionWidthB'].values + + records = [] + for j, E in enumerate(energy): + records.append([energy[j], l_value[j], spin[j], gn[j], + gg[j], gfa[j], gfb[j]]) + columns = ['energy', 'L', 'J', 'neutronWidth', + 'captureWidth', 'fissionWidthA', 'fissionWidthB'] + sample_params = pd.DataFrame.from_records(records, + columns=columns) + # Copy ResonanceRange object + res_range = copy.copy(self.file2res) + res_range.parameters = sample_params + samples.append(res_range) + + return samples + + +class MultiLevelBreitWignerCovariance(ResonanceCovarianceRange): + """Multi-level Breit-Wigner resolved resonance formalism covariance data. + Parameters + ---------- + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + + Attributes + ---------- + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + parameters : pandas.DataFrame + Resonance parameters + covariance : numpy.array + The covariance matrix contained within the ENDF evaluation + mpar : int + Number of parameters in covariance matrix for each individual resonance + lcomp : int + Flag indicating format of the covariance matrix within the ENDF file + file2res : openmc.data.ResonanceRange object + Corresponding resonance range with File 2 data. + formalism : str + String descriptor of formalism + + """ + + def __init__(self, energy_min, energy_max, parameters, covariance, mpar, + lcomp, file2res): + super().__init__(energy_min, energy_max) + self.parameters = parameters + self.covariance = covariance + self.mpar = mpar + self.lcomp = lcomp + self.file2res = copy.copy(file2res) + self.formalism = 'mlbw' + + @classmethod + def from_endf(cls, ev, file_obj, items, resonance): + """Create MLBW covariance data from an ENDF evaluation. + + Parameters + ---------- + ev : openmc.data.endf.Evaluation + ENDF evaluation + file_obj : file-like object + ENDF file positioned at the second record of a resonance range + subsection in MF=32, MT=151 + items : list + Items from the CONT record at the start of the resonance range + subsection + resonance : openmc.data.ResonanceRange object + Corresponding resonance range with File 2 data. + + Returns + ------- + openmc.data.MultiLevelBreitWignerCovariance + Multi-level Breit-Wigner resonance covariance parameters + + """ + + # Read energy-dependent scattering radius if present + energy_min, energy_max = items[0:2] + nro, naps = items[4:6] + if nro != 0: + params, ape = endf.get_tab1_record(file_obj) + + # Other scatter radius parameters + items = endf.get_cont_record(file_obj) + target_spin = items[0] + lcomp = items[3] # Flag for compatibility 0, 1, 2 - 2 is compact form + nls = items[4] # number of l-values + + # Build covariance matrix for General Resolved Resonance Formats + if lcomp == 1: + items = endf.get_cont_record(file_obj) + # Number of short range type resonance covariances + num_short_range = items[4] + # Number of long range type resonance covariances + num_long_range = items[5] + + # Read resonance widths, J values, etc + records = [] + for i in range(num_short_range): + items, values = endf.get_list_record(file_obj) + mpar = items[2] + num_res = items[5] + num_par_vals = num_res*6 + res_values = values[:num_par_vals] + cov_values = values[num_par_vals:] + + energy = res_values[0::6] + spin = res_values[1::6] + gt = res_values[2::6] + gn = res_values[3::6] + gg = res_values[4::6] + gf = res_values[5::6] + + for i, E in enumerate(energy): + records.append([energy[i], spin[i], gt[i], gn[i], + gg[i], gf[i]]) + + # Build the upper-triangular covariance matrix + cov_dim = mpar*num_res + cov = np.zeros([cov_dim, cov_dim]) + indices = np.triu_indices(cov_dim) + cov[indices] = cov_values + + # Compact format - Resonances and individual uncertainties followed by + # compact correlations + elif lcomp == 2: + items, values = endf.get_list_record(file_obj) + mean = items + num_res = items[5] + energy = values[0::12] + spin = values[1::12] + gt = values[2::12] + gn = values[3::12] + gg = values[4::12] + gf = values[5::12] + par_unc = [] + for i in range(num_res): + res_unc = values[i*12+6 : i*12+12] + # Delete 0 values (not provided, no fission width) + # DAJ/DGT always zero, DGF sometimes nonzero [1, 2, 5] + res_unc_nonzero = [] + for j in range(6): + if j in [1, 2, 5] and res_unc[j] != 0.0: + res_unc_nonzero.append(res_unc[j]) + elif j in [0, 3, 4]: + res_unc_nonzero.append(res_unc[j]) + par_unc.extend(res_unc_nonzero) + + records = [] + for i, E in enumerate(energy): + records.append([energy[i], spin[i], gt[i], gn[i], + gg[i], gf[i]]) + + corr = endf.get_intg_record(file_obj) + cov = np.diag(par_unc).dot(corr).dot(np.diag(par_unc)) + + # Compatible resolved resonance format + elif lcomp == 0: + cov = np.zeros([4, 4]) + records = [] + cov_index = 0 + for i in range(nls): + items, values = endf.get_list_record(file_obj) + num_res = items[5] + for j in range(num_res): + one_res = values[18*j:18*(j+1)] + res_values = one_res[:6] + cov_values = one_res[6:] + records.append(list(res_values)) + + # Populate the coviariance matrix for this resonance + # There are no covariances between resonances in lcomp=0 + cov[cov_index, cov_index] = cov_values[0] + cov[cov_index+1, cov_index+1 : cov_index+2] = cov_values[1:2] + cov[cov_index+1, cov_index+3] = cov_values[4] + cov[cov_index+2, cov_index+2] = cov_values[3] + cov[cov_index+2, cov_index+3] = cov_values[5] + cov[cov_index+3, cov_index+3] = cov_values[6] + + cov_index += 4 + if j < num_res-1: # Pad matrix for additional values + cov = np.pad(cov, ((0, 4), (0, 4)), 'constant', + constant_values=0) + + # Create pandas DataFrame with resonance data, currently + # redundant with data.IncidentNeutron.resonance + columns = ['energy', 'J', 'totalWidth', 'neutronWidth', + 'captureWidth', 'fissionWidth'] + parameters = pd.DataFrame.from_records(records, columns=columns) + # Determine mpar (number of parameters for each resonance in + # covariance matrix) + nparams, params = parameters.shape + covsize = cov.shape[0] + mpar = int(covsize/nparams) + # Add parameters from File 2 + parameters = _add_file2_contributions(parameters, + resonance.parameters) + # Create instance of class + mlbw = cls(energy_min, energy_max, parameters, cov, mpar, lcomp, + resonance) + return mlbw + + +class SingleLevelBreitWignerCovariance(MultiLevelBreitWignerCovariance): + """Single-level Breit-Wigner resolved resonance formalism covariance data. + Single-level Breit-Wigner resolved resonance data is is identified by LRF=1 + in the ENDF-6 format. + + Parameters + ---------- + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + + Attributes + ---------- + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + parameters : pandas.DataFrame + Resonance parameters + covariance : numpy.array + The covariance matrix contained within the ENDF evaluation + mpar : int + Number of parameters in covariance matrix for each individual resonance + formalism : str + String descriptor of formalism + lcomp : int + Flag indicating format of the covariance matrix within the ENDF file + file2res : openmc.data.ResonanceRange object + Corresponding resonance range with File 2 data. + """ + + def __init__(self, energy_min, energy_max, parameters, covariance, mpar, + lcomp, file2res): + super().__init__(energy_min, energy_max, parameters, covariance, mpar, + lcomp, file2res) + self.formalism = 'slbw' + + +class ReichMooreCovariance(ResonanceCovarianceRange): + """Reich-Moore resolved resonance formalism covariance data. + + Reich-Moore resolved resonance data is identified by LRF=3 in the ENDF-6 + format. + + Parameters + ---------- + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + + Attributes + ---------- + energy_min : float + Minimum energy of the resolved resonance range in eV + energy_max : float + Maximum energy of the resolved resonance range in eV + parameters : pandas.DataFrame + Resonance parameters + covariance : numpy.array + The covariance matrix contained within the ENDF evaluation + lcomp : int + Flag indicating format of the covariance matrix within the ENDF file + mpar : int + Number of parameters in covariance matrix for each individual resonance + file2res : openmc.data.ResonanceRange object + Corresponding resonance range with File 2 data. + formalism : str + String descriptor of formalism + """ + + def __init__(self, energy_min, energy_max, parameters, covariance, mpar, + lcomp, file2res): + super().__init__(energy_min, energy_max) + self.parameters = parameters + self.covariance = covariance + self.mpar = mpar + self.lcomp = lcomp + self.file2res = copy.copy(file2res) + self.formalism = 'rm' + + @classmethod + def from_endf(cls, ev, file_obj, items, resonance): + """Create Reich-Moore resonance covariance data from an ENDF + evaluation. Includes the resonance parameters contained separately in + File 32. + + Parameters + ---------- + ev : openmc.data.endf.Evaluation + ENDF evaluation + file_obj : file-like object + ENDF file positioned at the second record of a resonance range + subsection in MF=2, MT=151 + items : list + Items from the CONT record at the start of the resonance range + subsection + resonance : openmc.data.Resonance object + openmc.data.Resonanance object generated from the same evaluation + used to import values not contained in File 32 + + Returns + ------- + openmc.data.ReichMooreCovariance + Reich-Moore resonance covariance parameters + + """ + # Read energy-dependent scattering radius if present + energy_min, energy_max = items[0:2] + nro, naps = items[4:6] + if nro != 0: + params, ape = endf.get_tab1_record(file_obj) + + # Other scatter radius parameters + items = endf.get_cont_record(file_obj) + target_spin = items[0] + lcomp = items[3] # Flag for compatibility 0, 1, 2 - 2 is compact form + nls = items[4] # Number of l-values + + # Build covariance matrix for General Resolved Resonance Formats + if lcomp == 1: + items = endf.get_cont_record(file_obj) + # Number of short range type resonance covariances + num_short_range = items[4] + # Number of long range type resonance covariances + num_long_range = items[5] + # Read resonance widths, J values, etc + channel_radius = {} + scattering_radius = {} + records = [] + for i in range(num_short_range): + items, values = endf.get_list_record(file_obj) + mpar = items[2] + num_res = items[5] + num_par_vals = num_res*6 + res_values = values[:num_par_vals] + cov_values = values[num_par_vals:] + + energy = res_values[0::6] + spin = res_values[1::6] + gn = res_values[2::6] + gg = res_values[3::6] + gfa = res_values[4::6] + gfb = res_values[5::6] + + for i, E in enumerate(energy): + records.append([energy[i], spin[i], gn[i], gg[i], + gfa[i], gfb[i]]) + + # Build the upper-triangular covariance matrix + cov_dim = mpar*num_res + cov = np.zeros([cov_dim, cov_dim]) + indices = np.triu_indices(cov_dim) + cov[indices] = cov_values + + # Compact format - Resonances and individual uncertainties followed by + # compact correlations + elif lcomp == 2: + items, values = endf.get_list_record(file_obj) + num_res = items[5] + energy = values[0::12] + spin = values[1::12] + gn = values[2::12] + gg = values[3::12] + gfa = values[4::12] + gfb = values[5::12] + par_unc = [] + for i in range(num_res): + res_unc = values[i*12+6 : i*12+12] + # Delete 0 values (not provided in evaluation) + res_unc = [x for x in res_unc if x != 0.0] + par_unc.extend(res_unc) + + records = [] + for i, E in enumerate(energy): + records.append([energy[i], spin[i], gn[i], gg[i], + gfa[i], gfb[i]]) + + corr = endf.get_intg_record(file_obj) + cov = np.diag(par_unc).dot(corr).dot(np.diag(par_unc)) + + # Create pandas DataFrame with resonacne data + columns = ['energy', 'J', 'neutronWidth', 'captureWidth', + 'fissionWidthA', 'fissionWidthB'] + parameters = pd.DataFrame.from_records(records, columns=columns) + + # Determine mpar (number of parameters for each resonance in + # covariance matrix) + nparams, params = parameters.shape + covsize = cov.shape[0] + mpar = int(covsize/nparams) + + # Add parameters from File 2 + parameters = _add_file2_contributions(parameters, + resonance.parameters) + # Create instance of ReichMooreCovariance + rmc = cls(energy_min, energy_max, parameters, cov, mpar, lcomp, + resonance) + return rmc + + +_FORMALISMS = { + 0: ResonanceCovarianceRange, + 1: SingleLevelBreitWignerCovariance, + 2: MultiLevelBreitWignerCovariance, + 3: ReichMooreCovariance + # 7: RMatrixLimitedCovariance +} diff --git a/tests/unit_tests/test_data_neutron.py b/tests/unit_tests/test_data_neutron.py index 03746430d9..406ff0515f 100644 --- a/tests/unit_tests/test_data_neutron.py +++ b/tests/unit_tests/test_data_neutron.py @@ -37,9 +37,10 @@ def sm150(): @pytest.fixture(scope='module') def gd154(): - """Gd154 ENDF data (contains Reich Moore resonance range)""" + """Gd154 ENDF data (contains Reich Moore resonance range and reosnance + covariance with LCOMP=1).""" filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-064_Gd_154.endf') - return openmc.data.IncidentNeutron.from_endf(filename) + return openmc.data.IncidentNeutron.from_endf(filename, covariance=True) @pytest.fixture(scope='module') @@ -77,6 +78,13 @@ def na22(): return openmc.data.IncidentNeutron.from_endf(filename) +@pytest.fixture(scope='module') +def na23(): + """Na23 ENDF data (contains MLBW resonance covariance with LCOMP=0).""" + filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-011_Na_023.endf') + return openmc.data.IncidentNeutron.from_endf(filename, covariance=True) + + @pytest.fixture(scope='module') def be9(): """Be9 ENDF data (contains laboratory angle-energy distribution).""" @@ -97,6 +105,28 @@ def am244(): return openmc.data.IncidentNeutron.from_njoy(endf_file) +@pytest.fixture(scope='module') +def ti50(): + """Ti50 ENDF data (contains Multi-level Breit-Wigner resonance range and + resonance covariance with LCOMP=1).""" + filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-022_Ti_050.endf') + return openmc.data.IncidentNeutron.from_endf(filename, covariance=True) + + +@pytest.fixture(scope='module') +def cf252(): + """Cf252 ENDF data (contains RM resonance covariance with LCOMP=0).""" + filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-098_Cf_252.endf') + return openmc.data.IncidentNeutron.from_endf(filename, covariance=True) + + +@pytest.fixture(scope='module') +def th232(): + """Th232 ENDF data (contains RM resonance covariance with LCOMP=2).""" + filename = os.path.join(_ENDF_DATA, 'neutrons', 'n-090_Th_232.endf') + return openmc.data.IncidentNeutron.from_endf(filename, covariance=True) + + def test_attributes(pu239): assert pu239.name == 'Pu239' assert pu239.mass_number == 239 @@ -277,6 +307,101 @@ def test_rml(cl35): assert isinstance(group, openmc.data.SpinGroup) +def test_mlbw_cov_lcomp0(cf252): + # Testing on first range only + cov = cf252.resonance_covariance.ranges[0] + res = cf252.resonances.ranges[0] + assert cov.parameters['energy'][0] == pytest.approx(-3.5) + assert res.parameters['energy'][0] == cov.parameters['energy'][0] + assert isinstance(cov, openmc.data.resonance_covariance.MultiLevelBreitWignerCovariance) + assert cov.energy_min == pytest.approx(1e-5) + assert cov.energy_max == pytest.approx(1000.) + assert cov.covariance[0,0] == pytest.approx(1.225e-05) + + subset = cov.subset('energy', [0, 100]) + assert not subset.parameters.empty + assert (subset.file2res.parameters['energy'] < 100).all() + samples = cov.sample(1) + xs = samples[0].reconstruct([10., 100., 1000.]) + assert sorted(xs.keys()) == [2, 18, 102] + + +def test_mlbw_cov_lcomp1(ti50): + # Testing on first range only + cov = ti50.resonance_covariance.ranges[0] + res = ti50.resonances.ranges[0] + assert cov.parameters['energy'][0] == pytest.approx(-21020.) + assert res.parameters['energy'][0] == cov.parameters['energy'][0] + assert isinstance(cov, openmc.data.resonance_covariance.MultiLevelBreitWignerCovariance) + assert cov.energy_min == pytest.approx(1e-5) + assert cov.energy_max == pytest.approx(587000.) + assert cov.covariance[0,0] == pytest.approx(1.410177e5) + + subset = cov.subset('L', [1, 1]) + assert not subset.parameters.empty + assert (subset.file2res.parameters['L'] == 1).all() + samples = cov.sample(1) + xs = samples[0].reconstruct([10., 100., 1000.]) + assert sorted(xs.keys()) == [2, 18, 102] + + +def test_mlbw_cov_lcomp2(na23): + # Testing on first range only + cov = na23.resonance_covariance.ranges[0] + res = na23.resonances.ranges[0] + assert cov.parameters['energy'][0] == pytest.approx(2810.) + assert res.parameters['energy'][0] == cov.parameters['energy'][0] + assert isinstance(cov, openmc.data.resonance_covariance.MultiLevelBreitWignerCovariance) + assert cov.energy_min == pytest.approx(600) + assert cov.energy_max == pytest.approx(500000.) + assert cov.covariance[0,0] == pytest.approx(16.1064163584) + + subset = cov.subset('L', [1, 1]) + assert not subset.parameters.empty + assert (subset.file2res.parameters['L'] == 1).all() + samples = cov.sample(1) + xs = samples[0].reconstruct([10., 100., 1000.]) + assert sorted(xs.keys()) == [2, 18, 102] + + +def test_rmcov_lcomp1(gd154): + # Testing on first range only + cov = gd154.resonance_covariance.ranges[0] + res = gd154.resonances.ranges[0] + assert cov.parameters['energy'][0] == pytest.approx(-2.200001) + assert res.parameters['energy'][0] == cov.parameters['energy'][0] + assert isinstance(cov, openmc.data.resonance_covariance.ReichMooreCovariance) + assert cov.energy_min == pytest.approx(1e-5) + assert cov.energy_max == pytest.approx(2760.) + assert cov.covariance[0,0] == pytest.approx(0.8895997) + + subset = cov.subset('energy', [0, 100]) + assert not subset.parameters.empty + assert (subset.file2res.parameters['energy'] < 100).all() + samples = cov.sample(1) + xs = samples[0].reconstruct([10., 100., 1000.]) + assert sorted(xs.keys()) == [2, 18, 102] + + +def test_rmcov_lcomp2(th232): + # Testing on first range only + cov = th232.resonance_covariance.ranges[0] + res = th232.resonances.ranges[0] + assert cov.parameters['energy'][0] == pytest.approx(-2000) + assert res.parameters['energy'][0] == cov.parameters['energy'][0] + assert isinstance(cov, openmc.data.resonance_covariance.ReichMooreCovariance) + assert cov.energy_min == pytest.approx(1e-5) + assert cov.energy_max == pytest.approx(4000.) + assert cov.covariance[0,0] == pytest.approx(246.6043092496) + + subset = cov.subset('energy', [0, 100]) + assert not subset.parameters.empty + assert (subset.file2res.parameters['energy'] < 100).all() + samples = cov.sample(1) + xs = samples[0].reconstruct([10., 100., 1000.]) + assert sorted(xs.keys()) == [2, 18, 102] + + def test_madland_nix(am241): fission = am241.reactions[18] prompt_neutron = fission.products[0]