From e5f4a1c4df3dd9e19de5093390bbbf7fbc50a727 Mon Sep 17 00:00:00 2001 From: franciscovillaescusa Date: Wed, 28 Jun 2023 15:36:03 -0400 Subject: [PATCH] update --- docs/source/Examples/Density_fields.ipynb | 268 ++++++++++++++++++- docs/source/Examples/Reading_snapshots.ipynb | 1 + docs/source/tutorials.rst | 3 +- 3 files changed, 270 insertions(+), 2 deletions(-) create mode 100644 docs/source/Examples/Reading_snapshots.ipynb diff --git a/docs/source/Examples/Density_fields.ipynb b/docs/source/Examples/Density_fields.ipynb index ef0b0ef..d089fcc 100644 --- a/docs/source/Examples/Density_fields.ipynb +++ b/docs/source/Examples/Density_fields.ipynb @@ -1 +1,267 @@ -{"cells":[{"metadata":{},"id":"b63aadc6","cell_type":"markdown","source":"# Creating density fields from snapshots\n\n[![Binder](https://mybinder.org/badge_logo.svg)](https://binder.flatironinstitute.org/v2/user/fvillaescusa/Quijote?filepath=/Tutorials/Density_fields.ipynb)"},{"metadata":{"trusted":true},"id":"bfa73080","cell_type":"code","source":"import numpy as np\nimport readgadget\nimport MAS_library as MASL","execution_count":1,"outputs":[]},{"metadata":{},"id":"4ed62239","cell_type":"markdown","source":"Define the value of the parameters"},{"metadata":{"trusted":true},"id":"d999bae9","cell_type":"code","source":"snapshot = '/home/jovyan/Data/Snapshots/fiducial/0/snapdir_004/snap_004' #location of the snapshot\ngrid = 512 #the density field will have grid^3 voxels\nMAS = 'CIC' #Mass-assignment scheme:'NGP', 'CIC', 'TSC', 'PCS'\nverbose = True #whether to print information about the progress\nptype = [1] #[1](CDM), [2](neutrinos) or [1,2](CDM+neutrinos)","execution_count":2,"outputs":[]},{"metadata":{},"id":"9098a761","cell_type":"markdown","source":"Read the header and the particle positions"},{"metadata":{"trusted":true},"id":"8427b23a","cell_type":"code","source":"# read header\nheader = readgadget.header(snapshot)\nBoxSize = header.boxsize/1e3 #Mpc/h\nredshift = header.redshift #redshift of the snapshot\nMasses = header.massarr*1e10 #Masses of the particles in Msun/h\n\n# read positions, velocities and IDs of the particles\npos = readgadget.read_block(snapshot, \"POS \", ptype)/1e3 #positions in Mpc/h","execution_count":3,"outputs":[]},{"metadata":{},"id":"62506aef","cell_type":"markdown","source":"Print some information about the data"},{"metadata":{"trusted":true},"id":"d4808044","cell_type":"code","source":"print('BoxSize: %.3f Mpc/h'%BoxSize)\nprint('Redshift: %.3f'%redshift)\nprint('%.3f < X < %.3f'%(np.min(pos[:,0]), np.max(pos[:,0])))\nprint('%.3f < Y < %.3f'%(np.min(pos[:,1]), np.max(pos[:,1])))\nprint('%.3f < Z < %.3f'%(np.min(pos[:,2]), np.max(pos[:,2])))","execution_count":4,"outputs":[{"output_type":"stream","text":"BoxSize: 1000.000 Mpc/h\nRedshift: 0.000\n0.000 < X < 999.992\n0.000 < Y < 999.992\n0.000 < Z < 999.992\n","name":"stdout"}]},{"metadata":{},"id":"3d75aa92","cell_type":"markdown","source":"Define the matrix that will contain the value of the density / overdensity field"},{"metadata":{"trusted":true},"id":"951f674a","cell_type":"code","source":"delta = np.zeros((grid,grid,grid), dtype=np.float32)","execution_count":5,"outputs":[]},{"metadata":{},"id":"92a7fc19","cell_type":"markdown","source":"Now construct the 3D density field"},{"metadata":{"trusted":true},"id":"4c3b5a85","cell_type":"code","source":"# construct 3D density field\nMASL.MA(pos, delta, BoxSize, MAS, verbose=verbose)","execution_count":6,"outputs":[{"output_type":"stream","text":"\nUsing CIC mass assignment scheme\nTime taken = 5.884 seconds\n\n","name":"stdout"}]},{"metadata":{},"id":"05947788","cell_type":"markdown","source":"We can make some tests to make sure the density field has been computed properly"},{"metadata":{"trusted":true},"id":"ad8baf57","cell_type":"code","source":"# the sum of the values in all voxels should be equal to the number of particles\nprint('%.3f should be equal to\\n%.3f'%(np.sum(delta, dtype=np.float64), pos.shape[0]))","execution_count":7,"outputs":[{"output_type":"stream","text":"134217728.019 should be equal to\n134217728.000\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"As this point, delta contains the effective number of particles in each voxel. \nIf you want instead the effective mass in each voxel you can just do"},{"metadata":{"trusted":true},"cell_type":"code","source":"delta *= Masses[1]\n\n# now check that the mass in the density field is equal to the total mass in the simulation\nprint('%.3e should be equal to\\n%.3e'%(np.sum(delta, dtype=np.float64), pos.shape[0]*Masses[1]))","execution_count":8,"outputs":[{"output_type":"stream","text":"8.812e+19 should be equal to\n8.812e+19\n","name":"stdout"}]},{"metadata":{},"id":"77d940b3","cell_type":"markdown","source":"If needed, the overdensity is easy to calculate"},{"metadata":{"trusted":true},"id":"4593380c","cell_type":"code","source":"# at this point, delta contains the effective number of particles in each voxel\n# now compute overdensity and density constrast\ndelta /= np.mean(delta, dtype=np.float64); delta -= 1.0\n\nprint('%.3f < delta < %.3f'%(np.min(delta), np.max(delta)))\nprint(' = %.3f'%np.mean(delta))","execution_count":9,"outputs":[{"output_type":"stream","text":"-1.000 < delta < 1195.511\n = -0.000\n","name":"stdout"}]}],"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3 (ipykernel)","language":"python"},"language_info":{"name":"python","version":"3.7.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":5} \ No newline at end of file +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b63aadc6", + "metadata": {}, + "source": [ + "# Creating density fields from snapshots\n", + "\n", + "[![Binder](https://mybinder.org/badge_logo.svg)](https://binder.flatironinstitute.org/v2/user/fvillaescusa/Quijote?filepath=/Tutorials/Density_fields.ipynb)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "bfa73080", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import readgadget\n", + "import MAS_library as MASL" + ] + }, + { + "cell_type": "markdown", + "id": "4ed62239", + "metadata": {}, + "source": [ + "Define the value of the parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d999bae9", + "metadata": {}, + "outputs": [], + "source": [ + "snapshot = '/home/jovyan/Data/Snapshots/fiducial/0/snapdir_004/snap_004' #location of the snapshot\n", + "grid = 512 #the density field will have grid^3 voxels\n", + "MAS = 'CIC' #Mass-assignment scheme:'NGP', 'CIC', 'TSC', 'PCS'\n", + "verbose = True #whether to print information about the progress\n", + "ptype = [1] #[1](CDM), [2](neutrinos) or [1,2](CDM+neutrinos)" + ] + }, + { + "cell_type": "markdown", + "id": "9098a761", + "metadata": {}, + "source": [ + "Read the header and the particle positions" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "8427b23a", + "metadata": {}, + "outputs": [], + "source": [ + "# read header\n", + "header = readgadget.header(snapshot)\n", + "BoxSize = header.boxsize/1e3 #Mpc/h\n", + "redshift = header.redshift #redshift of the snapshot\n", + "Masses = header.massarr*1e10 #Masses of the particles in Msun/h\n", + "\n", + "# read positions, velocities and IDs of the particles\n", + "pos = readgadget.read_block(snapshot, \"POS \", ptype)/1e3 #positions in Mpc/h" + ] + }, + { + "cell_type": "markdown", + "id": "62506aef", + "metadata": {}, + "source": [ + "Print some information about the data" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "d4808044", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "BoxSize: 1000.000 Mpc/h\n", + "Redshift: 0.000\n", + "0.000 < X < 999.992\n", + "0.000 < Y < 999.992\n", + "0.000 < Z < 999.992\n" + ] + } + ], + "source": [ + "print('BoxSize: %.3f Mpc/h'%BoxSize)\n", + "print('Redshift: %.3f'%redshift)\n", + "print('%.3f < X < %.3f'%(np.min(pos[:,0]), np.max(pos[:,0])))\n", + "print('%.3f < Y < %.3f'%(np.min(pos[:,1]), np.max(pos[:,1])))\n", + "print('%.3f < Z < %.3f'%(np.min(pos[:,2]), np.max(pos[:,2])))" + ] + }, + { + "cell_type": "markdown", + "id": "3d75aa92", + "metadata": {}, + "source": [ + "Define the matrix that will contain the value of the density / overdensity field" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "951f674a", + "metadata": {}, + "outputs": [], + "source": [ + "delta = np.zeros((grid,grid,grid), dtype=np.float32)" + ] + }, + { + "cell_type": "markdown", + "id": "92a7fc19", + "metadata": {}, + "source": [ + "Now construct the 3D density field" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "4c3b5a85", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Using CIC mass assignment scheme\n", + "Time taken = 5.884 seconds\n", + "\n" + ] + } + ], + "source": [ + "# construct 3D density field\n", + "MASL.MA(pos, delta, BoxSize, MAS, verbose=verbose)" + ] + }, + { + "cell_type": "markdown", + "id": "05947788", + "metadata": {}, + "source": [ + "We can make some tests to make sure the density field has been computed properly" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "ad8baf57", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "134217728.019 should be equal to\n", + "134217728.000\n" + ] + } + ], + "source": [ + "# the sum of the values in all voxels should be equal to the number of particles\n", + "print('%.3f should be equal to\\n%.3f'%(np.sum(delta, dtype=np.float64), pos.shape[0]))" + ] + }, + { + "cell_type": "markdown", + "id": "a6deccd3", + "metadata": {}, + "source": [ + "As this point, delta contains the effective number of particles in each voxel. \n", + "If you want instead the effective mass in each voxel you can just do" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "1531aaa3", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "8.812e+19 should be equal to\n", + "8.812e+19\n" + ] + } + ], + "source": [ + "delta *= Masses[1]\n", + "\n", + "# now check that the mass in the density field is equal to the total mass in the simulation\n", + "print('%.3e should be equal to\\n%.3e'%(np.sum(delta, dtype=np.float64), pos.shape[0]*Masses[1]))" + ] + }, + { + "cell_type": "markdown", + "id": "77d940b3", + "metadata": {}, + "source": [ + "If needed, the overdensity is easy to calculate" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "4593380c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-1.000 < delta < 1195.511\n", + " = -0.000\n" + ] + } + ], + "source": [ + "# at this point, delta contains the effective number of particles in each voxel\n", + "# now compute overdensity and density constrast\n", + "delta /= np.mean(delta, dtype=np.float64); delta -= 1.0\n", + "\n", + "print('%.3f < delta < %.3f'%(np.min(delta), np.max(delta)))\n", + "print(' = %.3f'%np.mean(delta))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/source/Examples/Reading_snapshots.ipynb b/docs/source/Examples/Reading_snapshots.ipynb new file mode 100644 index 0000000..6531804 --- /dev/null +++ b/docs/source/Examples/Reading_snapshots.ipynb @@ -0,0 +1 @@ +{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Reading and manipulating snapshots\n\n[![Binder](https://mybinder.org/badge_logo.svg)](https://binder.flatironinstitute.org/v2/user/fvillaescusa/Quijote?filepath=/Tutorials/Snapshots.ipynb)"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport readgadget","execution_count":1,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"get the name of the snapshot"},{"metadata":{"trusted":true},"cell_type":"code","source":"snapshot = '/home/jovyan/Data/Snapshots/Om_p/32/snapdir_004/snap_004'","execution_count":2,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"read the header of the snapshot"},{"metadata":{"trusted":true},"cell_type":"code","source":"# read header\nheader = readgadget.header(snapshot)\nBoxSize = header.boxsize/1e3 #Mpc/h\nNall = header.nall #Total number of particles\nMasses = header.massarr*1e10 #Masses of the particles in Msun/h\nOmega_m = header.omega_m #value of Omega_m\nOmega_l = header.omega_l #value of Omega_l\nh = header.hubble #value of h\nredshift = header.redshift #redshift of the snapshot\nHubble = 100.0*np.sqrt(Omega_m*(1.0+redshift)**3+Omega_l)#Value of H(z) in km/s/(Mpc/h)\n\nprint('BoxSize = %.3f Mpc/h'%BoxSize)\nprint('Total number of particles:',Nall)\nprint('Masses of the particles:',Masses, 'Msun/h')\nprint('Omega_m = %.3f'%Omega_m)\nprint('Omega_L = %.3f'%Omega_l)\nprint('h = %.3f'%h)\nprint('redshift = %.3f'%redshift)\nprint('H(z=%.1f)=%.3f (km/s)/(Mpc/h)'%(redshift,Hubble))","execution_count":3,"outputs":[{"output_type":"stream","text":"BoxSize = 1000.000 Mpc/h\nTotal number of particles: [ 0 134217728 0 0 0 0]\nMasses of the particles: [0.00000000e+00 6.77240019e+11 0.00000000e+00 0.00000000e+00\n 0.00000000e+00 0.00000000e+00] Msun/h\nOmega_m = 0.328\nOmega_L = 0.672\nh = 0.671\nredshift = 0.000\nH(z=0.0)=100.000 (km/s)/(Mpc/h)\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"For N-body simulations, we only care about particle type 1 (and type 2 if neutrinos are included)"},{"metadata":{"trusted":true},"cell_type":"code","source":"mass_c = Masses[1]\nN_c = Nall[1]\nprint('Mass of a DM particle = %.3e Msun/h'%mass_c)\nprint('Number of DM particles = %d'%N_c)","execution_count":4,"outputs":[{"output_type":"stream","text":"Mass of a DM particle = 6.772e+11 Msun/h\nNumber of DM particles = 134217728\n","name":"stdout"}]},{"metadata":{"trusted":true},"cell_type":"code","source":"# we can check the value of Omega_m\nrho_crit = 2.775e11 #critical density at z=0 in (Msun/h)/(Mpc/h)^3\nestimated_Omega_m = N_c*mass_c/BoxSize**3/rho_crit\nprint('%.4f should be similar to\\n%.4f'%(estimated_Omega_m,Omega_m))","execution_count":5,"outputs":[{"output_type":"stream","text":"0.3276 should be similar to\n0.3275\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"Now lets read the positions, velocities, and IDs of the DM particles"},{"metadata":{"trusted":true},"cell_type":"code","source":"ptype = [1] #DM is 1, neutrinos is [2]\npos = readgadget.read_block(snapshot, \"POS \", ptype)/1e3 #positions in Mpc/h\nvel = readgadget.read_block(snapshot, \"VEL \", ptype) #peculiar velocities in km/s\nids = readgadget.read_block(snapshot, \"ID \", ptype)-1 #IDs starting from 0","execution_count":6,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets print some information about these quantities"},{"metadata":{"trusted":true},"cell_type":"code","source":"print('%.3f < X < %.3f Mpc/h'%(np.min(pos[:,0]), np.max(pos[:,0])))\nprint('%.3f < Y < %.3f Mpc/h'%(np.min(pos[:,1]), np.max(pos[:,1])))\nprint('%.3f < Z < %.3f Mpc/h'%(np.min(pos[:,2]), np.max(pos[:,2])))\nprint('%.3f < Vx < %.3f km/s'%(np.min(vel[:,0]), np.max(vel[:,0])))\nprint('%.3f < Vy < %.3f km/s'%(np.min(vel[:,1]), np.max(vel[:,1])))\nprint('%.3f < Vz < %.3f km/s'%(np.min(vel[:,2]), np.max(vel[:,2])))\nprint('%d < IDs < %d'%(np.min(ids), np.max(ids)))","execution_count":7,"outputs":[{"output_type":"stream","text":"0.000 < X < 999.992 Mpc/h\n0.000 < Y < 999.992 Mpc/h\n0.000 < Z < 999.992 Mpc/h\n-4777.000 < Vx < 5332.000 km/s\n-4387.000 < Vy < 4999.000 km/s\n-4977.000 < Vz < 4632.000 km/s\n0 < IDs < 134217727\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"You can get the position, velocity, and ID of a particle just by calling its index"},{"metadata":{"trusted":true},"cell_type":"code","source":"# lets consider the particle number 10\nprint('position =',pos[10],'Mpc/h')\nprint('velocity =',vel[10],'km/s')\nprint('ID =',ids[10])","execution_count":8,"outputs":[{"output_type":"stream","text":"position = [ 9.89725 996.024 15.1425 ] Mpc/h\nvelocity = [ 356.25 -327.125 -153. ] km/s\nID = 10\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"The particles IDs can be used to track particles across times. Lets take the particle with ID equal to 623 and find\nits position across redshifts"},{"metadata":{"trusted":true},"cell_type":"code","source":"part_ID = 620\nfor snapnum in [0,1,2,3,4]:\n snapshot = '/home/jovyan/Data/Snapshots/Om_p/32/snapdir_%03d/snap_%03d'%(snapnum,snapnum)\n \n # read header\n header = readgadget.header(snapshot)\n redshift = header.redshift #redshift of the snapshot\n \n # read positions and ids\n pos = readgadget.read_block(snapshot, \"POS \", [1])/1e3 #positions in Mpc/h\n ids = readgadget.read_block(snapshot, \"ID \", [1])-1 #IDs starting from 0\n \n index = np.where(ids==part_ID)[0]\n position = pos[index][0]\n print('z=%.1f -----> (X,Y,Z)=(%.2f, %.2f, %.2f) Mpc/h'%(redshift,position[0],position[1],position[2]))","execution_count":9,"outputs":[{"output_type":"stream","text":"z=3.0 -----> (X,Y,Z)=(2.14, 16.49, 86.31) Mpc/h\nz=2.0 -----> (X,Y,Z)=(2.87, 16.15, 86.48) Mpc/h\nz=1.0 -----> (X,Y,Z)=(4.18, 15.80, 86.84) Mpc/h\nz=0.5 -----> (X,Y,Z)=(4.97, 15.25, 87.17) Mpc/h\nz=0.0 -----> (X,Y,Z)=(6.31, 14.28, 87.94) Mpc/h\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"Keep in mind the simulations have periodic boundary conditions. For instance, this is the incorrect and correct way\nto compute the distance between them"},{"metadata":{"trusted":true},"cell_type":"code","source":"particle1 = pos[3]\nparticle2 = pos[4]\nprint('Position of particle 1: (%.3f, %.3f, %.3f) Mpc/h'%(particle1[0], particle1[1], particle1[2]))\nprint('Position of particle 2: (%.3f, %.3f, %.3f) Mpc/h'%(particle2[0], particle2[1], particle2[2]))","execution_count":10,"outputs":[{"output_type":"stream","text":"Position of particle 1: (4.534, 3.950, 998.616) Mpc/h\nPosition of particle 2: (4.922, 4.519, 2.621) Mpc/h\n","name":"stdout"}]},{"metadata":{"trusted":true},"cell_type":"code","source":"# this would be the incorrect way to compute the distance\nd = np.sqrt(np.sum((particle1-particle2)**2))\nprint('Incorrect distance = %.3f Mpc/h'%d)\n\n# this would be the correct way to compute the distance\nd = particle1-particle2\nindexes = np.where(d>BoxSize/2)\nd[indexes]-=BoxSize\nindexes = np.where(d<-BoxSize/2)\nd[indexes]+=BoxSize\nd = np.sqrt(np.sum(d**2))\nprint('Correct distance = %.3f Mpc/h'%d)","execution_count":11,"outputs":[{"output_type":"stream","text":"Incorrect distance = 995.995 Mpc/h\nCorrect distance = 4.064 Mpc/h\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"In simulations with massive neutrinos, you can read both dark matter and neutrino positions, velocities, and IDs"},{"metadata":{"trusted":true},"cell_type":"code","source":"# get the name of the snapshot\nsnapshot = '/home/jovyan/Data/Snapshots/Mnu_p/284/snapdir_002/snap_002'\n\n# read header\nheader = readgadget.header(snapshot)\nBoxSize = header.boxsize/1e3 #Mpc/h\nNall = header.nall #Total number of particles\nMasses = header.massarr*1e10 #Masses of the particles in Msun/h\nOmega_m = header.omega_m #value of Omega_m\nOmega_l = header.omega_l #value of Omega_l\nh = header.hubble #value of h\nredshift = header.redshift #redshift of the snapshot\nHubble = 100.0*np.sqrt(Omega_m*(1.0+redshift)**3+Omega_l)#Value of H(z) in km/s/(Mpc/h)\n\nprint('BoxSize = %.3f Mpc/h'%BoxSize)\nprint('Total number of particles:',Nall)\nprint('Masses of the particles:',Masses, 'Msun/h')\nprint('Omega_m = %.3f'%Omega_m)\nprint('Omega_L = %.3f'%Omega_l)\nprint('h = %.3f'%h)\nprint('redshift = %.3f'%redshift)\nprint('H(z=%.1f)=%.3f (km/s)/(Mpc/h)'%(redshift,Hubble))","execution_count":12,"outputs":[{"output_type":"stream","text":"BoxSize = 1000.000 Mpc/h\nTotal number of particles: [ 0 134217728 134217728 0 0 0]\nMasses of the particles: [0.00000000e+00 6.51631041e+11 4.92989376e+09 0.00000000e+00\n 0.00000000e+00 0.00000000e+00] Msun/h\nOmega_m = 0.318\nOmega_L = 0.682\nh = 0.671\nredshift = 1.000\nH(z=1.0)=179.513 (km/s)/(Mpc/h)\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"As can be seeing, particle type 2 (neutrinos) have millions particles and the masses are not zero"},{"metadata":{"trusted":true},"cell_type":"code","source":"mass_c = Masses[1]\nmass_n = Masses[2]\nN_c = Nall[1]\nN_n = Nall[2]\nprint('Mass of a DM particle = %.3e Msun/h'%mass_c)\nprint('Mass of a NU particle = %.3e Msun/h'%mass_n)\nprint('Number of DM particles = %d'%N_c)\nprint('Number of NU particles = %d'%N_n)\n\nOmega_m_estimated = (N_c*mass_c + N_n*mass_n)/BoxSize**3/rho_crit\nOmega_c_estimated = (N_c*mass_c)/BoxSize**3/rho_crit\nOmega_n_estimated = (N_n*mass_n)/BoxSize**3/rho_crit\nprint('Omega_cb = %.3f'%Omega_c_estimated)\nprint('Omega_nu = %.3e'%Omega_n_estimated)\nprint('Omega_m = %.3f'%Omega_m_estimated)","execution_count":13,"outputs":[{"output_type":"stream","text":"Mass of a DM particle = 6.516e+11 Msun/h\nMass of a NU particle = 4.930e+09 Msun/h\nNumber of DM particles = 134217728\nNumber of NU particles = 134217728\nOmega_cb = 0.315\nOmega_nu = 2.384e-03\nOmega_m = 0.318\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"Now lets read the positions, velocities, and IDs of the dark matter and neutrino particles"},{"metadata":{"trusted":true},"cell_type":"code","source":"pos_c = readgadget.read_block(snapshot, \"POS \", [1])/1e3 #positions in Mpc/h\nvel_c = readgadget.read_block(snapshot, \"VEL \", [1]) #peculiar velocities in km/s\nids_c = readgadget.read_block(snapshot, \"ID \", [1])-1 #IDs starting from 0\n\npos_n = readgadget.read_block(snapshot, \"POS \", [2])/1e3 #positions in Mpc/h\nvel_n = readgadget.read_block(snapshot, \"VEL \", [2]) #peculiar velocities in km/s\nids_n = readgadget.read_block(snapshot, \"ID \", [2])-1 #IDs starting from 0","execution_count":14,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Lets make a plot with the distribution of the dark matter and neutrino velocities"},{"metadata":{"trusted":true},"cell_type":"code","source":"# lets compute the modulus of the dark matter and neutrino velocities\nVc = np.sqrt(vel_c[:,0]**2 + vel_c[:,1]**2 + vel_c[:,2]**2)\nVn = np.sqrt(vel_n[:,0]**2 + vel_n[:,1]**2 + vel_n[:,2]**2)\nprint('%.3f < Vc < %.3f'%(np.min(Vc), np.max(Vc)))\nprint('%.3f < Vn < %.3f'%(np.min(Vn), np.max(Vn)))\n\nbins_histo = np.logspace(0,5,1000)\nhisto_Vc, edges = np.histogram(Vc, bins_histo)\nhisto_Vn, edges = np.histogram(Vn, bins_histo)","execution_count":15,"outputs":[{"output_type":"stream","text":"0.811 < Vc < 4902.378\n13.396 < Vn < 60315.773\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"As can be seen, neutrinos have, on average, larger velocities than dark matter"},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nplt.xscale('log')\nplt.yscale('log')\nplt.xlabel('V')\nplt.ylabel('Particle number')\nplt.plot(edges[1:], histo_Vc)\nplt.plot(edges[1:], histo_Vn)\nplt.legend(['Dark matter', 'Neutrinos'])\nplt.show()","execution_count":16,"outputs":[{"output_type":"display_data","data":{"text/plain":"
","image/png":"iVBORw0KGgoAAAANSUhEUgAAAjoAAAG1CAYAAADwRl5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAACDQElEQVR4nO3dd3xT9frA8U/SvQuUtrQUyoYyWmRfmVplXFDB/XMAKl61qFdcuEC9Ck4ujiou4Or1Kg4EJyp7yIayyqZAGV2UbrqS/P44bdrQlbZJTsbzfr36ysnJycnDoU2efMfz1RgMBgNCCCGEEE5Iq3YAQgghhBDWIomOEEIIIZyWJDpCCCGEcFqS6AghhBDCaUmiI4QQQginJYmOEEIIIZyWJDpCCCGEcFqS6AghhBDCabmrHYDa9Ho9586dIyAgAI1Go3Y4QgghhDCDwWAgPz+fiIgItNq6221cPtE5d+4cUVFRaochhBBCiCZITU2lbdu2dT7u8olOQEAAoFyowMBAlaMRQgghhDny8vKIiooyfo7XxeUTncruqsDAQEl0hBBCCAfT0LATGYwshBBCCKcliY4QQgghnJbLd12ZQ6/XU1paqnYYwk54eHjg5uamdhhCCCHM4LKJTmJiIomJieh0unqPKy0tJSUlBb1eb6PIhCMIDg4mPDxcShIIIYSd0xgMBoPaQagpLy+PoKAgcnNzawxGNhgMnD59mrKysgbn6QvXYDAYKCoqIiMjg+DgYNq0aaN2SEII4ZLq+/yuzmVbdMxRXl5OUVERERER+Pr6qh2OsBM+Pj4AZGRkEBoaKt1YQghhx6SJoh6V3Vqenp4qRyLsTWXiW1ZWpnIkQggh6iOJjhlkHIa4nPxOCCGEY5BERwghhBBOSxIdIYQQQjgtSXREk0yZMoUbbrhB7TCEEEKIekmi44SmTJmCRqNBo9Hg4eFBWFgY11xzDQsXLnTKekAajYZly5aZ7HvxxReJi4tTJR4hhBD2Q6aXO6kxY8awaNEidDod6enprFixgkcffZTvvvuOH3/8EXf3pv3X63Q6lxqIW1paKrPuXFBKViE7T11kT2oO53OLAfDy0OLn6UZsVDDDu7QmqqWUnBDCEUiLTiMYDAaKSstV+WlsXUcvLy/Cw8OJjIzkiiuu4Nlnn2X58uX89ttvLF682HjcvHnz6N27N35+fkRFRfHQQw9RUFBgfHzx4sUEBwfz448/EhMTg5eXF6dPn67xetu3b6d169a8/vrrtcZz8uRJNBoN33zzDcOGDcPHx4cBAwZw5MgRtm/fTv/+/fH392fs2LFkZmaanPeaa64hJCSEoKAgRowYwa5du4yPR0dHAzBx4kQ0Gg3R0dEsXryYl156iT179hhbtir/zTk5Odx33320bt2awMBArrrqKvbs2WM8X2VL0KeffkqHDh3w9vZu1HUXjquwpJxPN5xgzPz1jHprLU98u4cvtpxi5cF0Vh5M55e95/lmxxme+2E/w95Ywx2fbmHNoQy1wxZCNMApWnRSUlK45557SE9Px83NjS1btuDn52fx17lUpiNm1u8WP685kl8eja9n8/67rrrqKmJjY1m6dCn33XcfAFqtlnfffZcOHTpw4sQJHnroIZ566ik++OAD4/OKiop4/fXX+fTTT2nVqhWhoaEm5129ejWTJk3ijTfe4P777683htmzZzN//nzatWvHPffcw//93/8REBDAO++8g6+vL7fccguzZs3iww8/BCA/P5/Jkyfz3nvvYTAYePvttxk3bhxHjx4lICCA7du3ExoayqJFixgzZgxubm74+/uzf/9+VqxYwcqVKwEICgoC4Oabb8bHx4fffvuNoKAgPvroI66++mqOHDlCy5YtATh27Bjff/89S5culWKALqCotJzENcdYtOkkRaUVtbPctMRGBREXFUx0iB9ajYaSMh3ZhaVsPnGBnacusunYBTYdu8C43uH86/petPL3UvlfIoSojVMkOlOmTOGVV15h2LBhZGdn4+Ulbzh16d69O3v37jXe/+c//2ncjo6O5pVXXuGBBx4wSXTKysr44IMPiI2NrXG+H374gbvvvptPP/2UW2+9tcHXf+KJJxg9ejQAjz76KLfffjurVq3iyiuvBODee+81aXG66qqrTJ7/8ccfExwczLp16xg/fjytW7cGqtaequTv74+7u7vJvo0bN7Jt2zYyMjKMvyNvvfUWy5Yt47vvvjMmaaWlpXz++efGcwvnZDAY+P1AOv/6OZmzOZcAiG7ly73DOnJdnwiCfD3qfG5qdhGL/zrJ4r9O8uu+NLacyOa1Sb25tmd4nc8RQqjD4ROdAwcO4OHhwbBhwwCM38qtwcfDjeSXR1vt/A29tiUYDAaTMTYrV65k7ty5HDp0iLy8PMrLyykuLqaoqMhY/dfT05M+ffrUONfWrVv5+eef+e6778yegVX9PGFhYQD07t3bZF9GRlV3QHp6Os8//zxr164lIyMDnU5HUVFRrd1nDdmzZw8FBQW0atXKZP+lS5c4fvy48X779u0lyXFyFwpKmLl0H38mpwMQGezDC+NjuDYmDK224TFoUS19eWF8DBP7RvLEt3s4lJbPP/67k5eu68ndQ6KtHL1oFIMBclPh/F7IP1+1380TgqOgVRcIjARZy9BpqZ7orF+/njfffJOdO3dy/vx5fvjhhxofmomJibz55pukpaURGxvLe++9x8CBAwE4evQo/v7+TJgwgbNnz3LTTTfx7LPPWiVWjUbT7O4jtR08eJAOHToAyriZ8ePH8+CDD/Lqq6/SsmVLNm7cyL333ktpaakx0fHx8al1AHKnTp1o1aoVCxcu5O9//zseHnV/A65U/ZjKc16+r/rMsMmTJ3PhwgXeeecd2rdvj5eXF0OGDKG0tLTR//aCggLatGnD2rVrazwWHBxs3LZGt6ewH0fT85m6eDtnLl7Cw03DAyM68dDIzvh4Nv7LRK/IIH6cPpR//ZzMF1tOMWv5Abzctdw6oJ0VIhdmMRjg7C44+COcT4Lze+DSxfqf4+kP3cZB7K3QcRRopcvamaj+qV1YWEhsbCz33HMPkyZNqvH4kiVLmDFjBgsWLGDQoEHMnz+f0aNHc/jwYUJDQykvL2fDhg0kJSURGhrKmDFjGDBgANdcc40K/xr7tnr1avbt28djjz0GwM6dO9Hr9bz99tvGldm/+eYbs88XEhLC0qVLGTlyJLfccgvffPONWclOY2zatIkPPviAcePGAZCamkpWVpbJMR4eHsZ1ySp5enrW2HfFFVeQlpaGu7u7cRCzcC2bjmXxwH93kl9cTvtWvnx4Rz9iIupe9dgcnu5aXr6+Jz6ebny8/gQzl+7Dx9Od62IjLBS1MMulHNi5CHZ9AdnHTR/TukPrHtCiPWi0oNFA2SW4eBKyU6C0APZ9o/z4h8PAaTDoH+AVoMa/RFiY6onO2LFjGTt2bJ2Pz5s3j2nTpjF16lQAFixYwC+//MLChQuZOXMmkZGR9O/fn6ioKADGjRtHUlJSnYlOSUkJJSUlxvt5eXkW/NfYj5KSEtLS0kyml8+dO5fx48dz9913A9C5c2fKysp47733mDBhAps2bWLBggWNep3Q0FBWr17NqFGjuP322/n666+bPHW9Nl26dOGLL76gf//+5OXl8eSTTxpXD68UHR1tHOfj5eVFixYtiI6OJiUlhaSkJNq2bUtAQADx8fEMGTKEG264gTfeeIOuXbty7tw5fvnlFyZOnEj//v0tFrewP9/sSOXZpfso1xvo374FH9/dn5Z+likdoNFoeGZsd/KLy/lq22ke/yaJiCBv+kdbrytdVCjOhc2JsOVDKKl4P/fwhe5/h+hh0KYPhMaAex1jN3XlcG437F0C+7+DgjRY/S/YugDiX4TY/5NuLQdn1/97paWl7Ny5k/j4eOM+rVZLfHw8mzdvBmDAgAFkZGRw8eJF9Ho969evp0ePHnWec+7cuQQFBRl/KhMkZ7NixQratGlDdHQ0Y8aMYc2aNbz77rssX77cOJMoNjaWefPm8frrr9OrVy++/PJL5s6d2+jXCg8PN7YW3XHHHTVaUprjs88+4+LFi1xxxRXcddddPPLIIzVmfb399tv8+eefREVF0bdvXwBuvPFGxowZw6hRo2jdujVfffUVGo2GX3/9leHDhzN16lS6du3KbbfdxqlTp4zjhYTz0esNvPX7YZ76bi/legMTYiP4732DLJbkVNJoNLxyQy/G9gqnTGfgoS93kV3Y+C5WYSZdGfz1PrwTC+teV5Kc0Bi4PhGeOAo3fgr9JkNE37qTHAA3d4gaAH9/Cx4/ApM+gZYdoTATlifA59dBgZQRcGQaQ2MLtFiRRqMxGaNz7tw5IiMj+euvvxgyZIjxuKeeeop169axdetWAH777TeeeuopDAYD1157LfPmzavzNWpr0YmKiiI3N5fAQNMm7OLiYlJSUqSeiqhBfjccQ7lOz+Pf7mF50jkAHr6qM4/FdzVrwHFTFZWWc937mziWUcCYnuF8eOcVLlVk0ybSk2HZA8r4G4CQrjDqWehxvWVaX8pLYdtHsGYulBVCQATc8rmSEAm7kZeXR1BQUK2f39XZdYuOucaOHcu+ffvYv39/vUkOKIX0AgMD+eKLLxg8eDBXX321jaIUQtiSTm/giYokx12r4c2b+vD4td2smuQA+Hq6M//WONy1GlYcSOO7nWes+nouRVcO69+Cj4YrSY53MEx4Fx7cDD0nWq6Lyd0T/vYw3L9WSaLyz8GisbD9U2Wws3Aodp3ohISE4ObmRnp6usn+9PR0k/ooTZGQkEBycjLbt29v1nmEEPZHpzfw5Hd7WFaR5Hx4Zz9u7m+7bupekUE8dk1XAF788YCxTo9ohuwU+OwaZfyMvgy6joWErUr3lJuVhpu27grTVkOP65TX/OVx+ON5SXYcjF0nOp6envTr149Vq1YZ9+n1elatWmXSlSWEEJUMBgMvLN/P0l1ncdNqePf2vlwTY/sxWA+M6ES/9i0oLNXx0o8HbP76TuXI7/DxCDi3C7yDYOJHcPtXEGCDAo1eAUq3VfyLyv3N78OvT0qy40BUT3QKCgpISkoiKSkJwDhTprIg3IwZM/jkk0/4z3/+w8GDB3nwwQcpLCw0zsJqqsTERGJiYhgwQPpchXAmiWuO8b+tp9Fq4N+3xjGudxtV4nDTapgzsTfuWg1/JKez+lB6w08SpgwGWP8m/O8WZXZV2wFKN1XsbcoUcVvRaGDoYzDhHUAD2z+BVS/Z7vVFs6g+GHnt2rWMGjWqxv7JkycblwJ4//33jQUD4+LiePfddxk0aJBFXr++wUwy4FTURX437NMPu8/w2BJlgOrL19tHleI5vx7k4/Un6Brmz2+PDsfNymOEnEZZMfz4sFLbBqD/vTDmNWX8jJp2fQE/Tle2//42DLhP3XhcmLmDkVVPdNQmiY5oCvndsD9/Hcti8qJtlOkM3D+8I8+Oq7vMhC3lFpUx7I3V5BWXM//WOG7oG6l2SPavOBe++j84tVEp9jfuTeh/j9pRVVn/Jqx+BbQeMPU3mY2lEpeaddUU0nUlhPM4lpHPP77YSZnOwPg+bZg5prvaIRkF+XrwjxGdAJj35xHKdPoGnuHi8tNh0d+VJMczAO783r6SHIBhT0DM9coA5W/uhoJMtSMS9XDZREdmXQnhHApLynngv7vILylnQHQL3ro51upTyBtr6pXRhPh7cjq7iG93yHTzOuWchoWjIX0f+IXC1F+h40i1o6pJo1EKE1ZOPf9uqjL1Xdgll010hO1FR0czf/58tcMQTsRgMPDM0n0cyyggLNCLD+/sh7eH/S3I6OvpTsKozgC8u+ooxWWWqx7uNPLOw38mwMUUCG4P9/6uLN9gr7wC4Nb/KguCntygTHsXdkkSHSc0ZcoUNBoNr732msn+ZcuWWbRC68mTJ9FoNMYZcw3Zvn07999/v8VeX4j/bj3Nj3vO4abV8P7/XUGIfz2l/lX2f4PaERHkTVpeMf/dckrtcOxL4QX44gZlkc0W0XDPCmUZBnvXuhtc/76yvWk+HF+jajiidi6b6Dj7GB1vb29ef/11Ll68qHYolJYq6/20bt0aX19flaMRzmJPag7/+ikZgJljujPAzhfQ9HJ349H4LgB8suEEpeUyVgeA4jz47yTIPKQstXD3cgh0oJXfe06smnn106NQWqhuPKIGl010nH2MTnx8POHh4fUu0rlx40aGDRuGj48PUVFRPPLIIxQWVv2RajQali1bZvKc4OBg47T/Dh06ANC3b180Gg0jR44ElBalG264gVdffZWIiAi6desG1Oy60mg0fPrpp0ycOBFfX1+6dOnCjz/+aPJ669atY+DAgXh5edGmTRtmzpxJeXlVX/h3331H79698fHxoVWrVsTHx5v8G4RzulhYykNf7qJUp+famDDuG9ZB7ZDMMrFvW1oHeJGeV8Kv+86rHY76Sovgq9vgfBL4tlKSnBbRakfVePEvQmBbyDkFq19VOxpxGZdNdJrEYFCydTV+GlkFwM3NjTlz5vDee+9x5kzNwY/Hjx9nzJgx3Hjjjezdu5clS5awceNGpk+fbvZrbNu2DYCVK1dy/vx5li5danxs1apVHD58mD///JOff/65znO89NJL3HLLLezdu5dx48Zxxx13kJ2dDcDZs2cZN24cAwYMYM+ePXz44Yd89tlnvPLKKwCcP3+e22+/nXvuuYeDBw+ydu1aJk2ahItXTHB6er2BGd8kcTbnEu1b+fLmzbEOs2imp7uWuwe3B+CzjSmu/btaXgrf3AWnNoFXENz1g7LkgiPyCoAJ85XtrR/CmR2qhiNMWWmBECdVVgRzVGpSffYcePo16ikTJ04kLi6O2bNn89lnn5k8NnfuXO644w7++c9/AtClSxfeffddRowYwYcffmhWbZjWrVsD0KpVqxprj/n5+fHpp5/i6Vl/ca8pU6Zw++23AzBnzhzeffddtm3bxpgxY/jggw+Iiori/fffR6PR0L17d86dO8fTTz/NrFmzOH/+POXl5UyaNIn27ZUPj969e5t1bYTj+nDdcdYczsTLXcsHd1xBkI+H2iE1yh2D2/P+mmPsO5vL9pMXGdjBvrvcrEJXDt/fC8dWgocv3PENtIlVO6rm6XIN9LkN9n4NPz4CD2wArf0NjHdFLtui4+xjdCq9/vrrxuUzqtuzZw+LFy/G39/f+DN69Gj0ej0pKSnNft3evXs3mOQA9OlTNavCz8+PwMBAMjIyADh48CBDhgwx+bZ+5ZVXUlBQwJkzZ4iNjeXqq6+md+/e3HzzzXzyySd2MSZJWM9fx7J4+4/DAPzr+l70jAhSOaLGa+nnyaQr2gLw6YYTKkejAr0efnoEDv4Ibp5w25fQbrDaUVnGmLnKiuoZB2DPV2pHIyq4bItOQkICCQkJxsqKZvHwVVpW1ODRtEG8w4cPZ/To0TzzzDNMmTLFuL+goIB//OMfPPLIIzWe065dO0AZQ3N503pZWZlZr+vnZ17rk4eH6bdxjUaDXm/eIE03Nzf+/PNP/vrrL/744w/ee+89nnvuObZu3WocPyScR3peMY98vRu9AW7u15ZbBthuNXJLu3doNF9tO82fB9M5daGQ9q0a11rrsAwG+P0ZSPoSNG5w00LodJXaUVmOb0sY/iT88ZwyVqfnJPCUCRhqc9kWnSbRaJTuIzV+mjEG4bXXXuOnn35i8+bNxn1XXHEFycnJdO7cucZPZUtM69atOX++asDk0aNHKSoqMt6vPE6ns05NkB49erB582aTZGvTpk0EBATQtq3yjVij0XDllVfy0ksvsXv3bjw9Pfnhhx+sEo9QT7lOz8P/201WQSndwwN4+fpeaofULJ1DAxjRtTUGAyzadFLtcGxncyJsXaBs3/AB9JigbjzWMHAaBLVTCglu/VDtaASS6LiE3r17c8cdd/Duu+8a9z399NP89ddfTJ8+naSkJI4ePcry5ctNBiNfddVVvP/+++zevZsdO3bwwAMPmLTAhIaG4uPjw4oVK0hPTyc3N9eicT/00EOkpqby8MMPc+jQIZYvX87s2bOZMWMGWq2WrVu3MmfOHHbs2MHp06dZunQpmZmZ9OhhH2scCcv5aP0Jtp3Mxt/LnQ/v7IePp+OPfbh3qNLq+N3OMxSWuEBV3cO/wZ8vKNuj5ygrkDsjdy+4uuLfuXG+UiNIqEoSHRfx8ssvm3QJ9enTh3Xr1nHkyBGGDRtG3759mTVrFhERVYOt3377baKiohg2bBj/93//xxNPPGFSB8fd3Z13332Xjz76iIiICK6//nqLxhwZGcmvv/7Ktm3biI2N5YEHHuDee+/l+eefByAwMJD169czbtw4unbtyvPPP8/bb7/N2LFjLRqHUNeBc7nMX3kEgBev60mHEOfo5hnaOYToVr4UlJTz816VusRtJXU7fDsFDHq4YjIMfkjtiKyr100Q3gdK8mDjPLWjcXmyermsXi6aQH43bKO4TMf172/icHo+18aE8dFd/RxmKrk5Fqw7zmu/HSI2KpjlCVeqHY51FGTCR8OVrpwuo5XBx26ONVOuSY7+CV/eBB5+MOMA+LRQOyKnI6uXN8BVZl0J4cjm/XmEw+n5hPh7MndSb6dKcgBu6tcWDzcNe1JzOHDOsl2/dkFXrix4mX9OWQDzps9cI8kB6BwPoT2hrBC2f9bw8cJqXDbRcfbKyEI4uq0nLvBJxfTruZP60MqO17FqqhB/L67tqdSg+mrbaZWjsYLVLysLXnr6KwtgegWoHZHtaDRw5aPK9tYFUFasbjwuzGUTHSGE/covLuPxb/dgMMAt/dtyTUyY2iFZzf8NVMo5LNt9jqJSJxqUfPAn2PSOsn39+8oCmK6m1yQIioLCTKmroyJJdIQQdueVnw9y5uIl2rbw4YXxMWqHY1VDOraifcWg5J/2OMmg5Kxj8MODyvbgBGXhS1fk5gFDEpTtv94DvXVKcYj6SaIjhLArfyans2RHKhoNvHVzLAHezj2mQ6vVcNsApVXn+11nVY7GAkoLlTWsSvOh3RC45iW1I1JX37uUasnZx+FQ3ev+CeuRRMcMLj4xTdRCfies40JBCc8s3QvAfUM7MLhjK5Ujso3r4pSyDttPZnMu55LK0TSDwQA/PQoZyeAfBjcvdp3Bx3Xx8leKCIJSV0feO2xOEp16uLkpRclKS0tVjkTYm8oK0ZcvYSGa5/ll+8kqKKVrmD+PX+s6Yzoig30YGN0SgwHHrqmz7RPY962yvMPNiyEgvMGnuISB/wA3Lzi3S1mtXdiUy651lZiYSGJiYr3LF7i7u+Pr60tmZiYeHh5otZIXujqDwUBRUREZGRkEBwcbk2HRfCv2p/Hb/jTctRrm3RKHt4drXdsJcRFsO5nNj3vOcf/wTmqH03ip2+D3Z5Xta/8F7f+mbjz2xL819L0DdixUBmhHD1U7IpciBQMbKDhUWlpKSkqK2QtNCtcQHBxMeHi409V1UUtecRnXzFtHel4JD43sxFNjuqsdks1lF5Yy8NWVlOsNrHp8BJ1a+6sdkvmqFwXsORFuWtSs9fmc0oXj8F4/wAAP/gVhPdWOyOGZWzDQZVt0zOXp6UmXLl2k+0oYeXh4SEuOBRkMBl5cfoD0vBKiW/nyyNVd1A5JFS39PBnWJYQ1hzP5Mekcj13TVe2QzKPXwff3VBUFvO49SXJq06oTxFwHycuVujrXvad2RC5DEh0zaLVaKfMvhJWsP5rF0t1ncdNqeP3GPi7XZVXddXERSqKz5xz/jO/iGC2Gm+ZDynplqQNXKwrYWAOmKYnOgeUw9g3w8FE7Ipcgg06EEKq5VKpj9vL9AEz5WzSDXGSWVV2uiQnH20NLSlYh+8/mqR1Ow87uhDVzlO2/v+WaRQEbo/2VSgHBklxlNXdhE5LoCCFUM3/lEU5eKKJNkDf/jHfNLqvq/L3cubqHUgV6eZKd19QpKYDvp4G+XBmXE3u72hHZP60Wet+sbO9dom4sLkQSHSGEKvadyTWuZfXKDb2cvjCgua6PVWrq/LT3HDq9Hc8VWTFTKYIXGAnj/y3jcswVe5tye2ylMohbWJ0kOkIImyvT6Xnq+73oDTAhNsLYiiFgRLfWBHq7k55XwraUbLXDqV3yctj9BaCBiR+BTwu1I3IcrbtBRF+lJWz/92pH4xIk0RFC2NzH609w8HweLXw9mD3BudeyaiwvdzfG9moDwI/2uPZV7ln48RFle+g/ocMwVcNxSH0qWnX2fq1uHC7CZROdxMREYmJiGDBggNqhCOFSjmcW8M6qowDMmhBDiL+XyhHZn8olIX7dd57Scjuq4aXXw7IHoDhHaZUY+azaETmmXjeC1h3O7YbMw2pH4/RcNtFJSEggOTmZ7du3qx2KEC5Drzcw8/u9lJbrGdmtNTfERaodkl0a3LEVrQO8yL1UxvojdjSOY/N7FVPJfWHSp+DuqXZEjsm/NXSOV7b3SKuOtblsoiOEsL0vt55i+8mL+Hq68coNvRyjTowK3LQaxvexs+6rc0mw6l/K9pjXIKSzquE4vD63Krd7v1FayoTVSKIjhLCJ/OIy3vhdaaZ/ekx32rbwVTki+3Z9RWvXn8npFJWWqxtMaRF8fx/oy6D7eLjibnXjcQbdxoJXEOSdgVMb1Y7GqUmiI4SwiQ/XHie/uJzOof7cNbi92uHYvdi2QbRv5culMh1/JqerG8wfz8GFoxDQRpZ4sBQPH+h5vbK9R2rqWJMkOkIIq9t7JoeP1is1c564thtarXxQNkSj0XBdRU2dH5NU7L46+LOy6jbAxAXg21K9WJxN5eyr5OVKq5mwCkl0hBBWZTAYeH7ZfnR6AxNiIxjTK1ztkBzG9RWzr9YdyeRioQoLC1+6CD9OV7b/9jB0HGn7GJxZuyEQ3A5K8+Hwr2pH47Qk0RFCWNUXW06x90wuPh5uUjOnkTqHBtCjTSDlegMrDqTZPoDfn1OSndY94KpZtn99Z6fVQu9blO3kZaqG4swk0RFCWE16XjGv/3YIgMev7So1c5qgcvbVr/vO2/aFD/4MSV+CRgtjX5Op5NbSfZxye3wNlJeoG4uTkkRHCGE1r/12iMJSHX3bBXPPlR3UDschjeutJDp/Hb9gu+6rSxfhlxnK9pX/lC4ra2rTF/zDoLQATsrsK2uQREcIYRV/Hc/ih91n0Wjgpet6ygDkJuoQ4kePNoHo9Ab+SLZR99Xvz0NBOoR0hRFP2+Y1XZVWC11HK9tHVqgbi5NyikQnOjqaPn36EBcXx6hRo9QORwiXV1ym47kf9gNwx6B29GkbrG5ADu7vvZUB3L/ss0Gic2wVJP0X0MB174OHt/Vf09V1HaPcHlkBBjtesd5BOUWiA/DXX3+RlJTEmjVr1A5FCJf30boTpGQVEhrgxVNjuqsdjsMzdl8dyyKnyIrdVyUF8NM/le1B/4B2g6z3WqJKx5Hg5gU5pyHjoNrROB2nSXSEEPbhzMUiPlh7DIAXxscQ6O2hckSOr2Nrf7qHB1CuN/CHNYsHrnoJck8rU56vesF6ryNMefpBxxHK9pHf1I3FCame6Kxfv54JEyYQERGBRqNh2bJlNY5JTEwkOjoab29vBg0axLZt20we12g0jBgxggEDBvDll1/aKHIhRG3m/HqQknI9gzu2NM4YEs33995Wnn11ajNs+1jZnvAOePlb53VE7SrH6RyWcTqWpnqiU1hYSGxsLImJibU+vmTJEmbMmMHs2bPZtWsXsbGxjB49moyMDOMxGzduZOfOnfz444/MmTOHvXv32ip8IUQ1m45l8eu+NLQaePG6nrJopwWNrUh0Nh3LIreozLInL7tUVRiw753Q6SrLnl80rHKczpntUJilbixORvVEZ+zYsbzyyitMnDix1sfnzZvHtGnTmDp1KjExMSxYsABfX18WLlxoPCYyUln8rk2bNowbN45du3bV+XolJSXk5eWZ/Aghmq9Mp+elnw4AcNfg9nQPD1Q5IufSOdSfbmEBlOkM/HnQwt1Xa1+DC8fAPxyufdWy5xbmCWoL4b0BAxz9Q+1onIrqiU59SktL2blzJ/Hx8cZ9Wq2W+Ph4Nm/eDCgtQvn5+QAUFBSwevVqevbsWec5586dS1BQkPEnKirKuv8IIVzEF5tPcSS9gBa+Hsy4ppva4TilykHJv1my++rcbvjrPWV7/DzwCbbcuUXjdB2r3B6WcTqWZNeJTlZWFjqdjrCwMJP9YWFhpKUp0yzT09MZOnQosbGxDB48mLvvvpsBAwbUec5nnnmG3Nxc409qaqpV/w1CuIKsghL+vfIIAE+O7k6QrwxAtoZxFdPMNxzNIq/YAt1X5aWwfDoYdNBzEnT/e/PPKZqusvvq+Grl/0ZYhLvaATRXx44d2bNnj9nHe3l54eXlRWJiIomJieh0OitGJ4RreHPFYfKLy+kVGcitA6SV1Fq6hAXQOdSfYxkFrD6YwQ19I5t3wk3zIX0/+LSEcW9aJEbRDBEVVZIL0uHURhkrZSF23aITEhKCm5sb6emm/dHp6emEhzdvBeSEhASSk5PZvn17s84jhKvbk5rDNzuVltGXruuJm1RAtqpxFau/N3v2VcZBWPeGsj32DfALaWZkotm0WuhyrbIts68sxq4THU9PT/r168eqVauM+/R6PatWrWLIkCEqRiaEANDrDcz+8QAGA0zqG0m/9i3VDsnpjauYsr/2SCYFJeVNO4lep3RZ6cuU7pLeN1kwQtEsUiXZ4lRPdAoKCkhKSiIpKQmAlJQUkpKSOH36NAAzZszgk08+4T//+Q8HDx7kwQcfpLCwkKlTpzbrdRMTE4mJial3PI8Qon5Ld58lKTUHP083Zo6VCsi20C0sgI4hfpSW61l9KKPhJ9Rmy4dwdgd4BcL4f4OUAbAfnUZVVEk+BZmH1I7GKaie6OzYsYO+ffvSt29fQEls+vbty6xZswC49dZbeeutt5g1axZxcXEkJSWxYsWKGgOUG0u6roRonrziMl77TXkjfuTqLoQGyppItqDRaBhbMSi5SbOvLhyH1a8o29f+CwIjLBidaDZPP+gwXNmWRT4tQvVEZ+TIkRgMhho/ixcvNh4zffp0Tp06RUlJCVu3bmXQIFl/RQi1vbvyKFkFJXQM8WPqlR3UDseljO2ldF+tOZxBUWkjuq/0elieAOWXlA/TKyZbKULRLJ0rSqqkrFc3DieheqKjFum6EqLpjmXks/ivkwDMmhCDp7vLvpWoomdEIO1a+lJcpmfNoUzzn7jtIzi9GTz9lZXJpcvKPkUPVW5PbwGdhatguyCXfXeSrishmsZgMPDST8mU6w3E9whjZLdQtUNyORqNxlg88Nf9ZnZf5Z2v6rK65mVo0d5K0YlmC40BnxZQVqQUdBTN4rKJjhCiaf5ITmfD0Sw83bW8ML6H2uG4rMrigWsOZXCp1Ix6YH/OgtICaDsA+jVvMoewMq0W2l+pbJ/coG4sTkASHSGE2YrLdPzr52QA7h/Wkfat/FSOyHX1jgwiMtiHolId64400H11fDXs+wbQKDVztPLWb/eihym3KZLoNJfL/rbLGB0hGu/j9Sc4c/ESbYK8eWhUJ7XDcWlK95UZxQNLC+HHR5XtgdMg8gobRCearUNFopO6VZaDaCaXTXRkjI4QjXPmYhEfrD0GwLPjeuDr6fAryDi8ynE6qw6mU1xWR/fVujcg9zQERcHVs20YnWiW1j3At5UyTufsTrWjcWgum+gIIRpn7q+HKC7TM6hDS8ZXVOcV6oqLCiYiyJvCUh3ra+u+Sk+Gze8r2+PeBC9/2wYomk6rrZp9JeN0mkUSHSFEg7aeuMAv+86j1cDsCT3RyLRku6DRaBhdsfbVqoOXVUnW6+GXGaAvh+7jodtYFSIUzWIcpyP1dJrDZRMdGaMjhHl0emU6OcD/DWpHTESgyhGJ6iqn9689koFeX21tpD1fKTVzPHxhzGsqRSeapbJC8pntUFasbiwOzGUTHRmjI4R5vtmRSvL5PAK93ZlxTTe1wxGXGdShJQHe7qTnlbDpeJay81KOMp0cYMTTEBylWnyiGUK6gn8YlBcra5OJJnHZREcI0bD0vGLeWKGsZ/XYNV1p6eepckTict4ebkyIVdar+m1/mrJz3etQlAUh3WDwQypGJ5pFo6kapyPTzJtMEh0hRJ2e+2EfF4vK6BkRyJ2DpZKuvRrdUxmnszI5Hf2JDbB1gfLA2NfAXZJTh2Ycp7NO3TgcmCQ6QohaJZ/LY+XBDLQaeOe2ODzc5O3CXg3u2BI/Tzcy8y9x6bcXwKCH2P+DTlepHZporo4jldsz26EkX9VQHJW8cwkhatDrDbxW0WU1tncbOocGqByRqI+XuxsjurXmHrcV+GXuVgYgx0vNHKfQsgMEtVNmz0k9nSZx2URHZl0JUbc1hzNYfyQTd62GB0dIBWRHMLHdJZ5yX6LcufYVCAhXNyBhOVEVn1OpMnmmKVw20ZFZV0LUrrhMxxsrDgNw15D29IoMUjkiYY4R5xfjpSljva43pzvcpnY4wpLaDlRuz8jnVVO4bKIjhKjd/JVHOZyeT4i/J/cP76h2OMIcJ9bimfwtAPPLb+Sn+ta+Eo6nbUWLzpntYDDUf6yoQRIdIYTRhYISPtt4AoA5E3vTJshH5YhEg4pzYek/ADjW/lZ2GbqydNcZDPKB6DzCe4ObF1zKhuwTakfjcCTREUIYvfXHEcp0BnpHBnFtTxnj4RB2fQEFaRAQQfh1L+LlruV4ZiEHz8sMHafh7gltYpVt6b5qNEl0hBAAHM8s4OvtpwF45OouKkcjzJKfrqxODjDyafxbRTCsS2sA/kxOVzEwYXHVu69Eo7hsoiOzroQwNffXgxgMEN8jlGtiwtQOR5jj92egJBci+kLfuwC4tqfyf/dHcpqakQlLi5JEp6lcNtGRWVdCVFlzKIOVBzNw12qYOba72uEIcxxbBfu/B40Wxs8HrRsA8T3C0GrgwLk8zuZcUjdGYTmVLTpp+6G0SN1YHIzLJjpCCEVBSTkv/XQAgHuGdpDigI5Ar4cVzyjbA/8BEXHGh1r6edK/fUsAVh2U7iunERgJ/uFg0MH5PWpH41Ak0RHCxb2x4hAnLxQRGuDFw1d1VjscYY4tH0DWYfAKhFHP1Hg4PiYUkHE6TkWjgbb9lW1ZybxRJNERwoVl5BezdNdZAF6d2JsAbw+VIxINKsiENXOU7VHPgXfNgo7xPZRxOltOXCCvuMyW0Qlriuyn3J6RRKcxJNERwoW99ushCkrKad/Kl1HdWqsdjjDHny9AWaEyAHnQP2o9pGNrfzq19qNMZ2D9kUwbByisxtiiI2teNYYkOkK4qJNZhfy45xwA/741DndZndz+nVgLe75SBiCPnqN0Z9QhvmLm3ErpvnIeEX0BDeSmKqUFhFnknU0IF2QwGHj0692U6w0M6diKK9q1UDsk0RBdOfw5S9keMA3a/63ew6+p6L5afSiDMp3e2tEJW/AKgNAeyraM0zGbJDpCuKCNx7LYcyYXTzctcyb1VjscYY6/3lVm23gFwvAnGzy8b7sWtPLzJK+4nO0ns20QoLAJ4zgdKY1iLkl0hHAxabnFTP/fbgBu7NeWDiF+KkckGnRuN6x6SdmOfxH8Gx5P5abVcFV3ZfbVyuQMKwYnbKpynI4MSDabyyY6UhlZuKoF646Te6mMXpGBzJ4Qo3Y4whzbP1NuQ7pB3B1mP61ynM6fB9NkkU9nEVmR6JzbDXqdurE4CJdNdKQysnBF5To9P+9VBiA/fm03vD3cVI5INCjzMOz9RtkePw88vM1+6rAuIXi6a0nNvsTRjAIrBShsKrQHePhBaYHyuyEa5LKJjhCu6NONKWQVlNLC14OhnUPUDkc0RK+DpfeDrgQ6joT2Vzbq6b6e7sb/Zyke6CS0bhB5hbItA5LNIomOEC7iQkEJ8/48AsAz43rgIdPJ7d+Gt+F8EngGwKRP6p1OXpfK4oGS6DgRKRzYKPJOJ4SLeH3FIUrL9fRpG8TN/dqqHY5oSFE2bHpH2R73JviHNuk08T2U5yWl5pCRX2yp6ISaKlt0zu1WNw4HIYmOEC5g49EsvtlxBq0GHrumK5omtAwIGzIY4IcHlHEY4X0g9rYmnyo00JvYqGAAVh+U2VdOoU2ccptxEMpLVA3FEUiiI4STKynX8dqKgwDcNbg9o7o1rWVA2NCJNXD0d9B6wLi3mtRlVd01PWSRT6cS3A58WoC+DNIPqB2N3ZNERwgn982OM+w/m4e3h5Z7h3ZUOxzREIMB1r2hbA+4D9oNavYpK6eZbzyWRXGZTEl2eBpNVavO+SQ1I3EIkugI4cQMBgMfrz8OwGPxXWnXylfliESDTm6E05vBzQuufNQip+wWFkCbIG9KyvVSJdlZRMQpt+eS1IzCIUiiI4QTm7X8AKnZl/Bw03CTDEC2fxdPwo8PK9tX3A2BbSxyWo1Gw7AuyjTzDUezLHJOoTJp0TGbJDpCOKldpy/yxZZTADw4sjOt/L1Ujkg06MeH4WIKeAXBsBkWPfWwLsqyEWsPy4Bkp1DZopOeDOWlqoZi75wm0SkqKqJ9+/Y88cQTaocihOp0egOzlu8H4KZ+bZlxTVeVIxINyj4BKeuV7Vv+A4ERFj398C6t8XDTcCS9gMNp+RY9t1BBcHvwDlYGJGckqx2NXXOaROfVV19l8ODBaochhF34fpcyADnA252ZY7urHY4wx/q3ldtOV0OnURY/fZCvByO6Kq06Kw/K7CuHp9FUtepI91W9nCLROXr0KIcOHWLs2LFqhyKE6pQByCcAmD6qMyHSZWX/so5B0n+VbQsNQK7N8IpEZ6OM03EO4X2U2/N71Y3Dzqme6Kxfv54JEyYQERGBRqNh2bJlNY5JTEwkOjoab29vBg0axLZt20wef+KJJ5g7d66NIhbCvq09ksmxjAL8vdz5v0Ht1A5HNMRggNX/Ura7XAsdR1jtpSrXvdp56iKXSmWaucOrTHTS96sbh51rVKJTVlaGu7s7+/db7qIWFhYSGxtLYmJirY8vWbKEGTNmMHv2bHbt2kVsbCyjR48mI0MZULd8+XK6du1K164yBkGICwUlPP2d8u3utgFRBHh7qByRaNCGtyF5GWjdYeQzVn2pDiF+RAb7UKrTszXlglVfS9hAeC/lNv0A6PXqxmLH3BtzsIeHB+3atUOns9w3gbFjx9bb5TRv3jymTZvG1KlTAViwYAG//PILCxcuZObMmWzZsoWvv/6ab7/9loKCAsrKyggMDGTWrFm1nq+kpISSkqqS2Xl5eRb7twihts83nyIjv4QAb3fuGybFAe1edgqsrWiNjn+pag0jK9FoNAztHMKSHalsPJrFSKmS7dhadVHqLZUWQM5JaCl/87VpdNfVc889x7PPPkt2tvWLTpWWlrJz507i4+ON+7RaLfHx8WzevBmAuXPnkpqaysmTJ3nrrbeYNm1anUlO5fFBQUHGn6ioKKv/O4SwhTKdnq+3nwbg1Ym9CQ/yVjki0aB1r4O+HDpdBX+bbpOXHFpRT2fjMRmn4/Dc3CG0YrJBmnRf1aXRic7777/P+vXriYiIoFu3blxxxRUmP5aUlZWFTqcjLCzMZH9YWBhpaWlNOuczzzxDbm6u8Sc1NdUSoQqhul/3nSc9r4QQf0/G9AxXOxzRkMzDsHeJsn3V8zZ72Ss7h6DRwKG0fFnN3BmE91Zu0/apG4cda1TXFcANN9xghTAsY8qUKQ0e4+XlhZeXF4mJiSQmJlq0G04ItZy6UMizS5U3ulv6R+Hprvo8A9GQNXPAoIfu4yGyn81etqWfJz0jAtl/No9Nx7KY2FcqZju0sIpERwYk16nRic7s2bOtEUetQkJCcHNzIz3dtOZDeno64eHN+8aakJBAQkICeXl5BAUFNetcQqippFzH+Hc3UliqIy4qmIdGdVY7JNGQoyuVAchoYNSzNn/5oZ1bs/9sHhuOSqLj8IwDkiXRqUuTvvbl5OTw6aef8swzzxjH6uzatYuzZ89aNDhPT0/69evHqlWrjPv0ej2rVq1iyJAhFn0tIRzVV1tPk19SDsBj13TF36vR31+ELWWnwJI7lO2B0yCsp81DqFz3auPRLAwGg81fX1hQ64oxOjmpUFqkbix2qtHviHv37iU+Pp6goCBOnjzJtGnTaNmyJUuXLuX06dN8/vnnjTpfQUEBx44dM95PSUkhKSmJli1b0q5dO2bMmMHkyZPp378/AwcOZP78+RQWFhpnYTWVdF0JZ2AwGPgh6Zzx/t86tVIxGtEggwG+uQvKi8HDF4Y/qUoY/dq3wMtdS0Z+CUczCugaFqBKHMIC/ELAtxUUXYALR6FNrNoR2Z1Gt+jMmDGDKVOmcPToUby9q2Z1jBs3jvXr1zc6gB07dtC3b1/69u1rPH/fvn2NM6duvfVW3nrrLWbNmkVcXBxJSUmsWLGixgDlxkpISCA5OZnt27c36zxCqOnnvefZk5qDVgOrHh+Bh5uMzbFrZ3cpg0bdvGDqb+CvzvRubw83BnZoCchq5k6hslUn87C6cdipRrfobN++nY8++qjG/sjIyCbNhBo5cmSDTafTp09n+nTbTL0UwlGU6/R8ujEFgHuu7ECn1v4qRyQatOd/ym3MdVXrFKlkWJcQNhzNYuPRTO4d2kHVWEQzhXSFU5sg85DakdilRn/98/LyqrXI3pEjR2jdurVFgrKFxMREYmJiGDBggNqhCNEkH284wZ7UHPy93Ll3mHxQ2b1Tm2HHImU79nZ1Y0EZkAywNSWb0nKpquvQKlt0MiTRqU2jE53rrruOl19+mbKyMkCptHn69GmefvppbrzxRosHaC3SdSUcmV5v4IvNpwCYNT6GNkE+KkckGvTnLDDooNdNSoFAlXUPDyDE35OiUh27Tl9UOxzRHJUzr2QV81o1OtF5++23KSgoIDQ0lEuXLjFixAg6d+5MQEAAr776qjViFEJUYzAYePaHfZzPLSbAy53r4iLUDkk0JPMwnNkGaGD0q6DRqB0RWq2GKztXzb4SDqxNLKCBvLOQn97g4a6m0WN0goKC+PPPP9m4cSN79+6loKCAK664wmSZBiGE9fznr5N8vV2p6P38+B54e7ipHJGoV3EufP1/yna3cRBgP1Wrr+wcwvKkc2w4lsUTo7upHY5oKq8AaNVZmXWVvh8CmjdZx9k0ueDG0KFDGTp0qCVjsSmZXi4cUUm5jhd/SgaUqeS3DminckSiXno9fPA3yDsDPi1g9CtqR2Sisp7OvjM55BaVEeQrq907rNbdlEQn6wh0vlrtaOxKk+airlq1ivHjx9OpUyc6derE+PHjWblypaVjsyoZoyMc0cg31xq3Z02IUS8QYZ7UrUqSA3D9B3a3unSbIB86tfZDb4C/jkv3lUMzTjGXAcmXa3Si88EHHzBmzBgCAgJ49NFHefTRRwkMDGTcuHEkJiZaI0YhBHCxsJTzucoijNfEhNE9PFDliESDDixVbmNvh+7j1I2lDsO6KLOvNshq5o5NaunUqdFdV3PmzOHf//63SV2bRx55hCuvvJI5c+aQkJBg0QCFEIp//Zxs3H77Fql+avfO74XtnynbPa5TN5Z6XNk5hMV/nWR7SrbaoYjmaF0xxirzkFKB2w4GvNuLRrfo5OTkMGbMmBr7r732WnJzcy0SlC1IHR3hSFYfSmfpbmUtubdvjiXQW8ZS2DVdOax4RplO3mEEdB2tdkR16tNWWdT4eGYBl0plzKLDCukCaODSRSjMVDsau9KkOjo//PBDjf3Lly9n/PjxFgnKFmSMjnAk81ceNW5P7BupYiTCLFs/hFMbwcMPxv8btPY7My40wIsQf0/0BjiUVrMYrHAQHj7QIlrZzjioaij2xqyuq3fffde4HRMTw6uvvsratWuNK4hv2bKFTZs28fjjj1snSiFc2LGMfPaeUVpL507qjVYrTdJ2rSgb1sxVtse+Bq06qRtPAzQaDT0jglh3JJP95/Lo266F2iGJpgrtARdTlO6rjiPUjsZumJXo/Pvf/za536JFC5KTk0lOrhozEBwczMKFC3n++ectG6EQLqykXMe1/1YWy+3RJpDbBkSpHJFo0P7voawQwnpB3J1qR2OWnhGBrDuSyfaUbO4a3F7tcERThfaAw79Ki85lzEp0UlJSrB2HEKIWCzeeRF+x5u1j8V3QyABD+1ZaCBvmKdt97wStY6wmP7JbKB+sPc7Pe8/x/N97EBrorXZIoila91BuZYq5Ccf4KxTCBRkMBj7ZcMJ4f0B0SxWjEWb59SnIPwdB7aDfVLWjMdvADi3pEuqP3gD7zjrOpBJxmdDKxT0PKjOvBNCE6eUGg4HvvvuONWvWkJGRgV5vuurt0qVLLRacNUllZGHvfj+QRnZhKaDMtGrh56lyRKJeWUch6b/K9oR/g4djtYr0jgziaEYBB87lcXUPWULAIbXqAhotFOdAfhoEtlE7IrvQ6Badf/7zn9x1112kpKTg7+9PUFCQyY+jkFlXwp4ZDAYWbjoJwP3DO3Jjv7bqBiQatmORctt1DHR2vLX/YiKUApTJ52TmlcPy8IaWFYPfM2WcTqVGt+h88cUXLF26lHHj7LPKpxCOrqRcR89Zv1OuN+DprpXBoY7gyO+wdYGy3f8edWNpopg2FYnOeUl0HFpod2XNq4xD0OkqtaOxC41u0QkKCqJjR/tar0UIZ/LHgXTKK0Yg3ze0A1EtfVWOSNSrIAO+v6+qOGCXa9WOqEliIgLRaOB0dhFH0/PVDkc0lXFAsrToVGp0ovPiiy/y0ksvcenSJWvEI4RLKy7T8eovVW9Q/xhh3zVYXJ7BAN/fCyV54NsKbv/aYUvvB/t6El8xNmdZ0lmVoxFNZhyQLDOvKjW66+qWW27hq6++IjQ0lOjoaDw8TEvR79q1y2LBCeFqNh7NIi2vGG8PLeueHEWQjyz1YNcuHIeU9aB1hzu+BU/Hbn0b0bU1fyans/+sdF85rOpTzGXNK6AJic7kyZPZuXMnd955J2FhYQ5b10NmXQl79J/NJwG4bUA7wqSWiX3T62FtRQXktgMhsp+68VhAr0hlQsmu0xcpLCnHz6vRHxFCba06K4l3SR7knYUgmcjQ6N/iX375hd9//52hQ4daIx6bSUhIICEhgby8PIeaLSac16/7zrPhaBZuWg33XNlB7XBEQza8Bfu/U7ajHfv9sFKviEDatfTldHYRP+89x60D2qkdkmgsd08l2ck8pHRfSaLT+DE6UVFRBAYGWiMWIVzW5uMXeOhLpdt3RNfWtGvl2F0gTk+vh82JyrZvCAxJUDceC3F30zKut1J7JSlVCgc6rJYVE4ZyTqkbh51odKLz9ttv89RTT3Hy5EkrhCOEa1p9KN24/eBIGYBs9zIPKkXZPHzh8cPgE6x2RBbTp63Swr3vbI66gYimC6goFJh/Xt047ESju67uvPNOioqK6NSpE76+vjUGI2dnZ1ssOCFcwcXCUpZsTwXgrsHtZakHR7DhbeW27QBwc65xLF3D/AE4lVWEwWBw2HGYLi0wQrnNO6duHHai0X+h8+fPt0IYQriu+SuPkFdcTligF0+O6aZ2OKIhJ9YqK5QDXPW8qqFYQ9sWSrdpfkk5OUVlsvSII5JEx0STZl0JISxj+8ls/rNZ6Ud/9OquBHrLdHK7ptfBT48q212uhaiB6sZjBd4eboQGeJGRX8Lp7CJJdBxRi2jlVqaYA01IdE6fPl3v4+3aySh9Icy1YO1xAK7uHsptA6JUjkY06OifcPEkeAbAxI/UjsZquoT5k5Ffwo5TF4mNClY7HNFYbeJA46aM0clNhWDX/lxudKITHR1db5+t1KURwjy/7TvPqkMZADw9tjtarWt/67J7er0ypRyg/xTwdd6xVNf0CGPTsQv8tu889w6VUgcOx9MXQmMgfR+kH5BEp7FP2L17t8n9srIydu/ezbx583j11VctFpgQzsxgMPDu6mMA3NSvLV3DAlSOSDRoSyKc2Q4aLQy4T+1orGpMrza8+FMyO05dJD2vWIpXOqJWHZVEJztF7UhU1+hEJzY2tsa+/v37ExERwZtvvsmkSZMsEpi1SWVkoaadpy5y8Hwe3h5aXvh7jNrhiIac+qtqptWwx6vGQDip8CBvrmgXzK7TOfyZnM6dg9urHZJorMpaOhcl0Wl0HZ26dOvWje3bt1vqdFaXkJBAcnKyQ8UsnMemYxcAuDYmnCBfGYBs1wwG+HYqXLqoFAe88lG1I7KJqysW+PzreJbKkYgmaVHR5Zh9Qt047ECjW3Ty8kwXezMYDJw/f54XX3yRLl26WCwwIZzZb/uVQl69IqXKuF0zGOCr26AgTbk/bRV4uUY3Y0yE8rt5PKNQ5UhEk7SsTHSkRafRiU5wcHCNwcgGg4GoqCi+/vpriwUmhDPS6w089f1eDqXlA9BFxubYtzPb4cgKZbvLaKfvsqquU4hSODDlQiE6vQE3GSzvWIzLQJxWyiJo3dSNR0WNTnTWrFljcl+r1dK6dWs6d+6Mu7tzVQgVwtI+Wn+C73aeAcDTXcuQjq1UjkjU6+zOig0NXJ+oaii2FtnCB093LaXles7lXCKqpay/5lACIsDNC3QlyppXlYmPC2p0ZjJixAhrxCGEU9PrDfxn80leX3HIuG/9k6Pw9nDdb1l27+ifsGKmsj3iKfBvrW48Nuam1RDdypcj6QUczyyQRMfRaLUQFgPndsO5JEl0Guvo0aOsWbOGjIwM9Hq9yWOzZs2ySGBCOJOlu8/y0k/Jxvvbnrua0ACZsmu3DAZY9XLV/c7XqBeLijqG+HMkvYATmYWMlNVJHE9kPyXRObsTejnGjGhraHSi88knn/Dggw8SEhJCeHi4yXgdjUYjiY4QtVh1sGp18iEdW0mSY+/2fA1pe5Xtf2yANn3UjUclXcP8WXEAtqVkc48UDnQ8oT2U2wvH1Y1DZY1OdF555RVeffVVnn76aWvEI4TTKS3X89t+ZdbOP4Z35JlxPVSOSDToRMVYxLg7XTbJAaVw4Lurj7H6UAYl5Tq83KWr1aFUTjF38Vo6ja6jc/HiRW6++WZrxCKEU/p6e9X6cDf2a6tiJMIsqdtg7xJlu/vf1Y1FZT3aBBDk40GpTs/R9AK1wxGNVTnF/OJJZQkTF9XoROfmm2/mjz/+sEYsQjil3w8orTldQv1lqQd7ZzDA9/cq2xo3iIhTNRy1aTQaYtoo9XQOns9r4Ghhd4KiAA2UF0PRBbWjUU2ju646d+7MCy+8wJYtW+jduzceHqZVXR955BGLBWeOnJwc4uPjKS8vp7y8nEcffZRp06bZNAYh6pJdWGqsgvzhnVeoHI1oUNo+pe4IwN3LITBC3XjsQI82gWw+cYFkSXQcj5uHsvhs0QUozHC5mYOVGp3ofPzxx/j7+7Nu3TrWrVtn8phGo7F5ohMQEMD69evx9fWlsLCQXr16MWnSJFq1kvokQl3FZTqGzF0FKFN1O1YUYBN2KvMwfDRM2e50NXQYpm48dqKyQnLyOUl0HJJ/mJLoFKRDWE+1o1FFoxOdlBT7GtTk5uaGr69S36GkpASDwYDBYFA5KiHg5Z+TKSlX+sXvGtwerVSWtW/r3qja7jhStTDsTc+KRGff2VwZkOyI/EMhIxkKMtSORDUWW9SzqdavX8+ECROIiIhAo9GwbNmyGsckJiYSHR2Nt7c3gwYNYtu2bSaP5+TkEBsbS9u2bXnyyScJCQmxUfRC1O7MxSL+t1XpArmiXTAvXuea36QcSupW5TagDfSfqm4sdqRbWAChAV4UlerYlpKtdjiisfxCldv8NHXjUJHqiU5hYSGxsbEkJtZeXn3JkiXMmDGD2bNns2vXLmJjYxk9ejQZGVXZaXBwMHv27CElJYX//e9/pKen13ouIWzl8W/2GLdHdA1VMRJhlgM/QG4qoIGErS6zcKc5tFoNQzopQwH2nc1VORrRaC3aK7dZR9WNQ0WqJzpjx47llVdeYeLEibU+Pm/ePKZNm8bUqVOJiYlhwYIF+Pr6snDhwhrHhoWFERsby4YNG+p8vZKSEvLy8kx+hLCk91YdZWu1b773D3fd0usOY9M7ym2fW8A7SN1Y7FDlbMEjFYvRCgcS0Ve5PbdL3ThUpHqiU5/S0lJ27txJfHy8cZ9WqyU+Pp7NmzcDkJ6eTn6+8seXm5vL+vXr6dat7lrlc+fOJSgoyPgTFRVl3X+EcCl6vYEF66qqkO6ZdS0+njKmwa7t/VYpkw8w/El1Y7FTlYnOIUl0HE+bWOU28zDoytSNRSV2nehkZWWh0+kICwsz2R8WFkZamtLfeOrUKYYNG0ZsbCzDhg3j4Ycfpnfv3nWe85lnniE3N9f4k5qaatV/g3AtKRcKKSzVAbDrhWsI8vVo4BlCVUlfwdL7lO1OV0NIF3XjsVN92iqtXEfS88kvds0PS4cVEAHuPmDQVZVOcDFNWtRzw4YNfPTRRxw/fpzvvvuOyMhIvvjiCzp06MDQoUMtHWO9Bg4cSFJSktnHe3l54eXlRWJiIomJieh0OusFJ1zKmYtFzFq+H4C4qGBa+nmqHJGoV2kRLHug6v6A+9SLxc6FBXoTGezD2ZxL7Duby986yYQPh6HVKhWSM5IhOwVadVI7IptrdIvO999/z+jRo/Hx8WH37t2UlJQASrfRnDlzLBpcSEgIbm5uNQYXp6enEx4e3qxzJyQkkJyczPbt25t1HiEqvfxTsrE4YOWUXGHHTqyt2g5qB9FXqhaKI+gSptSBOn2hSOVIRKO1rBgnmHVY3ThU0uhE55VXXmHBggV88sknJlWRr7zySnbtsuxgJ09PT/r168eqVauM+/R6PatWrWLIkCEWfS0hmmvbyaoByO1a+qoYiTBL6hbltu+d8Ng+GYTcgKgWyu906kVJdBxOmzjl9qxrDkhudNfV4cOHGT58eI39QUFB5OTkNDqAgoICjh07ZryfkpJCUlISLVu2pF27dsyYMYPJkyfTv39/Bg4cyPz58yksLGTq1ObVuZCuK2FJaw5nkFOkjF1w12q4Lk6WDrBr+Wmw7RNlu+0AdWNxEJXJ+ylp0XE8kRUzr84nqRqGWhqd6ISHh3Ps2DGio6NN9m/cuJGOHRs/jXbHjh2MGjXKeH/GjBkATJ48mcWLF3PrrbeSmZnJrFmzSEtLIy4ujhUrVtQYoNxYCQkJJCQkkJeXR1CQfJMTTZd7qYz7P98BQGSwD5tmXqVyRKJe+WnwdsXMzJBu0PsWdeNxEJ1Dla4rWQrCAbWsGJeTe0ZZuFbjWlXaG53oTJs2jUcffZSFCxei0Wg4d+4cmzdv5oknnuCFF15odAAjR45scMmG6dOnM3369EafWwhbWHs4gzKdAU93LYumSuuAXSspqEpyAK6eBZ7SzWiOvu2CATiRVciFghJa+XupG5AwX+XitJWrmPu51mDyRic6M2fORK/Xc/XVV1NUVMTw4cPx8vLiiSee4OGHH7ZGjFYhXVfCEk5kFjBr+QEA7h/W0VhvRNipv941vd9trDpxOKBgX086hvhxIquQ5PN5DOvimithOyR3L2Vxz4J0pQK4iyU6jR6MrNFoeO6558jOzmb//v1s2bKFzMxM/vWvf1kjPquRWVfCEt764zC5l8oI8ffivmEd1A5HNGTvN1XbT58CrRRzbIwebZTZhAfPS/eVwwmMVG5zz6obhwqaXDDQ09OTmJgYBg4ciL+/vyVjEsIhHM8s4Nd9SuHKTyf3J9hX6ubYtbxzcDFF2b5/HfgEqxqOI+rRRmmx3HJCFvd0OEEViU6e6yU6ZnVdTZo0yewTLl26tMnBCOFIbvt4i3E7tq0MaLdrWUfh/f7KdptYiIhTNRxHNbZ3G9764whrD2dQUFKOv1eTas4KNQRVLHeU63qrAZj1W+qMs5JkjI5ojjKdnsx8pVjmyG6t0bjYLAaHs7pa1/oQxxlLaG86tfbH38udgpJy0vOK8W8trfkOw4W7rsxKdBYtWmTtOGxOppeL5vioYuFOL3ctCyfLTCu7lnUUkpcr2+P/DX1uVjceBxfi70lBSTlZ+SV0kkTHcbhw11Wjx+ikpKRw9OjRGvuPHj3KyZMnLRGTEHYtI6+Yd1YpfwMJozqj1Uprjt3KT6/qsuoyGvrfo248TiCkYlp5VkGpypGIRjF2XUmi06ApU6bw119/1di/detWpkyZYomYhLBbxWU67vpsG2U6A7Ftg3j4qs5qhyTqUpQNb3etut/7JvVicSKVic7awxkqRyIapbLrKv8c6MrVjcXGGp3o7N69myuvrLn43eDBgxu1irjaEhMTiYmJYcAA6XYQ5vvvllMcTs8H4Ob+UTI2x559Msr0fseRqoThbIJ8lDUOv915psFir8KO+IeC1h0MeihIUzsam2pSHZ38/Pwa+3Nzcx1qYK/U0RGNVVym45VfDgJw1+D23DGoncoRiTqd3AgXT1bd7zBCeaMXzXbn4PbG7aJSx3nPd3lat6oKyS7WfdXoRGf48OHMnTvXJKnR6XTMnTuXoUOHWjQ4IexFUWk53V9YYbz/9z5tpDXHXhkMVYOPKw2cpk4sTqhXZCBe7spHR3ahjNNxKIFtlVsXm2Le6CIIr7/+OsOHD6dbt24MGzYMgA0bNpCXl8fq1astHqAQ9mDFftOm3rioYHUCEQ1b/yZs+9h0n3ewKqE4I41GQwtfT9LyiskpKiOqpdoRCbO56MyrRrfoxMTEsHfvXm655RYyMjLIz8/n7rvv5tChQ/Tq1csaMQqhusqS9yH+XqTMHYe3hywdYJf0eljzas39Pi1sH4sTa+GnVAHPLpIWHYfiorV0mlTWMiIigjlz5lg6FpuSgoHCHGU6PTct2Mye1BwAnri2q3RZ2bOjv5ve7xwP7YZAuHwJs6SWfsqA5IvSdeVYWlSMr8o6rG4cNmZWorN371569eqFVqtl79699R7bp08fiwRmbVIwUJjjSHq+McmJCPJmTK9wdQMSdSsvga9uq7rv0wImfgx+rdSLyUmFBngD8M8lSfRr34Kolr4qRyTMEjVYuT29FcpLwd011uczK9GJi4sjLS2N0NBQ4uLi0Gg0tU4r1Gg00kIinMru0znG7d/+Odw4tVbYoTM7qrZ9WkDCNklyrGR0z3B+2K10f7z6y0EW3NVP5YiEWVp3V8arFedAxgGI6Kt2RDZhVqKTkpJC69atjdtCuILcS2U8v2w/ALf2j5Ikx54ZDLB4XNX9e/6Q6eRWdG1MGH/r1Iq/jl/gUFqe2uEIc2m10KYPpKyHtP0uk+iYNRi5ffv2xnEJp06dIjIykvbt25v8REZGcurUKasGK4Qtfbujagrm9X0jVIxENGjl7KrtuDuhdde6jxXNptVqmDOxNwAZ+SVSONCRhCn/b2QcVDcOG2r0rKtRo0aRnZ1dY39ubi6jRo2q5RlCOJ5f9503Fgd8ekx3/tYpROWIRL32LKnajp9d93HCYkIDlaUgikp1FJS41pICDi2gYpxhUZa6cdhQoxMdg8FQ66yTCxcu4OfnZ5GghFDbQ1/uMm7HRslgdbtmMFSVtJ/yq3RZ2YivpzsBXsroh/S8EpWjEWarLLVwKUfVMGzJ7OnlkyZNApQBx1OmTMHLy8v4mE6nY+/evfztb3+zfIRWItPLRV0uXVbWvkd4oEqRCLNUH4TcVtaus6XQQC/yM8vJyCumc6i/2uEIc/gEK7eXLqoahi2ZnehUTsE2GAwEBATg4+NjfMzT05PBgwczbZrjlFmX6eWiLsczC0zuVxZHE3bq8+urtl1kuqy9CAv05nhmIen5xWqHIsxlbNGRRKeGRYsWGQecvffee/j7S/YunIvBYOCF5fv575bTaocizJWyHsoKle2WHdWNxQWFBSr1dDKk68pxuGCi06gxOgaDgS+//JLz589bKx4hVHPm4qUaSc6iKdIVYrf0Ovj1qar7d/+oXiwuqnJAsozRcSCViU5xjjK+zQU0KtHRarV06dKFCxcuWCseIVSRXVjKnZ9tNdm3+ZmrGNVdBrbapUs58OnVkHkQ3Lzg8SMQHKV2VC6nskKydF05EN8Q0LqDvhzyzqkdjU00etbVa6+9xpNPPsn+/futEY8Qqvhq22lOXSgy3p/UN5I2QT71PEOo6sub4dxuZXvcmxAQpm48LiqsokUnI08SHYfh7gmtOivb6QfUjcVGGr2o5913301RURGxsbF4enqaDEoGaq2xI4S9O55RNQD55n5tefPmWBWjEfXS66uSnB7XQd+71I3HhVWO0ZGuKwcTGgOZh5QW0a7Xqh2N1TU60Zk/f74VwhBCPeU6PT/vVcad+Xu5kzCqs8oRiXoVZoK+DDRauGmhUtZeqCKsousqI7+4zhprwg4FRSq3BRnqxmEjjU50Jk+ebI04hFDN3N8OUarTA/D1/YOJDpHCl3arIBPerljeIaANuMn6Y2qqHIxcXKYnr7hc1oNzFH7K2pUUZqobh40066tQcXExeXl5Jj9COJKdp7L5bKOyUO2DIzvRM0KKA9qtQ7/CW9Va2wIj1YtFAODt4UbrACXZ+XmvawxsdQrGRMc1loFodKJTWFjI9OnTCQ0Nxc/PjxYtWpj8OIrExERiYmIYMECmD7uiwpJyvt2Ryo0fbjbue+SqLtL0bs++vt30vgxAtgsjuyofms/9sJ+Scqk07xCkRad+Tz31FKtXr+bDDz/Ey8uLTz/9lJdeeomIiAg+//xza8RoFQkJCSQnJ7N9+3a1QxEquP+LHTz53V7j/X7tW+Dj6aZiRKLRAmRFeXvw4MhOxu30XBmU7BB8Wym3kujU7qeffuKDDz7gxhtvxN3dnWHDhvH8888zZ84cvvzyS2vEKIRFGQwGNh2rqgV1fVwEn98zUMWIRIPy003v+4bAsBnqxCJMdGztT/tWvgCcz72kcjTCLC2ildv88y6xuGejE53s7Gw6dlRKrQcGBhqnkw8dOpT169dbNjohrKDoskU7n7i2G35ejR6XL2zpp0dN7z91HALC1YlF1BBeMc08TerpOAbfllXJTmWpBifW6ESnY8eOpKQogze7d+/ON998AygtPcHBwRYNTghruFhUatz+5h9DiGrpq2I0wixHfqvavn+tamGI2rUJqkh0ciXRcRhhvZTbC8fUjcMGGp3oTJ06lT179gAwc+ZMEhMT8fb25rHHHuPJJ5+0eIBCWNr7q5U/7CAfDwZ2aKlyNKJe53bDi0Gm+2S2ld1p6afMvLpYVKZyJMJslQOSi5x/SadGt9c/9thjxu34+HgOHTrEzp076dy5M3369LFocEJYWn5xGV9vTwUg95K8Kdu1nNPw8cia+30cZ3anqwjwVj5K8ovlb8phuNDMK7MTHb1ez5tvvsmPP/5IaWkpV199NbNnz6Z9+/a0b9/emjEKYTGbjrlG3QinsPu/NfdN+UWKBNqhwIpCgXnF5SpHIszmQomO2V1Xr776Ks8++yz+/v5ERkbyzjvvkJCQYM3YhLC4A+eqilp2bC0VkO1aca7p/eFPQfRQdWIR9ZIWHQfkF6LcukDRQLMTnc8//5wPPviA33//nWXLlvHTTz/x5ZdfotfrrRmfEBZz6kIh71WMz+kQ4sd/psqUcrul18GR3033dRurTiyiQYHeFS060h3sOKRFp6bTp08zbtw44/34+Hg0Gg3nzknZb2H/dp7KZsSba433F00ZILOt7NnpzXAxper+iKch8gr14hH1Cqxo0dl1OofCEum+cgiS6NRUXl6Ot7e3yT4PDw/KytTN4FNTUxk5ciQxMTH06dOHb7/9VtV4hH269z87jNsT+0bKwp32bv2bym338fDAJhj1rLrxiHoFeFeNm/pq22kVIxFmq0x0Ll0EnXO3xJk9GNlgMDBlyhS8vLyM+4qLi3nggQfw86v60Fi6dKllI2yAu7s78+fPJy4ujrS0NPr168e4ceNMYhKuq7hMx+8H0sipNu313qEdVIxINKgoG06sVbaHPQ7hvVQNRzQsIrjqS3BWQWk9Rwq74dMCNFow6JUp5k5cgNPsRGfy5Mk19t15550WDaYp2rRpQ5s2bQAIDw8nJCSE7OxsSXQEAI9/s4df9p032de2hY9K0Yh66XXw17tVg5B9Q6S7ykG08vdiaOcQNh7LorRcxm06BK1W+RsrzFC6ryTRgUWLFlklgPXr1/Pmm2+yc+dOzp8/zw8//MANN9xgckxiYiJvvvkmaWlpxMbG8t577zFwYM2BpDt37kSn0xEVFWWVWIVj+Xrb6RpJDiiFAoUd2rkYVr5Ydd/Jm9OdzfCuSqJTvfK4sHP+oUqik54M4b3VjsZqGl0Z2dIKCwuJjY0lMTGx1seXLFnCjBkzmD17Nrt27SI2NpbRo0eTkZFhclx2djZ33303H3/8sS3CFg5g5tJ9NfZ9+8AQNBqNCtGIBp3eYnq/JLf244RdauHrCUB2oSQ6DqPLNcrt/u/VjcPKVF/JcOzYsYwdW/e00Xnz5jFt2jSmTp0KwIIFC/jll19YuHAhM2fOBKCkpIQbbriBmTNn8re//a3e1yspKaGkpMR4Py8vr56jhSNKPpfHG78fqvWxnhGBNo5GmKW8FPZ9Y7pvpAxAdiQt/ZREJ0dadBxHZH/l9lK2unFYmeotOvUpLS1l586dxMfHG/dptVri4+PZvHkzUDVI+qqrruKuu+5q8Jxz584lKCjI+CPdXM7nH//dwdrDplMmN828inVPjsTXU/XcXtRm6X2m9zuOhKH/VCMS0UQh/spElT1ncnnkq92U6WSsjt3z8lduSwrUjcPK7DrRycrKQqfTERYWZrI/LCyMtLQ0ADZt2sSSJUtYtmwZcXFxxMXFsW9fzS6LSs888wy5ubnGn9TUVKv+G4TtpWZfMrk/bVgHIoN9aN9KBqjbpYIMSF5uui+yH7h71X68sEvtW1XVpfpxzznWHXb++iwOzzNAuS0tVDcOK3P4r7dDhw5tVHVmLy8vvLy8SExMJDExEZ1OZ8XohK0dPG/aFfnGjX24ZYC02tmtI78rg5Avp7Hr72CiFsEVY3QqyVA4B+BZ8eWvVFp0VBMSEoKbmxvp6ekm+9PT0wkPb95UuISEBJKTk9m+fXuzziPsy40f/mVy/4r2weoEIszzv1vg8K/KdsdRVfsl0XFI1ZMbnd6gXiDCPMZEx7lbdOz63cTT05N+/fqxatUq4z69Xs+qVasYMmSIipEJe1Su01NUatpCFy3dVfar/LJBqy2ioUVFMcfu420ejmi+92+vqnt0+d+isEOVY3R0JU5dzkH1rquCggKOHTtmvJ+SkkJSUhItW7akXbt2zJgxg8mTJ9O/f38GDhzI/PnzKSwsNM7CairpunI+h9PzTe5/cnd/3N3sOpd3bQWmLbW07gbX/gvy0yCkizoxiWb5e582/HulP8cyCigslTWv7J5HtS+CpQVKtWQnpHqis2PHDkaNqmqynjFjBqBUYl68eDG33normZmZzJo1i7S0NOLi4lixYkWNAcqNlZCQQEJCAnl5eQQFBTXrXMI+7DqdY3JfKiDbsfQDkHXUdF9oD/AKUH6Ew+oVEcixjAKKSuRLpN1z9wQ3T9CVKt1XkuhYx8iRIzEY6u/LnT59OtOnT7dRRMIRXSgo4YVl+433u4T60znUX8WIRJ0KMuDDWupdhfexfSzC4ny9lI8VadFxEJ5+cKnUqaeYu2y7fmJiIjExMQwYMEDtUIQFVG/N+d99g/j9n8PxkG4r+5SdYno/vA88tAV8W6oTj7AoP083QMboOIzASOX21EZ147Ail/0kkFlXziUtV6md4+mmZUinVmi1MrfVbpWb1jmi6xil20o4hcqinIUl0qLjECLilNtfHoeyS/Ue6qhcNtERzuOdlUd5YfkBAG4fGCVrWdmznNPw9R2m+yL6qhOLsAo/L6VFRxIdB9G12hJM+WnqxWFFqo/REaI5zuVc4t8rjxjv+8gSD/Yr/YDp2Bz/MBj+JHSre6074XjCAr0BOJvjnK0DTqfHeEADGKDYORfSddkWHRmj4/h+2XueWz/ebLJvTK/mFZIUVvTDP0zvtxsMA6dJCV0n06m1Mglg+8mLlMt6V46hdXflVhId5yJjdBzbsYwCEv63y7iu1XWxERx8eQxxUcHqBiZMpe2D83uV7cwjpo858SwPV9axdVVtlg/WHlcxEmE274oSK6e3qBuHlUg7v3BIz/1gunDr0C4h+FTM9hB2orwEFgxVtv3DlOqr1RVk2D4mYXW+nu7EtAkk+Xwe209mqx2OMEflBIG1c2DEU07XyuqyLTrCcWXkF7PtsjfQK9oFqxOMqFv1FpvLqyCD0y8k6Mpeu7E3APvP5jZYJ03YgYxDVdvlJXUf56BctkVHloBwTAaDgaPpBRgMShP549d0I7uolM6hUk3X7pQ1sFDg+Hm2iUPYXOU4nYtFZRSV6vDzctmPGsdQfRHd0kLw8FYvFitw2d8+WQLC8WxLyeYfX+ygT9tgANoEefP3Pm3UDUrUrbYVkdvEwuAE6DUJ3DxsH5OwCV9PNzzdtZSW68kuLJVEx97d+AksuVPZLi0Av1bqxmNh0nUlHMa0z3dwsaiMdUcyAQgNcK5vHU6ntKjmvv73QOytkuQ4OY1GQwtf5f84p8h5V8V2Gj0mVK1zVVbL362Dk0RHOIzcS6ZvmKEBXipFIsySc6rmvspy88LptfD1BCC7qFTlSIRZKhfTra0l1sFJoiMcVmtJdOzbd1Nr7pP1rFxGZaKTI4mOY/CoKAvghJMEXDbRkYKBjiX5XF6NfaGB0nVlt/JrmWUF4COJjqto6VfRolMoiY5D8KxMdKRFx2lIwUD7VqbTk3upjMKScorLdOxOvVjjmDBp0VFfcS7oLhuDsf97eLtr7cdXjgMQTi/QRxmAnF8sa145BGOi43xjdGQovLBL//fJFnadzkGnNxAe6M29QzvUOEZadFR2KQfmxUBwO0ioVlH1u3vqfo63zHB0FcZVzEsl0XEInkpJAOm6EsLKSsv15BaVsf3kRXR6pdBYWl5xrQMaZTCyys4nKbVyMg9C2n74742Qsr7+5zhZxVVRN7+KSuU5hWVSNNARVA5GLrqgbhxWIC06wm4UlpQz6q21XKwlqfkx6ZzJ/aiWPlKbQ22V3wABPrtWSXqOrTQ9ZuhjcPEkHPgB/GXBVVfiW/H3uWRHKgYMvHFTrMoRiXqF9lBu0/aqG4cVyCeFsBtrD2eSkV97+fGzOcpaLHcMasfonuH0jAi0ZWiiNuXFVdt1VUGOf1G5veoF8A+1ekjCfvhVW3vumx1nJNGxd5FXKLfn96gbhxVIoiPshjl9+e1b+TK8a2sbRCMaVD3RaUirTtaLQ9ilyjE6lYpKy2vsE3bEr+J9tSRf3TiswGXH6Mj0cvtTVNJwoiNvlHakvIFpw23lb8uV+Xm5mdw/ne18s3mciptSDqDBv2sH5LKJjkwvtz+FpQ0vsHr5m6dQUX0tOjctgsk/2S4WYXcu/1Jy6oIkOnbNvWJyh875Vi932URH2J8CadFxLLp6vvlF9gMPH9vFIuyO9rIZdqcl0bFvbpWJTik42Sw5SXSE3ci71PDifwEy08p+XN6iEz2satsn2KahCPvj42na+noq2/kq7joVd8+q7fq+xDggSXSE3Sgyo+uqbzuprGs3yi9r4h44rWrbM8C2sQi7c0W7YB69ugvDuoQAcPbiJZUjEvVyq1aX7PK/bQcnX4+F3bh0WaLTrqUv43q3wcfDjXK9nlHdQ2t8SxQquvzNsMNwmPAOeAeDVr5DuTqNRsNj13Tll73n2XA0y6wxeEJFbs7boiOJjrAbl8pM3wiviQlj5tjuKkUjGnR5ouMdDP2mqBGJsGM+nkrSW1wmiY5d02pB6wH6Mqdr0ZGvXUJVxzMLuFixuvHliU67lr5qhCRqc/Gk8gPKQMWzO6HkshXlZXkHUQtvd6UVVhIdB+CkM6+kRUeoZvvJbG5esJmB0S356K5+ZBeaNpfePrCdSpEJE+Wl8E5FVdvn0pRlHpbcaXrMXctsHpZwDN4V3c2Xf5ERdshJa+lIoiNU8+aKwwBsO5lN33/9afJY78ggPN2lwdEuVK+UmnsGdn1h+vjwJ6HTKNvGJBxGVYuOXuVIRIOctEXHZT9JpDKy+rIK6v5j0koviP0orzZbJj8NPP1MH3eTVeRF3SonEBTLYGT756QtOi6b6EhlZPWV6ev+hqeR8R72o6xaopO8DA4sNX08MMKm4QjH4u1RMRi5XBIduyctOkJYVrmu7uqbY3qF2zASUa+yahVtt39a8/Hu42wXi3A4lV1XZToD5TrpvrJrla2zTjbrSsboCNWUltf+pjd7Qgx3Dm5v42hEncrqKPR270qlmqqPFHEUdate+6q4XI+/m3y/tluV1ZGdrI6O/MYJVRgMBvKKa1/y4cZ+bfGQN0P7UVbLGkXj3oKoAdAm1vbxCIfiVW1SweVFQYWdqWzR+eo2KHWeJTvk00TYXFZBCWsPZ1JWR9eVnyzcaV8ub9HRukP/e9WJRTgcjUZTNU5HppjbN++gqu1jK9WLw8LkE0XYXP9X6v8DcpMpV/bl8kTHP1yWeBCN4uPhRnGZXmrp2Lu2/eDIb8q2h1/9xzoQebcSNlVUWl7v4yv+Oazex4WNbPkQdn8Jf7wA31drvYnsByNnqheXcEhBPh4AvLPqqMqRiHoNerBq2+A8A8elRUfY1JkGVjDuHh5oo0hEnbKOwYpakpmek+DmRbaPRzg8Py/lo+aXved5+boSWvlL7SW75OUPkf3h7A7Q1/+l1JFIi46wuPVHMrnu/Y0cOJdr3Fdarueuz7Zy7b/XqxiZqFNxHiwaB1s/gsLM2o9x87BtTMJpVB+bUyQDku2btqL9QxIdIeqW8L9d7D2Tyy0LNhv3HTiXy4ajWTWOvbV/lC1DE3XZ/QWc2gS/PVX3G5zMsBJNVH35BxmnY+eMiU7ts2IdkVMkOhMnTqRFixbcdNNNaocigJKK+jiF1b655RTV/kfz+k19aOnnaZO4RD00VbVOarzB9boRErbB4IdsG5NwGuXVqqBLi46dc6tMdJzn/8kpEp1HH32Uzz//XO0wRIV2LX1r7Mu5VHcBqtAA6a9XnV9I1XZ+mulj/uHQuhvIshyiid68qao1sKEJCUJl0nVln0aOHElAQIDaYYgKekNVfZxJH2wit6iMi4XO0wzqdI78YTqzatmDpo97yd+WaJ7hXVvTLUz5Pfq/T7ZyLKNA5YhEnSTRsbz169czYcIEIiIi0Gg0LFu2rMYxiYmJREdH4+3tzaBBg9i2bZvtAxVm0+mrEp1dp3P4aP1xci5JomO3/ndz/Y97y0w40XzBvlWD2Z/+fq+KkYh6SaJjeYWFhcTGxpKYmFjr40uWLGHGjBnMnj2bXbt2ERsby+jRo8nIyGjS65WUlJCXl2fyIyyreqIDcLGojJwi51o7xWEYDLDh7aoqp9kpsPIlKKj4+9m/tO7nVpIWHWEBvtXWvMqVLz72S1vx/yRjdCxn7NixvPLKK0ycOLHWx+fNm8e0adOYOnUqMTExLFiwAF9fXxYuXNik15s7dy5BQUHGn6gomfVjafrLEp1ynb7WN7Z+7WUxSKs7+iesehn+e6Nyf+Fo2DgPlt6v3P9uasPnkEU7hQX4VlvaxV2qn9svbUXLm7To2EZpaSk7d+4kPj7euE+r1RIfH8/mzZvreWbdnnnmGXJzc40/qamplgpXVNAZTBOd09lFZBeatujMvzWOzyb3t2VYrinvjOn9gnTl9vRmpbXHHK27WzYm4ZKqr2Iuy7zYscquK53ztLrZdWXkrKwsdDodYWFhJvvDwsI4dOiQ8X58fDx79uyhsLCQtm3b8u233zJkyJBaz+nl5YWXl8zysabLu662pmTXOOaGvpG2Cse1aev4E9e4md803bKj5eIRLqt615W7m11/x3ZtTjhGx64THXOtXNn4VVYTExNJTExEp3Oefkh7Ua43s6VAWF9diY7Wzfw3Mq1bw8cI0YDqZSQ8pEXHfskYHdsKCQnBzc2N9PR0k/3p6emEh4c369wJCQkkJyezffv2Zp1HKI6k53Oxonvq8hYdoaI6W3S0TvWNTdi/qGr1taTryo45YYuOXSc6np6e9OvXj1WrVhn36fV6Vq1aVWfXlLC945kFXPvv9Qx/cw0A5brGJTptW9QsMCgspHprjK7cdL8TvZEJ+9e+lZ9xWyvFJ+2XEyY6qnddFRQUcOzYMeP9lJQUkpKSaNmyJe3atWPGjBlMnjyZ/v37M3DgQObPn09hYSFTp5oxW6Qe0nVlOTtPXQQgv7gcnd5AqU5f57H92rdg6pXRJvvmTOyFp7uGOwe3t2aYrql6i46upGq7MWN0hLCA2LZBtG3hw5mLl2S9K3vm5nyzrlRPdHbs2MGoUaOM92fMmAHA5MmTWbx4MbfeeiuZmZnMmjWLtLQ04uLiWLFiRY0Byo2VkJBAQkICeXl5BAUFNetcrkBfkcB4eygtBAaDgZJy5X4L36q1qs5cLKrRdTW2Vzi/7VeWFfhscn+CfU3XtgoN9OaDO/pZ+V/gorTVVhwvK66238wWHY1dN/oKB6LRaHhtUh/u/GwrlyrWuyqr+FLkIYOT7YdxjI4kOhYzcuRIDA1Mc50+fTrTp0+3UUSiNnd8upXtJ7PZ8Xw8wb6ePPTlLn7bn8ZfM6+ipLzq29mhtPwaz/XxqOo+8fNS/VfOdS24smpbY2aiU9cYHyGaoHKKeVFZOXq9geFvrEFvMPDXzKtl3I69cMKuK5dNoxMTE4mJiWHAgAFqh+IQNp+4QLnewKqDSkXdyhaab3akGr+dAaRkFdZ47qQr2tLKz5PrYiPkm5utGap1EeSfr9rWmjkYWSMzroTltPJTWnMz80vIuVTG+dxi0vNKuFBY0sAzhc04YaLjsl/XpOuqacr1puNvdHqDSX/7qQtFNZ4ztEsIO1+4xuqxiVrU9WZVvUXHwxfKav6/ATK1XFhUm2BvNBooLtNzsdqyMObWrhQ24ISJjny9Fo1SetmMKp3eYNKic+qCaYuOh5s0RzfbmR1wekvTnlvXgOOLKXDwJ2XbzbP2Y0ASHWFRXu5uxno6p6t9KdJLpmM/JNERrk532YwqvQGTFp0zFy+ZPO6ulV+xZtGVwadXK2tUFTdhAVpD3TPgWDlbuS3OqfsYGaMjLCwi2AeA1ItViU5jS1IIKzImOs4zM85l38VkennTXF71WG8wbdG5fJVyadFpptJqLWQl+eAd2Ljnm/tmNW2N0sIT1hNyz1QlQTJGR1iYf8WEhIKSqhaDsnpKUggbc8IWHZdNdGSMTtOU6S5fmdxgMm4nr9j0j0MGHzdT9TebprSuGMxMdCKvUH4qVSY60qIjLMyz4j2h+heky99XhIqcMNGRTyFRp3KdHp3egL5aK47ussHIJeU6kzesy7lLi07z6EobPqa8pPZtaH7zs4zRERbm6a587BSWVE90pEXHblT+zTvR6uWS6Iha6fQG4uetY+w76026q8p0BpOCgF9uPU12Yd0fxtKi00zVE5faxttkHITX2sHKl2DrR/BqOBxfU+05zUx0PGR5DmFZlYnOwk0pxn2S6NgRJ6yM7LKfQlJHp37nci5x8kIRR9ILyCuuyuzL9foab0qrDmXUeR5JdJqpeotObUnL789BeTFsnAe/PaUkQ0unVT3e1BadiR9BUBTc+GnTni9EHTxreU+Qris7UlkN3Ylmwrnsp5CsXl6/6l1OhSaDBg01BiTXex6pdto81Vt0aktaSmsWaDR5g6pv1lV9Ym+Dx/ZDmz5Ne74Qdahs0amuXFp0hBW5bKIj6le9eyq/2gDjcp2BsnLz35TcpUWneRpq0SktqOVJ1RIdJ5oiKpxDbYlOfQsBC9Fc8ikkalW9Kbn6NNDauq7qI9PLm8lkjE4tLWklNdcWM23RkURH2JfaEh3puhLWJImOC3p/9VHm/naw3mOqNyUXVGvRKSnTU9aIrisZo9NMusu6rvLT4IO/wevRcC6p/hadta/Dn7NsEKQQ5vOq5T1Buq6ENbnsp5CrDkY2GAy89ccRPlp3gtTsOtY3wvQbVn5J1WDkojKdWV1XPdoohe3uGNSuGdEKyi/ruvrpn5BxAC5dhE9G1ZxODkqLTnkJrJ3T9DE6QliJdF0JW3PZamBSMBCK6ql/U70IYPUxOrmXyhrsuvrgjiu4qnsoxzIK6BnRyEq+wpTusunl2SdM72tq+65iqD0BEsIOSNeVsDWXTXREzZXIq6uezFye6DT07atrmD/eHm70inTNBNKiqrfo6HU1x9xoahkDZcC8QoNCqKD26eXSoiOsx2W7rgTUk+eYdl1VS3TyL5U1uACfl7tU07UYkxYdXc0iXtKiIxyMZy3vDzJGR1iTJDourL4WnXKTWVdVY3TM6bqqrWlaNJFJHR19LdlpbS06BtMESQg7UttMzFLpuhJWJJ9ILkxfT+XL5nRd1dY0LZrIpI6OvmaLTq1l2g2mXV5C2JHax+hIi46wHvlEchHf7zzD39/dwNmcS8Z91VttHv16Nw9/tRuDwcD8lUeYuriqYvTypHNVz9EbyLtU/2Jv0qJjAbv/CwuGwcWTVfsMtYzRqa2LqrQAPhhk1fCEaCovqYwsbMxlByMnJiaSmJiITucaBdUe/3YPAP/6Odm4T1fRopNfXGZMZp7/ew/mrzxa77myCyXRsbrlCcpt2t6qffpaxuhYootq0IPNP4cQZuoaFlBjX3GZJDrCelz2E8lV17qq/oZSucxD9fp/xWUNJ36Vi3wO7RxS6+OyvpWVGPTNX9Jh6m/wzBmIqtbiM/rV5p1TiEbo2Nq/xr76Sl0I0Vwum+i4quoDASsTnerrWpWYUQywsuvK17P22VWa2qY8i+Yz6Jqf6AS1Ba8AcPOs2lfrzC0hrMftsi9DRaW1jTUTwjLkHc7FVH+DqUxwqvePXzLjm1VuA4mOsJLa6ug0lmfNb9O11uIRwooun7BQKC06wook0XEx1VcTNyY61Vp0Cs34ZpVXMQvL18tlh3ipo7ZZV43l6WeZWIRohsvH8RWVSIuOsB5JdFyMe60tOlWJTlGJGWN0Klp0/KRFx7YsMUansstKWnGEii5PdMz5giVEU0mi42JMuq4MlS06VV1XRWYMRq7qupIWHZuyRNeVMcGRREeo5/KuKxmMLKxJEh0X46Gt2XVVfTByoRlNyMYWHS9p0bGp5iY5QtiJGi060nUlrEgSHQeXkV9M7qUyist0HE7LJyOvuN7jT14oNG5fLCwlNbuIPWdyjfv2nc2t7WkmTmQp55AWnSbKPQOlRVCcC/npyr6yYrh4qv7npW6zXAzSdSVUdHmLTk5RGanZRSpFI5ydfFI5sLziMga+ugqA8EBv0iqSnJ3Px9PK36vW52xNyTZuv/hTMi/+lGzy+P+2njb79f1lMHLjZR2F9/tDQATkV1ScfioFFo2DzINw32po26/25/71rgUDkURHqKdtCx8Op+cb718oLGXYG2v45h9DGNihpYqRCWfksi06iYmJxMTEMGDAALVDabLjGQXG7bRqLTn7z+XZ5PWHdQnh+wf/xpWdW7Hk/sGM6Nqad26Ls8lrO6xDvyi3+VXLapBxUElyAA4std5rT/y4altadISKXpnYi/BA7xr7v995RoVohLNz2a/kCQkJJCQkkJeXR1BQkNrhOJzwQG9a+XvRyt+LL+8bDMCgjq1UjsoB6BpYbLN6Ib+muuFDWHbZsg7XfwCxtzb/3EJYQJsgH7Y8ezWnLxQx/M01aocjnJzLtugI88gXfwurbRHO6txr73Jsthrnlf9YIYRrkERH1MtPBhxbVkOLcFqiRcec80oGK4RwEZLoiHr5SFFAyypvoOtKWnSEEMKiJNER9ZLqxxZmLy06QgjhIiTRcQEGg6Hhg+rgI11XlmWNFh1zVh+//LzSdSWEcBGS6LiA6ot2NtblFUxFM9XWolO94rFbExIdbzNmDdZo0ZFERwjhGuRTzAXompHoeGjlA9Giapt1Vb2Vx70JXUxNSnSEEMI1SKLjAsp0+oYPqoNWEh3Lqq2OTvVWHmu16EjXlRDCRUmi4wKa06IjeY6F1dqiU22ftgljoqTrSggh6uQUic7PP/9Mt27d6NKlC59++qna4did5ozRERZWa4tOtX1NaWmRFh0hhKiTw0+pKS8vZ8aMGaxZs4agoCD69evHxIkTadVKliOoVK6TRMduNNSi0xRmtehYqT6PEELYOYdv0dm2bRs9e/YkMjISf39/xo4dyx9//KF2WHalXN/0MToa6eKwrIZadJrCO7jhY2oMcpb/VyGEa1A90Vm/fj0TJkwgIiICjUbDsmXLahyTmJhIdHQ03t7eDBo0iG3bthkfO3fuHJGRkcb7kZGRnD171hahOwxp0bEj9tKiI11XQggXoXqiU1hYSGxsLImJibU+vmTJEmbMmMHs2bPZtWsXsbGxjB49moyMjCa9XklJCXl5eSY/1pKaXcST3+7haHq+cd/nm0/yzsqjbD5+gfv+s4NVB9MbPM/ypLO8/FMy+kaMtVl7OIOpi7bR9fnfmLp4e5Pidzp6Pax4BvZ9p9xPT4ZlCXDxlPnn2PUFrHwRDAY4thJ+ehRKi2DDPNj8ARTnwo+PQMqGitd7FrZ+DF9MhBeD4MLRmuf847mq7S9vgnVvNO7fZVai49G4cwqhgiU7Uhv1PieEOVQfozN27FjGjh1b5+Pz5s1j2rRpTJ06FYAFCxbwyy+/sHDhQmbOnElERIRJC87Zs2cZOHBgneebO3cuL730kuX+AfW49z/bOZJewK/7znPg5TEAzFp+wOSY9UczOfJK3f9+gEe/TgJgYIeWjOkVbtZrJ645xvaTFwFIySpsZORV7hzcns0nLuCu1TC0SwhrD2c2+VyqO/QzbPlA2e59E3w8Quk2Or8HHtxo3jl+nK7cdhsH/71R2dZoYcdCZTvrMOz6j/Jz+9ewpfYEvl5rXm3c8cHtGz7m8hacnhOV6xEYWfvxQqhkxYE0xvVuo3YYwomonujUp7S0lJ07d/LMM88Y92m1WuLj49m8eTMAAwcOZP/+/Zw9e5agoCB+++03XnjhhTrP+cwzzzBjxgzj/by8PKKioqwS/5H0AgAKS3V1HlNabv74mexC88dylDWyu8rX0w1fT3eyCpRulPBAbxZOGUBMRCDhQUNo19KP4jIdw95Y06jz2pXCy5K0yrEx6fsbf65LOVXbWdVaaS4cr/Z6WY0/b6XWPSDzYN2P//1tiOwHRReg41V1H9cmFm5eXHN/rxshoA2E9mh6jEJYwYWCZnblCnEZu050srKy0Ol0hIWFmewPCwvj0KFDALi7u/P2228zatQo9Ho9Tz31VL0zrry8vPDykhkol/tbpxDctPD7AaUrrV1LX2IiAgHo174lABn5xarFJy4T3gci+jZ8XFhvaNmx5n6NBqKvtHxcQghhZ+w60THXddddx3XXXdeo5yQmJpKYmIhOV3drixBCCCEcm+qDkesTEhKCm5sb6emmA3bT09MJDzdvrEpdEhISSE5OZvt2GagrhBBCOCu7TnQ8PT3p168fq1atMu7T6/WsWrWKIUOGqBiZEEIIIRyB6l1XBQUFHDt2zHg/JSWFpKQkWrZsSbt27ZgxYwaTJ0+mf//+DBw4kPnz51NYWGichdVU0nUlhBBCOD/VE50dO3YwatQo4/3KGVGTJ09m8eLF3HrrrWRmZjJr1izS0tKIi4tjxYoVNQYoN1ZCQgIJCQnk5eURFGRGHRIhhBBCOBzVE52RI0diMNQ/FXr69OlMnz7dRhEJIYQQwlnY9Rgda0pMTCQmJoYBAwaoHYoQQgghrMRlEx2ZdSWEEEI4P5dNdIQQQgjh/CTREUIIIYTTctlER8boCCGEEM7PZRMdGaMjhBBCOD+XTXSEEEII4fxUr6OjtsoaPnl5eRY/t76kyLhdef7q+y5/rKHzFBXmmxxbkJ9X6/kKC/Ipu1RY62N1Kb1UgJu26rXKLnnViCs/v9j4eHmxzirXzKoKLkFJRc2mvLyq7cr75qh8TkFh1falstq3q79eY10qr/+5+QWmMVc/tqCo6n5Rqfn/NiFsLD+/qMb71KXL3ueEjVW+bxWV2f17R+XvSUO1+DSGho5wUpVLQJSWlnL8+HG1wxFCCCFEE6SmptK2bds6H3fZRKeSXq/n3LlzBAQEMHDgwFrH7AwYMKDG/ob25eXlERUVRWpqKoGBgdb7B9QTj7Web86x9R1T12NynRt/rKWuc237XeU6m3N8Ux8353f68vvOfK3lvUOusyUZDAby8/OJiIhAq617JI7Ld11ptVpjJujm5lbrf0Rt+83dFxgYaJM/orpit8bzzTm2vmPkOtvfda5tv6tcZ3OOb+rj5vz+1vVcZ7zW8t4h19nSzFmrUgYjV5OQkGD2fnP32UpzX7sxzzfn2PqOketsuWMtdZ1r2+8q19mc45v6uDm/v2peZ0u8vr3+Tst7R9OPceTrXBuX77qylspV0XNzc23ybcFVyXW2DbnOtiPX2jbkOtuGPVxnadGxEi8vL2bPno2Xl5faoTg1uc62IdfZduRa24ZcZ9uwh+ssLTpCCCGEcFrSoiOEEEIIpyWJjhBCCCGcliQ6QgghhHBakugIIYQQwmlJoiOEEEIIpyWJjgp+/vlnunXrRpcuXfj000/VDsdpTZw4kRYtWnDTTTepHYpTS01NZeTIkcTExNCnTx++/fZbtUNySjk5OfTv35+4uDh69erFJ598onZITq2oqIj27dvzxBNPqB2K04qOjqZPnz7ExcUxatQoq72OTC+3sfLycmJiYlizZg1BQUH069ePv/76i1atWqkdmtNZu3Yt+fn5/Oc//+G7775TOxyndf78edLT04mLiyMtLY1+/fpx5MgR/Pz81A7Nqeh0OkpKSvD19aWwsJBevXqxY8cOee+wkueee45jx44RFRXFW2+9pXY4Tik6Opr9+/fj7+9v1deRFh0b27ZtGz179iQyMhJ/f3/Gjh3LH3/8oXZYTmnkyJEEBASoHYbTa9OmDXFxcQCEh4cTEhJCdna2ukE5ITc3N3x9fQEoKSnBYDAg31Ot4+jRoxw6dIixY8eqHYqwAEl0Gmn9+vVMmDCBiIgINBoNy5Ytq3FMYmIi0dHReHt7M2jQILZt22Z87Ny5c0RGRhrvR0ZGcvbsWVuE7lCae52F+Sx5rXfu3IlOpyMqKsrKUTseS1znnJwcYmNjadu2LU8++SQhISE2it5xWOI6P/HEE8ydO9dGETsmS1xnjUbDiBEjGDBgAF9++aXVYpVEp5EKCwuJjY0lMTGx1seXLFnCjBkzmD17Nrt27SI2NpbRo0eTkZFh40gdm1xn27HUtc7Ozubuu+/m448/tkXYDscS1zk4OJg9e/aQkpLC//73P9LT020VvsNo7nVevnw5Xbt2pWvXrrYM2+FY4vd548aN7Ny5kx9//JE5c+awd+9e6wRrEE0GGH744QeTfQMHDjQkJCQY7+t0OkNERIRh7ty5BoPBYNi0aZPhhhtuMD7+6KOPGr788kubxOuomnKdK61Zs8Zw44032iJMp9DUa11cXGwYNmyY4fPPP7dVqA6tOb/TlR588EHDt99+a80wHV5TrvPMmTMNbdu2NbRv397QqlUrQ2BgoOGll16yZdgOxxK/z0888YRh0aJFVolPWnQsqLS0lJ07dxIfH2/cp9VqiY+PZ/PmzQAMHDiQ/fv3c/bsWQoKCvjtt98YPXq0WiE7JHOus7AMc661wWBgypQpXHXVVdx1111qherQzLnO6enp5OfnA5Cbm8v69evp1q2bKvE6KnOu89y5c0lNTeXkyZO89dZbTJs2jVmzZqkVskMy5zoXFhYaf58LCgpYvXo1PXv2tEo87lY5q4vKyspCp9MRFhZmsj8sLIxDhw4B4O7uzttvv82oUaPQ6/U89dRTMmuikcy5zgDx8fHs2bOHwsJC2rZty7fffsuQIUNsHa5DM+dab9q0iSVLltCnTx9jP/0XX3xB7969bR2uwzLnOp86dYr777/fOAj54YcflmvcSOa+d4jmMec6p6enM3HiRECZUTht2jQGDBhglXgk0VHBddddx3XXXad2GE5v5cqVaofgEoYOHYper1c7DKc3cOBAkpKS1A7DpUyZMkXtEJxWx44d2bNnj01eS7quLCgkJAQ3N7caAwTT09MJDw9XKSrnI9fZduRa24ZcZ9uQ62wb9nadJdGxIE9PT/r168eqVauM+/R6PatWrZIuEwuS62w7cq1tQ66zbch1tg17u87SddVIBQUFHDt2zHg/JSWFpKQkWrZsSbt27ZgxYwaTJ0+mf//+DBw4kPnz51NYWMjUqVNVjNrxyHW2HbnWtiHX2TbkOtuGQ11nq8zlcmJr1qwxADV+Jk+ebDzmvffeM7Rr187g6elpGDhwoGHLli3qBeyg5Drbjlxr25DrbBtynW3Dka6zrHUlhBBCCKclY3SEEEII4bQk0RFCCCGE05JERwghhBBOSxIdIYQQQjgtSXSEEEII4bQk0RFCCCGE05JERwghhBBOSxIdIYQQQjgtSXSEEEII4bQk0RFCOI0JEyYwZsyYWh/bsGEDGo2GvXv32jgqIYSaJNERQjiNe++9lz///JMzZ87UeGzRokX079+fPn36qBCZEEItkugIIZzG+PHjad26NYsXLzbZX1BQwLfffsu9996rTmBCCNVIoiOEcBru7u7cfffdLF68mOrrFX/77bfodDpuv/12FaMTQqhBEh0hhFO55557OH78OOvWrTPuW7RoETfeeCNBQUEqRiaEUIPGUP1rjxBCOIErr7ySTp068fnnn3Ps2DG6dOnCmjVrGDlypNqhCSFsTFp0hBBO59577+X7778nPz+fRYsW0alTJ0aMGKF2WEIIFUiiI4RwOrfccgtarZb//e9/fP7559xzzz1oNBq1wxJCqEC6roQQTum+++5j6dKl5OXlcfr0aSIiItQOSQihAmnREUI4pXvvvZeLFy8yevRoSXKEcGHSoiOEEEIIpyUtOkIIIYRwWpLoCCGEEMJpSaIjhBBCCKcliY4QQgghnJYkOkIIIYRwWpLoCCGEEMJpSaIjhBBCCKcliY4QQgghnJYkOkIIIYRwWpLoCCGEEMJpSaIjhBBCCKcliY4QQgghnNb/A7+S368jDURuAAAAAElFTkSuQmCC\n"},"metadata":{}}]}],"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3 (ipykernel)","language":"python"},"language_info":{"name":"python","version":"3.7.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":5} \ No newline at end of file diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst index 1bac418..c66fae2 100644 --- a/docs/source/tutorials.rst +++ b/docs/source/tutorials.rst @@ -8,7 +8,8 @@ We provide multiple tutorials showing how to read and manipulate Quijote data. I .. nbgallery:: - + + Examples/Reading_snapshots Examples/Density_fields Examples/Density_fields2