forel_ule_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "__maintainer__ = \"[José M. Beltrán](<jobel.open.science@gmail.com>)\"\n",
    "__modified_from = {\"code\":\"matlab\",\"credits\":\"Shungudzemwoyo Garaba\",\"email\":\"<shungu.garaba@uni-oldenburg.de>\", \"dates\":{\"created\":20110224, \"modified\":20140910}}\n",
    "__credits__ = [\"José M. Beltrán, Tiago A. M. Silva\"]\n",
    "__license__ = \"GPL-3.0\"\n",
    "__status__ = []\n",
    "__creation_date__ = 20150319"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook to calculate the Forel-Ule Index\n",
    "\n",
    "It allows you to convert hyperspectral remote sensing data into the Forel-Ule color index.\n",
    "\n",
    "The two articles will provide you with some updated information on the algorithm and constant parameters (FUI_ATAN2.dat or FUI_ATAN210.dat).\n",
    "\n",
    "test.dat is the sample reflectance file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**References**:\n",
    "\n",
    "Wernand, M. R., & van der Woerd, H. J. (2010). Spectral analysis of the Forel-Ule ocean colour comparator scale. Journal of the European Optical Society, 5. doi:10.2971/jeos.2010.10014s\n",
    "\n",
    "Wernand, M. R., Hommersom, a., & Van Der Woerd, H. J. (2013). MERIS-based ocean colour classification with the discrete Forel-Ule scale. Ocean Science, 9, 477–487. doi:10.5194/os-9-477-2013"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "from scipy.interpolate import interp1d\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[**Notes from the colour and vision research laboratory**](https://www.cvrl.org/)\n",
    "\n",
    "The CIE 1931 2-deg CMFs (CIE, 1932), which form the basis for most practical colorimetry, are based on the chromaticity coordinates obtained by Guild (1931) and by Wright (1928). Chromaticity coordinates, however, provide only a relative measure of the ratios of the three primaries needed to match each spectrum color, whereas CMFs specify absolute energy values. In order to reconstruct the CMFs from the Wright and Guild data, it was assumed that the CIE1924 V(l) photopic luminosity function (CIE, 1926) is a linear combination of the three CMFs (see Wyszecki & Stiles, 1982), for a description of the reconstruction and for the tabulated values. \n",
    "\n",
    "It has long been clear that the CIE1924 V(λ) that was used to construct the CIE 1931 2-deg CMFs **seriously underestimates sensitivity at wavelengths below 460 nm**, so that these CMFs are seriously in error at short wavelengths. The Judd and Judd-Vos modifications are attempts to overcome this problem.\n",
    "\n",
    "The cmf used here are the CIE 1931 CMF with Judd and Vos corrections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fui = {\"discrete\":pd.DataFrame(),\"continuous\":[]}\n",
    "fui[\"continuous\"] = pd.read_csv(sep = \"\\t\", filepath_or_buffer = \"data/FUI_ATAN210.tsv\", names = [\"value\", \"alpha\"])\n",
    "# see Table 3 of Wernand et al 2013 to refer to the angles, alpha, of FU number in degrees and\n",
    "# the 20 boundary angles or colour transition angles, alpha_T, that are used in the discrete classification of ocean colour\n",
    "\n",
    "fui[\"discrete\"] = pd.read_csv(sep = \"\\t\", filepath_or_buffer = \"data/FUI_ATAN2.tsv\", names = [\"alpha\",\"alpha_T\",\"value\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>alpha</th>\n",
       "      <th>alpha_T</th>\n",
       "      <th>value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>229.943868</td>\n",
       "      <td>227.677431</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>225.410993</td>\n",
       "      <td>219.271054</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>213.131116</td>\n",
       "      <td>205.190884</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>197.250652</td>\n",
       "      <td>189.202645</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>181.154639</td>\n",
       "      <td>165.707958</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        alpha     alpha_T  value\n",
       "0  229.943868  227.677431      1\n",
       "1  225.410993  219.271054      2\n",
       "2  213.131116  205.190884      3\n",
       "3  197.250652  189.202645      4\n",
       "4  181.154639  165.707958      5"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Showing the first 5 rows in the dataframe\n",
    "fui['discrete'][:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Colour matching functions (cmf) \n",
    "\n",
    "$\\tilde{x}, \\tilde{y}, \\tilde{z}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmf = pd.read_csv(sep = \"\\t\", filepath_or_buffer = \"data/FUI_CIE1931_JV.tsv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>wavelength</th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "      <th>z</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>380</td>\n",
       "      <td>0.002690</td>\n",
       "      <td>2.000000e-04</td>\n",
       "      <td>1.226000e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>385</td>\n",
       "      <td>0.005311</td>\n",
       "      <td>3.955600e-04</td>\n",
       "      <td>2.422200e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>390</td>\n",
       "      <td>0.010781</td>\n",
       "      <td>8.000000e-04</td>\n",
       "      <td>4.925000e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>395</td>\n",
       "      <td>0.020792</td>\n",
       "      <td>1.545700e-03</td>\n",
       "      <td>9.513500e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>400</td>\n",
       "      <td>0.037981</td>\n",
       "      <td>2.800000e-03</td>\n",
       "      <td>1.740900e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>85</th>\n",
       "      <td>805</td>\n",
       "      <td>0.000007</td>\n",
       "      <td>2.607600e-06</td>\n",
       "      <td>5.862600e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>86</th>\n",
       "      <td>810</td>\n",
       "      <td>0.000005</td>\n",
       "      <td>1.836500e-06</td>\n",
       "      <td>4.213800e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>87</th>\n",
       "      <td>815</td>\n",
       "      <td>0.000003</td>\n",
       "      <td>1.295000e-06</td>\n",
       "      <td>3.031900e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>88</th>\n",
       "      <td>820</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>9.109200e-07</td>\n",
       "      <td>2.175300e-10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>89</th>\n",
       "      <td>825</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>6.356400e-07</td>\n",
       "      <td>1.547600e-10</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>90 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    wavelength         x             y             z\n",
       "0          380  0.002690  2.000000e-04  1.226000e-02\n",
       "1          385  0.005311  3.955600e-04  2.422200e-02\n",
       "2          390  0.010781  8.000000e-04  4.925000e-02\n",
       "3          395  0.020792  1.545700e-03  9.513500e-02\n",
       "4          400  0.037981  2.800000e-03  1.740900e-01\n",
       "..         ...       ...           ...           ...\n",
       "85         805  0.000007  2.607600e-06  5.862600e-10\n",
       "86         810  0.000005  1.836500e-06  4.213800e-10\n",
       "87         815  0.000003  1.295000e-06  3.031900e-10\n",
       "88         820  0.000002  9.109200e-07  2.175300e-10\n",
       "89         825  0.000002  6.356400e-07  1.547600e-10\n",
       "\n",
       "[90 rows x 4 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cmf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Colour matching functions')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEWCAYAAAB2X2wCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABPqklEQVR4nO3dd3zV1f348df7Zi8SMhkBslhhQ9igAlZFcWG1uOcXqdWvbW1r+6utX2vtUlttrVontQ4cpU4coEVkE/YMZEAI2QnZO/f8/vjcYAgZN8nnjtycp+Zxcz+f8znn5EPyzsn5nCFKKTRN07S+z+LqCmiapmnm0AFd0zTNQ+iArmma5iF0QNc0TfMQOqBrmqZ5CB3QNU3TPIQO6Fq7ROQCEclxdT0cTURWishvOzlfJSIJDihXRORVETktItvNzr+Lsj8VkVudWabmHDqgezgRuUFEUm2BKc/2wzzP1fVyJBE5LiIXmpGXUipYKZVpRl5tzAO+A8QqpWY4IH8AROT/ROT11seUUouVUv90VJma6+iA7sFE5MfAU8DvgBhgOPAscKUL6+TtqrLdzAjguFKq2tUV0TyIUkp/eOAHEApUAdd2ksYPI+Dn2j6eAvxs5y4AclqlHQusB8qAg8AVrc6tB+5q9f42YGOr9wr4AXAMyGqnHnG2NLcDJ4HTwApgOrDPVuYzrdInAl8BJUAx8AYQZjv3L8AK1Nq+/p/Zjs8DNtvyOgncZju+Evg78AlQCWwDEtvUPcnOtBcBaUA5xi/Or1vfl1bp7gTqgGZbHR9pe896UPY4YC1QChQA/w+4BGgAGm3l7G3774XRqHsIOAEUAq8BoW3+XW4Fsm33+petypwBpAIVtjL/7Orv+/7+4fIK6A8H/cMaP8xNgHcnaX4DbAWigShbwHvUdu4CbAEd8AHSbUHCF1hoCyqjbefPBAjb+7OCky0orAXCgYB26tESOJ4H/G2BsQ5431a3obZgc74tfRJGd4Wfrd4bgKda5XccuLDV++G2+l5v+1oigMm2cyttQXAG4I3xy2FVm7ondZUWiLQFtqW2c/fbAuk5Ab2De3TW+26WHQLkAQ/Y7l8IMNN27v+A19vke+bfC7jD9m+bAAQDq4F/tfl3eREIACYB9cBY2/ktwM22z4OBWa7+vu/vH7rLxXNFAMVKqaZO0twI/EYpVaiUKsJoKd7cTrpZGD+wf1BKNSilvgI+xgiQ9vq9UqpUKVXbSZpHlVJ1SqkvgGrgLVvdTgHfAFMAlFLpSqm1Sql6W73/DJzfxde5Tin1llKqUSlVopTa0+r8aqXUdtu9egOY3EleHaW9FDiolFptO/dXIL+TfHqio7KXAPlKqSdt969SKbXNzjxvxGhZZyqlqoBfAMvadI09opSqVUrtBfZiBHYwfmEliUikUqpKKbW1t1+g1js6oHuuEiCyiz7rIRh/arc4YTvWXrqTSilrm7RDu1Gfk3akKWj1eW0774MBRCRaRFaJyCkRqQBex2ghd2QYkNHJ+daBt6alnG6mHUKrr1EppQCzRwl1VHZXX19n2vse8MZ45tJVuXcCo4AjIrJDRJb0sA6aSXRA91xbMLotruokTS7Gw7kWw23H2ks3TEQsbdKesn1eDQS2OjeonTzMXNbz97b8JiqlBgA3AdJJWScx+t0dKQ+IbXkjItL6vR3Ouoci0t497EhnX19X972974Emzv5l2n7GSh1TSl2P0S32R+A9EQnqurqao+iA7qGUUuXAr4G/i8hVIhIoIj4islhE/mRL9hbwkIhEiUikLf3r7WS3DSPg/MyWxwXA5cAq2/k9wFJbGUkYLTdHCsF4yFcmIkOBn7Y5X4DRJ9ziDeBCEblORLxFJEJEJptcp0+ACbZ77Y3xELg7QXkvME5EJouIP0bft70+BgaJyA9FxE9EQkRkpu1cARDX5pdxa28BPxKReBEJxhgR9XYXXXUAiMhNIhJl+8utzHa4uRv11kymA7oHU0r9GfgxxiiGIoyW3L0YDxsBfosxSmEfsB/YZTvWNp8G4ApgMcZIh2eBW5RSR2xJ/oIxmqIA+CdGAHWkR4CpGKNJPsF4kNfa7zF+UZWJyE+UUtkYfdwPYDxY3MO3/cCmUEoVA9cCf8Lo7krGuLf1dl5/FOMh9TqM0UAbu1F2JcZD4ssxukeOAQtsp9+1vZaIyK52Ln8FY2TQBiAL46+6++ws+hLgoIhUAU8Dy5RSdfbWWzOfGF19mqaZydYizgFuVEr919X10foH3ULXNJOIyMUiEiYifhhDPAVjWKimOYUO6JpmntkYo02KMbo/rupimKammUp3uWiapnkI3ULXNE3zEC5bKCkyMlLFxcW5qnhN07Q+aefOncVKqaj2zrksoMfFxZGamuqq4jVN0/okETnR0Tnd5aJpmuYhdEDXNE3zEDqga5qmeQi9e4ymaf1OY2MjOTk51NW570oF/v7+xMbG4uPjY/c1OqBrmtbv5OTkEBISQlxcHMbCmO5FKUVJSQk5OTnEx8fbfZ3uctE0rd+pq6sjIiLCLYM5gIgQERHR7b8gdEDXNK1fctdg3qIn9esyoIvIKyJSKCIHOklzgYjsEZGDIvJ1t2vhYaobq/kg/QMamxtdXRVN0/oRe/rQVwLPYOwGfg4RCcNYH/sSpVS2iESbVrs+SCnFLzf+ki+zv6Sotoi7Jtzl6ippmtZPdNlCV0ptwNgUoCM3YGxem21LX2hS3fqkNw6/wZfZXxIdEM0L+16gsKZf3w5N05zIjD70UcBAEVkvIjtF5JaOEorIchFJFZHUoqIiE4p2L3uL9vJk6pMsGLaAVy95lSZrE0/vetrV1dI0zc3s2LGDiRMnUldXR3V1NePGjePAgQ57te1mxrBFb2AasAgIALaIyFbbllpnUUq9ALwAkJKS4lHr9pbXl/PTr39KTFAMj859lFC/UG5OvplXDrzCstHLmBA1wdVV1DStHY98dJBDuRWm5pk8ZAAPXz6uw/PTp0/niiuu4KGHHqK2tpabbrqJ8ePH97pcM1roOcBnSqlq276KGzB5v8a+4PEdj1NcW8yT5z9JqF8oAMsnLicyIJI/7PgDet15TdNa+/Wvf83atWtJTU3lZz/7mSl5mtFC/wB4xrbTuS8wE2PT4H7DqqxsyNnA4vjFjIv89rdykE8Q90+9n19t+hWfZH3CkoQlLqylpmnt6awl7UilpaVUVVXR2NhIXV0dQUFBvc7TnmGLbwFbgNEikiMid4rIChFZAaCUOgx8hrFz/HbgJaVU7zuD+pCMsgxO159m+qDp55y7IvEKEkMTee/oey6omaZp7mr58uU8+uij3HjjjTz44IOm5NllC10pdb0daR4HHjelRn1QaoGxrntKTMo55yxiYeHwhbxy4BUqGioY4DvA2dXTNM3NvPbaa3h7e3PDDTfQ3NzMnDlz+Oqrr1i4cGGv8tUzRU2wI38Hg4MGMzR4aLvn58fOp1k1syV3i5NrpmmaO7rllltYvXo1AF5eXmzbtq3XwRx0QO81pRQ7C3YyfdD0DqfqToicwADfAWzI2eDk2mma1p/ogN5LmeWZlNaVttvd0sLb4s3cIXPZeGojVmV1Yu00TetPdEDvpR35OwBIGdRxQAej26W0rpTDJYedUS1N0/ohHdB7aUf+DmICY4gNju003dyhcxGEDad0t4umaY6hA3ovKKVILUjttP+8Rbh/OBMiJ7AxZ6OTaqdpWn+jA3ovZJVnUVpX2u748/bMi53H/uL9lNZ1ttaZpmlaz+iA3gtn+s87eSDa2nlDz0Oh2HRqkyOrpWlaP6UDei/sKNhBdGA0w0KG2ZV+bMRYIvwj+CbnGwfXTNO0/kgH9B5SSpGab1//eQuLWJg3dB6bcjfRZG1ycA01TXNXv/rVr3j66W+X1v7lL3/JX//6117na8biXP1SdmU2JXUlTIuZ1q3r5sXO44OMD9hfvJ8p0VMcVDtN0+z26c8hf7+5eQ6aAIv/0OHpO++8k6VLl3L//fdjtVpZtWoV27dv73WxOqD30JHSIwAkRyR367rZg2fjJV58k/ONDuia1k/FxcURERHB7t27KSgoYMqUKURERPQ6Xx3QeyitNA0v8SIpLKlb14X6hTIpahIbT23kf6f+r4Nqp2ma3TppSTvSXXfdxcqVK8nPz+eOO+4wJU/dh95DR08fJW5AHH5eft2+dt7QeRwuPUxxbbEDaqZpWl9w9dVX89lnn7Fjxw4uvvhiU/LUAb2H0k6nMSp8VI+unR87H4CNp/QkI03rr3x9fVmwYAHXXXcdXl5epuSpA3oPlNeXk1+dz+iBo3t0/eiBo4kKiNIBXdP6MavVytatW7nzzjtNy9OeHYteEZFCEel0FyIRmS4izSLyXdNq56aOnjb2vx4d3rOALiLMHTqXzbmb9fBFTeuHDh06RFJSEosWLWLkyJGm5WtPC30lcElnCUTEC/gj8LkJdXJ7ZwJ6D1voAPOHzqeyoZL9xSYPl9I0ze0lJyeTmZnJk08+aWq+XQZ0pdQGoKvFR+4D/g0UmlEpd5dWmka4fziRAZE9zmPWkFlnhi9qmqaZodd96CIyFLgaeN6OtMtFJFVEUouKinpbtMuknU5j1MBRds8Qbc8A3wFnhi9qmqaZwYyHok8BDyqlmrtKqJR6QSmVopRKiYqKMqFo52uyNpF+Or1X3S0t5sfO18MXNU0zjRkBPQVYJSLHge8Cz4rIVSbk65ZOVJygwdrQ4weirc0fagxf1N0umqaZodcBXSkVr5SKU0rFAe8B9yil3u9tvu4qrTQNgFEDezYGvbVRA0cxJGgIa0+s7XVemqZp9gxbfAvYAowWkRwRuVNEVojICsdXz/2knU7D2+JNQmhCr/MSES6Ou5gtuVsory83oXaapvVn9oxyuV4pNVgp5aOUilVKvayUel4pdc5DUKXUbUqp9xxTVfeQdjqNhNAEfLx8TMnv4viLaVJNfJn9pSn5aZrm/p5//nkmT57M5MmTiY+PZ8GCBabkqxfn6qajpUeZNXiWafklhyczLGQYn2V9xtKRS03LV9M0+/xx+x/PrJ5qljHhY3hwxoMdnl+xYgUrVqygsbGRhQsX8uMf/9iUcvXU/24orSulqLbIlAeiLUSES+IuYXv+dr3XqKb1M/fffz8LFy7k8ssvNyU/3ULvBjMfiLZ2cdzFvLj/RdadWMd1o68zNW9N0zrXWUvakVauXMmJEyd45plnTMtTt9C7obdruHRk1MBRxIfG89nxz0zNV9M097Rz506eeOIJXn/9dSwW88KwDujdcLziOAP9BhLuH25qvi3dLqn5qRTV9N0ZtJqm2eeZZ56htLSUBQsWMHnyZO666y5T8tVdLt2QX53PoKBBDsn7krhLeG7vc3xx4gtuHHujQ8rQNM09vPrqqw7JV7fQu8GRAT0hLIGRA0fyadanDslf0zTPpwN6NxRUFzgsoAMsTVrK3qK97CzY6bAyNE3zXDqg26m6sZrKxkpiAmMcVsY1o64hwj+C5/Y+57AyNE0zKKVcXYVO9aR+OqDbqaC6AMChLfQA7wBuH3872/K2satgl8PK0bT+zt/fn5KSErcN6kopSkpK8Pf379Z1+qGonfKr8wHHBnSA60ZfxysHXuG5vc/x4kUvOrQsTeuvYmNjycnJwZ33ZfD39yc2NrZb1+iAbqf8GucE9ADvAO4YfwdPpD7B7sLdTIme4tDyNK0/8vHxIT4+3tXVMJ3ucrFTfnU+ghAdEO3wsq4ddS3h/uE8t0f3pWuaZj8d0O1UUFNARECEaassdibQJ5A7xt/BlrwtHCw+6PDyNE3zDDqg2ym/Op9BgY7tbmntqqSrANicu9lpZWqa1rfZs8HFKyJSKCIHOjh/o4jss31sFpFJ5lfT9Rw5qag9oX6hJIQmsLtwt9PK1DStb7Onhb4SuKST81nA+UqpicCjwAsm1MutKKWcHtABpkRPYU/RHqzK6tRyNU3rm+zZsWgD0OFC3UqpzUqp07a3W4HujbPpA6oaq6hpqnHopKL2TI6eTGVDJVnlWU4tV9O0vsnsPvQ7gQ4XIxGR5SKSKiKp7jz+sy1njUFvq2XIou520TTNHqYFdBFZgBHQO1wtXin1glIqRSmVEhUVZVbRDueqgD48ZDjh/uE6oGuaZhdTJhaJyETgJWCxUqrEjDzdibMmFbUlIkyKmsSewj1OLVfTtL6p1y10ERkOrAZuVkod7X2V3E9BdQEWsRAZEOn0sqdETyG7Mpvi2mKnl61pWt9iz7DFt4AtwGgRyRGRO0VkhYissCX5NRABPCsie0Qk1YH1dYn86nwiAyLxtjh/pYSWfvS9RXudXramaX1LlxFKKXV9F+fvAszZP8lN5dc4f8hii+SIZHwsPuwp3MOi4YtcUgdN0/oGPVPUDgXVBU6dJdqar5cv4yLG6QejmqZ1SQf0LiilKKgpICbIuWPQW5sSPYVDJYeob653WR00TXN/OqB3oaKhgtqmWpe10MGYYNRobdQLdWma1ikd0LvgqjHorU2OngzoCUaapnVOB/QuuENAD/cPJ25AnB7pomlap3RA70JLQHf2Oi5tjQ4fTXpZukvroGmae9MBvQsFNQV4i7dLJhW1lhCaQE5lDnVNdS6th6Zp7ksH9C7kV+cTFRiFl8XLpfVICEtAoThecdyl9dA0zX3pgN4FV04qai0xNBGAzLJMF9dE0zR3pQN6F5y99VxHRgwYgZd4kVGe4eqqaJrmpnRA74RSisKaQqIDo11dFXy9fBkWMky30DVN65DzV5vqQ6obq6lvrnf5A9EWCaEJuoVuovqmZj7am0dGURXZpTXknK5l2vCBPHDRKIL89I+G1vfo79pOnK4zdtYb6D/QxTUxJIYlsiFnA43Njfh4+bi6On1aQ5OVH7yxi3WHC/HxEoaGBRAd4s+rm7NYd7iAx787kZkJEa6upqZ1iw7onSipM/bqCPcPd3FNDPGh8TSpJrIrs0kMS3R1dfqshiYrP3jTCOa/uXIcN84cgZdFANiWWcJP39vH917Yyl3z4vnlZWMRERfXWNPso/vQO9HSQneXgN4SxDPKdLdLTzU2W7nvrV2sPVTAb64cxy2z484Ec4CZCRF89sP5XD9jOC9tzOLd1BwX1lbTuseeDS5eEZFCETnQwXkRkb+KSLqI7BORqeZX0zVK60oB9wno8aHxCEJmuX4w2lO/ev8Anx8s4OHLk7lldly7aQJ9vXnsqvHMiAvnsTWHKazUk7m0vsGeFvpK4JJOzi8GRto+lgPP9b5a7uF0vXv1oQd4BzAkeIge6dJDWzNLWLXjJHefn8Dtc+M7TWuxCL9bOoHahmYe+eiQk2qoab3TZUBXSm0ASjtJciXwmjJsBcJEZLBZFXSl0rpSAr0D8ff2d3VVzvD4kS4NNVBbZnq2jc1WfvX+AYaGBfDDRaPsuiYpOpj7Fibxyb481h0qML1OmmY2Mx6KDgVOtnqfYzuWZ0LeLlVaV+o2rfMWiWGJbMvbRrO12eXLEZiivhJ2/QtydkD+fijNAGWFAbEQMw4GT4SUO2DAkF4V88rGLI4VVvHSLSkE+Np/3+4+P5GP9+Xxqw8OMDMhnBB/PbpIc19mPBRtbwiAajehyHIRSRWR1KKiIhOKdqzTdafdpv+8RUJoAg3WBk5VnXJ1VXqnqR62Pg9PT4bPfwGnUiFyFJz3U7jw/2DEHCg/Cd/8Gf6WAhuegMae9WXnltXy1LpjXDg2hguTu7dqpq+3hd9fM4H8ijr+8bXu6tLcmxkt9BxgWKv3sUBuewmVUi8ALwCkpKS0G/TdSWldqcuXzW0rISwBMEa6DB8w3MW16aGsb+CDe6AsG+Lmw4WPQOy09tOWZsEXD8FXj8Luf8GSpyBxQbeK+81Hh1AoHr48uUfVnTp8IN8ZG8Mb205w78Ik/H084C8jzSOZ0UL/ELjFNtplFlCulOrz3S3gnl0uCaG2gN5X+9EPfwyvLwVvf7hpNdz6UcfBHCA8Hpa9Abd8AF5+8MZ34eB/7C5u54nTfHYwn3sXJDEsPLDH1b5jXjynaxr5z+4+/peR5tHsGbb4FrAFGC0iOSJyp4isEJEVtiRrgEwgHXgRuMdhtXUipRSldaVu1+US4htCdGB03xzpsvsNeOdmGDwJ7vgckhaBvZN2Ei6Au9ZC7HR47w7Y86Zdl/3j6wzCAn24Y17no1q6MjM+nOTBA3hlYxZKuf0fl1o/1WWXi1Lq+i7OK+AHptXITVQ1VtFkbXK7gA5GK73PjUXf9g/49GeQsAC+9zr4BXc/D/9QuOnfsOoGeP/70FgD0+/qMHl6YRVrDxdw34IkAn1717soItw5L54H3t3LxvRi5o+M6lV+muYIeqZoB9xtHZfWEsMSySzPxKqsrq6KfbI2wKcPwpglcMPbPQvmLXyD4Pq3YdRi+OQBOPxRh0lf+iYTXy8Lt8yJ63l5rSyZNJjIYD9e3phlSn6aZjYd0DvgbrNEW0sITaC2qfbMfqduraYUVt8NEYmw9AXw9ut9nj7+cN0/YchUeP8eKDn3eUJhRR2rd53i2pRYIoNNKBPw8/bi5lkjWJ9WRHphlSl5apqZdEDvQEtAd8cWenyo0R98vPy4ayvSFaXgw/ugugiuedloXZvF288I6mKBd26FxtqzTr+6+ThNVit3zUswr0zgxlnD8fW2sHKzbqVr7kcH9A60BPQIf/dbQjVuQByA++8vuuufcORjuPBhGDLZ/PzDhsPSF6FgP6z5yZnDlXWNvL71BIvHDyYu0sRfIkBksB9XTBrCf3adoq6x2dS8Na23dEDvgDv3oUcGRBLoHejeAb0kAz79ufEQdJYDn5mPusiYjLT7ddjzFgBv7zhJZV0Ty88zt3Xe4srJQ6huaGbDUfefHKf1Lzqgd6C0rpQgnyD8vMzpfzWTiBAXGseJihOurkrHvnwELF5w9fNgcfC32QW/gOGz4fNfoKqLeXN7NlOHhzFpWJhDipuVEEFYoA+fHugDzzC0fkUH9A6U1pUy0M/9WuctRgwY4b4BPXc3HPoAZt8LIU7YYNviBUv+AvWVlLz/CzKLqlk23XGzaH28LFyUHMO6QwXUN+luF8196IDeAXecVNRa/IB4cqtyqWtyw7W6v/wNBITDbCdOT4geC7N/QOSxd5jre4xLJzp2wc/FEwZTWd/EpvRih5ajad2hA3oH3HFhrtZGDBiBQnGy8mTXiZ0p6xvI+Arm/xj8Bzi16KpZD5CrInkiYCXB3o6dzTk3MZIQf28+2ae7XTT3oQN6B07XnXbLB6It4kLjADcb6aKU0XceMqTTGZyO8smRch5uvIXB9Vmw9VmHluXrbeE7yTGsPZRPQ1MfmeCleTwd0NuhlKK03r27XEYMGAHgXv3oaZ8a65pf8CD4BDi9+Ld3nCQz4nzUqEtg/R+hyrGjUC6bMJiKuiY2Z+huF8096IDejsrGSpqsTW7dQg/yCSIqIMp9JhcpBV//EcITYfJNTi8+vbCSXdllfG/6MOSi30JTLWz8s0PLnDcykmA/bz7dr7tdNPegA3o7Smvdd9p/a3Ghce7T5ZKTCnl7YNb3wcuMZfa75+0dJ/G2CEunxkLkSJh0A+x4Gcodt9ytn7cXF46N5vND+TQ2624XzfV0QG9Hy+bQ7h7Q3Wro4o6XwDcEJi1zetGNzVb+s/sUC8dEf7tuy/k/M7ay2/Anh5a9eMJgymoa2Z7V2ba7muYcOqC3o6WF7s5dLmAsAVBWX0ZZXZlrK1JdDAdXG8HcL8TpxW9KL6a4qoFrpsV+e3DgCJh2mzGDtNRxSw3PHxmJr5dFzxrV3IIO6O0ore8jXS7usqbLrtegucElI1sAPtyTS4i/NxeMbrNG+Xk/AYsPrP+Dw8oO9PVm6ogwNhzTD0Y117MroIvIJSKSJiLpIvLzds6HishHIrJXRA6KyO3mV9V5WtZxcfuA7g5DF63NkPqqsTdo9BinF1/b0MznB/O5dPxg/Lzb7PUZMghm/A/sewcKDzusDvNHRnE4r4KiynqHlaFp9rBnCzov4O/AYiAZuF5E2u62+wPgkFJqEnAB8KSI+JpcV6cprSsl2CcYXy/3/hKGBA/BW7xd249+9HMozzYCpwt8daSQ6oZmrpw8pP0E835kLNv7jeNGvMxLigTQwxc1l7OnhT4DSFdKZSqlGoBVwJVt0iggREQECAZKgSZTa+pE7rg5dHt8LD7EhsS6NqDveNGYSDT6MpcU/8GeU0SH+DEzoYNljgPDjb70A/+G08cdUofxQ0MJDfDhG93tormYPQF9KNB6fnmO7VhrzwBjgVxgP3C/UufujyYiy0UkVURSi4rc9yFSXwnoYPSjZ5W7aLOF0ixjmn/K7S4Zqlhe08j6tCKWTByCl6WTzaZn32ss4LXprw6ph5dFmJsUwcZjxXoDac2l7Ano7f2ktP2uvRjYAwwBJgPPiMg5C3kopV5QSqUopVKiotx3k113X8eltbjQOLIrsl2zv+ihD4xXFwxVBPjsYB4NzdaOu1taDBgMk643RrxUFjikLvOSosivqCOjSG9Np7mOPQE9BxjW6n0sRku8tduB1cqQDmQBzn9CZhJ3X2mxtREDRtBgbSCvOs/5hR9639jXM8xxS9V25oM9ucRFBDIxNrTrxHPvB2ujw9Z4mT/S6EfX3S6aK9kT0HcAI0Uk3vagcxnwYZs02cAiABGJAUYDjhv860BWZaWsrqzPBPSWoYsnyp3cj376hLHueXLbxynOUVBRx5bMEq6YPBTj0U0XIhIh+Spj9mhtmen1GRYeSFxEIBt1QNdcqMuArpRqAu4FPgcOA+8opQ6KyAoRWWFL9igwR0T2A18CDyql+uR3dmVDJU2qya03t2itZehiVoWT+9EP236nuyigf7IvD6XgiklddLe0Nu9H0FBpzGp1gHkjI9maWaKXAdBcxq4nWUqpNcCaNseeb/V5LnCRuVVzjZbNocMD+kYLPcI/gmCfYOcv0nXoAxg8CcLjnVuuzcf7chkzKISk6GD7Lxo8EZK+A9ueNx6U+vibWqd5SVG8vjWb3dllzIjvG98/mmfRM0XbODOpyK9v/ECKCAmhCWSWO7GHqzzHWCbXRa3zU2W17Mou4/LutM5bzLkPqotg/zum12t2YgQWgW+Oue8ILs2z6YDeRksLva8MWwRIDEskoyzDeQUesnW3jHVNQP90v/EA+NIJPdhmLv48GDQBNj8DVnO7RkIDfJg0LIyNels6zUV0QG/jTJdLH3koCkZAL6krcd4iXYc+gJjxEJnknPLa+HhfHuOGDCA+Mqj7F4vA7PugOA3S15let7mJkezLKaeyrtH0vDWtKzqgt1FYU4hFLEQEdDDz0A0lhCYAOKfbpSIXTm41Roy4wMnSGvacLOOy3mwCPX6pMbt1y9/Mq5jNnKQImq2KbZl6OV3N+XRAb6OwppAI/wi8Lc6f+dhTiWGJAKSXpTu+sMMfG68u6j//9IDR3bJkQg/6z1t4+cCsFZC1AfL2mVQzw9ThA/HztrBJr+uiuYAO6G0U1hYSFei+s1jbMzhoMIHegc5poR/9DCJGQtQox5fVjk/25TExNpThEYG9y2jqreAbDFueMadiNv4+XkyPC2dzeomp+WqaPXRAb6OopojogGhXV6NbWka6OPzBaGMdnNgESRc6tpwOnCytYW9OOZf15GFoWwFhMPUWY9Euk7epm5MUQVpBpV5OV3M6HdDbKKwpJDqwbwV0gISwBDLLHNxCz94MTXWQuNCx5XTg4329GN3SnpkrjG3qdrxoTn42cxP1crqaa+iA3kpDcwNl9WV9rssFjH70wtpCKhoqHFdIxlfGDkBxcx1XRic+2Z/LpGFhDAvvZXdLi4EjYMwSY4OOhhpz8sRYTjfE31t3u2hOpwN6K4U1hQDEBMa4uCbdlxhqPBh1aCs9478wfJaxYYSTnSip5sCpCpaY1TpvMeseqCuDfatMy9LLIsxKiNAPRjWn0wG9laJaY4ZfX22hA47rR6/Mh4IDkLTIMfl34RPbZKLFEwaZm/HwWTB4Mmx9ztSJRnMTI8g5XUt2iXktf03rig7orRTUGGtl98U+9CHBQ/D38iej3EEBPeO/xquL+s8/2ZfH5GFhxA40qbulhYjRSi8+CplfmZbtXL0tneYCOqC3UlRjtND72igXAItYiA+Nd1wLPeMrCIyEmAmOyb8Tx4urOZhbwZLeTCbqzLirIXgQbDFvrfSk6GCiQ/zYlKH70TXn0QG9laKaInwtvoT62bFhghty2JouVitk/hcSF4DF+d8y33a3OCige/vCjLsg40soPGJKliLCnMQItmTobek059EBvZWCmgKiAqPs2zDBDSWGJVJQU0BVg8nboBUcMFYoTHRN//ma/XlMGR7G0LAAxxUy7Xbw9jeW1jXJnKRIiqsaOJxXaVqemtYZuwK6iFwiImkiki4iP+8gzQUiskdEDorI1+ZW0zmKaov6ZP95izMjXcyeMZrxpa2ABebma4eW7hZTJhN1JigSJlwLe1dBjTnrsJw30ni4vkEvp6s5SZcBXUS8gL8Di4Fk4HoRSW6TJgx4FrhCKTUOuNb8qjpeX51U1MJhI10yvoLocRBi8ggTOzi8u6W1Wd+HplrY9Zop2Q0K9WfMoBC+TtMBXXMOe1roM4B0pVSmUqoBWAW0XZnpBoxNorMBlFKF5lbT8ZRSFNYUEhXQ94YsthgaPBRfi6+5Ab2hBrK3uqR1DsbolqmO7m5pETMO4ubD9hehucmULM8fFUXqiVKq6s3JT9M6Y09AHwqcbPU+x3astVHAQBFZLyI7ReQWsyroLNWN1dQ21fbJSUUtvCxexkgXM4cuZm+G5gaXBPSs4moO5VWYN9XfHrO+DxU5cOQjU7I7f3QUjc2KLXq0i+YE9gT09p4Qtn1s7w1MAy4DLgZ+JSLnLMcnIstFJFVEUouK3OvP0JZZon1xUlFriWGJ5s4WzVwPXr4wfI55edrp4725gIlrt9hj1CUQNgK2mvNwNGVEOIG+Xnx9tM/90ar1QfYE9BxgWKv3sUBuO2k+U0pVK6WKgQ3ApLYZKaVeUEqlKKVSoqLcK3AW1ho/cH25Dx1g5MCR5FbnmremS+Z6GDYTfE2e0GOHj/flkTJiIEOc0d3SwuIFM+82NvHI3d3r7Hy9LcxJjGR9WpEevqg5nD0BfQcwUkTiRcQXWAZ82CbNB8B8EfEWkUBgJnDY3Ko6VksLva8H9LHhYwE4XGLC7a8uhvz9kHBB7/PqpmMFlaQVVDpuMlFnptxkrJVuUiv9/NFR5JyuJau42pT8NK0jXQZ0pVQTcC/wOUaQfkcpdVBEVojICluaw8BnwD5gO/CSUuqA46ptvjNdLn34oShAcoQxAOlQyaHeZ5ZlG32a4Pz+84/35SHi5O6WFv6hMPkGY630yoJeZ3e+bfji10fdq5tR8zx2jUNXSq1RSo1SSiUqpR6zHXteKfV8qzSPK6WSlVLjlVJPOai+DlNYU0iITwiBPs7vWjDTQP+BDAkaYk5Az1wPfqEwZHLv8+oGpRQf78tlZnw40QP8nVr2GTNXgLUJdrzU66yGRwSSEBmkA7rmcHqmqE1RTVGffyDaIjkimYMlB3uXiVKQsR7i5xv9yk50JL+SjKJqlkzsxb6hvRWRaDwgTX3Z2Kmpl84bFcXWzBLqGptNqJymtU8HdJvC2r49qai15IhkTlae7N2D0dJMKM92Sf/5x/tysQgsHu/8iUxnmX0P1JTA/nd6ndX5o6Ooa7SyPcucWaia1h4d0G36+izR1lr60Xv1YDRzvfHq5P5zo7sljzmJkUQE+zm17HPEzYeY8cYqjL0coTIrPgJfbwv/TdPDFzXH0QEdsCorxTXFff6BaAtTHoxmrofQYUbXgxMdOFXBiZIa14xuaatlrfSiw9/+guuhAF8v5iRG8NWRQj18UXMYHdCB0rpSmlSTx7TQWx6M9rgf3doMWRsg4XwjqDnRx/ty8bYIl7i6u6XFhO9CUDRs7f1a6YvGxnCipIaMIpNXw9Q0Gx3QabWxhYcEdDBa6T1uoeftNfbZdHJ3i9Wq+GhvLvNHRhIW6OvUsjvk7QfT74JjX0DR0V5ltWiM8f217rDudtEcQwd0PGdSUWu9ejCaYduKLf48cyvVhe3HS8ktr+OqKW2XCnKxlDvAy6/XrfQhYQEkDx7Al4d7P7Zd09qjAzqeM+2/tXER44AePhhNXweDJ0Gwc+/HB3tOEejrxXeS3WyBtOAomLQM9rwJVb0bS37h2Gh2njjN6eoGkyqnad/SAR2jy0UQIgIiXF0V07Q8GO12P3rtaTi5DUZe5IBaday+qZlP9uVxUXIMgb7eTi3bLnPuM1ad3P5Cr7JZNDYGq0KPdtEcQgd0jC6XcP9wfCw+rq6KacL8w3o2YzTjv6CskPQdx1SsA1+nFVFR18SV7tbd0iJyJIy+FHa8CA09X5NlwtBQokL8+FL3o2sOoAM6njUGvbUePRhNXwf+YRCb4pA6deSDPblEBPkyPynSqeV2y9z/Nf6C2f1Gj7OwWIRFY6L5+mgRDU1WEyunaTqgA8bm0J4a0E9WnqS8vty+C6xWOLYWkhY5dbp/ZV0j6w4XsGTiYLy93PhbcvgsiJ0BW57p1Y5Gi8bGUFXfxLYsvemFZi43/ulxjmZrMycqThA3IM7VVTHdmQejpXY+GM3fB9WFTu9u+exAPvVNVvftbmltzn1QdgIOt11B2n7zkiLx87bobhfNdP0+oOdW5VLfXH9mg2VPMi7SCOh7CvfYd0H6WuM1aZFjKtSBD/bkMiIikCnDwpxabo+MuQzCE2DzX3u8HECArxdzkyJZd7hAzxrVTNXvA3rL/psJYQkuron5Qv1CGRs+lm152+y74Ng6GDLFqcMVCyrq2JxRzJWThiBOnpXaIxYvo5WeuxuOb+xxNt9JjiHndC2H8kzaWUrT0AGdjDJbQA/1vIAOMGvwLPYU7aGmsabzhDWlkLPd6d0tq3edwqpg6dRYp5bbK5Ouh8BI2PR0j7O4KDkGixjdTZpmFrsCuohcIiJpIpIuIj/vJN10EWkWke+aV0XHyizPJDowmhDfEFdXxSFmDp5Jk7WJ3YVd7I+ZaRuuONJ5AV0pxbs7TzI9biBxkUFOK7fXfAKMDTDS10JBz9bLiQj2Y0Z8OJ/qgK6ZqMuALiJewN+BxUAycL2IJHeQ7o8YW9X1GRllGSSGel7/eYsp0VPwsfiwNW9r5wmPrYOAgTB0mnMqBuzKLiOzqJprpw3rOrG7mX4n+ATC5r/1OItLJwwmvbCKYwWVJlZM68/saaHPANKVUplKqQZgFXBlO+nuA/4N9JlH91ZlJbM80yMfiLYI9AlkcvTkzgN6cxMc+xySLnTqcMX3duYQ4OPFpe6wVG53BYbD1Fth/7tQntOjLC4eZ6woqVvpmlnsCehDgZOt3ufYjp0hIkOBq4FOt0kXkeUikioiqUVFrt9fMb86n9qmWo98INrarMGzOFJ6hNK6DnbLyfra2Jkn+Sqn1am2oZmP9+ayeMIggv3ccKq/PWbfY4x02fpcjy6PGeDPtBEDdUDXTGNPQG9v6EHbsVZPAQ8qpTrdMFEp9YJSKkUplRIV5frNJFoeiHpylwsYAR1ge/729hMcXA1+A4wWupN8fjCfyvqmvtnd0iJsOIy/BnauNGaQ9sDi8YM4nFfBiZKeLyegaS3sCeg5QOufulggt02aFGCViBwHvgs8KyJXmVFBR8oszwQ8d4RLi+SIZEJ8Qtia2063S1MDHP7IWKfEx99pdXpvZw6xAwOYGR/utDIdYu7/QkMV7Hi5R5frbhfNTPYE9B3ASBGJFxFfYBlw1jQ5pVS8UipOKRUHvAfco5R63+zKmi2jLIMI/wjC/MNcXRWH8rZ4kzIopf1+9Mz/Ql05jLvaafU5VVbLpoxivjstFoulD4w978ygCZC4CLb9Axrrun35sPBAJsaG6oCumaLLgK6UagLuxRi9chh4Ryl1UERWiMgKR1fQkTLKMzz6gWhrswbP4lTVKU5Wnjz7xIHV4B8KiQudVpf3UnNQCq7pS2PPOzP3fmPJhL1v9ejyS8YPYu/JMk6V1ZpcMa2/sWsculJqjVJqlFIqUSn1mO3Y80qpcx6CKqVuU0q9Z3ZFzaaUIrMs0+O7W1rMGmL0o5/VSm+sg7Q1MOZy8HbOlm/NVsXbO7KZlxTJsPBAp5TpcPHnweDJxhBGa6ePkdq1eLwxymfNvjyTK6b1N/12pmhhTSFVjVX9poUePyCe6MDos/vRM76E+goY77zulq+PFpJbXscNM4c7rUyHE4F5P4TSDDjycbcvj48MYlJsKKt3nzK/blq/0m8DessaLmYF9MKKOj7dn8drW47z5Bdp/G7NYbJLuphu70Qiwryh89iUu4naJtuf9gdWQ0A4xJ/vtHq8uS2bqBA/99tmrrfGXgED42HjUz1atGvp1FgO51VwKFev7aL1XL8N6Jll5oxwUUqxelcOC5/8mu+/sYtff3CQv/83nVc2ZnHhX77mz2uPUtvQ/T/DHeGy+Muobqxm/cn10FADaZ9C8hXg5ZydmnLLavnqSCHXpcTi487rnvfEmUW7dsGJTd2+/PJJQ/DxElbv6tkkJU2DfhzQM8ozCPMLI9y/58Pmymsaue+t3fz4nb2MHRzCf+6Zw45fXsixxy5l44MLuWTcIP765TEu/PPXbE4vNrH2PZMyKIVBQYP4KOMjY6hiYzWMW+q08t/ecRIFLJvuQd0trU2+wVi065s/d/vS8CBfFoyO5v09uTQ1652MtJ7ptwG95YFoT5dsPV5czaV//YbPDuTz04tHs2r5bKYMH0hUiB9eFmFQqD9/vX4Kq5bPwt/Hwm0rd7DexRsDW8TCkoQlbM7dTPG2v0PkKOOBnhM0NVt5e8dJzhsZ5TkPQ9vyCYA59xrPJk52MImrE0unxlJcVc83bvDLX+ub+mVAV0qRXpbe4/7zvPJabnxpG7WNzbz3/Tn8YEESXh2Mp56VEMG/vz+HkdHBLH9tp8t3e7884XKaVTOfVmXCjOXGAz0n+G9aEfkVHvYwtD3T/wcCI2D977t96YIxUYQF+rB6l344qvVMvwzoJXUlVDRU9Cigl1TVc9NL26iobeS1O2Yw2Y5ddsICfXnjrpmMGhTM3S4O6glhCSRbAvkoZICxrreTvLntBDED/Fg0xvP2bj2LX7AxLj3jK8i2c2ORlku9vbhi0hC+OJhPRV2jgyqoebJ+GdCPlh4Fuj/CpaKukVtf3U7O6VpeujWF8UND7b42LNCX1+/8NqhvOOqixckqcrm8OI/Dvl6k1zhnduKJkmrWHy3ie9OHu/cm0GaZfhcERcH633X70qVTY6lvsuox6VqP9IOfrnPtLNyJl3gxIXJCt677vw8PciSvkudvmsbMhIhul9sS1BOjg1n+r1S2Zbpg1/fUV1hcWYWXWPgo8yOnFPnPzSfwEuEmT+9uaeEbZLTSM9fDiS3dunRSbCgJUUG8u1OPdtG6r38G9IKdjAkfQ5CP/bvkHMqt4D+7T3HX/AQW9KLbICzQl3/dOYOhYQHcsXIHu7N7tkpfjzTWQeqrRIy8mLlD5/FJ5idYlWNHVFTVN/Fu6kkumziY6AHOW/zL5VLuhKDobrfSRYTrpw9n54nTHDhV7qDKaZ6q3wX0huYG9hftZ1pM93bm+cNnRwgN8OH7F/R+IlJksB9v/s8sIkP8uPWV7c77wT24GmqKYeYKLk+4nIKaAjae6vlGx/ZYvSuHyvombpsT59By3I5vIMz7EWRtgPQvu3XpddOHEejrxSsbsxxUOc1T9buAfqD4AA3Whm4F9E3pxWw4WsS9C5IIDTBnEk7MAH/euGsmIf4+3PDiVse31Jsb4ZsnIToZ4s9j0fBFDAoaxMv7e7bsqz2sVsXKzceZNCyMKcMHOqwctzX9TmP26Ge/MO6/nUIDfLguZRgf7culsKL7Kzhq/Ve/C+g7C3YCMDV6ql3prVbFHz49wtCwAG6aNcLUusQODOTtu2cxMMiXm17a5tg+9dRXoCQdFj0MIvh4+XDbuNvYVbjrzD0x2zfpxWQWVXN7f2udt/D2g4t/B8VpsOOlbl1625w4mqyK17eecFDlNE/ULwN6UliS3Wugf7w/j/2nynngolH4+5i/32bswEDeXj6bQaH+3Prqdr455oDRL7VlsP4PxiSiURefObx05FLC/cN5cd+L5pcJvLopi6gQPy6d0Af3DDXL6MXG0sT//T1U2z9hKC4yiEVjYnh9WzZ1je6xdITm/vpVQG+yNrGnaI/d3S3NVsWTX6QxdvAArpo8tOsLemhQqD9v3z2buIggbn91B29sM7lV9s0TxhZpFz121kSiAO8Abk6+mU25mzhYctDUIjOLqlifVsSNM4fj692vvs3OJgKX/MFYZuHL33Tr0jvmxVFa3cCHe9puEKZp7bPrJ01ELhGRNBFJF5Gft3P+RhHZZ/vYLCKTzK9q76WdTqO6sdru7pYNx4o4UVLDvQuSHL6zTmSwH2/fPZt5IyP55X8O8ND7+2loMmEEyunjxm46k2+AwRPPOb1s9DJCfEJ4aV/3ugS68tLGLHy9LJ4/M9QeUaNhxt2w6zXI3WP3ZbMTIhgzKIRXNmWherCCo9b/dBnQRcQL+DuwGEgGrheR5DbJsoDzlVITgUeBF8yuqBl25tv6z2PsC+hvbz9JRJCv05Z6DQ3w4eVbp7Pi/ERe35rNTS9vo7y2lzMG1/0fWLxh4UPtng72Deb6sdezLnvdmU2ze6uosp73duZwzbShRIf0o6GKnbngQQiKhI9/CM1Ndl0iItw5L54j+ZWsd9VENK1PsaeFPgNIV0plKqUagFXAla0TKKU2K6VahmlsxdhI2u3sKtxFbHAsg4IGdZm2qLKedYcLuGZarFO7DLwsws8Xj+HpZZPZnX2a+1ftptnaw9ZZ2mdw8D/Gsq4DhnSY7KaxNxHgHcAL+8z5PfzaluM0Nlu5a37/2A3KLv6hcOnjkLsbNj9t92VXTh5K7MAAnvwiDWtPvw+0fsOeSDUUaL0RZY7tWEfuBD5t74SILBeRVBFJLSpybotDKcXOgp12t87/vSuHJqviupRhDq5Z+66cPJSHLx/H+rQinlp3tPsZVOTC+9+HQRNh/gOdJh3oP5AbxtzAmqw1HC453MMaG6rrm3htywkuSo4hMSq4V3l5nHFXQ/JVxgPqgkN2XeLrbeFHF47iwKkKPjuoN5LWOmdPQG+v87jdpoKILMAI6A+2d14p9YJSKkUplRIVFWV/LU2QWZ5JWX0ZKTEpXaZVSvH2jpPMiAsnKdp1QenGmcP5Xsow/vZVOp9354fZ2gz//h9oqofvvmoMn+vCHRPuINQvlL/s/EsvamyseV5e28jd5/ePrf267dInwC8EPrjH7q6Xq6YMZWR0ME98kabXStc6ZU9AzwFaN1NjgXMeu4vIROAl4EqllAsWKelcy1hre0a4bMsqJau4mu9Nd03rvIWI8MiV45gUG8oD7+wlvbDKvgs3PAEnNsJlT0Bkkl2XDPAdwPIJy9mSt4XNuZt7VN/GZisvb8xiRlw4U/vjRCJ7BEfBZU92q+vFyyI8cNFoMouq9b6jWqfsCeg7gJEiEi8ivsAy4MPWCURkOLAauFkp1YP+AcfblreNqIAohoV0HaRXbc8mxN/bLcZP+/t48dxN0/DztnDvm7u6HvmS9Q18/QeY+L1uL4+7bMwyhgYP5amdT/VojZc1+/M4VVbL8vN033mnxl0NyVcaY9NP2Tep6+JxMUyKDeXpdceob9Lj0rX2dRnQlVJNwL3A58Bh4B2l1EERWSEiK2zJfg1EAM+KyB4RSXVYjXugurGaDTkbWDh8YZc7FJXXNLLmQD5XTR5KgK/5E4l6YkhYAI9fO5Ej+ZX87atjHScsOARv3wjhiUYrsJubV/h6+XLvlHs5XHqYNVlrunWt1ap4bn0GSdHBLPT0Nc/NsOQpCBkE795mzBHogojw04vHcKqslte3Zju8elrfZNfwDaXUGqXUKKVUolLqMdux55VSz9s+v0spNVApNdn20XVHtRN9lf0Vdc11XJZwWZdp30k9SUOT1eXdLW0tHBPD0qlDeXZ9Bvtz2lnMqywbXl8KPoFw82qjn7YHLo2/lLHhY/nbrr9R12T/OiKfH8znSH4l9y10/Jh9jxAYDteuhIo8eP8esGOc+byRkcwfGclTa4+SX67XeNHO1S+m8H2S+QlDg4cyOWpyp+nqGpt58ZtMZidEdGvzCmd5eMk4IoN9+cm7e8/+s7u6BP61FBpq4KZ/Q1jPJ/NYxMIDKQ+QW53LywfsW7jLalU8/eUxEqKCWDKx4+GRWhuxKXDRo5C2BrY8Y9clj145noZmK7/64ICebKSdw+MDenFtMVvztnJp/KVddre8tzOHwsp67l1o34NEZwsN9OH3SyeQVlDJ375MNw7WnoY3rzVa6DesgphxvS5n5uCZXBp/KS/vf5ms8q6XcP3ikNE6/9+FIzvcW1XrwMwVMPYKWPswHN/UZfK4yCB+/J1RrD1UwJr9ehijdjaPD+ifH/+cZtXMpfGXdpqusdnK819nMHlYGHMSu78bkbMsHBPDNVNjee7rDA6mHYWVSyB/P1z3Txgxx7Ryfjr9p/h7+fPY1sc6bQlarYqn1h0jITKIyyfp1nm3icCVz0BEIqy6Hgq7ngdw57x4JgwN5eEPD1BW0+CESmp9hccH9DWZaxg9cDRJAztvdX+4J5ec07XcuyCpy5a8q/368mQmBpcTuuoKVGkW3PCOsaqfiSIDIrl/6v1sy9/GJ1mfdJjui0MFRt/5oiTdOu8p/1Cjq8w7AF6/Bso7337O28vCH66ZwOmaRn77Se8mgmmexaMDenZFNvuK93X5MNRqVTy7Pp0xg0L6xAiN0PI0Vnk9TIi1nFcT/wKJCxxSzrWjr2Vi5EQe3/E45fXnPog903ceGcTluu+8d8KGw03vQX2lEdS7GPkybkgod5+XwHs7c3g39WSnabX+w6MD+pqsNQjC4vjOW6+fHcwno6iae5ywqmKv7X0bXroQPy/hX2Oe49G9wWx10MYYFrHwq9m/oqy+jIc3P3zO2PT395zicF4F9y1KwtvLo7+VnGPQBFj2BpRmwhvXQk1pp8l/9J1RzE2K4Ber97M53f611jXP5bE/hUopPsn8hGkx0zpdjKu2oZk/rz1KfGQQl7nBRKIONTXAmp/Cf5bD0Klw9wbuuGYJw8MD+cm7e6mqt28aeXeNCR/DT1J+wpfZX561LEBlXSO/W3OEScPCuHKS49aK73fizzOWa8jbC68u7rT7xcfLwrM3TiM+Moi7X99JemGlEyuquSOPDeifn/ic4xXHuTLpyk7T/eqDA2QUVfHIFePctw+44CC8chFsfwFm3wu3fAAhMQT6evPktZM4VVbLQ//Z77BhbDeNvYllo5ex8uBK3kl7B4Cn1x2jpLqeR68c5/5/1fQ1Y5fATauNBdZevgiK0jpMGhrgwyu3TcfP24vbXt1BUWW9EyuquRuPDOhVDVX8afufGBs+lssTLu8w3TupJ3lvZw73LUjivFHOXSzMLk318NVv4R/nGcMSr3sNLn4MvL7dqDolLpwfXTiK9/fk8uqm4w6phojw4IwHmT90Pr/b9jvePrCWVzcfZ9n04UyMDXNImf1e/Hy4fQ1Ym+CViyF9XYdJh4UH8vKtKRRX1bP0uU0czqtwYkU1d+KRAf3Zvc9SXFvMQ7MewsvS/vT9w3kV/Or9A8xNiuD+C0c5uYZdaG6C/e/B8/Ngw+Mw/rvwgx3G+h/tuHdBEhclx/DYmsNsznBMX6q3xZvHz3+cpLAkHkt9kKDI7fzkIje7b55m0AS443MIGWw8KF37MDS3v+HJpGFhrFo+m4YmK0uf3cyn+/OcXFnNHXhcQE8rTePNw2/y3VHfZWLUuVuuAZTVNPCDN3YRGuDDU9+b4j5dLfWVsP1F+NtU+PedxnTwG/8NS/8BQR2PjbdYhCevm0RcRCD3vrmbnNM1DqlekE8QS4c8QmN1PCri3/xm+8/aHf2imSg8Hu76EqbdBpueglcvhdPt7zk7eVgYH907jzGDQ/j+G7v402dH9EJe/Yy4avpwSkqKSk01dw0vq7Jy66e3cqLiBB9d/RGhfudO388tq+WWV7aTXVLDv+6cwcwEF04iam6EwkOQuR6OrYXsrWBthNjpMPeHMPpSsNj/OzejqIqrntnEsPBAXr9rJuFBvqZW91hBJVc/u5mE6ECWnp/J07ufIsI/gp9N/xnfGfEdtx+/3+cd+Dd8eL/RDXPeAzDnf9td676+qZmHPzjIqh0nSYgK4rdXjWdOYqQLKqw5gojs7Gi9LI8K6C/tf4mndz3No3Mf5aqkq845f7Sgklte3k51fRMv3prCLGcFc6sVyrOh+JjxgKvoMOTtg6Ij0Gyb6RczHpIWwZglRkDvYXD8+mgRy19LZejAAF67YwaxAwNN+RLKaxq58u8bqapv5qP75jI4NICDJQd5aONDpJelMy5iHPdPvZ/ZQ2abUp7WgbKT8Pn/g8MfQngCXPx7GHVxu98v69MK+fUHB8kureHqKUP56cWjGRIW4IJKa2by+IBuVVb+svMvrDy4kotGXMTj5z+ORc5u2W48Vsw9b+wkwNeLlbfPYOzgAaaUfRaloOKUEawLDhoBuzgNitOhqfbbdIERxtZwgybA4EnGlP1O9vzsrh3HS7lz5Q78fbz45x29/1qbrYrbV+5gS0Yxb/3PLFLiwluda+bjzI/5+56/k1edx9ToqdycfDMXDLsAb4t3b78UrSPpX8KnD0LJMaMxMOsemPDdc1rsdY3NPPvfdJ7/OhOAa1Ni+f4Fiab9otecz6MDemNzIw9teog1WWtYNnoZP5/x87MehBZX1fO7NYdZvesUiVFB/NOsVmtTA5RmGOuo5O8zXvP2QW2rySBhwyFyNESNhshR33500h9ulrT8Sm59ZTvVDU388ZqJLB4/qEddIkopHvvkMC9tzOL3Sydw/Yz2V3JsaG7g3aPv8trB18itzmVw0GCuG30dl8VfxuBgNx7f35c1NcC+t2Hrs0bXXVAUTLjW6KobPhu8vv2FmnO6hue/zuCdHTlYlWLJxMFcNWUo85Ii9aSwPsYjA7pVWdl4aiMv7HuBvUV7uX/q/dw5/s4zQau6vonVu3J44ouj1DQ0cfd5ifxgQVL3N61QyhgyWHAA8g9A4UEoPGIEc6ttMo+XH8QkGy3uQRONj5hx4OfaTZJPldWy/LVUDuZWcP6oKH5z5ThGRATZfX1JVT0P/ns/6w4XcPOsETx61fgur2m2NrM+Zz1vHn6T7fnbAZgYOZGL4i7ivNjziBsQp/vazaaU8Rxm+wtGy725HvzDIHGhEdiHz4ToceDlTV55Lf/4OpPVu3KoqGsiIsiXSycMZm5SBNNGhBMV0vX+s5pr9Tqgi8glwNOAF/CSUuoPbc6L7fylQA1wm1JqV2d59jSgVzZU8n76+7x15C1OVp4kKiCKH6f8mCUJS6htaGZrZgn/2X2KtYcKqG1sZlZCOL+9ajxJ0XZs+NDcCKePG62d3N3fftS1jOQQY9RB1FiIHgNRY4w/dyNHnjU23J00NVt5bcsJ/rz2KA3NVm6fE8c102IZFdP5/dhwtIgH3t1LeU0jDy4ew+1z4ro9gSi7IpsvTnzBF8e/4HCpsYhUuH8402KmMSlqEiPDRpIQlkBMYIwO8mapr4LM/8KRNUaQr7Rt/+sbbHyvDhoPMeNpCB/FltJg3j3SyLq0IuoajWUd4iODGDdkACOjQ0iKDiYhKoghYQEM8PfW/0ZuolcBXUS8gKPAdzA2jN4BXK+UOtQqzaXAfRgBfSbwtFJqZmf59jSgv3vkP/xm269JCBnHtLDLibCkkFFYy8FTFRwrrMSqjNlzl00YxNUTo0kZ6oc01EBjDdRVQN1pqC0zFj+qKoSqfKgsgNNZxhoaLa1ui7fRyh4yxejnjplgtMJ97W/hupOCijoe++Qwn+zPo9mqGDMohMsmDCYuMojIYD8ig33JOV3L1swStmaWsDennJHRwTy9bArJQ3r/vOFkxUm2529nZ8FOdhXu4lTVt5sdB3oHMjhoMNGB0UQHRhMREMEA3wEM8BtAiG8Igd6B+Hv54+/tj5+XHz4WH3wsPnhbvLGIBS+LFxaxYMGCiCAiWLAYx9r56DeUgvKTkL0NTm4zugULDkJDqyUCvP1RocOo8gmn0DqAE3WBZNf4cqrOl3IVSJUKoAY/lHcgQSEDCAwMIjgwmKCgQIIDAwkMCCQwwI+gwEAC/PwI8PMm0NcLfx8v/Lwt+Hpb8PP2wttL8PWy4G0RvCyifzn0Qm8D+mzg/5RSF9ve/wJAKfX7Vmn+AaxXSr1le58GXKCU6nB2Q08D+jefvU596qOMamjGIgpB4SMKP4sVX4vCR5rxsjYizfZMgRYIjobgGFt/t62PO2qU8Seqj3+36+fuiqvqWbM/jw/25LLzxLkr+vl4CVOGDeS8UZHcNT8Bfx/H7KtaUltCZnkmWeVZZJZnkl+dT2FNIQU1BZTWltKkHLM2DYCXeCEIxv9yJtC3fC4iZ/2CaH1eRDjzX6ugdNY523Gh66DVWWCz5/ruMPJTxl+izQ3Gq7XR9toE1mbjVVmNdD2gbHX+9vXb0tvPsaPj357vrdb5O+LXSE/u1Cz/Sfz65jd6VF5nAd2eYQhDgdbrc+ZgtMK7SjMUOCugi8hyYDnA8OE92yZtbFw89aem4+fjjZ+vD/4+Xvj4+IB4Ga1qizd4+xr92t6+xh6bvkHgE2T0afuHQcBACAiDwMizHhz1B5HBftwyO45bZsdRXtNIQWUdxZX1FFXVExXsx5ThA52yOXZEQAQRARFMHzT9nHNKKWqbaqlsqKSioYK6pjrqmuuobaqlobmBRmsjjdZGmqxNNKtmrFYrzaoZhUIphVVZURivLR/NqhmllJHedl4pxZn/bOdajrVc19E1LfUEzsqj5dy3Lx3/uHd6zuRnW52V1U5iI9Bbm4xgr5rPDvbWZuNVWbFam2luNl6V1bhfytrqXikFymqU3vL+TCG2Y7bXcyvR+m1P74d7btM30N8xy3TbE83a+6XW9i7Zkwal1AvAC2C00O0o+xyRY+bCmLk9uVRrIzTQh9BAny77051NRAj0CSTQJ5CYoBhXV0fT+gx7OhRzgGGt3scCuT1Io2mapjmQPQF9BzBSROJFxBdYBnzYJs2HwC1imAWUd9Z/rmmappmvyy4XpVSTiNwLfI4xbPEVpdRBEVlhO/88sAZjhEs6xrDF2x1XZU3TNK09dj0RVEqtwQjarY893+pzBfzA3KppmqZp3dGPBuVqmqZ5Nh3QNU3TPIQO6JqmaR5CB3RN0zQP4bLVFkWkCGh/Ly3HiQQcs+lm36Tvx7n0PTmbvh/ncvU9GaGUandXe5cFdFcQkdSO1kDoj/T9OJe+J2fT9+Nc7nxPdJeLpmmah9ABXdM0zUP0t4D+gqsr4Gb0/TiXvidn0/fjXG57T/pVH7qmaZon628tdE3TNI+lA7qmaZqH8LiALiJeIrJbRD62vQ8XkbUicsz2OrBV2l+ISLqIpInIxa6rteOIyHER2S8ie0Qk1Xas394TEQkTkfdE5IiIHBaR2f38foy2fW+0fFSIyA/7+T35kYgcFJEDIvKWiPj3mfvRsk2Up3wAPwbeBD62vf8T8HPb5z8H/mj7PBnYC/gB8UAG4OXq+jvgfhwHItsc67f3BPgncJftc18grD/fjzb3xgvIB0b013uCsXVmFhBge/8OcFtfuR8e1UIXkVjgMuClVoevxPghxvZ6Vavjq5RS9UqpLIy13Gc4qaqu1i/viYgMAM4DXgZQSjUopcrop/ejHYuADKXUCfr3PfEGAkTEGwjE2H2tT9wPjwrowFPAzwBrq2MxyrZ7ku21ZXfWjja29jQK+EJEdto26Yb+e08SgCLgVVu33EsiEkT/vR9tLQPesn3eL++JUuoU8ASQjbHJfblS6gv6yP3wmIAuIkuAQqXUTnsvaeeYJ47hnKuUmgosBn4gIud1ktbT74k3MBV4Tik1BajG+PO5I55+P86wbS95BfBuV0nbOeYx98TWN34lRvfJECBIRG7q7JJ2jrnsfnhMQAfmAleIyHFgFbBQRF4HCkRkMIDttdCWvl9sbK2UyrW9FgL/wfhzsL/ekxwgRym1zfb+PYwA31/vR2uLgV1KqQLb+/56Ty4EspRSRUqpRmA1MIc+cj88JqArpX6hlIpVSsVh/On4lVLqJowNrG+1JbsV+MD2+YfAMhHxE5F4YCSw3cnVdigRCRKRkJbPgYuAA/TTe6KUygdOisho26FFwCH66f1o43q+7W6B/ntPsoFZIhIoIoLxPXKYvnI/XP1U2UFPqi/g21EuEcCXwDHba3irdL/EeCqdBix2db0dcB8SMJ7A7wUOAr/U94TJQCqwD3gfGNif74ftawwESoDQVsf67T0BHgGOYDR+/oUxgqVP3A899V/TNM1DeEyXi6ZpWn+nA7qmaZqH0AFd0zTNQ+iArmma5iF0QNc0TfMQOqBrmo2IrBcRUzf/ta3ueE+r9xe0rASqaWbTAV3THCsMuKerRJpmBh3QNbcnIj8Tkf+1ff4XEfnK9vkiEXldRJ4TkVTbGtaP2M4tFpF3WuVxgYh8ZPv8IhHZIiK7RORdEQlup8x204ixvvwjtuP7RWSM7XiUbZ3sXSLyDxE5ISKRwB+ARNta44/bsg+Wb9dkf8M2I1HTek0HdK0v2ADMt32eghEQfYB5wDcYM2BTgInA+SIyEViLMYU7yHbd94C3bUH2IeBCZSxaloqxhv4ZdqQpth1/DviJ7djDGMtNTMVYM2e47fjPMZaknayU+qnt2BTghxhraSdgrEOkab2mA7rWF+wEptnWpakHtmAE9vkYAf06EdkF7AbGAclKqSbgM+By27rWl2GsvzELI5BuEpE9GOtyjGhTXldpVreqV5zt83kYi8KhlPoMON3J17NdKZWjlLICe1rloWm94u3qCmhaV5RSjbZVNG8HNmOsw7IASARqMVrJ05VSp0VkJeBvu/Rt4AdAKbBDKVVp695Yq5S6vpMiu0pTb3tt5tufoe50m9S3+rx1HprWK7qFrvUVGzAC9waMVvkKjNbtAIx1zctFJAZjGdgW6zGWx/0fjOAOsBWYKyJJALZV9Ua1KcueNG1tBK6zpb8IY9EvgEogpDtfqKb1lA7oWl/xDTAY2KKMNbvrgG+UUnsxuloOAq8Am1ouUEo1Ax9jBPmPbceKMPaIfEtE9mEE7zGtC7InTTseAS6ydf0sxtjtplIpVYLRdXOg1UNRTXMIvdqipplARPyAZqVUk4jMxtgVabKLq6X1M7rvTtPMMRx4R0QsQANGN4+mOZVuoWuapnkI3YeuaZrmIXRA1zRN8xA6oGuapnkIHdA1TdM8hA7omqZpHuL/AzpYpRayeHxPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cmf.plot(x='wavelength', y=['x','y','z'], label=['x','y','z'])\n",
    "plt.legend()\n",
    "plt.title('Colour matching functions')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = pd.read_csv(sep = \"\\t\", filepath_or_buffer = \"data/test.csv\", header = None)\n",
    "# Renaming the columns of the test dataframe\n",
    "test.columns = [\"wavelength\", \"y1\", \"y2\", \"y3\", \"y4\"]\n",
    "#test['y1'] = np.cos(test['y1']/0.0175)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**From Wernand et al (2010):**\n",
    "\n",
    "`\"The chromaticity coordinates, based on transmission measurements, of the FU scale and the basic solutions. The white refers to the coordinates equal to 1/3\"`\n",
    "\n",
    "The transmission measurements used a TrIOS spectroradiometer with 4 nm band resolution.\n",
    "\n",
    "**From Wernand et al (2013):**\n",
    "\n",
    "`\"Because MERIS does not provide full-spectral range coverage, the reflection spectrum is first reconstructed by linear interpolation between band n=1 (412.5 nm) and band n=9 (708) with a resolution of 1 nm\"`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "#the interpolating part \n",
    "Delta_lambda = 5 # 1nm bins\n",
    "# adds +4 as it should include 720 resulting in 86 values, i.e. matching the len(cmf).\n",
    "iwavelength = np.array(range(380, 720 + Delta_lambda, Delta_lambda)) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# retrieving the reflectance values for lambda 380-720nm\n",
    "irho = interp1d(test[\"wavelength\"], test[\"y1\"], kind = 'linear')(iwavelength)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1d2ab260ac0>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAy5klEQVR4nO3deXzU5bn//9dFiBJXVHALKuhxBwwYdkWwKqIVLLYFfu1PtMcv9VTaU/stLba1WlpbWtqjxy5y3GoXK1AXRKV1qaAIciRIREWpgCgBiogGQQJmub5/zGfiJ8PMZJLMJLO8n4/HPDLz2eb+kDDX3Pd1L+buiIiIpEOnji6AiIjkDwUVERFJGwUVERFJGwUVERFJGwUVERFJGwUVERFJGwUVERFJGwUVKShmtiv0aDCzmtDrL7XieovM7JpMlLU1zGyEmVV1dDmkcHXu6AKItCd3Pyj63Mw2ANe4+zMdV6L2Z2ad3b2uo8sh+Uk1FRHAzDqZ2TQzW2dm281srpkdHuzrYmZ/DrZXm9lyMzvKzG4BzgV+E9R0fhPnunHPDfYtMrOfmdlLZrbDzB6Nvmewf7CZLQ3Oe8XMRoT2HW5mvzezzWb2oZnNM7MDgb8Bx4ZqX8ea2c1m9mBQjo+Aq8zsPjP7Seh6TWo4ZrbBzKaa2Soz+9jM7gnu+W9mttPMnjGzw9L+i5Ccp6AiEvEN4HLgPOBY4EPgt8G+ScChwHHAEcC1QI27fx9YDExx94PcfUqc68Y9N7T/SuArwXvWAbcDmFkp8ATwE+Bw4NvAQ2bWPTjvT8ABwJnAkcCt7v4xMBrYHJTnIHffHBw/FngQ6Arcn+K/yRXAhcApwGVEAtb3gG5EPju+keJ1pICo+Usk4qtEgkMVgJndDLxrZv8/UEskIPybu68CVrTgus2d+yd3fy14zxuBSjObBHwZWODuC4LjnjazCuASM3uKSPA4wt0/DPY/10w5XnT3ecHzGjNLpey/dvetQdkWA++5+8rg9SPAZ1K5iBQW1VREIk4AHgmamqqBN4B64CgitYIngdlBc9MvzKw4xes2d+7G0PN3gGIiNYETgC9EyxOU6RzgGCK1ng9CASUVG5s/ZB9bQ89r4rw+CJEYCioiERuB0e7eNfTo4u6b3L3W3X/k7mcAQ4HPEmm2Akg6zXcz50IkQEQdT6Rm835Qnj/FlOdAd58R7DvczLrGe8tERYl5/TGR5rOoo5Pdh0iqFFREImYBt5jZCQBm1t3MxgbPR5pZHzMrAj4i8sFfH5y3FTgx0UWbORfgy2Z2hpkdAEwHHnT3euDPwGVmNsrMioKE/wgz6+HuW4jkN35nZoeZWbGZDQ+V5wgzO7SZ+60k0pR2uJkdDXwzpX8lkWYoqIhE/DcwH3jKzHYCy4BBwb6jiSS5PyLSLPYckQ/96HmfD3pg3R7nusnOhUjz2H3Av4AuBMlvd99IJLn+PWAbkdrJVD79PxvN9bwJvEcQFNz9TeABYH3QbHZsgvv9E/AKsAF4CpiT+J9GJHWmRbpEOoaZLQL+7O53d3RZRNJFNRUREUkbBRUREUkbNX+JiEjaqKYiIiJpU9Aj6rt16+Y9e/bs6GKIiOSUFStWvO/u3ePtK+ig0rNnTyoqKjq6GCIiOcXM3km0T81fIiKSNgoqIiKSNgoqIiKSNgoqIiKSNgoqIiKSNgXd+0sk0+at3MTMJ9ewubqGY7uWMHXUqVzer7SjiyWSMQoqIhkyb+Umbnj4VWpqIzPdb6qu4YaHXwVQYJG8paAikmbR2smm6pp99tXU1jPzyTUKKpK3FFRE0ii2dhLP5jjBRiRfKKiIpEGy2kksB4bNeFb5FclLCioibZRK7SSW8iuSr9SlWKSNZj65pkUBJSqaXxHJJ6qpiLRSqk1eJcVFCYOO8iuSb1RTEWmFaJNXcwGltGsJPxvXh9KuJXH3dzKj17QnGDbjWeat3JSJooq0q4wGFTO72MzWmNlaM5sWZ7+Z2e3B/lVm1j+0714ze8/MXos5Z46ZVQaPDWZWGWzvaWY1oX2zMnlvUtiaa/IqKS7itvFlLJl2Ppf3K2XqqFMpKS7a57h6d5xPcywKLJLrMtb8ZWZFwG+BC4EqYLmZzXf31aHDRgMnB49BwB3BT4D7gN8Afwxf193Hh97jV8CO0O517l6W1hsRiSNZs1VpnJHz0efR0fWdzKiPWcpbY1gkH2QypzIQWOvu6wHMbDYwFggHlbHAH93dgWVm1tXMjnH3Le7+vJn1THRxMzPgi8D5GbsDkQSO7VoSt+mrtGsJS6bF/5O8vF9pY8DoNe2JuMcoxyK5LpPNX6XAxtDrqmBbS49J5Fxgq7u/FdrWy8xWmtlzZnZuvJPMbLKZVZhZxbZt21J8K5GIeSs3MWzGs2yqrsFi9pUUFzF11KkpXefYBDmWRNtFckUmg0rs/zmIjPtq6TGJTAQeCL3eAhzv7v2AbwF/MbND9rm4+53uXu7u5d27x11iWSSu2OS88+kfcDQhn2rTVbwcixHJrShpL7ksk81fVcBxodc9gM2tOGYfZtYZGAecHd3m7nuBvcHzFWa2DjgF0CL0khbxkvNO8iavRMI5lmitJ/ptSgMjJZdlsqayHDjZzHqZ2X7ABGB+zDHzgSuDXmCDgR3uviWFa18AvOnuVdENZtY96ByAmZ1IJPm/Ph03IgKJ8x2tzYNc3q+UJdPOp7RryT7Vcw2MlFyVsaDi7nXAFOBJ4A1grru/bmbXmtm1wWELiHzwrwXuAr4WPd/MHgBeBE41syoz+/fQ5SfQtOkLYDiwysxeAR4ErnX3DzJwa1KgMpUHSXewEulIGR1R7+4LiASO8LZZoecOXJfg3IlJrntVnG0PAQ+1tqwiiYRHzoebqaBlyflEEvUka2uw0gJh0hHMPdW8eP4pLy/3igqlXCSxeJNFRgNLvPEo2fIezV1z5GndWfjmNgUcaRUzW+Hu5fH2ae4vkSTSmZxPJJ1J+2TzkYWv+edl7zZuV8cASSfVVFRTkSR6TXsibh93A96ecWna3y86BiaecK0lXtMW0OIp+JO9h0giqqmItFKm8h2JJEvOb6qu4fo5lXxzTuU+tZnr51SmPMArGdVapK00S7FIEvEGKaYjOZ9Ic8HKY37Gbk8HdWeWtlBQEYkjOh3L9XMq2b9zJw47oBij5SPnWyrRbMZtFW/qimTUnVlaS81fIjFie05V19RSUlzErePLMt4kFJu0b6uS4iJ+Nq5P4zWjOZho769E76E5yKS1lKhXol5iJEqWp7PHVyridQtuiVSS7u3RZVryjxL1Ii2QLSPck3U1hk8//OMNyEy1iU5zkEm6KaciEiObpqWPzg+2Ycal3Dq+jNKuJY25nVvHl8Xd3tKcj+Ygk3RSTUUkxtRRp+7TJJTJHl+pCi/ylcr2lsqWGprkNtVURGJc3q+Un43r06Zv/7kom2pokrtUUxEJFPoEjNlaQ5PcopqKCE1XdXQ+TVIX0gqMsTW0riXFdCnuxPVzKrUapaRMQUWE+BNHFmKSOpq0v3V8GXvrGvhwd23BBllpHQUVEZSkjqUgK62loCKCktSxFGSltRRURGj/iSOznYKstJZ6f0nBiu3tdcXZpVoNMaCeYNJamvtLc38VpHhzXrVkepNCEA66h5YUYwbVu2sLPuBK8rm/Mtr8ZWYXm9kaM1trZtPi7Dczuz3Yv8rM+of23Wtm75nZazHn3Gxmm8ysMnhcEtp3Q3CtNWY2KpP3JrlNiejmqSeYtEbGgoqZFQG/BUYDZwATzeyMmMNGAycHj8nAHaF99wEXJ7j8re5eFjwWBO93BjABODM473dBGUT2oUR06hSApSUyWVMZCKx19/Xu/gkwGxgbc8xY4I8esQzoambHALj788AHLXi/scBsd9/r7m8Da4MyiOxDiejUJQq0m6prNChS9pHJoFIKbAy9rgq2tfSYeKYEzWX3mtlhLbmWmU02swozq9i2bVsKbyX5SL29Upcs0IabwqKrZfaa9oSCTQHLZFCJt4JpbK+AVI6JdQdwElAGbAF+1ZJrufud7l7u7uXdu3dv5q0k33TUMsG5rLkljmtq6/nmnEqun1PZZJqb6+dU0lMBpuBksktxFXBc6HUPYHMrjmnC3bdGn5vZXcDjrb2WFJaOXCY4l6W6xHHsNzgt9lWYMllTWQ6cbGa9zGw/Ikn0+THHzAeuDHqBDQZ2uPuWZBeN5lwCnwOivcPmAxPMbH8z60Uk+f9SOm5E8oMSzq0XXsirNWpq6/m/c19R01gByFhQcfc6YArwJPAGMNfdXzeza83s2uCwBcB6Ikn1u4CvRc83sweAF4FTzazKzP492PULM3vVzFYBI4Hrg/d7HZgLrAb+Dlzn7q1b3Fvyknp8tV1zTWHJ1LurS3IB0OBHDX4sGMNmPBu3+aa0awlLpp3fASXKTdFBkcmawlJVqoGUOanDBj+KZINocn5Tdc0+vTnU46vlok1ht40v26fWYjE/m6NaS/5RUJG8Fl58CyLJ4+gHnnp8tU28ZZdvHV/GhhmXcuv4ssbtRZY8xCivlV80oaTktXjJeUdNXulyeb/SuEE5vD3ePGuxlNfKHwoqkteUnO94qXRJ1kwG+UPNX5LXNB1LdkiWh1FeK78oqEhe03Qs2SVeHkZ5rfyiLsXqUpyXtBZIbohdKE2/m9yQrEuxciqSdzQdS26I/T1pOpf8oOYvyTuajiU36PeUnxRUJO+ox1du0O8pPymoSN5Rj6/coN9TflJQkbyjHl+5Qb+n/KREveSd8GA79SrKXvo95ScFFckb6p6aexJN8yK5S0FF8oK6p4pkB+VUJC+oe6pIdlBQkbyg7qki2UHNX5IXju1aEncGXHVPzR3KieUH1VQkL6h7am4LL6amdexzm2oqktNiJ47sUtxJE0fmoGQ5Mf0Oc4uCiuQsTRyZP5QTyx8Zbf4ys4vNbI2ZrTWzaXH2m5ndHuxfZWb9Q/vuNbP3zOy1mHNmmtmbwfGPmFnXYHtPM6sxs8rgMSuT9yYdTz2+8oembMkfGQsqZlYE/BYYDZwBTDSzM2IOGw2cHDwmA3eE9t0HXBzn0k8Dvd29L/BP4IbQvnXuXhY8rk3LjUjW0rfb/KGcWP7IZE1lILDW3de7+yfAbGBszDFjgT96xDKgq5kdA+DuzwMfxF7U3Z9y97rg5TKgR8buQLKavt3mD60ImT8ymVMpBTaGXlcBg1I4phTYkuJ7fAWYE3rdy8xWAh8BP3D3xbEnmNlkIrUijj/++BTfRrLR1FGnNsmpgL7d5jJN2ZIfMhlULM622LWLUzkm/sXNvg/UAfcHm7YAx7v7djM7G5hnZme6+0dNLu5+J3AnRJYTTuW9JLuox5dI9spkUKkCjgu97gFsbsUx+zCzScBngc+4uwO4+15gb/B8hZmtA04BtAh9HlGPL5HslsmcynLgZDPrZWb7AROA+THHzAeuDHqBDQZ2uHvSpi8zuxj4LjDG3XeHtncPOgdgZicSSf6vT9/tSDZQjy+R7Jaxmoq715nZFOBJoAi4191fN7Nrg/2zgAXAJcBaYDdwdfR8M3sAGAF0M7Mq4CZ3vwf4DbA/8LSZASwLenoNB6abWR1QD1zr7vsk+iW3qceXSHbL6OBHd19AJHCEt80KPXfgugTnTkyw/d8SbH8IeKjVhZWcoDm+CoPmActdGlEvOSH6IbOpugajaW8O9fjKL1obJ7dpQknJeuHJBiESUKLdBjWeIf8ob5bbVFORrBfvQ8aJBJQl087vmEJJxihvlttUU5Gspw+ZwqKZEnKbgopkPX3IFBbNA5bb1PwlWUvJ+cIUzY+p91duUlCRrBTbAyianI/mUvQhk980D1juUlCRrKTkvEhuUk5FspKS8yK5SUFFspKS8yK5Sc1fklWUnBfJbQoqkjWUnBfJfQoqkjWUnBfJfcqpSNZQcl4k9ymoSNZQcl4k9ymoSNbQ9BwSz7yVmxg241l6TXuCYTOeZd7KTR1dJElCORXpULGLMV1xdikL39ym6TkE0NoquUhBRTpMvA+Mh1Zs0voo0ijZ2ir6G8lOav6SDqPFmKQ56ryRexRUpMPoA0Oao84buUdBRTqMPjCkOeq8kXsyGlTM7GIzW2Nma81sWpz9Zma3B/tXmVn/0L57zew9M3st5pzDzexpM3sr+HlYaN8NwbXWmNmoTN6btJ0+MKQ5l/cr5Wfj+lDatQQjMhBWObfsZu7e/FGxJ5ld6O5PN3NMEfBP4EKgClgOTHT31aFjLgG+DlwCDAL+290HBfuGA7uAP7p779A5vwA+cPcZQaA6zN2/a2ZnAA8AA4FjgWeAU9y9aaN9SHl5uVdUVLT4/qVtwj2+Di0pxgyqd9eqt5dIjjCzFe5eHm9fa3t/3QMc38wxA4G17r4+KMRsYCywOnTMWCJBw4FlZtbVzI5x9y3u/ryZ9Yxz3bHAiOD5H4BFwHeD7bPdfS/wtpmtDcrwYivuTzIktsdXdU0tJcVF3Dq+TMFEJA8kDCpmNj/RLuCIFK5dCmwMva4iUhtp7phSYEuS6x7l7lsA3H2LmR0ZutayONdqWnizycBkgOOPby4uSrqpi6hIfktWUzkX+DKRJqgwI1IDaI7F2Rbb1pbKMalK6VrufidwJ0Sav1r5XtJK6vElkt+SBZVlwG53fy52h5mlMpCgCjgu9LoHsLkVx8TaGm0iM7NjgPfacC1pZ8d2LWFTnACiHl8i+SFh7y93H+3uC81sSriHVbBveArXXg6cbGa9zGw/YAIQ26Q2H7gy6AU2GNgRbdpKYj4wKXg+CXg0tH2Cme1vZr2Ak4GXUiintIPo/E3RxbfC1ONLJH+kkqg/GlhuZi8D9wJPegpdxty9zsymAE8CRcC97v66mV0b7J8FLCDS82stsBu4Onq+mT1AJCHfzcyqgJvc/R5gBjDXzP4deBf4QnC9181sLpGOAHXAdcl6fkn70eJbIoUjpS7FZmbARUQ+9MuBucA97r4us8XLLHUpbh/RGkosLb4lkpuSdSlOafBjUDP5V/CoAw4DHgzGjIgkpeS8SOFoNqiY2TfMbAXwC2AJ0Mfd/wM4G7giw+WTPKDpWEQKRyo5lW7AOHd/J7zR3RvM7LOZKZbkg+jI+WhyPtzQquS8SH5qNqi4+w+T7HsjvcWRfKHkvEhh0iJdkhHxRs5HA4qS89JasSuF6stJ9lFQkYxQcl7STUsL5watpyIZoeS8pJtWCs0NCiqSEVorRdJNtd/coOYvSavYtVK6FHfSWimSFpo3LjcoqEjaaK0UyaSpo05t8vcFqv1mIzV/SdqozVsySUsL5wbVVCRt1OYtmXZ5v1IFkSynmoqkjXp8iYhqKtJmmo5FRKIUVKRNNB2LiIQpqEibaDoWEQlTTkXaRMl5EQlTUJE2UXJeRMLU/CWtouS8iMSjoCItpuS8iCSS0eYvM7vYzNaY2VozmxZnv5nZ7cH+VWbWv7lzzWyOmVUGjw1mVhls72lmNaF9szJ5b4WsueS8AopI4cpYTcXMioDfAhcCVcByM5vv7qtDh40GTg4eg4A7gEHJznX38aH3+BWwI3S9de5elql7kggl50UkkUzWVAYCa919vbt/AswGxsYcMxb4o0csA7qa2TGpnGtmBnwReCCD9yBxKDkv2aKhoYGtW7fi7s0fLO0ik0GlFNgYel0VbEvlmFTOPRfY6u5vhbb1MrOVZvacmZ0br1BmNtnMKsysYtu2banfjTBv5SaGzXi2MTkfpuS8tLeGhgZGjhxJjx49GDFiBA0NDR1dJCGzQSX2cweadhJKdkwq506kaS1lC3C8u/cDvgX8xcwO2eci7ne6e7m7l3fv3j1h4aWpaHI+up5F+Jek2WKlI2zbto2lS5dSV1fH0qVL0ZfE7JDJ3l9VwHGh1z2AzSkes1+yc82sMzAOODu6zd33AnuD5yvMbB1wClDR1hsRjZyX7HPkkUcydOhQli5dytChQznyyCM7ukhCZoPKcuBkM+sFbAImAP9fzDHzgSlmNptIon6Hu28xs23NnHsB8Ka7V0U3mFl34AN3rzezE4kk/9dn6N4KjpLzkm3MjIULF7Jt2zaOPPJIImlW6WgZCyruXmdmU4AngSLgXnd/3cyuDfbPAhYAlwBrgd3A1cnODV1+Avsm6IcD082sDqgHrnX3DzJ1f4VGS7lKNurUqRNHHXVURxdDQqyQe02Ul5d7RYVax5JpbuS8cikihcfMVrh7ebx9GlEvCWnkvIi0lIKKJKTkvIi0lGYploSUnBeRllJQkYQ0cl5EWkrNX9JENDG/ubqGQ0uKKS4yaus/Tc9r5LyIJKOgIo1iE/PVNbUUdzIOO6CY6t21HKvkvIg0Q0FFGsVLzNc2OAfs15mVP7yog0olIrlEORVppMS8iLSVgoo0UmJeRNpKQUUaTR11KiXFRU22KTEvIi2hnIrs0+OrS3EnJeYl5zQ0NGhyySygoFLg4vX4Kiku4tbxZQomkjOiC3ZFp8FfuHAhnTqpIaYj6F+9wMXr8VVTW8/MJ9d0UIlEWk4LdmUPBZUCpx5fkg+iC3Z17txZC3Z1MDV/FahoHiXRwgfq8SW5RAt2ZQ8FlQIUm0eJpR5fkou0YFd2UFApQPHyKFFaJ0VE2kJBpQAlypcYaJ0UEWkTJeoLkEbOi0imqKZSQJpbb155FBFpKwWVAqH15kWkPWS0+cvMLjazNWa21symxdlvZnZ7sH+VmfVv7lwzu9nMNplZZfC4JLTvhuD4NWY2KpP3lmuaW29eAUVE0iFjQcXMioDfAqOBM4CJZnZGzGGjgZODx2TgjhTPvdXdy4LHguCcM4AJwJnAxcDvgusIGuQohaWhoYGtW7finmgklmRKJmsqA4G17r7e3T8BZgNjY44ZC/zRI5YBXc3smBTPjTUWmO3ue939bWBtcJ2CNm/lJobNeFaDHKVgROcB69GjByNGjKChoaGji1RQMhlUSoGNoddVwbZUjmnu3ClBc9m9ZnZYC94PM5tsZhVmVpHv8wNF8yibEtRGlJyXfKR5wDpWJoNKvHkSYr8wJzom2bl3ACcBZcAW4FcteD/c/U53L3f38u7du8c5JX80N8jxZ+P6KJcieUfzgHWsTPb+qgKOC73uAWxO8Zj9Ep3r7lujG83sLuDxFrxfQdEgRylEuTgPWHgtGHeP+zwX7gMyW1NZDpxsZr3MbD8iSfT5McfMB64MeoENBna4+5Zk5wY5l6jPAa+FrjXBzPY3s15Ekv8vZermcoEGOUqhis4Dls0fxNHOBPX19Y05oPPOOy/u81zKDWWspuLudWY2BXgSKALudffXzezaYP8sYAFwCZGk+m7g6mTnBpf+hZmVEWna2gB8NTjndTObC6wG6oDr3D1+20+e0yBHkewWXlSsvLycioqKxhyQmcV9vm3btpyYMNMKuctdeXm5V1RUdHQx0ireDMQa5CiFLluWGo6Ww9057rjjqKuro6ioiIEDB7J8+XKGDBmCmbF06dImz4cOHcqiRYuypuZlZivcvTzePo2ozzPNDXIUKTTZstRwuBxDhgxh6NChjWV69tlnef/995PmVNyd9957r8MDY3MUVPKMBjmKNBWvi3F7NiOFayfRcrz44ou8++67dOrUqTFIRMsU73m2BMZUZGeppMU0yFEkvtguxt26dWu30fbhgZhf/OIXm5Tj6KOPTrkzQS6NvVFNJQ9oJUeRxMJdjLt168b555+f8W/8qdZOUhUNjNFyZ/PYG9VU8oAGOYokF+1i/P777zf5xr9169a4tZZU5g4LHxP7PB21k7BoYKyqqsqqhH08Cip5oLlBjgooIhHhprAhQ4YwYcKExnEgdXV1+4wbCW8PB49EY0tGjBjB1q1bm9ROZs+enZZgkAtjb0DNXzktOh5FeRSR1ISbwsLdepcsWcLw4cNZvnx5k3Ej4e3hLr7NjS0JN1UdffTRaQ8E2dJFOh4FlRyVz3mU2tpaqqqq2LNnT0cXRTpYly5d6NGjB8XFxWm7ZvQbv7s3fviXl5ezfPly6urqWL58eeO4kfD2cPAIHxM7nuSoo47K6DQx2d4TTIMfc2zwY3i0fCK5Psjx7bff5uCDD+aII47Ium9h0n7cne3bt7Nz50569eqVkfeIfuPv3r17kw/q6LiR8PbY4JFobEmm/2a3bt1Kjx49qKuro3PnzlRVVbX7SHsNfswTzdVOID8mi9yzZw89e/ZUQClwZsYRRxyR0e6z0VoLsE/tIt722OARb2xJpmV7TzAFlSwRrYFsrq7h2AQ1jWS9vKLyJY+igCLQvn8H4QCTaHt7Bo9Esn0WZgWVLBBbA9lUXcMND7/auD8abJprqMzlPIqIpC5RAMwG2ZPdKWDxaiA1tfXcPP/1xpUbmwsoGo/SMS655BKqq6uTHvPDH/6QZ555plXXX7RoEZ/97Gf32V5ZWcmCBQtadc3q6mp+97vfpXTszTffzC9/+cukx8ybN4/Vq1e3qizSdqmMqWlPCipZINE4k+qa2mabu0qKi7htfJnGo7Sz6JiFBQsW0LVr16THTp8+nQsuuCCt799eQSUVCiodJzzQMlvWXFFQyQKtyYMYqp2Epfvb2n/913/Ru3dvevfuzW233QbAhg0bOP300/na175G//792bhxIz179uT9998H4Mc//jGnnXYaF154IRMnTmz8hn/VVVfx4IMPAtCzZ09uuukm+vfvT58+fXjzzTcBeOmllxg6dCj9+vVj6NChrFmzJmHZPvnkE374wx8yZ84cysrKmDNnDh9//DFf+cpXGDBgAP369ePRRx8F4PXXX2fgwIGUlZXRt29f3nrrLaZNm8a6desoKytj6tSp+1z/lltu4dRTT+WCCy5oUo677rqLAQMGcNZZZ3HFFVewe/duli5dyvz585k6dSplZWWsW7cu7nGSGdk4J5iCShaYOupUSoqLUj6+tGsJb8+4VLWTQLq/ra1YsYLf//73/O///i/Lli3jrrvuYuXKlQCsWbOGK6+8kpUrV3LCCSc0nlNRUcFDDz3EypUrefjhh0nWVb1bt268/PLL/Md//Edj4DnttNN4/vnnWblyJdOnT+d73/tewvP3228/pk+fzvjx46msrGT8+PHccsstnH/++SxfvpyFCxcydepUPv74Y2bNmsV//ud/UllZSUVFBT169GDGjBmcdNJJVFZWMnPmzH3uffbs2Y33sXz58sZ948aNY/ny5bzyyiucfvrp3HPPPQwdOpQxY8Ywc+ZMKisrOemkk+IeJ5kRO1lmNvQEU6I+C0QDQ3PjT0DJ+HjSPbX5Cy+8wOc+9zkOPPBAIPJhunjxYsaMGcMJJ5zA4MGD454zduxYSkoitc7LLrss4fXHjRsHwNlnn83DDz8MwI4dO5g0aRJvvfUWZkZtbW2LyvzUU08xf/78xiC1Z88e3n33XYYMGcItt9xCVVUV48aN4+STT056ncWLF/O5z32OAw44AIAxY8Y07nvttdf4wQ9+QHV1Nbt27WLUqFFxr5HqcdJ22dgTTEGlnYW7Dh9aUowZVO+ubexGfP2cyoRJ+Vwf1Jgp6e63n6wJLRpoWnJOrP333x+AoqIi6urqALjxxhsZOXIkjzzyCBs2bGDEiBGpFzh4/4ceeohTT236heP0009n0KBBPPHEE4waNYq7776bE088Mem1En0wXXXVVcybN4+zzjqL++67j0WLFrXpOEmPbOsJpuavdhTtOhztzVVdU8uHu2txPu1G3PWA+NNRRFduVEDZV7pncB0+fDjz5s1j9+7dfPzxxzzyyCOce+65Sc8555xzeOyxx9izZw+7du3iiSeeaNF77tixg9LSyO/2vvvua/b4gw8+mJ07dza+HjVqFL/+9a8bg1u0uW79+vWceOKJfOMb32DMmDGsWrVqn3PDhg8fziOPPEJNTQ07d+7ksccea9y3c+dOjjnmGGpra7n//vsTliXRcdI+Oro3mIJKO2pu8GJNbT3u7JNfUZNX89I5g2v//v256qqrGDhwIIMGDeKaa66hX79+Sc8ZMGAAY8aM4ayzzmLcuHGUl5dz6KGHpvye3/nOd7jhhhsYNmwY9fXJe/wBjBw5ktWrVzcm6m+88UZqa2vp27cvvXv35sYbbwRgzpw59O7dm7KyMt58802uvPJKjjjiCIYNG0bv3r33SdT379+f8ePHU1ZWxhVXXNEkmP74xz9m0KBBXHjhhZx22mmN2ydMmMDMmTPp168f69atS3icZF429AbL6NxfZnYx8N9AEXC3u8+I2W/B/kuA3cBV7v5ysnPNbCZwGfAJsA642t2rzawn8AYQ7a6yzN2vTVa+9p77q9e0J5odb2LArePLmh1dn8/eeOMNTj/99I4uRovt2rWLgw46iN27dzN8+HDuvPNO+vfv39HFynm5+vfQEdprXrAOmfvLzIqA3wIXAlXAcjOb7+7hDu2jgZODxyDgDmBQM+c+Ddzg7nVm9nPgBuC7wfXWuXtZpu6prY7tWtJsIv7YriVc3q+0oIJIvpg8eTKrV69mz549TJo0SQFF2l1sfjG6dHJ7JvEzmagfCKx19/UAZjYbGAuEg8pY4I8eqS4tM7OuZnYM0DPRue7+VOj8ZcDnM3gPaTV11Kl5O129wF/+8peOLoIUuI5YOjlWJoNKKbAx9LqKSG2kuWNKUzwX4CvAnNDrXma2EvgI+IG7L25d0dMrtsdXl+JOVO+ujdv7SzUUEWmLaH4xvAJldOnkTp06ZbzWksmgEq/UsSmFRMc0e66ZfR+oA6LdS7YAx7v7djM7G5hnZme6+0cx500GJgMcf/zxzd5EW8VOFlldU0tJcRG3ji9TABGRjAk3hUWXTo7WWv7xj3+wffv2jASYTNaFqoDjQq97AJtTPCbpuWY2Cfgs8KWg6Qx33+vu24PnK4gk8U+JLZS73+nu5e5e3r1791beWuoSTRY588nE03CIiLRVuKv93LlzG2st0SWSM9VDLJNBZTlwspn1MrP9gAnA/Jhj5gNXWsRgYIe7b0l2btAr7LvAGHdvnFTIzLoHCX7M7EQiyf/1Gby/lCSaLDLRdhGRdIk2hR111FGN07kMGDCgyRLJ6Z4vLGNBxd3rgCnAk0S6+s5199fN7Fozi3b1XUDkg38tcBfwtWTnBuf8BjgYeNrMKs1sVrB9OLDKzF4BHgSudfcPMnV/qWhoaOCog/eLuy9fFtPKBvNWbmLYjGfpNe0Jhs14lnkrN7X5mkOHDm32mNtuu61dJku87777mDJlStJjFi1axNKlS1t87fCEmGE//elPW3yttli8eDFnnnkmZWVl1NToC1e6hWstS5Ysyeh8YRmdpsXdFxAJHOFts0LPHbgu1XOD7f+W4PiHgIfaUt50ig5CWv1BEUdc/A286NOR8urllT7JFjhrS84qlQ/o2267jS9/+cuN82Slor6+nqKi1CcPTdWiRYs46KCDUgqGqfjpT3+adFLLdLv//vv59re/zdVXX91u71loki2dnNb3SevVpFF0ksOdry1k+99v5+iD99N09RmQqZzVQQcdBEQ+rEeMGMHnP/95TjvtNL70pS/h7tx+++1s3ryZkSNHMnLkSCAyqeOQIUPo378/X/jCF9i1axcQqQ1Mnz6dc845h7/+9a+MGDGCb37zmwwdOpTevXvz0ksvAfDBBx9w+eWX07dvXwYPHsyqVav2Kddjjz3GoEGD6NevHxdccAFbt25lw4YNzJo1i1tvvZWysjIWL17Mtm3buOKKKxgwYAADBgxgyZIlAGzfvp2LLrqIfv368dWvfjXuVB7Tpk2jpqaGsrIyvvSlLwHxlwKI9fe//53+/ftz1lln8ZnPfAaILPI1adIkLrroInr27MnDDz/Md77zHfr06cPFF19MbW0td999N3PnzmX69OmN7yeZlc4ZKPbh7gX7OPvssz3d6uvr/V//+pfX19f78OHDvXPnzj58+HBvaGhI+3vlq9WrV6d8bM/vPu4nxHn0/O7jbSrDgQce6O7uCxcu9EMOOcQ3btzo9fX1PnjwYF+8eLG7u59wwgm+bds2d3fftm2bn3vuub5r1y53d58xY4b/6Ec/ajzu5z//eeO1zzvvPL/mmmvc3f25557zM888093dp0yZ4jfffLO7u//jH//ws846y93df//73/t1113n7u4ffPBB49/SXXfd5d/61rfc3f2mm27ymTNnNr7HxIkTG8v5zjvv+Gmnnebu7l//+tcby/X444870HgP8e7f3b2iosJ79+7tu3bt8p07d/oZZ5zhL7/8cpPj33vvPe/Ro4evX7/e3d23b9/eWK5hw4b5J5984pWVlV5SUuILFixwd/fLL7/cH3nkEXd3nzRpkv/1r3+N+7toyd+DtA+gwhN8rmqW4lZqaGjYp/oYbfJqj257EpFoloJ05qwGDhxIjx49ACgrK2PDhg2cc845TY5ZtmwZq1evZtiwYUBkIa0hQ4Y07h8/fnyT4ydOnAhEJnD86KOPqK6u5oUXXuChhyItuOeffz7bt29nx44dTc6rqqpi/PjxbNmyhU8++YRevXrFLfMzzzzTZDXGjz76iJ07d/L88883Trd/6aWXcthhhzV7/4mWAgjPh7Zs2TKGDx/eWJ7DDz+8cd/o0aMpLi6mT58+1NfXc/HFFwPQp08fNmzY0Oz7S25R81crJJq0LXZdj+3bt2euiilA/AXO0p2zik5VD02nqw9zdy688EIqKyuprKxk9erVTRanip0yP/ZvwsziNkXFHvf1r3+dKVOm8Oqrr/I///M/7NmzJ26ZGxoaePHFFxvLs2nTJg4++OC412xOvHLFOybRdaP/fp06daK4uLjxuE6dOsX9t5TcpqDSCrHBY+vWrWzdupXu3btn3Sps+e7yfqX8bFwfSruWtHvOKjzl++DBg1myZAlr164FYPfu3fzzn/9MeO6cOZGJIF544QUOPfRQDj30UIYPH944VfyiRYvo1q0bhxxySJPzwlPk/+EPf4hbFoCLLrqI3/zmN42vKysrAZq8x9/+9jc+/PDDuOUrLi5uXCgslaUAhgwZwnPPPcfbb78NRPJDUpjU/NUKHTVSVeLrqAk4J0+ezOjRoznmmGNYuHAh9913HxMnTmTv3r0A/OQnP+GUU/YZfwvAYYcdxtChQ/noo4+49957gUhS++qrr6Zv374ccMABTYJG1M0338wXvvAFSktLGTx4cOOH+GWXXcbnP/95Hn30UX79619z++23c91119G3b1/q6uoYPnw4s2bN4qabbmLixIn079+f8847L+GsEpMnT6Zv377079+f+++/v3EpACDuUgDdu3fnzjvvZNy4cTQ0NHDkkUfy9NNPt+4fVnJaRqe+z3Ztmfo+mlNxd4477riMTzVdSPJ9qvMRI0bwy1/+kvLyuDOHS4x8/3vIRcmmvlfzVyvFG6mqJi8RKXRq/mqj8FTTavKSVGjNdslnqqmkQUYHEhWoQm6WlU/p7yD3KKhI1unSpQvbt2/XB0qBc3e2b99Oly5dOroo0gJq/pKs06NHD6qqqtI+e6rkni5dujQOPJXcoKAiWae4uDjhSHERyW5q/hIRkbRRUBERkbRRUBERkbQp6BH1ZrYNeKejy5FEN2DfZfkKQ6Heu+67sOTqfZ/g7t3j7SjooJLtzKwi0VQI+a5Q7133XVjy8b7V/CUiImmjoCIiImmjoJLd7uzoAnSgQr133Xdhybv7Vk5FRETSRjUVERFJGwUVERFJGwWVDmZmRWa20sweD14fbmZPm9lbwc/DQsfeYGZrzWyNmY3quFK3nZltMLNXzazSzCqCbXl/72bW1cweNLM3zewNMxuS7/dtZqcGv+fo4yMz+2a+3zeAmV1vZq+b2Wtm9oCZdcn7+3Z3PTrwAXwL+AvwePD6F8C04Pk04OfB8zOAV4D9gV7AOqCoo8vfhvveAHSL2Zb39w78AbgmeL4f0LUQ7jt0/0XAv4AT8v2+gVLgbaAkeD0XuCrf71s1lQ5kZj2AS4G7Q5vHEvngIfh5eWj7bHff6+5vA2uBge1U1PaS1/duZocAw4F7ANz9E3evJs/vO8ZngHXu/g6Fcd+dgRIz6wwcAGwmz+9bQaVj3QZ8B2gIbTvK3bcABD+ji96XAhtDx1UF23KVA0+Z2Qozmxxsy/d7PxHYBvw+aPK828wOJP/vO2wC8EDwPK/v2903Ab8E3gW2ADvc/Sny/L4VVDqImX0WeM/dV6R6SpxtudwffJi79wdGA9eZ2fAkx+bLvXcG+gN3uHs/4GMizR+J5Mt9A2Bm+wFjgL82d2icbTl330GuZCyRpqxjgQPN7MvJTomzLefuW0Gl4wwDxpjZBmA2cL6Z/RnYambHAAQ/3wuOrwKOC53fg0hVOie5++bg53vAI0Sq+fl+71VAlbv/b/D6QSJBJt/vO2o08LK7bw1e5/t9XwC87e7b3L0WeBgYSp7ft4JKB3H3G9y9h7v3JNIk8Ky7fxmYD0wKDpsEPBo8nw9MMLP9zawXcDLwUjsXOy3M7EAzOzj6HLgIeI08v3d3/xew0cxODTZ9BlhNnt93yEQ+bfqC/L/vd4HBZnaAmRmR3/cb5Pt9d3RPAT0cYASf9v46AvgH8Fbw8/DQcd8n0iNkDTC6o8vdhvs9kUgvl1eA14HvF9C9lwEVwCpgHnBYgdz3AcB24NDQtkK47x8BbxL50vQnIj278vq+NU2LiIikjZq/REQkbRRUREQkbRRUREQkbRRUREQkbRRUREQkbRRURLKImS0ys/I0X7OrmX0t9HpEdFZskXRTUBHJf12BrzV3kEg6KKiIpMDMvmNm3wie32pmzwbPP2NmfzazO8ysIlg740fBvtFmNjd0jRFm9ljw/CIze9HMXjazv5rZQXHeM+4xFlmL5kfB9lfN7LRge/dgfY6Xzex/zOwdM+sGzABOCtYymRlc/iD7dF2X+4MR3yJtpqAikprngXOD5+VEPpSLgXOAxURmBSgH+gLnmVlf4Gki03QcGJw3HpgTfND/ALjAI5NqVhBZV6dRCse8H2y/A/h2sO0mItP99Ccyn9rxwfZpRKabL3P3qcG2fsA3iazhcSKRuehE2kxBRSQ1K4CzgznL9gIvEgku5xIJKl80s5eBlcCZwBnuXgf8HbgsWE/jUiLzPA0m8mG+xMwqicz/dELM+zV3zMOhcvUMnp9DZHJS3P3vwIdJ7ucld69y9wagMnQNkTbp3NEFEMkF7l4bzCh9NbCUyNxdI4GTgBoitYUB7v6hmd0HdAlOnQNcB3wALHf3nUFT09PuPjHJWzZ3zN7gZz2f/j9uSRPW3tDz8DVE2kQ1FZHUPU8keDxPpHZyLZFv+YcQWRtlh5kdRWSK96hFRKa3/z9EAgzAMmCYmf0bQDCL7Skx75XKMbFeAL4YHH8RkckqAXYCB7fkRkVaS0FFJHWLgWOAFz2yJsgeYLG7v0Kk2et14F5gSfQEd68HHicSaB4Ptm0jslb5A2a2ikgAOS38RqkcE8ePgIuCZrjRRFYb3Onu24k0o70WStSLZIRmKRbJE2a2P1Dv7nVmNoTICpNlHVwsKTBqRxXJH8cDc82sE/AJkSY3kXalmoqIiKSNcioiIpI2CioiIpI2CioiIpI2CioiIpI2CioiIpI2/w+IHEyl7YFi4gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "test.plot(x='wavelength', y='y1', kind='scatter', marker ='.', c='k')\n",
    "plt.plot(iwavelength, irho, 'o')\n",
    "plt.title('Test spectrum')\n",
    "plt.legend(['original test data','interpolated to cmf'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wavelength overlap 380 720\n"
     ]
    }
   ],
   "source": [
    "# find overping wavelengths between cmf and multispectral image\n",
    "start = max([cmf['wavelength'].iloc[0], iwavelength.min()])\n",
    "end = min([cmf['wavelength'].iloc[-1], iwavelength.max()])\n",
    "print('wavelength overlap', start,end)\n",
    "cmf = cmf[ (cmf['wavelength'] >= start) & (cmf['wavelength'] <= end) ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating the tristimulus values for X, Y and Z by Riemann sum approximation of the integrals with $\\Delta\\lambda$ = 1 nm\n",
    "\n",
    "$$\n",
    "X = \\sum_{i=413}^{708} [\\rho_w]_N(\\lambda_i)\\tilde{x}(\\lambda)\\Delta\\lambda \\\\\n",
    "Y = \\sum_{i=413}^{708} [\\rho_w]_N(\\lambda_i)\\tilde{y}(\\lambda)\\Delta\\lambda \\\\\n",
    "Z = \\sum_{i=413}^{708} [\\rho_w]_N(\\lambda_i)\\tilde{z}(\\lambda)\\Delta\\lambda\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Creating a dictionary to hold the reflectance values aas tristimulus\n",
    "r = {\"X\": [], \"Y\": [], \"Z\": []}\n",
    "# ------- R_RS * cmf\n",
    "for i in range(0, len(cmf)):\n",
    "    r[\"X\"] = irho * cmf[\"x\"]\n",
    "    r[\"Y\"] = irho * cmf[\"y\"]\n",
    "    r[\"Z\"] = irho * cmf[\"z\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# ------ Sum\n",
    "s = {\"X\": [], \"Y\": [], \"Z\": []}\n",
    "\n",
    "s[\"X\"] = sum(r[\"X\"] * Delta_lambda)\n",
    "s[\"Y\"] = sum(r[\"Y\"] * Delta_lambda)\n",
    "s[\"Z\"] = sum(r[\"Z\"] * Delta_lambda)\n",
    "\n",
    "sum_xyz = s[\"X\"] + s[\"Y\"] + s[\"Z\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# ------ chromaticity\n",
    "chrom = {\"x\": [], \"y\": [], \"z\": []}\n",
    "chrom[\"x\"] = s[\"X\"] / sum_xyz\n",
    "chrom[\"y\"] = s[\"Y\"] / sum_xyz\n",
    "chrom[\"z\"] = s[\"Z\"] / sum_xyz\n",
    "\n",
    "sum_chrom_y = chrom[\"x\"] + chrom[\"y\"] + chrom[\"z\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ______ chromaticity - whiteness\n",
    "chrom_w = {\"x\": [], \"y\": []}\n",
    "chrom_w[\"x\"] = chrom[\"x\"] - (1 / 3)\n",
    "chrom_w[\"y\"] = chrom[\"y\"] - (1 / 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ______ calculate atan2\n",
    "# we use the average atan per scale\n",
    "# a_i, refers to the alpha M angle\n",
    "a_i = math.atan2(chrom_w[\"y\"], chrom_w[\"x\"]) * 180 / math.pi\n",
    "\n",
    "if a_i < 0:\n",
    "    a_i = a_i + 360\n",
    "else:\n",
    "    a_i = a_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- fui approximation\n",
    "\n",
    "if a_i >= fui[\"continuous\"][\"alpha\"][0]:  # FUI = 1 its > Average\n",
    "    fu_i = 1.0\n",
    "elif math.isnan(a_i):  # FUI = NAN = 0\n",
    "    fu_i = 0\n",
    "elif a_i <= fui[\"continuous\"][\"alpha\"][200]:  # FUI = 21\n",
    "    fu_i = 21.0\n",
    "else:\n",
    "    for c in range(0, 200):\n",
    "        if (fui[\"continuous\"][\"alpha\"][c] > a_i) and (a_i >= fui[\"continuous\"][\"alpha\"][c + 1]):\n",
    "            fu_i = fui[\"continuous\"][\"value\"][c + 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18.1"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The calculated FUI is\n",
    "fu_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}