MEANGB*key_units.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f82f0193-8002-4584-a7cf-67e4cfca2e15",
   "metadata": {},
   "source": [
    "## Units for the MEANGB\\* keywords"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2ee54d2-f91f-4459-9cf1-2de63f9f3e31",
   "metadata": {},
   "source": [
    "The MEANGB\\* keywords are in units of Gauss/pixel, not Gauss/Mm. To convert from Gauss/pixel to Gauss/Mm, multiply the values of MEANGB\\* by a constant factor of 2.74.\n",
    "\n",
    "This notebook describes the problem and derives the solution for converting these keys from Gauss/pixel to Gauss/Mm. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3daad82f-66a0-4740-a158-952a7e4cb99d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sunpy.map\n",
    "import scipy.ndimage\n",
    "import numpy as np\n",
    "import math\n",
    "import drms\n",
    "from datetime import datetime as dt_obj\n",
    "import matplotlib.pylab as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bf5d2325-4794-44dd-98db-a02196b75544",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define some constants\n",
    "radsindeg = np.pi/180.\n",
    "munaught  = 0.0000012566370614"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edf268df-a463-4e59-b1c4-f8a9715fec03",
   "metadata": {},
   "source": [
    "To query the SHARP data, register your e-mail address with [JSOC](http://jsoc.stanford.edu/ajax/register_email.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f323d7da-d8e8-4731-945e-ae1077d8120d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "http://jsoc.stanford.edu/cgi-bin/ajax/checkAddress.sh?address=mbobra%2540stanford.edu&checkonly=1\n"
     ]
    }
   ],
   "source": [
    "email = 'your@email.address'\n",
    "c = drms.Client(debug=True, verbose=True, email=email)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aeeaee7e-bb52-4e1f-ba4f-f51ac99a14bc",
   "metadata": {},
   "source": [
    "The MEANGB\\* keywords store the value of the derivative of the horizontal field (MEANGBH), vertical field (MEANGBZ), line-of-sight field (MEANGBL), and total field (MEANGBT). \n",
    "\n",
    "The formula to calculate the deriviative is $\\sqrt{\\left(\\frac{dB}{dx}\\right)^2 + \\left(\\frac{dB}{dy}\\right)^2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9ee15be-fd15-4656-aec0-f50b8b75c34f",
   "metadata": {},
   "source": [
    "The units of the magnetic field, or dB, are in Gauss. The units of length, i.e. dx or dy, are in pixels. Therefore the units of dB/dx or dB/dy in Gauss/pixel. The keys MEANGB\\* in the series hmi.sharp_720s and hmi.sharp_cea_720s are in Gauss/pixel.\n",
    "\n",
    "To convert from Gauss/pixel to Gauss/Mm, multiply the value of MEANGB\\* by a factor of (Gauss/pix)(pix/arcsec)(arsec/meter)(meter/Mm), which results in Gauss/Mm. Practically, this means multiplying the value of MEANGB\\* by the following factor: (1/cdelt1_arcsec)(RSUN_OBS/RSUN_REF)(1000000), where 1/cdelt1_arcsec is the angle subtended on the Sun in arcseconds as observed by the spacecraft (see [here](https://github.com/mbobra/SHARPs/blob/master/cdelt1_arcsec.pdf) for a derivation of 1/cdelt1_arcsec). This value of this factor is 2.74."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7df0a268-8e68-4256-95fd-ff33dc7ed703",
   "metadata": {},
   "source": [
    "First, we will look at a specific example for HARPNUM 86 at 2010.07.14_12:24:00 TAI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "98ea37e9-6550-4b2c-bdf9-9d09d4bb9019",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify a record\n",
    "t_rec = '2010.07.14_12:24:00_TAI'\n",
    "harpnum = '86'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8a1bfd2f-99de-4e8c-9b4b-7eee0c27ab05",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "http://jsoc.stanford.edu/cgi-bin/ajax/jsoc_info?op=rs_list&ds=hmi.sharp_cea_720s%5B86%5D%5B2010.07.14_12%3A24%3A00_TAI%5D%5B%3F+%28QUALITY%3C65536%29+%3F%5D&key=%2A%2AALL%2A%2A\n",
      "http://jsoc.stanford.edu/cgi-bin/ajax/jsoc_info?op=series_struct&ds=hmi.sharp_cea_720s\n"
     ]
    }
   ],
   "source": [
    "# Query the data\n",
    "hmi_query_string    = 'hmi.sharp_cea_720s['+harpnum+']['+t_rec+'][? (QUALITY<65536) ?]'\n",
    "keyhmi             = c.query(hmi_query_string, key=drms.const.all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "63d0f0d3-5b9a-43c5-b05b-123242a36dbc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.49428262675999274"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate cdelt_1 arcsec using the keyword metadata from this record\n",
    "cdelt1_arcsec = (math.atan((keyhmi['RSUN_REF'][0]*keyhmi['CDELT1'][0]*radsindeg)/(keyhmi['DSUN_OBS'][0])))*(1/radsindeg)*(3600.)\n",
    "cdelt1_arcsec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b0afc91c-482b-46d6-b50a-c896f5232512",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the factor using the keyword metadata from this record\n",
    "factor = (1/cdelt1_arcsec)*(keyhmi['RSUN_OBS'][0]/keyhmi['RSUN_REF'][0])*(1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "94eaed9e-9f0e-4e9f-b83e-a5221b1e04cb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For this record, we need to multiply the values of MEANGB* by a factor of 2.7440606028364725 .\n"
     ]
    }
   ],
   "source": [
    "print(\"For this record, we need to multiply the values of MEANGB* by a factor of\",factor,\".\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d647453b-e749-4459-a31d-f02a5888ab69",
   "metadata": {},
   "source": [
    "While this factor includes variable keywords, the variability cancels out to create a constant factor. The rest of this notebook will show why this factor is a constant.\n",
    "\n",
    "Calculating this factor depends on the value of four keywords in the hmi.sharp_cea_720s series:\n",
    "* `RSUN_REF`, a constant keyword that describes the reference radius of the Sun. The value of `RSUN_REF` is always 696,000,000.0 m.\n",
    "* `CDELT1`, a constant keyword that describes the map scale in the X direction (`CDELT2` describes the map scale in the Y direction, and both have the same value). The value of `CDELT1` is always 0.03 degrees for data in the hmi.sharp_cea_720s series.\n",
    "* `RSUN_OBS`, a variable keyword that describes the angular radius of the Sun in arcseconds.\n",
    "* `DSUN_OBS`, a variable keyword that describes the distance from SDO to the center of the Sun in meters.\n",
    "\n",
    "Now we will look at the maximum and minimum values of `RSUN_OBS` and `DSUN_OBS` and how they change over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0e46d8bf-3c56-4d4f-b339-55561d65d64a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "http://jsoc.stanford.edu/cgi-bin/ajax/jsoc_info?op=rs_list&ds=hmi.sharp_cea_720s%5B%5D%5B2019.01.01%2F800d%407h%5D&key=T_REC%2CCDELT1%2CRSUN_OBS%2CDSUN_OBS%2CMEANGBL\n"
     ]
    }
   ],
   "source": [
    "hmi_query_string    = 'hmi.sharp_cea_720s[][2019.01.01/800d@7h]'\n",
    "keyshmi             = c.query(hmi_query_string, key='T_REC,CDELT1,RSUN_OBS,DSUN_OBS,MEANGBL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5a0cba3e-fd89-486c-babc-8def7de7ddeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def parse_tai_string(tstr,datetime=True):\n",
    "    year   = int(tstr[:4])\n",
    "    month  = int(tstr[5:7])\n",
    "    day    = int(tstr[8:10])\n",
    "    hour   = int(tstr[11:13])\n",
    "    minute = int(tstr[14:16])\n",
    "    if datetime: return dt_obj(year,month,day,hour,minute)\n",
    "    else: return year,month,day,hour,minute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c04e8e44-653d-49d4-a708-96059e440d1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_rec_hmi = np.array([parse_tai_string(keyshmi.T_REC[i],datetime=True) for i in range(len(keyshmi.T_REC))])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d148951-557d-4827-8d36-ea9ef3ebd14f",
   "metadata": {},
   "source": [
    "Plot variation in RSUN_OBS over time (this is not smooth because it is overplotting many HARPs; otherwise it would be totally smooth)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "0b2540a4-22f1-40c6-b04d-c866ba9168ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Angular Radius of the Sun over time as observed by SDO')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnkAAAFYCAYAAAAvEeqvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABiV0lEQVR4nO3dd5hU5fXA8e/ZZZeyLH3pS+8dXMUGoiAgEHuLLRqj0Rhji0aj8UcsscXEGGOJJbZoxF5AQEVQQECQ3ttSl7J0WLaf3x/37uxs34WZuVPO53n22ffeuTNz7s7cnTNvFVXFGGOMMcZElzivAzDGGGOMMYFnSZ4xxhhjTBSyJM8YY4wxJgpZkmeMMcYYE4UsyTPGGGOMiUKW5BljjDHGRCFL8owph4hMF5FfeR1HeURkmIhs9dteLiLDQhzDIyKSKSI7qnn8eBF5O9hxRTMReVFE/uR1HIFQ+j0cCUSkg4ioiNSq4PZ0ERkR6riMqYwleSbsuQnXPhGp7XUsx8r9cDgiIodFZJuI/E1E4gPx2KraW1WnB+KxqkNEUoG7gF6q2rKc24P6AS4ibUXkQzfJPCAiS0Xk2mA9nxdE5FoRmem/T1VvUtWHvYrJhI6I/FFENrr/L7aKyHt+t00XkWwROSQiB0VkgYjcW/r/o4j0EpHP3GvkkIh8KyKnhv5sjJcsyTNhTUQ6AEMABc71NpqqVfQt39VfVesDZwCXAb8MTVQB1x7Yo6q7PHr+t4AtbhxNgWuAnR7FctyqeM+YYxSpf1cR+QVwNTDC/X+RBnxT6rDfqmoy0ArnC9flwCQREfcxOgOzgKVAR6A18DEwVUROCcmJmLBgSZ4Jd9cAc4DXgV/43yAir4vIv0RkovtNda77z63o9pEistr9Jvu8iMwoaoIt3XxYWVOMiHQWkWkissetPfqviDTyuz1dRP4gIkuAI1V9uKjqOpx/wAP8HuMfIrLF75v5EL/b6rrnuk9EVgAnlorP10zkHveI322lm3b/4NYkHnL/NsPLi1FEGorImyKyW0Q2icgDIhLnPs9XQGu3luH1UvdLAr70u/2wiLR2b050H/OQOE3MaX73a+3Wzu12azB+V8mf8ETgdVU9oqr5qrpQVb8s73zL+fuMF5EJFcVRzt/hVBH50X0P/VhUEyIil4vI/FLH3iEin7nl2iLyVxHZLCI7xWlqresfo/ta7AD+U+pxegIvAqe4f7/97n7fa+v3GPeIyC4RyRCR80VkjIisEZG9IvJHv8eME6e2Z737Pp4gIk0qOOfGIvKF+1rsc8tt/W6/VkQ2uH+/jSJyZQWPU1tEnhGR7e7PM1K2tumP7jWV7v847nmscJ9jm4j83u+2cSKySET2i8hsEennd1vpa/EBEfmg1HP+Q0SedcsNReRV9++3TZxuCPHubfHua5gpIhuAseWdZyknunHvE5H/iEgd97GWicjP/GJIcB93QHmPAUxR1fUAqrpDVf9d3pO518B0nC/Ap/jFOB74QVXvV9W9qnpIVZ/F+YL0RDXOw0QJS/JMuLsG+K/7M0pEWpS6/efAn4HGwDrgUQARaQZ8ANyHU9uzGjjWpgoBHsP5NtwTSMX5J1o6jrFAI1XNr/TBRHrg1E6u89v9I07S1wR4B3i/6AMC+D+gs/szilLJbrVPQqQ78FvgRLcWYBSQXsHh/wQaAp1wah6vAa5T1a+Bc4DtqlpfVa/1v5OqHil1e31V3e7efC7wP6AR8BnwnBtXHPA5sBhoAwwHbheRURXENgf4l5totavRH6GSOEpzk6CJwLM476G/ARNFpKl7v+4i0tXvLlfgvHbgfJB2w3lNu7jn9aDfsS1xXuv2wI3+z6uqK4GbcD6k66tqowrOoyVQx++xXwauAk7AeX89KCKd3GN/B5yP81q2BvYB/6rgceNwEs/2QDvgKMWvVZL79zjHfQ+dCiyq4HHuB052/wb9gZOAB0rF38yN/xfAv933KMCrwK/d5+gDTHOffxDwGvBrnNfkJeCzUsmj71rESWrGiEgD9/7xwKUUv05vAPk4r9FAYCRQ1Bf3BmCcuz8NuLiC8/R3Jc511Rnn9S863zdxXpsiY4AMVV1UzmPMAa4RkbtFJE2q0a1DVTcD83Fed4CzgffLOXQCcJqI1KvGuZhooKr2Yz9h+QOcDuQBzdztVcAdfre/Drzitz0GWOWWr8H5kCy6TXCa+H7lbo8H3va7vQNOk3Atd3t60bHlxHU+sNBvOx34ZRXnosBB4IhbfheoXcnx+3CadwE2AKP9brsR2Frq+Uf4/U0e8bttWNGxOB9ku4ARQEIlzx0P5OD0uSva92tgeunHrOD+ZW53/95f+233Ao665cHA5lLH3wf8p4LHbww8DiwHCnCSjBMreW7/v0+FcZTzPFcD80rt+wG41i2/DTzolrsCh4B67nvtCNDZ736nABv9YswF6lTyN7wWmFlqn++1dR/jKBDvbie776vBfscvAM53yyuB4X63tcK5tmpV4zocAOxzy0nAfuAioG4V91sPjPHbHgWk+8WfDyT53T4B+JNb3uy+5xqUeswXgIdL7VsNnFHRtQjMBK5xy2cD691yC5z3eV2/Y38OfOuWpwE3+d02Er//EeWcb3qp48f4PVdr9/3RwN3+ALinkr/dlcDX7vtoD3Cv323TKed/E84Xl5fdcj5+/zP8junhnkObql53+4mOH6vJM+HsF8BUVc10t9+hbC2W/+jOLKC+W26Nk9QBoM5/uGMaDCAizUXkf25zzkGcD/dmpQ7bUs5dSxvkxncZTmKT5Pccd4nISrdZcD9OLVrRc5Q4F2DTsZyHOs3Et+MkOrvcc2pdzqHNgMRSz7MJp8bleJR+reqI07TdHqd5d3/RD/BHnA/h8s5jn6req6q93WMWAZ+IOP2RjiOO0lpT9m/t/3d4BycpAKcW7xNVzQJScJK9BX7nM9ndX2S3qmZXM96K7FHVArd81P3t3zfxKMXXQ3vgY794VuIkyGX+xiJST0ReEqeZ/iDwHdBIROLVqam9DKemMUOcrhI9Koiv9N9vk7uvyD738cq7/SKcJGmTON0sivqRtQfuKvVeSS31uKWvxdKvU1EtXnsgwT2Posd6CWjuF39Nr7vSx7cGUKc2exZwkThdPc7BaZ0ol6r+V1VH4NRG3gQ8VEnNdpE2wF63nImTyJfWCijE+RJpYoAleSYsidN/6VLgDBHZIU7fpTuA/iLSvxoPkQH49yMS/22cb8j+TRZlRon6eQzn228/VW2A0+xSOqHQasSEOibg1Ag96MY2BPgDzvk2Vqd57oDfc2TgfJAVqayJstLzUtV3VPV0nA84pfz+OZk4tTztSz3ntsrOzf9pqnlckS04tVyN/H6SVXVMlU/kfAH4K86HaRNKnb/b1JVS/r2rtJ2SfwMo+XeYCjRz+1X9nOLkIRMnwertdz4N1elE7wu9iueu6d+wKltwmlj9/8Z1VLW81/QuoDtOrWADYKi7XwBUdYqqno2TMKzCaSYuT+m/Xzt3X5HGbvNvmdtV9UdVPQ8n4foEp5av6DweLXUe9VT1Xb/HKf23ex8Y5vYrvIDi12kLTk1eM7/HauB+eYCaXXdFSh/vf75v4PzvuASnlaHK60lV81T1fWAJTrN1ucQZ8X4C8L2762v3eUq71H3urKqe20QHS/JMuDofp6ahF05z0QCc/nDf4zTFVmUi0Feczui1gFsomfAsAoaKSDsRaYjTPFiRZOAwsF9E2gB31+REKvA4cKOItHQfPx/YDdQSkQeBBn7HTgDuE6dDfFvg1koedxFOH6Qm7mPfXnSDiHQXkbPc/kvZOIlIQekHcGuHJgCPikiyiLQH7sSpwayOnUBT9+9aHfOAg26H+bpuh/c+InJieQeLyBPu7bVEJBm4GVinqnuANTg1c2NFJAGnT9SxTr0zCegmIle4z3UZzvvxCwB1+l5+ADyFk2B+5e4vxEl8/i4izd2Y21SjJsbfTqCtiCQeY+ylvYjzerZ340kRkfMqODYZ572x3+2X+H9FN4hICxE5103OcnCuizLvIde7wAPuczXD+VJT+j30ZxFJdL/ojMPpi5ooIleKSENVzcPp5lD0HC8DN4nIYHEkua91ckUnrqq7cZo4/4PzZWKluz8DJ1F/WkQaiDM4pbOInOHedQLwO3Gm7GkM3FvRc/i5xT2+CU5t9Ht+t32CU5t/G04fvXKJM7BlrHvtxYnIOUBvYG45x9Zz4/0U5zqa5N70Z+BUEXnU/V+QLCK34vzv/EM1zsNECUvyTLj6BU6frM3qjC7boao7cDqAX1lB85qPW8NzCfAkTp+WXjgdk3Pc27/C+Qe8BKfv0heVPNyfcf45H8BJHj86nhNzn38pMAMnYZyCMyJ1DU4TTzYlm33+7O7fiPOh9FYlD/0WzgCGdPdY/w+Z2jjJZSZOk2VznA+i8tyKUyu2AadP0zs4Hd6rc26rcD7gN7jNYOU1CfsfXwD8DCeR3+jG9wpOk3V56uFMB7Hfja897vQ6qnoA+I17/23uORxTM72bNI7DqdnaA9wDjPPrPgDO32UE8L6WHHDzB5yBNXPcJs+vcWrHqmsaTp/DHSKSWdXB1fAPnMEiU0XkEE7n/sEVHPsMUBfndZiD09RcJA7n77Edp2nwDJy/d3kewbnmluBM5fGTu6/IDpxmw+04TZc3ue8dcPpDprt/u5twBy2o6nycARHPufddh9N/sSpFr9M7pfZfg9M1YYX7eB9Q3Mz5Ms61udiNvTrX/Ts4190G98d3vqp6FPgQZ0qTyh7rIM51uRnnPf4kcLOq+s+b+Jz7Ou7Eeb0+xOmDV+g+11qcPs39cf4XZOA0gY9S1VnVOA8TJcTpqmRMdBNnBOdW4EpV/dbreIwxscetpe+mqldVebAxAWA1eSZqicgoEWnkNk/+EadP0RyPwzLGxCC3Cfd6oNw574wJBkvyTDQ7BWcah0yc5sDz3SYTY4wJGRG5AacLxpeq+p3X8ZjYYc21xhhjjDFRyGryjDHGGGOikCV5xhhjjDFRqNJpKKJNs2bNtEOHDl6HYYwxxhhTpQULFmSq6rFO6B5bSV6HDh2YP3++12EYY4wxxlRJRI5pGcsi1lxrjDHGGBOFLMkzxhhjjIlCluQZY4wxxkQhS/KMMcYYY6KQJXnGGGOMMVHIkjxjjDHGmChkSZ4xxhhjTBSyJM8YY4wxJgpZkmeMMcYYE4ViasULY4zJOHCU7LxCjuTk07JhHZrVr+11SMYY4NGJK+jaIplGdRMY2bul1+FEBUvyAuzbVbuYs2EP783fQp/WDemcksT4c3sjIl6HZkzM+nzxdpomJXLFK3MrPGbBAyNokpRo16oxIaKqdLxvUoW3f3/PmaQ2qRfCiKKPqKrXMYRMWlqaBnPt2o2ZRzjzr9PL7J93/3CaJ9cJ2vMaY8r35dIMbv7vT9U+PiW5Nl/fcQYN6yUEMSpjTId7J1bruI9/cyoXPD+bB8b25JK0VBrWja1rU0QWqGrasd7f+uQFUGKt8v+cB7LyQhyJMbFt58FsOtw7sUYJHsDuQzk89MWKIEVljDlwNK/aCR7Abf9bBMAjE1fS/89Tmb0uk1iqnDpeluQFUOMKvv3f/cESNuw+HOJojIlNj01ayeC/fHPM9z9w1L6UGRMMM9dm0v/PU2t0n817s0psX/HKXDreN4ndh3ICGVrUsiQvgOomxJe7f2PmEcY+O5O35myybyDGBNEP6/fw0ncbqjyuQZ2KuyPXSyz/OjbGHLuD2Xlc9WrFfWJr6sRHv7bP02qwPnkBdiQnnzgRFOVwdj6NkxLZeySXuz9YwndrdnNGtxSevLgfLRpYHz1jAiU3v5Cb3l7AtFW7yr29cb0E/n1NGie0a0xcXPHAisJC5YMFW7nnwyUljt/42BgbgGFMgBw4mldpDV7/tg356DenEe9em/kFhazacYjXZm1k35Fcvl29u8L7pj8+NuDxhpPj7ZMXVkmeiNwG3AAI8LKqPiMi7wHd3UMaAftVdYCIdABWAqvd2+ao6k2VPX4okryKqCpvz93MoxNXUCchnkfP78vYfq0oKFS6PfAlrRrWYfrvh1Er3ipXjamp12dtZPznTl+6T285jfP+NQuAWnHC0vGjqFuN2rldh7I56VGnmfeKwe3487m9SbDr0Zjjkl9QSJf7vyz3tq/uGErXFsmV3l9VeW1WOg9X0lf2ljM7c/eoHscVZ7iKmiRPRPoA/wNOAnKBycDNqrrW75ingQOq+pCb5H2hqn2q+xxeJnlFNuw+zB0TFrN4y37OH9CaO8/uztCnvgWgc0oS39w1zNP4jIk0h3PyOePJb9lzJJcXrxrE6D6tyM0v5HBOPk2SEmv0WIWFypNTVvPijPUM6dqM564YRMO6CeQVFJKTX0j92jbrlDE1UdEgi5rWwFU13cqqh0dTp4IuU5EsmkbX9sSpjctS1XxgBnBB0Y3itJ1cCrzrUXwB0SmlPh/edAp3jOjG50syOPdfM323rd99pEajjowx8Mr3G9hzJJePf3Mqo/u0ApyR7jVN8ADi4oR7z+nBkxf344f1e7johdls3pPF7f9bRJ//mxLo0I2JaoFK8ABEhPTHx3Jq56bl3t7jT5NZteNgjR832oVTkrcMGCoiTUWkHjAGSPW7fQiw079mD+goIgtFZIaIDAllsMejVnwct43oyse/ObXcD6Lv11bc/8AYUyzzcA4vf7eBc/q0ZGC7xgF73EvTUnnr+sFkHs7hvH/NZOLSDACu/c88cvMLA/Y8xkSrzxdvL3f/8fahe+eGk1n18Ohybxv9zPc2GKOUsEnyVHUl8ATwFU5T7WIg3++Qn1OyFi8DaKeqA4E7gXdEpEHpxxWRG0VkvojM3707vJKnfm0bMfHWsrnp1a/OY92uQx5EZExkeW7aOrLzC/n9qO5VH1xDp3Ruyse/OY1G9Yq/iE1fvZtfvDYv4M9lTDS598Ml3PruwhL70to3DtggiToJ8cx/YES5t1XWpBuLwibJA1DVV1V1kKoOBfYCawFEpBZwIfCe37E5qrrHLS8A1gPdynnMf6tqmqqmpaSkhOI0aqRuYjw/3l/2zTrib995EI0xkWPL3iz+O3cTl6a1pXNK/aA8R8dmSXz8m1NL7Pthwx7e/CE9KM9nTDT4349byuz74OZTyzny2DWrX5vv7j6z3NvyC6y2vUhYJXki0tz93Q4nqSuquRsBrFLVrX7HpohIvFvuBHQFqp4gKwylJNfmwXG9yuyfuTbTg2iMiQxPT11NnAi3DS/z3S6gGtVL5O+X9S+x78FPl7M/Kzeoz2tMJCqvH97Gx8YE5bnaNa3Ht78fVmZ/l/u/5GhuQVCeM9KEVZIHfCgiK4DPgVtUdZ+7/3LKDrgYCiwRkcXAB8BNqro3dKEG1i9P78i1p3Yose+qV+eSk29vVGNKW7H9IJ8u3s51p3WkZcPgzzl5jjugw9+Ah76y/j/G+Hlpxvoy+4I952THZkl8f0/ZGr3Tn5gWtOeMJGGV5KnqEFXtpar9VfUbv/3XquqLpY79UFV7u8cOUtXPQx9xYI0/t3eZfWt32nJoxpT25JRVNKiTwM1ndA7J89VJKL9bxVtzNoXk+Y0Jd+M/W85jX64qse+TW04LyaTiqU3qldm354jVtEOYJXmGMp1JkytZfsmYWDRnwx6mr97Nb4Z1pmEF60UHQ0pybbq1KNn378FPl1NQaLV5xrw+O73EdqeUJAakNgrZ85c3qCM980jInj9cWZIXZprVr81pXYrnAWpQJ3QfYsaEO1Xl8S9X0aphHX5RqntDKEy5fSg9Wpacof+Kl+eEPA5jwknpfngdmyUxzYOJ/Rf+6ewS28P+Op19MV6jZ0leGHr7+sG+8uz1eziQlUfm4RwPIzImPExZvpNFW/Zz+4iunsxuLyK8f9MpJfbN3bjXBmGYmFVev9TJt3szbW3jpETWPnoOA9s18u0b+PBXnsQSLizJC0Miwvq/jKFL8/o88/UaBj48lbRHvvY6LGM8lV9QyJNTVtE5JYmLBrX1LI7kOgnMvvesEvs+X5LhUTTGeKv0vHRv/PIkatfybnmxhPg4RvZqWWJfLK+EYUlemIqPE+4Y0Y21uw5T1OVnxfbYfaMa88GCrWzYfYR7RvegVry3/7paN6pbYvujn7ZWcKQx0Wv1jrKT9p/Rzfv5aK8/vSMntC9eAeeTheWvvhELLMkLY+f0aVmi/8+YZ7+3SR5NTMrOK+CZr9cyqF0jRvZq4XU4ACz/8yjO6dOSa0/twMLN+9mw20bCm9gy6pmSk/Zv+Etw5sOrqcRacbx348m+7U8WbqMwRgdIWZIXxuLihLtGllyu6aXvInK+Z2OOy+uz09lxMJs/jO4RkikZqiOpdi1euOoEbh7WmTiBRyeutAnMTcyYva7ke/0Po3sQFxce1yY4a8T/6vSOAOw4mM289IidRve4WJIX5kb0bF5i+6kpqz2KxBhvHMjK4/lv13Fm9xQGd2pa9R1CrEWDOpzeNYVvVu3iqlfnkpWbX/WdjIlwV7wyt8T2zcNCM2dlTTwwrhdLx4+kXmI8ny7a5nU4nrAkL8yJCEvGjyyx75/frPUoGmNC7/kZ6ziUk889o3t4HUqFLhzYxleeaIMwTJRbX6prQun5XcNJcp0ERvduycQlGTG5gpQleRGgQZ0EkhKLRys9/dUaD6MxJnQyDhzl9VnpXDCgDT1bNfA6nAqN7F3cT3DC/LKLsxsTTYY/PaPEdrP6tT2KpHrOG9iGg9n5PPDxMjbtia0Jki3JixDL/jyKe88prsn4fu1uD6MxJjT+8fVaVOGOs7t5HUql6iXWYsKvT+HStLb8mL6vTE2HMdHi9+8vLrG99tFzPIqk+k7r3JRm9Wvz/oKtnPHUdK/DCSlL8iKEiHCT3zqdV786zyZgNVFt3a7DTJi/hStPblfu2pTh5qSOTfj9qO7ExwkTfrTaPBOdPlhQcrqgBI+nM6qOWvFxjO5TXNv+7epdHkYTWuH/6pgKvfmDLY5uotdfp6ymXmItfntmF69DqbbmyXU4q0dzPvxpK3k23ZGJMqVXXgqXKVOq42S/QVvX/edHDyMJLUvyIszqR0b7ytZka6LVT5v3MXn5Dm4Y0ommYd7fp7TL0lLJPJzLtFWxU1tgYkPplZfCacqUqgzp4v0kzV6wJC/C+C8X82P6Pg8jMSY4VJUnvlxFs/qJ/GpIR6/DqbFh3VNonlzbmmxNVDmYnVdiO/3xsR5Fcmwa1ktgqd9MFTe/vcDDaELHkrwIdOXgdr7yTW8tIDsv9oaFm+g1fc1u5m7cy++GdyWpdi2vw6mxWvFxXHRCW75dvYudB7O9DseYgPhoQeQv3ZeUWPz/5MtlO1i27YCH0YSGJXkR6NEL+vL8lYMAmLx8B9+tsWZbEx0KC5UnJ6+mXZN6XH5iu6rvEKYuTUulUMt2UjcmUo3/fIWvvPjBkZUcGb7i4oQFfnP6zY+BVTAsyYtQ/otA3/hWbFQ7m+j32eLtrMw4yF0ju5FYK3L/PXVslsTgjk14f/4WVGNzzUwTPU557JsS2w3rJXgUyfHz7+M7/vMVHM6J7hVqIve/aIwr3Yz10U9WY2AiW05+AX+duprerRvws36tvQ7nuF12Yirpe7KYuzH6awtMdMs4UNzt4LkrBnoYSWA0rFucpN774RIPIwk+S/Ii2Pf3nOkrvzXHplMxke2duZvZuu9o2C10fqzO6dOK5Nq1eM8GYJgIVnrAxbgo+AI2577h9G/bEIAvlmREdW27JXkRLLVJPXq3dpZ6Wrh5PwWF0ftGNdHtcE4+z01bx6mdmzKkazOvwwmIuonxnDugNZOWZnDgaF7VdzAmDPUbP9VX/t+NJ3sYSeDU9VsmFGD1zkMeRRJ8luRFuE4p9X3l8/81y8NIjDl2L3+3gT1HcvnD6B6IRH4tXpHLTkwlJ7+QzxZv9zoUY47b4I5NvA4hYAakNvKVRz/zvXeBBJkleRHu4fN6M6y7Mwhj6bYDUV3tbKJPdl4BE+Zv4ZXvNzCmb0v6+/3jjQZ92zSkR8tkmzPPRKRnvl5TYjuavoD9cWzPEtuvfL/Bo0iCy5K8CNeoXiLXnNLet/3RT9s8jMaYmnn5uw3c88ESjuQW8PuR3b0OJ+BEhMtOTGXptgOs2H7Q63CMqZFnvl7rK6+PoCXMqqN2rXg+ueU03/Y7czd7GE3wWJIXBfIKimvv7np/sYeRGFN9qsrTXzk1BWd2TynR9SCanD+gDYnxcUyYb7V5JnKUnmQ/PgoGQ5Xm32S7IfMI+VG43rQleVHgrB7NS2yXXkTamHA0a90eX/nf16R5GElwNU5KZGTvFny8cJutTmMiRo8/TfaVN0RZLV5Futz/pdchBJwleVEgIT6OJy/u59s+958zPYzGmOq56tW5vnJCfHT/K7r8xHYcOJrH1BU7vQ7FmBqLhimNKhJpa/DWVHT/Z40hl6alct85PQDYfsDWyzThLc+vWeSGIR09jCQ0Tu3clLaN6/Lej9HZ78dEl0Vb9vvKD53X27tAPDAvyiYvtyQviviPq129I3rn/TGR74kvV/nK94zu4WEkoREXJ1xyQiqz1u1hy94sr8MxplL+03Fdc0oH7wLxwKUv/eB1CAEVVkmeiNwmIstEZLmI3O7ue09EFrk/6SKyyO/4+0RknYisFpFRXsUdLi5NS/WVRz3znYeRGFO5V2ZuBKB+7VpR31Rb5OK0tojA+zYAw5iwsurh0SW2o2kqsrD57yoifYAbgJOA/sA4Eemqqpep6gBVHQB8CHzkHt8LuBzoDYwGnheR+HIfPEY0SUpk+u+HeR2GMZXaebC4O0G0zKBfHW0a1WVI1xTeX7DVVqcxYev0J6b5yj/cd5aHkYROnYSSqUPH+yZ5FEnghU2SB/QE5qhqlqrmAzOAC4puFGcWxkuBd91d5wH/U9UcVd0IrMNJEGNaUu1avvIt7/zkYSTGlG/E0zN85T5tGnoYSehdlpZKxoFsvl+72+tQjCkjr6CQrfuO+rZbNazrYTShNevekglttNTmhVOStwwYKiJNRaQeMAZI9bt9CLBTVYtmZ2wD+Ld7bHX3xbTkOsVJ3sQlGR5GYkxZuw5lcygnH4DHL+zrcTShN6JXcxrXS7A580xY+nRR7C6/16ZRyYT2prcXeBRJYIVNkqeqK4EngK+AycBiIN/vkJ9TXIsHUN6Y7jKpt4jcKCLzRWT+7t3R/+25TkI8p3cpXuDd5uUy4eT9+Vt95dF9WnoYiTdq14rngoFt+WrFTvbYfJYmzPzebzL9WJkbryJTlkfHdEdhk+QBqOqrqjpIVYcCe4G1ACJSC7gQeM/v8K2UrOlrC5T5GqKq/1bVNFVNS0lJCV7wYeScvsUfnl+vjI43qokOT01Z7Ss3qpfoYSTeuezEVPIKlI8X2hKEJnzklVrtIZrnxqtI6TnzcvMjfwWMsEryRKS5+7sdTlJXVHM3Alilqlv9Dv8MuFxEaotIR6ArMC+U8YarK05q5yt/sGBrJUca443Jtw/xOgTPdG+ZzIDURrz345ao6fdjIt+wp6b7yned3c27QMLIoxNXeB3CcQurJA/4UERWAJ8Dt6jqPnf/5ZRsqkVVlwMTgBU4zbu3qKq1TeIsin7FYCfR+27N7hKjGY3xysqMg75yj5YNPIzEe5edmMraXYdZ6DfprDFe2ra/eMDFrcO7ehiJt168apCv/MYPmzyMJDDCKslT1SGq2ktV+6vqN377r1XVF8s5/lFV7ayq3VU1+hadOw5/uaAv3/5+GIWKNQuZsDDm2e8BOLd/a48j8d64fq2omxDPhB9tAIbxntUoFxvdp1WJ7dnrMj2KJDDCKskzgdWxWRJp7RvzwYKtdhEbT+3PyqXoLXjDkE7eBhMGkuskMLZfKz5fvJ0jOflV38GYILrg+dm+crSv5VpTV7wyt+qDwpgleVHu4hPasm7XYRZvPeB1KCaGFdUmT/zd6fRtG1tz41XkshNTOZJbwMSlNtWR8dYi6zZQwpTbh5bYjuRKEkvyotyYfq2okxDHBwusWch4Q1V5d95m+rdtSO/WluAVSWvfmE4pSdZkazx1t9+0KcbRvWVyie3/zEr3JpAAsCQvyjWok8Do3i35bNF2mzPPeOKnzftYs/MwP/cb9W2cAVKXpqUyf9M+1u067HU4Jka97zcDQ+k1XGPZWT2a+8oPfRG5o2wtyYsBF5+QysHsfJszz3ji3XlbSEqM52c24KKMCwe1IT5ObAUM44mc/JJf/Euv4RrL/n7ZAK9DCAhL8mLAKZ2b0rphHZszz4TcgaN5fLFkO+cOaFNiXWXjaJ5ch7N6NOejn7aWmYzWmGA748npvvLS8SO9CyQMNaybUGI7UvvlWZIXA+LjhAsHtbU580zIfbpoG9l5hSUm6DYlXZaWSubhXL5ZucvrUEyM2eH3eZBcJ6GSI2PTUxf385Wv/c+PHkZy7CzJixEXndCWQnU+dI0Jle/W7KZTSpKNqK3EsO4pNE+ubU22JqRs6p6qXZJWvHLqjDW7PYzk2FmSFyM6NkuiWf1E0vdkeR2KiRHZeQV8vXIXjepaDUFlasXHcfEJbZm+ehc7DlhNuwmNi14onhtv8f9ZU21F7hhRvMTb375a42Ekx8aSvBiSeTiXd+ZuZs/hHK9DMTHg/H/NAuCnzfu9DSQCXJqWSqHChz9Zv1kTGqt2HPKVS/c/M8U6pST5ypOXRd6clpbkxaAVfmuIGhMs/h8ipnIdmiUxuGMTJszfQmFhZHbwNpEjK7e4qTZOPAwkAozrV7zM2ZqdkTfVkSV5MejqV+d5HYKJcv6Jyh9G9/Awkshx2YmpbNqTxZyNe7wOxUS5Xg9O8ZW//f0w7wKJACLCrWd18W0v2xZZq0dZkhdD7BubCZVrXy8eiXbzsM4eRhI5zunTiuTatWwFDBNS7ZsmVX1QjPP/HzbunzM9jKTmLMmLIRN/N8RXjtQ5f0xk+C5CR6J5qW5iPOcOaM2Xy3Zw4Gie1+GYKJVx4Kiv/McxVsteHfUSI3eOT0vyYkj7pvV85eXbrV+eCQ7//j6vX3eih5FEnstOTCUnv5DPbKojEySfL97uK9841GrZj4V/ohzuLMmLIf7fRsb9c2aZJW2MCYQ3f9jkKw/r3rySI01pfds0pEfLZN6zOfNMkPxl0iqvQ4hIk/xawk55bJqHkdSMJXkx5i8X9PWVL3tpjoeRmGj1+Jf2IXKsRIRL01JZtu0g6ZlHvA7HRJkDWcXdAL66Y6iHkUSeXq0blNiOlC5PluTFmMtPLJ7Be9GW/d4FYqLeb8/sUvVBpoxB7RsD8JHNmWcC7OIXiydA7toi2cNIIt8XSyJjzjxL8mJMXKkhtpHybcREhgK/qVNuG9HVw0giV9OkRACenbbO40hMtFm7K/LmeQsn/7pikK9867sLPYyk+izJi0F3nl28TMu5z83yMBITbfZl5QIw/me9SIi3fy/HIrlOcd9ZW+bMBEpeQaGv/Mktp3kYSeQa2buF1yHUmP0XjkE3nVE8ompphE3saMJbUVLSvEEdjyOJXEm1i5O8G9+a72EkJpo84ddXdkBqI+8CiWClv7h+5jdSOVxZkheDEmuVfNkPZtucXCYwpq/eBcDAdo28DSSC+X+QLNlqX8JMYLwycyMAQ7o28ziSyDbnvuG+8u8ioMnWkrwY9egFfXzlzXuyPIzERJOJS3dwQvvGtGpY1+tQItqD43r5yjbK1hyvfL+m2kvTUis50lSlZcPIaqWwJC9GXTm4va+8ea8leeb4bcw8wsqMg4zp26rqg02lrj21g6/8yswN3gViosK3q4tXoLHrM7ZUulaHiFxYjcfIVtVJAYrHeGDOhj124ZvjNmmpM6XAmL4tPY4k8vmPgn97zmYeOb9vJUcbU7kb3izu2xlvi5gft6tPbs9bczZVfWAYqGpBtpeBT4HK3hVDAUvyItDD5/fhT58s4/u1mV6HYqLAxCUZDGrXyJpqA+SmMzrz4oz1XodhIpz1uQ68h8/vw8Pn96n6wDBQVXPtl6r6S1W9rqIf4IdQBGoC76rB7aibEE9eQSGFhTZfnjl26ZlHWGFNtQF187DiUfA2MbI5VgP+PNVXnvfH4ZUcaaJRpUmeql5V1QNU5xgTnkSExy/qy9Z9R5mxZnfVdzCmApOWOU2151iSFzAN6ybQ211K6c4Jiz2OxkQq/+/vNrVR7KnWwAsRuUVEGvltNxaR3wQtKhMyY/q2okWD2lz3+o8csmp9c4wmLc1gQGoj2jSyptpA8p8z78BRuz5NzazdecjrEIzHqju69gZV3V+0oar7gBsCHYyI3CYiy0RkuYjc7rf/VhFZ7e5/0t3XQUSOisgi9+fFQMcTCxLi4zize3MAfv7yHI+jMZFo854slm07yFirxQu4py/p7ytvtKlUTA3d/cESX/nVX6R5GInxSlUDL4rEiYiou9CpiMQDiYEMRET64CSOJwG5wGQRmQi0Bc4D+qlqjog097vbelUdEMg4YlHdxHgAlm076HEkJhJNXFrUVGujagMttUk9X/n8f80i/fGxHkZjIs2iLft95eE9I29JLnP8qluTNwWYICLDReQs4F1gcoBj6QnMUdUsVc0HZgAXADcDj6tqDoCq7grw88Y8a2Izx2PS0gz6pzaibeN6VR9sauypi/t5HYIxJkJVN8n7AzANJ+G6BfgGuCfAsSwDhopIUxGpB4wBUoFuwBARmSsiM0TkRL/7dBSRhe7+IeU9qIjcKCLzRWT+7t02uKA8153W0Vd+64d07wIxEWfzniyWbjvAWKvFC5pe7uALgKzcfA8jMZFk/GfLfeW7zu7mYSTGS9VK8lS1EHgduF9VL1LVl1S1IJCBqOpK4AngK5xawsVAPk6TcmPgZOBunBpFATKAdqo6ELgTeEdEGpTzuP9W1TRVTUtJSQlkyFEjPk645IS2APzp0+VVHG1MMd+o2j7WHy9YerVqwGldmgLw06b93gZjIsbrs9N95VvO7OJdIMZT1R1dey6wCLeJVkQGiMhngQ5GVV9V1UGqOhTYC6wFtgIfqWMeUAg0U9UcVd3j3m8BsB6n1s8cgzv8vuk98/UaDyMxkWTS0gz6t21You+YCSwR4aWr04iPE+Zs2ON1OCYCxdkqFzGrus21/4czIGI/gKouAjoEOpiiQRUi0g64EKfv3yfAWe7+bjgDPjJFJMUdAIKIdAK6ArbI4zFq6Td/0jNfr/UwEhMptuzNYsnWAzYBcgjUr12LPm0aMnejJXmmahkHjvrKz1w2wLtAjOeqO7o2X1UPOK2kQfWhiDQF8oBbVHWfiLwGvCYiy3BG3f5CVVVEhgIPiUg+UADcpKp7gx1gtIqLE2rXiiMnvxCAwkK1b3+mUsVr1VqSFwond2rCazM3cjS3wDci3pjyzE/f5yuP7WfXZyyrbk3eMhG5AogXka4i8k9gdqCDUdUhqtpLVfur6jfuvlxVvUpV+7hNudPc/R+qam/32EGq+nmg44k1vxlW3G8j42C2h5GYSDBpaQb9rKk2ZE7u1JS8AuWnzfuqPtjEtP/9uNlXToiv7se8iUbVffVvBXoDOThNqAeB24MUk/HIb88qTvJue3ehh5GYcLdlbxaLrak2pNLaNyZOsH55plI5+QXMWue8R647rYO3wRjPVau5VlWzgPuB+91+cEmqalU9USber3l2/iarLTAV+9IdVTvGRtWGTHKdBPq2aWhJnqnUuGdn+sq3j7CxiLGuuqNr3xGRBiKSBCwHVovI3cENzXjhpjM6ex2CiQCTlu6gT5sGtGtqTbWhdHLnpizasp+Dts60qcDaXYd95YZ1EzyMxISD6jbX9lLVg8D5wCSgHXB1sIIy3hnUrpGv7K5iZ0wJ2/YfZdGW/dZU64FRvVuSV6BMXb7T61BMGPL/n21rSRuofpKXICIJOEnep6qaB1gGEIVG9i5euWDk37/zMBITrr50R9Xah0joDUxtRJtGdfliyXavQzFhaPb64qb8564Y6GEkJlxUN8l7CUgHkoDvRKQ9zuALE4WuPbUD4FT7Z+cFdGETEwUmLs2gd+sGtG+a5HUoMUdEGNevFTPXZrLvSK7X4Zgwc+Urc33lEEx5ZiJAdZc1e1ZV26jqGHXqgzcDZwY3NOOV8ef29pUv//ccDyMx4Wb7/qMs3GxNtV4a1681+YXKlOU7vA7FGBPmqjvw4i8i0shvVyNgfBDiMWFm0Zb9Xodgwsgka6r1XJ82DejQtB5fLMnwOhQTRrbtL17lYva9Z3kYiQkn1W2uPUdV9xdtqOo+YExQIjJhp6DQul8ax6SlGfRq1YAOzayp1itOk21rZq/PZPehHK/DMWHizdnpvnLrRnW9C8SEleomefEiUrtoQ0TqArUrOd5EuPE/6+UrPzJxhYeRmHCxff9Rftq835ZJCgPj+reiUGHyMqvNM46XvrOl201Z1U3y3ga+EZHrReSXwFfAG8ELy3jt2tM6+sofzN/qYSQmXHy5zOkDZv3xvNe9RTJdmtfnc2uyNaXMsqZa46fKJE+cITrvAo8APXGWN3tYVZ8McmwmTBzKyfc6BBMGJi3NoGerBnS0plrPiQg/69eaH9P3suOALT4U68Z/ttxXbt2wjoeRmHBTZZLnjqb9RFUnq+rvVfUuVZ0SgtiMx164cpDXIZgwkXHgKAs27WNs35ZVH2xCYlz/Vqg6U9qY2Pa62x+vTaO6NnWKKaG6zbVzROTEoEZiws7ZvVr4ylZbENu+XOo01Z5jTbVho3NKfXq2amATI8e4hZuL1xnv2qK+h5GYcFTdJO9M4AcRWS8iS0RkqYgsCWZgxnu14ovfHpOstiCmTVqaQY+WyXROsQ+RcDKuXysWbt7P1n1ZXodiPLJqxyFf+ZnLBngXiAlL1Z5CBegMnAX8DBjn/jZRbsEDIwD4YcOeKo400WrHgWzmb9pnAy7C0M/6tQZgog3AiFn3fbTUV25UL9HDSEw4qu6KF5tUdRNwFGfN2qIfE+Wa1q/NNae0Z+baTFviLEYVTdNhSV74ade0Hv3bNuRza7KNeaN7W39ZU1Z1V7w4V0TWAhuBGTjr2H4ZxLhMGDmrR3OO5hVYbV6MmrR0h2/KDhN+xvVrzbJtB0nPPOJ1KCbEJi8rXtruofN7V3KkiVXVba59GDgZWKOqHYHhwKygRWXCysmdmlIvMZ5vVu70OhQTYrsOZvPjpr1WixfGiiantgEYseemtxf4ys2TbeoUU1Z1k7w8Vd0DxIlInKp+CwwIXlgmnNRJiOf0Ls2YtnIXzow6JlZ8uWwHqjC2nzUFhavWjeqS1r6xrWUbY/z/F9usKaYi1U3y9otIfeA74L8i8g/AZsiNIcN7Nmf7gWxWZhyq+mATNSYuzaBbi/p0aZ7sdSimEuP6tWLVjkOs3WnXZ6zY5bdu8fz7R3gYiQln1U3yzgOygDuAycB6bHRtTDmzR3MApq2yJttYsetgNj+mW1NtJBjTtxUi2DJnMWTIk9/6yk3r21LypnzVHV17RFULVTVfVd9Q1Wfd5lsTI5on16F/24Z8vXKX16GYEJm83G2qtSQv7DVvUIeTOzbliyXbrUtFjMjNLwSgQ9N6HkdiwlmlSZ6IfFHVA1TnGBMdhvdsweKt+9nt10xgotfEJRl0bV6fri2sqTYSjO3Xig27j7B+t42yjXb//Gatr3xGtxQPIzHhrlYVt58uIp9VcrsAvQIYjwljZ/Vozt++WsO3q3dxaVqq1+GYINp1KJt56Xv53VldvQ7FVFPRaiTrdx+26W6i3NNfrfGV/3BODw8jMeGuqiTvvGo8Rm4gAjHhr3frBrRsUIf7P17K+QPakFirul06TaSZ4htVa021kSK5jvPv/NdvLSD98bEeR2NCpV5iVR/jJpZV+u5Q1RmhCsSEPxEhr6CQvALllZkb+M2wLl6HZIJk4tIMOqck0dVqhCJGUZJnoltBYXGfy1UPj/YwEhMJrCrG1MijF/QFYNa6TI8jMcGSnVfAvI17GdW7JWITcEWMlg2LJ8NdvcOmUolWh7OLZy+rkxDvYSQmEliSZ2rktC5NAZi1zgZXR6s1Ow9RqNCjVQOvQzE1ULtWPK3cRO+Jyas8jsYEy4bMw16HYCJIjZI8EUkQkYEi0jwYwYjIbSKyTESWi8jtfvtvFZHV7v4n/fbfJyLr3NtGBSMmU1KSX/8Pm6ohOp37nLNiYdOkRI8jMTVVNOBi2qpdHMmx+eqj0f0fL3N+j+npcSQmElQ1hcqLItLbLTcEFgNvAgtF5OeBDERE+gA3ACcB/YFxItJVRM7EGQDST1V7A391j+8FXA70BkYDz4uI1V0HWVxccfPdqzM3ehiJCQb/xN36eEWeP59bvEj9nsM2Ji4a5RY48+NddXJ7jyMxkaCqmrwhqrrcLV8HrFHVvsAJwD0BjqUnMEdVs1Q1H5gBXADcDDyuqjkAqlo0G+95wP9UNUdVNwLrcBJEE2R3nd0NgEcmrvQ4EhNoHy/c5isn1bYkL9J0SikeKDNxqa1+EW227sti3a7D3D+mJ3UTrU7DVK2qJM//q+DZwCcAqrojCLEsA4aKSFMRqQeMAVKBbsAQEZkrIjNE5ET3+DbAFr/7b3X3lSAiN4rIfBGZv3v37iCEHXv8P/wXbNrnYSQm0Pxfz84pNrI2Er11vfNd1/rlRZ9v3BWHhvcMSo8pE4WqSvL2i8g4ERkInIazbi0iUguoG8hAVHUl8ATwlfs8i4F8nGleGgMnA3cDE8QZ8lfesL8yncRU9d+qmqaqaSkpNjN4IFx5cjtfeceBbA8jMYH237mbAbhoUFuPIzHHqlXD4n/N9iUsuny9ciedmiWVqLE1pjJVJXm/Bn4L/Ae43a8GbzgwMdDBqOqrqjpIVYcCe4G1ODV0H6ljHlAINHP3+y+70BbYHuiYTFm1a8UzpGszAF6ducHjaEygHMrO85WfuKivh5GY49Heby3TV7636zNaHM7JZ+6GvVaLZ2qk0iRPVdeo6mhVHaCqr/vtn6KqdwU6mKJRuyLSDrgQeBenifgsd383IBHIBD4DLheR2iLSEegKzAt0TKZ8z1w2AICfNu/3NA4TOHdOWOwr14q32ZUiVUJ8nG+U7ZfLgtGzxnjh9VkbyS0o5KweLbwOxUSQqkbXNhOR/3OnMKkvIi+4U5x8KiLBWO7gQxFZAXwO3KKq+4DXgE4isgz4H/ALt1ZvOTABWIHTvHuLqhYEISZTjkb1iqfXeGN2uneBmID5asVOAHq0TPY4EnO8XrzqBF853x2NaSJXeuYR/jrVWa82rUNjj6MxkaSqr+vvALVxBj/MAzYAFwNfAK8EOhhVHaKqvVS1v6p+4+7LVdWrVLWP25Q7ze/4R1W1s6p2V9UvAx2PqVi831QqT09d7WEkJhD8p0651xY8j3hd/Jajs9q8yHftf4obqRKslt3UQFXvlhaq+kfgd0B9VX1KVVep6stAo6BHZ8La5Sc6XSIPZueXWE/RRJ4py4sTgTO62QClaDCoXSMAbn13obeBmOOWvicLgEvTbECUqZmqkrwCAHW+5pderNTaAGLc3aO6+8rb9x/1MBJzvB6dVDznoa1XGx0uTSsel3bYVr+IWP7N7Tbq3dRUVUleJxH5TEQ+9ysXbXcMQXwmjDWtX9tXvuiF2R5GYo7Xlr1Okn7+gNYeR2ICZXjP4g76ny+2iQci1e/fLx4QldahiYeRmEhU1ZT25/mV/1rqttLbJgbdPKwzL0xfz65DOV6HYo7R/qziOc+fvLi/h5GYQEpJLv4SdjTXxqRFqk8WFSfo/n2hjamOqqZQmaGqM4C5wB6cJtu5fvtNjBvbt5WvPH31rkqONOHqqlfn+sqJtaxTdzR585fO6hcPfbHC40jMscjOK07Oz+xufWVNzVU1hUotEXkSZ+LhN4C3gS0i8qSIJIQiQBPeOqUk+cr3fLDEw0jMsVq27SAAvx7ayeNITKCd2rmpr1xog6MizttzNvnK/74mzcNITKSq6mv7U0AToKOqnqCqA4HOOCNrrbnWUC+xuMX/wNG8So404e7kTk2rPshEFP9JrbPyrMk20jwysXhAlE2dYo5FVe+accANqnqoaIeqHgRuBsYEMzATOSb9bggAOfmFHLFRfBGlaO3ha05pz5k9bLmkaPTYhc4SdT/ZOrYRxb/m1UbVmmNVVZKn6j9LavHOAsDq/g0AvVo38JXfnbfZw0hMTX25LAOAa07p4G0gJmiGuX25pq2yPrORZOQz3/nKRWuFG1NTVSV5K0TkmtI7ReQqYFVwQjKRzL95wYS/SUsz6N4iucQKCSa6tGpYl7N7tWDq8h2U853dhKl1uw77yucPbONhJCaSVZXk3QLcIiLTReRpEfmriMzAWQHj5uCHZyLFG+4oPsA+SCLEzoPZzN+0jzF+I6RNdBrZqwXbD2SzfPtBr0Mx1WD/Q02gVDWFyjZVHQw8BKQDm4GHVPUkVd1WdJyI2IrJMe70LsXNCf5LZJnw9eXSDFRhbL+WXodigmx4zxbEiV2bkcI/GX/p6hM8jMREumoN11HVaar6T1V9VlW/KeeQ8vaZGOI/SefDX1iTbSSYtHQH3VrUp0vzZK9DMUHWJCmRkzo2YerynV6HYqph3D9n+sqjetuXMHPsAjUm26bhNtwwxFnpbputYxv2dh3M5sdNe62pNoaM7NWS1TsPkZ55xOtQTCWsqdYEUqCSPHtXGu47p6evvG7XoUqONF77ctkOp6nWkryYcXYvZy3br1ZYbV44+2H9Hq9DMFHEZlc0ARPn12R7zj++9zASU5WJSzPo2rw+XVtYU22sSG1Sj16tGli/vDC3+3DxOuBTbh/qYSQmGlhzrQmoB8Y6tXl5BWrNDmFq18Fsfkzfy9h+VosXa0b1bsmCzfvYfSin6oONJ2auzfSVu7e0L2Hm+AQqyRseoMcxEW5cv9a+8vQ1uz2MxFRk8nJrqo1VI3u3QBW+WWlNtuHq/QVbnd83neJxJCYaVJrkicghETno/hzy284SEd/6Vaq6N/ihmkjQsmEdX/njn7ZVcqTxysQl1lQbq3q0TCa1SV2mWr+8sDRxSYav3Kd1Qw8jMdGiqnnyklW1gfuTDLQGHgV2AP8IRYAm8ozq7XTw/mzxdo8jMaXtOpTNvHQbVRurRISRvVoyc20mh22d6bDzzNdrfOW6ifEeRmKiRbWaa0WkkYiMBxYDycCJqnpXMAMzkeuJi/r5yvkFhR5GYkqbUjSq1vrjxaxRvVuSW1DIjNXWnSKcqCpr3aXM0trb+gImMKpqrm0mIo8BPwH5wEBVfUBVbYy3qVCjeom+8uuz070LxJQxcWkGXZrXp5s11casE9o3pklSIlNX2CjbcHIkt8BXfvSCvh5GYqJJVTV5m4CfA28AWcD1InJn0U/QozMRq2iZs0cmrqSg0EbZhoPdh3KYt9GaamNdfJwwomdzpq3aRW6+1bSHi72Hc31lG1VrAqWqJO8p4D9uObmcH2PKdf3pHX3lxVv3exeI8Zm8fAeFNqrW4Kx+cSg7nzkbrFEmXHyyyBmo1qBOLY8jMdGk0neTqo4PURwmypzZo7mv/NFPWxnUzvqYeG3Skgw6pyTRrUV9r0MxHju9azPqJcYzdcUOhnZL8TocA/ztK2fQxavXnuhxJCaaVNUn79lSP/8QkT+JyOmhCtBErtQmdQF4e85mjyMxuw/lMHfjHsb2bYWIzV0e6+okxHNGtxS+WrGTQutO4bmNfusJt2lU18NITLSpqrl2Qamfn4DDwFMicntwQzOR7oGxvXzlRVv2exeI8TXVjrFRtcY1sncLdh7Mse4UYeCThdsQgW/uOoPWluSZAKqqufaN8vaLyIvAbOCZIMRkosTJHZv6yje+OZ9594/wMJrYNmlJBp1Skuhuo2qN66zuLagVJ0xdsZOB1p3CM6rKp4u2cWrnpnROsa4UJrCOaVkzVT0a6EAAROQ2EVkmIsuLagpFZLyIbBORRe7PGHd/BxE56rf/xWDEZI5dw3oJvvIuWyvTM5mHranWlNWwXgInd2rK1OU2lYqXFm3ZT/qeLM4b0MbrUEwUqnGSJyK1ROQ6YGsgAxGRPsANwElAf2CciHR1b/67qg5wfyb53W293/6bAhmPCYy/XdrfV953JLeSI02wTF7mNtXaqFpTysjeLVi/+wjr3El4Teh9snAbtWvFMbpPS69DMVGoJmvXHhSRQ8A24Bzg1wGOpScwR1WzVDUfmAFcEODnMCF24aC2vrJNvuqNSUsz6NQsiR4295Yp5exezhKEdm16I6+gkC+WZDCiZwsa1Emo+g7G1FBN1q5t4G63UNVLVTXQC5MuA4aKSFMRqQeMAVLd234rIktE5DUR8e880lFEForIDBEZEuB4TID8aZwzAGPiUvsgCbXMwznM2bCHMdZUa8rRqmFd+rdtyNTlO70OJSbNXJvJniO5nD/QmmpNcFRVk9deRBr6bZ/pTqNyp4gkVnbfmlLVlcATwFfAZJx1cvOBF4DOwAAgA3javUsG0E5VBwJ3Au+ISINyzuFGEZkvIvN377a1Gr1w/ekduemMzsxel8n+LGuyDaUpy62p1lRuZO+WLNqyn50Hs70OJeZ8vHAbjeolcIbNVWiCpKo+eROAJAARGQC8D2zG6TP3fKCDUdVXVXWQqg4F9gJrVXWnqhaoaiHwMk6fPVQ1p2gNXVVdAKwHupXzmP9W1TRVTUtJsQvJK2P7tiK/UK3GIMQmLc2gY7MkerayplpTvpFuk+1XK+zaDKXDOflMXbGDsX1bkVjrmMZAGlOlqt5Zdf2aZa8CXlPVp4HrcJOtQBKR5u7vdsCFwLsi4l8FcQFOsy4ikiIi8W65E9AV2BDomExg9GnTgNQmdZm4NMPrUGLGnsM5/LB+D2P6trSmWlOhLs3r07FZElNslG1ITV2+g+y8Qi6wploTRFUlef6fDGcB3wC4tWrB8KGIrAA+B25R1X3AkyKyVESWAGcCd7jHDgWWiMhi4APgJlXdG6S4zHESEcb0bcUsa7INmSnLd1pTramSiDCydwt+WL+HA0fzvA4nZnyyaDttG9flhPY2R6EJnqqSvGkiMkFE/gE0BqYBuLVrAf+kVtUhqtpLVfuralFCebWq9lXVfqp6rqpmuPs/VNXe7rGDVPXzQMdjAsuabENr0tIMOjStR69WZbqqGlPCyF4tyS9Upq/e5XUoMWHP4Rxmrt3N+QPaWC27CaqqkrzbgY+AdOB0VS36mtcSuD94YZlo1LdNQ2uyDZG9R3L5wUbVmmoamNqIZvVrM9X65YXE6p2HKFQY2K6R16GYKFfVFCqqqv9T1b+r6ja/m5YAzYIbmok21mQbOlOW76CgUK2p1lRLXJxwdq8WTF+1i+y8Aq/DiXpXvDwXgGSbG88EWVVTqDQQkftE5DkRGSmOW3EGOFwamhBNNPE12VqNQVBNWppB+6b16N3ammpN9Yzq3YIjuQX8sH6P16HEjPq1K10+3pjjVlVz7VtAd2Ap8CtgKnAxcJ6qnhfk2EwU6tumIW0b12XiEmuyDZa9R3KZvd6aak3NnNK5KQDXvf6j1bSHSLP6AZ1u1pgyqkryOqnqtar6EvBzIA0Yp6qLgh6ZiUoiwlhrsg2qqW5T7VhrqjU1ULtWvK/82qx07wKJcqoKQK04oXmDOh5HY6JdVUmebzy9qhYAG1X1UHBDMtFubD9rsg2miUszaNfEmmpNzfV0R2I/+81ajyOJXlm5Tp/Hu0d19zgSEwuqSvL6i8hB9+cQ0K+oLCIHQxGgiT5FTbaTbJRtwO2zplpzHP77q8G+clZuvoeRRCdV5akpqwFnEmpjgq2q0bXxqtrA/UlW1Vp+ZasmMMfEv8n2QJZNvhpIU1dYU605dk2SivuI9XpwioeRRKfXZqXz+ux0fnlaR87q0dzrcEwMsAXzjCfG9G1FXoEydYUtpRRIE5fuILVJXfq0se9g5tiM7WdfEIJh054jPPzFCro0r88DY3taTbsJCUvyjCf6tXVH2VqTbcDsz8pl9rpMa6o1x+Wak9v7yt+ushUwAuX+j5cBsHVfFnFxdn2a0LAkz3jCf2Jka7INjKnLd5JvTbXmOA3u1NRXnjB/i4eRRJdNe48AkJ0XrKXfjSnLkjzjmbHWZBtQE5dm0LZxXfq2aeh1KCZKLNy83+sQokJWbj4Z+7MBOKljE4+jMbHEkjzjmaImWxtle/z2Z+Uya10mY62p1gTAhzefCsCOg9nsOZzjcTSR78f0feQXKg+d15tXf5HmdTgmhliSZzxT1GQ705psj9vUFU5Tra1VawLhhPaNi8uPfO1hJNFh9vpMEuKFS05ItfVqTUhZkmc8ZaNsA2OS21Tbr6011ZrAs9q84zN73R4GtmtM3cT4qg82JoAsyTOe6t+2IW0aWZPt8TiQlccsG1VrAuzNX57kK//7uw0eRhLZ9mflsmz7AU7r3MzrUEwMsiTPeEpEGNvPmmyPx9QVO8grsKZaE1hDu6X4yi9ZknfM5mzYgyqc1qVp1QcbE2CW5BnPndEthbwCZem2A16HEpEmLc2gTaO69LemWhNgQ7oW1z5t23/Uw0gi1+z1e6iXGE//1EZeh2JikCV5xnOtG9UF4K056d4GEoH2HM7h29W7GdO3pTXVmoD76yX9fWWbGPnYzFqXyUkdm5AQbx+3JvTsXWc816JBbQCmLN/pcSSRZ+iT3wLQrUWyx5GYaNSiQR1f+ZGJKzyMJPJk5xXQ4d6JrN99xPrjGc9Ykmc8Vy+xlq88ZbmNsq2JI7kFADZqzwTNi1edADgrNWTnFXgcTeQ4/YlpvvKp1h/PeMSSPBMWmic7tXm/fmuBx5FEjm9XFzefHcrO9zASE81G92npK9sKGNWXeTjXV+7ZsoGHkZhYZkmeCQu/PqOzr7xi+0EPI4kcny7c5iuf2MGWSjLB878bT3Z+/7jZ40giw74juSW24+Ksv6zxhiV5Jixce2oHX/nGt+Z7F0gE+WTRdl+5S/P6HkZiot3JnZpyzSnt+XLpDvaWSmBMWY9/ucpXfurifh5GYmKdJXkmLMT7fdPdus+maqjKC9PX+8r//dVgDyMxseKKwe3ILSjkgwVbvA4lrBUUKu/NL/4b/ax/aw+jMbHOkjwTNu4e1d1XvviF2R5GEv6emFxcU3BqZ+vUbYKvR8sGnNC+Me/O24Kqeh1O2Fqydb+vPLZfK+ok2KAo4x1L8kzYuP70jr7y/E37PIwkstj8eCZUrjipHRszj/DD+j1ehxK2Lni++Avq037zDBrjBUvyTNgo/Y3XpmuoWBt3Aumpdwz1OBITS8b2a0XDugn8d54NwKgOq8UzXrMkz4SVpeNH+so9/jTZw0jCW7P6iQztlmKTIJuQqpMQz0WD2jJ1+Q4yD+d4HU7Ymb0+01d+cFwvDyMxxhFWSZ6I3CYiy0RkuYjc7u4bLyLbRGSR+zPG7/j7RGSdiKwWkVGeBW4CJrlOgtchhL3svAI2782iWf1Er0MxMeiKwankFSjvz9/qdShh54qX5/rK/jMGGOOVsEnyRKQPcANwEtAfGCciXd2b/66qA9yfSe7xvYDLgd7AaOB5EbG68SiQEF/cx+yblbbUWWnv/biFfVl5XDyordehmBjUpXkyJ3VswrvzNlNYaAMwipQejGJz45lwEDZJHtATmKOqWaqaD8wALqjk+POA/6lqjqpuBNbhJIgmwq14aLSvfP0bNmeev5z8Al6Yvp6TOjThFBtVazxy5eB2bN6bxcx1mVUfHCP+OW2dr/zS1Sd4GIkxxcIpyVsGDBWRpiJSDxgDpLq3/VZElojIayLS2N3XBvCfsGmru89EuIT4km9Lqy0oNuHHLew4mM1tI7raqFrjmdF9WtK4XgLvzLUBGEX+9tUaX3lU75aVHGlM6IRNkqeqK4EngK+AycBiIB94AegMDAAygKfdu5T3CVcmGxCRG0VkvojM3717dxAiN8Hww31n+cqd/jjJw0jCR05+Ac9PX09a+8Y2N57xVO1a8VySlspXK3ey62C21+GElT+O6eF1CMb4hE2SB6Cqr6rqIFUdCuwF1qrqTlUtUNVC4GWKm2S3UlzTB9AW2E4pqvpvVU1T1bSUlJRgn4IJkFYN63odQtiZMH8rGQesFs+Eh5+f1I6CQmXCfFsB47x/zfKVT+3czMNIjCkprJI8EWnu/m4HXAi8KyKt/A65AKdZF+Az4HIRqS0iHYGuwLxQxmuCq0WD2r7yA58s9TAS7+XkF/DCt+sY1K4Rp3exDxHjvY7Nkji1c1PenbeFghjuUpGbX8jiLfsBGNGzBX3aNPQ2IGP8hFWSB3woIiuAz4FbVHUf8KSILBWRJcCZwB0AqrocmACswGnevUVVbfbcKDL3jyN85bfnxHbfnw8WbGX7gWxuG9HNavFM2LhicDu27T/Kd2tityvM7e8t9JXTOjSu5EhjQq+W1wH4U9Uh5ey7upLjHwUeDWpQxngsN7+Q579dz4DURgztarV4JnyM7NWSZvUT+e/czZzZo7nX4YRcdl4Bk5bu8G1ffXJ7D6Mxpqxwq8kzpoR59w/3lV+budHDSLzz4U9b2bb/qPXFM2EnsVYcl6SlMm3VTjIOHPU6nJDbfah41Y9VD48mqXZY1ZsYY0meCW/Nk+v4yg99sYIjOfkeRhN6eQWF/OvbdfRv25Bh3WzgkAk/Pz+xHYXqTNIda579Zq2vbOvUmnBkSZ4Je3+/rL+v/NPmfR5GEnof/bSVrfusFs+Er3ZN6zGkazP+N28L+QWFXocTUu8vcJZ262uDLUyYsiTPhL0LBhYv33X1q7EzgDqvoJDnvl1Hv7YNObN77PV3MpHjysHt2XEwm29Xx84ADP+E9m+X9q/kSGO8Y0meiQh3j+rudQgh9/HCbWzZe5TfnWW1eCa8De/ZnObJtXl3XuyMgu9y/5e+cpOkRA8jMaZiluSZiHD96R195YUx0GSb7/bF69OmAcN7Wi2eCW8J8XGM7N2Caat2kZ0X/TNZ/emTZSW2m9avXcGRxnjLkjwTEeokxNOygTMI45UYGGX7yaLtbNqTZbV4JmIMSHXmiBv9zHceRxJchYXKW3M2+bYXPDCikqON8ZYleSZiTL97GCd1bMKUZTuierqG/IJCnpu2ll6tGnB2rxZeh2NMtYzs7bxX0/dksedwThVHR67PlxSvnlk3Id5q8UxYsyTPRIw6CfE8fUl/ClV564dNVd8hQn22eDvpe7JsRK2JKA3qJPjKczbs9TCS4CksVG773yLf9tLxI70LxphqsCTPRJTUJvU4u1cL3pm3maO50df3p6BQeW7aOnq2asBIq8UzESa5jjMZ8C3v/EReFE6n0uPBySW2a8XbR6gJb/YONRHnl6d1ZH9WHh8v3OZ1KAH3+eLtbMg8wm3Du1gtnok439x5hq/86MSVHkYSHLn5xYnrF7ee7mEkxlSPJXkm4pzUsQm9WzfgtVkbUVWvwwmYgkLl2Wlr6dEymZG9WnodjjE1lpJc3D/t9dnp3gUSBE9PXV1iu49NgGwigCV5JuKICNef3pF1uw7z/PT1XocTMF8s2c6G3Uf43fCuxMVZLZ6JPCLCW9ef5Nt+64d074IJsH9OW+crn9PHvoSZyGBJnolIY/u1AuCpKaujYqmzgkLln9PW0b1FMqN72weIiVxDuhavsfynT5d7GEngzF6XWWL7+SsHeRSJMTVjSZ6JSLVrFS8G/vni7ZUcGRkmLs1g3a7D3Dq8i9XimYh361ldfOWPF271MJLjl5NfwBWvzC2xz/rLmkhhSZ6JWE9f4qwX+fnijIheGL2wUPnnN2vp2rw+Y/q08jocY47bLWcWJ3l3vLeYgsLI7Tv7j6/Xlth+78aTPYrEmJqzJM9ErItOaMtLV59A5uEcvl+bWfUdwtSkZRms3XWYW60vnokSdRLiGf+zXr7t/Vm5HkZz7AoKtUy/38GdmnoUjTE1Z0meiWhndm9O43oJ3DlhUUSuaVtYqDz7zVo6pyQxtq/V4pno8YtTO/jK+7LyvAvkOHT+46QS2/NtCTMTYSzJMxEtsVYc5/Zvzb6sPC54frbX4dTY5OU7WLPzML8b3pV4q8UzUUREePi83gBMX73L42iO32MX9qWZLWFmIowleSbiXTiordchHJOiWrxOKUmM69fa63CMCbirT+nAKZ2a8tJ3G/h6xU5y8iNnlZoO904ssX35iakeRWLMsbMkz0S8fm2LJyWNpMmRe/3fZFbtOMStZ3WxWjwTtX43vCu7D+Xwqzfn0/2ByRGR6J306Ncltjc+NsZG1JqIZEmeiXgi4uvkff0b85mxZrfHEVXt88Xbyc5zRgT/zGrxTBQ7uVOTEtvn/ON7jyKpvl2HckpsW4JnIpUleSYqnD+wDfVr12Laql384rV5YT9B8q3vLvSVbZFzE81EhLtHdfdtb9h9hHW7DnsYUeVKN9P2adPAo0iMOX726WKiQqN6iVx/ekff9oXPzw7bptvZ60tO9xKucRoTKP7z5gGM+NsMjyKp3IJNe8vs+/g3p3kQiTGBYUmeiRq/PavkB0nH+yZVcKS3rnjZZs83sWf1I6NLbP/tqzUeRVK+7LwCLnrhhxL75v5xOAlW024imL17TdRIiI9jxt3DSux7M8wXSP/u7jO9DsGYkPBfihDg2W/WMmtd+Exi3uNPk0ts92rVgBYN6ngUjTGBYUmeiSrtmyaV2H7w0+XsOZxTwdHe+ssFfWnXtJ7XYRgTMmsfPafE9pWl1oT1ynPT1pbZN+m2IR5EYkxgWZJnos6HN59SYvuER76u4MjQ86+5GNffVrgwsSUhPo57Rncvse+0x6d5FI0jKzefv04t2XS86uHRFRxtTGSxJM9EnRPaN+G24V1L7Cs9Ys4Lb83ZxDWvzaNbi/rMuvcsGtRJ8DokY0LuN8NK9p3dtv+op7XtvR6cUmZfnYT4co40JvKEVZInIreJyDIRWS4it5e67fcioiLSzN3uICJHRWSR+/OiJ0GbsHTH2d3K7Fuz85AHkUB+QSH/9+ky/vTJMs7olsKHN59Km0Z1PYnFmHCw5pGSzbYnPPI1+QWFIY/jjvcWldmX/vjYkMdhTLCETZInIn2AG4CTgP7AOBHp6t6WCpwNbC51t/WqOsD9uSmkAZuwV7r/z8i/f8eBo6FdKP3A0Tyue/1H3vhhEzcM6cjL16SRbDV4JsYl1orjhSsHldj367cWhDSG3Ydy+HjhtpA+pzGhFjZJHtATmKOqWaqaD8wALnBv+ztwD2ATiplqS4iPo3uL5BL7zg7h/FzpmUe44PlZzNmwhycv6sf9Y3vZ8mXGuM7p24pRvVv4tr9ZtYuHv1gRsuc/8dGyfXWtFs9Em3BK8pYBQ0WkqYjUA8YAqSJyLrBNVReXc5+OIrJQRGaIiA2FMmVMuWMop3dp5tsuvVxRsMxen8l5/5rFviO5vH39YC61xc2NKeOlq9NKbL86cyPz08tOSBxIhYXKr96YX2b/hr+MCerzGuOFsEnyVHUl8ATwFTAZWAzkA/cDD5ZzlwygnaoOBO4E3hGRMuvPiMiNIjJfRObv3h3+a5qawHv7V4MZkNooZM/38cKtXPPqPJon1+bTW05ncKemIXtuYyJN6eTq4hd/YNv+o0F5LlXlmW/W8vXKnSX2//dXg4mzWnYThcImyQNQ1VdVdZCqDgX2AulAR2CxiKQDbYGfRKSlquao6h73fguA9UCZ3vaq+m9VTVPVtJSUlFCdigkzn9xyGsN7NAfgvo+Wsm5X4Adh5OQX8OCny7jjvcX0at2AD39zqs2DZ0wV4uKE1687scS+0x6fFpQRt8OfnsGz35SdE+80v9p+Y6JJWCV5ItLc/d0OuBB4U1Wbq2oHVe0AbAUGqeoOEUkRkXj3+E5AV2CDR6GbCPDUJf25YnA7PvppKyP+9h3X/Wces9ZlBmTt2Oy8An791gLe/GETAD1bNrApUoyppmHdmzP+Z71K7LvmtXkczA7cQKmCQmVD5pES+4b3aG798ExUC6skD/hQRFYAnwO3qOq+So4dCiwRkcXAB8BNqhrczhwmojVJSuQvF/Rl9r1ncefZ3Vi67QBXvjKXc/7xPR8s2Epu/rFN4ZCVm88vX/+RGWuKuwMUBiBxNCaW/OLUDlx1cjsAGtZNYM3OQ/zyPz+SlZt/XI+7cPM+7pqwmC17s8rc9sJVJxzXYxsT7iQQtRiRIi0tTefPL9vh1sSm7LwCPlu0nVdmbmDNzsOc1qUpj13Qr0ZNrIey8/jl6z+yYNM+/npJf+6c4IwP+vlJ7Xjswr7BCt2YqDV5WQandG7GzLWZ3PruT5zWpRmv/CKtzNq31fHYpJW89J3TwNOsfiKZh3N9t63/yxgb7W7CnogsUNW0qo8sX7jV5BkTMnUS4rn0xFSm3D6Uu87uxqx1exj61Ldc9595rK1g4uS8gkLW7jxE2iNf8/TU1Vz96jx+2ryfZ38+kAsHtaVJUiIA15/eIYRnYkz0GN2nFQ3rJjC2XyueuKgf36/N5NZ3FvomS1636zDfrdnNroPZqCp5BYWs23XI1+3ile83MHFJBqrqS/AAX4JXNyGeGXcPswTPxASryTMGZ2mlojU06yXGk5Vb4Lvt75f1p1/bRox7diZH8wpK3C8hXnjuikGM6t0ypPEaEyten7WR8Z+v4IKBbfjFqR04/1+zKjx29SOj6f7AZADO7J7Ct6uLu1B0bJZE28Z1ef7KQTYhuYkYx1uTVyuQwRgTqZrUS/SV/RM8gDveK2+KRse/r07jTHfUrjEm8K49rSNHcgt4aspqpq/eVemxr3y/EYBTOzdl7saSXbTvGtmNcf1aBy1OY8KRNdcaA9RNdPr7XHdah2rfp0fLZEvwjAmB3wzrzE1ndGZfVuWjbZ+ashqAq05uz1d3nuHb//09Z1qCZ2KS1eQZ49r4mDMp63/nbq7WSNvXrzsp2CEZYwAR4Q+ju5OdV8Drs9M5uVMT5mwoO5nC2b1acHavFpzVozl1EuJtehQT86xPnjGl7M/KZcBDX5XYJwJFl8rD5/VmcKemdCu1Lq4xJvgKCxURWL3zEPd8sIQlWw/4blvwwAia1q/tYXTGBNbx9smzJM+YKmQezqGZfXAYY4wJMZtCxZggswTPGGNMJLIkzxhjjDEmClmSZ4wxxhgThSzJM8YYY4yJQpbkGWOMMcZEIUvyjDHGGGOikCV5xhhjjDFRyJI8Y4wxxpgoZEmeMcYYY0wUsiTPGGOMMSYKWZJnjDHGGBOFYmrtWhHZDWwKwVM1AzJD8DxesnOMLrFyrnae0cXOM7rEwnnW9Bzbq2rKsT5ZTCV5oSIi849nQeFIYOcYXWLlXO08o4udZ3SJhfMM9Tlac60xxhhjTBSyJM8YY4wxJgpZkhcc//Y6gBCwc4wusXKudp7Rxc4zusTCeYb0HK1PnjHGGGNMFLKaPGOMMcaYKBTzSZ6IpIrItyKyUkSWi8ht7v4mIvKViKx1fzd29zd1jz8sIs+VeqzLRGSJ+zhPVvKcJ4jIUhFZJyLPioi4+4eKyE8iki8iF0fjOfrdfrGIqIgEbJRROJ2niPxdRBa5P2tEZH+gztPDc31URLaIyOFS+2uLyHvu32CuiHTw8DzPFpEF7muyQETO8nusSt+TVR0XrOsz3M7T7/aAXqPhdI5heH0G4lwj/fosN/5ynjPcPz+Ddo5+t1f/2lTVmP4BWgGD3HIysAboBTwJ3Ovuvxd4wi0nAacDNwHP+T1OU2AzkOJuvwEMr+A55wGnAAJ8CZzj7u8A9APeBC6OxnP0i+E7YA6QFq3n6XfMrcBrUfC+Pdl93sOl9v8GeNEtXw685+F5DgRau+U+wLaavFZeXJ/hdp5+MQT0Gg23cwyz6zMQ5xrp12e58Vf3NSV8Pj+Ddo5+MVT72gzYmzpafoBPgbOB1UArvxd5danjrqXkh+WJwNd+21cDz1fwhlnlt/1z4KVSx7weyDdpuJ0j8AwwDphenTdppJ6n3/7ZwNmR/L4t9RilP0SmAKe45Vo4E32Kl+fp7hdgD1C7Bq+V59dnOJxnKK5Rr8/Rb3/YXJ/Heq6l7h9x12dl8df0NQ329en1Odb02oz55lp/bjX2QGAu0EJVMwDc382ruPs6oIeIdBCRWsD5QGo5x7UBtvptb3X3hYTX5ygiA4FUVf3iOE6jSl6fp18c7YGOwLSan0X1hOhcK9MG2OI+Zz5wAKeGMKCO4TwvAhaqag7Vv+48vT7B+/MMxTXq9Tn6xRGO1+exnGtlIuH6rK5I+/wM6Dkey7VZqwZPHNVEpD7wIXC7qh6soNtDhVR1n4jcDLwHFOJ8O+xU3lOVd/cahntMvD5HEYkD/o5TmxQ0Xp9nqe3LgQ9UtaBGQVRTCM+10jDKe+gaPkblT1DD8xSR3sATwMgaxujZ9Qnen2corlGvz7HUdlhdn8dxrpU+bAAeo/InOP7zrPZTlbMvLD8/A32Ox3ptWk0eICIJOC/ef1X1I3f3ThFp5d7eCthV1eOo6ueqOlhVT8Gpyl0rIvFS3MH3IZysvK3f3doC2wN5PuUJk3NMxumjMF1E0nH6J3xWrc6j1RQm5+nvcuDd4zur8oX4XCuzFbf2z60NbAjsPbazKqum5ykibYGPgWtUdb1fjGVeq3C5Pt24w+E8g3qNhsk5+gub6/M4z7UykXB9VvTYYXF9hsk5Htu1GYw260j6wcma3wSeKbX/KUp2qnyy1O3X4te3yd3X3P3dGFgEdKvgOX90X6CiTpVjSt3+OoHtOBp25+geM53ADrwIq/MEugPpBKH/ixfn6nd86T4/t1CyY/cEr84TaAQsBi46lvdkdY4L9PUZrufpHhOwazTczjGcrs9AnKvf8RF5fVYUf03/HoG+PsPxHN1jqnVtBvSNHYk/OCMOFViC8wG3CBiD02fhG2Ct+7uJ333Scb4JHcbJunu5+98FVrg/l1fynGnAMmA98FzRPxmcTvBbgSM4nTWXR9s5HsubNFLPExgPPB5F79sn3fsVur/Hu/vrAO/j9O+bB3Ty6jyBB3Cun0V+P0VJbJXvSS+uz3A7z2Bdo+F2joTR9Rmgc43067Pc+L2+PsPpHI/l2rQVL4wxxhhjopD1yTPGGGOMiUKW5BljjDHGRCFL8owxxhhjopAlecYYY4wxUciSPGOMMcaYKGRJnjHGGGNMFLIkzxhjjDEmClmSZ4wxxhgThSzJM8YYY4yJQpbkGWOMMcZEIUvyjDHGGGOikCV5xhhjjDFRyJI8Y4wxxpgoZEmeMcYYY0wUsiTPGGOMMSYKWZJnjDHGGBOFLMkzxphyiEiBiCwSkeUislhE7hSRSv9nikgHEbkiVDEaY0xlLMkzxpjyHVXVAaraGzgbGAP8XxX36QBYkmeMCQuiql7HYIwxYUdEDqtqfb/tTsCPQDOgPfAWkOTe/FtVnS0ic4CewEbgDeBZ4HFgGFAb+JeqvhSykzDGxDRL8owxphylkzx33z6gB3AIKFTVbBHpCryrqmkiMgz4vaqOc4+/EWiuqo+ISG1gFnCJqm4M5bkYY2JTLa8DMMaYCCLu7wTgOREZABQA3So4fiTQT0QudrcbAl1xavqMMSaoLMkzxphqcJtrC4BdOH3zdgL9cfo2Z1d0N+BWVZ0SkiCNMcaPDbwwxpgqiEgK8CLwnDp9XBoCGapaCFwNxLuHHgKS/e46BbhZRBLcx+kmIkkYY0wIWE2eMcaUr66ILMJpms3HGWjxN/e254EPReQS4FvgiLt/CZAvIouB14F/4Iy4/UlEBNgNnB+a8I0xsc4GXhhjjDHGRCFrrjXGGGOMiUKW5BljjDHGRCFL8owxxhhjopAlecYYY4wxUciSPGOMMcaYKGRJnjHGGGNMFLIkzxhjjDEmClmSZ4wxxhgThf4fW8NwAttRfjEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,5))\n",
    "ax.plot(t_rec_hmi,keyshmi['RSUN_OBS'])\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('RSUN_OBS [arcsec]')\n",
    "ax.xaxis.labelpad=15\n",
    "ax.yaxis.labelpad=15\n",
    "ax.set_title(\"Angular Radius of the Sun over time as observed by SDO\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4be1e9f9-e45d-4bb7-9ca8-f555bb0db0c8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The keyword RSUN_OBS varies by 3.3429064757143063 percent.\n"
     ]
    }
   ],
   "source": [
    "print(\"The keyword RSUN_OBS varies by\",(np.nanmax(keyshmi['RSUN_OBS'])-np.nanmin(keyshmi['RSUN_OBS']))/np.nanmax(keyshmi['RSUN_OBS'])*100,\"percent.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7428dae1-5813-45dc-b441-d6ba067b5e32",
   "metadata": {},
   "source": [
    "Plot variation in DSUN_OBS over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "eec51568-a8ed-4163-9529-bcaf1bcde2e6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Distance from SDO to Sun')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnwAAAFYCAYAAADJOCHrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABXUElEQVR4nO3dd3zU9f0H8Nc7m5CQBAgQZth7SVgqQ5YGtGrVFq22WqvitrZaR93bVltn1f60jlbcmw2iqIDsEQgzrEBIwgokkP3+/XGXyyUkuYy7+3zve6/n45EHn7v75u71TXLkne9niaqCiIiIiOwrxHQAIiIiIvItFnxERERENseCj4iIiMjmWPARERER2RwLPiIiIiKbY8FHREREZHMs+IioSUTkNRF5wHSOxhKRG0UkW0TyRaSV6TxERL7Ago+IaiUiu0XklIicEJFjIrJURGaIiOv/DlWdoaqP1fO5Jvk2ccOISDiA5wFMUdUYVT3s59ePEJHnRCTTWXDuEpF/uD3u8evvPO5MEfnWeVyeiHwtIv3qeN2rReTHJuSOF5G3ROSg8zW3ichfGvt8ROR7LPiIyJMLVDUWQBcATwP4C4A3zUbymrYAogBsqulBEQnz8evfCyAFwAgAsQDOAbC22jF1fv1FZDSA+QC+BNAeQFcA6wH8JCLdfJT7HwBiAPQFEAfgFwB2+ui1iMgLWPARUb2oap6qfgXg1wB+JyIDAEBE3haRx53t1iLyjfNq1BER+UFEQkTkPQCdAXztvJJ1t/P4j51XifJEZImI9K94PefzviIis5xXkX4Wke5uj/cXkQXO18kWkfuc94eIyD0islNEDovIRyLSsvr5iEgvAFudN4+JyLfO+1VEbhaR7QC2O++7TkR2OF/rKxFp7/Y8KiI3ich2Z87HRKS7iCwTkePO14+o5cs6HMDnqnpAHXar6rsN+foDeBbAu6r6gqqeUNUjqvpXAMsBPFzDefcF8BqA0c7vxTHn/XEi8q6I5IrIHhH5a/UridVyv6+qR1W1XFW3qOonzudJdn5NXMWyiHwnIn9wtq8WkR9F5O8ictR5VTO1ltchIi9hwUdEDaKqKwBkAhhTw8N/cj6WCMfVs/scn6JXAdgLx9WqGFV91nn8HAA9AbQBsAbA/6o93+UAHgGQAGAHgCcAQERiASwEMBeOq1o9ACxyfs5tAC4CMM752FEAr9RwHtsAVBSY8ao6we3hiwCMBNBPRCYAeArArwAkAdgD4INqT3cegGEARgG4G8AbAH4DoBOAAc7zqMlyAHc6C8aBIiK1HOee2/X1F5FoAGcC+LiGQz8CMLmGz08HMAPAMuf3It750EtwXK3rBsfX7rcArqkj9xMico2I9PSUuQYj4Si2W8NRsL5Zn3MnosYL+oLPOQ4lR0TS6nHsWBFZIyKlInJptcfmOq9qfOO7tESWcQDAaVfNAJTAURR1UdUSVf1B69iwW1Xfcl6VKoLjatRgEYlzO+QzVV2hqqVwFINDnPefD+Cgqj6nqoXO5/jZ+dgNAO5X1Uy35720gd2zTzmvlJ2Co3B7S1XXOJ/vXjiujiW7Hf+Mqh5X1U0A0gDMV9UMVc2Do6gdWtvrAHjG+RqrAOwXkd/VI1/F178lHP+PZ9VwTBYcBZVHIhIKx5XDe51fy90AngNwVS2fcisc349bAGx2Xv1syFW6Par6b1UtA/AOHD8zbRvw+UTUQEFf8AF4G46/zutjL4CrAbxfw2N/Q+3/ORLZTQcAR2q4/29wXImbLyIZInJPbU8gIqEi8rSz6/U4gN3Oh9yLlINu7ZNwjBsDHFfOahsz1gXA584/wI4BSAdQhoYVFPvc2u3huKoHAFDVfACH4fgaVMh2a5+q4XYMaqCqZar6iqqeBSAejiuYbzm7XetS8fU/CqAcjoKpuiQAhzw8T4XWACLgdp7OdoeaDlbVU6r6pKoOA9AKjquJH9fUdV4L1/dVVU86mzV+jYjIO4K+4FPVJaj2i8s5/mauiKx2jkHq4zx2t6pugOM/2OrPswjACb+EJjJIRIbDUQicNsvTeXXoT6raDcAFcHRXTqx4uNrhVwC4EMAkOLoSkyteoh4x9gHoXsdjqaoa7/YRpar76/G8rlNxax+Ao4h0hBNpDkeR05Dn8/yCjiLqFTiKuLpm2Lq+/qpaAGAZgMtqOPRXqOzmPu3lqt0+BMfV2S5u93VGPc5RVY8DeBJAczgmjBQ4H4p2O6ydp+chIt8K+oKvFm8AuNX51+ufAbxqOA+RcSLSQkTOh2P82n9VdWMNx5wvIj2c47GOw3Flrcz5cDYc48MqxAIoguNqWTQcRUN9fQOgnYjcISKRIhIrIiOdj70Gx/iyLs5MiSJyYQOeu7r3AVwjIkNEJNKZ82dnt2eTOPOPF5FmIhLm7M6Nxekzdev6+t8DxySO25xfhwRxTKIZDcf4x5pkA+hYMZnE2bX6ERxft1jn1+5OAP+tJfcDIjJcHMvKRAG4HcAxAFtVNReOQvFK51Xc36P24pyI/IQFXzUiEgPnIGgRWQfgddTcXUIULL4WkRNwXDm7H45162obzN8TjskU+XBceXpVVb9zPvYUgL86u1r/DOBdOLoN9wPYDMdEgHpR1RNwTEi4AI7uwe1wLGkCAC8A+AqObuUTzucdWdPz1PO1FgF4AMCncIyL6w5gemOfr5pTcIyVOwjHVbabAVyiqhlux9T59VfVHwGcC+CXznx74BgzeLaqbq/ldb+FYymagyJS0e17KxxX5zLguHr7PoC3avl8BfAfZ+YDcHwvpjm7uwHgOgB3wVHM9wew1NMXgoh8S+oYTx00nIOvv1HVASLSAo6/Umst8kTkbefxn1S7fzyAP6vq+b5LS0RERNQwvMJXjXM8yi4RuQwAxGGw4VhEREREjRb0V/hEZCaA8XDMUssG8BAc3R3/gqMrNxzAB6r6qHOw9OdwrAlWCMeyEP2dz/MDgD5wzDQ7DOBaVZ3n37MhIiIiOl3QF3xEREREdscuXSIiIiKbY8FHREREZHMN2WrIVlq3bq3JycmmYxARERF5tHr16kOqmtjYzw/agi85ORmrVq0yHYOIiIjIIxHZ4/mo2rFLl4iIiMjmWPARERER2RwLPiIiIiKbY8FHREREZHMs+IiIiIhsjgUfERERkc2x4CMiIiKyORZ8RERERDbHgo+IiIjI5oJ2pw2ihlq/7xjKVdG1dXNER4QhIox/LxERUWBgwUdUh/JyRbf7Ztd5zEMX9MM1Z3X1UyIiqrAj5wRaNY/E8cISjPvbd1Uee/aSQfjV8E5mghFZEAs+ohoUFJWiTBXfbc31eOwjX29GWbniD2O6+SEZEQHAK4t34G/zttb6+N2fbsD+Y6cwpX9b9G8fBwA4drIY8dER/opIZCnskyKqwdnPfItBD8+Hqtbr+Mdnpdf7WCJqmpkr9tZZ7FV4YdF2/Hf5HhzKL8Kk57/HkEcX4OGvNuFkcakfUhJZi2UKPhF5S0RyRCStlsfHi0ieiKxzfjzovL+TiCwWkXQR2SQit/s3OdnR0ZMlAIDCkrJ6f84DX6ax6CPysblpWbj/840Y0im+zuO6tIpGl1bRWJ5xBJOf/x47cvIBAG8v3Y2r31rph6RE1mKZgg/A2wDO83DMD6o6xPnxqPO+UgB/UtW+AEYBuFlE+vkwJwWRFxftqHL7oQv64fu7xuPPU3qddux/l++tUvSpKgtAIi9auvMQbpu5DkM6xeO9a0dUeaz6JKrZt41BYkwkdh0qQJdWzas8tmL3EVz22lKs3H3E55mJrMIyY/hUdYmIJDfi87IAZDnbJ0QkHUAHAJu9m5CCxbKdh13txNhI/O2yQTize+sqx9wyoSdumdAThSVleOjLTQgJAVpEheP1JRkQCB69sD8en5WON3/chccvGoArR3Xx92kQ2Ura/jxc/+5qdGkVjbeuHo7YqHBsezwVi7fmYGzPRESFh0BEUFBUisiwEISFhuD2ST2x/+gpXJbSCd2rTb5aufsoLnttGd75/QiM65Vo6KyI/EesdAXCWfB9o6oDanhsPIBPAWQCOADgz6q6qYbPXwJggKoer+u1UlJSdNWqVV7JTYFPVbFi1xG8viQD327JQceEZnjtymHo374FRKTez/H0nC14fUkGrhrVBe8t3+N6bPfT03wVnci2NmQeQ2xUOADg0n8tRVR4KD65cTSS4po1+Lm+35aL3721osbH0h45FzGRlrn+QVQjEVmtqimN/fxA+glfA6CLquaLyFQAXwDoWfGgiMTAURDeUVuxJyLXA7geADp37uzzwGR9ZeWKBZsP4rXvM7Bu3zG0bB6BOyf3wlWjuiChecNm84kI7kntAwB4fUmGL+ISBY11+47hold+AgB0TGgGBfDutSMaVewBwLheiVj7wGTc/ekG3JvaBxOe+9712ICH5mHFfRPRpkWUN6ITWVLAXOGr4djdAFJU9ZCIhAP4BsA8VX2+Pq/FK3zBTVXx4cp9eH1JhnOMTzT+MKYbLhvWEVHhoU1+7oorfRV4hY+o/rLyTmH0U9+6bkdHhGLmdaMw2MNEjYZIvmfWafftempqva/oE/lbU6/wWWnSRp1EpJ0434kiMgKO7Ied970JIL2+xR7Rgs3ZuOezjYiNCsOrvzkD3/5pPK4a1aXJxR5QeaWvYmJH84hQ5JwobPLzEgWDYyeLqxR7APDGVSleLfYA4I+TTp941fXeuhdZJwpklin4RGQmgGUAeotIpohcKyIzRGSG85BLAaSJyHoALwKYro7Lk2cBuArABLclW6YaOQkKCD9sz8X1760GAPz7tymYOjAJoSHe/ateRHDLhJ748PpRKFfg4leWYlv2Ca++BpEdXfnmz1Vunz8oCWf3bF3L0Y13+6SeePN3p18s2XfkpNdfi8gKLNWl60/s0g1Oh/OLMOzxha7b259IRXiob//u2ZiZh9+/sxKFxWX415XDfPLLi8gObv7fGszamOW6ffHQDvjbpYMQ5sP36H2fb8T7P++tct/71408bWY+kWlB06VL1FQnCkuqFHsAfF7sAcDAjnH44uaz0CGhGa7+zwp8uHKv508iCjIH8wqrFHsA8PyvBvu02AOAJy8eiKS4qpM1rvj3z7UcTRS4WPBR0Bj48Pwqt28c391vr90hvhk+njEaZ/Zojb98uhHPzt2C8vLgvLpOVJNRTy2qcnvDw1P8NoFi9m1jTruvpkkdRIGMBR8FhUv+tfS0+/5yXh+/ZoiNCsebv0vBFSM749XvduLWD9Y2aOs2IruqXlz96zdnoIVz/T1/SGgeUePuOflF3HOX7IMFHwWF1XuOVrm96ykz83rCQ0PwxEUDcN/UPpi1IQtX/Hs5DucXGclCZAVp+/Oq3O4Q3wypA5P8nuOWCT1Pu2/AQ/P8noPIV1jwke1Vv3qw4r6JRtfaEhFcP7Y7/vWbM7DpwHEMe3whjhQUG8tDZNL5L/1Y5fZP90wwlATY9njqafcF68RGsh8WfGRrGbn5p91nldX0UwcmISU5AQBwx4frzIYhMuDf1XakWfzn8WaCOEWEheC2iVWv9HFtPrILFnxka+7bJwHW2/FC4LjSuGRbLscLUVDJOVGIJ2anV7mva+vmhtJUunPy6WP55qYdNJCEyLtY8JFtzam2xMMPd59jKEnt3Is8jheiYDLiiaqzck2Nq63J+oemVLk947+r8dGqfYbSEHkHCz6ypZPFpbjxf2tct0ckt0SnltEGE9XslnN6VLm9s4YuaCK7+XZLdpXbt0/saak9bOOahWP9g1WLvrs/2cCllCigseAjW/q/H3ZVuf3RjNGGktRtUr+2VW5z8gbZnari929X3eXojzV0o5oWF336sjDd7uN4PgpcLPjIdo6dLMbr3+903V5qcNZffTz1y4GudvOIMINJiHzvh+2Hqty2UldudWsfmGw6ApHXsOAj21m77xgKisvwxlXDsOHhKWgf38x0pDpNH97J1V6yPddgEiLfe2XxDle7S6toS3XlVpfQPAJv/q7q1qXzN3ECBwUmFnxkO9f8ZyUAoHubGL+u1t9YIoI1zisJn67O5Dghsq1Vu4/g511HADhm5H5neBmW+pjYt+qwi8dnpddyJJG1seAjW9l35KSr3c4i6+3VR8vmEXjx8qHYnpOPebyCQDb1yuIdaNk8ApsfPReL/zze0lf33GU8WdntvPfISbz90646jiayJhZ8ZCtfrT/gajePDKzxcNMGJqFr6+Z46dsdXN2fbOe/y/dg8dZcXD6iE6IDbKxqSEjVwvThrzcbSkLUeCz4yDZUFZ+v3Y9hXRKw9fHzTMdpsNAQwU3ju2Nz1nEs3ppjOg6R16gq/vpFGoDAuvLuLu2Rc6vcrr5lI5HVseAjW1BVfLc1Fzty8vHLMzogMizUdKRGuWhoB3SIb4YXF/EqH9nHd9sqJyOdLC4zmKTxYiLDMGNc9yr3uQ8hIbI6FnxkCx+s3Idr3nZM1pg2MMlwmsYLDw3BjeO7Y92+Y1i687DpOEReUTGRCgCuGNnZYJKmuevc3lVuj3l2saEkRA3Hgo9swX2iQ3x0hMEkTXfpsI5o2yISL3273XQUoib7ct1+V/vPU3ohNgBmztcmNETQu22s6RhEjcKCjwKeqmLbwRMAgJcuH2o4TdNFhYfi+rHdsTzjCFbuPmI6DlGT3P7BOlf7xvE9aj8wQMz749gqtxduzq7lSCJrYcFHAW/N3qM4kFeIxy8agAsGtzcdxysuH9EJrZpH4OVvd3g+mMii0vbnudpXjeqC0JDAWIalIf7w7irPBxFZAAs+Cnj/+Wk3YqPCcPHQDqajeE10RBiuHdMV32/LxYbMY6bjEDXK+S/96Go/emF/g0m8y31dPqJAwYKPAtrBvELMTTuIX6d0Crh19zy5alQXxDUL51U+CkglZeWudqeWzQJmkeX6qL4u3z2fbjCUhKj+WPBRQPvfz3tQporfjk42HcXrYqPCcc1ZyZi/ORvpWcdNxyFqkGfmbHG1Z982xmAS31jsti3cByv3mQtCVE8s+ChgFZWW4f2f92Jinzbo3CradByfuPrMZMREhlXZcJ4oEKzacxQAMLRzfEDPzK1N19bNq9w+cOyUoSRE9cOCjwJSWblixnurcbigGFef2dV0HJ+Jj47AVaO7YNbGLOzMzTcdh6je1u07BgB46pcDzQbxoV1PVY7lO/Ppb5F3ssRgGqK6seCjgPTc/K1YvNWxev9ZPVoZTuNb157dFZFhIXh18U7TUYjqZdehAle7T7sWBpP4lohUmXn88Wp27ZJ1seCjgPTqd47i57JhHW01GLwmrWMiccWILvhi3X5u5UQBYb3z6t6c2+03dq869/GJj89K55aIZFks+CjguBc9143tZjCJ/1w/thtCRfCv73mVj6zvmw1ZaB0TiR5tYkxH8bne7WLxC7f1P9/8cZfBNES1Y8FHAcd9/8ruifb/hQIA7eKicFlKR3yyKhNZeRwcTtaVfbwQi7fm4NJhHREeGhy/Yl68fCginOdaMXaRyGqC491ItnGisHJQdPu4KFuu3F+bGeO6o0wVbyzJMB2FqFYfr9qHsnLF9OGdTEfxq+TWjpUCvtmQhVK3NQiJrIIFHwWUr9dnudqvXjnMYBL/69QyGhcP7YCZK/Yi90SR6ThEpykvV3y4ah/O7N4KydWWLbG7xy+qnI18pKDYYBKimrHgo4By3+cbXe0hneLNBTHkpvHdUVRaznFCZEk/7TyEfUdOYfqIzqaj+N2Iri1d7bxTXJ6FrIcFHwWM4U8sdLVfu/IMg0nM6ZYYg/MHtcd7y3bjKK8ikMV8sGIfEqLDcW7/tqajGDX5H0uwes8R0zGIqmDBRwHDvRtzdLfWBpOYdfM53VFQXIb/LN1tOgqRy6H8IszffBC/PKMjIsNCTccx4u1rhrvaD321yWASotOx4KOAUFJtEHRctP22aqqvPu1aYEq/tnj7p104XsiuI7KGT1dnoqRMcfmI4Jqs4W5cr0RXOz3rhMEkRKdjwUcBoef9c1ztDQ9PMZjEGm6d0BPHC0vx3rI9pqMQQVXx4cp9SOmSgB5tYk3HMUZE8JuRjvGLZeWKmSv2Gk5EVIkFHwWcFjbciL2hBnaMw/jeiXjzx104WVxqOg4FuZ93HUHGoYKgnKxR3WMXDsBVo7oAAO53m2RGZBoLPrK8P7yzytWOjgjOsUE1uXVCDxwpKMb7P/MqApn1wYq9iI0Kw7SBSaajGBcSIq4Fp8sVOJzPJZTIGljwkeUtTM92tdc8MNlgEmsZ1qUlRndrhTeWZKCwpMx0HApSx04WY3baQVw8tAOa8Q8yAMDEvm1c7SdmpRtMQlTJMgWfiLwlIjkiklbL4+NFJE9E1jk/Hqzv55J9RIXzF4q7Wyf0QM6JIny8OtN0FApSn6/dj+LSckwfzu7cCmf1qFxF4LO1+w0mIapkmYIPwNsAzvNwzA+qOsT58WgDP5cCUPI9s1ztbY+nGkxiTaO7t8LQzvF4l0u0kAGqig9W7MPgjnHo176F6TiWVV6upiMQWafgU9UlABq1UmVTPpcCR0SYZX5cLUNEcEbnBBw4dsp0FApCa/Yew9bsE5ysUYNFfxrnat/z2QaDSYgcAu036GgRWS8ic0Skv+kw5Fs/bj/kavdsE2MwibXFNwtHQXEZ1u49ajoKBZkPVuxFdEQoLhjc3nQUy+meGIO2LSIBAB+tysQJrplJhgVSwbcGQBdVHQzgJQBfNPQJROR6EVklIqtyc3O9nY+87Mo3f3a1F9w5ro4jg1ups7vo4leXQpVdR+QfJwpL8M2GLPxicHvERIaZjmNJA9rHudrZxzlbl8wKmIJPVY+rar6zPRtAuIg0aH8tVX1DVVNUNSUxMdHzJxAFmG+35JiOQEHiy3UHcKqkjN25dXji4oGu9vkv/WAwCVEAFXwi0k5ExNkeAUf2w2ZTka+s3F05JPNPk3sZTGJ9vdpW7mywITPPYBIKJh+s3Is+7WIxuGOc54ODVLu4KFe7sKQcRwuKDaahYGeZgk9EZgJYBqC3iGSKyLUiMkNEZjgPuRRAmoisB/AigOnq7L+q6XNNnAN5z2WvLXO1b53Y02AS65s6sJ2r/cKi7QaTULBI25+HtP3HcfmIznD+HU61WHjnWFd76GMLDCahYGeZgReqermHx18G8HJjPpcCC7cKaxgRwRMXD8D9nzuWodyWfaLKVT8ibyovV/z7hwxEhoXgoqEdTMexvPbxzUxHIAJgoSt8RBV25OS72j/dM8FgksBxuduitze8t9pgErK7+79Iw5frDmBK/3aIa8Z9rT2Jjqh6XSX3BCdvkBks+MhyfvHyT652B/51XC8hIZXdaiVl5QaTkJ2VlStmrnDs3dyfCy3X23VjurrayzM49JzMYMFHlvW/P4w0HSGgnNffMZYvmvuZko88v2Crqz2+N1c6qK/7p/VztZfuZMFHZrDgI8sa1iXBdISA8uhFjrXIOyVEG05CdvXK4p0AgNYxkejTjlf4GuKHu8/BoI5xWLbzkOeDiXyABR9ZTnREKK49uyuiwnmlqiHaxEbhNyM7Y1nGYZwqLjMdh2wmv6hyMtVZPVoZTBKYOrWMxoVDOmD34ZPcCpGMYMFHlpKRm4+TxWVIbsWrVI0xdWASThaX4fttXICZvGvy89+72rdO6GEwSeAa3c1RKC9jty4ZwIKPLGVO2kEAwKR+bQ0nCUwju7ZEQnQ4Zm88aDoK2UxWXqGr3aMNl/1pjD7tYpEQHc5xfGQECz6ylNkbs3BG53gkxXF2bmOEhYbg3P7tsCg9G4Ul7NYl70jbX7mDS0QYf200VkiIYHT3VliecZj7XpPf8Z1LlrHncAE2HTiOqQOTTEcJaKkDk1BQXIYftnNwOHnHPxduc7VX3DfRYJLAN7pbK+w/dgp7j5w0HYWCDAs+soyKbshUFnxNcmb3VohrFo45G7NMRyGbWJheOSY0PjrCYJLAN7p7awBcnoX8jwUfWcactCwM7hTPxZabKDw0BFP6tcWC9GwUlbJbl5pmQ+YxV3v309PMBbGJ7onN0SY2Evd+thFvLNlpOg4FERZ8ZAn7jpzEhsw8TB3QznQUW5g6MAknCkvx0w5261LTuO98Q00n4hjHBwBPzt5iOA0FExZ8ZAlznbNzUwewO9cbzuzRCrFRYZytS2RBZ3bnOobkfyz4yBJmp2VhQIcW6Mz197wiMiwUk/u2xfxNB1Fcyr11qXHcZ3rPvG6UwST2MqIrCz7yPxZ8ZNyBY6ewdu8xXt3zstSBSTheWIpl3KydGumdpbtd7dG8KuU1XVs3x2TnWqMLNmcbTkPBggUfGVfZncvxe940pmdrxESGcbYuNdpTczjGzFf++eshGNQxDre8vwZr9x41HYeCAAs+Mm72xiz0aReLbokxpqPYSlR4KCb2bYN5mw6ipIzdutQwOScqd9b4929TDCaxp+aRYXjr6uGIiQzDWz/tNh2HggALPjLqYF4hVu05ysWWfSR1QBKOnizBzxlHTEehAPP5mv2u9sQ+bQwmsa/WMZHo2ro5DucXmY5CQYAFHxk1b5OjO5cFn2+M752I6IhQzE5jty41TMXeuX2TWiAkRAynsa8Dx05h6c7D2HQgz/PBRE3Ago+Mmr0xC73axqBHG3bn+kJUeCjO6dMG89IOoqyce3dS/f244xDCQwUf3sDZub50wFlYf7XugOEkZHcs+MiY3BNFWLH7CGfn+ti0gUk4XFCMFbvYrUv1szM3Hzty8nH/1L5oERVuOk5QCAvlVVTyLRZ8ZMy8TQehyu5cXxvfOxFR4SGYw25dqqeKpUIm9+fMeV979/cjAAAhwoKPfIsFHxkzJy0L3RKbo1dbduf6UnREGM7p3QZz2K1L9bRgczb6t2/Bfa39YGyvRMRGhuHoyWLTUcjmWPCREYfzi7A84wimDkiC8C9bn0sdmITcE0VYvYfrfVHdck4UYs3eo5jSj1f3/OWMLglYvCUX5fyDjHyIBR8ZMX9zNsrKFakD+UvFHyb0aYOIsBDM5iLM5MGi9ByoAlP6tzUdJWhcNLQ99h87hVX8g4x8iAUfGTF7Yxa6tIpGv6QWpqMEhZjIMIzrlYi5aQd5FYHqtGBzNjomNEOfdrGmowSNKf3aISo8BF+u2+/5YKJGYsFHfne0oBhLdx5GKrtz/WrawCQcPF6ItfuOmY5CFpVfVIofdxzClH7t+N70o+aRYZjcrx1mbcxCcSl3xSHfCKvrQRG5sx7PUaCqr3spDwWBBemO7typ7M71qwl92yAi1NGtO6xLguk4ZEFLtuWiuLQck/uxO9ffLhrSHl+vP4Al23IxiV9/8gFPV/juAhADILaOjz/5MiDZy1Nz0nH3JxvQMaEZBnaIMx0nqLSICseYnq0xZ2MWVNmtS6dbsDkb8dHhGJ7MPwj8bWyvRCREh+PL9VyAmXyjzit8AN5T1UfrOkBEmnsxD9nc699nAHD858YuI/9LHZiERVtysD4zD0M6xZuOQxZSUlaORenZmNyvHcJCOdrH38JDQzB1YBI+XZOJ/KJSxER6+vVM1DB1vqtV9W5PT1CfY4gAVNnpIXUAu3NNmNy3LcJDBXM4W5eqWbnrCI4XlrI716CLhnZAYUk55jv3GCfypnr9GSci8SJym4g8LyIvVnz4OhzZyw3vrXK1x/RMNJgkeMVFh+OsHq0xO43dulTV/M3ZiAwLwdherU1HCVrDOiegQ3wzfMl9dckH6nvdfjaAZAAbAax2+yCqt6MnS0xHIABTByRh35FT2HTguOkoZBGqivmbDmJMz0RER7Ar0ZSQEMGFQ9rjxx2HcCi/yHQcspn6FnxRqnqnqv5HVd+p+PBpMrKtu87tbTpCUJvcry1CQwSz2K1LTpsOHMeBvEJMYXeucRcO6YCycsUjX29GUWmZ6ThkI/Ut+N4TketEJElEWlZ8+DQZ2cpXbjPPbhrf3WASSmgegTO7t+JsXXKZvzkbIQJM7NvGdJSg17tdLESAr9cfwN/nbTUdh2ykvgVfMYC/AViGyu7cVXV+BpGb22audbU5O9e81AFJ2H34JNKzTpiOQhYwf9NBpHRpiVYxkaajEBzDLgDgh+2HDCchO6lvwXcngB6qmqyqXZ0f3XwZjOwj8+hJ0xGomin92yJEgDlp7NYNdvuOnMSWgyc4O9dCnrl0EABgy0H+QUbeU9+CbxMA/tamRsnKK3S11z80xWASqtA6JhKjurXCLHbrBr35m7MBgAWfhbivwXfsZLHBJGQn9S34ygCsE5HXuSwLNdQLC7cDAM7oHI+4ZuGG01CF1IFJyMgtwLbsfNNRyKD5mw6id9tYJLfmGvpWdNcnG0xHIJuob8H3BYAnACwFl2WhBlBV/LjDMQ7lhelDDachd+f2bwsRYDZn6watowXFWLn7CK/uWdDvz+oKwLHdHZE31Kvgc1+KxVfLsojIWyKSIyJptTw+XkTyRGSd8+NBt8fOE5GtIrJDRO7xZi5qmveW73G1W3NAuKW0iY3C8OSWHMcXxBZtyUG5OsZ0krUM61K5n/HqPUfqOJKofuos+ETkDU9PUJ9j6ultAOd5OOYHVR3i/HjU+fqhAF4BkAqgH4DLRaSflzJRE83ZWLlFULOIUINJqCZTB7TDtux87Mjh4PBgNH/TQbRrEYWBHeJMR6FqxvWu3I1o5e6jBpOQXXhaUv0iESms43EBcI43gqjqEhFJbsSnjgCwQ1UzAEBEPgBwIYDN3shFTbMs4zAA4IqRnQ0noZqkDkzCw19vxpyNB3HrxFjTcciPThWXYcn2XFw2rBOXSrKgmMgwRIaFoKi0HG/+uAszxnH9UmoaT126d6HqmL3qH6sA3O/LgNWMFpH1IjJHRPo77+sAYJ/bMZnO+04jIteLyCoRWZWbm+vrrEGvoKjU1X7wfF50taK2LaKQ0iUBs9O4WXuw+XHHIRSWlLM718IW3jkOAJB7ooiz6anJ6rzCZ7Ht09YA6KKq+SIyFY6JJD3huMpYXY3vDFV9A8AbAJCSksJ3j4/d//lGVzsqnN25VpU6MAmPfbMZuw4VoCtnagaN+ZsOIjYqDCO7tjIdhWrRMaGZqz3hue+x+M/jzYWhgFffWbrGqepxVc13tmcDCBeR1nBc0evkdmhHAAdqeArysy/WOb4ND/DqnqWdN6AdAM7WDSZl5YpFW3JwTu82iAgLmF8DQce9q33XoQKDScgOAuadLiLtxPnTLyIj4Mh+GMBKAD1FpKuIRACYDuArc0kJABZvzXG1+yZxbJiVdYhvhiGd4jlbN4is3nMURwqK2Z1LFEQaXPCJSIL4YISviMyEY6/e3iKSKSLXisgMEZnhPORSAGkish7AiwCmq0MpgFsAzAOQDuAjVd3k7XzUMD9nVC4jMLobu4ysburAdkjbfxx7D3NDnWAwf9NBRISGYFyvRM8Hk1F/v2ywq/3xqn11HElUN0/LsjwoIn2c7UgRWQxgJ4BsEZnkzSCqermqJqlquKp2VNU3VfU1VX3N+fjLqtpfVQer6ihVXer2ubNVtZeqdlfVJ7yZixrnk9WO/5gci/tyBqDVpTo3a+dVPvtTVSxIz8bo7q0QG8Wdb6zul0Mr5yBy1w1qCk9X+H4NYKuz/Tvnv4kAxgF40lehKLCdLC7FoXzH/o+vXHGG4TRUH51aRmNQxziO4wsC27LzsefwSXbnBoiQEP7BTN7hqeAr1sq54OcC+EBVy1Q1HZ7X8KMgdc1/VrraYaEBM0w06KUOSML6zDxkHmW3rp0t2OxYgmdyXxZ8geKOST1d7cP5RQaTUCDz9Nu4SEQGiEgiHAssz3d7LNp3sSiQ/bzLMX5vQIcWhpNQQ6Q6Z+vO5Zp8tlVWrnhx0Q4M6RSPNi2iTMeherrmzK6u9rDHFxpMQoHMU8F3B4BPAGwB8A9V3QUAznXw1vo2GgWi44UlrvbDF/Sv40iymuTWzdE3qQVeWbyjyqLZZB/pWcdRXFbOP8YCTFw0x1pS09VZ8KnqclXto6qtVPUxt/tnq+rlvo9HgebezyoXW05JbmkwCTXGoA5xOHqyBDe8t9p0FPKB81/6EQBwTu82hpNQU+Qcr2vHU6KaeZqlGyUivxORC8ThLyLyjYi84Fz0mKiKWRs46D+QtW0RCcCx7RbZF2fnBp7tT6S62nM3cdgFNZynLt13AUwBcC2A7wB0BvAygBMA3vZlMApsk/ryCkIgGtwp3tXOyM03F4S87pDbYP8urTgEO9CEu02Ae/BLLjVLDedppm0/VR0gImEAMlV1nPP+uc4FkIlc1uw96mq7LxZKgWOi28zNxVtz0S0xxmAa8qaXv93harflhI2AdG7/tpi3Kdt0DApQHpdlAQDnbhbV96ct80kiClh//HCdqx0fHWEuCHnFY99sNh2BvOjtpbsBAJ1aNjMbhBrtyYsHutpp+/MMJqFA5Kng6ygiL4rIS27titsdPHwuBZGycsUe57ZcV5+ZbDYMNcmLlw81HYG8rKi08u/zd64ZYTAJNUWrmEhX++r/rDCYhAKRpy7du9zaq6o9Vv02BbGNbn9t3jqhh8Ek1FTnuu3AcKSgGC2b82ptoHvHeXUPALvpbaJiNyOi+qqz4FPVdyraIhLjuEsLfJ6KAs6/l2S42u5/hVLgiQwLxYAOLZC2/zj+8ukG/Pu3KaYjURM9OXsLAGDawCTDSaip5tw+Bqkv/ADA0bMSyq3XqJ487nslIjeKyF4AewDsFZE9InKT76NRIJnl3IP1jM7xZoOQV3SMd8ziXLA5G4UlHK4byPJOVS6GPm0QC75A1zepctHstW4T5Yg88bQO318BXABgvHPx5VZwbLGW6nyMCJXbLQMvTOf4Lzu4YVw3V3vfEe6tG8i+Wl85324S98+1hZ/umQAAWJieYzgJBRJPV/iuAvBLVXX11znbvwLwW18Go8Dx0ap9rnanllzfyw6Gdk5wtbccPGEwCTXVA1+kudoRYR47dSgAdIhvhjO7t8KidC7RQvXn8d2vqqft4aKqpwCU+yQRBZy/fOrYTq17YnPDScib7k3tAwC4dSa3zbaDB8/vZzoCedHEvm2xPScfew/zCjzVj6eCL1NEJla/U0QmAOAeWoRTxZXju6b0b2cwCXnbVLcB/vlFpQaTUGMdL6wcv3flqC4Gk5C3VexmtJBX+aiePBV8twF4XUTeFpFbReQWEXkHwBsAbvF9PLK6nBOVF4D/PKW3wSTkba1iKpdjeXbuFoNJqLHu/9zRndsxoRm7c22mS6vm6J7YHN9u4Tg+qp86/wdQ1U0ABgBYAiAZQDdne4DzMQpyl722DAAwsEMclwewmeiIMISHOr6nBUWcqRtoCopK8bVzwsaYnq0NpyFfmNS3LX7ccQhLdx4yHYUCQL3G8KnqW6r6J1W9U1XfrD6uT0SW+S4iWVnOCceG7BXjvche7pjUCwDQOpaLLweazVnHXe0Lh3BjJDsa1zsRAHDFv382nIQCgbeu8XMn7iDXolm46QjkAzef0wOtYyJwrKDE88FkKe57rY7q1spgEvKV0W7f10P5RQaTUCDwVsGnng8hu3H/D2ZAhziDSciXuifGYFMWN2oPNI98vRkA8PhFAwwnIV8RqRxG88Xa/QaTUCDgKF5qtIo1oJqFhxpOQr40oU8bpO0/zgWYA1TqAM6eDwY/7eA4Pqqbtwo+jtYPQou35AIAlt972so9ZCMVy7NUbJ9HgSU+muMvg8HirbmmI5DFeavgu8pLz0MBoqi0DD9sz8UVIzsjLprj9+ysU8toDO4Yh1kbWPAFirJyxyibWyf04Ox5m/v8pjMr22szDSYhq/O0l+4uEclw+3C/vbPiOFVNq+t5yH6WZxxBQXGZa/FPsrdpg5KwcX8eV/UPEOnOGbrc6tD+3MdP//HD9QaTkNV5usKXAmC428cIAM/B0YW7zqfJyNK+Tc9GVHgIzuzO9b2CAbt1A8vXGw4gLEQwuW9b01HIx8JDORSf6sfTwsuHVfUwgKMAzgewGMBoANNU9RI/5CMLUlUsTM/B2T1aI4oTNoJCx4RoDO4Uj9ks+CxPVfHN+iyM6dkaCc05fi8YzL5tjKvN5VmoNp66dMNF5AYAmwGMAXChql6pqpv9ko4saWv2Cew/dgoTefUgqJw/kN26gWDtvmPYf+wULhjc3nQU8pN+7Vu42jf/b43BJGRlnq4F7wJwL4DXAMwGMFhEflnx4fN0ZEmL0h17N07sw/F7wSR1oGN5D3brWtvX6w8gIiwEk/vxD7Jg9POuI6YjkEV5KvgWwtGNOxjABdU+zvdtNLKqRenZGNQxDm1acIOVYNIxIRpDOsVj1sYDpqNQLcrKFbM2ZOGc3omIjeLs+WDyzu9HuNoFRaUGk5BVhdX1oKpe7accFCAO5Rdh7b5juH1iT9NRyIBpA5PwxOx07DlcgC6tmpuOQ9Ws3H0EOSeK2J0bhIYnJ7ja/R+ah91PTzOYhqyozoJPRO6sdpcCOATgR1Xd5bNUZFmLt+RAFZjE8XtBKXVgOzwxOx2zNmbhpvE9TMehar5efwDNwkMxgcMtgk50RJ2/zok8dunGVvtoAcdSLXNEZLqPs5HFnCouw12fbEBEWAj6uw0SpuDRMSEaQzvHcxFmCyopK8ectIOY1K8tf/kHqd+N7uJqs1uXqvO0LMsjNXzcDuBMAHf7JyJZRd8H5wIAikvLq2zaTcFl2sAkbDpwHLsPFZiOQm6W7jyMIwXFuGBQkukoZMgjFw5wtR+fxcU0qKpGrdioqkfA/XOD1nVjupqOQAalchFmS/pm/QHERoVhXO9E01HIAmau2Gc6AllMowo+EZkAx2LMFCROFld2D0zu185gEjKtQ3wzdutaTFFpGeZuOogp/dohMoyLoQezEcktTUcgi/K08PJGEdlQ7SMTwNMAbvJPRLKCfg/Oc7WT4rgcS7CbNjAJm7OOYxe7dS1hybZDOFFYigsGszs32L165Rmu9s7cfINJyGo8XeE7H6evvddbVUeo6hZfhyNr4obsVLG3Lrdas4av1x9AQnQ4zurBva2DXeuYSFf7N//+2WASshpPBV8ugAOqukdV9wBoBuB6EbnY20FE5C0RyRGRNA/HDReRMhG51O2+20UkTUQ2icgd3s4W7PI524uqaR/fDGewW9cSThWXYWF6Ns4bkITw0EaN0iGb+b/fpgAADh4vNJyErMTT/w5zASQDgIj0ALAMQDcAt4jIU17O8jaA8+o6QERCATwDYJ7bfQMAXAdgBBw7gpwvIlwV2Ite+na7q73x4SkGk5CVTBvUHpuzjiOD3UZGfbslByeLy9idSy4ju1WO41NVg0nISjwVfAmqWvHb/ncAZqrqrQBS4eWt1VR1CQBPmwDeCuBTADlu9/UFsFxVT6pqKYDvAXj9CmQwe/37DFeb2zVRhanOvXXZrWvW1+sPIDE2EiO7tjIdhSwiJrJyHcaK5bSIPBV87n8aTACwAABUtRhAua9C1UREOsBRyL1W7aE0AGNFpJWIRAOYCqCTP7MRBaOkuGYY1iUBszYeNB0laJ0oLMHirTmYNjAJoSFcKYsc3NdJLSzx669qsjBPBd8GEfm7iPwRQA8A8wFAROJ9HawG/wTwF1Utc79TVdPh6OZdAEcX9HoANQ46E5HrRWSViKzKzc31cVx7+G5r5cXUf/x6sMEkZEXTBiYhnd26xixMz0ZRaTm7c+k043pVrse49/BJg0nIKjwVfNfBsXduMoApqlrxU9MPwN99mKsmKQA+EJHdAC4F8KqIXAQAqvqmqp6hqmPh6BbeXtMTqOobqpqiqimJiVyctD6u/s9KV/vioR0NJiErSmW3rlFfr89yrIvYKcF0FLKYt68Z7mp/vy2njiMpWHjaWu2Uqj7t3E7tgIgkOu9fqqrv+SVhZZauqpqsqskAPgFwk6p+AQAi0sb5b2cAvwQw05/ZgsH71400HYEsKCmuGVK6JOAbztb1u2Mni7FkWy7OH5SEEHbnUjXu3boPfLnJYBKyCk8LL4uIPCQiuQC2ANgmIrki8qC3g4jITDhmAfcWkUwRuVZEZojIjHp8+qcishnA1wBuVlXuAuIF5eWVQzh7t401mISsbOrAJGw5eIKLvPrZ99tyUVqurq3uiKp75/cjXO0CLq8V9Dx16d4B4GwAI1S1laomABgJ4CznuD6vUdXLVTVJVcNVtaOzm/Y1Va0+SQOqerWqfuJ2e4yq9lPVwaq6yJu5glm3+2a72gnREQaTkJW5FmHmVT6/Wp5xBLGRYRjYIc50FLKokV0rl2eZ8d/VBpOQFXgq+H4L4HJV3VVxh6pmALjS+RjZVPW1m9hlRLVpFxeFlC4JmMVxfH61POMwRnRtydm5VKuo8FCc7dx95YfthwynIdM8FXzhqnraT4mq5gLggmw2tjyjcklE/kIhT6YNcnTr7shht64/HMwrxK5DBRjdnWvvUd3++4fK8deFJWV1HEl256ngK27kYxTg0vbnudqL/zTeXBAKCKkDkiDC2br+sjzjMABgVDcWfFR/z87dajoCGeSp4BssIsdr+DgBYKA/ApIZT8xOd7U7t4o2mIQCQUW3Lgs+/1iecRgtosLQN6mF6SgUQN76aZfng8i2PC3LEqqqLWr4iFVVduna1LGTlRdvf5XCtfeofqYNrOjWPWE6iu0tyziMEV1bcbgF1UuH+Gaudlk599YNVp6u8FEQ+mjVPlf78Yt4IZfqJ3Wgo1t31gZuteZLF778I/YcPolR3Vp6PpgIQP/2lVeCJzz3nbkgZBQLPqpCVfHk7C0AgP/7bQoiwvgjQvXTtkUUhndpyW5dH1uf6RhfywkbVF+vXTnM1d7DbdaCFn+bUxVHT5a42pP6tTWYhALRtEFJ2JrNbl1fWbbzsKvdpx3H71H9VF9Wq/qyWxQcWPBRFUdPcvI1NV7qgHbs1vWhy/+93NXm+D1qiIwnp7raq/dwM6pgxIKPqli92/EfQXgof5lQw7VpEYXhyS0xa+MB01FsJyvvlKt9w7huBpNQIHK/ynfpa8sMJiFTWPCRi6rirZ92oXfbWGx5LNV0HApQ0wYmYVt2PrZns1vXm9KzjrvaA9pzOzUiahgWfOTy047D2HLwBK4d05XdRdRorm5dTt7wqvX7KhdDv2Bwe4NJyA4WpWebjkB+xoKPXK5882cAwIVD+MuEGq9NiyiMSG6JWRtY8HlLebnihUXbAQCPXTTAcBoKVGmPnOtqX/vOKoNJyAQWfAQAOJRf5GpHhoUaTEJ2MG1QErbn5GMbu3W94vkF21ztCwYlGUxCgSwmMsx0BFtKvmcWnp27xXQMj1jwEQAg5fGFpiOQjZznmq3Lq3ze8Pna/a52iyhuckSNd+3ZXU1HsJV9RxzrGr763U7DSTxjwUdV9G4bazoC2UCbWGe3LsfxNVl5uWL/MccM3f7tW5y2phpRQ9w/ta+r/fr31i9SrG7Ms4tNR6g3FnyEj1ZWbqWW69a1S9QU5w9Kwg526zbZxv2VkzX+ftlgg0nIDtz/YHhqjvW7Icl7WPARHv1ms6s9g+t7kZec6+zWnZfGRZib4sJXfnK1+yZxdw1qurvO7e1qF5WWGUwS2AJt6SkWfIT8olJX+/qx3Q0mITtpExuFdi2isOcI9+4kspKbz+nhaq/cxV03GmvyP5a42oGwugULviC3eGuOq909sbnBJGRHzSPD8MnqTJwq5lWEpvp4xmjTEciGKpbjoqZ5LgCGW7DgC3LX/Gelq/3lLWcbTEJ2tCMnHwDw2dpMw0kCX592nFBF3nP9WA7faYrCkqp/xIaFWr+csn5C8huu0US+cv/naaYjBKQThSUAgD9N7oVYLsdCXnTn5F6u9tKdhwwmCUx9Hpjraj94fj+DSeqPBV8QW7X7iKs9sAP35iTvm3fHWNMRAtqcjY4JL6O7tzKchOwmKrxygf0r/s1u3aa4YmRn0xHqhQVfELv0tWWu9le3nGUwCdlVr7YxrvZcztZtsJkr96JHmxgM65JgOgrZ0D2pfVztnOOFBpMElrV7q050cS+erYwFH2FS37YQ4WKu5H3uP1cz/rvaYJLAs+XgcazdewzTh3fi+5N84ga3cXz/WLjdYJLA4n5FdOGdgdOLwYIvSGW7/TV3wWDuzUm+c8HgyuUKvt+WazBJYPlgxT5EhIbgkjM6mo5CNuX+h8TMFXsNJgksp9wmbPRoEziTqVjwBam/flE5iL4fF3MlH3px+hBX+/n5W80FCSCFJWX4bE0mzhvQDgnNI0zHIRub0q+tqz1/E4dd2BkLviB12G0LtZ7cP5d8SESQGBsJAFifmefhaAKAOWlZOF5YiukjOpmOQjb34uVDXe3r3+OwC09+83/LXe1Am5TGgi9IhYU4vvVXn5lsNggFhS9u5qSghpi5Yh+SW0VjdDfOziXfCpQJB1agqvhpx2HX7d4BtjYmC74g1SGhGQDgvql9DSehYNA+Lsp0hICxIycfK3YdwfQRnTlZg/ziXrfZunmnSgwmsbaTAb5jEAu+ILX7cAEGd4pHRBh/BMj3RMS1YXv1Feqpqg9X7kVYiHCyBvnNDeMq91Af/Mh8g0ms7ev1B1zt168aZjBJ4/C3fRDKyM3H2r3HkDqgnekoFEQqJgct3pLj4cjgVVRahk/X7Mfkfm1d4x6J/K2gqNR0BMtRVdzz2UbX7XP7B97vTxZ8Qejj1ZkIDRH8cmgH01EoiIztlYj2cVF4n8s/1Gr+pmwcKSjG5SMCY+V+so/z3AoY9z3WyeGiV5e62i9fMbSOI62LBV+QKS0rx6erM3FO70S0acFxVeQ/oSGC6SM644fth7D7UIHpOJb0wcq96JjQDGf3aG06CgWZf115hqu9wm3bTXJYv++Yq90rQFe2YMEXZJZsz0XOiSJclsLlHsj/fj28E0JDhIu81mDP4QL8tOMwfp3SCSEhnKxB/iUiGNsr0XX7kNvSXVQVCz4KCB+vykTrmAhM6NPGdBQKQm1bRGFy37b4eHUmiko5ecPdByv3IUTAP8bImNevrJyIkPL4QoNJrOVGm2wLyYIviBzOL8LC9GxcPLQDwkP5rSczfjOqM44UFGNuGlf1r1BSVo6PV2ViQp+2aMclbMiQZhFck68mc9z+r5r/x8BabNkdf+sHkS/WHUBJmfIKAhl1VvfW6NwyGrd/sA4/bOfeugCwKD0Hh/KLcDl31iDDKpZPAoCFm7MNJrGGp2anV7kdqN25AAu+oKGq+HjVPgzuFB/QP7AU+EJCBJcNc6wxd9WbKwynsYaZK/aiXYsojHMbQ0Vkws3n9HC1//DuKoNJrOH1JRmmI3gNC74gsHrPUVz86lJsOXgCv0rhYq5k3nS3ZUfyTgb3yv6ZR09iyfZc/Gp4J4RxqAVZDHfeqPTJjNGmIzSJZf53EZG3RCRHRNI8HDdcRMpE5FK3+/4oIptEJE1EZooIB8G4+dXry7DOOaX8gsHtzYYhAqosKvzfn/cYTGLeR6syAYB/jJFl/Pfaka52MO+8kXzPrCq3U5JbGkriHZYp+AC8DeC8ug4QkVAAzwCY53ZfBwC3AUhR1QEAQgFM913MwFJQVIqycgUADO4YhxZR4YYTETm0bB4BAPjbvK2Gk5hTWlaOj1ftw9ieieiYEG06DhEA4OyeXAeyuqd/OdB0hCazTMGnqksAeFrt8VYAnwKovjdTGIBmIhIGIBrAgeqfGKyy8k652n+a0ruOI4n8a9ZtZ7va6VnHDSYx5/ttucjKK+RkDbKcHm1iXO1F6cE3eePYyeIqt39lg8mOlin4PHFeybsYwGvu96vqfgB/B7AXQBaAPFWt8Rq0iFwvIqtEZFVubnDMDpz0/BJXm6v3k5W0c9vpJVi3cpq5Yh9ax0RiYt+2pqMQVbHwznGu9rXvBN/kjSGPLqhy2w6LoQdMwQfgnwD+oqpVVmsVkQQAFwLoCqA9gOYicmVNT6Cqb6hqiqqmJCYG32w4O/zAkn2ICP5yXh8AwMHjhUE3eeNgXiEWb83BZSkduS4mWV5xabnpCMbMvm2M6QheEUj/y6QA+EBEdgO4FMCrInIRgEkAdqlqrqqWAPgMwJnGUlrUggBeLJLs6/qx3VztF7/dbjCJ/328ah/KyhXThwd+VxHZ0ytXVO6ve+dH68wF8bPf/N/yKrf7tW9hKIl3BUzBp6pdVTVZVZMBfALgJlX9Ao6u3FEiEi0iAmAigPTanyk4qGqV7WDcZ0USWUWo21XnN3/cZTCJf5WXKz5ctQ9n9WiFLq2am45DVKNpg5Jc7W82ZBlM4l8/7Tjsar/5uxSDSbzLMgWfiMwEsAxAbxHJFJFrRWSGiMyo6/NU9Wc4CsA1ADbCcU5v+DywxeWeKHJtB/PwBf0QHx1hOBFRzebcXtldsu/ISYNJ/OfHHYeQefQUpg/v7PlgIoNuGFd5Fb5ixQc7yz5eWOX2mJ72Gf5lmYJPVS9X1SRVDVfVjqr6pqq+pqqv1XDs1ar6idvth1S1j6oOUNWrVLXIv+mt50RRqat99VldDSYhqlvfpMrukjHPLoaq/X+pzFyxFwnR4ZjSn5M1yNrucY6zBYDu9802mMQ/1uw56mqf178dIsIsUyY1mX3OhFy2Z5/AxOe+Nx2DqN7cx7HZfWX/3BNFWLA5G5ec0RGRYdysnqzNMVKq0qF8e19PufF/a1zt164aZjCJ97Hgs6ELX/nJ1f7zlF4GkxDVz8O/6O9qr3XuCmNXn67JRGm5VtlejsjK7p/a19W+9f21BpP41vPz7b0IPAs+GzpZXLlyzXVusyCJrCoqPBQ3ju8OwLEmX2FJmYfPCEyqig9W7MWI5JZVFrYlsjL33yPLMg7XcWRg+3TNflf7g+tHGUziGyz4bKakrHKtpEl927LLiALGb0d3cbXfXbbbXBAfWpZxGLsPn8R07qxBAWZcr8rJC1+s3V/HkYEpIzcf+49V7kw1qlsrg2l8gwWfjRSXlmPGe5VLsbzym6EG0xA1TFJcM8y5fQzG9GyNl7/dcdrWRnbwwYp9aBEVhqkDkzwfTGQhb18z3NW+48N15oL4gKpigtu49z7tYg2m8R0WfDayeGsOFm2p3GaYV/co0PRNaoG/TuuH/KJSvPTtDtNxvOpIQTHmph3EL8/oiKhwvjcpsFSfvLFws3321122s2o3tftSUXbCgs+mVtw/0XQEokbp3S4Wlw3rhHeX7cbew/ZZl++zNZkoLitndy4FrH/9pnLnjT+8a5/9dYvdhkLdObnXacWtXbDgs4msvFO4bWbl7Kk2sVF1HE1kbXdO6YWwkBA8O2+L6Sheoar4YOU+DOkUjz7t7LFNEwWf1GpDEdzHvAWqYyeL8c7S3a7b142x70RHFnw2cd27q1Dk3Nz6rnN7G05D1DRtW0ThurHd8M2GLKzde9TzJ1jc6j1HsSMnH1dwKRYKcLec08PVPuvpbw0m8Y6HvtqExVtzAQDhoYJmEfYdbsGCzwZUFWn7j7tu3+Rc3oIokN0wthtax0TiydnpAb/7xvsr9iImMgznD+ZkDQpsM2z2+2Wv23aOGx8+12AS32PBF+DKyrXKpvPnD0qy7fgDCi7NI8Pwx8k9sXL3UcwP4AHieadKMHtjFn4xpD2iI8JMxyFqkpjIMLSPqxwylHzPLINpmiZtfx7W7j0GAOjZJsb2k6lY8AW4l7/dgcdnpQNwdOU+e+kgw4mIvOfXKZ3Qo00Mnpmzpcoak4Hky3X7UVhSzu5cso2l91adFBioV+DPf+lHV/u2iT0NJvEPFnwBbtGWyisfN5/Tg1cQyFbCQkNwz3l9kHGoAB+s2Gs6ToOpKt7/eS8GdGiBAR3iTMch8oneD8w1HaHBMo9WduVOHdgOFwxubzCNf7DgC2CbDuRhQ2ae6RhEPjWxbxuM7NoS/1y4HScKS0zHaZD1mXnYcvAEpg/n1T2yl8cuGuBqF5cG3tX3s59Z7GpfP9Ze4xJrw4IvgE17sfJydKeWzQwmIfIdEcH90/ricEExXv8+w3ScBvlgxV40Cw/FhUPsf/WAgstVo7pUuR1IY/kOVFtOZlCQXH1nwRegysurjpn45lZ7rgxOBACDOsbjF4Pb4/9+zMDBvELTceolv6gUX60/gAsGJyE2Ktx0HCKvW/PA5Cq3A2Us3xfrKvcC/uiG0QgJCY6Jjiz4AlBW3il0u2+26/bDF/RDXDP+QiF7u+vc3igvB56bv9V0FI9eWrQdAx6ah5PFZZjOyRpkUy2bR1S53fXe2bUcaR2qis/WOAq+1jGRGNG1peFE/sOCLwBd9tqyKrerr35OZEedWkbjd2d2wSdrMpGeddzzJxhSWlaO5xZsAwD0bhuLoZ3izQYi8qH3/zCyyu28k9YeZ7txfx525OTjqV8OxKq/TjIdx69Y8AWYQ/lFyDxaOf7gq1vOQtsW3EaNgsMt5/REi6hwPDXHuluuLdme62pfMqwD18UkWzuzR+sqtwc/Ot9Qkvr5bM1+RISFYGoQXihhwRdAyssVKY8vrHLfoI7xZsIQGRAXHY5bJ/TAkm25WLIt1/MnGPDxqkwAQGxkGC4d1slwGiLfe+aSgVVu3/jf1YaS1K6gqBRPzU7Hf5fvweR+bYNyGBQLvgDy9YYDVW7fk9rHUBIic64a3QWdWjbDk7PTUVZunUHiac6uooXp2bj27K7Y+Mi5p41xIrKjX1dbdmhO2sHTJhaadvenG/D6kgyUlisuOaOD6ThGsOALEMWl5Xhu/rYq990wtpuhNETmRIaF4q5z+2DLwRP4bE2m6TgAgNV7juL8l37EpOe/R0mZ4tJhHU1HIvKr6mP53CcWmlZYUoZZG7Jct8f0TDSYxhwWfAHiw1X7qmzy/P1d4zk2iILWBYOSMLhjHJ6bvw2nistMx8El/1rqag/o0AJ9k1oYTEPkf9XH8gHA2r1HDSSpqrCkDH3cdgI5q0crhIcGZ+kTnGcdYI4WFOOBL9IwuGMcfrj7HNyb2gedW0abjkVkjIjgvql9cfB4Id76aZfRLNUXnL2M4/YoSO16amqV2xe/utT4H2THqs0a/u3oZDNBLIAFXwAY+tgCAEC/9i3QqWU0bhjXnVf3KOiN7NYKk/u1xb++24lD+UVGMtS00OwvgmBPTqKa1PR7qe+D5vbZ7X7fbIx6apHr9vDkBJzbv52xPKax4LM496sHXVs3N5iEyHruSe2DUyVleHHRdr+/9pGC4tMWmn30wv5I4EQNCmK7n5522n2mtl1zn9TVPbE53r9ulJEcVsGCz8IueOnHKrdDQ/jtInLXPTEGl4/ohPd/3ouM3Hy/va6q4g/vrDzt/mDuLiKq8ML0Iafd58/F0nfk5J9WZL577cigHbtXIbjP3sJUFRv351W574zO8WbCEFnY7RN7ITIsBM/M9d9izDe/vwZr9h6rct/3d4332+sTWdmFQzrg0xvPrHJf6gs/IL+o1C+vP+n570+7r0N8M7+8tpWx4LOgwpKyGvckHNo5wUAaImtLjI3EjHHdMW9TNlbuPuKX15y98WCV2+cPSkKXVhxyQVRhWJfTf18NeGiez1+3pu7jhy/o5/PXDQQs+CymvFzx8rc7Tru/+uwnIqr0hzHd0LZFJJ6cnV7jRApvUVX83w8Zp93/8hVn+Ow1iQLVtsdTT7vPV+P5VLXW5776rK4+ec1Aw4LPYm79YC1eXly14BuenMBZuUR1aBYRij9N7o21e4+ddvXNW4pLy3H/F2l4fFZ6lftrGqROREBEWAi+ufXs0+73dtGnqjX2igE1F53BigWfReSdKsGfPlpfZTXwCh/POLOGzyAid5cM64g+7WLxzNwtKC4t9+pzH84vwpVv/oz3f95b5f65d4zx6usQ2c2ADnFYeOe40+73ZtFXW7H32U1nIiKMZU4FfiUsYPHWHEz5x/f4Yt3+0x5Le+RcA4mIAk9oiOCe1D7Ye+Qk3lu+x2vPm551HL94+Ses33cML0wfgud/NRgA8PGM0ejTjjtqEHnSo00MQkNO76VqatG3YteROp/jDI57r4IFn0HHC0tw9yfrcc1/VqJFVDg+v6nqlbyMJ6ciJjLMUDqiwDOuVyLO7tEaj32zGU/P2eIaz5dzvLBRzzc3LQuX/GspysoVH88YjQuHdMDFQzvgh7vPwfDklt6MTmRrO5+seRx68j2zkHeqpMbH6nLxqz/hV68vq/Gxq0Z14VCLGogvBzhbWUpKiq5atcrY6y/Zlou/fLoB2ccLMWNcd9w+qSciw0Lxh3dWYWF6Nn78yznomMDt04gaatOBPEx70bGG5WMXDcCP23Mxb1M2AOD960bizO6n7/lZXWlZOXrcPwcAMKRTPN64ahjatIjyXWiiIFBeruh2X83dr2N7JeLd34/w+By7DxVg/N+/q/sYmxZ7IrJaVVMa/fks+HznRGEJSsoUhSVlOFlcilbNIxEWKnhydjpmrtiH7onN8dyvhmBIp3if5iAKNv/6bmed6/LdN7UPfp3SGXHR4VXu/+fCbfjnwqq7dmx57DxEhYf6JCdRMPLUlfvpjaPRIzEWcdHhOFVchl2HCnDT/1ZjTM/EOodrvHV1Cib0aevtuJbBgq+R/FHw9X9wLgqqbRzdIb4ZsvJO4bqx3fDHSb34i4TIR46dLMaQRxc0+XnserWAyCRvz9Rd/9AUxDUL93xgAGtqwccxfD5UvdgDHOP2Pp5xJu5N7ctij8iH4qMj8OmNo03HIKIabH/Ce8ulrHlgsu2LPW9gwednFw5pX+MK5ETkfcO6tOQVOiILCg8NafJ7c9Mj52L309PQsnmEl1LZm2UKPhF5S0RyRCTNw3HDRaRMRC513u4tIuvcPo6LyB1+Cd0ICdH8wSTyt91PT0PaI+diTE/PEzYqtImNxMaHp/gwFRHtfnpavd9nK+6fiEl92+DV35yBHU+kojlXsWgQy4zhE5GxAPIBvKuqA2o5JhTAAgCFAN5S1U9qeHw/gJGqWudCXP4Yw5d3qgQfrdwHEeDYyRK0bB6BK0Z2ZlcukUEHjp0CAHyyOhPPL9hW5bEQAcoVWP3XSWgVE2kiHlHQuvezjZi5Ym+Nj/3fb1MwqZ99J2TUh60mbYhIMoBv6ij47gBQAmC487jqBd8UAA+p6lmeXsv0sixEZA2qiiMFxSzwiCxGVVFUWs6LJE5BM2lDRDoAuBjAa3UcNh3ATP8kIiI7EBEWe0QWJCIs9rwoYAo+AP8E8BdVPX3qKwARiQDwCwAf1/YEInK9iKwSkVW5ubm+SUlERERkMYE04jEFwAciAgCtAUwVkVJV/cL5eCqANaqaXdsTqOobAN4AHF26vo1LREREZA0BU/CpateKtoi8DccYvi/cDrkc7M4lIiIiOo1lCj4RmQlgPIDWIpIJ4CEA4QCgqnWN24OIRAOYDOAGH8ckIiIiCjiWKfhU9fIGHHt1tdsnAbTydiYiIiIiOwikSRtERERE1Ags+IiIiIhsjgUfERERkc2x4CMiIiKyORZ8RERERDZnqb10/UlEcgHs8cNLtQZwyA+vY1IwnCPA87Qbnqe98DztJxjOtSHn2EVVExv7QkFb8PmLiKxqymbHgSAYzhHgedoNz9NeeJ72Ewzn6s9zZJcuERERkc2x4CMiIiKyORZ8vveG6QB+EAznCPA87YbnaS88T/sJhnP12zlyDB8RERGRzfEKHxEREZHNseCrRkQ6ichiEUkXkU0icrvz/pYiskBEtjv/TXDe38p5fL6IvFztuX4tIhucz/NsHa85TEQ2isgOEXlRRMR5/1gRWSMipSJyqR3P0e3xS0VERcQrs5WsdI4i8g8RWef82CYix7xxjobP9QkR2Sci+dXujxSRD51fg59FJNngeU4WkdXO78lqEZng9lx1/jx6Os5X702rnafb4159f1rtPH35HjV0noH+/qwxfw2vafXfnT47R7fH6//eVFV+uH0ASAJwhrMdC2AbgH4AngVwj/P+ewA842w3B3A2gBkAXnZ7nlYA9gJIdN5+B8DEWl5zBYDRAATAHACpzvuTAQwC8C6AS+14jm4ZlgBYDiDFjufodsytAN6ywc/sKOfr5le7/yYArznb0wF8aPA8hwJo72wPALC/Id8rE+9Nq52nr96fVjxPt2O8+h41dJ6B/v6sMX99v5+wzu9On52jW4Z6vze98kWw8weALwFMBrAVQJLbN31rteOuRtVfnsMBLHS7fRWAV2v5AdridvtyAK9XO+Ztb/7QWu0cAfwTwPkAvqvPD20gnqPb/UsBTA7kn9lqz1H9F8o8AKOd7TA4FhQVk+fpvF8AHAYQ2YDvlfH3phXO0x/vTyucp9v9Pn2P+vo8q31+wL0/68rf0O+nr9+fps+xoe9NdunWwXm5eyiAnwG0VdUsAHD+28bDp+8A0EdEkkUkDMBFADrVcFwHAJlutzOd9/mF6XMUkaEAOqnqN004jTqZPke3HF0AdAXwbcPPon78dK516QBgn/M1SwHkwXHl0KsacZ6XAFirqkWo/3vO6HsTMH+e/nh/Ol8nGRb4fvr6Peqn86xLILw/6yvQfnd69Rwb894Ma8ALBxURiQHwKYA7VPV4LUMlaqWqR0XkRgAfAiiH46/GbjW9VE2f3sC4jWL6HEUkBMA/4LjS5BOmz7Ha7ekAPlHVsgaFqCc/nmudMWp66gY+R90v0MDzFJH+AJ4BMKWBGY29NwHz5+mP9ydg/jyr3fbZe9SP51nn03rhOep+gaafZ71fqob7LPm709vn2Nj3Jq/w1UBEwuH4Zv5PVT9z3p0tIknOx5MA5Hh6HlX9WlVHqupoOC75bheRUKkcHPwoHBV7R7dP6wjggDfPpyYWOcdYOMY1fCciu+EY0/BVvQafBs45upsOYGbTzqpmfj7XumTCeVXQeZUwDsCRxp3V6Rp6niLSEcDnAH6rqjvdMp72vbLKe9OZ2wrn6dP3p4XO051P3qN+Ps+6BML7s7bntsT70yLn2Lj3pi/6tQP5A46K+l0A/6x2/99QdVDms9Uevxpu46Gc97Vx/psAYB2AXrW85krnN6xiUObUao+/De8OPLXcOTqP+Q7em7RhqXME0BvAbvhgvIyJc3U7vvoYoZtRdVD4R6bOE0A8gPUALmnMz2N9jvP2e9Oq5+k8xmvvTyuep6/eoybO0+34gHx/1pa/oV8Pb78/rXiOzmPq9d702g+1XT7gmL2oADbA8QtvHYCpcIxzWARgu/Pflm6fsxuOv5Ly4ajI+znvnwlgs/Njeh2vmQIgDcBOAC9X/IcDxyD6TAAFcAz23GS3c2zMD20gniOAhwE8baOf2Wedn1fu/Pdh5/1RAD6GYzzgCgDdTJ0ngL/C8d5Z5/ZRUdB6/Hk08d602nn66v1pxfOEj96jhs4z0N+fNeY3/f600jk25r3JnTaIiIiIbI5j+IiIiIhsjgUfERERkc2x4CMiIiKyORZ8RERERDbHgo+IiIjI5ljwEREREdkcCz4iIiIim2PBR0RERGRzLPiIiIiIbI4FHxEREZHNseAjIiIisjkWfEREREQ2x4KPiIiIyOZY8BERERHZHAs+IiIiIptjwUdERERkcyz4iIjqQUTKRGSdiGwSkfUicqeI1Pl/qIgki8gV/spIRFQbFnxERPVzSlWHqGp/AJMBTAXwkIfPSQbAgo+IjBNVNZ2BiMjyRCRfVWPcbncDsBJAawBdALwHoLnz4VtUdamILAfQF8AuAO8AeBHA0wDGA4gE8Iqqvu63kyCioMWCj4ioHqoXfM77jgLoA+AEgHJVLRSRngBmqmqKiIwH8GdVPd95/PUA2qjq4yISCeAnAJep6i5/ngsRBZ8w0wGIiAKYOP8NB/CyiAwBUAagVy3HTwEwSEQudd6OA9ATjiuAREQ+w4KPiKgRnF26ZQBy4BjLlw1gMBxjowtr+zQAt6rqPL+EJCJy4qQNIqIGEpFEAK8BeFkd42LiAGSpajmAqwCEOg89ASDW7VPnAbhRRMKdz9NLRJqDiMjHeIWPiKh+monIOji6b0vhmKTxvPOxVwF8KiKXAVgMoMB5/wYApSKyHsDbAF6AY+buGhERALkALvJPfCIKZpy0QURERGRz7NIlIiIisjkWfEREREQ2x4KPiIiIyOZY8BERERHZHAs+IiIiIptjwUdERERkcyz4iIiIiGyOBR8RERGRzf0/bQPJ3SUac5oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,5))\n",
    "ax.plot(t_rec_hmi,keyshmi['DSUN_OBS'])\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('DSUN_OBS [m]')\n",
    "ax.xaxis.labelpad=15\n",
    "ax.yaxis.labelpad=15\n",
    "ax.set_title(\"Distance from SDO to Sun\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "edda8d39-791f-473f-a729-2636fac48261",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The keyword RSUN_OBS varies by 3.342880394955568 percent.\n"
     ]
    }
   ],
   "source": [
    "print(\"The keyword RSUN_OBS varies by\",(np.nanmax(keyshmi['DSUN_OBS'])-np.nanmin(keyshmi['DSUN_OBS']))/np.nanmax(keyshmi['DSUN_OBS'])*100,\"percent.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91927260-e19e-409d-aaba-a5063ea24715",
   "metadata": {},
   "source": [
    "RSUN_OBS and DSUN_OBS vary inversely; that is, when DSUN_OBS is at its maximum value, RSUN_OBS is at its minimum value. So we will evaluate how much value of MEANGB\\* varies between two cases:\n",
    "\n",
    "1. When DSUN_OBS is at maximum and RSUN_OBS is at minimum, and \n",
    "2. When DSUN_OBS is at minimum and RSUN_OBS is at maximum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "8c2f4152-f9c3-4087-8629-b612102a726b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab data for when DSUN_OBS is at maximum and RSUN_OBS is at minimum\n",
    "dsun_obs_max = np.nanmax(keyshmi['DSUN_OBS'])\n",
    "rsun_obs_min = np.nanmin(keyshmi['RSUN_OBS'])\n",
    "cdelt1 = 0.03\n",
    "rsun_ref = 696000000.\n",
    "\n",
    "# Grab data for when DSUN_OBS is at minimum and RSUN_OBS is at maximum\n",
    "dsun_obs_min = np.nanmin(keyshmi['DSUN_OBS'])\n",
    "rsun_obs_max = np.nanmax(keyshmi['RSUN_OBS'])\n",
    "cdelt1 = 0.03\n",
    "rsun_ref = 696000000."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9158d340-ce23-48e0-ba51-375ce3d85385",
   "metadata": {},
   "source": [
    "Find the value of the factor for Case 1: when DSUN_OBS is at maximum and RSUN_OBS is at minimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "e2316312-1c20-43a0-9131-5ea579ad4cdf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.7440603945984874\n"
     ]
    }
   ],
   "source": [
    "cdelt1_arcsec_case1 = (math.atan((rsun_ref*cdelt1*radsindeg)/(dsun_obs_max)))*(1/radsindeg)*(3600.)\n",
    "factor_case1 = (1/cdelt1_arcsec_case1)*(rsun_obs_min/rsun_ref)*(1000000)\n",
    "print(factor_case1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb652a8b-e94b-45d5-9c1e-f9873f8ec86f",
   "metadata": {},
   "source": [
    "Find the value of the factor for Case 2: when DSUN_OBS is at minimum and RSUN_OBS is at maximum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ab5d22cf-3e04-489e-8fc2-7cd7cb9d0fb7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.744061135022291\n"
     ]
    }
   ],
   "source": [
    "cdelt1_arcsec_case2 = (math.atan((rsun_ref*cdelt1*radsindeg)/(dsun_obs_min)))*(1/radsindeg)*(3600.)\n",
    "factor_case2 = (1/cdelt1_arcsec_case2)*(rsun_obs_max/rsun_ref)*(1000000)\n",
    "print(factor_case2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}