{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7aa81dbc",
   "metadata": {},
   "source": [
    "# Structural modeling and experiment planning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b385d3c",
   "metadata": {},
   "source": [
    "What are the accessible structural details of a given system by phase measurements? How to perform data acquisition? Identifying proper experimental conditions is fundamental for successfully exploiting the X-ray dynamical diffraction (in its current conception).\n",
    "\n",
    "Our approach for identifying those suitable multiple-diffraction cases is by elaborating structural models and comparing their structure factor phases. So, this step-by-step tutorial will show you how to use **pyddt** to perform these tasks."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abe3434c",
   "metadata": {},
   "source": [
    "## CeFe$_{4}$P$_{12}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32611247",
   "metadata": {},
   "source": [
    "For illustrative purposes, let's consider the filled skutterudite CeFe$_{4}$P$_{12}$ and search for suitable experimental conditions to resolve the Ce oxidation state. Download its CIF on [Materials Project mp-16272](https://materialsproject.org/materials/mp-16272/).\n",
    "\n",
    "This tutorial aims to show the experiment planning performed by Morelhão, S. et al. during the investigation of vibrational dynamics of filled skutterudites. See [Phonon scattering mechanism in thermoelectric materials revised via resonant x-ray dynamical diffraction](https://link.springer.com/article/10.1557/mrc.2020.37) for a detailed reference."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "846f9dfa",
   "metadata": {},
   "source": [
    "### 0. Importing packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6d6c6d7c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5cdf1ef1f3a04ae1beb4d5364a74773b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rcParams['font.size'] = '14'\n",
    "\n",
    "import sys\n",
    "sys.path.append('path/to/pyddt')\n",
    "\n",
    "import pyddt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ec905e5",
   "metadata": {},
   "source": [
    "## 1. Structural modelling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fe9640c",
   "metadata": {},
   "source": [
    "The structural modeling within **pyddt** is based on the `Structure` class. \n",
    "\n",
    "So, in a few words, it follows what we will do in this section: after converting the CIF into a *.in* file, we will instantiate a structure. Then, carry out the desired changes in the unit cell (replace the atoms with ions, vary their B-factors) and save the resultant models."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42e87098",
   "metadata": {},
   "source": [
    "### 1.1 Converting CIF into .in file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "590fda76",
   "metadata": {},
   "outputs": [],
   "source": [
    "pyddt.to_in('Ce(FeP3)4.cif')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84ca1d94",
   "metadata": {},
   "source": [
    "### 1.2 Structure object"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05b00707",
   "metadata": {},
   "source": [
    "It's easy to generate a `structure` object from the *.in* file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5087ffe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfp = pyddt.Structure('Ce(FeP3)4.in')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c6fe15f",
   "metadata": {},
   "source": [
    "Let's check the current oxidation states and B-factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "895aef2f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Ce', 'Ce', 'Fe', 'Fe', 'Fe', 'Fe', 'Fe', 'Fe', 'Fe', 'Fe', 'P',\n",
       "       'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P',\n",
       "       'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P'], dtype='<U32')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cfp.atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fa671aa5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n",
       "       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cfp.bfactors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62fee2f2",
   "metadata": {},
   "source": [
    "Now, we will generate two simple models from this `structure` object: one presenting Ce$^{3+}$ and the other Ce$^{4+}$. Both including Fe$^{2+}$, P$^{1-}$, and all atomic displacements equal to 0.1 $\\mathring{A}$ $\\rightarrow$ $B = 0.7896\\;\\mathring{A}$$^{2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44de2705",
   "metadata": {},
   "source": [
    "#### 1.2.1 Varying the B-factors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a05dd1e",
   "metadata": {},
   "source": [
    "Let's start by changing all B-factors. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "524a19c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfp.replace_bfactor(':', 0.7896)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07f59e7d",
   "metadata": {},
   "source": [
    "The first argument is the indices of B-factors to be replaced, while the second is the new value. In this case, we replace simultaneously all B-factors using `':'`. As a sanity check, verify the current values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e965bfc0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896,\n",
       "       0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896,\n",
       "       0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896,\n",
       "       0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896, 0.7896,\n",
       "       0.7896, 0.7896])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cfp.bfactors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4ebd1fb",
   "metadata": {},
   "source": [
    "#### 1.2.2 Replacing Fe and P atoms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67b08adc",
   "metadata": {},
   "source": [
    "The procedure for replacing the atoms with the desired ions is very similar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "895e590d",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfp.replace_ion('Fe', 'Fe2+')\n",
    "cfp.replace_ion('P', 'P1-')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a5a432f",
   "metadata": {},
   "source": [
    "We said for **pyddt** to replace all `'Fe'` atoms with `'Fe2+'` ions, and the corresponding for phosphorous. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "82c0e4ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Ce', 'Ce', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+',\n",
       "       'Fe2+', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-',\n",
       "       'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-',\n",
       "       'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-'], dtype='<U32')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cfp.atoms  # Sanity check"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ba8ee1d",
   "metadata": {},
   "source": [
    "#### 1.2.3 Ce oxidation state"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df5a64a4",
   "metadata": {},
   "source": [
    "Now, we will change the Ce atom to Ce$^{3+}$ and save this structural model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9fbfc6b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfp.replace_ion('Ce', 'Ce3+')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d5e2d749",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Ce3+', 'Ce3+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+',\n",
       "       'Fe2+', 'Fe2+', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-',\n",
       "       'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-',\n",
       "       'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-'],\n",
       "      dtype='<U32')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cfp.atoms  # Sanity check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1dde362a",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfp.save_infile('ce3.in')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f78d03ef",
   "metadata": {},
   "source": [
    "The structural model `ce3.in` is saved in the current folder. Check this file. Are the ions and B-factors correct? If not, review this tutorial from section 1.2.\n",
    "\n",
    "By the end, replaces the Ce$^{3+}$ ion with Ce$^{4+}$ and save the other structural model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8bcf638f",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfp.replace_ion('Ce3+', 'Ce4+')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "71b7073a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Ce4+', 'Ce4+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+',\n",
       "       'Fe2+', 'Fe2+', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-',\n",
       "       'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-',\n",
       "       'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-'],\n",
       "      dtype='<U32')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cfp.atoms  # Sanity check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3b38f7a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfp.save_infile('ce4.in')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03d81fb1",
   "metadata": {},
   "source": [
    "Again, it's highly recommend that you check the saved model before proceeding. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42e8114f",
   "metadata": {},
   "source": [
    "## 2. Experiment planning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c002141",
   "metadata": {},
   "source": [
    "We have generated two structural models by changing the Ce valence from 3+ to 4+. Now, it is necessary to search for structure factors with phases susceptible to this change in the models. \n",
    "\n",
    "The function `pyddt.phase()` finds the susceptible cases from the structure factor list of both models. As we see in the previous tutorial, the structure factor list is obtained after creating `crystal` objects."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba359885",
   "metadata": {},
   "source": [
    "### 2.1 Crystal objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a47530f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ce3 = pyddt.Crystal('ce3.in')\n",
    "ce4 = pyddt.Crystal('ce4.in')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdc03d81",
   "metadata": {},
   "source": [
    "Check if the created objects have the expected valence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "00847a9c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Ce3+', 'Ce3+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+',\n",
       "       'Fe2+', 'Fe2+', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-',\n",
       "       'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-',\n",
       "       'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-'],\n",
       "      dtype='<U32')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ce3.structure.atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "1e6239f8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Ce4+', 'Ce4+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+', 'Fe2+',\n",
       "       'Fe2+', 'Fe2+', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-',\n",
       "       'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-',\n",
       "       'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-', 'P1-'],\n",
       "      dtype='<U32')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ce4.structure.atoms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d9ca37f",
   "metadata": {},
   "source": [
    "### 2.2 Phase comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68284e8c",
   "metadata": {},
   "source": [
    "For now, let's consider X-rays of 7100 eV and perform the comparison of structure factor phases. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e4dc834e",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "E = 7100  # eV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ca4a935",
   "metadata": {},
   "outputs": [],
   "source": [
    "pyddt.phase(ce3.diffraction(E), ce4.diffraction(E), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d62a9f4e",
   "metadata": {},
   "source": [
    "![pyddt_t2_1](https://user-images.githubusercontent.com/106104347/187934482-f32cb04b-e819-4296-94dc-8f1820d6e054.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c076fdce",
   "metadata": {},
   "source": [
    "The result is an interactive plot with the hover text presenting the Miller indices. All plotted information was saved in the current folder (filename displayed on the screen).\n",
    "\n",
    "The vertical axis presents the phase difference, so 002 is the best reflection for phase measurements."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f5ea2f5",
   "metadata": {},
   "source": [
    "### 2.2 Optimum resonance condition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3513f9db",
   "metadata": {},
   "source": [
    "In the previous tutorial, we calculated the Fe resonance amplitude near its absorption edge, and characteristic \"jumps\" are noted. How does it influence the structure factor phase? \n",
    "\n",
    "To check this, let's calculate the phase of 002 as a function of energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "f7eec856",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdQAAAFFCAYAAABVHRyKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABGYklEQVR4nO3dd3yV9d3/8dcne5BJBhAIEDayDSIigrhHa7XVulcddba1tb21vVv702q13t6tq632brG4W7V1ggNQqiiyp8gegYRAIJC9vr8/rhM4hAAHOCcn4/18PK7HOeda5/PlCnnnuq7vdV3mnENERESOTUS4CxAREWkPFKgiIiJBoEAVEREJAgWqiIhIEChQRUREgkCBKiIiEgRR4S6gtcrIyHC9evUKdxkiItKKzJs3b7tzLrO5aQrUg+jVqxdz584NdxkiItKKmNmGg03TIV8REZEgUKCKiIgEgQJVREQkCBSoIiIiQaBAFRERCQL18hURaUMaGhrYvHkz5eXl4S6l3YmOjiYrK4vk5OSjWl6BKiLShmzfvh0zY8CAAURE6CBjsDjnqKyspKCgAOCoQlVbQ0SkDdm1axfZ2dkK0yAzMxISEsjJyWHbtm1HtQ5tERGRNqS+vp7o6Ohwl9FuxcfHU1tbe1TL6pBvKG36EnashpTu3pCcA1Ex4a5KRNo4Mwt3Ce3WsfzbKlBDack/YM6f/UYYdMreF7Ap3SGlx/7vE9JB/1lERNocBWoIzev/I4oSLiA3agdd2E5qzTai9myG0s1QtAy+ngZ1lfsvFBV/iMD17eVGx4WnQSIiclAK1BB6c9kOnpu9C4gEsjHLpkvyaHJS4+meGU/3fvH07lRN7+gScmwHneu2EbWnwAvc0s2w6gMoKzxwxYlZ+wI2rRek9fReU3tBag+Iim3JZoqIBKSoqIgHH3yQt99+m82bN5ORkcGwYcO44447OPfccw+7/PLly7nttttYvnw5paWldOvWjUsvvZT77ruPmJjwn05ToIbQL84fzA3j89i8s5LNOyt8r5UU7Kpg7oadvLV4K/UNzjd3PNCTrKT+dE+Lp3taAt17xNMjOZK8uN10j9hBZn0xMWVboHSTF7jbVnh7ufXVft9qkNzNF7A99w/ctF7eIWcdUhaRFrZ+/XrGjRtHUlISDz30EMOHD6ehoYGPPvqI73//+2zcuPGw64iJieGaa65h5MiRpKamsmjRIm688Ubq6up45JFHml1m8uTJTJ48mZkzZwa5RQdSoIZQdGQEPdIT6JGeAHQ+YHpdfQOFu6so8AVtY/AW7Kpk4aZdvLtkK3V7AxcglYxO2XRPO8kL3b4J9EyPo19CGT0jiulcs5WIXRtg1wbYuR7WzoA9W/f/0oTO0GWobxjuvWb0g4jIUP5TiEgHd+uttwIwd+5cOnXqtHf8oEGDuPLKKwEoLS3l7rvv5l//+heVlZWMGjWK//mf/yE/Px+Avn370rdv373L9uzZk5kzZzJr1qwWbMnBKVDDKCoywtsTTUtgTDPT6xscRburKNjl28MtadzDrWRpQSnTlhVSW78vcGMi0+ie3o1enU8nNyOBngMS6J0aSV50CV0bConetR6KlkLhYvjiz1Bf4yskHrIHQ5dhXsB2HQ5ZgyEmoUX+HUSkfSspKWHq1Kk88MAD+4Vpo9TUVJxznHfeeaSkpPD222+Tnp7Oc889x6RJk1i5ciVdu3Y9YLnVq1czdepUvvnNb7ZEMw5LgdqKRUYY3VLj6ZYaz+he6QdMr29wbNlVycaSCjbsqGDDjnLvtaSCL9buoLymfu+8ZpCT2p/juuUztF8KQ09JYHhcEamlK72ALVwCy16HeX/zLRABnfv5AnbYvj3axAP3tEUkfH791jKWb9ndot85uFsyv/rGcQHPv3r1apxzDBo06KDzzJgxg4ULF1JcXEx8fDwA999/P2+99RZTpkzhpz/96d55TzrpJObPn091dTU33ngjDz744NE3JogUqG1YZITtPaQ8ru/+05xzbC+rYWOJL2R3VLC6uIxlBaVMW1a0d75uKVkMyfk2Q3tcz5ATkxmetJv03V95AVu4BDZ+Dkv/uW/FKT2g2wjoNtIbuo7wLvURETkI59xh55k3bx4VFRVkZmbuN76qqoo1a9bsN+6VV15hz549LFq0iLvvvpuHH36Ye+65B4BZs2Zxzjnn7J23rq6O2tra/faM7733Xu69995jaVKzFKjtlJmRmRRLZlIsx/fcP/BKK2tZtqWUpQWlLC3YzdKCUt5fvi9ku6YkMCTnTIbmXMzQkSkMTa8no2wlbF0MWxfClgWw4q19K0zrtS9gu430DhnHpbRMQ0U6uCPZUwyXfv36YWasWLGCCy+8sNl5GhoayM7ObvZ8aNP76vbo0QOAwYMHU19fzw033MDdd99NVFQU+fn5LFy4cO+8r7/+Oq+99hovvPDC3nHp6aHZCVCgdkAp8dGc1CeDk/pk7B23p6qWZVu8cF3iGz5cUUTjH5bZybEMzTmJ4d3PYcTwVIZlQMrOZV64blkABfNg2Rv7vqRz3/1DtsswiD3w3ImItH/p6emcddZZPPnkk9x5550HnEfdtWsXo0aNoqioiIiICPLy8gJed0NDA3V1ddTX1xMVFUV8fPx+HZeysrIOGBcqClQBICkumhPzOnNi3r5zpGXVdSzzhevSglIW+0K2UV5mIiO6n8qI3AsZPjaVQSm1xGxb5AvZhbDhM+9uUeCdk83oD7ljof/Z0Pc0iNT9SEU6iqeeeopx48aRn5/P/fffz7Bhw3DOMWPGDB566CE2bNjAuHHjuOCCC3jkkUcYOHAghYWFTJ06ldNPP53x48czZcoU4uLiGDp0KDExMcydO5d77rmH73znO8TGhv/6ewWqHFSn2CjG5HVmjF/I7q6qZfGmUhZt3sWCjbv4ZNV2Xl/gPe4oJjKCwd2SGdHjXEYMuJwRp6fSM3YPttUXsgXzvYCd9zfvetgRl8PIq6Bzn3A1UURaSF5eHvPnz+fBBx/kZz/7GQUFBXTu3Jnhw4fzzDPPYGa8++67/OIXv+DGG29k27ZtZGdnM27cOK6++moAoqKieOihh1i1ahXOOXr27Mltt93Gj370ozC3zmOBnCzuiPLz893cuXPDXUar55xjS2kVizbtYqFvWLK5lMpar4dxakI0w7unMrxHKiN6pHB8TidStnwC856DVdPANUCv8TDqGhj0Dd1WUeQwVqxYccjesnLsDvVvbGbznHP5zU3THqocEzMjJzWenNR4zh3qXSdWV9/Aqm1lLNy0a2/QPjl9FQ3O65k8pnc6l57wGGef3UDM0pdh/t/h9RsgLhWGX+qFa/bg8DZMROQIKVAl6KIiIxjUNZlBXZO57IRcAMqr61hSUMp/Vm3n34sKuPOlBWR0iuHi/PO5/Krv06N0rrfXOvev8MWfoPtoL1iHfFs3mBCRNkGHfA9Ch3xDp6HB8cmqYl74YiMfrSjCARP6Z3LlmJ6cmhtJ5JJXvHDdvtLbax15JeRfr3OtIuiQb0vQIV9pMyIijIkDspg4IIstuyp5+ctNvDxnIzf8fS7dUuK47IQz+e7V15NVMh++fNbbY539JPQ5DU64EfqdqXsPi0iroz3Ug9AeasuqrW/goxVFPP/5Rv6zejtREcaZx2VzzdhenJBRjc2f4vUO3rMVUnIh/zoYdTUkZhx+5SLtiPZQQ+9o91AjQlrVMTKzU8zsTTMrMDNnZtc2mW5mdp+ZbTGzSjObaWbHNZknzcymmFmpb5hiZqkt2Q45vOjICM4e0pXnbxjDjJ9M5Lpxvfh09Q6++8znnPe3Nbza6XKqblsIl/zdexzdR7+GxwbB6zfDpi9BfxiKSJi16kAFOgFLgR8Alc1M/ynwY+AOYDSwDfjAzJL85nkRGAWc7RtGAVNCWLMco94Zifz8vMF8fs9p/PaiodQ3OH76z8Wc9LtZ/G7TALZe+A+49Qs4/lr46h34v9Ph2Umw+B9QVxPu8kWkg2ozh3zNrAy43Tk32ffZgC3Ak8653/jGxeOF6k+cc382s0HAcuBk59ynvnlOBmYBA51zKw/2fTrk23o455i9dgeTP13PhyuKMDPOHtKF607qxfFdorDFr3iPo9uxCpK6wugb4Pjr9GQcaZd0yDf0OmKnpN5AF+D9xhHOuUoz+wQ4CfgzMBYoAz7zW+5ToNw3z0EDVVoPM9t77+FNJRVM+XwDL8/ZyDuLtzIkJ5lrTzqH82++lrgNM+Hzp2H6/fDJ72DYd+HEWyBLv3xEJPRa+yHfQ+niey1qMr7Ib1oXoNj57Yb73m/zm2cvM7vJzOaa2dzi4uIQlCzHqkd6AveeO4jP7z2N31w4hOraBn7yj0WMe3gmj63LpfhbL8Otn3s3iFj8Cjx9Ivz9Avh6GjQ0hLt8EWnH2nKgBp1z7hnnXL5zLr/pM/mkdUmIieKKMT15/0en8MINYxiZm8oTM1Yz7rfT+eknNawc/QDctQJO+yUUfw0vXgJP5sOcZ6G6LNzli3RIRUVF/OAHP6BPnz7ExsaSk5PDOeecw7vvvnvE61q1ahVJSUkHPLkmnNpyoBb6XrObjM/2m1YIZPrOtwJ7z71m+c0jbZiZMa5vBn+5ZjQf3TWBS0Z3581FWzjr959w9cur+ST7atwPFsG3/w/iU+Hdn8Bjg+Gj/wcVJeEuX6TDWL9+PaNGjWLatGk89NBDLF68mA8//JDzzjuP73//+0e0rpqaGi699FJOOeWUw847efJkJk6ceJRVH5m2HKjr8ELxjMYRZhYHjGffOdPZeD2Fx/otNxZIZP/zqtIO5GV24oFvDWX2f53G3WcNYMXW3Vz91zmc/cTnvFo9hurrPoDvfQh9JsKsx+D3w7w7MolIyN16660AzJ07l0suuYQBAwYwaNAgbr/9dhYvXgxAaWkpN910E1lZWSQlJTFhwgSa6xz6s5/9jGHDhnHxxRe3aBsOp1UHqpl1MrMRZjYCr9Zc3+dc37nQ3wM/M7OLzGwIMBmvE9KLAM65FcBU4M9mNtbMxuJ1Vnr7UD18pW1LS4zhtlP78p+fncqjFw/HDH76z8WM++0MHv86lZLz/gK3fAY5o+CtO+GDX+r8qkgIlZSUMHXqVG677bZmD9GmpqbinOO8886joKCAt99+mwULFnDKKacwadIktm7dunfed955h7fffpsnnniiJZsQkNbeyzcfmOH3+de+4TngWuARIB54CkgDvgDOdM7t8VvmcuAJYJrv85vA7SGtWlqF2KhIvnN8d749KofP1uzg2VlreeyDr3lqxmq+fXx3rj/r7/Sd+2v49A+weytc+Cfd0lDanvf+CwqXtOx3dhkK5/w24NlXr16Nc+6Ql/vMmDGDhQsXUlxcTHx8PAD3338/b731FlOmTOGnP/0pW7Zs4cYbb+SNN95oVedOG7XqQHXOzQTsENMdcJ9vONg8O4Erg1yatCGN51nH9c1gVdEe/vrpOv45bzMvz9nIwxf9iIsndYPpD0BEFFzwFES06gM3Im1OIPc7mDdvHhUVFTTtEFpVVcWaNWsAuOqqq7jlllsYM2bMQdcza9YszjnnnL2f6+rqqK2t3S+A7733Xu69994jbcZhtepAFQm2ftlJPHTRMH585gB+9MpC7n5tCRXfvJhrJjbAzAchNgnOfSTcZYoE7gj2FMOlX79+mBkrVqzgwgsvbHaehoYGsrOzmTVr1gHTkpOTAZg+fToff/wxv/71rwEvqBsaGoiKiuLpp5/mpptuIj8/n4ULF+5d9vXXX+e1117jhRde2DsuPT09iK3bR4EqHVJGp1j+ck0+t7+4gF+9uYzE71zOd04shc+fgqyB3uPiRCQo0tPTOeuss3jyySe58847Dzhcu2vXLkaNGkVRURERERHk5eU1u54lS/Y/tP3vf/+b3/zmN8yZM4ecnBwA4uPj6du37955srKyDhgXKjq2JR1WbFQkT14+kpP7ZvCz15cwPfd26HsGvHs3bJgd7vJE2pWnnnoK5xz5+fn84x//YOXKlXz11Vf88Y9/ZNiwYZx++umMGzeOCy64gPfee49169Yxe/ZsfvWrX+3dax0yZMh+Q05ODhEREQwZMoS0tLQwt1CBKh1cbFQkf7rqeAZ3Teb2lxezdsIfIKUHvH4TVJWGuzyRdiMvL4/58+dzxhln7L3sZdKkSbz55ps888wzmBnvvvsukyZN4sYbb2TAgAFccsklrFy5km7duoW7/IC0mZvjtzTdHL9jKSyt4vwnZpESH81bF8aSMOU8GHoxXPTncJcmsh/dHD/02uXzUEVaSpeUOB6/bCTrtpfz8y/jYfyPYfHLsGbG4RcWEUGBKrLXSX0yuH1SP95YUMD0zKsgrbd3PlXPWBWRAChQRfzcfmpfBnZJ4r/e/JrySb/xnrH6xZ/CXZaItAEKVBE/MVERPHrxcHaU1/Dw2p7Q93T4z/9C9Z7DLywiHZoCVaSJITkpXDEmlxe+2MjG4T+EyhL4XHup0nqoM2noHMu/rQJVpBk/PL0/CTGR/HJuLPQ/B2Y/ob1UaRUiIyOpra0NdxntVmVlJdHR0Ue1rAJVpBnpiTHcMakvM1cWs7zvjd41qQtfCndZIqSmplJUVESDnpAUVM45KioqKCgoICsr66jWoetQD0LXoUpFTR3jH57BkJwUnqu/Byp3wu1zdfN8CauGhgY2b95MeXl5uEtpd6Kjo8nKytp77+DmHOo6VN3LV+QgEmKi+N743jwydSUbzrmGnjPugNUfQP+zwl2adGARERHk5uaGuwxphv7UFjmEq07sSXJcFL/bNAASM2HBlHCXJCKtlAJV5BCS4qK59IRc3lu+g/IBF8HKqVBREu6yRKQVUqCKHMYVY3Kpb3C84SZAQy0s+We4SxKRVkiBKnIYPTsnMqF/Jk8si8VlD4VF6u0rIgdSoIoE4MoTe1K0u5pVWWfBlvmwa1O4SxKRVkaBKhKAiQMySU+M4cU9w70RX70T3oJEpNVRoIoEIDoygvOHdeWl1dHUZwyEFW+FuyQRaWUUqCIBumBEN6rrGliZNhE2fgbl28Ndkoi0IgpUkQCNyk2jR3o8L+4eCq4BVn8Y7pJEpBVRoIoEyMw4b2g3XtmUSkNCBqyZEe6SRKQVUaCKHIEzBmdR22BsSR8Da6aD7oUtIj4KVJEjMKJHGp0TY5hZNwTKt0HRsnCXJCKthAJV5AhERhiTBmbxt8Le3og108NbkIi0GgpUkSN0+uBs1lQlU5HSD9Z9HO5yRKSVUKCKHKHx/TKIiYzgq9ghsGkONNSHuyQRaQUUqCJHKCEmilE9U/moPA+qd8O25eEuSURaAQWqyFE4uW8G/9rhe8jzxs/DW4yItAoKVJGjcHK/TArIoDK+C2ycHe5yRKQVaNOBamb3mZlrMhT6TTffPFvMrNLMZprZceGsWdqHoTkpJMVF83XMcbBhtq5HFZG2Hag+K4GufsNQv2k/BX4M3AGMBrYBH5hZUksXKe1LZIRxUp/OfFSRB3u2wK6N4S5JRMKsPQRqnXOu0G8oBm/vFPgh8Fvn3GvOuaXANUAScHn4ypX24qQ+GXxY5rsedfOX4S1GRMKuPQRqnu+Q7joze9nM8nzjewNdgPcbZ3TOVQKfACeFoU5pZ0b3Sudr1536iFjYsiDc5YhImLX1QP0CuBY4G7gRL0A/M7POvvcARU2WKfKbth8zu8nM5prZ3OLi4tBULO3GgC5JxMfGURDXV4EqIm07UJ1z7znnXnXOLXbOfQicj9ema45yfc845/Kdc/mZmZlBrVXan8gIY1TPNObX9YItC3WDB5EOrk0HalPOuTJgGdAPaOztm91ktmy/aSLHZHSvND4p6wG15bB9VbjLEZEwaleBamZxwEBgK7AOLzjPaDJ9PPBZWAqUdmd0r3QWO99pex32FenQ2nSgmtmjZjbBzHqb2Rjgn0Ai8JxzzgG/B35mZheZ2RBgMlAGvBiumqV9Gd4jlU0ROdRExMOW+eEuR0TCKCrcBRyj7sBLQAZQDHwOnOic2+Cb/ggQDzwFpOF1YjrTObcnDLVKOxQXHclxOWms2tmH47SHKtKhtelAdc5depjpDrjPN4iExKjcNL4o7Mngwo+w+lqIjA53SSISBm36kK9IazC8RyqL6npidVXqmCTSgSlQRY7RiB6prHA9vQ+FS8JbjIiEjQJV5Bh1T4unNKEntRYDRQpUkY5KgSpyjMyMobkZrLMeULg03OWISJgoUEWCYHj3VBbU9KChcIke5SbSQSlQRYJgRG4qK1wuERXboazp7aNFpCNQoIoEwbDuqaxoaOyYpMO+Ih2RAlUkCFLio6noPND7ULg4vMWISFgoUEWCpF9ud7aQiSvSHqpIR6RAFQmSkT1SWVqfS90W7aGKdEQKVJEgGdEjjRWuJ1Ela6C2MtzliEgLU6CKBMmALkmssl4YDbBtebjLEZEWFtDN8c0sD5gA9MJ7eksxMB/41DlXFbLqRNqQmKgI6rOOgx14tyDMOT7cJYlICzpkoJrZFcAPgHygCNgCVALpwANAlZm9ADzs98g0kQ6ra8+BlG2PJ37rEiLDXYyItKiDHvI1swXAD4HngJ7Oua7OueOdcyc75wYDycAFvnXMNbOLW6JgkdZsRM90lrtcqjYtDHcpItLCDnUO9efOudHOuaecc5uaTnTOVTvnZjrnvg8MBtaHqkiRtmJkj1SWN/QkZvtyaGgIdzki0oIOGqjOuXcDXYlzrtg592VwShJpu7qnxbMhug/R9RWwc124yxGRFhRop6Tcg0xyQJVzrjh4JYm0XWYGXYbBVmDrIujcJ9wliUgLCfSymfXAumaG9UChme00s8fMLKCAFmnPMvKGU+siqd68MNyliEgLCjRQLwM2A78AzvANvwA2AtcD9wFXAf8d/BJF2pahPbNY5bpTsXFBuEsRkRYU6B7lLcCPnHOv+42bbmYrgR845yaY2Tbg18Cvgl2kSFsyvEcq7zf0pFfxsnCXIiItKNA91DHAkmbGLwVG+97PBroHoyiRtiwlPpqixP4k1O6APYXhLkdEWkiggboBuKmZ8TfiHfYFyARKglGUSJvXZRgAbuuiMBciIi0l0EO+PwZeM7NzgcbLY/KBPsC3fZ9HA68GtzyRtimj7/GwAXavm09K/7PCXY6ItICAAtU5946Z9cc7lzrAN/pN4E/OuY2+eZ4OTYkibc/g3t1Z35BNzMYFpIS7GBFpEQFf5uILzntCWItIuzGwSzLT6ckJ2/WwcZGOIuDHt5nZUDN70szeNbOuvnHfMrORoStPpG2KiYpgR9JA0qsLoGp3uMsRkRYQUKCa2Zl4505zgNPwHuEG3jlUXSYj0ozIbsMBqC5YHOZKRKQlBLqHej9wl3PuQqDGb/xM4IRgFyXSHnQdOAaAwhWfhbkSEWkJgQbqEKC5m+WX4D0bVUSaGDFoIAUug+r1c8Jdioi0gEADtQTvcG9To/BuSSgiTaQkRLM6ZiBpO3XIV6QjCDRQXwR+Z2bd8Z4wE2VmE4BHgb+HqrhgMrNbzWydmVWZ2TwzGx/umqT9q8waSWZ9EbW7toS7FBEJsUAD9Rd4T5fZAHQClgPTgf8AvwlNacFjZt8F/gA8CIwEPgPeO8Rj6USCIrnfSQBsXPJJmCsRkVALKFCdc7XOuSuA/sAlwOXAQOfcVc65+lAWGCR3AZOdc88651Y45+7Ae2LlLWGuS9q5fsPHUeMi2b1KHZNE2rsjen6pc24NsCZEtYSEmcUAx+Mdnvb3PnBSy1ckHUlmWgorIvOIL5of7lJEJMQOGqhm9tdAV+Kcuz445YREBhAJFDUZXwSc3vLlSEezM304I4vfpLqmmtiY2HCXIyIhcqhDvplNhm8DFwJ9fcO3gIvwAqtdMLObzGyumc0tLi4OdznSTiT1PYl4q+GrBZ+GuxQRCaGDBqpz7huNA14nnmlAd+fcKc65U4AewFTgi5Yp9ahtB+qB7Cbjs4H9HlbpnHvGOZfvnMvPzMxsqfqknesz2nvazPalH4W5EhEJpUB7+d4J3OecK28c4Xt/P3BHKAoLFudcDTAPOKPJpDPw/lAQCamEzt3ZHJVL0tbZ4S5FREIo0EDtBHRrZnxXICF45YTMY8C1ZnaDmQ0ysz/gtedPYa5LOojS7DEcV7uULTt0o3yR9irQQH0N+JuZXWpmvXzDpcD/Aa+HrrzgcM69AvwQ73rahcDJwLnOuQ1hLEs6kNTjTifRqlk2d2a4SxGREAk0UG8B3gIm4102swZ4DngHuDUklQWZc+5p51wv51ysc+5455yutJcW022416G87KvpYa5EREIloOtQnXOVwK1mdjfeI9sA1vifUxWRg7PEDArj+9F1xxfsqaolKS463CWJSJAF/IBx8DoiOecW+waFqcgRcH1P53hbySeLV4e7FBEJgYMGqpn9xcx6Hm4F5rnCzK4Kbmki7Ut2/reItnq2zH0r3KWISAgc6pDvJmCxmX2Bd/50LrAFqALSgMF4nXu+C6wHbg5ppSJtXESP0ZRHpdKlcCbl1XeRGHtEd/4UkVbuUDd2+DXeHZFmAjfhPVlmLV6oLgP+CuQC1znnxjnnloa8WpG2LCKSyp6nMd4WMnWxHiMs0t4c8hyqc67YOfegc24o3i0GRwHjgAFAmnPuO86591ugTpF2ofPxF5Bq5Sz/7J1wlyIiQRZwpyTn3E7n3CLn3OfOudXOORfKwkTaI+t3JjWRifQvfp+1xWXhLkdEguiIevmKyDGKjqd+wPmcEzmH179sU09CFJHDUKCKtLD4Ud8l2SoonPsWVbX14S5HRIJEgSrS0npPoCYug7Nqp/OvBQXhrkZEgkSBKtLSIqOIzr+GSZELePOTOTQ0qDuCSHtwRIFqZhlmNsbMYkNVkEhHYPnXEgGM3fUWM1ZuC3c5IhIEAQWqmSWZ2avANrxniOb4xv/JzO4LXXki7VRqLq7/WVwRNYPHpy3VXqpIOxDoHurDeCE6Cqj0G/82cGGwixLpCCLG3ko6pQwtfot3l24NdzkicowCDdRvAj90zi0E/P+UXgHkBbsokQ6h13hc9zHcEfM2j767hIJdlYdfRkRarUADNQ3Y0cz4JED9/kWOhhk24W6yXTETKj/i/Mdn6XyqSBsWaKB+ibeX2qhxL/VmvHOqInI0+p4OOfn8d8LrDO20m+v+9iWPTlupc6oibVCggXovcL+ZPYv3hJq7zGw6cBXwi1AVJ9LumcE3fk9UQw3P1d/LXcdV8OSM1dzx8gLd9EGkjQkoUJ1znwEnATHAGuA0vKfOjHXOzQ9deSIdQJeh8L1pWFQMd2y8k6fG7OSdxVu5+q9z2F1VG+7qRCRAR3Jz/CXOuWucc0Occ4Odc1c655aEsjiRDiNrENzwIZbeh/OW3sVLkypYsHEn1/x1DmXVdeGuTkQCEOh1qIPNbIDf5zPM7Hkzu8fMIkNXnkgHktQFrnkTMvozds7tvHBaNYs3l3L9376kokahKtLaBbqH+ldgJICZ9QD+DaQDtwEPhKY0kQ4oIR2u/jek53HCF3fyf+ckMHdDCT94eaE6Kom0coEG6kCg8Vzpd4AvnHPn4nVKuiwUhYl0WImd4fJXITqeiXPv4KEzsvhgeREPT/sq3JWJyCEEGqiRQI3v/WnAu773a4DsYBcl0uGl9oDLX4GK7Vyy/r+5ZkwOf/54rZ5OI9KKBRqoS4FbzGw8XqBO9Y3PAbaHojCRDq/bSDj/f7GNs/lV8tuc0Cudn7+xhA07ysNdmYg0I9BA/RlwIzATeMmvd+83gTkhqEtEAIZfCiOuIGLWozw1royoyAjufGkBNXUN4a5MRJoI9DrUT4BMIMM5d73fpD8Dt4SiMBHxOfd30LkPmR/9mEcv6MOizaU8PXN1uKsSkSaO5DrUeufczibj1jvndPNRkVCKSYRvPA67NnBG0V/5xvBuPD1jDau3lYW7MhHxE3CgmtmpZvaMmU01s+n+QygLFBGg1zg4/jr4/Gn+3+g64qIjuPeNJTinS2lEWotAb+xwLfAe3tNlJgLFeE+gGQUsD1FtIuLvjF9DfDpps+7jv84eyJx1JbyzRM9RFWktAt1D/Qlwu3PuMqAWuMc5NxJ4HtBxJ5GWEJcCp94DG/7Dd1OWMLBLEo9MXakOSiKtRKCBmgd86HtfDXTyvX8SuDbINYnIwYy6FjIGEPnBL7nnrD5sLKng+c83hLsqESHwQN2Bd7gXoAAY4nvfGYgPdlGBMrOZZuaaDC83mSfNzKaYWalvmGJmqWEqWeTYREbB6fdByRpOqZrJuL6deXLGat3rV6QVCDRQZwFn+t6/CjxuZn8DXgI+CEVhR+BvQFe/4eYm01/EO9d7tm8YBUxpyQJFgmrAOZA9FPvPY9x1Wh9Kymt4ec6mcFcl0uEFGqi344UnwEPA7/D2Tl8FbghBXUeiwjlX6DeUNk4ws0F4IXqTc262c242XuCe7//0HJE2xQxO+QnsWM3xZR9zQu90nvlkrc6lioRZoDd2KHHObfG9b3DOPeyc+6Zz7ifOuV0hrfDwLjWz7Wa2zMweNbMkv2lj8TpNfeY37lOgHO+B6SJt06BvQuZAmPUYt03sQ+HuKt5YsDncVYl0aAFfhwpgZt3MbISZjfIfQlVcAF4ErgBOBe4Hvg285je9C1Ds/C7W873f5pu2HzO7yczmmtnc4uLikBYuckwiImDsbbBtGafErOS4bsn833/W6bpUkTAK9DrUkWa2DNiE9xi3uX7Dl8EsyMweaKajUdNhIoBz7hnn3DTn3BLn3MvAd4EzjjbkfevLd87lZ2ZmBq9RIqEw9GKIT8O+fJZrxvbi66Iy5qwrCXdVIh1WoHuoz+CF6Xi8S2h6+w15Qa7p98CgwwwHuyH/XKAe6Of7XAhkmpk1zuB7n+WbJtJ2RcfDyKtgxdt8s3cDyXFRTNElNCJhExXgfIOBkc65r0NZDIBzbjtH/0i4oXjPbm28fcxsvGtmx7LvPOpYIJH9z6uKtE2jvwefPUHc4ilcnH8Rz322nm17qshKigt3ZSIdTqB7qEto5pxjOJlZHzP7pZnlm1kvMzsXeBlYgNfxCOfcCrxnt/7ZzMaa2Vi8J+S87ZxbGbbiRYIlrRf0PR0WvsQVJ3SnrsHxj7nqnCQSDgcNVDNLbxyAe4FHzOx0M8v2n+abHg41eA87nwasBB4H3gdOd87V+813ObDIN9803/urWrZUkRAacRns3kxe2XxG90rj9fmb1TlJJAwOdch3O+D/v9LwAqvpOId3mLVFOec2ARMCmG8ncGXoKxIJkwHnQWwKLHyJC0f+nHvfWMKSglKGdU8Nd2UiHcqhAvXUFqtCRI5edBwMuRAWv8p5pz7EfZERvD6/QIEq0sIOGqjOuY9bshAROQbDL4N5k0nZMJXTBg3grUVb+Pl5g4iOPKJLzUXkGBzyf5uZJZjZU2ZWYGbbzOxFM8toqeJEJEA9xkByd1j+by4cmcOO8ho+XX20neVF5Ggc7s/XX+M9nu0dvB60ZwB/DHFNInKkzGDwBbBmOqfkxpIQE8m0ZUXhrkqkQzlcoF4EfM85d5Nz7k7gPOBbZtbinZBE5DAGXwD1NcSt+4BTB2bxwfJC6hvU21ekpRwuUHvgPboNAOfcHKAO6BbKokTkKHQfDUndYPm/Ofu4Lmwvq2Hehp3hrkqkwzhcoEbiXe/pr47A77AkIi0lIgIGfxNWfcCpveOJiYxg6lLdYVOkpRwuUA143szebByAOODZJuNEpDUYeD7UV9Np8yzG98tg2rJC3eRBpIUcLlCfA7YAO/yG5/FulO8/TkRag9wTvZs8rJrGWUO6ULCrkmVbdoe7KpEO4ZCHbp1z17VUISISBJHR0OdU77DvqY8CMHPlNobkpIS5MJH2T1d9i7Q3/c+CsiIyy1YyNCeFmSuLw12RSIegQBVpb/qeARisep+JAzKZv3Enuyqa9i0UkWBToIq0N50yIWcUfD2NiQOyaHAwa5XumiQSagpUkfao31lQMI8R6XWkJkTrsK9IC1CgirRH/c8EHJFrpzO+XyYff11Mg+6aJBJSClSR9qjLcIhPg7UzOXVAJtvLqnX5jEiIKVBF2qOICOg9AdbOZHxf7wFRs1brsK9IKClQRdqrvImwZwuZ1RsYkJ3EZ6t1DxaRUFKgirRXfU71XtfOZFzfDL5cX0JVbX14axJpxxSoIu1VWi9vWDuTcX07U13XwHw9fUYkZBSoIu1Z3kRYN4sTcpOIjDA+XaPrUUVCRYEq0p7lnQo1e0jasYQRPVL5j86jioSMAlWkPet9CmCwdgbj+nRmyeZdlFbWhrsqkXZJgSrSniWkQ9fhsHYmJ/XNoMHBF2u1lyoSCgpUkfYubyJs/pKR2ZHERUfw6WqdRxUJBQWqSHuXNxEa6ogt+IITenfm0zXaQxUJBQWqSHuXeyJExnqXz/TpzOptZRTtrgp3VSLtjgJVpL2LjvdCde3HjPPdhvA/epybSNApUEU6gryJsG0Zg5MqSUuI1vWoIiGgQBXpCPImAhCxfhYn9c3gs9U7cE6PcxMJJgWqSEfQdTjEpcDamZzcN4PC3VWsKS4Pd1Ui7YoCVaQjiIj0bvKwdiYn9+kMoMtnRIKs1Qaqmd1kZjPMbJeZOTPr1cw8aWY2xcxKfcMUM0ttMs9QM/vYzCrNrMDMfmlm1lLtEGk18ibC7s30oJDc9AT+o0AVCapWG6hAAvA+cN8h5nkRGAWc7RtGAVMaJ5pZMvABUASMBn4A3A3cFZKKRVqzvMbHuc1gXN8MPl+zg7r6hvDWJNKOtNpAdc793jn3EPCf5qab2SC8EL3JOTfbOTcbuBk438wG+Ga7Ai+Yr3HOLXXO/RN4GLhLe6nS4aTnQUqPvedR91TXsbigNNxVibQbrTZQAzAWKAM+8xv3KVAOnOQ3zyznXKXfPNOAbkCvFqhRpPUwg94TYN0sxvZOxQw+1fWoIkHTlgO1C1Ds/Pr++95v801rnKeoyXJFftP24ztvO9fM5hYXF4egZJEwy5sIVbtI372C47ol6zyqSBC1aKCa2QO+DkaHGia2ZE3+nHPPOOfynXP5mZmZ4SpDJHTyJniva2cyrm8G8zfupKKmLrw1ibQTLb2H+ntg0GGGOQGuqxDI9D8X6nuf5ZvWOE92k+Wy/aaJdCydsiDruL3nUWvrHXPWlYS7KpF2oUUD1Tm33Tn31WGGigBXNxvohHeetNFYIJF951VnA+PNLM5vnjOALcD6Y2uNSBuVNxE2fs7onHhiovQ4N5FgabXnUM2si5mNAPr7Rg02sxFmlg7gnFsBTAX+bGZjzWws8GfgbefcSt8yLwIVwGQzG2JmFwH/BTzmdN816ajyJkB9NXGFc8nvmcanq/U4N5FgaLWBCnwfWAC84Pv8ju/zN/3muRxYhNdzd5rv/VWNE51zpXh7pN2AucBTwP8Aj4W4dpHWq+dJEBG19zzq8q272VFWHe6qRNq8Vhuozrn7nHPWzDDZb56dzrkrnXPJvuFK59yuJutZ4pw7xTkX55zr6pz7tfZOpUOLTYLuo2H1R4zv5z3O7ZNV6tUucqxabaCKSAj1OxMKFzOkUzkZnWKZ/pUCVeRYKVBFOqL+ZwMQseYDJg3M5OOV26jVbQhFjokCVaQjyhoEqbmwciqTBmazu6qOeRt2hrsqkTZNgSrSEZl5e6lrZ3Jyr0RiIiOY/tW2cFcl0qYpUEU6qv5nQV0lnbbMZkxeOh+taHqXThE5EgpUkY6q13iIToSv3+O0gVmsKS5n/fbycFcl0mYpUEU6qqhY6HMqrJzKaQO9e1d/sFx7qSJHS4Eq0pENvgD2bKFH+VKG5CTzzpKt4a5IpM1SoIp0ZAPOgchYWPYG5w7tysJNuyjYVXn45UTkAApUkY4sNgn6nQHL/sW5g7MAeE97qSJHRYEq0tEddyGUFdKrYgmDuybz1mIFqsjRUKCKdHT9z4aoeFj2Ot85vjuLNu3i8Y9WoVteixwZBapIRxfbybsmddkbXH1CVy4cmcNjH3zNj19dRHVdfbirE2kzFKgiAiOugIodRK2exmOXDOfHZ/Tn9QUFXPmXLygprwl3dSJtggJVRKDvaZDUDRY8j5lxx2n9ePLykSzeXMq3nvqUFVt3h7tCkVZPgSoiEBEJIy6H1R/C7i0AnD+sGy/fdCJVtfVc+PSnvLFgc5iLFGndFKgi4hl5JbgGWPDCvlG5abx958kM757Kj15ZxC/+tUTnVUUOQoEqIp703tBnEnz5LOzZdwvCrKQ4XrhhDDdPyOP5zzdyyZ9ms2GH7vkr0pQCVUT2mXgPVJXCH8fCV+/sHR0VGcE95wziT1cez9rt5Zz7h1m8OneTLq0R8aNAFZF9epwAN38CyTnw8uXw79uhes/eyWcP6cLUH57C0O4p/PSfi7nl+fnsVC9gEUCBKiJNZQ6AGz6Ck++CBc/Dn06GdZ/snZyTGs+LN5zIPecM5KOvijjr95/w/rLCMBYs0jooUEXkQFExcPqv4Lp3vc/PfQP+fRtUlAAQEWHcPKEP/7ptHOmJMdw0ZR63PD+Pbburwli0SHgpUEXk4HqeBLfMhnE/hIUvwVMnwJJ/gu/c6XHdUnjrjpO5+6wBfPTVNk577GNe/GIjDQ06tyodj6lTQfPy8/Pd3Llzw12GSOtRuATevBO2zIde4+Hsh6DL0L2T1xaXcc/rS/hiXQnDe6Tyq28MZlRuWhgLFgk+M5vnnMtvbpr2UEUkMF2Gwg0fwrmPQtEy+NN4L2DLigHIy+zESzeeyKMXD2frrkouevozfvTKQgpLdRhYOgbtoR6E9lBFDqFyJ3z8CMx5BqIT4OQfwgk3ezfaB8qr63h65mqenbWOSDNuHN+bG07JIzkuOrx1ixyjQ+2hKlAPQoEqEoDtq+D9/4av34OEDBh/F+R/D6LjANhUUsFv3/uKd5ZsJSU+mpsn5HHtSb1IiIkKc+EiR0eBehQUqCJHYNOXMOMBWDsTkrrC+B97tzKMjgdgaUEp//P+SmasLCajUwy3TOzLZSf0ULBKm6NAPQoKVJGjsG4WTH8ANn0OiZkw5mYYfQPEe52T5m0o4dFpXzN77Q7SEqK5emwvrjmpF+mJMWEuXCQwCtSjoEAVOUrOwfr/wKd/gNUfQHQiHH8tnPh9SM0FYO76Ev708Ro+XLGNuOgILh2dy/XjepPbOSG8tYschgL1KChQRYKgcCl89rh37SoO+p3l7bH2mQQREXxdtIdnPlnLvxYUUO8cE/tnctXYnkzon0VkhIW7epEDtMlANbObgMuAkUAK0Ns5t77JPOuBnk0Wfdg5919+8+QCTwGTgErgReAnzrlD3oBUgSoSRLs2wbzJMP85KC+GtF6Qfz2MuAISMygsreLFORt5ac5GivdU0yM9nivG9OSiUTlkJcWFu3qRvdpqoP4QiMcLwf/l4IH6HPBHv9Flzrky3/RIYCGwA7gL6Oyb/zXn3B2H+n4FqkgI1NXAijdh7l9hw6cQEQX9zoThl0L/s6khmveXFzJl9ga+WFeCGQzNSWFUbhojc1MZlZtG97R4zLT3KuHRJgO1kZnlA19y8EB90jn36EGWPQd4B+jpnNvkG3cl8Bcgyzm3+2Dfq0AVCbFtK2Dhi7D4VSgrhLgUOO4iGPZd6DGGVcXlvLukkM/WbGfx5lIqa70Hm2cmxTIqN5WRuWmMyk1jWPcU4qIjw9wY6Sjae6DGAdHAJuAfwO8aD+ea2f8Dvu2cO85vmUxgGzDJOTfjYN+rQBVpIQ313uU2i16GFW9BXSV0yoYB58Kgb0Cvk6mzaL4q3MOCjTuZv3EX8zfuZMOOCgCiIoz+2UkMzUlhSE4yQ3JSGNQ1WSErIXGoQG3rF4E9DizAO6R7AvBboDdwg296F6CoyTLbgXrfNBEJt4hI6HuaN1Tvga+neYeFF78K8/4GEVFEZfRnSJehDMkewlX5Q+H8oWx3SSz0heuSglLeX17IK3M3ARAZYfTL6sRx3VIY6heyibFt/VeetGYtuodqZg8APz/MbKc652b6LXPQPdRm1n8J8AqQ4ZzbYWbPAH2dc5P85jGgFrjKOfdSk+VvAm4CyM3NPX7Dhg2BNk1Egq22EtZ+DJvneL2FC5fAni37pid1hewh3j2GuwzBZQ6kIKIbS4uqWballCUFpSwtKGV72b7+hz3S4xmQncyALp0Y0CWZgV2S6J2RSHSkbmsugWk1h3zNLAPIOMxsG51zFX7LHEmg9gTWAyc6577QIV+RdqZ8BxQt2RewRUuh+CtoqPOmWySk53kPSc8ciMscwI74PBZVZrK8uIavivbwdeEe1m4vp973iLnoSKNPZif6ZycxoEsS/bOTyMtMpEdaAjFRClrZX6s55Ouc2453yDVURvhet/peZwO/MLPuzrnNvnFnANXAvBDWISKhkNgZ8iZ6Q6O6atj+NRSv9MK1+Cvv/cr3MFdPBnAaxmlpvSBzIAzpS21qbwoicviqJpOFu+L5els58zbs5M1F+/aAIyOMHmnx5GV2ondGIr0zEsnLSKRHegKZSbE6RysHaLUnFMysC955zv6+UYPNLBVvD7bEzMYCJwIzgFJgNN7lNW865zb6lnkfWAb83cx+jHfZzO+AZw/Vw1dE2pCoWN9h36H7j6+rgZI1+wK28XXNdKLrq+kF9ALOjoqH9N7QK4/qlN4URnVjg+vCiupMFpcmsHZHBZ+t2U5VbcN+q09LiCY7OY6s5Diyk2LJTo4jOznW++x7n9EpVoeTO5BW28vXzO4DftXMpOucc5PNbBTwNDAQiAU2AC8DjzQ5ZJzrm6/xxg4vAHc756oP9f065CvSTjU0wO7NULIWdqzxe10DO9dDvd89X6LiISUHl9ydyoQubI/Iosg6s7UhnU01nVhflcCaini27qlj257qvYeRG5lB58RYspP3BW5mp1gykrywzegUS+dOMWR0iiU5LkrX17YBreYcaluiQBXpgBrqoXSzF647fAFbutkbdhfAnkKgmd+Zcam4xExq4zpTGZPGnsg0dloqxQ1JbKlLYnN1Iuuq4llVFs+68iicOzA4YyIj9oZrhu+1s+99ZlIsnRNjyUjyxqclxOjWjGHSas6hioi0ahGRkNbTG/pMOnB6XQ3s2Qq7t0DFdu82imXFUF6MlRcTU76dmD1rSSkvpntlSbNf4eKjaUjoTG2sF75lESnstk6UuE5sr0tkW10CBSXxbCqIY15lLMX1ndhDPLAvQCMMkuOjSYmPJjnOe02JjyY5Pprk+CiS46JJiovyhtjG995rclw0neKiFMghoEAVEQlUVMy+wD2c+lqoKIHybV7wlm/fG7yR5cVElm8nrryYtIoCqCyBqtID1xHtDS4iitqYFKqjU6mITGZPRBKlJLPLJVJSl8jOXbHs2B7D9ppYVtREU1ofSxnxlLl4yoinglgc+5/LTYyJ3Buy/oGbFBdNcjPjGudrDOvE2CidH25CgSoiEgqR0ZCU7Q2BqK+Dql1QudML4sqSva9WUUJMZQkxFSUkVe4ku7IEKlZ789RV7b+eKA74ze4w6qMSqIlMpCYygaqIBCotnnISKGuIY3d5HLv3xLGrPpaSulh21Mawvj6OMuL2hnK5i/OFcxwNvnCOiYwgITaSxJgoEmIiSYiNIjEmkoSYKBJjfa/7jY8kLjqS+JhI4qO9Ic7vfXzj9OhIoiOtzZ1TVqCKiLQGkVGQmOENR6K2EqrLoGaP91q9B2p8r77BasqIqi4jqno3CXunlUHNTqjevW+5hlpfLb7hIOoiYqmNiKMmIo4ai6OKWKpqYqmoiaPCxVDuYihviGV3fQx7GqLZXRfNLmKpJJYqF0MVvsH/PdFUuRiqiaGaaOoiYoiPjvIFcMS+AA4wkBuXifObNrhrMlEh3KtWoIqItGXR8d5A5rGvq656XxD7B2/17n2fayqIqi0nqqaC+NoKqK2AmsbXcqgtgdryfeNqK7xD10fIYdRZDLUulpqaWGpqvbCtIoZqF02li6bSxVDhoilviKa8IYpKF0MZMexwXijvH9rRPP6Lu4lKiD/2f6eDUKCKiIgnKtYbjnQv+VCc8/ai9wZupfcAhNoq73B1XZVv3P6vVldFdG0l0XXVJPjPv9+8e3zvq3C+eayu8qClNETcFbx2NUOBKiIioWMGMQneEMygbvo1jW+c8/a094Zwpfe5tpKImISQfT8oUEVEpD0xg+g4bwjd0d1mqc+ziIhIEChQRUREgkCBKiIiEgQKVBERkSBQoIqIiASBAlVERCQIFKgiIiJBoEAVEREJAgWqiIhIEChQRUREgsCcc+GuoVUys2JgQxBWlQFsD8J62oKO1FboWO1VW9sntfXI9XTONftoHwVqiJnZXOdcfrjraAkdqa3QsdqrtrZPamtw6ZCviIhIEChQRUREgkCBGnrPhLuAFtSR2godq71qa/uktgaRzqGKiIgEgfZQRUREgkCBKiIiEgQK1MMws/Vm5poZ3vGb51YzW2dmVWY2z8zGN1lHrJk9YWbbzazczN40s+5N5sk1s7d807eb2eNmFtNS7fTVEIy2zmxm+ZebzJNmZlPMrNQ3TDGz1BZqZmMNh2yrmZ3i204FvvHXNrMOM7P7zGyLmVX62n5ck3naS1snN7P8503mOezPeagF0NZ7zOxLM9ttZsW+/3NDmqyjvWzXQNraJrarr47Dtfc2M1vsa+9uM5ttZuc1WUdot61zTsMhBiAT6OI3jAQagGt8078L1AI3AoOAJ4AyINdvHX8EtgBnAKOAmcBCINI3PRJY4hs/yjffFuCJNtjWmcBfm6wnpcn3vAcsA8b6hmXAW62srecCDwLfASqAa5tZx8+APcC3gSHAq77tltQO2zoZ+KDJetKbzHPIn/NW0tZpwHW+7TUUeAMo9G9LO9qugbS1TWzXANt7AXAO0BfoD/wG7/fVsJbati32j9FeBuDnwC4g3vf5C+DZJvOsAh7yvU8BaoAr/Kb38P0gnOX7fI7vcw+/ea4EqoDkttJW3+eZwJOHWOcgwAHj/Mad7Bs3oLW0tcm0MpqEDGDAVuDnfuPiff9Zb25PbfWNnwy8fYh1HvbnvLW11Te9E1APfKM9b9fm2tqWt2sg7fXNU+K33UK+bXXI9wiYmQHfA553zlWad0j2eOD9JrO+D5zke388EO0/j3NuE7DCb56xwArf+EbTgFjf8i3uKNva6FLf4aFlZvaomSX5TRuL90v7M79xnwLlzaynRTRta4CL9cb7K9l/u1YCn7D/dm0PbW10spltM7OvzexZM8vymxbIz3mLCrCtSXinvnb6Prfn7dq0rY3a1HaFw7fXzCLN7FK8PyIat1PIt23UEbekYzsDb6M86/ucgXe4tqjJfEXA6b73XfD+Kmx6D8ki37TGeZquY7tvuS6Ex9G0FeBFvHsgbwGOAx4ChgFn+qZ3AYqd708/AOecM7NttJ62BqKx1ub+PXL85mkPbQWYCrwOrAN6AQ8A083seOdcNYH9nLe0QNr6B7zDl7N9n9vzdm3aVmib2xUO0l4zG4rXvji8YLzQObfENznk21aBemRuBL50zi0KdyEt4Kja6pzzv3h6iZmtBb4ws1HOuflBrTB4tF0Pwznn37FsiZnNw/vD6Ty8X8it0SHbamaP4R3OO9k5V9+ilQXfUbW1jW5XOHh7VwIj8A5Vfwd4zswmOueWtkRROuQbIN9hkAvY/y+ixr3I7CazZ+Od/Mf3Gom3h3eoeZquo3GPsJAWdgxtbc5c33L9fJ8LgUzfIZvG7zMg6zDrCYmDtDUQjbUebtu3h7YewDm3BdjM/tv1cD/nLeZwbTWz/wUuAyY559b6TWp32/UQbT1Aa9+ucOj2OudqnHOrnXPznHP34O2R/8g3OeTbVoEauGuBauClxhHOuRpgHt7hB39nsO8Y/Dy8nmZ75/F1OR/kN89sYFCTruhn+L5vXtBaELhrObq2Nmco3n/Irb7Ps/HOa4z1m2cskHiY9YTKtTRpa4DW4f0H89+uccB49t+u7aGtBzCzDLzDZI3bNZCf85Z0LQdpq5n9gX0B81WTye1qux6mrQdoA9sVjuznOAKvLwq0xLYNZy+ttjLg9Q77miY9XH3TvovXC+4GvB+yP+Adu+/pN88f8f7qOx2vq/cMmr9sZrpv+ulAAS182cyxthXoA/wSyMc7H3MuXueF+fh1scfrlr6Efd3Sl9DClxwE0NZOeIeORuBdSvJL33v/S4R+BpQCF+F1wX+Z5rvgt+m2+qY/6qu/FzAR7xfP5iZtPeTPeStp61PAbmAS+1+C0akdbtdDtrWtbdcA2vtbvHDshfeH/EN4vZHPaalt26L/GG11AE7F6zZ9wkGm3wqsZ98e5SlNpsfiXbO5w/cL6y38LpHxzZMLvO2bvgN4HIhtS23F607/sa/+amA1Xug2va4tDXje9599t+99amtqq++Xi2tmmOw3jwH34f01X+Vr+5D21la8SwumAdvw/qDagHe5RdOf4cP+nLeCtjbXTgfc1w636yHb2ta2awDtnexrQ7WvTR/S5NKeUG9b3RxfREQkCHQOVUREJAgUqCIiIkGgQBUREQkCBaqIiEgQKFBFRESCQIEqIiISBApUEQkZM/ubmf0yBOv9nZk9Eez1ihwLBapImJnZZDNzzQyfh7u2Y+F78se3gN8HOP8oX7vHH2T6K2bWePu3R4BrzCwvGLWKBIMCVaR1+BDo2mQ4N5RfaGZR/jcBD4E7gNecc7sDmdl5TyNaCFzfdJqZdcYL57/45i3Ge67lLUGqVeSYKVBFWodq51xhk6GkcaJvz+0mM/uHmZWb2Vozu9J/BWaWY2Yvm9lO3/COmfXzm36fmS01s2vNbA3eLdoSzay/mX1sZlVmttLMzjWzMjO71rfcdDN7ssl3JZtZhZld1FxjzCwSuATvNnX+42PM7GEz2+xb/kszO8tvlr8AF5tZpyarvNJX7yt+497Eu/G7SKugQBVpO34J/BsYjhcsfzWzXAAzS8C7aXkVMAHvpt5bgQ990xr1Bi4HLvatpwZ4A6gDTsR7ksev2PeEDvAek3W5mfmPuwzvwQj7BaafYXjPpJzbZPzffPVdjndz8ueAt8xsuG/6C3gPi/huk+W+B7zinCv3GzcHyDGzPgepQaRFKVBFWoezfXuF/sPDTeaZ4px73jm3GvhvvBA8xTftUrwbf1/nnFvsvEd13Yz3RJHz/dYRA1zlnJvvvIcunwoMAK52zi10zs3Ge35klN8yr+M9teNCv3HXA393ztUepD098W5i3vgYMHzBdxlwiXPuE+fcWufck8C7vlpxzu0CXsML0MblRuM9PeQvTb5ji++110FqEGlRUYefRURawCfATU3G7WryeXHjG+dcnZkV4z34GOB4vL3PPU1OiybgPVav0WbnXJHf54HAFudcgd+4L/ECtPG7qs1sCl6IvmxmxwEn4O3NHkw8UOuca/AbNwov9Jc3qTEW79GFjf4CzDCzgb4/DK4HljrnvmjyHZV+3yUSdgpUkdahwrfneShN9wYd+44yReB16Lm0meVK/N6XNzM9EH8BFvsOMV8PzHbOrTjE/NuBGDNLcM5V+NXogNEc2JZKv/cf4z3673oz+xXeXu19zXxHuu+1+EgaIhIqClSR9mE+XvBs9x02DdRXQDcz6+acazyEmk+T00HOuWVm9gVwI14HoZ8fZr0Lfa+D2XcedQHeHmoX59yMgy3onHNm9lfgB7764oEpzcw6BC+YlxymFpEWoXOoIq1DrJl1aTJkHsHyLwBFwL/NbIKZ9TazU8zsf/x7+jbjA2Al8JyZDTezE4HH8M7PNn1Y8rPAT4FE9u9tewDfZS3zgZP9xn3tq3OymX3HzPLMLN/MftJMb+HngAzgUeBfzrkdzXzNeGCW3x6wSFgpUEVah9PxOvD4DwsCXdgXKqcAa4F/4O3ZPQekATsPsVxjZ6NYvF6zzwG/wQvTqiazv4LXK/hV59yeAMp6Briiybjr8Hr6PuKr8W1f3Rua1LUFr7NSGgd2Rmp0GV7Ii7QK5lzTP0JFpCPzXcKyEMh3zs3zG98N2AhMcM59GsB64vBC8yrn3Kwg13ge8DtgmHOuLpjrFjlaOocq0sGZ2YV4nZVW4V2C8hiwCO+QLWYWDXQGHgQWBBKmAM65KjO7mn2dh4IpEe8SIYWptBoKVBFJAh4GeuAdHp4J/MjtO3w1Du+mEavw7n4UMOfcJ8Erc7/1vhqK9YocCx3yFRERCQJ1ShIREQkCBaqIiEgQKFBFRESCQIEqIiISBApUERGRIFCgioiIBMH/BxtKMxANeph+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "E = np.linspace(7000, 7300, 2000)  # eV\n",
    "\n",
    "fig = plt.figure(figsize=(7, 5))\n",
    "\n",
    "plt.plot(E, np.angle(ce3.Fhkl(E, [0, 0, 2]))*180/np.pi, label='Ce3+')\n",
    "plt.plot(E, np.angle(ce4.Fhkl(E, [0, 0, 2]))*180/np.pi, label='Ce4+')\n",
    "\n",
    "plt.xlabel('Energy (eV)')\n",
    "plt.ylabel('Phase (deg)')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71bf6c43",
   "metadata": {},
   "source": [
    "Although the phase difference between the structural models is less than 10 degrees, the phase shift can be at 180° for just 100 eV of variation in the X-ray energy. This phenomenon allowed a multi-wavelength use of X-ray dynamical diffraction exploited by [Morelhão, S. et al](https://link.springer.com/article/10.1557/mrc.2020.37)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87a1ab33",
   "metadata": {},
   "source": [
    "### 2.3 Phase triplets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a2a5a3",
   "metadata": {},
   "source": [
    "Until now, we selected the 002 as the primary reflection and might acquire data using X-rays in the range of 7100 and 7200 eV. \n",
    "\n",
    "Lastly, check the three-beam cases predicted to show opposite profile asymmetries on each model structure for a fixed energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27f9a7b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "E = 7150\n",
    "pyddt.triplet(ce3.diffraction(E), ce4.diffraction(E), [0, 0, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cd2f7ca",
   "metadata": {},
   "source": [
    "![pyddt_t2_2](https://user-images.githubusercontent.com/106104347/187935676-81303ffe-6583-427a-b885-d1ad713afd0f.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b68031e2",
   "metadata": {},
   "source": [
    "Again, the plotted information was saved in the current folder (filenames displayed on the screen)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a13e52d",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d265484e",
   "metadata": {},
   "source": [
    "## Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c660497c",
   "metadata": {},
   "source": [
    "The first part of this tutorial has shown some tools for structural modeling. However, the functions for appending or deleting atoms and changing their occupancy numbers weren't covered, as well the `visualizer_in` method. Please, check the [API documentation](../api.rst) for more details. On the other hand, the second part covered all functions usually used for experiment planning. If you develop new tools using **pyddt** classes, please contribute to the project on [GitHub](https://github.com/rafaela-felix/pyddt).\n",
    "\n",
    "In the last tutorial, we will analyze two Renninger scans and indicate how compatibility analysis is carried out using **pyddt**."
   ]
  }
 ],
 "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.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}