From c06b26b12a847981a99d3c36f0808adb380fdaf9 Mon Sep 17 00:00:00 2001 From: Gao Date: Thu, 18 May 2023 16:09:20 +0200 Subject: [PATCH] prepare to add fft --- 20230509_Data_Analysis.ipynb | 2491 ++++++++++++++++++++++++++++++++++ DataContainer/ReadData.py | 138 +- 2 files changed, 2628 insertions(+), 1 deletion(-) create mode 100644 20230509_Data_Analysis.ipynb diff --git a/20230509_Data_Analysis.ipynb b/20230509_Data_Analysis.ipynb new file mode 100644 index 0000000..eef41a0 --- /dev/null +++ b/20230509_Data_Analysis.ipynb @@ -0,0 +1,2491 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Import supporting package" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "import numpy as np\n", + "import copy\n", + "\n", + "from uncertainties import ufloat\n", + "from uncertainties import unumpy as unp\n", + "from uncertainties import umath\n", + "\n", + "import matplotlib.pyplot as plt\n", + "plt.rcParams['font.size'] = 18\n", + "\n", + "from DataContainer.ReadData import read_hdf5_file, read_hdf5_global, read_hdf5_run_time\n", + "from Analyser.ImagingAnalyser import ImageAnalyser\n", + "from Analyser.FitAnalyser import FitAnalyser\n", + "from ToolFunction.ToolFunction import *\n", + "\n", + "from ToolFunction.HomeMadeXarrayFunction import errorbar, dataarray_plot_errorbar\n", + "xr.plot.dataarray_plot.errorbar = errorbar\n", + "xr.plot.accessor.DataArrayPlotAccessor.errorbar = dataarray_plot_errorbar\n", + "\n", + "imageAnalyser = ImageAnalyser()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Start a client for parallel computing" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "
\n", + "
\n", + "

Client

\n", + "

Client-1ad26782-f4a3-11ed-8014-80e82ce2fa8e

\n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "
Connection method: Cluster objectCluster type: distributed.LocalCluster
\n", + " Dashboard: http://127.0.0.1:8787/status\n", + "
\n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "

Cluster Info

\n", + "
\n", + "
\n", + "
\n", + "
\n", + "

LocalCluster

\n", + "

39ebed45

\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + "\n", + " \n", + "
\n", + " Dashboard: http://127.0.0.1:8787/status\n", + " \n", + " Workers: 14\n", + "
\n", + " Total threads: 56\n", + " \n", + " Total memory: 782.31 GiB\n", + "
Status: runningUsing processes: True
\n", + "\n", + "
\n", + " \n", + "

Scheduler Info

\n", + "
\n", + "\n", + "
\n", + "
\n", + "
\n", + "
\n", + "

Scheduler

\n", + "

Scheduler-5010ee37-2d62-418f-908a-3d869cf85668

\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
\n", + " Comm: tcp://127.0.0.1:54873\n", + " \n", + " Workers: 14\n", + "
\n", + " Dashboard: http://127.0.0.1:8787/status\n", + " \n", + " Total threads: 56\n", + "
\n", + " Started: Just now\n", + " \n", + " Total memory: 782.31 GiB\n", + "
\n", + "
\n", + "
\n", + "\n", + "
\n", + " \n", + "

Workers

\n", + "
\n", + "\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 0

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54920\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54937/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54877\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-6k5w8crs\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 1

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54966\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54967/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54878\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-_7zysh1d\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 2

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54952\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54955/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54879\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-pivw5um_\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 3

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54940\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54943/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54880\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-4mhsz13r\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 4

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54972\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54973/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54881\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-7_6k4mgb\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 5

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54961\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54964/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54882\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-_vtk580o\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 6

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54960\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54962/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54883\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-ofu5zf80\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 7

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54939\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54941/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54884\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-u45d30hc\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 8

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54953\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54958/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54885\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-old86hcm\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 9

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54975\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54976/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54886\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-0xsdxrw4\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 10

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54946\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54947/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54887\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-hnelt761\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 11

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54969\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54970/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54888\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-ywnxkrum\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 12

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54945\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54948/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54889\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-07tq1zpz\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "

Worker: 13

