SRIM_handler.ipynb 85.9 KB
Newer Older
Parish, Chad's avatar
Parish, Chad committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rule #1: generate figures programmatically\n",
    "<br>\n",
    "This parses the annoyingly formated .txt files output by SRIM (http://www.srim.org/), uses the parsed data and minimal user input, and generates well-formatted figures of the damage in dpa (displacements per atom) and the ion implantation level in the material identified by parsing the input files.\n",
    "<br>\n",
    "<br>\n",
    "The principle illustrated by this notebook is that a small investment in time can result in an easily reusable code that will make the subsequent generation of figures fast and efficient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
Parish, Chad's avatar
'minor'    
Parish, Chad committed
25
26
27
28
29
30
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T22:44:43.359243Z",
     "start_time": "2020-05-12T22:44:43.047052Z"
    }
   },
Parish, Chad's avatar
Parish, Chad committed
31
32
33
34
35
36
37
38
39
   "outputs": [],
   "source": [
    "import tkinter\n",
    "from tkinter import filedialog\n",
    "import os\n",
    "import glob\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import gridspec\n",
Parish, Chad's avatar
'minor'    
Parish, Chad committed
40
41
42
    "from warnings import warn\n",
    "\n",
    "plt.style.use({'font.sans-serif': 'arial', 'font.size': 16, 'font.weight': 'bold', 'svg.fonttype': 'none'})"
Parish, Chad's avatar
Parish, Chad committed
43
44
45
46
47
48
49
50
51
52
53
54
55
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find and read in the needed text files\n",
    "In this case, PHONON.txt, RANGE.txt, and TDATA.txt are needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
Parish, Chad's avatar
'minor'    
Parish, Chad committed
56
57
58
59
60
61
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T22:44:48.958171Z",
     "start_time": "2020-05-12T22:44:43.360177Z"
    }
   },
Parish, Chad's avatar
Parish, Chad committed
62
63
64
65
66
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
Parish, Chad's avatar
'minor'    
Parish, Chad committed
67
      "Path: C:/work/Collaborations/isotope_production/30MeV_He_on_Bi\n",
Parish, Chad's avatar
Parish, Chad committed
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
      "Lines in PHONON.txt: 124\n",
      "Lines in TDATA.txt: 39\n",
      "Lines in RANGE.txt: 137\n"
     ]
    }
   ],
   "source": [
    "# get the data directory\n",
    "root = tkinter.Tk()\n",
    "path_dir = filedialog.askdirectory(\n",
    "    initialdir=\"/\",\n",
    "    title=\"Find SRIM data:\")\n",
    "print(f'Path: {path_dir}')\n",
    "root.destroy()  # remove the GUI panel\n",
    "os.chdir(path_dir)\n",
    "\n",
    "# find PHONON.txt\n",
    "fn = glob.glob('PHONON.txt')\n",
    "assert len(fn) == 1, 'Did not find exactly 1 PHONON.txt!'\n",
    "file_lines_phonon = []\n",
    "with open(fn[0], mode='r') as fh:\n",
    "    for line in fh:\n",
    "        file_lines_phonon.append(line)\n",
    "    else:\n",
    "        pass\n",
    "print(f'Lines in PHONON.txt: {len(file_lines_phonon)}')\n",
    "\n",
    "# find TDATA.txt\n",
    "fn = glob.glob('TDATA.txt')\n",
    "assert len(fn) == 1, 'Did not find exactly 1 TDATA.txt!'\n",
    "file_lines_tdata = []\n",
    "with open(fn[0], mode='r') as fh:\n",
    "    for line in fh:\n",
    "        file_lines_tdata.append(line)\n",
    "    else:\n",
    "        pass\n",
    "print(f'Lines in TDATA.txt: {len(file_lines_tdata)}')\n",
    "\n",
    "# find RANGE.txt\n",
    "fn = glob.glob('RANGE.txt')\n",
    "assert len(fn) == 1, 'Did not find exactly 1 RANGE.txt!'\n",
    "file_lines_range = []\n",
    "with open(fn[0], mode='r') as fh:\n",
    "    for line in fh:\n",
    "        file_lines_range.append(line)\n",
    "    else:\n",
    "        pass\n",
    "print(f'Lines in RANGE.txt: {len(file_lines_range)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parse out the important information from PHONON.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
Parish, Chad's avatar
'minor'    
Parish, Chad committed
128
129
130
131
132
133
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T22:44:48.978112Z",
     "start_time": "2020-05-12T22:44:48.961155Z"
    }
   },
Parish, Chad's avatar
Parish, Chad committed
134
135
136
137
138
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
Parish, Chad's avatar
'minor'    
Parish, Chad committed
139
      "Density: 2.82e+22 atoms/cm^3\n"
Parish, Chad's avatar
Parish, Chad committed
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
     ]
    }
   ],
   "source": [
    "#Find the important lines:\n",
    "density_found = 0\n",
    "for i, l in enumerate(file_lines_phonon):\n",
    "    if 'DEPTH' in l:  # line that defines the data columns\n",
    "        i_start = i + 3  # find the first line of data\n",
    "    if 'atoms/cm3' in l:  # find the density\n",
    "        i_density = i\n",
    "        density_found += 1\n",
    "        assert density_found <=1, 'Can currently only handle single layers'\n",
    "\n",
    "# Calculate the length of the data. SRIM is always 100\n",
    "length = len(file_lines_phonon) - i_start\n",
    "data_phonon = np.zeros((length, 3), dtype='f4')  # hold the data here\n",
    "\n",
    "# parse out the density (atoms/cm^3)\n",
    "density_line = file_lines_phonon[i_density]\n",
    "start_phrase = 'Density ='\n",
    "start_phrase_length = len(start_phrase)\n",
    "end_phrase = 'atoms/cm3'\n",
    "st = density_line.find(start_phrase) + start_phrase_length\n",
    "en = density_line.find(end_phrase)\n",
    "density = np.array(density_line[st:en]).astype('f4')\n",
    "print(f'Density: {density:3.3} atoms/cm^3')\n",
    "\n",
    "# read the data\n",
    "for i in range(i_start, len(file_lines_phonon),1):\n",
    "    data_phonon[i-i_start,:] = np.array(file_lines_phonon[i].split()).astype('f4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parse out the important information from TDATA.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
Parish, Chad's avatar
'minor'    
Parish, Chad committed
183
184
185
186
187
188
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T22:44:48.998056Z",
     "start_time": "2020-05-12T22:44:48.981118Z"
    }
   },
Parish, Chad's avatar
Parish, Chad committed
189
190
191
192
193
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
Parish, Chad's avatar
'minor'    
Parish, Chad committed
194
195
196
197
198
      "Ion string: \"He\"\n",
      "Energy: 30000.0 keV\n",
      "Target: \" Layer 1\"\n",
      "Displacement: 23.0 eV\n",
      "Number of ions: 200000\n"
Parish, Chad's avatar
Parish, Chad committed
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
     ]
    }
   ],
   "source": [
    "kp_flag = False  # Check for Kinchin-Pease\n",
    "\n",
    "# loop over the TDATA.txt lines\n",
    "for i,l in enumerate(file_lines_tdata):\n",
    "    # find the projectile\n",
    "    if 'Ion = ' in l:\n",
    "        st = l.find('Ion = ')\n",
    "        en = l.find('(')\n",
    "        ion_string = l[st+6:en-1]\n",
    "        print(f'Ion string: \"{ion_string}\"')\n",
    "\n",
    "    #find the displacement energy\n",
    "    if 'Displacement = ' in l:\n",
    "        st = l.find('Displacement = ')\n",
    "        en = l.find(',')\n",
    "        displacement = float(l[st+14:en-2])\n",
    "        print(f'Displacement: {displacement} eV')\n",
    "    \n",
    "    # Check for Kincihn-Pease model\n",
    "    if 'Kinchin-Pease Estimates' in l:\n",
    "        kp_flag = True\n",
    "        \n",
    "    # Find the ion (projectile) energy\n",
    "    if 'Energy  = ' in l:\n",
    "        st = l.find('Energy = ')\n",
    "        en = l.find('keV')  # if keV not found, find() returns -1:\n",
    "        if en > -1:\n",
    "            energy = float(l[st+10:en-1])\n",
    "        else:\n",
    "            energy = float(input('Energy not parsable -- input ion energy in keV:'))\n",
    "        print(f'Energy: {energy} keV')\n",
    "    \n",
    "    # Find the identity of layer 1\n",
    "    if 'TARGET MATERIAL' in l:\n",
    "        temp_l = file_lines_tdata[i+1]\n",
    "        st = temp_l.find('Layer # 1 - ')\n",
    "        en = temp_l.find('\\n')\n",
    "        target = temp_l[st+11:en]\n",
    "        print(f'Target: \"{target}\"')\n",
    "    \n",
    "    # Check for multilayer targets\n",
    "    if 'Layer #2' in l:\n",
    "        warn('Currently this code only parses single layers properly')\n",
    "    \n",
    "    # find # ions calculated\n",
    "    if 'Total Ions calculated =' in l:\n",
    "        st = l.find('Total Ions calculated =')\n",
    "        en = l.find('\\n')\n",
    "        num_ions = int(l[st+23:en])\n",
    "        print(f'Number of ions: {num_ions}')\n",
    "\n",
    "# If Kinchin-Pease not found...\n",
    "if kp_flag is False:\n",
    "    warn('Warning: this code is intended to use Kinchin-Pease data.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parse out the important information from RANGE.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
Parish, Chad's avatar
'minor'    
Parish, Chad committed
269
270
271
272
273
274
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T22:44:49.013029Z",
     "start_time": "2020-05-12T22:44:49.000051Z"
    }
   },
Parish, Chad's avatar
Parish, Chad committed
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
   "outputs": [],
   "source": [
    "#Find the important lines:\n",
    "for i, l in enumerate(file_lines_range):\n",
    "    if 'DEPTH' in l:  # line that defines the data columns\n",
    "        i_start = i + 3  # find the first line of data\n",
    "\n",
    "# Calculate the length of the data. SRIM is always 100\n",
    "length = len(file_lines_range) - i_start\n",
    "data_range = np.zeros((length, 3), dtype='f4')  # hold the data here\n",
    "\n",
    "# read the data\n",
    "for i in range(i_start, len(file_lines_range),1):\n",
    "    data_range[i-i_start,:] = np.array(file_lines_range[i].split()).astype('f4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform the calculation\n",
    "The details here are from:<br>\n",
Parish, Chad's avatar
Parish, Chad committed
297
298
299
300
    "Stoller, R.E., Toloczko, M.B., Was, G.S., Certain, A.G., Dwaraknath, S. and Garner, F.A., 2013. On the use of SRIM for computing radiation damage exposure. <i>Nuclear instruments and methods in physics research section B: beam interactions with materials and atoms</i>, <b>310</b>, pp.75-80.\n",
    "\n",
    "<span style=\"color:firebrick\">Addendum 2021-07-23:</span>\n",
    "Please see the recent paper S. Agarwal, Y. Lin, C. Li, R. E. Stoller, and S. J. Zinkle, \"On the use of SRIM for calculating vacancy production: Quick calculation and full-cascade options,\" <i>Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms</i>, <b>503</b>, pp.11-29, prior to using results output from this code in order to be sure you understand the proper range of application and internal calculation options of SRIM."
Parish, Chad's avatar
Parish, Chad committed
301
302
303
304
305
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
Parish, Chad's avatar
'minor'    
Parish, Chad committed
306
307
308
309
310
311
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T22:44:57.342276Z",
     "start_time": "2020-05-12T22:44:49.018025Z"
    }
   },
Parish, Chad's avatar
Parish, Chad committed
312
313
314
315
316
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
Parish, Chad's avatar
'minor'    
Parish, Chad committed
317
      "Please enter the fluence (ions/cm^2): 1.3e18\n"
Parish, Chad's avatar
Parish, Chad committed
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
     ]
    }
   ],
   "source": [
    "fluence = float(input('Please enter the fluence (ions/cm^2): '))\n",
    "\n",
    "# Calculate the NRT damage in dpa\n",
    "depth_in_nm = data_phonon[:, 0] / 10  # always angstroms\n",
    "ion_contribution = data_phonon[:, 1]\n",
    "recoil_contribution = data_phonon[:, 2]\n",
    "P = ion_contribution + recoil_contribution  # (P)honons\n",
    "T_damage = P * fluence / density * 1e8  # 1e8 is angstrom --> cm\n",
    "NRT = T_damage * 0.8 / 2.0 / displacement  # NRT displacements per atom. See citation.\n",
    "\n",
    "# Calculate the ion implantation in atomic fraction\n",
    "atom_frac = data_range[:, 1] * fluence / density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
Parish, Chad's avatar
'minor'    
Parish, Chad committed
339
340
341
342
343
344
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T22:44:57.560670Z",
     "start_time": "2020-05-12T22:44:57.343219Z"
    }
   },
