{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysing lipid membrane data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This Jupyter notebook demonstrates the utility of the *refnx* for:\n",
    "\n",
    " - the co-refinement of three contrast variation datasets of a DMPC (1,2-dimyristoyl-sn-glycero-3-phosphocholine) bilayer measured at the solid-liquid interface with a common model\n",
    " - the use of the `LipidLeaflet` component to parameterise the model in terms of physically relevant parameters\n",
    " - the use of Bayesian Markov Chain Monte Carlo (MCMC) to investigate the Posterior distribution of the curvefitting system.\n",
    " - the intrinsic usefulness of Jupyter notebooks to facilitate reproducible research in scientific data analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step in most Python scripts is to import modules and functions that are going to be used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use matplotlib for plotting\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os.path\n",
    "\n",
    "import refnx, scipy\n",
    "\n",
    "# the analysis module contains the curvefitting engine\n",
    "from refnx.analysis import CurveFitter, Objective, Parameter, GlobalObjective, process_chain\n",
    "\n",
    "# the reflect module contains functionality relevant to reflectometry\n",
    "from refnx.reflect import SLD, ReflectModel, Structure, LipidLeaflet\n",
    "\n",
    "# the ReflectDataset object will contain the data\n",
    "from refnx.dataset import ReflectDataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order for the analysis to be exactly reproducible the same package versions must be used. The *conda* packaging manager, and *pip*, can be used to ensure this is the case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('0.1.23.dev0+1b59a5a', '1.7.1')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# version numbers used in this analysis\n",
    "refnx.version.version, scipy.version.version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `ReflectDataset` class is used to represent a dataset. They can be constructed by supplying a filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "pth = os.path.join(os.path.dirname(refnx.__file__), 'analysis', 'test')\n",
    "\n",
    "data_d2o = ReflectDataset(os.path.join(pth, 'c_PLP0016596.dat'))\n",
    "data_d2o.name = \"d2o\"\n",
    "\n",
    "data_hdmix = ReflectDataset(os.path.join(pth, 'c_PLP0016601.dat'))\n",
    "data_hdmix.name = \"hdmix\"\n",
    "\n",
    "data_h2o = ReflectDataset(os.path.join(pth, 'c_PLP0016607.dat'))\n",
    "data_h2o.name = \"h2o\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `SLD` object is used to represent the Scattering Length Density of a material. It has `real` and `imag` attributes because the SLD is a complex number, with the imaginary part accounting for absorption. The units of SLD are $10^{-6} \\mathring{A}^{-2}$\n",
    "\n",
    "The `real` and `imag` attributes are `Parameter` objects. These `Parameter` objects contain the: parameter value, whether it allowed to vary, any interparameter constraints, and bounds applied to the parameter. The bounds applied to a parameter are probability distributions which encode the log-prior probability of the parameter having a certain value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "si = SLD(2.07 + 0j)\n",
    "sio2 = SLD(3.47 + 0j)\n",
    "\n",
    "# the following represent the solvent contrasts used in the experiment\n",
    "d2o = SLD(6.36 + 0j)\n",
    "h2o = SLD(-0.56 + 0j)\n",
    "hdmix = SLD(2.07 + 0j)\n",
    "\n",
    "# We want the `real` attribute parameter to vary in the analysis, and we want to apply\n",
    "# uniform bounds. The `setp` method of a Parameter is a way of changing many aspects of\n",
    "# Parameter behaviour at once.\n",
    "d2o.real.setp(vary=True, bounds=(6.1, 6.36))\n",
    "d2o.real.name='d2o SLD'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `LipidLeaflet` class is used to describe a single lipid leaflet in our interfacial model. A leaflet consists of a head and tail group region. Since we are studying a bilayer then inner and outer `LipidLeaflet`'s are required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameter for the area per molecule each DMPC molecule occupies at the surface. We\n",
    "# use the same area per molecule for the inner and outer leaflets.\n",
    "apm = Parameter(56, 'area per molecule', vary=True, bounds=(52, 65))\n",
    "\n",
    "# the sum of scattering lengths for the lipid head and tail in Angstrom.\n",
    "b_heads = Parameter(6.01e-4, 'b_heads')\n",
    "b_tails = Parameter(-2.92e-4, 'b_tails')\n",
    "\n",
    "# the volume occupied by the head and tail groups in cubic Angstrom.\n",
    "v_heads = Parameter(319, 'v_heads')\n",
    "v_tails = Parameter(782, 'v_tails')\n",
    "\n",
    "# the head and tail group thicknesses.\n",
    "inner_head_thickness = Parameter(9, 'inner_head_thickness', vary=True, bounds=(4, 11))\n",
    "outer_head_thickness = Parameter(9, 'outer_head_thickness', vary=True, bounds=(4, 11))\n",
    "tail_thickness = Parameter(14, 'tail_thickness', vary=True, bounds=(10, 17))\n",
    "\n",
    "# finally construct a `LipidLeaflet` object for the inner and outer leaflets.\n",
    "# Note that here the inner and outer leaflets use the same area per molecule,\n",
    "# same tail thickness, etc, but this is not necessary if the inner and outer\n",
    "# leaflets are different.\n",
    "inner_leaflet = LipidLeaflet(apm,\n",
    "                             b_heads, v_heads, inner_head_thickness,\n",
    "                             b_tails, v_tails, tail_thickness,\n",
    "                             3, 3)\n",
    "\n",
    "# we reverse the monolayer for the outer leaflet because the tail groups face upwards\n",
    "outer_leaflet = LipidLeaflet(apm,\n",
    "                             b_heads, v_heads, outer_head_thickness,\n",
    "                             b_tails, v_tails, tail_thickness,\n",
    "                             3, 0, reverse_monolayer=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Slab` Component represents a layer of uniform scattering length density of a given thickness in our interfacial model. Here we make `Slabs` from `SLD` objects, but other approaches are possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Slab constructed from SLD object.\n",
    "sio2_slab = sio2(15, 3)\n",
    "sio2_slab.thick.setp(vary=True, bounds=(2, 30))\n",
    "sio2_slab.thick.name = 'sio2 thickness'\n",
    "sio2_slab.rough.setp(vary=True, bounds=(0, 7))\n",
    "sio2_slab.rough.name = name='sio2 roughness'\n",
    "sio2_slab.vfsolv.setp(0.1, vary=True, bounds=(0., 0.5))\n",
    "sio2_slab.vfsolv.name = 'sio2 solvation'\n",
    "\n",
    "solv_roughness = Parameter(3, 'bilayer/solvent roughness')\n",
    "solv_roughness.setp(vary=True, bounds=(0, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once all the `Component`s have been constructed we can chain them together to compose a `Structure` object. The `Structure` object represents the interfacial structure of our system. We create different `Structure`s for each contrast. It is important to note that each of the `Structure`s share many components, such as the `LipidLeaflet` objects. This means that parameters used to construct those components are shared between all the `Structure`s, which enables co-refinement of multiple datasets. An alternate way to carry this out would be to apply constraints to underlying parameters, but this way is clearer. Note that the final component for each structure is a `Slab` created from the solvent `SLD`s, we give those slabs a zero thickness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_d2o = si | sio2_slab | inner_leaflet | outer_leaflet | d2o(0, solv_roughness)\n",
    "s_hdmix = si | sio2_slab | inner_leaflet | outer_leaflet | hdmix(0, solv_roughness)\n",
    "s_h2o = si | sio2_slab | inner_leaflet | outer_leaflet | h2o(0, solv_roughness)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Structure`s created in the previous step describe the interfacial structure, these structures are used to create `ReflectModel` objects that know how to apply resolution smearing, scaling factors and background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_d2o = ReflectModel(s_d2o)\n",
    "model_hdmix = ReflectModel(s_hdmix)\n",
    "model_h2o = ReflectModel(s_h2o)\n",
    "\n",
    "model_d2o.scale.setp(vary=True, bounds=(0.9, 1.1))\n",
    "\n",
    "model_d2o.bkg.setp(vary=True, bounds=(-1e-6, 1e-6))\n",
    "model_hdmix.bkg.setp(vary=True, bounds=(-1e-6, 1e-6))\n",
    "model_h2o.bkg.setp(vary=True, bounds=(-1e-6, 1e-6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An `Objective` is constructed from a `ReflectDataset` and `ReflectModel`. Amongst other things `Objective`s can calculate chi-squared, log-likelihood probability, log-prior probability, etc. We then combine all the individual `Objective`s into a `GlobalObjective`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "objective_d2o = Objective(model_d2o, data_d2o)\n",
    "objective_hdmix = Objective(model_hdmix, data_hdmix)\n",
    "objective_h2o = Objective(model_h2o, data_h2o)\n",
    "\n",
    "global_objective = GlobalObjective([objective_d2o, objective_hdmix, objective_h2o])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `CurveFitter` object can perform least squares fitting, or MCMC sampling on the `Objective` used to construct it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitter = CurveFitter(global_objective)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll just do a normal least squares fit here. MCMC sampling is left as an exercise for the reader."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "53it [00:26,  1.67it/s]/Users/andrew/miniconda3/envs/dev3/lib/python3.8/site-packages/scipy/optimize/_numdiff.py:557: RuntimeWarning: invalid value encountered in subtract\n",
      "  df = fun(x) - f0\n",
      "53it [00:26,  1.98it/s]\n"
     ]
    }
   ],
   "source": [
    "fitter.fit('differential_evolution');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAELCAYAAADHksFtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABb50lEQVR4nO3dd3zN1//A8de5N7lZiAhiJBIjdggJQoWWWq1RqrU6KLrQXdXqUK1Snb6/6rC1NVqjiipaVaOqJMTeI8QKEbGy7/n9cRJupox7M67zfDzuw/187ud+PudefN73nPcZQkqJpmmapuXEUNwF0DRN00o2HSg0TdO0XOlAoWmapuVKBwpN0zQtVzpQaJqmablyKO4C2ELFihWln59fcRdD0zStVAkPD78kpayUeb9dBgo/Pz/CwsKKuxiapmmlihAiMrv9dtX0JIToIYSYFhcXV9xF0TRNsxt2FSiklCuklE+7u7sXd1E0TdPshl0FCk3TNM367DJHoWmalp3k5GSioqJISEgo7qIUK2dnZ7y9vXF0dMzT8TpQaJp214iKiqJs2bL4+fkhhCju4hQLKSUxMTFERUVRs2bNPL3HrpqedDJb07TcJCQk4OnpedcGCQAhBJ6envmqVdlVoNDJbE3T7uRuDhLp8vsd2FWg0DRNsza/Mb/hN+a34i5GsdI5ikzS/0E8264WYx5oUMyl0TTNno0bN44yZcpw4cIFVqxYgclkonbt2syePZvy5csXd/Fu0TUKC5NWHbj1/NuNxzNsa5qm2UqnTp3Yu3cvu3fvpm7dukycOLG4i5RBiQ8UQgg3IcRcIcR0IcQgW17r243HM2yv3nfelpfTNK0UCY+Mtcp5JkyYQN26dWnbti2HDh0CoHPnzjg4qAaekJAQoqKiAJV8HzJkCAEBATRr1oz169dbpQz5VSxNT0KIWUB3IFpK2dhif1dgCmAEZkgpJwF9gMVSyhVCiJ+AebYo0/z/TmXZd/l6IsEf/sG1hBTMUuJoNNC4Wjne6NaAIF8PWxRD07QSxDI4DJqxlXnDQgr1fz88PJyFCxcSERFBSkoKzZs3JygoKMMxs2bNol+/fgBMnToVIQR79uzh4MGDdO7cmcOHD+Ps7FzgMhREceUo5gBfAd+n7xBCGIGpQCcgCtguhFgOeAN70g5LtVWBft97jrf+mskDhzbneMx1kyuvdn+FR09d4ednWutgoWl2buvxmFvPk1PMbD0eU6j/95s2baJ37964uroC0LNnzwyvT5gwAQcHBwYNUo0nmzdvZtSoUQDUr18fX19fDh8+TJMmTQpchoIolkAhpdwohPDLtLslcFRKeRxACLEQ6IUKGt5ABLk0lQkhngaeBqhRo0a+y9StcVUiPH34t0bTHI95ZO+fDIz4nbFeIwv9D0bTtJIvpJbnreeODoYM29Y2Z84cVq5cybp160pcF96S1OupOnDaYjsKaAX8D/hKCPEgsCKnN0sppwkhzgE9TCZTUE7H5WRgqxq81bQzPzftnHMBr0bT6MIxADxcTfm9hKZppYzlj8HCNjsBtGvXjsGDB/Pmm2+SkpLCihUreOaZZ1i9ejWTJ09mw4YNt2obAKGhocybN48OHTpw+PBhTp06Rb169QpVhoIoSYEiW1LKG8CQPB67AlgRHBw8PL/XuVOiSgDHK1Sn26F/MAiIvZmU30tomlaKWaMFoXnz5vTr14+mTZtSuXJlWrRoAcDIkSNJTEykU6dOgEpof/vttzz//PM899xzBAQE4ODgwJw5c3Bycip0OfKrJAWKM4CPxbZ32r48E0L0AHrUqVMn3xe3bIsEcHU0gBBUcHXk+fv8qVelLNs2/0D5hOs4GbFpFVTTNPs1duxYxo4dm2Hfa6+9lu2xzs7OzJ49uyiKlauSFCi2A/5CiJqoANEfGFhUF7e88Ts7Gvghm2pm5fsCMP7zEwsebUigzk9o2l3h5KQHi7sIxa5YxlEIIRYA/wL1hBBRQoihUsoUYCSwBjgA/Cyl3Jef8xZmrqe8tEX6+KsKz97dx63Wp1rTNK2kK65eTwNy2L8KWFXQ8xam6Qnu/MvhSKoT/sCSP3Zx4FCSVZJbmqZpJV2JH5mdH7aePXbXTfV1ud+8eqtPtaZpmr2zq0Bh6/Uo6jdSi3x4JlyzeZ9qTdO0ksKuAoWtaxSNm9QC4KEazrrZSdPuFuPc1eMuZleBwubKlweDgdAKBh0kNE3Lt5MnT9K4ceNCH5PZu+++y59//lmYouWqJHWPLbTCJrPvyGCAChUgRucmNE0rOcaPH2/T89tVjaJIlkL19NSBQtPuRqe3WeU0qampDB8+nEaNGtG5c2fi4+MJDw+nadOmNG3alKlTp946ds6cOTz00EN06tQJPz8/vvrqKz7//HOaNWtGSEgIly9fBmDw4MEsXryYuLg46tWrd2v68gEDBjB9+vRCl9muAkVRuF7GndNHo/Q4Ck27G1gGh7k9rRIsjhw5wogRI9i3bx/ly5dnyZIlDBkyhP/7v/9j165dWY7fu3cvS5cuZfv27YwdOxZXV1d27txJ69at+f777zMc6+7uzldffcXgwYNZuHAhsbGxDB+e7xmNstCBIh/CI2P5Nw6uRp1n0IytOlhomr07uen289SkjNsFVLNmTQIDAwEICgri5MmTXLlyhXbt2gHw+OOPZzj+vvvuo2zZslSqVAl3d3d69OgBQEBAACdPnsxy/k6dOhEQEMCIESOYMWNGocsLdhYobN09duvxGOJcylE+/poeR6FpdwO/0NvPjaaM2wVkOamf0Wjk0qVLeT7eYDDc2jYYDKSkpGQ53mw2c+DAAVxdXYmNtc6PWbsKFLbOUYTU8uSqWzkqxF/V4yg07W7g0/L28yeXZ9y2kvLly1O+fHk2b1aLps2bV7hFPL/44gsaNGjA/PnzGTJkCMnJyYUuo10FClsL8vWgebM6uKQk8v79tXUXWU27m9ggSKSbPXs2I0aMIDAwECllgc9z6NAhZsyYwWeffUZoaCjt2rXjww8/LHT5RGEKVVIFBwfLsLAwq583PDKWZc+P44NV/+PeUXP57NUeOlhoWily4MABGjRokL83pQ+2G2ebJu3ikt13IYQIl1IGZz5Wj6PIh63HYzhVRjU3Vb58Xi+Hqml3AzsLEAVhV01PRZGjOF+xGgC+V87r5VA1Tbsr2FWgsIaAuQEEzA3I9rUgXw+GDLgXMwKfK+cZv3Kf7iKraZrd04Ein2KS4Wy5inhfOa+7yGqadlfQgSIHEdER2e4PqeVJVPkq1LhyHqNB6C6ymqbZvRIfKIQQtYQQM4UQi219rV1R22m5/zre0UkMXzs8x2Bx2qMqNa6cByFsXSRN04pZbs3RdwubBgohxCwhRLQQYm+m/V2FEIeEEEeFEGNyO4eU8riUcqgty5luV+QWpn9yku5brpCYmsiKYyuyHLP1eAwnylel8o1YXG5c001PmqblWU5TiL/++uvUr1+fJk2a0Lt3b65cuVL0hcuFrWsUc4CuljuEEEZgKtANaAgMEEI0FEIECCFWZnpUtnH5Mmji354Dvs6E7L+ORLLs6LIstYqQWp4cr+IHQMPYU7rpSdO0QuvUqRN79+5l9+7d1K1bl4kTJxZ3kTKwaaCQUm4ELmfa3RI4mlZTSAIWAr2klHuklN0zPaJtWb7MAisH8l/DMjQ5Fo9LoplkczJhFzIO3Avy9WDEqD4A9HfW/as17W6RU1N0fmU3zXjnzp1xcFDD2kJCQoiKigIgISGBIUOGEBAQQLNmzVi/fr1VypBfxZGjqA6cttiOStuXLSGEpxDiW6CZEOLNXI57WggRJoQIu3jxYoELt7WhG46pkuaHbiCRXE28muWYpOreXDW5Erdtp55FVtPsmGVwyC1vmR/ZTTNuadasWXTr1g2AqVOnIoRgz549LFiwgCeffJKEhIRClyG/SnwyW0oZI6V8VkpZW0qZY31MSjkNeB/YYTIVbCBcRHQEO/3dSHIQhOy/DpClRgGw9cRlDlXyo97FkyQl6y6ymmavLP//Z9fCUBDZTTOebsKECTg4ODBo0CAANm/ezGOPPQZA/fr18fX15fDhw4UuQ34VR6A4A/hYbHun7St2YRfCSHAysNPfldb7bgBw4PKBLL8iPFxNHKrkS4OLJ5HSrEdoa5qdCva6Pe2Ro8Exw3ZBZZ5mPH2q8Dlz5rBy5UrmzZuHKGE9KosjUGwH/IUQNYUQJqA/sNwaJy7sFB7uJvW+fxuVod7pBDyvJGM2m7P8ioi9mcSeqv6US7xBrctn2XtW5yo0zR4FVg689Xx65+kZtq1p9erVTJ48meXLl+Pq6nprf2ho6K1pxw8fPsypU6eoV6+eTcqQG1t3j10A/AvUE0JECSGGSilTgJHAGuAA8LOUcp+VrleohYsOXj4IwJbGZQBove86BmHI8isipJYnu33UrIuBZw+xOFwvjapp9s5WQQJg5MiRXLt2jU6dOhEYGMizzz4LwPPPP4/ZbCYgIIB+/foxZ86cDDWSomLT2WOllANy2L8KWGWD660AVgQHBxdokViJmnL9YA1nYsoauWfvDZq/+kWWfyBBvh4EdW7N1bluND97gF8COuqZZDVNuyM/Pz/27r09rOy1114DYNy4cdke7+zszOzZs4uiaLkq8cns/ChsjaJn7Z4ASINga+OydDiQTFxCbLY9HRp6exBRtS7Nzh7CDFyLL/wqUpqmlTx7ntzDnif3FHcxipVdBYrC5igsaw7XOrbF9coN1i39ONtucbE3k9hZrT71LkbilniTGZtP6OYnTdPskl0FisLWKED9eng35F2meh3FLCB01zUSUhOyJLRDanmyy7s+Rmmm6bnDmKXU3WQ1rRSwx1U98yu/34FdBQprLVz056k/uVLWgT21XAjdrcZTZB54F+TrQddnHiZFGLjn1G5MDgY9nYemlXDOzs7ExMTc1cFCSklMTAzOzs55fo9eCjUb99e4ny1nt7ChaVleWBqN55XsB9o82qER14Na0PfyQUKGhehktqaVcN7e3kRFRVGY2RvsgbOzM97e3nk+3q4CRWF7PaV7pN4jhF8IZ0PgYl5YGk27XddYUUENvMvcA+pq23upOmUy52Ivgw4UmlaiOTo6UrNmzeIuRqljV01P1jSp3SS82zzAuQqO3BdxDSllllpFeGQsIy9WREjJ9PdnMGnVgWIqraZpmu3oQJGLIQFPsbF5eUL2XceUkJIlT7H1eAy7qtXlqsmVe07s5NuNx3Ww0DTN7thVoLBGr6fM/mzmhkuSJGTfNWbvm82iQ4tuvRZSyxNpNLLVtwntTuwEKZm26bjuJqtpml2xq0BhrV5P6cIuhBFWz404NyMdw1Vt4pejv9x6PcjXg6dDa/Fn7ZZ4X42mUfRxpER3k9U0za7YVaCwtmCvYHB0ZEPTstwbcQ2HFJllNtkxDzSg+hP9SBUGuh7agqNR6G6ymqbZFR0ochFYOZDZXWezp50/7jdSaXHwerazybZt04D/fBrT7dA/UMKmB9Y0TSssHSjuILByIPX7v8gNZwOdt1/FjPnWdOTpth6PYU29NtS5HIXv+ZMs2RFVTKXVNE2zPrsKFLZIZgPsv3mcDYFl6Rh+FYcUyeYzmzO8HlLLk3X12wDQ9dA/etpxTdPsil0FCmsns2+dF8maFu54XE+l1f7rbIjakCFPEeTrQfv7Atnu3ZBe+zeQnJyqaxWaptkNuwoUttKzdk/+beLONRcD3f6Ly3bwXZ/m3vzS5H7qXI4i8OwhFoWd1rUKTdPsgg4UeRBYOZDXQ9/hryB3Ouy4imtK1lXvgnw9uNazDzcdnXhkzx8kp0qW6lqFpml2QAeKPPL38Oe3EHfKxpsJjcg+B1K2cgV+qxdKjwMbcUlKIPpaYhGXUtM0zfpKRaAQQjwkhJguhPhJCNG5OMoQdiGM7Q3LctHdga7/xLDi2Iosxzzc3JulTe+nbFI83Q7/w18Ho5n/36liKK2maZr12DxQCCFmCSGihRB7M+3vKoQ4JIQ4KoQYk9s5pJTLpJTDgWeBfrYsLwDj3NXDQrBXMAYHR35r7U7o7mus2/FThuk8QDU/1erTjeMe1Xh8xypSzZJ3f92rcxWappVqRVGjmAN0tdwhhDACU4FuQENggBCioRAiQAixMtOjssVb3057X5ELrBxIrzq9WHGPB46pks5bL/PB1g+yBIs+QT78ENyTZucO0TzqAKlmvfKdpmmlm80DhZRyI3A50+6WwFEp5XEpZRKwEOglpdwjpeye6REtlI+B36WUO7K7jhDiaSFEmBAizGqLkpzelmGzZ+2eHPNxY7+vMw9tikUimfDfhCxdZcs8O4wrzmUYtv0XJHAtPtk65dE0TSsGxZWjqA6cttiOStuXk1HA/UBfIcSz2R0gpZwGvA/sMJlMBS/Zqf/gRAqcToGZnSBszq2XAisH0t6nPb+29aBhZAJ1T8WTKlOzdJV19nBnfmBXuhzZis+V88zYfEI3P2maVmqVimS2lPJ/UsogKeWzUspvczmu8APuTm6CFfGwMUltr3o1Q81iSOMhrG7jSaKDoM/GWAwYskzpEVLLkx+De5IqDDwV9iupZqkH4GmaVmoVV6A4A/hYbHun7SsUq0zh4VoR6jvC8RRIkGA2q+CRJrByICM7vsu64HL02HIFx6QUJm2blKX5aeRj7Vne6F4G7FpD5WuX+Gn7ad0DStO0Uqm4AsV2wF8IUVMIYQL6A8uLqSwZnd8FDR3ADBxKBgH4hWY4JC4pjkXtPSh300zn7XEkm5OzND8NbFWDk8+/gsFsZtSWn3QPKE3TSq2i6B67APgXqCeEiBJCDJVSpgAjgTXAAeBnKeW+wl7LOnM9SahuBHcBe1NAmlWuwqL5KdgrmN0NPThRxcSj6y9jEFmbnwDu6xbCz4Fd6Ld7LT5XzuseUJqmlUpF0etpgJSyqpTSUUrpLaWcmbZ/lZSyrpSytpRygjWuZZWmp6YD1ZoSjR3hWApcN6v9c3veChaBlQOZ2XUWEb1bEXg0nronbzJ5++QMzU+gmqBiXx5NisGBlzbP0z2gNE0rlUpFMjuvrFKj8Gmp/gx0BAnsTruxpyRmyVVcG/Qw8SbBo+suZdv8BGCoXo25Qd3pve9vGp0/qntAaZpW6thVoLDqehQVjeBthB3JICVgBpeMS5w28W/P6jaePLj1Ch43JGevn81Sqwip5cl3bR4lxtWd8X98izlVT0GuaVrpYleBwurrUQQ7QowZTqQCAuIz5hcCKwfS4L2vcUmS9Pj7AosOL2L42uFZekC93i+ET+4dTNDZg/TZ+5de2EjTtFLFrgKF1TVyBFcB/yWB0ZSl9xNAQkN/tjZwY+CfMTikSJJSk7LtAWV6ajDh1erz1vpZuMfF6FqFpmmlhl0FCqsvheogVK3iSAq0/ep2/sLCimMr+KGLJ1Uup9B5exxCiCxrVQD0Dq7B2w++iFvSTcav/Yaftp3S4yo0TSsV7CpQ2GQp1JYmMAKvPJZl7idQy6RualKWE1VMPLn6EoEVmxJYOTDLcUG+HjTvdg9ftH2Mboe38OC+v3l72R4dLDRNK/HsKlDYhJsBgkywKxk+6QgLB2UIGD1r98RodGRu14o0jEzAeeOWLAntdH2aezM7pA/h1eozYc1UfC6f1cFC07QSL0+BQggxSgjhYevCFJbVm57S3ZNWq/g7EQ6uhFldb00WGFg5kN7+vVnZxoOL7g48sfJ8tosagapVvNe7KS/2Gk2qwcjUXz/GMTlJBwtN00q0vNYovIDtQoif0xYcErYsVEHZpOkJoKwBWplgTzKcSwWZmmGywJ61eyKdnfixsydt9l3n0Np5OdYqBraqwfNPduS1B1+m8YVjfLjma8xmqYOFpmklVp4ChZTybcAfmAkMBo4IIT4SQtS2YdlKlrZOqgfUmgQ1rsKcCrvmA7cXNfqpQwWuuhp4avm5bAffpRvYqgYdXh/GlHsG8MjeP3n2vyWYJXouKE3TSqQ85yiklBI4n/ZIATyAxUKIyTYqW/EZl03TlbOADk4QmQp7UgAJO+dnqFWklnFjfqeKdNhxFZ8TuTd/DWxVg0qfTmRFg3aM2TCHBw9s0tORa5pWIuU1R/GiECIcmAz8AwRIKZ8DgoCHbVi+4jMu7vZj6B9Q/0Fo5gjVDapWcdMM5hRVq9j0GYGJSYxuMZofOlXgurMB8eEHvPjXizk2QQEMDPHl5jfT2e7dkC9XfkqHo//p6cg1TStx8lqjqAD0kVJ2kVIuklImA0gpzUB3m5Uun2yWzPZpCf3nw/A/4Y3+ap2K1QkqVxE2G9aNh7k9iYvey9UyDszv5Enn7XGc3rKKoWuG5hos+rWvy++TZrC/ci2+XjaRe46G6XyFpmklSl4DRS0pZaTlDiHEDwBSygNWL1UB2SyZnc6nJTz0IrRzUs1Pey1mgk1NIjghAQfhwNwunlxzMfD8smiSzEk59oJK92Db+jzV/wOOefowY8kHdNu/SQcLTdNKjLwGikaWG0III6rZ6e5zchO0Nak1K1bGQ2zaNOQGBwLrP8xbrd7iRhknvu9SkfvDr9LwRDzLji7LtVYR5OvBq/1bM2DgRCKq1eX/lk/msbAVOlhomlYi5BoohBBvCiGuAU2EEFfTHteAaODXIilhSeMXCkYBD7uo1e8W3YRkCc0Ggk9LHqn3CHO6zWHHwHbEljEyaskFUmVqrr2gQCW33+jfmif7fcC6Oi0Z/+d3fPj7/zF+8Q4dLDRNK1YOub0opZwITBRCTJRSvllEZSrZfFpC9ymw8mV4yAUWxsPKBHCYBQhoOoBAn5aMaj+WOT238vL8KFruv457qzs3hw1sVQOA54xv8fLGHxmxdRH1L57kxRtvsO9sMH2aexPkW+LHPWqaZmdyDRRCiPpSyoPAIiFE88yvSyl32Kxkt8vQAHgRqAisk1J+Y+tr3lHwYPBqqHo8nf9Wjdj2NICYBTvnweCVBPq05Ngbkzi3eggv/HyWJ+qrRfweqfdIrqdODxZvG4zs9arNJ79P4beZI3nn9HM8+t+9dGxYhWfa19YBQ9O0IpNroABeAZ4GPsvmNQl0yO3NQohZqF5R0VLKxhb7uwJTUBNjzJBSTsrpHGnJ8meFEAbge6D4AwXcnkm23Sy4bIb1iVBOQCAqgPi0JFbE89XDXkyYHsX9/8XwkeEj/D38s5000FJ6sHhXCB6oUofPVn7OlJWf0fnIVt6/9jT9DkXTL9hH1zA0TSsSd2p6ejrtz/sKeP45wFeoGzxwKxE+FegERKGmBlmOChoTM73/KSlltBCiJ/Ac8EMBy2EbJzep9bV7Oqu1tZcngEmA+B6qBBJcI5iv23jy2JpLvLjoAn81d2fFsRV3DBSggkW9KmVZsiOKgeUnMfzfxbz4zwLandjB56GP8X1ydxZuO0XHBl66hqFpmk3ldcDd7rTEdr6m7JBSbgQuZ9rdEjgqpTwupUwCFgK9pJR7pJTdMz2i086zXErZDRiUSxmfFkKECSHCLl68mJ9iFpxfKDi4gNEI/dzU0qlL4mF/Aqx6lcDEJN6sP5DP+1ehekwyj629dMceUJaCfD34qHcA7/duyndtHqXT0K8Jr96Q99ZNZ/Wskdx/aAtr953n0e/+1QlvTdNs5k5NT+l6AP2An4UQZuAn4GcpZUHuTtWB0xbbUUCrnA4WQtwL9AGcgFU5HSelnCaEOAf0MJlMRdN116clPLlc1SxcPIFX4IersCgekoFfR/LIpUMc8Pbgr2ZlGb7yIitDPQm7EJanWkU6y9rF/DZN+Wn5r7z291y+++Ujdlatx+T2T/LuMglA7M0kQmp56hqGpmlWk6dAkTbYbjIwWQjhD7wDfIxqLrIpKeXfwN95PHYFsCI4OHi4LcuUgU/L2/mK2BMgv4CFN+GXeLixG0JM9Lx+g/cf9aLtO8d44edzrGu8h4joiHwFiyBfj1s3//n1vehWpyUP7V7Hy5vnsWDhW+yoVo/vDj3MWv8QnEwOzBsWooOFpmlWkdcaBUIIX1Stoh+QCowu4DXPAD4W295p+wpNCNED6FGnTh1rnC7/nMuBk4BBripQrE2Ey2YCu8LA8vH80NmToasu8fP6lQw9s5mZXWbmK1ikS69hfLuhGh0b38vDu/5g2PZf+O6XjzhWoTozWvRmanU37m/mp2sYmqYVWp4ChRDiP8ARWAQ8IqU8Xohrbgf8hRA1UQGiPzCwEOcrOfxCweAADinQ1wX+TIQtSXBJkvhaM6b1vEj3f6/w1g9nGfiea76boCwF+Xow/YlgwiNr8+Wf1egQ2JWuh7bwzLYlTFzzFTEbv2dxwP0satqFLyp6M6xtTcq6OOqgoWlavgk1e/gdDhKinpTyUL5PLsQC4F7UGIgLwHtSyplCiAeAL1FNV7OklBPye+7cBAcHy7Cw3EdC20zYHLWokTkVhAGifGDObpJcDTz3fA08rqby6denmfhYdeq8/cUdx1XkRXhkLINmbCUp2QxIWkXu4bGdv9H5yFYczals9m3K/MBurKvbCulo4pFgHxpVc9e1DU3TMhBChEspg7Pszy1QCCEek1L+KIR4JbvXpZSfW7GMhWbR9DT8yJEjxVeQ09tUgtsvVP057334+ToyTvJ7V3cqnEyh0Yl4+k6sxz0NO9Cz6bAC1yzShUfGsvV4DB6uJsav3EdSshnP65d5dPcfDNi1Bu+r0VxydWd5g/YsCejIvsq1QAicHQ06n6FpGlDwQPGMlPI7IcR72bwspZTjrVlIaynWGkVmp7fB3J5wIx5+i4c9yRz0caLWuSTWNyvL68/74GQ0Mb3LrEIHi3TpQeNafDIzNp9ApqQQemIH/Xb/QYdj23BKTeFAJT+WNO7A8ob3Ua+ZPy/dX1cHC027yxUoUFi8+R4p5T932lfcSkyNIrPT29Ro7fDvOXkghfIrb1ImPhUHM7wwyoeNQe6M9GjOsGYjbvegshLLmsbes3H8sfkA3fb8TZ+96wg8d5hUYWBTzWYsa3I/7v0epl5NL90kpWl3qcIGih1SyuZ32ldSlKgahaWwOfDby+xNNGBcGU+DPfGkGOClF3xoXz2VR24mwAOfqbmkbMQycET88S9+q5bw0N71VLt2iatObqys35ZljTsQUaMRj7SooacJ0bS7SEGbnloDbYCXgC8sXioH9JZSNrVyOa2ixAYKsEh2pxARLmm68hoCWNuyHNXbG2nkBgQ9CU0HWL12kVl6Ejw5MZmQU3vos+8vuh7agltyApHlq/BLo/tY0aQjIZ1a6uS3pt0FChoo2qN6LT0LfGvx0jVghZSyBLXvlOCmp8zSmqJmHFuGy8prDPrzMokOqguYQ0sTtHUCN0eb1y5ABYslO6JYHB5FSooZ56R4uh3eQu+9f9EmcjcGJNurN2Rp4w6satCWxDLleLd7Ix00NM0OFbbpyTfzUqglWYmuUViI2DufEVsmMu/dwzgnmdlR141u2+IQJqCNE4S4QJMHoExlm9cwMucyFodHUfHyBXrt/5s+e//CP+Y0iUZH/vAPYVlABzb6NQNHRx7Rs9hqmt0obKD4AzXQ7kratgewUErZxdoFtYbSEigAxv87niO/zWHuRyf4JdSDfaFleHHJBdwPJoGLgNYmaGkCV2cYvNLmzVHpMtQ0klNpdP4oD+/7i577N1Ah/ioXXcuzvGF7ljbuwJGqtXU+Q9PsQGEDxU4pZbM77SspSlOgiIiOYPja4Ty/IJIhqy4y4iVftjRx45Od0dy/+gocSQFnoIUJ+nWEVo9CfIwao1EEQcOyphF7MwlPR9jwfz/QM+LPW11tD1b05ZfGHVjRrBMj+92jm6U0rZQqbKAIRyWvT6Vt+wK/lLReT6UmR5FJRHQE07d/xYvP/ojn1RT6fOjP1XKOzD53kcDj12BzIhxIUROuNHWE1s7g5apmri2iGoal9NrGH5sP0HnvBvrs/YugswdJMjjwR73WLGjahfDagbzTI0AHDU0rRQobKLoC04ANgABCgaellGusXVBrKE01inQR0RF8OKM/8947zNZGZXjhJT9G1e7NsLg4OLcbdoXBvwmwK1lNyVjHAR5qAY8Mg8TYIqthWLJsnvKNjmRAxBp671mHR8I1TnpU5afArixu3JG4sh46l6FppUChAkXaCSoCIWmbW6WUl6xYPqsqjYECYNGhRZz44GVGzzvD5Me8SXh+OD1r9yQwMUmN7k5JhBspEJYM25PghlRrdQc7QnM3aPNEkXSrzcyyeWrSsp103LeJARGraXl6H0kGB9bWbc28wK7srB3Iuz0a61qGppVQha1RCNTqcrWklOOFEDWAKlLKbdYvauGV1kABEHFhJ659+uO37QiD3qnNUb8yvNXqLR5x9b29QNKBX+HwetiXBNuS4EyqapZq7AitXOHBwVA1sEhzGeksg8aPc37n4bDfeXjvOsonXOdQxRp8H9yTXxrdi9nFVXez1bQSprCB4hvADHSQUjZI6/W0VkrZwvpFLbjSmqPI7McNX9DpodHEOxnoP642ia5OzO46+/ZcUOnzR6UkAmY4Z4btibA3Wa2sV8UAzU3QxAlcHaHZY8VW01iyI4oVW4/Rbe8GBocvp+GF41xxLsNPTbswL6g7UeUqYXLQExNqWklglSk8LHs6CSF26ZHZthERHcH/ffkI0yYdZW0Ld8Y8V4NRQS8wLGDY7YPSZ6h18YTzEbBzPtxIhD1JsCMJzptVLaOhIzRzhFou0Pzx4m2acnFk5TeLGLTtV7oc2oIAfq/bhhkhD9OwV0eql3fRtQtNK0aFDRT/oaby2J4WMCqhahS6e6yNLDq0iPNvv8ioxef45HEfOn26PPfZZdMnHtw5H1KS4FwK7EhWgSMJ8BAQaIKgstDvk2JploLbQaP6tUtEf/w5/cNXUS7xBv/6NuG7Vg+z1T+YecNb62ChacWgsIFiEGoJ1ObAXKAv8LaUcpG1C2oN9hAoACLO78Cj/1P4/LOX32a+RkKrIOKS4gj2Cs45aGSuafw3D/begIgkOJmqjqlthKYmaFwWhq0oli62oIJG+J6TeC+eR7Olc6h6LYYDlfzYPfAZLj/Yi5Z1q+iAoWlFyBq9nuoDHVHdY9dJKQ9Yt4jWYy+BAmDPkU24h3bCJSGV/uNqE+3hiLPRmemdp+dt/QrLmsalRNiVBBGJECfBCWjjBy+8Db2eAiFs/GmyFx4Zy+DvNtFt1188/d9S6sScJrJ8Fb5r25+AMSO5nCR1k5SmFYGCTgpYIbeTSikvW6FsdySEcEON4RgnpVx5p+PtKVDM2DODlSs+Zt4Hxzle1Ykhb9Yk2eRASLUQnmv6XN4XO7Ksaax6A46l1TIOpCXAa3hBl+Yw9Dmo5nV7hb4inDJk6/EYzl2+QfT8xYz6ZwEBF45xqnwVvmrTj1VNOzL3mbY6WGiaDRU0UJwAJKoWQdpz0rallLLWHS46C+gOREspG1vs7wpMQU2YOkNKOekO5xkPXAf2322BIn2Kj3u2XeTL/4tkVYg7bzzjjUEYMRlNea9ZWDq9Df6eCMf/hoRU2J+iahqRqepvtqZRjQBvXA56flzkU4YMmrGV5ORUOh7bzqhN8wi4cIzI8lXZ/tjzRPfqSyt/Lx0wNM0GChoo2kopNwshnKWUCQW4aDvUDf779EAhhDACh4FOQBSwHRiAChoTM53iKaAp4Ima8ejS3RYoQAWLsAthtJq1joAvFzL1ocp8+1BljMLIyGYjM/aGyqv0LrapSarJSZrhcooa+b0rCa5IMAGNHCHQEWqWgW5FEzQyrP+9Yi+hB7fy4j/zaXz+GCc8qvJlhyE88emrBPnlWuHVNC2fChoowqWUQYVZzU4I4QestAgUrVFNSF3Stt8EkFJmDhLp758AuAENgXjUnFPm3K5pb4HiFim5PKAXFX5awdinfVgb6lWwGkU6y+ao1WNU0DAY1YiZkwmqaWpfkmqaqmCAZk7QxAQVnItsnqn0oHE29ibR8xfz2oa51Lt0ilN1mxA24k18e3XRtQtNs5KcAoXDHd6XLISYBngLIf6X+UUp5QsFKEt14LTFdhTQKqeDpZRjAYQQg1E1imyDhBDiaeBpgBo1ahSgWKWAEFSYu4hrZ9oxfmYYDRreV7jz+bS8fbP3ang7LwG3A8ivo2HPDVXLWBcPf8VD7Rtw/X3o1A7q3mfTgBHk60GQr4dqktrZmr9rB9N331+8uOEH+rw4gD++bs3eqV/QuGOO/4Q0TSukO9UoKgL3Ax8D72Z+XUo5944XyFqj6At0lVIOS9t+HGglpRxZkA+Q6Vp2MTL7TnYf3YRTh854X0zi2TH++HceqOaEKmjNIjeWtY75r0F4WhL8qlmtlxHoAi+8CdVdiqxJ6uyVeH755zCDty/nua2LcEtJImbAE6zo/TRNg+vpGoamFVBhx1E0lVLuKuCF/ShE01M+r3VXBIoZe2awcN3nzPnoKGVumhn8Vk3O1ChfuGaovEgPGpdPweJZsCMBDqaopipvIwS7wSuTges2DRq3Et4pZionXmNO1BpqLfmBeEdnptz7BHXeeY3LSWbdpVbT8qmwgaIu8A3gJaVsLIRoAvSUUn6Yh/f6kTFQOKCS2R2BM6hk9kAp5b58fJ5c2W2OIk16T6iK564x96NjAAwbU5uqzdvnr8tsQVkmwm9K2JUA4UlwyawS4AGO0KosjPldHW+DrrbptYuQWp5sPR7D0gXreO+P72h3cif7vWrxdufn2e/bUM8hpWn5UNhAsQF4HfjOYq6nvZZdXnN43wLgXqAicAF4T0o5UwjxAPAlqqfTLCnlhPx9nByvd1fUKEAFixXHVrBrwwKmfXSYFKNg6Bs1OVu9LG2rt8XTxdN2zVGQNRGekghRqRCWAPuSIQWoVwUa3oQGBnBxgq6TbNJryrJL7QOH/uGtP6dT9XoMPzXpzPX3P2DoQ8Uz8lzTSpvCBortUsoWmSYFjJBSBlq/qIVn7zUKSxHRESz7dQKjXluGWcAzr/lxxMcZAJPBxMwuM4umhmEZNK4nwB4z7HGE05dVx+ZAJ2jhqHpPObhYvdeUZZfaz5Zs59kN8xkc9iuibFnOjH6X31o+QCv/yrp2oWm5KGivp3SXhBC1SRtwl5aQPmfF8lmFRY2iuItSZAIrB0KvsTx7+QBfTTrEnInHGfGyLxH+biSbkwm7EGb7QJFd76kRoSAlvN8N/rsO2xJha6KaZ6qVGY5vUMdbqVkqvXcUQL0qZdnarSmHU16lxjuj8R37CiHVvuOdHi8z/s1HdbDQtHzKa42iFmop1DZALHACGCSljLRt8QrmbqpRpIuIjmDj5u/p9fzXVI5NYswzPmwMrkBv/940qNDgzpMJ2kp6beO6AT55B7bfhOtSTRkScBOaGMHVyWbjMqb+dYSjX3zHu39OwzU5nh1DXsQ05g22norTyW5Ny6TQkwKmncQNMAA3gf5SynnWK2Lh3U05ipzs3fcXHo88QdWDZ/jfw17MfLDircn+LCcTTB/tXaTB4/Q2OLIe9ifBtzNhX6RKfgc5wcuvQkgvqye+0/MX5eJieP/P7+h2YDN7vGrz6oMvc6paLZ3s1jQLBR2ZXQ4YgRok9yvwZ9r2q8BuKWUv2xS3cO7GGkUG8fEc69Oe2qu383tLd94fUo0bLsZbU34EewXz+O+PA+RvJlprOr0NJnSDLdfV6G9hgPoGaO0EvmWsWsOw7CEVPWceLSa/Q9nE63za/knKvzWakDqVbr2ug4Z2NytojuIHVFPTv8BwYCxq2rjeUsoIaxeysO7GHEW2XFy4Nuc7vhzZjVFLLtAwMp7Rz/tytGZZzl4/y/Lry28dWmR5jMx8WsLY31UNwlQXpnwOSzfAvhvglwhO06CfhMjNha5hWOYvwgcPomecJ++vmMLYv2ZyPm4//ds8yyk3T70kq6bl4E41ij1SyoC050ZUArtGQSYILEp3fY0iTUR0BKdXzafjmJk4XbrCtO6eTOtRCYPJmSRzElCMNYrMTm+Dad1h+3X4N1GN/PZygDYmaFoGnrLeAkvhkbFsPXaJB8N+p/p7bxIvBW90HcXaBm15pXM9Rtx3l//Q0O5aBW16yjAZYGEmByxKOlBkEhPDkcHd8V+5lSPeTnz8mDe+PZ+kapmqecpRFFk+Iz3xXS0EZn0F03+Gi2bwMMCwPvDBjxC9y6p5jD0bd2AeMICmZw8zL7g7DX6cRvN6Va3wYTSt9ClooEgFbqRvAi6oRHb6ehTlbFDWQtOBIquI6AimT+zF2B/OUS0mmbiu9+L+xTdEVEjINQhEREcUTz7j9DaY3QMO3IBNiXAmBapWhubxEGgEZ+v1lNpx5AKGsW8RuGgWNxs1YemYz2kQ2lw3QWl3nZwChSG3N0kpjVLKcmmPslJKB4vnJS5ICCF6CCGmxcXFFXdRSpzAyoEMf/NX1vz+JWfHjsL9n3Bkw4ZEP9ieP5ZMYvja4Sw6tIgZe2YQER1x631hF24H3PR8RpHwaQlDVsDz78OWzbBmDVRyg9+uwf/i4N+rcHi9VS7V3N+LwJ9ncnTGfJKOHafX0B7Mev1zwiNjrXJ+TSvt8tU9trTQNYo8uHiRnWMHU/vH1ZSLN7O3pgtLQz1Y3cqd5HJujG4xmrikONxN7kzePplkczKOBsfizWec3gbjusKf1+BUilqydeRgaOUB/oWf7nzq+qMsWLSJr5ZNJPDcEbY/PoIWc/4Hhlx/T2ma3bDKOIrSQgeKvImIjuDZpQPptfkKD2+IpW5UAslGwbYGbqxrXo7NTcoQ61XuVtAolgF7mZ3eBic2QpQzfD4NwvdAeQH3l4VP14BvSIFPnT7mQiQkMGHt1/TZ/Sd0787Oj75iy6Vk3X1Ws3s6UGjZSk9UuzuWY8mCN+my7Sodd16lxgXVK+qspyM3Wgfh32UQNGvG7qqw7cYBgr3Uv6Wc8htFkgDf+ClMex/W3YTzZqhTDf43AxpVKHC32ltjLmpWIOi3BciXXuJ4OS+GPfwu57x8dPdZza7pQKHdkWXQWPTreJrtvULw4Zvce8KA48UYAMxCBY/TVZw5XsWRU15OXKrsRpd7n+FcRRNN67QDKJoEePp058mJcMAM/5WFU2fAzxE6OUMN10InvJd+OZ/2bz2HQUqe6fsO7Z/qrbvPanbrrggUegoP68lQI6jUFM6cYe0vH3N4/c/UOpeI37lEfC8k4ZaQcWXa684G4qt4crhMPNEeDlz0MNGwSSfatuwHPj7g7U1E6mnCosOtU9tI71LrFwpegfBGf5i+DG5Itb73u6Ph4Q8KfPrwyFjenLyU7xa8S7WrFznz5TfUGjm0cGXWtBLqrggU6XSNwjbSF0xKNidjFEaQknJXE6l+MZkqlxKoeimJqrGpNEuujPnUSSpdTsIzLgVjpn9iCY6C8xUcuVDRiTpNOuBZvznUqgW1a0Pt2kTIMwUPJKe3wXfdYcM12JoADk7w3BPQ1Q8a3F+g2kV4ZCw7dx6l/4RRlAn7j6g33uPXzoMIqV1RN0NpdkUHCs0qLGsaoHIU2fWMSn+tvKEM3/wxDq/LyXjHCToY6nLu4HaqXk6iWkwKdeIccb2UsTvzdWcDJ6s6caqaC81CB5BQtybbK8ZTN7gbgVWD7pz/SK9lGGvDxKmw8k8oK6BbOfhkDdRoVbAPn5DA5UcGUmHlL8wO6sHkrs/w4/A2OlhodqOw61FoGqDGY1jenNOf+3v4Z7l5B1YOZMaeGUR7OBLt4ch+YaSsf2uWH7t4K6iMbjGaSRvHUf1SMrViJF1T/YnZ+x81zyXS7MA1qv4zFYCaQLzpPWLq1+KoRwxnfV341N+d14bMJbB6i4yFtFwfY/QxqPwf/H4Dfo6DY4/Dj8uhfv38f3hnZxa8NAnnC2aGbv+VMknx/HdvHR0oNLunA4VmFZkDSLpgr2Ccjc63AkOP2j3oUbvHraASdiGMJJOBE9WcOFXdiKd/a5YH3Q4kfap2JmLTT9SNSqD+6USCL9ygU9hV+m6IBc6S8uE93AhowIl6lXDr0JWa3Z+AypVv1zo8qhJY0xWGG2FnKmw4D02awCuvwJBuEB2Wr95RIXUqMajz01w3ufLiPwu4POVNwv2+Zevpq7r7rGa3SnzTkxDiXuADYB+wUEr5953eo5ueSpbcmoos8x6Zm63Sm7csXx/dYjSTt31MpegbBJ5I4skbDbm55W/qn7yBS5L6t5zgX5MV1a+wvb4bOwM8+aTdaAJjz6mA4FwT3ngDZs8GdwM84AIN8zeteXoX2l5/zsf7o/f4y78Vz/d6A5yddfdZrVQrlhyFEGIW0B2IllI2ttjfFZgCGIEZUspJuZyjPTAGuAB8KKU8eqfr6kBRutwp55D5dcvtsAthTNkxBYcUMwGRSYy81pQqOw5Tfvs+ysWbMQuIaeBLpYcGQZcu0Lo1ODrC16Pg/a8h2gyNHWHCW9BzXL7L/vdL73PvlHFs9GvGM33fYeQDAbr7rFZqFVegaAdcB75PDxRp05UfBjoBUcB2YAAqaEzMdIqngEtSSrMQwgv4XEo56E7X1YHi7pFTjeSZ1cPwP3aVtvtu8NjpypQJ3wOpqVC+PHTvDqGN4eRnsPUabEyAcu4w5f+gfd18DdYLj4xl8agPmLjiC/7yb4X7quUE1als40+tabZRbL2ehBB+wEqLQNEaGCel7JK2/SaAlDJzkMh8HhMwX0rZN4fXnwaeBqhRo0ZQZGSJXM5bs4HsaiRZ9sXFwbp1sHw5rFgBly+rGWiD60DovbBuO2zbBv4m6O4CFZzz3BwVHhnL9S/+R/sp46B/f/jxRzAabfmRNc0mSlKvp+rAaYvtKCDH/opCiD5AF6A88FVOx0kppwkhzgE9TCZTkHWKqpUG2SXSs+xzdyeibS3C/NsS/NFIAg9ehWXLSFryM6bNU0l1dcFY3xuORMHUJOieAu035SlQBPl6wJfvQVVnGDMGypSBadNurVWuaaVdie/1JKVcCizN47ErgBXBwcHDbVsqrbTJbl0NGj/FE83+otkRV3ptvUbPsMs4pAISWHoDYpfCkqfAI4/J6TfegGvXYMIELkhHFg98WQ/K0+xCcQSKM4CPxbZ32r5C02tmaznJaV0NaRDsqOfGrvrluPLJMzy14xpX5s6g3L9HMazfAlWrwuuvwzvvEHFl/50nOvzgAy6cuYjXzG+IO5bAoLaP6J5QWqlXHBPtbwf8hRA10/IO/YHlxVAO7S6SPp7DKIw4GhwJ9grOsq+5T2sienUmdLgTHb6ox8/3V8Kcmgoffoi5nBu7nurCT+s+Z/ja4RkWd8pACBY/9ior64cy5u85tDvwL1uPxxTpZ9U0a7N1r6cFwL1ARVT31veklDOFEA8AX6J6Os2SUk6w5nV1ryctO3lJes/YM4MpO6YAYBRGXvTqxpAXPoHDqUggxQirWnuQ/OpL9H3o3WyvEx4Zy9Bv/ub7uW9QO+Y0kcvW0PCBdkX0KTWt4O6KuZ707LFaYWXpbuvVgcB/vuXEEUnlJddxSFULxptSJNwfCg81gZ5PZEl6h0fGsnv7AQaN6ovJ0UH1qKpSpXg+lKbl0V0RKNLpGoVWGBlqGYlJas2L1CT23zDh9Qt4Hr0ALZrCvt1wU0INR5j8JfR7PuvJdu6Etm254V+fHybNpUWD6jpfoZVYd0Wg0DUKzSYs17yoGAAvvgjTp4OPgbP1nSm3PZEyV1KhUyf46CMIzvj/7Nh3P1D72Sf4tUF73ugzmnnDW+tgoZVIOQUKu1o1Xkq5Qkr5tLu7e3EXRbMnPi0h9FX1p4uLGiPxxXuknpc4bE/iuRE+fDGgOinh26FFC+jbFw4cuPX21XVb80no4/Q6sIHeO1br5LZW6thVoNC0IvPSOH79bjgpRsG0TyI54+nAvJXjYdw4WLsWGjeGYcMgOpqQWp7MatePzX6BvPvHNO41Xyru0mtavthVoBBC9BBCTIuLi7vzwZpWSLUefI4h7zfkcA0XPv0qks6/HoR334Xjx1Xz1PffQ716BP22gB+fasWRT6biUK4MjV55BhISirv4mpZndpWjSKeT2VpRiYiOYOepf+k1aSUVlqyCgQOh5goi3Jw46v0IXaaFU3bzNmjeHN4bAQe3wBszYdQo/Fy7ALDkOb1KnlYy3BXJ7HQ6UGhFTkqYOBHGjuWGjwO9X6nFhfKOmAyOLEroj9+7n8GFi9DMEYwOEBbP0IffYV2dVjg7GvToba1EuCuS2brpSSs2QsBbb8EHT2K6kMr0iSepcimZZJnKn609Ye6L0NoEu5JhXwIJFd35ZNUUvK5dIjnFrBPcWolmV4FC93rSit2Tz3PiqXJ4XE/h+4+O438uRa3U17ATdHaGZ9zAw4jzpTjKJVznk9+n4GgUhNTyBNRAvanrjxIeGVvMH0TTbrOrQKFpxc6nJXXfWsvZL56krHBmwSfnCDyVorrWDv0DBrwPmzbAO+9gRNLuxE7+TNhCkK8H4ZGxPPzNFj5Zc4hBM7bqYFGMAuYGEDA3IOc5ve4yOlBomrX5tKT+07Nw/W8nDuXKw333wd9/3x6PUasNjB+P2LABTCa8P5sAr7zCtoNnb51CN0fdln7TLqpzWQaHXCeAvIvoQKFptlKnDmzeDDVqQNeuamU9S6Ght/d98QVPvtSPWjFRADg6GG41R2mKNW/YuZ0rpynp72Z2FSh0MlsrcapXh40boUkT6N2b0c/5ZPxF27kzDB4MRiOuF87yx0+v85X7Od0LKk1ef91bs6YQ7HW700/6lPR3O7sKFDqZrZVInp5qve62bfloWhQdwq9mvElVXAyOqVC/PsZateg+9mmCfpquutyWYtZoMsrvr3tr1BQsF6Wa3nl6zotU3UXsKlBoWolVtiy7Z01gn58Ln3xzmhmfPXr7puZmgI7O8O+/vB10Uc0V9cYb8NhjEB9frMW2hsI0GeXl170tagp7ntzDnif3WD1IWDPfUpR0oNC0IrLtxgGef9WP41Wd+OTLo5xeNV/NTAvQ3BFzNSOvLb6gJh388EOYPx9ql+H+L+qXupuLtRLCefl1b881hXwFlnHu6mEDxbFmtqbdlYK9gpniZuSZ1/2YO/EkD4z6Gj5PUi8aBId6l6Xe11e48O4reP1vFgQEwCO9WDjuGC+NqkFEdESx3dz8xvwG5H26kZ2RW2h0/CZesSl43LhCUviHYKoNly5BTAxcvQpGI5svbCXZKLiv5v3g4ACOjuDqCpUr33q03H+dy+UcGPrBPmJdTnF8cs8M18pvTaEkKM6/y4Io8YFCCGEAPgDKAWFSyrnFXCRNK5D0G8Plcg7cWLUMY+9nYfRsGJBKhI8Lw4IqM76lmfbT57JnxJME9OzJwefK47rgJjM/PsEbCY8wePSiIr/BWI7nGDRja9ZE+/XrsGsXhIfDjh2wYweD9+9nSGqqxVmiwMlJ5Ws8PcHdHRIScL+eitEsIeUIpKSox7VrKqCkvX/mrXM8RqLRgYQlfjjX84fataFOHQJr16bm2QTOVDQxvVfJrSlkrmWVlloN2DhQCCFmAd2BaCllY4v9XYEpqDWzZ0gpJ+Vyml6ANxADRNmwuJpmcxl+0f75p+oi+8N5Do8qQ2JVwYzulXngv6skf/kZvOXCZt8yfP+OF998FsknXx7jb+9p8NLX+bpmem3g5KQHC1Rmy/Ec6eM7giqa4Jdf4Mcf4Y8/wGxWB1SuDEFBiF69eOnKXM5WNPFe9//RqH57VVMQIsO5B6Y1q2T5pW82Q2wsREdDdDTH9h1nzrJtVL96kZpx5wmNjMJ182YVVIDlQKoQJE96FJo3VbWxxo3Vn7VqgdFYoM8O+a9N5SS7JrK8Bop81UBOb8uyNG9h2TpHMQfoarlDCGEEpgLdgIbAACFEQyFEgBBiZaZHZaAesEVK+QrwnI3Lq2lFp3ZtFSzM8NA3F/GMS+GIjzPrg8rTZMHfsO8PghMSiHcz8uzrfhyo6cL9r02DhQsLdDnLmoHfmN9u3QDvJH08h8GcSvtTEfT/6m3w8oLHHyfq351806IPLF8OZ87A+fOwahV88AHrgt054OdCo6AHwM0tS5CwlCWHYTComkeDBtC+PasbhPJD8+5MuncIzz/0JrO/WgpxcRAdzcElq3mp+6tMDXmUTY4VSdgertYFefhhqFsXypZVC0oFmaC7C4SFQWJivr+zwo6Wz2+323zledJzXaCW7rXctgKb1iiklBuFEH6ZdrcEjkopjwMIIRYCvaSUE1G1jwyEEFFAWkMuqZlftzjuaeBpgBo1ahS+8JpWFBo2hI3bMYWGMuX/TvHUGzWpNnEqDp0HwfozBLok8e35aLa5OMOP/0OM/h4GDVI3uiefvOPp79RsFB4Zm2E7u9pHULUyDNu2lOHbl+F1/bJqNho4EB57jNDf4pDCQMsmbQiqlvHXdurNGqTcrJXlGukiJlaCalWAOzfFWA4+vDUYUQioVIl1HrVZ1igFAKOAVzrXY0TLqrB/P+zdC3v2wO7d8E8y7EiG31qoXEjjxhAUpB7BwWqsi8mU4brZ1qYKWKvIbzI9XzWQk5tuP09NUttWrFUUR46iOnDaYjsKaJXL8UuB/xNChAIbczpISjlNCHEO6GEymYKsUlJNKwpBQfDDDzTt25dxs89Q758B0O1HmLUIhks2ebgQ5uzMc62Gwu8D4KGH1CC9hAQCnL8Cck7SZnejs5Rd8DjpPBDGAePiYP16GDmSt/fvV2uCP/00dO8Ozs6ER8YixZZszxMeGcvNyOdzvAZAmLPz7bKl3Qgf+vyMKkOmZjLL92Y+V7ZBxM1N1SJatLh9kvfKwRUJobNVPiU8HJYuhRkz1OtOTmrdkJAQ9WjVipCaFbKe2wry0oyUrxqIX+jt50ZTxm0rKPHJbCnlTWBoHo9dAawIDg4ebttSaZqVPfwwfPABPd55ByZNgrffhnvugXAnZg606PLo6qqaefr2hWefZdDAKszrXDHH02Z3E73Tr+SpKT0Jvb6LJgMGqGaumjXVNXv0yHDu3M6Tl1/iwRar/N2+EZ7J+AHSu3uOi8sxx5JbEMlACPAQ6rvr21ftkxIiI2H7dvjvP9i6Fb75Br74Qp27ShWmlfVjZ7V69BrZj/peLtmf2wbyVQOxrD08udzqOYriCBRnAB+LbW+y/OsoGCFED6BHnTp1rHE6TStaY8eq5pK33lKJ4vvugy1/s6dFJH4pC24f5+ysfgkHujFm/nkcUyQ8manZKO0GGzTu9nQ22d1ELX8lh0fG4pCawpVtjtT+5whmDmN47z01+M8l6w1y6sEXQQwDacryazvbX/mZBCYm3Xp++0Z45lZZ8tPEU9BEPUKAn596PPKI2pecrJqq0gJHnZXr6HxkK2yYq5qmgoNVEL/nHvjnMXA1qNpXHhS0e26+ekdZOUhA8Qy42w74CyFqCiFMQH9UpwVNu7sJATNnQsuWalT2gAFwXfLOtkFAxnwDJhP0dWFVK3de/fkCzL1zr/H0G29Ov8DD9p1m1uL3Gbt+Nlt9Apg38zeVFM4mSAAYXU/hWmNGlvPkdo3Mhl6Jo2lCIoGVA62aOM7VnRK9jo6qOfD55+H776l1+QxcvAi//qrWQgeYMkU1AX5yHb66DsOGqR5gZ6zym7fEsWmgEEIsAP4F6gkhooQQQ6WUKcBIYA1wAPhZSrnPGtfTcz1ppZ6LCyxbBh4e8MEH3GgSyLNbl+CYmpz15mkQjB1ena0N3WDYMFpH7sr2lCcnPZjxF/c4d5WHwOKGfvUqA8c9wz2Ru3ij2yhG9B1Dw9Dmdyyu0fVUxvNkI9vX0m7WI2PjmH4+Gk5vu2M+pTD6vPnF7Y2C9AqqWBF69oTJk+Gff1SPq02boKMTeBpgyRJ4/HHw9lY9rZ55BhYsgHPnrPYZipNNA4WUcoCUsqqU0lFK6S2lnJm2f5WUsq6UsraUcoK1rqdnj9XsQtWqKicQE0Ni7FWqX7vIQ/vWZ3vzTHEw8MqIGsTX9uW7Xz6izqVTef4lXrbBGDU9xOXLcP/9lN0Zxgs9XqdF82OsM71GkOFIns6TU7NPlgCV4UXVS8cBcJQSTm7KU3NVQYUYDtzeSO8VVBjOztC2LbR1ggGuaoDgjh3w2WdQr57K7QwcCNWqqS6+zz2nBiWWUnY115OuUWh2o1kz+OEHKpw+TpyTG4MiVme5eUY4mTBKyTVXA32fcSHB0cicReN48ctVhJv983SZCldTVC5k1y5YupRzDb3o67CZ6oYYm/THvyWtV04KkCwE+IXmPSldAFvNDW5vZOoVlJ8xJTkyGtXf2SuvqDVGYmJUgnzyZNUZ4Mcf4ZM2+Z6LKX1ywuJW4ns95YdOZmt2pU8fePdd3MePJ/DcYZbeU5aG6TfP09sIc3ZGAgjB6QqCZ57qybxvFvH1wvf59/EgglxzrxFUik1mxuSTcMUAK1dCp06EbH369gE26I9/S9o5HQCHp9ZmuUaWIGE52vjqVTh+/Pbjp7fgihnKBKg8T/rDYLj1fOzpOC4ZjJQVN3Gq3QhOz4Qaf4CvLy1Pn+FMucpq+hAHK90SHRxU0js4GF5/XSXIP7DxQlR5TKgXhF0FCt09VrM7777Lf7OW0DJqHw1//Ba6tgVg8rcz6FwmAZMsRzLgIIzsqtCaUT1rMm3pBGr8Gg39clnP4qqZOd+cwPNqCqxZD+3aATD62WEw8yd1jA3642crp0B06j84mwpHUmBJOzDWhqgL6te6JRcB5QXU91aBwWxW3V7THnE3ErlhSuKgrIEBM00vJ+D2yy8qQQ38nH6eacPUQlM1aoCvr5r6IzgYWrWCKlUK9xkdHcGQ88j0ks6uAoWm2R2jkZG9xrDp26dwXrAAvvwSKlZkq7kBoxN/ujVqu/X9H/PQAXfW1fHlzPuTqPHuaFhjguHZzPsjJSyPp2KcZNhoP+anBQmgQP3xrTkTqlNKkkrKP7sCliyASzdAAB4GqJ2iurDWqnX7UbMmfOmr3jxuRbbn/HH9UT5ZcwgAI6m80qUhI+6ro9b6OHWKx95dRPWr0XzcorwaU3HqFPz7L/z0062JCalRQwWM9Efz5mpMSxq/Mb8VvItuKWBXgUI3PWn26GIZDz4LfYyxf89WXTI3bWKHrAvAJte0UduNBwKqnb3G4PawygRbk+DZ++HbPzPe8Hcmw7FUPn+8Kntqu2a9YLpcgoRVZ0KNjobffoPlyzm0di3cvAllykBoC3D6F/wdoJwbPPlDgZrB0vM6RlJxJOV2nsfFBerVY3PNZgA8OjTTpH83b8LOnWo8Rfpj0SL1mtGopvwwxIO3EZ8q5wv22UsJuwoUuulJs0cnJz0I5m5QYanqmjl9OmomHJhZPpvk6MlN0NkJ4syw6hr8NB1eS7vBnj4NaxPAz8jP91XI+t48KsxMqH5jfgMpOXk2BTYlwvgqqpbj7a2mJunZE+69V02pkZ78LcRo4/Sb/ysOiwgxHCDIt8+t13KdC8vV9fbAunTnz8O2bbcDx+ZkCE9mE8Ng/WR44AHo1k015VlMUVLa2VWvJ02zWwYDvPyyej5qFA0vHM/5WL9Q1Vbf2wW8HODDn1XSV0oYPhzMQE8XZCHazPM7E6qloKj9LFzwJvx4UyWh33tP/XI/dQqmToUuXVSQsGSFhPoIh+VZuvzme+xGlSoqkE2YALM/gjfKwnNuJHVxg2qeavqPLl3UzLc9eqjtkycLXfbipgOFppUWTz2lAoDJxNRfJ0JiDsnq9Juqo4DFP4PBCL17q5vWmjVwvzN4GGiakLeptrNToGVFw8OhWzeWzButRjt3c4aRZVSgCAzMdRpyWynU2I2Tm1SZKxsxtHKE8Q+rMSm//QZDhsC+fWp0d82aapbgNQlwPAWSku587hLGrpqedI5Cs2s+PtC1K2zfjk/MBVhhgHoy9xts296wwFU1h7z4IoQ0gxZHAZh+PprhVSoXulg5BomkJNXUtWGDqins2AFOTkS6e2GUEtbEwuoE+MxF9QpKXwo1/eHgAHHXoYyAqGG3k9fpj8qVCx1cCjV2w6JHWDIOOPiFquaqBx5QDynh8GH4/Xe1RsdfB1Te6NdK8OCDqvtz164qH1PC2VWg0DkKze4NGwa//84vjTvyyN51DPjNzILuavW2HHvddO4MderAkSPgX+HWzdVRygwzuBZaUpKa1mLFCnVjPJLNOA5PT6JFeaLcK+PtcUX1aGo9So0zSElRf6Y/UlIg4iRck2qcx4ULGc/l6qom80u+qabRqLdAzZNVq1aBAki+B/hZNIcNSnqLpZmbx4RQo7Tr1YOXXoK3ysGJFCjzqJqmZcEClcfo3FkFjR49oELB80a2ZFeBQtPsXvfuJFesRLmE66z1b8UbS7YTVbEZ5LaG0ezZ6qYdGAjz1sNjzlDTgXBnJ35zc+Pe/HZvTU8w16xB+Wsp8MMPKjisWaMGw5lM6ubo6Khuln36wGuvqWVJTSYeSRsF3TttvinGTc7lWmkjpsedhxs3VHv/iRMZH/8eVDfggWnnq1BBBYz0h+WaFDYQbvZnh6x75xlvTQLqOcK46aoZ8J9/1CzAv/yipmwxGlUSv3dv1butenWbljs/dKDQtNLEZGJPx150/HkW9w3/Dq9rMXw6cxcMDVODwzI7d05NK9G+vZr9tE0bWHyAfSPKM9yvMlKI/HdvvW6GiGTmzj9O06M3QT6h5qfq109NB7J0KSxerG56P/ygejNZg5sbNGqkHpbGuYNZwsObVI+k9MeHH95ay/uknx+43YTqRjVdSUCA6iBQSOFmfwYlvQXkvEBTthwc1N9J+/ZqbEz6IkpLl8LIkeoREqKCRrduajW+YsjhpLOrZLaeFFC7G5ieHo6DNNPz4EaG9X2T2LKOatW5yMisB7/0EiQkqC617u7q12uqxGvBdUzJKhme3r31jqKi4IUXYMp1WJeIc5KZ73pWUnMaRUWpBG76WhoTJqj1wK0VJO7EIKBpU9Wra/p0FQyuXoWNG+HTT1XNIioV1iaqmpWXFzz6KHz7LTUvn1H5hALYam5AUtrv7QLPeCuECvIffQQHD6o1SSZMUM1vb7yhxmt4e6vvd8ECNQFhEbOrGoXOUWh3g8YdWrKhZnNe3LqAdR1O8fyr1fl18kWVQP3nHyhfXh14NAV+/hnGjwf/tEkC69aF3i5UXBjPO3PP8vaw6jga79C99cQJtere7NnqhtrYEdqa6BdcE4DnmzWDiRNV7yUfH9i8Wf0aLm5ubhAaqh4A41arsSXNpsJff8G6dbBoEeuBc2U84dwD0LEjdOigPkcehBgOYCKFBIzWm/G2QQP1eOstFYDXrlXNer/+CnPmqMASFKRyG126QOvWqpnPhuyqRqFpd4sxXUdhMqTw6dxtnPIyqSaLI0fUkqpJSZAsOb02mRNVTDB6dMY313OE9iZ6/XOFV366wIx7v86+2enwYfUr1t9f3aCGDlXX6OUCniqB7nU5Ge6/Xy3d+uijEBGR9yAxLs6mE9lle70vrsETT6jPc+oUHD7MW11GEO7dUPVOGjxYTdfhaVRrSvz0kxo5noMgwxEMvjMxVVpt9RlvAVWTeOopVY6LF9Ugv/ffV+NMPv5YNV15ekKvXvD113DsmHWvn0YHCk0rhc6VqwQ9XGh8Ip6RS6PVr+AZM9Qv5eHDYUMCPheTGf9ktayD1wDaO/HzvR4MWX2Jpr2eVU0e6aJTYclN9at24ULVXn78uErA+vndOqxD+FUWv3NUNT3NmQPz5qnmrdJCCPD3Z35gN0b2ekP1qtq1C7o4QUWD+uz9+6tmqoAA1Yy3fDlcuZLhNEbXUzhV/Nv6QSIzo1E1ob3zjqq1xcSoHwgDB6qlW0eMUL3bVq+2+qXtqulJ0+4qDR1Z3N6ZIb9fUs0oTzyhegW99x4I+PWe8oQ1yKGPvhB8MLg6/wSUZcr8SDXJ3QsvwKFDsOwGOMKsrhX5vktF/n7hy9vvM5vhcDJsTmLK6avs93Wm/Nodqkkrn/K7LrbNGQwqHxDiBCHA2zFq7Mdff6nHtGlqCVSDQTX9dOgAZ1JwqWYm3qkYfnO7u6tkd+/eqknw6FHVTNWmjdUvVeIDhRAiFBiEKmtDKaX1vwVNK6X6tkmBC/XVMpy7dsHYsWpwW3Q0sWWM1IuMVzf3HHr4/BVUDjpMUbWQjz9WvXFCHCHUiS8aWUytnZysfmF//DHsiwd3wcRBVVl0rwc78hEkcp1bqaRxcLjdxXbMGEhMVE0/6YHj888hOZl/5+3neHUnWD9YJaWDglRi3WJ2WZtLqx3dykVZma3XzJ4lhIgWQuzNtL+rEOKQEOKoEGJMbueQUm6SUj4LrATuvIK8pt1NTELdwGNiVA5hzhzVpu5lYPCaGBa/d0y1YffsqZbpDAtTXUmlJGTfdWZNPK5+kYIa8OXgALtT4JSaXtsl0Qz/+59q0njiibQ5pJxhVBnmd/Ik2TF/txBbrottc05OarK/ceNUb6rYWBjkyrSelThXwVE1+YwapX7RlyunaidDhqjAvXUrpPUys9mqgTZk6xrFHOAr4Pv0HUIIIzAV6AREAduFEMsBIzAx0/ufklKmZ5IGAkNtXF5NK32aNlVLbr70krpZtW0LHXZxv3tVWhy6yURzVzWNxoq09RqcgLIGpl+6xoXyDqof//Dh6hfwgQMwaBD8tJOvDkbS5NhNuP6iOufXX6ueVe+XL3BRbbkudpFzc4M6Dnxd0wuAPU/shjNn1JiI8HAVlH/7TQVvUKPQKxpgXigEd4W6gWpQXfrD21sF9WIcL5ETmwYKKeVGIYRfpt0tgaNSyuMAQoiFQC8p5USge3bnEULUAOKklNdsWV5NK038EubfnrbjhRdU+/Qff8C338KiNlzwNLGyjYmJT05Xx5w9mzau4AmIMTP+wWosa1ueHcNevH3SBg3Ur9/33iPk04/Z0rgM9331e8aptgvBlutiZ1gutTgIoW723t6qFxKo3EFUFPzwLqyeD9FmuJYKf/wNP63IOn7DyQmqVbsdPDw81Eh3Jyf1p+Xz7PZ5edlNjqI6cNpiOwpodYf3DAVm53aAEOJp4GmAGjVqFKZ8mlb6CKEGup05oybMW5TNMdWqqV48B58BYFHNHOYVMplg4kSC660AIdhjpSCRmVWChGUzztyehVq3wiaEUGMyHn8OkpaqfQ4uqpxVmqn1Lc6cyfiIirpdM7l6VXV3TkxUf6avuJeT9u3h77+t/jFKfDIbQEr5Xh6OmSaEOAf0MJlMQUVQLE0rWUwmFSSspQQ2gWRxctPt56lJarskBYp0OS0x6+OT58F9gAoUloHD8nlios0WSyqOQHEGsPxmvNP2aZqm5Y/FVN8YTRm386jI17ouTCAzGtUSri4u1itPHhTHgLvtgL8QoqYQwgT0B5Zb48RSyhVSyqfdS9OgH03TCi6nX+qaVdm6e+wC4F+gnhAiSggxVEqZAowE1gAHgJ+llPusdD09KaCm3a10kLAZW/d6GpDD/lXAKhtcT08KqGmaZmV2NdeTrlFomqZZn10FCp2j0LQSYpz77ZXwtFLPrgKFrlFomqZZn10FCl2j0DRNsz67ChS6RqFpmmZ9QhZwrdiSTAhxEchmAeEsKgJFvwCtdZTWspfWcoMue3HRZS86vlLKSpl32mWgyCshRJiUMpfFgkuu0lr20lpu0GUvLrrsxc+ump40TdM069OBQtM0TcvV3R4ophV3AQqhtJa9tJYbdNmLiy57MburcxSapmnand3tNQpN0zTtDnSg0DRN03Jll4FCCNFVCHFICHFUCDEmm9edhBA/pb3+n+W63kKIN9P2HxJCdCnSglPwsgsh/IQQ8UKIiLTHtyWw7O2EEDuEEClCiL6ZXntSCHEk7fFk0ZX61vULU/ZUi+/dKmur5Eceyv6KEGK/EGK3EGKdEMLX4rWS/r3nVvaS/r0/K4TYk1a+zUKIhhavFet9Jt+klHb1AIzAMaAWYAJ2AQ0zHfM88G3a8/7AT2nPG6Yd7wTUTDuPsZSU3Q/YW8K/dz+gCfA90NdifwXgeNqfHmnPPUpD2dNeu17Cv/f7ANe0589Z/JspDd97tmUvJd97OYvnPYHVac+L9T5TkIc91ihaAkellMellEnAQqBXpmN6AXPTni8GOgohRNr+hVLKRCnlCeBo2vmKSmHKXtzuWHYp5Ukp5W7AnOm9XYA/pJSXpZSxwB9A16IodJrClL245aXs66WUN9M2t6KWH4bS8b3nVPbilpeyX7XYdAPSew4V930m3+wxUFQHTltsR6Xty/YYqVbciwM88/heWypM2QFqCiF2CiE2CCHyv3hw4RTmuysN33tunIUQYUKIrUKIh6xasjvLb9mHAr8X8L3WVpiyQyn43oUQI4QQx4DJwAv5eW9JYtMV7rQidQ6oIaWMEUIEAcuEEI0y/arRbMNXSnlGCFEL+EsIsUdKeay4C5WZEOIxIBhoX9xlya8cyl7iv3cp5VRgqhBiIPA2UOR5IGuwxxrFGcDHYts7bV+2xwghHAB3ICaP77WlApc9rRobAyClDEe1e9a1eYmzKVea/Hx3peF7z5GU8kzan8eBv4Fm1izcHeSp7EKI+4GxQE8pZWJ+3mtDhSl7qfjeLSwEHirge4tfcSdJrP1A1ZKOo5JE6UmmRpmOGUHGhPDPac8bkTHJdJyiTWYXpuyV0suKSrCdASqUpLJbHDuHrMnsE6iEqkfa89JSdg/AKe15ReAImZKaxV121A30GOCfaX+J/95zKXtp+N79LZ73AMLSnhfrfaZAn7e4C2Cjv8QHgMNp/8DGpu0bj/pFAuAMLEIlkbYBtSzeOzbtfYeAbqWl7MDDwD4gAtgB9CiBZW+Bao+9garB7bN471Npn+koMKS0lB1oA+xJ+4+/BxhaAsv+J3Ah7d9GBLC8FH3v2Za9lHzvUyz+T67HIpAU930mvw89hYemaZqWK3vMUWiapmlWpAOFpmmalisdKDRN07Rc6UChaZqm5UoHCk3TNC1XOlBomqZpudKBQtM0TcuVDhSaVkyEEP+XtsZFizwcW0sIMVMIsbgoyqZplnSg0LRiIIRwAyoDzwDd73S8VNNZD7V5wTQtG3r2WE3LByGENzAVtfiMEVgFvCotJqvLdPy3wA9Syn8s90spbwghqqIms6thcXwAMDHTaZ6SUkZb7UNoWj7pGoWm5VHaAlFLgWVSSn/AH3BBrTWQkxDUgjuZz+UJuALXgJT0/VLKPVLK7pkeOkhoxUoHCk3Luw5AgpRyNoCUMhV4GXhCCFEm88FCiAbA4bTjMnsb+BQ1aVyjO11YCOGZVjtpJoR4sxCfQdPyTTc9aVreNQLCLXdIKa8KIU4CdVCzhFrqBqzOfBIhhB9q9tNXgLZp592S24WlWmvk2YIVW9MKR9coNM12upBNoAA+BMZLNXXzAfJQo9C04qRrFJqWd/uBvpY7hBDlgCqodQUs97sC5aWUZzPtDwT6AG2FEFNR64vssWGZNa3QdI1C0/JuHeAqhHgCQAhhBD4DvpJSxmc69j7UYjWZfYxa2MZPSukHNEXXKLQSTgcKTcujtKai3kBfIcQR1Ep3ZinlhGwOz5KfEEJ0AFyllH9anPMCUEYIUcF2Jde0wtEr3GlaAQkh2gALgN5Syh2ZXtsBtJJSJhdL4TTNinSg0DRN03Klm540TdO0XOlAoWmapuVKBwpN0zQtVzpQaJqmabnSgULTNE3LlQ4UmqZpWq50oNA0TdNy9f+qZeBel7qhBwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "global_objective.plot()\n",
    "plt.yscale('log')\n",
    "plt.xlabel('Q / $\\AA^{-1}$')\n",
    "plt.ylabel('Reflectivity')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can display out what the fit parameters are by printing out an objective:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_______________________________________________________________________________\n",
      "\n",
      "--Global Objective--\n",
      "________________________________________________________________________________\n",
      "Objective - 140312336727488\n",
      "Dataset = d2o\n",
      "datapoints = 137\n",
      "chi2 = 417.78640890990016\n",
      "Weighted = True\n",
      "Transform = None\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "________________________________________________________________________________\n",
      "Parameters: 'instrument parameters'\n",
      "<Parameter:    'scale'    , value=1.03184 +/- 0.00368, bounds=[0.9, 1.1]>\n",
      "<Parameter:     'bkg'     , value=-6.60703e-08 +/- 1.05e-07, bounds=[-1e-06, 1e-06]>\n",
      "<Parameter:'dq - resolution', value=5 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:  'q_offset'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "________________________________________________________________________________\n",
      "Parameters: 'Structure - ' \n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:  ' - thick'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - sld'    , value=2.07 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - isld'   , 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:       ''       \n",
      "<Parameter:'sio2 thickness', value=11.9604 +/- 0.23 , bounds=[2.0, 30.0]>\n",
      "<Parameter:   ' - sld'    , value=3.47 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - isld'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'sio2 roughness', value=4.64994 +/- 0.704, bounds=[0.0, 7.0]>\n",
      "<Parameter:'sio2 solvation', value=0.12718 +/- 0.00815, bounds=[0.0, 0.5]>\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:'area per molecule', value=57.2217 +/- 0.246, bounds=[52.0, 65.0]>\n",
      "<Parameter:   'b_heads'   , value=0.000601 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_heads_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_heads'   , value=319 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'inner_head_thickness', value=9.7611 +/- 0.161, bounds=[4.0, 11.0]>\n",
      "<Parameter:   'b_tails'   , value=-0.000292 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_tails_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_tails'   , value=782 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'tail_thickness', value=13.6668 +/- 0.156, bounds=[10.0, 17.0]>\n",
      "<Parameter:' - rough_head_tail', value=3 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - rough_fronting_mono', value=3 (fixed)  , bounds=[-inf, inf]>\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:'area per molecule', value=57.2217 +/- 0.246, bounds=[52.0, 65.0]>\n",
      "<Parameter:   'b_heads'   , value=0.000601 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_heads_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_heads'   , value=319 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'outer_head_thickness', value=6.06557 +/- 2.12 , bounds=[4.0, 11.0]>\n",
      "<Parameter:   'b_tails'   , value=-0.000292 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_tails_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_tails'   , value=782 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'tail_thickness', value=13.6668 +/- 0.156, bounds=[10.0, 17.0]>\n",
      "<Parameter:' - rough_head_tail', value=3 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - rough_fronting_mono', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:  ' - thick'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'd2o SLD'   , value=6.17663 +/- 0.00475, bounds=[6.1, 6.36]>\n",
      "<Parameter:   ' - isld'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'bilayer/solvent roughness', value=1.24812 +/- 5.02 , bounds=[0.0, 5.0]>\n",
      "<Parameter:' - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n",
      "\n",
      "\n",
      "________________________________________________________________________________\n",
      "Objective - 140312336636128\n",
      "Dataset = hdmix\n",
      "datapoints = 97\n",
      "chi2 = 116.09164939809303\n",
      "Weighted = True\n",
      "Transform = None\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "________________________________________________________________________________\n",
      "Parameters: 'instrument parameters'\n",
      "<Parameter:    'scale'    , value=1 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:     'bkg'     , value=5.5744e-08 +/- 1e-07, bounds=[-1e-06, 1e-06]>\n",
      "<Parameter:'dq - resolution', value=5 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:  'q_offset'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "________________________________________________________________________________\n",
      "Parameters: 'Structure - ' \n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:  ' - thick'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - sld'    , value=2.07 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - isld'   , 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:       ''       \n",
      "<Parameter:'sio2 thickness', value=11.9604 +/- 0.23 , bounds=[2.0, 30.0]>\n",
      "<Parameter:   ' - sld'    , value=3.47 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - isld'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'sio2 roughness', value=4.64994 +/- 0.704, bounds=[0.0, 7.0]>\n",
      "<Parameter:'sio2 solvation', value=0.12718 +/- 0.00815, bounds=[0.0, 0.5]>\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:'area per molecule', value=57.2217 +/- 0.246, bounds=[52.0, 65.0]>\n",
      "<Parameter:   'b_heads'   , value=0.000601 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_heads_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_heads'   , value=319 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'inner_head_thickness', value=9.7611 +/- 0.161, bounds=[4.0, 11.0]>\n",
      "<Parameter:   'b_tails'   , value=-0.000292 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_tails_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_tails'   , value=782 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'tail_thickness', value=13.6668 +/- 0.156, bounds=[10.0, 17.0]>\n",
      "<Parameter:' - rough_head_tail', value=3 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - rough_fronting_mono', value=3 (fixed)  , bounds=[-inf, inf]>\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:'area per molecule', value=57.2217 +/- 0.246, bounds=[52.0, 65.0]>\n",
      "<Parameter:   'b_heads'   , value=0.000601 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_heads_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_heads'   , value=319 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'outer_head_thickness', value=6.06557 +/- 2.12 , bounds=[4.0, 11.0]>\n",
      "<Parameter:   'b_tails'   , value=-0.000292 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_tails_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_tails'   , value=782 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'tail_thickness', value=13.6668 +/- 0.156, bounds=[10.0, 17.0]>\n",
      "<Parameter:' - rough_head_tail', value=3 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - rough_fronting_mono', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:  ' - thick'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - sld'    , value=2.07 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - isld'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'bilayer/solvent roughness', value=1.24812 +/- 5.02 , bounds=[0.0, 5.0]>\n",
      "<Parameter:' - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n",
      "\n",
      "\n",
      "________________________________________________________________________________\n",
      "Objective - 140312336633968\n",
      "Dataset = h2o\n",
      "datapoints = 104\n",
      "chi2 = 251.58933827912017\n",
      "Weighted = True\n",
      "Transform = None\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "________________________________________________________________________________\n",
      "Parameters: 'instrument parameters'\n",
      "<Parameter:    'scale'    , value=1 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:     'bkg'     , value=7.71044e-08 +/- 1.13e-07, bounds=[-1e-06, 1e-06]>\n",
      "<Parameter:'dq - resolution', value=5 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:  'q_offset'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "________________________________________________________________________________\n",
      "Parameters: 'Structure - ' \n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:  ' - thick'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - sld'    , value=2.07 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - isld'   , 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:       ''       \n",
      "<Parameter:'sio2 thickness', value=11.9604 +/- 0.23 , bounds=[2.0, 30.0]>\n",
      "<Parameter:   ' - sld'    , value=3.47 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - isld'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'sio2 roughness', value=4.64994 +/- 0.704, bounds=[0.0, 7.0]>\n",
      "<Parameter:'sio2 solvation', value=0.12718 +/- 0.00815, bounds=[0.0, 0.5]>\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:'area per molecule', value=57.2217 +/- 0.246, bounds=[52.0, 65.0]>\n",
      "<Parameter:   'b_heads'   , value=0.000601 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_heads_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_heads'   , value=319 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'inner_head_thickness', value=9.7611 +/- 0.161, bounds=[4.0, 11.0]>\n",
      "<Parameter:   'b_tails'   , value=-0.000292 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_tails_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_tails'   , value=782 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'tail_thickness', value=13.6668 +/- 0.156, bounds=[10.0, 17.0]>\n",
      "<Parameter:' - rough_head_tail', value=3 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - rough_fronting_mono', value=3 (fixed)  , bounds=[-inf, inf]>\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:'area per molecule', value=57.2217 +/- 0.246, bounds=[52.0, 65.0]>\n",
      "<Parameter:   'b_heads'   , value=0.000601 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_heads_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_heads'   , value=319 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'outer_head_thickness', value=6.06557 +/- 2.12 , bounds=[4.0, 11.0]>\n",
      "<Parameter:   'b_tails'   , value=-0.000292 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - b_tails_imag', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   'v_tails'   , value=782 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'tail_thickness', value=13.6668 +/- 0.156, bounds=[10.0, 17.0]>\n",
      "<Parameter:' - rough_head_tail', value=3 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:' - rough_fronting_mono', value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:  ' - thick'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - sld'    , value=-0.56 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - isld'   , value=0 (fixed)  , bounds=[-inf, inf]>\n",
      "<Parameter:'bilayer/solvent roughness', value=1.24812 +/- 5.02 , bounds=[0.0, 5.0]>\n",
      "<Parameter:' - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(global_objective)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's example the scattering length density profile for each of the systems:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAELCAYAAADeNe2OAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABC+ElEQVR4nO3dd3hUVfrA8e+ZmfReSSAkQVoogRCqhSaCoCJ2RGAFCyq6VpBV97frurqWde1YUMFGWRVEQGRVRAQL0mtIKFICSQghCQnpM+f3xwxICZAymZvMvJ/HeUjuPXPuO/OMb86ce4rSWiOEEMJzmIwOQAghhGtJ4hdCCA8jiV8IITyMJH4hhPAwkviFEMLDWIwOoCYiIyN1YmKi0WEIIUSTsXbt2sNa66jqzjWJxJ+YmMiaNWuMDkMIIZoMpdTes52Trh4hhPAwkviFEMLDSOIXQggP0yT6+KtTWVlJZmYmZWVlRodiKF9fX+Li4vDy8jI6FCFEE9FkE39mZiZBQUEkJiailDI6HENorcnLyyMzM5NWrVoZHY4Qoolosl09ZWVlREREeGzSB1BKERER4fHfeoQQtdNkEz/g0Un/OHkPhBC11WS7eoQQojHIP1bB1oNHySos5WhZFVVWG1U2TaXVRnWr3p91IfxqCvv7WLi7f2unxguS+J3mySefJDAwkJycHBYuXIi3tzetW7dmxowZhIaGGh2eEMLJMnKK+NfiNJZn5Fab4Ovi9C/wkYE+kvibgsGDB/Pss89isViYMmUKzz77LM8//7zRYQkhnGjZ9kPc9clafC0m/jywDX1aR9Ai1I9QP28sZoXFrPAymTCZGmdXrCF9/EqpUKXU50qp7UqpNKXUhUbEUV/PPPMM7dq145JLLiE9PR2AIUOGYLHY/5726dOHzMxMwH4zevz48SQnJ9OtWzeWLVtmWNxCiLrbl1fCfbPW0a5ZIN9PGsDDQ9pzUetIEiICCPH3IsDHgo/F3GiTPhjX4n8VWKK1vkEp5Q3416eyfyzcyraDR50TmUPH5sH8fXins55fu3Ytc+bMYcOGDVRVVZGamkr37t1PKTN9+nRGjhwJwNSpU1FKsXnzZrZv386QIUPIyMjA19fXqXELIRrWU4u2ATBtbA8iA30MjqZuXN7iV0qFAP2A9wG01hVa6wJXx1FfK1as4Nprr8Xf35/g4GCuvvrqU84/88wzWCwWRo8eDcDKlSsZM2YMAElJSSQkJJCRkeHyuIUQdbc9+yjfpeVwd//WNA/1MzqcOjOixd8KyAVmKKW6AmuBB7TWx04upJSaAEwAiI+PP2eF52qZG+GDDz5g0aJFLF26VIZbCuFGZv66Dx+LiT9dmGh0KPViRB+/BUgF3tJadwOOAX85vZDWeprWuofWukdUVLVLShuqX79+zJ8/n9LSUoqKili4cCEAS5Ys4YUXXmDBggX4+//Rg9W3b19mzpwJQEZGBvv27aN9+/aGxC6EqL0qq41Fmw4ypFMMIf5Ne4kUI1r8mUCm1nqV4/fPqSbxN3apqamMHDmSrl27Eh0dTc+ePQG47777KC8vZ/DgwYD9Bu/bb7/NxIkTueeee0hOTsZisfDBBx/g49M0+weF8ESbDhSSX1LJ0E4xRodSby5P/FrrbKXUfqVUe611OjAI2ObqOJzhiSee4Iknnjjl2KRJk6ot6+vry4wZM1wRlhCiAfy88zAAF7aOMDiS+jNqVM+fgZmOET27gfEGxSGEEDXy8648OsQGEx7gbXQo9WZI4tdabwB6GHFtIYSorbJKK2v25vOnPglGh+IUTXqRNiGEcIWN+wuoqLK5RTcPSOIXQojz2pZlnyDaJS7U2ECcRBK/EEKcx/asIiICvIkKco+ReJL4hRDiPLZnHyUpNsjoMJxGEn8d7dmzh86dO9e7zOn+9re/8d1339UnNCGEE1ltmvScIpJigo0OxWlkWeZG5qmnnjI6BCHESfbmHaOs0kZSjLT4BWC1Wrnzzjvp1KkTQ4YMobS0lLVr19K1a1e6du3K1KlTT5T94IMPuOaaaxg8eDCJiYm88cYbvPTSS3Tr1o0+ffpw5MgRAMaNG8fnn39OYWEh7du3P7Hc86hRo3j33XcNeZ1CeLLt2UUAdIiVFn/j8vVfIHuzc+uMSYZhz52zyI4dO5g9ezbvvvsuN910E3PnzuWFF17gjTfeoF+/fkyePPmU8lu2bGH9+vWUlZXRpk0bnn/+edavX89DDz3ERx99xIMPPniibEhICG+88Qbjxo3jgQceID8/nzvvvNO5r1EIcV7p2UUoBW2iA40OxWmkxV8PrVq1IiUlBYDu3buzZ88eCgoK6NevHwBjx449pfzAgQMJCgoiKiqKkJAQhg8fDkBycjJ79uw5o/7BgweTnJzMvffey3vvvdegr0UIUb39R0qIDfbF18tsdChO4x4t/vO0zBvKyYusmc1msrKyalzeZDKd+N1kMlFVVXVGeZvNRlpaGv7+/uTn5xMXF+ekyIUQNbX3SAnxEfXaK6rRkRa/E4WGhhIaGsrKlSsBTizDXFcvv/wyHTp0YNasWYwfP57KykpnhCmEqIV9R0qID5fEL85hxowZ3HvvvaSkpKC1rnM96enpvPfee/znP/+hb9++9OvXj6efftqJkQohzqe0wkpuUbnbJX5Vn+TkKj169NBr1qw55VhaWhodOnQwKKLGRd4LIRpGenYRl7/yI6+N6sbVXZsbHU6tKKXWaq2rXQxTWvxCCHEW+46UALhdi18SvxBCnMXxxJ8giV8IITzD/iMlBPlYCG3ie+yeThK/EEKcxb4jJcSF+6OUMjoUp5LEL4QQZ5FVWEaLUF+jw3A6SfxCCHEW2YWlxIRI4hcOZ1tyefLkySQlJdGlSxeuvfZaCgoKXB+cEKLeyiqt5JdUEhMsiV+cx+DBg9myZQubNm2iXbt2PPvss0aHJISog+zCMgBiQvwMjsT5DEn8Sqk9SqnNSqkNSqk1539G41TdssxDhgzBYrEvgdSnTx8yMzMBKCsrY/z48SQnJ9OtWzeWLVtmZOhCiPPIPmpP/LFu2NVj5CJtA7XWh51R0fO/Pc/2I9udUdUJSeFJTOk15ZxlqluWecyYMSfOT58+nZEjRwIwdepUlFJs3ryZ7du3M2TIEDIyMvD1db8PlRDu4I8Wv/v9PypdPfVQ3bLMxz3zzDNYLBZGjx4NwMqVK0/8UUhKSiIhIYGMjAxXhyyEqKGs44nfDfv4jWrxa+AbpZQG3tFaTzu9gFJqAjABID4+/pyVna9l3lBOX5a5tLQUsO+2tWjRIpYuXep243+F8BTZhaUE+VoI8HGP1etPZlSL/xKtdSowDLhXKdXv9AJa62la6x5a6x5RUVGuj7COlixZwgsvvMCCBQvw9/9jmnffvn1PLNOckZHBvn37aN++vVFhCiHOI/tomVv274NBiV9rfcDx7yHgC6CXEXE0hPvuu4+ioiIGDx5MSkoKd999NwATJ07EZrORnJzMyJEj+eCDD075xiCEaFyyC8vcckQPGNDVo5QKAExa6yLHz0OAp1wdR30lJiayZcuWE79PmjQJgCeffLLa8r6+vsyYMcMVoQkhnCCrsIz2MUFGh9EgjOi8agZ84ej7tgCztNZLDIhDCCGqZbVpDheX08wNb+yCAYlfa70b6Orq6wohRE3ll1Rg0xAZ6J7dsU16OGdT2D2socl7IITz5RaVAxAVJIm/UfH19SUvL8+jE5/Wmry8PJkEJoSTHS5278TfZAeoxsXFkZmZSW5urtGhGMrX15e4uDijwxDCrRxv8btrV0+TTfxeXl60atXK6DCEEG5IunqEEMLDHC4ux8/LTIC32ehQGoQkfiGEOE1uUTmRQd5uu+SKJH4hhDhNbnE5UW7avw+S+IUQ4gyHiyrc9sYuSOIXQogz5BaXu+2NXZDEL4QQp6i02jhyrEISvxBCeIojxyoA9x3DD5L4hRDiFO4+hh8k8QshxCly3Xy5BpDEL4QQpzjR4peuHiGE8Azuvk4PSOIXQohTHC4uJ8jHgp+bLtcAkviFEOIU9uUa3Le1D5L4hRDiFLlF7r1cA0jiF0KIUxwuti/Q5s4k8QshxEmkxd+AlFJmpdR6pdQio2IQQoiTlVdZOVpW5dZj+MHYFv8DQJqB1xdCiFMcLnb/5RrAoMSvlIoDrgTeM+L6QghRHU9YrgGMa/G/AjwK2M5WQCk1QSm1Rim1xtM3VBdCuMZhSfwNQyl1FXBIa732XOW01tO01j201j2ioqJcFJ0QwpMdX6dHunqc72LgaqXUHmAOcKlS6hMD4hBCiFN4wnINYEDi11o/prWO01onAjcD32utx7g6DiGEOF1uUTlh/l54W9x7pLt7vzohhKiF3CL33nLxOIuRF9da/wD8YGQMQghxnLvvtXucoYlfiOoUllSyM7eInYeKOVhQxpFjFRw5VsGxiiqsNo3WYDIpArzNBPhYCPSxEBHgTbMQX2JDfIkN8SMhwh8vs3yhFbWTW1ROanyo0WE0OEn8wnDlVVZWZBzmh4xDrNp9hB2Hik85H+rvRXiAN4E+FkxKYTYpqmyagwWllJRXUVReRVFZ1SnP8TabaBMdSIfYYC5qHUHfdpFEB/m68mWJJkZrLV09xymlBgM3AVO11huUUhO01tMaPjTh7rYeLOTDn/fw9ZZsisqqCPA20yMxnGu6tSApJog20YG0CPXDUoOWe1mllZyjZWQXlpGZX0pGThHbs4v4If0Qc9dlAtA9IYzRveO5skssPhb3XWtd1M2xCiullVZJ/A63AfcAf1VKhQMpDRqRcHu/7MrjtaU7+GV3Hv7eZq5IjuXKLrFc0iayzt0zvl5mEiICSIgIoPdJx202zbasoyzPyOXztZk8/OlGXv4ug8eGdWBY5xiUUs55UaLJ85RZu1CzxF+ktS4AJimlngN6NmxIwl3tyi3m2cVpfJd2iJhgXx6/IomRPeMJ8fNqsGuaTIrOLULo3CKEiQNa80NGLs8t3s7EmesY2imG56/vQoh/w11fNB1/7LXr/l2CNUn8Xx3/QWv9F6XUnxswHuGGqqw2pq3YzSvf7sDbYuLRoe257eJW+Hq5trtFKcXA9tH0bRPJ+yt/59//S+f6t3/mo9t60TzUz6WxiMbHk1r85/1erbX+EkApFen4/fWGDkq4j315Jdzw9i+8sCSdQR2iWTZpABMHtHF50j+ZxWzirv6t+eSO3uQUlnHTO79w6GiZYfGIxuFQkf0zIIn/VNMbLArhln7MyGX4GyvZnVvMqzen8Obo1Eb1P1WfCyL45I7eHDlWwa0zVlNaYTU6JGGg3KJyLCZFaAN2PTYWtUn8chdM1IjWmnd/3M24Gb8RG+LLwj9fwoiUFo3yRmrXlqFMHZ3K9uyj/O3LLUaHIwyUW1ROZKAPJlPj+5w6W20Sv26wKITbsNk0z3yVxjOL0xjaOYZ5Ey8iISLA6LDOaWD7aO4b2IbP1mby9eYso8MRBvGUWbsgLX7hRFVWG5M+28h7K39n3EWJvDEqFX/vpjFH8IFBbekYG8zfF2zlaFml0eEIA3jK5C2oXeJ/rMGiEE2e1aZ55LONzFt/gIcHt+Pvwzs2qa/MFrOJ565P5nBxOa8v3WF0OMIAnrDJ+nE1Tvxaa+kAFdWy2TSPz9vMlxsO8ujQ9tw/qG2j7M8/ny5xoVyXGseHv+zlQEGp0eEIF7LaNHnHKqTFXxtKqVBn1COaHq01/1i4lf+u2c/9l7Zh4oA2RodULw8NbgcgrX4Pk19SgdWmJfEfp5TqrpT6u1IqTCkVpJTqo5S6XSn1klLqf0qpA8Cehg9VNEavLd3Jh7/sZUK/C04kzaasRagfI3u0ZN66A+TI2H6Pceio50zegpq1+N8BFgH7gHTgn9jX69kJJAPdtNahDRSfaMTmrcvk5e8yuD41jseGJTXJ7p3q3NG3FVU2GzN+2mN0KMJFjv+Rbxbs/ss1QM0S/8/AZGAdcAB4V2v9Z631m0C51vpQQwYoGqdfduUxZe4mLmodwbPXJbtN0gdIiAhgWOdYZq3aK5O6PES2I/HHhkjiB0BrfT9wm9a6P3A50Ecp9YtSahgytt8j7cot5q6P15AYEcBbY7q75f6kYy9M4GhZFV/JuH6PkFVYhkl5TldPjQZZa61LHP8eAR5WSiUATwPNlFIDtdbLGjBG0YgUlVUy4aM1eJlNTB/Xs3Yra1aWwsH1cGAdFGVBaQF4B0BgFMT1tD+8G8dkr96twrkgKoBZq/ZyQ/c4o8MRDSy7sJSoIB+P2bWtTrNrtNZ7gbFKqf8AzymlnnR8IxBuzGbTPPLpRvbklfDJ7b1pGe5//idpDTu/gw2zIGMJVJbYj1v8wC8MKo9BWaH9mFcAdLkJ+kyEKGNvFCuluKVXPE9/lcbOQ0W0iQ4yNB7RsLIKy4jxkP59qOfWi1rrDcBQpdTAmj5HKeUL/Aj4OK7/udb67/WJQ7jG1GU7+WZbDn+7qiMXto44d2GbFTZ/Bj+9Coe2gX8EdL0Z2gy2t+wDo/4oW1oAmWtg6xewcTas/9ie/PtPAZ/ABn1N5zK8a3OeWZzGwo1ZPDRYEr87yy4s44KoxvFt0xVqsvXiWCATmAhUAT9qrd86uUwtu3rKgUu11sVKKS9gpVLqa631r7WoQ7jYsu2HeOm7DK5Jac74ixPPXXjvz/D1o5C9GaI7wrXvQOfrwXyWbiG/UGh7mf0x+B/w3ZPw82uQ8T+4eRZEGjM3oFmwL71bhbNw00EevKxpTkoTNZNdWMbFbSKNDsNlatKh1RO4Umt9o9Z6FJBUnwtqu+O7aXs5HnKTuBHbc/gY989ZT4eYYJ69rsvZE2BpAXxxN8wYBiX5cMN0uOdne0v/bEn/dAGRMOIN+NMCKDkM714Kv//otNdSW8O7Nmd37jHSsooMi0E0rKKySorKq4jxkBE9ULPEfxSIVErdqZS6Aaj39yGllFkptQE4BHyrtV5VTZkJSqk1Sqk1ubm59b2kqKNj5VXc9fFazCbFO2O74+d9lg1Udv8Ab10Mmz6FfpPhvtX2Vn5dW8kX9IcJP0Bwc5h5E+wyZvzAsM6xmE2KhZsOGnJ90fByPGwoJ9Qs8f8fMB8IB7yBem+9qLW2aq1TgDigl1KqczVlpmmte2ite0RFRZ1Rh2h4WmsenbuJHYeKeH1Ut+pv5lqr4Ju/wkcjwMsP7vgWLv0reNfgxu/5hMbDuEUQfgHMvhn2/lL/OmspPMCbS9pEsnDjQbSWL6buKLvQPmvXUyZvQc3G8Wut9Xzgfa31LK2101avcmzivgwY6qw6hfNM+3E3X23KYsrQJPq2reaP77HD8PE18PPr0ON2uOtHaNHduUEERMKtCyEkDuaMgrxdzq2/BoZ2jiEzv5SMnOLzFxZNzsFCe0rzpFE9Lt96USkVdXxRN6WUHzAY2O6MuoXzrNiRy/NLtnNll1gm9LvgzAIH18M7/SFzNVzzNlz1knNa+dUJiIDRn4Eywcwb7PcSXGhg+2gAvt8uk9TdUWZ+KUpBbKgk/uo4a0hDLLBMKbUJWI29j3+Rk+oWTrAvr4T7Zq2nXbMgXri+mpu5mz+H9y+399/ftgRSRtWoXq01Nm3Dpm21Dyr8AvsIn4J9sOA++/wAF4kJ8aVT82C+357jsmsK18k8UkJssC8+lrPcv3JDtRnH75T/07TWm4BuzqhLOF9JRRUTPl4DwDtjuxPgc9JHRGv7uPzv/g7xF8HIj+1dMWeRVZzFigMrWJ29moz8DA4UH6Dcau9PbR7QnLZhbekX14+hrYYS7B18/uDi+8Bl/4BvnoBVb0Ofe+r1Wmvj0qRopi7bSUFJBaH+3i67rmh4+46U1GwyohupTeKXQcxuTmvNo59vIj2niA/G9zp1r1yb1T42f/V70Ok6uOYt8Drzq3FZVRmLf1/Mgl0LWJuzFoBo/2g6R3Smb4u+BHgHUGWrYn/RfrYe3sryzOW8tPYlbu14K7cl34aP+TxrpVx4r32ewDf/B3G9IM7J9xTOYmBSNK9/v5PlGbmMSGnhkmsK19h3pIT+7TxrAEltEr9svejmpv24m0Wbsnh0aPtT/0eoKIG5t0P6Yrj4ARj0JJhO7SUsLC9kVtosZm+fTX55PonBidyXch+XJ15OQnBCtWP/tdZsO7KN9ze/z5sb3+R/e/7HywNfplVIq7MHqRRcMxXe7gfz7oS7VzbcvYWTdI0LJSLAm++3H5LE70ZKK6wcKionXlr81dNab1FKJQEjgOOf/APAAq11WkMEJ1znh/RD9pu5ybHc07/1HyeKc2H2SPvN3CtehF53nvK80qpSZqbNZPqW6RRVFNE/rj+3drqVHs16nHemq1KKThGdeGnAS6w8sJLHVzzOLV/dwluXvUVKdMrZn+gXZk/+Hw63dztd8e96vPKaMZsUfdtGsnLnYbTWMovXTWTm29eOio/wrMRf45u7SqkpwBzsXT6/OR4KmK2U+kvDhCdcYXv2Ue6btZ72McG8cMNJN3MP74T3L4OcbTBy5ilJX2vNwl0LuWreVby67lVSo1OZe/Vc3hj0Bj1jetY6MV7S4hLmXDWHSL9IJnw7gY25G8/9hFb9oPc98Ns0++QxF7ioTSSHiytkWKcb2XfEnvg9rY+/NqN6bgd6aq2f01p/4ng8B/RynBNN0KGjZdw2YzUBPmamj+vxx83cfavg/cFQXmyfRJV0xYnn7Du6jzu/vZPHVz5OTEAMHw79kDcGvUG7sPqtqNk8sDnTL59OpF8k939/P/uL9p/7CZf9HSLbwfyJLhnieZFjYbqfdh5u8GsJ19jvSPye1tVTm8RvA5pXczzWcU40MSUVVdz+4RoKSit5/9aexIb42U9s+xI+utrepXLHtxDXA4BKWyXvbX6P6xZcx9bDW/lr77/y8RUfk9os1WkxRflH8eagN6myVfHID49Qaa08e2EvP7j2bSjKhiUN/6UzLsyfhAh/ft6V1+DXEq6x70gp/t5mIgI8a6RWbRL/g8BSpdTXSqlpjscSYKnjnGhCyqus3PXxWrYeLOT1Ud3o3CLEfuKXN+HTWyGmC9z+rX38PLDh0AZuWngTr657lX5x/fjymi8ZmTQSk3L+xhWJIYk8dfFTpB1J4/X1r5+7cIvu0PcR+3LOaQudHsvpLmodyardeVRZpa3jDnYfLiYhIsDj7tnU5ubuEqVUO+xdOyff3F2ttZaNSZuQSquNP89az4odh3nhhi4M6tAMbDb7+Phf34QOw+G6d8HLj6KKIl5d9yqfpn9Ks4BmvH7p6wxoOaDBYxwUP4gb293IjK0z6BvXl54xPc9euN9k+yYvCx+Eln1OXevfyS5qHcHs3/ax5eBRUlqGNth1hGtkZBfRq1W40WG4XK2aa1prm9b6V631XMfjV621VSk1vqECFM5ls2kmf7aRb7bl8OTwjtzUo6V9B6zZN9uTfu974MYP0RZfvt37LSPmj+CzjM8Y3WE080fMd0nSP25yz8m0CGzB078+fe4uH4s3XDcNyotg0YMNOqv3QunndxtFZZUcLCyjbTPP22THWd/T/+GkekQDqrLamPTZRuZvOMjky9sz7uJWcHgHvDsIdi21D9cc9hzZpbncv+x+Hv7hYSL8Iph1xSym9JpCgJdrdyjys/jxl15/YXfhbmamzTx34egO9lVBty+CjXMaLKbIQB/aNQvkt9+PNNg1hGscH53V3gMTf427ehxr61R7CmjmnHBEQymrtPLn2ev5dlsODw9ux70D20DGN/aJWWZv+NMCquJ7M3Prh0zdMBWAR7o/wpiOY7CY6rVDZ70MaDmA/nH9eWvjW1x5wZVE+Z+jG+fCeyH9a/sM41Z97St6NoAeieEs3HAQq01jNnlW37A72ZFj31ynnQcm/tq0+JsBfwKGV/OQYQ6NWP6xCsbN+I1vHd079/dPsK+hP+tGCEuECT+wMSCQmxfdzItrXqRHsx7Mu3oe4zqPMzTpHzel5xQqrBVM2zTt3AVNZrjmTfvyEvMn2u9bNICeiWEUlVeRni27cjVl6TlF+HmZiQvzMzoUl6tN4l8EBGqt95722AP80CDRiXpLzy5ixNSfWLe3gJdHdmVce6t9UpZjDf0jt8zhqbQZjF08lvzyfF4e8DJTB00lLqhhWst10TK4Jde2vZbPd3zOgeID5y4c3goufwZ+Xw6r3jp32TrqkWC/Gbh2r3T3NGU7copp2ywQkwd+a6tx4tda3661XnmWc7c4LyThDFpr5q8/wHVv/kRppZU5d3bn2pK58E5fKNhH+U0f8n5CJ65cdANzd8xldIfRLLhmAZclXNYoh7ZN6DIBEybe2lCDZN59HCRdBd/+Dfb/5vRY4sL8aBbsw+o9+U6vW7iG1pptWUc9sn8fardIm2giCkoqeGL+Fr7alEWPhDDeGWgjYvE1cGgrVe2uYHGXobyx7U2yjmXRP64/D3V/iNahrc9br5FiAmIYmTSSmWkzuavrXbQMann2wkrBiKnwTj/4bBzctcK+mYuTKKXokRjOmj3S4m+q9uSVcORYBakJYUaHYgjnz74RhrHZNP9dvY9B/1nO/7Zk83RfXz6NmEbEnCupKCvgs4EPMNyngCfWvkioTyjvD3mfNwa90eiT/nHjOo3DrMx8uPXD8xf2C4WbPoRjufDFXU7v7++ZEMbBwjIOFDhtJ1LhQuv22r+tpcZL4hdNlNaa77fnMGLqT0yZu5n+ITms7vQ5Y9bcSM7u73i961AubxHNU3u+INQnlNcGvsacq+bQK7aX0aHXSrR/NMNbD2f+zvkcKatBa7t5Nxj6HOz8FpY/59RYeiTa+/ml1d80rd2XT5CPhbbRgUaHYojaDOe8EPhVaxfueSfOqaSiisWbs3l/5e/sy8rh5qCNTI9bhW/ebyyzhrKkXQo/VeShj6bRN64vozuM5sLYCxtlH35N3drpVubtmMfs7bO5N+Xe8z+hx21wYC0sfx4i2kCXm5wSR1JMEIE+FtbsyZf1+ZugdXvzSYkP9cgbu1C7Pv4/AVOVUhnAEmCJ1jq7YcISZ1NSUcVPO/P4Zms2v23eRveqDUz230xM6EbWeCv+7hvGr4kJVGIj1svC7e1v54Z2N9A8sLr19ZqeC0IuYGDLgczePpvxncbj73WeVRWVgqtegfy98OW9EBpv38KxnixmE93iQ1ktLf4m52hZJek5RVzeKcboUAxTm7V67gFwbMYyDPhAKRUCLMP+h+AnWbPHuSqtNjLzS0nPKmT77r3k7tuGKXc90ZaddPLNpHnEUdJ8vPm7jy8FJnvXQ0JQAje37MfQxKEkRyY36db92dzW+TaWfb2ML3Z+wegOo8//BIu3fX/g9y6zL00x/mv7TN966pEQzitLMzhaVkmwr1e96xOusSLjMFrDJW3Pvl+0u6v1qB6t9XZgO/CyUsoPGAjcCLwE9HBuePVTUVGOTVvRWmPTgLahtQ2btqHRoDU2rdE2K1qDDfv54+W1toHWaK3R2NA2x79ao202bGh7GcCmbWDTaDQ2+8XQ2LDZbI5d6jVWq5WySiuVlRWUlxVTUlZIWVkRpeVFFJcVU1JeREV5IRWV+ZRWFVKmj4G5FKu5nDwLHAiwkB9sdrw6hYUw2gbFMzA6hdSY7vSJ7UNMgPu3YlKiU+gS1YU52+cwKmlUzVYI9Q+HMXNh+lD46Bq4bYl9zH89dIsPRWvYtL/Qo5NIU7N0ew6h/l50M3qRvdN7zc/Wi25y/q3Yeg3n1FqXAosdjxpRSrUEPsI+E1gD07TWr9YnjrO5eGYqZQ3wpjUoH8cDMGsI1iaCzWHE+oYzKKgFLSI70iIiiYSQRNqGtsXL7JktzVFJo3hsxWP8cvAXLm5xcc2eFN4Kxn4BH1wBH42A2/4HwbF1jqGrI3Fs2J8vib+JsNk0y9Nz6d8uCou5FrlBa/v2o3t/hpytULgfjh2G0nywVoCtCqyVYKu0/+wsAdEweYfz6nMwYhx/FfCI1nqdUioIWKuU+lZrvc3ZFxpm6oTVVoVCYf9PwUk/K0AfP+roEjm9jH0pIk4qc/yZjtJKndSdojilluPlj9etTFhMJsxmM2aLL37eAfj6BBDoE0xoYDDB/iH4BkQSEBBDqG8YAV6et054TV2ecDkvrn6R2dtn1zzxAzTrCKPn2jea+fgaGLe4zmP8Q/y8aB0VwPp9BXV6vnC99fsLyDtWwaVJ0TV7QtlRWP0urJlhT/YAgc0grBVEtLZvVmTxAZMXmC32f00W+72lM1Rz7Ixyp/3u3TA7g9Up8SulTIBJa13rP21a6ywgy/FzkVIqDfv6/k5P/E/d+qmzqxSNhJfZixva3cC0TdPYX7T/3BO6ThfXHUbNgU+uh5k3wK0LwKduMzi7xYexbPsh2YC9iViw4QDeFhMD2p8n8WsN6z+2z/4uzYfWl8KAx6DNZRDU9NekrHU/iFLqPiAH2KuU2qSUuqOuF1dKJQLdgFXVnJuglFqjlFqTm5tb10sIN3ZT+5swKzNzttdhGeZWfeHGDyBrI8wZDVXldYohpWUoeccq2H9EJnI1dhVVNhZuymJwx2aE+J2ji/RYHsy8ERb8GaI7wZ3L7F2E3Ua7RdKHuk3gegRI1lq3AC4HLlZKPVnbSpRSgcBc4EGt9dHTz2utp2mte2ite0RFNdyOSqLpivaPZlDCIL7Y+QUllSW1ryDpCvvSDr8vty9Pba1932y3+FAA1u+XdXsau+UZuRw5VsH1qeeYd3EoDd4dCL//CMNegFsXQgvn7SndWNQl8RcDh+BEt83twHW1qUAp5YU96c/UWs+rQwxCAHBL0i0UVRTx1e9f1a2ClFFw+bP2/XrrsHtX+2ZB+HmZpZ+/CZi5ai9RQT70bXuWhuTO7+C9wVBVBuMXQ++7GmRETWNQl1f1FvCZUqqN4/d4oMbNLWXvCH0fSNNav1SH6wtxQrfobiSFJzF7+2zqPKn8won2fXvXfwwr/lOrp1rMJpLjQtiwv6Bu1xYusSu3mB/ScxnTOwGv6kbzbPsSZt0M4Yn2rp24RjUy3elqnfi11m8CM4H3lFL5wE4gXSl1o1KqbQ2quBgYC1yqlNrgeFxR2ziEAPvIqVFJo9iRv4M1OWvqXtHAJ6DzDfD90/ZdvGqhW8tQth08SnmVzF9srD76eQ/eZhO39I4/8+TGOfZVXFukwrivIMT9l+Co0/cYrfU8rfUAIApIBb4HLgLeqcFzV2qtlda6i9Y6xfGo8TwAIU53RasrCPEJYVbarLpXohRc/TrEdoG5d8KR32v81G7xoVRYbWw7eMatKtEIHC2r5PO1mVzVNZaoIJ9TT6790L56a2Jf+w1c3xBjgnSxenVgaa2rtNabtNYfaq0f0lpf6qzAhKgpX4sv17e9nu/3f09WcVbdK/L2h5GfgDLBvDvtE3JqIKWlfWlf6edvnD5dvZ9jFVbGX3TaTO2t82HhA/Yhmrd8Ct4BhsRnBPe8cyE8zs3tbwZgTnodhnaeLDQehr8Mmath+Qs1ekpMiC+xIb7Sz98IWW2aj37ZS4+EMJLjTmrN715u/+Peshfc9DF4+RoXpAEk8Qu3EBsYy6D4QczdMZeyqrL6Vdb5eug6yn6jN3tLjZ6S0jJUhnQ2Qt9vP8S+IyWMv/ik1n5uhn3uRnhr+0S+Bpod25hJ4hduY1TSKArLC1n8uxNuGV3+L/suXoserNHuXd3iQ9l/pJTDxXWbCCYaxke/7CE2xJchnRwTr8qL4b9j7MssjPncvnifB5LEL9xGj2Y9aBfWjplpM+s+tPM4/3AY8oy9y2fdB+ctfryff4P08zcavx8+xoodhxnVK/6PIZwLH4C8HXDDdAiJMzZAA0niF25DKcXoDqPJyM9gbc7a+lfY9WZIuMQ+xLPs3CN2kluEYDYp6edvRGb+uheLSXFzT8c6Tlu/gC2fw4DH4YL+xgZnMEn8wq0cH9o5M21m/StTCob8E0ry4OfXz1nUz9tMUkyQ9PM3EmWVVj5bm8nlnWOIDva1L6H81SP2fZgvecjo8AwniV+4FV+LLze1u4ml+5bye2HNx+KfVYtU6HQd/PIGFOWcs2i3+FA27i/EapNtqY22cONBCksrGdM7wX5gyWP2b20j3rQvn+zhJPELtzO6w2i8zd7M2DLDORVe+lf7ZhsrXjxnsZSWYRSXV7Ert9g51xV19smve2kTHUifC8Jh/2+w+VO4+AH7fgxCEr9wPxF+EVzb5loW7l5I9rFsJ1TY2j68c91HUHzorMWOr9S5bq909xhpR04RGzMLGdUr3r6tyZLHIDBGunhOIolfuKVxncehtebDrR86p8JLHrK3+n9546xFLogMINTfi3X7JPEbacHGg5gUDO8aC1vmwoE1MOj/wCcQq83KlsNb2F242+gwDSWJX7ilFoEtuKLVFczdMZf8Mick4ojW0OlaWP0+lByptohSitT4MNbJkE7DaK35csNBLm4TSXSAF/zwLDTrDF1vYXfBbq5dcC2jvhrFiPkjeHzF41Q5c3/cJkQSv3BbdyTfQVlVGe9vft85FV7yMFQUw9qz3ztIjQ9l56FiCktqts6PcK4N+wvYd6SE4V2bw5Z5kLcT+j9KesEORi8eTWF5Ic9c8gy3d76dhbsX8uaGN40O2RCS+IXbuiD0Aoa3Hs7s7bOd09cf0xla9YfV08+6W1dqvGPBNhnWaYgFGw/ibTExtGMU/PgCRHckP/ES7v/+fvwt/sy5cg5Xt76aB7s/yPALhjNj6wwyizKNDtvlJPELtzYxZSIazdsb33ZOhb3vgqOZkF79jl9dW4ZiUkh3jwGsNs2iTVkMbB9F8O7FcDgD+k3m6d/+RW5pLq9d+hqxgbEnyt+fej9omLW9Hst5N1GS+IVbaxHYgpHtRzJ/53wy8jPqX2G7ofYVPFdNq/Z0gI+F9jHBrJcbvC736+48covKGZHSAn59E8JasTQwiG/2fsM9Xe+hU2SnU8rHBMRwWcJlfLnzSypruAS3u5DEL9zeXV3uItA7kH+t+lf91/AxmaHnnbB35VlX7uyeEMr6fQUykcvFvtxwgEAfC4OC9kLmasp73ckLa16kXVg7xnUeV+1zhrceztGKo/x88GfXBmswSfzC7YX6hvJA6gOszVlb903ZT9ZtDFj84LfqN5xLjbdP5NpxqKj+1xI1Ul5l5est2Qzp1AyfNdPAJ5iZfiYOHjvI5J6T8TJ5Vfu8C2MvJMArgOWZy10csbEk8QuPcF2b6+gU0YkXV79Y/+Gd/uGQfANs/rzaxduO3+Bdt7egftcRNfZDei5FZVXc2M4E277kSMpI3t32Ef3i+tEnts9Zn+dl9qJns56sylrlwmiNJ4lfeASzycw/LvoHhRWFPP3r0/Xv8uk+DipLYOu8M04lRPgTHuAtE7lcaMGGg0QEeNM7dx5oGzOCAiipKuGR7o+c97m9Y3uzr2gfB4sPuiDSxsHliV8pNV0pdUgpVbOtjYRwkvbh7bk35V6+2ftN/TdradEdojrAuo/POPXHRC5J/K5QVFbJd2k5XN05AtP6j8hvN4T/7v2aYa2GcUHoBed9fu/Y3gAe1eo3osX/ATDUgOsKwbhO40iJSuGpX56q37R9pSB1rH05gENpZ5xOTQhld+4x8o9V1CNaURPfbsuhvMrG2NAtUJLHJ9FxlFWVcWfynTV6fpvQNkT4RvBr1q8NHGnj4fLEr7X+Eah+zrsQDcxisvDv/v/G1+LLg8sepLiiHitpdrkZTF7VtvplIpfrfLnhIHFhfrTa9zlHw+KZlfMzlyVcRuvQ1jV6vlKKbtHd2JS7qYEjbTwabR+/UmqCUmqNUmpNbm6u0eEINxITEMOL/V9k39F9TP5xMpW2Oo7hDoiApCth42yoOnWv3ZSWoXiZFat+lzZOQ8orLmflzsOMbW9D/b6c2YkpFFcWM6HLhFrV0zmyM5nFmRSUFTRMoI1Mo038WutpWuseWuseUVFRRocj3EzPmJ78X5//Y+WBlfztp79h0+ffUL1aqWOh9Aikn3rPwNfLTErLUH7dLYm/IS3enIXVprlefU+5ycTMkt/pF9ePpPCkWtXTObIzAFvztjZEmI1Oo038QjS069tdzwOpD7Bo9yL+vfrfdRvpc8FACI6zr9V/mt6tIthyoJDics9cAdIVvtxwkI7RfkTu+JzFF/Qkv6KQWzveWut6OkbYN2jZctgzxpxI4hce7fbOtzO241g+SfuEF1a/UPvkbzJDt9GwaxkU7D/lVJ8LIrDaNGv2SKu/IWTml7Bmbz73tdyFLs7hEx9oF9aOnjE9a11XkHcQicGJbMmTxN8glFKzgV+A9kqpTKXU7a6OQYjjlFJM7jGZMR3G8EnaJzz969O17/ZJGQ1oe1//SVITQrGYpJ+/oSzcmAXAwGNfszq8ORklWYzuMBqlVJ3q6xzZmW2HtzkzxEbLiFE9o7TWsVprL611nNbaSYulC1E3Sike7fkot3W+jU8zPuXvP/8dq81a8wrCEuzLNa//GGx//NHw97bQtWUov+7Oa4CoxZcbDnBZi0r89i7jk2bxhPmEcUWrK+pcX1J4EodKDzln455GTrp6hMCe/B9MfZCJXScyf+d8Hl9Zy92ZUv8EBftgz4pTDvduFc6mzEKOST+/U2XkFLE9u4h7Q35hv1nxQ1kWN7S7AV+Lb53rbBvaFoAd+TucFWajJYlfCAelFPek3MMDqQ+w+PfFPLbisZon/6QrwTcE1n9yyuHj/fxrZQN2p1qw4SAWZaNL7kJmteyIWZm5OenmetXZLrwdADsKJPEL4XHuSL6Dh7s/zJI9S3j0x0drNs7fyw+Sb4S0BVBacOJw94QwzCYl3T1OpLXmy40HuKvFHsqKDzLfXMbgxMFE+0fXq94I3wjCfMKcs29DIyeJX4hqjO88nkk9JvHt3m+ZvHxyzTbq6DYWqspgy+cnDgX4WOgSF8LPuyTxO8v6/QXsP1LKKMsy5oc3o9haztgOY+tdr1KKdmHtpKtHCE92a6dbmdJzCkv3LeX/fv6/84/2ie0KzZLPWMKhb5tINmUWUFAi6/Y4w4INB2lhKSQ25wdmhYXRJaoLyVHJTqm7bVhbdhbsrPuEviZCEr8Q5zCm4xju73Y/X+3+iv+s+c+5Cytl36QlawNkbz5xuH/7aGwaVuw43LDBeoBKq40FGw8ypdkaVvp6sc96jDEdxjit/rZhbSmtKnX7Ddgl8QtxHnck38EtSbfw0baP+HDrh+cu3OUmMHvD+pknDqW0DCXEz4vlGbLmVH0tT88l/1gZg8v+x8cx8UT7R3NZwmVOq79dmOMGr5t390jiF+I8lFJM6TWFIQlD+M+a/7B8/zm26fMPt4/w2TTnxMJtZpOib9tIlmfkYpN9eOtl3vpMrvTfzr6KHFZRxqikUWfdVrEuWoe2RqHc/gavJH4hasCkTDx9ydMkhScxZcUUdhecYy3/bmOhNP+UhdsGtI8mt6ictOwzt2oUNVNYUsl32w5xb9AKPomIws/sy43tbnTqNfwsfsQHx7v9kE5J/ELUkJ/Fj9cufQ0fsw9//v7PHK04SxK/YIB94baTxvT3axcJ2PeGFXXz1eYsQqx5RBX/zFd+3lzdZgQhPiFOv07b0LbS1SOE+ENMQAwvD3iZg8UHefLnJ6tf1M1khpRbYOdSKLTfJIwO8iW5RQjfbstxccTuY966TO4O+YVPg/ypxMboDqMb5Dptw9qy9+heSqtKG6T+xkASvxC1lNoslftT7+fbvd8yJ31O9YW6OZLS2g9OHBqWHMOG/QUcKHDfhNJQdh4qZv3ew4zgW/4bFkb/uP60CmnVINdqG9YWjT53d14TJ4lfiDq4tdOt9G3Rl3+v/jdpeWfuuUtYIrQfBmumQ2UZAMM6xwLw9eYsF0bqHmat2seVltWs8CrmCDbGdqz/hK2zOb5mjzvf4JXEL0QdmJSJZy55hjDfMCYtn1T93r197oGSPNj8GQCtIgPoEBvM11uyXRxt01ZWaWXu2v08GPg/PgyPoF1YO3rF9Gqw67UMaomv2detb/BK4heijsJ8w3ih3wscKD7AU788dWZ/f2JfaNYZVr0NjnNXdI5h7d58sgvLDIi4afpqUxaty7ex07SPXWb7vIq6rrlfE2aTmdahrd36Bq8kfiHqoXuz7tybci9f7/maeTvmnXpSKeh9F+RsObFc87Bke3fPV9LdUyNaaz7+dS/3+/+PaeHhJAbFMyRhSINft22Ye4/skcQvRD3dnnw7fWL78Nxvz7Ezf+epJ5NvBP8I+OVNANpEB9IlLoRPV++v2x6/HmbN3nwOZ2Zg89lIupeZO7vehdlkbvDrtg1tS15ZHkfK3HP3NEn8QtSTSZl4tu+z+Hv5M2n5pFOHAXr5Qa8JkPE1ZG0C4Oae8aTnFLF+f4ExATch7yzfzcM+XzItNIQW/jEMazXMJddtG+bem7JI4hfCCSL9Inmu73PsLtzNc789d+rJ3neDTwgsfx6A4V1j8fMy89/f9ldTkzhu56Fi0rdvIsT/N7b6eHFH17ucujzDuRxP/O46skcSvxBOcmHzC7kj+Q7m7ZjH4t1/LNeAXyj0uRu2L4LsLQT5ejG8aywLNx2kWLZkPKu3l+/iPq/5vBIeQmJgHCPajHDZtSP9Ign3DSf9SLrLrulKhiR+pdRQpVS6UmqnUuovRsQgREOYmDKRbtHd+Mcv/2Df0X1/nOhzD/gEw/f/BGBUr3hKKqx8ulpa/dXZkVPExvWr0CHr2eNl4aGek13W2j+uY0RHth3Z5tJruorLE79SygxMBYYBHYFRSqmOro5DiIZgMVl4od8LeJm9mLR8EhVWx+YrfmHQ92HIWAK7vqdbfBg9E8N4f+XvVFrde9OPuvj3ku1M9vmIN8OC6R7ZhYEtB7o8hk4RndhVsMstl24wosXfC9iptd6tta4A5gCu+w4nRAOLCYjh6YufJu1IGv/89Z9/jN7pM9E+o3fJ42Ct4u7+rTlQUMqXGw4aGm9js3ZvPrb0r/kp4gBHzWam9Plrg47bP5tOEZ2waZtbdvcYkfhbACd/v810HDuFUmqCUmqNUmpNbq6saCialgEtB3B317uZv3M+H2z9wH7Q4gNDnoHcNFj1NgPbR9OpeTAvf5tBWaXV0Hgbi0qrjX/O+41rg2cyLyiQP3X8Ex0iOhgSS8cIe0fE1rythly/ITXam7ta62la6x5a6x5RUVFGhyNErd3T9R4uT7ycl9e+zNJ9S+0Hk66EdsNg6VOYDqfz+BUdOFBQyse/7DU22EZi2o+7uTp/Kq9HmmjpG8U93e41LJZo/2gi/SLZlud+/fxGJP4DQMuTfo9zHBPCrZiUiacvfprOkZ2Z8uMUfsv6zT6b9+rXwCcQvpjAxa1CGNg+ile+y/D4VTt35BSx4fs5bIzZSo6XF/8a+BJ+Fj/D4lFK0TmiM5tyNxkWQ0MxIvGvBtoqpVoppbyBm4EFBsQhRIPztfgyddBUWga15L7v72P9ofUQGA1XvQJZG2HxZJ66uhMaeGzeZo/dmvFYeRXPfLSA1PDpfB/gz0OpD5ISnWJ0WHRv1p09R/eQW+Je3c0uT/xa6yrgPuB/QBrwqdba/TrRhHAI8w3j3SHv0sy/GRO+mWDfs7fj1XDJw7B2Bi3TZ/CXYUn8mJHL2z/uMjpcl6uy2nhi5nIu4Z+8GR7A5c37MrbzeKPDAqBHTA8A1uasNTgS5zKkj19rvVhr3U5r3Vpr/YwRMQjhSpF+kXww9ANah7bmgWUP8Mm2T9AD/wodR8A3f2Ws6Ruu6hLLi/9LZ8kWz1nAzWrTPD1nGR0PP8DrURYuDOvIvy59xZBRPNVJCk8iwCuA1dmrjQ7FqRrtzV0h3E2EXwTTL59O37i+PL/6eR7+cRJHhj0P7a9AfT2Zl6MWkRoXxJ9nr2eJB6zZf6y8iqdmzMP70H283kzRLaQNrwybgbfZ2+jQTrCYLKRGp7I6RxK/EKKO/L38eW3ga0zqMYkf9v/A1V/dyKcpV1OZcgteP/2HOb7PMrBZGffMXMur3+2gyk0nd63dk8dLrz3I3srHmB1uYUTsJUy7+lP8vfyNDu0MvWN783vh72QWZRoditOoprA0bI8ePfSaNWuMDkMIp9pVsIunf32aNTlraB7QnLEhHbly7WeEVlbwY/BwpmT1Jyy2FY8MbsegDtGNpvujPtIP5rN48VQyS+bwQxB4KxOTUh/iuuRxjfb1ZRZlMmzeMB7u/jDjG8m9h5pQSq3VWveo9pwkfiGMo7VmxYEVvLvpXTbkbsCizPQ0B9Erdy/dS8soqmrP0rJuZAankpzSi75JMXRqHoKvV8OvSe8MZZVW0nan89v6z9id+xP7vPezxdeCRWuujEjl/oEvEB0YY3SY53XzoptRKGZfNdvoUGpMEr8QTUD6kXS+2v0VKw6sYGfBHxu6NKuqokVVFWFVGi+rDzabP8rkh9kShLfFH1+zDxaLNxaLN8rshUmZUcqESSlMJoVJmRytaQVa2f/B8dDa8bMNtA17NrA5toq0oTVobTte2lFO23+y2bBqGzar/d8Kayll1mOUWkso06Uco4QjlipyLCaqHK35BKsXl7foz00XTaZZUHPXvsH1MH3LdF5e+zKLr11My+CW539CIyCJX4gmJq80jw25G9hdsJvduZvJKdhNbukR8qqOUYINa+PsFcGkNf42TYBWBNq8iDQFEesXQ2pCX/p2uZ7IJpTsT5Z9LJth84ZxQ9sbeKLPE0aHUyPnSvwWVwcjhDi/CL8IBsUPYlD8oGrPV9mqKLeWU1ZxjPKKo5SXHaOktJSqyjKs1iqs2orVasNqs1Fls55ozWtsmAB14mECNEpZUEqhlAkc3xZw/KyUGZPjZ1AnvlGYzWa8vSz4WCx4WyyEBkXi7xfeaPvq6yMmIIYRrUcwb8c8JnSZQJR/015GRhK/EE2QxWTBYrIQ4BUAAdFGh+MRbu98O/N3zueVda/wzCVNe/qRDOcUQogaaBncktuTb2fBrgW8vfFto8OpF2nxCyFEDU3sOpGs4iymbpjK5sObGdNhDF2jujbK+QfnIolfCCFqyGwy88+L/0nbsLZM2zSNHzN/BCDIO4gArwDMyozFZMGszj/cVnP+gTVhPmF8OOzDesd9Okn8QghRC2aTmfGdxzOy/Uh+y/6NnQU7yTmWQ5m1jCpbFVablSpdheL8N7nPdyM80CvQWWGfQhK/EELUgb+XPwNaDmBAywFGh1JrcnNXCCE8jCR+IYTwMJL4hRDCw0jiF0IIDyOJXwghPIwkfiGE8DCS+IUQwsNI4hdCCA/TJNbjV0rlAntddLlI4LCLrtVUyHtyKnk/ziTvyZmMfk8StNbVrh/dJBK/Kyml1pxt8wJPJe/JqeT9OJO8J2dqzO+JdPUIIYSHkcQvhBAeRhL/maYZHUAjJO/JqeT9OJO8J2dqtO+J9PELIYSHkRa/EEJ4GEn8QgjhYSTxOyil/q2U2q6U2qSU+kIpFXrSuceUUjuVUulKqcsNDNOllFJDHa95p1LqL0bHYwSlVEul1DKl1Dal1Fal1AOO4+FKqW+VUjsc/4YZHasrKaXMSqn1SqlFjt9bKaVWOT4r/1VKeRsdoysppUKVUp87ckiaUurCxvwZkcT/h2+BzlrrLkAG8BiAUqojcDPQCRgKvKlUDTbUbOIcr3EqMAzoCIxyvBeepgp4RGvdEegD3Ot4H/4CLNVatwWWOn73JA8AaSf9/jzwsta6DZAP3G5IVMZ5FViitU4CumJ/bxrtZ0QSv4PW+hutdZXj11+BOMfPI4A5WutyrfXvwE6glxExulgvYKfWerfWugKYg/298Cha6yyt9TrHz0XY/4dugf29OL4L9ofANYYEaAClVBxwJfCe43cFXAp87ijiae9HCNAPeB9Aa12htS6gEX9GJPFX7zbga8fPLYD9J53LdBxzd576us9KKZUIdANWAc201lmOU9lAM6PiMsArwKOAzfF7BFBwUsPJ0z4rrYBcYIaj++s9pVQAjfgz4lGJXyn1nVJqSzWPESeVeQL71/uZxkUqGhulVCAwF3hQa3305HPaPibaI8ZFK6WuAg5prdcaHUsjYgFSgbe01t2AY5zWrdPYPiMWowNwJa31Zec6r5QaB1wFDNJ/THA4ALQ8qVic45i789TXfQallBf2pD9Taz3PcThHKRWrtc5SSsUCh4yL0KUuBq5WSl0B+ALB2Pu3Q5VSFker39M+K5lAptZ6leP3z7En/kb7GfGoFv+5KKWGYv/6erXWuuSkUwuAm5VSPkqpVkBb4DcjYnSx1UBbx2gNb+w3uBcYHJPLOfqv3wfStNYvnXRqAXCr4+dbgS9dHZsRtNaPaa3jtNaJ2D8T32utRwPLgBscxTzm/QDQWmcD+5VS7R2HBgHbaMSfEZm566CU2gn4AHmOQ79qre92nHsCe79/Ffav+l9XX4t7cbTqXgHMwHSt9TPGRuR6SqlLgBXAZv7o034cez//p0A89iXDb9JaHzEkSIMopQYAk7TWVymlLsA+ACAcWA+M0VqXGxieSymlUrDf7PYGdgPjsTesG+VnRBK/EEJ4GOnqEUIIDyOJXwghPIwkfiGE8DCS+IUQwsNI4hdCCA8jiV8IITyMJH4hhPAwkviFaABKqdeVUuuUUj2NjkWI00niF8LJHCszRgN3YV/7SYhGRRK/EHWklHpbKXXx6ce11seAWOAH4DVXxyXE+UjiF6Lu+mDftOcUSqkIwB8owr6+kxCNiiR+IU6jlLpbKbXB8fhdKbWsmjIdgAyttbWaKv4KvAhsxb5lpxCNiiR+IU6jtX5ba50C9MS+1vpL1RQbBiw5/aBjl66LgP9i36ZREr9odCTxC3F2r2Jfb35hNecup5rEDzwNPOXYyEcSv2iUPGoHLiFqyrEbWwJwXzXn/IFQrfXB046nANcBlyilpmLfoWpzgwcrRC1J4hfiNEqp7sAkoK/W2lZNkYHYd5w63fPYd3D7zlFPM+ybkgjRqEhXjxBnug/7TlLLHDd43zvt/Bn9+0qpSwH/40kfQGudAwQqpcIbOmAhakN24BKilpRS64DeWutKo2MRoi4k8QshhIeRrh4hhPAwkviFEMLDSOIXQggPI4lfCCE8jCR+IYTwMJL4hRDCw0jiF0IID/P/DwA9yPezd9wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(*s_d2o.sld_profile(), label='d2o')\n",
    "ax.plot(*s_hdmix.sld_profile(), label='hdmix')\n",
    "ax.plot(*s_h2o.sld_profile(), label='h2o')\n",
    "\n",
    "ax.set_ylabel(\"$\\\\rho$ / $10^{-6} \\AA^{-2}$\")\n",
    "ax.set_xlabel(\"z / $\\AA$\")\n",
    "ax.legend();"
   ]
  }
 ],
 "metadata": {
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}