\n", + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n", + " \n", + "\n", + " \n", + "\n", + "
\n", + " Comm: tcp://127.0.0.1:54951\n", + " \n", + " Total threads: 4\n", + "
\n", + " Dashboard: http://127.0.0.1:54954/status\n", + " \n", + " Memory: 55.88 GiB\n", + "
\n", + " Nanny: tcp://127.0.0.1:54890\n", + "
\n", + " Local directory: C:\\Users\\data\\AppData\\Local\\Temp\\dask-worker-space\\worker-lfu5si46\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "\n", + "
\n", + "
\n", + "\n", + "
\n", + "
\n", + "
\n", + "
\n", + " \n", + "\n", + "
\n", + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from dask.distributed import Client\n", + "client = Client(n_workers=14, threads_per_worker=4, processes=True, memory_limit='60GB')\n", + "client" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Set global path for experiment" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "groupList = [\n", + " \"images/MOT_3D_Camera/in_situ_absorption\",\n", + " \"images/ODT_1_Axis_Camera/in_situ_absorption\",\n", + " \"images/ODT_2_Axis_Camera/in_situ_absorption\",\n", + "]\n", + "\n", + "dskey = {\n", + " \"images/MOT_3D_Camera/in_situ_absorption\": \"camera_0\",\n", + " \"images/ODT_1_Axis_Camera/in_situ_absorption\": \"camera_1\",\n", + " \"images/ODT_2_Axis_Camera/in_situ_absorption\": \"camera_2\",\n", + "}\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "img_dir = '//DyLabNAS/Data/'\n", + "SequenceName = \"Evaporative_Cooling\" + \"/\"\n", + "folderPath = img_dir + SequenceName + \"2023/05/09\" # get_date()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Check the stability of our BEC" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The detected scaning axes and values are: \n", + "\n", + "{'runs': array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.,\n", + " 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.,\n", + " 22., 23., 24., 25., 26., 27., 28., 29., 30., 31., 32.,\n", + " 33., 34., 35., 36., 37., 38., 39., 40., 41., 42., 43.,\n", + " 44., 45., 46., 47., 48., 49., 50., 51., 52., 53., 54.,\n", + " 55., 56., 57., 58., 59., 60., 61., 62., 63., 64., 65.,\n", + " 66., 67., 68., 69., 70., 71., 72., 73., 74., 75., 76.,\n", + " 77., 78., 79., 80., 81., 82., 83., 84., 85., 86., 87.,\n", + " 88., 89., 90., 91., 92., 93., 94., 95., 96., 97., 98.,\n", + " 99., 100., 101., 102., 103., 104., 105., 106., 107., 108., 109.,\n", + " 110., 111., 112., 113., 114., 115., 116., 117., 118., 119., 120.,\n", + " 121., 122., 123., 124., 125., 126., 127., 128., 129., 130., 131.,\n", + " 132., 133., 134., 135., 136., 137., 138., 139., 140., 141., 142.,\n", + " 143., 144., 145., 146., 147., 148., 149., 150., 151., 152., 153.,\n", + " 154., 155., 156., 157., 158., 159., 160., 161., 162., 163., 164.,\n", + " 165., 166., 167., 168., 169., 170., 171., 172., 173., 174., 175.,\n", + " 176., 177., 178., 179., 180., 181., 182., 183., 184., 185., 186.,\n", + " 187., 188., 189., 190., 191., 192., 193., 194., 195., 196., 197.,\n", + " 198., 199., 200., 201., 202., 203., 204., 205., 206., 207., 208.,\n", + " 209., 210., 211., 212., 213., 214., 215., 216., 217., 218., 219.,\n", + " 220., 221., 222., 223., 224., 225., 226., 227., 228., 229., 230.,\n", + " 231., 232., 233., 234., 235., 236., 237., 238., 239., 240., 241.,\n", + " 242., 243., 244., 245., 246., 247., 248., 249., 250., 251., 252.,\n", + " 253., 254., 255., 256., 257., 258., 259., 260., 261., 262., 263.,\n", + " 264., 265., 266., 267., 268., 269., 270., 271., 272., 273., 274.,\n", + " 275., 276., 277., 278., 279., 280., 281., 282., 283., 284., 285.,\n", + " 286., 287., 288., 289., 290., 291., 292., 293., 294., 295., 296.,\n", + " 297., 298., 299., 300., 301., 302., 303., 304., 305., 306., 307.,\n", + " 308., 309., 310., 311., 312., 313., 314., 315., 316., 317., 318.,\n", + " 319., 320., 321., 322., 323., 324., 325., 326., 327., 328., 329.,\n", + " 330., 331., 332., 333., 334., 335., 336., 337., 338., 339., 340.,\n", + " 341., 342., 343., 344., 345., 346., 347., 348., 349., 350., 351.,\n", + " 352., 353., 354., 355., 356., 357., 358., 359., 360., 361., 362.,\n", + " 363., 364., 365., 366., 367., 368., 369., 370., 371., 372., 373.,\n", + " 374., 375., 376., 377., 378., 379., 380., 381., 382., 383., 384.,\n", + " 385., 386., 387., 388., 389., 390., 391., 392., 393., 394., 395.,\n", + " 396., 397., 398., 399., 400., 401., 402., 403., 404., 405., 406.,\n", + " 407., 408., 409., 410., 411., 412., 413., 414., 415., 416., 417.,\n", + " 418., 419., 420., 421., 422., 423., 424., 425., 426., 427., 428.,\n", + " 429., 430., 431., 432., 433., 434., 435., 436., 437., 438., 439.,\n", + " 440., 441., 442., 443., 444., 445., 446., 447., 448., 449., 450.,\n", + " 451., 452., 453., 454., 455., 456., 457., 458., 459., 460., 461.,\n", + " 462., 463., 464., 465., 466., 467., 468., 469., 470., 471., 472.,\n", + " 473., 474., 475., 476., 477., 478., 479., 480., 481., 482., 483.,\n", + " 484., 485., 486., 487., 488., 489., 490., 491., 492., 493., 494.,\n", + " 495., 496., 497., 498., 499., 500., 501., 502., 503., 504., 505.,\n", + " 506., 507., 508., 509., 510., 511., 512., 513., 514., 515., 516.,\n", + " 517., 518., 519., 520., 521., 522., 523., 524., 525., 526., 527.,\n", + " 528., 529., 530., 531., 532., 533., 534., 535., 536., 537., 538.,\n", + " 539., 540., 541., 542., 543., 544., 545., 546., 547., 548., 549.])}\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAG+CAYAAABCjQqZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAACRpUlEQVR4nO29e3wV1bn//5lsQkgIgSREwexotFUu6ulF23ohFX6nRsV60AAKSq3aiqU9rVxErVIgtlWpKOjx1uLtUAS5JF7qldADR9BqxfqtR1SqgggREQLkTi47z++PcXb2ZS5rzW3Pzn7er9e8kj0za2bNmjVrPetZz/MshYgIDMMwDMMwjOdkpToDDMMwDMMwmQILXgzDMAzDMD7BghfDMAzDMIxPsODFMAzDMAzjEyx4MQzDMAzD+AQLXgzDMAzDMD7BghfDMAzDMIxPsODFMAzDMAzjE/1SnQHGGT09Pfj8888xaNAgKIqS6uwwDMMwTMZBRGhubsYxxxyDrCxznRYLXmnO559/jrKyslRng2EYhmEynt27dyMcDpuew4JXmjNo0CAA6ssuKChw7bpdXV1Yv349KisrkZ2d7dp1GUYErn9MKuH6x8jS1NSEsrKyaJ9sBgteaY42vVhQUOC64JWXl4eCggJueBjf4frHpBKuf4xdREx+2LieYRiGYRjGJ1jwYhiGYRiG8QkWvBiGYRiGYXyCBS+GYRiGYRifYMGLYRiGYRjGJ1jwYhiGYRiG8QkWvBiGYRiGYXyCBS+GYRiGYRif4ACqDJMCIhFg82Zg715g+HCgogIIhVKdK4ZhGMZrWPBiGJ+prQWuvx7Ys6d3XzgM3HsvUFWVunwxDMMw3sNTjQzjI7W1wKRJ8UIXANTXq/tra1OTL4ZhGMYfWPBiGJ+IRFRNF1HyMW3fzJnqeQzDMEzfhAUvhvGJzZuTNV2xEAG7d6vnMQzDMH0TFrwYxif27nX3PIZhGCb9YMGLYXxi+HB3z2MYhmHSDxa8GMYnKipU70VF0T+uKEBZmXoewzAM0zdhwYthfCIUUkNGAMnCl/Z76VKO58UwDNOXYcGLYXykqgpYtw4oLY3fHw6r+zmOF8MwTN+GA6gyjM9UVQETJnDkeoZhmEyEBS+GSQGhEDB2bKpzwTAMw/gNTzUyDMMwDMP4BAteDMMwDMMwPsGCF8MwDMMwjE+w4MUwDMMwDOMTLHgxDMMwDMP4BHs1MgzD9EEiEQ5ZwjBBhAUvhmGYPkZtLXD99cCePb37wmF15QQO0sswqYWnGhmGYfoQtbXApEnxQhcA1Ner+2trU5MvhmFUWPBiGIbpI0QiqqaLKPmYtm/mTPU8hmFSAwteDMMwfYTNm5M1XbEQAbt3q+cxDJMaWPBiGIbpI+zd6+55DMO4DwteDMMwfYThw909j2EY92HBi2EYpo9QUaF6LyqK/nFFAcrK1PMYhkkNLHgxDMP0EUIhNWQEkCx8ab+XLuV4XgyTSljwYhiG6UNUVQHr1gGlpfH7w2F1f7rF8YpEgE2bgFWr1L/skcmkOxxAlWEYpo9RVQVMmJD+kes5ECzTF2HBi2EYpg8SCgFjx6Y6F/bRAsEmxiTTAsGmo/aOYQCeamQYhmECBgeCZfoyLHgxDMMwgYIDwTJ9GRa8GIZhmEDBgWCZvgwLXgzDMEyg4ECwTF+GBS+GYRgmUHAgWKYvw4IXwzAMEyg4ECzTl2HBi2EYhgkcfS0QLMNocBwvhmEYJpD0lUCwDBNL2mi82tra8NJLL+F3v/sdqqqqcNxxx0FRFCiKgoULFwpdY9++fZgzZw5GjBiB3NxcFBUVoaKiAo888ghIL2BMAp988gmuu+46HH/88RgwYACOOuoonHfeeaipqRG6/z/+8Q9MmzYN4XAYOTk5GD58OC655BL8z//8j1B6hmGYTEMLBDt1qvqXhS4m3Ukbjdff//53jB8/3nb6t99+G+eddx4aGhoAAPn5+WhubsaWLVuwZcsWrF27Fs899xxycnJ007/44ouYPHky2traAAAFBQVoaGjA+vXrsX79elx99dV49NFHoRhYgz7yyCOYMWMGuru7AQCDBw/Gvn378Mwzz+CZZ57BggULhAVIhmEYhmHSk7TReAFAYWEh/v3f/x1z587FqlWrMGzYMKF0jY2N+OEPf4iGhgaMHDkSb731Fpqbm9Ha2or7778f2dnZWL9+PWbNmqWbfufOnbj00kvR1taGs88+G9u3b0djYyMaGxsxf/58AMDjjz+Ou+66Szf93/72N/zsZz9Dd3c3Lr74YuzevRuHDx/G/v37cd111wEAqqursWbNGhulwjAMwzBM2kBpQnd3d9K+4447jgDQggULTNPOmzePAFBubi7t2LEj6fjtt99OACgUCtH27duTjk+bNo0A0LBhw+jQoUNJx6dPn04AqKCggA4ePJh0fMyYMQSATj31VOrs7Ew6ft555xEAOu6443Sf04zGxkYCQI2NjVLprOjs7KRnnnlGN78M4zVc/5hUwvWPkUWmL04bjVfIwcT+8uXLAQBTpkzB8ccfn3T8l7/8JfLz8xGJRPDkk0/GHWttbY3acM2YMQNDhgxJSv/rX/8aANDU1IRnnnkm7tiOHTuwZcsWAMANN9yA7Oxsw/S7du3Cq6++KvdwDMMwDMOkDWkjeNll+/bt+OyzzwAAF1xwge45+fn5qPgqEt/69evjjm3ZsgXt7e2m6cvLyzFq1Cjd9HV1ddH/zz//fN30Y8aMwaBBg3TTMwzDMAzTd+jzgtd7770X/f+UU04xPE879v777xumP/nkky3Tb9u2TTf9UUcdhaOOOko3bSgUwsiRI3XTMwzDMAzTd+jzgtfnn38e/b80MRJfDNqxpqYmtLS0JKUvLCxEXl6eZfrY+8X+Nru3WXqGYRiGYfoOaRNOwi7Nzc3R/80Ep9hjzc3NyM/Pj0tvljb2eOz93EifSEdHBzo6OqK/m5qaAABdXV3o6uoyTSuDdi03r8kwonD9s0ckAmzZokSDjY4ZQxz3ygZc/xhZZOpKnxe8+hp33HEHqqurk/avX7/eUrizQ6yNGsP4Ddc/cf72t+F45JFT0dCQG91XXNyOn/70/3DmmXtTmDNrIhHg/feLcejQABQWHsHo0Q2BEBi5/jGiaDE+RejzgpdmtA6oBVNQUKB7XmyhxabR/rcqVO14bFo30ify61//GrNnz47+bmpqQllZGSorKw2fzQ5dXV2oq6vDueeeq+uJyTBewvVPjqefVvCHP4SQuADHwYMD8Ic/fAdPPRXBJZdYr86RCp5+WsHs2SHU1/cGny4tJdxzT+ryzPWPkUWbfRKhzwtexxxzTPT/+vp6Q+Gkvr4egBqRXptmjE1/6NAhtLW1GWqVtPSx94v9rR03wih9Ijk5ObrR9bOzsz1pILy6LsOIwPXPmkgEmDMHSUIXABApUBTghhv6YeLE4C23U1sLTJmSnPfPP1cwZUq/lC+GzfWPEUWmnvR54/pYT8ZYD8VEtGOjR482TG/mcailT/R81NJ/+eWX2L9/v27aSCSCDz/8UDd9phKJAJs2AatWqX8jkVTniGGCyebNwJ49xseJgN271fOCRCQCXH+9kcCo/p05k799pu/R5wWvESNG4NhjjwUAvPzyy7rntLa2YvNXrVJlZWXcsTFjxiA3N9c0/a5du/DBBx/opj/33HOj/xulf+2116JG9YnpM5HaWqC8HBg3Drj8cvVvebm6n3EfFnLTm72C5lui5/lFugqMDOOUPi94AcCVV14JAHjqqafw6aefJh1/4IEH0NLSglAohCuuuCLu2MCBAzFx4kQAwEMPPYTGxsak9IsWLQKg2mddfPHFccdOOOEEjBkzBgBw991363o+3HnnnQCA4447Dt///vflHq6PUVsLTJqU3CDX16v7WfhyFxZy05/hw909zy/SVWBkGKekleB16NAhHDhwILr19PQAUA3TY/fHxuEC1KV6hg0bhra2Nlx44YV4++23AQCdnZ146KGH8Jvf/AYAMH36dJx00klJ973tttswcOBA7N27FxdddBE++ugjAKqm7LbbbsPDDz8MAJg3bx4KCwuT0v/hD39AKBTCP//5T0yZMiVqz3Xw4EH8/Oc/x0svvRR3XqZiNfVApB5njYw7sJDbN6ioAMJhQFH0jysKUFamnhck0lVgZBjH+LB2pGtoi2JbbT/+8Y+T0m7dupWKi4uj5wwaNIiys7OjvysrK+nIkSOG937hhRcoLy8vev7gwYMpFApFf1911VXU09NjmH7ZsmXUr1+/6PlDhgwhRVGiv60W+jaiLy2SvXGjJl6Zb9XVvmWpz9LdTRQOG5exohCVlannpQJepFiOmhr1nSlK8ntUFPV40NDqYGKeg1AHg17/urvV9nLlSvVvqr5Tppc+uUi2U0477TRs27YNs2bNwoknnoiuri4MHDgQY8aMwbJly/DSSy/pegtqjB8/Hu+++y6uvfZalJeXo729HUOGDMG5556LdevW4fHHH4diNOQE8NOf/hRvvvkmLr/8cpSWlqKtrQ1HHXUULr74Yvz1r3/FwoULPXjq9EJ0SmHBAtbGOIXta/oWVVXAunVA4gIZ4TBS7hloRCgE3Huv+n9i06n9Xro0eJ6YqcbIPGDtWrbVTBcUIr2JHSZdaGpqwuDBg9HY2Oh6HK8XX3wR48ePt3STjUTUDlqLll1RYa+x3LRJbUREKCsDdu7kRtkuq1apjbYVK1cCU6d6n59EZOof04tb36Kf1NaqJgSxA4GyMlXoSpXAGNT6p5kHiPTa4TBwzz1ASUl61Qcv8OO7kOmL+3wcL8Zb9BrNcFgdyco2mpqtipkmRkPTxowdK3cPRoXta/omoVD6fRNVVcCECeknMPqNmQ2sHnv2AJdeGr/PbtuczrjZR7lFxkw1Mu7jtnF27NSDCOztZJ90Nchm+iaawDh1qvpXE7o41EkvVuYBImSa40xQHYhY8GJs4VXww6oqQGcpSl1YG2Mftq9hgo6TUCfpJrCJ5NeNgWYmBaYNcoBeFrwYW3hpnH3rrao2xgjWxrhDOhpkM5mBE01FusWmE82vWwPNTHGcCbIDEQtejC28DH6oaWMUxb42Jt1GvKmiqgr49FNg40bVkH7jRtVpoa8IXVwP0g8nmoqgTi0ZIZNfK/MAWfq6qUaQA/Sy4MXYwmvjbCfamHQb8aYaI/uadIfrQXpiV1MR5KklPWTza2YeYIe+bqoRZAciFrwYW/hhnG1HG5NuI17GG7gepC92NRVBnlrSw05+jQakMmSKqUaQHYhY8GJs4Zdxtow2Jt1GvIw3cD1Ib+xqKoI8taSH3fzqDUjXrjW3i9XIJMeZIDsQseDF2CZoxtnpNuJlvIHrQXojYstUUqJqL2Pt9oI8taSHk/wmDkgnTRITxjLNcSZofZQGB1BlHBGk4IfpNuJlvIHrQXqjaSomTVKFLz3N5f79wLRp6v9aMMwJE9T/6+v10yiKejwoU2yagOlWfvWC515ySTDaZhG8ii4fpD5KgwUvxjFBiZadbiNexhtk6kGqlthJx6V9/ETTVCRGHNdDs9tbt85YYEv11JIeZgKmW/lNddssWs+9ji6f6nJIhKcamT5DkI0pGf8QrQcHDqTG65G9LcWItWVasUKdXtQj1m5vwoRgTi0ZEdSpMBGsQrWI1vNMdIRhwYvpM4RCqs2D2Vpmbo94OU5U8BAxqp0yRV3Hzu/GPhM7GSdomorSUnV60YhYu710i02XbvkFrIUq0XqeqY4wLHgxrhAEAaS2Fli82Pj4DTe425ix5iK4mGkSVq9W66nfjX2mdjJuIGu3l26x6dIpv1ZC1dq15vWcCPjZz4DOzsx1hGHBi3FMEAQQs05N46mn3OvUWHMRfIw0CSUlqWnsM7WTcQO23wwGIoOHX/zC2i5v/351EPTss2L37WuOMCx4MY4wEkD27AEmTgRuu01M2HGqMbPq1AD3OjXWXKQPepqEVHk9Zrq3pZNvnO03g4HI4MFsSjiW/ftV0w8R+ppAzYIXY4pZYymiZVqwwFr75YbGzM9OjTUX6U2qtCeZrLWx843Htj2bNwNLlqj7gxYMM5PwYlAQCmWeQM2CF2PI3/42HF//ej/DxlJEywSo52jTb4mC3Lp17kzZ+dmpZbrmIt1JlfYkU7U2dqbl9QS1WbNUO8109ADsK4i2nyUl4utJRiLqYDWTBGoWvBhdnn5awaJF30F9ffz+2MZSVrCYPj25MZ0yxZ0pO7vRru2QyZqLvkCqlhIJ8hImXmFnWt5MUFu8GLjnnvTyAOxLiA4eHnxQ7rozZ2aYQE1MWtPY2EgAqLGx0bVrdncTlZb2ENBDvX4ovZuiEJWVEa1fn3zMi23jRrF819SoeVMU62uGw+r5dssnHDa+j1Y+3d32rs8QdXZ20jPPPEOdnZ2e3aOmRn2Pse+urMx+vQj6fVPBxo1y37j2bRmd59e35Uf9Sze6u9X3NHNm77tIfDeK0luPa2qIhg4Vf//a9Veu7P2dTsj0xRy5nkli82agvt5YdaTZMPk1EhHVrNmNdi37HH5EnGa8J1VLiQRxCROvkJ2Wl7GfDFIk8r6OXmT5rKx4TWU4rLZ7WntaVQX88IfqfiOD+9hlkYIWXd5LWPBikhBtLFtavM2HhsyUXWynVl+v2oXoffSaTYEW7Vq20zMS8hIbHybYpKqxz5RORnZanu0n7eHlElTa1G/idLEmdGltqN49+/cHHn5YTQ+4N0hN9yW32MaLSSJItkl2FrW1E+3aDukYcVqGIATFZdIbWYcCtp+Ux8s4ilae64oC1NSYCz5WwYyLiuTamCDEjXQKa7yYJCoqgNJS+sqwXtA1xSOuvdb+SMaP0bNfmgu/R3h6UwtDhwLTpvWOboH0HnUy3iM7La8JavX1+p197NQUY6yNcmJKEYtbU7960+sHDqgzEjILY3v9vH7BGi8miVAIuOcedeihKCZBugQZOhQoLraX9sQT7d+3r4ye/R7hGXmVHTigdpLjxgFHH61u6TzqZPxBZiHoTPT8tIsfgZzdHLzGBjM+eFB+rdS+FLiaBS9Gl0suIdx001uuCCWXXw786U/20jq5v8g0RzisfqhBnU6zuzSR3WlCkaC4ANDQoG4yeWLsk+7TvjLT8jKCWibjRyBnLwavdgUomecN/Pfig5cl4yFehJMg6nWnfuWVLldCQtTUED31lPj5brmN19SY36e4WD7MhF9uz3Zd6/XCFYiGzxB1//f6vWWKO79IXXLyPtOZVIYXSIf6t3Kl2De5cqX9e3gROkc2xIjs886cmZrvRaYvZo0XY8q+fc6voXkPDh0ql86PKYVErY22xuTatfrn+zntZ2dE63TxbqfeYm6MsjMFkbqUyYux662zKYqexiPwWhBJ/DCl8GLq1+70pehzLF0a/O+FBS/GFDemGrXOeNMmsfOLi92ZUtBU2naYOlXNQyx+d4KyDZQbNhBu2bv1JXd/LzpskbrUl2xa/ERPoO2LNomaKYURoktQWdVvt6d+7QqMIquTGAmAgftevFW+MV7j9VRje3unqapZZps3T+y8DRvceQan02ZAr3o6FRG1ZVXyMucbTeNYTS2IbtXVzp49KFM9XkzzWdUlgKioiGjxYntTMumA29OIiVHVRbbESOuxBKX+mVFTk2wqIfJsidcQrd9678zOe3QyfWm0OolMe+XV9yLTF8ObLDB+4bXg1dnZGa3sToWYDRu8W2pHrwEQtQkw27T82LFLcNq5yDZQbtlAuPG+RRp9M4LQ8RmVg2inZoQbA4LYzYkNTypwW5jVu55MPdVrc4JQ/8yw+kaLi8WELif128l7NBOgrO5ttOSWqNDt1ffCglcG4YfgReRe4+bkgzPCqAGornanY9OEJ5mP2q3ORaa8nHToeuus2X3fZh2aKKnu+LzUcLoxIEisn+mC28KsW4PCxDJMdf0zQ0RjGg6b102n9duN9+hkzVK9Qa1do323YMErg/BL8CKKV+eXlMRXZk3lLSIguLlIsFkDoOXLacOsfdyiH7UXnYtIeYlME4ZC4o2t2fv2WihIdcdntxEX0XK6pfFKt8XY3RZmRQQQmW88llTXPzPcEDCcXMPN9+jmlLMXHpgysOCVQfgpeMWi98HICFRufHAiDYAmeDkRvrT8iXzUHR3eaEpEy8srG4iODnvCl121fqo7Pjuu+qJaTjfs6JxOd6YCN2wQ7VzPTp1Pdf0zw40wEk6ukWrNkhlezKiIwuEkGM/Rc/WWCZLoxFVcQyTcQkMDsHBhskeOKMXFvUvhiLhVv/66N0ENRcvLzANp5kyxe9XXJ3s5vf66+bqXRgR9VQAjZD2vZDxeY+uSXdIxmKiop+uzz4qFbHHDc1bU8y9IuBFGwsk1gryQeboE3+W1GhlX8WvtQkD8wz7xRFUg1NYJ++gjVRgjsk7b0KB2BFVVvR914hqG4bAqdFVVqcKKm3m3g966aBUV6u+lS63Tz5gBNDf3/g6HVeFBFk1odYPEtSrPOksVBvfuBY46Sj3nyy/dWzNSZs1Aq7APWhy7CRN681VVpS4QPHWquHv7kiVqSAS7z+j3ep+JyMRhSkRvLT6nQn26LkHkxnqWTq4hWu4ffSR2nhF266tR+xeod+yd4o3xg1RNNfqB1XSDXZW3mRu23pRO4tSgWb6CrIa3O8Vld0rMSUiJRK/axCk8M1s1t6JUi05b2H3nMlNlTm1TghD93gsbRNH6XFyc/M2b2ZV61f65ZdNkZVJQXe3MLMFsWq67m6i01LrMrQz8rZ4v1fVVFrbxyiD6quAl8uHZMaa06wUlKiil2sDTCqPGVmQLhcTTFRc7e0at/q1e3WVLUHTLnkPEbtGuvYyMd6OTZ/EqLIaTvLhlg2hVn2fOtBdzyov2T68ulZaKCUmi19MTMM0EFruOTqIe43YGmEGqrzKw4JVBBMm43i1kPjyZUZsTLygZI/FUGniKsHatfS9F0U7S6TOqmq5nqLS0x3Ye3RJwvdK8iqZzojlMReBfK9yOw+Sml7SG24KX6IBPVqsTWzerq+0JLHbacq/WiQxifRWFBa8MIhWCl5dqYDsfnmjD68QLSnbk5kVn4AZO43PpBV/14hk7Ozvpt7/dbDufdt+bHexqOUWmypxM1xAFd+rb7ThMbg8E3RS8ZAZ8dgdnfgssXtWroNZXEWT6YjauZ6TQvLeI4vfrGb/aQWZhaM2IX9SY0o5Bu4ihqh5BNPA0encyTJgALF7c+1yaYfsXX6hejyUlQFGRahjr9FkPHRrg7ALwx7MqFFIN3ydPTj5mZsCteTdOmqSeF/tetHT33uusHIPqgabnhOPE4NtPpx5ZrNq0WIj0nTGc3kOv3XSCGwb+egS1vroNC16MMHa8t0SuGSuc1NeLpUv88EQaXlkvKKdeT0HpDLRFcK+91pnQpbndJz5XbS1w883Jnp733utMCC8sPGI/8Vf4Ec6ithaYNUv/WKzHqx4inrJOcCP0gF+ICKLp5oEIyAsJsUKS5o1sNXgTvYcWKsboeqKehF69q3Sqr47wQQPHeIifU41uq4H1pr2GDvVO1Szr1ZeKqUFtymTFCqIlS9S/Tr2f3IjubTT94ZUhbKyNlx1HAL9sQaxsd9asEbuOVzaTQXf20CMI0/RuTjXaNXGwWlPVzj0SbTsT12iVNSFx+12lY33VYBuvDMJPwcuOQaVRh2LXu9Dph+eGG7ZXmAlJTtZ5dCp0GS2466VdSaJXo8xz+OXEkC6GwEF39tDDS+cdEbyw8XLjW9Te2Zo18eWjrZhh1wN47lz7Ayi335Vb9dXvOsSCVwYRZI2X0QhqzRoxLYxXHUUQRtR6ebJqNGUbHaearqIiVRA1arC8NIR1EsfLr3eZTobAemVYUhIfboHpxSuvRlHByKx+6x0Ph3uFJ7uhYszanaB4vop+16mIA8aCVwbh9yLZompgq8WrRbZEtbibHWqqR9SJeREVREUbQNnpDUVR87Bhg3iZeOVSTpRc/xLfV0dH7+8NG+TybYZMvfDy+b3AbNHzVAenDNL3SORfHC8n7WNiOk1zJRocWnYTGUC4vei1nWulKg4YC14mbNmyhS677DIKh8PUv39/Gjp0KI0bN46eeOIJ6unpMU37xRdf0OzZs+mkk06iAQMGUGFhIY0ZM4aWLVtmmZaI6OOPP6bp06dTeXk55eTkUElJCVVWVtK6detsP4/f4SRE1MBuaFuAXvsmow8vaI21XWSFJJEGUCY4p90GSTTfS5bIv6NUBPCVHSW7veizHwQxOGUQo5T7Ebm+utpZPDO9dxgOi0WVt7NZDSCC8B5TOf3PgpcBt956KwGIbkOGDKH+/ftHf59//vl05MgR3bRbt26l4uLi6Ln5+fnUr1+/6O/KykrDtEREL7zwAuXl5UXPLygooKysrOjvq6++Wkh4SyQocbxitVFO4mWJChhB+MjdQkZIEmkAieSXo7FTbnaWgBF9R2b1r7tb1XDNm6duGzbEa8DsCDZ2BBJRDfDatc7rqqzgpnd+EG3SgigIEvXWv/b2Tk8FZifxzPzerNrjILzHVE7/s+Clw7Jly0gTcqZMmUK7d+8mIqKOjg5asWIFDRo0iADQtddem5T28OHDNGzYMAJAI0eOpLfeeiua9v7776fs7GwCQDNmzNC9944dO2jgwIEEgM4++2zavn07ERE1NzfT/Pnzo/latGiR9HMFMXK9rCCht5WUGHv0BeUjl8GN9R1lGg0RoaioSBVanK4BKGNXIvqOzDSuelMpWVn2BRsnAomVBtiJ0XLsPWQEN6PzvVzmxQ5BFAQ12ts7acqU96moqEe43N3CTWN80c1sOTCr9yDyHvXMGLzQAqdy+p8FrwS6u7ujgtO3v/1tXc3S448/TgAoKyuL3n333bhj8+bNIwCUm5tLO3bsSEp7++23EwAKhUJRoSqWadOmEQAaNmwYHTp0KOn49OnTSdOCHTx4UOrZgrhWo4wgIbuMRpAbayP0OsJYw3UvbLy0+/rhzSZr/C7yHHr1r6ZGrl6JPqPTUbKRBtjKiUSkHGQHGW7YVvplkxZU54SaGkoSuLz6dszyYNdQXmZLHCDYaSvsaOhk15QURcb8we0+ggWvBN544w3StEpPPvmk7jk9PT109NFHEwCaM2dO3LFjjz2WtOlAPZqbmyk/P58A0Pz58+OOtbS0UG5uLgGgaoNF13bu3BnN32OPPSb1bEEUvJxMwVg1dkFtrI2w8lbUQjW47dUYe38/vNliR69Lljh/R3rG9bK2K6KCqugoeeZM42ePnfp85ZXe307KQXaQ4ZZtpZffTmw9ES0fP50TeoV74zVC/RrcyQ5o7GyxZgZ2PQndmOHQytWpUCujLXRbe8mCVwKrV68mTbB55513DM/7/ve/TwDo5JNPju778MMPo2nXmERDvOCCCwgAnXHGGXH7X3755Wj6v//974bpR40aRdo0qAxBFLyIxLUtsQFDzRZu1hq7FSuC11gbIaPJ0oQvo/O1BtCOet5vbzY31P2J9c/LdTZlrq2nYXLaOa5cqS+8LV4s93xObYO8FijsBvN1WxA0+oZkBVc/8pW4b+1a43ZVNN8TJ7rrsOSmTZobdVBUW+i29pIFrwRiBa+tW7cannf22WcTAOrXrx91dHQQEdG6deuiad9//33DtHPnziVtujCWxYsXR9O3trYapp88eTIBoFNPPVXq2YIqeBHJjaBkVMSpaBTtIGvgHtvQ6kWud+JQIDtl5cT+wg2tZGL9ExW49TYjAU9U6DfqFNwKUFtd7SwEgPZ8sp6sXnZCidgpKy8EQbNvSFaAcHNwJ/NtG7WrXtjwibQDXtikbdhgsyBNysjrOsaCVwJvvvkmacLPE088oXtOV1cXDR06NHrenj17iIjovvvui+4zK9ClS5dGz2tubo7unz17NgGgwsJC0zzOnDmTAFBxcbHUswVZ8CIS78BFO40VK9JnSQlZFbxXXkOyU1Z6jVZpqXhUf5l4b0bE1r+1a4kGD3bWiCfWQbvaF+09uTGtpyjuxFyS1XhVVydP23ppNG6nrLwQBK2+IdkwDm4N7qyE0tggxmYDMxEByMxpSS9fMsKgmzZpRUXG7160T+nu9negLtMXZ8Qi2aeddhqGDRuGL774AosWLcIVV1yBfv3iH/2Pf/wjDhw4EP3d1NSE0tJSNDc3R/fl5eUZ3iP2WHNzM/Lz86P/W6WNPR57Pz06OjrQ0dERl08A6OrqQldXl2laGbRruXHNs8/u/b+nR90SKSlRILJm+9FHd+Puu4EpU0JfLc6qRI8pCgEAFi+OoKeHdO/jJ6LPpLF7dze6uihpfyQC/OpX/UAEAErcMSL1ua+/Hhg/vlt3Udr//V8Fe/YY54NIXZR348ZuHDyolm3iverrCQsW9P4uLSXcc08El1ySnF8AuPtuxdE70urdjTcC6phG0T/RBEUhFBUBP/4xUF/fm76oiHDwYPSsmBRi99m9uxsbN8K0TK0hEAFHouuA23u+0lLgjDO60dUFnHEGUFraD59/Hl/mieefdFIEQCjunkSE7u6Ibv1zilX906O0lHD33RFcdBHBjWZN5BtasSL5mBGhEOGLL5yXl1m+NBYsAP70J8Jll/Vg9eqsuLqsfYc9PWo+jL47tW4D+/crmDYtPq3eN/z004phOzBpEvDUU/HpLroIeOopBbNnh+LyZ5eDB/Xv8/TTyfcwe47iYrE22KjtlUGmr8wIwSsUCmHhwoX42c9+hg8++AAXXnghbr/9dpx66qk4dOgQVq5ciV//+tfIzs6OFl5WVlaKc63PHXfcgerq6qT969evtxTu7FBXV2d6PBIB3n+/GIcODUBh4RGMHt0gvSK9dp3i4ko0NAyAfgNEGDq0HU1NdcjJAW68cTgeeeRUNDTkRs8oLm7HT37yHnJy9uLFF+Xz4DbWzxTPrl1v4MUXG5L2/9//FaO+foxhOiIFe/YAixe/iVNPTU7/6qulAE63vP8LL/w//PnPJ4MovlNWif9dXw9cdlkIN930Fs48c2/Stdx4R6+9NhxLl4o2UYlCkyrYNCQXh4HQpfdbn1273sD/+38DIFKmRmRlEXp6stDaavcK6vNdccVbeOWV3vKfNm04Fi36DozK4/TTP8bUqV9PuprV+5Qltl3YvTsfwEjLNJMnb0dZWXNcO+LWdyzyDR04ABQUdKCpqT+s6kIkAkyd6ry8rPKlUV8P3HNPcp+U+N6Mvjuja+q980gE+PnPK3XbAVWYI/ziF53o168urq3PyQHuu6/3vQ8efAT33XeacPsXjwKi+Pu89tpw3HXXd4SfAwB27SoGYF2+Rm2vDG1tbeInO1ewpQ8333wzAdDdjj76aFqwYEH09759+4goeFONR44cocbGxui2e/duAkAHDhygzs5O17bW1lZ65plnqLW11fCc1au7qLQ03vuntLSHVq/usnVPdUHkHlKU+Gtq+xKv297eSXV1XbR8eRfV1XVRe7t7z+/WtVev7iLVQ8rMS6qHwuEew3ssX94lpC5fvly/3OvqxNIvXtwtNR1glW8n5djU1EoFBUeE85KVlVwP1ZAAxuUuu8U+r2iZJm5VVd2W9UFkC4eNvzO97zIc7qFVq7T9RqESrN+naJ1PvL/IVldnr90Q2US/oV/+slu3DfKqvETzpW7W702rm0880UWLF3fTY4910dChcu9ctG6LvC+jNl3mG6ir66JVq7qSvnGRd9He3kmlpcbv060639nZSQcOHCArOUEDlmf0Md544w36yU9+QqeeeiqVlZXRaaedRrfccgvt378/Go+rsLAwGuuLjev1bby8CmIahAWs7Rixm9kdGAX9FC0vp8bqojZXdg3YvXBkkBFsbrklOXL9hg32nsW4YU/2xpUxKNaCSLoR7kEkBpGTiOh2jK81m6Prr7dXtl7aZcra+tTUkJTg6KT+u+kRqLcE0dCh8s/gdhBSvfa0uJgoP1/sPjK2d3rvwq94hmxcb5Pzzz+fAND48eOj+zichP6SLV4GMU3lunZWhq568a9EBLXubrVhLCqKP09EqHTDWF2k8bHbCXgRukNGE6DX2DqNLSSyQLus27qo15nZFgqpQqaGWWiExP1+dKiymx/elDJ5XLtWTdfe3kmTJ3/gef1PRZR6q2fwUkCPrY+igyMRj2Ord+HHgJ4FLxvs2rWLQqEQAaDVq1fHHdMCqF5zzTW6aVtaWoQCqN5222266T/99FPShLN0CKCabkFMRZHxvtIEK7/CNLgxarNqfOx2AqnUeJWU6JehXSFSE2JF134UieOllbFbgSa18jYS+OfO9X7JILdCaXipzbaTR20A09nZSb/97WahNHpes17n0+0t9p27MdATQeQ+MkKXVd31ekDPgpcknZ2ddN555xEAOuWUU5KEDW3JoLy8PNq5c2dS+kWLFhFgvWTQ8OHD6fDhw0nHZ8yYQQBo0KBBabFkUCrXw/ISO0sdmYUDcHsKxY1Rm1XjI+MW7uUUUXt7JxUXt5GVHYimodB7TlkhUkaIjS3HDRvi16EzEtrcmlaaN08+9EFsfXXaoToNpTFvnvfabCd53LhRbf9qap75yj7IuLyKi50vgE6kTW+6Uz9k64XeO/dzuTGz+8jUc6O1Z82+VTfrHwteOnzyySd066230ttvv03t7e1EpK7huGnTJhozZgwBoPz8fPrHP/6RlDZ2kezRo0dHg7B2dHTQgw8+SP379ydAbJHsiooK+te//kVEqjasurqaFNXHPrCLZCd21qIq4nTTeLmlkfCyHGRGbU60a1adltdTRJ2dnXTTTW+aGjnPnWt+jbVrjfMOJAvNokKs3UC2qZ5W0gQFpx2qUwHSzvcgW5ed5HHlyt72TzUO1y8vt78NzRxB9p06qQ9m+TRbbN1NwcVsQGnnPcZ+i1ZtmZvx61jw0uGdd94hbTpPURQqLCykfv36Rfcdc8wx9Prrrxum37p1KxUXF0fPHzRoEGVnZ0d/V1ZW0pEjRwzTv/DCC5SXlxc9f/DgwdGpTQB01VVX6S7ebYXXgtfq1V26wTTdGDkHDTcNXRMbcrcR0VzJCAeJ14vV2ugZ7Xrt8GBW/0pK1AWozVi71niaQlvAOnaJHr2Rsh5OnUrMRviA+Xfl1rZggZgNmxF2Byh22wU7gq6TQdSSJarGVRt4Gt3fK223kSCiN4UsE7HezjuPbRf02gFZwUXGJlHbb1dzrS38LXKuG20ZC146HDp0iObPn0/f//736ZhjjqH+/ftTcXExnXnmmbRo0aK4EBBGfPHFFzRr1iw68cQTacCAATRkyBAaM2YMLVu2jCKRiGX6jz/+mK699loqLy+P3v/cc8+ldevW2X4uLwUvI41DbGU2GjmvWZM6A3kzzAQWN6KR621ua/6MOgKtzI3U80aNjKhzgJvv0+p6ZhpXq3vPnWv+PubMcaaxctrZmo3w3Y4ArrclOngMHWotyMZiZ4Bidxq3utqeoOt0EFVa2kM33fRm1NTCb62/jIAiapMlareohxte7Ha9xe06pYiumeqWooAFrwzCK8HLysbGyL7BaGTm5ZIkooh8+DU17nZy4bD7a805nV4QWW/Qy6lEs/egdSzLl3fRb3+7mVpaOqU6izVr7JeLlTOEm8uPWIUf8WIAIPvsZnmXzZ/MmqKi1zbrMJ1O66oDTjVWmh5u27k6tUPy0ibL7oDDjgBtpWHzcnM6QGbBK4PwSvAS9SpL9OhZu9b/jlwEUQGju1s8vgxgbDMTe9xu/C+9c91qiLR7eRkWRPY9aOUVuz8xaKJZB97R4WxNR0CdjtHCNdgVgNyYWo6tF/Pmed/pyBjWm2lVjbbSUvFpWDt12ex69gcqanBNJ16zoh6ibtgheRUyQfRZlyyJ7wdkBWiZNJMmuf8NOP1uWfDKILwSvETjKMVW1lR05CLI5Et2imLmTHVkZidAqozqvbubaPFi9xoZbaFctzqPxLya2WvI5TVxFQPj6VLRYJFWW0mJmH2IW+VlhVe2h7J516uvWVli1zX7DrSYTjIDntjNrMM0EkbWrHGmwXQr5IKMsCk7TStiQyWCV05Hdjc7YSac1n0RWPDKIFKt8YqtrDIjIz+FLxkBw24jY9YBJTbCInYLsTY3ZpHv7W4lJeJaC61jE2m8zYRJtwQI0elSvzerztZu5+eV7aHZu9Z7r2ZlPHCgvfJxY1rVqsM0Knen04VOp/dk32txsXvG+jKmH34K/jJbSYl5nQyFxNoFtvFipPHaxst4javkyiojtPhp8yXTwHrZyGzcKB6vJxRSVe9u25zFvj/ZfIvYx5lN58pOT4nky0+hxKo8ZV3zS0vFXfO9qgd6ZZqIiMY4HBbXyGr3cENgNordJIIbGl/R6T094c9OW1NdLfeMbthwdnR4o2Fyus2caS74alprL0J/6MGCVwbhh1ej6IhONgCpXzZfMg2sl7GWtIZCtmPxsvEyGxXG2l5YNd4inbPbjbfXgrLMZmZLIypgWA1G3FhyyOz9hMP6Rt2iZSxqi7ZypfuduZ2BnPW3bmzjlXgdO6Fd7AxEZLRebph++O3oIbMZDQhjv0Wr/LsZFocFrwwiFXG8jCqrrNDil82XrD2GV1NXQRw1xpZB4m8tLIhI4y2z7ppbZetkatiNLdaY2Gx6UcbI2Gww0t3tnh2b3rtPnM7WhBnRMhYVvKqr3X8OuwM54+lCc69G2esblbndei+CU41eUKbw9cou0XTDTPB16jEqCgteGUQqItdbjZBkPYncNkYmSv7YFiyQa7BratzVNokaIadimznTeeRo0U7XaHrAbsMrq/EqKyNavdq592NimBCjb0Q2f1aDEbenawFrxxBRTduGDdYDHK8DxNoJ36KvkYqP42UHEaFb1A4pdhP1vnNiwyaSd6ffkJMt1aGJ9JDpi7PgkFdffRVvvPGG8Pl///vf8eqrrzq9LeMjoRBQUQEMHw7s3Qts3gxEIvrnVlUB69YBpaXi19+71518atTWAuXlwLhxwOWXAz/4AVBdrX9uOKzmt6oqfn9VFbBmjXt56ulx71p6LFigPkssJSViaSdMAD79FNi4EVi5Uv27c6daBm6/mwkT9OtHcbFRCor7pSjq36VLe+tlONy73wxFAaZMAebMARobxfJrdN32duDZZ9X/E+vbuHHq79pa+fIjAnbvVr8xPSZMELuOyLu//npgwwYgN9c4LwCwbJl5GSsKUFYGjB0L3Htv777Ec7TrUfwrdZU9e4Df/z55fyQCbNoErFql/o1tv6qq1Pq/YQMwb566PfpoBN/9rrPKv3mzmh8zIhG1PETqr8bw4e6et29fcnsukvfGRrWeieZdO6+6GvjP/xRLo4dxW5FGOJXyFEWhY445Rvj88vJyCoVCTm/LfIUfGi+7EYfdDDgpiox6vLraOPjixo1qyAU3p8a82rRRfqLWpaPDuY2HqMbmlVfkpnO1EALz5qkxeYymnhKdO/SmuUW0rFpgX1EvJ80410obZHRNWW1R4ma0kLTotLmRXZ62aetcir5fLQCmiL2nXntRUmKsdfZiiw3GO3Nm8jR/Yvult8RUcXGbo6lGUY2TnsZZ5BuyQsb0I7E8ZPIuqsGO/Xad2Gb6aR8sg69TjYqi0PDhw4XPLy8vp6ysLKe3Zb7CDxsvu14xbsW6EUXWnkbv3kEyJtWmZkQ6GT3MwlDIurxbNazhsLEXkWjnHL/FC1xFRcaCstH1tJAZIkKoXgchOlVkVS5OnTWMvEetytmoTGJDlchMR8kE6NQTZOxM3Ysu+5K46a2qYVRWxktM9ZCi9Nju4GWdeuws92WFrDesrGCkGbjrtTPaN6tnW6V9j2bfhUyIniAQaMGrpKSEcnJynN6W+QovBa+ammeotFQ/nERsw+vEDd7NUYudUVSsti1IxqSxDa1Rw2YWFX/tWvPrW0XUj0VEoxSrAbLqnO2Us6igb2SLaCfGnFsek0baIrvPrsWBSxRkRMMYxCJrgC0ay83JdzR4cK/ArGnuvPguFUVEGOyx3cHbHXi6GYFeVvDSvCZl8m4VssJoCTmzgZpofr2wD7ZLYAWvDz/8kBRFobKyMqe3Zb7CS8Hrt7/d7Ljy+yl42fFwiw0OaqXdsBqBid4zeUojWbBKbGhjp+bmzTOPXbRmjZgmRmbKorqaqLDQuiOzWozXSdwtJ6NcGc2OJlz853/ay2fiNnOmc01q7BRi4nWstIFW79ZNrbTTuGqxSzZpBEELbaeDNwuULOLBajfSfOw17JTbLbeYC72xebd7D6uBmmxw5yAg0xf3k7UJe/bZZ/GsZlX6FY2NjbjmmmvM7Mhw+PBhbN68GYqioKKiQva2TApoaBggdF59vf7+SEQ14DVCUYCZM1WD4VBIPn+JiBqT6qURMSY1M5APh4G77wZmz1bLgyj5HEVRz/v4Y+D111XD6+HDVQNxLQ+x+2LLJBQC/v3f1c2M2lrg0kvNzwHUZ928WTWItrre9ddblw2gPvPu3eqzGV1XpJytri+S70RE68ZHH6mG8XbzqMeTTwKLF6v1XHvHH32kGq2L3kd79smTk48dOgQsXAicckqyk4gVoZBqED9pUrwBPJDsyCDCpk32yk6718MPA/37xx+rqoovu/XrgSeekL+HE2SdJKy+m3BYLVej96U5jmjPvHlzcptgRiQC/Nd/2XsXt9+ubuEwcMMNqkNC7HVi8273fROp7/ypp4BPPkluDzdvVu9hhd53HYmYt6WBQFaqW7hwISmKQllZWaQoivR21FFH0fbt2+0IlIwOXmq8rrnmXaFRx5Il+tfwai1AI2TjiMVGvXYSDyp2lG5sL5I8XeR0RGv0/KL5NnIj1/KlTZHJlofZKNSNuFt2Rrkimh07oQ5Ew4QYRYR3ayFspzYvbkxvyYRgEZkmNSIVqxUYtVF637HVVGt1tblW2Oh9OFks28n2q1+pbby2vqvd1Upkvw23pmn9WiXFU43XN7/5Tfz4xz+O/v7v//5v5Obm4lKTYXZWVhYKCgpwyimnYOLEiRgyZIgNEZHxm4KCDqHzjFzXRUeJboUsMBu963HwoBpqIhwGrr3W/n3371dHbAcPqpoNI264QR0l6o2Gw2E177Iai1hktUmJo0UZ7ZbMdUWPyV5fZmQbWzcS0epKV5d1nYlNAwDjxwPPP299vl4dD4V6NXebNgG/+53YvfUgsq8NBJK1SrKagtpatWxFy2/NGvXadu6V+C5F72mXkpJerXQset9LaSlw5IhxnhQFuO++ZG1n7PdvVJb19ep+vfA3sXmSeQ8i3HdffB5j35Mb37PRtyGjiY1E1DAiCxYkX0uk3HzHqZQna1zPuEuQbbz81nhp2B3xOQnuuGKF9T1Fl9+xi8zoU2+xYqcGzCJaFzeWZNLWsJQd2Ro5KeTnExUUyOVBNsCsyELObixVlQqbFze8ie1gVAdEvmPtuLnXsBq9/vrrkzU9bjviaN+/6EoRdsLHuJXHRNtTp/XW6VqYIuvf+uEF6atx/aZNm+j11193ehnGJkH2avQ7nETivbVYXCLLk8RONdlpRERjlpktGyRaHm5ESddrPN1ulI0QW93AuN4ZdaxmeXCro6ysVN+1NrXsZh03CxMhmr9UeHnJrtHq5rSP2VSfWblpnbd5vYivg5pg79VUp6KILytWXZ2cBy+WkxKpz3ZWK5H5Nswci2S/ay+/D14yKINIVRwv0UbU6qP0eu5dNiSAXoNm5iGoNR4rVrjXuMmOABM7BLOGKBRSNUZOykhvs2MPpBdjauJE9zsGLzrKWO2aaFwtu+VSVkY0Z45YPUxFXCNRTatMCBOnWMV1SxQckutHDyUKXtr79HKxcr82p0uYJbZRRvVWJrafzLsMh621g3qblxphFrwyiFRFrpfpaO3EoXILWePP2HAC2ihaxK3arZhPZo2DVbyc2BG8kfAVGzzTbhlpW0mJvsGtKHraCrcWvY7tGLzsKLXydDP2klH9M8uH25okGUTr/oYN/uZLxoElcbUKI42rori7hmsqNk1j62R1DjPHnBUr4o3x9UKgiH4bZm2enWevrpavR6L4Lnj19PTQ448/TpWVlTRs2DDq378/ZWVlGW68ZJB7+CF4ETnzwhMRGLxCViAy0jZZdawiU06i0whmXj5GaWI1HrJCgJ2FnL16b24JsFrHIBtAUnaL1SCmylNVT4vpJ6k0KXAbNwdQdjevlikz8wSUnSqUaSc17ZTst+GVptqreuir4HXkyBH693//d6nwErxkkHv4JXjJkjyCTE2jLPrxihqGW7mAm2nGNNW4nQ7KiwjjsmWkbaKj1dg8JC4ZYm0XaG5bKFIOfoYd8GrwkCoHFVncnG51k0QtzPLlxqERiMQ1rvn55m2J1ZJFZt+/mYbdbv20eg+iDklmbZTbA2yvhGCvvhVfBa/q6uqoQDVx4kRavnw5bdiwgTZt2mS6Me4QRMHLjlehlx2HiAGm084hdq21REEzVlCx20HJRF6XxWgqOLax19ZcEx2t6q3VF7uZeSCqZaRvYwOYe63Fdgx+ai+8Gjx4+d7dxs3pVq/yY1UP3aozms2lTHy22O/fqCxFp80TDe1F3kNsGxb7vYm0UTIaeVHcMjvw61vxVfAaNWoUZWVl0cKFC51eirFB0AQvu95jXnccZo2w085B79pDh+ob8hqdb5UHGTsas+V6ZIM92rHDMwsiK9KIExGtXt1FxcVtumUk6rDhpOG2U4e9GDyki8ZLw4vpVjvItEOx9VBE42q1HJe23iGR+PsrKUn+FvTKUjR8Q2mp/GApFr2BkxvmCjL1lDVeJgwYMIBCoRA1Nzc7vRRjgyAJXk6mdtyMXm8leCQaf1pFkTbDrnpdtoMSjbxu5PGoJ+yVllrFMZK3iVizRvydm42CtXAmdXVdumUk4rAh2nBfdpm+ICzrNeXF4KEv2U/5hWw7lFiGVhpXmfZMRFDSW5/SDBGbLCdTvEYDST3HHA0vNLNuxbbz61vxVfAqKSmhoqIip5dhbBIkwcvOCMXNj0EmqKbZ1KDoEhNeqNetnk/W9sONRktUKO7uthdLSO/6VvVPROAVabg1wdJIEJYxzvdqJC0yPR0UTVMQsKspiX1/RhpX0cWbtUGd2dJbToUjL4KG2h1IeqWZtRsjzM2yFsVXwWv8+PEUCoWooaHB6aUYGwRJ8JKd2nHyMYi63Ovdw8r2QzRfoo3NkiXudYhGwqWV5srJJjJK7e4WDyIrcn2z+ifr4SliU2cmuKxZk/oYWmbT03YGHLL1MZ0EO7tTzLH10EjjKjN1GPtbL26W03A6GzaI5UVm4GR3IOmlZtaOzbCTtUDt4qvgtWHDBsrKyqJ58+Y5vRRjgyAJXrIjTbsfg96HKNoxitp+iDQUdhp4K5W9CImdoGgDbHezarjtNIyJ1098pvZ24/onO7q2sqkzmoaNtZFZvdq4nvjltSdjoyc64BDR7tpNlyrc0HgZtX9uT38Z2YGK4Pb0nlOtlRuDHCO0NKLOCma2rl7hexyvu+66i/r160e33XYbtba2unFJRpAgCV6yjZIdAcSu8b72MbrpbenE+POyy9xrDLzy/hERPp28j1j3+WTBp4duuulN3fpnp8Mxm0oUyX84rDoNBMlrz47mT+8cqzADbk+TeY0dW9PEYMBm7Z9b01+J9curcAuiGi83BDk7gxyZZw+yzaOvgte4ceNo3LhxVFhYSFlZWZSbm0unn356dL/e9v/9f/+f09syXxEkwYvIXY1SIk7jMsm4dos0Mk5Hv2ZTDTKjQi+8f0Q6VifvQ7u+tqRI8nHVuHn16i7bz2vV4cjkX8uvnUCQXiHj6WpnCslvG0Y3cTIgCIdVGy8rG0O9JYm8/N4ScSqE2NWci3xXMoMcWY1YUGPG+Sp4iQZN5QCq3hA0wYtIrfiiRtaJH7HZR+dUwLAjeIlMszkZ/eo1FLKjQqcCkJ43pIgmx8n7EPMa7KFwuMdQIHA66pXNf9AEDVENhWi9T6zrMjaMQSmTWOxOgavfs7HGVSOxrXK6Xqud+mUWusVKk2nk5eyFNklUiNfTfuu1fUaatVQOjGT64n5wyIIFC5xeguljVFUB7e3AtGnW5+7d2/t/bS1w/fXAnj29+8Jh4N571WvGniuDoqjXGTsW+N3v5NJUVJifV1UFrFuXnG9RiICZM4EJE4BQSC2DSZPU/bHU16v7161T7xlLKATccw9w6aVy91YU9e+f/qTef/NmtYyHD1efOxQyTy/zPoYOBX7xC2DEiN7rb95sVWYK9uxRzxs7tndvKKTWiUmT1GeILSvtmZYudTf/gHqf3buT85Mqhg9393qJ5SFaPrNmAXff3fud6hGJ9Navo45S9335pXhds0NVVW+9rq8H9u8HiouBhgb175w56r5EiNR69Oijp2DhQiA7W//6oVB8Pdi0yVl+ZetXbS2weLHx8Rtu0H8fRm3M55/37nPyXelh9a1rzz55cvIxvbYv9t1qbdaBA2pdNOo/AoUPgiDjIUHUeBGJR1iONYC2UkPbDVchE14gMY0oTrz6tHLo7jb3TjQbcdopG6f2SU41Ik5tSpxGSrersZONReTVCFxU82d3CkmmfMy+GSvNk5uG+qLlLfpsdXXJU91m93bDu1jUi9hKm2d3+tiuBtwMp3aoVtq2INgi+m5cz6SOIApeIrGPYj8kUTV0R4f81EFigyEyNWi3kXHSuGjxfkTO1Zv+lJl2cksIcDrl5+R5Y/NgV7Cxa6MnYqzc3a0+X6Jbu1tCRmwcOq2sjTocu+9Jtnz0riNiayXTOZq9b5kpetHvZfly/wUvkfpl185Rxi7QyYDBK8/refP0A2MHwRaRBa8MImiCl6i9UWxDK9OIiCxHo21GmhYjw1gnrt0yz2HU0CV20kab3ojYqwCGVtg1dO3utg4AaWTj5Uf+nTTgRlH1RcpFNM9W4VREBxyiXo0ywmli1HaRNKIetGYrM8hoPLzQeLnh5CK6UoRdbbEfa3/asR+zU052+g8vYcErgwia4CX6EVRX96YRbQxWrJDTeK1YYZxPL6aA7GpQSkrkRoR6DUgq3aztTPmJ1pP58+Mz7MV7EzHCFhWYZLW9dvJqVr/MBg92p2ZljdS1TtuOIGIVI0qvLAH5KXrr76WHhg5tpfZ27wJI622x7WJsXhPrvB8aLzPseC7q/a/3W2SL/R6DspC8r4JXVlaW9BYKhZzelvmKoAledj4CGVshmY9Tb+FZr7GjIVi7VrzcYhfgFb23H3YOsgKRnakeLwN5xua/utqegCLrXSo7AndjSsWu4Cpjw6g9lx1BRK9zdBpGxqi8zb8Xa6/GRNzQeCU+v1Gd1zyCvZo+NvuurPJkVj+N7Mdk10WNfUa3I/jbhcNJZBBBE7zsjMREtTWy7tp+GlbGIqMhmDtXrtz0RsRW905lkE8jZKd63IoBJIqd68h2vLIjcD+mVMyeW1ar6pbGy604dXrlbfS9WMXxMio7p1HtY5/fqs5rMfC8mD42uoaIRstqM7IfW7PGXpndfXcwgqr6Knht2rTJdHv22WfplltuoZKSEho6dCitXr2aNm3a5PS2zFcETfCyO+Uloq2x69Xo5ken1zFZ7dPToJSUxEfuF2m0zbRdVnkMGtbP20ODBh2hV17psnSq0N6xaAwgr5DV8MgKSF5PqYhoFGW0qjKCiNl36tbKDEblrfe9OA0gLSt8JT6/k7hXotPHVjalsnkS3YzqpxMBW5tqToW2XyOQNl779u2j0aNH08iRI6mlpcWv2/Z5giZ4ETkz5DVrRJx4DbmhZtbLX3Fxcp70OnsRYciq0Q6a1sopop2UaDBeo87Dr4ZXpuOwMxjwUuMlY5wuo1V1w6vRqcbLzuDLafsnI6DoPb/Mu7Yz0BKxRUysT25pHo3qpxMB28jWz09tfyAFLyKijRs3kqIodOutt/p52z5NEAUvIvtTXlbu4nY/TKeGlTJLkDjp7NNlqtAt7EYX97rjtYNMnDg779MrBwo7tmOinb2Zh6do/RZ5bs1jzi2Nh9P2r6NDfAkhvef3Urspk7fYe3gdi8sNATscVqcygx65Hj7kJ0pPTw/l5ubSiBEj/LxtnyaogheRu1NeTtXcTu1e7Bp+2nnmdJgqdJPubrlwGqmoA6JYafHM1udMRK8eeOFA4ZUmzYkHptG1zJ7bzUGL0/ZPxmFI7/m9fCeyGmQ7Gi879dMNGzm/vnM9ZPriLD+j5CuKgqysLHz22Wd+3pZJEdqSGlOnqn+dLAtivbyMPooClJVZL/3j9r2Jepf/kMXNcksHQiF1O3jQ2/vYXXJKBm0JqdLS+P1FRUB1tbosS1ERsGqVusRMJKJ/ndpaoLwcGDcOuPxy9W95uXpM7/rhsP5yUiKIlotM+UUi6jJaRPrHFQWoqRFfLsioXGOfu6oK+PRTYONGYOVK9e/OnalZLka0rI4+Wv/5KyrUZ9OW6tEjFFKXyRFFWypIJk1RkfouIxHrPGlt7dq19uqnthSYdi27+PGdO8VXwevtt99GW1sb8vLy/Lwt0wew8zE5XV/Myb1l0kYiaids1Rn3ZfxoLN1e29AIPQHgyy+BU04Bvva1ZGGqtjY+vdZBJgr72pp1gLsChmi5yJSf6Np8MgMTEcEqFFIFhOHD1Tq1ebP+9+T1N+e0TGOFECMiEXVtw7Vrre9jJQgbcfAg8IMfqPX02WeNBaPYtnbSJPv100zAvuoqsTz79Z07wgcNHBERvfnmmzRq1CjKysqi//iP//Drtn2eIE81uomd+X+ZaQaz6T0ntgdWam8vY1NZEaQpTSdlXFQUDHdyM2TCYfi9/IkXtmNObZTs1k1Rz0yrc5y2f26V6Zo1yasSJG6hkOrdaIYb9lNuTulavV8jT/Egf+e+2niNGzfOdDvzzDMpHA5TVlYWKYpCOTk59Pe//93pbZmvyBTBS+Sjs2tYadUQ27E9EGkErGxg3BAa7T6z3zi176iudt/2yY1n2rhRjT9nZswcW1fSbeknI5w8h1mATrN6LiLcigrAbrR/bpSpjMBkdj03QnLE1lMngzYnbU8qg0RbEdgAquXl5fTKK684vSUTQ6YIXkTefHSiDbFMfB4ZQ1Kz64jE7bLTiMmEDvATuzGQALUDCJJHqB1vTa0jE31eP/Jst/ycxvMTKQO9wZHZNxkOi2sT3Wr/nJapjMDkpcdgYj1NRMbL1Unb092tv/h8EDy/ZfpihYjIyVRldXW16fF+/fqhsLAQ3/jGN3DWWWdBcWI1xyTR1NSEwYMHo7GxEQUFBa5dt6urCy+++CLGjx+P7Oxs167rlNpa1VYh1n6krEy1LZC1cYlEVNsFI1sURVFtC3buVG0u9O5dXKz+bWiQy8+mTaqNjxXV1cD8+frHNFugxC9Y+8T0jFlln9lv9MpYhI0bVUeESES169m7V7X1EDXedhOj92LFypVqnkXqhfa8buNm+WnlAOiXRU1NfP20qpuJxNbzoiKxchNh40bg7LN727+srGxHZeKkTEXbidi869ULrWzr6+XrZSLz5gELF/Y+g943Gw6r9mAy79eq7dG7T1GRuu/WW1PvhCTVF3suBgaM9evX0+TJk+nYY4+lnJwcGjBgAB1//PF0+eWXW0bU/+KLL2j27Nl00kkn0YABA6iwsJDGjBlDy5Yto56eHst7f/zxxzR9+nQqLy+nnJwcKikpocrKSlq3bp3t5wmSxssvmyG37mN3eSORyPVWOF2bUURjpjcCTtVUlgzd3epSQbNnv0Uvv9wVaLuORJyEPdHqTjo9rxVGcbz0wmrYXZkiHCa65Rb3NDorV/a2f6tXd6V0Sl42aLSZJtSJRjlx08pARoPldPo5iFr6WAIbxyuV9PT00HXXXUcAotuAAQMoNzc3bt+sWbN002/dupWKi4uj5+Xn51O/fv2ivysrK+nIkSOG93/hhRcoLy8ven5BQQFlZWVFf1999dVCwlsiQRG8gmYzJEIqp3VkOhkn69clpk3VM8sKp7H1L8h2HYnYndIpKVFtwTZuVA2lzTrI2KWmgo5Mh+nW0kBOt40b1fp3001vkqL0GJ7nx3uQDRptx5nHzqa9UzOhMHGQYLftEXU46egwX2fUa6UAC146PPbYY6QJOZMmTaJ//etf0WMffvghTZgwIXq8trY2Lu3hw4dp2LBhBIBGjhxJb731FhERdXR00P3330/Z2dkEgGbMmKF77x07dtDAgQMJAJ199tm0fft2IiJqbm6m+fPnR++7aNEi6ecKguCVDqMRPVKp/enuFg8YqicEiTZiM2f6+8xGQT9lhfLE+hck+y0z3BAewmF1AWSjziboAxoNWQ9NN+2Q9O6l2XiJaBPb2zupuLiNAGPBS8Sb0MvyMytLq+tu3Eg0b5535a3Xlthte0TTJTqwxGrl/FAKpETw6unpoZqaGpo0aRKVl5dTXl4eDRw4kMrLy2ny5Mn09NNP29LouMXYsWMJAH3961+nrq6upOOdnZ10wgknEACaMmVK3LF58+YRAMrNzaUdO3Ykpb399tsJAIVCoahQFcu0adMIAA0bNowOHTqUdHz69OmkacEOHjwo9VypFrxS4f7uFqLTOmYjKSdUV9triIjkGqPEkZ9XU1lGa1ka3cdMKNerf0EKf2GEG8KDVjZz5tgru6Ag29G6FbncrLxEtad1dV3C1/fqPchGipfNh1flnbhpA0e7bY/dwYzZc3nxDfkueH3xxRc0ZswYysrKioaNiN20/RUVFbR37143binNiBEjCABNnDjR8JyqqioCQD/84Q/j9h977LGkTQfq0dzcTPn5+QSA5s+fH3espaUlOp1ZXV2tm37nzp2kab0ee+wxqedKteCVDjZDZlg1xHqaB7dGS1b2G2ZCUHe3+NIfiWXvp3eoVcNo9HxB9KoVQaQzKykhWr7cer08s/hNemXnlWBq97p2ppbctEPStkTNqIj2dPlyccHLq4GljB2o3fbIi/I2a3/stD1eaULdVgr4Knh1dHTQN77xjajA9b3vfY/mzZtHDz/8MD388MM0b948OuOMM6IC2De/+U3q6OhweltpLrjgAhLVeC1cuDC6/8MPPyRNKFpjMqmvXf+MM86I2//yyy9H05vFLxs1ahTpadusSLXglUo7KStkXJz1GuK5c72fQnUyTfurX4mV/YoV4s/sJHSAGw2zRroKXkRinYtbnYlWdl5Npzi5rt1BmVt2SIA6nWY0cDFrG2Q0XkZ12Cmi5bdhg7P7GLUFohp5o80sZIhM2+O1Zs6td+er4LV06VJSFIUGDx5Mf/nLXwzPe+GFF2jw4MGUlZVF9957r9PbSvPcc89RrI3XRx99FD324Ycf0sUXX0wA6Gtf+1pcwa1bty6a7v333ze8/ty5c0mbLoxl8eLF0fStra2G6SdPnkwA6NRTT5V6Li8Fr5qaZ6iurstUcPFL4yU76pbtMBKv39Hh3xSqXSFoyRKxsl+yRP95N2xQO6Z589T/vfYONdr0hPJ0FryIrN+pW4bkWvwyLwYIMtH29b5NJ9Paiddcu9aeMGa33RGx8YrdEm0p3cBPD1e7keKLi+1pz+20515p5txSCvgqeJ111lmUlZVFK/SG1QmsWLGCFEWhs846y+ltbbFkyRLq378/aYJQbm5udBpwyJAhNGPGDGpoaIhLc99990XPNyvQpUuXRs9rbm6O7p89ezYBoMLCQtO8zZw5kwBQcXGx6XlHjhyhxsbG6LZ7924CQAcOHKDOzk7XtpUrj3zV8PRW0NLSHlq9uivuvPb2Tiot7TH0/lGUHgqHe6i9Xe7+7e2dVFfXRcuXd9H8+d1UWtpjmRdtW72666v89CTlRVGM08VuoiPeujrra8k+b11dl1B5/ed/dgvl8Ykn4vO4enWXVHlabTLTMqJl2NraSs888wy1tra6Wq/d3szem9kxWY2K0fbyy9q7dP/7E7nuqlXmdWnVqq6vrmH/W9QrT6+eO7b+3XTTm4bXT9xKSuzfy2zT2rLE9tVO+Xl1f732JBx2P2969ykoEHs/frThBw4cIFHBq595lC9rPvjgA2RnZ+Oyyy6zPPeyyy7DT37yE3zwwQdOb2uLmTNn4sQTT8Q111yDL7/8Eu3t7dFjHR0daG5uRmNjI4qKiqL7m5ubo/+bLe4de6y5uRn5+flx6a0WBteOx95PjzvuuEM3aO369etdW3z8b38bjkWLvpO0v74euOyyEG666S2ceWbvqsbTpmnnE4DYALkEIuCKK97CK6+Ir4L8t78NxyOPnIqGhty4a4nkJRIBfv7zShCFEvICECkACL/4RSf69aszDbj36qulAE63zOtLL/0/tLbWWz+UIAUFQGsr8Mor6rO8/34xDh0agMLCIxg9uiGa57/9bTjuvz/5Hemxe/cbePHFhmg6mXebiF6edu0qBjBG+lkBwtCh7WhqqsOLL+qfUVdXZ+O6/qBXT4uL2/HTn/5fXBnGvlONSAQoLq5EQ8MAJNZTFUJWFqGnRzE8PnRoO9588x+orzcueyIFe/YAixe/iVNPbTA8L5H/+79ioetOnZr8EWl16eKLP8bmzWEAyd1McXE7fvKT95CTs9fw3RtRUABs3DgcjY3fAJCjlztb7U4iZ54JzJ37FhYvPh1EWabn7t+vSJexCDk5wI036tczu+Xn9v0B4L779NsqN/OWkxN/n8GDj+Dee08DYPQNAb39hvE3ZNb+yNDW1iZ+sqVoZoEWSFSUwsJCGjBggNPbStPa2kqXXnopAaDTTz+d1q9fTwcOHKD9+/fT+vXr6fTTTycANHToUPrnP/8ZTff73/+eoL49XdswjT/96U/R8z7//PPo/muvvZYAUGlpqWn+brnlFgJA/fv3Nz3Pa42X6Eg3cXTn1qjHSFslmhe3NFV+a7xERneaJsHqHfVu8eVj991a5UnTehjHPJLXeARd4+WGVtVKmzB7dreltkFU27h8uVw9FddiWr3zxOPqvlWrkvMjqvW1aiOKi51rW2Lrn6hmWbaMZTY7GvEg39+N64lrjf3RGMpovBwLXscddxxlZWXRrl27LM/duXMnKYpCxx13nNPbSvPzn/+cANBJJ51EbW1tScfb2tropJNOIgA0ZsyY6P6gTTUm4raNlxObLadeVXaNtGPz4paxfyojiFvZ1sgYvfoVOVpzRNCz9VA7w/j9VjZsnZ3BtfHq6BBf9NoKK1swq+Ne2Vh6HVMrsXxEbTJF2ohw2Pl3GVv/0t1zW7Zd9jpsi1uOIDKxDP2IAeirjdePfvQjUhSFqqqqTON09fT00CWXXEJZWVl05ZVXOr2tFE1NTdEo8/fdd5/hebFC1r59+4go84zr0yWau1Fe3GwkUxExXSQummjg1USDX68jR+sZQGsNnJPI9V5gt3OpqbEfxsNuXsyOezVA8CPGU6xHpqhzgF9CUGz9S+dlnGSFHK+DjbrpCCJTF/pc5Pp33nknGqfrnHPOoQ0bNsQ1lqpKsI7OOeccUhSFQqEQvfPOO05vK8XWrVtJE35eeOEFw/Neeuml6HlvvPEGEWVeOIlUju5WrHDWgBO53xH5HTHdTU1DYrmIekE6SedWA+el4GW3c5GNVebG4ESkPL32ajTSYjrdVq6UD8Ds18Awsf6l07JVGrL1wusVSES0lUVF4h7WQROIfQ+gumTJEooNlNq/f3865phjqLS0lPr37x8XVHWJnm+7x/zjH/8gTfh58MEHDc974oknoudt27Ytul8LoHrNNdfopmtpaREKoHrbbbfppv/000+j9011ANXeymxkB+RNZZbRJMRuiZHZtWs5bSRjO7wNG9TNj4jpoh1LUZF4gyMaF8luOjcFDQ2vBC+7nYudaXCngxNRAVFmIWo38lBWpq5T6FQjtnGj/EAvFRovIvX9V1cna5udDMK81MLICrR+rEAiM6gU1bIFSSBOyZJBf/nLX2j06NGUGLVe204++WR6/vnn3bqdFG1tbVHh59vf/raukXx3dzedddZZpNljdcfUMG3JoLy8PNq5c2dS2kWLFhFgvWTQ8OHD6fDhw0nHZ8yYQQBo0KBBgVgyaO5cIjPD7epqdxsLO1HPtc0ofo4TTVUqF/wWbZyqq8UaHNGytZvOzc4uFi8ELye2WbLLtzjtpEQFRKv3ZBVrSwSjtHZjK8WWj6wGyy8tR2z902sPiorUb9DufbxuY4Io0MrErpMRnPyelTAipYtkv/vuu/TYY4/RnXfeSXfeeSc99thj9O6777p9G2l++ctfkqZVOv/88+ndd9+lSCRCkUiE/vnPf1JlZWX0eOLSPrGLZI8ePZq2bt1KRGrU/gcffDAaG0xkkeyKioroAt0tLS1UXV1NiqIQEIxFsnsbcn3BKz/f3cbCatkcJw2Bnc7Ga3W7FTIdi1WDI6OlsZvOzc4uFrcFL6e2WbKdhtNvQkT7IBrkV8/2zq1O3mrVB6uBgZ0O3w8th1b/VA9K++2BXhtkp42RbctkBVo/pnBlzShk2hU/bLisSKngFVTa2tro/PPPJ024AkA5OTmUk5MTt2/q1Klx2i6NrVu3UnFxcfS8QYMGUXZ2dvR3ZWUlHTlyxPD+L7zwAuXl5UXPHzx4MIVCoejvq666ytYi4m4KXnamU5w2dk6WpXC7s/dD3S6CTMdi1uCINnRLlsSnk3knXgmkbgpebthmiZZlSYnzspB5b3a/HTffm5lGzEoTYVeD5ZaWwyjvqqbrmaTwKTLtgV4eS0vl12e1ox0LosbLrsOG03v6JZCx4GVAT08PrV27liZMmEDhcJj69+9POTk5VFZWRhMnTrScCv3iiy9o1qxZdOKJJ9KAAQNoyJAhNGbMGFq2bBlFIhHL+3/88cd07bXXUnl5OfXv35+Ki4vp3HPPpXXr1tl+JjcFL7uG3U48p0Q99Pzo7IPkNu5Gx2J3kWKZd+GVSt8twcst2yyRTqOkRNVCOUX0vf3nf9r7Xp1+tzLIOAfIarCcdqpmAk1nZyf99rebbdcXJ+YTsdd0apMoKtDaOd+uV7Ds9LRdLZvfJiOeC1633347fetb36Jp06YJnd/T00PTpk2jb33rW3T33XfbuSVjgJuCl9P14/yKFeRVZx+0Bb+ddiyygqSskJKoKXPzedwSvNy0zfLLkNcPjZeT79YL/LbTsRJoVq/uotmz37LVHtgR9vWu6VQDL1tfRc93KtDIOu3YqZ+pMBnxVPD64osvKDc3l7Kzs+M8/6zYtm0bZWdnU35+ftJ6iIx9gqDxMmqArJDx4PPDqzBIGi83kB3Fyrx/WU2JbGPtluDlxDbLyD7HawFB9L1pNl5ONCt2vluv8GtaSESgCYd7qLransbLjZAwdjw+9ZCtr1bnuyXQdHerbbrZjIeTmZRUmIx4Knjdc889pCiKrSCoV199NWVlZdH9998vnZbRxwsbL7sNuVcarwRfB88IWlwYN5AZ9coIKTKChp3G2m+NV6Jtlpmg6IeAIKN9MDrPq+/WiCAYOIsg3u5s/mopLLn2wMnMgROPTyNk34vR+V4INF5okVM1gPZU8Bo/fjxlZWVRXV2ddMb++te/kqIodOGFF0qnZfTxyqvReN099wQSEUGvuNjfBjxIcWHcQnTU64UgbLexdtvGS8Y2y4tpCrsetiLvzW6sLbPvVja/qQzBIouoQDN79ltRr0aZ9sCJraxTj08v8So/bmuRU2Uy4qngFQ6HKSsry9SDz4gjR46QoihUVlYmnZbRx4s4XjU1lOTNo3niuC2QWBlbpqLhDkpcGDPcGsUmnmMlpMiug2e3sfbCq1HUS9SLUb1doUT0PcvG2jL7bmXzm+oQLLKI1snf/nazYRwvPe9Mrfw3bLAWeIuLvfP49AovBRo3taV9UuOVm5trueCzGYWFhZSbm2s7PROPF4IXEVF7u+rVs3x5l+f2LUEUdII8beKldsFtjZ/dxtqLOF5uav1EG+0gCCUy35dsfoMSgkUGEYEmHO6hmpr4yPVG7YFe+YoMVL30+PSCoGngjEiVwOqp4DVw4EAaOHCgrYy5kZ6JxyvBy6jj80ogCbKgEyT86MjdFISDoPHSsKpj3d1E8+aJ5XfePDEtY1CEEhmNp0x+06UzTsRKoFm9ukuo/pl9j7ECmJPvKCgD06Bp4MxIhcDqqeB13HHHUVZWFu3fv186Y/v37ydFUei4446TTsvo47fgxfiP1mmuWGF/yRu793QqCNttrJ3UP7fsqaw2t4NYpho74UdEBVW37WncqJ9mAo1I/RPzjnTHIzsoA9MgaeCs8Ftg9VTwuuiiiygrK4tWrFghnbE///nPpCgK/fCHP5ROy+jDglffxo5AEJSOXMNOY223/tmZhrUb7NKqswlaXDgrZPLrRywmI9ycajcSaETqX7oJ1m4RFA2cCH4KrDJ9cT9IUllZieeffx533nknLr30UmRnZwul6+rqwp133glFUVBZWSl7W4bpc0QiwObNwN69wPDhQEUFEAr1Hq+tBSZNUps2GfbudTefTqmqAtatA66/Htizp3d/OAwsXaoedwOj8qqvV/evW5d8r0hEzZdsGQNqGkUBZs4EJkyIf3eA+k5FED3Pa0Tz8dFHwMKFYmWmKOp7rqhwlLUodt6xGaEQMHasvbyIfmdW58W2A0cdpe778kv9NiEIVFWp9d2s7QoKTt6vp8hKdU1NTVRUVERZWVk0ZcoU6hBYI6Ojo4OmTJlCiqJQUVERNTU1yd6WMYA1XumJ1ajdSfTroI6wvYxcb9eeyo1gl0Zlnk42MUSiRufi9dLt6Sc/bebc1HiZhV+x0hwGNSQHk4xMX5wlK6gNGjQIf/jDH0BEWLNmDU477TSsWrUKLS0tSee2tLRg5cqVOO2007BmzRooioJFixZh0KBBLoiMDBMsIhFg0yZg1Sr1bySif542ao/V/gC9o/baWnU0mXjcCkUBysrc0y64jTb6nDpVzePmzdZlJYpVeREBu3er58XilnZQ7zqhEHDvver/ihJ/TPu9dGlwNAUi+b32WvF6GQ7La6DMsPuOvaKiQn3GxLJKZOFC9ZtOxKgdiCW2TUglom0bI4hd6e6WW24hRVEoKyuLsrKyqF+/fjRixAg688wz6ayzzqIRI0ZQv379KCsrixRFIUVR6Oabb7Z7O8YA1ngFA1G7E9FR+4oVchqXIBq3GiFSVrL1z649lZcaL7PnDapNDJF5fkXLWcTrUxY/beZE619Njdi3maiJk9Fop1ozmk7BcVOJ54tkazz11FN07LHHRgWrWEEsdl9ZWRmtWrXKya0YA1jwSj0yIR68WgA5yB15LKJlpRdHzgy7hs4iQWNDIefThUHxShPFKL+pNCj3895m7V9i2SxcKJ8vOwJ/KkwIghCHLl3wTfAiUivo2rVr6Wc/+xlVVFTQ6NGjadSoUVRRUUE/+9nPaO3atdx5ewgLXqlF1u5EdNS+YoXYUjcrVqRHR04kXlZr1yavnGA1whZdGkivvKy8LufOTR8Xeq9Jpd2an/c2av/0tD/5+WLfdKwmzs56jn57vwYpDl064KvgxaQWFrxSi+woXOZ8K4FgzRrnWhQ/NTFyo/we3WcWCQshEhoiUZCzmg5Mt+lCL0llLCe/7q3X/tkNO5LYBhClh8YrU8Nl2IUFrwyCBa/UImt3IjtqN+rw5851bnfht+2GnVG+WdmIPpPRtRI7apHI9uk0XeglqRRE/bh3YvvnxMsYICoo0LfxEhHkUqVZSrc4dKnG0zher776qitG/d///vdduQ7DpBLZWE2a59ikSao3FFHvOXqebnoxcw4cAC69ND4tIBfLyO14SCI4jVdF1Ou1ZhSbJ7a86uuBWbOA/fv1r5UYg8sq5k9gYwJ5gFWMuVTGckrFve14GcfS1AQ8+2zvN2XWDsSSSu/XdItDl1bISnWxBvR2t1AoZEuiZJLxQuPV3U1UV9dFs2e/RXV1XRk9srfCrt2J3VG7G3YXqbLdcKo1kB1h81SJPdiLLVnj5VRbq33fIu1AKrSIerhlU5cpmmJP43h9Jaw53phgUlsLlJcD557bD/fcczrOPbcfystTH0cmqNiN1VRVBXz6KbBxI7Bypfp3505rLZMbsYxSFQ8ptqycIDrCdiuyeCYhEmMuE3FDq6P3TSW2Axs2qJtMm+AVbsSh0/qTceOAyy9X/3J/AkhPNe7cuVP6Jl9++SV+97vf4fnnn2ehK8CkYvopaFhNsehhd0kcO1NXbggTqRRIqqrU6b2lS+XTyi4/w1Mlcpgtn6Q3NZtJaMFS6+uNpwVFMAq0G9QpbCfLfXF/YoLr+rYYWltb6bbbbqPBgwdHY3uNGjWKnn76aS9vm1G4NdXIrsPOp1j8UKm7MX2W6ik4Ox5ddrzW0m3JnlSTinoR1GkoL7wa03laW/Y9ZWJ/knKvxu7ubrr//vtp2LBhUYGrrKyMHn30UYpEIl7cMmNxS/BKdWecatIlUKAbwkSqBRKxoKXx4STs2rqkMvRBuuG3F1uQbcmMvLrnzrUncPVFQcOMTOxPPLfxMuOpp57CyJEj8atf/Qr79u3DkCFDsGjRInz00Ue45pprkJXl+i0ZF8hkexirKRZAnWIJwvpkbthdpHoNQav7KwqwYkUEv/3tFixf3u3I1kWbKiktjd/v9jqCfQHRKdejjnK+bl862pJFIuozW5EO63J6TSb3J0K4Je298sor9O1vfzuq4crLy6ObbrqJDh8+7NYtGB1Y4+WcdHx2N2IZpTooqNn93Y4jF9QprSAhogktLnaupUqHaSi9+ifaTpSUpO6bCgrp2KY6xdM4Xols3boVN998MzZu3AgiQigUwtVXX42FCxfimGOOcS4ZMr5gZTwqa9icTtgdndkxxHcLN2IZpTIWk9X9u7rcvVeQDZiDglWMOSKgoSE5nayxtIxXbZDemWg7sWSJqmFNxTcVFDK5PxHBtuD18ccf45ZbbkFNTQ3oq5K95JJLcPvtt2PEiBGuZZDxB6uAfkR9V1Vux/uttlbf0+fee/2bvnJDmEi1QJLq+zPxGHmxDRmiDjSampLTEMl5PKbrNJRoO1FaynVaNlB0piFtcPXFF19gxowZOPnkk7Fu3ToQEc455xy88cYbqKmpYaErjdEa3aKi5GPFxf7nxy+00VmibYaGogBlZb2jMz/tUyIR5/Y0DCODFluqurq3LTh0SF/o0ojVUlnx0Udi+ZAN8+H1tyLbTmQ6bF9pguw85sCBA6N2XN/85jfppZdekp8MZVzD7cj1vd598osUpzOi3m9+2qcE2evLa3it0NRiN3SClcdjdzdRaan1dcJhuW/I7W+lvb2TfvvbzbR8eVecTaDVQuwzZ7INYSJ27CvT0SbT03ASsUsGHXfccXT88cdLbyeccIKtB2OScVPwSgejVy8RMTb3y2g0XcJbeEVfFryC3qk4WdrJqt6Lfj/V1eL5tfutGL2Hmhqi0tL4gWesEKfXToRCmTlA8oJ0HXB6Lng53bKysmw9GJOMm4JXJnqiJGLVKfoR6yjTBWCivit4pUOnYjfArUiddPv7sfutGL2HuXPFNP5aOzFzpvF9M2GA5DbpPOD01KtxwYIF7s51MoEhXY1e3cTK2NvLZWg0L8m//jU9vb4Yc9JlCRXZ71vGWNrt78eOh6TRe9izB7jrLu1XvCEXUbIDQUUF8KMfGd83k5dYsoOdJatS6VnuBBa8mCi8tp01XrlJ63lJWtGXBeC+Rjqtgyj7fYus26fh9vcjO1g0ew9WJApx6RoWI6jIlmcQPMvtwmHkmSjstWONF1HfjbwkrchkATjdkOlUUo1VOwAAJSXAihWQXlXA7e9HdrBo9R5E0IQ4niEwxo6HqUx5puPKB7Gw4MVEiW8U44eEQYi9EpTQCm66SdsZgbMAnH6kUyctsqTTww8DV1yhah5k2wM3vx/ZwaIb5asJcTxDoE9tLVBeDowbB1x+ufq3vNxaGJJZsipdlngzxAebM8ZD3A4nQaTv1ZPqZS+sjJJT4Snmxj1lDZnTwcjUDfqacX06Oq54vaSUW9+szELodhwHEp9fy2eqF5sPIk6M40XLc8OGYH5Lnno1MsHCC8GLyDiOTSqw+pjnzg2+p5gRol5eQRGARXHaqfY1wStdO+mgh77QEBUSrd6D1bZ2bfJ9RYW+vo4b3tgi5emHZ7kdWPDKILwSvILS8dmNKZQuDZ/oCHzevGB3fLG4ETIhKPXPTbiT9o7ublUTMm+eum3YYPytWAVBldWipHqx+aDgllbXqjyDqj1mwSuD6OuCl5OpgaBqEWJJV02IESJTDSJalKDUP7fhTtp97Aj6emmcaFHSRTPoJW5qoszKM6htpqdxvBjGT5wYwxIF3527Ly0mKxIyYfr09HUBd4OqKjVkRDrGHgoidmOjxb6H+nrgF78AGhut72dkAM6LvbvrbGBWnn2hzWSvRibQuOERFARPMTP6ymKyIiETGhrS1wXcLbROZepUe16BjIqVoA+Ye7dp7+GKK4Bly6zvx57E5vgZjijd20wWvJhAIxJTyIp9+1IfgsKKqirg00/VuEgrV8rHRwoC9fX20ol0kkz64FfYFzdjo02eDMyda3xcUYKvRUk1XsQ4NCOd20yeamQCjZlaWTT9rFm9v4M8rZXO0xW1targZBen08LpunRIX8PPaOJux0b7wx+A734X+PnPgf37e/eXlYlH5s90NE2UXh3wogzTtc1kjRcTeIzUymVl6ihVC+qoR+JoO9OmtfxAs7M5cMD5texMC9sN2Mi4i9/RxL0IYDppkloH6+q6MXv2VtTVdaeNFiUopLMmyi9Y48WkBZox7KZN6gaoI52xY4EzzkgeYYVC+lMcmpF3UNbFS3ecrH2nh6xNX7osPN3XScValF6tmxoKAeecQ2htrcc553yD2wgbpKsmyi9Y48X4jl0bkGefBa66Cvjd79TtBz9QNRtA/AhryRLza8rYfjDmiK59N3QoUFzsruGtU+Nqxj1SsRal3zZFDOMWLHgxvmJ3WshqGuPZZ3s9xY4+WiwvQfd2TAdEy3DpUuBPf1L/d6uTTKeFp/s6qVqLMt2925jMJCMEL0VRhLdx48YZXmffvn2YM2cORowYgdzcXBQVFaGiogKPPPIISGCu5ZNPPsF1112H448/HgMGDMBRRx2F8847DzU1NW4+bmCxawMiq9ngxWv9Q7QMS0vd7yTTaeHpvk4qvzm2KWLSjYyw8TraQgXS1dWFgwcPAgC+853v6J7z9ttv47zzzkNDQwMAID8/H83NzdiyZQu2bNmCtWvX4rnnnkNOTo5u+hdffBGTJ09GW1sbAKCgoAANDQ1Yv3491q9fj6uvvhqPPvooFCdxEwKMExsQGc3G2LHe2X4wvWhehPX1QEmJalgvUtZuBhBlAdsd3PAITfU3xzZFTFrhfSD94LN48WICQADoww8/TDp++PBhGjZsGAGgkSNH0ltvvUVERB0dHXT//fdTdnY2AaAZM2boXn/Hjh00cOBAAkBnn302bd++nYiImpubaf78+dF7L1q0SDrv6bJkkJP1tewsRcHr4nmH6HIrTspapP4FdekQLW9BX0Kmu5uoupqoqEhuuR0j+tI311eXrGK8g9dqlGTUqFEEgMaMGaN7fN68eQSAcnNzaceOHUnHb7/9dgJAoVAoKlTFMm3aNAJAw4YNo0OHDiUdnz59OgGggoICOnjwoFTe00XwcrKOl12hjdfFcx+jtRj1NidlLVr/gtjZu7FIuNfU1BAVF7svMPeVb44FL0YWFrwkeO2116IapyeeeEL3nGOPPZYA0NVXX617vLm5mfLz8wkAzZ8/P+5YS0sL5ebmEgCqrq7WTb9z585oHh577DGp/KeL4OVE4+VEs5EOmod0QXsPZu+vpIRoxQrnZS1T/4LU2YssEp5qamrEtJV2tYV94ZtjwYuRRaYvzgjjejMeffRRAKrN1eTJk5OOb9++HZ999hkA4IILLtC9Rn5+Piq+Ml5Yv3593LEtW7agvb3dNH15eTlGjRqlm76v4GQdLydu47wunnuIhI7Yv181nvezrINiXJ0O4S20PFpBZN8jlL85hjEnowWvlpYWrFmzBgBw+eWXIy8vL+mc9957L/r/KaecYngt7dj7779vmP7kk0+2TL9t2zaBnKcfTmPusNt46gmyF2EQOvt0CG8hGndNgz1CGcZ9MsKr0YinnnoKLS0tAICf/vSnuud8/vnn0f9LE3v9GLRjTU1NaGlpQX5+flz6wsJCXcEuMX3s/foaTtfxctMjjpHHay/CWO+6khIl7QKfBlkwtXtv9ghlGPfJaMHrkUceAQB84xvfwGmnnaZ7TnNzc/R/M8Ep9lhzc3NU8NLSm6WNPR57Pz06OjrQ0dER/d3U1ARADYnR1dVlmlYG7VpuXhMALroIGD8e2LJFiQpPY8YQQiFA9FZnn937f0+PujHec8YZQGlpP3z+OUCUPGesKITSUuCMM7qF36XG008rmD07hPp67br9UFxcif/6rx5MmuRuHXSDSCS5DpeUKBBpUktKutHVpTMf6QOieQQI4bC9d9kX8Kr9Y/ouMnUlYwWvbdu24c033wRgrO0KInfccQeqq6uT9q9fv95SuLNDXV2d69fUKCgAWluBV17x7BaMy0ybNhyLFn0Hqi9IrPBFIAKuuOItvPKKnFrlb3/TrhlPQ8MAXH458M9/voUzzwzOnNff/jYcjzxyKhoacqP7iovbcc01/4fi4lPR0DAA8WWjQRg6tB1NTXV48UXfshtHJAIUF1ea5BFQ3629d9nX8LL9Y/oWWoxOETJW8NK0XQMGDMAVV1xheN6gQYOi/7e1taGgoED3vNhCj02j/W/1UrTjsWn1+PWvf43Zs2dHfzc1NaGsrAyVlZWGebNDV1cX6urqcO655yI7O9u16zLpzfjxwLe/HflKO9W7PxwG7r47gksu+RaAbxmmT9QUnXkm4Re/0JqhREFAgaIQnnzyO1i4sDsQU8pPP63gD38IJRnQHzw4AIsXfwezZvVgyRJAFUR7n0dR1AQPPNAfF1003r8M6/DggwqmTAES86hRXAw8+GDyu9TT8gXhnXgBt3+MLNrskwgZKXh1dnZixYoVAICJEyeisLDQ8Nxjjjkm+n99fb2hcFP/VS9UUFAQnWaMTX/o0CG0tbUZaqW09LH30yMnJ0c3On52drYnDYRX12XSl0svBSZOTLS1UxAKmTcntbXJ9n1Dh6pR740gUrBnD/DGG9kpj0weiQBz5hh5LaoCzNq1IaxZA8yalWjHqHxlx5iaJjcxOr1eHouK1Pdz663J79Lo3U2bptpd9lVbS27/GFFk6klGCl7PPvssDnzV2ltNM8Z6Mr733nvRsA+JaN6Lo0ePNky/bds2wyWJtPRmno8MExRkl2jR1ulMFFrMhK5YguBdJ+IRuHu3KpB8+mlwnED0hKZwGLjnHnW5J6s8mr27pUvVLRxWvZbZu5hhrMnIcBLaNOPXv/51nHPOOabnjhgxAsceeywA4OWXX9Y9p7W1FZu/8hGvrKyMOzZmzBjk5uaapt+1axc++OAD3fQMk+6YxbcSJQjedaLC37PPBiO8BWC+MP1llwEHD5rnUfTdWS10zzBMLxkneH322WfYsGEDAOCaa64RWpT6yiuvBKCGn/j000+Tjj/wwANoaWlBKBRKshcbOHAgJk6cCAB46KGH0NjYmJR+0aJFAFT7rosvvljmcRgm8MjGjopFUcgwsK7fiAp/Tz6Z2iCpGm4EdBV9d0EJEMsw6UDGCV6PPfYYenp60K9fP1x11VVCaW644QYMGzYMbW1tuPDCC/H2228DUG3FHnroIfzmN78BAEyfPh0nnXRSUvrbbrsNAwcOxN69e3HRRRfho48+AqBqym677TY8/PDDAIB58+aZ2psxTDpif5pQ7c3NAuv6SUWFOo1oxf79qQ2SquFGQFeZdxeEALEMkw5klI1XT08PnnjiCQDA+PHjMVxwCDt48GA8//zzOO+88/D+++/j9NNPx6BBg3DkyJFo7I7KykosUd2Zkjj++OOxZs0aTJ48GZs3b8ZJJ52EwYMHo6WlBZGvhodXXXUV5s6d6/whGSZgiGqKSkpUoUVj6NB2PPBA/5QZpCcSCqnG5EuXWp8bBJs0NwK62pniDcKzM0yQySiN14YNG7Br1y4A8rG7TjvtNGzbtg2zZs3CiSeeiK6uLgwcOBBjxozBsmXL8NJLL+l6G2qMHz8e7777Lq699lqUl5ejvb0dQ4YMwbnnnot169bh8ccfF5r2ZJh0Q3Sdzj17etdbrKvrxh//WIdLLklNoFEjJkwQOy8INmlurDRg9e6c3JdhMhWFyInJK5NqmpqaMHjwYDQ2Nroex+vFF1/E+PHj2Z2acYxm5A3E2xxpHXrieptBrX+RCFBerhqT67WciqIKKjt3pn561K28Gnk12r1eOhDU+scEF5m+OKM0XgzDpIa+ssi508Xe/cStvGrvLhw2Pidoz84wQYYFL4ZhfKGqSo1vpU0nbtyoakfSRejSSCch0q28xr67mTNVezwn12OYTCYYVqsMw2QEsoFXg0pVlWrvFZQgqWa4lVft3Y0dCyxenB7PzjBBhAUvJi1IXPKEG3om1aSTEOl2XtPp2RkmaLDgxQQeoyVPeIkShmEYJt1gwYsJNEYeVdoSJWxXEnxYWxlc+N0wjP+wcT0TWNxY8oRJLbW1akiDceOAyy9X/5aX85p+QYDfDcOkBha8mMDixpInTOowW6CZF1ROLfxuGCZ1sODFBBY3ljxhUgNrK4OL1bshUo/zu2EYb2DBiwksbix5wqQG1lYGF6t3A6jHf/97f/LDMJkGC15MYBFd46+iwt98MdawtjK4iJb5ggU85cgkE4kAmzYBq1apf1kzKg8LXkxgSaflWZh4WFsZXGTKXJsO5s6WAdghwy1Y8GICTTotz8L0wtrK4KK9GxF271anHLmzZdghwz1Y8GICT19Z4y+TYG1lcIl9NyIsWMCdbabDzjLuwoIXkxZoS5RMnar+5Q47+LC2MrhUVQHV1fbTc2ebWbCzjLuw4MUwjGewtjK43Hqr+JSjHtzZZg7sLOMuvGQQwzCewgsqBxNtynHSJPV37DSSouhPK+nBnW3fh51l3IU1XgzDMBmK2XSw6FQkd7Z9H3aWcRcWvBiGYTIYo+ngm28GSkqM03Fnmzmws4y7sODFMAyT4SQ6rzz7LPC1rwH79+ufz51t5sHOMu7BNl4MwzBMFC1ek5mNVzisCl3c2WYWVVXAhAmqQ8Xeveo0c0UFC9+ysODFMAzDADCP16RRUgJ8/DHQv79/+WKCAzvLOIenGhmGYRgAYgto798PvP66P/lhmL4IC14MwzAMAI7XxDB+wIIXwzAMA4DjNTGMH7DgxTAMwwDgeE0M4wcseDEMwzAAOF4Tw/gBC14MwzBMFI7XxDDewuEkGIZhmCiRCFBUBNx5p+rBWFKiCmEcr4lh3IEFL4ZhGAaAGjz1+uvjQ0qEw+r0IwtdDOMOPNXIMAzDRCPWJ8bxqq9X99fWpiZfDNPXYMGLYRgmwzGLWK/tmzlTPY9hGGew4MUwDJPhWEWsJwJ271bPYxjGGSx4MQzDZDgcsZ5h/IMFL4ZhmAyHI9YzjH+w4MUwDJPhnHWWGjbCCI5YzzDuwYIXwzBMBlNbC3zta2rMLj04Yn0ykQiwaROwapX6l50OGBk4jhfDMEyGooWQ0PNm1AiHVaGLI9armMU64zJiRGDBi2EYJgMxCyGhUVICfPwx0L+/f/kKMkaCqhbrjJdUYkTgqUaGYZgMxCqEBKBOP77+uj/5CToc64xxCxa8GIZhMhAOISEHxzpj3IKnGpm0JxJRG7u9e1V3d17Ml2Gs4RAScrCgyrgFa7yYtKa2FigvB8aNAy6/XP1bXs7ryjGMFRUVqlG45rWYCIeQiIcFVcYtWPBi0hZe1Jdh7BMKqZ54QLLwxSEkkmFBlXELFryYtIQNXRnGOVVVqideaWn8/nCYPfQSYUGVcQsWvJi0hA1dGcYdqqqATz8FNm4EVq5U/+7cyUKXHiyoMm7AxvVMWsKGrgzjHqEQMHZsqnORHlRVARMmsEMPY5+M1Hg1NTVh0aJFOOuss1BSUoKcnByEw2GMGzcOCxcuxOHDh3XT7du3D3PmzMGIESOQm5uLoqIiVFRU4JFHHgGZRSH8ik8++QTXXXcdjj/+eAwYMABHHXUUzjvvPNTU1Lj8hH0fNnRlGCZVaILq1KnqXxa6GBkyTuO1ceNGTJ06Ffv27QMA9OvXD/n5+aivr0d9fT02bdqEiy++GN/85jfj0r399ts477zz0NDQAADIz89Hc3MztmzZgi1btmDt2rV47rnnkJOTo3vfF198EZMnT0ZbWxsAoKCgAA0NDVi/fj3Wr1+Pq6++Go8++igUI8tNJg7N0LW+Xt/OS1HU42zoyjAMwwSJjNJ4vfbaa7jwwguxb98+/OAHP8CWLVvQ0dGBQ4cOoa2tDVu3bsWtt96KwYMHx6VrbGzED3/4QzQ0NGDkyJF466230NzcjNbWVtx///3Izs7G+vXrMWvWLN377ty5E5deeina2tpw9tlnY/v27WhsbERjYyPmz58PAHj88cdx1113eV4GfQU2dGUYhmHSEsoQWltb6YQTTiAANHHiRIpEIsJp582bRwAoNzeXduzYkXT89ttvJwAUCoVo+/btScenTZtGAGjYsGF06NChpOPTp08nAFRQUEAHDx6Ueq7GxkYCQI2NjVLprOjs7KRnnnmGOjs7Xb2u29TUEIXDRKreS93KytT9TPqSLvWP6Ztw/WNkkemLM0bj9ec//xk7duxAbm4uHn74YWRliT/68uXLAQBTpkzB8ccfn3T8l7/8JfLz8xGJRPDkk0/GHWttbY3acM2YMQNDhgxJSv/rX/8agGp79swzzwjni2GPLIZhGCa9yBjBSxOeJkyYgKFDhwqn2759Oz777DMAwAUXXKB7Tn5+Piq+MiZav3593LEtW7agvb3dNH15eTlGjRqlm56xhg1dGYZhmHQhIwSvjo4ObN26FQBwzjnnYMeOHfjJT36CcDiMnJwcDBs2DBMmTMBLL72UlPa9996L/n/KKacY3kM79v777xumP/nkky3Tb9u2TeCJGIZhGIZJRzLCq/HTTz9FZ2cnAGDPnj34t3/7N7S2tqJ///7Iy8vDvn378Nxzz+G5557Dz372Mzz00EPRtJ9//nn0/9LEqHkxaMeamprQ0tKC/Pz8uPSFhYXIy8uzTB97Pz06OjrQ0dER/d3U1AQA6OrqQldXl2laGbRruXlNhhGF6x+TSrj+MbLI1JWMELwOHToU/f+OO+5AQUEBVq1ahYkTJyI7Oxu7d+/GjTfeiKeeegoPP/wwRo4cieuvvx4A0NzcHE1rJjjFHmtubo4KXlp6s7Sxx2Pvp8cdd9yB6urqpP3r16+3vIcd6urqXL8mw4jC9Y9JJVz/GFG0UFEiZITg1dPTE/f/ww8/jMsuuyy6r6ysDE8++SS2b9+Od955B7/73e/wi1/8Av36Ba94fv3rX2P27NnR301NTSgrK0NlZSUKCgpcu09XVxfq6upw7rnnIjs727XrMowIXP+YVML1j5FFm30SIXiShQcMGjQo+n9ZWVmc0KWRlZWFOXPmYNq0aThw4ADefvttfO9734tL29bWZijcxEq7sWm0/62kYe14bFo9cnJydIO0Zmdne9JAeHVdhhGB6x+TSrj+MaLI1JOMMK6Ptc0aOXKk4XmaZyEA7Nq1CwBwzDHHRPfV19cbptWOFRQURKcZY9NrQVqt0sfej2EYhmGYvkVGCF5FRUVR4ctsSR6KWXtGOy/WkzHWQzER7djo0aPj9semN/NY1NKbeT4yDMMwDJPeZITgBQCVlZUAgA8++MBwQesPPvgg+r8WKHXEiBE49thjAQAvv/yybrrW1lZs3rw57j4aY8aMQW5urmn6Xbt2Re+dmJ5hGIZhmL5DxgheV199NQBg9+7dWL16ddLxnp4e3HPPPQDUqclvf/vb0WNXXnklAOCpp57Cp59+mpT2gQceQEtLC0KhEK644oq4YwMHDsTEiRMBAA899BAaGxuT0i9atAiAat918cUXyz8cwzAMwzBpQcYIXhUVFZg0aRIAdeme1atXR+Nu7N69G1dccQXeeecdAMDvf//7uCWFbrjhBgwbNgxtbW248MIL8fbbbwMAOjs78dBDD+E3v/kNAGD69Ok46aSTku592223YeDAgdi7dy8uuugifPTRRwBUTdltt92Ghx9+GAAwb948FBYWelQCDMMwDMOkmozwatR44okn8OWXX+LVV1/FlClTkJOTg7y8vLg4X/Pnz8ePf/zjuHSDBw/G888/j/POOw/vv/8+Tj/9dAwaNAhHjhyJCm+VlZVYsmSJ7n2PP/54rFmzBpMnT8bmzZtx0kknYfDgwWhpaUEkEgEAXHXVVZg7d65HT84wDMMwTBDIGI0XoE77bdy4EcuWLcP3v/99DBw4EC0tLSgtLcWUKVPw2muv6QYnBYDTTjsN27Ztw6xZs3DiiSeiq6sLAwcOxJgxY7Bs2TK89NJLumEeNMaPH493330X1157LcrLy9He3o4hQ4bg3HPPxbp16/D444+bGv4zDMMwDJP+ZJTGC1Djdf30pz/FT3/6U+m0Rx99NO65556oLZgsX/va1/CnP/3JVlqGYRiGYdKfjNJ4MQzDMAzDpBIWvBiGYRiGYXyCBS+GYRiGYRifYMGLYRiGYRjGJ1jwYhiGYRiG8QkWvBiGYRiGYXyCBS+GYRiGYRifYMGLYRiGYRjGJ1jwYhiGYRiG8QkWvBiGYRiGYXyCBS+GYRiGYRifYMGLYRiGYRjGJzJukWyGYRg/iUSAzZuBvXuB4cOBigogFEp1rhiGSRUseDEMw3hEbS1w/fXAnj29+8Jh4N57gaqq1OWLYZjUwVONDMMwHlBbC0yaFC90AUB9vbq/tjY1+WIYJrWw4MUwDOMykYiq6SJKPqbtmzlTPY9hmMyCBS+GYRiX2bw5WdMVCxGwe7d6HsMwmQULXgzDMC6zd6+75zEM03dgwYthGMZlhg939zyGYfoOLHgxDMO4TEWF6r2oKPrHFQUoK1PPYxgms2DBi2EYxmVCITVkBJAsfGm/ly7leF4Mk4mw4MUwDOMBVVXAunVAaWn8/nBY3c9xvBgmM+EAqgzDMB5RVQVMmMCR6xmG6YUFL4ZhGA8JhYCxY1OdC4ZhggJPNTIMwzAMw/gEC14MwzAMwzA+wYIXwzAMwzCMT7DgxTAMwzAM4xMseDEMwzAMw/gEC14MwzAMwzA+wYIXwzAMwzCMT7DgxTAMwzAM4xMseDEMwzAMw/gER65Pc4gIANDU1OTqdbu6utDW1oampiZkZ2e7em2GsYLrH5NKuP4xsmh9sNYnm8GCV5rT3NwMACgrK0txThiGYRgms2lubsbgwYNNz1FIRDxjAktPTw8+//xzDBo0CIqiuHbdpqYmlJWVYffu3SgoKHDtugwjAtc/JpVw/WNkISI0NzfjmGOOQVaWuRUXa7zSnKysLITDYc+uX1BQwA0PkzK4/jGphOsfI4OVpkuDjesZhmEYhmF8ggUvhmEYhmEYn2DBi9ElJycHCxYsQE5OTqqzwmQgXP+YVML1j/ESNq5nGIZhGIbxCdZ4MQzDMAzD+AQLXgzDMAzDMD7BghfDMAzDMIxPsODFRGlubsbChQtx6qmnIj8/H4MHD8Z3vvMd3H333ejs7Ex19pgA09DQgMcffxzTpk3D6NGjMXDgQOTk5CAcDuPiiy/G008/bXmNffv2Yc6cORgxYgRyc3NRVFSEiooKPPLII0LLcHzyySe47rrrcPzxx2PAgAE46qijcN5556GmpsaNR2TSjDvvvBOKokQ3M7juMb5CDENEn376KZWXlxMAAkB5eXmUk5MT/f2tb32LDh48mOpsMgGlX79+0boCgAYMGEADBw6M23fBBRdQa2urbvqtW7dScXFx9Nz8/Py4a1ZWVtKRI0cM7//CCy9QXl5e9PyCggLKysqK/r766qupp6fHq8dnAsaHH35IAwYMiKt/RnDdY/yGBS+Guru76dRTTyUANHz4cKqrqyMiokgkQk899RQNGjQo2nEyjB4A6Lvf/S49+OCD9Mknn0T379y5k37yk59EO6Fp06YlpT18+DANGzaMANDIkSPprbfeIiKijo4Ouv/++yk7O5sA0IwZM3TvvWPHjqiQd/bZZ9P27duJiKi5uZnmz58fvfeiRYs8eHImaEQiETr77LMJAJ155pmmghfXPSYVsODF0COPPBJtIF5//fWk4ytXrowe37BhQwpyyASd//mf/zE9ft1110Xr0GeffRZ3bN68eQSAcnNzaceOHUlpb7/9dgJAoVAo2rHFMm3aNAJAw4YNo0OHDiUdnz59elQTwVrbvs/SpUsJAF1xxRW0YMECU8GL6x6TCljwYqiiooIA0Lhx43SP9/T00PHHH08A6Morr/Q5d0xf4O9//3u0A6ytrY07duyxx0anZPRobm6m/Px8AkDz58+PO9bS0kK5ubkEgKqrq3XT79y5M3rvxx57zJ0HYgKJpoEqLi6mL7/80lLw4rrHpAI2rs9w2tra8NprrwEALrjgAt1zFEXB+eefDwBYv369b3lj+g4DBgyI/h+JRKL/b9++HZ999hkA4/qXn5+PiooKAMn1b8uWLWhvbzdNX15ejlGjRummZ/oW1157LVpbW3HPPfegpKTE9Fyue0yqYMErw/nggw/Q09MDADjllFMMz9OOffHFFzh48KAveWP6Dps2bYr+f+qpp0b/f++996L/i9S/999/P25/bPqTTz7ZMv22bdvEMsykHcuWLcNf//pX/OAHP8CVV15peT7XPSZVsOCV4Xz++efR/0tLSw3Piz0Wm4ZhrDh8+DDuuOMOAEBFRQVGjBgRPSZb/5qamtDS0pKUvrCwEHl5eZbpue72Terr6zF37lzk5ubij3/8o1AarntMqmDBK8Npbm6O/m/WeMQei03DMGb09PTgRz/6Efbu3YucnBz813/9V9xxp/VP+98sbexxrrt9k+uuuw6NjY1YuHAhTjjhBKE0XPeYVMGCF8MwnnH99dfj+eefBwA8+OCD+MY3vpHiHDF9jRUrVuCFF17AN7/5TcyePTvV2WEYS1jwynAGDRoU/b+trc3wvNhjsWkYxogbbrgB999/PwBgyZIluOaaa5LOcVr/tP/N0sYe57rbt/jyyy8xc+ZMhEIhLFu2DP369RNOy3WPSRUseGU4xxxzTPT/+vp6w/Nij8WmYRg9brzxRtx9990AgLvuugszZ87UPU+2/hUUFCA/Pz8p/aFDh0w7QC09192+xU033YSGhgZMnz4dI0eOREtLS9wWu9RZ4j6ue0yqYMErwxk1ahSystRqEOulk4h2bNiwYSgqKvIlb0x6MnfuXNx1110AgD/84Q+44YYbDM+N9SYTqX+jR482TG/mNaalN/M+Y9KPnTt3AgAeeughDBo0KGnTnDoARPfdeOONALjuMamDBa8MJy8vD2effTYA4OWXX9Y9h4jwyiuvAAAqKyt9yxuTftxwww1YvHgxAFXomjt3run5I0aMwLHHHgvAuP61trZi8+bNAJLr35gxY5Cbm2uafteuXfjggw900zOZC9c9JmWkOoIrk3q0JYMURaE33ngj6fjq1at5ySDGkjlz5kTryeLFi4XTacu25OXl0c6dO5OOL1q0SGjZluHDh9Phw4eTjs+YMYMA0KBBg3jZlgxDdMkgrnuMn7DgxVBXV1d0kezS0tKocBWJRGjNmjVUUFDAi2Qzptx4443RDu6ee+6RShu7UPHo0aNp69atRKQuVPzggw9S//79hRcqrqiooH/9619EpC7pUl1dTYqi8ELFGYqV4MV1j0kFLHgxRKSuKVZeXh5tpPLy8mjAgAHR39/61rd4xMbosmvXrmg9ycrKoqOPPtp0u+uuu5KusXXrViouLo5eZ9CgQZSdnR39XVlZSUeOHDHMwwsvvEB5eXnR8wcPHkyhUCj6+6qrrqKenh4vi4EJIFaCFxHXPcZ/2MaLAaCuKfbuu+9i/vz5OOWUU6AoCrKzs3Haaadh8eLFeOONN1BYWJjqbDIBRFtySvt/3759plts9G+N0047Ddu2bcOsWbNw4oknoqurCwMHDsSYMWOwbNkyvPTSS8jJyTHMw/jx4/Huu+/i2muvRXl5Odrb2zFkyBCce+65WLduHR5//HEoiuLJ8zPpDdc9xm8UIqJUZ4JhGIZhGCYTYI0XwzAMwzCMT7DgxTAMwzAM4xMseDEMwzAMw/gEC14MwzAMwzA+wYIXwzAMwzCMT7DgxTAMwzAM4xMseDEMwzAMw/gEC14MwzAMwzA+wYIXwzAMwzCMT7DgxTAMwzAM4xMseDEMwzAMw/gEC14MwzAMwzA+wYIXwzAMwzCMT7DgxTAMwzAM4xMseDEMwzAMw/gEC14MwzAmjB07FoqiYOHChejq6sLdd9+N008/HUOGDIGiKNi0aRMWLlwIRVEwduxYw+ts2rQJiqJAUZSkY4np//rXv+LCCy9ESUkJBgwYgFGjRqG6uhpHjhwxvP4rr7yCqqoqhMNh9O/fHwUFBTjhhBNQWVmJxYsX4+DBg06LgmEYF+iX6gwwDMOkA0eOHMHYsWPx+uuvo1+/fhg0aJAn97nrrrtw0003AQAGDx6Mzs5OfPjhh1i4cCH+93//F3V1dQiFQnFpbrvtNixYsCD6Oy8vD0SEnTt3YufOnairq8Ppp59uKhgyDOMPrPFiGIYR4IEHHsC7776Lxx9/HE1NTTh48CAOHDiAf/u3f3PtHv/85z9x88034+abb8aXX36JQ4cO4fDhw5g/fz4AYOPGjfjv//7vuDS7du1CdXU1AGD27Nmor69Ha2srmpubcfjwYWzevBk///nPPRMUGYaRgzVeDMMwArS0tOC5557DRRddFN1XXFzs6j0OHz6MBQsWYOHChdF9BQUFqK6uxnvvvYfa2lqsWrUK11xzTfT4m2++iZ6eHpx00km4++674643ePBgjBkzBmPGjHE1nwzD2Ic1XgzDMAKcfPLJcUKXF+Tk5OCGG27QPTZhwgQAwLvvvhu3f8iQIQCA5uZmtLa2epo/hmGcw4IXwzCMAGeffbbn9zj55JORn5+ve+yYY44BgCQj+e9+97sYOnQo9u7di+9973u4//778eGHH4KIPM8vwzDysODFMAwjwFFHHeX5PczssPr1Uy1Duru74/YPGTIEq1atQklJCbZt24Zf/vKXGDVqFAoLC/Ef//EfWLFiBbq6ujzNN8Mw4rDgxTAMI0CiJ2GQ+MEPfoCdO3di+fLl+PGPf4wTTzwRjY2N+Mtf/oIf/ehH+Na3voX6+vpUZ5NhGLDgxTAM4xhNG2UWZ6uxsdHTPAwcOBA/+tGP8MQTT+Bf//oX9uzZg0WLFmHAgAFRTRjDMKmHBS+GYRiHFBYWAgB2795teM6bb77pV3YAAKWlpbjxxhsxZ84cAEBdXZ2v92cYRh8WvBiGYRzyjW98AwDw+eef44033kg6/uWXX2LZsmWe3Lujo8P0eG5uLoBgT5UyTCbBghfDMIxDzjrrLBx33HEAgKuuugpbt24FEaGnpwebNm3C2LFj0dPT48m9Fy1ahAsuuAB//vOfsWfPnuj+jo4OrFmzBnfddRcAYPz48Z7cn2EYOTiAKsMwjEOysrLwxz/+ERdddBG2b9+O73znO8jLy0NPTw+OHDmCE088EQ888ACmTp3q+r17enrw8ssv4+WXXwagarhyc3Nx6NChaEiJUaNG4Z577nH93gzDyMMaL4ZhGBc477zzsHnzZvzwhz9EYWEhIpEIysrKcPPNN+Ptt9/GsGHDPLnv9OnT8ac//QlTp07FKaecgry8PDQ1NaGwsBAVFRVYunQp/vGPf3h2f4Zh5FCIo+wxDMMwDMP4Amu8GIZhGIZhfIIFL4ZhGIZhGJ9gwYthGIZhGMYnWPBiGIZhGIbxCRa8GIZhGIZhfIIFL4ZhGIZhGJ9gwYthGIZhGMYnWPBiGIZhGIbxCRa8GIZhGIZhfIIFL4ZhGIZhGJ9gwYthGIZhGMYnWPBiGIZhGIbxCRa8GIZhGIZhfIIFL4ZhGIZhGJ/4/wGJWoqjT/FUOQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "shotNum = \"0007\"\n", + "filePath = folderPath + \"/\" + shotNum + \"/*.h5\"\n", + "\n", + "dataSetOfGlobalDict = {\n", + " dskey[groupList[i]]: read_hdf5_global(filePath, groupList[i])\n", + " for i in [0]\n", + "}\n", + "\n", + "dataSetDict = {\n", + " dskey[groupList[i]]: read_hdf5_file(filePath, groupList[i], datesetOfGlobal=dataSetOfGlobalDict[dskey[groupList[i]]])\n", + " for i in [0]\n", + "}\n", + "\n", + "dataSet = dataSetDict[\"camera_0\"]\n", + "\n", + "print_scanAxis(dataSet)\n", + "\n", + "scanAxis = get_scanAxis(dataSet)\n", + "\n", + "dataSet = auto_rechunk(dataSet)\n", + "dataSet = swap_xy(dataSet)\n", + "\n", + "dataSet = imageAnalyser.get_absorption_images(dataSet)\n", + "\n", + "imageAnalyser.center = (959, 876)\n", + "imageAnalyser.span = (100, 100)\n", + "imageAnalyser.fraction = (0.1, 0.1)\n", + "\n", + "dataSet_cropOD = imageAnalyser.crop_image(dataSet.OD)\n", + "dataSet_cropOD = imageAnalyser.substract_offset(dataSet_cropOD).load()\n", + "\n", + "Ncount = imageAnalyser.get_Ncount(dataSet_cropOD).load()\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "Ncount.plot.errorbar(ax=ax, fmt='ob')\n", + "\n", + "plt.ylabel('NCount')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dataSet_cropOD = auto_rechunk(dataSet_cropOD)\n", + "\n", + "fitAnalyser = FitAnalyser(\"Two Gaussian-2D\", fitDim=2)\n", + "params = fitAnalyser.guess(dataSet_cropOD, dask=\"parallelized\")\n", + "fitResult = fitAnalyser.fit(dataSet_cropOD, params, dask=\"parallelized\").load()\n", + "\n", + "fitValue = fitAnalyser.get_fit_value(fitResult)\n", + "fitStd = fitAnalyser.get_fit_std(fitResult)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "BEC_Ncount_val = fitValue['A_amplitude']\n", + "BEC_Ncount_std = fitStd['A_amplitude']\n", + "\n", + "thermal_Ncount_val = fitValue['B_amplitude']\n", + "thermal_Ncount_std = fitStd['B_amplitude']\n", + "\n", + "BEC_width_x_val = fitValue['A_sigmax']\n", + "BEC_width_x_std = fitStd['A_sigmax']\n", + "BEC_width_y_val = fitValue['A_sigmay']\n", + "BEC_width_y_std = fitStd['A_sigmay']\n", + "\n", + "thermal_width_x_val = fitValue['B_sigmax']\n", + "thermal_width_x_std = fitStd['B_sigmax']\n", + "thermal_width_y_val = fitValue['B_sigmay']\n", + "thermal_width_y_std = fitStd['B_sigmay']\n", + "\n", + "BEC_center_x_val = fitValue['A_centerx']\n", + "BEC_center_x_std = fitStd['A_centerx']\n", + "BEC_center_y_val = fitValue['A_centery']\n", + "BEC_center_y_std = fitStd['A_centery']\n", + "\n", + "thermal_center_x_val = fitValue['B_centerx']\n", + "thermal_center_x_std = fitStd['B_centerx']\n", + "thermal_center_y_val = fitValue['B_centery']\n", + "thermal_center_y_std = fitStd['B_centery']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "total_Ncount_val = BEC_Ncount_val + thermal_Ncount_val\n", + "total_Ncount_std = BEC_Ncount_std + thermal_Ncount_std\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "total_Ncount_val.plot.errorbar(ax=ax, yerr=total_Ncount_std, fmt='ob')\n", + "# plt.ylim([0, 1100])\n", + "plt.ylabel('Ncount from fit')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()\n", + "\n", + "print(total_Ncount_val.mean())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "BEC_Ncount_val.plot.errorbar(ax=ax, yerr=BEC_Ncount_std, fmt='ob')\n", + "plt.ylim([0, 750])\n", + "plt.ylabel('Ncount of BEC part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_Ncount_val.plot.errorbar(ax=ax, yerr=thermal_Ncount_std, fmt='or')\n", + "plt.ylim([0, 500])\n", + "plt.ylabel('Ncount of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "BEC_width_x_val.plot.errorbar(ax=ax, yerr=BEC_width_x_std, fmt='ob')\n", + "\n", + "plt.ylabel('X-axis Width of BEC part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "BEC_width_y_val.plot.errorbar(ax=ax, yerr=BEC_width_y_std, fmt='ob')\n", + "\n", + "plt.ylabel('Y-axis Width of BEC part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_width_x_val.plot.errorbar(ax=ax, yerr=thermal_width_x_std, fmt='or')\n", + "\n", + "plt.ylabel('X-axis width of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_width_y_val.plot.errorbar(ax=ax, yerr=thermal_width_y_std, fmt='or')\n", + "\n", + "plt.ylabel('Y-axis width of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "BEC_center_x_val.plot.errorbar(ax=ax, yerr=BEC_center_x_std, fmt='ob')\n", + "\n", + "plt.ylabel('X-axis center of BEC part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "BEC_center_y_val.plot.errorbar(ax=ax, yerr=BEC_center_y_std, fmt='ob')\n", + "\n", + "plt.ylabel('Y-axis center of BEC part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_center_x_val.plot.errorbar(ax=ax, yerr=thermal_center_x_std, fmt='or')\n", + "\n", + "plt.ylabel('X-axis center of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_center_y_val.plot.errorbar(ax=ax, yerr=thermal_center_y_std, fmt='or')\n", + "\n", + "plt.ylabel('Y-axis center of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fitFullResult = fitAnalyser.get_fit_full_result(fitResult)\n", + "condensateFraction = fitFullResult['A_amplitude'] / (fitFullResult['A_amplitude'] + fitFullResult['B_amplitude'])\n", + "condensateFraction_value, condensateFraction_std = seperate_uncertainty(condensateFraction)\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "condensateFraction_value.plot.errorbar(ax=ax, yerr=condensateFraction_std, fmt='ob')\n", + "\n", + "plt.ylabel('Condensate Fraction')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "val = Ncount.mean().item()\n", + "std = Ncount.std().item()\n", + "print(f'The total Ncount is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = total_Ncount_val.mean().item()\n", + "std = total_Ncount_val.std().item()\n", + "print(f'The total Ncount from fit is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = BEC_Ncount_val.mean().item()\n", + "std = BEC_Ncount_val.std().item()\n", + "print(f'The Ncount of the BEC part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_Ncount_val.mean().item()\n", + "std = thermal_Ncount_val.std().item()\n", + "print(f'The Ncount of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = BEC_width_x_val.mean().item()\n", + "std = BEC_width_x_val.std().item()\n", + "print(f'The x-axis width of the BEC part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = BEC_width_y_val.mean().item()\n", + "std = BEC_width_y_val.std().item()\n", + "print(f'The y-axis width of the BEC part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_width_x_val.mean().item()\n", + "std = thermal_width_x_val.std().item()\n", + "print(f'The x-axis width of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_width_y_val.mean().item()\n", + "std = thermal_width_y_val.std().item()\n", + "print(f'The y-axis width of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = BEC_center_x_val.mean().item()\n", + "std = BEC_center_x_val.std().item()\n", + "print(f'The x-axis center of the BEC part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = BEC_center_y_val.mean().item()\n", + "std = BEC_center_y_val.std().item()\n", + "print(f'The y-axis center of the BEC part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_center_x_val.mean().item()\n", + "std = thermal_center_x_val.std().item()\n", + "print(f'The x-axis center of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_center_y_val.mean().item()\n", + "std = thermal_center_y_val.std().item()\n", + "print(f'The y-axis center of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = condensateFraction_value.mean().item()\n", + "std = condensateFraction_value.std().item()\n", + "print(f'The condensate fraction is: {val: .4f} \\u00B1 {std: .4f}')" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "i=0\n", + "da_test = read_hdf5_run_time(filePath, datesetOfGlobal=dataSetOfGlobalDict[dskey[groupList[i]]])" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:  (runs: 550)\n",
+       "Coordinates:\n",
+       "  * runs     (runs) float64 0.0 1.0 2.0 3.0 4.0 ... 546.0 547.0 548.0 549.0\n",
+       "Data variables:\n",
+       "    runTine  (runs) datetime64[ns] 2023-05-09T14:30:03 ... 2023-05-09T15:56:53\n",
+       "Attributes: (12/101)\n",
+       "    TOF_free:                          0.02\n",
+       "    abs_img_freq:                      110.858\n",
+       "    absorption_imaging_flag:           True\n",
+       "    backup_data:                       True\n",
+       "    blink_off_time:                    nan\n",
+       "    blink_on_time:                     nan\n",
+       "    ...                                ...\n",
+       "    y_offset_img:                      0\n",
+       "    z_offset:                          0.189\n",
+       "    z_offset_img:                      0.189\n",
+       "    runs:                              [  0.   1.   2.   3.   4.   5.   6.   ...\n",
+       "    scanAxis:                          ['runs']\n",
+       "    scanAxisLength:                    [550.]
" + ], + "text/plain": [ + "\n", + "Dimensions: (runs: 550)\n", + "Coordinates:\n", + " * runs (runs) float64 0.0 1.0 2.0 3.0 4.0 ... 546.0 547.0 548.0 549.0\n", + "Data variables:\n", + " runTine (runs) datetime64[ns] 2023-05-09T14:30:03 ... 2023-05-09T15:56:53\n", + "Attributes: (12/101)\n", + " TOF_free: 0.02\n", + " abs_img_freq: 110.858\n", + " absorption_imaging_flag: True\n", + " backup_data: True\n", + " blink_off_time: nan\n", + " blink_on_time: nan\n", + " ... ...\n", + " y_offset_img: 0\n", + " z_offset: 0.189\n", + " z_offset_img: 0.189\n", + " runs: [ 0. 1. 2. 3. 4. 5. 6. ...\n", + " scanAxis: ['runs']\n", + " scanAxisLength: [550.]" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "da_test" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "filepath = '//DyLabNAS/Data/Evaporative_Cooling/2023/05/09/0007/2023-05-09_0007_Evaporative_Cooling_549.h5'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Close to the BEC transition point, in evaporative cooling 2 with truncation value = 0.77" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "shotNum = \"0015\"\n", + "filePath = folderPath + \"/\" + shotNum + \"/*.h5\"\n", + "\n", + "dataSetDict = {\n", + " dskey[groupList[i]]: read_hdf5_file(filePath, groupList[i])\n", + " for i in [0]\n", + "}\n", + "\n", + "dataSet = dataSetDict[\"camera_0\"]\n", + "\n", + "print_scanAxis(dataSet)\n", + "\n", + "scanAxis = get_scanAxis(dataSet)\n", + "\n", + "dataSet = auto_rechunk(dataSet)\n", + "\n", + "dataSet = imageAnalyser.get_absorption_images(dataSet)\n", + "\n", + "imageAnalyser.center = (879, 956)\n", + "imageAnalyser.span = (200, 200)\n", + "imageAnalyser.fraction = (0.1, 0.1)\n", + "\n", + "dataSet_cropOD = imageAnalyser.crop_image(dataSet.OD)\n", + "dataSet_cropOD = imageAnalyser.substract_offset(dataSet_cropOD).load()\n", + "\n", + "Ncount = imageAnalyser.get_Ncount(dataSet_cropOD).load()\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "Ncount.plot.errorbar(ax=ax, fmt='ob')\n", + "\n", + "plt.ylabel('NCount')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "Ncount.plot.errorbar(ax=ax, fmt='ob')\n", + "plt.ylim([0, 3000])\n", + "plt.ylabel('NCount')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dataSet_cropOD = auto_rechunk(dataSet_cropOD)\n", + "\n", + "fitAnalyser = FitAnalyser(\"Gaussian-2D\", fitDim=2)\n", + "params = fitAnalyser.guess(dataSet_cropOD, dask=\"parallelized\")\n", + "fitResult = fitAnalyser.fit(dataSet_cropOD, params, dask=\"parallelized\").load()\n", + "\n", + "fitValue = fitAnalyser.get_fit_value(fitResult)\n", + "fitStd = fitAnalyser.get_fit_std(fitResult)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "thermal_Ncount_val = fitValue['amplitude']\n", + "thermal_Ncount_std = fitStd['amplitude']\n", + "\n", + "thermal_width_x_val = fitValue['sigmax']\n", + "thermal_width_x_std = fitStd['sigmax']\n", + "thermal_width_y_val = fitValue['sigmay']\n", + "thermal_width_y_std = fitStd['sigmay']\n", + "\n", + "thermal_center_x_val = fitValue['centerx']\n", + "thermal_center_x_std = fitStd['centerx']\n", + "thermal_center_y_val = fitValue['centery']\n", + "thermal_center_y_std = fitStd['centery']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "total_Ncount_val = thermal_Ncount_val\n", + "total_Ncount_std = thermal_Ncount_std\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "total_Ncount_val.plot.errorbar(ax=ax, yerr=total_Ncount_std, fmt='ob')\n", + "plt.ylim([0, 3000])\n", + "plt.ylabel('Ncount from fit')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_width_x_val.plot.errorbar(ax=ax, yerr=thermal_width_x_std, fmt='or')\n", + "\n", + "plt.ylabel('X-axis width of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_width_y_val.plot.errorbar(ax=ax, yerr=thermal_width_y_std, fmt='or')\n", + "\n", + "plt.ylabel('Y-axis width of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_center_x_val.plot.errorbar(ax=ax, yerr=thermal_center_x_std, fmt='or')\n", + "\n", + "plt.ylabel('X-axis center of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_center_y_val.plot.errorbar(ax=ax, yerr=thermal_center_y_std, fmt='or')\n", + "\n", + "plt.ylabel('Y-axis center of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "val = Ncount.mean().item()\n", + "std = Ncount.std().item()\n", + "print(f'The total Ncount is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = total_Ncount_val.mean().item()\n", + "std = total_Ncount_val.std().item()\n", + "print(f'The total Ncount from fit is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_width_x_val.mean().item()\n", + "std = thermal_width_x_val.std().item()\n", + "print(f'The x-axis width of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_width_y_val.mean().item()\n", + "std = thermal_width_y_val.std().item()\n", + "print(f'The y-axis width of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_center_x_val.mean().item()\n", + "std = thermal_center_x_val.std().item()\n", + "print(f'The x-axis center of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_center_y_val.mean().item()\n", + "std = thermal_center_y_val.std().item()\n", + "print(f'The y-axis center of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## At the end of ODT loading" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "shotNum = \"0020\"\n", + "filePath = folderPath + \"/\" + shotNum + \"/*.h5\"\n", + "\n", + "dataSetDict = {\n", + " dskey[groupList[i]]: read_hdf5_file(filePath, groupList[i])\n", + " for i in [0]\n", + "}\n", + "\n", + "dataSet = dataSetDict[\"camera_0\"]\n", + "\n", + "print_scanAxis(dataSet)\n", + "\n", + "scanAxis = get_scanAxis(dataSet)\n", + "\n", + "dataSet = auto_rechunk(dataSet)\n", + "\n", + "dataSet = imageAnalyser.get_absorption_images(dataSet)\n", + "\n", + "imageAnalyser.center = (550, 800)\n", + "imageAnalyser.span = (900, 1600)\n", + "imageAnalyser.fraction = (0.1, 0.1)\n", + "\n", + "dataSet_cropOD = imageAnalyser.crop_image(dataSet.OD)\n", + "dataSet_cropOD = imageAnalyser.substract_offset(dataSet_cropOD).load()\n", + "\n", + "Ncount = imageAnalyser.get_Ncount(dataSet_cropOD).load()\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "Ncount.plot.errorbar(ax=ax, fmt='ob')\n", + "\n", + "plt.ylabel('NCount')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dataSet" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "Ncount.plot.errorbar(ax=ax, fmt='ob')\n", + "plt.ylim([0, 150000])\n", + "plt.ylabel('NCount')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dataSet_cropOD = dataSet_cropOD.chunk((1, 900, 1600))\n", + "dataSet_cropOD" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "# dataSet_cropOD = auto_rechunk(dataSet_cropOD)\n", + "\n", + "fitAnalyser = FitAnalyser(\"Gaussian-2D\", fitDim=2)\n", + "params = fitAnalyser.guess(dataSet_cropOD, dask=\"parallelized\")\n", + "fitResult = fitAnalyser.fit(dataSet_cropOD, params, dask=\"parallelized\").load()\n", + "\n", + "fitValue = fitAnalyser.get_fit_value(fitResult)\n", + "fitStd = fitAnalyser.get_fit_std(fitResult)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "thermal_Ncount_val = fitValue['amplitude']\n", + "thermal_Ncount_std = fitStd['amplitude']\n", + "\n", + "thermal_width_x_val = fitValue['sigmax']\n", + "thermal_width_x_std = fitStd['sigmax']\n", + "thermal_width_y_val = fitValue['sigmay']\n", + "thermal_width_y_std = fitStd['sigmay']\n", + "\n", + "thermal_center_x_val = fitValue['centerx']\n", + "thermal_center_x_std = fitStd['centerx']\n", + "thermal_center_y_val = fitValue['centery']\n", + "thermal_center_y_std = fitStd['centery']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "total_Ncount_val = thermal_Ncount_val\n", + "total_Ncount_std = thermal_Ncount_std\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "total_Ncount_val.plot.errorbar(ax=ax, yerr=total_Ncount_std, fmt='ob')\n", + "plt.ylim([0, 160000])\n", + "plt.ylabel('Ncount from fit')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_width_x_val.plot.errorbar(ax=ax, yerr=thermal_width_x_std, fmt='or')\n", + "\n", + "plt.ylabel('Y-axis width of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_width_y_val.plot.errorbar(ax=ax, yerr=thermal_width_y_std, fmt='or')\n", + "\n", + "plt.ylabel('X-axis width of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_center_x_val.plot.errorbar(ax=ax, yerr=thermal_center_x_std, fmt='or')\n", + "\n", + "plt.ylabel('Y-axis center of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "thermal_center_y_val.plot.errorbar(ax=ax, yerr=thermal_center_y_std, fmt='or')\n", + "\n", + "plt.ylabel('X-axis center of thermal part')\n", + "plt.tight_layout()\n", + "plt.grid(visible=1)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "val = Ncount.mean().item()\n", + "std = Ncount.std().item()\n", + "print(f'The total Ncount is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = total_Ncount_val.mean().item()\n", + "std = total_Ncount_val.std().item()\n", + "print(f'The total Ncount from fit is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_width_x_val.mean().item()\n", + "std = thermal_width_x_val.std().item()\n", + "print(f'The y-axis width of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_width_y_val.mean().item()\n", + "std = thermal_width_y_val.std().item()\n", + "print(f'The x-axis width of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_center_x_val.mean().item()\n", + "std = thermal_center_x_val.std().item()\n", + "print(f'The y-axis center of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')\n", + "\n", + "val = thermal_center_y_val.mean().item()\n", + "std = thermal_center_y_val.std().item()\n", + "print(f'The x-axis center of the thermal part is: {val: .2f} \\u00B1 {std: .2f}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "l = list(np.arange(0.001, 0.025, 0.0005))\n", + "# l = np.logspace(np.log10(100e-3), np.log10(20), num=20)\n", + "\n", + "l = [round(item, 7) for item in l]\n", + "#random.shuffle(l)\n", + "\n", + "print(l)\n", + "print(len(l))\n", + "np.mean(l)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## ODT 1 Calibration" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "v_high = 2.7\n", + "\"\"\"High Power\"\"\"\n", + "P_arm1_high = 5.776 * v_high - 0.683\n", + "\n", + "v_mid = 0.2076\n", + "\"\"\"Intermediate Power\"\"\"\n", + "P_arm1_mid = 5.815 * v_mid - 0.03651\n", + "\n", + "v_low = 0.0587\n", + "\"\"\"Low Power\"\"\"\n", + "P_arm1_low = 5271 * v_low - 27.5\n", + "\n", + "print(round(P_arm1_high, 3))\n", + "print(round(P_arm1_mid, 3))\n", + "print(round(P_arm1_low, 3))" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## ODT 2 Power Calibration" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "v = 0.7607\n", + "P_arm2 = 2.302 * v - 0.06452\n", + "print(round(P_arm2, 3))" + ] + } + ], + "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.13" + }, + "vscode": { + "interpreter": { + "hash": "c05913ad4f24fdc6b2418069394dc5835b1981849b107c9ba6df693aafd66650" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/DataContainer/ReadData.py b/DataContainer/ReadData.py index 79e8fc3..3e587be 100644 --- a/DataContainer/ReadData.py +++ b/DataContainer/ReadData.py @@ -5,6 +5,7 @@ from functools import partial import copy import glob import os +from datetime import datetime def _read_globals_attrs(variable_attrs, context=None): @@ -185,4 +186,139 @@ def read_hdf5_file(filePath, group=None, datesetOfGlobal=None, preprocess=None, ds.attrs = copy.deepcopy(datesetOfGlobal.attrs) - return ds \ No newline at end of file + return ds + +def _assign_scan_axis_partial_and_remove_everything(x, datesetOfGlobal, fullFilePath): + scanAxis = datesetOfGlobal.scanAxis + filePath = x.encoding["source"].replace("\\", "/") + shotNum = np.where(fullFilePath==filePath) + shotNum = np.squeeze(shotNum) + runTime = _read_run_time_from_hdf5(x) + x = xr.Dataset(data_vars={'runTine':runTime}) + x = x.expand_dims(list(scanAxis)) + + return x.assign_coords( + { + key: np.atleast_1d(np.atleast_1d(datesetOfGlobal.attrs[key])[int(shotNum)]) + for key in scanAxis + } + ) + +def _read_run_time_from_hdf5(x): + runTime = datetime.strptime(x.attrs['run time'], '%Y%m%dT%H%M%S') + return runTime + +def read_hdf5_run_time(filePath, group=None, datesetOfGlobal=None, preprocess=None, join="outer", parallel=True, engine="h5netcdf", phony_dims="access", excludeAxis=[], maxFileNum=None, **kwargs): + + filePath = np.sort(np.atleast_1d(filePath)) + + filePathAbs = [] + + for i in range(len(filePath)): + filePathAbs.append(os.path.abspath(filePath[i]).replace("\\", "/")) + + fullFilePath = [] + for i in range(len(filePathAbs)): + fullFilePath.append(list(np.sort(glob.glob(filePathAbs[i])))) + fullFilePath = np.array(fullFilePath).flatten() + + for i in range(len(fullFilePath)): + fullFilePath[i] = fullFilePath[i].replace("\\", "/") + + if not maxFileNum is None: + fullFilePath = fullFilePath[0:int(maxFileNum)] + + kwargs.update( + { + 'join': join, + 'parallel': parallel, + 'engine': engine, + 'phony_dims': phony_dims, + 'group': group + } + ) + + if datesetOfGlobal is None: + datesetOfGlobal = xr.open_mfdataset( + fullFilePath, + group="globals", + concat_dim="fileNum", + combine="nested", + preprocess=_read_shot_number_from_hdf5, + engine="h5netcdf", + phony_dims="access", + combine_attrs=_read_globals_attrs, + parallel=True, ) + + datesetOfGlobal.attrs['scanAxis'] = np.setdiff1d(datesetOfGlobal.attrs['scanAxis'], excludeAxis) + + _assgin_scan_axis = partial(_assign_scan_axis_partial_and_remove_everything, datesetOfGlobal=datesetOfGlobal, fullFilePath=fullFilePath) + + if preprocess is None: + kwargs.update({'preprocess':_assgin_scan_axis}) + else: + kwargs.update({'preprocess':preprocess}) + + ds = xr.open_mfdataset(fullFilePath, **kwargs) + + newDimKey = np.append(['x', 'y', 'z'], [ chr(i) for i in range(97, 97+23)]) + + oldDimKey = np.sort( + [ + key + for key in ds.dims + if not key in datesetOfGlobal.scanAxis + ] + ) + + renameDict = { + oldDimKey[j]: newDimKey[j] + for j in range(len(oldDimKey)) + } + + ds = ds.rename_dims(renameDict) + + ds.attrs = copy.deepcopy(datesetOfGlobal.attrs) + + return ds + +def read_hdf5_global(filePath, preprocess=None, join="outer", combine="nested", parallel=True, engine="h5netcdf", phony_dims="access", excludeAxis=[], maxFileNum=None, **kwargs): + + filePath = np.sort(np.atleast_1d(filePath)) + + filePathAbs = [] + + for i in range(len(filePath)): + filePathAbs.append(os.path.abspath(filePath[i]).replace("\\", "/")) + + fullFilePath = [] + for i in range(len(filePathAbs)): + fullFilePath.append(list(np.sort(glob.glob(filePathAbs[i])))) + fullFilePath = np.array(fullFilePath).flatten() + + for i in range(len(fullFilePath)): + fullFilePath[i] = fullFilePath[i].replace("\\", "/") + + if not maxFileNum is None: + fullFilePath = fullFilePath[0:int(maxFileNum)] + + kwargs.update( + { + 'join': join, + 'parallel': parallel, + 'engine': engine, + 'phony_dims': phony_dims, + 'group': "globals", + 'preprocess': _read_shot_number_from_hdf5, + 'combine_attrs': _read_globals_attrs, + 'combine':combine, + 'concat_dim': "fileNum", + } + ) + + datesetOfGlobal = xr.open_mfdataset(fullFilePath, **kwargs) + + datesetOfGlobal.attrs['scanAxis'] = np.setdiff1d(datesetOfGlobal.attrs['scanAxis'], excludeAxis) + + return datesetOfGlobal + \ No newline at end of file