Parish, Chad's avatar
Parish, Chad committed
345
346
347
   "outputs": [
    {
     "data": {
Parish, Chad's avatar
'minor'    
Parish, Chad committed
348
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsYAAAGaCAYAAAAMzIWkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde3yU5Z3///ecZ5KZSUISzocQkIMSMCioSFXUKm3drWV72tbawx5a/fZg3e6q60qr669d99fWVXuwdl2/7W7r2q7ttrZWW1EUS1WUg4AICiEchYSc5pA5398/JjMJmkBC5r5nwryej0cegeSe+/7MXAHeXPO5r8tmGIYhAAAAoMzZi10AAAAAUAoIxgAAAIAIxgAAAIAkgjEAAAAgiWAMAAAASJKcxS4gp60tZOr5a2oq1NkZNfUaKA2MdflgrMsHY10+GOvyUayxrq8PDPm9spkxdjodxS4BFmGsywdjXT4Y6/LBWJePUhzrsgnGAAAAwIkQjAEAAAARjAEAAABJBGMAAABAEsEYAAAAkEQwBgAAACQRjAEAAABJBGMAAABAEsEYAAAAkEQwBgAAACQRjAEAAEpO585tevSiuTq2dWOxSykrBGMAAIASc/iPzyjWfkRtm18sdillhWAMAABQYsL7WyRJyXCoyJWUF4IxAABAiQnt2yOJYGw1gjEAAECJYca4OAjGAAAAJSSdSChyaL8kKRkhGFuJYAwAAFBCIof2y8hkJDFjbDWCMQAAQAnJtVFIBGOrEYwBAABKSO7GO0lKRnqKWEn5IRgDAACUEGaMi4dgDAAAUEJyM8beugkEY4sRjAEAAEpIaN9euYNV8k+ZrlQkLMMwil1S2SAYAwAAlAgjk1H4wF75p82UKxBUJpVUJhEvdlllw1nsAgAAAJAVPXJImURcgemNMtJpSdk+Y4fHW+TKygMzxgAAACUid+Odf9pMufwBSdyAZyVmjAEAAEpEaF82GAemz1Q6FpXEkm1WYsYYAACgROSCsX/6TLkqmTG2GsEYAACgRORaKQK0UhQFwRgAAKBEhPbtkcPrk69+olz+oCSCsZUIxgAAACXAMAyF9++Vf2qDbHa7nLkZ4wjB2CoEYwAAgBIQ7zymZLhHgekzJYke4yIgGAMAAJSAgUu1SaLHuAgIxgAAACUg1LpHkhSY0ShpQDCmlcIyBGMAAIASEBqwIoU0sJWCdYytQjAGAAAoAe9spWBVCqsRjAEAAEpAaF+LbE6nKidPkyS5Kv2SaKWwEsEYAACgBIT3t6hy0jTZnU5Jkt3plMNXoWQkXOTKygfBGAAAoMiSkZBix9ryN97luCoDtFJYiGAMAABQZKF9x/cX57j8BGMrEYwBAACKLJxfkaLhuK+7Kv30GFuIYAwAAFBkuRnjwPS3tVL4A0r3RpVJpYpRVtkhGAMAABTZ25dqy8mvZcwNeJYgGAMAABRZaF921zv/21sp/GzyYSWCMQAAQJGF9rXIN2GynF7fcV/Pb/JBn7ElCMYAAABFlE7EFX3roALTZ77je/3bQhOMrUAwBgAAKKLwgVbJMBSYNkgw9hOMrUQwBgAAKKKhbryTBgRjWiksYUowTiQS+tKXvqTm5mZdeeWVeuWVV8y4DAAAwJiXu/Hu7Uu1Sf2tFClmjC1hSjB+7rnn9NRTT+mRRx7Rueeeq7vuusuMywAAAIx54dyud4P1GLMqhaWcZpx09uzZqqio0Pjx4xUMBuVyuU76mJqaCjmdDjPKyauvD5h6fpQOxrp8MNblg7EuH+U21vG39kmSGpqb5K06/rnHpk6QJLmMxGn5upTaczIlGAeDQU2fPl3Lly9XJpPRj370o5M+prMzakYpefX1AbW18TZEOWCsywdjXT4Y6/JRjmN9bPeb8lSPUyjhUOhtzz2ayk4adh9tP+1el2KN9YnCuCmtFA8++KBaW1v105/+VH/5l3+pL33pS4rH42ZcCgAAYMzKpNOKHGiVf5D+YklyVebWMWbnOyuYEowrKyvldrvl8Xjk9/vV3d2tFHt8AwAAHCf61gFlUslBl2qTWK7Naqa0UnzqU5/S66+/ro9+9KMKBoP653/+Z1VWVppxKQAAgDErvG+vJMk/vWHQ77Ncm7VMCcYVFRW69957zTg1AADAaeNES7VJkt3tkd3pYsbYImzwAQAAUCSRwwckSZWTpw36fZvNJmeln2BsEYIxAABAkSS6OyVJnuraIY9x+QMEY4sQjAEAAIok0dMlSfJU1wx5jMsfUDLCBh9WIBgDAAAUSbxvxtgdrB7yGFdlUMlwSIZhWFVW2SIYAwAAFEmiu1MOr08Oj3fIY1z+gGQYSkUjFlZWngjGAAAARZLo6ZK7aug2Cokl26xEMAYAACiSRHenPFVDt1FIkquSTT6sQjAGAAAogkw6rURP9wn7iyV2v7MSwRgAAKAIkqFuSRp+KwXB2HQEYwAAgCJIdPct1TbsHmOWbDMbwRgAAKAI8ku1nSwY02NsGYIxAABAESSGsYax1D9jnIqETa+p3BGMAQAAimDYM8b+oCRmjK1AMAYAACiC/HbQw22lYB1j0xGMAQAAiiAx7BljeoytQjAGAAAogv5gzDrGpYJgDAAAUAS5HuOT7nzHcm2WIRgDAAAUQa7H+GStFE5fpSRmjK1AMAYAACiCRHeXbA5HftWJodjsdjkrAwRjCxCMAQAAiiDe1Sl3oEo2m+2kx7r8AValsADBGAAAoAgSPV0nbaPIcTFjbAmCMQAAgMUMw1Ciu3P4wdhPMLYCwRgAAMBi6d6oMsnESbeDznH5A8okE0on4iZXVt4IxgAAABbrX6pt+DPGEitTmI1gDAAAYLH+pdqGOWNcSTC2AsEYAADAYolTnjFmkw8zEYwBAAAsFh/mdtA5+WAcDZtWEwjGAAAAlkt0D2/XuxxXZXYTEFopzEUwBgAAsNiIgzE331mCYAwAAGAxVqUoTQRjAAAAi+Vuvhv2Osa5VSnYFtpUBGMAAACL9S/XxoxxKSEYAwAAWGzEM8Z9wTjFcm2mIhgDAABYLN7dKaevUg63e1jHuyr9kmilMBvBGAAAwGKJ7s5ht1FIksvPcm1WIBgDAABYLNHTLc8wN/eQJGdF34wxwdhUBGMAAAALZZJJJcM9clcPf8bY4XbL4fHSSmEygjEAAICFEqFuScO/8S7H5Q8wY2wygjEAAICF8itSjKDHWMquZUwwNhfBGAAAwEK5NYyHu+tdDjPG5iMYAwAAWCjedYozxv6AUr0RZdJpM8qCCMYAAACWGunmHjm5baFT3IBnGoIxAACAheJ9wXikrRTO3LbQkXDBa0IWwRgAAMBCuR5j9wjWMZbY5MMKBGMAAAALjbaVgrWMzWNqMP7BD36gZcuWafny5Xr88cfNvBQAAMCYED/V5dpyrRTMGJvGadaJ161bp/vvv18PP/ywtmzZooMHD5p1KQAAgDEj0ZPd4MMzgp3vJIKxFUwLxuvXr1dVVZW+9rWvKRQK6atf/eoJj6+pqZDT6TCrHElSfX3A1POjdDDW5YOxLh+Mdfk43cfaiPbI7nRq0oxJstlsw35c+6R6SZLXljhtXqNSex6mBeOenh61tbXp61//un7zm9/oxhtv1Lp164b8AejsjJpViqTsC9/Wxv+wygFjXT4Y6/LBWJePchjrSHu7XMFqtbePbHWJ3oxLktTxVttp8RoVa6xPFMZN6zH2+/2qq6vTsmXLdNlll6mtrU2dnZ1mXQ4AAGBMiHd3jnipNklyVfol0UphJtOC8ZIlS3T06FG99NJLevnll1VXV6eampH/EAAAAJwuDMNQoqdrxDfeSSzXZgXTWikuv/xy/c3f/I2+8IUvKBgM6pvf/OaI+mgAAABON6loWEYqdYrBmOXazGZaMJakG2+8UTfeeKOZlwAAABgz4l25NYyrRvzY/DrGzBibhg0+AAAALJI4xe2gJZZrswLBGAAAwCL920GPPBg7vD7ZHA5aKUxEMAYAALBIojsbjE9lxthms8lVGWDG2EQEYwAAAIvkt4MOVp/S413+gJLhnkKWhAEIxgAAABYZTSuFlA3GqcjINgbB8BGMAQAALJK/+a76VINxUMlISIZhFLIs9CEYAwAAWGTUrRSVARnptNKx3kKWhT4EYwAAAIvkZoxH00ohsWSbWQjGAAAAFsn3GI/i5juJ3e/MQjAGAACwSKKrUy5/UHbnqW0+zO535iIYAwAAWCTe3XVK20Hn9LdSsGSbGQjGAAAAFkl0d55yf7FEj7HZCMYAAAAWSCcSSvVGTmnXuxxnpV8SPcZmIRgDAABYYLSbe0iSqzIoiRljsxCMAQAALJAY5RrGEq0UZiMYAwAAWCA+yjWMJW6+MxvBGAAAwAK5VorR9Bjnl2uLhAtSE45HMAYAALBA/653p95K4Q7QY2wmgjEAAIAFChGM+1spugtSE45HMAYAALBAonv0rRTOCr9kszFjbBKCMQAAgAUKcfOdzW6XqzJAMDYJwRgAAMACuRnj0SzXJmXbKViVwhwEYwAAAAvkZoxH00ohiRljExGMAQAALJDo6ZLd5ZbDVzGq87gCQSXDIRmGUaDKkEMwBgAAsECiu1PuqhrZbLZRncflDyiTSiodjxWoMuQQjAEAACyQ6O4cdX+xNGCTD9opCo5gDAAAYDIjk1Gip0ue6tH1F0vZVgqJYGwGgjEAAIDJkpGQjExmVJt75PTPGLMyRaERjAEAAEyW3/WuEK0UfoKxWQjGAAAAJosXYNe7nP5gTCtFoRGMAQAATJbf3KMAwdjt7+sxjhCMC41gDAAAYLJEAbaDznHlgnGIVopCIxgDAACYLG5KjzEzxoVGMAYAADBZokDbQUsDZoxppSg4gjEAAIDJEj2F6zHOzRgnaKUoOIIxAACAyQo7Y9zXSsGMccERjAEAAEwWPXJIkuStqx/1ufKtFPQYFxzBGAAAwGSh1t3y1k3I71o3Gs6KSslmY1UKExCMAQAATJROJBQ5tF+BGY0FOZ/NZpPLH1AyQjAuNIIxAACAicIH9srIZBScMatg53T5g7RSmIBgDAAAYKJQ625JUqDQwZhWioIjGAMAAJioPxgXppVCUl8rRUiGYRTsnCAYAwAAmKpnb+FnjN3+oIx0WulYb8HOCYIxAACAqfIzxtMLO2MsSckw7RSFRDAGAAAwUah1j3wTJsvpqyjYOfuDMTfgFZKpwfjVV1/VggUL9Itf/MLMywAAAJSkVKxX0bcOKFjA/mKJTT7MYlow7u7u1pe//GUlk0mzLgEAAFDSwvtbJBW2v1jqnzFOsDJFQZkWjG+++WZdeeWVZp0eAACg5PX3F5sTjNnko7CcZpz0P/7jP3Ts2DHde++9evDBB4f1mJqaCjmdDjPKyauvH/02jBgbGOvywViXD8a6fJxOY93aflCSNGXhWQV9Xkcn1kuSvLbkmH69Sq12U4Lxf/3Xf+nYsWO64IILJEm33367pk+frnPPPXfIx3R2Rs0oJa++PqC2NvpwygFjXT4Y6/LBWJeP022sD7/2miTJqJ5c0OcVk1uS1HH46Jh9vYo11icK46YE45/+9KdKpVKSpMsuu0xf/OIX1dTUZMalAAAASlaodbdksykwfWZBz+uq5OY7M5gSjCdOnHjc72tqauTxeMy4FAAAQMkKte5RxcQpcni8BT0vy7WZw5RgPNDOnTvNvgQAAEDJSUUj6j16WBPPv7jg53YHmDE2Axt8AAAAmCC0b4+kwi/VJrHznVkIxgAAACYItfYF4wYTgnElrRRmIBgDAACYoCe3hrEJM8YOX4VsDoeSEYJxIRGMAQAATJDb3CNoQjC22WxyVQbY+a7ACMYAAAAmCLXukc1uV+WUGaac3+UP0GNcYARjAAAAE4Rad6tyynQ53G5Tzu8KBJWilaKgCMYAAAAFlgz3KHbsqCn9xTmuyqCS4ZAMwzDtGuWGYAwAAFBgPbkVKaY3mnYNlz8gI5NRKhox7RrlhmAMAABQYCETV6TIceU2+aCdomAIxgAAAAVmSTDOrWXMyhQFQzAGAAAoMDOXasvp3/2OGeNCIRgDAAAUWKh1j2xOpyqnTDftGm4/rRSFRjAGAAAosFDrbvmnzJDd6TTtGvkZY1opCoZgDAAAUEDx7i7FuzpM7S+WBrZSEIwLhWAMAABQQP033pm3VJs0IBjTSlEwBGMAAIACsmJFCkly9fUYJ2ilKBiCMQAAQAFZsSKF1B+MWZWicAjGAAAABRTK7XpHK8WYQzAGAAAooNC+3bK73KqYNM3U6zBjXHgEYwAAgAIxDCO7VNu0mbI7HKZei1UpCo9gDAAAUCDxrg4leroVNLmNQpIcHq9sTiczxgVEMAYAACgQq1akkCSbzSa3P0gwLiCCMQAAQIFYGYylbDsFrRSFQzAGAAAoEKs298hxVQaYMS4ggjEAAECB9C/VZtGMcSCoZCQkI5Ox5HqnO4IxAABAARiGoY7Xtsjhq1DFhMmWXNPlD0iGoVQ0bMn1TncEYwAAgALo2L5ZodbdmnLxlbLZrYlYrkrWMi4kgjEAAEABtPzqYUnSzD//iGXXdAUIxoVEMAYAABilTDKpvb/7hTzj6jT5wsssu25uk48EK1MUBMEYAABglA79cY3iHe1qeM8q2V0uy67rqsztfseMcSEQjAEAAEap5dePSJJm/vlHLb0urRSFRTAGAAAYhUSoWweeflzBmWdo3IJmS6+da6Vgk4/CIBgDAACMwr7f/0qZRFwz/+wjstlsll67v5WCYFwIBGMAAIBRyLVRNPzZhy2/dv+MMa0UhUAwBgAAOEXhg/t0dMMfNX7JhfJPmW759d35HmNmjAuBYAwAAHCK9v7m55KsXbt4IDb4KCyCMQAAwCkwDEMtv/5vOTxeTb/i/UWpgVaKwiIYAwAAnIKObZvU0/KGpqx4j9yBqqLUkFuujQ0+CoNgDAAAcApafv3fkorXRiFJDrdHdpebGeMCIRgDAACMULG2gB6Myx8gGBcIwRgAAGCEirUF9GBc/qCSEYJxIRCMAQAARqhYW0APxuUPKBmix7gQCMYAAAAj0LP3TR146jcKNs61fAvowbj8AaWiYWXS6WKXMuYRjAEAAIbJMAy9/PWblEkltehLt1q+BfRg3P7syhSpaLjIlYx9BGMAAIBhOvD04zr8/BpNXLZC0y7/s2KXI0ly5tYypp1i1AjGAAAAw5CK9eqVf7lFNqdT5/7jXSUxWyxlb76TxA14BWBKMA6Hw7ruuuvU3Nysyy+/XGvXrjXjMgAAAJZ57cF7FDm4T/OuvV5VjXOKXU5erpWCJdtGz5Rg/OMf/1hbt27VL3/5Sy1ZskQ333yzGZcBAACwRPjAXr327/8mX/1ENV3398Uu5zi5baETtFKMmtOMk37iE5/Q+9//fk2ZMkW1tbVKD+MuyZqaCjmdDjPKyauvD5h6fpQOxrp8MNblg7EuH6U41i985atKx2O6+I5/1eSGycUu5zg1E+slST57siRfuxMptXpNCcaBQECBQEBPPPGEHnroIX3uc5876WM6O6NmlJJXXx9QWxtvMZQDxrp8MNblg7EuH6U41oeeX6Pdj/9K9Ysv0LiL/qzk6ovLLUnqONxWcrWdSLHG+kRh3LSb7x577DHdeOONWrlypa6//nqzLgMAAGCadCKhl79+k2x2u5b807+WzA13A9FKUTimBOPNmzfr5ptv1ooVK3TbbbcpEomYcRkAAABTvf6f31do75s646N/pZp5TcUuZ1C5YJyMEIxHy5Rg/MADDyiVSumpp57SeeedpyVLligej5txKQAAAFO89cKz2vq9u+SpqdXCL9xa7HKG5KrsC8asSjFqpvQYf+973zPjtAAAAJZo+fUjeuGfPi9JOu9ffyhPVXWRKxqaK1AliQ0+CsGUYAwAADAWGYah7T+8W1v+7Q65AkFdfN9PNGHpu4pd1gn1t1IwYzxaBGMAAABJmVRKG+78e735s4dUMXGqVvzg56o+Y36xyzopWikKh2AMAADKXioa0fNf+YwOrn1SNXMX6JIf/FwV4ycVu6xhcbjdcni8SoZppRgtgjEAAChrsY52rb3uwzq2daMmLluhi/7tR3L1bbM8Vrj8AWaMC4BgDAAAylb06GGt+czV6tmzU41Xf0zn3X6P7C5XscsaMVclwbgQTNvgAwAAoJSFD7bqD594j3r27NS8T16v8/+/747JUCxJrkCQYFwAzBgDAICy09Pyhtb81dWKvnVQTdf9g5o+f0tJ7mo3XC5/QKneiDKplOxO4t2pYsYYAACUlc6d2/SHa9+r6FsH1fx3d2jhF/5xTIdiacDKFJFwkSsZ2wjGAACgbLS/+oqe+uRVih1r05Lbvqkz/+qLxS6pIFyB7M2CrEwxOsy1AwCA0174YKv2/O/D2vF/v6t0b0QXfP17arz6Y8Uuq2D61zImGI8GwRgAAJyWUr1R7f/DY9r9y5/oyIvPSZKcFX4t/9ZDmn7l+4tcXWHllpfjBrzRIRgDAIDTSvfunXr9x99X6+9+kZ9BrV98gWat+rimX/n+/Ozq6aS/lYJgPBoEYwAAcFqIHjmkV7/7L9rzi/+SkcnIN2Gy5nzsb9R49ccUbJhV7PJMRStFYRCMAQDAmJbo6dL2f79HO//z+0rHYwo2ztWiL92qqZe+T3aHo9jlWcLlzwVjZoxHg2AMAADGpHQ8pl0P/7u2/eBbSnR3yjd+khZ+4RY1vv9jZbeWr5tWioIor58aAAAw5nW+vlW7H/1Ptfzm50p0d8oVCOrsG7+muR//Wzl9FcUuryhyrRQJWilGhWAMAABKXqKnS3t/+z/a/Yv/Usf2zZIkb229zvqbGzX/01+Qp7qmyBUWV66VIkUwHhWCMQAAKCmZZFI9e99Q167X1Llzu7p2bdeRF59TOh6TzeHQlBUrNWvVJzTloitkd7mKXW5JyPcYR2ilGA2CMQAAsFT0yCGF9+9VrKNNsY5jine0KdbRrlj7UUUPtqhj5w5lUsnjHhOceYYaP/BxNb7/o/LVTyxS5aXL5a+SJMW7O4tcydhGMAYAAJYIH2zVq9/5hlp+/YhkGIMe46yoUM38JlXPOUvVc87s+3yWvDW1Flc7trirquWtrVfHa6/KMAzZbLZilzQmEYwBAICpetuPavsPvqk3HnlImVRS1XPO1JSLV8pbWyfPuHp5x9XLW1sn77g6TZ3bqGMd0WKXPObYbDbVN5+n/U/9RpFD++WfMr3YJY1JBGMAAGCKRKhbOx66T6//6PtK9UbknzpDC79wqxre90HZ7PZBH1Mu6w6bIReM2za9SDA+RQRjAABwSnrbj6pr12vqadmlWEeb4p0dine0K951TLGOY4oc2q9UNCxv7Xg1f+V2zfqLa+Vwu4td9mmrfvEFkqS2TS9q5lUfKnI1YxPBGAAAHCeTTCoZCSkZ7lEyHFIyElYyHFKs/Yi63ngtu1rEru2Kd7QPeQ53sFoVEyZp5p9/VPM+8Tk5KyotfAblqWb+Qjk8XrVvfKHYpYxZBGMAAMpcOpFQ2yt/0sHnntTBZ3+v0N43T/oY/9QZql+0RFVzzlL17Hny1U+Qp6ZWnpo6eapqWEatCBxut2qbFuvoK39SItQtd6Cq2CWNOQRjAABOY8lwj+JdHTLSGWXSKRnptIxMWkY6pY4dr+rQs7/X4fVrlYqGJUlOX6UmLF0ud7BaLn9Azgq/XP6gXP6APNXjVDV7nqrPmJ/faQ2lpb75fB19eb3at7ysycsvK3Y5Yw7BGACA04hhGOratV0Hn/29Dj33e7VvfklGJnPCxwSmN2ryxddoysVXavy5y+RweyyqFoVWv/h8SVLbphcIxqeAYAwAQJFlUiklujsV7+pQvPOY4l0dSnR3Zvt8I+F8j28qElYqFpXTVyl3sEouf1DuQJVcgaAcHo/aNr6oQ+v+oOhbByVJNrtdtYuWKDhjlmwOh2wOp2wOh+x2u2wOhyonT9Pki65QsGF2kV8BFErd2UslSW0bXyxyJWMTwRgAgAJJJxIyMmk5PN5BN1hIRsLq2rVdna9v7fvYplDrbiV6ugpWg7uqRg3v+5AmX/RuTX7X5fJUjyvYuVH6PFXVqpo9X+2vvqxMMkmv9wgRjAEAkGRkMooeOaRQ626FWveoZ++bihxslc3ukLOisu/DL1dFpZwVFUqEetR79LCiRw6r98ghRdveyq/SYLPb5azwy1nZd3ylX8meboX2txy345vd6ZJ/eqOq556VvXGtalzfDWzj5Kmq6evxDchV6ZerMns+p69CqWhEiXCPkj3d+c/JaFjVc85S3cJzZXfyz3s5q28+T91v7lDnzm2qXdBc7HLGFP7kAABOC6loRJFD+xU+tF+x9iNKRsJKRSN9n7O/TkUjyiQTSieTyiQT+Y9UNKrwwValY72ndG2nr1K+CZNUfcZ82V2e7PUiYSWjYSV6uhQ5fEAOt0cTli5XzdwFqpnXpJp5TQo2zj2ldX091ePE4mcYSv055+vNn/9ftW16gWA8QgRjAIBlUrFe9ezZ2bfG7ZQh3+ZNxXrVtXObOl7boo4dr8oejyiezMhmt8tmd8hmt8lmtysZDilyaL8ih/Yr3tUx8oJsNjncHjk8HgUbZivQMFuBGbMUnNGowIxZ8k+bKZvNpmQ0cly4TvVG5KwIqGLCJFVMmCSXPzjKVwYonPrmvhvwNr6oeZ+4rsjVjC0EYwAoQ9kbukKDf9MwlEmlZGQyMvqW98qkU0pFwtm2gaOHFH3rkKJHDyv61iGlomH5pzUo2HCGgjNnKzjzDAVmzJY7WKVQ6261b9mg9ldf0bFXX1bnru0yUilJ2XYD34TJ8k+Zrsop01Uxcap6jxzUse1b1LNnp4x0eljPxeH1qXLyNI1b0KzKSdNUOXmafOMn9rUfBI5rg3D6KuRwu2V3uWV3e4a9/bC3tn5YxwGlwD91hrx1E9S28QUZhjFovzsGRzAGgBJgGIYiB1t1bNtmdWzfpNixNo1b0Kzxi89X1RlnDjvADSYV61Xn61t1bOtGdWzfpGPbNqmn5Y3jel1Plc1ul93jVefrW9/xPbvbo0wi3v97l1vjzjxb485cpGQkpMjB/Yoc2qejr/xJenl9/jiHr0J1i5Zo3PyFGnfW2ao5c5Gmzp2l9rYeycj0rcObkWFk5F8hii4AACAASURBVPRWyDOujn/4gQFsNpvqF5+n/b//tSIHW+Wf2lDsksYMgjEADGBkMooc2qfu3TsV2tcim8Mhp8cru8crp9crh8cnh9crd1WNfLXj5a4e947QahiGEt1dCh/Yq/CBVoUP7FUyEpbd4ZDN6ZTd4ZTN4ZTd6VCs41g+rCa6O487z57//akkyeUPqu7spapffJ7GnXm2ZGSUikaV6o0o1Rvt+4hke1rDoexHNPvrRE+XQq27j5t9dVYGNGHJhfLWTxzydbAPWNrLZnfI7nTK4fWpYsLkvvaBKaqYOFne2vGyORyKtR9RT8sb6ml5Uz1731DP3t3qPXpYVbPnqW7huapbeK6q5y4YtJ82nUgo+tZBRQ7vl69uggINs9/xmlbUBeQzvCMeT6BcjV98vvb//tdq2/giwXgECMZAEUQO7dfBZ5+UJE2+6Ar5p0wvckWnJp2IK97VIYfbK4fHk12iym7Pf98wDKVjvUqGe5QI9fSFth4lQz1KhLuVDOV+3aNUNCLf+Il9b8dn35IfuLOWYRiKtR/tWzFgt0L7WuR125RyVshdVSN3VbU8wWq5q2pkd7uVjsWUjvfmP6diMWWS2aW0lDGUSafzs4/JaEQ9LbvU/ebr6t6zS+ne6PBfBJtNnppaecfVyzuuVolQTzYIh3pG9Fr6p83UpGUrNO6sZtU2NctbU6f2V19W28YXdHTjCzr8/FM6/PxTIzqnzemUqzKg2qZzVLugWeMWNKt2wWIFG2YfN06F4KufKF/9RE1Y+q4RP9bhdiswfaYC02cWtCagnOX7jDe9qJl//pEiVzN2EIxx2kuGe3Rs6yYlIz0KzJilwPRGOTzWzjwZhqHOHa/qwNOP68DTj7/jbeeaeU2aeul7NfXS96pm/sL828KGYSje0Z6ddTzYqt4jh5QIdefDZDLUo2S4R+l4XBWTpiowvbEvYDQqMKNR3roJSkVCCvfdnJT7iB4+IJvdIVcg2LdBQDD/a1dFpRwer+wez3GBN9Hdpe49O9WzZ1ff5zcU3t/yjh217E6XHF5f9saoaDjfT3oqfPUTFZjRqGQkrNC+FqWG6oktELvLrWDjGaqaPV9Vs+YqOCO76UE6EcsG7L7PqVhU8a5OxY4dVexYm+Id7eo9ekjdb+6Qw+OVf2qD/OdckP08rUH+qQ1yB6pkZNLKpLJb8WbSKRmplFyVAdXMXyhPdc076qmaPU+zVl0jSYoda1PbphfV9cZrsrvc2Z5ZX4Wcvr7PFRV9PbXZrXudlf4h19IFcPqrmdckh69CRze+UOxSxhSbYRSgyawA2trM/Qevvj5g+jVKUao3qq43dgxYTH6revbslMNbIe+4Onlr6+UZV9c321Un+2DLBhmG0olE393Y4fwSSKloRIZhZEOVPxus3H2fnRWV2beL828bZ99CttnsyqSSyhy3VFJSmVRywE0+6b4ewuzn3G5N2fM58ud1+ir6Q50/KFegSu5AUJ5Et95Y+6yOvfqy2re8rO49O4/rpbTZ7aqcMj07M9l4hirGT1Iq1jvgefU9x95oNsDkbkJKpfKBJpNMZWsf+FxSKTk8nmw4qfDLWZlde9Th8+nY1o2KHj4gKRscJ5x/kaZe+l5JNh14+nEdefE5ZZIJSVLFxKmqnnOmIocPKHKgVaneyEnH2eZ0DhpA7U6XMqnkqf3wnISnepyCjXPlGz8xu/xVPK5MIqZ0PK50PKZMOi2XPyBXZaB/jPyBvo/g8T83gaCcvkpFDh9Uz943FGp5Qz1731RPyxvZZa483mzgnzGr76NRgemNqptYqyOth7I7hnV3KtHTpUR3pzLJpBwerxx9rQ9Or1d2t1cOtzvfGtC/uoFddo9HwYYz5J86Y1Trv6YTCdldLsKoCcr17/ByxFgXzlOfukpHNvxRH/pTi9zB6mKX8w7FGuv6+sCQ3yvrYHzkpXXa+/gvskEm9w92Zfazw+PNzgz1RpSK9SoVjSodiyoV61UmET9+DcxENiA5vL78HdD5xdgr/DIyGaXjsQEf8exbu73R49fYjISVjEZkpNNyeH35mbrsR3b2zuZ0Zv8hHxA6ZbMpHes7X6w3W2dvrxI9XQrv23PcjJ7N4VBgxixlkgnFjrUrFQ2b+roXm9NXqdqmZtUuPFee6lqFWt9U955dCu19U7FjbcM+Tz7g2x3ZcO5yy+50ye5yZX/tcsvucCidiGfbBfrGNBfIXYGgplx0paZe+l5Nftdl71jaKRnu0aF1a3Tg6cd18LknlQz1yFkZUGDqDFVOnSH/1BnyT21QxcQpcldV5/8D4vJXyeUPyGazqbftLYX27cl+tLYovG+PIof3y101TpWTp2Xv/J+cvWO/ctJUGYbR3+LQN/OcCHUr3RvN/owmYkonEsrEY0on4nL6KhRsnKuqxjkKNs6Rd1xdQcdqKOlEXHana9C3/vkHtHww1uWDsS6cLffeqW33f1OX3P9zTbno3cUu5x1KMRiXdSvF3t/+j978+Y+KXUae05ddUsjmcCjR05UP0qdy57jd5Zar0q/6xef3LyY/v0lVs+Yd10aQ6o0q1tGueEe7Yh3tQ84uOtyefOB3Vvizv67ILi8/8C39XMhK9UbySzwZ6b63jlPZWeBcmHS43NnZNacrGzKdfTf72B39s3oOe18/aKp/9rbvcyoaUTLU/Y7rV02cIP+cRapdeI6qZs8f8m7+RE+XelreVG/7ETl9uedW2fc8A3L6fLL1BbJTmQE0Mpn8f3481eNOuC2nyx/UjPd8QDPe8wFlkkklI2G5q6pHdN3sTVGTNWHJ8hHXWsocbk+xSwCAMWlgn3EpBuNSVNbBeMnqb2veJ/9P/13ckZ78r9OJuJxen5xenxy+imwPn7dCDp+vfw3Mvo/s27POvhAUHrDbUkipaEQ2u6Nv9nfgLLCnP2BW+uX0VQ46I2YYRt/b1NmZZiOTlpFK9bcbpLNv8zt9Pjm8FdnPHt+w3w52+irknzJ9VDd/OSsqpfGTTvnxhTbc/4G6g9WqW3SuaXXY7Pb8OwcjYXe5Bu03BQBgJOrOXiLZbGrfRJ/xcJV1MLY7HKpqnFPsMk7IltuVye2Rhp75BwAAOI47UKXqM+ar/dVXlEkmT/jOJbIKu14PAAAASkb94guUjvWqY8erxS5lTCAYAwAAnKbqm8+TJLXRTjEsBGMAAIDTVP3ibDA++Mzv1Nt2pMjVlD7TgnEmk9Hq1at1zjnnaNWqVWppaTHrUgAAABhE5eTpCkxv1JGXntcvLpmnP3zyfdr50x8Skodg2jrGTz75pL7yla/o4Ycf1t133y2Xy6X7779/yOPZ4AOFwliXD8a6fDDW5YOxLrze9qNq/d2j2vfkr9SW2wnPZtP4cy7Q+CUX9q0V37dUqs0um8OeXSnL5M2Kmq5epWRggqnXGExRNvj4xje+ofXr1+uxxx7Tgw8+qAceeEAvvvjikMenUmk5nYOvNwsAAIDRCx8+qDce+4Xe+NXPdeilPxa1lvkfvkZXfrd09pOQTFyuLRQKyevNbiTh9XoVCp34f3+dnVGzSpHE/0DLCWNdPhjr8sFYlw/G2mTOoKZ+4FOa+oFPKXrkkHpa3pCRMbL7JKTTMoyMjHRaMjInP9eo2DT/ineXz853fr9fsVhMkhSLxRQIsAgvAABAqcjtmFosvtqAwiX2nyDTbr5btGiR9u7dqx07dmj9+vVqbm4261IAAADAqJk2Y7xy5Uq99NJLuuaaazRjxgx961vfMutSAAAAwKiZFowdDoduv/123X777WZdAgAAACgYNvgAAAAARDAGAAAAJBGMAQAAAEkEYwAAAEASwRgAAACQRDAGAAAAJBGMAQAAAEmSzTAMo9hFAAAAAMXGjDEAAAAggjEAAAAgiWAMAAAASCIYAwAAAJIIxgAAAIAkgjEAAAAgiWAMAAAASCqDYJzJZLR69Wqdc845WrVqlVpaWopdEkbo4x//uObOnau5c+eqqalJLS0tWrVqlc455xytXr1amUxGknTfffdp6dKlWrlypTZv3ixJIzoWxfPjH/9Yc+fOVTweN218hzoW1ho41vv378//2Z47d65Wr14tSXrkkUd0wQUX6JJLLtHTTz8tSWpvb9e1116r5uZmffGLX1Rvb++Ij4U1wuGwrrvuOjU3N+vyyy/X2rVrtXnzZq1cuVJLly7VfffdJ2nof59HeyysMdg4/+lPfzruz/QDDzwgaYz9/W2c5p544gljwYIFxtatW43PfOYzxmc/+9lil4QRSKfTxtlnn2389re/Nbq7u42enh7js5/9rPGZz3zG2L59u3HWWWcZTzzxhLF161Zjzpw5xrp164xbbrnFuOqqqwzDMEZ0LKzX29tr/Mu//Isxb948Y86cOUYsFjNtfAc7FtYZbKwff/xxY+nSpUZ3d7fR3d1t9Pb2Gm1tbcaZZ55p/OxnPzPuueceY+nSpUYikTDuuOMO46qrrjLefPNNY9myZcaDDz44omNhne9+97vGhRdeaLS0tBg333yzcd555xlXXXWVccsttxjr1q0z5syZY7z66qtD/vs82mNhjcHG+Yc//KFx9dVX5/9Mx+PxMff392k/Y7xx40Y1NDRowYIFWrZsmTZt2lTskjACLS0tikajuueee/ShD31I69at08aNG3XhhRfqzDPPVENDgzZt2qSNGzeqoqJCy5cv10UXXaRdu3YpHA6P6FhYr62tTS0tLbr++uvzXzNrfAc7FtYZbKy3bdumRCKhVatW6Ytf/KI6Ojq0ZcsWpVIpXX755VqxYoW6urq0Z88ebdy4UUuXLtWsWbO0YMECbdq0aUTHwjqf+MQn9Mgjj6ihoUG1tbVKp9PatWuXLr74Yi1fvlwVFRXavHnzoP8+h8PhUR8Laww2ztu2bdOhQ4d09dVX69Zbb1UsFhtzf3+f9sE4FArJ6/VKkrxer0KhUJErwkgYhqEPf/jD+sY3vqH3vve9uummm9Td3f2OMQ2FQvL5fJKU/xwOhwcd/6GOhfWmTZum+++/X1OmTMl/bSRjNtpjYZ3BxnrmzJn69Kc/re985zvq7u7WnXfemR8Xn8+XH7/cuA421sM9FtYJBAKaMmWKnnjiCT300EN697vfLUmDjt9QYzqaY2GNt4/ztddeq7POOkt/+7d/q3vuuUdbt27VfffdN+b+/naaduYS4ff7FYvFJEmxWEyBQKDIFWEkZs+erZtuukl+v19VVVX63ve+J0mKx+OS+sd04Djn+gn9fr/8fv+wj0VpGMmYjfZYFNcHPvABJRIJ+Xw+XXjhhXriiSe0atUqSdkxyo1fblwHG+vhHgtrPfbYY7rpppv0nve8R7fccoseffTRE/5ZffuYjuZYWGfgOF9//fXKZDLKZDLyeDxqbm7Wrl27dOmll46pv79P+xnjRYsWae/evdqxY4fWr1+v5ubmYpeEEfjlL3+ppUuXavPmzVqzZo18Pp+WLVum9evXa/v27WptbVVzc7MWLVqkSCSi9evX69lnn9WcOXPk9/u1cOHCYR+L0jCSMRvtsSiuT33qU/rMZz6jgwcP6sUXX1RTU5OamprkcDj09NNP65lnnlF1dbUaGxu1cOFCbdiwQbt379bWrVvV3Nw8omNhnc2bN+vmm2/WihUrdNttt8kwDM2cOVPPPfec1q1bp2g0mv+z+vZ/nwOBwKiPhTXePs6RSERXXnmlbrrpJu3bt09btmxRU1PT2Pv727Tu5RKRSqWM1atXG4sXLzY+8IEPGHv27Cl2SRiBZDJp3HrrrUZzc7NxxRVXGGvXrjX27NljrFq1ymhubjZuu+02I51OG4ZhGPfee6+xZMkS44orrjA2bdpkGIYxomNRPI8++mj+hiyzxneoY2GtgWO9a9cuY9WqVcbZZ59t/PVf/7Vx7NgxwzAM47//+7+NZcuWGRdffLGxZs0awzAMo62tzfjkJz9pnH322cbnP/95IxqNjvhYWOO6664z5syZc9zHpk2bjJUrVxrnnnuuce+99xqGMfS/z6M9FtYYbJxfeuklY+XKlcbixYuNG2+80YhEIoZhjK2/v22GYRjmxW4AAABgbDjtWykAAACA4SAYAwAAACIYAwAAAJIIxgAAAIAkgjEAAAAgiWAMAAAASCIYAwAAAJIIxgAAAIAkgjEAAAAgiWAMAAAASCIYAwAAAJIIxgAAAIAkgjEAAAAgiWAMAAAASCIYAwAAAJIIxgAAAIAkgjEAAAAgiWAMAAAASCIYAwAAAJIkZ7ELyGlrC5l6/pqaCnV2Rk29BkoDY10+GOvywViXD8a6fBRrrOvrA0N+r2xmjJ1OR7FLgEUY6/LBWJcPxrp8MNb90ol4sUswVSmOddkEYwAAgLHi6Ct/0iOLJ+ngs08Wu5SyQjAGAAAoMQee+Z2MTEZ7fvnTYpdSVgjGAAAAJaZ9ywZJ0qHn15z2LRWlhGAMAABQQjLJpDq2b5YkpaJhHXlxXZErKh8EYwAAgBLSuXOb0rFeVc2eLynbVgFrEIwBAABKSK6NYt4nr5M7WJ3tNzaMIldVHgjGAAAAJaR9czYYjz9nmSZffIV6jxxSx2tbilxVeSAYAwAAlJD2LS/JXVWjwIxZmnbp+yRJB59+vMhVlQeCMQAAQInobT+q8IFW1S1aIpvNpknLL5Xd5dYBgrElCMYAAAAlItdfXHf2EkmSqzKgCeddpM6d2xQ+uK+YpZUFgjEAAECJyAXj+kVL81+beul7JEkH1z5RlJrKCcEYAACgRLRvfkmy2VS7cHH+a1MvWSlJtFNYgGAMAABQAjKplI5t26TqM+bLVRnIf71i4hSNO+tsHdnwvBKh7iJWePojGAMAAJSArl3ZjT3qBrRR5Exd8R4ZqZQOrXuqCJWVD4IxAABACWjbfPyNdwNNvfS9ktgFz2wnDcaZTEarV6/WOeeco1WrVqmlpeW47z/99NO65JJLtGzZMv3sZz/Lf/2SSy7R3LlzNXfuXK1cubLwlQMAAJxG8itSLHpnMK6eu0CVk6fp0HO/VyaZtLq0snHSYPyHP/xBv/zlL/WjH/1INTU1uuuuu/LfSyaTuuWWW/QXf/EX+vKXv6zbb79dbW1t6ujo0OHDh/XII49ow4YN+p//+R9TnwQAAMBY1755g9zBagUbZr/jezabTVNXvEfJUI+Ovry+CNWVh5MG440bN6qhoUELFizQsmXLtGnTpvz3du/era6uLq1YsUKXXXaZUqmUtm7dqq1bt0qSbr75Zn384x/Xtm3bzHsGAAAAY1zsWJvC+1uyG3vYB49nU1bk2ilYncIszpMdEAqF5PV6JUler1ehUCj/vXA4LEny+Xzy+Xz546dMmaKPfexj+uAHP6if/OQnuuGGG/T888/L6Rz6cjU1FXI6HaN6MidTXx84+UE4LTDW5YOxLh+Mdfkox7He/cozkqTpyy4c8vmPe+9KPR+s0uHnnlTdt78rm81mZYmmKLWxPmkw9vv9isVikqRYLKZAIHDc93Jf7+3tlSQFAgGde+65mjdvnvx+v1auXKlHH31UbW1tmjRp0pDX6eyMjuqJnEx9fUBtbaGTH4gxj7EuH4x1+WCsy0e5jvWe59ZJkirOWHTC5z9p+eVqffxRvfHHF1Qzd4FV5ZmiWGN9ojB+0laKRYsWae/evdqxY4fWr1+v5uZmRSIR9fb2qrGxUYFAQM8884zWrFkjp9OphQsX6jvf+Y4uueQS7d69W2vXrtWECRM0fvz4gj4pAACA00X7lr6NPZoWn/C4qSv6dsF7hl3wzHDSGeOVK1fqpZde0jXXXKMZM2boW9/6lj73uc+prq5Od999t+666y7deeedisViWr16terq6vTpT39aO3bs0Ac/+EFNmzZN3/nOd+RwmNsmAQAAMBZlUim1b92oqlnz5A5UnfDYCUvfJUk6tm2jFaWVHZthGEaxi5Bk+lR6ub41U44Y6/LBWJcPxrp8lONYd+x4Vb/7i4s064PX6vw77j3p8Y++a44cHo+ufmqrBdWZZ0y2UgAAAMA8ufWL6wdZv3gwNfObFDm0X/HuLjPLKksEYwAAgCJq3/ySJKnu7HduBT2YmnkLJUldO8f2jHEpIhgDAAAUUfuWl+UOVik484xhHV8zv0lStgUDhUUwBgAAKJJY5zGFWnerduG5Q27s8Xbj5mWDcefrzBgXGsEYAACgSHL9xXXD7C+WJP/0Rjl8FQRjExCMAQAAiqR9U19/8QiCsd3hUM2cs9S9e6fSibhZpZUlgjEAAECRHHn5j7I5HKpvHt6Ndzk185pkpFLqfvN1kyorTwRjAACAIkhFI+rYulHjzlwkV+XQa+sOpoY+Y1MQjAEAAIqgbfMGZVJJjV+yfMSPJRibg2AMAABQBEdffl6SNGHJhSN+bPWcM2Wz2wnGBUYwBgAAKIIjG/4om92u+sXnj/ixTl+FAg2z1fn6NhmGYUJ15YlgDAAAYLFUrFfHXn1FNfOa5A5UndI5auY1KRnuUeRga4GrK18EYwAAAIu1b9mgTDJxSv3FOTXzs1tD005ROARjAAAAix3d8EdJp9ZfnJPbAY+toQuHYAwAAGCxoy//UbLZVH/OslM+R/XcBZKYMS4kgjEAAICF0vGY2jZvUM2cs+Spqj7l8/jqxstXP1GdO7YVsLryRjAGAACwUPvWV5RJxDV+6an3F+fUzGtS9K0Dind1FKAyEIwBAAAsVIj+4hw2+igsgjEAAICFjryU3dhjNP3FOQTjwiIYAwAAWCSdSKh9ywZVnXGmvDW1oz5fzXyCcSERjAEAACxybNtGpWO9BWmjkCT/tJly+ioJxgVCMAYAALDI0b42igmj2NhjILvDoeq5Z6l7zy6l47GCnLOcEYwBAAAscuTl7I13488dfX9xTs28JhmplLp3v16wc5YrgjEAAIAFMsmk2je9pGDjXHlr6wt23vwNeDtopxgtgjEAAIAFOl7brFRvRBMKsH7xQOPmL8yen62hR41gDAAAYIHcMm3jC3TjXU7VGfNls9vV+To74I0WwRgAAMACR3Ibe5xb2GDs9PoUnDlHnTu3ychkCnruckMwBgAAMFkmlVLbxhcUaJgtX/2Egp+/Zl6TUpGQwgdaC37uckIwBgAAMFnHa1uUioYLtkzb2/XvgEef8WgQjAEAAEx2tK+NotD9xTnsgFcYBGMAAAATRQ4f0I4ffVc2u73gK1Lk5GaM9z/1W6UTCVOuUQ4IxgAAACZJRsJ69vqPKtZ+RM1//8+qGD/JlOt4x9Vp9oc+qe43d2j7D79tyjXKAcEYAADABJl0Wn/8+79S585tmv3hT2vetdeber3Ff//Pqpg4Vdt+8E1aKk4RwRgAAMAEm755mw6ufVITl63Qklv/VTabzdTrufxBnXfHPTJSKf3pH69XJpk09XqnI4IxAABAgb3xyEN6/UffU7Bxrt717Ydkd7ksue7k5Zdp1qpr1Pn6Vm3/97stuebphGAMAABQQIfXP6MNd35FnppaXfL9R+QOVlt6/cX/cKd8EyZr2/f/f3XuZDe8kSAYAwAAFEjnzm1a9+VPymZ36KL7fqLAtAbLa3AHq3Xe1/5NmVRSL9z6f2ipGAFnsQsAAAAYq9KJhNo3v6hDz6/RoXV/UNfO7ZKkZXc9oPGLzy9aXVMuvkKNV39Me/73p3rtwXu04HNfKVotY8lJg3Emk9HXvvY1/fa3v9WMGTP0rW99SzNnzsx//+mnn9Ydd9yhRCKhG264QR/+8IfV29urf/iHf9Dzzz+vhQsX6tvf/rZqa2tNfSIAAACFZhiGEj3dinceU7zrWPZzR7tiHe1q37JBb73wnFLRsCTJ7vZo0oWXqvHqj6nhfR8scuXS4pu+rsPrn9HW792lqZe9V9VnnFnskkqezTAM40QHPPnkk/rKV76ihx9+WHfffbdcLpfuv/9+SVIymdTy5ct1zTXXaOLEifra176mtWvX6rHHHtODDz6oH//4x7rhhhu0dOlS3XbbbScspK0tVLhnNYj6+sA7rpGK9aprhL03J3m5TDfYHa1W1FSoO2mtqLW6ukJdXdHhP2C4NY3iNRjq9Su1sTuunpHUZvKd1kMZ8VhrGH+GivxnvBQU+++5nIFjlRvrYdc21HEDvm5o6HPZNMjPyQmOt9QIf17zr5lhSMbbvvY2A19zwzAkI9P32Rj0MSc6fqjH9Ndt9J870/dZhgJ+t3q6e6VMRoaR6ft+Jn8+I5Pp+33ftTIZGZl039cNGUZGRjqdf7yRyQz4tZE9NpVSJpVUJpWSkU4pk8z+OhnpUTKU/UiE+34d7smeYwiBGbM0+V2Xa9LyyzVhyYVy+ipOPCAWO/jsk1p73UdkczhUOXma/FMbsh/TGuSfMkPeunrJZsv+zNuyH7aBn010xoUXqDOcMvUag6mvDwz5vZPOGG/cuFENDQ1asGCBli1bpgceeCD/vd27d6urq0srVqzQ5MmT9U//9E/aunWrNm7cqKamJs2aNUtLlizRpk2bCvNMCmzDHTdqz/8+XOwyAABACXFWBuQOBOWbMFlVs+fJU1UjT01t30edPDXj5KmpVdWseUXpIR6JKRdfqXNu/rpan/yVwvv36q0/rS12SXn7PnKtFn/13mKXcZyTBuNQKCSv1ytJ8nq9CoX6Z13D4exbBz6fTz6fL398KBRSTU3NoI8ZSk1NhZxOx8ifwQi8/X8I53/+Ro2bPn3kM0Qj+R+UYRRuNu1EdZr5v7pCzx4VotZTeV1P8JiT/a/4HTOKI/0ZOJFSG7sB9QxntqBUZheHZZh/hsyeJRnIMIwRz+oXor6TnqdYf88NPGf+l2+rdZjXGvL5nWSsT/TnfTh/V1jy8zPC1yNf03GzgG9/3Ntfc3t+5nDgLOKAg44/3t5/vM1mk81uP0ltxx+Xf6zd1neuvq8P+L7N4ch+z27PPr7vmjaHo/97Nnv+mMG+JrtdDpdLNocz+9nplMOZ/ewOBOX2B2R3mJtHrFb/dzfpXX93kyQpGYmoS1AwnQAAB3NJREFUZ/9edbfuUffeFsU6jx3/rsCAz6ay2TTnzz94wtnbYjhpMPb7/YrFYpKkWCymQCBw3PdyX+/t7ZUkBQKBEz5mKJ2dI3s7dKQGa6WwTT5Dcz93i6nXhfUGG2ucnhjr8sFYlw8rxzrXIJHu+xxPSOowN4+UhNrpCtROV2DxJUUto1h/rk8Uxk+6XNuiRYu0d+9e7dixQ+vXr1dzc7MikYh6e3vV2NioQCCgZ555RmvWrJHT6dTChQu1aNEibdu2TXv27NGGDRvU3Nxc0CcEAAAAFNpJb75Lp9O644479Jvf/Ca/KsXq1atVV1enu+++W2vWrNGdd96pWCymG264QR/5yEcUiUR0yy23aN26dWpqatLdd9/NqhQAAAAoaScNxgAAAEA5YOc7AAAAQARjAAAAQBLB+P+1d38hTbZ9HMC/Lz2Qk0lBhxJoxRa6qbfWok0yw2yFBzqok8hqRKKEB53MIe4ggjrVVUQRQUcZSAdR7ESzBkYF22xKMMr1hzpZRWtT53L7PQfRnp7e+eLN8F76fj8gg/FFLvhyXdePcY8REREREQHgYExEREREBICDMRERERERAA7GREREREQA/g8G42w2C4/Hg4aGBjgcDkSj0WIviVQ6evQojEYjjEYjzGYzotEoHA4HGhoa4PF4kM3++O0ir9cLi8UCu92OUCgEAKqyVDy3bt2C0WjEwsLCivW7VJa09WvX79+/z+1to9EIj8cDABgeHsbu3buxd+9ejI2NAQA+ffqEzs5OKIqC3t7e3K+tqsmSNpLJJLq7u6EoClpaWjA+Po5QKAS73Q6LxQKv1wtg6fu50CxpI1/PT548+deevnbtGoBVdn7LGufz+cRkMkk4HBan0yldXV3FXhKpkMlkpK6uTu7fvy/xeFy+ffsmXV1d4nQ6ZXp6Wqqrq8Xn80k4HBaDwSB+v1/cbre0tbWJiKjKkvbm5+fl4sWLsn37djEYDJJKpVas33xZ0k6+rh88eCAWi0Xi8bjE43GZn5+XWCwmVVVVcufOHRkcHBSLxSLpdFrOnTsnbW1t8urVK7FarXLjxg1VWdLO5cuXxWazSTQalb6+Ptm1a5e0tbWJ2+0Wv98vBoNBXrx4seT9XGiWtJGv5+vXr0t7e3tuTy8sLKy683vNf2IcCARQUVEBk8kEq9WKYDBY7CWRCtFoFHNzcxgcHMThw4fh9/sRCARgs9lQVVWFiooKBINBBAIBlJaWorGxEXv27EEkEkEymVSVJe3FYjFEo1H09PTk3lupfvNlSTv5up6amkI6nYbD4UBvby++fPmCyclJLC4uoqWlBc3Nzfj69StmZmYQCARgsViwdetWmEwmBINBVVnSzrFjxzA8PIyKigps2rQJmUwGkUgETU1NaGxsRGlpKUKhUN77OZlMFpwlbeTreWpqCh8/fkR7ezv6+/uRSqVW3fm95gfjRCKBkpISAEBJSQkSiUSRV0RqiAiOHDmCCxcu4NChQ3C5XIjH4//VaSKRgE6nA4DcazKZzNv/UlnS3ubNm3H16lWUl5fn3lPTWaFZ0k6+risrK3Hy5ElcunQJ8Xgc58+fz/Wi0+ly/f3sNV/Xy82SdsrKylBeXg6fz4ebN29i//79AJC3v6U6LSRL2vi9587OTlRXV+P06dMYHBxEOByG1+tddef3Xyv2n/8Qer0eqVQKAJBKpVBWVlbkFZEa27Ztg8vlgl6vx4YNG3DlyhUAwMLCAoB/Ov2155/PE+r1euj1+mVn6c+gprNCs1RcHR0dSKfT0Ol0sNls8Pl8cDgcAH509LO/n73m63q5WdLWvXv34HK5cPDgQbjdboyMjPzPvfp7p4VkSTu/9tzT04NsNotsNov169dDURREIhHs27dvVZ3fa/4T49raWrx58wYvX77ExMQEFEUp9pJIhbt378JisSAUCmF0dBQ6nQ5WqxUTExOYnp7G27dvoSgKamtrMTs7i4mJCTx69AgGgwF6vR41NTXLztKfQU1nhWapuE6cOAGn04kPHz7g6dOnMJvNMJvNWLduHcbGxvDw4UNs3LgRW7ZsQU1NDZ4/f47Xr18jHA5DURRVWdJOKBRCX18fmpubMTAwABFBZWUlHj9+DL/fj7m5udxe/f1+LisrKzhL2vi959nZWRw4cAAulwvv3r3D5OQkzGbz6ju/V+zp5T/E4uKieDweqa+vl46ODpmZmSn2kkiF79+/S39/vyiKIq2trTI+Pi4zMzPicDhEURQZGBiQTCYjIiJDQ0Oyc+dOaW1tlWAwKCKiKkvFMzIykvtC1kr1u1SWtPVr15FIRBwOh9TV1cmpU6fk8+fPIiJy+/ZtsVqt0tTUJKOjoyIiEovF5Pjx41JXVydnzpyRubk51VnSRnd3txgMhn/9BYNBsdvtsmPHDhkaGhKRpe/nQrOkjXw9P3v2TOx2u9TX18vZs2dldnZWRFbX+f0fEZGVG7uJiIiIiFaHNf8oBRERERHRcnAwJiIiIiICB2MiIiIiIgAcjImIiIiIAHAwJiIiIiICwMGYiIiIiAgAB2MiIiIiIgDA31juVJTENjWkAAAAAElFTkSuQmCC\n",
Parish, Chad's avatar
Parish, Chad committed
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
      "text/plain": [
       "<Figure size 864x504 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('seaborn'):\n",
    "    fig = plt.figure(figsize=(12, 7)) \n",
    "    gs = gridspec.GridSpec(2, 1, height_ratios=[4, 1]) \n",
    "    ax0 = plt.subplot(gs[0])\n",
    "\n",
    "    ax0.plot(depth_in_nm, NRT, color='xkcd:brick red')\n",
    "    ax1 = plt.subplot(gs[1])\n",
    "    ax1.plot(depth_in_nm, atom_frac, color='xkcd:brick red')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
Parish, Chad's avatar
'minor'    
Parish, Chad committed
371
372
373
374
375
376
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T22:44:59.424641Z",
     "start_time": "2020-05-12T22:44:57.562631Z"
    }
   },
Parish, Chad's avatar
Parish, Chad committed
377
378
379
380
381
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
Parish, Chad's avatar
'minor'    
Parish, Chad committed
382
383
384
      "Uncertainty in fluence (for \"+/-20%\" input \"0.2\"): 0\n",
      "Uncertainty in displacement energy (for \"+/-20%\" input \"0.2\"): 0\n",
      "Combined error: 0.0\n"
Parish, Chad's avatar
Parish, Chad committed
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
     ]
    }
   ],
   "source": [
    "try:\n",
    "    err_fluence = float(input('Uncertainty in fluence (for \"+/-20%\" input \"0.2\"): '))\n",
    "except ValueError:\n",
    "    err_fluence = 0\n",
    "\n",
    "try:\n",
    "    err_displace = float(input('Uncertainty in displacement energy (for \"+/-20%\" input \"0.2\"): '))\n",
    "except ValueError:\n",
    "    err_displace = 0\n",
    "err_both = np.sqrt( err_fluence * err_fluence + err_displace * err_displace )\n",
    "print(f'Combined error: {err_both}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now, generate the figure\n",
    "This does require some hand-tweaking of parameters on the first attempt, but all subsequent uses will take only a few seconds, where most of the time is spent waiting for the user to click on the new data directory, or to type in the fluence or error."
   ]
  },
  {
   "cell_type": "code",
Parish, Chad's avatar
'minor'    
Parish, Chad committed
412
413
414
415
416
417
418
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T22:45:04.497537Z",
     "start_time": "2020-05-12T22:44:59.426655Z"
    }
   },
