{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting started"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we will briefly demonstrate loading in ellipsometry datasets, creating a model and optimising that model to the dataset. For a more detailed description of model creation, please see the [getting started](https://refnx.readthedocs.io/en/latest/getting_started.html) tutorial on *refnx*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting an ellipsometry dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by importing all of the relevant packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "from os.path import join as pjoin\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "\n",
    "import refnx\n",
    "from refnx.analysis import CurveFitter\n",
    "from refnx.reflect import Slab\n",
    "\n",
    "import refellips\n",
    "from refellips.dataSE import DataSE, open_EP4file\n",
    "from refellips.reflect_modelSE import ReflectModelSE\n",
    "from refellips.objectiveSE import ObjectiveSE\n",
    "from refellips.dispersion import RI, Cauchy, load_material"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For reproducibility, it is important to note the versions of software that you're using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "refellips: 0.0.3.dev0+92f1c50\n",
      "refnx: 0.1.30\n",
      "scipy: 1.8.0\n",
      "numpy: 1.22.3\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"refellips: {refellips.version.version}\\n\"\n",
    "    f\"refnx: {refnx.version.version}\\n\"\n",
    "    f\"scipy: {scipy.version.version}\\n\"\n",
    "    f\"numpy: {np.version.version}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading a dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*refellips* has the capability of loading data directly from output files of both Accurion EP3 and EP4 ellipsometers, as well as Horiba ellipsometers using `open_EP4file` and `open_HORIBAfile` respectively.\n",
    "\n",
    "Alternatively, other datasets can be imported using `DataSE`. The file must be formatted to contain four columns: wavelength, angle of incidence, psi and delta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "pth = os.path.dirname(refellips.__file__)\n",
    "dname = \"testData1_11nm_PNIPAM_on_Si_EP4.dat\"\n",
    "file_path = pjoin(pth, \"../\", \"demos\", dname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now use `DataSE` to import our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = DataSE(data=file_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating a model for our interface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with *refnx*, `ComponentSE` objects are assembled into a `StructureSE` object which describes the interface. The simplest of these `ComponentSE` objects is a `SlabSE`, which is what we will use here.\n",
    "\n",
    "We begin by loading in dispersion curves which describe the refractive index for each layer within our `StructureSE` (i.e., for each `ComponentSE`). *refellips* offers multiple ways to prescribe the refractive index of a layer. Here we will demonstrate the different ways to prescribe refractive indices. \n",
    "- `RI()`: Users can provide a refractive index ($n$) and extinction coefficient ($k$) for a given wavelength by `RI([n, k])` or alternatively, for spectroscopic analysis users load in their own dispersion curve of $n$,$k$ as a function of wavelength by providing a path to the desired file as `RI('my_materials\\material.csv')`.\n",
    "- `Cauchy()`: Creates a dispersion curve of a given material using the provided $a$, $b$ and $c$ values.\n",
    "- `load_material()`: Searches the *refellips* materials database for the provided material.\n",
    "\n",
    "Note, dispersion curves provided in the *refellips* materials database are downloaded as a .csv file from [refractiveindex.info](https://refractiveindex.info/). When providing a dispersion curve, files must contain at least two columns, assumed to be wavelength (in microns) and refractive index. If three columns are provided the third is loaded as the extinction coefficient. Further details on modelling material optical properties are provided on the [FAQ](https://refellips.readthedocs.io/en/latest/faq.html#where-do-i-find-dispersion-curves-for-a-material) page."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "si = load_material(\"silicon\")\n",
    "sio2 = RI([1.4563, 0])\n",
    "PNIPAM = Cauchy(1.47, 0.00495)\n",
    "air = RI(pjoin(pth, \"materials/air.csv\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have defined the refractive indices of our layers, we can create a `Slab` object for each interfacial layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is a 20 Angstrom layer\n",
    "silica_layer = sio2(20)\n",
    "\n",
    "polymer_layer = PNIPAM(200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each `Slab` has an associated thickness (as defined above), roughness, and volume fraction of solvent. As this is a dry film we will leave `vfsolv` as 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "silica_layer.name = \"Silica\"\n",
    "silica_layer.thick.setp(vary=True, bounds=(1, 30))\n",
    "silica_layer.vfsolv.setp(vary=False, value=0)\n",
    "\n",
    "polymer_layer.name = \"PNIPAM\"\n",
    "polymer_layer.thick.setp(vary=True, bounds=(100, 500))\n",
    "polymer_layer.vfsolv.setp(vary=False, value=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now create the `Structure` by assembling the `Component` objects. Our structure is defined from fronting to backing, where the thickness of the fronting and backing are defined to be infinite (i.e., `np.inf`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "structure = air() | polymer_layer | silica_layer | si()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can create our model. A wavelength must be provided here, however, if your ellipsometry dataset contains a wavelength that will be automatically used. We have the option to define the `delta_offset` parameter here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = ReflectModelSE(structure)\n",
    "\n",
    "model.delta_offset.setp(value=0, vary=False, bounds=(-10, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now have a quick preview of how our model compares to our dataset prior to fitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAAEGCAYAAAAubTHtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA20ElEQVR4nO3dd3hU1dbA4d8KoQVRxAJKSVApAkJAUFS6AtaLvWEviIICIqhwVVC5VwSxXkv8ROGKihULIk2qUi5iUJqIQDCIBQREkRJY3x/7xAwxDTJnzplkvc8zT2b2nDOzKJM1u5y1RVUxxhhjwiIh6ACMMcaYSJaYjDHGhIolJmOMMaFiickYY0yoWGIyxhgTKolBB1AUCQkJWrFixaDDMMaYuLJ9+3ZV1bjrgMRFYqpYsSJ//PFH0GEYY0xcEZE/g47hQMRdJjXGGFOyWWIyxhgTKpaYjDHGhEpczDEZY0xR7N69m8zMTHbs2BF0KDFVoUIFatasSdmyZYMOJSosMRljSozMzEwqV65MSkoKIhJ0ODGhqmzatInMzEzq1KkTdDhR4dtQnojUEpHpIrJMRJaKSG+vPVVE5olIuogsFJGT/Hj/sWMhJQUSEtzPsWP9eBdjTJjs2LGDww47rNQkJQAR4bDDDitRvUQ/55iygH6q2hBoBfQUkYbAo8AQVU0F7vceR9XYsdC9O2RkgKr72b27JSdjSoPSlJSyHeifWURGicjPIrIkom2c13FIF5G1IpIe8dy9IrJKRL4RkS7FjzxvviUmVd2gqou8+9uA5UANQIGDvcMOAX6I9nsPGgTbt+/btn27azfGGPOXV4AzIxtU9TJVTfU6D+8A7wJ4HYvLgUbeOc+KSBk/gorJqjwRSQGaAfOBPsBwEfkeGAHcm8853b2hvoVZWVn79X7r1uXXbntPGWP8VaZMGVJTU2ncuDGXXHIJ23N/S45w9tlns2XLltgFl4uqzgJ+zes5cd2wS4HXvaauwBuqulNV1wCrAF+mYnxPTCJyEC7r9lHV34Bbgb6qWgvoC7yU13mqmqaqLVS1RWLi/q3RqF07n3bJhKeegj/j8mJoY0wcqFixIunp6SxZsoRy5crx/PPP53vsxx9/TJUqVfwMJzH7C753674f57YBflLVb73HNYDvI57P9NqiztfEJCJlcUlprKq+6zVfi9c1BN7Ch4w79Ow5JLFvCaMk/mBo9aegd2849lhLUMYY31dJtWnThlWrVrFhwwbatm37V09q9uzZAKSkpLBx48aovmcuWdlf8L1b2n6cewU5vaWY8nNVnuB6Q8tVdWTEUz8A7bz7HYFvc59bXN0+voo0biaZtQh7SWYtadxMt7JvwfTpUL++JShjSjufV0llZWUxceJETjjhBF577TW6dOlCeno6ixcvJjU1NSrv4RcRSQQuBMZFNK8HakU8rum1RZ+q+nIDWuMWOnwFpHu3s732L4DFuDmnEwt7raSkJN0vIqruv9q+N5GcY6ZPV23f3rUfdZTqk0+qukq8xpg4tWzZsqIfnJyc9++J5ORixZCQkKBNmzbVpk2baq9evXTnzp06c+ZMPfbYY/WBBx7QL7/8MiKEZP3ll1+K9X7Z8vqzA39o4b+rU4AludrOBGbmamvk/d4uD9QBVgNlCnv9A7n5lpiiedvvxLQ//+EiE1T16qpPPGEJypg4tV+JqShfYA9ApUqV8mxfv369pqWladOmTXX06NGqGnxiwg3VbQB24+aMbvTaXwF65HH8IOA74BvgrIJeuzi3klkrb+hQSEraty0pybXn1r69G96bPh0aNIA+feCYY+DJJ22Iz5iSLN9VUvm0F0NGRgbVqlXj5ptv5qabbmLRokVRf48DoapXqOpRqlpWVWuq6kte+3Wq+rdVG6o6VFWPVdX6qjrRr7hKZmLq1g3S0iA5GUTcz7Q0156f7AQ1YwYcf7wlKGNKuv35AltMM2bMoGnTpjRr1oxx48bRu3fvqL9HSSJe9yzUKlWqpDHfKHDmTBgyxCWr6tXhnnvcxKjtpGtMaC1fvpzjjz++6CeMHeuuvF+3zvWUhg4t+AtsiOX1ZxeR7apaKaCQDljJ7DFFQ7t28Omn1oMypiTr1g3WroW9e93POE1KJY0lpsJYgjLGmJiyxFRU+SWoJ56wBGWMMVFkiWl/5U5QfftagjLGmCiyxHSgLEEZY4wvLDEVV3aCmjkTGja0BGWMMcVkiSla2raFadMsQRXCdhY2JV32theNGjWiadOmPPbYY+zdu7fAc9auXUvjxo0BSE9P5+OPP45FqKFliSnaLEHly3YWNqVB9rYXS5cuZcqUKUycOJEhQ4YU+XxLTJaY/GMJ6m9sZ2ETNn734I888kjS0tJ45plnUFX27NlD//79admyJU2aNOGFF17Y5/hdu3Zx//33M27cOFJTUxk3bhwLFizglFNOoVmzZpx66ql888030Q0yjPwqwhfN234XcQ2jmTNVO3Z0RSKrVVMdOVL1jz+CjiqmfKqZacxf9qeI66uvqiYl7ft/MSnJtRdHXkVcDznkEP3xxx/1hRde0IceekhVVXfs2KEnnniirl69WtesWaONGjVSVdWXX35Ze/bs+de5W7du1d27d6uq6pQpU/TCCy/M830PtLp4GG/WY4qVyB5Uo0Zw552uB/X443/vRpRQMayZaUyhgujBT548mTFjxpCamsrJJ5/Mpk2b+Pbbgrek27p1K5dccgmNGzemb9++LF261L8AQ8ISU6xlJ6hZs6Bx41KVoGJYM9OYQq1bt3/tB2r16tWUKVOGI488ElXl6aefJj09nfT0dNasWUPnzp0LPP++++6jQ4cOLFmyhA8//JAdO3ZEN8AQ8nMH21oiMl1ElonIUhHpHfHc7SKywmt/1K8YQq1NG5g6tVQlqAMp+m6MX2LRg//ll1/o0aMHvXr1QkTo0qULzz33HLt37wZg5cqV5C5QXblyZbZt2/bX461bt1KjRg0AXnnllegFF2Z+jRECRwHNvfuVgZVAQ6ADMBUo7z13ZGGvVSLmmAoza5bq6aeX6jkoY4orDHNM2TvYNmzYUJs0aaLDhw/XPXv2qKrqnj179N5779XGjRtro0aNtH379rply5Z95pg2bdqkLVq00KZNm+obb7yhn3/+udatW1dTU1N10KBBmpzPDrslaY4pZtteiMj7wDPAzUCaqk4t6rmBbHsRlNmz3XYb06ZBtWpw991wyy1/HwMzxvzN/m57UYJ2vbBtL/aXiKQAzYD5QD2gjYjMF5GZItIyn3O6i8hCEVmYlZUVizDDIb8hvpEjS+wQnzFBsV0vwsn3xCQiBwHvAH1U9TcgEagKtAL6A2+KiOQ+T1XTVLWFqrZITEz0O8zwyZ2g+vWDOnUsQRljSjxfE5OIlMUlpbGq+q7XnAm86w2BLgD2Aof7GUdci0xQJ5xgCcqYQsRqeiJMStqf2c9VeQK8BCxX1ZERT43HLYBAROoB5YCNfsVRYmQnqNmzoUkTS1DG5KFChQps2rSpxP2iLoiqsmnTJipUqBB0KFHj2+IHEWkNzAa+xvWKAAbiVuSNAlKBXcBdqvppQa91oIsftm+HChVcuZESZ84ct0hi6lQ48ki3SKJHD1skYUq13bt3k5mZWSqu9YlUoUIFatasSdmyZfdpj9fFDzFblVccB5qYBgyA0aOhUyfo3Nn9POooHwIMkiUoY0w+LDH56EAT08cfw2uvwZQp8PPPru2EE1yS6tzZjY5VrBjlYINiCcoYk4slJh8V9zqmvXvhq69g8mR3mz0bdu2C8uVdhaAuXVyiatzYVSSIa7kT1IABLkFVirv/m8aYYrLE5KNoX2C7fbtb5DZ5MkyaBMuWufbq1XN6U2ec4a5vjVuWoIwp9Swx+cjvyg+ZmW64b/Jk93PTJteempqTqFq3dj2suGMJyphSyxKTj2JZkmjvXvjyy5xhv88+g9273VxU+/Y5ier44+Ns2O+zz1yCmjLFEpQxpYQlJh8FWStv2za3hVJ2osrePLJGjX2H/Q6Pl0uELUEZU2pYYvJRmIq4ZmTkDPtNnQqbN7ue04kn5iSqU06BcuWCjrQQlqCMKfEsMfkoTIkp0p49sHBhTm9q7lzXVqkSdOiQk6jq1QvxsF9kgjriCJegbr3VEpQxJYAlJh+FNTHl9ttvMH16TqJatcq1167tElSXLtCxI1StGmycebIEZUyJY4nJR/GSmHJbvTonSU2b5hJXQgK0bJnTmzr5ZMhVRSRYlqCMKTEsMfkoXhNTpKwsWLAgJ1HNn+9WAFau7HpR2YnquOOCjtRjCcqYuGeJyUclITHltmULfPqpu8B30iS3qALcnoDZSapjRzjkkEDDhM8/dwlq8mRLUMbEGUtMPiqJiSmSqpuPyu5Nffop/P47lCnjhvqyE1XLlhDYnomWoIyJO5aYfFTSE1Nuu3fDvHk5JZMWLnTJq0oVOP30nESVkhJAcLkTVP/+cNttlqCMCSFLTD4qbYkpt02b3OKJ7ESVmena69bNSVIdOrj5qpixBGVM6BWWmERkFHAu8LOqNo5ovx3oCewBJqjqAK/9XuBGr/0OVZ3kS9w+bhRYCxgDVAMUSFPVJyOe7weMAI5Q1QJ3sC3tiSmSqqs+kT3sN326K0qbmOgu7M2ulN68uRsK9J0lKGNCqwiJqS3wOzAmOzGJSAdgEHCOqu4UkSNV9WcRaQi8DpwEHI3b9LWequ6Jdtx+7u2aBfRT1YZAK6Cn9wfLTlqdgXU+vn+JJAINGsAdd8BHH8Gvv7o5qbvugj/+gH/+E046yRVzuOwyeOkl+P57HwM69VTXjfvsM2jWzM091akDw4e7gEzojR3rhoUTEtzPsWODjsjEiqrOAn7N1Xwr8Iiq7vSO8XazoyvwhqruVNU1wCpckoo63xKTqm5Q1UXe/W3AcqCG9/TjwABcT8oUQ/nybhjv3/+GL75wGyK+9hr84x+usPhNN7kLfI8/Hnr3hgkT3MKKqLMEFZfGjoXu3d2qUFX3s3t3S04lSKKILIy4dS/COfWANiIyX0RmikhLr70GEPk1N5Oc3+lRFZM5JhFJAWYBjYEOQEdV7S0ia4EWNpTnD1W319SkSW6kbeZM2LHDXdB72mk5w36pqe7bclTNneuG+CZNchVuBwywIb4QSknJuVQhUnIyrF0b62hMtBVl8YP3+/mjiKG8JcB04A6gJTAOOAZ4Gpinqq96x70ETFTVt6Mdt59DeQCIyEHAO0Af3PDeQOD+IpzXPTvLZ2Vl+RtkCSUCjRrBnXfCJ5+4grNTpkCfPu7+vfe64rPVqsGVV8Irr8D69VF681NOcW/6+efuTQYMcL8FrQcVKuvyGUzPr92UCpnAu+osAPYChwPrgVoRx9X02qLO1x6TiJQFPgImqepIETkBmAZs9w6pCfwAnKSqP+b3OtZj8sePP+ZUSp882Q0Dgktm2bX92rSBpKQovFnuHtQDD7hq5oFdmGXAekwl3QH2mHoAR6vq/SJSD/c7uzbQEHiNnMUP04C6fix+8HNVngCjgV9VtU8+x6zFhvJCYe9e+PrrnCQ1ezbs3OnmsNq0yVmW3qRJMSulz50Lgwa55YSNG8NTT7lJMhOI7Dmm7dtz2pKSIC0NunULLi4THUVYlfc60B7XI/oJeAD4LzAKSAV2AXep6qfe8YOAG3CjX31UdaIvcfuYmFoDs4GvcV1BgIGq+nHEMWuxxBRK27e75JSdqJYsce3VqkGnTi5JdeoE1asfwIurwnvvQb9+7mv5RRfBiBEBXTFsxo513xXWrXMLZYYOtaRUUtgFtj6yxBS89etzhv2mTIGN3leJpk1zelOtW0OFCvvxon/+CY895pYU7t3rroG6+25bIGFMlFhi8pElpnDZuxfS03N6U3PmuDJKFSpAu3Y5iapRoyIO+33/vUtIr78ONWu6BRKXXRbi3RWNiQ+WmHxkiSncfv/dLUXPTlQrVrj2o4/OSVJnnOEKQxQ4bjRnjrty+Msv3cTWk0+6a6KMMQfEEpOPLDHFl3Xr9h3227zZtTdP2UTnzJfpnDWBtsyiDHv/PtO+Zw+MGgUDB7oigTffDA8/7GU1Y8z+sMTkI0tM8WvPHli0yLvI9+H5zN3ZnCzKchzf0o/HuJbRVEyu9ve1yZs3u+XlzzzjqtMOHuwu0A3Vdr/GhJslJh9ZYiohEhL4TQ/iY85mJHfyP07iCH7mdp7hto0PcthheZyzbJm7InjKFGjYEJ54wi0HNMYUKl4Tk++VH4z5S+3aHMw2Lmcc8zmZGbTjJBZwPw9SuzbcfjusWZPrnIYNXXdr/HhXT6lzZzj/fFi9OoA/gDEmFiwxmdgZOvSvMhICtGMWHyVdxpJHPuLSS+GFF+C44+Dyy11B2r+IQNeusHQp/OtfMHWqq0o7cKBPFWmNMUGyxGRip1s3t9AhOdklm+RkSEuj0d3n8vLLrrd0110wcSK0aAEdO7r7f402V6jgCvx98w1ceqm7/ql+fbfSLw6GpI0xRWNzTCZ0tm6FF19000nr17vKRf37u55UuXIRB86d68b/vvjCbbvx1FOuYKwxBrA5JmOi5pBDXM9p9WpX8Rzg2mvhmGNc5aLffvMOPOUUWLDA7Ya4ahW0bOk2oMquRmuMiUuWmExolSvnEtJXX8HHH0O9eq7nVKuWKxTxww+4jaRuuAFWrnT7e4weDXXrwsiRsGtX0H8EY8wBsKE8E1cWLnQVi95+G8qUcdNWd93lyh8Bbv6pb183OVW/vhsPPPPMIEM2JjA2lGdMDLRoAePGwbffwi23uPuNG8O557qySFqvvuteffSRK+p31llw3nnuBGNMXLDEZOLSMcfA00+78kdDhsD8+dC+PZx8sutN7TnzHLfB1LBhMGOG61LdfTds2xZ06MaYQlhiMnHt8MPh/vtdgnruOfj1V7jkEjeK9+xL5dnea4DrLXXrBo8+6iaqRo92vSljTChZYjIlQsWKbqf2b75xPabDD4eePd2lUkNeqM7G4S/DvHmuovl117nl5QsWBB22MSYPviUmEaklItNFZJmILBWR3l77cBFZISJfich7IlLFrxhM6VOmjNsQd+5cmDXLrSgfPNjlo55jTua7V+e6NegZGW7c77rrYMOGgKM2xkTyc2v1o4CjVHWRiFQGvgDOB2oCn6pqlogMA1DVuwt6LVuVZ4pj2TK3Ue5//+uqnV90EfS/7Q9aTnwQHn8cypd344G9e+e6gteY+Gar8nJR1Q2qusi7vw1YDtRQ1cmqmuUdNg+XqIzxTcOG7hrctWvddVCTJ8NJHSrRfv4wJjyzhr3tOsCAAW5534QJQYdrTKkXkzkmEUkBmgHzcz11AzAxn3O6i8hCEVmYlZWV1yHG7Jejj4ZHHnELJR57DL77Ds69pQZN1n7AK3d+xS4p79adn322m6wyxgTC9wtsReQgYCYwVFXfjWgfBLQALtRCgrChPOOH3bvhjTfcBbtffw1HH630bjabW2Z145A/f3TbvN9/v6uRZEwcsqG8PIhIWeAdYGyupHQdcC7QrbCkZIxfypaFq6+GxYvhk0+gQQPh7gltqaUZ9D/+QzJHvumWl48aZcvLjYkhPxc/CDAa+FVV+0S0nwmMBNqp6i9FeS3rMZlYWbTI9aDefBMSRLmy6kTu+uVuTmhRwVUvP+WUoEM0psjitcfkZ2JqDcwGvgayv24OBJ4CygObvLZ5qtqjoNeyxGRibc0at2DvpZeU7duFs8pPo//Oh2nfrSby6DA3YWVMyFli8pElJhOUTZtcRYmnn9rLz78kcKJ8Qf9yT3HRP48n8a4+bvNCY0LKEpOPLDGZoP35p7sOasS/d/Ht2nLUYTV3HjaG6//TgkqXnuN25DUmZCwx+cgSkwmLPXvggw9g+D+3MHdZFaqyiZ51JtJrdEuObFM/6PCM2YclJh9ZYjJh9NnMLB7tlcEHS46lAn9y3QmL6PfKCRzX/OCgQzMGiN/EZEVcjTlAp7VL5P2vj2X5Z79yVYMvGPV1C+qdeBAXn7ia+Z/vCTo8Y+KWJSZjiqnBqVV5cXlr1k5ayT01xzJt0aG0Oq0MbVO3/rVfoTGm6Gwoz5hoUmXbK+/wf32X8vjW6/me2hx/3G7uurcs3bq5erHGxEq8DuVZYjLGD3/8we5/DefN4RkMz+rLYm3CUdWV3n2EW26BKlWCDtCUBpaYfGSJycSttWvRfncx9d2tPFrhAabuaM1BBynduwt9+kCtWkEHaEoyS0w+ssRk4t6nn0Lv3ny5JJER1YYzbuPpiAiXX+624mjSJOgATUkUr4nJFj8YEwsdO8KXX9LsmZsYu+tSvtt7DD0bTue995SmTeHMM2HaNIiD74nG+M4SkzGxkpgIPXvCt9+SfOvZPLHkDNaVPY6h/5hPerpyxhnQooXbisO2IDMlgYjUFZG3RWSZiKzOvhV2niUmY2LtsMPgP/+BL7+kamptBn7QirVHnETaXd/w++9wxRVQt64rZm4j2MZPIjJKRH4WkSURbYNFZL2IpHu3syOeu1dEVonINyLSpQhv8TLwHJAFdADGAK8WdpIlJmOC0qSJm3t66y0q/PYzN49owPImlzE+7Wdq1IDevd3iiPvug59+CjpYU0K9ApyZR/vjqprq3T4GEJGGwOVAI++cZ0WkTCGvX1FVp+HWM2So6mDgnMKCssRkTJBE4OKLYcUKGDKEhAkf0vWOZOacMZjPpu2gXTsYOhSSk+GWW2DlyqADNiWJqs4Cfi3i4V2BN1R1p6quAVYBJxVyzk4RSQC+FZFeInIBcFBhb2SJyZgwqFjRbeO+YgV07QpDhnDq9fV574o3Wb5MueYaGD0aGjSACy+EuXODDtjEiUQRWRhx617E83qJyFfeUN+hXlsN4PuIYzK9toL0BpKAO4ATgauAawp7c98Sk4jUEpHp3qTXUhHp7bVXFZEpIvKt9/PQwl7LmFKjdm23+mHGDDj0ULjsMur36EBaz8VkZMDAge6pU0+F1q3h/fet5JEpUJaqtoi4pRXhnOeAY4FUYAPwWDHeP0VVf1fVTFW9XlUvAmoXdlKBiUlEDvZ+Vs3rVshrZwH9VLUh0Aro6Y1R3gNMU9W6wDTvsTEmUrt28MUXbpfCJUugeXOqDb6Vh/tsZN06eOIJyMyE88+Hhg3h//4PduwIOmhTEqjqT6q6R1X3Ai+SM1y3Hoi8JLym11aQe4vYto/CekyveT+/ABZ6P7+IeJwvVd2gqou8+9uA5bhuX1dgtHfYaOD8woI0plQqUwZ69IBvv3XLzF98EerW5aCXn6Z3zyxWrYLXXoOkJLj5ZkhJgX/9CzZvDjpwE89E5KiIhxcA2Sv2PgAuF5HyIlIHqAssyOc1zhKRp4EaIvJUxO0VXKel4BhiUflBRFKAWUBjYJ2qVvHaBdic/TjXOd2B7gDlypU7cefOnb7HaUyoLVkCffq4K3EbNYInn4TTT0fVNQ0fDpMnQ6VKLlH16eMWTZjSq7DKDyLyOtAeOBz4CXjAe5wKKLAWuEVVN3jHDwJuwCWXPqo6MZ/XbQo0A4YA90c8tQ2YrqoFfn0qUmISkdOAdFX9Q0SuApoDT6jquiKcexAwExiqqu+KyJbIRCQim1W1wHkmK0lkjEcVxo+Hfv1gzRq44AJ47DGoUweAxYthxAh4/XV3+GWXuZJHqamBRWwCFHRJIhFJVNX9vly8qInpK6Ap0AS37v3/gEtVtV0h55UFPgImqepIr+0boL2qbvC6jDNUtcA9qS0xGZPLjh0uIf3rX26/9/794Z57XHcJ/pqHevFF+P136NTJHXLGGW6FuikdgkpMIvI1rseVJ1UtsDpkURPTIlVtLiL3A+tV9aXstgLOEdwc0q+q2ieifTiwSVUfEZF7gKqqOqCg97fEZEw+MjPh7rvdZFPNmvDoo3D55X9ln82b4YUX3Kjfjz+6nlP//nDJJVC2bLChG/8FmJgKHERW1YwCzy9iYpoJfAJcD7QFfgYWq+oJBZzTGpgNfA1kL2gdCMwH3sQtGczA9bwKvMDLEpMxhfjsM7jjDli0CE47zdUzap7zvXHnTnj1VTfMt2KFW5Xety/cdBMcVOjljiZeBT2U58WQDNRV1akiUhFI9BbE5X9OERNTdeBK4H+qOltEauOG48ZEI/DCWGIypgj27IGXX3YXO23c6LLO0KFwxBF/HbJ3L0yY4DpWc+a4S6VuvRVuvx2qVw8wduOLoBOTiNyMW8RWVVWPFZG6wPOqenqB5xWUmESkAtADOA7X83npQCayissSkzH7YcsWePBBePppN+c0eLBbbp5r7G7ePLeS77333FN33AGPPOJWqZuSIQSJKR13HdR8VW3mtX1d0GgbFH4d02igBS4pnUXxrgA2xsRClSowciR89RW0auXG7Jo2dWvJI7RqBe+844b2rrzSDfNdcYUb9jMmSnaq6q7sByKSSAGLIrIVlpgaqupVqvoCcDHQpngxGmNi5vjjYeJE+OAD2LULunRxdfi++26fw+rVcyOAI0bAW2/Beee5lXzGRMFMERkIVBSRTsBbwIeFnVRYYtqdfSeIITxjTDGJuEyzdKkbp/v0U1fD6N57/5Z9+vWDUaPcxbqdOsGvRa05bUz+7gF+wY263QJ8DPyzsJMKm2PaA2RP7ghQEdju3VdVPbh4MReNzTEZEyU//OCS0pgxcNRRMGwYdOsGCTnfUcePdxfm1q0LkyZBjcLqR5vQCnqOyYvhCABV/aWo5xTYY1LVMqp6sHerrKqJEfdjkpSMMVF09NFu/4y5c911T9dc45aX/+9/fx1y/vnwySeQkeEqmK9aFVy4Jj6JM1hENgLfAN+IyC/etbCFsv2YjCmNWrVyy/JGjXKljU46CW680S0zBzp0gOnT3Whf69au1JEx+6EvcBrQUlWrqmpV4GTgNBHpW9jJMSniWlw2lGeMj377DR56yNUwOvhgNxd1442QkMCKFW6+ads2+PBDaGPLn+JKgJUfvgQ6qerGXO1HAJOzl47nx3pMxpR2Bx/sLmhKT4fGjaF7dze8l55OgwauqET16tC5s7s415giKJs7KcFf80yFFsOyxGSMcRo1ctvjjhnjlpSfeCL06UPtKr8xe7Z7+vzzYezYoAM1cWDXAT4H2FCeMSYvmzfDoEHw/POuuzRyJL+ddRldzxdmzHBFJXr1CjpIU5gAh/IiV3Tv8xRQQVUL7DVZj8kY83eHHgrPPgvz57uVfFdcwcEXdWLikyvp2tXV1hsyxG0PZUxuuVZ0H5xrdbcN5RljiqFlS5ecnnkGFi6kQssTeLvRA1x3VRaDB0Pv3q4wrDHRZEN5xpii+fFHt5nTq6+yN+UY+qdOYeT4Y7jqKrfq3PZ3Cp8wXGB7IKzHZIwpmurV4b//henTSahYnhHjj2Xo8a/y6qtw4YXw559BB2hKCt8Sk4iMEpGfRWRJRFuqiMwTkXQRWSgiJ/n1/sYYn7RvD+npyCOPMHBtd54rewcTJihdOu9ly5aggys5xo6FlBRXLSolpXSthvRtKE9E2gK/A2NUtbHXNhl4XFUnisjZwABVbV/Ya9lQnjEhlZEBvXsz7v3yXM2rNDpuB5/MqUy1akEHFt/GjnWXk23fntOWlARpaa60YVHZUF4uqjoLyF2fWIHsGnuHAD/49f7GmBhITobx47nsg6v4sNpNrFyVQOu6P7L2f0Wu12nyMGjQvkkJ3ONBg4KJJ9ZiPcfUBxguIt8DI4B78ztQRLp7w30Ls7Jsxw1jfBONMaPzzqPL6ueYeu2rbNxWntNO3s3S+95w272b/bZu3f61lzSxTky3An1VtRauyN9L+R2oqmmq2kJVWyQmJsYsQGNKlewxo4wMd1FSRoZ7fCDJKSmJU165hVnvb0HLlqPtw52Y3/hGWLgw+nGXcLVr7197SePrcnERSQE+iphj2gpUUVUVEQG2FmX7DJtjMsYnKSkuGeWWnAxr1x7wy67+Tul82u/8+HMC7+kFdLqtLgwd6rZ9N4WyOabY+gFo593vCHwb4/c3xkTyaczomGOFOemVObZhBc5J+Ji3n/sF6td3y83j4NrJoHXr5pJQcrLbhDg5ef+TUjzzc1Xe60B74HDgJ+AB3IZRTwKJwA7gNlX9orDXsh6TMT7xqceUbfNmt7P73LnK87X/zc1rB0G7dq7cUcOGxX59U7B47TFZ5QdjSrNojRkVYPt2uPhimDgRHrlgPnfPOMtt8NSvH9x3H1SKu9+bcSNeE5NVfjCmNIvBmFFSEowfD1dcAfe8dzIDrsxEr7oahg1zvabx4214z+zDekzGmJjYu9dVJX/2WbjhBnjh6jkk3n4rLFkC554LTz0FdeoEHWaJYj0mY4wpQEKCK1J+//2u6OtlT7dmx+eLYMQImD7d9Z6GDoWdO4MO1QTMekzGmJh78kno0wc6dnQjeZW3ZkLfvvD221CvnutWnX560GHGPesxGWNMEfXu7XZwnznT5Z+NFWrCW2+5FRJ79sAZZ8CVV8KGDUGHagJgickYE4irr4b33oOvvoK2bSEzEzjzTDfn9MAD8O677tqnp54CK0tWqlhiMsYE5rzzYNIkWL8eTjsNVq4EKlSAwYPh66/hlFNc96plS5g3L+hwTYxYYjLGBKpdO5gxw2002Lo1LFrkPVG3LnzyiRvi++UXl6S6d4dNm4IM18SAJSZjTOCaNYM5c6BiRbcP4cyZ3hMi7urc5cvhzjvdcr4GDeDll936c1MiWWIyxoRCvXrw2WdQsyZ06QIffBDxZOXK8NhjrjtVr567EKptWzdBZQ5YXjuNRzzXT0RURA73HouIPCUiq0TkKxFp7ldclpiMMaFRsybMmgVNmsCFF7qVe/to0gRmz3Y9pxUroHlzV9po27ZA4i0BXgHOzN0oIrWAzkBkNd+zgLrerTvwnF9BWWIyxoTK4YfDtGluSO/aa+GJJ3IdkJAA118P33zjek4jR8Lxx7troOLguswwyWencYDHgQG4XcezdQXGqDMPqCIiR/kRlyUmY0zoVK4MEya4XlPfvq7W699yzmGHubp+c+fCEUfAJZfAWWfBqlWBxBxSidk7gXu37oWdICJdgfWqujjXUzWA7yMeZ3ptUWeJyRgTSuXLw5tvwo03wsMPQ69e+ax3aNUK/vc/V07i88+hcWO33HzHjliHHEZZ2TuBe7e0gg4WkSRgIHB/bMLLmyUmY0xolSkDL74IAwa4KkXdusGuXXkcmJgId9zh5p0uuACGDHEJ6pNPYh5znDsWqAMsFpG1QE1gkYhUB9YDtSKOrem1RZ0lJmNMqIm4HTKGDYM33oDzz993+6h9HH00vP46TJnistpZZ7khvszMWIYct1T1a1U9UlVTVDUFN1zXXFV/BD4ArvFW57UCtqqqLzWjfEtM+S1DFJHbRWSFiCwVkUf9en9jTMkyYIDrPU2aBJ06ud1x83XGGW4p+cMPw0cfucURI0fC7t0xizceeDuNzwXqi0imiNxYwOEfA6uBVcCLwG2+xeXj1uptgd9xqzgae20dgEHAOaq6U0SOVNWfC3stqy5ujMn2zjuuvmv9+i5JHVXYurA1a9xGUBMmwAknuDHB1q1jEmvQrLp4LvksQ7wVeERVd3rHFJqUjDEm0kUXuRyzerXLL6tXF3JCnTrw4YeuYuyWLdCmjVtm/ssvsQjXHIBYzzHVA9qIyHwRmSkiLfM7UES6Zy9xzLLKwsaYCGec4a512rLFFX/9+utCThBxk1PLl8Pdd8N//+u6XGlpVtoohGKdmBKBqkAroD/wpohIXgeqalr2EsfExMRYxmiMiQMnn+yqRCQkuOpEn39ehJMqVYJHHoHFi10ViVtugVNPhS+/9D1eU3SxTkyZwLvelcMLgL3A4TGOwRhTQjRq5OrrHXGEWxAxaVIRT2zY0G3nPmaMm4Nq0cJtr7F1q6/xmqKJdWIaD3QAEJF6QDlgY4xjMMaUICkprnxevXpuf6dx44p4oojbrXDFCujRA55+2lUuf/11K20UMD+Xi+e1DHEUcIy3hPwN4Fr1a1mgMabUqFbN7enUqhVccQU8//x+nHzoofCf/8CCBVCjhlvy16mTq8VnAuHbcvFosuXixpii+PNPuPRSd+nSww/DwIGuY1Rke/bACy+4E7dvdxdPDRwISUm+xewnWy5ujDEBq1gR3n0XrroK/vlPtyPGfi26K1MGbrvN9ZYuvxyGDnUTWR995FvM5u8sMRljSpSyZWH0aFc67/HH3SVL+33FSbVqbmHEjBku2513nqvBl5HhR8gmF0tMxpgSJyHB7eP04IMuSV188QEWG2/XDtLTXaG+yZPdar5hw/KpJGuixeaYjDEl2n/+47bMaN8e3n8fDj74AF8oIwP69IHx413tvWefdS8aYjbHZIwxIdSzJ4wdC3PmQMeOxahElJzsyhp9+KFbZdGhg1tu/tNPUY3XWGIyxpQCV17pekvLlrn6euvWFePFzj0Xli51qyvGjXOljZ591q3oM1FhickYUyqcfbabJvrpJ1dfb8WKYrxYUhI89JAr0teiheuWnXyy20nXFJslJmNMqdG6Ncyc6bZlat0aFi4s5gvWr+82JXz9dVi/3iWn224rZLMoUxhLTMaYUqVpUzffVLmymyb69NNivqCIu+ZpxQq379MLL7iENWaMlTY6QJaYjDGlznHHueKvyclu9/Xx46PwooccAk8+6bphxxwD117rVu0tXRqFFy9dLDEZY0qlo49222Y0b+42H3z55Si9cLNmbg+OtDQ3B5Wa6vaAskteiswSkzGm1KpaFaZOdRsP3nADPPZYlF44IQFuvtmVNrrmGnj0UXft0/jxNrxXBJaYjDGlWqVK8MEHcMklcNddrmZr1HLHEUfASy+5Sa0qVVxZo/POK8J+8KWbJSZjTKlXvrxbWHfLLfDvf7vtmaJ6WdJpp8EXX7gu2cyZrjDsww/Dzp1RfJOSwxKTMcbgCos/95zrMaWluX2doloSr2xZuPNOWL7c9Zruu89t7z51ahTfpGTwc6PAUSLys7cpYO7n+omIiohtq26MCQ0Rt9PFiBHw1lsuf/z+e5TfpGZNePNN+OQTtydHp05uufkPP0T5jeKXnz2mV4AzczeKSC2gM1CcoiDGGOObfv1g1CjXmenUCX791Yc36dLFrdobPNgtimjQwC033+89Okoe3xKTqs4C8vrnfBwYANjSFGNMaF1/PbzzDixaBG3b+tShqVABHngAliyBU0911ctbtIC5c314s/gR0zkmEekKrFfVxUU4truILBSRhVn2DcIYE4Dzz3cjbhkZbv3CqlU+vdFxx8HEifD227Bxo0tSN98Mmzb59Ibh5ut+TCKSAnykqo1FJAmYDnRW1a0ishZooaobC3sd24/JGBOkhQtdhYgyZWDSJFfWyDfbtsGQIW6nwypV3GRXhw4H9FK2H1PhjgXqAIu9pFQTWCQi1WMYgzHG7LcWLWD2bLewrl07d1mSbypXdqsvvvwSWrZ0dfdKmZj1mPJ4bi3WYzLGxJF166BzZ/fz7bfdVhphZj2mXETkdWAuUF9EMkXkRr/eyxhjYqF2bddzatgQunaF114LOqKSydceU7RYj8kYEya//eYS08yZ8PTTbp/AMLIekzHGlBIHH+wW0f3jH9CrFzz4oNVmjSZLTMYYcwAqVHDzTNdd5y5F6tPHFXIwxZcYdADGGBOvEhNd8fBDD4XHH3cVIkaNcqv3zIGzxGSMMcWQkOCKhh9+OAwaBFu2uFJ4FSsGHVn8sqE8Y4wpJhFXlfy552DCBFcGb+vWoKOKX5aYjDEmSnr0cPs6zZsH7dvDTz8FHVF8ssRkjDFRdNll8OGHsHIltG4Na9cGHVH+8tqeSEQeEpGvRCRdRCaLyNFeu4jIUyKyynu+uV9xWWIyxpgo69LFbZmxcaMr/rpsWdAR5esV/r490XBVbaKqqcBHwP1e+1lAXe/WHXjOr6AsMRljjA9OOQVmzXLXN7VpA/PnBx3R3+W1PZGq/hbxsBI5WxR1BcaoMw+oIiJH+RGXJSZjjPHJCSe4gq+HHgqnnx7ILuqJ2dsHebfuRTlJRIaKyPdAN3J6TDWA7yMOy/Taos4SkzHG+OiYY1xyOvZYOOcct/lgDGWpaouIW1pRTlLVQapaCxgL9PI3xL+zxGSMMT6rXh1mzHC7WFx6Kbz4YtARFdlY4CLv/nqgVsRzNb22qLPEZIwxMXDooTB5slsY0b07DBsWdER5E5G6EQ+7Aiu8+x8A13ir81oBW1V1gx8xWOUHY4yJkaQkGD/e1de75x63c/qwYe4C3SB42xO1Bw4XkUzgAeBsEakP7AUygB7e4R8DZwOrgO3A9b7FZdteGGNMbO3dC7ffDs8+CzfeCM8/7+ruRVu8bnvhW49JREYB5wI/Z+9gKyLDgfOAXcB3wPWqusWvGIwxJowSEuCZZ+Cww+Chh2DzZrfpYPnyQUcWDn7OMb3C3y/cmgI0VtUmwErgXh/f3xhjQkvE7eP0xBPw7rtuxd62bUFHFQ6+JaZ8LtyarKpZ3sN5uFUdxhhTavXuDWPGuFV7p5/uqkWUdkGuyrsBmJjfkyLSPfuisKysrPwOM8aYuHf11fDee/DVV9C2LWRmBh1RsAJJTCIyCMjCrZHPk6qmZV8UlujHrKAxxoTIeefBpEmwfr2rr7dyxAeQkuImpFJSYGy+vy5LnJgnJhG5DrcoopvGw5JAY4yJkXbt3JDen1t20Lr/KSzKqOqK7WVkuIufSklyimliEpEzgQHAP1R1eyzf2xhj4kGzZjDnoLOoyHY6MJ3PONU9sX272yK3FPAtMXkXbs0F6otIpojcCDwDVAameHt9PO/X+xtjTLyqt2Emn3EaJzOfmkRMOK1bF1xQMWQX2BpjTNikpLjhu9ySk/dr58F4vcDWauUZY0zYDB3q6hdFSkpy7aWAJSZjjAmbbt0gLc31kETcz7Q0114K2FCeMcaUUDaUZ4wxxkSBJSZjjDGhYonJGGNMqFhiMsYYEyqWmIwxxoRKXKzKE5G9wJ9Bx+FJxBWgDTOLsfjCHh+EP8awxwclP8aKqhp3HZC4SExhIiILVbVF0HEUxGIsvrDHB+GPMezxgcUYVnGXSY0xxpRslpiMMcaEiiWm/ZcWdABFYDEWX9jjg/DHGPb4wGIMJZtjMsYYEyrWYzLGGBMqlpiMMcaEiiWmQojIWhH52ttxd2FE++0iskJElorIo2GLUURSRWRedpuInBRgfFVE5G3v72u5iJwiIlVFZIqIfOv9PDSo+AqIcbj3+CsReU9EqoQpvojn+omIisjhQcVXUIxh+azk828cps9JfS+O7NtvItInbJ+VmFBVuxVwA9YCh+dq6wBMBcp7j48MYYyTgbO8+2cDMwKMbzRwk3e/HFAFeBS4x2u7BxgW8N9hXjF2BhK9tmFBxphXfN79WsAkICP3/4EwxBimz0o+8YXmc5Ir1jLAj0By2D4rsbhZj+nA3Ao8oqo7AVT154DjyYsCB3v3DwF+CCIIETkEaAu8BKCqu1R1C9AV94sC7+f5QcQH+ceoqpNVNfuK+3lAzTDF5z39ODAA9+8dmAJiDMVnpYD4QvE5ycPpwHeqmkGIPiuxYompcApMFpEvRKS711YPaCMi80Vkpoi0DDA+yDvGPsBwEfkeGAHcG1BsdYBfgJdF5EsR+T8RqQRUU9UN3jE/AtUCig/yjzHSDcDE2IcG5BOfiHQF1qvq4oDiipTf32FYPiv5xdeHcHxOcrsceN27H6bPSkxYYipca1VtDpwF9BSRtrjaVVWBVkB/4E0RkZDFeCvQV1VrAX3xvikGIBFoDjynqs2AP3DDEX9RN0YR5Df+AmMUkUG4WmVjgwkvz/gGAwOB+wOKKbf8/g7D8lnJL76wfE7+IiLlgH8Ab+V+LgSflZiwxFQIVV3v/fwZeA84CcgE3lVnAbAXCGziOZ8YrwXe9Q55y2sLQiaQqarzvcdv435B/CQiRwF4P4McDs0vRkTkOuBcoJv3SyEI+cVXB1gsImtxw4yLRKR6MCHmG2NYPiv5xReWz0mks4BFqvqT9zhMn5WYsMRUAG+4pHL2fdxk+BJgPG5SFxGph5tI3RiyGH8A2nmHdQS+DSI+Vf0R+F5E6ntNpwPLgA9wvxTwfr4fQHhA/jGKyJm4+Zt/qOr2kMW3SFWPVNUUVU3B/eJt7h0blhiXEZLPSgHxheJzkssV5AzjQYg+K7FilR8KICLH4Hog4IYCXlPVoV5XexSQCuwC7lLVT0MWY2vgSa9tB3Cbqn4RUIypwP/hfimtBq7HfSl6E6iNW1F2qar+GkR8BcT4P6A8sMk7bJ6q9ghLfKq6OeL5tUALVQ3kC5IXQyp//zv8g/B8VvKKrxEh+Zx4MVYC1gHHqOpWr+0wQvRZiQVLTMYYY0LFhvKMMcaEiiUmY4wxoWKJyRhjTKhYYjLGGBMqlpiMMcaEiiUmYyKIyPlepe4Gudpbi8gCrzL1iojST4jIYBG5q5DXTRCRMSLyuYg08it+Y0oCS0zG7OsKYI73EwCvmsJrQA9VbQC0Bm4RkXP243U7A/OBC4B+0QvXmJLHEpMxHhE5CJd0bsQV0czWE3hFVRcBeBexDiBXzb9ClMGV49kLBFlX0ZjQs8RkTI6uwCequhLYJCIneu2NgNzVABZ67UU1CVf65gNgZHEDNaYkSww6AGNC5ApceRqAN7zHUSlP4+3rdHmhBxpjLDEZAyAiVXFFPE8QEcUNvamI9McV+zyRfYtnnggsjXmgxpQCNpRnjHMx8F9VTfYqdtcC1gBtgP8A13lFQLOLag7DbXm9DxHpJSK9Yhe2MSWPJSZjnCvIqdKe7R3gCm/30KuAF0VkBfA5MEpVP8zjdRqQU43cGHMArLq4MVEkIh8BF6rqrqBjMSZeWWIyxhgTKjaUZ4wxJlQsMRljjAkVS0zGGGNCxRKTMcaYULHEZIwxJlQsMRljjAmV/wfgmwfuMKTxyAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "axt = ax.twinx()\n",
    "\n",
    "aois = np.linspace(50, 75, 100)\n",
    "\n",
    "for dat in data.unique_wavelength_data():\n",
    "    wavelength, aois, psi_d, delta_d = dat\n",
    "    wavelength_aois = np.c_[np.ones_like(aois) * wavelength, aois]\n",
    "\n",
    "    psi, delta = model(wavelength_aois)\n",
    "    ax.plot(aois, psi, color=\"r\")\n",
    "    p = ax.scatter(data.aoi, data.psi, color=\"r\")\n",
    "\n",
    "    axt.plot(aois, delta, color=\"b\")\n",
    "    d = axt.scatter(data.aoi, data.delta, color=\"b\")\n",
    "\n",
    "ax.legend(handles=[p, d], labels=[\"Psi\", \"Delta\"])\n",
    "ax.set(ylabel=\"Psi\", xlabel=\"AOI, °\")\n",
    "axt.set(ylabel=\"Delta\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating an objective"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now create an objective. The `Objective` object is made by combining the model and the data, and is used to calculate statistics during the fitting process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "objective = ObjectiveSE(model, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimisation of our `Objective` is performed by *refnx*'s `CurveFitter`. Data can be fit using a local optimisation such as `least_squares`, or a more global optimisation technique such as `differential_evolution`. For more information on the available fitting methods, see [*refnx*](https://refnx.readthedocs.io/en/latest/refnx.analysis.html#refnx.analysis.CurveFitter.fit)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitter = CurveFitter(objective)\n",
    "fitter.fit(method=\"least_squares\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimised model and data post fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now view our optimised objective, including our fit parameters. Users can plot this using the above method, or alternatively use the `objectiveSE.plot()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEGCAYAAADBr1rTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7+ElEQVR4nO3deZzP9fbA8dcxYw3JkhaaoagQY0nupaJFKkXaUzdSk9K+uElFuSq07ynSItUtrb8WRDtKsitKlilFukmJDOf3x/lMpjG7+Xw/3+93zvPx+D5m5vPdznTvd473+33e5y2qinPOORdPKkQdgHPOOZeXJyfnnHNxx5OTc865uOPJyTnnXNzx5OSccy7upEYdQHFUqFBBq1atGnUYzjmXUDZu3KiqmpCDkIRITlWrVuX333+POgznnEsoIvJH1DGUVkJmVOecc8nNk5Nzzrm4E1pyEpGGIjJNRBaJyEIRuTy4niEiM0RkjojMEpH2YcXgnHMuMYW55pQNXK2qs0WkBvC5iEwGRgI3q+pbInJc8HPnEONwzpUTW7ZsISsri02bNkUdSkxVqVKFBg0aULFixahDKTOhJSdVXQ2sDr7fICKLgb0BBWoGD9sV+D6sGJxz5UtWVhY1atQgPT0dEYk6nJhQVdatW0dWVhaNGjWKOpwyE5M1JxFJB1oDM4ErgFEisgq4AxgUxnuOHw/p6VChgn0dPz6Md3HOxZNNmzZRp06dcpOYAESEOnXqJN1oMfTkJCLVgZeAK1T1V+Ai4EpVbQhcCYwp4HmZwZrUrOzs7BK95/jxkHleNitWgCqsWGE/e4JyLvmVp8SUIxl/ZwnzyAwRqQi8AbyjqncF19YDtVRVxf6LrlfVmoW9zi677KIl2eeUXvc3VqyrvsP1tDq/sfynHa8755LD4sWLOfDAA6MOIxL5/e4islFVdynoOSIyFugOrFHVFsG154H9g4fUAn5R1YzgvkFAP2ArcJmqvlPGv8ZfwqzWE2xUtDgnMQW+Bw4Pvj8CWFrW771yXbUSXXfOubKSkpJCRkYGLVq04NRTT2Xjxo0FPva4447jl19+iV1wOxoHdMt9QVVPV9WMICG9BEwEEJFmwBlA8+A5D4lISliBhTmt1xE4BzgiKBufE1TnXQDcKSJzgVuBzLJ+431Yme/13fgfcy4ezdZ1v5T1WzrnHGAdbebMmcOCBQuoVKkSjzzySIGPffPNN6lVq1bsgstDVT8Afs7vvmCAcRowIbjUA3hOVTer6rfA10BoW4FCS06q+pGqiqq2zMnCqvpmcL2tqrZS1UNU9fOyfu/hde6iGn+fBhS28TN1aP1wJrXqpnB0o68ZOnAjU6bAb7+VdQTOuYQQcuXUoYceytdff83q1as57LDD/hpRffjhhwCkp6fz008/lel75pGas3Yf3EoyGDgU+FFVc2a39gZW5bo/K7gWioTorVdSve89BPpewuAtQ1jJPuzDSv6TejOdRvXgkz9a88nYxXz09R7cMqoKOgpSUpRWrYROnfjrtueeUf8WzrlQjR8PmZmQM+22YoX9DNC7906/fHZ2Nm+99RbdunXj2Wef5ZhjjmHw4MFs3bq10Km+Mpatqu1K+dwz2T5qirmkTE707k1voPfgzrByJeyzDwwfDr17kg6cNSgNFi9m/ZALmfFiFh/LYXy05hQeG70v991ng8n99oNDD91+23dfSMKCGOfKr8GDtyemHBs32vWdSE5//PEHGRkZgI2c+vXrx4wZMzjvvPPYsmULPXv2/Ov+eCUiqUAvoG2uy98BDXP93CC4Fg5VjftbtWrVNDRLlqj26aOakqJ/VtpFZ548Qu+44X/ao4dq7dqqVoyuuueeqqefrvrQQ6oLF6pu2xZeSM650lm0aFHxHyyy/QOe+yayUzHssssu+V7/7rvvdPTo0dqqVSt98sknVVU1LS1N165du1PvlyO/3x34XYv4+wqkAwvyXOsGvJ/nWnNgLlAZaAQsA1KKev3S3rzxa5Mm8MQTsGQJFf91Ju1fHczVI3bnlT0vYu3nK1mwAB5+GDp3ho8+gosvhubNoX59OOUUuO8+mDcPtm2L+hdxzpXIPvuU7PpOWLFiBfXr1+eCCy7g/PPPZ/bs2WX+HqUhIhOA6cD+IpIlIv2Cu84gz5Seqi4EXgAWAW8DA1R1a2jBhZX1yvIW6sgpr+XLVS+8ULViRbtlZto1tdHS11+rjhmjeu65qunp2/+xVaeOaq9eqvfdpzp/vurWrbEL2TlnSjRyeuYZ1WrV/j5qqlbNru+E/EZO48aN0+bNm2tGRoZ26tRJly1bpqrxMXKK11uom3DLSkk34ZaJlSvh9tthzBgbFvXpY3PR6el/e9iKFfD++zBtGrz3Hixfbtfr1oXDD4cjj4QjjoCmTX3NyrmwlXgT7vjx9rn+29r0zhdDRKE0m3DjmSenoqxaBSNGwGOPWZLq2xeuv36HJJVj+XJLUtOmwdSpkJVl1/fe25JUzi2EmQPnyj3vEOHJKaYiTU45vvsObrvNkpSqJanBgwvNMqrw9deWpHJuOVsamjSBo46yW5cusNtuMfo9nEtinpw8OcVUXCSnHFlZlqQef9yyT79+NpJq2LDIp27bBgsXwpQp8O67NsL6/Xfb/9e2rSWqo4+Gf/4TKlcO/1dxLtl4cvLkFFNxlZxyrFq1PUmJwAUXwKBBNn9XTH/+CZ9+aslqyhSYMQO2boVddrHqwK5d7bb//r5e5VxxeHLy5BRTcZmccqxcCbfeaoUTKSnQvz9cdx3ssUeJX+rXX200NWmS3ZYGTUMaNIBjjoFu3Wx0FWErLufimicnT04xFdfJKce338J//gNPPgmVKsGAATBwINSrt1MvOXmyJaopU2D9est/hxxiiapbN5sOrOC71ZwDPDklU3LyP2tlpVEjGz19+SWceircdZddGzwYfs636W+xXjIzE1580QopPvrIZg7//BOGDIH27W0z8DnnwIQJsG5dGf9OzrkSyzkyo3nz5rRq1Yo777yTbUXs0l++fDktWrQAYM6cObz55puxCDWueXIqa/vtZ6OnhQvhhBNsyq9RI7jlFpu3K6XUVOjYEYYNg88+gx9/hGeesdHTW2/BWWfB7rvbY4YPhy++sHoN51xs5RyZsXDhQiZPnsxbb73FzTffXOzne3IKRL0LuDi3mHaIKGvz5qmedJLtPq9dW3XECNXffivTt8jOVp0+XfXGG1Xbtt2+2X3vva3BxWuvqf7+e5m+pXNxqUQdItSaQaSlWTu9tLSdbg6hqjt2iPjmm2+0du3aum3bNs3OztZrrrlG27VrpwcddJA+8sgjqqr67bffavPmzXXz5s3asGFDrVu3rrZq1Uqfe+45nTlzpnbo0EEzMjL0H//4h3755Zf5vm+ydYiIPIDi3BI6OeX47DPVY4+1/+T166vee6/qpk2hvNXq1apjx6qefLJq9er2llWqqB53nDWuXbUqlLd1LnIlSU4hdS/Kt33Rrrvuqj/88IM++uijOmzYMFVV3bRpk7Zt21aXLVv2V3JSVX3iiSd0wIABfz13/fr1umXLFlVVnTx5svbq1Svf90225BTmMe0NRWSaiCwSkYUicnmu+y4VkS+D6yPDiiGutGsHb74JH38MBx4Il19uPY3GjIHs7DJ9qz32sD3COWtVkyfDhRfactjFF9uWrHbtbKZx7lyf/nPlU2EnZoRl0qRJPPXUU2RkZHDIIYewbt06luaU5RZg/fr1nHrqqbRo0YIrr7yShQsXhhdgHAlzzSkbuFpVmwEdgAEi0kxEumDH/bZS1ebAHSHGEH/++U9rFTF5smWR88+HZs2soiGE1uaVK1v5+T33WLeKhQtte1alSjB0KGRk2JLYZZfZxuAtW8o8BOfi0sqVJbteWsuWLSMlJYXdd98dVeX+++9nzpw5zJkzh2+//ZauXbsW+vwbb7yRLl26sGDBAl5//XU2bdpUtgHGqTCPaV+tqrOD7zcAi7EjfS8CblfVzcF9a8KKIW6JWMaYMQNefRWqVLGKhowMeP310IYyIpYHr7sOPvkEvv/e9hC3bGldmY46yqr/zj0XXn55x39VOpdMYnFixtq1a+nfvz+XXHIJIsIxxxzDww8/zJbgX4FLliwh7zaZGjVqsGHDhr9+Xr9+PXsHm/vHjRtXdsHFuZhU64lIOtAamAk0BQ4VkZki8r6IHFzAczJzzr3PLuNpr7ghAieeCHPm2Mhp0yb7uWNH240bsj32sO5Lr71m038TJ0L37vZzr17WWf2kk+Dpp+GXX0IPx7mYGj4cqlX7+7Vq1ez6zsg5Cbd58+YcddRRdO3alSFDhgBw/vnn06xZM9q0aUOLFi248MILyfv3rUuXLixatIiMjAyef/55Bg4cyKBBg2jduvUOj01qYS9qAdWBz4Fewc8LgPsBAdoD3xJsBi7olhQFEcXx55+qo0dbmR2odu2qOmtWJGFMmaI6YMD2UCpWVO3WTfXxx1XL6PgZ58pcPFTrRSXZCiJC7RAhIhWBN4B3VPWu4NrbwAhVnRb8/A3QQVXXFvQ6CdEhoiz98Qc89JAtDq1bZ0fu/uc/1mQvxrZts/5/L71kt2+/tS4Vhx9uYfXqZVOBzsUD7xDhHSKKJCICjAEW5ySmwCtAl+AxTYFKwE9hxZGQqlaFq6+GZcvgppvg7bftbPjMTDu6I4YqVIAOHWDUKPjmG/j8c/j3vy2Miy+GvfayAxVHj4a1Bf7zwjnnSibMNaeOwDnAESIyJ7gdB4wFGovIAuA54FwNc/iWyGrWhJtvtqxwySUwbpx1oPj3v0vdEmlniECbNjYnv3gxzJtnp4VkZVmp+p57Whf1xx8vOLzx4+2cxgoV7Ov48bH8DVx5UB7/nCTj7+yNXxPJ8uXWVO/pp2HXXS1JXXbZjqu6MaZqieq//4Xnn7eS9dRU66R+5pnQowdUr26JKDPz71WA1arZqCtBT8Z2cebbb7+lRo0a1KlTBykn58yoKuvWrWPDhg00atTob/cl8rSeJ6dENH++7RR8/XUbrtx8s+26TU2NOjJUra/fhAnw3HM2qqpa1doMTpuW/9RfWprlXed21pYtW8jKyio3e4FyVKlShQYNGlCxYsW/XffkFDJPTgX46CMbPX3yiRVL3HYb9OwZNycTbttmoU2YYKOqgtakRELZf+xcuefJKWSenAqhaht5Bw2y/kT/+AeMHAmdOkUd2d9kZ9shwWvy2XK95562Idg5V7YSOTn5kRmJTsRGS/PnW5uHFSvg0ENtM++iRVFH95fUVDviKr/lsdWrrdDi3nvzT17OufLHk1OySE21Pn1Ll9oZUu+/DwcdZBUIcTIs6d3bih/S0iynpqXBgw9aUqpQAa64wkZXJ55oe6o2b446YudcVHxaL1n99JNt3H3oIahY0fZNXXst1KgRdWQFWrjQzml85hkbTdWubdV+551nIyvnXMkk8rSeJ6dk9803Vtn3/PNQr561Ir/gAktYcSo72zqkP/mk9fvbvBlatbIk1bs31KkTdYTOJQZPTiHz5FQGPv3URk4ffGDnSI0cafNncVLZV5D//c9K0seOhVmz7KiPHj2sYe3RR9t0oHMuf4mcnPyjXV60b2+dzl97zf6i9+wJnTvDZ59FHFjhdtsNLrrIwpwzx76fOhW6dYPGjW3mMsYdnZxLGiIyVkTWBB17cl/f4UBYEUkXkT9ydfx5JNTYfORUDmVnW4+hIUOsPO7MM62IIj096siKZfNmq54fPdqm/ypUsKM+MjMtaaWkRB2hc/GhqJGTiBwG/AY8paotgmtdgMHA8aq6WUR2V9U1wdFHb+Q8Lmw+ciqPUlOhf3/rM3TDDfDKK7aJd+DAhDi4qXJlOO00mDLFfoWBA2HmTEtQ6ekwbJgVVDjnCqeqHwB5O2HGxYGwnpzKsxo17C/5kiV2Eu8dd1hj2QcfTJjz2vfd1xpjrFpl5ecHHmiN3PfZB0491VomJcDkgHNhSc05tDW4ZRbjOYUdCNtIRL4Irh8aUsyAT+u53L74wkrOp02zkdQdd8Dxx8d90UReS5fCo4/CE09Yd/QDDrCBYp8+1i/XufKiOAUReafrgvWnacBlwMHA80Bj7Hij6qq6TkTaYscfNVfVX8OI3UdObrvWrW0R57XX7OcTToCjjrJKhATSpInl1awsK0ffddftG3wvvjiuGmc4F4+ygInBYbqfAtuAuqq6WVXXAajq58A32CgrFJ6c3N+JWFKaPx8eeADmzrUdsOefDz/8EHV0JVK1KvzrXzBjhh2SeOqpVpLevLnl3Fdfha1bo47SubjzCvkcCCsi9UQkJbjeGGgCLAsrCE9OLn8VK8KAAVZxcNVV8NRTNiS59VY7Rj7BtGlj03yrVtlhiV99ZdX0++0Hd94J69dHHaFzsSciE4DpwP4ikiUi/Sj4QNjDgHkiMgd4EeivqqGdehrampOINASeAuoDCoxW1Xtz3X81cAdQT1ULPabd15ziQE5Z3MsvW7XBiBFw+ukJtx6VIzvbihTvuw8+/NAOQ+zXDy6/HPKc1+ZcwvJNuPnLBq5W1WZAB2CAiDSDvxJXV2BliO/vytJ++1kvoWnTtje969gx7jfxFiQ1FU45xRpmzJplXScefNB+zVNOgY8/9io/56IUWnJS1dWqOjv4fgOwGNg7uPtuYCA2onKJpHNn+2s+ZgwsW2adJ849N246n5dG27bWbHb5chscTp1qx2F16GDl6b4u5VzsxWTNKShVbA3MFJEewHeqOreI52Tm1OZnZ2fHIkxXXCkp1oV16VK47jprfte0qS3mJOB6VI69996+Z+rBB2HdOhtFHXAAPPJIQv9qziWc0Pc5iUh14H1gOPA2Vj/fVVXXi8hyoJ2vOSW4ZcusqezEiXZI06hR9lc9QdejcmzdaktsI0fa7GW9enDZZVaOXrt21NE5VzRfcyqAiFQEXgLGq+pEYF+gETA3SEwNgNkiskeYcbiQNW5s819Tp9qmotNOgy5drAw9gaWkWI6dOdOW2g4+GG680epBrr46oWcynYt7YVbrCfAk8LOqXlHAY5bjI6fksnWrNZUdPNjOu7jgAmsdXrdu1JGVifnzbSQ1YYIlr759bZ2qceOoI3NuRz5yyl9H4BzgiFwt1o8L8f1cPEhJgQsvtPWoSy6xRNWkidVsJ0i/vsIcdBA8/bS1I+zb1/ZONW0K55xjJ/k658qG99Zz4Vq40HoHTZkCzZrBvfdae4Yk8f33ton3kUdg40abBrzpJktizkXNR07OFaR5c5g0yXa8btpkx9fusosVS6Snw/jxUUe4U/bay5LTihU2k/nOO9CypSWpefOijs65xOUjJxc7TzxhU365p/eqVoXHHoPevaOLqwz9/DPcfbcNEDdsgF69bCTVqlXUkbnyKJFHTp6cXOykp9sQI6+6de1E3gQvPc/t55/hnnssSf36q42kbrnFzptyLlYSOTn5tJ6LnZUFdKv66Sfo2hUWL45tPCGqXduS0fLldtjw229DixbWTGNZaH2cnUsenpxc7OyzT/7Xa9e2lkgtW9pm3g0bYhtXiHbbzQ4bXrYMrrwSXnjBznHs3x+++y7q6JyLX56cXOwMHw7Vqv39WrVqVma+ZIkNK+64w/56P/tsUnVerVfPfrVvvoHMTDtXat994ZprrE2Sc+7vPDm52OndG0aPthZHIvZ19Gi7Xq+e7YmaOdOa3PXuDYcfnnQlb3vtZX37liyxxu533WVJasQI793nXG5eEOHiz9atNrQYNAh++QUuvRSGDrXWSElmwQL7Nd94w3LyLbfYADIlJerIXDLwggjnylJKirU9WrLEvt57r7UGHz8+qab6wIokXn8d3nsPGjSwAw9btrRklWS/qnMl4snJxa/ateHhh22qr0EDOPtsayibhH2CDj8cpk+HF1+0bWAnnGD7lRO8d65zpebJycW/gw+GGTPg0Uet82pGhlUSJFFVH9gy3MknW+697z744gto3doGjz/8EHV0zsWWJyeXGFJSrMztq6+gTx/rGXTggXZUR5LNf1WsaMtsX39tbQnHjbPeubfe6kUTrvzw5OQSS9261u5o+nSr8DvlFDjuOPtLnmR2282q+RYtgiOPtN59SZqPnduBJyeXmDp0sONp770XPvnEKgtuvtmayyaZJk2sb+7UqVCzpuXjo4+2pOVcsvLk5BJXaqqdm/7ll7ZYM3SoJalJk6KOLBRdusDs2XD//fD559ZM9qqrYP36qCNzruyFlpxEpKGITBORRSKyUEQuD66PEpEvRWSeiLwsIrXCisGVE3vuaWXm775ra1PHHANnnAGrV0cdWZlLTbUzHJcssaW3e+6xhhrjxsG2bREH51wZCnPklA1crarNgA7AABFpBkwGWqhqS2AJMCjEGFx5csQR1lHilltsHuyAA+CBB2xTb5KpV8+W3j791Jq99+0LnTsnZZW9K6dCS06qulpVZwffbwAWA3ur6iRVzQ4eNgNoEFYMrhyqXBluvNFaL3ToYGVvhxxi82BJqF07W3J7/HFLTBkZ1nFi48aoI3Nu58RkzUlE0oHWwMw8d50HvFXAczJFZJaIzMrOzs7vIc4VbL/97JyK556z9t/t28Pll9vhSkmmQgXrLPHVV7ZP+fbb7QDi//u/qCNzrvRCT04iUh14CbhCVX/NdX0wNvWX7zndqjpaVduparvU1NSww3TJSAROP90KJi66yCoJmjWDl19OylrsunXtsOH33rMDhrt3tzoRP5rDJaJQk5OIVMQS03hVnZjreh+gO9BbE6HzrEtsu+5qa08zZthf8F69oGfPgg8/THCHHw5z5tim3TfftHz82GNJmY9dEguzWk+AMcBiVb0r1/VuwEDgRFX1mXEXO+3b26GGd9wBU6bYX+2774YknDauVMnWnubPhzZtrLnGkUcm5V5ltxNEZKyIrBGRBXmuXxpUVS8UkZG5rg8Ska9F5CsROSbU2MIauIhIJ+BDYD6QU+R6PXAfUBnIOWJthqr2L+y1/MgMV+ZWrIABA2xhpnVrG1q0bRt1VKFQtYKJa66xprLDhtnym8+WJ7+ijswQkcOA34CnVLVFcK0LMBg4XlU3i8juqromqLaeALQH9gKmAE1VNZRyWD/PyZVfqtYL6NJLYc0a+4t9yy1QvXrUkYXiu+/g4ovhtdesym/sWDjooKijcmEqznlOQcHaG7mS0wvAaFWdkudxgwBU9bbg53eAoao6PYzYvUOEK79ErBfQ4sU273X33Vbm9uabUUcWir33tu1fzz9vA8e2ba2y76mnbK9UhQr2dXy+JUouQaXmVD0Ht8xiPKcpcKiIzBSR90Xk4OD63sCqXI/LCq6FwpOTc7Vq2blRH35oo6bjj7cOE0l4ToUInHaa9eXr0cPWpfr2tWSlal8zMz1BJZHsnKrn4Da6GM9JBWpjzROuBV4IaghiypOTczk6dbJDlIYNs3LzAw+0ua8EmPouqbp14YUX7GvetkcbN1oHdFduZQET1XyK1QzUBb4DGuZ6XIPgWig8OTmXW6VKcMMN1gapZUvb3ZqkZW4isG5d/vclaZW9K55XgC4AItIUqAT8BLwGnCEilUWkEdAE+DSsIDw5OZef/feHadNg9GhrBX7QQTBihJW7JZF99sn/eu3aSTlgdHmIyARgOrC/iGSJSD9gLNA4KC9/Djg3GEUtBF4AFgFvAwPCqtQDr9Zzrmjff28VfRMnWvO6xx9PmrLz8eNtjSl3L74KFWyq78wzbSlu112ji8/tnOJU68UrHzk5V5S99rKS84kT4ccfbTPvwIFJ0V21d28bHKal2TRfWpodvzFsmK1JtWoFH30UdZSuPPKRk3Ml8csv8O9/21/0/fazzbudO0cdVShmzoSzzoLly61A4qabfONuovGRk3PlRa1a8Oijdma6qh1Pe+GFSXkc7SGHWI++c86xkdRhh1mpuXOx4MnJudLo0sUq+q65xtagmjWz1gtJpkYNm+Z77jk7L6p1a3j99aijcuWBJyfnSqtaNRg1yua/6tSxXa1nnQVr10YdWZk7/XQrWmzUCE48cXufPufC4snJuZ3Vrp11O7/5ZnjxRRtFPf980tVi77uvnbo7YADceadN8/l+KBcWT07OlYVKlaxiIGd4ccYZdm7U6tVRR1amKle2o7FeeGH7sfBvvBF1VC5eiUgTEXlRRBaJyLKcW3Ge68nJubLUooUNL0aNsmPimzWDJ59MulHUqadaHk5PhxNOgOuvh62hbcd0CewJ4GHs1PMuwFPAM8V5oicn58paaqotysyda50l+vSB446DVauKfGoi2W8/y8MXXAC33Wa/YkHtkFy5VVVV38W2La1Q1aHA8cV5oicn58LStCm89x7cfz988IGNqsaMSapRVJUqtuVr9Gj7Vdu1s/Jz5wKbRaQCsFRELhGRk4BiHZgW5jHtDUVkWjDXuFBELg+u1xaRySKyNPi6W1gxOBe5ChXgkku2n5d+/vnQrVvSVRJccIHl3y1b4J//hGeKNXHjyoHLgWrAZUBb4GzgX8V5Ypgjp2zgalVthp0LMiA45vc64F1VbQK8G/zsXHJr3BjefRceegg+/tgONXz00aQaRR1yCHz+uXV3OuccO1jYy83LvXRV/U1Vs1S1r6qeDBTQbvjvCm1fJCI1VfVXEamd3/2q+nNxIxSRV4EHgltnVV0tInsC76nq/oU919sXuaSyfLkdxTF1qh3HMWaMNbVLElu2WOvBe+6xzk4vvmjbwFzsRd2+SERmq2qboq7l+9wiktMbqtpdRL4FFMh9GqKqauNiBpgOfAC0AFaqaq3gugD/y/k5z3MygUyASpUqtd28eXNx3sq5xKBqCzXXXGMdV++806b8Yn/gaGieecZycMOGVm5+wAFRR1T+RJWcRORY4DjgNOD5XHfVBJqpavsiXyPsxq8iUh14HxiuqhNF5JfcyUhE/qeqha47+cjJJa3co6iuXa0VUsOGRT4tUUyfDj17wubNtjeqa9eoIypfIkxOrYDWwM3ATbnu2gBMU9X/FfkaxUlOItIRmKOqv4vI2UAb4B5VLXRVV0QqAm8A76jqXcG1r/BpPee227YNHnnE5sJSUmw+rE+fpBlFrVhhLY8WLLBf7ZJLkuZXi3txMK2XqqrZpXlucQsiHgY2BtnwauAb4OkighJgDLA4JzEFXgPODb4/F3i1RBE7l2wqVICLL7ZGsq1bw3nnQffudshhEkhLsxqQ7t3hssvsV/VCieQmIvNFZB4wW0Tm5b0V6zWKOXKaraptROQm4DtVHVPUopaIdAI+BOYD24LL1wMzsaN+9wFWAKcVVVjhIydXbmzbBg8+aGdGValivYLOPDMphhpbt1oniZEj4aijrFDCT9kNV4TTeoVW+KhqkYevFDc5vY+dGd8XOAxYA8xV1YOKF+rO8eTkyp2lS+Hcc23RplcvOy99992jjqpMPPGEHQ3frBm8+SbsvXfUESWvqKf1ghjSgCaqOkVEqgKpqrqhqOcVd1rvdGAz0E9VfwAaAKNKHa1zrnBNmsCHH9ow4403rLvExIlRR1Um+vaF//s/WLYMOnSwtSiXnETkAuBF4NHgUgPgleI8t9DkJCJVROQK4Abgd2A6gKquVNWnShmvc644UlLg2mutw2rDhnDyyXD22fC/Igud4l7XrpZ7t26FTp2s9ZFLSgOAjsCvAKq6FCjWFEBRI6cngXbYutGxwJ2lj9E5VyrNm8OMGTB0qJ0TddBBMGlS1FHttIwMm7Xcay845hiYMCHqiFwINqvqnzk/iEgqtme2SEUlp2aqeraqPgqcAhxa+hidc6VWsSIMGWJJatdd7a/5xRdDgq/F5lTydehghwiP8sWCZPO+iFwPVBWRo4H/Aq8X54lFJae/Cj5LW6vunCtDbdtaA7urrrK9URkZdm5FAtttN3jnHTjtNNvqdd11SdVysLy7DliLzb5dCLyJLRMVqaj2RVuxtSaw1kVVgY3B96qqNUsfc/F5tZ5z+Xj/favoW7XKSs+HDrUTeRPU1q22QfeRR+DCC62iPiUl6qgSW5xU69UDUNW1JXleoSMnVU1R1ZrBrYaqpub6PiaJyTlXgMMPt427ffvaaX/t2yd06VtKijVtHzTIGrb37g1//ln081z8ETNURH4CvgK+EpG1wV7ZYvHDBp1LZDVrWj++116D1att2u+uu2wzbwISgVtvhREjrPajZ0/YuDHqqFwpXIlV6R2sqrVVtTZwCNBRRK4szguE3vi1LPi0nnPFsGaN7W599VXo0gXGjYN9inV0TlwaPRr697dS89df924SpRFhh4gvgKNV9ac81+sBk1S1dVGv4SMn55LF7rvDyy/b+VCffWYl508/nbDVBZmZVl4+fbrl2p9+Kvo5rmREZKyIrBGRBbmuDRWR70RkTnA7LrieLiJ/5Lr+SCEvXTFvYoK/1p0qFic2T07OJRMRaxw7bx60bAn/+peVwa1bF3VkpXL66TZjuXgxHHEErC3RkrorhnFAt3yu362qGcHtzVzXv8l1vX8hr1vYamGxVhI9OTmXjBo1srYLt99u03wtW8LkyVFHVSrHHmvTekuXWoJasybqiJKHqn4AFPtE8xJoJSK/5nPbABSrJ6snJ+eSVUqKlZjPnGkLNl27whVXwB9/RB1ZiR11lLUY/OYbT1AllCois3LdMov5vEuC4y3Gikjuw2AbicgXIvK+iBTYlCFPpXfNPFXfxZrW84II58qDP/6wRHX//dYOafx4aNUq6qhKbNo0OP54GxhOnQr160cdUXwrTkGEiKQDb6hqi+Dn+sBPWJuhYcCeqnqeiFQGqqvqOhFpizVwba6qv4YRu4+cnCsPqlaF++6Dt9+29af27eGOOxKu5LxLFztmY/ly6NzZqudd2VLVH1V1q6puAx4D2gfXN6vquuD7z7FDZ5uGFUdoyamAKpAMEZkRVHrMEpH2Yb2/cy4fxxwD8+fb8OPaa+HooyErK+qoSqRzZ3jrLWuM0aWLJ6iyJiJ75vrxJGBBcL2eiKQE3zcGmgDLwoojzJHTOHasAhkJ3KyqGcBNwc/OuViqWxdeesk2786YYcUSL70UdVQlcthhNgjMyrL86mXmpSMiE7CjkPYXkSwR6QeMzHXMehdsQy3YQbPzRGQOdkZT/6JOMd+p2MJcc8pnLvMdYKyqPi8iZwInqOpZRb2Orzk5F5KlS61P0GefWQn6vfdC9epRR1Vs06bBccfZMtq77/pG3bziobdeacV6zekKYJSIrALuAAYV9EARycypMMnO9obozoWiSRMYMMDaII0day3Cb7456qiKrUsXePFFmDsXundP+BNEXC6xTk4XAVeqakNsqDimoAeq6mhVbaeq7VJTU2MWoHPlyvjxdi7Ur0HBVXa2dTc/7TRrE54Ajj8enn3WTg456STYvDnqiFxZiPW03nqglqqqiAiwvjjdzX1az7mQpKfDihX533fYYfDMM3ZEfAIYN84atPfsCS+8YOczlnc+rVd83wOHB98fASyN8fs753JbubLg+2bPtmKJF1+MXTw7oU8f28b1yiuWpBKsSt7lEWYpeX5VIBcAd4rIXOBWoLi7lZ1zYSioa3laGnzxha1JnXoqnH9+QizoXHKJHbkxfrx9nwA9BlwBvEOEc+XZ+PHW/jv3oUnVqtl5Fb17w5YttgZ1222WqCZMgDZtIgu3uP79bxg5EoYPh+uvjzqa6Pi0nnMuMfXubYkoLc06mqelbU9MYAs3w4dbr6Dff4cOHeDuu+N+zuy22+Dss2HwYFuLconHR07OueJZt86m9155Bbp1s7/6cdzc7s8/rbx86lTran7ssVFHFHs+cnLOJb86dWDiRHjoITuOo1UrmDQp6qgKVKmSNb5o2dKWzT77LOqIXEl4cnLOFZ8IXHSR/aWvW9d69V1zjQ1T4lCNGtYotl492w/1zTdRR+SKy5OTc67kWrSwBHXxxXDnndCxI3z9ddRR5WuPPeCdd2yZ7Jhj/CyoROHJyTlXOlWrwoMPwssv25CkTRtr1RCHmja1wwq//x5OOCEhz1ssdzw5Oed2Ts+eMGeOLe707m0NZOOwgKlDB6uEz+lxmwC1YOWaJyfn3M7bZx8rkrjxRqvia9vWElac6dEDbr8dnnsOhg2LOhpXGE9OzrmykZoKt9xiZ1ds2ACHHAIPPBB3Q5Rrr4Vzz4UhQ6wHn4tPvs/JOVf21q61ZndvvmmtwseMseM44sTmzXDkkfD55/DBB3DwwVFHFA7f5+Scc7nVq2c7X++80762bg3Tp0cd1V8qV7Y6jj32sKm+BDupvlzw5OScC0eFCnDVVfDxx/b9oYdaw7s4aX2Ukz83bLAE5ZMz8cWTk3MuXO3b2/EbJ51kHVmPOy5uNhu1aGHFEV98YetQcZI3HZ6cnHOxUKuWVR88/LBV9WVk2Nc4cPzxMGqUtToaMSLqaFwOT07OudgQgf794dNPra/QkUdadV8cHAd/1VVwxhlwww0weXLU0Tjwaj3nXBQ2bLAefePHW5J65hmrTojQb7/ZRt0ffrAqvrS0SMMpE16tlw8RGSsia0RkQZ7rl4rIlyKyUERGhvX+zrk4VqMGPP00PP64FUxkZNj+qAhVr25N17dsgVNOgU2bIg2n3AtzWm8c0C33BRHpAvQAWqlqc+COEN/fORfPRKBfP+sntNtucPTRtjM2wmm+pk3hqadg1iy49NLIwnCEmJxU9QPg5zyXLwJuV9XNwWPio2THORedFi0sG5xzjq1BHX20za1FpEcPGDTIBnWPPx5ZGOVerAsimgKHishMEXlfRArcly0imSIyS0RmZWdnxzBE51zM7bILPPkkjB0LM2bYNN/UqZGFM2wYHHUUXHKJ5U0Xe7FOTqlAbaADcC3wgohIfg9U1dGq2k5V26WmpsYyRudcVPr2tWq+nGm+iKr5UlKsg3n9+nDyyXZCvYutWCenLGCimk+BbUDdGMfgnItnOQcZnnmmrUF16xbJpt26dW3v0+rVljMToLA5qcQ6Ob0CdAEQkaZAJeCnGMfgnIt31atbNd9jj8FHH9k03wcfxDyMdu1sY+7rr1uDdRc7YZaSTwCmA/uLSJaI9APGAo2D8vLngHM1ETZaOediTwTOPx9mzrTS8yOOsEwR4x5DV1xhXSSuuSYuj6jaKflt+RGRoSLynYjMCW7H5bpvkIh8LSJficgxocaWCLnBN+E6V879+itccIG1QOre3YonateO2duvXQutWkHNmrZBd5cE2dZa1CZcETkM+A14SlVbBNeGAr+p6h15HtsMmAC0B/YCpgBNVTWURUFvX+Sci381a1qH1vvvh3fegTZtrHAiRurVs2YWS5Yk1/6nArb8FKQH8JyqblbVb4GvsUQVCk9OzrnEIGK13R9/bD936mTJKkazP126wODB8MQTVsmXIFJztuQEt8xiPu8SEZkXTPvlnBK5N7Aq12Oygmuh8OTknEssBx9sR3Accwxcdpl1bN2wISZvPWQIdOwIF14I33wTk7fcWdk5W3KC2+hiPOdhYF8gA1gN3BlmgAXx5OScSzy1a8Orr8Ltt8OLL1rCWrCg6OftpNRUm95LSbFK9z//DP0tY05Vf1TVraq6DXiM7VN33wENcz20QXAtFJ6cnHOJqUIFO7xw6lT45Rc71PDpp0N/27Q0GDPGtmINHRr628WciOyZ68eTgJys/xpwhohUFpFGQBMgtIU/T07OucR2+OF2lG379vCvf0FmZugtxXv1sp61I0ZYt6VEVcCWn5EiMl9E5mH7Uq8EUNWFwAvAIuBtYEBYlXrgpeTOuWSRnQ033mhTfa1b23Rf48ahvd2vv0LLllC5suXGatVCe6tS8/OcnHMuaqmpcNtt1s7h22+hbVv7PiQ1a1rl3pIl1sXclS1PTs655NK9u1XzNW4MJ54I119vo6oQdOliBYP33RdpE/Wk5NN6zrnktGmTZY7HHrPWRxMmwO67l/nbbNxos4ibNsH8+Taiihc+reecc/GmShUYPdrm3j75xDJIzgbeMlStmnVTysqCK68s85cvtzw5OeeSW58+VlJXtSp07gz33FPmXSU6dLCq9rFj4Y03yvSlyy2f1nPOlQ/r11uieuUVOP10O4O9evUye/nNm62afc0a2w9cp06ZvXSp+bSec87Fu113hYkTbXPSf/9rmeTLL8vs5StXhqeeslNzk6k5bFQ8OTnnyg8RGDgQJk+Gn36ytkf//W+ZvXyrVtYcdsIEeOutMnvZcinMwwZ3OMQq131Xi4iKiB/R7pyLvSOOsHLzFi3gtNPgqqtgy5YyeenrroMDD4SLLgJfjSi9MEdO44BueS+KSEOgK7AyxPd2zrnCNWgA779vc3B33w1HHgmrV+/0y1aubEWCK1bATTeVQZzlVGjJqZBDrO4GBgLxX4nhnEtulSrZDtrx4+2I2zZt4KOPdvplO3WyYzXuucde1pVcTNecRKQH8J2qzo3l+zrnXKHOOsvKzatXt7YP99670+Xmt99ue34vuCC0BhVJLWbJSUSqAdcDxRroikhmzumN2f6/rHMubAcdBLNmwfHHwxVX2IFNv/1W6perVcsO6v3iC8t1rmRC3eckIunAG6raQkQOAt4FNgZ3NwC+B9qr6g+FvY7vc3LOxcy2bVZufsMNcMABVn6+//6leilV6NkTpkyxvU+NGpVtqEXxfU7FoKrzVXV3VU1X1XTs/Pk2RSUm55yLqQoVrM34O+/Ajz9aufnLL5fqpUTggQfsJS+6qMwbUyS1MEvJ8zvEyjnnEsNRR1m5+f772+mC118PW0t+tl7DhnDrrZbrJkwIIc4k5e2LnHOuMLm7mx99NDz7LNQt2RbNrVuhY0dYtszOf6pVK5xQ8/JpPeecS1Y53c0ffxw++MAOMZw1q0QvkZICDz1kTSluvjmkOJOMJyfnnCuOfv2274Hq1AnGjCnR09u0sbLy+++HRYtCiC/JeHJyzrniatfOdtUeeiicfz5kZlo78mIaPhxq1LBZwgRYUYmUJyfnnCuJunXh7betou+xx+Cww2DVqmI/ddgwePfdUhcAlhteEOGcc6X18stw7rm2LvX889ZdogjZ2TbF9+uvNr1XrVp44XlBhHPOlUcnnQSffmonCx51FNxxR5Hzdamptu60YgWMGhWjOBOQj5ycc25nbdgAffvCSy/Bqafaee1FnLJ7xhnw6quweDGkp4cTlo+cnHOuPKtRww4tHDHCEtQhh9iGpkKMGmWdI66+OkYxJhhPTs45VxZyTtmdNGl726NXXy3w4Q0bWtOJiROt9577O5/Wc865srZiBZx8spWd33ADDB1qO3Hz2LQJmje3AwrnzbP1qLLk03rOOee2S0uzDbvnnQf/+Q907w4/73j2apUqcNddtu5Uwj29Sc+Tk3POhaFKFWt59OijtrGpXTuYM2eHh514ojWcGDJkp46PKhURGSsia0RkQT73XS0iKiJ1g587i8h6EZkT3EI9hN6Tk3POhUXEukh88AH8+Sf885/wzDM7PGTkSFumuuuumEc4DuiW96KINAS6Aivz3PWhqmYEt1vCDMyTk3POha1DB1t/at8ezjnH+hdt2fLX3f/4hy1R5SSpWFHVD4Ad5xvhbmAgEFlRgicn55yLhfr1YfJkuOoq24V7xBHww/azVm+91Qokbgl1PFI0EekBfKeqc/O5+x8iMldE3hKR5mHG4cnJOedipWJFuPNOO3Vw9mzrY/TJJwA0bWozgKNHF7lFqiRSRWRWrltmYQ8WkWrA9UB+60mzgTRVbQXcD7xSZlHmI8yTcHdYaBORUSLypYjME5GXRaRWWO/vnHNx64wzYMYMa6zXubMd9qTKkCFWVn799WX2Ttmq2i7XbXQRj98XaATMFZHlQANgtojsoaq/qupvAKr6JlAxp1giDGGOnMax40LbZKCFqrYElgCDQnx/55yLXwcdZIcWdu0KAwZA377Ur/kH115rTSZmzIh9SKo6X1V3V9V0VU0HsoA2qvqDiOwhIgIgIu2x/LEurFhCS075LbSp6iRVzQ5+nIFlZeecK59q1YLXXrNNuk8+CR07cvXJy6lfH669Nvwzn0RkAjAd2F9EskSkXyEPPwVYICJzgfuAMzTELg6hdogQkXTgDVVtkc99rwPPq+ozOzzR7s8EMgEqVarUdnMJDvRyzrmE88YbcPbZkJLCI70/5KL7m/Hqq7YPqrQSuUNEJMlJRAYD7YBexcm83r7IOVcufP019OrFlvlfsk+VH/lpU3WySSEt5XuGZy6n90OdSvRyiZycYl6tJyJ9gO5A7zCHhM45l3D22w+mT+eFJtfz86aqZFMRqMCKrQ3IfLg14y/+KOoIYyamIycR6QbcBRyuqmuL+zo+cnLOlSfpqVms2LrjknxaShbLs4u/VJ/II6fQklOw0NYZqAv8CAzBqvMqs73CY4aq9i/qtTw5OefKkwqyDc1nYkvYxjYt/oSXJ6eQeXJyzpUnPnLyDhHOORd3hmcupxp//wd5NX5neObyaAKKgCcn55yLM70f6sToi74gLSULYRtpKVmMvuiLElfrJTKf1nPOuSTl03rOOedcGfLk5JxzLu54cnLOORd3PDk555yLO56cnHPOxZ2EqNYTkW3AH1HHEUgFsot8VLQ8xp0X7/FB/McY7/FB8sdYVbUELSXiSEIkp3giIrNUtV3UcRTGY9x58R4fxH+M8R4feIzxLCEzqnPOueTmyck551zc8eRUcqOjDqAYPMadF+/xQfzHGO/xgccYt3zNyTnnXNzxkZNzzrm448nJOedc3PHkVAQRWS4i80VkjojMynX9UhH5UkQWisjIeItRRDJEZEbONRFpH2F8tUTkxeC/12IR+YeI1BaRySKyNPi6W1TxFRLjqODneSLysojUiqf4ct13tYioiNSNKr7CYoyXz0oB/xvH0+dk/yCOnNuvInJFvH1WYkZV/VbIDVgO1M1zrQswBagc/Lx7HMY4CTg2+P444L0I43sSOD/4vhJQCxgJXBdcuw4YEfF/w/xi7AqkBtdGRBljfvEF3zcE3gFW5P3/QDzEGE+flQLii5vPSZ5YU4AfgLR4+6zE6uYjp9K5CLhdVTcDqOqaiOPJjwI1g+93Bb6PIggR2RU4DBgDoKp/quovQA/sjwXB155RxAcFx6iqk1Q1Z2f+DKD452PHIL7g7ruBgdj/3pEpJMa4+KwUEl9cfE7ycSTwjaquII4+K7HkyaloCkwSkc9FJDO41hQ4VERmisj7InJwhPFB/jFeAYwSkVXAHcCgiGJrBKwFnhCRL0TkcRHZBaivqquDx/wA1I8oPig4xtzOA96KfWhAAfGJSA/gO1WdG1FcuRX03zBePisFxXcF8fE5yesMYELwfTx9VmLGk1PROqlqG+BYYICIHIb1uqoNdACuBV4QEYmzGC8CrlTVhsCVBP9ijEAq0AZ4WFVbA79jUxN/UZuviPJf/oXGKCKDsd5m46MJL9/4hgLXAzdFFFNeBf03jJfPSkHxxcvn5C8iUgk4Efhv3vvi4LMSM56ciqCq3wVf1wAvA+2BLGCimk+BbUBki9EFxHguMDF4yH+Da1HIArJUdWbw84vYH4kfRWRPgOBrlFOjBcWIiPQBugO9gz8MUSgovkbAXBFZjk05zhaRPaIJscAY4+WzUlB88fI5ye1YYLaq/hj8HE+flZjx5FSIYOqkRs732AL5AuAVbKEXEWmKLa7+FGcxfg8cHjzsCGBpFPGp6g/AKhHZP7h0JLAIeA37w0Dw9dUIwgMKjlFEumHrOSeq6sY4i2+2qu6uqumqmo798W0TPDZeYlxEnHxWCokvLj4neZzJ9ik9iKPPSix5h4hCiEhjbCQCNi3wrKoOD4bdY4EM4E/gGlWdGmcxdgLuDa5tAi5W1c8jijEDeBz7w7QM6Iv9w+gFYB+s0uw0Vf05ivgKifEzoDKwLnjYDFXtHy/xqer/ct2/HGinqpH8IymIIYMd/xv+Tvx8VvKLrzlx8jkJYtwFWAk0VtX1wbU6xNFnJVY8OTnnnIs7Pq3nnHMu7nhycs45F3c8OTnnnIs7npycc87FHU9Ozjnn4o4nJ+cAEekZdPY+IM/1TiLyadDJ+stc7aEQkaEick0Br3eFiPwr+H4vEZkqIq+KSPXg2iUicl6Yv5NzicyTk3PmTOCj4CsAQbeFZ4H+qnoA0Am4UESOL+yFRCQV68X3bHDpMuBSbI/N2cG1scE151w+PDm5ci8YzXQC+mENN3MMAMap6myAYIPrQPL0BszHEVgHh5yO5ilY255tgASvtRFYHuX5Qc7FM09OztmRBG+r6hJgnYi0Da43B/J2C5gVXC9MxzzPewB4FOgPPJPntQ4tbdDOJTNPTs7ZVN5zwffPkWtqr5T2xI5nAEBVV6jqYap6gqpuyPW4NcBeO/leziWl1KgDcC5KIlIbm4Y7SEQUm4JTEbkWawzalr832mwLLCziZf8AqhTj7asEj3XO5eEjJ1fenQI8rappQYfvhsC32HTbg0CfoGFoTgPOEdix2YVZDOxXjPduinWQd87l4cnJlXdnsr2re46XgDOD00fPBh4TkS+BT4Cxqvp6Ea/5FnYkeFE6ApNLGK9z5YJ3JXcuBCLyMjBQVfM9H0hEWgNXqeo5sY3MucTgycm5EASH2tVX1Q8KuP9oYKmqLo9pYM4lCE9Ozjnn4o6vOTnnnIs7npycc87FHU9Ozjnn4o4nJ+ecc3HHk5Nzzrm48//GBRDE/gmftQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = objective.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "________________________________________________________________________________\n",
      "Parameters: 'instrument parameters'\n",
      "<Parameter:'delta offset' , value=0 (fixed)  , bounds=[-10.0, 10.0]>\n",
      "________________________________________________________________________________\n",
      "Parameters: 'Structure - ' \n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:  ' - thick'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:  ' - rough'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n",
      "________________________________________________________________________________\n",
      "Parameters:    'PNIPAM'    \n",
      "<Parameter:  ' - thick'   , value=132.339  +/- 982 , bounds=[100.0, 500.0]>\n",
      "<Parameter: ' - cauchy A' , value=1.47 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter: ' - cauchy B' , value=0.00495 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter: ' - cauchy C' , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:  ' - rough'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n",
      "________________________________________________________________________________\n",
      "Parameters:    'Silica'    \n",
      "<Parameter:  ' - thick'   , value=1 +/- 1.01e+03, bounds=[1.0, 30.0]>\n",
      "<Parameter:  ' - rough'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:  ' - thick'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:  ' - rough'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n"
     ]
    }
   ],
   "source": [
    "for i, x in enumerate(objective.model.parameters):\n",
    "    print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also view the resultant refractive index profile of the interface as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(110, 3.7, '658.0 nm')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAX8UlEQVR4nO3dfXBV1bnH8e+TF6KoF4sJ4R20CiKUF42iYGe8QpVa0etbxaGCta3Tqq0yeu+Vyui0joMdb+1VxHoZBLR1VLQIFrEVKxatiII3viIKll5FkGBFS21Rkuf+cXZiEhNykpxkn6z1+8ycyTl7r5zzsDW/rKy99l7m7oiISNdXkHYBIiKSGwp0EZFAKNBFRAKhQBcRCYQCXUQkEEVpfXBpaakPHjw4rY8XEemS1q9fv9Pdy5ral1qgDx48mHXr1qX18SIiXZKZ/aW5fRpyEREJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUCkNg+9rTZu/xuPvvxeVm3NjHOO7s/AQ7p3cFUiIunrcoG+acdu5qzalFVbd/i0uob/nHRkB1clIpIH3D2VxzHHHOMdbeisFX7jo693+OeISG689957Pm3aNC8tLfWSkhIfNmyYP/XUU3X7p0+f7kCDx9ixYxu8x7Zt2/xb3/qWl5eX+/777+8jR470X//61y1+9kMPPeTDhg3zbt26+bBhw3zJkiU5//flArDOm8nV4MfQXSsyiXQJu3btYvz48bg7jz76KBs2bGDOnDn06tWrQbuJEyeybdu2useKFSsa7J82bRobNmxg2bJlvPrqq0ybNo0LL7yQ1atXN/vZa9as4fzzz2fq1KlUVlYydepUzjvvPNauXdsh/9YO01zSd/SjM3roR856zG/47Wsd/jki0n4zZ870cePG7bPN9OnT/Rvf+MY+2xxwwAG+YMGCBtsGDhzoN998c7Pf881vftMnTpzYYNuECRN8ypQpzX7PwoUL/YADDvAnnnjChw8f7t27d/eTTjrJ33777bo2119/vQ8fPtwXLVrkgwYN8u7du/tFF13ke/bs8blz53r//v29Z8+ePmPGDK+urt7nv6sWsfbQzTJ/k4lI/lu6dCljx47l/PPPp1evXowePZrbb7/9C39lP/PMM/Tq1YshQ4bwve99jx07djTYf+KJJ7J48WI++OADampqWLZsGVVVVUycOLHZz16zZg2nnHJKg22nnnoqzz777D5r3rNnD7Nnz2bBggWsWbOGXbt28f3vf79Bmy1btrBs2TKWL1/OkiVLePDBBznjjDN44YUXePzxx5k/fz5z5szh4YcfzuYw7VOXOynaGkbmxKiI5L+3336bO+64gxkzZnDNNddQWVnJD3/4QwAuv/xyACZNmsTZZ5/NoYceypYtW5g1axYnn3wy69evp6SkBIDFixczZcoUSktLKSoqoqSkhPvuu4/Ro0c3+9nbt2+nvLy8wbby8nK2b9++z5r37t3L3LlzGTp0KABXX301F198Me6OmQFQXV3NwoUL6dGjByNGjGDSpEn88Y9/ZOvWrXTr1o1hw4Yxfvx4Vq1axTnnnNOmY1cr7EA3w9VHF+kSampqqKioYPbs2QCMGTOGt956i7lz59YF+pQpU+raf+UrX+GYY45h0KBBPProo5x99tkAzJo1i507d/LEE09QWlrK0qVLmTZtGqtXr2bUqFE5rbmkpKQuzAH69u3Lp59+yocffkjPnj0BGDhwID169KhrU15ezpAhQ+jWrVuDbY3/0miLsAMd9dBFuoo+ffpw1FFHNdg2bNgwbr311ma/p2/fvvTv35+33noLgM2bNzNnzhwqKyvrwnvUqFE8/fTTzJkzh/nz5zf5Pr179+b9999vsO3999+nd+/e+6y5qKhhhNb2ymtqauq2FRcXf6FNU9uqq6v3+VnZCHoMHUu7ABHJ1vjx49m4cWODbW+++SaDBg1q9nt27tzJ1q1b6dOnDwCffPIJAIWFhQ3aFRYWNgjZxk444QRWrlzZYNvKlSsZN25cq/4NaQs60DM9dHXRRbqCGTNm8Nxzz3HjjTeyadMmHnzwQW677TYuu+wyAHbv3s3VV1/NmjVr2LJlC0899RSTJ0+mV69enHXWWQAceeSRHH744Vx66aU8//zzbN68mZ///OesXLmyrg3AhAkTmDlzZt3rK664gieffJKbbrqJN954g9mzZ7Nq1SquvPLKTj0G7RV2oJtpBF2kizj22GNZunQpixcvZsSIEVx77bXccMMNXHrppUCml/3KK69w5plnMmTIEKZPn87QoUNZs2YNBx10EJAZ3lixYgVlZWVMnjyZkSNHcs8997Bw4UImT55c91mbN29m27Ztda/HjRvH/fffz6JFi+q+54EHHmDs2LGdexDaydLqwVZUVHhHryk6+qePM3lkX274txEd+jkiIp3FzNa7e0VT+8LuoaddgIhIJwo70DVtUUQiEnago2mLIhKPsANdl/6LSESCDnQw9dBFJBpBB3rmoi0luojEIexAR2PoIhKPsAPdFOgiEo+gAx3QtEURiUaLgW5m+5nZ82b2kpm9ZmY/aaLNRWZWZWaVyeO7HVNu65hOiopIRLK5fe4e4GR3321mxcAzZvaYuz/XqN0D7n557ktsO01bFJGYtBjoyRp2u5OXxcmjS+SkToqKSEyyGkM3s0IzqwR2ACvdvamlsM8xs5fN7CEzG9DM+1xiZuvMbF1VVVXbq86SLv0XkZhkFejuXu3uo4H+wHFm1vj2hb8FBrv7SGAlcHcz7zPP3SvcvaKsrKwdZbeC8lxEItGqWS7uvgtYBUxqtP0Dd9+TvJwPHJOT6tpJY+giEpNsZrmUmdnByfP9ga8BbzRq06feyzOADTmssc0y89AV6SISh2xmufQB7jazQjK/ABa7+3Iz+ymwzt0fAX5kZmcAe4G/Ahd1VMGtYbojuohEJJtZLi8DY5rYfl295zOBmY3bpE1DLiISk6CvFNW0RRGJSdiBrkWiRSQiYQc6OikqIvEIOtDRGLqIRCToQNf6FiISk7ADXZf+i0hEwg50NMtFROIRdqBrxSIRiUjQgQ5asUhE4hF0oGvFIhGJSdiBrmmLIhKRoAMdNIYuIvEIOtDNDPXRRSQWYQc66qGLSDzCDnSNoYtIRIIPdBGRWIQd6Jjutigi0Qg70DXkIiIRCTvQ0UlREYlH0IGOViwSkYgEHehasUhEYhJ2oGuWi4hEJOxAR2PoIhKPsANdKxaJSETCDnTUQxeReIQd6FqxSEQi0mKgm9l+Zva8mb1kZq+Z2U+aaFNiZg+Y2SYzW2tmgzuk2jbQkIuIxCKbHvoe4GR3HwWMBiaZ2fGN2nwH+NDdDwd+Afwsp1W2kVYsEpGYtBjonrE7eVmcPBrH5JnA3cnzh4AJZnkwaVCX/otIRLIaQzezQjOrBHYAK919baMm/YB3ANx9L/ARcEgT73OJma0zs3VVVVXtKjwbBkp0EYlGVoHu7tXuPhroDxxnZiPa8mHuPs/dK9y9oqysrC1v0SqZm3Mp0UUkDq2a5eLuu4BVwKRGu7YCAwDMrAjoAXyQg/raxUh/1EdEpLNkM8ulzMwOTp7vD3wNeKNRs0eA6cnzc4EnPQ9uoqJpiyISk6Is2vQB7jazQjK/ABa7+3Iz+ymwzt0fAe4CfmVmm4C/AlM6rOJW0P3QRSQmLQa6u78MjGli+3X1nv8TOC+3pbWfViwSkZiEf6Vo2kWIiHSSoAMdNIYuIvEIOtBNKxaJSETCDnRQF11EohF2oGsMXUQiEnagow66iMQj7EDXikUiEpGwAx310EUkHmEHui79F5GIBB3ooJOiIhKPwANdl/6LSDyCDvQ8WDNJRKTThB3oaRcgItKJwg50nRQVkYiEHehoHrqIxCPsQFcPXUQiEn6gp12EiEgnCTvQNW1RRCISdKCjHrqIRCToQM/cDz3tKkREOkfYga4Vi0QkImEHOmgMXUSiEXagawxdRCISdqCjeegiEo+wA10rFolIRFoMdDMbYGarzOx1M3vNzK5oos1JZvaRmVUmj+s6ptzWUQ9dRGJSlEWbvcBV7v6imR0ErDezle7+eqN2T7v76bkvsR106b+IRKTFHrq7b3P3F5PnfwM2AP06ujAREWmdVo2hm9lgYAywtondJ5jZS2b2mJkNb+b7LzGzdWa2rqqqqvXVtpLpjugiEpGsA93MDgR+A1zp7h832v0iMMjdRwFzgKVNvYe7z3P3CnevKCsra2PJ2cvcbVFjLiISh6wC3cyKyYT5ve6+pPF+d//Y3Xcnz1cAxWZWmtNK28DQPHQRiUc2s1wMuAvY4O63NNOmd9IOMzsued8PclloW+h+6CISk2xmuYwHLgReMbPKZNuPgYEA7n4ncC7wAzPbC/wDmOJ5MNahFYtEJCYtBrq7P0ML6y27++3A7bkqKlfUQxeRmAR+pajG0EUkHkEHOph66CISjaAD3bTChYhEJOxAR2PoIhKPsANdY+giEpGwAx3TlaIiEo2wA109dBGJSNiBjsbQRSQeYQe6achFROIRdKCDhlxEJB5BB7rpdugiEpGgAx1QF11EohF0oGfutigiEoewA10rFolIRMIOdDTiIiLxCDvQdT90EYlI4IGuFYtEJB5hBzrqoYtIPIIOdHQvFxGJSNCBbkp0EYlI2IFuaAxdRKIRdqCjMXQRiUfYga4RFxGJSNiBrhWLRCQiYQe6eugiEpGwAx2NoYtIPFoMdDMbYGarzOx1M3vNzK5ooo2Z2W1mtsnMXjazozum3FbSDdFFJCJFWbTZC1zl7i+a2UHAejNb6e6v12vzdeCI5DEW+GXyNVWKcxGJSYs9dHff5u4vJs//BmwA+jVqdiZwj2c8BxxsZn1yXm0r1XbQdWJURGLQqjF0MxsMjAHWNtrVD3in3ut3+WLoY2aXmNk6M1tXVVXVylLbTnkuIjHIOtDN7EDgN8CV7v5xWz7M3ee5e4W7V5SVlbXlLVrFkkEX5bmIxCCrQDezYjJhfq+7L2miyVZgQL3X/ZNtqdKQi4jEJJtZLgbcBWxw91uaafYIMC2Z7XI88JG7b8thnW1Se1JUcS4iMchmlst44ELgFTOrTLb9GBgI4O53AiuA04BNwCfAt3NeaRt83kNPtw4Rkc7QYqC7+zO0MAPQM2Mal+WqqFwxqx1DV6KLSPiCvlK0lnroIhKDoANdF4qKSEzCDvTaaYvqoYtIBMIO9NqTohpDF5EIhB3oyVf10EUkBmEHel0PXUQkfGEHet0YuiJdRMIXdqCrhy4iEQk60EVEYhJ0oNddKaouuohEIOxAr32iQBeRCIQd6JqHLiIRCTrQa2nIRURiEHSg637oIhKTsAPdNA9dROIReKBnvirORSQGYQd68lUddBGJQdCBjlYsEpGIBB3omocuIjEJO9A1hi4iEQk70LVikYhEJOxA15WiIhKRsAM9+aoeuojEIOxA1xi6iEQk7ED/fJ6LiEjwgg702jzXpf8iEoMWA93MFpjZDjN7tZn9J5nZR2ZWmTyuy32ZbaMxdBGJSVEWbRYBtwP37KPN0+5+ek4qyqHam3OJiMSgxR66u68G/toJteSceugiEpNseujZOMHMXgLeA65299eaamRmlwCXAAwcODBHH9282g76R//4jBWvbKNGyS4iHcAwTjy8lB7di1OtIxeB/iIwyN13m9lpwFLgiKYauvs8YB5ARUVFp6XrTb/bwJ82fdBZHyciEbrsX7/Mv596ZKo1tDvQ3f3jes9XmNkdZlbq7jvb+97tVdtD/9OmDzh1eDlXnzI03YJEJEhn3/Esf99TnXYZ7Q90M+sNvO/ubmbHkRmXz4vu8BG9DqJvj/2odue7Xz2MI8oPSrskEQlQt6ICPquuSbuMlgPdzO4DTgJKzexd4HqgGMDd7wTOBX5gZnuBfwBTPE8mfo/o14NnZ05IuwwRCVxRoVFdk37stRjo7n5BC/tvJzOtUUQkSkUFBXxWnX6gh32lqIhIJygqNPbWpD/kokAXEWmnogJjbx4MuSjQRUTaqbiwgL15cFJUgS4i0k6FBcZejaGLiHR9RYUFGnIREQlBcYFOioqIBKGwwDRtUUQkBMWFBXlxYZECXUSknYoKTbNcRERCUKQhFxGRMBQVaMhFRCQIRYXGZ5rlIiLS9RXpwiIRkTAUaZaLiEgYigstLxa4UKCLiLRToe62KCIShqIC3W1RRCQIxYXqoYuIBKGwoECzXEREQlCsJehERMJQVFBAjUNNysMuCnQRkXYqKjSA1K8WVaCLiLRTUUEm0NO+uEiBLiLSTkWFmShN+46LCnQRkXaq7aGnPRddgS4i0k61Y+h5P+RiZgvMbIeZvdrMfjOz28xsk5m9bGZH575MEZH8VVyQDLnke6ADi4BJ+9j/deCI5HEJ8Mv2lyUi0nUU5smQS1FLDdx9tZkN3keTM4F73N2B58zsYDPr4+7bclWkiEg+qx1ymbbgeboVttxPPv/YAXz3q4flvo4cvEc/4J16r99Ntn0h0M3sEjK9eAYOHJiDjxYRSd/xhx3C2WP68c+91Vm1Lz2wpEPqyEWgZ83d5wHzACoqKtK/8YGISA6U/8t+3HL+6LTLyMksl63AgHqv+yfbRESkE+Ui0B8BpiWzXY4HPtL4uYhI52txyMXM7gNOAkrN7F3geqAYwN3vBFYApwGbgE+Ab3dUsSIi0rxsZrlc0MJ+By7LWUUiItImulJURCQQCnQRkUAo0EVEAqFAFxEJhGXOaabwwWZVwF864aNKgZ2d8Dm50FVqVZ2511Vq7Sp1QteptbV1DnL3sqZ2pBboncXM1rl7Rdp1ZKOr1Ko6c6+r1NpV6oSuU2su69SQi4hIIBToIiKBiCHQ56VdQCt0lVpVZ+51lVq7Sp3QdWrNWZ3Bj6GLiMQihh66iEgUFOgiIoEIOtDNbJKZbUwWsL4m7XpqmdkAM1tlZq+b2WtmdkWyvaeZrTSzt5KvX0q7VgAzKzSz/zWz5cnrQ81sbXJcHzCzbmnXCJAsf/iQmb1hZhvM7IR8PKZmNiP57/6qmd1nZvvlyzFtalH45o5hmgvEN1Pnzcl/+5fN7GEzO7jevplJnRvN7NTOqrO5Wuvtu8rM3MxKk9ftOqbBBrqZFQJzySxifRRwgZkdlW5VdfYCV7n7UcDxwGVJbdcAf3D3I4A/JK/zwRXAhnqvfwb8wt0PBz4EvpNKVV90K/A7dz8SGEWm5rw6pmbWD/gRUOHuI4BCYAr5c0wX8cVF4Zs7hmkuEL+IL9a5Ehjh7iOBN4GZAMnP1hRgePI9dyT50FkW8cVaMbMBwCnA/9Xb3L5j6u5BPoATgN/Xez0TmJl2Xc3Uugz4GrAR6JNs6wNszIPa+pP5IT4ZWA4Ymavaipo6zinW2QP4M8mJ/nrb8+qY8vkavD3J3L56OXBqPh1TYDDwakvHEPgf4IKm2qVRZ6N9ZwH3Js8b/OwDvwdOSPOYJtseItPx2AKU5uKYBttDp/nFq/OKmQ0GxgBrgXL/fLWn7UB5WnXV89/AfwA1yetDgF3uvjd5nS/H9VCgCliYDA/NN7MDyLNj6u5bgf8i0yvbBnwErCc/j2mt5o5hPv+MXQw8ljzPuzrN7Exgq7u/1GhXu2oNOdDznpkdCPwGuNLdP66/zzO/nlOdU2pmpwM73H19mnVkqQg4Gvilu48B/k6j4ZU8OaZfAs4k8wuoL3AATfw5nq/y4Ri2xMyuJTOseW/atTTFzLoDPwauy/V7hxzoeb14tZkVkwnze919SbL5fTPrk+zvA+xIq77EeOAMM9sC3E9m2OVW4GAzq13tKl+O67vAu+6+Nnn9EJmAz7djOhH4s7tXuftnwBIyxzkfj2mt5o5h3v2MmdlFwOnA1OSXD+RfnV8m8wv9peRnqz/wopn1pp21hhzoLwBHJLMHupE5KfJIyjUBmTPZwF3ABne/pd6uR4DpyfPpZMbWU+PuM929v7sPJnP8nnT3qcAq4NykWep1Arj7duAdMxuabJoAvE6eHVMyQy3Hm1n35P+D2jrz7pjW09wxzKsF4s1sEpnhwTPc/ZN6ux4BpphZiZkdSuaE4/Np1Ajg7q+4ey93H5z8bL0LHJ38P9y+Y9qZJwY6+0Fm8eo3gc3AtWnXU6+uE8n82foyUJk8TiMzPv0H4C3gCaBn2rXWq/kkYHny/DAyPxCbgAeBkrTrS+oaDaxLjutS4Ev5eEyBnwBvAK8CvwJK8uWYAveRGdv/LAma7zR3DMmcIJ+b/Hy9QmbmTpp1biIz/lz7M3VnvfbXJnVuBL6e9jFttH8Ln58Ubdcx1aX/IiKBCHnIRUQkKgp0EZFAKNBFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRALx/8p514+OSqzzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "structure.reverse_structure = True\n",
    "plt.plot(*structure.ri_profile())\n",
    "plt.text(110, 3.7, f\"{structure.wavelength} nm\", fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using the plotting tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.insert(1, \"../tools\")\n",
    "from plottools import plot_ellipsdata, plot_structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAEYCAYAAABBfQDEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABPk0lEQVR4nO3de7yVc/r/8de127uTSnSSohLjUDopOktC5RCJEJVBjOSUcZicZ8yXYRg0DpEJE3IWyk9yNqNUKipGQ6UGJZKUtOv6/fG5t1a7vat9WOtea+338/G4H2ute933va57rb3XutZnXZ/Px9wdEREREREJcuIOQEREREQknShBFhERERFJoARZRERERCSBEmQRERERkQRKkEVEREREEuTGHcCOyMnJ8WrVqsUdhohUQGvXrnV3V2NCktStW9ebNm0adxgiUgHNnDnzW3evV9R9GZEgV6tWjZ9++inuMESkAjKzdXHHkM2aNm3KjBkz4g5DRCogM1tc3H1qFRERERERSaAEWUREREQkgRJkEREREZEEGVGDLCKps2HDBpYuXcrPP/8cdygpVbVqVRo3bkxeXl7coYiISMyUIIvIFpYuXUrNmjVp2rQpZhZ3OCnh7qxcuZKlS5fSrFmzuMMREZGYKUEWkS38/PPPFSo5BjAz6tSpw4oVK+IORZJk06ZNXHPNNaxevZr27dszZMgQ+vTpQ5MmTahRowa33XYbM2fO5E9/+hM1atSgV69eDBky5Nf933jjDcaNG0d+fj633noru+++e4xnIyLJpgRZRLZSkZLjAhXxnNPVxRdfzOzZs0u1b5s2bfjb3/621foXXniBpUuXUqdOHRo3bgxA9erV2bRpEw0aNADg/fffZ+TIkXTs2JEhQ4ZskSDfd999PP7448yfP5+xY8dyzTXX/HrfgQceyNChQ5k5cyZjx47llltu4fvvv+f777+nVatWv/468eijj5bqnEQk9dRJT0REst6nn35K586duf3227n33nsBeOqppxgzZgxfffUVc+fOpW/fvvz+97+nR48eDB06dIv93Z2cnByaNGnC0qVLt7ivcePGjBw5kk6dOv2a2A8cOJC7776bf/3rX/zf//0f1apV4/vvv0/FqYpIOcjKFuRrD32L+ctq0ajnvjTauzqNGkHjxtCoUVg0KZ9IeqtUqRIHHngg+fn57L///jz88MNUr169yG379u3LY489Ru3atVMbpCRNUS3AZdW4cWMqV64MhL8vgJyc0EZUv3591qxZw5gxY5gwYQJ77rknJ510EkccccSv++fk5LBp0yaWLFnyawt0gZ122gmAvLw81q9fD0CtWrWoUqUKtWrVAqBy5cq/3iciO27cuHFMmDBhu9t17dqVUaNGldvjZmWCvGaNM++/VXn1v/n8WMT9deuGhLlg2WOPLZfGjaFKlZSHLSKRatWq/doSN2jQIO677z4uvfTSIredNGlSCiOTTNW/f39GjBjBO++8Q/fu3QEYMmQI1atXJz8/n8svv5yff/6Zyy+/nJo1a9KhQwcAzjjjDB599FGGDRvG2WefzYYNG7jlllviPBWRCuUf//gHH374Ifvvv/82t/vxx6IyvtIzdy/XAybDTjvt5CWeanr+fLj8cla//DbLdu/AsiGjWLp3D5Z9lcPSpfDll/x6+d13W+/eoAE0aQJ77rl5adJk87LLLqCSRclGCxYs2O4b0RbGj4dRo2DJkvCPctNNMGhQmWKoUaMGa9asAULt59y5c7nmmmsYOHAgq1evJj8/n3vvvZdu3br9OlVx3bp1y/SYUPS5m9lad9+pzAeXIrVv39411bSIFKd79+7k5uby+uuvl/uxzWymu7cv6r6sbEEG4IAD4KWXqDV1KrUuu4z9/+9w6NAB7rgDunTZYtO1a0OiXLAsWbL58qOP4OWXYd26LQ9fo0ZIlJs23XJp1ixc7rqrEmipAMaPh2HDwj8RwOLF4TaUOUkGyM/PZ/LkyfTu3ZvHHnuMo446ilGjRrFx40bWFjymiIhkrbgacrM3QS5w+OEwcyb885/whz9A165w8slwyy0hkwWqV4d99w1LUdxh5crw2V+wLFq0+fK992DVqi33qVkzJMt77RWWguvNm4eH3dESjiQ0zomUn1GjNifHBdauDevL8Ie6bt062rRpA0C3bt0466yzeP/99/ntb3/Lhg0bOP7443+9X0REslscowxlf4IMkJMDgwfDiSfCbbfBX/4CL7wAl1wCV10FUSeK4piFuuW6deGgg4reZtWqkCwvWgRffLF5+c9/4P/9vy1boM1CnXPz5puXffaBvfcO12vWDNsluXFOpOyWLCnZ+h2UWINcoHv37rz99tu8/PLLDB06lEsvvZTBgweX6XFERCS9ubsS5KTbaSe47jo4++zQmnzzzfDQQ6FZ9swzIerZXBq1a0ObNmEpzB2++QY+/xz++98tlxdfhOXLt9x+t91Csvzhh0lpnBMpP3vuGb65FbW+nC1evJjGjRtzzjnnsH79embNmqUEWUQkyylBTqVGjeDhh2HECLj4YjjnHLjnHrjzTujWrdwfziwkvbvtBp07b33/jz+GZHnhQvjss82XxfVLXLwYLroI9tsvlIXstx80bKiaZ4nBTTdt+TMHhJqlm24q94d68803ufXWW8nLy6NGjRo88sgj5f4YIiKSXlSDHIf27eGdd+DJJ+H3v4fu3WHgwFCCkYQWsOLUrFl063PTpkU3zlWuDA8+uGVOUqsW7L//lkuLFuEYOZoORpKl4KeMci6ULxjBIlHhmc0KLFq0qEyPJSIi6UstyHExC0nxscfCrbeGznsvvABXXAGXXx5aw2JSXOPcmDFw2mmwbBl88klYFiwIyyuvwLhxm7evVm1zstyiBbRsCQceGMZ7VouzlItBg1TzIyIiSZNVCbKZVQXeBqpEj/O0u19nZs2AJ4A6wEzgDHf/JVlx7LDq1UN98plnhsT4hhvgH/+A22+H/v1jySa31zhXMNFJr15b7rdqVUiW583bvEydCo8+unmbWrU2J8utWkHr1uH6dvorioiIiKRMNrYgrwd6uvsaM8sD3jWzycClwB3u/oSZ3QecBdybxDhKZs894Ykn4He/CzXKAwaEoeLuuiuMrZxipWmcq10bOnUKS6Lvvw/J8kcfwccfh8sJE+D++zdv06xZSJYLSj7atlVrs4iIiMQj62qQPZxRQSFhXrQ40BM4LVr/MHA96ZQgFzj0UJg1K2SPV18dmllHjIDrr4edd447ulLZZZcwDHTXrpvXuYcZBefMCcvcueHyhRfCfQB16mxOlg86KJRuN2+upFlERESSKxtbkDGzSoQyir2BvwP/BVa5e360yVKgUTH7DgOGAVSuXDmZYRYvNxeGDw81yldfHUa5ePzxMJbyoEFZkSGahRbiPfaAY47ZvP6nn0Ky/OGHm5e77oJfomKY2rWhXbuQLLdvD4ccopZmERERKX9xJMhJHd/A3Te6exugMXAwsF8J9h3j7u3dvX1ubsx9CevWhfvug+nTw/zSZ5wRWpg/+ijeuJJop51Cicb558MDD8CMGWE4ulmzwu2BA2H1avjb38LEhE2ahGHsjj0W/vjHMDlK4dkFC4wfv3l0jaZNw20RERGRwrKyBbmAu68yszeATkBtM8uNWpEbA8tSEUO5aN8e/v3vMLnIlVeGmoMRI0KHvgrQu61y5XDKbduGuVYgtCjPnRu+O0yfDtOmwUsvhfvMwggaBfXQnTqFWb/PO0+zA8q2VapUiQMPPJANGzaQm5vL4MGDueSSS8jZxpiFixYt4phjjuHjjz9m9uzZ/O9//6Nv374pjFpERMpbXDXISWtBNrN6ZlY7ul4NOAJYALwBDIg2GwK8kKwYkiInJ2SHn34aLu+8M8zUMWHC5qLdCqRy5fC94fzzw/ByCxaEluMpU8L3hiZN4Nlnw1PVogUMHVr87IAiBQqmmp43bx5Tpkxh8uTJ3HDDDTu8/+zZs5k0aVISIxQRkVSIqwU5mSUWDYE3zGwu8AEwxd1fAq4ALjWzhYSh3sYmMYbkqVMnlF1Mmwa77w6nnAJHHRWmwKvgdt45DD13zTUwaRJ8+21InMeOhU2bit5n8WJ44w1Yty61sUrZJbtkpn79+owZM4bRo0fj7mzcuJHf//73dOjQgVatWnF/4jAswC+//MK1117LhAkTaNOmDRMmTGD69Ol06tSJtm3b0rlzZz799NPyDVJERJImq0os3H0u0LaI9Z8T6pGzQ4cOIUm+997QDNqyZSi/uOoqqFo17ujSQk5OaGTfbz+48caiZwcE6NkztEgfcgj06BFG1+vYEapUSWm4UgLjx285mU2ySmb22msvNm7cyPLly3nhhRfYeeed+eCDD1i/fj1dunThyCOP/PUNtHLlytx4443MmDGD0aNHA7B69WreeecdcnNzee211/jDH/7AM888U34BiohIUmRjC3LFUakSXHBBmNLuxBNDFtiyZagzkC3cdNPWkxMWzA740ktw4YXw889hux49YNddoXfvMMnhhx8W3wIt8Rg1KvUlM6+++iqPPPIIbdq04ZBDDmHlypV8tp1fbn744QdOOukkWrZsySWXXMK8efOSF6CIiJSbrBsHuUJq2BAeewx++9tQlHvkkaEZ7fbboX79uKNLC9ubHfDoo8PlqlXw1lthBsCpU8PkhgD16oWntXfvcKmnNV5LlpRsfWl9/vnnVKpUifr16+Pu3H333Rx11FFbbLNo0aJi97/mmms47LDDeO6551i0aBE9evQo3wBFRCQp1IKcTXr1CkM7XHstPPlkqC148EE1f0YGDYJFi8LTsWhR0T/F164N/fqFsZfnzYNly+CRR0KZ96uvhpH2GjQIE5dcfTW8/z5s3JjiExH23LNk60tjxYoVnHfeeVxwwQWYGUcddRT33nsvGzZsAOA///kPP/300xb71KxZkx9//PHX2z/88AONGoUh18eNG1d+wYmISNIpQc4mVauGYRzmzIEDD4RzzgljJ8+fH3dkGWn33UNS/Oij8PXXYVzmgnKNm28OQ8g1bBhGyXj66TBGsyRfcSUzN91UtuOuW7eONm3a0KJFC3r16sWRRx7JddddB8DZZ5/NAQccQLt27WjZsiXnnnsu+fn5W+x/2GGHMX/+/F876V1++eVcddVVtG3bdqttRUQkfcXVgmxx1XaUxE477eSFW4gyinsYA23kSFizBv7wh9CJT73PysV334WJSV56CSZPhu+/h7y80MnvhBNCS3SDBnFHmTkWLFjA/vvvv8Pbjx9ffMlMpinq3M1srbvvFFNIacPMqgJvA1UI5XlPu/t1hbYZCtzK5vHtR7v7g9s6bvv27X3GjBnlH7CIZIVWrVrRvHlznnvuuXI/tpnNdPf2Rd2nFuRUMIMzzwyd+AYMCC3LbdvCe+/FHVlW2HVXOPXUkKgtXw5vvx06+/3nP3DuuaFluWvXUApe3AgaUno7UjIjWWE90NPdWwNtgN5m1rGI7Sa4e5to2WZyLCKyPapBrgjq1w+d+CZNgp9+Clnb734HP/wQd2RZIzcXunWD226DhQtDhct114WG+5Ejwzi9nTqFKbKXLo07WpHM4cGa6GZetKT/T5AikvGUIFcUffqEnmcXXxzGN2vRYvP8zFJuzKBVq5Agz54dEuY//zkMI3fJJbDHHiGZvuceWLky7mjTSyaUXpW3injOJWVmlcxsNrCcMPnTtCI2O9HM5prZ02a2RzHHGWZmM8xsxooVK5IZsohkOLUgVzQ1asAdd4ThF3bZBY49Fk47DfRhkTTNm4fS7w8/DNUuN94Y6pWHDw9lGMcfD888A+vXxx1pvKpWrcrKlSsrVMLo7qxcuZKqmtxnm9x9o7u3ARoDB5tZy0KbvAg0dfdWwBTg4WKOM8bd27t7+3r16iU1ZhHJbOqktw0Z30lve375JQzF8Kc/hXma77orTF0dwx9ERTRnThgdY/z4MEJG7dowcCCcdRa0b1/xXoYNGzawdOlSfv7557hDSamqVavSuHFj8vLytlivTnpFM7NrgbXuflsx91cCvnP3nbd1HHXSE5FtOeCAA2jRogVPPfVUuR97W530NFFIOqhcOYyZ3L9/yMpOOy3UKt9/fxjfTJKqdeuw3HJLmJTkkUfCcv/9oUTjnHNCx7Nddok70tTIy8ujWbNmcYchacbM6gEb3H2VmVUDjgBuKbRNQ3f/Krp5HLAgxWGKSBZSiUVF17Il/Otf8Ne/hkytRYswPFwGtPJng0qVwux8//wnfPUV3HdfGC5uxIhQgnH66fDuu3o5ZEtm9pCZLTezjwutH2Fmn5jZPDP7S8L6q8xsoZl9amZHbX3EtNUQeMPM5gIfEGqQXzKzG83suGibC6PznQNcCAyNKVYRyRKqQZagUiW49NLNE4yceWaYecEsDMEwfnzcEVYIO+8choibMQNmzYKzzw79KLt1g3btYOxYWLs2bDt+fHhpcnL0ElVQ44DeiSvM7DCgH9Da3VsAt0XrDwBOAVpE+9wTlSKkPXef6+5t3b2Vu7d09xuj9de6+8To+lXu3sLdW7v7Ye7+SbxRi0imU4IsW9pnn/Dbfl5eGHYBwiC+55yjDCzF2raF0aPDdNdjxoQprc8+Gxo3hqOPDtcXLw4ty4sXw7BheokqEnd/G/iu0OrfATe7+/pom+XR+n7AE+6+3t2/ABYCB6csWBGRDBNXXzklyOnsmmtgw4Yt161bB1dcEU88FdxOO4XvJ3PmwFtvQa9eYUjrwn3Z1q4NM8tJ1sgtGJIsWobtwD6/AbqZ2TQze8vMOkTrGwFfJmy3NFonIiLFiKMFWZ300tmSJUWvX7YsDLtw+ukVb4iFNGAG3buHpbinXzP2ZZX84no5b0MusCvQEegAPGlme5V7ZCIiWU4lFrK1Pfcsen2VKjB4MJxwQhiXTGLTpEnR6/Py4Mknw/TLUiEtBZ6NZp+bDmwC6gLLgMTJMxpH60REpAhKkGVrN90UOuglql4dHnggzKX8yith5IskjA0oO6aol6hyZahXL4yl3KoVTJgQ6palQnkeOAzAzH4DVAa+BSYCp5hZFTNrBuwDTI8rSBGRdKcaZNnaoEGhV1iTJuG3/CZNwu0zzoCRI8OUcHvtBSefHLb9/vu4I65winqJHnooVMc88UTouHfKKWFAkieeUKKcjczsceDfwL5mttTMzgIeAvaKhn57AhgStSbPA54E5gOvAMPdXX8VIiLboJn0ipH1M+mVRX4+/PnPYd7k3XaDf/wDjjgi7qgksmkTPP003HADzJ8fEuVbboHevVU+nik0k15yaSY9EdmW5s2b07lzZx599NFyP/a2ZtJTC3Kmy80Ns/C9/z7UrBlmuhgxYvMgvRKrnJzQwP/RR/D44/DTT9C3Lxx+eBhjWURERIqnGmQpm/btw4wWF18cBu1t21YZWBrJyQmlFgsWwF13hYS5Qwc49VT4/PO4oxMREUlPqkGWsqtWDe64I0xTvXYtdOoUepGp8DVtVK4cGvj/+98wVvILL8D++4chr9XoLyIisjW1IEv56NkT5s6FE0+Eq6+GQw+FL76IOypJUKsW/OlP8NlnMGBAuH7AAfDcc6Fjn4iIiKjEQsrbLruEotd//jP8nt+6NTz8sLKvNNOoUZiW+q23QtLcv3/owPef/8QdmYiISPyUIEv5MwvjkM2dG2qShw4NhbAaDi7tdO8eSsjvvDP0t2zZEv74x61nGhcREalIVIMsydOkCbz+ehgO7tlnQ2vyO+/EHZUUkpsLF14YWo9PPDEMTtKhQxjuWkREpKJSC7IkT6VKcNVV8N57oadYjx4hA8vPjzsyKaRBg1Ad8/zz8M03IUm++mpYvz7uyERERFJLJRaSGgcfHJokzzgj/IbfrZvGGUtT/fqFyUVOPz0MRtKuHXzwQdxRiYiIpI4SZEmdmjVh3LjQTLlgQahPnjAh7qikCLvsEl6qSZNg9Wro3BluvTXM0CciIpLtVIMsqXfKKTB7dhhf7JRT4JxzNBhvmurTB+bMgeOOg8svDyNdfP113FGJiIgkl1qQJR5Nm8Lbb4f65LFjw4x8H30Ud1RShF13haefhvvuC30sW7eGV16JOyoREZHkUoIs8cjLCyNcvPoqfPdd6BV2770aMzkNmcG554ZZxOvXDy3Lv/+9+lqKiEh2UguyxK9Xr/A7fo8ecP75MHAg/PBD3FFJEVq0gOnT4bzz4LbbQsnFypVxRyUiIlK+VIMs6aFBg9Aj7Oabw5jJ7drBzJlxRyVFqFYtNPQ/9FAouTj4YPj447ijEhERKT9Z14JsZnuY2RtmNt/M5pnZRdH6681smZnNjpa+yYpBSiknB664Isx//MsvYeiE0aNVcpGmzjwzvFTr1kHHjuF7jYiISLbIqgQZyAdGuvsBQEdguJkdEN13h7u3iZZJSYxByqJLlzDKxRFHwIgRMGAArFoVd1RShI4dQ11yy5ZhFr7rrtNQcCIikvmyrgXZ3b9y91nR9R+BBUCjZD2eJEmdOjBxYhh8d+JEOOggzX2cpnbfHd58E4YOhRtvhCFDYMOGuKMSEREpvayuQTazpkBbYFq06gIzm2tmD5nZLsXsM8zMZpjZjHx10Y9XTg5cdln4HX/9eujUCcaMUclFGqpaNdQk/+lP8M9/wrHHwpo1cUclIiJSOlnXglzAzGoAzwAXu/tq4F6gOdAG+Ar4a1H7ufsYd2/v7u1zc3OTHabsiM6dQ+vxoYeGscYGD4affoo7KinEDEaNggcegClT4PDD4dtv445KRESkdLIuQTazPEJyPN7dnwVw92/cfaO7bwIeAA5OZgxSzurVC6Nc3HADjB8fhk5YsCDuqKQIZ58dOuzNnQtdu8LixXFHJCIiUjJZ14Js4WzGAgvc/faE9Q0TNjsB0MBUmaZSJbj22jCxyIoVIUl+6qm4o5Ii9OsXXqZvvgk/AGgYOBERySTZWIPcBTgD6FloSLe/mNlHZjYXOAy4JIkxSDL16hVKLg48EE4+GS69VL3C0lC3bmGcZICePWH+/HjjkcxkZlXNbLqZzYmG7ryhiG2qmNkEM1toZtOi/iciIqUWVwty0op73f1doKgz0rBu2aRRozB0wmWXwR13hLHGJkyAhg23u6ukTsuW4WU69NCQJL/5Juy3X9xRSYZZD/R09zVR+dy7ZjbZ3d9P2OYs4Ht339vMTgFuAQbGEayIZI+sKrGQCqRyZbjrrlCTPHNmmH2voMlS0sY++8Drr4fBR3r2hM8+izsiySQeFIyJkhcthX/77Ac8HF1/Gjjc4vhkE5GskXU1yFIBnXYaTJsGNWuGDOyuuzQUXJrZb7+QJG/YEF6izz+POyLJJGZWycxmA8uBKe4+rdAmjYAvAdw9H/gBqFPEcX4dxnPFihVJjlpEMlk21iBLRdSyJXzwAfTtCxddFGarWLs27qgkQYsW8Npr4WU57DCNbiE7LhqBqA3QGDjYzFqW8ji/DuNZr169co1RRLKLWpAle+y8Mzz3XJjO7Z//DFNWL1oUd1SSoHXrMEby6tWhJXn58rgjkkzi7quAN4Dehe5aBuwBYGa5wM7AypQGJyJZRwmyZI+cHLjmGnjxRfjiizBF9ZQpcUclCdq1g8mT4X//C8PBrVsXd0SSzsysnpnVjq5XA44APim02URgSHR9APC6x/X7qIhkBbUgS3Y6+uhQctGwIfTuDbfeqrrkNNKxY+hbOW0anHEGbNoUd0SSxhoCb0RDdH5AqEF+ycxuNLPjom3GAnXMbCFwKXBlTLGKSJZQDbJkr332gfffh/794fLL4fTT1VyZRvr3h9tug2eegSuVzmQ9M/tjVP5QcLuWmf1je/u5+1x3b+vurdy9pbvfGK2/1t0nRtd/dveT3H1vdz/Y3dUNVETKRC3Ikt1q1IAnn4SbboLHH4fatcEMmjYNTZgSq0sugfPPDw38990XdzSSZLnANDNrZWZHEFqDZ8Yck4hIsbJqohCRrZhBkyZh3OT168O6xYth2LBwfdCg+GKr4MzgzjtDX8rhw8PL1KdP3FFJMrj7VWb2GjAN+B7o7u4LYw5LRKRIakGWimHUqM3JcYG1a8N6iVVubpgEsVWrUHax++6hr6Ua+bOLmXUH7gJuBN4E7jaz3WMNSkSkGKpBlophyZKi1y9eHGavkFjVqAG//S38/DN89VXoT1nQyK8kOWvcBpzk7v/n7qcBDwCvxxyTiEiR1IIsFcOeexZ/X58+8N13qYtFivTXv269To38WaWTu88vuOHuzwJdYoxHRGSblCBL9rvpJqhefct11auHJsp33oFDDoEFC+KJTYDiG/mLWy8Zp7mZTTWzjwHMrBXwu5hjEhEpklqQpWIYNAjGjAm9wAo67Y0ZA/ffD2+8EaZ269gxzGAhsSiukX9bjf+SUR4ArgI2QBi+DTgl1ohERIqhGmSpOAYNCsMlbNoULgtGr+jcOUwq0qwZHHNMGFZBk4qkXFGN/GYqscgi1d19eqF1+bFEIiKyHWpBFoHQTPnee3DccXDxxWFwXnXeS6nCjfwNGoTRLCZP1veVLPGtmTUHHMDMBgBfxRuSiEjxlCCLAOy0U5jW7YorwqwVffvCqlVxR1WhJDbyf/01/OUv8NxzcPfdcUcm5WA4cD+wn5ktAy5GNcgikqbUgiySKCcHbr4ZHnoI3noLOnWChZrLIC6XXBKqXq64Aj77LO5opCzc/XN37wXUA/Zz967uvijmsEREihRXDbJm0pP0duaZsNdeYeaKQw6B55+Hbt3ijqrCMQv9KFu0gLPOgjffDN9hJHOY2aXFrAfA3W9PaUAiIjtALcgixTn0UJg2DerVg169NGNFTHbfHf72tzAa39//Hnc06cPMHjKz5QXDphW6b6SZuZnVjW6bmd1lZgvNbK6ZtUthqDWjpT2hpKJRtJwHpDIOEZESUYIsUpy994Z//SuUWpx+Otx4o3qMxWDw4DCfy5VXwn//G3c0aWMc0LvwSjPbAzgSSBxBug+wT7QMA+5NQXwAuPsN7n4D0Bho5+4j3X0kcBCgQfxEJC2pBVlke3bdFV59Fc44A667LpRf/PJL3FFVKGZhhIvcXDj77NCJr6Jz97eBoqaAvAO4nGi0iEg/4BEP3gdqm1nDFISZqAGQ+I/zS7ROREQiSpAls1SuDA8/DDfcEC6POgq+/z7uqCqUxo3DdNRvvhnqkiuAXDObkbAM294OZtYPWObucwrd1Qj4MuH20mhdKj0CTDez683semAa8HCKYxAR2a6CDnpxtCCrk55kHjO49trQee+ss8IEI5MnQ9OmcUdWYZx1Fjz5JFx+eRiFr0mTuCNKqnx3b7+jG5tZdeAPhPKKtOPuN5nZK0DXaNWZ7v5hnDGJiGyLSixESuL000PJxddfh+mpZ86MO6IKwwwefDBcP+cclYMX0hxoBswxs0WEmt9ZZrYbsAzYI2HbxtG6VJsNPAU8B6w0M9Ugi0jaibMFWQmyZLZDDw2d96pWhe7dYdKkuCOqMPbcMwxVPWUKPPts3NGkD3f/yN3ru3tTd29KKKNo5+5fAxOBwdFoFh2BH9w9pbPYmdkI4BtgCvAS8HJ0KSKSVuIaAxmUIEs22H9/+Pe/Yb/9whTVY8bEHVGFce650LIlXHYZ/Pxz3NHEw8weB/4N7GtmS83srG1sPgn4HFgIPACcn4IQC7sI2NfdW7h7K3c/0N1bxRCHiMg2qQVZpKwaNgwz7h15ZMjaRo3S7/4pkJsbxkZetAjuuCPuaOLh7qe6e0N3z3P3xu4+ttD9Td392+i6u/twd28eJaYzYgj5S+CHGB5XRKRE1ElPpDzUqAETJ8L558Of/wzLlsEDD0BeXtyRZbXDD4d+/cJTPnRo+K4iae1z4E0zexlYX7BSM+mJSLpSC7JIWeXmhrHHCoaB69cP1qyJO6qsd9ttsH59aLiXtLeEUH9cmc2z69WMNSIRkSLEWYOsFmTJPgXDwO2+eyi3OOwwePllqF8/7siy1t57w8UXh0R5+HA46KC4I5LiRLPpiYikPdUgiyTD2WfD88/DvHnQpQt8/nncEWW1q6+GevXgootU/p2OzOxv0eWLZjax8BJzeCIiW1ENskiyHHssvP46HHMMdOoUJhRp1y7uqLJSrVpw001hXOQnn4SBA+OOSAp5NLq8LdYoRERKSC3IIsnQsSO89x5UqwY9eoSEWZLizDOhdesww966dXFHI4ncfWZ0+VZRS9zxiYgUphILkWTbd98woUiTJtCnDzzzTNwRZaVKleDOO2HJEvj73+OORkREMllWThRiZnuY2RtmNt/M5pnZRdH6Xc1sipl9Fl3ukqwYRLaw++7w9tvQoQOcdFIY7ULK3aGHQq9ecOutsHZt3NFkPjPbx8yejt5LPy9YYoijyPf0Qtv0MLMfzGx2tFyb6jhFJHtkawtyPjDS3Q8AOgLDzewA4EpgqrvvA0yNboukxi67wKuvwtFHw3nnwR//qB5lSXDddbB8Odx3X9yRZIV/APcS3lMPAx4B/lnWg5pZ9RLuUtx7emHvuHubaLmxrHGKiGRVguzuX7n7rOj6j8ACoBHQD3g42uxh4PhkxSBSpOrV4dlnYciQMBzchRfCpk1xR5VVunaFnj3hL39RLXI5qObuUwFz98Xufj1wdGkPZmadzWw+8El0u7WZ3bO9/bbxni4ikhTZ2oL8KzNrCrQFpgEN3P2r6K6vgQbF7DPMzGaY2Yz8/PxUhCkVSV4e/OMfcNllMHo0DB4MGzbEHVVWue46+OYbGDMm7kgy3nozywE+M7MLzOwEoEYZjncHcBSwEsDd5wDdS3KAQu/phXUyszlmNtnMWhSz/6/v7ytWrChR8CJScWRlDXIBM6sBPANc7O6rE+/zcOZFnr27j3H39u7ePjdXo9FJEpiFJs4//xnGj4f+/dXcWY66dw+DhtxyC/z8c9zRZLSLgOrAhcBBwOnA4LIc0N2/LLRq447uu633dGAW0MTdWwN3A88X8/i/vr/Xq1dvxwMXkQolbVuQzaxWdLlrUcv2Dm5meYQ30vHu/my0+hszaxjd3xBYXrZTECkDM7jqKrj33jDbXp8+sLrwZ76U1rXXwldfwQMPxB1JRmvq7mvcfam7n+nuJwJ7luF4X5pZZ8DNLM/MLiOUS2xXMe/pv3L31e6+Jro+Ccgzs7pliFVEJP0SZOCx6HImMCO6nJlwu1gWzmYssMDdb0+4ayIwJLo+BHihhDGLlL/zzgutyO+9F6am1s++5aJHD+jWDW6+Wa3IZXDVDq7bUecBwwn1w8uANtHtbdrGe3riNrtF22FmBxM+Y1aWIVYRqcDSdiY9dz8mumxWimN3Ac4APjKz2dG6PwA3A0+a2VnAYuDkUhxbpPydemqYDm7AgJDVvfYaNG4cd1QZzSzUIvfqBWPHwvDtpmFSwMz6AH2BRmZ2V8JdtQgjSpT60O4+qBT7FfeevieAu98HDAB+Z2b5wDrgFI+ziFBEMlqcbx87VNxrZl2A2e7+k5mdDrQD/ubuS4rbx93fBYpL+Q8vcaQiqXD00WEYuGOO2ZwkN28ed1QZrWdP6NIltCKffTZUqRJ3RBnjf4Rf646LLgv8CFxShuO+Z2aLgAnAM+6+akd22s57esE2o4HRZYhNRORXaVuDnOBeYK2ZtQZGAv8FHk1aVCJx6tYtTEe9enW4Pn9+3BFltIJW5KVL4aGH4o4mc7j7HHcfBzR394cTlmfd/fsyHPc3wNVAC2CWmb0UNXyIiKSldE6Q86OfyfoBo93970DN5IUlErODDoK33gqTiBx6KMyaFXdEGa1XLzjkELj9dg05vaPM7CMzm0tIYucWXspybHef7u6XAgcD37F5bHoRkbSRtjXICX40s6sIwwt1j8bkzEteWCJpoGVLeOcdOPzw0HFv8mTo3DnuqDKSGVx0EZx2GrzyCvTtG3dEGeGYZBw0Gp3oBOAUoDnwHCFRFhFJK5kwDvJAYD1wlrt/DTQGbk1aVCLpYu+94d13oUEDOOIImDo17ogy1oknQsOGcPfdcUeSGaJZ8xa7++Jo1T7R9eWEVt/SmkMYueJGd/+Nu1/h7jO3s4+ISMqlbQ2ymVU1s4sJ9Wo/Af8GcPcl7v5I8sMTSQN77BFakps3D534Jk2KO6KMVLlyGE3vlVfg00/jjiZzmNk5wNPA/dGqxhQzAccO2svdL3H3f5c1NhGRVEi7BJlQl9Ye+AjoA/w16RGJpKMGDeCNN6BFCzj+eHj++bgjykjnnhtm+f773+OOJKMMJwyxthrA3T8D6pf0IGb2t+jqRDPbaim3aEVEykk61yAf4O4HApjZWGB68kMSSVN16oQSi969w1jJ48fDwIFxR5VRGjQIT9k//gF/+lMYdlq2a727/1LwAWFmuUBpCvMKRh66rbwCExFJpnSuQd5QcMXdyzIwvUh2qF0bpkwJnfVOOw0eUaVRSV14IaxZAw9r3IQd9ZaZ/QGoZmZHAE8BL5b0IAl1xm3c/a3EhVCTLCKSVtK2BhlobWaro+VHoFXBdTNbnYoARdJOzZphRIuePWHoUBgzJu6IMkqHDmHIt9GjNeTbDroSWEEodTsXmEToF1JaQ4pYN7QMxxMRSaq0K7Fw90qpCkQko+y0E7z4Yhia4dxzIT8fzj8/7qgyxoUXwqBBYdLC3r3jjia9ufsmM3seeN7dV5T2OGZ2KnAa0KxQzXFNyjYqhohIUqRzC7KIFKdqVXj2WWjbFoYPD4P9Nm0aapNlmwYMgN1205Bv22LB9Wb2LfAp8KmZrTCza0t5yH8ROlp/El0WLCOBo8ojZhGR8pTONcgisi1PPw2ffLL59uLFMGyYkuTtKBjybdIk+OyzuKNJW5cQRq/o4O67uvuuwCFAFzO7pKQHi8ZUfhMYBExLqD9eQBg6TkQkragFWSRTjRoF69ZtuW7t2rBetklDvm3XGcCp7v5FwQp3/5wwo+ngMhz3SSCx+nsjoeOfiEhaUoIskmmWLCl6/eLFRa+XX+22G5x0Ejz0UPhOIVvJc/dvC6+M6pDzynDcXHf/JeF4vwCVy3A8EZGkUAuySKbac8/i7/vzn1MXR4YaNgx+/BGeey7uSNLSL6W8b3tWmNlxBTfMrB+wVSIuIhI31SCLZKqbboLq1bdcV60adOkSyiz+7//iiStDdOsW+jVqOOkiJQ6zubrQkJsHluG45wF/MLMlZvYlcAVh+DgRkbSiFmSRTDVoUBgHuUmTMIpFkybwwAPw1lvhvj/8AW6+Oe4o01ZODpxxBrz2GixbFnc06cXdK7l7rSKWmu5e6hILd/+vu3cEDgD2d/fO7r6w/CIXESlfaTcOsojsgEGDwlLYww+DO1x1VUier7gi9bFlgDPOgD/+MQz8cfnlcUdTMZjZ0UALoGrBB4+73xhrUCIihagFWSQbVaoUkuTTToMrr4Rbbok7orS0zz5h5u6C7xOSXGZ2HzAQGAEYcBLQJNagRESKoBpkkWyVmxsyv1NPDUnyrbfGHVFaGjwY5s+HWbPijqRC6Ozug4Hv3f0GoBPwm5hjEhHZilqQRbJZbm7ohXbKKaGG4I474o4o7Zx8MlSpEr5LSNL9HF2uNbPdgQ1AwxjjERHZJiXIItkqNxcefTTMsXzppTB6dNwRpZVddoHjjoPHH4dfyjKAmeyIF82sNnArMAtYBDwWZ0AiIkVRC7JIRZCbC489BscfDyNGwP33xx1RWhkyBL79FiZPjjuS7GVmOcBUd1/l7s8Qao/3c/drYw5NRGQrqkEWqSjy8mDCBDjmGDjvPBg7Nu6I0saRR0L9+hoTOZncfRPw94Tb6939hxhDEhEpllqQRSqSypXh6aehd2845xwV3kby8sJoeS++CCtXxh1NVptqZidaHJ84IiIloARZpKKpUgWefRYOPxx++1t44om4I0oLgwfDhg2hkV2S5lzgKWB9wcx8ZrY67qBERIqjBFmkIqlWDV54Abp2hdNPh+eeizui2LVpA61aZVajupk9ZGbLzezjhHW3mtknZjbXzJ6LOsUV3HeVmS00s0/N7KgUxtklulrP3XPcvXLCzHy1UhWHiMiOUg2ySEVVvTq89BJ06AADB8KkSXFHFLvBg2H6dPjkk7gj2WHjgN6F1k0BWrp7K+A/wFUAZnYAcAphFrvewD1mVilFcd4VXf4rRY8nIlImKrEQqchq1gxDNxx4IPTvD1Onxh1RrE47LczMnSllFu7+NvBdoXWvunt+dPN9oHF0vR/wRNQ57gtgIXBwikLdYGZjgMZmdlfhJUUxiIjsMCXIIhVd7drw6qvwm9+EAYHfeSfuiGLTsCF06RJKtNNErpnNSFiGlXD/3wIFg9c1Ar5MuG9ptC4VjgFeB9YBM4tYtsnM9jCzN8xsvpnNM7OLitjGooR7YVRe0q6cz0FEKqA4EuTclD+iiBStTh2YMgV69ICjj4bXXoODU9W4mF769w/zqSxcCHvvHXc05Lt7+9LsaGajgHxgfPmGVHLu/i3whJktcPc5pThEPjDS3WeZWU1gpplNcff5Cdv0AfaJlkOAe6NLEZESUw2yiAQNGoTEuF69MAzc3LlxRxSL/v3DZSb3WzSzoYRW20G++V1+GbBHwmaNo3WptM7MphZ0KjSzVmZ29fZ2cvev3H1WdP1HYAFbt373Ax7x4H2gtplpGmsRKRWVWIjIZo0ahTrk6tXhiCPg00/jjijlmjSBgw6CZ56JO5LSMbPewOXAce6+NuGuicApZlbFzJoRWlqnpzi8BwidBjcAuPtcQsfBHWZmTYG2wLRCd+1QCYmZDSsoWVmxYkVJHlpEKhAlyCKypaZNN3fW69ULFi2KM5pY9O8P06bB0qVxR7JtZvY48G9gXzNbamZnAaOBmsAUM5ttZvcBuPs84ElgPvAKMNzdN6Y45OruXjgpzy9yyyKYWQ3gGeBidy/V+MnuPsbd27t7+3r16pXmECJSgShBFpHN9t03dNxbsyZMKPK//8UdUUqdeGK4fP75WMPYLnc/1d0bunueuzd297Huvre77+HubaLlvITtb3L35u6+r7tP3taxk+RbM2sOOICZDQC+2pEdzSyPkByPd/eiulGmQwmJiGSJrKxBLmbw/OvNbFnUojLbzPom6/FFskLr1vDKK7B8eWhJrkA/R++7LxxwQFqNZpEthgP3A/uZ2TLgYuC8be5BGKECGAsscPfbi9lsIjA4Gs2iI/CDu+9Q8i0iUli2lliMY+vB8wHuSGhV0awIIttzyCFhMpEvvggd9374Ie6IUqZ/f3jrrQr1vSDp3P1zd+8F1AP2Aw4Fuu7Arl2AM4CeiY0cZnaemRUk2JOAzwnjOz8AnF/+ZyAiFUVWJshFDZ4vIqV06KGhx9rcuXDMMbB27fb3yQInngibNsHEiXFHkvnMrFY0zfVoMzsCWAsMISSzJ29vf3d/193N3VslNnK4+33uXlBj7e4+PCohOdDdZyT3rESkIsiqBHkbLogGkH/IzHYpbqPEXs75+Tvcf0Qke/XtC+PHw3vvhczxl1/ijijpWreGZs1UZlFOHgX2BT4CzgHeAE4CTnD3fnEGJiJSlKysQS7GvUBzoA2hU8hfi9swsZdzbq7mMxEB4OSTYcyYUJd8+umwMdUDIKSWWSizeO21ClVZkix7uftQd78fOBU4ADjK3WfHG5aISNGyssSiKO7+jbtvdPdNhPq0ijlNmEhZnH023HYbPPUUnHsuxPgNOxX69w+N5S+/HHckGW9DwZVoaLml7v5zjPGIiGxTnAlySptmzaxhQo/mE4CPt7W9iBRj5MjQpPrHP0KtWvDXv4bm1izUsSM0bBjKLE47Le5oMlprMysYt9iAatFtI5QP14ovNBGR4mVVghwNnt8DqGtmS4HrgB5m1oYw/uYi4NxkPb5I1rvhBli1Cu64A+rUgVGj4o4oKXJy4IQTYNy40DexevW4I8pM7l4p7hhEREoizhrkpCXI7n5qEavHJuvxRCocM/jb30KSfPXVsMsucH52jqrVvz/cc0+YN+X44+OORkREUqHC1CCLSDnLyYGxY+HYY+GCC+Cxx+KOKCkOPRR23VWjWYiIVCRKkEWk9PLyYMIE6N4dhgzJyt5sublhlLvJk8O4yCIiUnEoQRaR0qlWLcym0aoVDBgA77wTd0Tlrk8f+PZbmDkz7khERCQVKtI4yCKSLLVqhfGRmzQJs+19+GHcEZWrI48MZdeTJ8cdiYiIpIJKLESkfNSrB1OmwM47Q+/e8NlncUdUburWhQ4dwncAERHJfkqQRaT87LFHGO5h0yY44ghYtizuiMpN794wbRp8913ckYiISKooQRaR8rHffqEWYeVKOOqorMko+/QJef+UKXFHIiIiyaYaZBEpf+3bwwsvhDKLo4+Gn36KO6Iy69AhDPemOmQRkeynEgsRSY6ePeHxx2H6dDjxRPjll7gjKpNKlUJnvVde0XBvIiLZTgmyiCRP//4wZgz8v/8XxknO8MyyTx/45huYMyfuSEREJBXiSJCTNtW0iKSRs84KgwhfeSXUqQN33x3GTMtARx0VLidPhrZt441FRESSRzXIIpJ8l18OI0fC3/8ON94YdzSl1qABtGun4d5ERLKdSixEJPnM4NZbQ5nF9dfDPffEHVGp9e4N//oXrFoVdyQiIpIsSpBFJDXM4MEH4dhj4YIL4Mkn446oVPr0gY0b4bXX4o5ERESSTQmyiCRfbi5MmABdu8Lpp4dJRTJMx45hskCVWYiIZC+1IItIalWrBhMnwgEHhFEupk2LO6ISyc0NkwS+8grE2IdDRESSSJ30RCT1atcOGWaDBtC3LyxYEHdEJdKnT5hF+6OP4o5ERESSQS3IIhKP3XYLJRZ5eWEGji+/jDuiHVYw3JvKLEREspsSZBFJvebNQ5a5enVIkr/9Nu6IdkijRtCqlaadFhHJVmpBFpF4tWkDL74IX3wBRx8Na9bEHdEO6d0b3n0Xfvwx7khERKS8qQZZROLXvXsY3WLGDOjUCZo0gZwcaNoUxo+PO7oiHXEE5OeHJFlERLKLWpBFJD306xempf74Y1iyJAwRsXgxDBuWlkly586hfPqNN+KOREREypsSZBFJH0WNi7x2LYwalfpYtqN69dDY/frrcUeS/czsITNbbmYfF3N/DzP7wcxmR8u1qY5RRLKTEmQRid+SJSVbH7OePWHWLPj++7gjyXrjgN7b2eYdd28TLTemICYRyWKqQRaR9LHnniVbH7OePUMlyFtvxR1JdnP3t4Hv4o5DRCqOOEssclP+iCKS3m66KdQcr1275fre22s8jMfBB4eJAd94A44/Pu5oKrxOZjYH+B9wmbvPK2ojMxsGDAPYM0lfvNavX09+fn5Sji0iqaEEWUTSx6BB4XLUqFBWscceUKsWPPhgmHHvuOPija+QKlWga1fVIaeBWUATd19jZn2B54F9itrQ3ccAYwDat29f7r+hTp8+nS5duihBFslwLVq0AJQgi0i6GDRoc6IMYVzkww+Hk08OM3Mcdlh8sRWhZ0+46ipYvhzq1487morJ3VcnXJ9kZveYWV13T/nMM4sWLSI/P5+RI0fSoEGDVD+8iJSDZ555hmnTpsX2+EqQRWT7atSASZPg0ENDC/LUqaG2IU0U5OtvvhlyeEk9M9sN+Mbd3cwOJvRxWRlHLBs3bgTgnHPOYd99940jBBEpo/r16/+aIFevXj3lj68EWUR2TJ06YQi4bt2gTx94+22Ifv6K20EHQc2aocxCCXJymNnjQA+grpktBa4D8gDc/T5gAPA7M8sH1gGneExd0AtKK3Jz9REnkqlOO+00dt55Z6pUqULr1q1T/vh69xCRHbf77jBlSij6PeKIMIXdXnvFHRW5uaFxW3XIyePup27n/tHA6BSFs01KkEUyX15eHsfH2PNaw7yJSMnstVdIktevh1694H//izsiINQhf/YZLF0adyQSt4ISi0qVKsUciYhkKiXIIlJyLVqEznorVoSW5G9T3g9rKwV1yJp2WtSCLCJlpQRZRErn4INh4kT4/HM48khYtSrWcFq1gl13VZmFKEEWkbJTgiwipXfYYfDcc/Dxx6Hj3o8/xhZKTk4I5/XXw8x6qWJmD5nZcjP7OGHdrmY2xcw+iy53idabmd1lZgvNbK6ZtUtdpBVHQYKsEgsRKa2kJcgl+dAQkQzWuzc8+SR88AEce+zWM/Cl0GGHhblNvvgipQ87Dig8zeCVwFR33weYGt0G6EOYPGMfwkxy96YoxgqloAZZLcgiUlrJbEEex45/aIhIJjv+ePjnP8PQbyecEDrwxaBnz3CZyjILd38b+K7Q6n7Aw9H1h4HjE9Y/4sH7QG0za5iSQCsQlViISFklLUEu4YeGiGS6U06BsWPDWMknnwwbNqQ8hP32g912K/cEOdfMZiQsw3Zgnwbu/lV0/WugYDq3RsCXCdstjdZJOVKCLCJllep3j+I+NLYSfQgNA6hcuXIKQhORMjvzTFi3DoYPh4EDYcIEyMtL2cObhVbkqVNDHbJZuRw2393bl3bnaGa5WCbMqKg0zJuIlFVsnfSiGZaK/dBw9zHu3t7d26sVQCSDnH8+3Hln6Lw3cCD88ktKH/6ww+Cbb+CTT1L6sIV9U1A6EV0uj9YvA/ZI2K5xtE7KUX5+PmZGTo76oYtI6aT63aO4Dw0RySYXXrg5ST7llJQmyXHUIRdhIjAkuj4EeCFh/eBoNIuOwA8Jv6pJOcnPz1d5hYiUSaoT5OI+NEQk21x4Idx1V8pbkps1g0aN4J13UvJwmNnjwL+Bfc1sqZmdBdwMHGFmnwG9otsAk4DPgYXAA8D5qYmyYsnPz1d5hYiUSdK+YkcfGj2Auma2FLiO8CHxZPQBshg4OVmPLyJpYMSIcHnhhZtrkpPcp8AMuneHt94q1zrkYrn7qcXcdXgR2zowPLkRycaNG9WCLCJlkrR3kJJ8aIhIFhsxImSpI0bAgAFhzOSqVZP6kN26weOPh0n+mjdP6kNJGlKJhYiUlXowiEjyXXAB3HMPvPRSmHFv9eqkPlz37uEyVWUWkl6UIItIWSlBFpHU+N3vwmQi77wDhx8O336btIfaf3/YdVclyBXVxo0bVYMsImWiBFlEUue00+D55+Hjj0Mz77LkjHCWkwNdu4aJ/aTiUQuyiJSVEmQRSa1jjoFXXoGlS0MWu3BhUh6me/dw6K80iFqFowRZRMpKCbKIpN6hh8Ibb8CPP4Ykefr0cn+Ibt3CpcosKh4N8yYiZaUEWUTicdBB8O67UL16SJiffLJcD9+2Ley0kxLkikjDvIlIWSlBFpH47LcfTJsWkuWBA+FPfwqDF5eDvDzo1El1yBWRSixEpKyUIItIvOrVg6lT4fTT4ZprYPBgWL++XA7dvTt89BF8/325HE4yhBJkESkrJcgiEr8qVeCRR+CPfwxDwR1+OHzzTZkP261baJB+771yiFEyhoZ5E5GyUoIsIunBDK6+OkxHPXMmtGoVRrsog0MOCaUWqkOuWNSCLCJlpQRZRNLLySfDBx9A/fph1r1LLy11yUW1atChgxLkikYJsoiUlRJkEUk/LVuGod+GD4c77gi97T79tFSH6t495Ntr15ZzjJK2NMybiJSVEmQRSU/VqsHo0fDCC7BkCbRrB3//O+Tnl+gw3bqFXaZNS1KcknY0zJuIlJUSZBFJb8cdB3PmQOfOcMEF0Lp1iWqTu3QJ5c0a7q3iUImFiJSVEmQRSX+NGsGrr8Jzz8Evv4Ta5N694eOPt7vrzjuHnFp1yBWHEmQRKSslyCKSGczg+ONh3jy4/fZQM9G6NZxzTigy3sYEI927w7//DRs2pC5ciY+GeRORslKCLCKZpXJluOQSWLgwlFw8+igcfDDsv3+YiW/Roq126dYtdNKbNSv14WYLM3vIzJabWZHN9hbcZWYLzWyumbVLdYwF1IIsImWlBFlEMlOdOnDnnfD11/DAA9CgQZiJr1mzkBFfckm4/4UX6LbrPEBlFmU0Dui9jfv7APtEyzDg3hTEVCQlyCJSVkqQRSSz1a4NZ58Nb70FX3wRWpHXrIExY+Dii+H442lweEt+w6e8fcu/4442Y7n728B329ikH/CIB+8Dtc2sYXnH0bVrV3Jycra5zJs3j8qVK5f3Q4tIBaKv2CKSPZo2hVGjwuIOK1eGkotFi/j9hEVU2aVa3BFms0bAlwm3l0brviq8oZkNI7Qys+eee5boQYYMGULPnj23u92JJ55YouOKiCRSgiwi2ckM6tYNS/v2nD0g7oCkgLuPAcYAtG/fvvjelUU455xzkhKTiEgilViIiEh5WAbskXC7cbRORCTjKEEWEZHyMBEYHI1m0RH4wd23Kq8QEckEKrEQEZHtMrPHgR5AXTNbClwH5AG4+33AJKAvsBBYC5wZT6QiImWnBFlERLbL3U/dzv0ODE9ROCIiSaUSCxERERGRBEqQRUREREQSKEEWEREREUmgBFlEREREJIGFfhXpzcw2AetKuFsukJ+EcOKSbecD2XdO2XY+kH3nVJrzqebuakxIEjNbASxOwqHrAt8m4bhlkW4xpVs8kH4xpVs8kH4xpVs8sOMxNXH3ekXdkREJcmmY2Qx3bx93HOUl284Hsu+csu18IPvOKdvOR4qXjq91usWUbvFA+sWUbvFA+sWUbvFA+cSkVhERERERkQRKkEVEREREEmRzgjwm7gDKWbadD2TfOWXb+UD2nVO2nY8ULx1f63SLKd3igfSLKd3igfSLKd3igXKIKWtrkEVERERESiObW5BFREREREpMCbKIiIiISIKsSZDNbJGZfWRms81sRrRuVzObYmafRZe7xB3njirmfK43s2XRutlm1jfuOHeUmdU2s6fN7BMzW2BmnTL59YFizykjXyMz2zch5tlmttrMLs7U12gb55ORr4+UjJn1NrNPzWyhmV0Zw+PvYWZvmNl8M5tnZhdF62P9fzKzSmb2oZm9FN1uZmbToudpgplVTnE8afe5YGaXRK/Zx2b2uJlVTfXzZGYPmdlyM/s4YV2Rz4sFd0WxzTWzdimK59bodZtrZs+ZWe2E+66K4vnUzI4q73iKiynhvpFm5mZWN7pdqucoaxLkyGHu3iZh7Lsrganuvg8wNbqdSQqfD8Ad0bo27j4ptshK7k7gFXffD2gNLCDzX5+izgky8DVy908LYgYOAtYCz5Ghr9E2zgcy8PWRHWdmlYC/A32AA4BTzeyAFIeRD4x09wOAjsDwKIa4/58uYvP7FMAthP+HvYHvgbNSHE9afS6YWSPgQqC9u7cEKgGnkPrnaRzQu9C64p6XPsA+0TIMuDdF8UwBWrp7K+A/wFUA0d/5KUCLaJ97ov/JVMSEme0BHAksSVhdquco2xLkwvoBD0fXHwaOjy+UisvMdga6A2MB3P0Xd19FBr8+2zinbHA48F93X0wGv0YJEs9Hst/BwEJ3/9zdfwGeIPwdp4y7f+Xus6LrPxISv0bE+P9kZo2Bo4EHo9sG9ASejimedP1cyAWqmVkuUB34ihQ/T+7+NvBdodXFPS/9gEc8eB+obWYNkx2Pu7/q7gWzkr4PNE6I5wl3X+/uXwALCf+T5aqY5wjgDuByIHEEilI9R9mUIDvwqpnNNLNh0boG7v5VdP1roEE8oZVKUecDcEH0E8FDmfJzN9AMWAH8I/p570Ez24nMfn2KOyfIzNco0SnA49H1TH6NCiSeD2T+6yPb1gj4MuH20mhdLMysKdAWmEa8/09/IyQOm6LbdYBVCUlOqp+ntPtccPdlwG2E1sevgB+AmcT7PBUo7nlJh7/33wKT447HzPoBy9x9TqG7ShVTNiXIXd29HaEpfbiZdU+808N4dpk0pl1R53Mv0BxoQ/jn/Wt84ZVILtAOuNfd2wI/Uehnswx8fYo7p0x9jQCIauuOA54qfF8GvkZFnU9Gvz6SWcysBvAMcLG7r068L5X/T2Z2DLDc3Wem4vF2UNp9LkRfmPsRkvfdgZ0o4mf8uKXTe7GZjSKUFI2POY7qwB+Aa8vrmFmTIEff/HD35YRaw4OBbwqa0aPL5fFFWDJFnY+7f+PuG919E/AASfjZIkmWAkvdfVp0+2nCG2PGvj4Uc04Z/BoV6APMcvdvotuZ/BpBofPJgtdHtm8ZsEfC7cbRupQyszxCcjze3Z+NVsf1/9QFOM7MFhFKTnoS6n9rR6UEkPrnKR0/F3oBX7j7CnffADxLeO7ifJ4KFPe8xPb3bmZDgWOAQb55Uo244mlO+GIzJ/o7bwzMMrPdShtTViTIZraTmdUsuE4o0P4YmAgMiTYbArwQT4QlU9z5FKqZOYFwjmnP3b8GvjSzfaNVhwPzydDXB4o/p0x9jRKcypblCBn7GkW2OJ8seH1k+z4A9olGHqhMKLGZmMoAovrescACd7894a5Y/p/c/Sp3b+zuTQnPx+vuPgh4AxiQ6niimNLxc2EJ0NHMqkevYUFMsT1PCYp7XiYCg6ORGjoCPySUYiSNmfUmlOwc5+5rC8V5iplVMbNmhI5x05Mdj7t/5O713b1p9He+lNBo9TWlfY7cPeMXYC9gTrTMA0ZF6+sQent+BrwG7Bp3rGU8n0eBj4C50QveMO5YS3BObYAZUezPA7tk6uuznXPK5NdoJ2AlsHPCuox9jYo5n4x9fbSU6LXvS+hZ/9+C988UP35Xwk/gc4HZ0dI3Hf6fgB7AS9H1vQjJy0JCGVKVFMeSdp8LwA3AJ4Qvz48CVVL9PBG+1H8FbCAkemcV97wARhi15b/Re1v7FMWzkFDXW/D3fV/C9qOieD4F+qTqOSp0/yKgblmeI001LSIiIiKSICtKLEREREREyosSZBERERGRBEqQRUREREQSKEEWEREREUmgBFlEREQkA5jZLmb2mplNiabrliTRKBYiIiIiGcDMBhBm+TPgS988CY2UM7UgS0Yws+PNzM1sv0Lru5rZdDP7JFqGJdx3vZldVszxLjazwdH13c3sdTN7IZoaFjO7wMx+m8xzEhEpKTPbaGazzWyemc0xs5FmlhPd197M7trGvk3N7LTURbvV419oZgvMLNZpiROZWQ8ze6kE2w81s91L8TjHm9kBpdwvcfrkt4BrouXNhO30mVXOlCBLpjgVeDe6BCCaQvIx4Dx3348wOP+5Znb0tg4UTRn622hfgAuBEcCDwOnRuoeidSIi6WSdu7dx9xbAEYTp1K8DcPcZ7n7hNvZtCsSWIAPnA0d4mMXvVwnTOGeCoYQW3K2YWaVt7Hc8UOIEmTBb3T0JtzcC3wOrCK3IBfSZVc6UIEvai1p1uxJm7zkl4a7hwDh3nwXg7t8S3kyu3M4hewKz3D0/ul0J2BQtFh1rLbDIzA4ur/MQESlP7r4cGAZcEE2j+2trqJkdGrU0zzazD82sJnAz0C1ad0nUovyOmc2Kls7Rvj3M7E0zezr6ZW58NPUyZtbBzP4VtV5PN7OaZlbJzG41sw/MbK6ZnVs4VjO7jzAj3eTosa83s0fN7D3g0SiW16P9p5rZntF+48zsXjN738w+j2J7KGqJHlfU82JmN5vZ/OhYtyUcZ0DCNmsSdqllZi+b2admdp+Z5UTnNM7MPjazj6KYBwDtgfHRc1jNzBaZ2S1mNgs4yczOiZ6HOWb2jIVpqzsDxwG3Rvs1j5ZXzGxm9Bps8etoFONvgPXRZ1uBk4FngWeAgQl/C/rMKmeZ9K1NKq5+wCvu/h8zW2lmB7n7TKAF8HChbWdE67elCzAz4fZowpSiP7Bl68oMoBspmEdeRKQ03P3zqOWyfqG7LgOGu/t7USPDz4TGg8vc/RgAM6tOaNH92cz2IUzf2z7avy3hvfR/wHtAFzObDkwABrr7B2ZWC1hHaLz4wd07mFkV4D0ze9Xdv0iI8zwz6w0c5u7fmtn1hBbVru6+zsxeBB5294ejUoG7CK2uEKag7kRIMicS3sPPBj4wszbuPrvgccysDnACsJ+7u5nV3oGn8eAolsXAK0B/4Augkbu3jI5b291XmdkF0XM4I1oPsNLd2xU8vrs/EF3/E2EK5LvNbCJhiu+no/umEn79/MzMDiG0EvcsFFcXYFahdacTGocKpk9ObF3WZ1Y5UoIsmeBU4M7o+hPR7ZnFb75dDYEFBTfcfTHQvYjtlgNbfasXEckA7wG3W6j3fdbdl0bJXKI8YLSZtSH8dP+bhPumu/tSADObTSjP+AH4yt0/AHD31dH9RwKtElpodwb2ISSZ2zLR3ddF1zsRElMIDRZ/SdjuxSjZ/Qj4xt0/ih53XhTX7IRtfyB8GRgbtabvSH3xdHf/PDrm44RfLKcCe5nZ3cDLwKvb2H9CwvWWUWJcG6gB/L/CG0dfWDoDTyW8JlWKOG5DYEXCfk2B2u4+J7q9i5k1S/gios+scqQEWdKame1K+FZ9oJk5oRzCzez3wHzgIOCFhF0OAuZt57DrgKo78PBVo21FRNKSme1FSG6XA/sXrHf3m83sZaAvoUX3qCJ2vwT4BmhNKLn8OeG+9QnXN7LtfMGAEe6+VTK4HT/t4HYFsWwqFNemwnG5e35UZnA4MAC4gPAZkk9UVmqhU2PlxN0KPZ67+/dm1ho4CjiPUNpQXCe4xPMYBxzv7nPMbCjQo4jtc4BV7t6mmOMVWEf4slFgENDIzBZFt3cm/Op5U3Rbn1nlSDXIku4GAI+6exN3b+ruexBaJboRfl4aGrV+FPy0dgtbtjwUZQGw9w489m+Aj0sbuIhIMplZPeA+YLQXGrPVzJq7+0fufgvwAaFl8UegZsJmOxNahDcBZxAaILblU6ChmXWIHqOmhQ52/w/4nZnlRet/Y2Y7lfB0/sXmPiaDgHdKuD/RY9cAdnb3SYQvAK2juxYRGlAglGrkJex2sJk1ixLngcC7ZlYXyHH3Z4CrgXbRtoWfw8JqAl9Fz0ViZ8Rf94ta3r8ws5OimC1Kxgsr/Fk1CDg4+ixsGp1P4mPoM6scKUGWdHcq8Fyhdc8Ap7r7V4R6rAfM7BPCG+xD7v7ido45maJLKgrrAkwpYbwiIslULeroNQ94jfDT/w1FbHdx1MFsLrCB8L43F9gYdSC7hFC/OsTM5hAS6G226Lr7L4QE8u5onymEVssHCb/ozTKzj4H7Kfkv1COAM6N4zwAuKuH+BWoCL0XHeRe4NFr/AHBoFHcntjzXDwh9URYQGmCeAxoBb0blJf8Eroq2HQfcV9BJr4jHvwaYRihx+SRh/RPA7y10mGxOSGzPiuKZR+hrU9jbQNsogW4HbHL3zwrujMpCNkT3gT6zypUmCpEKycyeAy5PfLMpdH9b4FJ3PyO1kYmIiARmdiehBvu17Wynz6xyphZkqaiuJHSAKE5dQkuAiIhIXP4MVN+B7fSZVc7UgiwiIiIikkAtyCIiIiIiCZQgi4iIiIgkUIIsIiIiIpJACbKIiIiISAIlyCIiIiIiCZQgi4iIiIgk+P+1o7Y/G+T4VQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n",
    "\n",
    "plot_ellipsdata(ax[0], data=data, model=model, xaxis=\"aoi\")\n",
    "plot_structure(ax[1], objective=objective)\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Saving the objective"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you would like to save the `Objective` or model to a file, this is best done through serialisation to a Python pickle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "pickle.dump(objective, open(\"my_objective.pkl\", \"wb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can then simply reload your objective."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "objective = pickle.load(open(\"my_objective.pkl\", \"rb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you would like to save your objective as a `.csv` file, this can be done as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"my_objective.csv\", \"wb\") as fh:\n",
    "    data = objective.data\n",
    "    wav = data.wavelength\n",
    "    aoi = data.aoi\n",
    "    psi_d = data.psi\n",
    "    delta_d = data.delta\n",
    "    for dat in data.unique_wavelength_data():\n",
    "        wavelength, aois, psi_d, delta_d = dat\n",
    "        wavelength_aois = np.c_[np.ones_like(aois) * wavelength, aois]\n",
    "        psi_m, delta_m = objective.model(wavelength_aois)\n",
    "    save_arr = np.array([wav, aoi, psi_d, delta_d, psi_m, delta_m])\n",
    "\n",
    "    np.savetxt(\n",
    "        fh,\n",
    "        save_arr.T,\n",
    "        delimiter=\",\",\n",
    "        header=\"Wavelength, AOI, Measured Psi, Measured Delta, Modelled Psi, Modelled Delta\",\n",
    "    )"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "a1ddefe8f63dd6da83e48fd4468e36c4710169d87cae4380bd80228d76a38af9"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}