{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ac92d7c1-e5af-4f1a-b3df-1fd361073edc",
   "metadata": {},
   "source": [
    "# Variations of the IP, the PG and the RMM with the MJO phase\n",
    "Here we plot variations of various parameters with the MJO. In particular, we plot the ionospheric potential (IP), the fair-weather potential gradient (PG) and the two components of the RMM index. We also plot the sunspot number and outgoing long-wave radiation (OLR)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce2dfdec-3074-422d-b615-4f0c770c671f",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "26804089-74df-41d6-8c8e-b60af3173175",
   "metadata": {},
   "outputs": [],
   "source": [
    "# data processing\n",
    "import numpy as np\n",
    "import csv\n",
    "\n",
    "# plotting the data\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.path import Path\n",
    "from matplotlib.patches import PathPatch\n",
    "\n",
    "# dates\n",
    "import datetime as dt\n",
    "\n",
    "# testing correlation significance\n",
    "import scipy.stats as st"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7e68cb6-0540-4cca-b1b9-32e52321b9c4",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Loading the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "53dda0e6-8893-4ae2-82f3-c57c8550d2d8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "map_contrib = np.load('./data/DAILY-IP-MAP-V4.3.npy')\n",
    "# original data with the shape (number of days, number of latitudes, number of longitudes)\n",
    "# contains IP values (not normalised) depending on (d, lat, lon)\n",
    "# d (axis 0) is the number of a day starting with 0 and ending with 4991\n",
    "# every third day is taken, 0 corresponds to 1 Jan 1980 and 4991 corresponds to 29 Dec 2020\n",
    "# lat (axis 1) describes the latitude (an integer in [0, 180])\n",
    "# lon (axis 2) describes the longitude (an integer in [0, 360])\n",
    "\n",
    "map_contrib /= np.mean(np.sum(map_contrib, axis=(1, 2)))\n",
    "map_contrib *= 240e3\n",
    "# normalisation of contributions to the IP to the global mean of 240 kV\n",
    "\n",
    "ip = np.sum(map_contrib, axis=(1, 2))  # total IP values for different days"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d663a94f-942b-4583-9179-08699c48bb4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "map_olr_source = np.load('./data/OLR_41years_NOAA.npy')\n",
    "# original data with the shape (number of days, number of latitudes, number of longitudes)\n",
    "# contains OLR values depending on (d, lat, lon)\n",
    "# d (axis 0) is the number of a day starting with 0 and ending with 14976\n",
    "# every third day is taken, 0 corresponds to 1 Jan 1980 and 14976 corresponds to 31 Dec 2020\n",
    "# lat (axis 1) describes the latitude (an integer in [0, 180])\n",
    "# lon (axis 2) describes the longitude (an integer in [0, 360])\n",
    "# the data downloaded from https://psl.noaa.gov/mddb2/showDataset.html?datasetID=37\n",
    "\n",
    "olr_days_ind = np.flatnonzero(np.amin(map_olr_source, axis=(1, 2)) >= 0)\n",
    "# indices of the days for which the OLR data are available\n",
    "map_olr = np.take(map_olr_source, olr_days_ind, axis=0)\n",
    "\n",
    "olr = np.sum(\n",
    "    np.mean(\n",
    "        map_olr[:, 75:105, :], axis=2\n",
    "    ) * (\n",
    "        np.cos(np.arange(75, 105) * np.pi / 180) -\n",
    "        np.cos(np.arange(76, 106) * np.pi / 180)\n",
    "    ), axis=1\n",
    ") / np.sum(\n",
    "    np.cos(np.arange(75, 105) * np.pi / 180) -\n",
    "    np.cos(np.arange(76, 106) * np.pi / 180)\n",
    ")\n",
    "# global means of OLR values for 15°S–15°N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9fd8402a-3763-47d5-b510-413eb03ee3aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "d_start = dt.date(1980, 1, 1)  # first date\n",
    "d_end = dt.date(2021, 1, 1)  # last date\n",
    "\n",
    "pg = np.full(((d_end - d_start).days, 24), np.nan)\n",
    "\n",
    "# hourly data obtained from 10-s data\n",
    "with open('./data/'\n",
    "          'vostok_hourly_from_10_s_without_calibration_and_empty.tsv') as file:\n",
    "    tsv_file = csv.reader(file, delimiter='\\t')\n",
    "    next(tsv_file)  # skip the header\n",
    "    for line in tsv_file:\n",
    "        timestamp = dt.datetime.strptime(line[0], '%Y-%m-%d %H:%M:%S')\n",
    "        pg_value = float(line[1]) / 3  # hourly PG value in V/m (3 is the form factor)\n",
    "        if 0 < pg_value <= 300:\n",
    "            day = (timestamp.date() - d_start).days\n",
    "            hour = timestamp.hour\n",
    "            pg[day, hour] = pg_value\n",
    "            # negative PG values and PG values greater than 300 V/m are omitted\n",
    "\n",
    "# hourly data obtained from 5-min data\n",
    "with open('./data/'\n",
    "          'vostok_hourly_from_5_min_without_calibration_and_empty.tsv') as file:\n",
    "    tsv_file = csv.reader(file, delimiter='\\t')\n",
    "    next(tsv_file)  # skip the header\n",
    "    for line in tsv_file:\n",
    "        timestamp = dt.datetime.strptime(line[0], '%Y-%m-%d %H:%M:%S')\n",
    "        pg_value = float(line[1]) / 3  # hourly PG value in V/m (3 is the form factor)\n",
    "        day = (timestamp.date() - d_start).days\n",
    "        hour = timestamp.hour\n",
    "        if np.isnan(pg[day, hour]) and 0 < pg_value <= 300:\n",
    "            pg[day, hour] = pg_value\n",
    "            # negative PG values and PG values greater than 300 V/m are omitted\n",
    "            # 5-min data are used only if there are no 10-s data available\n",
    "\n",
    "diff = (np.amax(pg, axis=1) - np.amin(pg, axis=1)) / \\\n",
    "       np.maximum(np.mean(pg, axis=1), 1)\n",
    "# the elementwise maximum in the denominator is taken in order to avoid division by zero\n",
    "\n",
    "fw_days_ind = np.flatnonzero(diff <= 1.5)\n",
    "# indices of fair-weather days according to the criteria used\n",
    "pg_fw = np.take(np.mean(pg, axis=1), fw_days_ind)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "dabe9bd5-2e61-48ff-b63f-cd2d0010326c",
   "metadata": {},
   "outputs": [],
   "source": [
    "sn = np.full(((d_end - d_start).days), np.nan)  # sunspot number\n",
    "\n",
    "with open('./data/sunspot_number_data.csv') as file:\n",
    "    # the file downloaded from https://www.sidc.be/silso/datafiles\n",
    "    csv_file = csv.reader(file, delimiter=';')\n",
    "    for line in csv_file:\n",
    "        timestamp = dt.date(*list(map(int, line[:3])))\n",
    "        day = (timestamp - d_start).days\n",
    "        if day >= 0 and line[4] != -1 and timestamp.year < 2021:\n",
    "            sn[day] = int(line[4])\n",
    "            # only the available data for 1980–2020 are used\n",
    "\n",
    "sn_days_ind = np.flatnonzero(np.logical_not(np.isnan(sn)))\n",
    "# indices of the days for which the sunspot number is available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c45b721f-6f09-4c8e-85ed-528bc15e02d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert sn.shape == sn_days_ind.shape\n",
    "# all days are available for the sunspot number in 1980–2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c61cb9c0-5ead-4abf-a725-726b35ea055f",
   "metadata": {},
   "outputs": [],
   "source": [
    "sn_fw = np.take(sn, fw_days_ind)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa57840a-f79c-4c0a-9727-ee3233bd0594",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Averaging over MJO phases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e9d4adde-6492-40ee-b11b-5c47cc0e9dd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "rmm_source = np.genfromtxt('./data/rmm.txt')\n",
    "\n",
    "def phase_separation(rmm, val):\n",
    "    \"\"\"\n",
    "    Compute average values of a variable for each MJO phase.\n",
    "    :param rmm: an array of RMM values, the shape: (number of days, 2)\n",
    "    :param val: an array of values, the shape: (number of days)\n",
    "    :return: average values for MJO phases, the shape: (8)\n",
    "    :return: average square values for MJO phases, the shape: (8)\n",
    "    :return: number of days for each phase, the shape: (8)\n",
    "    \"\"\"\n",
    "\n",
    "    assert rmm.shape[0] == len(val)\n",
    "\n",
    "    angle_rmm = np.arctan2(rmm[:, 1], rmm[:, 0])  # phase angles of RMM values\n",
    "    phase_rmm = np.floor((angle_rmm / np.pi + 1) * 4).astype(int)  # phase numbers\n",
    "\n",
    "    phase_avg_val = np.zeros((8), dtype=float)\n",
    "    phase_avg_sqr = np.zeros((8), dtype=float)\n",
    "    counter = np.zeros((8), dtype=int)\n",
    "\n",
    "    for i in range(len(phase_rmm)):  # summing over each phase\n",
    "        counter[phase_rmm[i]] += 1\n",
    "        phase_avg_val[phase_rmm[i]] += val[i]\n",
    "        phase_avg_sqr[phase_rmm[i]] += val[i]**2\n",
    "\n",
    "    phase_avg_val /= counter\n",
    "    phase_avg_sqr /= counter\n",
    "    # averaging over each phase\n",
    "\n",
    "    return phase_avg_val, phase_avg_sqr, counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8f99af7a-43bf-4bd9-8535-1374a4534dc7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# phase separation for the RMM1 and RMM2\n",
    "rmm = rmm_source[:, [3, 4]]  # RMM1 and RMM2\n",
    "\n",
    "phase_avg_rmm1 = phase_separation(rmm, rmm[:, 0])[0]\n",
    "phase_avg_rmm2 = phase_separation(rmm, rmm[:, 1])[0]\n",
    "\n",
    "# phase separation for the IP for the whole period of 1980—2020\n",
    "rmm = rmm_source[::3, [3, 4]]  # RMM1 and RMM2\n",
    "# the array should look like the IP data (with every third day taken)\n",
    "\n",
    "phase_avg_ip, phase_avg_sqr_ip, counter_ip = phase_separation(rmm, ip)\n",
    "\n",
    "# phase separation for the IP for 2006–2020\n",
    "# (when the PG data are available)\n",
    "d_start = dt.date(1980, 1, 1)  # first date (for the IP)\n",
    "d_pg = dt.date(2006, 1, 1)  # first date for the PG\n",
    "ip_pg = ip[((d_pg-d_start).days + 2)//3:]\n",
    "# the first index corresponds to the first day in 2006 for which\n",
    "# the IP data are available\n",
    "rmm = rmm_source[((d_pg-d_start).days + 2)//3*3::3, [3, 4]]  # RMM1 and RMM2\n",
    "# the array should look like the IP data\n",
    "# (with every third day taken, starting from 2006)\n",
    "\n",
    "phase_avg_ip_pg, phase_avg_sqr_ip_pg, counter_ip_pg = phase_separation(rmm, ip_pg)\n",
    "\n",
    "# phase separation for the PG\n",
    "rmm = rmm_source[np.ix_(fw_days_ind, [3, 4])]  # RMM1 and RMM2\n",
    "# the array should look like the PG data\n",
    "# (with only fair-weather days retained)\n",
    "\n",
    "phase_avg_pg, phase_avg_sqr_pg, counter_pg = phase_separation(rmm, pg_fw)\n",
    "\n",
    "# phase separation for the sunspot number for 2006–2020\n",
    "# (when the PG data are available)\n",
    "rmm = rmm_source[np.ix_(fw_days_ind, [3, 4])]  # RMM1 and RMM2\n",
    "# the array should look like the PG data\n",
    "# (with only fair-weather days retained)\n",
    "\n",
    "phase_avg_sn, phase_avg_sqr_sn, counter_sn = phase_separation(rmm, sn_fw)\n",
    "\n",
    "# phase separation for OLR for 1980—2020\n",
    "rmm = rmm_source[np.ix_(olr_days_ind, [3, 4])]  # RMM1 and RMM2\n",
    "# the array should look like the OLR data\n",
    "\n",
    "phase_avg_olr, phase_avg_sqr_olr, counter_olr = phase_separation(rmm, olr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c323b584-c1f9-4299-9fba-14009f743ad6",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Averaging of partial data sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "befd58ab-d5e8-4ebd-88e5-6e909da21585",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# phase separation for the IP for 1981–1990, 1991–2000, 2001–2010, 2011–2020\n",
    "phase_avg_ip_part = np.zeros((4, 8), dtype=float)\n",
    "phase_avg_sqr_ip_part = np.zeros((4, 8), dtype=float)\n",
    "counter_ip_part = np.zeros((4, 8), dtype=int)\n",
    "for i in range(4):\n",
    "    d_start = dt.date(1980, 1, 1)  # first date (for the whole IP dataset)\n",
    "    d_first = dt.date(1981 + 10*i, 1, 1)  # first date\n",
    "    d_last = dt.date(1991 + 10*i, 1, 1)  # last date\n",
    "    idx_first = ((d_first-d_start).days + 2)//3\n",
    "    idx_last = ((d_last-d_start).days + 2)//3\n",
    "    ip_part = ip[idx_first:idx_last]\n",
    "    rmm = rmm_source[idx_first*3:idx_last*3:3, [3, 4]]\n",
    "    # the array should look like the IP data\n",
    "\n",
    "    phase_avg_ip_part[i], phase_avg_sqr_ip_part[i], counter_ip_part[i] = \\\n",
    "        phase_separation(rmm, ip_part)\n",
    "\n",
    "# phase separation for the PG for 2006–2015 and 2011–2020\n",
    "phase_avg_pg_part = np.zeros((2, 8), dtype=float)\n",
    "phase_avg_sqr_pg_part = np.zeros((2, 8), dtype=float)\n",
    "counter_pg_part = np.zeros((2, 8), dtype=int)\n",
    "for i in range(2):\n",
    "    d_start = dt.date(1980, 1, 1)  # first date (for the whole IP dataset)\n",
    "    if i == 0:\n",
    "        d_first = dt.date(2006, 1, 1)  # first date\n",
    "        d_last = dt.date(2016, 1, 1)  # last date\n",
    "    else:\n",
    "        d_first = dt.date(2011, 1, 1)  # first date\n",
    "        d_last = dt.date(2021, 1, 1)  # last date\n",
    "    idx_first = (d_first-d_start).days\n",
    "    idx_last = (d_last-d_start).days\n",
    "    fw_idx_first = np.sum(np.where(fw_days_ind < idx_first, 1, 0))\n",
    "    fw_idx_last = np.sum(np.where(fw_days_ind < idx_last, 1, 0))\n",
    "    pg_fw_part = pg_fw[fw_idx_first:fw_idx_last]\n",
    "    rmm = rmm_source[np.ix_(fw_days_ind[fw_idx_first:fw_idx_last], [3, 4])]\n",
    "    # the array should look like the PG data\n",
    "\n",
    "    phase_avg_pg_part[i], phase_avg_sqr_pg_part[i], counter_pg_part[i] = \\\n",
    "        phase_separation(rmm, pg_fw_part)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de93dfa2-280f-4b7e-98b8-7da504bf7cec",
   "metadata": {},
   "source": [
    "## Estimating standard errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3e2b96da-d6db-429d-9f8c-8bc844718c9c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def std_error(avg_val, avg_sqr, counter):\n",
    "    \"\"\"\n",
    "    Estimate the standard error from the average value\n",
    "    and the average value of the square.\n",
    "    :param avg_val: the average value\n",
    "    :param avg_sqr: the average square value\n",
    "    :param counter: the size of the sample\n",
    "    :return: the standard error\n",
    "    \"\"\"\n",
    "\n",
    "    return np.sqrt((avg_sqr - avg_val**2) / (counter - 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b9634e8-fe13-4295-b76b-e1f00be37561",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Estimating correlation coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "cbfb26db-80ae-4bd7-bb7e-8f16015adf32",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Correlation coefficient between the IP and RMM1: 0.33046323421940305.\n",
      "Correlation coefficient between the IP and RMM2: -0.9300119158570614.\n"
     ]
    }
   ],
   "source": [
    "print('Correlation coefficient between the IP and RMM1: '\n",
    "      f'{np.corrcoef(phase_avg_ip, phase_avg_rmm1)[0, 1]}.')\n",
    "print('Correlation coefficient between the IP and RMM2: '\n",
    "      f'{np.corrcoef(phase_avg_ip, phase_avg_rmm2)[0, 1]}.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a9648327-d770-4534-823b-733ffe35b5e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Correlation coefficient between the PG and RMM1: -0.6786301924889234.\n",
      "Correlation coefficient between the PG and RMM2: -0.6688733027123144.\n"
     ]
    }
   ],
   "source": [
    "print('Correlation coefficient between the PG and RMM1: '\n",
    "      f'{np.corrcoef(phase_avg_pg, phase_avg_rmm1)[0, 1]}.')\n",
    "print('Correlation coefficient between the PG and RMM2: '\n",
    "      f'{np.corrcoef(phase_avg_pg, phase_avg_rmm2)[0, 1]}.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "814e19ad-7a71-4aef-95c9-6a471db693ec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Correlation coefficient between OLR and RMM1: -0.4962456609250617.\n",
      "Correlation coefficient between OLR and RMM2: 0.8422317866250696.\n"
     ]
    }
   ],
   "source": [
    "print('Correlation coefficient between OLR and RMM1: '\n",
    "      f'{np.corrcoef(phase_avg_olr, phase_avg_rmm1)[0, 1]}.')\n",
    "print('Correlation coefficient between OLR and RMM2: '\n",
    "      f'{np.corrcoef(phase_avg_olr, phase_avg_rmm2)[0, 1]}.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "4ef049d7-9d36-461f-a837-01c00423836b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Correlation coefficient between the IP and the PG: 0.5036181549956392.\n",
      "Correlation coefficient between the IP and the shifted PG: 0.9019059036035604.\n"
     ]
    }
   ],
   "source": [
    "print('Correlation coefficient between the IP and the PG: '\n",
    "      f'{np.corrcoef(phase_avg_ip_pg, phase_avg_pg)[0, 1]}.')\n",
    "print('Correlation coefficient between the IP and the shifted PG: '\n",
    "      f'{np.corrcoef(phase_avg_ip_pg, np.roll(phase_avg_pg, 1))[0, 1]}.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a046c8d5-6f7c-45ba-82c2-fca9d4456c57",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Significance level: 0.01.\n",
      "Number of points: 8.\n",
      "Critical value (((P − 2)r^2/(1 − r^2))^(1/2)): 3.707428021324907.\n",
      "Threshold correlation coefficient: 0.834341625597055.\n"
     ]
    }
   ],
   "source": [
    "P = 8  # the number of points\n",
    "a = 0.01  # significance level\n",
    "\n",
    "q = st.t.ppf(1 - a / 2, P - 2)\n",
    "r = q / np.sqrt(q**2 + P - 2)\n",
    "\n",
    "print(f'Significance level: {a}.')\n",
    "print(f'Number of points: {P}.')\n",
    "print(f'Critical value (((P − 2)r^2/(1 − r^2))^(1/2)): {q}.')\n",
    "print(f'Threshold correlation coefficient: {r}.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d9108a37-3893-439d-a59a-c645e6ee653b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum correlation coefficient between the IP and RMM is 0.987443102662834 at the angle 290°.\n",
      "Maximum correlation coefficient between the PG and RMM is 0.9520895875523926 at the angle 225°.\n",
      "Minimum correlation coefficient between OLR and RMM is -0.9781956110458861 at the angle 300°.\n"
     ]
    }
   ],
   "source": [
    "ip_rmm_corr = np.corrcoef(phase_avg_ip[np.newaxis, :],\n",
    "                          phase_avg_rmm1[np.newaxis, :] *\n",
    "                          np.cos(np.reshape(np.arange(360), (360, 1)) * np.pi/180) +\n",
    "                          phase_avg_rmm2[np.newaxis, :] *\n",
    "                          np.sin(np.reshape(np.arange(360), (360, 1)) * np.pi/180))[0, 1:]\n",
    "pg_rmm_corr = np.corrcoef(phase_avg_pg[np.newaxis, :],\n",
    "                          phase_avg_rmm1[np.newaxis, :] *\n",
    "                          np.cos(np.reshape(np.arange(360), (360, 1)) * np.pi/180) +\n",
    "                          phase_avg_rmm2[np.newaxis, :] *\n",
    "                          np.sin(np.reshape(np.arange(360), (360, 1)) * np.pi/180))[0, 1:]\n",
    "olr_rmm_corr = np.corrcoef(phase_avg_olr[np.newaxis, :],\n",
    "                          phase_avg_rmm1[np.newaxis, :] *\n",
    "                          np.cos(np.reshape(np.arange(360), (360, 1)) * np.pi/180) +\n",
    "                          phase_avg_rmm2[np.newaxis, :] *\n",
    "                          np.sin(np.reshape(np.arange(360), (360, 1)) * np.pi/180))[0, 1:]\n",
    "ip_rmm_angle = np.argmax(ip_rmm_corr)\n",
    "pg_rmm_angle = np.argmax(pg_rmm_corr)\n",
    "olr_rmm_angle = np.argmin(olr_rmm_corr)\n",
    "print('Maximum correlation coefficient between the IP and RMM is '\n",
    "      f'{ip_rmm_corr[ip_rmm_angle]} at the angle {ip_rmm_angle}°.')\n",
    "print('Maximum correlation coefficient between the PG and RMM is '\n",
    "      f'{pg_rmm_corr[pg_rmm_angle]} at the angle {pg_rmm_angle}°.')\n",
    "print('Minimum correlation coefficient between OLR and RMM is '\n",
    "      f'{olr_rmm_corr[olr_rmm_angle]} at the angle {olr_rmm_angle}°.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35ab8aa1-635b-4e52-ab8e-a14b4c2f29c1",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "92ca748a-eec7-4731-8d67-2c3d138cc0c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAFGCAYAAAAM69Z+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7UElEQVR4nO3de3gU5fn/8fdDEiDkRMBwCiERECpy0oKiIEGLBNRS/VZQtAVK9acgFQ/IKRGCoFULIq3wRREFFVprPaRA/KJYI6CIAnJUEQ9AOBkEEhJCyGHv3x9JtiGEZJPs7szu3q/rmkt2Z/aZeyYzH3d2Zp4xIoJSSqkLa2B1AUopZXcalEopVQMNSqWUqoEGpVJK1UCDUimlaqBBqZRSNbBFUBpjGhtjPjfGbDfG7DbGzLS6JqWUKmfscB2lMcYAYSKSZ4wJATYAE0TkM4tLU0opgq0uAEBK0zqv7GVI2WB9giulFDY59AYwxgQZY7YBWcAHIrLJ4pKUUgqwyTdKABEpAXoaY5oC7xhjuorIrvLxxpgBwIDKnxs4cOCM7OxsbrjhBho2bOilapWv2LZtGz179rS6DGUzR48eZffu3WzYsGGJiNxd4wdExHYDMB2Y6Mq0M2bMkKysLBk3bpzk5uaKUhXNmDHD6hKUzWzbtk2mTp0qJSUlAqSKCzlji0NvY0xM2TdJjDGhwA3AN65+PiYmhtTUVCZPnkxeXl7NH1BKBaTt27fzxhtvMHv2bBo0cD3+bBGUQGvgI2PMDuALSn+jXFWbBjQslVLVqWtIgk2CUkR2iMjlItJdRLqKyON1aUfDUilVlfqEJNgkKN1Jw1IpVVF9QxL8MChBw1IpVcodIQl+GpSgYalUoHNXSIIfByVoWCoVqNwZkuDnQQkalkoFGneHJARAUIKGpVKBwhMhCQESlKBhqZS/81RIQgAFJWhYKuWvPBmSEGBBCRqWSvkbT4ckBGBQgoalUv7CGyEJARqUoGGplK/zVkhCAAclaFgq5au8GZIQ4EEJGpZK+RpvhyRoUAIalkr5CitCEjQonTQslbI3q0ISNCjPoWGplD1ZGZKgQXkeDUul7MXqkAQNyippWCplD3YISdCgvCANS6WsZZeQBA3KamlYKmUNO4UkaFDWSMNSKe+yW0iCBqVLNCyV8g47hiRoULpMw1Ipz7JrSIIGZa1oWCrlGXYOSdCgrDV/CssBAwbw0ksvWV2GX9B1WXd2D0nQoKwTXwrLhIQEQkNDCQ8Pp2XLlowePdo2NR84cIDw8PBzBmMMc+fOtbq0Ktl5XZabP38+F198MWFhYVx66aV8++23VpdULV8ISdCgrDNfCsuVK1eSl5fH1q1b2bx5M7Nnz7a6JADatWtHXl6ec9i5cycNGjTgt7/9rdWlXZBd1yXASy+9xJIlS1i9ejV5eXmsWrWKiy66yOqyLshXQhI0KOvFl8ISIDY2liFDhrBr1y7ne/v376dv375EREQwaNAgfv75Z+e4YcOG0apVK6Kioujfvz+7d+92jktPT6dLly5EREQQGxvLnDlznONWrVpFz549adq0Kddccw07duxwqb5XX32V/v37k5CQUP+F9TC7rUuHw8HMmTOZN28eXbp0wRhDhw4daNasmQeWvv58KSQBEBGfHmbMmCFWy8rKknHjxklubq7VpZwnPj5ePvjgAxEROXDggHTp0kVSUlJERCQxMVHat28ve/bskfz8fElMTJTJkyc7P7tkyRI5deqUFBQUyIQJE6RHjx7Oca1atZJ169aJiMiJEydky5YtIiKydetWiYmJkc8++0yKi4tl6dKlEh8fLwUFBdXW6XA4pH379vLKK6+4celF3Ll92Hld7t+/XwB57rnnpG3btpKQkCDTp0+XkpISty2/u2zbtk2mTp1qi9qAVHEhZywPuvoOdghKEfuGZXx8vISFhUlUVJS0a9dOxo4dK/n5+SJSunPPmjXLOe2CBQskKSmpynZOnjwpgGRnZ4uISFxcnCxatEhycnLOme6+++5zhke5Tp06SUZGRrV1rlu3TsLCwty+/twdlHZdl5988okAcuONN8rJkyflxx9/lEsuuURefPHFei2zu9kpJEVcD0of+M7rG+x8GP7uu++SnZ3N/v37WbhwIaGhoc5xrVq1cv67SZMmztpLSkqYMmUKHTp0IDIy0nk4XH44+dZbb5Genk58fDyJiYls3LgRKD38nDt3Lk2bNnUOmZmZHD58uNoaly1bxm9/+1vCw8PduehuZ9d1WV7HpEmTaNq0KQkJCdx7772kp6d7ZD3Uhc8dblfgW9XanJ3DsrZWrFhBWloaa9euJScnh3379gGlRyAAvXv3Ji0tjaysLG655RaGDx8OQFxcHMnJyWRnZzuH/Px8RowYccF5nTlzhjfffJNRo0Z5fLms4I112blzZxo2bIgxxvlexX9bzZdDEjQo3c5fwjI3N5dGjRrRvHlz8vPzmTZtmnNcYWEhy5cvJycnh5CQECIjI50b/z333MOiRYvYtGkTIsLp06dZvXo1ubm5F5zXO++8Q3R0NNddd53Hl8sK3liXTZo04fbbb+eZZ54hNzeXgwcP8uKLL3LzzTd7bTkvxNdDEjQoPcIfwnLkyJHEx8cTGxtLly5d6NOnzznjX3vtNRISEoiMjGTRokUsX74cgF69erF48WLGjx9PdHQ0HTt2ZOnSpdXOa9myZfz+97+31Tcgd/LWunz++ecJDw+nTZs2XH311dx5552MGTPGk4tWI38ISUBP5niSXU/wBBI7bx/+zm4nbqqCnsyxnj98s1SqLvzmm2QZWyyBMSbOGPORMeYrY8xuY8wEq2tyFw1LFWj8LSTBJkEJFAOPiEgXoA9wvzGmi8U1uY2GpQoU/hiSYJOgFJEjIrK17N+5wNdArLVVuZeGpfJ3/hqSYJOgrMgYkwBcDmyyuBS307BU/sqfQxIg2OoCKjLGhANvAQ+KyKlK4wYAAyp/xtcuUq4Ylk8//bTt70RRqib+HpJgo2+UxpgQSkNyuYi8XXm8iGSISGrlwRd6mqlMv1kqfxEIIQk2CUpTeqXxEuBrEXnW6nq8QcNS+bpACUmwSVACfYHfA9cbY7aVDTdaXZSnaVgqXxVIIQk2CUoR2SAiRkS6i0jPssE+3Z54kIal8jWBFpJgk6AMdBqWylcEYkiCBqVtaFgquwvUkAQNSlvRsFR2FcghCRqUtqNhqewm0EMSNChtScNS2YWGZKnAXXKb07BUVtOQ/K/AXnqb07BUVtGQPJeuAZvTsFTepiF5Pl0LPkDDUnmLhmTVdE34CA1L5Wkakhema8OHaFgqT9GQrJ6uER+jYancTUOyZrpWfJCGpXIXDUnX6JrxURqWqr40JF2na8eHaViqutKQrB1dQz5Ow1LVloZk7ela8gMalspVGpJ1o2vKT2hYqppoSNadri0/omGpLkRDsn50jfkZDUtVmYZk/ela80MalqqchqR76JrzUxqWSkPSfXTt+TENy8ClIeleugb9nIZl4NGQdD9diwFAwzJwaEh6hq7JAKFh6f80JD1H12YA0bD0XxqSnqVrNMBoWPofDUnP07UagDQs/YeGpHfomg1QGpa+T0PSe3TtBjANS9+lIelduoYDnIal79GQ9D5dy0rD0odoSFpD17QCNCx9gYakdXRtKycNS/vSkLSWrnF1Dg1L+9GQtJ6udXUeDUv70JC0B9useWPMy8aYLGPMLqtrURqWdqAhaR92WvtLgcFWFzFgwABeeuklq8uwBQ1L69g5JANxH7HNX0BE1gEnvDGvhIQEQkNDCQ8Pp2XLlowePdpWQbBt2zauvfZaoqKiaNu2LbNmzbKsFl8Ly+LiYg4fPszhw4f5+eefOXv2LGfPnkVErC7NZXYISbvvIxXrCw8PZ9CgQR6dX7BHW3cjY8wAYEDl90eNGlWn9lauXMnAgQM5dOgQSUlJzJ49m6eeeqpeNbrLnXfeya233kpGRgb79u2jX79+9OjRg6FDh1pST8WwfPrppwkPD7ekjnIiwtGjR9m5cyc7d+4kJycHYwzGGIKCgmjRogXGGM6ePcuWLVt49tlnKSwsdIalw+FwttO2bVt69uzJ5ZdfTsOGDa1cLMAeIVnOzvsI/Lc+rxAR2wxAArCrNp+ZMWOG1FZ8fLx88MEHztcTJ06Um266SUREEhMTJSUlRa655hoJDw+XG264QY4dO+ac9rbbbpOWLVtKZGSkXHvttbJr1y7nuNWrV8ull14q4eHh0qZNG/nLX/7iHLdy5Urp0aOHREVFydVXXy3bt2+/YH2hoaGye/fuc+b55JNP1no53S0rK0vGjRsnubm5Xp+3w+GQTZs2SWpqqkyZMkUWLlwo69evl5MnT1b7ueq2D4fDIQcPHpSVK1fK448/LikpKTJz5kxZv369OBwO9y6AC7Zt2yZTp06VkpISr8+7MrvvI5XrqysgVVzJJlcm8tZgRVAeOHBAunTpIikpKSJSuhG0b99e9uzZI/n5+ZKYmCiTJ092fnbJkiVy6tQpKSgokAkTJkiPHj2c41q1aiXr1q0TEZETJ07Ili1bRERk69atEhMTI5999pkUFxfL0qVLJT4+XgoKCqqsb+rUqTJ58mQpLCyUb775RmJjY+Xzzz+v9XJ6gjfDsqSkRD799FNJTU2VqVOnyttvvy35+fm1aqO228fZs2clLS1NJk2aJM8++6wcOXKkVp+vKzuFpIj995H4+Hhp0aKFXHTRRXLDDTfItm3b6rScGpTViI+Pl7CwMImKipJ27drJ2LFjnTtgYmKizJo1yzntggULJCkpqcp2Tp48KYBkZ2eLiEhcXJwsWrRIcnJyzpnuvvvuc25k5Tp16iQZGRlVtvvJJ59Ihw4dJCgoSACZPn16rZfRkzwdlpmZmfL444/L1KlTJS0tTc6cOVPntuqyfZQ7fPiwzJ07VyZNmiRpaWlSVFRU57aqY7eQFLH/PrJhwwbJz8+X06dPy5NPPiktW7as8eiiKj4XlMDfgSNAEXAQ+KMrn3PHoXdFiYmJsnjxYufrV155Rfr27SsiIsXFxTJ58mRp3769RERESFRUlADy3XffiYjI559/LkOHDpWmTZtK//795dNPPxURkSFDhkhoaKhERUU5h9DQUFmxYsV58z9+/LhERETIsmXLpKioSDIzM+Wqq66SBQsW1Ho5PckTYXngwAF57LHHZM6cOeftSHVVn6As53A4ZP369TJt2jR55pln3LrMdgxJEXvvI1Xp3Lmz/Pvf/671croalLY5mSMiI6yuoSYrVqwgLS2NtWvXkpCQQE5ODtHR0eVBT+/evUlLS6OoqIjnn3+e4cOHk5mZSVxcHMnJySQnJ9c4jx9++IGgoCBGjhwJQNu2bbnjjjtIT09n3LhxHl2+2nDnCZ6TJ0/y17/+lfDwcB599FEiIiLcWGn9GWPo168f/fr14+DBg8yePZvOnTszcuRIgoKC6tyunU7cuIs39pGqGGOc8/AE//jreElubi6NGjWiefPm5OfnM23aNOe4wsJCli9fTk5ODiEhIURGRjo3/nvuuYdFixaxadMmRITTp0+zevVqcnNzz5tHp06dEBFWrFiBw+Hg6NGjvPHGG3Tv3t1ry+mq+l46VFxczAsvvMDcuXMZN24cjzzyiO1CsrK2bdvy1FNP0b17dx5++GE2bdpUp3b8MSTBO/vIgQMH+OSTTygsLKSgoIC//OUv/Pzzz/Tt29djy+U/fyEvGDlyJPHx8cTGxtKlSxf69OlzzvjXXnuNhIQEIiMjWbRoEcuXLwegV69eLF68mPHjxxMdHU3Hjh1ZunRplfOIjIzk7bffZt68eURHR9OzZ0+6du1KSkqKpxevTuoalt9//z0TJkzg6quvZvbs2cTExHiwSvf75S9/ybx58/jqq69ITU2t1bL7a0iCd/aR3Nxcxo4dS3R0NLGxsfzf//0f7733Hs2bN/fcgrlyfG7nwR2/Qan6q81vlitXrpTp06fL2bNnPV6XN7aPw4cPy8MPP+z8va06dv1NMlDh4m+U/vW/M2UZV75ZlpSU8PTTT5OXl8fMmTNtcYG3O7Ru3Zo5c+awadMm/vWvf11wOn/+Junv9K+l3Ka6sDx+/DgPPfQQN954I3fccYdFFXqOMYYHH3wQh8PB/PnzzzuxoCHp2/QvptyqqrDcsmULs2fPZtasWXTr1s3iCj1r+PDh9OrVi+TkZIqKigANSX9gm8uDlP+oGJa33HILn3/+OXPnzg2YkOjbty8tW7bk4Ycf5vbbbyc9PV1D0sdpUCqPiImJ4Te/+Q0pKSmsXbs24EKiY8eO3HbbbUyYMIHVq1cH3PL7G/3rKY/YuXMn69evZ+XKlUyZMsVWXXR5w/bt21mzZg0fffQRs2bNIj8/3+qSVD1oUCq3O3ToEEuXLmXmzJm0aNHCp/qzdIeKv0lGRkaSkpJCcnIyJSUlVpem6kiDUrlVXl4eTz755Dm/yfla57/1UdWJm9atW3P33Xfz5JNPWlydqisNSuU2DoeDlJQUZsyYQWho6DnjAiEsqzu7fdlll9GjRw/efPNNi6pT9aFBqdzm5Zdf5q677qJFixZVjvfnsHTlEqChQ4eyd+9evv76ay9Xp+pLg1K5xeHDhzl48CC9e/eudjp/DMvaXCc5adIkFixYQEFBgZeqU+6gQanq5MEHH2TJkiXO5888++yzTJw40aXP+lNY1vZi8uDgYB566CEWLVrkheqUu2hQqjrZuXMn999/P1deeSX33nsv1113Xa36pPSHsKzrHTcdOnQgNzeXn376yYPVKXfSoFR1Vv6UwxdffJHU1FRef/31WnWe6sthWd/bEh944AGef/55D1SmPEGDUtVJ5UDcvHkzd999N8OHD69VO74Ylu64dzsqKoq4uDh27tzp5uqUJ2hQKrcIDQ1l4MCBLFmypNaf9aWwdGcHF2PGjOHll1/26CMMlHtoUKo6qbhzt2nThscee4yVK1cSGRlZp/Z8ISzd3QtQcHAwgwcPZvXq1W6oTnmSBqWqk7NnzwLQvXt30tLSmDp1KsaYerVp57D0VFdpSUlJ/Oc//9FvlTanQanq5MCBA9x444188skn9OrVy23t2jEsPd2f5PXXX8+GDRvc3q5yHw1KVScffvghq1evrtdjai/ETmHpjU53Bw8eTHp6ukfaVu6hQanqpHPnzh5t3w5h6a2eyYODgwkPDyc7O9tj81D1o0GpbMvKsPT24xuGDRvG22+/7fH5qLpxeQswxkRf4P227itHqXNZEZZWPOOmU6dO7N271yvzUrVX41ZgjOlkjPkaOG6MOWSMqXxF8VeeKU2pUt4MSysfBBYVFUVOTo5X56lc48qWMB94E2gO3A88a4yZUmF8/a4JUcoF3ghLq5+WePPNN7Nq1Sqvz1fVzJWtoTcwU0ROisi7QB/g98aYJzxamVKVeDIsrQ5JKO3c95tvvrFk3qp6rmwRDiCi/IWIHAQGAEOMMX/1UF1KVckTYWmHkATqfcG+8hxXtopPgVsrviEix4DrKf122cQDdSl1Qe4MS7uEZDljjN6lY0OubBmPArsqvyki2cCvgD+4uSalauSOsLRbSAK0atWKo0ePWl2GqqTGrUNE9orIFxcYlysir7q/LKVqVp+wtGNIAnTr1o1du877XqIs5srlQS/XNHijUKWqUpewtGtIAnTt2lX7qLQhV7aS0cDVwDHg0AUGpSxTm7C0c0gCREdH662MNuTKlnIrpReVjwQuAtJF5LGKg0crVMoFroSl3UNS2Zcrv1GmichvgS7ADkovON9jjEkxxjR1VyHGmMFl7X5X6YJ2pVxSXVj6Ukjavb5A5PJfpOyC8/8FkoA0IBW43B1FGGOCgAXAEEoDeYQxpos72laBpaqw9KWQVPbk0lZjjGlgjBlijPkHsAdoBvxKRD5yUx1XAt+JyA8iUgj8A/iNm9pWAaY8LCdNmkpmZqaGpKq34JomMMbMBYYDO4FXgdEiUuDmOmKBzAqvDwJXVapjAKV3BJ2jb9++bi5F+YNiRzD7gm7ju+Yfc73J5fHHH7e6JJd99JG7vn8od6kxKIGHgO8pvY1xHDCu8q1WItLf/aWdS0QygIzK76emps5ITU319OyVD8kvcDBlwTEKTuRzOtTBe4Sw7aGHiI6Ksro0l+k27R0zZ850aTpXgtIbd94cAuIqvG6LXnak6qA8JL/68SwzxrRk979/4n7Tit/Mepz01JkeeXSF8n81BqWILPNCHV8AlxhjLqY0IO8A7vTCfJUfqRiSj425iMQrmnDk340YKLl80awZD0yZwl+fesr2Yan3etuPK79RjqxpmvrexigixcaY8cAaIAh4WUR216dNFViqCkmAvUQyr+FBehZFwHXXMXnyZJ5++mlbh6UGpf24cui9FPgOOErVnfQKpSd56kVE0gF9FJ2qk8PHijlwtOickAQooQFdu/Vm7JffsHD3bv7zpz/ZOiw1JO3J1R7OmwC5wPPAQBG5tsLg8RM5StWkY1xDlj/e5pyQdOoxgJmOQzQNDuLTY8csf7pjdY4cOULr1q2tLkNV4sqdOQ8B8cBC4H+AfcaYxcaYfp4uTqnaCAu9wOYc9wua9R3K1yNuZeq119riUbgXsmvXLrp162Z1GaoSl67AFZESEVktIrcDvwBOAhnGmOs8Wp1S7mAM3DCSFhf/92Yvu4blzp07ueyyy6wuQ1VSm8fVRhlj7gXeo7SjjFnANg/VpZTH2TEss7OziY6u8snQykKu9Ef5a2PMm8A3lN7b/aiIXCIiM0XkpMcrVKqC/AIHu74/67b27BiWyn5c+UaZBnQFXgeygCRjzOMVB49WqFSZ/AIHUxccY/LzWeTklbitXbuEZUlJid6PblOu/FVeBT6jtC/KuCqGth6rTqky5SG5+8ezTPp9c6LCg9zavh3CcuPGjVx55ZWWzFtVz5Wz3qNF5A9VDcA8KjzKVilPqBiSla+TrI/Xd+zg13//u/PaRavDcu3atQwcONDr81U1c+U3yibGmFnGmJXGmGeNMZHGmPbGmHeATyg9HFfKI9wekgX58OFy2LebopISVn37LSsqPKPGqrAUEQoLC2nUqJHX5qlc58qh9wLg15Q+DmIg8BbwMbAbuFhE7vdceSqQlTiEaQvd/E0yOAS++D/Y8j6jevbkl61bM2ntWvIKC52TWBGWmzZtonfv3l6Zl6o9V4IyCRgkIpOBGyl9lvedIpIiIj97tDoV0IIaGJKuDnPr4TbBIdC1L3z9GQ0KC5g/eDCHc3N5esOGcybzdliuWrWKG2+80ePzUXXjSlCGi0gWgIgcBPJEZL1ny1J29+2333plPkOuDndfSJbrPgCKC+HrjfRt144RXbsyZ+NG9lV6+qG3wvL06dOEhIToYbeNuRKUwcaY64wx1xtjrgeo+Lr8PRVY1qxZw48//mh1GXUT1xmiW8H2jwF4euBADDDpgw/Om9QbYfnmm28ybNgwj7St3MOVoMwCXgaWlA3HK71+yWPVKdsaPXo0y5Z5o6tSDzAGuifCjzsh52fioqKY0q8fn2Rmcjw//7zJPRmWJSUl7Nixgy5d9Fl6dubK5UEJInJxNUN7bxSq7CUiIoKgoCCyKx2u+ozuiYDAznUAPHrNNewZP57mTao+zPdUWC5btozf/e53bmtPeYbeBqDqbMyYMSxcuLDe7eQXOHj8pZ85lFXkhqpc1Lw1XHUTXFR6v0RoSAjhDRtW+xF3h2Vubi579+7liiuuqHdbyrM0KFWdxcbGEhUVxc4K1yHWVnnP5Ou25bPviBeDEmDI3fCL2t0J486w/Nvf/saf/vSnerWhvEODUtXLfffdx+LFiykpqf2915Uf39C3h5vPbnuIO8Jy3759NGzYkDZt2ri5OuUJGpSqXoKCgvjDH/7ASy/V7pzehZ5x4yvqG5YLFy5k3LhxHqhMeYIGpaq3yy+/nGPHjnHw4EGXpvf1kCxX17D89NNP6dmzJ00ucOJI2Y8GpXKLhx56iHnz5rk07Xuf5tk+JLNOn+bG5ctJ37u32ulqG5Y5OTn84x//4I477nBXqcoLNCiVW4SFhTFo0CD+9a9/1Tjt/1wXwcJJrWwbkgBNGzfmuxMneHjNGopq+P3V1bAUEVJTU5kxY4b2O+lj9K+l3CYpKYk9e/awffv2aqczxtCpXfWX4nhVSQls+w98t835VsOgIJ5NSmLP8eMs+OKLGptwJSyff/557rzzTpo3b+6uypWXaFAqt5o6dSpLly7l8OHDVpfiugYNYN2/4JO3z3n7pksuIalDB1IzMjh2+nSNzVQXlqtXryYqKkp7CPJRGpTKrRo0aMATTzzBE088QW5urtXluMZ5S+MuyPm5wtuGeUlJ5BUW8thHH7nUVFVhuWXLFnbt2sXIkSM9Ur7yPA1K5XZNmjQhNTWVRydP5+cTPvLAru79AYEdH5/z9qUxMYy/8kpe3LKFbUePutRUxbD86quv+Oc//8mkSZM8ULTyFg1K5RFhEc0h7mFGp2yjoKDA6nJq1qw1xP2iNCjLHg1RbkZiIs2bNGHNd9+53FxMTAzDhg1jxIgRTJo0CWOMuytWXqRBqdyu/DrJ7w4Jf7ilHY8++iinTp2yuqya9RgAxzLhyA/nvB0dGso399/P5H79XG5qzZo1bNiwgffff5/p06fro3B9nAalcqvKF5P/dlA7pk+fTnJyMl9//bXV5VXvsr4QFHze4TdwwV6FqrJkyRIOHTpESkoKLVu2tPzpjqr+NCiV21zojpuYmBiee+453nnnHd58802Lq6xGaDj8Mgki63b5jsPhYPbs2cTGxjJmzBjn+1Y/3VHVnwZlJQMGDKj1fcuq5tsSg4KCmDZtGiEhIcycOZOiIi/3FOSqG++Ga35T64+dOXOGiRMnMnToUAYPHnzeeH8Ky0DcRwIyKBMSEggNDSU8PJyWLVsyevRoW228jz32GN26dSM4OJjU1FSry3GJAUKCqfG2xFtuuYURI0bw4IMPctTFs8h2t27dOqZMmcLEiRPp3r37BafzpbC0+z5S7uOPP8YYQ0pKikfnE5BBCbBy5Ury8vLYunUrmzdvZvbs2VaX5NSxY0eeeeYZbrrpJqtLcVlo4wbMeaCFS7cldurUiWeeeYa//vWvvPbaazgcDi9U6H4nT54kOTmZzMxMnnvuOZe6TPOlsLTzPgJQVFTEhAkTuOqqqzw+r4ANynKxsbEMGTKEXbt2Od/bv38/ffv2JSIigkGDBvHzz/+9CHnYsGG0atWKqKgo+vfvz+7du53j0tPT6dKlCxEREcTGxjJnzhznuFWrVtGzZ0+aNm3KNddcw44dOy5Y06hRoxgyZAgRERFuXlrPqs0lMGFhYTz55JNccsklPPTQQ3z66acerMy9iouLGTN/Pr987jkemDCBu+66q1bL7kthCfbcRwDmzp3LoEGD+MUvfuHGpa1awAdlZmYm6enpXH755c73VqxYwSuvvEJWVhaFhYXn/DGHDBnC3r17ycrK4oorruCuu+5yjvvjH//ICy+8QG5uLrt27eL660sfUPnll18yZswYXnjhBY4fP869997L0KFDOXv2rPcW1Kb69OnDvHnz2LdvHxMnTmTbtm1Wl3RBIsKaNWuYOHEiF198MT82aMCaI0fq1JYvhaUd95H9+/fz8ssvM336dA8tdSUi4tPDjBkzpLbi4+MlLCxMoqKipF27djJ27FjJz88XEZHExESZNWuWc9oFCxZIUlJSle2cPHlSAMnOzhYRkbi4OFm0aJHk5OScM919990nKSkp57zXqVMnycjIqLbOu+66S+qyfL6qsLBQXnnlFXnkkUdk8+bN4nA46t2mO9bfsWPHZP78+fLoo4/KqlWrxOFwSInDIVcuXiyt58yR3LNn69x2VlaWjBs3TnJzc+tdpzvZfR8ZOnSo/OMf/xARkVGjRklycnKdlhNIFRdyxvJvlMaYYcaY3cYYhzGml7fm++6775Kdnc3+/ftZuHAhoaGhznGtWrVy/rtJkybO/+OXlJQwZcoUOnToQGRkJAkJCQDOw4633nqL9PR04uPjSUxMZOPGjUDp//3mzp1L06ZNnUNmZqZvdRxRJr/AwZsfnsLhkJonrqWQkBBGjx7Nn//8Z7799lumTZvGE088wY4dOxBx//yqU1JSwsZ/vsKSqeN55ZVXuPXWW52/GxtjaGAM8wcP5kheHn9ev77O87HzN0u77iMrV64kNzeX22+/3VOLfp5gr83pwnYB/wO8YHUhNVmxYgVpaWmsXbuWhIQEcnJyiI6Odu7EvXv3Ji0tjaKiIp5//nmGDx9OZmYmcXFxJCcnk5ycbPES1E/FS4C6d2xE5/hGHplPSEgII0aMYMSIEZw6dYrVq1fz97//nbCwMG6++WZ69OjhkVsCy5+x/d5775Gbm8vkmBz6xAZh7p9Y2nFGJX3atuV33bszd+NG/njFFbSPjq7TfCuG5dNPP014eHh9F8Uy3thHPvzwQzZv3uwM65ycHIKCgti5cydpaWkeWS7Lg1JEvobanQiwSm5uLo0aNaJ58+bk5+czbdo057jCwkLefPNNbr75ZqKiooiMjHR2znrPPfdw6623MnDgQK688kry8/PJyMigf//+VZ6wKSoqoqSkBIfDQXFxMQUFBYSEhBAUFOS1Za2s8nWSngrJyiIjI52hmZeXx+rVq/nnP/9JcHAwzZo1o1u3bnTr1o2YmJhabUMlJSX8+OOP7N69m927d3PmzBkaNGhAt27dGDt2LNHR0bB5DaxaVHpLY5sOVbbz1K9+xdtff82jH3zAW8OH13k5/SUsvbGPzJo1iylTpjhfT5gwgTZt2vDYY495bLksD0pfMnLkSNasWUNsbCzNmjVj1qxZ/O///q9z/Guvvcb48eMpKSmhc+fOLF++HIBevXqxePFixo8fz969ewkNDaVfv37079+/yvncc889LFu2zPn6iSee4JVXXmH06NEeXb4LscszbsLDw7n99tudh1zHjx9n586dvPHGG85DO2MMIuIMzYyMjPOuRXU4HDRo0ID27dtz2WWXMWjQoHMOK50u6wvvvQTbMy4YlLGRkUzr148XtmzheH5+rW51rMwfwtIb+0hERMQ54RkaGkpYWBjNmjXz2HIZb/z2Y4xZC7SqYlSyiKSVTZMBTBSRzRdoYwAwoPL7o0aNmrF06VI3Vaoqs0tI1lVqamr9Ltp/4xk48BU8/FLpfeBVKCguxiFCk5CQus+ngmPHjpGamuqzYelLjDEzRSS1pum8cjJHRAaKSNcqBpd/UBCRDBFJrTyU/1isPGPuihM+G5Ju0SMRTufA99suOEnj4GC3hSTY+wRPoLL8rLeytz8ObcrMewI0JAE6XgGhEbD9/B6FPEnD0l4sD0pjzK3GmIPA1cBqY8waq2tS/9XmomD69gjQkAQIDoGe10Nw8Hkd+nqahqV9WB6UIvKOiLQVkUYi0lJEkqyuSalzJI2GWydUeYmQp2lY2oPlQamUPyq/o8MdNCytp0GpACjxwJ02gSrr9GmSXn+d9L173damhqW1NCgV+QUOHn4ui9Wf6A7oDk0bN+ZATg4PrVlDYUmJ29rVsLSOBmWAK79OcvcPZwkP1c3BHRoGBTEvKYm9J07wt02b3Nq2hqU1dM8IYL5+MbmdDbnkEm685BIeX7eOn9wcaBqW3qdBGaA0JOvo1An4ofoOZcs9O2gQ+UVFpPznP24vQ8PSuzQoA5CGZD383xL411woKa5x0s4XXcQDV17Jki+/5Ms6dvBbHQ1L79GgDED7jhTxw6FCDcm66J4I+aeqvaWxoscSE4kJC2PtDz94pBwNS+/QoAxAXS5uxPLH22hI1kXHy2t1S2PTxo3ZM348j/bt67GSNCw9T4MyQEWFW9e3pU8LDoGu/WDP51Bw2qWPNG3c2MNFaVh6mgalUrXVYwAUF8JXG62u5Bwalp6jQalUbcVeAs1al3boazMalp6hQenH8gscfHug0Ooy/I8x0OM6KDwDRfZbvxqW7qdB6afKLwF69K9Z5OU7rC7H/1z7P3DvXAhpWKePe/rJAhqW7qVB6YcqXif50J3NCG+if2a3a1D3k2FLtm7lhtdew6Fh6TN0D/IzFUMyZcxFDNBLgGyncXAwH/74I8u2bfP4vDQs3UOD0o9oSPqGO7t1o0/btkz98ENOnT3r8flpWNafBqWfKCkRpmpI+gRjDPMHD+an06d5cv16r8xTw7J+NCj9RFCQ4drLQzUkfcSVsbGM6tGDeZ99xncnTnhlnhqWdadB6Uduuz5SQ9IKxUVwqPa9mf/5V7+iYVAQE99/3wNFVU3Dsm40KJWqrw+Xw8vJcMa1WxrLtY6IIPnaa9n+008cz8/3UHHn07CsPQ1Kpeqra18oKYKvPq31Rx/q04evxo2jeRPvHgloWNaOBqVS9dWmIzSPhR0Ztf5oo+BgQkNC3F+TCzQsXadB6WPyCxw8tew4x7Jr7jhWeYkx0CMR9n8FJ7OsrqZWNCxdo0HpQ8qvk1z7xWm+yyyyuhxVUbfE0v/udK2fSjvRsKyZBqWPqPz4hqu7hVpdkqoougXEX1baoa+Hb030BA3L6mlQ+gB9xo2P6DEA8nPg1PF6NSMiHu80oyoalhemQWlzGpI+pHsiPPIyRF1U5yZ+ysvj+ldfZeW337qxMNdpWFZNg9Lm3snI1ZD0FcEhpUM9NAsN5ae8PB5es4azxdacsNOwPJ8Gpc3dMSiS+Q+31JAMECFBQcxLSuL7kyeZv2mTZXVoWJ5Lg9LmghoYLmvfyOoylBcldezIrzt1Yta6dRy1MKQ0LP9Lg1IpG5o7aBBni4uZ9uGHltahYVlKg1IpG7qkeXMmXHUVS7dtY/Phw5bWomGpQWkbZwocFBb53vV3qhonf6rXx1P696dleDjr9u93U0F1F+hhqUFpA+WXAKUuPmbJ9XPKAza/D/Pvq9ctjVGNG7Nn/HgevvpqNxZWd4EclhqUFssvcDB1wTF2/3iWpD7hGGOsLkm5Q4eepf/dUb9bGiMb2etEXqCGpQalhSqGpF4n6WectzR+5JO3NFYnEMPS8qA0xvzFGPONMWaHMeYdY0xTq2vyBg3JANBjAJw4Uqfez+0u0MLS8qAEPgC6ikh34FtgqsX1eJyGZIDocjUEN4TtGVZX4hGBFJaWB6WIvC8i5fdqfQa0tbIebyhxQHGJaEj6u8Zh0Lk37NpQ+lwdN7HTCb9ACUvLg7KSMcB7VY0wxgwwxqRWHvbt2+fdCt0gokkD/jZRb0sMCD0GQMFpOPydW5p7ccsWBixbRonD4Zb23CEQwtIrQWmMWWuM2VXF8JsK0yQDxcDyqtoQkQwRSa08JCQkeGMR3K5BAz27HRA69ISHXoR2l7qluchGjVi3fz+vbNvmlvbcxd/D0itBKSIDRaRrFUMagDFmNHAzcJfY6bhCqfoKCobI5m5r7vbLLqNvXBzJ//kPOQUFbmvXHfw5LC0/9DbGDAYmAUNFxHvP7FTKBxlj2DBmDD9NnEhU48ZWl3Mefw1Ly4MSeB6IAD4wxmwzxiyyuiB3yS9wsHJ9rq1+fK9owIABvPTSS1aX4Rd0Xf6XP4al5UEpIh1FJE5EepYN91ldkzuU35Y4/42T/HDIugeBJSQkEBoaSnh4OC1btmT06NG22XizsrIYMWIEbdq0ISoqir59+7LJwj4Ya2LndQlw3XXXERMTQ2RkJD169CAtLc2yWvwtLC0PSn9U+fENHdo2tLSelStXkpeXx9atW9m8eTOzZ8+2tJ5yeXl59O7dmy1btnDixAlGjRrFTTfdZOsdy67rEmD+/PkcOXKEU6dO8eKLL/K73/2OI0eOWFaPP4WlBqWb2fkZN7GxsQwZMoRdu3Y539u/fz99+/YlIiKCQYMG8fPPPzvHDRs2jFatWhEVFUX//v3ZvXu3c1x6ejpdunQhIiKC2NhY5syZ4xy3atUqevbsSdOmTbnmmmvYsWNHlfW0b9+ehx9+mNatWxMUFMT/+3//j8LCQvbs2eOBpXcvu61LgO7duxMcHAyU/pZZVFREZmamOxe71vwmLMuf+Oarw4wZM8QuTp8pkT/NOSq/un+/ZGw5bXU5IiISHx8vH3zwgYiIHDhwQLp06SIpKSkiIpKYmCjt27eXPXv2SH5+viQmJsrkyZOdn12yZImcOnVKCgoKZMKECdKjRw/nuFatWsm6detEROTEiROyZcsWERHZunWrxMTEyGeffSbFxcWydOlSiY+Pl4KCghpr/fLLL6VRo0aSnZ3trsUXd24fvrAub7rpJmnUqJEAkpSUJCUlJW5b/vrIysqScePGSW5urtWlnANIFRdyxvKgq+9gp6Cc/kKWrUJSpHTnDgsLk6ioKGnXrp2MHTtW8vPzRaR05541a5Zz2gULFkhSUlKV7Zw8eVIAZ4jFxcXJokWLJCcn55zp7rvvPmd4lOvUqZNkZGRUW2dOTo507dpVnnzyyVovY3XcHZS+sC4LCwslPT1d5s6dW+tl9CQ7hqWrQamH3m40+uYoZtxtr8NtgHfffZfs7Gz279/PwoULCQ0NdY5r1aqV899NmjRxHh6VlJQwZcoUOnToQGRkJOUX9pcfTr711lukp6cTHx9PYmIiGzduBEoPP+fOnUvTpk2dQ2ZmJoer6aX7zJkz/PrXv6ZPnz5MnWrvW/3tvi4BQkJCGDJkCO+//z7//ve/3bn49eLLh+EalG50cZuGXNvTXiFZVytWrCAtLY21a9eSk5ND+a2iUnapU+/evUlLSyMrK4tbbrmF4cOHAxAXF0dycjLZ2dnOIT8/nxEjRlQ5n7Nnz3LLLbfQtm1bXnjhBa8sm7d5a11WVlxczPfff++RZaorXw1LDUpVpdzcXBo1akTz5s3Jz89n2rRpznGFhYUsX76cnJwcQkJCiIyMpEGD0k3pnnvuYdGiRWzatAkR4fTp06xevZrc3Nzz5lFUVMRtt91GaGgoy5Ytc7bhb7yxLr/55hvee+89zpw5Q1FREa+//jrr1q0jMTHRa8vpKl8MS//cMlW9jRw5kvj4eGJjY+nSpQt9+vQ5Z/xrr71GQkICkZGRLFq0iOXLS2/R79WrF4sXL2b8+PFER0fTsWNHli5dWuU8Pv30U1atWsX7779P06ZNCQ8PJzw8nPXr13t68bzKG+tSREhNTaVFixbExMQwf/583njjDa644gpPL16d+FxYuvJDpp0HK07mlJQ4vD5PVTd2Otmnzmf1CR70ZI5n5Bc4eGR+FhlbTltdilI+z1e+WWpQ1kL5xeQ7vz+rDwFTyk18ISw1KF1k5ztulKro+xMnaDR7NqPffdfqUlxm97DUoHSBhqTyJR2aNeOhPn1Ytn07nx86ZHU5LrNzWGpQ1kBDUvmi5GuvJTYigk8tvte7tuwalhqUNfhmfyF7DxRqSCqfEtGoEXvGj+fBSpci+QI7hqUGZQ2u6NyY5Y+30ZBUnuNwgKPE7c2GNbS2e7/6sFtYalC6oFlUkNUlKH91/AjMvw/2bLa6EtuxU1hqUCplpaYtSp/5vSPD6kpsyS5hqUGplJWCgqDbtaXfKPPPv4db2SMsNSjL5Bc42HfEumfbqADWPREcxbD7E4/PyiH2fNBdTawOSw1K/nsJ0CPzf+LMWYfV5ahA07o9xMTBjo89Opv//eIL+r38MiUO39zGrQzLgA/KitdJPjC8GaGNAn6VKG8zBnoMgMxv4ITnHgbWvEkTNh48yEtbt3psHp5mVVgGdCroxeTKNrr1BwzsWOexWQzr0oVr27Uj5aOPyC4o8Nh8PM2KsAzYoNSQVLYSdRHcmQx9fu2xWRhjmD94MMfz83n8Y88e5nuat8MyIIOysEg0JJX9dPolNPbstnh569bcfcUV/O3zz/mmwuN0fZE3wzIgg7JhiKHXpY01JFVAmn399TQJCeHhNWusLqXevBWWARmUACNvjNKQVAGpRVgY0/v354eTJzmen291OfXmjbAM2KBUKpA9cNVV7Bg7luZN/OPLgqfDUoNSqQAUEhREwyD/6sPAk2GpQamU8hueCku/Dsr8Agdzlx8nJ8/9XVgp5TEF+bBsBmx+3+pKfJInwtJvg7L8Osn3Np7mm/2FVpejlOsahULeSdj+kdWV+Cx3h6VfBmXli8mvuizU6pKUcl3FWxqPe+6WxsocZc+w9hfuDEu/C0q940b5Bectjd65g+ZoXh7XLFnCW19/7ZX5eYu7wtKvglJDUvmNqIvg4m6lQemFb3kxTZpwpriYie+/z5ki/+pu0B1h6VdB+Y8PTmlIKv/RIxFOHoXMPR6fVVCDBjyXlMT+nBye3bjR4/PztvqGpeVBaYyZZYzZYYzZZox53xjTpq5t/X5IFHMntNCQVP7h0qshuKHXTupcd/HF/PbSS3lywwYOnTrllXl6U33C0vKgBP4iIt1FpCewCphe14ZCgg09LmnstsKUslSjULjlAY/2KFTZX264gRKHgykffui1eXpTXcPS8qAUkYr/6woD/Oe0m1L11bUvxLT12uwujo5m4jXX8PqOHWw6eNBr8/WmimHpKmOHywGMMU8AI4Ec4DoROVbFNAOAAVV8PAnwxW5QEoB9FtdQWwlozd6SgFV1h4U1ZODAXqxZ8zkFBcW1+GQCvreue4vITTVN5JWgNMasBVpVMSpZRNIqTDcVaCwiM2rRdqqIpNa/Su/yxbq1Zu/xxbr9ueZgL9SCiAx0cdLlQDrgclAqpZSnWf4bpTHmkgovfwN8Y1UtSilVFa98o6zBU8aYzoAD2A/cZ3E9Sil1DsuDUkR+a3UNSilVHcsPvd0gw+oC6ijD6gLqIMPqAuogw+oC6ijD6gLqIMPqAuogw5WJbHF5kFJK2Zk/fKNUSimP0qBUSqka+EVQurNjDW8xxvzFGPNNWd3vGGOaWl2TK4wxw4wxu40xDmNML6vrqY4xZrAxZo8x5jtjzBSr63GFMeZlY0yWMWaX1bW4yhgTZ4z5yBjzVdm2McHqmmpijGlsjPncGLO9rOaZ1U7vD79RGmMiy+8ZN8Y8AHQREVtfZmSMGQT8R0SKjTFPA4iI6zefWsQYcymll3K9AEwUkc0Wl1QlY0wQ8C1wA3AQ+AIYISJfWVpYDYwx/YE84FUR6Wp1Pa4wxrQGWovIVmNMBLAFuMXO69oYY4AwEckzxoQAG4AJIvJZVdP7xTdKX+xYQ0TeF5Hy+2g/A7zX80E9iMjXIuL5DhLr70rgOxH5QUQKgX9QekODrYnIOuCE1XXUhogcEZGtZf/OBb4GYq2tqnpSqrz7oJCy4YK54RdBCaUdaxhjMoG7qEdXbRYZA7xndRF+JhbIrPD6IDbfef2BMSYBuBzYZHEpNTLGBBljtgFZwAcicsGafSYojTFrjTG7qhh+AyAiySISR+n94uOtrbZUTTWXTZMMFFNaty24UrdSlRljwoG3gAcrHeXZkoiUlPWD2xa40hhzwZ86LL8zx1W+2LFGTTUbY0YDNwO/Ehv9WFyLdW1nh4C4Cq/blr2nPKDsd763gOUi8rbV9dSGiGQbYz4CBgNVnkTzmW+U1fHFjjWMMYOBScBQEcm3uh4/9AVwiTHmYmNMQ+AO4N8W1+SXyk6MLAG+FpFnra7HFcaYmPIrTYwxoZSe9LtgbvjLWe+3gHM61hARW397MMZ8BzQCjpe99Zndz9QDGGNuBf4GxADZwDYRSbK0qAswxtwIPAcEAS+LyBPWVlQzY8zfKe2g+iLgJ2CGiCyxtKgaGGP6AeuBnZTugwDTRCTduqqqZ4zpDiyjdNtoAPxTRB6/4PT+EJRKKeVJfnHorZRSnqRBqZRSNdCgVEqpGmhQKqVUDTQolVKqBhqUSilVAw1KZVvGmH3GmDPGmDxjzFFjzNKy2+Qo+7dUvq3SGDOv7P3RZa9Hl72eV2m635S9v7TCey+WdcvmKP+8UqBBqezv1yISDvSktLOFqRXGfQuMLH9hjAkGhgPfV2rje2B42fhyo8o+X9F2YByw1S2VK7+hQal8gogcBdZQGpjlVgL9jDHRZa8HAzuAo5U+fpTSu0aSAIwxzYBrqHRLo4gsEJEPgQJ31698mwal8gnGmLbAEOC7Cm8XAGmU3scNpd8uX71AE6/y32+fd5R97qz7K1X+SINS2d27xphcSvuWzOL8XqFeBUaWdXCQCLx7gXbeAQYYY6KoPlCVOo8GpbK7W0QkgtKOIn5BaWcRTiKygdIOOpKBVSJypqpGyt5fDaQAzUXkE08WrfyLBqXyCSLyMbAUmFPF6NeBR6j5W+KrZdO97tbilN/zmY57laK0y7R9xpgeld7/K6XdfK2r4fMfU9rv4JdVjSzrt7IBYIAQY0xjoFBEHFVNrwKHfqNUPkNEjlH6rXB6pfdPiMiHNfUSX/ZAqQ9F5EIP73ofOEPpGfEXy/7dv/6VK1+n/VEqpVQN9BulUkrVQINSKaVqoEGplFI10KBUSqkaaFAqpVQNNCiVUqoGGpRKKVUDDUqllKqBBqVSStXg/wNZDMHvHK7pwgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5),\n",
    "                       constrained_layout=False)\n",
    "\n",
    "for axis in ['top', 'bottom', 'left', 'right']:\n",
    "    ax.spines[axis].set_linewidth(0.5)\n",
    "ax.tick_params(length=6, width=0.5)\n",
    "\n",
    "ax.set_xlim(-3., 3.)\n",
    "ax.set_xticks(np.arange(-3, 4))\n",
    "ax.set_xticklabels(f'−{-x:d}' if x < 0 else f'{x:d}'\n",
    "                   for x in np.arange(-3, 4))\n",
    "ax.set_xlabel('RMM1', fontsize='large')\n",
    "ax.set_ylim(-3., 3.)\n",
    "ax.set_yticks(np.arange(-3, 4))\n",
    "ax.set_yticklabels(f'−{-y:d}' if y < 0 else f'{y:d}'\n",
    "                   for y in np.arange(-3, 4))\n",
    "ax.set_ylabel('RMM2', fontsize='large')\n",
    "\n",
    "ax.set_aspect(1)\n",
    "\n",
    "ax.axhline(color='0.', linewidth=0.5)\n",
    "ax.axvline(color='0.', linewidth=0.5)\n",
    "ax.plot([0, 3], [0, 3], linewidth=0.5, color='0.')\n",
    "ax.plot([-3, 3], [3, -3], linewidth=0.5, color='0.')\n",
    "\n",
    "ax.plot([0, -3 / np.tan(np.pi * ip_rmm_angle/180)], [0, -3], linewidth=1.5,\n",
    "        linestyle=(0, (5, 5)), color='coral', clip_on=False, zorder = 4)\n",
    "ax.plot([0, -3 / np.tan(np.pi * pg_rmm_angle/180)], [0, -3], linewidth=1.5,\n",
    "        linestyle=(0, (5, 5)), color='royalblue', clip_on=False, zorder = 4)\n",
    "ax.plot([0, -3 / np.tan(np.pi * olr_rmm_angle/180)], [0, -3], linewidth=1.5,\n",
    "        linestyle=(0, (5, 5)), color='teal', clip_on=False, zorder = 4)\n",
    "\n",
    "for i in range(8):\n",
    "    ax.text(0.5 + 0.4 * np.cos(np.pi * (1 + i/4 + 1/8)),\n",
    "            0.5 + 0.4 * np.sin(np.pi * (1 + i/4 + 1/8)),\n",
    "            f'Phase {i + 1}', fontsize='large', ha='center', va='center',\n",
    "            transform=ax.transAxes, zorder = 5,\n",
    "            bbox=dict(edgecolor='1.', facecolor='1.', boxstyle='square, pad=0.'))\n",
    "\n",
    "ang = np.pi * 7/8\n",
    "d = 0.1\n",
    "coef = 1.5\n",
    "th = np.arcsin(coef * d / 2)\n",
    "P = Path.arc(-ang * 180/np.pi, ang * 180/np.pi)\n",
    "ax.add_patch(PathPatch(P, fill=False, linewidth=0.5))\n",
    "ax.arrow(np.cos(ang) + d * np.sin(ang - th),\n",
    "         np.sin(ang) - d * np.cos(ang - th),\n",
    "         -d * np.sin(ang - th),\n",
    "         d * np.cos(ang - th),\n",
    "         linewidth=0.5, color='0.', head_width=coef*d,\n",
    "         head_length=coef*d, length_includes_head=True)\n",
    "\n",
    "fig.savefig('figures/rmm_diagram.eps', bbox_inches='tight')\n",
    "fig.savefig('figures/rmm_diagram.png', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "306b70f0-69aa-4f0a-9b7a-7cceacd0b57e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqoAAALZCAYAAAB/MwMFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAADct0lEQVR4nOzdd3xUVfrH8c9DL6GLFKmioCJgbygExS5iVyyAupa1rL2jRMXe17Y2BFzUtbMo6NoCllV/6oKCCKiAQADpEHrg+f1xb+IkpMLcmWHyfb9e98XMuffc58wYT07OPcXcHRERERGRVFMl2QUQERERESmOGqoiIiIikpLUUBURERGRlKSGqoiIiIikJDVURURERCQlqaEqIiIiIilJDVURkZCZZZuZm1l2sssiIiJqqIqIiIhIilJDVURERERSUrkaqmZ2rZn9z8yWmNkGM1tsZmPMbJ+oCygikgxmNsDMppvZWjP7ysz2TnaZREQqGyvPFqpm9i7QC5gNrAV2BWoAK4Gd3X1BlIUUEUmEcGxqT2AVUBX4DegIVAOWENR3S5JWQBGRSqa8j/5vABq5+y7uvgewe5heDzguioKJiCRRLeAAd+8MnBCmNQYuS1qJREQqofI2VHcAPjSzZWa2CZgWc65l/IslIpJUP7r7RAB3fw9YGqZ3TV6RREQqn2plXWBm7YF/AzWBXOB7IA/YP7ykamSlExFJjpLGRFlCSyEiUsmVp0d1L4JGKsCR7r4PcHV0RRIRSbquZtYFwMyOBhqF6ROTVyQRkcqnzB5VYBKwkaDn9H0zmwlsH2WhRESSbC3wjZn9SjCZCmAZ8GTSSiQiUgmV2aPq7lOBAcAvBDP91wDHRlwuEZFk+ha4mGBSlQPfEDxRWpzUUomIVDLlWp5KRERERCTRtDOViIiIiKQkNVRFREREJCWpoSoiIiIiKUkNVRERERFJSWqoioiIiEhKUkNVRERERFKSGqoiIiIikpLUUBURERGRlJRyDVUzq2lmL5jZLDNbaWYTwr2288/XMbOnzGyRmS03s/Ex58zM7jOzxeFxn5lZcj6JiEh0yqorRUTSQbVkF6AY1YDZQE/gd+AY4DUz6+LuM4Fnw2t2BZYAe8TkvRA4AehGsO3hh8AM4B+JKbqISMKUVVeKiGzztoktVM3sB+B2YDLBntut3H1FMdd9CQxz92fD9+cDF7j7AYksr4hIMuTXle7+ZrLLIiISD6nYo1qImTUDOhI0UvcDZgG3m9k5wDwgK6ZS7gxMjMk+MUwr7r6ZQGbR9N69ew/u2rVrnEovIsnw8MMPv+Duf0l2ORKpSF0Zm56J6jqRtJXu9V1K96iaWXVgLPCru19kZjcDdxH0rt4NHAi8B+zr7lPMbCPQ2d1/DvPvDEwDqng5P2hWVpZnZWXF/8OISMKY2e3unpXsciRK0bqyPHlU14mkh3Sv71JuMlU+M6sCvASsBy4Lk9cAG4Ah7r7e3ccBnwJHhOdzgfoxt6kP5Ja3kSoisq0poa4UEUkLKdlQDWfqvwA0A0529w3hqR+KuTy2ETqZYCJVvm4UeQwmIpIuSqkrRUTSQko2VIGnCWb193H3NTHp4wlmt95kZtXMrDvQC/ggPD8CuNrMdjCzlsA1wLDEFVtEJKFKqitFRNJCyjVUzawtcBHBslPzzSw3PM4Kewv6EizDshx4DuifPyYVeAYYDfwITCIYv/pMPMr1xBNPsM8++1CzZk0GDhxY6Nzzzz/PTjvtREZGBkcddRQ5OTkF59atW8fFF19Ms2bNaNy4MX369GHu3LkF55csWcKJJ55I3bp1adu2LS+//HKZZXnvvfc4+OCDadiwIc2bN+cvf/kLK1euLBTzvPPOo379+jRv3pyHH3644NxXX33F4YcfTuPGjWnatCmnnnoq8+bNKzjv7txwww00adKEJk2acMMNN6CREyKpp7S6cmvuq7pORFJJyjVU3X2Wu5u713L3jJhjZHh+srsf6O513X03d387Jq+7+/Xu3jg8ro/X+NSWLVsyaNAgzjvvvELp2dnZ3HzzzYwaNYolS5bQvn17+vXrV3D+scce47///S8//PADOTk5NGrUiMsvv7zg/KWXXkqNGjVYsGABI0eO5K9//SuTJ5c+WmH58uUMGjSInJwcpkyZwty5c7nuuusKzmdlZTF9+nRmzZrFp59+yv3338/7778PwNKlS7nwwguZOXMms2bNol69epx77rkFeZ999lneeecdJk6cyA8//MDo0aN55pm4tPVFJI7Kqiu3lOo6EUkp7q4j5hg8eLCX5pZbbvEBAwYUvL/mmmv8kksuKXg/d+5cB/yXX35xd/eLL77Yr7vuuoLz7777rnfs2NHd3XNzc7169eo+derUgvNnn32233DDDaWWoag333zTd99994L3LVq08A8++KDg/aBBg/z0008vNu93333nGRkZBe8PPPBAf+aZZwreP//8877//vtXqDwiyUawbF3S65NUPlTXqa6T9JDu9V3K9ahui4Kfk8KvJ02aBMD555/PF198QU5ODqtXr2bkyJEcfXSwy+G0adOoVq0aHTt2LMjfrVu3MnsZiho/fjydOwfLxS5dupR58+bRrdufc8pKu2dsXoDJkyeXO6+IVC6q60Qk0VJ+wf9Ud9RRR3HGGWdw8cUXs/POO3PHHXdgZqxevRqAnXfemdatW7PDDjtQtWpVunTpwhNPPAFAbm4u9evXL3S/Bg0aFBqDVZYPP/yQ4cOH8/XXXxfcM/8+Zd3zhx9+4I477mDUqFEFabm5uZvlzc3Nxd0JJhiLSGWkuk5EkkE9qlupd+/e3H777Zx88sm0a9eOdu3aUa9ePVq1agUE47LWrVvH4sWLWbVqFSeddFJBL0NGRgYrVhTeCXbFihXUq1cPgA4dOpCRkVFwtGvXrtC1X331FWeeeSZvvPFGQU9FRkZGwX2Ku2e+X375haOPPprHHnuMQw45pCC9aJlWrFhBRkaGKm6RSk51nYgkgxqqcXDppZcyffp0FixYwMknn0xeXh677747ABMmTGDgwIE0btyYmjVrcvnll/PNN9+waNEiOnbsSF5eHtOnTy+418SJEwseT/3666/k5uYWHDNnziy47n//+x/HH388Q4cO5bDDDitIb9SoES1atGDixInF3hNg1qxZ9O7dm1tvvZVzzjmn0Gfp3LlzqXlFpPJSXSciiaaGajnl5eWxdu1aNm7cyMaNG1m7dm1B2qRJk3B3fv/9dy688EKuuOIKGjVqBMC+++7LiBEjWL58ORs2bOCpp56iZcuWbLfddtStW5eTTjqJ2267jVWrVvHFF18watSozSrUoiZNmsRRRx3F448/Tp8+fTY7379/f4YMGcLSpUv5+eefee655wqWmZk7dy6HHnool112GRdffHGxeR9++GHmzp1LTk4ODz300GZL1IhI+lJdJyIpJdmzuVLtKGkm7ODBg51gF6yCY/Dgwb506VLv0qWL16lTx5s1a+Y33nij5+XlFeRbtGiRn3nmmd60aVNv0KCBd+/e3b/++uuC84sXL/a+fft6nTp1vHXr1j5y5Mhi48caOHCgm5nXrVu34Nhtt90Kzq9du9bPPfdcr1evnm+//fb+0EMPFZzLyspyoFDeunXrFpzftGmTX3fddd6oUSNv1KiRX3fddb5p06YyyySSSkjzWbDxOFTXqa6T9JDu9Z0Fn1HyZWVleVZWVrKLISJbwcxud/esZJcjlamuE0kP6V7f6dG/iIiIiKQkNVRFREREJCWpoSoiIiIiKUkNVRERERFJSWqoioiIiEhKimtD1czuN7Pd4nlPEREREamc4t2jujPwvZl9Z2aXm9l2cb6/iIiIiFQScW2ouvuJQAvgBaAfMMfMRpnZSWZWPZ6xRERERCS9xX2Mqrsvdfen3P0goAswEXgAmGdmT8Q7noiIiIikp0gnU7n7dGAIcBOwEth8w2URERERkWJE1lA1s4PM7B/APOAu4EVgp6jiiYiIiEh6qRbPm5lZe+AcoD/QBHgD6Ovun8czjoiIiIikv7g2VIFpwEfArcDb7r42zvcXERERkUoi3g3Vtu6eY2Yti2ukmtke7j4hzjFFREREJA3Fe3mqnPDlf8yscew5M9sXGBvPeCIiIiKSvqKaTPUsQWM1A4KJVcC/gfMjiiciIiIiaSbej/4BcPe/m1lDYIyZ3QMMA85y94+iiCfJkZmZCUB2dnZSyyEiIiLpKZKGKoC732FmDYB/Ace5+/ioYomIiIhI+olbQ9XMZgNeJLlKePzTzABw9zbxiikiIiIi6SuePapnx/FeIiIiIlLJxa2h6u7j4nUvEREREZHIxqhKGsg6sfTzMyeV77qst+NTHhEREalUolqeSkRERERkq6ihKiIiIiIpSY/+RUREKoHv786s0PV73ZwdSTlSlV1gFbrenyu60JFEIWENVTMbCnwODHf3jYmKKyIiItsONRglViJ7VA04E7gG6JzAuBKR7IGHJLsIIiIiksYS1lB193MBzKx6omKKiIikEj1+F6mYhE+mcvcNZV1jZjXN7AUzm2VmK81sgpkdHZ5rZ2ZuZrkxx61F8g41sxVmNt/Mro7y84iIJIOZXWZm35rZOjMbluzyiIhEIeotVDdTzi1UqwGzgZ7A78AxwGtm1iXmmobunldM3ixgZ6At0Bz41Mx+cvf3yxFXRKRMZtYceBE4CJgMXOvuX8acX+Hu9SMuRg4wBDgSqB1xLBGRpEjJLVTdfRVBgzPfu2Y2A9gb+K6M7AOAge6+FFhqZs8BAwE1VEUkXh4jaChmEvxBPdrMLnf3l8PzFZsNsgXc/S0AM9sHaBV1PBGRZNgmtlA1s2ZAR4Kei3yzzMyBD4Hr3H2RmTUCWgATY66bCJwQVdlEpFLKBNq6+1rgf2b2CTDGzOq6+3OU4+mSiIiULbLJVGa2B3AIsB0xvQvuflsF71MdGEmwrNXPZpYB7AtMAJoAT4bnjwQywmzLY26xHKhXzH0zCX7ZFHL66aeTk5NTkSKmr4494nMffZ+SfqoRU3+6+w9hnfKhmW1W3yST6rrU4h2OqND18fxvlMzYFXHhPhdW6Pp4lTNZcaV0kTRUzexC4BHgP8DRwFjgCGBUBe9TBXgJWA9cBuDuucC34SULzOwyYF74yyE3TK8PrI15vbLovd09G8gump6VlTW4ZcuWFSlm+po2Pj73OfOq+NxHJHV8R/DH8Zv5Ce7+S9go/Biom5xibU51XWqZ/+t/KnR9y9NvTovYFfHst89W6PpnLnpmm44rpYtq1v/1wFHufiKwJvz3FKDMGf/5zMyAF4BmwMmlrBaQ/4itSjgudR7QLeZ8NwoPGZBtXGZmJpmZmckuhlRutwKriia6+yyCMau3J7xEIiJpKKpH/9u7+2fh601mVsXdx5rZyArc42lgV6C3u6/JTzSz/YFlwHSgEfB3INvd8x/3jwAGmdm3BI3cC4Bzt+rTiIjEcPevSzk3F7gj6jKYWf7wg6pAVTOrBeSVsBqKiMg2KaqG6hwza+fuM4FpQF8zW0TwCL9MZtYWuAhYB8wPOlchTNsE3A1sD6wgmEzVLyb7YIJG7ixgDXCflqYSkXgyszLH2rt71I3VQQT1Xb6zCXpysyKOKyKSMFE1VO8n6A2dSdCz8AZQA/hbeTKHj89KW97llVLyrgPOC4+0l/8IPDs7O6nlEKlksoCpwP9RfF0V+ax/d89CjVIRSXORNFTdfVjM67HhslE1wolQIiLbuquA/gRrO48A/hk+8hcRkTiK585U5u4evi46SSsPyAvHqm6KV0xJU1knln5+5qTyXZf1dnzKI1KEuz8GPGZmuxFsMvKlmU0DhgOvh092RERkK8Vz1n/s2qV5BDP8Y4/8NBGRtODuP7n7DUAH4H/AMKB7UgslIpJG4vnov3PM6/ZxvG/lpt5FkZRlZrsS9KieAfwKnA98mdRCiYikkXhuoTo75u2p7v5g0WvM7Grg4XjFFBFJBjO7nGCMah2CTUkOKVIHiohIHEQ16/82YLOGKsFyKmqoisi27jGCWf/fArsBd8UsoweAu/dPQrlERNJKXBuqZnZo+LKqmfWi8LItO1LMVqYiFZU98JBkF0HkDhKwBJWISGUX7x7VF8J/awFDY9IdmA9cHud4IiIJF65hKiIiEYtrQ9Xd2wOY2Qg99koM9S6KJJ6Z9SjrGncfn4iyiIiks6gW/C9opBZdU1XrqIpIGsgG/iDYFrqknanaJLJAIiLpKJ7rqBYws73M7L9mtgqtoyoi6WcUQWN0NHCKu7cucqiRKiISB1HN+h9OUIGfB6yOKIaISFK4+4lm1hjoB/zdzBoQLFM1QstUiYjET1QN1bbALflbqoqIpBt3XwI8CTxpZp2Ay4DfzKynu2vRfxGROIjk0T/wNnBERPcWEUkJFjgSGAycCvwT+C25pRIRSR9R9ajWAt42s88JlqUqoNUARGRbZ2ZdCHamOh34CRgBnO/ua5JaMBGRNBNVQ/Wn8BARSUcTCXamehbIIfjjvF/s7lTuPrT4rCIiUl5RLU91exT3FRFJEeMJZv0fVsJ5p/CmJyIisgWi6lHFzA4HzgC2d/c+ZrYPUN/dP4kqpohIIrh7ZrLLICJSGUTSUDWzy4ErgOeBU8LkNcDfgYOiiCkikgrMbHvgOne/LtllkZJ9f3dmha7f6+bsSMqRyuyC4vayKJk/V/xCPxX9riujeH3X6SiqWf9XAr3d/V4gfyeqn4FOEcUTEUmYcLb/+Wb2uJldYWbVzayhmT1CMOt/t2SXUUQkHUT16L8ekL/odX6zvzrBdoMiItu6Bwlm/H8BnAzsB+wPfAUc4O6Tklg2EZG0EVVDdTxwI3BXTNrfgE8jiicikkinAT3c/Tcz24VglZPT3P2NJJdLRCStRNVQvRwYbWYXAPXMbCqwEjguongiIonUwN1/A3D3n81stRqpIiLxF9XyVPPMbF9gX4LtVGcD37j7ptJziohsE8zM2gP5MyDyirwnvyErIiJbLqpZ/6PcvS/wTXjkp7/l7idFEVNEJIHqAr8Q0zAFfo157UDVhJZIRCQNRfXov1cJ6ZkRxRMRSRh3j2rFFBERiRHXhqqZ3RG+rBHzOt+OwKx4xhMRERGR9BXvHtXW4b9VYl5D8BhsNpAV53giIiIikqbi2lB193MBzOxLd38unvcWERERkcolqln/z5nZzkA/YAdgLvCKu0+PIp6IiIiIpJ9IJgSYWR/gO2AXYAnB1qnfmtnxUcQTERERkfQT1az/u4G+7l6wE5WZZQJPAP+OKKaISOTMbDZ/bg1dIndvk4DiiIiktagaqq2Az4qkfR6mi4hsy86Oeb0vMAD4O8GqJm2By4ARSSiXiEjaiaqhOgG4BrgvJu3qMF1EZJvl7uPyX5vZk8CR7j43Jm0s8D7wUBKKJyKSVqJqqP4VGG1mVxAsS9UGWAX0iSieiEgytARyi6TlEkwiFRGRrRTVrP+fzWxX4ACCijwH+NrdN0QRT0QkSf4N/NvMhgBzCNaPvgmNxRcRiYsotwH0IsemCGOJiCTDxcB/gX8A34f/fh2mi4jIVopqeaquwHTgNeA64HVgupl1K2f+mmb2gpnNMrOVZjbBzI4u5rrbzMzNrHeRvEPNbIWZzTezq+P0sURECnH3te5+o7t3cPfa7r5j+H5NIuKbWWMze9vMVoX15ZmJiCsikihRjVEdCjwJPOzubmYGXBWm713Ocs0GegK/A8cAr5lZF3efCWBmHYBTgXlF8mYBOxPMvm0OfGpmP7n7+1v7oUREijKzTkA3ICM23d2HJiD8k8B6oBmwB/CemU1098kJiC0iErmoGqodgUfd3QHCxupjBI3IMrn7qiLXvmtmMwgauTPDtCeBG4CnimQfAAx096XAUjN7DhhIMAtXZItlZmYCkJ2dndRySOows5uB24CJwOqYU07wh3mUsesCJwO7u3su8LmZ/Rs4B7gxytgiIokSVUN1DHA88HZMWh/gvS25mZk1I2j8Tg7fnwqsc/cxQWdtwXWNgBYEvzTyTQROKOaemUBm0fTTTz+dnJycLSlmNDr2iM99tuQzJSt2Mj9zKdavXx/eNoV+PiTZrgT2c/cfkhC7I5Dn7tNi0iYSPIkqUJ66bt5nwyoU+Nllsyp0PcDgPoOLTY86dklxmw98uUL3ueiZi7aJuPGMPbfsSwopqW5MVtyK/mzNvb1ikUv7XZDM2OkmqoZqVeBVM/uO4BF+a4Le0FFmVrAQtrv3L+tGZlYdGAkMD1cTqEew89XhxVye/+hteUzacqBe0QvdPRvILpqelZU1uGXLlmUVK3GmjY/Pfc68atuJnczPXIoaNWoAkFI/H5Jsa4CfkxQ7A1hRJG2z+q48dd38X/9TocDPzhhX9kVFPHPRM8WmRx27pLgV9ey3z24TceMZe1tX0Z+tlqffnBax001UDdVJ4ZHvJ+CDit7EzKoALxGMwbosTM4CXsofq1pE/nqG9YG1Ma9XVjR2RemxsEildCvwuJllAQtiT7h71Cud5BLUb7ESUt+JiCRKVOuo3r619wgnYL1AMEngmJg1WA8DWpnZJeH7pgQTre5z9/vMbB7BxIYPw/PdCIcMiIjE2bDw37/EpBnBGNWqEceeBlQzs53dfXqYpvpORNJKVD2q8fA0sCvQu8hSL4cB1WPe/x/B9qxjw/cjgEFm9i1BI/cC4NzoiyvbvKwTSz8/c1L5rst6u/Tzkk7aJyuwu68ys7eAO8zsLwSz/vsCByWrTCIi8ZaSDVUzawtcBKwD5sdMmLrI3UcWuXYjsDSc9QowmKCRO4tg/Nh9WppKRKLg7hWfVRRflxCsLvAHsBj4q5amEpF0kpIN1bDytzIvDK5tV+T9OuC88BARiZSZHU8w0347Yuqt8kwW3VruvoRiVjUREUkXKdlQTUl6LCwiRZjZYILtUl8l2IDkGeBM4F/JLJeISLqIagvVv5vZQUXSDjKzR6OIJyKSJOcBh7v7VcD68N8+QLuklkpEJE1E1aPaD7i2SNp3wDsEC2SLbHOyBx6S7CJI6mno7vlL8a03s+ru/o2Z9Sw1l4iIlEtUDVVn897aqsWkiYhsy341s87hBKZJwF/NbCmwNMnlEhFJC1E1VD8DhpjZ9e6+KVy4PytMT0vqbROplAYBTcLXNwIvE+wYdUmJOUREpNyiaqheAbwLzDOzWUAbYB7B2C0RkbTg7mNiXn8D7JTE4oiIpJ2odqaaY2Z7AfsDrYDZwDcJ2FJQRERERNJEZMtThY3S/0Z1fxERERFJb3FrqJrZFHffNXw9m2BC1WbcvU28YoqIiIhI+opnj+oFMa/PjuN9RURERKQSiltD1d0/BzCzqgSLYF8YbmcqIlJpmJkBh7j7+GSXRURkWxf3dU3dfSNwBKCJUyJSGdUAPk12IURE0kFUk6keAW43s8HuviGiGCIiSWFm/Us5XSNhBRERSXNRNVQvB5oDV5vZQmImVmkylYikgRcJtoUubniTJbgsIiJpK6qGqiZTicgWy8zMBCA7Ozup5SjFdOAGd9/sEb+Z1QJWJ75Iko78uWIX0BGpNKJa8H9cFPcVEUkR44BdKH4s6sbwvJTDXjdnVyzDBeqwFqlMImmomllN4DagH9DE3RuY2RFAR3d/IoqYIiKJ4u4XlXJuA9ArgcVJOPXyybagwn8ESUqK+6z/0CPA7sBZ/Dk+dTLw14jiiYiIiEiaiWqM6onATu6+ysw2Abj7XDPbIaJ4IiIJY2ZDy7rG3c9LRFlERNJZVA3V9UXvbWZNgcURxRORbUnWiaWfnzmpfNdlvR2f8lTcQGAq8G+C+k5ERCIQVUP1dWC4mV0FYGYtgEeBVyOKJyKSSCcC/cPjHWCEu/83qSUSEUlDUY1RvRmYAfwINCRYyiUHuCOieCIiCePuo9z9ZGA34AfgYTObamaDzKxhcksnIpI+Immouvt6d7/K3TOAZkC98H1xi2OLiGyT3H2puz8NHAmMArKAPZNaKBGRNBLVo3/MrAHQCcgI3wPg7p9EFVNEJFHMrApBA3UA0BN4DzhM60iLiMRPVOuoDgSeBHIpvEOLAztGEVNE0kf2wEOSXYRSmdlDwGkEw5tGAAPdfW1ySyUikn6i6lG9CzjF3cdGdH8RkWS6CvgVqAdcAlyS/9Qon7v3SEK5RETSSlQN1WrAfyK6t4hIsp2b7AKIiFQGUTVU7wMGmdmd7r4pohgiIknh7sOTXQYRkcogbg1VM5vNn9ulGtAcuN7MCi3y7+5t4hVTRCQZzKx/Wde4+4hElEVEJJ3Fs0f17DjeS0QklQ0DfgHmE/xhXpQTTLISEZGtELeGqpZkEZFK5DHgVGAlQYP0Ha0TLSISf5Es+G9mNczsDjObbmarwn/vNLNaUcQTEUkkd78KaAs8BZwEzDSz58zs4OSWTEQkvUS1herTwKHA34B9w38zCSp1EZFtnrtvdPf33P10YBdgKZBtZr2SXDQRkbQR1az/E4AO7r4sfP+TmX1NMKbrvIhiiogkVLgD3xkEu1M1Be4EJiSzTCIi6SSqhup8oA6wLCatNjAvongiIgljZn2A/sDBwCjgOnf/IrmlEhFJP1E1VF8C3jezx4E5QGvgUmCEmR2af5G7fxJRfJG0kpmZCUB2dnZSyyEFRgFTgX8Ca4AjzezI2Avc/bZkFExEJJ1E1VC9KPz35iLpF4cHBMu37Fg0o5nVJBjL2htoTLBN4U3uPtbMdiOYYdshvPw74G/u/lOY14B7gb+E558HbnR3R0QkfkYQ1GHbJSO4mV0GDAS6AK+4+8BklENEJGqRNFTdvf1WZK8GzAZ6Ar8DxwCvmVkXIAc4BZhFMBHsUuBVoGuY90KC8bHdCH6JfAjMAP6xFeURESkkBRqGOcAQ4EiCYVUiImkpqln/hZhZLzPrUZ5r3X2Vu2e5+0x33+Tu7xI0Nvd292VhuhMssr0R2Ckm+wDgIXef4+5zgYcIeh1ERBLCzLqY2etRxnD3t9z9HWBxWdeKiGzLolpHdZyZdQ9f30DQ6/mKmRUdClCeezUDOgKTY9KWAWuBx4G7Yy7vDEyMeT8xTBMRiRszqxOuDT3azB42s/pmtqOZvQ38F/gj2WUUEUkHUY1R3R34Knx9AdCLYAeXLyjcsCyVmVUHRgLD3f3n/HR3b2hmdQl6UGfFZMkAlse8Xw5kmJkVHadqZpkEa7sWcvrpp5OTk7N5YTqWq0O4bMXduzTJipvM2JUtLkD2v0o9vf6POcGtX36k9Ptknl7x2MmQzO86Pp4E9gQ+AI4mGCu6CzAcuMDdFyWrYEWVp67zDkdU6J7F1pEJcuE+F1bo+niVNVlxZdu0Lf0/leqiaqhWAdzMOgAWM9mpUXlvYGZVCFYPWA9cVvS8u68ys38AC81sV3f/A8gF6sdcVh/ILW4ylbtnA9lF07Oysga3bNly8wJNG1/eopfuzKsqdn2y4iYzdmWLW47YNVYvA6BlWWXcktjJkMzvOj6OBPZw9z/C1U1+B3q6+2dbe2MzyyYYo1+cL9y9Qrtflaeum//rfypUxpanV/jhWNw8++2zFbr+mYue2abjyrZpW/p/KtVF1VD9HHgCaAG8DRA2WsvVyxDO3n8BaAYc4+4bSri0CsF6rTsQPGqbTDCR6pvwfDdihgyIiMRJRvjHMe4+x8xy49FIDe+XGY/7iIikg6gaqgOBa4CFwANh2i7AY+XM/zSwK9Db3dfkJ5rZ4QSN3R+AugSzXpcCU8JLRgBXm9kYgln/1xCMYxWRLaD1W0tULdwq1fITir6Pcp1oM6tGUH9XBaqaWS0gz93zooopIpIMUS1PtZgia6i6+3vlyWtmbQnWYV0HzA86VyFMW0/Q8GxFsMj2N8BR7r42vOYZgrVZfwzfPx+miYjE0x/A0Jj3i4u8L3ad6DgaBAyOeX82cDuQFWFMEZGEi6ShGk6CGgScA7QkWPPvJeAud19fWl53n0VMr0QxSlz2JRyLen14iKSN7IGHJLsIEsPd2yU5fhaVtFHqz2n/FpHKJKpH//cD+xHsQjULaAvcSjC5aRuZ7SEiUjnsdXN2souQ8tRAFkmOqBqqpwLdwiEAAFPN7HuCdU3VUBUREZG0pT/+4ieqhmpJj+5Le6QvIomWdWLp52dOKt91WW/HpzwiIiIxotpC9XVgtJkdaWa7mtlRwDvAaxHFExEREZE0E1WP6vUEk6me5M/JVK8QLCclIiIiIlKmqJanWg/cFh4iIiIiIhUWVY8qZtaJYGeojNh0dx9afA4RERERkT9FtY7qzQS9qROB1TGnnMKLYotICtP6rSIikkxR9aheCezn7j9EdH8RERERSXNRzfpfA/wc0b1FREREpBKIqqF6K/C4mbUwsyqxR0TxRERERCTNRPXof1j4719i0oxgjGrViGKKiIiISBqJqqHaPqL7ioiIpLwvvviCMWPGUK1aNcy0KaNsGXcnLy+PY445hu7duye7OEkR1TqqswDCR/3NgAXuvimKWCIiIqnkxRdfxMwYMmSIGqmy1dyd4cOHM23aNM4999xkFyfhIhkzamb1zWwEsBaYC6wxs+Fm1iCKeCIiIqlixowZDBw4UI1UiQszY+DAgcyYMSPZRUmKqCY3/R2oC+wO1Aa6AHXCdJGtsnLdBlau21DpYovItqFKFc0blvirrD9XUY1RPQrY0d3zF/ufZmbnAr9GFC+h3vhpLnu3aEj7RnVZumY9l4/9gbG/LACgT8fmPHZUVxrUqh5J7JnLVtGuYd2C9/+aNIc3puTgDifs0oKzu7aOJC7Af379g2ETfmfywhWsXJdHvZrV6Ny0Pufu0YbDO2wfSczHv/6VIzpsT6ft6jF7+Wr6vfktX85eghkc3KYJL5+0DzvUr512sd/5OYdflqzizC6taFCzOrd9OoXflq2md/umXLrfjpHELE3epk1cOHoCQ/vuFcn93b1Q79N/Zy/hrSk5ONC3U3MOabtdJHFFRCS1RdU8Xws0LZK2HbAuongJdfUHk2gYNkSveP9H1uZt4sNzuvPhOd1Zk7eJK96Pbp+Drk9/WvD6H9/O4KoPJrFPi4bst0MjbvxoMk9+81skcR/57y8MeOc7dmpcl9t6dOLZPnswuOcu7Ny4LgPe+Z7Hvormb5C7P59Gq7AxeNnYH9i7RUMWXX8MC687hr2aN+SSMRMjiZvM2Pd+Po1Bn0xhxMTZ9Hjxcx788he2r1uTfVo05K7PpnHf59MiiVuajZuc4RN/j+z+De59r+D1W1NyOPKfX5K7Po/VGzZy3Ctf8frkuZHFFhGR1BVVj+rzwIdm9jAwC2gLXAU8G1G8hFqyZn1Bj+nHMxYy/fLe1KkefJUv9t2T9o99GFlsj3n95P/N4M3T9uPA1o0ByGy3HQPf+T6SHrcHvvyFTwcczC7b1SuUftKuLenXpRW9hn/OFQd0iHvc3PUbqVkt+HvqqzlLef3U/ahRNXh//+Gdaf7Q2LjHTHbsZ76byZfn9WCTO60f+YATdmlBt+bB8O7DO2xP/7e/44aDO8Y97qHDPy/x3Eb3Es/FQ+zd7/l8Gv86ZV+O3rkZEDwpuOGjyZzaeYdIyyASpe/vzozLffa6Obvc17Zr144FCxZQtWpVMjIyOOqoo3jiiSfIyMhg4MCBDB8+nHfeeYe+ffsW5Lnqqqt49NFHefHFFxk4cCDDhg3j3HPP5corr+SRRx4puG7UqFGccMIJDBgwgGHDhgFw4YUXMm7cOKZPn87QoUMZOHDgFn9OuyA+43v9ufLXXYn8vqZNm8Z1113Hl19+ycaNG9l33335+9//TqdOneLyudNJVD2qdwH3AqcAD4X/3h+mb/M6bZfBR78tBKBJ7RrMW/lnR/GC3HV4hL/UY//XnbdyLQe0alTwfr8dGjFnxZpI4q7akEfLerWKPdc8oyarN2yMJO5eLRrwxk85AOzUuC6T/lhRcG7KopXUrhbdsrzJir1kzXpa1KvFDvVrU7dGtYJGKgT/jXNWro0k7tdzl3Jo+6ac1aXVZke/3VtFEjNf7M/1jKWrOXKnP4eSHL5jU2YsXRVpfJF0NXr0aHJzc5kwYQL/+9//uOeeewrOdezYkREjRhS8z8vL47XXXqNDh8KdDh06dOC1114jLy+vIG348OF07Fj4D+Zu3brx1FNPsdde0QwRSoREfV/Lli3j+OOPZ+rUqSxYsID99tuvUANY/hTV8lQODA2PtHNf786c9dZ3XLx3O07YpQWHv/QF5+3ZFoAXJ8zi0n2jG0O4Nm8j/d/+Dgh6uRasWkfzjKABuWzt+oIev3g7edeW9HnlK27rsQtdm9WnQa3qrFi3gYnzVzDks6mcsmvLSOI+cmQXjh75X8b+soB9Wzak94gvOHGXFpgZb/+cw92H7hZJ3GTGblCzOms2bKR29aoMOqTwL4Io/xvv0bwBu2yXwSm7bd5zuS5vI5e8F90wi3UbN3Hbp1MK3i9fu4FGtWsAsHrDRs2eFtlKzZs358gjj2TChAkFaX369OGll15i6dKlNGrUiPfff5+uXbuycuXKzfJmZGTwwQcfcOyxx7JkyRK+/PJLzjnnHBYuXFhw3aWXXgpArVrFd2psS6L+vvbbbz/222+/gjxXXXUVQ4YMYfHixTRp0iQhn3FbEdffeGZ2aFlHPOMlS+8dt+fTAd35ZekqRv44m7kr1/Lgl9P5z69/MOiQTgzO3CWy2Lcc0okOjerSoVFdrty/A8vW/jkDffysxRwR0aSmfxy3Bwe1akz/d76j2YNjqX3Xv2n24FgGvPMdB+zQiKeP2yOSuHu1aMhPlxxGq3q1+fGPFTStW5Pv5i1j/cZNvHHqfly0T3R7SyQrdr/dWxX0jBd9xP/a5Bz23aFhJHGv3L8DjcPGYVHVq1bhxYgmUgH0230HZq9Yw+wVa+jTqTmzlq8uOPf+Lwvosn39yGKLVAZz5sxh7Nix7LTTTgVptWrVom/fvrz66qsAjBgxgv79+xebv3///gW9ia+++ip9+/alZs2a0Rc8SRL9fY0fP57mzZurkVqMePeovlDGeQcSP2U5Ars1rc/Ik/ZJeNzSGsHHd2rB8Z1aRBK3RtUq3NO7M/f07syytevJXb+RjBpVaVir+IZNPDWpU4O7Douu5zTVYt93eOcSz12wV1su2KttJHFLGwNaxYwBe7SJJC7AsBP2LvHcybvtwMnF9PKKSNlOOOEEzIzc3FwOPfRQbr/99kLn+/fvz3XXXUe/fv0YN24cw4cP58knn9zsPieeeCJXXXUVy5cvZ8SIETz00EOMHRvdHIFkScb3NWfOHC699FIefvjhSD7Tti6uParu3r6MIy0aqcVZtT4vIetr/m/esshjlGZt3kZmL1/D3BVrWLU+mnGpqWrqopU8/vWvPP71r0xbnBt5vGmLc3l7Sg4jJv7O21NymLY4FzOL/DF4cXGjtmTN+shjiFRG77zzDitXriQ7O5uff/6ZRYsWFTp/8MEHs3DhQu666y6OO+44atcufsm92rVrc+yxxxY8nk7X7TwT/X0tXLiQI444gksuuYR+/frF/fOkg6hm/ae1u8ZP5ZYewcy8xavXc9Zb3/KfX//ADHq1a8rLJ+/D9nWjeSSy97PZdGhcl/5dW9O/WxvaNqwTSZyictfncdmYifxr8lzWb/xzN9wd6tVmcM9OnL9Xu0jizs9dy7mjvufL2Uvo3LQ+Dx7RmYNa//lopP4977LipuMiiZ057DMeObILe7ZoyFtTchjwzvf0arcdZjDo0ym8dOLekfRg/758Nae/8X9MnL+CDo3r0qBmNVasy+PXpavo1qwBr56yD20axP+/e7LiAmz/wNhg1Yo92nDSri0KVtEQkfjo2bMnAwcO5Nprr+Wdd94pdO7ss8/mjjvu4NNPPy0+c6h///4ceuihDB48OMKSpoZEfF9Lly7liCOO4Pjjj+eWW26JV9HTjn4bbIH7vphe0FC97sNJ1KtZjXnXHAUE66pe/+GkUh9lbo1ggk0nhk/8nTvGT+WQNk0YuEcbTt61JXVrRPef88LRE1i9IY9P+h/MJnfu/3I6h7Zryq5N63HNfyaRu35jJMtTXfH+j7TMqEX2gIMZN2sRfV75iseP7sqZXYKNDaJcNGnighXsEc64v33cz7x35gH0CBee/+L3xVz07oRIGqrnjvqeQ9o04eP+3Qs12Fatz+OOcVMZ+M73fDLg4LSJC8HQkkPaNGFw9hQueW8ip+zWkoF7tCn4vkVk61155ZW0a9eOiRMLT4z829/+xiGHHEKPHj1Kzd+zZ08+/PBD9txzz2LPr1+/nk2bNuHubNiwgbVr11KjRo1tdkelKL+vFStWcOSRR9K9e3fuvffeuJY73WybPz1JFts4+ui3hTx1TDeaZdSiWUYtnjymGx/+trDEvFurajhW8JMBB/PL5YfTq912DBk/leYPvc/Ad74je2Y0scdMn8/LJ+/Dga0b071NE0acsDcP/fcXjuiwPf86ZR8e/TqaBf+zZy7iyWO7sWeLhlx5wE58OuBgrv9wMs99NxMovKxRvFWrYiwNJ6vlrFzLQeF6tQAHtGrMrOXRLAX29ZylDDl0t816FevWqMYdvXbh67lL0youBN/14Mxd+PVvR/DumQdQxYw+r3zFjo/9h6zsKfym5alEtlrTpk3p378/d9xxR6H0xo0bc9hhh5U5rMjMOOyww2jcuHGx54844ghq167Nl19+yYUXXkjt2rUZP3583MqfaFF+X2+//Tb/93//x4svvkhGRkbB8fvv0W2ssq1Sj+oWMIItHzd50GhtUufPCUWNa1dnxbq8EvPGU9uGdbi15y7c2nMXvvh9McMn/s5J//qGJTccG/dYDWpVZ9naDQWNmGVrN1C1SvA/6W5N67NwVTSbjuVt2kTepj//NOjarAHZAw/m8Je+ZOX6aL/nU3fbgcvG/MCzffbgnK6tueezaQzq0Qkn2D0qqpnorRvU5t1p8zmpmCW/xkxfQJsG0Wzbmqy4RfVoux092m7H40d35a0pOYz4YTZ3fzaN9bdqjUHZdlVkof54mTlz5mZpTz/9dJn5Pv/8z80/Bg4cWOLC/UOGDCn0Pjs7uyLFK1VFFuqPl0R+XwMGDGDAgAEVLmNlpIbqFshdn0e1O0fhDmYwYf5y9mrREIDpS1bRtE50M+FL+l+3e5smdG/ThMeP7hpJ3Av3akfvEV9y/p5tcOD572cxsFswC/y3pasKthqNt71bNOSDXxYUmvW9U+MMsgcczKEjPmfVhugaqw8dsTsXjP4frR5+nzYN6jDpjxXc9VmwfWmbBrUZdcYBkcR94uhunPza1zz831/o1qxBwZq1E+YvZ/LClbx52n5l32QbigvF/1zXrl6Vs7q25qyurclZGU3vtYiIpDY1VLfAjCuOKPS+ad0/G6bL1m7g7giXM3rmuG6lnq8Z0W5Jt/ToRKv6tXl3+nwArjtop4IJVI1qVWd0v2gabXf22rXg8Xustg3r8Nm5h/D897MiiQtBQ+mfJ+3D1EUr+WbuUuasWEvt6lXo1qwBh7RtQrWIxl0dtmNTfv3bEbw1JYfJC1fwx6p11KtRjYF7tOGEXVqwXZ1oJuolKy7AzWVsCduyXmJ6c0VEJLWooboF5qxYw46N6tKiXi3W5W1kyPipjPllAQDH7dy8YKJVFNo2qMO8lWuLjd2nY3NuPqRTZDsXDdijTcFamjOWruLJb37DgaN22p6dm2REEjNvk9OtWTChqaTPG5Uvfl/Mjo3q0mm7erRtWIch46fy+k8LMJsd6XedH/eCvduxNv8zT1/A+N8XMWfFmrSLC3BTkR24ZixdxZjpCwp+vnZqHM3Pl4iIpDZNptoCA0d9X/D62v9M4tOZi7jlkE4MOqQTn85cxA0fTk5K7E9mRBd71yc/Kng9buYiuv3jU96dPp8x0xew5zPZfDIjmklcJX3eWyL+vEVjX//hZLJnLmJQj+i/68oWF5L38yUiIqlNPapbIL9HE2DU1PlMuLhXwdaTme22o/NTn/DIUV3SKvacFWsLXg/6dApPHNOV/uEY1ZE/zGbQJ1P48vymcY9b2uftlcDv+p2f5yXsu65scSF5P18iUdi0aVPZF4lUUGX9uVKP6hZo06AO/xcu1VOzWpVCs9LzNjlr8qLbsSlZsWMX4Zi6aCX9dm9V8P6M3Vvx86KVkcStjN91ZYsLyfv5EolC+/btGTZsGO6Jn7ku6cfdefHFF2nfvn2yi5IU6lHdArf17MRpb/wft/XoxF/2bMtxL/+Xv+0fLHb/+De/cWbML9l0ib1h0yZe/N8snGBtuPUbN1E9HK+Yt2kTGyOqkCvjd13Z4kLyfr5EonDuuefyxRdfMGjQIKpVqxb5tseSvtydvLw8jj32WA466KBkFycp1FDdAmfs3oomtWuQNe5nvs1ZxoaNm+j/9ne0ql+bc/dow609o5vgk6zY++/QiBE/zAZgt6b1+GnhSvbdoREA42YtplNEk6kq43dd2eJC8n6+RKLSvXv3Evd3F5HyU0N1Cx3eYXsO77A9m9xZkLuO2tWr0LBWdOunJjt29sBDSjy3/w6NGHPWgZHFrmzfdWWMm8yfLxERSV0pN0bVzGqa2QtmNsvMVprZBDM7Ojx3gJl9aGZLzGyhmb1uZi1i8pqZ3Wdmi8PjPov4mUsVM1rUq5WwhlOqxI7VoFb1SNfYzFcZv+vKFrc4ifr52laUVkeKiKSblGuoEvTyzgZ6Ag2AQcBrZtYOaAQ8C7QD2gIrgRdj8l4InAB0A7oCfYCLElNsEZGEKK2OFBFJKyn36N/dVwFZMUnvmtkMYG93fzP2WjN7AhgXkzQAeMjd54TnHwIuAP4RaaFFRBKktDoSmJmMMomIRMVSffkMM2sGzAL2cPefi5y7EjjD3Q8I3y8HjnD3r8P3+wCfunu9Yu6bCWQWE3Jf4P+2sLjtSM4vimTFTWbsyhY3mbG3xbit3P0v8StK6iqtjgzPZ5I+dV0yY1e2uMmMXdnibm3s9K7v3D1lD6A68BHwTDHnugJLgENi0jYCu8S83xmCFW8SVN6sJH1PSYlbGT+zvuv0j7stHaXVken636ay/Tzqu07/uMmOnepHKo5RBcDMqgAvAeuBy4qc2wkYC1zh7p/FnMoF6se8rw/kevhTICKS6sws28y8hOPzmOtKrCNFRNJFyo1RhWD2PvAC0Aw4xt03xJxrS9CDcKe7v1Qk62SCiVTfhO+7hWkiItsEd88s65rS6kgRkXSSqj2qTwO7An3cfU1+opntAHwCPOHuxU2QGgFcbWY7mFlL4BpgWALKKyKSSMXWkSIi6SblelTDHtOLgHXA/JhlUC8CdgJ2BLLMLCv/hLvnb1vzTHj+x/D982GaiEhaKK2OdPeRSSuYiEgEUq6h6u6zgNIW6b+9lLwOXB8eyZBdyeImM3Zli5vM2JUtbkorRx2ZCNmVMHZli5vM2JUtbrJjp7SUX55KRERERCqnVB2jKiIiIiKVnBqqIiIiIpKS1FAVERERkZSkhmocmNllZvatma0zs2EJilnTzF4ws1lmttLMJpjZ0YmIHcb/p5nNM7MVZjbNzBK6fZuZ7Wxma83snwmKlx3Gyw2PqYmIGxP/DDObYmarzOxXMzsk4ni5RY6NZvZ4lDFjYrczszFmttTM5pvZE2aWchM/K6Nk1HVh3KTVd5WtrgtjJq2+S3RdF8ZUfZfC1FCNjxxgCDA0gTGrAbOBnkADYBDwmpm1S1D8e4B27l4fOB4YYmZ7Jyg2wJNs+T7lW+oyd88Ij06JCmpmhwP3AecC9YAewG9Rxoz5nBlAc2AN8HqUMWM8BfwBtAD2IPgZvyRBsaV0yajrILn1XWWs6yAJ9V0y6jpQfZfq1FCNA3d/y93fARYnMOYqd89y95nuvsnd3wVmAAmpQN19sruvy38bHh0SEdvMzgCWAR8nIl4KuB24w92/Cv9bz3X3uQmMfzJBRfpZWRfGSXvgNXdf6+7zgfeBzgmKLaVIRl0Xxk1afae6LqGSXdeB6ruUo4ZqmjCzZkBHErhlrJk9ZWargZ+BecCYBMSsD9wBXB11rGLcY2aLzOwLM8tMREAzqwrsAzQ1s1/MbE74aKh2IuKHBgAjPHFr2T0KnGFmdcLd6I4mqLxFgMTXd5WwroME13cpUteB6ruUo4ZqGjCz6sBIYLi7/5youO5+CcHjmUOAtwh2yonancAL7j4nAbFi3UCw69kOwLPAaDNLRK9KM6A6cArB97wHsCfBo8/Ihbsg9QSGJyJeaDxBj8IKYA7wLfBOAuNLCktGfVfJ6jpITn2X1LoOVN+lKjVUt3FmVgV4CVgPXJbo+O6+0d0/B1oBf40ylpntAfQGHokyTnHc/Wt3X+nu69x9OPAFcEwCQufv4/64u89z90XAwwmKDXAO8Lm7z0hEsPDn+X2CxkBdYDugEcG4NankklnfVZa6DpJW3yW7rgPVdylJDdVtmJkZ8ALBX6Inu/uGJBanGtGP28oE2gG/m9l84FrgZDP7PuK4xXESsI2luy8l+Cs79jFUIreT609iexcaA22AJ8JfkouBF0nsLytJQSlU31W2ug4SUN+lQF0Hqu9SkhqqcWBm1cysFlAVqGpmtRK0vMTTwK5AH3dfU9bF8WJm24dLiGSYWVUzOxLoR/QD/p8l+AWxR3j8A3gPODLKoGbW0MyOzP/vamZnEcxGTdQ4oheBy8PvvRFwFfBu1EHN7CCCR3+Jmv1K2IsyA/hr+F03JBgz9kOiyiAlS2JdB0mo7ypbXQdJr++SUteB6ruU5u46tvIAsvhzNmj+kRVxzLZhnLVAbsxxVgI+b1NgHMFs1BXAj8AFSfre/5mgz/t/wMrwM38FHJ7Az1mdYAmTZcB84O9ArQTEfQZ4KQn/XfcAsoGlwCLgNaBZosuho9j/Ngmv68K4SanvKltdF/OZk1LfJauuC2OrvkvRw8IvSkREREQkpejRv4iIiIikJDVURURERCQlqaEqIiIiIilJDVURkVKYWZaZuZlpQL+ISIKpoSoiIiIiKUkNVRERERFJSeVuqJrZQDP7r5nlmtkaM5sSLgQsIpIWzKyemQ0zs5VmtsjMHiJY21FERJKgXDuKmNkDBFu4QbAI8G9Ae2BfYGQ0RRMRSbhngTPC1wsJtlSsnbziiIhUbmUu+G9mbQi2+KoCfA4c6+4rwm30dnL3SdEXU0QkWma2I/ALwZ7m/3T3c8ysHvAdsDOAu0e637mIiBRWnkf/+8Vc95i7rwBw97VqpIpIGtmdoJEKwTaGuPtKErTXuIiIbE6TqUREREQkJZWnofoNsCl8fXn4KAwzq2FmnSMrmYhIYk0C8sdCnQJgZhnAsUkrkYhIJVdmQ9XdfwceDt/2AOaY2USCiQYXRFg2EZGEcfffgNfDt/3N7FeC8fltk1cqEZHKrVyP/t39OuBc4GugKtARyAG+ja5oIiIJdwHwEpALNCZouD6azAKJiFRmZc76FxERERFJBk2mEhEREZGUpIaqiIiIiKQkNVRFREREJCWpoSoiIiIiKUkNVRERERFJSWqoioiIiEhKUkNVRERERFJSSjZUzewyM/vWzNaZ2bAi5+qY2VNmtsjMlpvZ+JhzZmb3mdni8LjPzCzhH0BEpBzMrKaZvWBms8xspZlNMLOjw3NnmVluzLHazNzM9g7Pq74TkQopqX2VyvVNSjZUCXa9GgIMLebcswQ7xuwa/ntVzLkLgROAbkBXoA9wUZQFFRHZCtWA2UBPoAEwCHjNzNq5+0h3z8g/gEuA34Dvw7yq70SkooptX6VyfZOSDVV3f8vd3wEWx6ab2S7A8cCF7r7Q3Te6+3cxlwwAHnL3Oe4+F3gIGJigYouIVIi7r3L3LHef6e6b3P1dYAawdzGXDwBG+J/bCaq+E5EKKal9VYyUqW9SsqFaiv2AWcDt4aP/H83s5JjznYGJMe8nhmkiIinPzJoBHYHJRdLbAj2AETHJqu9EJO5Srb6ploggcdQK2B14E2gJHAi8Z2Y/ufsUIANYHnP9ciDDzCzmrwIAzCwTyCwaoHfv3oO7du0aSeFFJDEefvjhF9z9L8kuR0WYWXVgJDDc3X8ucro/8Jm7z4hJK1d9V1Jd17Vr18Ht27ePT+FFJGlGjRo1m+KHSma7e/YW3HKL65sobGsN1TXABmCIu+cB48zsU+AIYAqQC9SPub4+kFvclxj+x8sump6VlTU4Kysr7gUXkcR5+OGH5yS7DBVhZlWAl4D1wGXFXNIfuLtIWrnqu5LquszMzMHvvPPOFpdZRFKDmf3m7llxvOUW1zdR2NYe/f9QTFrslzSZYKBvvm4UeYQmIpJKwpmzLwDNgJPdfUOR890JniC9USSr6jsRiatUrG9SsqFqZtXMrBZQFahqZrXMrBowHvgduCm8pjvQC/ggzDoCuNrMdjCzlsA1wLDEfwIRkXJ7mmAVkz7uvqaY8wOAN919ZZF01XciUiGltK/ypVx9k5INVYIlWtYANwJnh68HhT0NfYFjCMZHPAf0jxnP9QwwGvgRmAS8F6aJiKSccNLCRcAewPyYNQzPCs/XAk4DhheTXfWdiFRUse0rSN36JiUbquFyLVbkyArPTXb3A929rrvv5u5vx+Rzd7/e3RuHx/XxGj+xbt06zj//fNq2bUu9evXYY489GDt2bMH5jz/+mF122YU6derQq1cvZs2aVSjveeedR/369WnevDkPP/xwoXuvXr2aSy65hO22244GDRrQo0ePMsvzwAMPsPvuu1OvXj3at2/PAw88UOj8zJkz6dWrF3Xq1GGXXXbho48+Kjg3fPhw9t57b+rXr0+rVq24/vrrycvLKzi/ZMkSTjzxROrWrUvbtm15+eWXK/x9iUjZ3H1WWL/Vil3D0N1HhufXuntDd/+4mLyR1Heq60TSVxntq4TXN+WxrU2mSpq8vDxat27NuHHjaNOmDWPGjOG0007jxx9/JCMjg5NOOonnn3+ePn36cOutt3L66afz1VdfAZCVlcX06dOZNWsW8+fPp1evXuy2224cddRRAFx44YXk5eUxZcoUGjduzIQJE8osj7szYsQIunbtyq+//soRRxxB69atOeOMMwDo168fBx54IGPGjGHMmDGccsopTJ8+naZNm7J69WoeffRR9t9/fxYuXMjxxx/Pgw8+yI033gjApZdeSo0aNViwYAETJkzg2GOPpVu3bnTurJVvRNKd6rpo6rovvviCMWPGUK1aNbSBmJSXu5OXl8cxxxxD9+7dk12c5HB3HTHH4MGDvby6dOnib7zxhj/zzDN+4IEHFqTn5uZ6rVq1fMqUKe7u3qJFC//ggw8Kzg8aNMhPP/10d3efMmWK16tXz5cvX17uuMW5/PLL/bLLLnN396lTp3qNGjV8xYoVBecPPvhgf/rpp4vN+9BDD/lxxx1XUPbq1av71KlTC86fffbZfsMNN2xV+UQSCcjyFKhPUvno2bNnub9P1XVbZ+jQof7iiy/6pk2b4n5vSX+bNm3yF1980YcOHVrseYJlqJJep0R1pOSj/23BggULmDZtGp07d2by5Ml06/bnZLi6devSoUMHJk+ezNKlS5k3b16h8926dWPy5GCy3DfffEPbtm0ZPHgw2223HV26dOHNN9+sUFncnc8++6ygF2Dy5MnsuOOO1KtXr9iYRY0fP74g77Rp06hWrRodO3YsV14RSW+q67bejBkzGDhwoHpSZYuYGQMHDmTGjBllX5yG1FDdAhs2bOCss85iwIAB7LLLLuTm5tKgQYNC1zRo0ICVK1eSm5tb8L7oOYA5c+YwadIkGjRoQE5ODk888QQDBgxgypQp5S5PVlYWmzZt4txzzwUotTxFDR06lG+//ZZrr722IG/9+vXLlVdE0pvquvioUkW/amXrVdafo8r5qbfCpk2bOOecc6hRowZPPPEEABkZGaxYsaLQdStWrKBevXpkZGQUvC96DqB27dpUr16dQYMGUaNGDXr27EmvXr34z3/+w7Rp08jIyCh0nHPOOYXiPPHEE4wYMYL33nuPmjVrllmeWO+88w433XQTY8eOZbvttqtQXhFJb6rrRCQVqKFaAe7O+eefz4IFC3jzzTepXr06AJ07d2bixD+3wF21ahW//vornTt3plGjRrRo0aLQ+YkTJxY8fipuu9b8x0MdO3YkNze30PHSSy8VXDd06FDuvfdePv74Y1q1alWQ3rlzZ3777bdCPQOxMQHef/99LrjgAkaPHk2XLl0K0jt27EheXh7Tp08vMa+IpDfVdZJo2dnZhf7bppOLL76YO++8s1zXZmZm8vzzz0dcom1MsgfJptpR2mSqiy66yPfff39fuXJlofQ//vjD69ev72+88YavWbPGr7/+et9///0Lzt9www3eo0cPX7JkiU+ZMsWbN2/uY8eOdXf39evXe4cOHfyOO+7wDRs2+Oeff+4ZGRkFkxNK8s9//tObNWvmP/30U7Hn999/f7/mmmt8zZo1/tZbb3mDBg38jz/+cHf3jz/+2Bs3buzjxo0rNu/pp5/uZ5xxhufm5vrnn3/u9evX90mTJpVaHpFUgiZTlXmUNplKdV18Fft7ZfAJ8TkqoG3btv7hhx/G50PF2aeffuo77LBDwuPOmDHDAd+wYUNc7vfiiy969+7dtzh/z549/bnnniv2XEntEzSZSgBmzZrFM888w4QJE2jevHnB46mRI0fStGlT3nzzTW655RYaNWrE119/zauvvlqQ9/bbb6dDhw60bduWnj17ct111xUs11K9enVGjRrFmDFjaNCgARdccAEjRoxgl112KbU8gwYNYvHixey7774FZbn44osLzr/66qt8++23NGrUiBtvvJE33niDpk2bAnDnnXeyfPlyjjnmmIK8Rx99dEHep556ijVr1rD99tvTr18/nn76afUyiFQSqutU14mklGS3lFPtqMjyVCKSmlCPaplHRZankq2Tij2qGzdu9DvvvNPbtGnjTZs29XPOOceXLVvm7n/2Mg4bNsxbt27tTZo08SFDhhTcZ+3atX7FFVd4ixYtvEWLFn7FFVf42rVr3d194cKFfuyxx3qDBg28UaNGfvDBB/vGjRsL4t99992+6667esOGDX3gwIG+Zs0ad/+zR/XBBx/0pk2bevPmzUtcjqk4gD/22GPevn17b9KkiV977bUFcUv7rK1bt3bA69at63Xr1vUvv/zS3d1feOEF32WXXbxhw4Z+xBFH+MyZMwvFevrpp32nnXbyBg0a+CWXXOKbNm3yn376yWvWrOlVqlTxunXreoMGDdzdfcCAAX7LLbe4u/uSJUv82GOP9e22284bNmzoxx57rM+ePbvg3upR3fxQj6qIiEglM2zYMIYNG8ann37Kb7/9Rm5uLpdddlmhaz7//HOmTp3Kxx9/zB133FGwQsNdd93FV199xYQJE5g4cSLffPMNQ4YMAeChhx6iVatWLFy4kAULFnD33XcXWpZr5MiRfPDBB/z6669MmzatIB/A/PnzWb58OXPnzuWFF17g0ksvZenSpeX+TG+//Tbffvst33//PaNGjWLo0KFlftbx48cDsGzZMnJzcznwwAMZNWoUd999N2+99RYLFy7kkEMOoV+/foVivfvuu/zf//0fP/zwA6+99hoffPABu+66K//4xz848MADyc3NZdmyZZuVMX/VilmzZvH7779Tu3btzb53KUwNVRERkUpm5MiRXH311ey4445kZGRwzz338OqrrxbaYnbw4MHUrl2bbt260a1bt4KJciNHjuS2225j++23p2nTpgwePLhg8lv16tWZN28es2bNonr16hxyyCGFGqqXXXYZrVu3pnHjxtxyyy288sorBeeqV6/ObbfdRvXq1QuGa0ydOrXcn+mGG26gcePGtGnThiuvvLLg3uX5rLH+8Y9/cNNNN7HrrrtSrVo1br75ZiZMmFBou+Abb7yRhg0b0qZNG3r16lWuXdYAmjRpwsknn0ydOnWoV68et9xyC+PGjSv3Z6yM1FAVERGpZHJycmjbtm3B+7Zt25KXl8eCBQsK0po3b17wuk6dOgVr5RaXNycnB4DrrruOnXbaiSOOOIIdd9yRe++9t1Dc1q1bF5sPgkZctWp/7uweG7M8Srp3eT5rrFmzZnHFFVfQsGFDGjZsSOPGjXF35s6dW3BNSd9NWVavXs1FF11E27ZtqV+/Pj169GDZsmVs3Lix3J+zslFDVUREpJJp2bJloR7C33//nWrVqtGsWbMtytuyZUsA6tWrx0MPPcRvv/3Gv//9bx5++GE+/vjjgmtnz55dbL54KOnepX3W4nYLa926Nc888wzLli0rONasWcNBBx1UZhnK2n3soYceYurUqXz99desWLGiYOhBMNRUiqOGqoiISCXTr18/HnnkEWbMmEFubi4333wzp59+eqEezdLyDhkyhIULF7Jo0SLuuOMOzj77bCAYu/nLL7/g7jRo0ICqVasW2lHpySefZM6cOSxZsoS77rqL008/PW6f6YEHHmDp0qXMnj2bxx57rODepX3Wpk2bUqVKFX777beC+1x88cXcc889BdvpLl++nNdff71cZWjWrBlz5sxh/fr1xZ5fuXIltWvXpmHDhixZsoTbb799Kz91+lNDVUSkgszsfjPbLdnlENlS5513Hueccw49evSgffv21KpVi8cff7xceQcNGsQ+++xD165d6dKlC3vttReDBg0CYPr06fTu3ZuMjAwOPPBALrnkEnr16lWQ98wzzywYFtChQ4eCfGW5++67Cy0tVpy+ffuy9957s8cee3Dsscdy/vnnl/lZ69Spwy233EL37t1p2LAhX331FSeeeCI33HADZ5xxBvXr12f33Xdn7Nix5SrnoYceSufOnWnevHnBLmixrrzyStasWcN2223HAQccULB8m5TM1N1cWFZWlmdlZSW7GCKyFczsdnfPivD+bwNHA5OBYcAr7r4oqnhRyMzM9Ozs7GQXo1LIyspCv1egXbt2PP/88/Tu3Tvu9zYzpk+fzk477RT3e6eKkn6OzGycu2cmvEAJoh5VEZEKcvcTgRbAC0A/YI6ZjTKzk8ysenJLJyKSPtRQ3cZlZmaSmZmZ7GKIVDruvtTdn3L3g4AuwETgAWCemT2R3NKJiKSHskdNi4hIqdx9upkNASYB9wEXA1rFWyTGzJkzI7u3hjGmr0gbqmZWH8iITXP3nBIuFxHZ5pjZQUB/4FRgCfAiMCKphRIRSRORNFTNrDfwLNAWiF1UzIGqUcQUEUkUM2sPnEPQQG0CvAH0dffPk1owEZEUY2a/UbgtuNkl7t6upJNR9ai+ANwJvAqsiSiGiEiyTAM+Am4F3nb3tUkuj4hIqvpLzGsHRgEnxLx/r7TMUTVUawEvurv2BBORdNTW3XPMrGVxjVQz28PdJyShXCIiKcXdP4l9b2Z5sWlmllda/qhm/T8CXG9l7SUmIrINihlr/x8zaxx7zsz2Bcq3OriISOWz3MxaA5hZC2BVaRdH1VB9E7ggLMxvsUdE8UREkuFZgsZqBhRMrPo3cH5SSyWyDRk5ciRHHHFEua7Nysoq2K41Gdq1a8dHH30EBLtl/eUvfykjhxTjXWCsmd0P/Ad4v7SLo3r0/wbwGfA6GqO6VQ695PdSz0+cvq5c133yVJu4lUlEAu7+dzNrCIwxs3sIdqk6y90/SmrBJOWVVWeXV0Xq9nbt2rFgwQKqVq1K3bp1Ofroo3niiSfIyMgoM1+8dpSaOXMm7du3Z8OGDVSrFjRBzjrrLM4666ytvnei3XzzzXG5T3HfSZq7FrgJ2AsYDQwp7eKovpH2wJ7uvimi+4uIpAR3v8PMGgD/Ao5z9/HJLpNsmfzNU9J5a9nRo0fTu3dv5s6dy5FHHsmQIUO49957k12spNi4cSNVq2ohokQws47AaUBLIAcYWd5trqN69D8KODSie4uIJI2ZzTaz32MP4HSC+vSfMWkiKWuHHXbg6KOPZtKkSQD8+9//pnPnzjRs2JDMzEymTJkCwDnnnMPvv/9Onz59yMjI4P777wfgq6++4qCDDqJhw4Z069atUOM+MzOTW2+9le7du1OvXj2OOOIIFi1aBECPHj0AaNiwIRkZGfz3v/9l2LBhHHzwwQX5r7jiClq3bk39+vXZe++9+eyzz8r9ue6//35atGhBy5Ytef755zEzfvnlFwAGDhzIX//6V4455hjq1q3Lp59+ynvvvceee+5J/fr1ad26NVlZWYXu99JLL9G2bVuaNGnCXXfdVehc0WEI8fxO0omZ9QW+A3YHVob//s/Mji9P/qgaqjWBf5vZB2Y2IvaIKJ6ISKKcTbCGauxxFnBskbQymVlNM3vBzGaZ2Uozm2BmR8ecr2NmT5nZIjNbbmbjY86Zmd1nZovD4z5NYJXymj17NmPGjGHPPfdk2rRp9OvXj0cffZSFCxdyzDHH0KdPH9avX89LL71EmzZtGD16NLm5uVx//fXMnTuXY489lkGDBrFkyRIefPBBTj75ZBYuXFhw/5dffpkXX3yRP/74g/Xr1/Pggw8CMH588CO8bNkycnNzOfDAAzcr27777suECRNYsmQJZ555Jqeeeipr15a9Atz777/Pww8/zEcffcQvv/xSbM/4yy+/zC233MLKlSs5+OCDqVu3LiNGjGDZsmW89957PP3007zzzjsA/PTTT/z1r3/lpZdeIicnh8WLFzNnzpxiY0f9ncSLmV1mZt+a2TozG1bkXFT1zd3A8e5+hrvf4O5nAMeH6WWKqqE6mWAbwS+BX4scIiLbLHcfV56jnLerBswGegINgEHAa2bWLjz/LNAY2DX896qYvBcSrEXYDegK9AEu2sqPJ2nuhBNOoGHDhhx88MH07NmTm2++mX/9618ce+yxHH744VSvXp1rr72WNWvW8OWXXxZ7j3/+858cc8wxHHPMMVSpUoXDDz+cffbZhzFjxhRcc+6559KxY0dq167NaaedxoQJE8pdxrPPPpsmTZpQrVo1rrnmGtatW8fUqVPLzPfaa69x7rnn0rlzZ+rUqbNZ7yhA37596d69O1WqVKFWrVpkZmbSpUsXqlSpQteuXenXrx/jxgX/+77xxhscd9xx9OjRg5o1a3LnnXdSpUrxzaaov5M4yiEYEzq0mHNR1TctCeYtxRoXppcprmNUwx2pxrn77fG8r4hIOnL3VUBWTNK7ZjYD2NvMahH0OrRy9xXh+e9irh0APOTucwDM7CGC1Vb+EXnBZZv1zjvvbDYpKicnh7Zt2xa8r1KlCq1bt2bu3LnF3mPWrFm8/vrrjB49uiBtw4YN9OrVq+B98+bNC17XqVOH3NzccpfxwQcf5IUXXiAnJwczY8WKFQWPyUuTk5PDPvvsU/C+devWm11TNO3rr7/mxhtvZNKkSaxfv55169Zx6qmnFtwv9vq6devSpEmTYmNH/Z3Ei7u/BWBm+wCt8tPNbBeiq28mAFcCD4Z5Dbg6TC9TvCdTXQu8YmZfEOw0MMbdi/9JFxGRQsysGdCR4KnUfsAs4HYzOweYB2S5+5vh5Z2BiTHZJ4ZpRe+ZCWQWTc/MzCxXL1VlcthhhwHE/XtZvXp1XO8Xby1btuTHH38seO/uzJ49mx122AGAok94W7duzTnnnMNzzz1X4VhlPS3+7LPPuP/++/n444/p3LkzVapUoVGjRrh7mfdu0aJFoUfzs2fPLjP+mWeeyWWXXcbYsWOpVasWV155ZUGjuEWLFgVjdSH477h48eJiY0f5ncTGL+Fns6GZZRWTnu3u2eUsRlzqmxJcAow2sysJniC1JlgR6rjyZI5rQ9XdjzKzOsBhwDHAIDNbBowhaLh+qZUAREQ2Z2bVgZHAcHf/2cxOIph08CbBI7IDgffM7Cd3nwJkAMtjbrEcyDAz85jf6uEvquyi8TIzMwd36tQpqo+T+rJO3Czp42HB08lbN35fzPVvb3GoOnXqbHHeRDjttNO49957+fjjj+nRowePPfYYNWvW5KCDDgKgWbNm/Pbbn8ugn3322ey777588MEH9O7dmw0bNvDVV1+x00470apVq5LCANC0aVOqVKnCb7/9RseOHTc7v3LlSqpVq0bTpk3Jy8vj3nvvZcWKFcXcqfjPcd5553HOOefQtm1b7rzzzjLzrFy5ksaNG1OrVi2++eYbXn755YI1XU855RT2339/Pv/8c/bbbz9uu+02Nm0qvgkT5XeSr06dOpTw/+yy8s6gL0Ur4lDfFMfdp4Q9tgeGceYStAdL3ZEqX9zHqLr7ancf7e5/dfe2BJMMlhKMiZhnZq+a2f7xjltZ7dn3X+zZ91/JLoaIbAUzqwK8BKwHLguT1wAbgCHuvj4c9/opkL8yei5QP+Y29YHcsn5piBTVqVMn/vnPf3L55Zez3XbbMXr0aEaPHk2NGjUAuOmmmxgyZAgNGzbkwQcfpHXr1owaNYq7776bpk2b0rp1ax544IESG3Gx6tSpwy233EL37t1p2LAhX331VaHzRx55JEcddRQdO3akbdu21KpVq9hH+MU5+uij+dvf/kavXr3YaaedOOCAAwCoWbNmiXmeeuopbrvtNurVq8cdd9zBaaedVnCuc+fOPPnkk5x55pm0aNGCRo0aldjojPI7SZBI6xt3z3P3z9z9FXcf7+555Z1gb/Gs08ysobsvK+V8A4IPvTRVF8TOysry4gZglyWq9feSsSi0yLbOzG6PQw/DlsYeCnxO0DO6sRzXG8HEhnbAMe6+Jkw/jGAr1jr5PQ9m9m/gY3d/zMy+BF509+fCc+cBF7r7AeUpZ2ZmpqfzeqFlKqZHNTPsUc0eeEgx1295j2pWVlaxE3skWlOmTGH33Xdn3bp1abGQfkk/R2Y2zt0zK3IvMxtCMB51YPg+svomXIbqRqA5EDvOoQ3wO4C7ty8pf7x7VOeb2etmdpyZbbaKrrsvd/fXU7WRKiISBwacCfxQzuufJphl2ye/kRoaT1CJ32Rm1cysO9AL+CA8PwK42sx2MLOWwDUEO2OJVFpvv/0269atY+nSpdxwww306dMnLRqp8RLWJbWAqkBVM6tlZtWItr55nOCP8b8QbC99PnAesDbmfYni3VDdH5gJPAPkmNkjZrZHnGOIiKQsdz/X3XsDe5R1rZm1JVjiZQ+CP/Rzw+Msd98A9CUY778ceA7o7+4/h9mfIdh+8EdgEsE8gGfi/HFEtinPPPMM22+/PR06dKBq1ao8/fTTyS5SqhlE8Jj/RoI1odcAgyKub+q5+/Pu/knM8SmwMf99aZnjPZlqIjDRzG4ADidY9PpzM/sNGE6wZdb8su5jZpcBA4EuwCsxXdPtgBnAqpjL73P3O8PzNQl6J04BVgP3u/vDcflwIiIVEFb8ZV0zi8KPwoqen0wwAaG4cw5cHx4iQrDov5QsHBKVVcK5qOqbISWk31GezJH0h4cz+z8APjCzusDJwO0EuxCUPKr5T/kL0h4J1C7mfMMSZotlATsDbQnGQnwazljTT66IxIWZzQbKHNzv7hoYLiKVnrs/bGatCHpwWxMsg/VPd3+wPPkjHbgRjmM4G+gPNAFeLk++khakLYcBwEB3XwosNbPnCHpm1VAVkXg5u+xLZFtU7CSqONBCDBIP2+rPkZntC4wiGDrQH3gN+MHMjnP3Mpc4iHtD1cxqE/SgDgB6EGybdR/wprvHa9XjWWbmwIfAde6+yMwaAS3YfEHaE+IUU0SECmyPKgLApk2bWLduXanLJImUZt26dWzcWOYiIqnqAYJOxP+Y2anuPsDMjgQeArqXlTneW6gOB04kWMz1JeA8d998a4gttwjYl2DbrSbAkwQLZB9JsBgtbL4gbb0SyppJMbu1nH766eTk5Gx2/fD3lm+WFuv3eUEb/J7nppR63YBjG5R6vqjD91xV9kXlUNxnEpGtF04YPQTYjpjxpu5+W7LKJKnlggsu4Prrr+f+++9XY1UqbN26dVx//fVcc801yS7KluoCFFrtyd0/MLNXy5M53j2qq4EjytOVuyXcPRf4Nny7IJx0Nc/M6hEsRgvBIrRrY16vLOFe2RSzW0tWVtbgli1bbnb9h/8rfQOFpblVw+vqlnrdTRdsfu/SlBW3vCoaV0TKZmYXAo8A/wGOJliH8AiCx1yyhaJalzpZ2rRpwzXXXMOdd95J1apVy71lpoi7s3HjRq655hratNlmh73nATUI2mYGBeu2zixP5njP+v9rPO9XnpDhv1XcfamZzQO6EQwJIHw9OcFlEpHK43rgKHf/zMyWuvuJZnY0cEayCyappU2bNgwZUtLkZ5G09g3B8qXjgGpm9j6wD8Ew0TKl5Cq44eKz1YhZkJagRb43sAyYDjQC/g5ku3v+c/kRwCAz+xZoBlwAnJvY0otIJbK9u38Wvt5kZlXcfayZjUxqqUREUsdfgfy9ZO8A5gNnufvi8mROyYYqwYK0g2Pen02wvNVUgiWutgdWEPSc9ou5bjDBOqqzCBaxvU9LU4lIhOaYWTt3nwlMA/qa2SJgfXKLJSKSGtx9TszrByqaPyUbqqUtSAu8Ukq+dQTbcp0X/1KJiGzmfoLtT2cS9BS8QTAW629JLJOISMows08oYWMTd+9VVv6ENVTNrA0w19232fUVSrNn338luwgikmDuPizm9dhwmbwa4cRPERGBfxZ535xgWGaJHY+xEtmjOhOYaWb3uPtzCYwrIhI3ZmbhdoKYWZUip/OAvHCs6qbNc4uIVC7uPrRompn9CxgOlLmMXyIbqu0Jts7qmcCYIiLxtpxg6TsIGqZFt4uxMK1qIgu1zco6cfO0mZNKPpf1drTlEZFEmAN0Lc+FkTRUzWx/d/86Ns3dZ5lZM3e/K4qYIiIJ0jnmdfuklUIKpNu6qyLpxMwGFEmqA5wGfFme/FH1qH7Inz0Osd4HGkcUU0QkckV22zvV3R8seo2ZXQ08nLhSiYikrPOLvM8F/gtsVncWJ95bqFYheOxlFmy9ETvLqwPBYzIRkXRxG8VXtoNQQ1VEBHfvsTX5492jGjteq2ijdBOgx/4iUqZUf5RrZoeGL6uaWS8K/1G+IyVs3SwiUhmZWW2gKYXryvEE85bc3WeVlDfeDdX2YSHGAbEtaAcWuvuaOMcTEUmGF8J/awGxM1qdYNeVyxNeIhGRFGRmWcAtQNFVUgB+JWg3FncOiHNDNaZF3Dae9xURSSXu3h7AzEa4e/9kl0dEJIVdAuxG0CjN5wQ7jDZk85VTColq1n9j4FpgDyAj9tzWjlUQEUkVsY3Uomuqah3VLZc98JASzx16ye+bpU2cvq7Ec5881SZ+BRORLVHN3acXTQyXpC5zE6ioZv2/DNQEXgNWRxRDRCSpzGwv4EmC9QBr5SejdVSlAlJ9TLbIVjqqhPTDy5M5qobqQUBTd18X0f1FRFLBcGA0cB76o1xEpDi7mtmuJZz7GsDMznT3l4u7IKqG6g9AKwqPRxARSTdtgVvyt1QVEZHNFF1HNZ8R/LEP0J/gafxmomqofgK8b2YvEsyALVDcnq8iUrkUN5YwVmljDmOlwPjDt4EjgA+SXRARkVRUnrlJ7l7S8IDIGqqHEOzjWnT8gVN4KRcRkW1ZLeBtM/uczf8o12oAIlLphRNNjwRy3f2ziuaPpKHq7r2iuK+ISIr5KTxERKR4zxE0VKub2QPu/qCZnQns5+5XlpU5qh5VzKwJcAzQ3N0fMLOWQBV3nxNVTBGRRHL325NdhlSRzJnre/b9V8JjbrGsEzdPmzmp5HNZb0dbHpHonQJ0BJoAbxBsOz0GuBe4sqzMJe4EsDXMrCcwFTiLYC9sgJ2Bp6OIJyKSLGZ2uJm9YGajw/f7xGyxKiJS2eUCS939J6BFmLYCaFCezJE0VIFHgdPDwbF5YdrXwH4RxRMRSTgzu5zgD/Dp/Llt9BpgSDnz1wwbubPMbKWZTTCzo8Nz7czMzSw35ri1SN6hZrbCzOab2dVx/ngikmbM7DIz+9bM1pnZsJj0KOubR4Abi6RdAUwrT+aoHv23c/ePw9f5y7asjzCeiEgyXAkc5u4zzeyGMO1noFM581cDZgM9gd8Jhku9ZmZdYq5p6O55xeTNInhS1RZoDnxqZj+5+/sV/hQiUlnkEPwhfSRQu5jzUdQ3lwCtzexKoJ6ZLQHWAaeWp8BRNRx/MrMj3T12yZbewI8RxRORNLINjTmsR9DQhD//KK9O8Id5mdx9FcEvgHzvmtkMYG/guzKyDwAGuvtSYKmZPQcMBNRQFZFiuftbEAxRIljvvry2pr75S2wRgEXAVHcvVz0ZVUP1GoIK9z2gtpk9A/QB+kYUT0QkGcYTPNK6Kybtb8CnW3IzM2tGMOlgckzyLDNz4EPgOndfZGaNCMZ6TYy5biJwwpbEFREJxb2+cfdPtqZAUS1P9ZWZdSOYTDWUoMdhP834F5E0czkw2swuIHikNRVYCRxX0RuZWXVgJDDc3X82swxgX2ACwWzZJ8PzRwIZYbblMbdYTtDDW/S+mUBm0fTMzEymTp1a0WIC8Nb4TZul/T4v2EH2nuemFJvnpB4lTIk48JwKxT7891UVur7Ez1jBuGzhd1WeuM+EacVGiFdcSWcNzSyrmPRsd88u5z0WEYf6pjhmNrg815W0ikpkY0bdfS5wf1T3FxFJNnefZ2b7ElTwbQn+KP/G3TdvyZUiXBD7JYIhA5eF984Fvg0vWWBmlwHzzKwewSxagPrA2pjXK4spYzaQXTQ9MzNzcKdO5R1KW9iHj22+Y9jS3KrBuf/VLTbPTReUsIvYK0XnWJQR+4/HKnR9vOJyZJyWiUpWXElny9w9a2tuEK/6pgQdynGNlXQibg1VM3vW3S8MX7/En+O1CtFuLSKSLsxslLv3Bb4Jj/z0t9z9pHLew4AXgGbAMe6+oYRL8+vUKu6+1MzmAd0IHtERvp5cbE4RkYqJW32zte2+ePaozoh5/Usc7ysikqpK2oUvswL3eBrYFejt7mvyE81sf2AZwdJXjYC/EzzKy3/8NgIYZGbfEjRyLwDOrUjhRaRyMbNqBG2/qkBVM6tFsIzo3qRofRO3hqq73xPzWru1iEjaMrM7wpc1Yl7n2xGYVc77tAUuIliqZX7QuQph2ibgbmB7gsWxPwT6xWQfTNDInUWwdut9WppKRMowiKDuyHc2cDvBEOmUrG/i+ei/XDuxbO3sLxGRFNA6/LdKzGsIHpfNpvCSUyVy91mUMjYLeKWUvOuA88Ij6bahJcVEKq1wLGtWCadTsr6J56P/F4q834Gg0l5MMIPMgDkEvQ0iItssdz8XwMy+dPfnkl0eEZF0Fc9H/+3zX5vZzQSN01vdfbWZ1QHuIGi0ioikBXd/zsx2JnhEtgMwF3jF3acnt2QiIqnDzI4CTiNYjzUHeL28QwdKWNhuq10F3OjuqwHCf28CtBd1GsnMzCQzMzPZxRBJGjPrQ7CD1C7AEoKtU781s+OTWjARkSQys7FmdkL4+gaCp+45wChgPjDUzK4vz72iWkd1FbAf8EVM2r7A6ojiiYgkw91AX3cv2IkqXGD/CeDfSSqTiKSB/I6g7OzspJZjCx3In7uRXgn0cvef80+Gy5h+QjnW24+qoXor8L6ZjSaYWNCaYKeWSyOKJyKSDK2Az4qkfU7F9tAWEUk3eUBdgk1MNgJFdwmZSbCySZkiefTv7i8B+wNTCHYv+Bk4IEwXEUkXE4BriqRdHaaLiFRW2cAwM9sFuBN4ysx2MrMaZtYBeI7Cy2SVKMotVH8Cforq/iIiKeCvwGgzu4Lg6VEbgqFPfZJaqjjYxh87ikhyXQo8BfzIn52i5xCsAJW/69VZbL5i1GYia6iGkwl6AtsRs06gtlAVkXTh7j+b2a7AAUBLgskCX5eyDaqISNpz9wXAyWZWhWCmf80tvVckj/7NbDDwTHj/UwmWpTqSYHsuEZF04kWOco27EhFJd+6+yd3nuvtvxR3luUdUParnAYe7+yQzO9fdrzKzVwi27iqTmV0GDAS6EKxJOLCYa24j2PbrcHf/KEyrSbDF1ykEKwzc7+4Px+HzVEqHXlJ07HNhE6evK9d1nzzVJm5lEkklZtYVeIegt2AuwSSqtWZ2ortPTGbZRGTbUdzv0dJ+x25Lv1fN7MWSzuVvnlKaqNZRbejuk8LX682surt/QzAUoDxygCHA0OJOhgNxTwXmFTmVBewMtAV6AdeHi8yKiERhKPAk0Mrd9yNY9P8JSqi7REQqodkxxxyC4aAnAGvKkzmqHtVfzayzu08GJgF/NbOlwNLyZHb3twDMbB+KX+blSeAGgoG6sQYAA919KbDUzJ4j6Jkt1+4HIiIV1BF41N0dwN3dzB6j5L20U1PWiZunzZxU8jkei7Q4kjiaNCdRc/fbiqaZ2QEES5mWKaqG6iCCLVQBbgReBjKAS7b2xmZ2KrDO3ceYWWx6I4IBu7GP2yYStNqLu08mkFk0/fTTTycnJ2ez6w/fc9WWFzpGcfcuTbLilif2zI83luu6LYkt6S2ZP9dxNgY4Hng7Jq0P8F5yiiMikvrc/Ssz61Gea+PeUA1neK0FvgoL8w2wU5zuXY9gJ5jDizmdEf67PCZtOVCvuHu5ezbBOl+FZGVlDW7ZsuVm13/4v7yKFbYEN12w+b1Lk6y45Ym9NLdqeF3duMeW9JbMn+s4qwq8ambf8efmJnsDo8xsRP5FWu1ERORPZlYD6GtmVdy91AmocW+ouvsmMxvl7sU2ELdSFvCSu88s5lxu+G99goZy/uuVEZRDRASCoU2TYt7/BHyQpLKISBrZs++/kl2EuDCzkwiGauYCZ4YdmAcAtcpqpEJ0j/7Hm9kB7v5VnO97GNDKzPKHEDQFXjOz+9z9PjObB3QDPgzPdwMmx7kMIiIAuPvtyS6DiEiKexi4DmgE3AscSrBs6aMEw6dKFVVDdRYw1sxGETwOy9+FoNhBtUWZWbWwbFWBqmZWi2Df2MOA6jGX/h/BdoVjw/cjgEFm9i3QDLgAKHPpA9ky6fLXnoiIiESmkbu/ZGbVCZcpdffJ4QpOZYqqoVqbYG1BKH7WflkGUXgP2LOB2909K/YiM9sILHX3/Mf+gwnWUZ1FsOzBfe6uGf8iIhWUPfCQZBdB4q2iqztkvb15mkjFjTOz/dz9m3BsKma2O+XcBCqShmp5FnAtI38W5Vjexd3bFXm/jmCzgfO2Jr6IiIhIPFXipcBmAqPN7GWghpk9ApwB/KM8mSNpqJrZjiWcWgfMK8/gWRERERHZ5jUkmGTaBHidYJvpi9z93+XJHNWj/1/4c1yqxbwG2GRm/wYucfcFEcUXEYmcmf0deNXdv4xJOwg4zd2vTFrBRERSxNYuzxfVFqoXECzy3xGoBXQCXiJY8L8LQQP5yYhii4gkSj/g2yJp3wFnJqEsIiIpycx2M7PbzOxpM7vVzHYrb96oelRvB3Zy9/z1TH8Jl5Sa5u7PmNlAYHpEsSXNVeJxPpJ6nM3/4K9aTJpIStKkOYmamZ1D0Dn5LsFKULsCX5vZJe7+Uln5o2qoVgHaAT/HpLUhqMABVkUYW0QkUT4DhpjZ9eFmJ1UIJoJ+ltxiSSLoj2aRchkCHFVkiFR3YCTB0/ZSRdVYfBT4xMxeJGg9tyJYz/TR8PwxwH8jii0ikihXEPQSzDOzWQR/kM8D+iS1VCIiqaMOwbr3sb4BSt9/PRTV8lT3m9kPwKnAXgQV9/n5a5q6+zv8uc6qiMg2yd3nmNlewP4Ef5DPBr7RyiYilVxF16zlsUiLk2QPA/eY2W3uvtrM6gB3Ao+UJ3Nkj9/DRqkW2xeRtBY2SvWESESkeBcQ/CF/uZktJlimyoC5ZnZB/kXu3r64zBrwLyJSAWY2Jeb1bDP7vbijnPeqaWYvmNksM1tpZhPM7OhirrvNzNzMehfJO9TMVpjZfDO7Oj6fUETSlZldZmbfmtk6MxtWwjXxrm/+AhwFHE2w0+jR4fvzixzF0oQmEZGKuSDm9dlbea9qBMMFegK/E4zff83Murj7TIBwP+xTCYZQxcoCdgbaAs2BT83sJ20bHY1DL9n8b4+J09eVeO6Tp9pEXiaRLZBDMLnpSILt7guJor5x90+2psBqqIqIVIC7fw5gZlUJtmu+MNy+eUvutYrC20W/a2YzgL0Jth2EYFmXG4CnimQfAAx096XAUjN7DhiIhlyJSAnc/S0AM9uH4HF8UXGvb8xscCnluT285iZ3v6e4a9RQlZRTXO9ErNJ6MWKpR0Oi5O4bzewIgu0A48LMmhFslDI5fH8qsM7dx5hZ7HWNgBbAxJjsE4ET4lUWEalcIqxvOpQUMuZ1seNTIY4NVTN7icJbpRZra7fSEhFJIY8At5vZYHffsDU3MrPqBOsKDnf3n82sHnA3cHgxl2eE/y6PSVsO1CvmvplAZtH0zMxMpk6dGrw58JwKlfXw31dV6Hrgz1hFRRw7yrgzP94YnNtz83PxiktJ99kSyYxd2RTzXT8TphX3rW7lz3VDM8sq5rJsd88uz/3iVd8UpzztPne/sKRz8exR/SWO9xIR2RZcTjBe62ozW0jMH+vuXu4u/XCjgJeA9cBlYXIW8FL+WNUicsN/6wNrY16vLHph+Isqu2h6Zmbm4E6dOgVvXrmxvEUF4MM/Kr6Uzk0XlPB1RBw7yrhLc4M9bD783+bLQcYrLke+XbHrS5PM2JVNYn+ul7l7VoVusLks4lDflMTMWhGM6W8NzAJGuvvc8uSNW0M1f5yBiEglsrWTqbDgGdsLQDPgmJie2cOAVuH20wBNCSZa3efu95nZPKAb8GF4vhvhkAERkQqKrL4xs32BUcBooD/wGjDRzI5z96/Kyh/ZGFUzqwF0ArYjZhzC1s7+EhFJFe4+Lg63eZpg7+ve7r4mJv0woHrM+/8DrgbGhu9HAIPM7FuCRu4FBDsASoLs2fdfyS6CSIWYWTWCtl9VoKqZ1QLyiLa+eYBgItZ/zOxUdx9gZkcCDwHdy8ocSUPVzA4GXgdqEnQPryAYyzAb2DGKmCIiiWZmNYHbgH5AE3dvEE6w6ujuT5Qjf1vgImAdMD9mAsNF7j6yyLUbgaXunv8YbjBBI3cWsAa4T0tTiUgZBhHUHfnOBm4vOnQgzvVNF+Cj2AR3/8DMXi1P5qh6VB8B7nf3R8xsqbs3NrPbgNURxRMRSYZHgB2As/iz52FymF5mQ9XdZ1F45mtp17Yr8n4dwfJY55W/uCJSmYUN0qxyXNeuyPutqW/ygBoE41sNwMwO488l+EoVVUO1I5tvXHsvMAN4MKKYUknocZukkBOBndx9lZltAnD3uWa2Q5LLJSKSKr4B9gfGAdXM7H1gH+Dk8mSOqqG6nOCR/zJgnpntBizmzyUORETSwXqK1KNm1pSgvhMREfgrf643fQcwHzjL3ctVT0bVUH2LYCvAl4GhwKfABuCNiOKJiCTD68BwM7sKwMxaAI8C5Rp7JbIlMjMzAcjOzk5qOUTKw93nxLx+oKL5I2mouvuVMa8fNLOvCXpTP4ginohIktwM3Af8CNQBpgPPEfQaiIhUemb2CSWMxXf3XmXlT8gWqu7+WSLiiKQr9aCkJndfD1wFXBU+8l/k7mXu0CciUon8s8j75gRLW71Snszx3EL1fXc/Knz9GSVsp+ruPeIVU0Qk2cysAcGa0Rnhe0BrRouIALj70KJpZvYvYDjB8n6limeP6oiY18/H8b4iIinJzAYCTxJsMRi7/J6jNaNFREoyB+hangvjuYXqyzFvf3b3r4teY2b7xSueiEgKuAs4xd3HlnmliEglZGYDiiTVAU4DvixP/qjGqH5IsDxVUe8DjSOKKSJxprGxZaoG/CfZhZD0deglv2+WNnH6uhLPAXzyVJtIyyRSQecXeZ8L/Jdyrqsf14aqmVUhmNllFgzUip3l1YFgdwIRkXRxH8H+13e6+6YyrxYRqWS2dm5SvHtU8/hzElXRRukmgsdkIlJEST0j+crqQcmnnpTomdls/qznjGAG6/VmVmjxanfXfwwREcDMagNNKdyBOR7oCXi4nXSx4t1QbR8WYhwQ24J2YKG7r4lzPBGRRDs72QUQEdlWmFkWcAtQpZjTvxK0G4s7B8S5oRrTIm4bz/uKiKQKdx+X7DKIbOs0/r1SuQTYjaBRms+BFUBDSljONF8kk6nMrDFwLbAH4dqCBSXTOqoikibMrAYwCOgHtARyCLZPvcvd1yazbCKyOTWQk6Kau08vmmhm7u4by8wcTZl4GagJvEbhtQVFRNLJ0wSL/f8NmEXwNOlmYAfgvCSWS9LYnn3/lewiiFTEUSWkH16ezFE1VA8Cmrr7uojuLyJxoElcW+0EoIO7Lwvf/2RmXwO/oIaqiFRyZnYLMM7M2hadMFXcevvFiaqh+gPQisLjEURkC6kHJWXNJ1i8ellMWm1gXlJKIyKSWv4g2Lm0HaVMmCpNVA3VT4D3zexFgoq8QHF7vopsCzS2SYrxEkFd9zjBloCtgUuBEWZ2aP5F7v5JksonkhZU/26b3P05MxtKMBxqi0TVUD2EoNIuOv7AgTIbqmZ2GTAQ6AK84u4Dw/TdCFrmHcJLvwP+5u4/hecNuBf4S3j+eeBGdy91RpmIyBa6KPz35iLpF4cHBPXejgkrkUiqyTpx87SZk0o+l/V2tOWRhAonTJU+fqwUkTRU3b3XVt4iBxgCHEnwGC02/RSCSQtVCHouXgW6hucvJBgz1o3gl8OHwAzgH1tZHhGRzbh7+2SXQUQknUXVo4qZNQGOAZq7+wNm1hKo4u5zysrr7m+F99iHYKxrfvoywrFgYe/pRmCnmKwDgIfyY5jZQ8AFqKEqIglgZr2Aje4+PtllEan04tCTqyEHybdFA1vLYmY9ganAWcBtYfLOBEu5xOP+y4C1wOPA3TGnOgMTY95PDNNEROLOzMaZWffw9Q0ET3heMbOiQwFERGQLRNWj+ihwurt/bGZLw7Svgf3icXN3b2hmdQl6UGOXO8gAlse8Xw5kWLiqbOw9zCwTyCx679NPP52cnJzNYh6+56qtLjdQ7L1Lk6y4yYydzM9cmvXr10dyX0jd73rmxxvLdV26/DfeArsDX4WvLwB6ASuBLyj8R7SIlFNxy+GVtlReJV4eb5tnZgcAA9z9ryVdE1VDtZ27fxy+zm8gro9nPHdfZWb/ABaa2a7u/geQC9SPuaw+kFvcZCp3zwayi6ZnZWUNbtmy5WbxPvxfXlzKfdMFm9+7NMmKm8zYyfzMpalRowYAxf18bK1U/a7bHfZGeF1i45ZXvP8bb4EqgJtZB8BiJnY2Sm6xRERSk5m1Bs4Jjx2AUmfPRdVQ/cnMjnT3D2LSegM/xjlOFYI1DHcgWKtrMsFEqm/C893CNBGRKHwOPAG0IKxsw0bromQWSiTVZQ88JNlFkAQyszrAyQRPwg8BxgN3AW+6+5rS8kbVUL0GeNfM3gNqm9kzQB+gb3kym1m1sGxVgapmVgvII3istohgQ4G6BCsDLAWmhFlHAFeb2RiCntxrCMaxipRJuzTJFhhIUM8sBB4I03YBHktWgUREUkm4pv7JBEtUjQD6u3u5x21FtTzVV2bWFTibYN3U2cB+5ZnxHxoEDI55fzZwO0Hv6OMEKwGsIeg5Pcrd14bXPUOwXmF+z+3zYZqISNy5+2KKrKHq7u8lqTgishU0NjYyZxIsL3o/QQ9qhSYpRNJQNbNr3f3BsFCx6Ve7+8Nl5Xf3LCCrhNOvl5LPgevDQ0QkUmZWneAP63OAlgSV8UvAXe6+PpllE0kn8dpGWkMOkqItwSpQ1wJPmdk7BPXkh+6+qazMkSxPxZ9LUhU1KKJ4IiLJcD/B+PuLCcbEXwwcCtxXnsxmVtPMXjCzWWa20swmmNnR4bndzOxbM1saHh+Fu/Pl5zUzu8/MFofHfeH60iKRyczMLFhbVLY9ZnZZWK+sM7NhMemR1TfuPt/dH3L3rsDBBEOlRgCzzewBM+tWWv649qjG7G1dNVz4OvZD7EiwbIuISLo4FegWDgEAmGpm3xOs4XxVOfJXIxga1ZNg/NYxwGtm1gXtxCci8ZfUnT/dfQIwwcyuBY4mmFz1NVCrpDzxfvT/QvhvLYKxqQVlA+YDl8c5nohIMpXUo1DenoZVFB7m9K6ZzQD2dvc30U58IhJHidz508x6A42A/7h77Br3uPtG4F2COq9haWWOa0M1f99rMxvh7v3jeW8RkRT0OjDazG4n6BFtSzDE6bUtuZmZNQM6ErOsXrgTXwZBL0fssKpy7cRX0uYmmZmZTJ06NXhz4DkVKufhv1d8w4aCWEVFHLuyxY1nbIq5z+rVq0uPEYe4qfRdl7bpSanfQWI/c0Mzyyrmsuxwzfhy29r6poiVwF+AJ8xsCjAGGOPuk2IvChvJJYpq1n//cImpgwjWOJ0D/Nfd47PKt0gSxGswv6SV6wkapk/y52SqVwgerVVIODFrJDDc3X/OT9/anfhK2twkMzNzcKdOnYI3r9xYobJ++EfFV9+66YISZkhHHLuyxY1nbI7cfB32OnXqAFDws1Me2/B3XdqmJyXGjVPs0hSJvSychL7V4rHzZ8y9vgbOMLMqBO3BY4GRYQ/qmPD42N1Xl1amqGb9dyLo0q1NMP6qNbDWzPq4+5RSM4uIbCPCmf23UfIE0nIJK/KXCHbwu6yYOFu1E5/IltByTZVTvOubcGb/5+Fxk5ntQDAe/3xgqJl95+5HlZQ/qgX/nwaeBR7M/xDhwNmnCBbtFxFJC+Ef5t0IehwKuPvQ4nNslt8Ixvc3A45x9w0lXKqd+EQkUSKrb9x9LvAc8Fz4JKlnaddH1VDdAzi8SEv7UeCWiOKJiCScmd1M0Js6EYh9fOUUnlBamqeBXYHesVsJmtnhaCc+EYmjZOz8aWaDSzrn7rcDH5WWP6qGag5BC/mTmLRDwnQRkXRxJcGuez9sSWYzawtcBKwD5scsS3gRwTAA7cQnIvGUjJ0/OxR534ygTfhueTJH1VC9Gfi3mb1LMBi3LcEg2rMjiicikgxrgJ/LvKoE7j6L0pey0k58IhI3ydj5s7hVoMKlqwaUJ38kO1O5+7+BPYFJQL3w373dfVQU8UREkuRW4HEza2FmVWKPZBdMRCSFfQL0Lc+F8d6Zqg5Bt/LuwPfAPe6+Lp4xRERSyLDw37/EpBnBGK6qCS+NSMS0TJ9UlJkdViSpNsGufjPKkz/ej/6fBPYBxhJsxdUE7UYlIumrfbILICKS4p4r8n47golbR5Ync7wbqkcBe7n7PDN7HBiPGqoikqbCMab566A2AxaEawaKiAjg7jvGvg+X5LsROBm4t6z88R5HVdfd54UFmw00iPP9RURShpnVN7MRwFpgLrDGzIabmeo+EZFihBOz7geuKs/18e5RrWZmvfhzFmvR97j7J8XmFBHZ9vydYM3B3flzhZO7wvRyzWgVEUlnZlZ0vH5t4AyCP/DLFO+G6h8UXuR6cZH3TrAOV9qqWgXuv3x7rnnsj0ju36xxVRYs2VjwPnPvOvTcM9h7+YsfVvPRN6VumbvFnru5Odnfr+bDr1fxx9KNZWeIk1o1jYHHNqBti+p8M3kNo8bncsu5Tdhvt9r8Nnc9941YQs6ivEhim8FJmfXYoWk1Rn+ey9IVG7myX2NabFeN739eywv/XkZeBF9Fsv4bAxyxf12OPqgu7VpUp1YNY+Gyjfw8cz0vjV3O7AXRfM+QvJ+vODgK2DFmr+ppZnYu8GsSyyQikkrWU3gZPif4w/788mSO66N/d2/n7u1LOdK6kQpQpQp03blmZPd//pYWBa/7HJLBJSc3ZOqsdfw8ax0X9G1I3x4ZpeTecq2bVWevXWox8o6WPHTF9hy+Xx1qVi9t+cf4uOqMxuy4Q3X+++Maunetwz2XNGXlqk1c8fACfsvZwOWnNYos9kUnNuTALrVp1awaD/xte/r2rEf296t5/eMV7Ne5FuceF83T3WT9N+5/TH3OObo+3/y0ljc+WcmiZRv55NvVzF+Sx6NXN2OvTtH9XCfr5ysO1gJNi6RtR7CAv/x/e/cdJ1V1/nH88+zSWXpvgiKIBcWosURlNfYkGksUO/YeNbHGwtqiiT879hJRNMbesYtdrKAiWEBAei+7bN/n98edXYbdZeudubOz3/frdV/M3HPvfc4Uzp459xQRERhAsJBA+dbF3Tdx9xpXpCqXqAn/09r4q/tuMC0jwX9bLe76B+2eRc79S/n+lyIApvxYyEXHdeOF93NDj1tc4vz91sX06prJPju259gDOvHXI7ryweS1vPZpHt/8lJi/y9tv0YZjrpxPfqHzzud5PPfv/vzx3rkUFjt3P7OS/1234c+isfbYrh0nXruAzAzj2X/1481JeRWttz/NKeK6M3ty/wurQo8b1Wd84G4dOOvGhRWtue9+uZarT+3Oydct5LufCznt4C6cdsPC0ONCdN+vEDwAvGlmN7Pu1v/5wH2R5kpEJEW4e6NWJVVFtQE6tsvgnmdXsGBZ1VuhLTON686s3MASHvd1j7t1yqyowABMn11Ej86Jnbpx0fJSHp2wmkcnrGarwa3ZZ8f2XHtaD9asLePoK8NfITfD1r1mBzAoi+3w+DcjAdq2ySAv3wFnbaGv18Vg9sISOmclZk73yD5jg+Wr1912X766lE5ZQazPpxUw5uTEFxfJ/n6F4DqCpaGPAvrGHv+b9bs8iYhIA6mi2gA//VpEYbHz9Q9VW3latqh5PcTGatXSuOT4bgBkZBhdOmawYnUwG077tkZJaWIqb1bNi/puRiHfzSjkjieXs+s27RIS96sfCsg5pTsffZPPzsPb8vUPBZxxSBde+jCXA36XxdSZiWtpW51bSrs2xtoC59b/Ll8vrVNWBgWFiXmvo/qMP5uazyXHd+OJN1ZjZhy5T0e++iHo657VNoO8gsTNuhTV96uxYqNXHyJNK6Z3fT6TM3eof4+tAb1aNKpP8/L8Irq2bVXn49u0Ngb0bMH8pSWxH5cNV9Nr3qRfS2bOK27U9cu5O5/NW8G8NcH/sX4d2vDbfl2w6v4zVKNj+wxW5zX+/2R93+uw5BaV8OOyXDbp0o7ObWqPH+ZnXJfXPHSjVvTskklZGfy6uDihffTLhfn9SieqqDbAI6+uoqCo+v8oxSXwtwQNpAJ47LV1t5qfeaeYrLbrKjHbbNqGL6bVaRBdvX3784YrhMUlwW3iRLjl8eWcdFBndhneluffX8NPc4oYc3J37rigF9NnFXHj+GUJiQvw3MRcOrbPYG1B0Fcz3i5bt2XS1PyExI3qM77zqRWc9Zcu/OvsngB8/n0BdzwZVNCz2hp3PrUiIXEhuu9XQ5nZnrUd09RmOLn5k5/Xe+7uXP/hTxSUBN+9v+28aZ2v9e9zenLk5XVrAf9ozjJOfulrMsx46MBtufzdacxckUdRqfPkYTuw84CuVc45d1QXbnsi+D5uNbg1l53QjflLS+jXowW3PL6cSVPr9n+kvq/53kt7s2BpCe9+sZZ3vshj9sKGVV7emLGYM1+ZwpBu7enXoS0Ac1fn8/PyPO76wzbsM7jnesf/dos2nDuqK0tXlnDHkyv4x+hutGpptGxh3PDIsmobTarTkPd6QK8WnHlYF9xh7JMrOGb/juy6TTvmLi7mhkeWMaeO78GZr0zhrj9sA8CHc5Zx1DNfMLhre35ense9f9yGA4b0Xu/4sD7ja9//gct33wyA75es5s9PTKK4zHF3/nfYDuzYf/3XvPWQ1pxxSGdy852hG7XiuxmFdGiXQUmpc/24ZSyp42DPbxet4pSXJjNvTQH7b9qTf+21JV1ilePf3j+Rz07JrnJOWN+vdKOKagNMqaW/XCL70z3y6uoNpn38bT4ff5uYytOldy1JyHVrk5vvFYVVuXNvTtwPgXjPvLtmg2kTPs5jwsd5CYkb1WecV+D8+9Hl1aYtWFbKgmWJiQvRfb8a4cFa0pvcDCdjJk7ngCG92LJHh4ruJ6XurCms/o/l2X+pfiCjWdACX1fnv/4tTx62A7lFpfzh8U95ftSO7LpRN75asJJzJnzDRyfuXuWcLQatG9h3wp86ceW9S/jp12L6dMvkypO717kSU9/XPHNeMdc/vIw9d2jHtaf3oKDIeSdWqYifqaM25772DW8dtwuDOrdfb/8vK/I44PFPmHbWXuvtP/mgzlx612Ky2mZw41978o+7ljBtVhEb9W7BZaO717nveEPe678f1ZX/vbWGtq2Nm87tyX3Pr+Tfjy5n5+Ft+evhXbng9rqVx5/OXVe2XPHuNJ4ftSO/6dOZmSvyOPypz6tUVMP6jJ+dNr+ionrhm1O5bb+t2X9ILz6bt4LzXvuWj09a/zWfdVgXLrpjMatyy+jdLZMzD+3CX29axHbD2nDhMV256I66lVVnvDKFnOxh7NS/Cw98NZtd//MBL47aicFd21NcVn1DV1jfr3Sjimoj9O6WycZ9W9G6lbFkRQmz5heTV5DYfpOb9m/Jz3OTf2sgqrjl+vdswaA+LWnbJoP8gjJ+mV/MvCXJ/7XZsX0G3TtnMmdhcUKmpooXxfcr3oCeLdh+izZA0Lo6d3Hi3u+ov1/15e5pt3Tq1DP35O9vfEdeUSljsjejXcsWjJsyhzHZw6o9fr+d23P3MysoruZrsef2de+qUVzmDO8VzKDRo30rdt0o6Pbymz6dyS+u/T9Z+zYZ/PRr8N1ZsKyUjDreOof6v2YcZi0o5qEXV/HQi6sYNrAVe2zfjtv/3ovFK0o55/8W1SluSZnTv2PbKvv7dWxLcTVde9ypaLksLHamzQr6rc9ZWFJtt5kNach73bZNBp/Efhyf8MdOFXc3Pvk2n9F/aNjMJ6sLi/lNn84AbNKlfcW4gw1pzGccb/6aAvYf0guA3/brQn41hXhGBqzKDVrUFy8vpVfXoJr05fQCzjqsc51jrSkqYb9Ng1gX7DKE7fp0Zr/HPubRg7fbcPfAkL5f6UYV1Qbo2jGDy0/sztabBr/43GFtoZOZAc+8s4b/vBz+SPBy91zSm/lLS3hjUh5vTkrer6yo4vbskskVJ3VncP+WzF9SQl5+Ge3aZtC3ewtmzivmmgeXJmzezY16t+CKE7vTt3sLxr++mplzi7jk+G60bmWsySvj4rFLmLUg/MpVVN+vm8/ryV1Pr+DnucXsNqItFx/Xjck/FuAOJ/ypMzc8vCxhrblRfb9knY06teOpv/yWF6YvYO9HP+b8nQbXePwPs4uYtaCYqTOLqqQdX48KTHwl5frfb7FeWlFp9X0wB/Ruwf2X9caA3t1akNXWyM13zKBFPf6q1fc1V65hTJ9dxPTZRdz9zMp6TUt44oiB7HD/e4zash8DOgUV1l9X5fPE1HmctO3AKsfn5pfxx12zaN/GWLO2jMP27MDEL9fym2FtyC+sez/VhrzXGXGN40+9s/5dpvq819OX5rL13e/gwKyVa1mRX0SXtq0oc682dlif8cwVeRz4309xnLmr81lbXEK7lsEFiquJ++PsIi44pitf/1DALlu3ZfJPQctt65ZGRj2n9VlVUEynNi0B2GPjHjxz+G859MnPWJ5f9f8MENr3K92ootoAFx7TjdkLi7n2oaVkmHHsAR1ZsLSEDybn87eju3Linzrx0EuJqUwUFjvjJ6xm353ac9wBnfj250Je+zSPD75eu8F+s0057kXHduPbGYVccNtiCovXxWnTyjj2gE5cfFy3hC2ucPZfuvDqx7mUlQWPb39yBQdfNI+MDDjz0C6c8ufOXHZ3+Leso/p+De7fqqJV87gDOvGPu5bwTazv6JabtOJvR3VNWEU1qu+XVHXQsD7stUkPciZOr7bVr1zO/UspKq7+s6nPDA3X7LF5ReXhz8PWTTc3Y3kex22zUbXnnHD1gvWe58cGNnZol8HDDfghV9fX/L83N9wtpz5dvi7dbSgHDevNiz8s5JNfg1vi/Tq25bFDtmOLHh2rHH/DI8s4Zr+OOHDRHYv5/fbt+dc5PVi0vJSbHq++u051GvJev/BeLm1aGwWFzotx0+L17REsfFJX0876/XrPs1oF1Y/l+UVcvcfmVY6v/BmXlwP1/YxfGLXTes9LY7fdF+UWcMb2VW+M3Pz4cv6waxZbbtyar6YXVHTxcuDisXX/W3Px74YwbekadorrA7t1r068fdzvuOa9H6o9J6zvV7pRRbUBthrcmsvuWUJZ7MfYnU+v5LGr+/LEm2v417hljL2wV8IqqmVl8MakPN6YlEevrpnsvWN7jtmvI+ce0YX3v17L65/m1dqHtinF3XxQKy65c3GV2+wFRc7DL6/k+Rv7hx6z3JABrbjojiVBxfSwLrz+aVBglZXBQy+t5JGcxMzhGtX3q7TU6dAugzVry+jWKZPv4mZUmPZLUcUtsESI6vsl1WvfqgU37rNVjcesWRvOLBAHbtan2v2Du7bnot8NqTZtQy3uq/PK+GByw35M1eU1Vx5U2Rhb9OhYbaW0OktWlHLLf9f11X/6nTU8/c6G+9BvSEPe65c/rH7O5vlLSrjz6ZV1jj2wc/XdQbq3a80hm1ctS8P6jEcO6l7t/l5ZbTjrt1W7kZeWsV6FHNbNsFCfOz1HDR9Q7f6NOrXj/gO3rTYtzO9XOlFFtQFWrCmlX49106/07d6CtbGpexavKKVdm8TMr1nZouWljJ+wmvETVrPlJq3Yd6csrjq1O3++cF7axF28spSdh7ettmDaccu2CV1us7wbVFlZUDGObz3KL/SErZwU1ffrva/X8tcjunDT48t587M8jtq3I+MnrMYMjty3I7/MT04f0qi+183daz8vquhTt7KgiAvemMpn81awVc+O3LLvVvTKarPe8W1aG6P27sjuI9rSvUsLSkqCuYZf+iC34kddXSzMLeCqidPJMOPqPTbnjs9m8My0BWzePYvb9tuaPh3aVDnn3kt688HktbzzxdqELaG8/2MfM+HoXRISd3VhMdd/8CNz1+RzwKa9OXL4uh/c8aPjy2VkwAG7ZPG7bdrSPTaP8tKVpXw0JZ8JH+eygbv2VXwxfwUXvjmVfh3acP3vt+DEF7/ms3krGNoti/v/NIIRvTtXOad1S+PP2Vng8OzENey5fXt2G9GWOYtKghlwQpimr7r3Oqzv16qCYq7/8Eeen76AxXmFmEHP9q05aLM+XLLrkCpTY51yUCf+99YaVueVMXSjVlx5cne8zGmRaVz/yLI6t2zWNy4k53vdFKmi2gBPvLGaW87vxXtfrcWAkdu1Y1zsVsSgPi1ZkMAv2Ib6kE+dWcTUmcu548n0inv7/5Zz1Sk9OGzPImbMKyYvv4z2bTIY3L8lg/q2Iue+xI0Wn7+0hJ5dMlm8opSDLpi7Xtom/VqydFViKslRfb/ufmYlfz+6K09e149FK0rYuE9Ljtkv6Gu4aHkJV9ybuPc6qu+XrPOPt79fN/jjjan0zmrNS0fuxLPT5nPay5N5vtIt1MtGd+PDKflcNHYJ2b9pR5vWxrtfrOWY/TvSv2cLHnyxbq3+o5//ij8M6UVecSl7jPuQo4f359WjduL56Qs4/ZXJVW7dAmS1yyCrXQY3ndeTFauD6ePe/XIty+r5f/KrBSur3e8OkxdWzX9YcU944SuGdM3i0M378tDXc3h62jweP2R7WrfIXG90fLl/HN+N3PwyHnllFUtWBrF6dM5k353a848TunHNg3Wbpu/MV6ZwVfbmrCwoZpeHPuCWfbfizWN/x9szl3DGK1P45KSRVc65+LiuLF5RSutWxvVn9mTOwmL+9+Yadtm6LeeP6sr14+oWu77vdVjfr8Of/pw9B3Vn4uhd6R37sbUwt4Bxk+dw+FOf88axv1vv+B23alux4uBph3TmmgeX8sPsIvr3bMFlJ3TjjH/VbUBTfeNCeN+vdKOKagO8+nEe85aUsNNWQV+max9aWjGP3bwlxZx3c+JG5t1cS3+k6kbgNuW4X/9QyDFj5rPbiLYM6tOSzh0yKSgs4/VP8/hwytJQJrzekBvGLSM3v/rrd2qfwX9eWpmQuFF9v4qKnesfXsaAni0YNqgV3bu0oKjImTmviCk/F1Z0RUiEqL5fUr0v5q9g8unBVLHn77wp46bMqXJM724tKlq2nn5nDXdd3IvxE1bz70eX858r+tS5IrEor5BzdgwGMd31+S9cvOtQAM7ZcTAPfj272nNy15Zxz7MruefZlQwf3Jo9d2jHvZf0ZvaiYt75PI9XPqpbi9sO909k5MDuVNcmuLKg6h2EsOLOWJ7HM4fvCMCfh/Xluvd/YM9HPuLFairlAEM2asXxV63fZ3PpylKmzSpi3Jjqb+dXp7jMK0a9X/zWVA7boh8Av9+kBxe8Wf1/8P49W3J1rCL89PX9Kqaj+nZGIfdf1rvac6pT3/c6rO/XrJVrK75TFdfOasPFuw7loclVv9eZmUZGRnAnrXVL44fZwcCnuYtLaNmi7nfR6hsXwvt+pRtVVBtowdISps4sZNaC4vWm7Skugd1GtE1YX5N3vlhLzy6ZDN2oVZXYEEwLk4jYUcWFYOT/mrVlPPDCKopLnYN2z2JQn5YsW13KpO8SM/k9BFO/bNq/JX17tGDSdwUUlzgH7p5F3+4t+OqHAt7/OnHzikb1/erZJZOBfVoybVYRcz9bP0YiP+Mov18SWJxXyM2f/Iy7s7qwBHevWCWpumkfC4qcrQa35rsZhewyvC1rYj8a3TfcQl6d+JHox22zfr++uizC9u2MQr6dUcgd/1vBdpu3YY/t2tX5D/rm3Ttw7x9HMKRbVpW0Abe8nrC4haVllLlXTLN02e6b0a9jG3Z/+ANyi6r+KluztoyR27bl/cn5FfO9msHI37Qjtx59hdu0yOCNGYtZVVCMGTw/fT5/HtaX92YtJbMOH1qVRU7qcde/vu91WN+vgZ3a8u+PfuL4bQZUdF9ZlFvAw5PnMKCagXMvvr+GG87qwX9fX83n3xdw1l+68MHktWw7tA0z6jGFXn3jVtaY71e6UUW1AXbYog1XntSdhctK6Ncz+NV3x/9WVBTm5x/ZNWF/VKOKHVXc/Xdpz4l/6oy7c+Q+pXwwOZ+eXTLJzDCuOLE7Y59awWufJOY/b02xL09g7Ob2GUcdWwKn/GZQxUT3x2+zEUvXFtGjfWsW5hYwonfV6aZu+e9yLji6K/16tmTW/CL+b3zQKt4pK4Pn36v7QJ+DNutDblEJWa1acO2e66ZM+nl5LptVU6mBYEnLyso8mO/38+/r/uM1J3vYBufwvGP/4QmL+6ehvXnnlyXstcm6FahGjxhI76w2nDPhmyrHX/PgUk49uDN/HdW1omKa1S6DyT8WcM1DS+sc954/jOCit6aSYfD6Mbtw9+e/MPr5r+jXsS33/2lEtef8MKeoYtT/jePX3fmI7ztfF/V9r8P6fv3vsB244cOfGPnwhyzKK8SAXlmtOXBoH578yw5Vjn9uYi4z5xVz4O5Z9O/ZkswM6N+jBR9OWcv4CXUfxFrfuBDe9yvdqKLaACcf2Jlr/7OUSd8V0KVDBpeO7sY1p/dgzH1LKCmt36+9phI7qrhH7NWR829ZhAH/ubIPY59aUTFv43tfr+W0gzsnrKIaVezm9hlHHVsCY7KHMX3pGuatLmDH/l0qpg/qndWGo4ZXnV1j5rxgCc3unTL5flZRxaCaVbll9Voc4uo9Nq827qZdszj5N1XnFAW49qFlDOjVgu6dM5kWFxuCHz11/aN+2Bb9mL50DW/PXLJebIA2LTITFvffe29Vbdz9Nu3F7fttXeX4RctLeertNTz11hrmLy1hQK8WbLlJa2YvKGbhsrr3X9ymdydu3Xc489cU0L9jW27bf2tu2z+I99rP1Xcpuumx5Qwb2AonmDt3YO8W7LBFW35dXFyvVQLr+16H9f3q0rYVJ2y7EXsP7sFO/buuFzd+AGG85atLefnD3EZ9xg2JG9b3K90kZ3h6mukTuxUMsGJNGZfcuYT8wjKuP6tnwkaCRx07qrhdO2Uyd3EJvy4uobDY15tc/MvpBQmdMimq2M3tM446tgTumDSDg56YxB2fzWCru97mhenr+kT+4+3vqxx/cHYWV5/agz9nd+Chy/uwy9brbmeefGDnhMUtj33t6T04uJGxb2/Aaw4jbk2v+bJ3qsY97oCOnP2XLpw7qguH79WBE/7YmVYtjSP36cjR+9VtiisIXu/B/6vfe33cAR05+/AunDeqCycf1Im/HtGVNq2NI/euf+wovl/lccd+NjOpn3F944YZO92oRbUBcteW0aNLJktiUyOVlQW/hC48piv/d27P9VbySJfYUcUtKCwjMyOY267ylCStW1qtS+81xdjN7TOOOrYE7vtqNl+emk1WqxbMWpnHYU9+zqyVazl3p8HVdkX8w++yOP1fCykodHp1zSTnlO707taCZ99dU2WFnTDjlsc+7YbGx76/Aa85jLj1fc27b9uOU/+5kJYtjKdv6McRl81jbYHz5FtruPOiXjz22oYniq/8er+o53sdZuwovl9Rfcb1jRtm7HSj4r8BvpxewH47ta+y/8bxy5kxr4hW9RgZ2FRiRxX36x8K6dcz+D11+/9WrJe20/C2zJyXuLk9o4rd3D7jqGNLoMy94vbkoM7tmTh6Vyb8vIi/vf4tXs2Psgyj4tbkouWlnH/rYn67RRvOOLRzvf6m1jdulLGjiltaFvRVLCwO5hJdWxAcU1TseD1m42jIex1V7Ob2GYcZO92ootoAtz2xnCffqr4z963/XcFRV9R9+cCmEjuquNePW8achdX3R/r6hwJy7q/7QIKmEru5fcZRx5ZAr/atmbxwZcXzrFYtePmonVi6tohvF1dtNVuxpozB/VtWPC8odP5x9xI6ZWWycb+WVY4PK26UsaOKW1KyboGR029YWLG/fRujPjd2GvJeRxW7uX3GYcZONyl569/MzgZGA8OB/7r76Nj+nYBrgO2AUmAi8Fd3XxBLN+AG4OTYpR4ALvEN/XxpoJJSKKlh3pRErpYUVewoX/OGrMpN4MSeEcZujp9xKn6/ksHMWgN3AXsBXYEZwKXuPiHZ5d0jB29Hi4z1221aZGTwyMHbcdp2g6ocf/24ZRXrppcrKwvmH37pg7qPyq5v3ChjRxX3vFsWVcwlHP/pZmYaNzxStwn3GxI3ytjN7TMOM3ZNUr1+VZ2UrKgC84FrgX2B+AnHugD3Aa8DJcBY4D/AfrH0U4E/A9sQzPD2JvALcE8yMi0iUk8tgF+BkcAc4ADgSTMbTpLLu/41zO34u426Vdm3dOWGfzzEDzwMO26UsaOKu6EFL1bnldVr0ZOGvNdRxW5un3GYsWvR5OpXKVlRdfdnAcxse6B/3P4J8ceZ2VjgvbhdxwM3ufvcWPpNwCmooioiKcjd84CcuF0vm9kvwHbu/kz8sSrvRKSxmmL9ypLQattgZnYt0L+8abqa9POAUe6+U+z5KmAfd58Ue7498K67d6jm3Gwgu5rL7gB83sAsDwJmNfDcxogqbpSxm1vcKGM3xbj93f3k2g9LLWbWC5gNjHD36ZXSzqMB5V0NZd2JwMwGZrUzsLKB5zZWVLGbW9woYze3uI2NvQnwUDX7J7r7xOpOSGT9KnTunrIbQfP0wxtI2xpYDuwWt68UGBb3fAhBE7UlKb85Eb1PkcRtjq9Z73X6x43w820JvAXcW01aSpV3+n+Q/nGb42tuTu91U6pfNclR/2a2KTABONfdP4hLygXiZyDuCOR67F0VEUlFZpYBPAoUAWdXSlN5JyJJkYrlTZOrqJrZQIJWh2vc/dFKyVMJOvqW2ya2T0QkJcVG0z4I9AIOdffiuDSVdyKSFKla3qRkRdXMWphZGyATyDSzNrF9/YB3gLHuXl0H3keAv5lZPzPrC/wdeDhpGRcRqb+7gc2BP7l7fvlOlXciEramWL9KyVH/wOXAmLjnxwBXEfSH2ATIMbOc8kR3z4o9vDeW/m3s+QOxfSIiKSfWgnEaUAgsDBpXIbZvU1TeiUi4mlz9KqVH/Tc1ZpbtGxhhl45xo4zd3OJGGbu5xZXa6f9B+seNMnZzixt17FSniqqIiIiIpKSU7KMqIiIiIqKKqoiIiIikJFVURURERCQlqaIaAjM728y+MLNCM3s4STFbm9mDZjbbzNaY2WQz2z8ZsWPxx5vZAjNbbWY/mllSl6s0syFmVmBm45MUb2IsXm5s+yEZcePijzKzaWaWZ2YzzGy3BMfLrbSVmtkdiYwZF3uQmb1qZivMbKGZjTWzVJ2hpFmJoqyLxY2svGtuZV0sZmTlXbLLulhMlXcpTBXVcMwnWI6surV2E6UF8CswEuhEMOXEk2Y2KEnxrwcGuXtH4EDgWjPbLkmxAe4EPk9iPICz3T0rtm2WrKBmtjfwL+AEoAOwOw1fo71O4l5nFtAbyAeeSmTMOHcBi4E+wAiC7/iZSYotNYuirINoy7vmWNZBBOVdFGUdqLxLdaqohsDdn3X354FlSYyZ5+457j7L3cvc/WXgFyApBai7T3X3wvKnsW1wMmKb2ShgJfB2MuKlgKuAq93909hnPc/d5yUx/qEEBekHtR0Yko2BJ929wN0XAq8BWyYpttQgirIuFjey8k5lXVJFXdaByruUo4pqmjCzXsBQkriEopndZWZrgenAAuDVJMTsCFwN/C3RsapxvZktNbOPzCw7GQHNLBPYHuhhZj+b2dzYraG2yYgfczzwSBLXkL8VGGVm7WKrpexPUHiLAMkv75phWQdJLu9SpKwDlXcpRxXVNGBmLYHHgHHuPj1Zcd39TILbM7sBzxKsrpNo1wAPuvvcJMSKdzHBqhz9gPuAl8wsGa0qvYCWwGEE7/MIYFuCW58JZ8HKSSOBccmIF/M+QYvCamAu8AXwfBLjSwqLorxrZmUdRFPeRVrWgcq7VKWKahNnZhnAo0ARcHay47t7qbt/CPQHzkhkLDMbAewF3JLIONVx90nuvsbdC919HPARcEASQpev/X6Huy9w96XAzUmKDXAs8KG7/5KMYLHv82sElYH2QHegC0G/NWnmoizvmktZB5GVd1GXdaDyLiWpotqEmZkBDxL8Ej3U3YsjzE4LEt9vKxsYBMwxs4XABcChZvZVguNWxwGr9ajGBnFfQfArO/42VDKXkzuO5LYudAU2AsbG/kguA/5Dcv9YSQpKofKuuZV1kITyLgXKOlB5l5JUUQ2BmbUwszZAJpBpZm2SNL3E3cDmwJ/cPb+2g8NiZj1jU4hkmVmmme0LHEniO/zfR/AHYkRsuwd4Bdg3kUHNrLOZ7Vv+uZrZ0QSjUZPVj+g/wDmx970LcD7wcqKDmtkuBLf+kjX6lVgryi/AGbH3ujNBn7FvkpUH2bAIyzqIoLxrbmUdRF7eRVLWgcq7lObu2hq5ATmsGw1avuUkOObAWJwCIDduOzoJr7cH8B7BaNTVwLfAKRG97+OT9Ho/B9bEXvOnwN5JfJ0tCaYwWQksBG4H2iQh7r3AoxF8riOAicAKYCnwJNAr2fnQVu1nk/SyLhY3kvKuuZV1ca85kvIuqrIuFlvlXYpuFnujRERERERSim79i4iIiEhKUkVVRERERFKSKqoiIiIikpJUURURiWNmV5vZHDMrMTM3s9FR50lEpLlK1rQiIiIpz8wOBK6IPZ0JLIltIiISAVVURUTWGR73eAt3T8ZSmSIisgF1vvVvZieb2RdmttbM8szsfTPbO5GZExFJFjObCFwbt6sgdut/UDQ5EhGROlVUzexm4H5gO2AewaS0uwGvm9lBicueiEjSfE9QvpWbFNvUqioiEpFaJ/w3s4EEfbUyCFYguSq25vJzwEHAj+6+WcJzKiKSYGaWA4wBcPeErm0uIiK1q0uL6m/jjssxMwfKCCqpAEPNrFsiMiciIiIizVd9B1N9T7D+b2UtQ8iLiIiIiEiFulRUPwccMGACcKHH+guY2cbACHdfmLgsioiIiEhzVOutf3efBdwRe/p3YJ6ZfW1mCwn6rp6buOyJiIiISHNV11v/5wHfAScDWwGdgfnAu8D4RGRMRERERJq3Wkf9i4iIiIhEoc4T/ouIiIiIJJMqqiIiIiKSklRRFREREZGUpIqqiIiIiKQkVVRFREREJCWpoioiIiIiKUkVVRERERFJSSlXUTWz1mb2oJnNNrM1ZjbZzPaPpQ0yMzez3Ljtirhz/8/MfoqdN93MjovulYiINJ6ZTTSzgrgy74e4tKNiZWWemT1vZl2jzKuISNhSrqJKsFrWr8BIoBNwOfCkmQ2KO6azu2fFtmvi9ucBf4qddzxwm5ntkpxsi4gkzNlxZd5mAGa2JXAvcCzQC1gL3BVhHkVEQhdqRdXM/mRmT5nZtrHn/67vNdw9z91z3H2Wu5e5+8vAL8B2dTh3jLtPj503CfgA2Lm+eRARaQKOBl5y9/fdPRe4AjjEzDpEnC8RkdCE3aL6D+BfwDVmlg0Ma+wFzawXMBSYGrd7tpnNNbP/mFn3DZzXFtih0nkiIk3R9Wa21Mw+ipWtAFsCU8oPcPcZQBFBeSkikhZahHy92e7+hZkdCjwObNqYi5lZS+AxYJy7TzezLILK52SgG3BnLH3fak6/h6AQf30D184Gsivv32uvvcZsvfXWjcm2iETs5ptvftDdT446HyG5GPieoBI6CnjJzEYAWcCqSseuAtZrUVVZJ5Le0qy8q8LcPbyLmV3g7v8Xe5wF3O3uxzbwWhkEld2OwEHuXlzNMb2BBUBHd18Tt/9GYE9gD3dfXZ+4OTk5npOT05Asi0iKMLOr3D0n6nwkgpm9BrwC7AV85O7/jktbA2S7+5e1XUdlnUh6SOfyDkK+9V9eSY09zm1EJdWABwkGCBxaXSW1PEzs34rXYWZXAfsD+9S3kioiUldmlhkbkd86yaEdMIJuTdvE5WcToDXwY5LzIyKSMGHf+gfAzDoBfwW2Jbg9VcHd96nDJe4GNgf2cvf8uOvuCKwEfgK6ALcDE919VSz9UuAoYDd3X9b4VyIiUj13LzWzjUng7Clm1hnYEXgPKAGOAHYHzgVaAp+Y2W7AV8DVwLPxd5dERJq6hFRUgaeATOA5IL+WY9djZgOB04BCYGHQuAqxfWXAP4GewGrgTeDIuNP/SdCP6+e48/7p7v9s0KsQEanZVcDdZjYGmMu6uzy4e1kI128JXEswMLUUmA782d1/BDCz0wn66XcD3gJOCCGmSGSys7MBmDhxYqT5kNSRqJaAnYD93X2suz8Yv9V2orvPdndz9zZx8wZmuftj7v5fd9/Y3du7ex93P87dF8ada+7eutJ5oVRSx44dy/bbb0/r1q0ZPXr0emkPPPAAm266KVlZWey3337Mnz+/Iq2wsJDTTz+dXr160bVrV/70pz8xb968ivTly5dz8MEH0759ewYOHMjjjz9ea15eeeUVdt11Vzp37kzv3r05+eSTWbNmXSNKYWEhJ554Ih07dqR3797cfPPNFWmffvope++9N127dqVHjx785S9/YcGCBRXp7s7FF19Mt27d6NatGxdffDFh9mMWSTMPAMcBMwl+JBcTtHxuqLtSvbj7Enffwd07uHtnd9/J3d+MS3/c3TeKlYkHufvyxsYsLCzkpJNOYuDAgXTo0IERI0YwYcKEivS3336bYcOG0a5dO/bYYw9mz5693rkbKnsA1q5dy5lnnkn37t3p1KkTu+++e635ufHGG9lqq63o0KEDG2+8MTfeeON66bNmzWKPPfagXbt2DBs2jLfeeqsibdy4cWy33XZ07NiR/v37c9FFF1FSUlKR3pDyV0SSK1EV1Q8JYWqqVNK3b18uv/xyTjzxxPX2T5w4kX/84x+88MILLF++nI033pgjj1zXyHvbbbfxySef8M033zB//ny6dOnCOeecU5F+1lln0apVKxYtWsRjjz3GGWecwdSpNc+otWrVKi6//HLmz5/PtGnTmDdvHhdeeGFFek5ODj/99BOzZ8/m3Xff5d///jevvfYaACtWrODUU09l1qxZzJ49mw4dOnDCCesaYe677z6ef/55pkyZwjfffMNLL73Evffe26j3TiSNbRzbNonbyp83SSUlJQwYMID33nuPVatWce2113L44Ycza9Ysli5dyiGHHMI111zD8uXL2X777TniiCMqzq2p7AE49dRTWb58OdOmTWP58uXccssttebH3XnkkUdYsWIFr732GmPHjuWJJ56oSD/yyCPZdtttWbZsGddddx2HHXYYS5YsAYKK8a233srSpUuZNGkSb7/9Nv/3fxVDKRpU/opIkrl76BvBrfkvCKaPujJ+S0S8MLcxY8Z4TS677DI//vjjK57//e9/9zPPPLPi+bx58xzwn3/+2d3dTz/9dL/wwgsr0l9++WUfOnSou7vn5uZ6y5Yt/YcffqhIP+aYY/ziiy+uMQ+VPfPMM77VVltVPO/Tp4+//vrrFc8vv/xyP+KII6o998svv/SsrKyK5zvvvLPfe++9Fc8feOAB33HHHeuVH5GoATmexHKD4Ed/n2TGbOxWW1kXb/jw4f7000/7vffe6zvvvHPF/tzcXG/Tpo1PmzbN3Wsue6ZNm+YdOnTwVatW1Tludc455xw/++yz3d39hx9+8FatWvnq1asr0nfddVe/++67qz33pptu8j/+8Y8VeQ+j/JVwjRw50keOHBl1NpqUZJd3yd4S1aJ6HTCAYNT+kLitUfOqpqrge7L+4++++w6Ak046iY8++oj58+ezdu1aHnvsMfbff38AfvzxR1q0aMHQoevm595mm23q/Yv+/fffZ8sttwSCFtMFCxawzTYVg4FrvGb8uQBTp06t87kizZ2ZdTazx4EC4OfYvgPN7NpocxaeRYsW8eOPP7LllltWKR/at2/P4MGDmTp1aq1lz2effcbAgQMZM2YM3bt3Z/jw4TzzzDP1you788EHH1SUWVOnTmWTTTahQ4d1U8fWtbwLq/wVkcRK1GCqUcBQd19Q65FN3H777ceoUaM4/fTTGTJkCFdffTVmxtq1awEYMmQIAwYMoF+/fmRmZjJ8+HDGjh0LQG5uLh07dlzvep06dVqvv2lt3nzzTcaNG8ekSZMqrll+ndqu+c0333D11VfzwgsvVOzLzc2tcm5ubi7uTtwANREJ3AOsAAYSTMoP8AlwE3B5VJkKS3FxMUcffTTHH388w4YNIzc3lx49eqx3THn5UlvZM3fuXL777jsOPfRQ5s+fzyeffMIf/vAHtthiCzbffPM65ScnJ4eysrKK7kqVy6vymPHjAMo99NBDfPHFFzzwwAMV5za2/BWRxEtUi+pMQhpMkOr22msvrrrqKg499FAGDRrEoEGD6NChA/379weCPlCFhYUsW7aMvLw8DjnkkIoW1aysLFavXn+q19WrV1e0DgwePJisrKyKbdCgQesd++mnn3LUUUfx9NNPV7QKZGVlVVynumuW+/nnn9l///257bbb2G233Sr2V87T6tWrycrKUiVVpHq/B/4a+1HuEAyAIuj+1KSVlZVx7LHH0qpVq4of1zWVWbWVPW3btqVly5ZcfvnltGrVipEjR7LHHnvwxhtv8OOPP65X1mVlZXHssetPwz127FgeeeQRXnnlFVq3bl1rfuI9//zzXHrppUyYMIHu3bvX61wRiVaiKqqPAi+a2ZFmtmf8lqB4kTrrrLP46aefWLRoEYceeiglJSVstdVWAEyePJnRo0fTtWtXWrduzTnnnMNnn33G0qVLGTp0KCUlJfz0008V15oyZUrFrakZM2aQm5tbsc2aNaviuK+//poDDzyQhx56iN///vcV+7t06UKfPn2YMmVKtdcEmD17NnvttRdXXHFFlT8GW265ZY3nish6VgHd43eY2UYEK+Y1We7OSSedxKJFi3jmmWdo2bIlULV8yMvLY8aMGWy55Za1lj3VLdda/gN46NCh65V1ubm5PProoxXHPfTQQ9xwww28/fbbFY0A5fmZOXPmeq2glcus1157jVNOOYWXXnqJ4cOHV+yvrfwVkRSRiI6vwC8b2GZG3Sm3tm1DAwyKi4s9Pz/fL7nkEj/mmGM8Pz+/Yt+3337rZWVlPnv2bB85cqRfeumlFeeNHj3aDznkEF+5cqUXFRX5dddd53379q1IP+KII3zUqFGem5vrH374oXfs2NG/++67avNQ7ttvv/WePXv6E088UW36xRdf7LvvvrsvX77cp02b5r179/YJEya4u/vcuXN9k0028RtvvLHac++++24fNmyYz5071+fNm+dbbLHFBgcmiKQqkjS4ALgE+BjYg2Axkp2Bd4HzkhG/MVtNg6lOO+0033HHHX3NmjXr7V+8eLF37NjRn376ac/Pz/eLLrpovcGWNZU9RUVFPnjwYL/66qu9uLjYP/zwQ8/KyqoYiLUh48eP9169evn3339fbfqOO+7of//73z0/P9+fffZZ79Spky9evNjd3d9++23v2rWrv/fee9We25DyVxJLg6nqL1nlXVRb5BlItW1DhfeYMWOc4NZexTZmzBhfsWKFDx8+3Nu1a+e9evXySy65xEtKSirOW7p0qR911FHeo0cP79Spk//ud7/zSZMmVaQvW7bMDzroIG/Xrp0PGDDAH3vssWrjxxs9erSbmbdv375i22KLLSrSCwoK/IQTTvAOHTp4z549/aabbqpIy8nJcWC9c9u3b1+RXlZW5hdeeKF36dLFu3Tp4hdeeKGXlZXVmieRVJLEiqoRrBL1PZAHTAPOAywZ8RuzbaismzVrlgPeunXr9cqI8ePHu7v7m2++6Ztttpm3adPGR44c6b/88kvFuTWVPe7u3333ne+0007erl0733zzzf3ZZ5/dwCe4zqBBg7xFixbr5eW0006rSP/ll1985MiR3qZNGx86dKi/+eabFWnZ2dmemZm53rn77bdfRXpDyl9JLFVU6y/dK6oWvMZwmNmvwATgVeBNd88L7eJJkpOT4zk5OVFnQ0QawcyucvecqPORylTW1Z1WS0oevdf1l+7lXdh9VH8LTAKOBWaZ2Ztmdr6ZbRZyHBGRlBDrf3+/mb0S+/f3tZ8lIiJ1EWpF1d0XeLBU6qFAH+B6oB/wnJn9bGa3m9l+ZtY6zLgiIlEws78DTwDLgVeAZcDjsf0iItJIiZpHFXcvAd6JbReY2SDgD8A5wFbA/234bBGRJuFvwJ7u/l35DjN7FHiTYC5VERFphFArqmZ2EfAe8IW7l8anufssgiVV7wwzpohIxH6u9HwmsTlVRUSkccLuo7oX8BawyszeMrMrzGykbvWLSLows4zyDcgBHjSzIWbW1syGAvcBYyLNpIhImgi1RdXd9zGzTGA7YLfYdi7Q3sy+AN5z9ya/rKCINGslrGsxLV+y7chK+44CHkhyvkRE0k7ofVRjt/w/i203mVln4FSCvly7kAbrX4tIs7Zx1BkQEWkuQq+omlk3YPe4rQ/wCfBv4IOw44mIJJO7z446DyIizUXYg6m+B9oQDKj6ELjH3X8IM4aISKows07AX4Ftgaz4NHffJ5JMiYikkbBbVOcBWwMDgdlAfzP71d3XhhxHRCQVPAVkAs8B+RHnRUQk7YQ9mGrv2GCq3xDc9j8HeMzMZhPc9v/A3V8IM6aISIR2Arq7e1HUGREJk5YylVQR9vRUuHupu3/u7je5+5+BYcAzwDHAs2HHExGJ0IcEZZyIiCRAMgZTbQ3MB94m6LsqIpIuRgOvmtkkYFF8grtfHUmORETSSCIGU20G/AK8D9wOvO/uv4QZR0QkRVwHDABmAR3j9mtlKmk03X4XCb9F9WqCiun8kK8rIpKKRgFD3X1B1BkREUlHYQ+meiLM64mIpLiZQHHUmRBpKl6/PqvG9OVz8ut03L6X5oaWJ0ltofdRFRFpRh4FXjSzO6jaR/WdaLIkIpI+VFEVEWm4s2L//rPSfgc2SXJeRETSjiqqIiIN5O4bR50HEZF0Fvo8qhtiZhvFFgMQEREREalVMltUZwGzzOx6d78/iXFFRBLCzH5lA1NRuftGSc6OiEjaSWZFdWNgI4JFAERE0sExlZ73Ac4FNAOKiEgIklZRdffZwGzgg2TFFJHG0YTjNXP3KqvtmdlE4DXgtqRnSEQkzYS9MtVxtR3j7o+EGVNEJMUUEtxBEhGRRgq7RfWUDex3YHOgK6CKqkg9qWUzNZnZ1ZV2tQMOACZEkB0RkbQT9spUu1XeZ2ZbA9fEnl4SZjwRkYgNqPQ8D7iZYCEAERFppIT1UTWzIcDVwL7ArcAx7r4mUfFEpP7sqqtqPmDWrDod52PGhJSjpsXdT4g6DyIi6Sz0iqqZbQTkAIcC9wCbuvvysOOIiKQCM9sM2AZYb3Fyd38omhyJiKSPsAdTjQWOBR4Ghrj74jCvL5Ku1LLZNJnZP4ArgSnA2rgkB0KrqMbuUH0LPO3ux5hZNvBOpZhnufu4sGI2J69fn1Vj+vI5+XU6bt9Lc0PLk4gEwm5RPZOgj9YhwMFmVuUATYItImnkPOC37v5NguPcCXxead98d++f4LgiIpEKu6K6R8jXE5EonRBNF8wmNMtBPjA9kQHMbBSwEvgY2DSRsUREUk3Yo/6rTH4tIpLGrgDuMLMcYFF8gruXNfbiZtaRYFDqnsDJlZJ7mtkigtv/zwOXu3teY2OKiKSSZC6hKiINFVHLptTq4di/8ZVII+ijmhnC9a8BHnT3uZW6Uk0HRsT+HQiMI5gW67TKF4j1Z82uvP+II45g/vz5IWSx6Ws3eHSN6Zltn4sdd3CNx4X9fhYVFSXkulHGTtX3WlKXKqoiIg2XsBWozGwEsBewbeU0d18ILIw9/cXMLgJeppqKqrtPBCZW3p+TkzOmb9++4WW4Cft2xsM1ppfmB4Op1tZyXN/Dx4aUo0CrVq2C60bwOSUqdqq+15K6EjE9VRvA3D0/7GuLSHpIl1kO3H12Ai+fDQwC5sRaU7OATDPbwt1/UzkrQEYC8yIiEolQCzYzO55gypSJZla5P5WIiNTdfcBgglv8IwjmpX4F2NfM9jCzgRYYANwAvBBVRkVEEiXsFtVrgc2AlsD3wAMhX19EpFlw97XEzZNqZrlAgbsvMbNtgfFAF2AZ8BxwWSQZlSZJc8dKUxF2RfUn4CCCQQQzQ762iEiz5e45cY9vJhg8JSKS1sKuqB5OMAF2BvCXkK8tIpJSzOx2d/9rNftvdffzIsiSNCFq1RSpXdjzqC4FLg/zmiLSDDWd6bhGA1UqqgRLSZ+X1JyIiKShhE5PZWYZwOYE8wpOc/fSRMYTEUkGMzsx9rBF3ONymwBLk5wlkbRw4ylto86CpJhQK6pmdp+7nxp7PJBgXr+NCaZOmWtmf3T3GWHGFBGJwLGxf1vFPYagrFsEHJ/0HImIpKGwW1RHAafGHt9EMMn0NgQtqjcBNwKHhBxTRCSp3H0PADO71t3V3UlEJEHCrqjGr/G3CzCsfL1rM7sc+DnkeCIikXH3y82sC/AnoB8wD3jZ3ZdHmzMRkfQQ+komZraxmQ0GyoibAzD2uH0dzm9tZg+a2WwzW2Nmk81s/1jaIDNzM8uN266IO/dwM/vYzNaa2cSQX5qIyHrMbGdgBnA6sDXBEqY/x/aLiNQqOzub7OzsqLORssJuUW1P0Gpa3rK6E/Bh7PFwgtaGuuTpV2AkMAc4AHjSzIbHHdPZ3UuqOXc5cCswDNizvpkXEamnW4Ez3f2J8h1mdgRwO7BDVJkSEUkXYU9PVVMLbTFwRh2ukQfkxO162cx+AbYDvqzl3LcAtHyriCTJUODJSvueJljuVEREGimh01PFc/fvCZZVrRcz60Xwx2Bq3O7ZZubAm8CFsflb63vdbCC78v4jjjiC+fPn1/dyIo1y6tChoVynvt/d5hY3AX4iGET6eNy+vxB0BxARkUZKWkXVzFoA97l75TkHazqnJfAYMM7dp5tZFsHttMlAN+DOWPq+9c2Pu08kmJVgPTk5OWP69u1b38tV9C+ZOLHKJUVqdd+PP4ZynXuPPFJxk+s8grs+fwVmA4OAIcAfI8yTiEjaSFpFFcgkmFuwThXV2GIBjwJFwNkA7p4LfBE7ZJGZnQ0sMLMO7r4m/CyLiGyYu38cGzz6B6Av8BLwqkb9i4iEI+wJ/9+pITmzHtcx4EGgF3CAuxdv4FCP/Rv67AUiInXh7ivM7F2Ciup8VVJFRMITdovqjsD1wIJq0loCu9bxOncTLL26l7vnl+80sx2BlQT9wroQjKyd6O6rYumZsTgtgAwzawOU1lDRrTO76qqaD5g1q07H+Zgxjc2KiKQIM9uIoPvRTsAKoIuZfQoc4+6zI82ciEgaCLuiOhmY7u5PV04ws9bAXbVdILb06mlAIbAwaFyF2L4y4J9AT2A1wWCq+E5qxwL/iXueD4wDRtfvZYiI1Mk4gtlI9nP3vFg/+mti+7OjzJiISDoIu6J6K8FcptUpBk6o7QKxVgir4ZD/1nDuw8DDtcUQEQnJdsA+5Xdt3D3XzC4GlkWbLRGR9BD2PKpP1ZBWRtDKICKSLj4Ffgt8FLdve+CTaLIjIpJeEj7q38yucfcraj+yiTuh1sZiEUk/M4BXzewVghX1BhCspve4mV1dfpC7XxlR/kREmrRkTE91LpD+FVURaY7aAM/GHvck6Fv/HNCWoNIK62YnkSbqxlPaRp0FkWYrGRXVmvqbiog0We6uWykiIgmUjIqqWhNEJC2Z2RbAMndfFBvxfyHB7CQ3uvvaaHMnItL0JWOifLWoiki6+i/QOfb4/4DdCeZUvTeqDImIpJNktKj+OQkxRESiMMjdf4itpncIsAXB/M2/RJstEZH0kIwW1QwzO9/MdklCLBGRZCowsw4EU1TNcfelBAOq2kSbLRGR9BBqi6qZ/Rd4290fiD2/iGCVlm+Aa83sdHd/NMyYIiIRehx4B+gAjI3t+w1p3KKanZ0NwMSJEyPNhySWZjqQVBH2rf/fEUxHhZllEAwsOMrdnzGz/YEbAFVURSQtuPv5ZrYPUOzu78Z2lwHnR5gtEZG0EXZFtbO7L4493pbg9tfzseevUcPypyIiTZG7v1Hp+RdR5UXSi1o1k0d3ClJX2BXVpWY2yN1nAXsAn7h7aSytPVC6wTNFRJoYM9sYuA4YAWTFp7n7RlHkSUQknYRdUX0AeMXMXgeOA86JS9sdmBZyPBGRKD1OsIzq3wHNmyoiErJQK6ru/k8zmwdsD5zr7vG3+nsAN4UZT0QkYlsCv3P3sqgzIiKSjkKfR9XdxwHjNrBfRCSdvE/QH//LqDMiIpKOkjHhv4hIupoFvGZmzwEL4xPc/cpIciQiKeX167NqTF8+J79Ox+17aW5oeWpKVFEVEWm49sDLQEtgQMR5ERFJO6FXVM2sDWDunh/2tUVEUom7nxB1HkRE0lnYK1MdD5wGZJrZ/eUrVImIpCMz22RDae4+M5l5ERFJR2G3qF4LbEZwG+x7gumqRETS1c+AAxa3z2P/ZiY/OyIi6SXsiupPwEEEBbRaE0Qkrbl7RvxzM+sNjAE+iCZHIiLpJaP2Q+rlcIJ5BbcA/hLytUVEUpq7LwTOA64P87pmNsTMCsxsfNy+o8xstpnlmdnzZtY1zJgiIqkg1Iqquy9198vd/R+xAlskbWRnZ1esBy1Sg82AdiFf807g8/InZrYlcC9wLNCLYFWsu0KOKSISuYROT2VmGcDmBP23prl7aSLjiYgkk5l9wLo+qRBUULcErg4xxihgJfAxsGls99HAS+7+fuyYK4BpZtbB3deEFVtEJGphj/q/z91PjT0eSDC/4MYEBflcM/uju88IM6aISIQqDxjNA6a4+09hXNzMOhJUevcETo5L2pKg4gqAu88wsyJgKJVWyTKzbCC78rWPOOII5s+fX+88FRUVATTo3FTVbvDoUK5T3/ckqrhRxo7yNdckkd/r2l5zZtvnYscdXONx6fR/rj7CblEdBZwae3wTMBHYhqBF9SbgRuCQkGOKiEQiCUtDXwM86O5zzeInFiALWFXp2FVAh8oXcPeJBGXxenJycsb07du33hlq1aoVAA05N1V9O+PhUK7T9/CxTSJulLGjfM01SeT3urbXXJofTDu/tpbjwn7NTUXYFdX4knQXYJi7lwGY2eUEU7mIpCS76qqaD5g1q07H+ZgxIeVIUp2Z/Q14x90nm9lOwJNAKXC0u39c89m1XnsEsBewbTXJuUDHSvs6ArrtLyJpJRErU21MMEirjKCDf7m1BMsNioiki/OBB2OPrwduJqgs3gLs2MhrZwODgDmx1tQsgsVUtgBeI7hbBVQsPNAa+LGRMUVEUkrYFdX2BK2m5S2rOwEfxh4PB+aFHE9EJEqd3H2VmXUgqDju5e6lZnZTCNe+D3gi7vkFBBXXM4CewCdmthvwFUE/1mc1kEpE0k2oFdXKk19XUkxQwIo0TSdoWXep4lcz24VgcNP7sUpqR4Lb/43i7muJuytlZrlAgbsvAZaY2enAY0A34C1AX1ARSTsJnZ4qnrt/T7CsqohIurgQeBooAg6N7fsj8FnYgdw9p9Lzx4HHw44jIpJKklZRNbMWwH3ufmKyYoqIJJK7vwpUHib8VGxrkl6/PqvG9OVz8ut0HMC+l+aGkicRab7CXkK1JpnA8UmMJyKSdO5e7O7FUedDRCQdhD3h/zs1JGeGGUtERESkLsK6U6C7BMkX9q3/HQmmaFlQTVpLYNeQ44mIiIhImgq7ojoZmO7uT1dOMLPWwF0hxxMRERGRNBV2H9VbgeUbSCtG06eISBoxs8lmdp6Z9Yo6LyIi6SjUiqq7P+Xu1fZTdfeyJKyLLSKSTFcDuwMzzWyCmR1lZm2izlS6ys7OJjs7O+psiEgSJXPUv4hIWnH3Z939EGAA8AJwJrDQzB4ysz2jzZ2ISNOniqqISCO5+3JgHHAPMIdg8v/7zOxHM9sr0syJiDRhSZvwX0Qk3ZiZAfsAxxKsSPUJcAPwnLvnm9mhwHigd3S5FJFUduMpbaPOQkoLvUXVzDLMbE8zaxX2tUVEUswC4CbgG2ALd9/f3R9393wAd38GmBZlBkVEmrLQW1TdvczMXnD3DmFfW0QkxfzR3b+o6QB33yNZmRERSTeJuvX/vpnt5O6fJuj6IiKRMLNN4p4ur/S8grvPTFKWRETSVqIqqrOBCWb2AvAr4OUJ7n5lgmKKiCTDzwRlmtVwjKNlo0VEGi1RFdW2wPOxx/3j9nvVQ0VEmg5312wpIiJJkpCKqrtrBaokKZ/8euLEiZHmQ0RERCRsCZueysyGAEcC/YB5wH/d/adExRMRSQYze83d94s9/oAN3Cly992TmjERkTSUkIqqmf0JeAx4maC/6mbAF2Z2rLu/mIiYIiJJ8kjc4wciy4WISDOQqBbVfwIHufu75TvMLBsYC6iiWg921VU1HzBrVp2O8zFjQsqRSPPm7o/HPR4XZV5ERNJdoiqq/YEPKu37kPUHVomINHlm1gv4LdCduJkA3P2hyDIlIpImElVRnQz8HfhX3L6/xfaLiKQFM/szwRKpPwFbAlOBrQh+mKuiKiLSSImaZuVM4GQzm29mk8xsPnAqcEZdTjaz1mb2oJnNNrM1ZjbZzPaPpQ0yMzez3LjtikrnPmRmq81soZn9LSGvUEQErgVOcPdtgbzYv6cCX0abLRGR9BBai6qZbePuUwDcfZqZbQ7sBPQF5gOT3L24Hvn6FRgJzAEOAJ40s+Fxx3R295Jqzs0BhgADgd7Au2b2vbu/1oCXlfpO0ExgIhHayN2fqrRvHLAQuCCC/CTcjae0jToLItKMhHnr/wOgI4CZ/eTuQwhuf9Wbu+cRVDjLvWxmvwDbUXtLxfHAaHdfAawws/uB0UB6VlRFJEqLzayXuy8CZpnZzsBStCqViEgowqyorjSzPwLfA33MbGOqWWKwIetfxwYrDCXo/1Vutpk58CZwobsvNbMuQB9gStxxU4A/V3PNbCC78v4jjjiC+fPnV8nDqUOH1jfb1aru2jWJKm5zFOV73dy+X2n0vb4f2BV4BrgFeBcoA26KMlMiIukizIrqucCtBLfcM4AZ1RxT7/WvzawlwZys49x9upllATsQDMzqBtwZS98XyIqdtiruEquADlUy4j4RmFh5f05Ozpi+fftWycd9P/5Yn2xv0L1HHlmv46OK2xxF+V43t+9Xunyv3f1fcY8fMbOJQHt3nxZdrkRE0kdoFVV3fw54DsDM1rh7lcphfZlZBvAoUAScHYuTC3wRO2SRmZ0NLDCzDkBubH9HoCDu8ZrG5kVSh5aNlVTl7nOizoOI1J/6XqeuRE1P1a2xFzAzAx4EegEH1DAQq3z5wgx3X2FmC4BtCLoEEHs8tdozRUTqycx+ZQPLpsZz942SkB0RkbSWkIqquxeFcJm7gc2Bvdw9v3ynme0IrCSYt7ALcDsw0d3Lb/c/AlxuZl8QVHJPATQ0XkTCckzc4x0IBnDeTrBc9ECCuz+PVHOe1OL167NqTF8+J79Ox+17aW6N6SLSdCSqRbVRzGwgcBpQCCwMGlchtq+MYInWnsBqgpbT+I5qYwgqubOBfOBfaTs1lYgknbu/V/7YzO4E9nX3eXH7JhDMMqIBVSIijZSSFVV3n001MwbE+W8N5xYCJ8Y2aYLsqqtqPmDWrDod52PGhJQjkQ3qy7q+8eVygX4R5EVEJO0kamWqapnZH5IZT0QkwV4EXjSzvc1sczPbh2BQ6YsR50tEJC2EXlE1syFmdqiZbRO370Az+xL4T9jxREQidDrwCXAP8BVBt6NJsf0iItJIod76N7PRBBNgLwe6mdnfgD2BrQn6az0UZjxpprRsrKQIdy8ALoltIiISsrBbVC8GDnT3XsAhBJXTGcBm7n5n/Oh9EREREZGahF1R7evuE2KPXwJKgUtqmANVRERERKRaYVdUK0bqu7sDa1VJFRFpGDMbb2YLzGy1mf1oZifH9g8yMzez3LjtiqjzKyIStrCnp2pvZvFLCHaq9FyrtYiI1N31wEnuXmhmw4CJZvY1sCyW3tndS6LLnohIYoVdUd0z5OuJiKSs2FLPJxMsOtLd3bc2s92B3u7+ZGOv7+7xyz97bBvMuoqqiEhaC7WiGr9iS2VmlkmwatQGjxERaWKuBvYGbiWYogpgLnAL0OiKKoCZ3QWMBtoCXwOvAt1jybPNzAlW6LvQ3ZeGEVNEJFUkc2WqFsBlwJVJjCkikkijgW3dfamZ3R3b9wuwSVgB3P1MMzsH2BnIJlhaeimwAzAZ6AbcCTwG7Fv5fDPLjp23niOOOIL58+dXiddu8OhwMg7VXr8mtcXObPtc7LiDkxq3rppK3ChjN7e4UcdOB8leQrWmZVFFRJqaTNYtoeqxf7Oouqxqo7h7KfChmR0DnOHutwNfxJIXmdnZwAIz6+DuayqdOxGYWPmaOTk5Y/r27Vsl1rczHg4t330PH1uv42uLXZofzHC4tpbjwo5bV00lbpSxm1vcqGOng6Quocq6glxEJB28CtxsZq2hos/qNQTT8yVCC4I+qpWVl63JLtNFRBIq7JWpahpM1SrMWCIiKeBvwDhgFdCSoCX1DeC4xl7YzHoSDFB9GcgH9iIYtHWkme0IrAR+AroAtwMT3X1VY+OKiKSSsG/9P1hL+pxa0kVEmgx3Xw0cbGa9gI2AX919YViXB84gGKSVAcwGznP3F83sSOCfQE9gNcFgqiNDiisikjLCHvW/cZjXExFJZWZWfqt9SWzDzDLcvayx13b3JcDIDaT9F/hvY2OIiKS6hAymMrMsghGq3QkK70mVO/iLiKSBEqrpe29mJcB84FlgjLuHOrhKRKS5CL3jvZmdBywgNsgAmADMN7O/hR1LopWdnU12dnbU2RCJ0jnAO8A+wOYE00O9DVxEcNt+F4I5VkVEpAHCHkw1GrgEOAl4xt1LYxP9HwbcZmYr3P0/YcYUEYnQ34DfxA1i+tHMvgC+dPfBZvYt8GV02RMRadrCblE9Hzje3Z+MzfuHu5e6+/+A4wkKdRGRdNERaFdpXzugU+zxQoIVpUREpAHC7qM6GHhrA2lvE+JqLSIiKeAR4E0zuw34FegPnEswZRUEXQJ+iChvIiJNXtgV1TVAP6qfhqpfLF2aCLvqqpoPmDWrTsf5mDEh5Ugk5VxIMJfpKKAvQf/8O4H7Y+nvUs2qUNIwN56ixmmR5ibsiurzwF1mdpi7F5TvNLO2BIX3cyHHExGJTGwaqntiW3XpBdXtFxGRugm7onoJwS3+WWY2gaB1oQ+wPzAXOCbkeCIikYpN9v9bgun4rHy/uz8UWaZERNJE2BP+rzKznQmWD9wL2AFYClwGPOruRWHGExGJkpn9GRhPcPt/S2AqsBXwIaCKqohII4U+4b+7FxMspVrbcqrS1J1wQtQ5EInatcAJ7v5UbPq9bc3sBIJKq4iINFLoE/6LiDQjG7n7U5X2jSO4qyQiIo2kiqqISMMtjvVRhaBv/s4E0/RlRpgnEZG0oYqqiEjD3Q/sGnt8C8F0VFOAuyLLkYhIGgm9j6qISDNyY2yKKtz9ETObCLR392nRZktEJD0kpKJqZl2BC4ARQFZ8mrvvnoiYIiLJZGaZQK6ZdXb3QgB3r26xExERaaBEtag+DrQGngTWJiiGiEhk3L3UzH4EugHzo86PiEg6SlRFdRegR3krg4hImnoMeNnMbiNY1MTLE9z9nchyJSKSJhJVUf0G6A/MSND1RURSwRmxf3Mq7Xdgk+RmRUQk/SSqovoO8JqZ/QdYGJ+gZQVFJF24+8ZR50FEJJ0lqqK6G8FtsL0r7Xe0rKCIpBEzawnsBPR19/+ZWXsAd8+LNmciIk1fQiqq7r5HIq4rIpJKzGw48CJQSNDd6X/ASOB44IgIsyYikhYSPuG/BTLKt0THExFJoruBK919GFAc2/ce6xYBEBGRRkhIxdHM+pnZc2a2DCghKMDLNxGRdLElMD722KHiln/byHIkIpJGEtXCeQ9QBPweyAV+Q3B77PQExRMRicIsYLv4HWb2W+DnSHIjIpJmEjmP6kbunmdm7u5TzOwk4GOCtbFFRNLBFcArZnYP0MrMLiX4QX5KtNkSEUkPiWpRLSW45Q+w0sx6AHlAvwTFExFJOnd/GdgP6EHQN3UgcIi7vxFpxkRE0kSiWlQnAQcAzwGvE4yEzQe+SFA8EZGkM7Pu7v41cGbUeRERSUeJalE9lqB1AeA8ggUAvgOOSlA8EZEozDGzV83s6PL5U0VEJDyJmkd1ZdzjfODaRMQREYnYRsDhBEup3mNmLwOPAxPcvaTGM0VEpFaJmp6qtZldZ2YzzWxVbN8+ZnZ2IuKJiETB3Ze6+13uviuwFTAFuA5YEG3ORETSQ6Ju/d9CUGgfTWxuQWAqQauDiEg66gn0AroDK6PNiohIekjUYKqDgU1j01OVAbj7PDPTqH8RSRtmtgVwZGxrCzwJ/NndP4s0YyIiaSJRFdWiyteOTVG1LEHxRESi8BHwDHAa8K67lwGYWUb5YxERabhEVVSfAsaZ2fkAZtYHuBV4IkHxRESi0Mvdi8qfmNlw4HiCGU76RpYrEZE0kag+qv8AfgG+BToDPwHzgasSFE9EJOncvcjMepjZuWb2FTAZ2B44N9qciYikh0RNT1UEnA+cH7vlv9TdvZbTRESaBDNrCRwIjAb2BX4G/kuwMtXh7r44utyJiKSPUCuqZrbRBpIGmBkA7j4nzJgiIhFYBJQBDwNj3P0rADPTClUiIiEK+9b/LIJb/r/EHlfefqntArE5WB80s9lmtsbMJpvZ/tUcd6WZuZntFbevn5m9YGbLzWyumZ3e+JckIlLFNwTdmnYEdjCzLokIYmbjzWyBma02sx/N7OS4tN+b2XQzW2tm75rZwETkQUQkSmFXVKcQ9Ee9nOAWWMtKW6s6XKMF8CswEugUu9aTZjao/AAzGwz8haqTao8nqAz3Av4A/NPM9mj4yxERqcrds4HBwBvABcBCM3sJaE9Q1oXlemCQu3ck6GpwrZltZ2bdgWeBK4CuwBfA/0KMKyKSEkKtqLr7tsBhBAXnR8CrwCiglbuXuntpHa6R5+457j7L3cvc/WWCyud2cYfdCVxMMA0WAGaWBWQD17l7sbtPAZ4GTgzn1YmIrOPus939GncfAvye4IdzGTDFzP4dUoyp7l5Y/jS2DQYOAaa6+1PuXgDkANuY2bAw4oqIpApL1BgnM8sA9iYYbLA/sGd5P656XqcXMBsY4e7TzewvwDHufpCZzQJOdve3zKwDsJpgupjFsXPvB7aPVaArXzeboGJb2Q7A5/XNZ8wggi4OyRZV3ChjN7e4UcZuinH7u/vJtR8WLjNrQ7DgyXHuXqXLUgOveRdBOdoW+BrYnWCZ1lbufkbccd8R9Jd9ptL52aRPWRdl7OYWN8rYzS1uY2NHUt4ljbsnZAM2A/5J8Ma/C2zcgGu0BN4C7o0970DQtWBQ7PksYK+44z8E7gDaAL8BlgM/JOo1VpPfnGTFSoW4zfE1671O/7ipuAGZwK4EXaFaAg8CN1Q65iNgdLp/Ns3t+6j3Ov3jRh071bdQb/2bWVczO8vMPgOeB3KB3d19D3evdSBVpWtlAI8S3N4/O7Y7B3jU3Wdt4LSjgY0J+rjeTdBndW49X4aISErxoOvUh0B/4AyCsrVjpcM6AmuSnTcRkUQKex7V+QT9SR8FPo3t29TMNi0/wN3fqe0iFsxl9SDBoKgD3L04lvR7oH/cFDA9CAZa/cvd/+Xus4E/xl3ncUBrbotIumhB0Ed1KsEKWACYWfu4/SIiaSPsiupCgtvup8S2yhzYpA7XuRvYnOC2fn7c/t+z/ojaz4G/ARMAzGxzghbUQuBwYJ/YdUREmhQz6wnsCbwM5AN7AUfGtk+AG83sUOAV4ErgG3efHlF2RUQSItSKqrsPauw1YnMBnkZQ2VxYvlAAcJq7P1bp2FJghbvnxnbtC1wGtCMYdLCfuy9pbJ5ERCLgBLf57yGYoWU2cJ67vwgQq6SOJejiNIlghhURkbSSkCVUGyN2+95qPZCqFWN3vxW4NfRM1d3EZhY3ytjNLW6UsZtb3JQQ+5E9sob0t4CopqOaGFHcKGM3t7hRxm5ucaOOndISNj2ViIiIiEhjhL0ylYiIiIhIKFRRFREREZGUpIqqiIiIiKQkVVRDYGZnm9kXZlZoZg8nKWZrM3vQzGab2Rozm2xmoSzZWMf4481sgZmtNrMfzSypy7eZ2RAzKzCz8UmKNzEWLze2/ZCMuHHxR5nZNDPLM7MZZrZbguPlVtpKzeyORMaMiz3IzF41sxVmttDMxppZyg38bI6iKOticSMr75pbWReLGVl5l+yyLhZT5V0KU0U1HPOBa4GHkhizBcEKXCOBTgRLKz5pZoOSFP96gqVsOwIHAtea2XZJig1wJw1fp7yhznb3rNi2WbKCmtnewL+AEwiWEd4dmJnImHGvMwvoTTCP51OJjBnnLmAx0AcYQfAdP7OmEyRpoijrINryrjmWdRBBeRdFWQcq71KdKqohcPdn3f15YFkSY+a5e467z3L3Mnd/mWBVsKQUoO4+1d0Ly5/GtsHJiG1mo4CVwNvJiJcCrgKudvdPY5/1PHefl8T4hxIUpB8kKd7GwJPuXuDuC4HXgC2TFFtqEEVZF4sbWXmnsi6poi7rQOVdylFFNU2YWS9gKElcQtHM7jKztcB0YAHwahJidgSuJliRLNmuN7OlZvaRmWUnI6CZZQLbAz3M7Gczmxu7NdQ2GfFjjgce8eTNZXcrMMrM2plZP2B/gsJbBEh+edcMyzpIcnmXImUdqLxLOaqopgEzawk8BoxL5hKK7n4mwe2Z3YBnCVYTS7RrgAfdfW4SYsW7mGD5337AfcBLZpaMVpVeBMsGH0bwPo8AtiW49ZlwsZXiRgLjkhEv5n2CFoXVBEsifwE8n8T4ksKiKO+aWVkH0ZR3kZZ1oPIuVami2sSZWQbwKFAEnJ3s+O5e6u4fAv0JlntMGDMbQbDe+S2JjFMdd5/k7mvcvdDdxwEfAQckIXR+7N873H2Buy8Fbk5SbIBjgQ/d/ZdkBIt9n18jqAy0B7oDXQj6rUkzF2V511zKOoisvIu6rAOVdylJFdUmzMwMeJDgl+ih7l4cYXZakPh+W9nAIGCOmS0ELgAONbOvEhy3Ok4dl/ptVBD3FQS/suNvQyVzObnjSG7rQldgI2Bs7I/kMuA/JPePlaSgFCrvmltZB0ko71KgrAOVdylJFdUQmFkLM2sDZAKZZtYmSdNL3A1sDvzJ3fNrOzgsZtYzNoVIlpllmtm+wJEkvsP/fQR/IEbEtnuAV4B9ExnUzDqb2b7ln6uZHU0wGjVZ/Yj+A5wTe9+7AOcDLyc6qJntQnDrL1mjX4m1ovwCnBF7rzsT9Bn7Jll5kA2LsKyDCMq75lbWQeTlXSRlHai8S2nurq2RG5DDutGg5VtOgmMOjMUpAHLjtqOT8Hp7AO8RjEZdDXwLnBLR+z4+Sa/3c2BN7DV/CuydxNfZkmAKk5XAQuB2oE0S4t4LPBrB5zoCmAisAJYCTwK9kp0PbdV+Nkkv62JxIynvmltZF/eaIynvoirrYrFV3qXoZrE3SkREREQkpejWv4iIiIikJFVURURERCQlqaIqIiIiIilJFVURERERSUmqqIqIiIhISlJFVURERERSkiqq0qyZ2SAz8yROWi4iEgmVd9IUqaIqkTCzWWZWZGbdK+3/OlaQDoo9zzGz8XHpZmYXmtlPZpZvZnPM7Hoza53klyAiUicq70QaThVVidIvBMsRAmBmw4F2tZxzO3AqwZrMHYD9gd8TrOYhIpKqVN6JNIAqqhKlRwkK4HLHA49s6GAzGwKcSbBs4ifuXuLuU4FDgf3MbM8NnDcx1grxmZmtNrMXzKxrpcOOjrVWLDWzy+LO/a2ZfWJmK81sgZmNNbNWsTQzs1vMbHHsut+a2VaxtNZm9n+xay4ys3vMrG1D3iQRSQsq70QaQBVVidKnQEcz29zMMoFRwPgajv89MNfdP4vf6e6/xq61dw3nHgecCPQBSghaKuLtCmwWi3GlmW0e218KnA90B3aOpZ8ZS9sH2B0YCnQCDgeWxdJuiO0fAWwK9AOurCF/IpLeVN6JNIAqqhK18laGvYFpwLwaju0OLNhA2oJY+gbjuPt37p4HXAEcHvtjUe4qd8939ynAFGAbAHf/0t0/jbVmzALuBUbGzikmuB03DDB3n+buC8zMCG7Xne/uy919DfBPgj9MItJ8qbwTqSeN/JOoPQq8D2xMDbfBYpYStBBUpw9BH7AN+TXu8WygJesX9AvjHq8FsgDMbChwM7A9QX+yFsCXAO7+jpmNBe4EBprZs8AFQJvYsV8GZTgABsT/oRCR5kflnUg9qUVVIuXuswkK3AOAZ2s5/B1ggJn9Nn6nmQ0AdgLeruHcAXGPNyJoHVhahyzeDUwHhrh7R+AfBIVwef5vd/ftgC0Ibn1dGLtuPrClu3eObZ3cPasO8UQkTam8E6k/VVQlFZwE7Bm7TbVB7v4jcA/wmJntZGaZZrYl8Azwlru/VcPpx5jZFmbWDrgaeNrdS+uQtw7AaiDXzIYBZ5QnmNkOZrajmbUE8oACoMzdy4D7gVvMrGfs2H5mtm8d4olIelN5J1IPqqhK5Nx9hrt/UdMhcY/PBh4gGISQC7wGTCQYCVuTR4GHCW55tQH+WsfsXQAcBawhKIz/F5fWMbZvBcHttWXAjbG0i4GfgU/NbDXwFsHgBRFpxlTeidSPuXvtR4lExMxuBjLc/bxGXGMiMN7dHwgrXyIiYVN5J1KVWlQlZZlZZ2BfoKbWBxGRJk/lnUj1VFGVlGRmfwRmAJPQKiwiksZU3olsmG79i4iIiEhKUouqiIiIiKQkVVRFREREJCWpoioiIiIiKUkVVRERERFJSaqoioiIiEhKUkVVRERERFLS/wPDWbnn0JhWYQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x864 with 7 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(10, 12),\n",
    "                       constrained_layout=False)\n",
    "\n",
    "for i in range(3):\n",
    "    for j in range(2):\n",
    "        for axis in ['top', 'bottom', 'left', 'right']:\n",
    "            ax[i, j].spines[axis].set_linewidth(0.5)\n",
    "        ax[i, j].tick_params(length=6, width=0.5)\n",
    "\n",
    "        ax[i, j].set_xlim(0.5, 8.5)\n",
    "        ax[i, j].set_xticks(np.arange(1, 9))\n",
    "        ax[i, j].set_xticklabels(np.arange(1, 9).astype(int),\n",
    "                                 fontsize='large')\n",
    "\n",
    "        if i == 2:\n",
    "            ax[i, j].set_xlabel('MJO phase', fontsize='large')\n",
    "\n",
    "        ax[i, j].grid(color='0.8', linewidth=0.5, axis='y')\n",
    "\n",
    "ax1 = ax[1, 1].twinx()\n",
    "for axis in ['top', 'bottom', 'left', 'right']:\n",
    "    ax1.spines[axis].set_linewidth(0.5)\n",
    "ax1.tick_params(length=6, width=0.5)\n",
    "ax1.grid(color='0.8', linewidth=0.5, axis='y')\n",
    "\n",
    "ax[0, 0].set_ylim(220, 260)\n",
    "ax[0, 0].set_yticks(np.arange(220, 261, 10))\n",
    "ax[0, 0].set_yticklabels([f'{y:d}' for y in np.arange(220, 261, 10)],\n",
    "                         fontsize='large')\n",
    "ax[0, 0].set_ylabel('Ionospheric potential, kV', fontsize='large')\n",
    "\n",
    "ax[0, 1].set_ylim(-2, 2)\n",
    "ax[0, 1].set_yticks(np.arange(-2, 3, 1))\n",
    "ax[0, 1].set_yticklabels([f'−{-y:d}' if y < 0 else f'{y:d}'\n",
    "                          for y in np.arange(-2, 3, 1)],\n",
    "                         fontsize='large')\n",
    "ax[0, 1].set_ylabel('RMM1 and RMM2', fontsize='large')\n",
    "\n",
    "ax[1, 0].set_ylim(120, 160)\n",
    "ax[1, 0].set_yticks(np.arange(120, 161, 10))\n",
    "ax[1, 0].set_yticklabels([f'{y:d}' for y in np.arange(120, 161, 10)],\n",
    "                         fontsize='large')\n",
    "ax[1, 0].set_ylabel('Potential gradient, V/m', fontsize='large')\n",
    "\n",
    "ax[1, 1].set_ylim(220, 270)\n",
    "ax[1, 1].set_yticks(np.arange(220, 271, 10))\n",
    "ax[1, 1].set_yticklabels([f'{y:d}' for y in np.arange(220, 271, 10)],\n",
    "                         fontsize='large')\n",
    "ax[1, 1].set_ylabel('Ionospheric potential, kV', fontsize='large')\n",
    "ax1.set_ylim(120, 170)\n",
    "ax1.set_yticks(np.arange(120, 171, 10))\n",
    "ax1.set_yticklabels([f'{y:d}' for y in np.arange(120, 171, 10)],\n",
    "                      fontsize='large')\n",
    "ax1.set_ylabel('Potential gradient, V/m', fontsize='large',\n",
    "              rotation=270, va='bottom')\n",
    "\n",
    "ax[2, 0].set_ylim(249, 252)\n",
    "ax[2, 0].set_yticks(np.arange(249, 253, 1))\n",
    "ax[2, 0].set_yticklabels([f'{y:d}' for y in np.arange(249, 253, 1)],\n",
    "                         fontsize='large')\n",
    "ax[2, 0].set_ylabel('Mean OLR for 15° S–15° N, W/m²', fontsize='large')\n",
    "# thin spaces after '°'\n",
    "\n",
    "ax[2, 1].set_ylim(30, 50)\n",
    "ax[2, 1].set_yticks(np.arange(30, 51, 5))\n",
    "ax[2, 1].set_yticklabels([f'{y:d}' for y in np.arange(30, 51, 5)],\n",
    "                         fontsize='large')\n",
    "ax[2, 1].set_ylabel('Average daily sunspot number', fontsize='large')\n",
    "\n",
    "for j in range(2):\n",
    "    fig.align_ylabels([ax[i, j] for i in range(3)])\n",
    "\n",
    "for j in range(2):\n",
    "    ax[0, j].text(0.02, 0.98, '1980–2020', ha='left', va='top',\n",
    "                  transform=ax[0, j].transAxes, fontsize='large')\n",
    "    ax[1, j].text(0.02, 0.98, '2006–2020', ha='left', va='top',\n",
    "                  transform=ax[1, j].transAxes, fontsize='large')\n",
    "ax[2, 0].text(0.02, 0.98, '1980–2020', ha='left', va='top',\n",
    "                  transform=ax[2, 0].transAxes, fontsize='large')\n",
    "ax[2, 1].text(0.02, 0.98, '2006–2020', ha='left', va='top',\n",
    "                  transform=ax[2, 1].transAxes, fontsize='large')\n",
    "                    \n",
    "ax[0, 0].bar(np.arange(1, 9), phase_avg_ip / 1e3,\n",
    "             yerr=std_error(phase_avg_ip, phase_avg_sqr_ip, counter_ip) / 1e3,\n",
    "             width=0.6, color='coral')\n",
    "\n",
    "ax[0, 1].bar(np.arange(1, 9) - 0.2, phase_avg_rmm1,\n",
    "             width=0.4, color='peru', label='RMM1')\n",
    "ax[0, 1].bar(np.arange(1, 9) + 0.2, phase_avg_rmm2,\n",
    "             width=0.4, color='darkgreen', label='RMM2')\n",
    "\n",
    "ax[1, 0].bar(np.arange(1, 9), phase_avg_pg,\n",
    "             yerr=std_error(phase_avg_pg, phase_avg_sqr_pg, counter_pg),\n",
    "             width=0.6, color='royalblue')\n",
    "\n",
    "ax[1, 1].bar(np.arange(1, 9) - 0.2, phase_avg_ip_pg / 1e3,\n",
    "             yerr=std_error(phase_avg_ip_pg, phase_avg_sqr_ip_pg, counter_ip_pg) / 1e3,\n",
    "             width=0.4, color='coral', label='Ionosph. potential')\n",
    "ax1.bar(np.arange(1, 9) + 0.2, phase_avg_pg,\n",
    "        yerr=std_error(phase_avg_pg, phase_avg_sqr_pg, counter_pg),\n",
    "        width=0.4, color='royalblue', label='Potential gradient')\n",
    "\n",
    "ax[2, 0].bar(np.arange(1, 9), phase_avg_olr,\n",
    "             yerr=std_error(phase_avg_olr, phase_avg_sqr_olr, counter_olr),\n",
    "             width=0.6, color='teal')\n",
    "\n",
    "ax[2, 1].bar(np.arange(1, 9), phase_avg_sn,\n",
    "             yerr=std_error(phase_avg_sn, phase_avg_sqr_sn, counter_sn),\n",
    "             width=0.6, color='darkgoldenrod')\n",
    "\n",
    "leg = ax[0, 1].legend(fontsize='large', framealpha=1, edgecolor='0.',\n",
    "                      ncol=2, columnspacing=1., loc='lower right', handlelength=1.5)\n",
    "leg.get_frame().set_linewidth(0.5)\n",
    "\n",
    "handles, labels = [(a + b) for a, b in\n",
    "                   zip(ax[1, 1].get_legend_handles_labels(),\n",
    "                       ax1.get_legend_handles_labels())]\n",
    "leg = ax1.legend(handles, labels,\n",
    "                 fontsize='large', framealpha=1, edgecolor='0.', handlelength=1.5)\n",
    "leg.get_frame().set_linewidth(0.5)\n",
    "\n",
    "for p in range(1, 9):\n",
    "    ax[0, 0].annotate(str(counter_ip[p-1]), xy = (p, 221),\n",
    "                      ha='center', va='bottom', rotation=90,\n",
    "                      fontsize='large', color='0.')\n",
    "    ax[1, 0].annotate(str(counter_pg[p-1]), xy = (p, 121),\n",
    "                      ha='center', va='bottom', rotation=90,\n",
    "                      fontsize='large', color='1.')\n",
    "    ax[1, 1].annotate(str(counter_ip_pg[p-1]), xy = (p - 0.2, 221),\n",
    "                      ha='center', va='bottom', rotation=90,\n",
    "                      fontsize='medium', color='0.')\n",
    "    ax1.annotate(str(counter_pg[p-1]), xy = (p + 0.2, 121),\n",
    "                 ha='center', va='bottom', rotation=90,\n",
    "                 fontsize='medium', color='1.')\n",
    "\n",
    "for i in range(3):\n",
    "    for j in range(2):\n",
    "        ax[i, j].text(-0.25, 1.05, chr(ord('a') + 2 * i + j), fontsize='x-large',\n",
    "                      fontweight='semibold', ha='left', va='bottom',\n",
    "                      transform=ax[i, j].transAxes)\n",
    "\n",
    "fig.subplots_adjust(hspace=0.25, wspace=0.3)\n",
    "\n",
    "fig.savefig('figures/ip_pg_rmm_variation.eps', bbox_inches='tight')\n",
    "fig.savefig('figures/ip_pg_rmm_variation.png', dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c169197e-495c-475c-9a97-3bc3ca7494b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Largest positive IP deviation from the long-term mean: 12251.515085229592.\n",
      "Largest negative IP deviation from the long-term mean: -11782.171890650236.\n"
     ]
    }
   ],
   "source": [
    "print('Largest positive IP deviation from the long-term mean: '\n",
    "      f'{np.amax(phase_avg_ip) - 240e3}.')\n",
    "print('Largest negative IP deviation from the long-term mean: '\n",
    "      f'{np.amin(phase_avg_ip) - 240e3}.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63796d95-d604-49d4-b4de-1174787e526b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(10, 10.5),\n",
    "                       constrained_layout=False)\n",
    "\n",
    "for i in range(3):\n",
    "    for j in range(2):\n",
    "        for axis in ['top', 'bottom', 'left', 'right']:\n",
    "            ax[i, j].spines[axis].set_linewidth(0.5)\n",
    "        ax[i, j].tick_params(length=6, width=0.5)\n",
    "\n",
    "        ax[i, j].set_xlim(0.5, 8.5)\n",
    "        ax[i, j].set_xticks(np.arange(1, 9))\n",
    "        ax[i, j].set_xticklabels(np.arange(1, 9).astype(int), fontsize='large')\n",
    "\n",
    "        if i == 2:\n",
    "            ax[i, j].set_xlabel('MJO phase', fontsize='large')\n",
    "\n",
    "        ax[i, j].grid(color='0.8', linewidth=0.5, axis='y')\n",
    "\n",
    "ax[0, 0].set_ylim(210, 260)\n",
    "ax[0, 0].set_yticks(np.arange(210, 261, 10))\n",
    "ax[0, 0].set_yticklabels([f'{y:d}' for y in np.arange(210, 261, 10)],\n",
    "                         fontsize='large')\n",
    "\n",
    "ax[0, 1].set_ylim(200, 250)\n",
    "ax[0, 1].set_yticks(np.arange(200, 251, 10))\n",
    "ax[0, 1].set_yticklabels([f'{y:d}' for y in np.arange(200, 251, 10)],\n",
    "                         fontsize='large')\n",
    "\n",
    "ax[1, 0].set_ylim(220, 270)\n",
    "ax[1, 0].set_yticks(np.arange(220, 271, 10))\n",
    "ax[1, 0].set_yticklabels([f'{y:d}' for y in np.arange(220, 271, 10)],\n",
    "                         fontsize='large')\n",
    "\n",
    "ax[1, 1].set_ylim(220, 270)\n",
    "ax[1, 1].set_yticks(np.arange(220, 271, 10))\n",
    "ax[1, 1].set_yticklabels([f'{y:d}' for y in np.arange(220, 271, 10)],\n",
    "                         fontsize='large')\n",
    "\n",
    "ax[2, 0].set_ylim(120, 170)\n",
    "ax[2, 0].set_yticks(np.arange(120, 171, 10))\n",
    "ax[2, 0].set_yticklabels([f'{y:d}' for y in np.arange(120, 171, 10)],\n",
    "                         fontsize='large')\n",
    "\n",
    "ax[2, 1].set_ylim(110, 160)\n",
    "ax[2, 1].set_yticks(np.arange(110, 161, 10))\n",
    "ax[2, 1].set_yticklabels([f'{y:d}' for y in np.arange(110, 161, 10)],\n",
    "                         fontsize='large')\n",
    "\n",
    "for i in range(2):\n",
    "    ax[i, 0].set_ylabel('Ionospheric potential, kV', fontsize='large')\n",
    "ax[2, 0].set_ylabel('Potential gradient, V/m', fontsize='large')\n",
    "\n",
    "for j in range(2):\n",
    "    fig.align_ylabels([ax[i, j] for i in range(3)])\n",
    "\n",
    "for k in range(4):\n",
    "    ax[k // 2, k % 2].text(0.02, 0.98, f'{1981 + 10*k}–{1990 + 10*k}',\n",
    "                           ha='left', va='top',\n",
    "                           transform=ax[k // 2, k % 2].transAxes,\n",
    "                           fontsize='large')\n",
    "ax[2, 0].text(0.02, 0.98, '2006–2015', ha='left', va='top',\n",
    "              transform=ax[2, 0].transAxes, fontsize='large')\n",
    "ax[2, 1].text(0.02, 0.98, '2011–2020', ha='left', va='top',\n",
    "              transform=ax[2, 1].transAxes, fontsize='large')\n",
    "\n",
    "for k in range(4):\n",
    "    ax[k // 2, k % 2].bar(np.arange(1, 9), phase_avg_ip_part[k] / 1e3,\n",
    "                          yerr=std_error(phase_avg_ip_part[k],\n",
    "                                         phase_avg_sqr_ip_part[k],\n",
    "                                         counter_ip_part[k]) / 1e3,\n",
    "                          width=0.6, color='coral')\n",
    "\n",
    "for j in range(2):\n",
    "    ax[2, j].bar(np.arange(1, 9), phase_avg_pg_part[j],\n",
    "                          yerr=std_error(phase_avg_pg_part[j],\n",
    "                                         phase_avg_sqr_pg_part[j],\n",
    "                                         counter_pg_part[j]),\n",
    "                          width=0.6, color='royalblue')\n",
    "\n",
    "for p in range(1, 9):\n",
    "    ax[0, 0].annotate(str(counter_ip_part[0, p-1]), xy = (p, 211),\n",
    "                      ha='center', va='bottom', rotation=90,\n",
    "                      fontsize='large', color='0.')\n",
    "    ax[0, 1].annotate(str(counter_ip_part[1, p-1]), xy = (p, 201),\n",
    "                      ha='center', va='bottom', rotation=90,\n",
    "                      fontsize='large', color='0.')\n",
    "    ax[1, 0].annotate(str(counter_ip_part[2, p-1]), xy = (p, 221),\n",
    "                      ha='center', va='bottom', rotation=90,\n",
    "                      fontsize='large', color='0.')\n",
    "    ax[1, 1].annotate(str(counter_ip_part[3, p-1]), xy = (p, 221),\n",
    "                      ha='center', va='bottom', rotation=90,\n",
    "                      fontsize='large', color='0.')\n",
    "    ax[2, 0].annotate(str(counter_pg_part[0, p-1]), xy = (p, 121),\n",
    "                      ha='center', va='bottom', rotation=90,\n",
    "                      fontsize='large', color='1.')\n",
    "    ax[2, 1].annotate(str(counter_pg_part[1, p-1]), xy = (p, 111),\n",
    "                      ha='center', va='bottom', rotation=90,\n",
    "                      fontsize='large', color='1.')\n",
    "\n",
    "for i in range(3):\n",
    "    ax[i, 0].text(-0.25, 1.05, chr(ord('a') + 2*i), fontsize='x-large',\n",
    "                  fontweight='semibold', ha='left', va='bottom',\n",
    "                  transform=ax[i, 0].transAxes)\n",
    "    ax[i, 1].text(-0.18, 1.05, chr(ord('a') + 2*i + 1), fontsize='x-large',\n",
    "                  fontweight='semibold', ha='left', va='bottom',\n",
    "                  transform=ax[i, 1].transAxes)\n",
    "\n",
    "fig.subplots_adjust(hspace=0.25, wspace=0.25)\n",
    "\n",
    "fig.savefig('figures/ip_pg_variation_partial.eps', bbox_inches='tight')\n",
    "fig.savefig('figures/ip_pg_variation_partial.png', dpi=300, bbox_inches='tight')"
   ]
  }
 ],
 "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}