Parish, Chad's avatar
Parish, Chad committed
419
420
421
422
423
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
Parish, Chad's avatar
'minor'    
Parish, Chad committed
424
      "File name (enter to skip):5hours\n"
Parish, Chad's avatar
Parish, Chad committed
425
426
427
428
     ]
    },
    {
     "data": {
Parish, Chad's avatar
'minor'    
Parish, Chad committed
429
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAu4AAAG7CAYAAACcthATAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXxU1fnH8c/JMpMQMhMSdhARWQKiqCwiuFEriiuL1QLa4lKwVq1aa7W2bq3+ftX+urgUxRVRcGVR3AqCgLgCAiIKYQkIIkIyIQEySUjO74+ZDAkkJCEzcyfD9/16zSsz95x77jPRP54cnvtcY61FRERERERiW4LTAYiIiIiISN2UuIuIiIiINAFK3EVEREREmgAl7iIiIiIiTYASdxERERGRJiDJ6QBi1XnnnWffe+89p8MQERERkfhn6jNJO+612Llzp9MhiIiIiIiEKHEXEREREWkClLiLiIiIiDQBStxFRERERJoAJe4iIiIiIk2AEncRERERkSZAibuIiIiISBOgxF1ERERiUk5ODhdddBEtWrSgQ4cO3HjjjezZsweAPXv2cOWVV9K8eXPatGnDgw8+WO3cSI7Xde6BjDGkpKQcdDwlJQVj6tW+u0HrNnWdO3fGGMMPP/zgdCh1WrBgAcYYrrvuuuhc0FqrVw2vvn37WhEREXFGcXGx7datmwXswIEDbefOnS1gr776amuttddee60FbM+ePW2bNm0sYF988cXQ+ZEcr+vcAwHW7XYfdNztdttAKnZ4LrnkEvuzn/3ssM+PVUcffbQF7LZt25wO5ZByc3NDsU6YMKGxy9UrP3U8QY7VlxJ3ERER58ydO9cC9uKLL7bWWltQUGBTUlJsSkqK3b17t3W73bZDhw62tLTUfvXVVxawp512mrXWRnS8rnNrEqnEPV41hcT9qaeespmZmRaIauKuUhkRERGJOT169ODFF1/k9ttvB8Dj8ZCamkppaSnLly+npKSE/v37k5ycTO/evcnIyGDJkiVUVFREdLyucxvj+++/Z/jw4aSlpdGqVStuuukmSkpKap1/YKlMcXExt956K+3btyclJYWBAwfy4YcfVpvfu3dvXn75ZY455hiysrIYN25ctWu88MIL9OrVi9TUVNq2bcv48eMpLi4GYMvHC1nx7MQGfafdu3czadKkBp1TH4sXL2bgwIGkpaXh8XgYOnQo69ev55xzzsEYw9tvvx2ae8MNN2CM4ZlnngEO/XvOzc3FGMPQoUMZPnw4Ho+Hxx9//KDr//Wvf8XtdkevRCZIibuIiIjEnI4dOzJ27FgGDx4MwKuvvorP5yM7O5vvv/8egMzMzND8zMxM/H4/+fn5ER2v69zalJWVMXz48GqvsrKy0Li1lhEjRjBr1ix69+5Nu3btePTRR7n++uvr/TsbNWoU//znP0lMTOSUU07hiy++4JxzzmHx4sWhOZs2bWLChAl07tyZoqIiJk+ezAsvvADAypUrGTduHNu2bWPIkCGkpaXx1FNPhf54Wj7pET7537sPum5paSm5ubk1vh5//HEmTJjAH/7wh3p/j7rs3r2biy66iKVLl9KvXz/at2/PnDlzuOWWWxgzZgwAM2fODM1/6623SE5OZsSIEfX+Pc+dO5dPP/2U7t27c8oppxwUw5133snq1atrHIukpKheTURERKSBvvjiC6699loAbr75Zvx+PwDJycmhOUlJgZSmuLg4ouN1nVubiooKZs2aVev4/Pnz+fzzzxk9ejRTp04FYPDgwbzwwgs8/PDD1f5QqMnixYt599136d69O8uWLSMtLY1JkyYxYcIE7r77bj744AMgkPTOnj2bCy64gLvuuosHH3yQpUuX8qtf/Yrc3FystVxyySVMnDiRkpISnnnmGQYMGACAb0MOJYW7sNZWu6l21apVDBo0qNbY3G43Dz30ENZaHnrooUN+j/ooKSnhvvvuo3Xr1lx++eXk5+eTlZVFTk4OI0eO5Prrr+ett96ioqKClStXsnnzZs4//3wyMzOZN2/eIX/Play1LF68mGOPPbbGGCZMmNDo73E4tOMuIiIiMevrr79m2LBh7N69m5EjR3LttdeGykPKy8tD8yp3r1NTUyM6Xte5tXG73QfVK7vd7mrfE2DatGkYYzDG8PHHH7Nv3z6WL19e16+JTz/9FIARI0aQlpYGwJVXXlltrNK5554LQM+ePQFCZSJDhgyhd+/eTJ48mczMTEaNGkVSUhJ9+/bFVlRQsHEdtrycfXv3Vlvv5JNPxu/3H/J1wQUX8PDDD4e6AjVGVlYWo0aNYtOmTVx44YV06dIFCPxR5fV6Of/889m+fTuffPJJ6I+lyy+/HKj/77l169a1Ju1O0o67iIiIxKStW7dy7rnnkpeXx9lnn83UqVMxxtC2bVsAfD5faK7P5yM1NZWsrKyIjtd17uGqTP6zs7Pp0aNHtbFD/UFQKSGh9r3YqrvjLpcr9C8ElT+ttQCkp6fz+eef88Ybb/DOO++wYMEC5s2bxwsvvMD8WTMoD/5rQ2lRIcnBPw4AlixZQv/+/euM8U9/+lPoj4rGyM3N5dRTTwXgpptu4q677mLQoEGh7zFmzBimT5/OzJkzmTdvHm63m0suuQSo/+/Z4/E0Os5I0I67iIiIxBxrLWPGjGHr1q307duXWbNmhXao+/TpQ1JSEp999hllZWWsXr2agoIC+vXrhzEmouN1nXu4evXqBUDXrl2ZOXMmM2bMoEePHpx22mkcf/zxdZ5/0kknATBjxozQrvaUKVMAQkkucMgY3333Xa6++moSEhKYOnUq3333He3atWPZsmVsXLYkNK+kaFe180488US2bdtW4+sf//gHEEja//KXv9TnV1Gn119/nR9++IGxY8dy5513kpiYWG38wgsvxOv18uKLL7Js2TLOO+88vF4vUP/f84FrxgrtuIuIiEjMmTNnDgsXLgx9Hjt2bOj9888/z+jRo5kyZQp9+vQJ3RRa2eHD4/FEbLyucw/X0KFDyc7OZvbs2fTp04eKigpWrVpF//79ue222+o8/4wzzuDMM89kwYIFZGdn06VLFz766COSk5O5//776xVDZmYmr7/+Om+88QbPPvssPp+Pbdu2BZLdndtD80oKqyfuSUlJoX+JONCECRPweDxcc8019Yqhqp/+9KehfxWo9Morr3D00UcD8Oijj7JixYpQKdDu3buBQFnSiBEjeP7554H9ZTJQ9+95586dDY4zqurbN/JIe6mPu4iIiHNuvPHGUI/sA1/btm2zRUVF9he/+IVNS0uzrVq1sg888EC18yM5Xte5B6Kefdw3btxoL7nkEtu8eXPr9XrtqFGjDtnL/MB1i4qK7G9/+1vbvn1763a77cCBA+2CBQtqnT9t2jQL2F/+8pehY7Nnz7b9+/e3aWlptmXLlnb48OF23bp1dv4fb7b/bJlg/9kywebOe/+Q37exKvu41/T68ssvbUVFhb3xxhttRkaGbdWqlb3mmmvsoEGDLGDXrl1rrbX2/ffft4BNTU21RUVF1dY/1O9548aNFrA9evSoV6zPPfdcVPu4GxusB5Lq+vXrZ5csWVL3RBEREZE4N3P0hXz30XzK/X7Of/plul/yM6dDOqR//etf3HLLLYwaNYrXX3/d6XDqo151VqpxFxEREZFDKlifQ6veJwJQekCpTCx57733GDFiBHfccQfgXNvGSFHiLiIiIiK1Ki8tZdfmjbTp0xeAkqJChyOqnTGG//73v7Ro0YL/+Z//4ZxzznE6pLDSzakiIiIiUqtdmzdiy8tp3edkMCamd9zPPffcsPSKj1XacRcRERGRWhVsyAEgs1s2rubpMb3jHu+UuIuIiIhIrQrWBxL3Fsd2x+3xxvSOe7xT4i4iIiIitfKtX0tKi0xSWmTi8ngPegCTRI8SdxERERGplW99DhnHdgfA1dxDaaFKZZyixF1EREREalWwfi0tunQDwO3xHPTkVIkeJe4iIiIiUqOyPXvYvW0rLSp33FXj7igl7iIiIiJSo4KN6wDIOLZyx93L5rw8hgwZQlZWFunp6VxwwQVs3boVgJycnFrH4oWT31GJu4iIiIjUyLd+LQAZwVIZl8fLzsIi9u3bx3333ccVV1zBO++8w8033wzA1q1bax2LF05+Rz2ASURERERqVNkKMuOYrgC40710smX8Y+5cktxuAKZMmcLXX38NwKBBg1i0aFHo/KpjNTHG4Ha78fv9kfoKNRo7diw5OTl8/vnnDT63od8xnLTjLiIiItIICxcu5MILL6RNmzYYYzDGcO+999br3Keffpp+/fqRmZmJy+Wiffv2XHDBBdUSQyfifPnllzn55JM5/fa7ucsHY8aNY926dbg9XpKqPD31o48+Ys+ePQwePBgAl8sVWuPAsZpccsklXHzxxeH7kvX08ccfHzKuQ2nodwwnJe4iIiIijbBs2TLee+89MjMzG3zu4sWL+e677+jUqRPZ2dns2LGDd955h6FDh5Kbm1vjOe+99x7ffPPNQcenT5/Opk2bGh3nM888w+jRo/nyyy/JcCVhjeGNN95g8ODB7CovB6CkcBdff/01o0aNokePHvztb3+rtsahxqqaOXMmr7766iHjCbcffviB3NxcBg0a1Kh16vsdw0mJu4iIiEgjXHnllRQWFvLFF180+NyJEyeyfft2li9fzsqVK3niiScA8Pv9LF269KD5BQUF/PznP+fss88mJycndPy1117jsssu45prrmlUnKWlpdxxxx0AjBo1ivuO8vDs1aNJT0/nxx9/5Jm33wdgxZdLGTJkCOnp6cydO7faHwOrVq2qdexAxhhSUlJCn4uLi7n11ltp3749KSkpDBw4kA8//PCgc3r37s3LL7/MMcccQ1ZWFuPGjaOkpASAF154gV69epGamkrbtm0ZP348xcXFofMXL14MBEpecnNzMcYwfPhw7r//frKysmjbti0vvfQS7733Hl27diUtLY1x48ZRWlp6WN8xrKy1etXw6tu3rxURERGpr6KiIgtYwN5zzz31Pm/BggX2lFNOsb1797bJyckWsCkpKTY3N7fG+bNnz7Yul8t27NjRbtiwwU6fPt0mJSXZTp061XpOfeP86KOPQmPPT3rS/rNlgl3y+P/Zc845xwK2y1FH2btbGJuVkWETExPtQw89ZKdNm2anTZtmrbV28+bNtlWrVjWO1QSwbrc79HnYsGEWsB07drRnnHGGTUhIsElJSfajjz6qdk7z5s2tx+OxZ511Vuh3NmnSJLtixQprjLEZGRl22LBhtkuXLhawN9xwQ+j8W2+91R599NHWWms3btwY+n1nZWXZfv36hWJq1qyZPeuss2yzZs1C6x/Od6yneuWnjifIsfpS4i4iIiINcbiJ+4wZM0LnAbZ169Z24cKFdZ6TlJRk27dvb5OTk23Hjh3t+vXrGx3ntGnTQmOvPP6I/WfLBLv+vTftFVdcEUhw3W77G4+pFm/ly1pr58+fX+tYTaom7pV/NHTv3t3u3r3bWmvtk08+aQH7k5/8pNo5gJ09e7a11to//vGPFrATJkyws2bNsoD95S9/affu3Wt9Pp/9+9//Xu33eeqpp9oxY8ZYa/cn7omJiTYnJ8dWVFTYY489ttrv5i9/+YsF7E033XRY37Ge6pWfqlRGRERExEHDhw+noqKCbdu2ccMNN/Djjz8yZswYNm/efMhzbrvtNr7//nvKysqYPHkyXbp0aXQsgbw4YPe2QG/yjGO7Vzve1WVYNfW5g3eDgbPOOqvmneJ6+PTTTwEYMWIEaWlpQKC8p+pYVeeeey4APXv2BKCkpIQhQ4bQu3dvJk+eTGZmJqNGjSIpKYm+ffuG5ixbtuyg+vaOHTvStWtXjDG0b98egDPOOAMg9LmyFKcx37GxlLiLiIiIOMwYQ9u2bXnggQcA2LJlS6jevSZz5szhX//6F1lZWSQkJDB+/PiwPASoU6dOofff5azFJCbi7XQMP/74IxBIcAFKCwsbfa0DJSTUnpYaY6p9drlcJCUFuppX/rTWkp6ezueff86UKVMYMWIE3377LTfffDOnn346AEuXLqWkpOSgLjBV6+wr42jWrFmdcUVb7EQiIiIiEqfuvPNOsrOzOfvss0PH9u7dy1NPPVXtxsm33nor9H7Pnj01rjVv3jwuueQSPB4PCxcu5LnnnmPjxo385Cc/Ydu2bY2Ks3///mRlZQHw38++wHNUZ7bv3Bna8T7vvPMAKCna1ajr1OSkk04CYMaMGaHvPmXKFABOPfXUanMPTOQrvfvuu1x99dUkJCQwdepUvvvuO9q1a8eyZcvIy8vj448/pnnz5hx//PFhjz8alLiLiIiINML06dPp2rUrffr0CR175JFH6Nq1K2PHjgVg27ZtrFmzhvXr14fmlJaWMn78eFq0aMFxxx1H165dueKKKwBITk4OnVtVQUEBI0eOpFmzZsydO5devXrxi1/8gqeffpqcnJxQacnhxulyuXjwwQcB+Dh3C3d+s4WePXtSVFREy5YtufOPfyS5WVqoj3s4nXHGGZx55pmsXbuW7OxszjzzTH7961+TnJzM/fffX681MjMzef311xk3bhw//elP6d+/P9u2baNXr15kZWXx8ccfM3DgQBITE8MefzTEdeJujOlhjHnFGLPdGFNmjNlqjHnGGNPe6dhEREQkPhQWFrJ+/Xo2bNgQOubz+Vi/fv0hy1dSUlK44oorOOqoo8jNzWXTpk20a9eOkSNHsmjRIgYMGHDQORkZGTz33HPMnTu32q7xVVddxeTJk5k4cWKj4xw/fjxTpkyhY3IC+f5SjDGMGDGCxYsX0759e1zpHkoikLgnJCQwe/Zsfvvb31JRUcFnn33GgAEDmDt37kE77rU55ZRTmDlzJieeeCKffvopmzdvZvjw4bz55psAfPLJJ43u3+4kE61i+mgzxvQAPgO8NQx/B/S11u6o7fx+/frZJUuWRCo8ERERkZi1+4fvefr4oxjyt8foc/Wvq41NHnQcLbOP44Jno/vgpDhXc+3PAeJ5x/1u9iftzwFXAG8GPx8F3OxEUCIiIiKxzrd+LQAtju120Jjb441IjbvULcnpACKob/DnWmvt1QDGmNeAPKA5cLJTgYmIiIjEsoJg4p7R5eDEPVKlMlK3eN5xz6vhmGH/P0X8GMVYRERERJoM3/ocEt1u0jscddCY2+OltCj87SClbvGcuE8K/uxujHnaGDMGeBVIA0qBxx2LTERERCSGFaxfS0aXbpgaepi7Pd6IdJWRusVtqYy1drIxpjnwb+Ca4AugABhlrf3cseBEREREYphvwzqyuvesccylGnfHxO2OuzGmJTAWOLBRZwZwZzCpP/Cc8caYJcaYJTt21NpwRkRERCRuVezbx67c9TXWt0Ogxn3f3r2Ul5VFOTKJ28QdeB44FSgBLgY8wHXBsZ8Cjxx4grV2krW2n7W2X6tWraIVp4iIiEjMKPxuExVlZTV2lAFwpwea9qnOPfriMnE3xrQDzg9+fM1a+5a1tsha+ySwOHh8rDEmbkuFRERERA5HZSvIjGO71zju8ngA+GbFci666CJatGhBhw4duPHGG9mzZw8AX331FWeddRbNmjWjR48ePP3004e85p49e7jyyitp3rw5bdq0CT29taG6d++OMYZXXnml2vGSkhLS09MxxvDNN98c1tqxIOYTd2NMkjGmjTHmqgac1pn93WO2HTD2ffCnC9C2uoiIiEgVBRvWAdCilsTd7fFSZi0/G3c1s2fPJjs7G5fLxWOPPcZNN92E3+9n2LBhLFy4kAEDBrB9+3Z+9atf8cYbb9R6zZtvvpkXX3yRTp06YYzhrrvu4qWXXmpw7KNGjQJg5syZ1Y7Pnz+f3bt307NnT3r2rLl2H+CEE07AGFPj64MPPmhwPOHmeOJujEk0xjxkjNlijCk1xpRXfREodfkeOPSfatVVLVA/vcq1koB+wY8l1NwyUkREROSIVbB+LS6Pl9SWNe9vuj1eNpbBhs2bufjii/nkk09Yvnw5KSkpTJ06lc8++4z8/Hyuv/56PvzwQ5588kkAXn/99RrX27NnD1OmTKFDhw6sWLGCuXPnAvDEE080OPZLL70UgHfffZeyKjX4b731VrXxmvj9/tBufIcOHTjllFNCr4EDBzJgwIAGxxNusVAqcjtwG4FEOgfoRiDx3gUcDaQS6Lle738zsdauM8asBE4ABhpjXgLeA34GHBOcNttaWxquLyEiIiISD3wbcmjRpRvGmBrHXR4vrRPh/35/C6dcEtjh9ng8pKamsmvXLgYMGEBRURF+vx+A7du3A5CZmVnjesuXL6ekpIT+/fuTnJxM7969ycjIYMmSJVRUVJBQQ0vK77//nuuvv545c+bQrFkzRo8ezcMPP0zfvn3p3Lkzubm5zJ8/n6FDhwIwe/ZsYP+OfE1WrFjBvn37APjTn/7EddddV+tcpzi+406g80s+0NVaexywAJhvre0FZBHYaW8FrG7gutcCRcH3Y4AXgIuCn7cBtzYybhEREZG4s3fHj6S1aVvruDvdS0aiYVi/kxg8eDAAr776Kj6fj+zsbFJTU0lMTCQtLY3LLruMW2+9lc6dO3PHHXfUuN733weqmKsm9pmZmfj9fvLz8w+ab61lxIgRzJo1i969e9OuXTseffRRrr/+emB/cj5r1iwgkJBv3ryZrl270qdPn1q/19KlS0PvjzvuuFrnOSkWEvdjgOnW2q3Bz18AZwBYa/0EOsGsBW5pyKLW2i+Akwkk7D8A+wgk7M8A/a21m8MSvYiIiEgc8fvySGmRVeu4yxPoKlNSGOgq88UXX3DttdcCgVr1ShUVFcycOZPy8nJatWpFUVHRwYtBaGc+OTk5dCwpKVAUUlxcfND8+fPn8/nnnzN69Gg+++wzVq5cyaBBg3jhhRfIz88PlcPMmjULa229ymSgeuJ+xhlnVKtvX7BgAQD33XcfxhgSEhJIS0ujb9++TJ8+/ZDrhlMslMpA9Zr0tUAHY0yGtbbAWlthjHkfGN7QRa2164BfhitIERERkXjn9+WRknmIxD090FWmtHAXX3/9NcOGDWP37t2MHDkylMADGGPYuXMnixYt4qKLLmLkyJF8++23B62XkpICQHl5eehYZX16amrqQfO//vprAKZNm8a0adOqjS1fvpwhQ4bQsWNHtmzZwpIlS3jzzTeBQ5fJwP7E3ePxVLuB1RhDv36BWyRXrVoFwCOPPEJRUREPPfQQl156KYsWLQr960MkxULivgmoetvyuuDP44FFwff7UAcYERERkYjaV1zMvuLiQ+64J7ndJLrdbNm6hZvOPZe8vDzOPvtspk6dGqqL37FjBy6XC6/XywUXXEDHjh1Zs2YNO3fupGXLltXWa9s2UJbj8/lCx3w+H6mpqWRlHRxHZVKfnZ1Njx49qo2lpqZijGHkyJE88sgjPPnkkyxZsoTOnTuHku+alJSUsHp1oCp7+PDhTJ48ucZ5X331FS6Xi1//+tckJibi9Xr5zW9+wyuvvBKVxD0WSmXeBi42xlxljEkEvgSKgd8AGGMyCOy2b3EuRBEREZH4V+wLNNxLreVG0kqudA/3vjKTrVu30rdvX2bNmoXb7QZg0qRJtG7dmnvuuQcI1LBv376dtLS0Gm9Q7dOnD0lJSXz22WeUlZWxevVqCgoK6NevX403yPbq1QuArl27MnPmTGbMmEGPHj047bTTOP7444H9ZTHPPvss1to6d9tXrlwZ+oPgxBNPrHFOSUkJ69ato3v37iQmJgLQrVvgIVXbth3YfTwyYmHH/X+BkQRuQk2y1j5ljJkE3GyMORNIIfDU03scjFFEREQk7vnzA4n7oXbcAXJsMl9t3b+nOnbs2ND7SZMm8ec//5l///vfLFu2jJycHEpLS/nDH/4Q6hAzfPhwXC4Xr776Kh6Ph9GjRzNlyhT69OkTuiG1tq4uQ4cOJTs7m9mzZ9OnTx8qKipYtWoV/fv357bbbgNg8ODBtG3blh9++AGof5kMwMSJE6s9wOmYY45h2rRprF69mvLycrKzs0NjGzZsAODYY4895Prh4njibq3NM8acBPyawI2pAHcCbmA0gd33pwgk+CIiIiISIX5fZeJ+6B33VXv390hfunRptcT3iSeeYM6cOdx66618/vnntGrVigcffJDbb789NKfqDj3Af/7zH4wxvPHGGzRr1owHHniAMWPG1HjthIQE3n33XW6++WY++OADEhMTGTVqFI899li1OSNGjGDixIl06NCBgQMHHvL7VI0/Jyen2liXLl0C3zlY315Z/15SUsJjjz2GMYaf//znh1w/XIy1NioXamr69etnlyxZ4nQYIiIiIlGT8+brvH3N5Yz98EtaHXdCrfPeGDWUfcXFXP7OolrnxJvbb7+dhx9+mOuuu44ePXowefJkli9fzu9+9zv+/ve/N3b5mpvmH8DxHfdKJlDENBDoA3iBncCn1tqvHQ1MRERE5Aixv8b90KUy7nQPe3/8IRohxYzKHfennnqK9PR0srOzefbZZ7nqqquiFkNMJO7GmJ8CE4EulYeCP60xZgkw3lq7wpHgRERERI4Q/mB9eV017i6Pl5LCXdEIKWa88847TofgfOJujBkEzCbQ4eYl4HMCD0zKAAYReOrpPGPMIGvtGscCFREREYlzfl8eSc2akRTsrV4bV7qH0iMscY8FjifuwH0E+rQPCT7ttKqnjTFPAx8A9wOXRzs4ERERkSNFcX4eqXXstgO4PV5KdxdRUV5OQrA1okReLPRx7w+8VkPSDoC19mPgdeAnUY1KRERE5Ajj9+XVWSYDgcQdoGzP7kiHJFXEQuJeChTVMaeQ2IhVREREJG75ffmk1HFjKgRKZYAjrs7dabGQDE8DxhpjOtc0aIxpS+ABTS9HMSYRERGRI44/P6/OHu6wf8ddde7RFQs17k8Bg4EvjTGPAAuBrUAqgTKa3wOJwAJjzMVVT7TWvhnlWEVERETiVrEvr85WkBDoKgPacY+2WEjcVwKWQAvIPwffV1XZGnLaAccsgYReRERERBrJVlRQUuBrUI27dtyjKxYS9/s5OFkXERERkSgq2VWAraioV6lMqMa9SIl7NDmeuFtr73U6BhEREZEjnd9X+XfCBoYAACAASURBVPCletS4p6tUxgmxcHOqiIiIiDis2JcHUK8+7q5QqUxhRGOS6qK+426MefYwT7XW2mvCGoyIiIiIAIGOMkC92kEmpaaSkJREaZES92hyolRmXA3HKmvcTS1jlTejKnEXERERiQB/cMe9PjenGmNwebwqlYkyJxL3kw74nAm8BOQTuFH14+D75gTaQd4DZAAXIyIiIiIRUZwfqHGvTztICHSWKdXNqVEV9cTdWrui6mdjzDPAPuBMa21elaG9wNvGmEXAlwSS+suiFqiIiIjIEcTvy8MkJOD2ZtRrvitdO+7RFgs3p44AZh2QtIdYawuBt4BzoxqViIiIyBHE78vDndECk1C/9NDt8erm1CiLhcTdEiiFOZSOgD8KsYiIiIgckYrz8+rVUaaSKz1dfdyjLBYS94+Anxljzqxp0BgzEhgO/DeqUYmIiIgcQfy+/Hp1lKnk1s2pUef4A5iAPwFDgDnGmHeBpUAR4AUGAz8BtgN3ORahiIiISJzz5+fRvH2Hes93pXvVDjLKHE/crbVfGWPOAP4NXBR8hYaB94AbrLWbnYhPRERE5Ejg9+XRstfx9Z4f6CpTiLUWY2rq6C3h5njiDmCtXQ6caYxpB/QBWgA+4Etr7XZHgxMRERE5AhT78hpUKuPyeLHl5ZTt2YOrefMIRiaVYiJxr2St3QZsczoOERERkSPJPr+ffXv3NujmVLfHA0Bp0S4l7lESCzenioiIiIiD/AWBhy+lZGbW+xxXuhdAN6hGkRJ3ERERkSOcPz/wOJ2UBrWDDOy4K3GPHiXuIiIiIke4Yl8gcU9tYDtIgFIl7lGjxF1ERETkCHc4O+6hxF0tIaNGibuIiIjIEc4f3HFvaFcZQE9PjSIl7iIiIiJHuOL8wM2pDesqU1kqox33aGkSibsx5pfGmHnGmHeNMf9rjElvwLmnB88rCL4+NMYMiWS8IiIiIk2J35dHUmoqSamp9T4nOa05JiFBN6dGUZNI3IHOwFnARKAT8GJ9TjLGjAMWAOcB3uDrTGCuMWZoBOIUERERaXL8+XkNqm8HMMbgSvfo5tQoiqkHMB3C88CH1toFwJvGmGPrOsEY04VAom+AfOBewAX8FUgB/gH0jlC8IiIiIk2G39fwxB0Cde7acY+eJpG4W2s3AZuqfF5fj9N+TSBBBxhjrX0fwBiTRmDXfasxprm1dne44xURERFpSorz80lpUf+HL1Vyp3t0c2oUxVziboxpB2Raa782xiRZa/cd5lLnBn/uBuZUHrTW3t/YGEVERETiib8gn5bZxzX4PFe6VzenRlFM1LgbY1KDN51uB7YAK4JDvwvelNqjgeslApXnrAV+YYzZYIwpNsYsMsYMCl/0IiIiIk2b35fXoFaQldwej/q4R5HjibsxpjmwCLgd8AMbCNSlAzQjcFPqImNM5wYsm0Ggnh3gWOA54BgCpTOnAfONMWc2MnQRERGRJs9ai9+Xrxr3JsDxxB34E3AycBOB7jEvVQ5Ya+8BfglkAn9uwJrNqrz3EqiPH0/gxtQyAkn9U8aYat/fGDPeGLPEGLNkx44dDf8mIiIiIk1MSeEubHk5qYe14+6lVDXuURMLNe6XAe9Zax8DMMbYqoPW2inGmFFAQ3qvlxzweYy19uPg+mnALUA34CRgaZVrTQImAfTr188iIiIiEuf8+cGnph7Gzamu9MCOu7UWY0zdJ0ijxMKOe3tgeR1z1gDtGrBmAVARfF8OfFZlbFGV910bsKaIiIhI3PH7gon7Ye64V5SVUe73hzssqUEsJO47gJ51zOkdnFcv1tpSYGPwYyL7690hUCpTqby+a4qIiIjEo+LQjvth1LinewBU5x4lsZC4vwVcZIwZVtNgsExmGPBOA9f9sMr7C6u871fl/bcNXFNEREQkrlTuuB9ejXswcVede1TEQo37fcBFwFvGmHeALABjzL0EkuxhwI/AXxq47hPA1QQ61DxljOlI4Pv+Pji+wlq7qtHRi4iIiDRh/vx84DB33D1eALWEjBLHE3dr7fZgX/UngAvY3wry7uDPhcB4a+3WBq67xBjzl+A6XuAfVYaLgGsaFbiIiIhIHPD78sAY3N6MBp/rTg8m7iqViQrHE3cAa+13wAXGmLYEWkNmEHji6UprbW4j1r3HGPMNcCtwPFAMzAP+bK39ptGBi4iIiDRxxfl5uL0ZJCQmNvjcyh131bhHR0wk7pWstT/Q8Fr2utZ8GXg5nGuKiIiIxAt/QT6ph1EmA4GuMqAd92hxPHE3xtxd9yzKgb3Ad8BHwQRfRERERBrJ78s7rFaQUGXHXTXuUeF44g7cC1Q+7Kimzv0HjpUZY+621v4t0oGJiIiIxDt/fj7NWrc5rHNdzdMBlcpESyy0gzyeQM/1fOAu4CygO4GOMjcAm4GdwOUEbij9CnjQGDPCiWBFRERE4kmxL++wWkECJCQmkpzWXKUyURILO+43Euj6cpK1dssBY8uMMW8Cy4DTrLW/Nca8DHwN3ATMiG6oIiIiIvHFn593WK0gK7k9XrWDjJJY2HEfCUyrIWkHIHh8OjA6+LkYmA2cELUIRUREROJQeWkpZXt2Nypxd3m8KpWJklhI3F1AszrmJB4wZ2/wPBERERE5TP78wFNTUzIzD3uNwI67EvdoiIXEfTkwwhjTpaZBY0xnArvyK6scHgBsinhkIiIiInGs2BdI3A+3HSSAK91DSaFKZaIhFhL3+wEP8IUx5h5jzPnGmJONMWcZY24HPiJQA/8AgDHmGeBM4DXHIhYRERGJA/t33BtZ465Smahw/OZUa+08Y8zlwH+Ae9jf/hECLSDzgbHW2reNMe2Bq4DPgH9EPVgRERGROOIP7rirxr1pcDxxB7DWTjfGvAsMI9AGsiVQCCwFZllr9wanFgH9gWXWWlvjYiIiIiJSL8WhxL0RNe7pHkpU4x4VMZG4Q6hbzPTgq7Y5RQSSeRERERFpJH9+PtDYGncv5X4/5aWlJLrUOySSYiJxN8a4gdOB1gQ6yFQ+JdUAyUAWMMxae5YjAYqIiIjEIX9BPoluN0nN6mrwVzu3xwsEnp7arGWrcIUmNXA8cTfGHAPMB46qY2pFFMIREREROWJUPnzJGFP35Fq4vYHEvVSJe8TFQleZvwCdgHeBO4CdwFzgTmAKsA/4EejgVIAiIiIi8ajYl0dqIzrKQODmVEB17lHg+I478BNgqbX2QgBjzIlAO2vtQ8HPU4D3gauB/3UsShEREZE448/Pb1RHGQB3+v4dd4msWNhxzwIWVPm8Auhb+cFa+wGBxP3SKMclIiIiEtf8vrxGdZQBcHszACjZVRCOkOQQYiFx38P+m1EBNgDNg7XvlVYCnaMZlIiIiEi88/vyGvXwJahSKqMd94iLhcR9BTDE7L8r4hsCiXy/KnPaEeg2IyIiIiJhYK3F78tvVCtIqNpVpjAcYckhxEKN+3PA88AHxpibga+AdcBDxpgioC1wObDEsQhFRERE4kxpUSEV+/Y1fsc93RNcTzvukeZ44m6tfcEYcwJwM9DLWrvSGHMH8CrwdnDaPuAep2IUERERiTf+/MY/NRUgITERV/N01bhHgeOJO4C19jZjzD8Af/DzdGPMYAI77X7gFWvtCidjFBEREYknxb5A4t7YdpAQqHNXjXvkxUTiDmCt/f6Az58BnzkUjoiIiEhcC+24ZzRuxx0Cde5qBxl5MZO4G2N6EOgc465tjrX2zagFJCIiIhLH/AX5AI2ucQdwpXspKdLNqZHmeOJujDkaeB04+VDTAIs6y4iIiIiEhT8/mLg3sqsMgNvrZe+OHxu9jhya44k78AiBBy59SKA0xu9oNCIiIiJHgMoa95SMFo1ey+3xUrA+p9HryKHFQuJ+GvCetfZ8pwMREREROVL48/NwezNISGp8OujyeClRO8iIi4UHMFUAXzsdhIiIiMiRxO/LC0uZDIDbk0HJrgKstWFZT2oWC4n7DOBcY0ws7P6LiIiIHBGK8/PCcmMqBEplKsrKKPer4jmSYiFZ/j0wB5hvjHkU2AiU1DTRWrsymoGJiIiIxCu/L59mLVuFZS2XJ/D01JLCXSSlpoZlTTlYLCTuycBu4CxgUB1z1VVGREREJAz8vjwyu2eHZS23xwtASdEu0tq0DcuacrBYSNz/QyBp3wx8DuxxNBoRERGROFdeVsbubVtp3q5jWNZzezIAKNlVEJb1pGaxkLifA3wCnGGtLXc6GBEREZF4V7h5IxVlZWR26xGW9Sp33PX01MiKhZtTARYpaRcRERGJDt+6NQC06BqexN1VWSqjxD2iYiFxnwOc4XQQIiIiIkeK/JzwJu7acY+OWEjcbwOOMsa8Zow53RjTzhjjqenV2AsZY943xlhjTG7jwxYRERFpmnw539KsVeuwPDUVwO0N1rgrcY+oWKhxn0ugs8woYOQh5lkaEa8x5pfA0MM9X0RERCRe+NatpUXX8HSUAUhOa45JSFDiHmGxkLh/D2wFvonUBYwxrYD/i9T6IiIiIk1J/rpv6XrBiLCtZ4zBle5RqUyEOZ64W2vPisJlHgHC82gwERERkSasOG8n/vw8MruFb8cdAjeolhSqHWQkxUKNe0QZY84Hfk6g1GaHw+GIiIiIOCrUUebY7mFd1+3JUKlMhDm+4w5gjGkDXAS0JvB0VFM5RKD+PQs411rbpYHrNgcmBj8+BfQAzgxHzCIiIiJNUShxD/OOu9vjpaSwMKxrSnWOJ+7GmD7AAiCdQKJuK4eCP23wfd5hLP8/QCfgB+APwMxGBSsiIiLSxOWvW0Oiy4WnU+ewruv2eCjauiWsa0p1sVAqcy/gAZ4ALge2EEiwfw7cD+wCtgNdG7KoMWYgcH3w403W2jqLrowx440xS4wxS3bsUFWNiIiIxB/fujVkdOlGQmJiWNd1ezNU4x5hsZC4DwYWWGt/Y619DfgAaGetfdVaey9wFpAB3FHfBY0xLuBpAt9vdnDdOllrJ1lr+1lr+7Vq1aqBX0NEREQk9vly1oTtwUtVudK9qnGPsFhI3DOAz6t8XgX0McYYAGvtSmA2MKwBa94CHBd8P8cYM9wYMxxoGTzWLHjs5MaFLiIiItJ0lJeVsWvThogk7m6Pl9KiQqy1dU+Ww+J4jTtQALirfF4PpADdgTXBYznATxuwZs8q7/9dw3grYAYwGRjXgHVFREREmqxdueup2LePzG4R2HH3eLHl5ZTt2YOrefOwry+xseO+FDjfGJMS/LyawM2og6vMORbYF+3AREREROJJqKNMJHbcvV4A1blHUCwk7o8TuPF0mTFmsLU2B/gS+Jsx5jpjzL3ACAIJfr1Ya8dZa82BLwLdawA2BY+NC+9XEREREYld+TkRTNzTA4m7np4aOY4n7tba2cBNQHugXfDwLUAqgaT+bmA3cKcjAYqIiIjECd+6NTRr3Ra3xxv2tV2eyh13Je6REgs17lhrHzPGTCLw8CWstQuNMb2A4YCfQGeY752MUURERKSp861bQ2aYH7xUqfKPgZJdKpWJlJhI3AGstaUHfN4MPBLma5wVzvVEREREmgprLfk539L9kp9FZH23NwOAkiLtuEdK1BN3Y8wZh3uutXZhOGMREREROVIU5+2kpMBHi64R2nFXjXvEObHj/iFwuA0+w/uILxEREZEjxP6OMt0jsv7+GvfCiKwvziTuj3Bw4n450AZ4H/gYyAeaA/2Bi4FNBG5UFREREZHD4Mv5FiBiNe5JqakkJCerxj2Cop64W2tvrvrZGDOewAORLrLWvnPgfGPM6cAcIDk6EYqIiIjEH9/6tSS63aR37BSR9Y0xwaenqlQmUhxvBwncBkyvKWkHsNYuAl4HbohqVCIiIiJxJD/nW1oc252ExMhVHrvSvWoHGUGxkLh3AH6oY04B0DIKsYiIiIjEJd+6tRF58FJVbo9XN6dGUCwk7jnARcaY9JoGjTFtgJHAV1GNSkRERCROlJeWsmvThsgn7t4MSnYpcY+UWEjcHwE6A/ONMSOMMZ2MMS2MMZ2NMWOBBQRuXP0fJ4MUERERaap25a7HlpeT2S2yibvL41Ef9why/AFM1tpnjTHHAH8gUMt+oBLgRmvtm9GNTERERCQ+5Ac7ykR8xz1dpTKR5HjiDmCt/bMxZjJwGXAC0ALwAUuBV4JPURURERGRw+BbtxaIfOLu8njVDjKCYiJxB7DWrgMedDoOERERkXiTn/MtaW3b42pe4y2FYeP2ZlC6u4iK8vKIdq85UsVCjbuIiIiIRJBv3ZqIPXipKnfw6amlu4sifq0jkRJ3ERERkThmrcW3bg0tunaP+LXc6R4A1blHiBJ3ERERkTi2d8ePlOwqoEXXKOy4ezMAVOceIUrcRUREROKYb90agKjsuLuCpTJ6empkKHEXERERiWO+dYFWkFGtcVcv94iIeuJujJlnjPlFtK8rIiIiciTyrVtLUmoq6R2Oivi1XOnacY8kJ3bczyLwpFQRERERibD8nG/J6NIdkxD5tG9/jbsS90hQqYyIiIhIHPOtWxuV+naoUiqjHfeIUOIuIiIiEqdKd++mcPNGMiP8xNRKiS4XiSkplKjGPSKcStytQ9cVEREROWLkzHoVW1HB0T85N2rXdKd71Q4yQpxK3O81xpQ38LXPoVhFREREmqRVLz1LZveetOt/atSu6fZm6ObUCEly6LqFgP4UExEREYmQ/LXfsO2LTzj9vocxxkTtui6PVzXuEeJU4v5Pa+39Dl1bREREJO6teulZEpKT6fmzK6J6XbfHQ0lRYVSveaTQzakiIiIicaa8tJRvXnmBLudeRLNWraN6bdW4R44SdxEREZE4s+H9tyjO28lxY6+O+rXd3gyVykSIEncRERGROPP1S8/SvH1Hjh4yNOrXdnm8agcZIU4k7pOB5Q5cV0RERCTuFW39jtx573Pc6HEkJCZG/fpuj5d9e/dSXlYW9WvHu6gn7tbaq6y1bzb0PGPM6EjEIyIiIhJPVr88Gayl1+hxjlzf7ckA9PTUSHCsVMYY080YM8EYc5Mx5sxDzDvaGPMu8GIUwxMRERFpcmxFBV9PfY6jzjgb79HHOBKDy+MBUC/3CIh6O0hjTALwKDABMFWOLwQutdbmBT8b4HfAPUAaoP/6IiIiIofw3aJ5FG7OZfCfHnAsBrfHC2jHPRKc2HH/DfBroASYAvwfsAo4E3gSwBjTBvgQ+BuBpH060MuBWEVERESajFUvPYs7owXHDhvuWAyuYOJeUqiWkOHmxAOYxgB+oL+1djWAMeYPBJL4nxtjegNvAN2A74DrrbVvOxCniIiISJNRnJ/H+rdncPwvxpOUkuJYHJU17iqVCT8ndty7AzMrk3YAa60FHgzGM4NA0j4Z6K2kXURERKRu374+lfLSUo674hpH43CHdtyVuIebE4m7B9hQw/H1wZ9dgBuD3WeKGnMhY0x7Y8xEY8wmY0ypMWaHMeY1Y8xxjVlXREREJJZYa/n6pWdoc2I/Wh13gqOx7K9xL3Q0jnjkROKeCJQeeNBaWxJ8O99a+3hjL2KMOQZYBlwHdAKSgZbApcDnxpgBjb2GiIiISCzYvGAuO1d/xXFjr3I6FFzplV1lVOMebrH45NSPw7TOg0Cb4PsXgSuAZ4OfmwFPhOk6IiIiIo7ZV1zMvN//howu3ej183FOh0NCUhLJac1VKhMBTtycWpd9jV3AGJMIXBz8uNpae2Xw/UvGmBOAfsBJxph21tptjb2eiIiIiFM++7+/sit3PaNmzHX0ptSq3B6v2kFGgFOJuz3MsfpKJdBmsgOw9oCxXAKJO4AXUOIuIiIiTdKOr1ey9PG/02v0OI46bYjT4YS4PF7tuEeAU4n7LcaYmoqw7KHGrLXH1mdxa+1u4O4DjxtjXMCg4Md9wJZ6xisiIiISUyrKy/ngd9fhzmjB6fc+5HQ41bg9GUrcI8CpxD0j+GroWGPdB7QPvn8/mOCLiIiINDkrn3uCH5Z+xnlPTCE1M8vpcKpxezwU5+10Ooy4E/XE3VrryA2xxpibgTuCH8uBe2qYMx4YD9CpU6foBSciIiLSAEVbv2PxX//I0UOG0mPkaKfDOYjL46Vg4/q6J0qDxGJXmbAzxvwW+GeVQ3dba5ceOM9aO8la289a269Vq1bRC1BERESknqy1zL/jRmxFOT95+D8YY5wO6SBuj5fSIpXKhFvcJ+7GmOuBf1U59Li19kGn4hERERFpjPVvz2DDe28x8PZ78R59jNPh1MjtyaBkl/q4h1vUS2WMMfMO81RrrT27gde6FHisyqH/ADce5vVFREREHOUv8DH/jpto1ftETr7uZqfDqZXL46W8tJR9fn/MtKiMB07cnHpWA+dbwNDANpHGmB7A88FzAZ6w1v6mgdcWERERiQkFG9bx5pXDKc7bwUVTZpCQFIuP4wlwe7wAlBTuUuIeRk78F29Rz3ndCTzd9CSgFPjfBl7nYSAt+D4f+MAYM/yAOfOttSrAEhERkZj23aJ5zL76MowxjHjtPdqe1N/pkA7J7fEAUFJYQFrrNnXMlvpyoqvMIRPl4FNPbwf+DKQAHwO/stZ+U99rGGM6AxdWOZQJvFbD1JOA5fVdV0RERCTaVjw7kQ//+FtadO3BxVNmknFMvR5r4yi3J9DZW09PDa+Y+jcWY8wA4CmgN1AE/M5aO/Ewljqd/SUyIiIiIk1OeVkZH/7xt3z1/JMcM/QCznviRdzpHqfDqhdXlVIZCZ+YSNyNMWnAA8BvgERgFvAba+33h7OetXYKMCV8EYqIiIhET3F+Hm9ffRlbFn9IvxtvZ9BdfyUhMdHpsOqtssZdO+7h5Xjibow5n0C3l07ANuBGa+10Z6MSERERccauTRuZ+fMLKPwul3P/M5meP7vC6ZAaLLTjrpaQYeVY4m6MaQU8AlwWPDQJuN1aW+hUTCIiIiJO+nHll8wcfSHlpSWMfP2/dBh4mtMhHRa3N1DjXlKktC6cHHkAkzHmKuAb4HJgLXCmtfY6Je0iIiJypNr04Rxeu/gsEpOTuWz2wiabtAO40pqDMSqVCTMnHsD0Aft7uS8D/gZkGmMurutca+2bEQxNRERExBHfvPYic266hswevRg+bTbN23VwOqRGMQkJuNI9ujk1zJwolRlS5f3JwMv1OKfyAUxN564MERERkTpYa1n62N/56P476HjaEC6a/Eboxs6mzu3NUI17mDmRuN/nwDVFREREYsa+4mI2/Hc2q1+eTO7cd+k+/DKGPvY8SW6306GFjdvjpaRIO+7h5MQDmJS4i4iIyBGnorycLR/N59vXp7Ju9nRKdxeR1qYdp95xPwNuuROT4MithxHjTvdSWqjbF8PJ8XaQIiIiIvGsOD+PZRP/yeppz7Nn+zZczdPpeuFIsi8dQ8fThjSp/uwN4fJ42P39VqfDiCtK3EVEREQiwO/LZ9nEf/LlpEco27uHLudeSPalY+ky9EKSUlOdDi/i3N4M8tasdjqMuKLEXURERCSM/LsK+PKJf/Hlk/+mtKiQbpf8jIG/v5usHr2cDi2qAqUyqnEPJyXuIiIiImFQsHE9q1+ezPKnH6O0cBddLxzJwN/fTctexzsdmiNcHi8lhbuw1mKMcTqcuKDEXUREROQwFefnsXbWq3z72lS2ffExAF2GXcypt99Lq959HI7OWW6PF1teTtmePbiaN3c6nLigxF1ERESknir27WPXpg38uHIZa2e+ysY571BRVkZW9nEM/vP/kD1qNOkdjnI6zJjg9gb60ZcW7VLiHiZK3EVERESCbEUF/gIfxfk7Kd65g70//kB+zrfkrVlN/prV+Natoby0FIBmrdty4rU3kH3pWFodf6LKQQ7gSg8k7iW7Cpr8k2BjhRJ3EREROWKVFBWy4qnHWDPzFfbu2I4/Pw9bUXHQPE+nzmT16MXRZ59HVveeZGUfR6veJ5KQpFSqNp5OnQHIX/sNWdnHORtMnND/bSIiInLEKSncxfKnHmXZE/+ipMBHx9OG0H7AIFKzWpGa2ZLUrKzA+6xWZHTpplKPw9D6hJNJbpbGlsUL6HbxpU6HExeUuIuIiMgRw7+rgOWTHuHLJ/9Nya4Cupx7Iafc9mfanNjP6dDiTmJyMu0GDGLLJwudDiVuKHEXERGRuFFSuIu8Nasp2rIZvy+f4vyd+PPz8PvyKM7PY9uSTykt3EWXYRcz8Hd/pnWfk50OOa51HHQmHz/4J4rzdpKa1dLpcJo8Je4iIiIS8yrKyyndXUTZ7iJKg6+y3UUUbs4lb83q0M2ju7dtPehcV7qHlMwsUltkccw559P3+t/R+oSTHPgWR56Og84AYOuni+h6wQiHo2n6lLiLiIhIzCn6fgubP5zDpg/n8N2ieRTv3FHr3KRmzcjs1pOjTh9CZo/jyOrRE+/RXUjJbElKi0wSk5P/n737jo/jrvM//vqod1mWe+/G6T2QkATScAhwQELgSAjwo4VyhHYHhCPAATk4CAfkwtHhINQAoZMCSUhPSLUTO3G3Yzu2LEtyUdfu5/fHzErjza60WpVded/Px2MeszPzmdmv5quRPjv7ne93HEsuUdOPP5mSykq233eXEvdRoMRdRERExk2sp4e9a5+kr7sLj8WIx2N4Xx/xWIy+jnZ2PHAPW++8jZZn1gBBl4sLzn4Z9QuXUFZTS1l1DaU1NZRV11JaU0vtrDnUzVuAFRXl+CeTVIrLyph50ovYfu/fc12Uw4ISdxERERkz7k7rxnVsu+M2tt55K9vvuZPejva08cXl5cx+0Rkc+ca3Mv8l59G44ij1jz7BzTn9LO7/4qfpamulYlJDroszoSlxFxERkX4de5o4uGsn3W0tdLbspautha6WFrraWujZv4+e9oP0th+kt709aGfefhCPxSirrQumujrKa+soq62nr7ODbXf9jQPbtwEwaeESVrz+cuacdiZldfUUFRdjxcVYUTFFJSUUlZYy5QVHUVJZmeOzIKNp9ovOBHd2PHA3krMvIwAAIABJREFUi1e+KtfFmdCUuIuIiEwg7j7kHeiegwdo3/0c7bueo3PvHorKyiitqqa0uiZoalJdQ3F5Bfuf3cLeNatpfvop9q59kua1q9O2JS+pqqK8tp6ymlpKw2NUTZ1G6YJFFBUX03NgPz0HgodFg9f7AZh92pmcfOVHmf+S86lfsGjUz4fkvxknnEJxeTk77rtLifsIKXEXERHJkXhfH31dXcT7eon19BDv7Qnmfb10teylbctG9m3ZxL4tG2kL553NeyiprOxPxEurqimpqqa4pJSO5ibadz9Hb/vBYZWjpKqKxuVHsuj8V9C44mjq5s6jYtJkKiY3UtEwmYpJkympqBijsyCHu5KKCmac+EK236f+3EdKibuIiMgwuDudzXuIx/qwomKsqKh/AujYs5sDO57l4I5nObBze/B653a69rUFTUw62sOmJgeJdXcP/YZm1M6eS/2CxSxe+Sqqps0g1t1FT/tB+jragyYr7QeJ9/Yy9ejjWHjey6mePpPq6TOonj6TyinT8L6+SBOXsAwdHdTMms2UFUdTP3+hHu6UMTXntDN56Cufp3v/Psrr6nNdnAlLibuIiEx4Ho/TubeZg7t20r5rJyUVFdQvXELtrDkpE1KPx2nduI6mJx5h9+OPcHDXTmpnzaF29jxq58ylds48amfPo6y2jpZn1rDnqSdofmoVe558gj1PPUF3W2vGZauaNoPaWXOoaJhM7aw5/XfJS2vCu+XlFRSXlVNUVkpxaRnFZWUUlZZRVltH/YJF1M1dQEl5+WieLpFxN+e0M3nwy59l5wP3sPD8C3NdnAlLibuIiAwq3tdH68Z1lFZVUztn3qj28NHV2kLT6sdoWvUYLevW4rG+1GWIxfBw6u8+sK+PztYW2nftpKNpF/G+5+9bXF5O/YLFTFq0hEkLl4A7ux9/hD2rH6Pn4AEASiorqZkxm823/pG+zs60ZS2prKRxxdEsfdXFNC4/kpKKcuKxGMQdj8eDyeNUNk6hdtZcamfPpXrmbCXdIsCME19IcVkZ2++/S4n7CChxFxGZIOJ9fXTv3xfcoa2ooKi4OKvj9HV30/7cjuAhx+KwN4+iYqykBDOjbfMGmlY/xp7Vj7Nn9eM0r13d36SjatoMZp50KjNOfCEzTzyV6cedRGl1NRA0IYl1d9PX2RE2xWg/ZKTL3vaD9Bw8QMeeJvY8+ThNqx7r720EoHr6TIpTtaN2D3odCXsgKSopCZqoFBdT0dDA5KUvoGbGLKpnzKJmxkyqps+kr7ODts0baNu0oX++9fZbwIypRx3HikvexLTjTmT6cScxeekLKCopwd3pam3hwPZtHNixjQPbn6WrrYWGJcuZevRxTFq4JOtzLlLoSquqmH78KWrnPkLm7rkuQ1466aST/OGHH851MUTGTTwWo2XdWurnLexPxAqFx+NBkhn2hNEdzj0ep37+IurmLRh05MV4X1+Y7D1LSWUl5fWTKK9voLx+Usr9+hPc7i68ry+4UxuL4R7ctY3HYnTs2U3rhmdoXf8MLRuepnX9M7Rt3kC8t7f/OEUlJRSXV1BSWUlJZRWVk6dQ2TiFqilTqWycSsXkKZTV1NK+ayf7n93C/me3sn/bFtp3P5fReSmf1MC0o49n6tHHMfWoY+k5cIBdjzzAcw8/SNum9QBYcTGVkxvp7eigr7MDj8eHPrAZDYuWMvWY45l29PFMO+Z4ph59PJWTGzMqV7YSd8WLSnTPSiQX7rvmk/zj61/k3Rv2UlZTm+vi5JuMvspU4p6GEncZS+5Ox56moFu1qqqclaNzbzNb77iFzX/9C1tvv4Wu1haKy8uZc9pZLDj35Sw89wImLVqSct++zk7279hGZ/OeIOE9uL+/O7ieA/spKi2lbt4C6ucvon7+IqqmTT+kiUU8FqOjaRf7t2/j4I5n6WzZS2l1NWW1iT6gw6mmjuKKirAdcNkh7ZV7Dh6kbeM6Wjeto3Xj+vD1BmLd3ZRUVFAc7pN47X19QTkPHuxP0nsO7Ken/SAM8rfQioupm7uASQsXM2nRUqqnzeDAc9vZt2Uj+7Zs5sD2rSmbaUDYhV7dJDweJ9bdRV93F7GurozrqKikhPqFS5i8dDkNi5dTNW0G8d6e/uP0dXcR6+6mt/0AnXub6dzbTMfePXTu3UNfR0f/MWrnzKN2znzq5s6nbt4CamfNoaikhHhfX9D8JNZHvC9oilI3bz5Tjz6e2tlz0zaL6WzZy65HHmTXIw/S0dwUtNWurKKksipov11ZSUlVVfA7XlMbjHhZU0tZdS3l9ZPUT7dIAdp6523c9LqVvPrnf2LBOStzXZx8o8R9JJS4jx135+DO7TSvWc3ep5+iee1qOpp2UdHQSOXkKVQ0TqGqcWrQDdmkhrQPlvV2dtDX0R72rNDRPy+uqDg08QunkvKKga/Zi4spKi4Jvh6PxyNdsfUS6w3m8b7eMKEJ29bGg9e4h8cooaj/eMHAIWU1tZGEs7b/a/X23bvY/cTD7H78EZoef5jdTzxKR9MuAGpmzmbSoqVBG9xwXlJRGen9oaP/dV9Pd//Q4N7Xd0hb33hfL7HeXuK9iZ+jFzxOSVV1f7/Lif6b4319bLvrb+x65EFwp3LKVBacvZI5p51J89NPseWvf6F1wzMATFq0lAXnrKSotJQDz25l//ZtHNi+lY49Tekr2ex5iXBJZSX18xdRVlvPwee2075rZ9pkdzBFpaVBXZaUPO8Bwdo585i0aCmllZX0dXcT6+4m1tNNX1cXse4uikpK+usmUVeJpLK8rp6yuuCDQuLDg8fj7Nu66ZDmFm2b1tNz8AAVDZOpX7CY+vkLgzbUCxZRM3su8Z4euva10r1vH9372uje10r3/n1YcXF/E5eS8uCDRElFBUUlpUGPJMVBryRFRcVQVERFw2QmL30BdfMWDnq3fzC9HR30HNhP5ZSpauIhIjnX297O/y6ZzAnv+TAv/uQ1uS5OvlHiPhK5StwfvPbzlFRUUFpTM5Bc1NT1f6XU2xl8HZ2Y93V2BslJou/f/nkvjvcnaqXVAwNvFJdXEO8J79h1d4VJTTd9XZ0D3ZSF855wdLyikpL0SUeijWzY7jRxF6+vqzMoY1dn8DV6Vyftu56jee1qevbv6/+Za2bNoWbWHLpbW+hsaaartSXr81dUWnpIU4JcK62qpri8fOBnMmPyshVMP/ZEph59PH0d7bRuWh8mhuvTDnySUFxeTlFxCVZSMtDeN/wAUlRaSlFpGcWlpRSVlFJUVoaZHdL1XOLDDWZMP/5kFp57AQvOvYDpx574vA9IbZs3suVvN7Plr3/m2XvvBKBu9jxq586jbs586uYuoHbuPKqmTn/eXfLSqmpiPT3sf3ZL0Af11s3s27qJ/ds2071vHzWzZge9d8yeEzzEN2celY1T6evs6G+m0nNwPz3799F9YD+x7p6Bu9Xd3cS6u4j19FAzczYNi5fSsGQ59QsWj8u3F+5OX2dnTr8pERGZqH7x8hcD8Po/35PjkuQdJe4jkYvE3eNxvj6zPLM2ooMxo7isDCCzPoKTHDKwR3UNJZVVeCwWSfIHvu4Pvl4PvmpPVY7Syqqg/W1FMFU2TqXxiKOYsuJopqw4isYVR1ExqeGQ3eJ9fXS1ttC5dw9dabpcM7PgK/nqgQFISquqgw8NsViQpEbbK+/fR6y3Z+DudOJOdW8vFBUF3a+VlB4yt8QDe/0JcjAsN2bB3fe+voG737E+Yj099Bw8ECadB/rfv7ezg4bFy5l+3IlMPeo4ympq0p777v37aNu8gVhPT/ChKzrASmXlqPSzHI8FZR9OTxex3t7gg9ko9iYiIiKF557PXsWj37iWd29oKbjnqYagxN3MyoGrgMuBWcAO4EfANe7eM9i+ubrj3v+QXNgLQ8+B/f09MvQnwhXBg2glVVWUVlZRXFZOcXk5xaVlFJWVHfKVeLyvL+nu+UH6OjspLi8/5M55cVk5xRUVQfKbxVfq7h42OQmSWSspCZJfJXoiIiIS2vK3m/ntGy7ktb+6hXlnnZvr4uSTjBKmw/3R+huAiyPLC4FPASuA1+ekREOwoqL+JjKjoaikhPK6+jEfpczM+u9Koz6LRUREJIVZp56OFRez/d6/K3HPwmE7vrGZnc9A0n43cGk4B7jEzM7LScFEREREClRZTS3TjjlB/bln6bBN3AmaxwA48EZ3/ylwWbgc3S4iIiIi42TO6Wex+7GHBh2pWFI7nBP308P5RnffDuDu24ANSdtFREREZJzMedGZxHp6uPOqK9nxwD1BV8uSkcMycTezMmB+uLgjaXNieX4YJyIiIiLjZO6Z57D0n17H2l/+mBtfeRbfOWo2t135djbd8gfdhR/C4fpwaj0DT+ceSNrWHs6LgDqgebwKJSIiIlLoSioquPC7P6f7wH62/u1mNt78e9b/8Tc89dMfUFJZSdXU6cHAdEXF4QB14XwUukQeyvlf+x7Tjj1hzN8nW4dr4l4ReZ38/Ut0+ZAxt83sncA7AebNmzc2JRMRERERymvrWPbqS1j26kuI9fSw/b6/s+Wvf6GrtQX3eP/o5cTjxOOx543IPRaKKyqGDsqhw7IfdzObBuwOF//o7q+MbPsjcGG4ONXdU95xz1U/7iIiIiJScDLqx/2wbOMO7Geg95jkcckTw3TFwzgRERERkbx3WCbu7t4F7AwXZyZtnhXOtw41eqqIiIiISL44LBP30IPhfLmZzQYws5nAknD9fTkplYiIiIhIFg7nxP2n4bwI+ImZvRH4OQM/8w05KZWIiIiISBYO115lcPdfm9kfgFcCZ4VTwm/d/ebclExEREREZPgO5zvuABcDVwObgB5gC/A54A05LJOIiIiIyLAdtnfcAcKHTz8bTiIiIiIiE9bhfsddREREROSwoMRdRERERGQCUOIuIiIiIjIBmLsPHVWAzGwPsHWM32YK0DzG7yH5QXVdOFTXhUN1XThU14Uhl/Xc7O4rhwpS4p5DZvawu5+U63LI2FNdFw7VdeFQXRcO1XVhmAj1rKYyIiIiIiITgBJ3EREREZEJQIl7bn071wWQcaO6Lhyq68Khui4cquvCkPf1rDbuIiIiIiITgO64i4iIiIhMAErcRUREREQmACXu48zMys3sM2a22cy6zWyTmX3azMpyXbZCZWZXmJmnmdoicWZmV5rZ2rDutpvZV82sLs1xLzOzx8ysy8yazOwHZjYjTezLzOw+M+sws1Yz+7WZLUkTe7KZ3WZmB8LpVjPL6+6rcsnMvptcl0nbc15PZrbMzH4THrPdzO41s/PTxM40sx+a2R4z6zSzR83sjZmej8PZYHVtZi8Y5Dp3MzsuKV51nWfMbJaZ/a+ZbTWznvC83GhmR6aI1XU9gWVa1wV5Xbu7pnGcgBsBTzH9ItdlK9QJ+EaaOnGgLRL35TQx9wLFScd8b5rY9UBdUuyrgFiK2CZgXlLsyUBHith24MRcn8t8m4BzgN7kusynegLmA3tSxMaAVyTF1gEb05T53bk+33le15cMcp07cJzqOn8nYCGwK835aAdOicTqup7A0zDruuCu65xXUCFNwPmRCroLeGM4T6w7L9dlLMQJuCc8/48Ar06aLgxjVkQu+NXAZcBvInX3jsjxGoED4fpngbcSPKmeiP18JLYE2BaubwPeDXw+EvuTpLLeF67vAf4V+AgDycq9uT6X+TQBF4V/ZBPnsi1pe17UE/DT6HuGx94XLm8FSiKx10Rivx2W+dlweT8wOdfnPR/rOoz5XLitO8V1/mqgXnWdvxPws8j5+DFwKfC9yLpHwzhd1xN8yrSuw9iCu65zXkGFNAE3hJUTB+aE6+aFyw78ONdlLLQJsMgFdv0gcZ+LXGynheuqwgvNgbsjsW+PxL4x8j7rw3XbIrHnRmKviqz/a7iuE6gJ1y2OXvCR2O9G1i/O9TnN9QQsAX4VOSfpEvec1xNQE+7rwK2R2H+PxJ4TWZ/4A/9MZN1lkdi35fr852Ndh7F/CLc9NcQxVdd5NgHFDHwweypp2z8i52SmruuJPQ2nrsN1BXddq437+Do9nG909+0A7r4N2JC0XcbPAoKvrgDWm9kpZnZ+iraQibrpBh4EcPcO4IFw/SlmVpIUC/D3MNaBO8J1c81sXrrY0O3hvAI4IXz94iFik49XqF5BcAcW4CZgd5q4fKinE8N9h4w1s/nAnHDdXUMct1BkWtcAx4Tzp8N2sSvNbHmKONV1/qkErgW+D/woaduWyOt6dF1PdMOpayjA61qJ+zix4OHT+eHijqTNieX5podUx9uxkdfXECTltwA7zOzrkWR8WTjf5e6xyD6JuovWbyI2DjyXIhZgaVJs8vaRxha63cD73f21QFeamHyop9GIfY7gZ4jGFpIh69rMJhF8uwlwAbAW+AvBP/u/mdm0SLjqOs+4+0F3v9rd3+buX0ysD/9fnhYu9gHb0XU9oQ2nrgv1ulbiPn7qCb6ug6D9XVR7OC9i4O6vjI9jIq8rI6+LgH8BrguXJ4XzdHUH0JAU2+Hu8Qxjk4890thC9htgvrtfN0RcPtTTiGPDu4mdSbGFItO6TnedA5wN3Bq5aaK6njg+A8wKX9/i7gfRdX24SlXXBXldK3EfPxWR17GkbdHl5F8+GVuVBG3cmwkeZKki6J2iOdz+LjNbxkD9ZVJ32cQmbx9pbMFy923u3p1BaD7U02jERpcLqv6HUdeNBL1UxIB/I/gHuhx4ONx+LHB5+Fp1PQGY2QeAj4WLMeBT4Wtd14eZQeq6IK9rJe7jJ/rPpThpW3S5Exk37v5xd58EzHX337l7p7vfTtD1IwTfkqxkoP4yqbtsYpO3jzRWhpYP9TQasdFl1X8K7n6Tu88kePjsS+6+z93XAe+LhL08nKuu85yZXQn8d2TV1e7+SPha1/VhZLC6LtTrWon7+En0PgLBXd2o6nAeD+NknLl7ctvY1ZHXswnuykP6ugNoDeeJ2Eozswxjk4890lgZWj7U04hjw7JXJsVKChlc56C6zmtm9h7gq5FV17v7NZFlXdeHiQzqGii861qJ+zgJf7F2hoszkzYn2m1tdfee8StVYbPAsWZ2oZmtTNpcE3l9kGDABIAZSf8MEnXXS9BfK5HYYmBqilgIuiaLxsKhvxcjjZWh5UM9jUbsDAb+lqv+UzCzRWZ2rpm9KWlT8nUOquu8ZWYXA/8TWfUNgmeRonRdHwYyqetCva6VuI+vB8P5cjObDcHQtwR9EUPQ4b+Mk/CBkN8DfwRuMLP6yOZoIn8/A3VXBZwKYGblBKOrATzk7n3h6wcj+54deZ3oYmo7QV+ug8WeEc67gEczjAX9Dg1HPtTTIwQ9JGQSu4lgdD6Alw4RK4f6FHAb8CMzi3b1lnydg+o6L4Vd/P2QgU4evunu7w3/jkfpup7ghlHXhXldj2VH+pqeNwDARQx0sn8nwcipf4+sW5nrMhbaBHwpcv4fIHiQ5WsMDIq1muAD7omRuFUEAyb8OrLuisgxZzMwIMNW4C3AtyKxX4jEVhJ0Z+cEX5FdwaGDPf08qbyrwvU9wIcJRnPrSZQ/1+czHyeCvn+d5w/AlBf1xMAAIg58Njx2W7j8LIeOunddJPabBKPuJUYDPAg05vp852ldXxg5bzuB94R1khh8rZNwyHPVdX5OBDdZEudjL3AxKUbJ1HU98adh1HVBXtc5r6BCm5J+IaPTTbkuWyFO4cX8RJo6aQGOjMRelybuoeiFGcZ+OE3sJmBSUuzrCJ4oT47dCyxMij2dgX9K0akTODXX5zMfJ9Ikc/lSTwR9++5NERsDXpMUOzUsW6oyX5nrc53raYi6/k6a8xYDLlVd5+9EMFBePM25iE7HhfG6rifolEVdF9x1nfNKKrSJYKCeTxK0deoGNhN8QivPddkKdSLoQurasC56CLqX+hFB/9DROCN4Wn1NWHc7CZL5+jTHfQPB12hdBF+P/QiYlSb2POBuoIPgk/pNwPI0sScCNxM8yHyAYMhmJe3p63cLaZK5fKknguZyNxLcCeog+Lr0gjSx04HvAU1hmR8HLs/1ec6HabC6ZmBshsfD89ZGMNjamarr/J6ANzF0ItefzIX76LqegNNw67oQr2sLDyIiIiIiInlMD6eKiIiIiEwAStxFRERERCYAJe4iIiIiIhOAEncRERERkQlAibuIiIiIyASgxF1EREREZAJQ4i4iIiIiMgEocRcRERERmQCUuIuIiIiITABK3EVEREREJoCSXBcgX61cudJvvvnmXBdDRERERA5/lkmQ7rin0dzcnOsiiIiIiIj0U+IuIiIiIjIBKHEXERERkSG1blzH1jtuzXUxCtqQbdzN7PJsD+7uP8p2XxERERHJH/d/4VNsuuUPvHtDC8VlZbkuTkHK5OHUHwIeWbYUywnR9QBK3EVEREQOA3ueWkVfZye7H/sHs049PdfFKUiZJO4fTFouBj4C1AP/B9wHtAA1wMnA24Bm4KrRK6aIiIiI5EpfVxdtm9YDsP3+u5W458iQibu7fy26bGZXAXXAGe7+aFL4L83se8ADwEnAr0aroCIiIiKSGy3r1uKxGAA77vs7fOBjOS5RYcrm4dQrgF+nSNoBcPengRuBN4+kYCIiIiKSH5rXrgZg7hlns/Oh+4j39eW4RIUpm8R9MtCeQVxVFscWERERkTzTvGY1xeXlHHnpW+ltP0jT6sdyXaSClE3ivhp4jZnNTLXRzJYCFwP/GEnBRERERCQ/NK9ZzeRlRzD3xWcDsOO+u3JcosKUTeL+BWA6cL+ZfdDMzjSz483sJWb2CeAegrvtnx7FcoqIiIhIjuxd+yRTVhxF9fQZNCxexvb7/p7rIhWkTHqVOYS7/87M3gF8GbiW53cNuQd4nbvfk02BzKwYeBmwODz2OuA2d0/ualJERERExljn3mbadz/HlCOPAWD2aWey7nc3Eo/FKCouznHpCsuwE3cAd/+emf0KuBA4BmgAWoFHgD+7eyZt4J/HzI4Gfg0sSdq0xsxe4+7rszmuiIiIiGQn8WDqlBVHATD7RWfw5I+/S/Oa1Uw7+rhcFq3gZJW4A7j7PuCn4TRavgfsAC4iuNNeAbyIYBCobwMvHcX3EhEREZEhNK95EoApK44GYM5pZwFBO3cl7uMrmzbuAJjZSjP7hZk9Y2ZN4bpLzexqM0vbo0x4Vz3VegNOBK5x99Xu3u3u+9z9ZoIPBy/MtqwiIiIikp3mNauobJxC1fQZANTOnkvd/IVsv1/t3MdbVom7mX0T+BPwOmAR0BhuOongodS/mVlNmt0fN7Ofhr3P9AvbsG8F3mpmdZH3Wgy8AlAzGREREZFx1rz2SRpXHE1wjzUw+4VnsOP+u9EjiONr2Im7mb0LeCdBW/SlwOcjm/8D+D5wKvDhNId4OUEb9qfM7HtmNi+y7X3Aq4BmM9sZ3slfB8wA3jvcsoqIiIhI9jweZ+/TT/a3b0+Yc9qZdLXspeWZNTkqWWHKduTUVe5+ibtvJNKrjLu3uvvbCfpwvyTVzu5+i7ufArweOBlYZ2bXm9ksd/8zMB/4CPCrcPoQsMzd786irCIiIiKSpX1bNtHX0cGUIw5t6Zxo575d/bmPq2wS9+XAzUPE3AksGCzA3W8CjgXeBpwHbDCzL4fbvu7u73f397j719x9VxblxMyuMDNPM7Vlc0wRERGRQtHfo0xS4l43fyE1s+aoP/dxlk3i3gFMGyJmVhg3KA/8BFgBXEnQZn6zmX3OzCZlUbZkx4zCMUREREQKUvPaJ8GMxuVHHrLezJj9IrVzH2/ZJO73AK81s7mpNoYPnb4GuHewg5hZg5mdED586u7+HYI2858kuAu/2cw+YWbVWZQxIZG4PxqWKTpdOoLjioiIiBz2mp9axaQFiymtfn46Nue0M+lo2kXbJvUfMl6ySdz/AygHHjKzDxE0ncHMzjKzjwD3AaXAf6ba2czqzewXQDNBW/h1wDYze5W797j71whGTf0vgvbtm83sw2ZWMZxCht1LJr7XecDdf5s0/Wm4P7iIiIhIIdm79kkajzgq5bbZ/e3c1VxmvAw7cXf3R4HXAsXAl4E3AAbcTpBslwKXuvuDaQ7xNeCVwFUEPcz8M7AB+Hmii0h373D3/wQWAt8CPhXGDMcCINGt5HozO8XMzjezGcM8joiIiEjB6e3ooG3zBqYckbrlccPiZVRNnc6O+9V/yHjJauRUd/+Lmc0H/gk4AZgEHARWATeFo6qm8wrge+7+xcQKM7sdaALOJ9Jfu7vvBz5pZl8HPjrMYh4beX0NUBm+jpvZ9cCH3L1vmMcUERERKQgt69bg8fjzuoJM6G/nft9duPsh/bzL2Bh24m5mZwJb3H0b8PNwSo5ZAZzi7v+X4hAxgqYwUYnBmPanek9330PQReRwRD8eVkZeFwH/QvDNwLujO5jZOwn6qGfevGj38iIiIiKFpXlN6h5louacdhbrf/8r9m/bQv38heNVtIKVTRv3O4A3DxHzFuD6NNuuA1aa2ZPhCKq/B/4KbARuyqI86VQC+wja0r8aqALOCZcB3mVmy6I7uPu33f0kdz9p6tSpo1gUERERkYmlee2TlFRWUr8g+X7rgNmnnQmonft4GfKOu5ldDJwWXQVcYGYNaXYpIxhcqT3VRnf/nJmtAf4fQXOWVuB/gS+6+8FhlH1Q7v5x4ONmVuHuXeHq28O+4r8Q/hwrCR6OFREREZGI5jWraVx+JEXFxWljGpcfQUXDZHbcfzdH/vNbxq9wBSqTpjJPADcQJOQQjJT6wnAazCfSbXD33wC/yaSAIxVJ2hNWR17PHo8yiIiIiEw0e9euZsG5Lx80xoqKwv7cNYLqeBgycXf39WZ2CtDAQO8xPwRStV93oBfYEbaBz4mwK8hjgDlAzN2jI73WRF6P2h1+ERERkcNFe9NuOvY0pe1RJmrWqS9m459/R+feZiobp4xD6QpXRg+nuvuqxGsz+wxwh7vn7UefZHTWAAAgAElEQVQrd/ew7fw8YK+ZLY70dLMyEnr/+JdOREREJL/tXZt4MDV1jzJRk5e+AIDWjeuUuI+xbPpx/0w+J+0RvwznjcAtZna5mX2N4MFZgCcJvj0QERERkYjmtU8CZHTHvWHJcgBaNzwzpmWSLPtxN7PTgXcAywhGUU3Vcae7+4kjKNtIXU3QL/wxwKnhlNAKvMHd47komIiIiEg+a16zmqqp06maMnQve3Vz51NUWkrrxvVDxsrIZNOP+2sJ7mYPdbfesyrRKHH3TjM7C/gkwUivs4EW4Fbgk+6+NZflExEREclXzWtWZ9RMBqCopIRJC5fQulF33MdaNv24f5zgAdRLgQZ3L0ozpe87aJy4e5u7f9jdF7p7mbvPcPfLlbSLiIiIpBaPxdj7zFMZNZNJmLR4KW264z7msmkqcxRwg7v/bDQLYmYVwHyCpjcpRR+SFREREZHRt2/LRmJdXTSuyOyOOwTt3LfefgvxWGzQft9lZLJJ3NtIM7hSNsysEfg2weimQ9FvgoiIiMgYan4quE865YijM96nYdFSYt3dHNi+jfr5C8eqaAUvm8T9t8CrzOxjKQY3ysZXgdcAG4BHgNE4poiIiIhkoXntk1hREY3Ljsh4n2jPMkrcx042ifvHgZOAO8zsOmA90J0qMMOmLecD9wFnqpcXERERkdxqXruaSYuWUlJZmfE+DYuXAdC6cT0Lzlk5RLRkK5vEvYWgxxgDThkiNpOmLeXAvUraRURERHJvz+rHmX78ycPap3LKVMrq6tWX+xjLJnH/EaPb1eMtwBmjeDwRERERyUJny172b9vCMW+5Ylj7mRkNS5bTtmndGJVMIIvE3d3fMspl+CBwr5n9DPhvYDPpm97sH+X3FhEREZFQ06pHAZh27AnD3rdh8VJ23HfXaBdJIrLpxz0jZlaVYWgr8A/gEuB+YFe4LnlqGYNiioiIiEioadVjAEw7+vhh79uweDkHdjxLb0fHaBdLQtk0lcHMjiFItKcRtGO3xCagFGgEXgzUZnC4/wYuAjqBNYxiV5MiIiIikrmmJx6hbv5CKhomD3vfhsVLAWjbtJ6pRx072kUTskjczewlBO3SSwgS9cSDqgmJ9u9PZHjIi4CngDPcfd9wyyMiIiIio6Np1aNMP2b4zWQg0iXkxnVK3MdINk1lriJI2j8GvJCgO8ifhK//H7ANaAYuyPB4FcBflLSLiIiI5E5XWyv7tmxiapaJ+6SFSwBo3aAHVMdKNon7ScCf3f1L7v4QcAdwpLs/5O4/BM4CKoFPZni8ewF9LBMRERHJoT2rg/bt07N4MBWgtLqa2tlz1bPMGMomca8GnowsrwGOMLMSAHffBvwOeEmGx/sI8CIzu9bM5mZRHhEREREZod1PhD3KZHnHHWDS4mW06I77mMkmcd/LoQ+dbiR4IPUFkXXPAvMyPN5XgD3AB4AtZtZlZi0ppr1ZlFVEREREMtC06lFq58yjsnFK1sdoWLyMto3rcB/NIX8kIZvE/QHg1WaWqNUnCR5OPTcScyyZ9w6zjKBnmm3h9BywL8WkPtxFRERExkjTE49m1X97VMPiZXTva6Ozec8olUqisukO8lrgTuApM3uzu99sZn8HPm9mM4AZwErgN5kczN0XZFEGERERERkl3Qf207ZpPSsuedOIjtOwZBkQ9CxTNXXaaBRNIoZ9x93d7wUuJhgQqTxc/X7gAPBvwOXAVuCjo1RGERERERlDexIDL43CHXeA1g3PjLhM8nxZDcDk7r8DfmdmFi6vNrMlwNlAF3CPuw9r2CwzuxR4G0EzmyqCtvRPAf/n7j/NppwiIiIiMrSmVSN/MBWgdu58isvKaN20fjSKJUmyStwTPPLkgbsfBH4/3GOEyf9PgNcTtJVvI3jgtQE4DzjXzF7u7peNpKwiIiIiklrTqkepmTmb6mnTR3ScouJi6hcu0R33MTJk4m5m78/24O7+9QzCrgDeANwOvN/d10TeewlwPfDPZvY3d/9BtmURERERkdR2j8KDqQkNS5bRul6J+1jI5I77VwEnuBs+HA5kkri/neAO+yvdvfOQA7hvMLPXAKuAdwFK3EVERERGUc/BA7RueIblr3n9qByvYdEyNt/6J+J9fRSVjKhxhyTJ5Gy+dYzLsAL4XnLSnuDuHWb2F+DNY1wOERERkYKz58knwH1U77jHe3vZ/+xWJi1cPCrHlMCQibu7/98Yl6EXqBkipgaIj3E5RERERApOUzhi6vRjThyV40V7llHiPrqy/v7CzOqACwl6gakHmgkGZ7rN3XuGcaiHgH8ys/nuvjXF+ywAXg38I9uyioiIiEhqTaseoWraDKpnzByV4zUsWQ5A68b1LDxvVA4poawSdzN7M/AVYBKHtn13YLuZvcPdb83wcP8F3ALcYWafAe4mGCl1FnAGcBVQB3wpm7KKiIiISHpNqx5j+ig1kwGomNxI+aQGWjfqAdXRNuzE3cwuBL5PkFx/nuCO+S6CJP404F8I+ng/092HvEvu7reZ2ZUEI7J+P/ntgD7gA8P4ICAiIiIiGehtb6dl3VqWXPiaUTummdGwZDltG9WX+2jL5o77JwhGTT3V3Tclbfurmf0ceBD4NEFTmiG5+3Vm9kfgMoKmN3UEI7E+Dtzg7puzKKeIiIiIDGLPU0/g8fioPZia0LBoKc/ec8eoHlOyS9yPJkimk5N2ANx9nZn9Grgok4OZ2ZnAljA5/2yamBXAKePwoKyIiIhIwWha9RgA044dnQdTExqWLGftL39Mb3s7pdXVo3rsQlaUxT77gOIhYhzoyvB4dzB0V49vIRiISURERERGSdMTj1A5ZSo1M2eP6nEbFi8FoHWTmsuMpmzuuH8b+KiZfSdVG3YzWw68Dvhaqp3N7GKCtvD9q4ALzKwhzfuVAa8H2rMoq4iIiIik0bTqUaYdcwJmwx1nc3D9PctseIZpRx83qscuZNkk7g8Aa4F7zexnwF3ADqASOBl4J9ADtJrZ+6M7uvvXgSeAGwgScgjuzr8wnAbziSzKKiIiIiIp9HV2sveZNSx62StH/diTFi4BM9p0x31UZZO43xx5/aZw8nA5+nHtWp7fVeTX3X29mZ0CNITbbwd+CKRqv+4EAzTtcPdtWZRVRERERFLYs2YVHosx7ZjjR/3YJZWV1M6ZR+sGdQk5mrJJ3N860jd191WJ12Hf7Xe4+10jPa6IiIiIZCYxYupoP5ia0LB4Ka3qEnJUDTtxH+2eXdz9M5nEmdlCdQspIiIiMjqaVj1KRcNkaufMG5PjNywOepZx91FvQ1+osho5dbSZ2cuBNwLTCHqsSdSuAaVAI7CMoXuzEREREZEhuDu7H/sH044d/QdTExoWL6PnwH469jRRPW36mLxHoclm5NQi4L0EifYCoDxNqLt7YwbHey1wI4e2h0/WDvxueCUVERERkVSevvEGmtes5iVveseYvcf0408CYN1Nv+D4d71/iGjJRDb9uH8S+Cpwari8L820P8PjfQjoAy4BZgCPAd8JX58NPELwkOpHsyiriIiIiEQc2LmdOz5+JbNOOZ1j3nrFmL3PzJNeyNwzz+Ghr15Dz8EDY/Y+hSSbxP3NwDZgobvPdPeF6aYMj3c08Ft3/5W7NwH3AC929yZ3vxN4GdCNuoMUERERGRF3568feAfxvl7Ov+77FBWPbSvk0//983Q27+HRb351TN+nUGSTuE8DfunuW0epDBXAhsjy08AyMysHcPcW4LcM3c+7iIiIiAziyR9/h6133MoZn/ovJi1aMubvN+P4k1ly4Wt49Ppr6WjeM+bvd7jLJnF/FBjNmt4NTI0sbyQo15GRdc3AnFF8TxEREZGCsm/rZu66+l+Ze+Y5HPOWd43b+5521Wfp7WjnH1/7wri95+Eqm8T948DLzewKG53HkP8OXGRmy8LlJ8L5P0ViTgdaRuG9RERERAqOx+Pc+v63YWac97XvYkXZpIDZmbxsBUe84c2s+v432L9d42mOxLBrzd3vBb4FXA/sM7OnzOzRFNMjGR7yC0AlsNrMLnb33cAfgKvM7BdmdidB4n7bcMtqZuVm9hkz22xm3Wa2ycw+bWZlwz2WiIiIyET1+HevZ8d9f+esz32FujHqt30wp/7r1WDGA/+V0fA9koa5+/B2MPsg8GUG774Rgu4gM3riwcxOAT4DfMXdbzOzOcBfGGgu8xDwqvDh1eGU9Ubg4hSbfunurx9s35NOOskffvjh4bydiIiISN5p3biOn7z0BOa++KW86ie/z9lgSHdd/REe+9bXuOyuJ2hcfkROypDHMqqUbL4neT+wl6C3lxp3L0ozZfyYsrs/5O4XuPtt4fJ2dz8aOA54AfCiLJL28xlI2u8GLg3nAJeY2XnDOZ6IiIjIRBKPxXj2njv4yzsvpbi8gnO/8u2cjmB68pUfo7Sqmvv+85M5K8NEl83IqdOBbyaS7LHk7qtGsPvlicMAb3T37WZ2D7CF4FPN5WTR/EZEREQkX7k7Tase5Zlf/4xnbvoF7bt2Ulpdw8u+8X9Uz5iZ07JVNk7hxPd9hPu/8Cmee+RBZp546tA7ySGyaSrzCPCUu18+ZHAOmdlmgpFdN7j70sj6dcBSYLO7L0q3v5rKiIiISL5yd7paW2jf/Vww7dpJ26YNrP/9r2jduI7isjIWnHsBy1/7zyw870JKq6pyXWQAeg4e5AcnL6Fx+ZFcdNNfc/oNQJ7J6ERkc8f9c8BPzexGd/9DFvuPufDh0/nh4o6kzTsIEvf5Zlbm7j3jWjgRERHJW+6Ox+P9E4nXHj9kvcfjeCx26PpYLHwd638djyVeDyyTWJ9Y19dHvLeHeG8vsd4eYj09xHt66OvqorN1L10te+lsaaardS9dLS10tjTT0bSLWE9SCmPGnNNfwonv/TBLXnkRFZMacnMSB1FWU8OpH/p37rzqSn523qk0vuBIGpcfQePyI5i8/Ajq5s4f1x5vJppsEvcVwFrgt2a2hWDwpPYUce7uF42gbCNRz8Anl+QxdhNlLQLqCPqIzxvfO2ERxOO5Lobk0HC/BUsYjbsW2bz3WN8tyfZ8TBTpzt+o/dzZHGe4dTqadZTqvSPHH+y8jPm5TPMeGR0/Tcwh+/a/ThVriQJkGB/Z59B3TF/GEUj9cwwR7/68+SHCnzVxzvtj4/Fwn/gh72tmA/uQ2OfQRLv/vdIXLJMfd9yV10+ioqGRismNVE2dRuMLjqBq+kxqps+ievoMqqfPpHrGLKqnz8ybO+uDOfrN76SjuYnnHn6AbXfextpf/Kh/W0llJWU1tUBQn4l6taKicbk7f+H3f8mME04Z8/fJVrZ33BMWhlMqufztr4i8jiVtiy5XRjeY2TuBdwLMmzf+XSUBzDvj7MM+UZGhDfeP01gnJuPxvoMZbplG64/7cI+VTfxgRu2f1HCOk22djkZZB3vvyPGHTJ7dn1eeUf9gm/QeGR0/XUyq46RI0PvfP/LeKeMj+yT2O6R8Q5U1xfnLxKHnYOj9+xOyRMKdeE3SzxopkxUVHZrEJaZI4h/d14qLMSt6/n4pf2wPticmK4Iw3oqLD1nff7ziIoqi7xHGFBUXh/uE+xUXB+siy/3ri4LY4rIyikrLKC4tDeZlZRSXV1AxqYGikmzStfxVXFbGaR//j/7lrrZWWtatpeWZNbSsf5rejvaBD3Phhy0fp5ua5fWTxuV9spVNG/f5Q0cF3H3rsEs0CsxsGsGIrAB/dPdXRrb9EbgwXJzq7invuKuNu4iIiIiMk7Fp456rZHyY9hPc8Tcg+Tuj6nAeD+NERERERPLekIm7md2e5bHd3c/Jct8RcfcuM9sJzAaS+z6aFc636sFUEREREZkohmwqY2bZNirKeOTUsWBmvwZeS3BnfZ677zCzmcB2ggdTf+Lulw2y/x5grL9dmEKePRwrY0Z1XThU14VDdV04VNeFIZf13OzuK4cKyqSpTLqHT/PdTwkS9yLgJ2b2beBdDIwWe8NgO7v71LEtHpjZw+5+0li/j+Se6rpwqK4Lh+q6cKiuC8NEqOchE/cJ0qb9edz912b2B+CVwFnhlPBbd785NyUTERERERm+w72H+4uBq4FNQA+whaA7yzfksEwiIiIiIsN2eHUMmiR8+PSz4ZSPvp3rAsi4UV0XDtV14VBdFw7VdWHI+3oedj/uIiIiIiIy/g73pjIiIiIiIocFJe4iIiIiIhOAEvdxZmblZvYZM9tsZt1mtsnMPm1mZbkuW6EysyvMzNNMbZE4M7MrzWxtWHfbzeyrZlaX5riXmdljZtZlZk1m9gMzm5Em9mVmdp+ZdZhZq5n92syWpIk92cxuM7MD4XSrmeV191W5ZGbfTa7LpO05ryczW2ZmvwmP2W5m95rZ+WliZ5rZD81sj5l1mtmjZvbGTM/H4WywujazFwxynbuZHZcUr7rOM2Y2y8z+18y2mllPeF5uNLMjU8Tqup7AMq3rgryu3V3TOE7AjYCnmH6R67IV6gR8I02dONAWiftymph7geKkY743Tex6oC4p9lVALEVsE8HgYdHYk4GOFLHtwIm5Ppf5NgHnAL3JdZlP9QTMB/akiI0Br0iKrQM2pinzu3N9vvO8ri8Z5Dp34DjVdf5OBGPK7EpzPtqBUyKxuq4n8DTMui646zrnFVRIE3B+pILuAt4YzhPrzst1GQtxAu4Jz/8jwKuTpgvDmBWRC341cBnwm0jdvSNyvEbgQLj+WeCtBE+qJ2I/H4ktAbaF69uAdwOfj8T+JKms94Xre4B/BT7CQLJyb67PZT5NwEXhH9nEuWxL2p4X9UQwWFz/e4bH3hcubwVKIrHXRGK/HZb52XB5PzA51+c9H+s6jPlcuK07xXX+aqBedZ2/E/CzyPn4MXAp8L3IukfDOF3XE3zKtK7D2IK7rnNeQYU0EYzW6kAcmBOumxcuO/DjXJex0CbAIhfY9YPEfS5ysZ0WrqsKLzQH7o7Evj0S+8bI+6wP122LxJ4bib0qsv6v4bpOoCZctzh6wUdivxtZvzjX5zTXE7AE+FXknKRL3HNeT0BNuK8Dt0Zi/z0Se05kfeIP/DORdZdFYt+W6/Ofj3Udxv4h3PbUEMdUXefZBBQz8MHsqaRt/4ick5m6rif2NJy6DtcV3HWtNu7j6/RwvtHdtwO4+zZgQ9J2GT8LCL66AlhvZqeY2fkp2kIm6qYbeBDA3TuAB8L1p5hZSVIswN/DWAfuCNfNNbN56WJDt4fzCuCE8PWLh4hNPl6hegXBHViAm4DdaeLyoZ5ODPcdMtbM5gNzwnV3DXHcQpFpXQMcE86fDtvFrjSz5SniVNf5pxK4Fvg+8KOkbVsir+vRdT3RDaeuoQCvayXu48SCh0/nh4s7kjYnluebHlIdb8dGXl9DkJTfAuwws69HkvFl4XyXu8ci+yTqLlq/idg48FyKWIClSbHJ20caW+h2A+9399cCXWli8qGeRiP2OYKfIRpbSIasazObRPDtJsAFwFrgLwT/7P9mZtMi4arrPOPuB939and/m7t/MbE+/H95WrjYB2xH1/WENpy6LtTrWon7+Kkn+LoOgvZ3Ue3hvIiBu78yPo6JvK6MvC4C/gW4LlyeFM7T1R1AQ1Jsh7vHM4xNPvZIYwvZb4D57n7dEHH5UE8jjg3vJnYmxRaKTOs63XUOcDZwa+Smiep64vgMMCt8fYu7H0TX9eEqVV0X5HWtxH38VERex5K2RZeTf/lkbFUStHFvJniQpYqgd4rmcPu7zGwZA/WXSd1lE5u8faSxBcvdt7l7dwah+VBPoxEbXS6o+h9GXTcS9FIRA/6N4B/ocuDhcPuxwOXha9X1BGBmHwA+Fi7GgE+Fr3VdH2YGqeuCvK6VuI+f6D+X4qRt0eVOZNy4+8fdfRIw191/5+6d7n47QdePEHxLspKB+suk7rKJTd4+0lgZWj7U02jERpdV/ym4+03uPpPg4bMvufs+d18HvC8S9vJwrrrOc2Z2JfDfkVVXu/sj4Wtd14eRweq6UK9rJe7jJ9H7CAR3daOqw3k8jJNx5u7JbWNXR17PJrgrD+nrDqA1nCdiK83MMoxNPvZIY2Vo+VBPI44Ny16ZFCspZHCdg+o6r5nZe4CvRlZd7+7XRJZ1XR8mMqhroPCuayXu4yT8xdoZLs5M2pxot7XV3XvGr1SFzQLHmtmFZrYyaXNN5PVBggETAGYk/TNI1F0vQX+tRGKLgakpYiHomiwaC4f+Xow0VoaWD/U0GrEzGPhbrvpPwcwWmdm5ZvampE3J1zmorvOWmV0M/E9k1TcInkWK0nV9GMikrgv1ulbiPr4eDOfLzWw2BEPfEvRFDEGH/zJOwgdCfg/8EbjBzOojm6OJ/P0M1F0VcCqAmZUTjK4G8JC794WvH4zse3bkdaKLqe0EfbkOFntGOO8CHs0wFvQ7NBz5UE+PEPSQkEnsJoLR+QBeOkSsHOpTwG3Aj8ws2tVb8nUOquu8FHbx90MGOnn4pru/N/w7HqXreoIbRl0X5nU9lh3pa3reAAAXMdDJ/p0EI6f+PbJuZa7LWGgT8KXI+X+A4EGWrzEwKNZqgg+4J0biVhEMmPDryLorIseczcCADFuBtwDfisR+IRJbSdCdnRN8RXYFhw729POk8q4K1/cAHyYYza0nUf5cn898nAj6/nWePwBTXtQTAwOIOPDZ8Nht4fKzHDrq3nWR2G8SjLqXGA3wINCY6/Odp3V9YeS87QTeE9ZJYvC1TsIhz1XX+TkR3GRJnI+9wMWkGCVT1/XEn4ZR1wV5Xee8ggptSvqFjE435bpshTiFF/MTaeqkBTgyEntdmriHohdmGPvhNLGbgElJsa8jeKI8OXYvsDAp9nQG/ilFp07g1Fyfz3ycSJPM5Us9EfTtuzdFbAx4TVLs1LBsqcp8Za7Pda6nIer6O2nOWwy4VHWdvxPBQHnxNOciOh0Xxuu6nqBTFnVdcNd1ziup0CaCgXo+SdDWqRvYTPAJrTzXZSvUiaALqWvDuugh6F7qRwT9Q0fjjOBp9TVh3e0kSObr0xz3DQRfo3URfD32I2BWmtjzgLuBDoJP6jcBy9PEngjcTPAg8wGCIZuVtKev3y2kSebypZ4ImsvdSHAnqIPg69IL0sROB74HNIVlfhy4PNfnOR+mweqagbEZHg/PWxvBYGtnqq7zewLexNCJXH8yF+6j63oCTsOt60K8ri08iIiIiIiI5DE9nCoiIiIiMgEocRcRERERmQCUuIuIiIiITABK3EVEREREJgAl7iIiIiIiE4ASdxERERGRCUCJu4iIiIjIBKDEXUQkT5nZp83Mk6a4mbWb2Xoz+7aZvSAH5fpnM1sUWX5LWLYPjHdZREQKSUmuCyAiIkP6HcHoehDccKkDjgXeAVxmZq9z9z+NR0HM7IvAvwHHj8f7iYjIACXuIiL577fu/sPklWb2coIhu39hZse5+4ZxKMv0cXgPERFJQU1lREQmKHf/M/BJoDqci4jIYUyJu4jIxPY/QBdwkZkd8i2qmZ1tZreZ2b6wXfz9ZnZx8gHC9uk3mNlLzewhM+s0s81m9jkzq4jEbQHeHC4+Fi5HFZnZh8zsGTPrNrNNZvbvyeXKlJktCMv2aTN7VaRsTWb2HTObkuLn+K6ZnWVmd5tZh5k9Z2bXmFmxmR1hZjeb2UEz22Fm15lZVTZlExHJBSXuIiITmLt3AI8S3HU/LrHezN4O/BU4BvgF8C1gGnCjmV2V4lDHAjcDHcD1QCvwCeBPZpb4X/FV4Inw9bfC5aiPAp8B7gH+FygFPgv814h+SHglQZOg54CvAzuAtwM/TxH7QuBWYE9Yhm7g42F57wWKgW+EP9/7gM+PsGwiIuPG3D3XZRARkRTM7NPAp4C3pmrjHon7JfA64FXu/gczmwNsADYBZ7j73jCukiCZfyFwrLs/Ga5P/CO43t3fF64rAX4JvCb6/mb2Q4K77se7++PhurcAPwAOAie5+zPh+llhObqBRnePD/PnXwBsDhcvcfcbw/WlwGPAkcASd9+Y9HN80N2/Gq5bDjwdrr/W3T8Srq8DngW63F3t9kVkQtAddxGRia87nNeF88uAcuDqRNIO4O6dBB8Eihho8pJwkEg7eXfvA/41XLw0w3L8MpG0h8fYCTwCTAIaMjxGKpsSSXt43F6CDyAAS5NiuwnuqCdinwGaw8UvR9bvB9YC08IPNCIieU+9yoiITHy14fxgOD8xnJ9jZkclxdaE8+OS1q9y99boCnffaGYtBM1oMrE+xbrEB4eayOvhWpdi3b5wXp60/ll370la1w7UuPuupPVdkWN0Zlk2EZFxo8RdRGTiWxDON4XzSeH8ikH2mZy0vCNN3C5gSYbl6Bpkm2V4jFS6U6xLNItJPm77MI4hIjKhKHEXEZnAzKyBoK13G7AmXJ24877Y3Tel3PH50jUXmcRAUxMREckhtXEXEZnY3kVwE+YX7h4L160K5yclB5vZUjP7spm9MmnTiZHeYxKx84FZwIOR1erRQEQkR5S4i4hMUGZ2NnA1wR32/4xsugGIAZ83sxmR+BLgOuDDQGPS4WYy8DBqoueWr4SL34/E9YbzslH4EUREZBjUVEZEJP+9OuwaEYI23fXACcAZBA9VvsHdtyaC3X29mf0bcC3wlJn9jqDf8guAFcAfCZL7qAPAZ8MPA2uAc4CjgR+7+x8jcYm28Nea2V/d/TPD/WHC7iMXAD909y3D3V9EpFApcRcRyX//FE4JHcAWglFTv5roxzzK3b9iZk8T3F2/iGDgoY3h8vVhd49Rm4APEgyqdCZB/+kfIBjwKOp64HSCDw1HmNm1Wfw8bwHOAu4Mfw4REcmABmASESlw4cBFT7h7cheRY/mejwHvcff7x+s9RUQmOrVxFxGRcWVmRwLLGBjRVEREMqDEXURExtvbgLclD/gkIiKDU1MZEZECl4umMiIiMnxK3EVEREREJgA1lRERERH5/+3WAQkAAACAoP+v2wXAzCQAAAAdSURBVBHoCmFA3AEAYEDcAQBgQNwBAGBA3AEAYCAV5OjBjwxWEQAAAABJRU5ErkJggg==\n",
Parish, Chad's avatar
Parish, Chad committed
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
      "text/plain": [
       "<Figure size 864x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# the 'with' block sets temporary parameters:\n",
    "with plt.style.context({'font.sans-serif': 'Arial', 'svg.fonttype': 'none'}):\n",
    "    #Create a figure with unequal-sized subplots\n",
    "    fig = plt.figure(figsize=(12, 7)) \n",
    "    gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1]) \n",
    "    \n",
    "    # create the first subplot\n",
    "    ax0 = plt.subplot(gs[0])\n",
    "    ax0.plot(depth_in_nm, NRT, color='xkcd:brick red', label='Damage')  # plot the NRT damage\n",
    "    if err_both > 0:  # if the user input error bounds, plot the error bounds\n",
    "        ax0.fill_between(depth_in_nm, y1=(NRT * (1-err_both)), y2=(NRT * (1+err_both)), \n",
    "                     color='xkcd:merlot', alpha=0.15, label='Error bounds')\n",
    "        ax0.legend(loc=3, prop={'size': 14})  # only show legend if error bounds present\n",
    "    ax0.set_ylabel('NRT damage, dpa', fontsize=20)\n",
    "    \n",
    "    # Make annotations\n",
    "    pre, post = str(fluence * 1e4).split('e+')  # split fluence into A x 10^B, with conversion from cm^2 to m^2\n",
    "    ion_val_string = str(num_ions) + ' ' + ion_string + ' ions $\\Rightarrow$' + target\n",
    "    fluence_string = '$' + pre + r'\\times 10^{' + post + '}$ ions/m$^2$'\n",
    "    disp_string = str(displacement) +' eV $E_D$'\n",
    "    if err_both >0:\n",
    "        flu_err_string = r'Fluence $\\pm$ ' + str(err_fluence*100) + '%'\n",
    "        disp_err_string = r'$E_D \\pm$ ' + str(err_displace*100) + '%'\n",
    "    else:\n",
    "        flu_err_string = ''\n",
    "        disp_err_string = ''\n",
    "    ax0.text(0.98, 0.98, f'{ion_val_string}\\n{fluence_string}\\n{disp_string}\\n{flu_err_string}\\n{disp_err_string}',\n",
    "                 transform=ax0.transAxes, fontsize=14, verticalalignment='top',\n",
    "                 horizontalalignment='right')\n",
    "    \n",
    "    ax0.tick_params(axis='both', labelsize=20)\n",
    "    ax0.spines['right'].set_visible(False)\n",
    "    ax0.spines['top'].set_visible(False)\n",
    "    \n",
    "    # create the smaller plot\n",
    "    ax1 = plt.subplot(gs[1])\n",
    "    if atom_frac.max() < 0.01:  # Plot atom fraction if <1%\n",
    "        atom_string = 'Implanted\\natom fraction'\n",
    "        atom_mult = 1\n",
    "    else:\n",
    "        atom_string = 'Implanted\\natom %'\n",
    "        atom_mult = 100\n",
    "\n",
    "    ax1.plot(depth_in_nm, atom_frac*atom_mult, color='xkcd:brick red', label='Implantation')\n",
    "    if err_both > 0:  # if user input error bounds\n",
    "        ax1.fill_between(depth_in_nm, y1=(atom_frac * (1-err_fluence))*atom_mult,\n",
    "                         y2=(atom_frac * (1+err_fluence))*atom_mult,\n",
    "            color='xkcd:merlot', alpha=0.15, label='Error bounds')\n",
    "        ax1.legend(prop={'size': 14})   # only show legend if error bounds present\n",
    "    ax1.set_ylabel(atom_string, fontsize=20)\n",
    "    ax1.set_xlabel('Depth, nm', fontsize=20)\n",
    "\n",
    "    ax1.tick_params(axis='both', labelsize=20)\n",
    "    ax1.spines['right'].set_visible(False)\n",
    "    ax1.spines['top'].set_visible(False)\n",
    "    \n",
    "fn = input('File name (enter to skip):')\n",
    "if len(fn) > 0:\n",
    "    fn = '..\\\\' + fn  # go up one directory  \n",
    "    fig.savefig(fn + '.png', dpi=300)\n",
    "    fig.savefig(fn + '.svg', dpi=300)\n",
    "    fig.savefig(fn + '.eps', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
Parish, Chad's avatar
'minor'    
Parish, Chad committed
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-12T22:45:04.503487Z",
     "start_time": "2020-05-12T22:45:04.498499Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.45283368"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "NRT[50]"
   ]
Parish, Chad's avatar
Parish, Chad committed
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
Parish, Chad's avatar
Parish, Chad committed
554
   "version": "3.7.9"
Parish, Chad's avatar
Parish, Chad committed
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "position": {
    "height": "409.4px",
    "left": "1166px",
    "right": "20px",
    "top": "120px",
    "width": "351.4px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}