diff --git a/examples/accretion/Makefile b/examples/accretion/Makefile new file mode 100644 index 000000000..cd123bddb --- /dev/null +++ b/examples/accretion/Makefile @@ -0,0 +1,35 @@ +export OPENGL=0# Set this to 1 to enable OpenGL +export SERVER=1# Set this to 1 to enable the visualization web server +include ../../src/Makefile.defs + +# CCPROBLEM is defined in Makefile.defs to allow for +# a compact cross platform Makefile +.PHONY: all librebound +all: problem.c librebound + @echo "Compiling $< ..." + $(CCPROBLEM) + @echo "" + @echo "Compilation successful. To run REBOUND, execute the file '$(EXEREBOUND)'." + @echo "" + +librebound: + @echo "Compiling shared library $(LIBREBOUND) ..." + $(MAKE) -C ../../src/ + @-$(RM) $(LIBREBOUND) + @$(LINKORCOPYLIBREBOUND) + @echo "" + +clean: + @echo "Cleaning up shared library $(LIBREBOUND) ..." + $(MAKE) -C ../../src/ clean + @echo "Cleaning up local directory ..." + @-$(RM) $(LIBREBOUND) + @-$(RM) $(EXEREBOUND) + +rebound_webgl.html: problem.c + @echo "Compiling problem.c with emscripten (WebGL enabled)..." + emcc -O3 -I../../src/ ../../src/*.c problem.c -DSERVERHIDEWARNING -DOPENGL=1 -sSTACK_SIZE=655360 -s USE_GLFW=3 -s FULL_ES3=1 -sASYNCIFY -sSINGLE_FILE -sEXPORTED_RUNTIME_METHODS="callMain" --shell-file ../../web_client/shell_rebound_webgl.html -o rebound_webgl.html + +rebound_console.html: problem.c + @echo "Compiling problem.c with emscripten (WebGL disabled)..." + emcc -O3 -I../../src/ ../../src/*.c problem.c -DSERVERHIDEWARNING -sSTACK_SIZE=655360 -sASYNCIFY -sSINGLE_FILE -sEXPORTED_RUNTIME_METHODS="callMain" --shell-file ../../web_client/shell_rebound_console.html -o rebound_console.html diff --git a/examples/accretion/problem.c b/examples/accretion/problem.c new file mode 100644 index 000000000..2d88e4128 --- /dev/null +++ b/examples/accretion/problem.c @@ -0,0 +1,189 @@ +#include +#include +//#include +#include +#include +#include "rebound.h" + +void heartbeat(struct reb_simulation* r); + +double e_init; // initial energy +int Nparticles = 1000; +double tmax = 6*M_PI; + +char title[100] = "trace_nc_ns"; +char title_remove[100] = "rm -rf trace_nc_ns"; + +char remove_snapshots[100] = "rm -rf *trace_snapshot_*"; +char snapshot_1[100] = "trace_snapshot_1"; +char snapshot_2[100] = "trace_snapshot_2"; +char snapshot_3[100] = "trace_snapshot_3"; +char snapshot_4[100] = "trace_final_orbit_snapshot_4"; + +double snap1_time = 0.0; +double snap2_time = 10.0 * 2 * M_PI; +double snap3_time = 500. * 2 * M_PI; +int snap1=1; +int snap2=1; +int snap3=1; + +clock_t begin; +// accretion of the moon +int main(int argc, char* argv[]){ + struct reb_simulation* r = reb_simulation_create(); + r->integrator = REB_INTEGRATOR_TRACE; + //r->ri_ias15.adaptive_mode=2; + + r->dt = 0.2;//0.003 / 10.56789; + //r->integrator = REB_INTEGRATOR_IAS15; + //r->softening = 3e-8; + //r->collision = REB_COLLISION_DIRECT; + //r->collision_resolve = reb_collision_resolve_merge; + + // Initialize masses + struct reb_particle earth = {0}; + earth.m = 1; + double earth_r = 1./2.9;//0.00592; // in units of Roche radius + earth.r = earth_r; + reb_simulation_add(r, earth); + + r->rand_seed = 1; + r->heartbeat = heartbeat; + + double lunar_mass = 0.0123; + + double mtot = 0.; + double rad_tot = 0; + // Test particles + for (unsigned int i = 0; i < Nparticles; i++){ + double m = reb_random_powerlaw(r, 3.2e-7, 3.2e-4, -1.); + mtot += m; + double rad = pow(m, 1./3.) * 1.185 * earth_r; + rad_tot += rad; + double a = reb_random_powerlaw(r, earth_r, 1.5, -1.); + double e = reb_random_uniform(r, 0., 0.95); + double inc = reb_random_uniform(r, 0, 50. * M_PI / 180.); + double Omega = reb_random_uniform(r, 0, 2 * M_PI); + double omega = reb_random_uniform(r, 0, 2 * M_PI); + double f = reb_random_uniform(r, 0, 2 * M_PI); + reb_simulation_add_fmt(r, "primary m r a e inc Omega omega f", earth, m, rad, a, e, inc, Omega, omega, f); + } + //printf("%f\n", mtot/lunar_mass); + //printf("%f\n", rad_tot/Nparticles); + system(title_remove); + //system(remove_snapshots); + //exit(1); + + reb_simulation_move_to_com(r); // This makes sure the planetary systems stays within the computational domain and doesn't drift. + +/* + FILE* f = fopen(snapshot_1,"a"); + for (unsigned int i = 1; i < r->N; i++){ + struct reb_particle* p = &r->particles[i]; + + double dx = p->x - e->x; + double dy = p->y - e->y; + double dz = p->z - e->z; + + double r = sqrt(dx*dx+dy*dy); + fprintf(f, "%f,%f,%f\n",p->m,r,dz); + } + fclose(f); +*/ + + e_init = reb_simulation_energy(r); + //printf("Initial Energy: %e\n", e_init); + begin = clock(); + reb_simulation_integrate(r, tmax); + //printf("Final Energy: %e\n", reb_tools_energy(r)); + //printf("Conservation: %f\n", (reb_tools_energy(r) - e_init) / e_init); + //printf("Final N: %d\n", r->N); + //printf("Time Spent: %f\n", time_spent); +/* + FILE* f4 = fopen(snapshot_4,"a"); + struct reb_particle* e = &r->particles[0]; + for (unsigned int i = 1; i < r->N; i++){ + struct reb_particle* p = &r->particles[i]; + struct reb_orbit o = reb_tools_particle_to_orbit(r->G, *p, *e); + + double dx = p->x - e->x; + double dy = p->y - e->y; + double dz = p->z - e->z; + + double r = sqrt(dx*dx+dy*dy); + fprintf(f4, "%f,%f,%f,%f,%f,%f\n",p->m,r,dz,o.a,o.e,o.inc); + } + fclose(f4); +*/ + reb_simulation_free(r); +} + +void heartbeat(struct reb_simulation* r){ + if (reb_simulation_output_check(r, 0.01*2.*M_PI)){ + reb_simulation_output_timing(r, tmax); + } +/* + if (snap2 && r->t > snap2_time){ + struct reb_particle* e = &r->particles[0]; + snap2 = 0; + FILE* f = fopen(snapshot_2,"a"); + for (unsigned int i = 1; i < r->N; i++){ + struct reb_particle* p = &r->particles[i]; + + double dx = p->x - e->x; + double dy = p->y - e->y; + double dz = p->z - e->z; + + double r = sqrt(dx*dx+dy*dy); + fprintf(f, "%f,%f,%f\n",p->m,r,dz); + } + fclose(f); + } + + if (snap3 && r->t > snap3_time){ + struct reb_particle* e = &r->particles[0]; + snap3 = 0; + FILE* f = fopen(snapshot_3,"a"); + for (unsigned int i = 1; i < r->N; i++){ + struct reb_particle* p = &r->particles[i]; + + double dx = p->x - e->x; + double dy = p->y - e->y; + double dz = p->z - e->z; + + double r = sqrt(dx*dx+dy*dy); + fprintf(f, "%f,%f,%f\n",p->m,r,dz); + } + fclose(f); + } +*/ + //if (r->t < 10. * 2 * M_PI){ + + if (reb_simulation_output_check(r, 0.1)){ + clock_t end = clock(); + double time_spent = (double)(end - begin) / CLOCKS_PER_SEC; + FILE* f = fopen(title,"a"); + fprintf(f, "%f,%e,%d,%e\n", r->t, fabs((reb_simulation_energy(r) - e_init) / e_init), r->N-1, time_spent); + //fprintf(f, "%f,%d\n", r->t, r->N-1); + fclose(f); + } + + //} + /* + if (r->t > 10. * 2 * M_PI && r->t < 10. * 2 * M_PI){ + if (reb_output_check(r, 10. * 2.*M_PI)){ + FILE* f = fopen(title,"a"); + fprintf(f, "%f,%e,%d\n", r->t, fabs((reb_tools_energy(r) - e_init) / e_init), r->N-1); + fclose(f); + } + } + + if (r->t > 100. * 2 * M_PI){ + if (reb_output_check(r, 100. * 2.*M_PI)){ + FILE* f = fopen(title,"a"); + fprintf(f, "%f,%e,%d\n", r->t, fabs((reb_tools_energy(r) - e_init) / e_init), r->N-1); + fclose(f); + } + } + */ +} diff --git a/examples/ms_close/Makefile b/examples/ms_close/Makefile new file mode 100644 index 000000000..cd123bddb --- /dev/null +++ b/examples/ms_close/Makefile @@ -0,0 +1,35 @@ +export OPENGL=0# Set this to 1 to enable OpenGL +export SERVER=1# Set this to 1 to enable the visualization web server +include ../../src/Makefile.defs + +# CCPROBLEM is defined in Makefile.defs to allow for +# a compact cross platform Makefile +.PHONY: all librebound +all: problem.c librebound + @echo "Compiling $< ..." + $(CCPROBLEM) + @echo "" + @echo "Compilation successful. To run REBOUND, execute the file '$(EXEREBOUND)'." + @echo "" + +librebound: + @echo "Compiling shared library $(LIBREBOUND) ..." + $(MAKE) -C ../../src/ + @-$(RM) $(LIBREBOUND) + @$(LINKORCOPYLIBREBOUND) + @echo "" + +clean: + @echo "Cleaning up shared library $(LIBREBOUND) ..." + $(MAKE) -C ../../src/ clean + @echo "Cleaning up local directory ..." + @-$(RM) $(LIBREBOUND) + @-$(RM) $(EXEREBOUND) + +rebound_webgl.html: problem.c + @echo "Compiling problem.c with emscripten (WebGL enabled)..." + emcc -O3 -I../../src/ ../../src/*.c problem.c -DSERVERHIDEWARNING -DOPENGL=1 -sSTACK_SIZE=655360 -s USE_GLFW=3 -s FULL_ES3=1 -sASYNCIFY -sSINGLE_FILE -sEXPORTED_RUNTIME_METHODS="callMain" --shell-file ../../web_client/shell_rebound_webgl.html -o rebound_webgl.html + +rebound_console.html: problem.c + @echo "Compiling problem.c with emscripten (WebGL disabled)..." + emcc -O3 -I../../src/ ../../src/*.c problem.c -DSERVERHIDEWARNING -sSTACK_SIZE=655360 -sASYNCIFY -sSINGLE_FILE -sEXPORTED_RUNTIME_METHODS="callMain" --shell-file ../../web_client/shell_rebound_console.html -o rebound_console.html diff --git a/examples/ms_close/problem.c b/examples/ms_close/problem.c new file mode 100644 index 000000000..686f81d93 --- /dev/null +++ b/examples/ms_close/problem.c @@ -0,0 +1,188 @@ +#include +#include +#include +#include "rebound.h" +#include + +void heartbeat(struct reb_simulation* r); + +double e_init; // initial energy +double err; +double emax = 0.0; +double tmax = 3000. * 2. * M_PI; + +char title_stats[100] = "peri_stats"; + +const double k = 1.; // Constants for the Harmonic Oscillator +const double m = 10000000.; + +void derivatives(struct reb_ode* const ode, double* const yDot, const double* const y, const double t){ + const double omega = sqrt(k/m); + //struct reb_orbit o = reb_orbit_from_particle(ode->r->G, ode->r->particles[1], ode->r->particles[0]); + //double forcing = sin(o.f); + yDot[0] = y[1]; + yDot[1] = -omega*omega*y[0];// + forcing; +} + +int main(int argc, char* argv[]){ + + // Initialize masses + struct reb_simulation* r = reb_simulation_create(); + struct reb_particle star = {0}; + star.m = 1; + + // Jupiter + double jm = 9.55e-4; + double ja = 5.2; + double je = 0.05; + double jmse = 1.-je; + + double peri = 0.1; + int integrator = 0; + if (argc == 3){ + peri = pow(10, -1. * atof(argv[2])/2.); + integrator = atoi(argv[1]); + } + + // Saturn + double sm = 2.857e-4; + double sa = 9.58; + double se = (sa - peri) / sa; + double omse = 1. - se; + double si = M_PI / 2.; + + char grators[10][10] = {"TRACE", "WHFAST", "MERCURIUS", "WHFASTr", "MERCURIUSr", "TRACEf"}; + switch(integrator){ + case 0: + //printf("Using TRACE\n"); + r->integrator = REB_INTEGRATOR_TRACE; + r->dt = 0.15 * 2 * M_PI; + r->ri_trace.peri_crit_distance=2.; + break; + case 1: + //printf("Using WHFAST %f\n", se); + r->integrator = REB_INTEGRATOR_WHFAST; + r->dt = 0.15 * 2 * M_PI; + break; + case 2: + r->integrator = REB_INTEGRATOR_MERCURIUS; + r->dt = 0.15 * 2 * M_PI; + break; + case 3: + r->integrator = REB_INTEGRATOR_WHFAST; + r->dt = 2*M_PI*sqrt((1-se)*(1-se)*(1-se)*sa*sa*sa/((1+se)*r->G*star.m))/50.; + break; + case 4: + r->integrator = REB_INTEGRATOR_MERCURIUS; + r->dt = 2*M_PI*sqrt((1-se)*(1-se)*(1-se)*sa*sa*sa/((1+se)*r->G*star.m))/15.; + break; + case 5: + r->integrator = REB_INTEGRATOR_TRACE; + r->dt = 0.15 * 2 * M_PI; + r->ri_trace.peri_crit_fdot=17.; + break; + } + + //r->integrator = REB_INTEGRATOR_TRACE; + //r->dt = 2*M_PI*sqrt((1-se)*(1-se)*(1-se)*sa*sa*sa/((1+se)*r->G*star.m))/15.;//0.15 * 2 * M_PI; + //r->dt = 0.15 * 2 * M_PI; + //r->ri_trace.r_crit_hill = 0.; + //r->ri_ias15.adaptive_mode=2; + //r->ri_trace.peri_crit_fdot=0.; + + //struct reb_ode* ho = reb_ode_create(r,2); // Add an ODE with 2 dimensions + //ho->derivatives = derivatives; // Right hand side of the ODE + //ho->y[0] = 1; // Initial conditions + //ho->y[1] = 0; + reb_simulation_add(r, star); + reb_simulation_add_fmt(r, "m a e", jm, ja, je); + + //r->dt = orb.P / 15.12345; + reb_simulation_add_fmt(r, "primary m a e inc omega", star, sm, sa, se, si, M_PI/2.); + + + struct reb_orbit orb = reb_orbit_from_particle(r->G, r->particles[2], r->particles[0]); + + reb_simulation_move_to_com(r); // This makes sure the planetary systems stays within the computational domain and doesn't drift. + r->heartbeat = heartbeat; + r->exact_finish_time=0; + e_init = reb_simulation_energy(r); + //system("rm -rf ho.txt"); + //system("rm -rf energy_wh.txt"); + + //FILE* f = fopen("energy_trace_best.txt","a"); + //fclose(f); + + clock_t begin = clock(); + //while(r->tt, ho->y[0]); + //} + clock_t end = clock(); + double time_spent = (double)(end - begin) / CLOCKS_PER_SEC; + //printf("\n%e\n", fabs((reb_simulation_energy(r) - e_init) / e_init)); + printf("\n%s %e %e %e\n", grators[integrator], peri, emax, time_spent); + + FILE* f = fopen(title_stats,"a"); + fprintf(f, "%s,%e,%e,%e\n", grators[integrator], peri, emax, time_spent); + fclose(f); + reb_simulation_free(r); +} + + +void heartbeat(struct reb_simulation* r){ + //fprintf(f,"%e %e\n",r->t, (e - e_init) / e_init); + //if (reb_simulation_output_check(r, 10.*2.*M_PI)){ + // reb_simulation_output_timing(r, tmax); + //} + + double e_curr = fabs((reb_simulation_energy(r) - e_init) / e_init); + //printf("%f %e\n", r->t, e_curr); + if (emax < e_curr){ + emax = e_curr; + } +/* + if (reb_simulation_output_check(r, 0.15 * 2.*M_PI)){ + // Once per 4 days, output the relative energy error to a text file + + double err = fabs((reb_simulation_energy(r) - e_init) / e_init); + if (err > emax){ + emax = err; + } + FILE* f = fopen("energy_wh.txt","a"); + + // rotate whole simulation to rotating frame + //reb_simulation_irotate(r, r1); + + struct reb_particle* sun = &r->particles[0]; + struct reb_particle* jup = &r->particles[1]; + struct reb_particle* sat = &r->particles[2]; + + struct reb_orbit o = reb_orbit_from_particle(r->G, *sat, *sun); + struct reb_orbit oj = reb_orbit_from_particle(r->G, *jup, *sun); + + const double dvx = r->particles[0].vx - r->particles[2].vx; + const double dvy = r->particles[0].vy - r->particles[2].vy; + const double dvz = r->particles[0].vz - r->particles[2].vz; + const double dv = sqrt(dvx*dvx + dvy*dvy + dvz*dvz); + + const double dvxj = r->particles[0].vx - r->particles[1].vx; + const double dvyj = r->particles[0].vy - r->particles[1].vy; + const double dvzj = r->particles[0].vz - r->particles[1].vz; + const double dvj = sqrt(dvxj*dvxj + dvyj*dvyj + dvzj*dvzj); + + const double vcirc = sqrt(r->G * r->particles[0].m / o.a); + const double vcircj = sqrt(r->G * r->particles[0].m / oj.a); + double e = (reb_simulation_energy(r) - e_init) / e_init; + fprintf(f,"%e %e %e %e %e %e %e %e %e %e %e %f %f %f %f\n",r->t, e, r->particles[0].x, r->particles[0].y, r->particles[0].z, r->particles[1].x, r->particles[1].y, r->particles[1].z, r->particles[2].x, r->particles[2].y, r->particles[2].z, o.e, o.inc, dvx/vcirc, dvxj/vcircj); + + //double e = reb_tools_energy(r); + //fprintf(f,"%e %e\n",r->t, (e - e_init) / e_init); + fclose(f); + + //reb_integrator_synchronize(r); + //printf("\n Jacobi: %e\n", e); + } +*/ + +} diff --git a/ipython_examples/HybridIntegrationsWithTrace.ipynb b/ipython_examples/HybridIntegrationsWithTrace.ipynb new file mode 100644 index 000000000..854208753 --- /dev/null +++ b/ipython_examples/HybridIntegrationsWithTrace.ipynb @@ -0,0 +1,212 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Hybrid integrations with TRACE\n", + "REBOUND comes with several integrators, each of which has its own advantages and disadvantages. TRACE is a hybrie time-reversible integrator that is based on the algorithm described by Hernandez & Dehnen (2023). It uses a symplectic Wisdom-Holman integrator when particles are far apart from each other and switches over to a high order integrator during close encounters. Specifically, TRACE uses the efficient WHFast and BS integrators internally." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's start out by showcasing the problem with traditional fixed timestep integrators such as WHFast. We setup a simulation of the outer solar system and increase the masses of the planets by a factor of 50. " + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbIAAAHDCAYAAABBDZ94AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB5f0lEQVR4nO3dd3yT5fo/8E920pXupnvQAbLLatkbAfFwQMAtyEFUXMBRwQHiguNCRQ84AQWBAwoqIFDKKtIyWsruoHTvkSbpyH5+f/hrvlQKdCR5Mq7365UXNkmfXIaST+/7uZ775jAMw4AQQgixU1y2CyCEEEK6goKMEEKIXaMgI4QQYtcoyAghhNg1CjJCCCF2jYKMEEKIXaMgI4QQYtcoyAghhNg1CjJCCCF2jYKMEEKIXbOrICstLcWjjz4KHx8fSCQS9O7dG+fOnTM9zjAMVqxYgcDAQEgkEowfPx65ubksVkwIIcTS+GwX0F5yuRzDhg3DmDFj8Mcff8DPzw+5ubnw8vIyPeeDDz7A559/js2bNyMyMhJvvvkmJk2ahKtXr0IsFt/1NYxGI8rKyuDu7g4Oh2PJ/x1CCCF3wDAMVCoVgoKCwOXeZczF2IlXX32VGT58+G0fNxqNjEwmYz788EPTffX19YxIJGK2bdvWrtcoLi5mANCNbnSjG91s5FZcXHzXz267GZH99ttvmDRpEmbNmoXjx48jODgYzz77LBYsWAAAyM/PR0VFBcaPH2/6HqlUiiFDhiA1NRUPPvjgLcfUaDTQaDSmr5n/vxFAcXExPDw8LPx/RAgh5HaUSiVCQ0Ph7u5+1+faTZDduHED69evx5IlS/Daa6/h7NmzeOGFFyAUCvHEE0+goqICABAQENDq+wICAkyP/d3q1auxatWqW+738PCgICOEEBvQntM8dtPsYTQaER8fj/fffx/9+/fHU089hQULFmDDhg2dPuby5cuhUChMt+LiYjNWTAghxBrsJsgCAwNxzz33tLqvR48eKCoqAgDIZDIAQGVlZavnVFZWmh77O5FIZBp90SiMEELsk90E2bBhw5Cdnd3qvpycHISHhwMAIiMjIZPJkJycbHpcqVTi9OnTSExMtGqthBBCrMduzpEtXrwYQ4cOxfvvv4/Zs2fjzJkz+Prrr/H1118D+Gse9aWXXsK7776LmJgYU/t9UFAQpk+fzm7xhBBCLMZugmzQoEHYvXs3li9fjrfffhuRkZH49NNP8cgjj5ie88orr6CxsRFPPfUU6uvrMXz4cBw4cKBd15ARQgixTxympeecQKlUQiqVQqFQ0PkyQghhUUc+j+3mHBkhhBDSFgoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF3js10AIcT6GIaBVqsFAHA4HAAAl8sFn08fCcT+0E8tIXZOpVKhsrISCoWi3TeVSgU/Pz9UVVWZjtO3b1/k5OTAzc2tzZurqyuCgoIgk8kQEhKC0NBQSKVSUxASwhYKMkLshFarRV5eHrKyslrdAMDb2xs3btwwPVcikcDT0xMeHh6QSqWQSqUICAgw/beHhwdcXFxMIcQwDDgcDpRKJRoaGtDY2AiVSoXGxkY0NDSgqKgIDQ0NqK2tRX19vel13NzcEBoaagq2lj8jIyPh7+9PIUesgoKMEBvDMAyqqqpuCay8vDwYDAYAQGhoKOLi4vDggw8iNjYWgYGB8PT0NIWUQCCwWH0NDQ0oLi5GcXExSkpKTH8ePHgQ5eXlMBqNcHFxgVgsRkJCAiIjIzF48GD06tULQqHQYnUR58VhGIZhuwhboVQqIZVKoVAo4OHhwXY5xEkwDIPc3FykpaXh1KlTUKlUSEtLAwC4uLige/fu6NGjB7p3747u3bsjNjYWbm5uLFfdNr1ej7KyMuTn5yMrKwvl5eU4evQolEolRCIR+vfvj8GDB2Pw4MHo0aMHuFzqNyNt68jnMQXZTSjIiLWUl5cjNTUVR48eRWpqKmprayEQCNC/f3+MHj0aMTExiIuLQ3BwsN1/2BuNRuTk5ODs2bM4ffo0MjIy0NzcDDc3NwwYMACDBw/GoEGDEB0dTVORxISCrJMoyIiltIy6Dh48iEOHDuHixYuIj48Hh8PB0KFDkZiYiPj4eEgkErZLtTi9Xo8rV67gzJkzOHv2LDIzM6HT6eDl5YXBgwdj9OjRGDp0KP0bdHIUZJ1EQUbMiWEYXLlyBb/++isOHTqEgoICuLq6YvTo0Zg4cSJGjx4NT09PtstknUajwYULF3DmzBmcOXMG169fh16vx+DBgzFt2jSMGDECLi4ubJdJrIyCrJMoyIg51NXV4ZdffsG2bdtQVVUFkUiE8ePHY9KkSRg6dChEIhHbJdq02tpaHDlyBElJSaipqUFtbS0mTpyIf/7zn+jRowdNPzoJCrJOoiAjnWU0GvHnn39i27ZtOHDgABiGweTJk/Hggw9i2LBh4PF4bJdolyoqKvDbb7/ht99+Q2VlJWJiYjB9+nRMnjwZ7u7ubJdHLIiCrJMoyEhHlZWV4X//+x927NiBkpISxMTE4KGHHsLMmTPh7e3NdnkOw2g0Ii0tDbt370ZKSgp4PB7Gjx+P6dOno1+/fjRKc0AUZJ1EQUbaw2Aw4NixY9i8eTOOHTsGiUSC+++/Hw8++KCpgYNYTm1tLX7//Xf8+uuvKCkpQUJCAu6//36MGzfO7js8yf+hIOskCjJyJxqNBnv27MFnn30Gg8EAf39/PPzww5g2bZrNXtflyIxGI86ePYtt27YhNTUV4eHhmDt3LiZNmkRTuQ6AgqyTKMhIW3Q6Hf73v//h448/hru7O3r27ImFCxeib9++bJdG/r+rV6/i+++/R0pKCoKDg/HEE09gypQpFl3hhFgWBVknUZCRmxmNRvz666/48MMPkZ+fj+nTp+Pll19GVFQU26WR28jJycH333+Po0ePQiaT4bHHHsP9999PS2PZIQqyTqIgI8Bf138lJSXhP//5D65evYoJEybg1VdfRc+ePdkujbTTjRs3sHHjRiQlJSEwMBBPPPEEpk2bRlOOdoSCrJMoyMiff/6J1atXIz09HUOHDsXy5csxcOBAtssinVRYWIhdu3Zh165diIiIwPPPP4+hQ4eyXRZpBwqyTqIgc14ZGRlYs2YNUlJS0K9fPyxfvhwjRoygDkQHkZ2djc8//xwqlQohISFYunQpfHx82C6L3EFHPo+pV5U4tZycHMybNw9Tp05FVVUVNm7ciP3792PkyJEUYg4kLi4OX3zxBebOnYuLFy/i4Ycfxt69e0G/xzsGGpHdhEZkzkOlUuGjjz7Cnj174ObmhiVLlmD69Ol0DsUJKBQKrFu3Dn/88QcGDRqEV155BUFBQWyXRf6GphY7iYLMORw7dgyLFy9GeHg4JkyYgAULFlBXmxM6ffo0/vOf/0ChUOCpp57CrFmz6IJqG0JTi4S0QaVSYcmSJZg1axYiIyOxbt06LFq0iELMSQ0ZMgRbt27FtGnTsG7dOixcuBD5+flsl0U6gYKMOIWjR49ixIgR+OWXX/Dhhx/i559/RmhoKNtlEZZJJBK89NJLWL9+PRobGzF37lxs3LgROp2O7dJIB1CQEYfWMgqbPXs2oqKikJKSgrlz51IjB2mld+/e2LRpEx599FFs3LgRTz75JK5cucJ2WaSd6BzZTegcmWM5evQoFi9ejPr6eqxatQqPP/44BRi5q+vXr+P9998HAIwZMwaPPPIInTtjAZ0jI05NpVJh8eLFmD17Nrp164aUlBQ88cQTFGKkXaKjo/H1118jISEB33zzDV5++WXI5XK2yyJ3QEFGHMrZs2fx9NNPY/fu3fjoo4+wa9cuOhdGOozP5+Opp57Cxx9/jOrqajz//PO4fv0622WR26AgIw6BYRhs2LAB9913H+RyOY4dO0ajMNJlgwYNwkcffQSJRIJFixYhJSWF7ZJIGyjIiN1TKpWYO3cuXn/9dSxcuBC//vorIiIi2C6LOAh/f3+sW7cOCQkJeOONN7BlyxZaEcTG8NkugJCuyMrKwvLly5GZmYkff/wRU6ZMYbsk4oDEYjFWrlyJ8PBwfPPNNygoKMArr7xC1yDaCBqREbt1+PBhTJgwATqdDkePHqUQIxbF5XLx5JNPYuXKlTh+/DhefPFF1NXVsV0WAQUZsUMt58PmzJmD4cOHY/v27TSVSKxm7NixWLduHSorK7Fw4UJqArEBdhtka9asAYfDwUsvvWS6T61WY9GiRfDx8YGbmxtmzpyJyspK9ookZqfT6bB06VIsX74czz33HLZs2QI3Nze2yyJOpnv37vj666/h5eWFjz76CBkZGWyX5NTsMsjOnj2Lr776Cn369Gl1/+LFi/H7779j586dOH78OMrKyjBjxgyWqiTmVl9fj1mzZmHLli1Yt24dVq1aRavVE9b4+vris88+g4eHB5YvX44zZ86wXZLTsrsga2howCOPPIJvvvkGXl5epvsVCgW+++47fPLJJxg7diwGDBiAjRs34tSpU0hLS2OxYmIOBQUFmDBhAi5evIjdu3fj0UcfZbskQiCRSPDee+9h4MCBeP3113Hq1Cm2S3JKdhdkixYtwtSpUzF+/PhW96enp0On07W6v3v37ggLC0Nqamqbx9JoNFAqla1uxPZcuXIFEyZMQExMDJKTkzFs2DC2SyLERCAQYNWqVRg6dChWrFiBEydOsF2S07GrINu+fTsyMjKwevXqWx6rqKiAUCiEp6dnq/sDAgJQUVHR5vFWr14NqVRqutEKELbn3LlzmDJlCmQyGT7//HNERkayXRIht+Dz+VixYgVGjRqFVatW4ciRI2yX5FTsJsiKi4vx4osvYuvWrRCLxWY55vLly6FQKEy34uJisxyXmMfJkydx//33o3v37ti7dy98fX3ZLomQ2+LxeHj99dcxfvx4vPvuuzh06BDbJTkNu7kgOj09HVVVVYiPjzfdZzAYcOLECXzxxRc4ePAgtFot6uvrW43KKisrIZPJ2jymSCSCSCSydOmkEw4ePIjHH38ciYmJ+Omnn+Di4sJ2SYTcFZfLxbJlyyAQCLBmzRrodDpMnTqV7bIcnt0E2bhx43Dp0qVW982bNw/du3fHq6++itDQUAgEAiQnJ2PmzJkAgOzsbBQVFSExMZGNkkkn/fLLL1iwYAHuvfdefP/99/TLBrErHA4HS5cuhUAgwLZt2yASiW45p0/My26CzN3dHb169Wp1n6urK3x8fEz3z58/H0uWLIG3tzc8PDzw/PPPIzExEQkJCWyUTDph586dWLZsGWbOnIn//ve/4PPt5keUEBMOh4MXXngBH330Ef7zn//Azc2NPocsyG7OkbXH2rVrcd9992HmzJkYOXIkZDIZfvnlF7bLIu30888/41//+hceeOABbNiwgUKM2DUOh4MlS5YgMTERH330Ea5du8Z2SQ6Ldoi+Ce0QzZ59+/bh0UcfNYUYXehMHIVWq8WKFSuQn5+PL7/8kpqW2ol2iCZ25fDhw3j88ccxdepUrF+/nkKMOBShUIhXX30VPB4Pb731FrRaLdslORwKMsKqM2fOYPXq1Rg7diy+//57mk4kDsnLywurVq1Cfn4+Pv30U9rPzMzoU4OwJj8/H7NmzUJcXBx++OEHp97biWEY1NbWoqysDFVVVZDL5aYtQiorK6HRaKDVaqHRaG773xKJBFqtFtHR0airq4O7uzvc3Nzg7u5uurm5ucHDw6PV156envD29qaRsIXFxMRg6dKlWL16NaKjo2kdWDOiICOsqK+vx4wZM+Dp6Ynt27dDIpGwXZLFGY1GVFdXo6SkBKWlpa1uZWVlUKvVAAAXFxdwOBx4eXmhW7duaG5uhlAohLu7O3x9fSEUCiEWiyEUCk3XQgqFQvB4PBiNRgB/vb8qlQoNDQ1QqVQoLCw0fa1UKluNCAIDAyESiSAQCBAREYGIiAiEh4cjIiICvr6+4HA4rLxfjmjs2LHIy8vDhg0bEBER0eq6WNJ51OxxE2r2sA6tVov7778fly9fxrFjxxAdHc12SWan1+tx48YN5OTkIDs7Gzk5OdBqtbhx4waAvy6clclkCAkJQVBQEIKDgxESEoLg4GD4+/tbNNgZhkFTU5Mp2ORyOSorK5Gfn4/CwkIUFBSgubkZwF+h2hJqNwecu7u7xepzdEajEa+//jqysrLw5ZdfIigoiO2SbFJHPo8pyG5CQWZ5DMNg4cKF2LlzJ/bt24ehQ4eyXZJZ6HQ6ZGVlITMzExkZGbh48SKCg4NRWFiI8PBwxMbGIi4uDuHh4QgODkZAQIDNng9kGAbV1dWmUCsoKEBhYSGKi4uh1+sBADKZDP3790dUVBTi4+Nvu3oOaVtDQwMWLVoEoVCIzz77jFauaQMFWSdRkFnemjVr8M4772Djxo2YPXs22+V0mk6nw9WrV5GRkYHMzExcunQJGo0GLi4u6NOnD/r374++ffsiJibGbGuDsk2v16O0tNQUagUFBUhPT4fBYEBQUBDi4+MRHx+PXr160Wos7VBYWIjnn38e8fHxWLlyJU3h/g0FWSdRkFnWjh078OSTT2LFihV49dVX2S6nwxQKBU6dOoWTJ08iLS0N7u7uaG5uRt++fdGvXz/Ex8cjJibGqZommpqacOnSJWRkZCAjIwNVVVUQCATo2bOnKdhCQkLoQ/o2UlNTsWLFCjzxxBO0x97fUJB1EgWZ5fz555+47777MGvWLHz11Vd288FWV1eHpKQkpKSk4Pz58zAajejRowdGjBiBoUOHIiYmBlwuXcUC/DUlWVZWZgq1y5cvQ6vVwtfXF/Hx8UhISECfPn0gEAjYLtWmbN26FZs2bcJ7772HwYMHs12OzaAg6yQKMssoLCzEiy++CLVajd9++83m2+ybm5tx7NgxHDx4EKdPnwaXy8XEiRPRr18/DBs2jFZmaCetVosrV66Ygq2pqQkajQYTJ05EQkIC4uLi7OYXGktiGAZr167FyZMn8dVXX8HPz4/tkmwCBVknUZCZX3NzM8aOHYv6+nqkpKTYbAgYDAacPn0aBw8exLFjx6BWq9G3b19MnjwZY8eOhVQqZbtEu1daWoqUlBSkpaWhuLgYoaGhmDBhAkaNGgU3Nze2y2OVSqXC008/jaCgIPznP/+hUT4oyDqNgsz8XnzxRfzwww84duwY+vbty3Y5t6ipqcHu3btx5MgR3LhxAxEREZg8eTImTZqEwMBAtstzSAzD4OLFi0hKSsKZM2fA5XIxbNgwTJw4EbGxsU47Srtw4QL++9//YsqUKfjHP/7Bdjms68jnsW32/xKH8PPPP+OXX37BF198YVMhxjAMzp8/j507d+Lo0aMQCoWYPHky3nzzTfTo0cNpP0ithcPhoG/fvujbty/q6+tx9OhRJCUl4dixYwgLC8OECRMwcuRIpxul9e3bF71798amTZuQmJgIf39/tkuyGzQiuwmNyMynrKwM8fHxGDt2LLZu3WoT4aDX63HgwAEcOXIEKSkpCA8PxwMPPID77rvP6T40bU3LKO3QoUM4e/YsuFwupkyZgrFjxyIkJITt8qymsbERCxcuRFRUFFatWmUT/27YQlOLnURBZh5GoxH33Xef6Torb29vVuvRaDT47bffsHnzZpSXl+Pee+/F9OnTMXDgQKf+oLBV9fX1SE5ORnJyMqqrqzFo0CBMnz4dsbGxbJdmFampqXj77bexbNkyjBo1iu1yWENB1kkUZObxxRdfYOnSpdi/fz/GjRvHWh2NjY3YtWsXtmzZgvr6ekycOBFz585FTEwMazWR9tPr9UhJScGePXug1WoRGhqKhx56CJGRkWyXZnHvvfceLl26hK+//tppP4soyDqJgqzrrl69ioSEBPzrX//CJ598wkoNzc3N2LZtG3bs2IH6+npMmzYNTzzxBEJDQ1mph3QNwzA4ffo0duzYgbKyMgwdOhSzZ8926GYcuVyOBQsWYOjQoViyZAnb5bCCgqyTKMi6RqvVYvjw4dBqtUhNTbX6ivZ6vR579uzBhg0boFQqMXfuXMycORMBAQFWrYNYhsFgwPHjx7Fr1y7I5XKMHTsWDzzwALy8vNguzSIOHjyITz/9FO+//z769+/PdjlWR0HWSRRkXfP666/js88+w8mTJ9GvXz+rvS7DMEhOTsYXX3yBoqIiTJkyBc8++yytKu6gtFotDh06hN27d0Oj0Zja1V1dXdkuzawYhsGyZctQVVWFDRs2ON36lRRknURB1nkpKSmYMGEC3n33Xfz73/+22uteuXIFP/74Iw4ePIihQ4fixRdfdJqmAGfX1NSEvXv3Yu/evfD398f48eMxYcIEh1rrsrS0FM888wxmzJiBuXPnsl2OVVGQdRIFWcckJSVhxYoVKCkpQXV1Nbp374709HSrfJAolUqsW7cOO3fuRM+ePfHCCy9gyJAhFn9dYnsUCgX27NmDgwcPIiwsDPPnz3eohp5du3bhf//7Hz777DOHPi/4dxRknURB1n5PPvkkNm7ceMv98+fPx7fffmux1zUajfjtt9+wdu1a6HQ6LFq0CA8++KBD/RZOOicvLw/fffcd8vPzMWbMGDz00EMOsQGoVqvFv/71L/Ts2dMud43oLAqyTqIga5+kpCRMnDjxto8fPnzYIm33eXl5ePvtt5GZmYkpU6Zg6dKlNrt2I2GH0WhEcnIytm/fDplMhqlTpyIxMdHurxc8dOgQPvvsM3z66acONdq8EwqyTqIga5/ExESkpaXd9vGEhASkpqaa7fUMBgO+++47bNq0CVFRUXjhhRdouwtyR0qlElu3bsXJkycxePBgzJs3z67/TRsMBjz33HPw8vLCe++9Z/fB3B4d+TymJZZJh5WVlXXp8Y64ceMGHn30UXz55ZeYM2cOvv/+ewoxclceHh545pln8MILL+Dq1at49dVXce7cObbL6jQej4e5c+fiwoULOH/+PNvl2BwKMtJhd2trN0fbu9FoxKZNmzBr1iw0NDTgxx9/xIsvvmjze5kR2zJkyBB88MEHiImJwdq1a7F+/Xo0NjayXVanDB48GD179sT3338PmkhrjYKMdNjbb799x8fffffdLh2/vLwcy5YtwyeffIIHH3wQu3btQp8+fbp0TOK8pFIpFi9ejGeeeQYZGRn4+OOPkZOTw3ZZHcbhcDBv3jzk5+fj2LFjbJdjUyjISIdNmDAB8+fPb/Ox+fPnd6nR4/jx45gxYwZycnKwceNGvPzyy053ISgxPw6Hg+HDh2PNmjXgcDh47733sH//frsb2fTo0QOJiYn47bffoNfr2S7HZlCQkU759ttvcfjwYSQkJCAsLAwJCQk4fPhwp1vv9Xo9PvroIzz77LMYMGAANm/ejAEDBpi5auLsfHx88Nprr2Hy5MnYtm0bPv/8czQ1NbFdVoc89thjuHHjBo4ePcp2KTaDuhZvQl2L7CgvL8fSpUtx+fJl/Pvf/8Zjjz3mFF1ZhF3p6en4+uuv4ebmhhdeeAHh4eFsl9Ru77//PgoKCrBhwwZwuY45HqGuRWI3UlNT8eKLL6K6uhpbtmzB448/TiFGrGLAgAF45513IJFIsGrVKhw/fpztktpt1qxZKC8vx8mTJ9kuxSZQkBHW/PTTT1iwYAECAwOpoYOwwt/fHytWrMCwYcOwd+9ebN26FUajke2y7iomJgb9+/e3y/N8lkBBRqxOr9fj7bffxjvvvINHHnkEn3zyCaRSKdtlESclFAoxf/58TJo0CUlJSfj000+hVqvZLuuuZs6ciaysLFy7do3tUlhHQUasSqFQYMGCBfjf//6Hd955B8uXL6d1EolNGD9+PJYuXYrs7Gy8++67qKurY7ukO+rTpw9kMhkOHjzIdimsoyAjVlNYWIhFixYhKysLGzduxAMPPMB2SYS00rt3b7z55ptobm7G+vXrUVpaynZJt8XhcDBx4kT8+eefaGhoYLscVlGQkbtiGKbL8/BXr17FQw89BJ1Oh//9738YNGiQmaojxLxCQkKwYsUKqNVqrFmzBgUFBWyXdFvjxo2D0Wh0+lZ8CjJyV10NsTNnzuCxxx5DYGAg1q9fj9DQUDNVRohlSKVSLFu2DP7+/vjggw+Qm5vLdkltkkqlSEhIwIEDB5y66YOCjNxRyz+OzrbEJycnY/78+ejduzd++OEHeHt7m7M8QizG1dUVL7/8MsLCwvDRRx/h6tWrbJfUpkmTJqGkpMSpmz4oyMgdMQzT6RDbs2cPnnvuOYwZMwZff/01XF1dzVwdIZYlFouxZMkSxMbGYu3atbhw4QLbJd2ipenj0KFDbJfCGgoycltdGY3t2LEDn3/+OR588EGsXbuWVq0ndksoFOLFF19Enz598Pnnn+Ps2bNsl9RKS9PHyZMnnbbpg4KM3FZnR2O7du3CG2+8gVGjRmHFihXUXk/sHp/Px7PPPovBgwdj/fr1NreihrM3fVCQkTZ19sTxL7/8gtdeew0PP/wwVqxYQctNEYfB4/GwYMECDB8+HL///vsdd0m3Nk9PTyQkJODgwYNO2fRBQUbuqCNBtHv3bixbtgxz5szBypUrKcSIw+FyuZg3bx7i4uLw3Xff2dRuzZMmTUJxcbFTNn1QkJE2dXRa8cCBA/jkk08wa9YsrFq1ymFX5CaEw+Fg7ty5iI+Px+7du3H9+nW2SwLg3E0f9GlDbtHRJo+0tDS89NJLGDZsGN555x0KMeLwuFwuFixYAHd3d3z55Zeorq5muyRwOBzce++9KCwsRHNzM9vlWBV94pBbdGQ0lpWVhYULF2LIkCEUYsSp8Pl8PPPMM3BxccG6detsIjyGDRsGrVaLzMxMtkuxKvrUIa10ZDRWWlqKefPmITw8HP/9738hEAgsXR4hNsXNzQ3PP/88FAoFvvrqK9a3gPH394dYLMaZM2dYrcPaKMhIK+0djdXX12PevHkQiUT47rvv6GJn4rRkMhmefvppXLt2DTt27GC7HPTv3x+ZmZlO1b1IQUZMjEYjDAbDXYNMq9Xi7bffRn19PTZt2gQ/Pz8rVUiIberRowcefvhhHDlyBMeOHWO1lvj4eCiVSuTl5bFahzXx2S6A2I72TIswDIMVK1Zg79692LZtGyIiIixfGCF2YNSoUSgvL8e2bdvg5+eHnj17slJHXFwcXFxccP78eURHR7NSg7XRiIyYGI3GuzZrbNmyBTt27MB7772HAQMGWKkyQuzD7Nmz0bNnT3z11VcoLy9npQYej4e+ffsiIyODlddnAwUZAfB/o7E7BVlaWhpWrVqFuXPnYtasWdYqjRC70dKW7+3tjXXr1rG29mF8fDxyc3OhVCpZeX1royAjAP4KMg6Hc9vzY6WlpXj22WcxZMgQvP7661aujhD7IZFI8Pzzz8NgMGDr1q2sdDL269cPDMPY5Gr9lkBBRgDceVpRr9fjhRdeQJ8+ffDFF1+Az6dTq4TciY+PD5588kmUlJSwstKGt7c3IiIinGZ6kYKM3HVa8cMPP0RmZiZeeukleHl5WbM0QuxWXFwcBg4ciL1796KkpMTqr+9MbfgUZOSO04rHjx/Hhg0b8Morr6Bfv37WL44QOzZ16lTIZDJs3rwZer3eqq89YMAAp2nDpyAjt51WrKqqwuLFizFq1CgsWLCAhcoIsW98Ph9z585FWVkZ/vjjD6u+dmxsrKkN39FRkDm5200rGo1GLF68GFwuF2vXrqU1FAnppJCQEEyZMgUHDhxAYWGh1V6Xx+MhISEBRUVFVntNttCnk5O73bTili1bcPLkSXz66afw8fFhqTpCHMO9996LkJAQbNq0CTqdzmqvGx4ejvT0dBgMBqu9JhsoyJxcW9OK2dnZWLlyJRYvXozhw4ezVBkhjoPH42Hu3LmoqanB77//brXXjY6Ohk6ns+pIkA0UZE6srWlFvV6PJUuWICIiAs888wxbpRHicAIDAzFt2jQcPnwYN27csMprRkVFgcfj2czmn5ZCQebEdDqdaWqxxVdffYXLly/j448/hkgkYrE6QhzP+PHjERERgc2bN0Or1Vr89YRCIcLDw5Gbm2vx12ITBZkTMxqN4PF4pq9zc3Px8ccf46mnnkJ8fDyLlRHimLhcLp544gnI5XLs2bPHKq8ZHR1NQUYck9FoBMMwpmlFhmHw/vvvY9iwYVi6dCnL1RHiuAICAjB9+nQcPXrUKtd4xcTEQKvVoqmpyeKvxRa7CbLVq1dj0KBBcHd3h7+/P6ZPn47s7OxWz1Gr1Vi0aBF8fHzg5uaGmTNnorKykqWKbdvfz4/t3r0bhw4dwsKFCyEWi9ksjRCHN2bMGPTu3Rt79uyxeEdhZGQkamtrWVldxFrsJsiOHz+ORYsWIS0tDUlJSdDpdJg4cSIaGxtNz1m8eDF+//137Ny5E8ePH0dZWRlmzJjBYtW2q6VbkcPhQKVS4Z133sF9991HXYqEWAGHw8G0adNQWFiItLQ0i75WUFAQOBwOiouLLfo6bLKb1V8PHDjQ6utNmzbB398f6enpGDlyJBQKBb777jv89NNPGDt2LABg48aN6NGjB9LS0pCQkMBG2TbLYDCYzo99+umnUKlUWLFiBctVEeI8QkNDMWDAAPzxxx8YNGgQhEKhRV5HIBBAJpPRiMwWKRQKAH+t8gwA6enp0Ol0GD9+vOk53bt3R1hYGFJTU9s8hkajgVKpbHVzBgzDgGEY8Hg85Obm4ttvv8ULL7yA4OBgtksjxKlMnToVKpUKJ06csOjrBAcHU5DZGqPRiJdeegnDhg1Dr169AAAVFRUQCoXw9PRs9dyAgABUVFS0eZzVq1dDKpWabqGhoZYu3Sa0zMlzOBy88cYbCAkJwdNPP81yVYQ4H19fXwwfPhyHDh2yaDNGaGgoSktLLXZ8ttllkC1atAiXL1/G9u3bu3Sc5cuXQ6FQmG6OPId8s5bzY0eOHMHp06fxzjvvWGxagxByZ5MmTYLBYMDhw4ct9hohISGora112M5Fuwuy5557Dnv37sXRo0cREhJiul8mk0Gr1aK+vr7V8ysrKyGTydo8lkgkgoeHR6ubMzAajTAYDFi5ciXGjRtnOqdICLE+Dw8PjBkzBseOHTOdMjG34OBguLu7o7y83CLHZ5vdBBnDMHjuueewe/duHDlyBJGRka0eHzBgAAQCAZKTk033ZWdno6ioCImJidYu12YxDAOj0Yiff/4ZN27cwJIlS9guiRCnN27cOAiFQott9RIQEIDGxkZUV1db5Phss5uuxUWLFuGnn37Cr7/+Cnd3d9N5L6lUColEAqlUivnz52PJkiXw9vaGh4cHnn/+eSQmJlLH4k2MRiOam5vxySefYMaMGejZsyfbJRHi9CQSCSZOnIhff/0V48aNg5+fn9mPLxQKIZfLzXpcW2E3I7L169dDoVBg9OjRCAwMNN127Nhhes7atWtx3333YebMmRg5ciRkMhl++eUXFqu2PQaDAT///DNqa2vx8ssvs10OIeT/GzFiBKRSKfbu3Wv2Y3M4HHh5eTlskNnNiIxhmLs+RywW48svv8SXX35phYrsU0NDA9auXYu5c+ciIiKC7XIIIf+fQCDAlClTsHXrVowfP97sXdSenp4OG2R2MyIj5rF9+3bI5XLMnz+f7VLIHTAMA71eD41GA7VajebmZjQ1NaGxsbHVrbm5GWq1GhqNBlqtFnq93rSOJrE/gwcPRnR0tEU6GL28vG5phnMUdjMiI12n0+nwzTffYNq0aTQasyF6vd4UQgaDAQaDwbQWJvDXpow8Hs+03U7Ln1wu13RNYMtF7jfv9s3j8cDlck1//n0DVWJ7uFwuEhIS8NNPP6GyshIBAQFmO7aXl5fDXmJEQeZEdu/ejZKSEmzatIntUpya0WiERqMx3VrCSCQSQSAQQCgUmsKrJYTaq6UrtSXYjEYj9Hq96fGWQLs5GIltGTBgAPbu3Yvjx49j9uzZZjsunSMjdo9hGKxbtw6jR49Gnz592C7H6RiNRjQ1NaG5uRkajQYcDgdCodDUTSYSicwSLBwOp9Uecy1aQq0l2Fq28GkZwVGo2Q4+n48xY8bg4MGDuO++++Di4mKW43p5eaGpqQkajcbhNs2luQYncfjwYWRnZ+Ppp5+mDy0rUqvVqK2tRWlpKerq6sAwDLy8vBAYGAg/Pz94eHhALBZb/O+kJeAEAgHEYrFpRNZycbzBYGjzvFpSUhISExMRHh6OxMREJCUlWbRO8pdBgwZBr9fj3LlzZjtmy/J9jniejILMSXz++ecYMGAAXVNnBQzDoKGhAWVlZaioqIBarYZUKkVwcDD8/f3h6urK+vmqlilGPp9vGsH9PcyefPJJTJw4EWlpaSgqKkJaWhomTpyIf/3rX2yV7TTc3d3Rp08fpKammq1xx8vLCwAccnqRgswJnDt3DkajEQsXLmxz2omYj0qlQnFxMWpqasDj8RAQEIDg4GB4eHjY7HvfMlrj8/mmkWFSUhI2btzY5vO/++67VivoEMtITExERUUFCgoKzHI8CjJi1zZu3IjKykqMHj2a9ZGAo2psbERRURGqq6shEokQFBSEgIAASCQStkvrlLvtTffGG29YqRLnFRMTAx8fH5w6dcosx5NIJHBzc2u1GbGjoE81B1dbW4s9e/bgiSeeoBZsC9BqtSgrK0N5eTl4PB6Cg4MREBBg97sJlJWVdelx0nUcDgeJiYnIzMw0y6r1HA4HXC7XIVfAp081B/frr79CIBBgzpw5AEBBZiYMw6CmpgaFhYXQ6XQIDAxEcHCww3SDBQUFdelxYh6DBw+G0Wg0W9NHcHCwQzZ70aeaA2MYBt988w0mT55s6liiIOs6jUaDgoIC1NXVwdvbG2FhYXB1dWW7LLN6++237/j4Sy+9ZJ1CnJy5mz7kcjk0Go0ZKrMt9KnmwC5duoScnBzMnDkTRqPRIX8Tsza5XG46+R4eHg4fHx+HfF8nTJhw22XMXF1dsXTpUpw4ccLKVTmnlqaPwsJCtkuxWRRkDmznzp3w9fXF6NGjYTAYHPID11oMBgOKi4tRWVkJT09PREREOMw04u18++23OHz4MBISEhAWFoaEhAQcPnwY+fn5uOeee3Dvvfdi7dq1tK6jhcXExCA4ONgs04uO+ndFK3s4qJbtWv75z3+Cz+ejubmZphU7SaPRoKioCAaDASEhIXBzc2O7JKsZN24cxo0bd8v9+/fvx8qVK/H777+jpqYGq1atAp9PHyeWwOFwEBcXh3PnzmHGjBn077gN9I44qJMnT6KystK0VhtNLXZOU1MTbty4AQCIiopyqhC7Ez6fj/feew8LFy7Exo0bMX/+fOh0OrbLclj9+vVDQ0OD6WexKxzxc4B+hXJQ//vf/xAZGYn+/fsDgGlldNJ+CoUCJSUlcHFxQVhYmM1e0MymOXPmQCKR4PHHH4der8fGjRvt/tIDWxQSEgJvb29kZmYiOjqa7XJsDo3IHFBzczP27t2L2bNnm8KLgqxjamtrUVJSAqlUioiICAqxO7j//vvx008/Yf/+/XjsscccsiuObRwOB3379kVOTk6rLX46ylE/ByjIHNCBAwfQ0NCAWbNmme5z1B9gS2gJMU9PT4SEhND71g5TpkzBjh07cPjwYbz22ms0zWgB99xzD+rq6hx2T7GuoCBzQKdOncK0adMQGRkJ4P86legD+e7kcjmKi4vh4+OD4OBgtsuxKxMnTsT27duRmZmJN998k+1yHE5ERAQkEgmys7M7fQxH7VqkIHMwarUaW7ZsMZ0bAyjI2kupVKKoqAje3t4ICQlhuxy7NGHCBDz66KNYv349Nm/ezHY5DoXL5SI2NhY5OTlsl2JzKMgczIkTJ6BWqzFlyhTTfRRkd9fQ0ID8/Hx4eHggNDSU7XLs2rx58zB//nwsXboUaWlpbJfjUGJjY1FUVNSl9RId8XOAgszBHDp0CEOGDEFMTIzpvpaTw474A2wOWq0WeXl5kEgkiIiIoPfJDNasWYPBgwfjlVdeQUVFBdvlOIy4uDgwDIPc3NxOfT9NLRK7kJKSgtjY2Fb3tfzw0oWUt2IYBnl5eeBwOIiKiqIQMxOhUIhNmzahoaEBixYt6lKnHfk/np6eCAgI6NJ5MkdEn2wOpKqqCllZWRgxYkSr+xmGoa7F22iZpunWrRutTGFm/v7++PDDD3HkyBH88MMPbJfjMOLi4pCdnd2p0ZWrqytcXFwsUBW7KMgcyMmTJwEAw4cPb3V/S5CR1mpra1FdXe2Qq9fbinHjxuHll1/GqlWrUFRUxHY5DiE2NhY8Hg81NTUd/t6qqiqHvCaSgsyBpKSkIDo6GoGBgWyXYvM0Gg0KCwvh4+MDPz8/tstxaC+88ALc3Nzw8ssv0y9UZhAeHm7qsO0ItVoNhmHg7u5uocrYQ0HmQFJSUm6ZViRta1mzLiwsjOVKHJ+bmxv+85//QC6X49ChQ2yXY/fEYjH8/Pw6HGT19fUAYNqb0JFQkDmI8vJy5ObmYuTIkbc8RlOLrdXU1EChUCAqKsohp1ls0eTJk+Hq6op3330XBoOB7XLsXlhYWIeDTC6XAwC8vLwsURKrKMgcxO3Oj7WgIPuLXq9HQUEBfH19HfI3U1vF4XDw5ptvIisrCzt37mS7HLsXFhaGioqKDi0F1jIik0qlFqqKPRRkDiIlJQXdu3eHv7//LY/dvHCws2vZZTc8PJzlSpxPfHw8pk2bhjVr1tDCwl0UFhYGo9GIkpKSdn9PfX09JBKJQ24IS0HmIC5fvkznx+6ioaEB1dXVCAkJgUAgYLscp7R8+XLIZDLs2LGD7VLsWkBAAIRCYYemF+vr6x1yWhGgIHMINTU1OHv2LIYMGdLm4zQi+0tRUREEAgECAgLYLsVpxcTEoF+/fvj++++d/uexK7hcLrp164aqqqp2f099fb3DTqd3OMieeOIJnDhxwhK1kE66evUqEhISWi0UfDO6EBpQqVSoq6tDWFgYvR8smzZtGrKzs2kdxi7y9PREWVlZu58vl8spyFooFAqMHz8eMTExeP/991FaWmqJukgHXLlyBenp6aZtW27HmX8DLiwshIuLC3x9fdkuxekNHToUUVFR2LJlC9ul2LWAgABUV1e3e/kvGpHdZM+ePSgtLcUzzzyDHTt2ICIiApMnT8auXbtoMz2W5OTkoFu3brdtJXf2EYhKpYJSqaTRmI3gcDh49NFH8fvvv6Ouro7tcuyWv78/DAYDamtr7/pchmFoRPZ3fn5+WLJkCS5cuIDTp08jOjoajz32GIKCgrB48eJOr8xMOic7O/uWhYJvxuFwYDQanXbh1tLSUgiFQpsejRmNRuj1etPNYDA49Ah6zpw5AECt+F3Q0qHcnvNkTU1N0Ov11OzRlvLyciQlJSEpKQk8Hg9TpkzBpUuXcM8992Dt2rXmqpHcRW5u7l2DDHDOqUWdTofKykqbavDQaDRQKBSoqqpCaWkpiouLUVJSgoqKClRVVaGqqgrV1dWmm1wuR0NDA9RqNfR6Pdvlm4W3tzceeughnD59mu1S7JabmxtcXV3bNSKrq6sDh8Nx2BFZh5f71ul0+O2337Bx40YcOnQIffr0wUsvvYSHH34YHh4eAIDdu3fjySefxOLFi81eMGmtoaEBZWVldwyylu1bnHFEVlFRAQ6HA5lMxmodTU1NUKlUaGxshNFoBI/Hg1AohFgshkAgAJ/PB4/HA4fDMe1UwDAMDAaDaTTd1NQEHo8HLpdr+h57XrE/MTERzz//vM39omEvOBwOvL292xVkNTU14HA4DruuaIf/FQQGBsJoNOKhhx7CmTNn0K9fv1ueM2bMGIdNflvTMo3bnhGZMwZZeXk5fH19WblujGEY1NfXo66uDjqdDgKBAB4eHnBzc4NYLO7U+TqDwQCdTmeaftTr9XYbaKNGjQKHw8Hx48cxe/ZstsuxS1KpFEql8q7PKy8vh0QiMQ02HE2Hf/rXrl2LWbNmQSwW3/Y5np6eyM/P71JhpH1ycnIA3DnIgL9GZc42tahQKKBWqxEdHW3115bL5aiqqoJer4eHhwcCAwPNsg8Uj8czNfW0BJlOp4PRaASfz7erzVO9vLzQp08fHDt2jIKskzw8PNr1WVtWVobAwECHbXbq8E/9Y489dscQI9aVk5MDmUx2160ZuFyu0y3WWl1dDR6PZ9UT3Gq1Grm5uSgpKYFYLEZ0dDRCQ0Mtspkhj8eDSCSCWCwGl8s1NYq0/MKSlJSExMREhIeHIzExEUlJSWavoatGjx6NEydOON3Pprl4eHi0e0QWFBRkhYrYYT+/vpE25eTkICYmpl3PdbYPi8rKSvj4+Fjtt9DKykpkZ2dDr9cjMjIS4eHhVlnXjsPhmKYXWzpUn3zySUycOBFpaWkoKipCWloaJk6ciH/9618Wr6cjRo4cidDQUGRnZ7Ndil3y8PBAc3PzHS99YhgG5eXlDr1PIQWZncvJyUFcXNxdn8fj8ZwqyJqamtDc3GyVlnuDwYC8vDwUFxfDx8cHcXFxcHNzs/jr/h2XywWPx0NycjI2btzY5nO+++47JCcnW7my2+vduzeys7ORmZnJdil2qWUl+zuNympqauDm5oaQkBBrlWV1FGR2zGg0wt/fH927d7/rc50tyKqrq8Hlci0+rajT6ZCVlQWlUono6GiEhISwfp5q5cqVd3z8jTfesFIldyeRSDBs2DBaIaiTWpo37hRkhYWFkMvlDr3jAwWZHWtqasKhQ4fg4+Nz1+fyeDyn6lqsq6uDl5eXRTfO1Gq1uHLlCtRqNWJjY22mU/du6+91ZH0+a/D19cWZM2fYLsMutSfI8vPz4ePjc9fz6PaMgsyOKRQKAO3bupzH44FhGKcIM2ssx6PX63Ht2jUYjUb07NkTrq6uFnutjrrbSX1bO+nfrVs35OXlsV2GXRIKhfDy8kJjY+Ntn1NQUHDXdVjtHQWZHevIjq8tQeYoK0PcSWNjIzQajcV2wjUajbh27Rqam5vRo0cPm+viffvtt+/4+LvvvmulStonOjoaQUFBkMvlbJdil7hc7m2DzGg0orCwEBEREdYtysooyOxYS5C1Z+TB5/PBMIxTLOysVCrB4XAsdvHnjRs3oFQqcc8990AikVjkNbpiwoQJmD9/fpuPjR07FuPGjbNyRXcmk8lw7do1VFZWsl2KXRKLxVCr1W0+Vl5eDo1GQyMyYrs6MrXI4XDA4/GcYkSmVCohkUgsstpFTU0NSkpKEBkZ2amg1Ov1aG5uRkNDAxoaGtDY2Ijm5mZotVqzXrD+7bff4vDhw0hISEBYWBgSEhIwdepU+Pj42FxjRUtnaXV1NcuV2CeRSASNRtPmY/n5+eBwOAgLC7NyVdZlf+vaEJOOTC0CgEAgcJoRmSXOWel0OmRnZ8PX17fd55mMRiMaGhqgUqmgVqvBMAy4XK6pVZ7D4Zi+5nK5EAqFpltXg3jcuHGtRl/19fUYOHAgXn31VZvaC6wlyGpqaliuxD65ubndtlO2oKAAgYGBNjf9bW4UZHZMoVCYVnZoDz6ff9spCEfS1NRkkUWCr1+/DoZh7rocGPDXtWXV1dWmBV1dXFwglUohFoshFotNIQbAtDhwy/qJWq3WtIaiSCQyW+elp6cn1qxZg3nz5uHIkSMYO3asWY7bVWKxGPHx8U7xs2kJRqMRzc3NbT5WU1ODHj16WLki66OpRTsml8s71NAgEAhaLWHkiBiGgVarNfu5K5VKhYqKCoSFhUEoFN7xuTU1Nbhy5Qqqq6vh5eWFmJgYREZGwtfXF25ubqYVOFq0rMwhFovh6uoKT09PuLi4gMvlQqPRmKYdzbHk1MyZM/Hkk0/igw8+sKkO1rKyMhqRmZlKpcK1a9ccvtEDoCCzax3durxlqsqRz5O1fPCbO8jy8/PB5/MRGhp62+cYDAbk5uYiLy8P7u7u6N69O4KCgu4afG1pCTaRSAQOh2O2Jac4HA7mzJmD06dP4+DBgx2uy1Isfc2fM8rKygKAdi2YYO8oyOyYQqHoUJC1bGXiyEGmVqtNe32ZS2NjI8rKyhAZGXnbdRvVajUuXrwIhUKBuLg4REZGmmXrGB6Ph2PHjmHTpk1tPt6ZJaeGDh2Ke++9F7///nuX6zOXlvO9xHyuXbuGoKAgm7lQ35IoyOxYfX19h6YWWxoKtFqtBatil06ng06nM2uQFRUVgc/n37bBQ61W49KlSzAajejVq5fZPzhWrFhxx8c7s+TUo48+iitXrpi2AWKbTqezyz3VbEFbpwoYhsG1a9ec4vwYQEFm19zd3Ts8/83n8x06yFpGm+b6UDQajaioqEBwcHCbnWF6vR5Xr14FwzDo3bu3RbrDLLHk1Pjx41FSUoKff/65s2WZlV6vZ2XzU0dVVVWFuro6CjJi+1r2o+oIoVDo8EHW0jxhDrW1tVCr1QgODm7z8aysLKhUKvTq1cuso8CbWWLJKZFIhMcee8xm1jikEZl5Xbt2DVwut10dto6AgsxOJSUlYceOHdiwYUOHOtiEQiF0Op3DniczGAwwGAxmW4G+oqICfD6/zYufKyoqTO3Nltg4s4WllpwaNGgQtFot64sIGwwG9OrVq12LX5P2uXbtGrp162aV/fBsAQWZHWrpYFMoFGhqaupQB1vL1JcjX7PD4XDMdolBXV0d/P39b7lfr9fj+vXr8PLyavNxc7rTklP9+/fv9JJTgwcPRkZGBjIyMrpSXpdVVVXhwoULDr06uzXp9XpUV1ejd+/ebJdiNRRkdiYpKalLmya2rB7hqEHWsu+aOa6R0ul0qK2tbbN5o7i4GA0NDe3enbur2lpyauHChfD09ERxcXGnjunn54eYmBhcuHDBzNV2TMuI0NZW5bdXubm5KC8vR8+ePdkuxWooyOyMOTrY7rTIqL1ruRbJHJuItuzx9PfNORmGQUlJCQIDAy06pfh348aNQ2pqKgoLC5GamooPP/wQAG7bmt8effv2veNeVtZAQdY1fD6/1c9hZmYmfHx8bnte1xFRkNkZc3SwicVi6PV6hzxP1hJk5vh/UyqV4PF4cHNza3V/bW0t5HI56wuxuru7Y9y4cbh+/XqnR6ABAQGsLyJcWloKb29vm9xJwB40NjaaptIZhkFmZib69+9/22seHREFmZ0xRwebI58nEwgEZuvMbG5uNi3qe7Py8nK4urreMlJjw9ixY3HgwAGcP3++U98vlUpx9epVM1fVMeXl5QgMDGS1BnvG5XJN5xfz8/OhVCrRt29flquyLocMsi+//BIREREQi8UYMmSIzbQYm4M5Otgc+TyZWCxGU1PTbRdR7QitVttm11dVVRW8vb27fHyj0Yht27Z16RiDBg3CuHHjcOnSpU59v1QqRUNDA6vrb+bl5TlNm7glVFdXm345vXTpEqRSKaKioliuyrocLsh27NiBJUuWYOXKlcjIyEDfvn0xadIkVFVVsV2aWdypg23+/Pnt7mBz1PNkEokEHA4HTU1NXT5Wc3PzLaMxg8EApVJplk07165di4cffhgff/xxp4/RcoH7uXPnOvX9Li4u6N69O2sLCOv1ety4ccOpOuzMyWAwoKmpCR4eHjAajUhLS8PAgQPNdvmJvXC4/9tPPvkECxYswLx583DPPfdgw4YNcHFxwffff892aWbTVgfb4cOH8e2337b7GI56nozD4cDFxeW2W7935ng3a2hogNFoNEuQfffddwDQ5Z/NQYMGdXqfuebmZuTk5LC2YG9WVhYaGxudZgUKc1OpVAD+Ol969epVKJVKDBkyhOWqrM+hLqXXarVIT0/H8uXLTfdxuVyMHz8eqamptzxfo9G02lmV7e6tjvj7pokd1TJlplarb2lmsHeurq5oaGjo8nFu3jOshVqthl6v7/RSVE8++STKy8sB/N/q5NeuXcPkyZMBAIGBgR0ONldXV1RUVHSqnsbGRotsQtpe58+fh0gkoiDrpJbPLA8PDyQnJyM4OPiOOzQ4KocKspqaGhgMBgQEBLS6PyAgwPShcbPVq1dj1apV1irPpvB4PAgEAjQ2NjpckHl4eJhlbysOh3NLIOr1enC53E4tp9TQ0IBNmzbdcj6KYRgcOHDA9Jqff/55h/5OXF1dOz2y1uv1GDp0aKe+1xwyMjLQu3dvWmexk1qCjMvl4tKlS/jnP//JckXscLipxY5Yvnw5FAqF6dbZC0vtlUQiadW66yg8PDxQW1vb5WlTiURyy5Qdl8uFWq3u1DklNzc3nDlz5rY7FkilUpw7d67Dv1gwDNPpX0auXbtmlvOJnWE0GlFdXY2RI0ey8vqOQKVSQSQS4eLFi2AYBoMGDWK7JFY4VJD5+vqCx+OhsrKy1f2VlZWQyWS3PF8kEsHDw6PVzd505SS9m5vbHbdJt1c+Pj5gGAZ1dXVdOo6Liwu0Wm2rVn6BQACRSNTp92zgwIG3nQasqKhAfHx8h49ZVVXV6XOC+fn5iIyM7NT3dtXFixdx5coV9O/fn5XXdwRKpRLu7u5ITU1F7969HW52pb0cKsiEQiEGDBjQapkmo9GI5ORkJCYmsliZZXV2RCUUCiEQCMxyPsmWeHp6gs/no7a2tsvHUavVkMvlpvtazr915T3bsmVLm/dv3bq1U8fT6XSdCoPm5mZ4eHigX79+nXrdrkpNTcWgQYPQq1cvVl7fEahUKuj1epSUlDj0Z9zdOFSQAcCSJUvwzTffYPPmzbh27RqeeeYZNDY2Yt68eWyXZpPc3NwcbnqRw+EgICCgy5dceHp6QiKRtDqORCKBl5dXl87B/fDDDwCA0NBQnDp1CiEhIQCAzZs3d+p4J0+e7NRv4ufPn8elS5dwzz33dOp1u8JoNOL3339HVFSU07WKm5NSqURJSQnc3d2dam3Fv3OoZg8AmDNnDqqrq7FixQpUVFSgX79+OHDgwC0NIOQvbm5ukMvlaG5utuq6gZbm4+ODq1evwmg0dvqDksvlwtvbGxUVFa2uc3J1dUVpaSkGDhzYqeMuWrQICQkJWLNmDbhcLgoLC7Fs2TIMGDCgw8eqqKiAv78/EhISOvy9aWlp8PX1ZWVq8cKFC6isrDR1a5LOkcvlKCoqwpw5c5z6FwKHCzIAeO655/Dcc8+xXYZVdHU9NZFIZJpedKQgCwoKQnp6Oqqqqto8P9peMpkMV65cgV6vN3UqhoeHIy0tDc3NzZ1aH3DOnDmYM2eO6Wsul4sPPvigU/UdO3YMV69e7fC5NYZhkJubi+nTp7OyJt+JEyfQr18/uhC6CxiGQV5eHgQCQad+kXEkzhvhDqSr04KOOL3o7e0NNzc3FBUVdek44eHhUCgUrY4THh4OrVaL7OzsrpbZJQzD4ODBg5g1a1aH9/K6dOkSTp06xcp5lYaGBuzatQuJiYlOPYroKpVKhbKyMoSFhTn9WpX0U0Tg5uYGg8HgcEtWRUREoLi4uEsB7enpCZlMhtzcXNN9YrEYMTExOH/+vFm2i+ms9PR0yOVyjB49usPfm5ycjO7du2Pw4MHmL+wu9u3bB51Oh3/84x9Wf21HcvHiRdTU1GDChAlsl8I6CjJiml5sWe7GUYSHh5t+a+2KuLg45ObmQqFQmO7r06cPjEYja5tSGgwGfPnllwgMDOxwGNXW1iIjIwNTp061+ojIaDTi2LFjmD59Ovz8/Kz62o7m6NGjcHd3x4gRI9guhXUUZHaOw+GYZUrQzc0NKpXKoaYXfX19IZVKkZeX16XjxMbGQiqVIjMz03Sft7c3IiMjkZqaysrSZnv27IFQKMSCBQs6HEabN29GUVERK6tAHDlyBBkZGbj//vut/tqOpKamxrRAMFvrZNoSCjIHwDBMl1cvd3d3h8FgcLhrymJjY3Hjxo0uTZvy+Xz06tULFy5caHVN2bBhw+Dp6Yl9+/ZZdfX469evY/v27YiOju5wk0dZWRmOHDmCBQsWdPi8WlcZjUYcPHgQQ4cOZaXl35EkJSVBp9Nh1KhRbJdiEyjIHIA5us6EQiEkEkmr6TNHEB0dDZFI1OXNI1tWTThx4oTpPpFIhFGjRkGhUODw4cNdLbVdamtr8frrr0Mmk3W4M5dhGHzyyScIDw/HjBkzLFTh7SUlJeHkyZN0TWcX1dXVISUlBWFhYQgPD2e7HJtAQeYgzDEl6OHhgaamJrPsrmwrRCIRIiIicOXKlU5vdQL8NSobOXIksrOzce3aNdP9ISEhGDZsGM6dO4cjR45YdGq2rq4O//73vxEWFoY33nijw63/SUlJyMnJwezZs61+qYVOp8NXX32FkSNHok+fPlZ9bUdz6NAhGAwGhIWFtWtHeGdAQeYAzHWezMPDAxwOB/X19V0vyob06dMHOp2uzR0QOiI6OhqxsbE4dOhQq8aY3r17Y/LkyTh9+jT27NnTamsgc7l+/TreeOMNcLlcvPDCCx2+wL+kpATfffcdBgwYwMp01O7du1FZWYlnn33W6q/tSOrq6nD69GlER0fDx8eH1S14bAkFmYMwGo1dDjMOhwNPT0/U19eztmOwJbi6uiI6Otp0YXNXTJo0CVKpFL/88kur1vv+/ftj+vTpKC8vx/r161FQUNDFqv9iMBiwe/durF27FlwuF6tXr+7wflONjY14//33IRQK8e9//9ssdXVEbW0tvv32Wzz66KM0FdZFhw8fhkQigY+PD43GbkJB5gBautbMMSrz9PSE0Wh0uHNl/fr1g1Kp7HK7vEQiwYQJE1BWVmbaQ6xFjx49MGfOHPj5+eGHH37Azp07b9mJob2MRiNOnTqFZcuWYfv27RgyZAg+/vhj+Pv7d+g4Op0Ob731FsrLy/H222+z8hv8N998g5CQEDz88MNWf21HIpfLkZaWhtGjR6O2tpaC7CYOuUSVMzLX9KJAIIC7uzvkcjm8vLzMUJlt8PDwQJ8+fXD+/HnExcV1abuL4OBgTJ8+HX/88Qd4PB4mTZpkarjx8fHBQw89hIyMDKSnp+Pjjz9Gr169EBMTg+7du9/xPTUajSgoKMCFCxeQnJwMoVAIT09PrFmzplPrIer1evznP/9BeXk53njjDVbWVExOTsa+ffvw1ltv3XYfNtI+hw8fNu2mfe7cOQqym1CQOQgOhwOj0WiWa0q8vLxQVFSExsZGh5qDj4+PR3Z2Nk6fPo1x48Z16Vg9evRAU1MTkpKSoNfrMWXKFNPImMvlYuDAgejfvz8uX76M8+fP49ChQ/j111/h7e0NsViMgIAAGI1G03YzSqUSpaWlqK+vh0wmQ69evTB+/HhER0d3qr7m5masXLkSeXl5ePnllzu1IHFX1dbW4pNPPsHYsWO7/H47u/r6eqSmpmLy5Mmora0Fh8Pp0hqijoaCzEFwudwun/9p4eLiApFIhJqaGocKMqFQiCFDhiA1NRWlpaUIDg7u0vEGDBgADoeDI0eOoK6uDrNmzWrVScjj8dC3b1/07dsXjY2NuH79OiorK1FYWAiVSoWKigp4e3ujqakJISEhppFbVFRUl34hKS8vx4cffgitVos33niDlY0rGYbBxx9/DD6fjyVLllj99R1Ny2hsxIgROHLkCPz8/CAQCNguy2ZQkDmIlqmtrmxbcjMfHx+UlJR0eoV3W9Wy3FRKSgpmzJgBoVDYpePFx8fD3d0dBw4cwH//+1/Mnj27zWYMV1dX9O3bt0uv1R4pKSn4+OOPERYWhqVLl7K2+/Mvv/yCK1euYPny5TSl2EVyuRxZWVkYPXo0xGIxamtrTXvYkb9Qs4eDuDnIzMHDwwNCobBLG0jaIg6HY7qIOSUlxSzHjImJwSOPPAIvLy9s3boV+/fvt/q1eNXV1Xj//ffx6aefYty4cVi9ejVrIXbu3Dls2LABM2fOdPrtRczhwIED0Gg0GDVqFBobG1FfX4+wsDC2y7IpNCJzMOZsm/fz80NpaSnUajXEYrHZjss2Dw8PjBw5EidOnEBISAji4uK6fExvb2888cQTSElJQVJSEvLy8jBo0CAMGjTIomvhqVQq7N+/HydPnkRtbS2WLFmCoUOHsrLHGAAUFRXhrbfeQnx8PB555BFWanAkN27cQHp6OubMmQOxWIzMzEyo1WrWfkmxVRRkDoTH43Vp9Yq/8/DwQHV1NWpqahxuKqNHjx4oLi5GcnIyfH194ePj0+Vj8ng8jB49Gn369MGRI0ewa9cuJCUlYejQoRgyZEiXOiX/rqX2/fv3w8XFBWPHjsXMmTNZPaepUqnw2muvwdfXFytWrKDFbLvIYDBg165dCAsLM+1GnpOTg9DQUIea7jcHCjIH0nJuzFznyTgcDnx9fVFWVgaNRgORSNTlY9qSMWPGoLq6GklJSZg+fbrZRp3e3t544IEHMHz4cPz555/Yt28f9u3bh/j4eHTr1u2ubfhtYRgGhYWFuHjxIjIzM5GRkYGwsDDMmjULEydOZP08lFarxaefforGxkb897//dagmIbYcO3YM1dXVWLJkCTgcDtRqNYqKijB27Fi2S7M5HMaR9u3oIqVSCalUCoVCAQ8PD7bL6RSNRgMul2u2jqaW7dSFQqFDzssrFAps27YNMpkM9913H/h88/9u19jYiHPnziEnJwcXLlyAVCqFh4cHPDw8EBISAg8PD7i7u0MoFILL5cJgMKCxsRFyuRyVlZVQKpUoLi5GQUEBIiIiEBUVhb59+2LQoEE20bmm0WiwcuVKqFQqLFq0iFa2N4Pa2lp88MEHGDFiBO677z4AwJUrV7Bv3z48++yzZh3d26qOfB5TkN3EEYJMp9PBaDSadfQkl8tRWlqKbt26OeSURllZmWkK5/7777foZpMtbfj5+fkoKCiAXC5HbW0txGIx6uvrodPpEBAQgKamJjQ0NCAkJAQhISGQyWSIiYlBbGysTYRXC41GgzfffBNXr17Fu+++i379+rFdkt1jGAbffPMNqqqq8Morr5g6a3/55Rc0Nzc7zbnHjnwe09Sig+HxeNDr9WAYxmwn/D09PVFTU4OKigqHPMkcFBSE+++/H3v27MEff/yBKVOmWKxZoqUN/+ZWfKPRiKamJqjVajAMAy6XC4lEAolEwlrTRns0NzfjzTffRHZ2Nt5//31a1d5MMjMzkZ2djfnz55tCTKvVIj8/HyNHjmS5OttE7fcOhsvlgmEYs10cDcC0ioBKpWJlN2RriIiIwNSpU3H9+nXs37/fqosmc7lcuLm5wdfXF35+fvDx8YGLi4tNh5hKpcJHH32EnJwcrF69mkLMTJqbm/Hrr7+iT58+raZob9y4AYPBgNjYWBars10UZA6Ix+O1WpndHNzd3eHu7o6ysjKHWhn/ZjExMbj33ntx6dIl7Ny506wdoI6koKAAL774IuRyOT788EP06tWL7ZIcRss1iNOnT291f05ODgICAlhv6rFVFGQOiM/nm2Vbl78LCgqCTqdDdXW1WY9rS+Li4jBnzhwUFRXhp59+QnNzM9sl2ZRTp05h8eLFEIlEeOWVV8xyDR75S2FhIVJTUzFlypRWgaXX65GXl0ejsTugIHNALc0K5h6ViUQi+Pn5oaqqyqF2kf67yMhIPProo2hqasLGjRtRVVXFdkmsMxqN2LJlC9555x0MHDgQn3zyCS1aa0YGgwE7d+5ESEgIhg4d2uqxgoIC6HQ6CrI7oCBzQBwOx+wXR7fw9/cHn89HWVmZ2Y9tSwIDA/Hggw/CxcUF33//PS5evMh2SaypqqrC+++/j1OnTuHxxx/Ha6+95pDdq2w6ceIEKioq8MADD9zSNZudnQ0fHx+zXLTvqCjIHFTLeTJzn8/icrkIDAyEUqmESqUy67FtjZeXFx555BF0794de/bswb59+6DRaNguy2oYhsH+/fuxaNEi5OXl4emnn8ZDDz1k000o9qiiogJnz57FuHHjbllBR6fTITc3l85D3gW13zsoPp8PDocDnU5n9hU5PD09UVtbi7KyMsTGxjr0B5tAIMD06dMRFRWFvXv3IicnB9OmTev0PmH2ory8HF988QUuXryIe++9F/PmzYOLiwvbZTkcrVaLLVu2QCAQtLlnW1ZWFvR6PXr06MFCdfaDgsyB8fl86PV6iywtFRwcjNzcXFRVVSEgIMDsx7c1ffr0QVhYGH799Vfs3r0bsbGxGDNmjN1eOH87TU1N2L17N06fPg2NRoN3333XKtvPOCOGYfDzzz9DoVDghRdeaHNLoYsXL6J79+4O93NmbhRkDkwgEECn00Gv15t96SWxWAwfHx9UVFRAKpU61Or4t+Pp6YnHHnsMmZmZOHnyJD799FMMGzYMw4cPt/t1KHU6Hfbv349du3ZBq9XigQcewLRp05zi75Utp0+fRmZmJh555BH4+fnd8nhpaSkqKysxbNgwFqqzLxRkDozL5Zp2jrbEGoIymQz19fUoLCx0+CnGFlwuF/Hx8bjnnnuQkpKCkydP4uzZsxgxYgQGDRrU5Y06rU2j0eD48ePYs2cPKisrMWHCBMyePRve3t5sl+bQSkpK8Ntvv2Ho0KG3vZg8PT0d3t7eiIiIsG5xdoiCzMEJBAJoNBqzLlnVgsvlIjw8HDk5OaioqEBgYKBZj2/LxGIxJkyYgMGDB+PEiRPYt28fDh8+jISEBCQmJtr8VFBlZSUOHTqE5ORkeHt7o3///pg6dSqCgoLYLs3hNTc3Y8uWLQgMDMTUqVPbfE59fT2uX7+OCRMmWLk6+0RB5uD4fD7UajW0Wq1Fpr9cXV0hk8lQUVEBDw8Pp9u+QyqVYtq0aRg5ciT+/PNPpKam4vTp04iKisKAAQMQGxtrM/tyabVaZGRk4MyZMzh58iRcXFwwbtw4TJgwga4JsxKGYbBjxw6o1WosXLjwtjMl58+fh0QioSaPdqIgc3AcDgd8Pt8i3YstZDIZFAoFCgoK0KNHD4uuHm+rpFIppkyZgrFjx5r2C/v+++9NiwTfc8896Natm0WmeO+ksbER58+fx9mzZ3HhwgV4eHggODgYTz/9NIYNG2b35/bszYkTJ3Dt2jXMmzfvtnvSaTQaXL58GQMGDLD6z4u9onfJCQiFQjQ2NkKn01lkCxAOh4PIyEhcu3YNJSUlDrlvWXuJxWIkJCRgyJAhKCsrQ2ZmJgoKCnDy5Enw+XzExMQgKioKYWFhCA4ONvuFxXK5HHl5ecjLy8Ply5ehUChQXV2Nbt264Z///CcGDx7sVFPAtiQ/Px8HDhzAmDFj0L1799s+79KlSzAajdQt2gEUZE6Ax+OBz+dDo9FYbC8rkUiEkJAQFBUVQSqVOv3iphwOB8HBwQgODgbDMKisrERWVhaysrJw+vRp/PrrrwAAHx8fREVFwcXFBV5eXpBKpfD09DRttMnj8Uw3g8EAtVqNxsZGNDY2QqlUora2FgqFAvn5+SgpKUFDQwMMBgO8vb3Ro0cPjB07Fn379qXmDZapVCps3boVERERmDhx4m2fZzQacf78ecTFxTndNH1XUJA5CYFAgMbGRot1MAKAr68vFAoFCgsL0b17d7vr4LOUlm1wZDIZRo8eDaPRiKqqKpSUlKC4uBj19fW4cOECFAqFaVmxgIAAlJaWmo7h6uoKlUoFFxcXNDQ0APhruTC5XI5u3brBzc0N48aNQ1hYGKKiom47bUWsz2g0Ytu2bWAYBg8//PAdp95zc3OhUqkwYMAAK1Zo/yjInIRAIACPx4NarbboNunh4eG4evUq8vLy0L17d6doye8oLpdrCraBAwea7mcYBk1NTVAoFGhsbIRarYbBYIDBYDDtZCASieDq6gpXV1e4ubnZ/L5lzq5lma+6ujo88sgjcHd3v+PzL168iNDQ0DavKyO3R0HmRMRiscVHZXw+H926dUNWVhaKi4ud+nxZR3E4HFNIEcdw9OhRpKam4h//+AeioqLu+Nz8/HyUlpZixowZVqrOcThfe5kTu3lUZkmurq4ICwtDVVUVamtrLfpahNiqtLQ0JCcnm643vJvU1FQEBQXRL3+dQEHmZMRiMfR6PfR6vUVfx8/PD76+vigsLERjY6NFX4sQW3Px4kXs3bsXw4YNw6hRo+76/Ly8PFRVVSExMdEK1TkeCjInY61RGQCEhYVBJBIhNzfXoTfiJORm2dnZ2LlzJ/r374/Jkyff9RymwWDAn3/+iaioKISGhlqpSsdCQeaErDUq43K5pl1tc3Nzzb43GiG2prCwENu2bUNcXBz++c9/tqsRp6VjlRYH7jwKMidkzVGZQCBAbGws1Go1cnNzTd13hDiaiooK/PDDDwgJCcGDDz7YrhVumpubcebMGfTq1Yt2gO4CCjInJRKJoNVqTdctWZKLiwu6deuG+vp63Lhxw+KvR4i11dXVYePGjfD29sZjjz3W7q7g1NRUAEBCQoIly3N4FGROqmXViJaLay3N09MT3bp1Q3V1NQoLC63ymoRYg1KpxPfffw+xWIy5c+e2e/3KmpoaXLlyBUOGDDH7UmXOhoLMibm6usJgMKC5udkqr+fr64uIiAiUl5e3WrWCEHvV1NSETZs2wWg0Yt68eR26BvDEiROQSqW33Y+MtB8FmRPj8/kQi8VoamqyWiOGTCZDSEgIiouLUVlZaZXXJMQSGhoaTEtPzZs3D56enu3+3ry8PJSUlGDkyJE2s82PPaOVPZycq6srtFotGhoarLYZZEhICPR6PfLz88Hn8+kkN7E7tbW1+PHHH6HX6/HYY491aEmplnb7sLAw2v3ZTGhE5uRalkXSarVWvdYrIiICvr6+uH79Ourq6qz2uoR0VWlpKb799lvweDz861//QkBAQIe+/9y5c1Cr1Rg5cqSFKnQ+FGQEIpHItDq+Ndvju3XrBk9PT2RlZdFSVsQu5ObmYuPGjfDx8cH8+fM7NJ0I/DWSS09PR79+/WhrHTOiICMAADc3NxiNRqs1fgB/jQZjY2Ph7e2Na9eu0TkzYtPOnz+Pn376CVFRUXjiiSfg4uLSoe83Go1ITk6Gp6cnbdNiZnSOjAD4a/NNiUSCpqYmiEQiq52A5nA4iIuLA5/PR25uLgwGA4KCgqzy2oS0B8MwSElJQXJyMgYOHIipU6e262Lnv7tw4QJqamowY8YMavAwMwoyYiKRSKDRaNDQ0GDVHZ45HA6io6PB5/Nx48YN6PV6WgGc2ASj0Yg//vgDZ86cwZgxYzBq1KhO7f9WXV2N06dPY8CAAZDJZBao1LlRkBETDocDNzc31NfXo6mpqcNTJ10VEREBPp+PgoIC6PX6u+7fRIgl6fV6/Pzzz7h27Rruv//+Tk8H6nQ6JCUlwcfHB4MGDTJzlQSgICN/IxAIwOfzIZfLIRAIIBAIrPr6ISEh4PP5uH79OnQ6HWJjY2kHZGJ1TU1N2LlzJ4qKivDQQw8hLi6u08f6888/0dDQgNmzZ3dqSpLcHQUZuYWnpye0Wi3q6urg7+9v9SCRyWTg8XjIycmBRqPBPffcY7EdrQn5u5KSEuzbtw8ajQZPPPFEl6a58/LycPXqVYwZM6bDHY6k/ejXA3ILDocDb29v6PV61NfXs1KDn58fevXqhcbGRpw/f96q3ZTEOTEMg9TUVGzatAl8Pr/LIdbQ0IBjx46hW7du6NGjhxkrJX9HQUbaJBAI4OnpicbGRjQ1NbFSg1QqRb9+/cDhcJCRkUEXThOLaWpqwvbt25GUlISEhATMmzevSw1PDMPg8OHDEAgEGD16tPkKJW2i+RpyW66urtBoNJDL5RAKhaxM70kkEvTr1w9ZWVm4ePEiIiMjERYWRufNiNkUFhbil19+gcFgwMMPP4zo6OguH/PcuXOoqanBlClT2r0aPuk8CjJyRzefL/Pz82MlQPh8Pnr27ImCggLk5+dDpVKhe/fudN6MdEnL9WHHjx9HWFgYZsyYAXd39y4ft6CgAOnp6Rg8eDBdE2kl9ElA7ojL5cLHxwdVVVVQKBSsnbDmcDiIjIyEh4cHrl27hnPnzqFHjx5Wvd6NOI6Ghgbs3r0b+fn5GDVqFEaOHGmWX9LkcjmSk5MRGRmJ+Ph4M1RK2oPOkZG7EggEkEqlaGhoYL3pwsfHBwMHDoRAIEBGRgZu3Lhh1fUhif3Ly8vDV199herqajz++OOdvsj577RaLQ4cOAB3d3eMHTvWDJWS9rKLICsoKMD8+fMRGRkJiUSCbt26YeXKlbes1n7x4kWMGDECYrEYoaGh+OCDD1iq2PG4ublBIpGgtrYWer2e1VrEYjHi4+MRGRmJwsJCnDt3jrWGFGI/tFot/vjjD/z2228IDAzEwoULzbaNitFoxPHjx6FWqzFp0iSrX3/p7OxiajErKwtGoxFfffUVoqOjcfnyZSxYsACNjY346KOPAPy13fjEiRMxfvx4bNiwAZcuXcKTTz4JT09PPPXUUyz/HzgGLy8vaLVaVFVVQSaTsXpxJ4fDQUREBLy9vXH16lWcOXMGMTExCA4OZq0mYruys7Nx+PBhaDQajBgxAgMHDjTr+d6TJ0+isLAQkydPpuluFnAYO52X+fDDD7F+/XrcuHEDALB+/Xq8/vrrqKiogFAoBAAsW7YMe/bsQVZWVruOqVQqIZVKoVAorLbJpL3R6/UoLy+HQCBAQECATXQPGgwGXL9+HaWlpfDx8UGPHj1MPwPEucnlchw8eBA3btxAdHQ0Jk6caPbzvBkZGTh37hxGjx6N2NhYsx7bmXXk89guphbbolAoWu3nk5qaipEjR7b6AJs0aRKys7Mhl8vbPIZGo4FSqWx1I3fG5/Ph7+8PjUaD6upqtssB8NfK/XFxcejbty9UKhX+/PNPFBUV0bkzJ6bX63HixAl8/fXXqK2txQMPPIDZs2ebPcSys7Nx7tw5DBw4kEKMRXYZZNevX8e6deuwcOFC030VFRW37NTa8nVFRUWbx1m9ejWkUqnpFhoaarmiHYhIJIKfnx+amppsakNMHx8fDBkyBD4+Prh8+TJOnTrF2sokhD15eXn4+uuvcerUKQwZMgRPPfWURUKmuLgYJ06cQI8ePahDkWWsBtmyZcvA4XDuePv7tGBpaSnuvfdezJo1CwsWLOjS6y9fvhwKhcJ0Ky4u7tLxnImLiwt8fX2hUqmgUCjYLsdEIBCgT58+SExMBMMwOHXqFC5dunRLYxBxPAqFArt27cKOHTvg5eWFBQsWYPTo0RZpvKipqUFSUhLCwsIwbNgwsx+fdAyrzR5Lly7F3Llz7/icm7fyKCsrw5gxYzB06FB8/fXXrZ4nk8lu2WG45evb7f8jEonoqvsucHNzg16vh1wuB4/Hg5ubG9slmXh5eWHYsGEoLi5GdnY2KioqEBsbS6uCOCC1Wo20tDRcu3YNOp0O06dPxz333GOx16uvr0dycjJ8fX0xbtw4WtHeBrAaZH5+fvDz82vXc0tLSzFmzBgMGDAAGzduvOWHJzExEa+//jp0Op3pN7CkpCTExcXBy8vL7LWTv3h6esJgMKCmpgZcLtfqe5jdCYfDQVhYGGQyGbKzs3HlyhWUlJQgLi4Ovr6+bJdHukiv1yMjIwOnTp2CwWDAoEGDkJCQYNFGn/r6euzbtw8SiQSTJk2i1WVshF10LZaWlmL06NEIDw/H5s2bW20T3jLaUigUiIuLw8SJE/Hqq6/i8uXLePLJJ7F27dp2t99T12LnVVVVobm5GQEBARCLxWyX06b6+npkZWWhpqYGvr6+6N69O22tYYe0Wi0yMzNx5swZ8Hg8REZGYvjw4RafEVAoFNi3bx9EIhGmTp1qsz/njqIjn8d2EWSbNm3CvHnz2nzs5vIvXryIRYsW4ezZs/D19cXzzz+PV199td2vQ0HWeQzDoKKiAmq1GoGBgTb9j7yiogJZWVloaGhAYGAgoqOj6dofO6DRaJCeno5z585Bq9Wid+/eGDJkiFV+GVEqldi3bx8EAgGmTp0KiURi8dd0dg4XZNZCQdY1LWGm0Wggk8lsOswYhkFpaSlyc3OhUqkQEBCA2NhYmoa2Qc3NzTh37hzS09NhNBrRp08fDBkyxCwL/LaHQqFAcnIyDAYDpk6dalPT546MgqyTKMi6riXMmpqaEBgYaPP/6P8eaH5+foiIiIBMJqOmEJbV1NTg/PnzKCgoQENDA/r3749BgwbB1dXVqjUcPHgQEokEkydPppGYFVGQdRIFmXncHGYBAQE21c14Oy01X79+HbW1tZBIJIiMjER4eDh1tlqR0WhEXl4eMjIyUFxcDFdXV/Tv3x99+/a1eoiUl5fj8OHD8PLywoQJE+jnwMooyDqJgsx8GIZBVVWVaZRjT+eg6uvrkZ+fj5KSEjAMg+DgYERGRrZaSYaYV319Pa5evYrc3FxUV1cjODgY/fv3R0xMTKvmLmspKCjAsWPHEBgYiHHjxlF3IgsoyDqJgsz8ampqUF9fD29vb7sLAq1Wi6KiIuTn56OxsRFSqRQREREIDg6mtRzNQKvVIicnB1euXEFpaSmEQiF69uyJnj17wt/fn7W6cnJy8OeffyIiIgKjRo2i68RYQkHWSRRkllFXV4e6ujp4enra5fVbLaPLGzduoKqqCnq9Hv7+/ggJCaFQ6yCNRoP8/Hxcv34dlZWVUCgUCA8PR8+ePREdHc3qyIdhGJw/fx7Z2dmIiIhAQkICnSdlEQVZJ1GQWY5CoUBNTQ0kEgnrW8B0hVqtRmlpKUpKSlBTUwMAFGp30dTUhLy8PFy/fh3FxcUwGo0ICAhAXFwcYmJirNZ9eCc6nQ4pKSkoKirCwIED0atXL7ZLcnoUZJ1EQWZZDQ0NqKiogEAgQHBwsN2fd7g51Kqrq8HhcBAQEICAgADIZDKb+IBmA8MwqKmpQWFhIW7cuIGysjJwOByEhISgW7du6Natm029N42NjUhOToZKpcLIkSNp8XAbQUHWSRRklqfRaFBaWgoACAoKsulrzTqiJdQqKipQXl4Oo9EIFxcXyGQy+Pr6ws/Pz6pt49ak0+lQXV2NiooKVFRUoLa2FnK5HHw+H5GRkQgLC0NUVJRNtq5XVVXhyJEj4PP5GDduHF1HaEMoyDqJgsw6DAYDSktLoVarERAQYFcdje2h1+tNH+w1NTWmrW4kEgn8/Pzg4+MDT09PeHp62mVLd319vSm0KisrUVNTA4ZhTJutBgUFISgoCIGBgax0HLZXVlYWLly4AA8PD4wZM8ZhfqlyFBRknURBZj0Mw6CyshL19fWQSqUOfQGyVqtFdXU1ampqUF1dDblcDp1OB+Cv7XBaQs3T09O0Nx7b5xANBsMtm86qVCoUFxdDo9EAALy9vSGTyUxTqd7e3nbxd6jVapGamorCwkL06NEDAwYMYP39JreiIOskCjLrUygUqKiogFAodJpmCYZhoFKpIJfLUV9fb7o1NjaCy+XCYDBAKBTC1dUVEomk1c3FxcX032KxuMMfwAaDAXq9HjqdDmq1Gs3NzVCr1aagUigUUCqVaGxsNH0Pj8eDh4cHfH194eHhAZlMBn9/f7scTdbW1uLEiRPQaDQYOnQowsLC2C6J3AYFWSdRkLGj5byZXq+HTCZz2vdeq9Wivr4eSqUSDQ0NaG5uRnNzM5qamkyBc/M/V4PBAB6PB4lEAp1OZ9qMtiXcOBwO+Hw+mpubodPpoNfrTd8vFovR1NQEAODz+eByufDw8LjlJpVK4eLiYhcjrbvJysrCuXPn4OXlhZEjR9pUwwm5VUc+j+27bYw4BJFIhIiICJSXl6OsrAwqlQoymcymz69YglAohL+//20vBmYYBmq12hRsGo3GFFAGgwEMw8BoNJqeazQawePxTIF2800oFEIkEkEsFkMsFtt9B+mdNDc3Iy0tDdXV1YiNjcWAAQOc7mfL0TnuTy+xK1wuF8HBwaZGgqamJgQFBTlsp19ncDgc07QiaZ/8/HycPXsWPB4Pw4YNQ3BwMNslEQugICM2xdPTE66urigrK0NhYSG8vb3h7+9PJ+NJh6jVapw+fRrFxcWIiIjAoEGD7PKcHmkfCjJicwQCAcLDw1FXV4fKyko0NDQgKCjI5reEIbYhPz8f586dAwCMHDmSGjqcAAUZsVne3t5wc3NDaWkpbty4AW9vbwQEBND5DdImpVKJc+fOoaKiAhEREYiPj6drw5wEBRmxaUKhEBEREabRmVKpRGBgoMNdRE06T6/X48qVK7h27RpcXFwwevRoBAUFsV0WsSIKMmLzOBwOfHx84OHhgfLychQXF6OmpsYudqAmllVSUoL09HQ0NzejZ8+euOeee2jE7oQoyIjdEAgECAsLg0qlQnl5Oa5fvw5PT0/IZDKnuJCa/J/a2lpcuHDBtNfduHHj7GIncmIZFGTE7ri7u8PNzQ1yuRyVlZXIzs6Gr68v/P396bdxB6dSqXDx4kUUFxdDKpUiISGBphEJBRmxTxwOB97e3vD09ER1dTWqq6tRV1eHgIAA+Pj4OMRKFOT/NDc348qVK8jLy4NEIsGQIUMQERFBf88EAAUZsXNcLhcBAQHw9vZGZWUlysrKUFNTA39/f3h5edEHnZ1rampCdnY2CgoKAAB9+/ZFTEwMjbxJKxRkxCEIBAKEhITA19cXFRUVKCoqQmVlJXx8fODj40MffHamoaEB165dQ2FhIfh8PmJjYxEbGwuBQMB2acQGUZARhyIWixEREYHm5mZUVlaitLQU5eXl8PHxgb+/PzWF2Di5XI7s7GwUFxdDJBKhd+/e6Natm0OvBUm6jn46iEOSSCSIiIhAcHCw6RxaVVUVPD094e/vTx1uNsRoNKK0tBTXr19HTU0N3N3dER8fj4iICBpJk3ahICMOTSAQICgoCDKZDHV1daiqqkJOTg5cXFzoPBrLmpubcePGDeTn50OtVsPf3x9Dhw5FUFAQ/Z2QDqEgI06By+XC19cXvr6+UCqVqKqqQkFBAUpKSuDt7Q0fHx9aVd4KDAYDysvLUVRUhKqqKgBAeHg4unXr5rT70JGuoyAjTqdl00i1Wo2qqirU1NSgrKwMEokE3t7e8Pb2plAzI4ZhUFtbi6KiIpSWlkKn08Hb2xt9+/ZFSEgINXCQLqMgI05LLBYjLCwMISEhUCgUqKurQ0VFBUpLSynUuohhGMjlcpSVlaGkpARNTU1wdXVFt27dEBYWRucoiVlRkBGnx+Vy4eXlBS8vLxiNRigUCsjlclOoicViU6jR2o63ZzAYUF1djfLyclRUVECj0UAkEiEwMBChoaHw8fFhu0TioCjICLnJzaHGMIxppNZysbVYLIa7uzs8PDzg7u7u1O38Le9PTU0NampqUF9fD7VaDTc3N4SGhiIwMBDe3t7UuEEsjoKMkNvgcDjw9PSEp6cnGIaBUqlEXV0d5HI5ysvLAeCWYHPkERvDMFCpVKitrUV1dTVqa2uh0+nA4/Hg7e2NmJgY+Pv7w93dne1SiZOhICOkHTgcDqRSKaRSKSIjI6HVaqFSqaBUKqFSqVBTUwOGYcDn8+Hu7m4KN1dXV7u8FophGDQ1NUGhUKC+vh719fVQKBQA/ppC9PLyQlRUFPz8/ODp6Qkul8tyxcSZUZAR0glCodC0/BXw14d7Q0ODKdhKS0tRVFQEg8EAFxcXiEQiSCQSiMViiMVi03+zHQAGgwGNjY1oaGi45U8OhwOdTgcXFxdIpVLExsaaRqj2GM7EcVGQEWIGPB7PNGID/m9E09DQgKamJjQ1NUEul0OtVsNoNJq+TygUmoJNJBJBKBRCIBCAz+eDx+OBy+Xe8uedws9gMECv10On00Gv17e6abVaaDQaqNVqqNVqaDQaNDU1mb5XJBLB1dUVnp6eCA4ONo0qnfk8ILEPFGSEWACHw4GrqytcXV1veUyr1UKtVqO5udkUKs3NzVAqlWhqagLDMOBwODAYDLc9dkuw6fV6MAwDAK3+uy0ikQg8Hs90Xs/Pzw8SicRUJ13PRewVBRkhViYUCiEUCttcyYJhGNOIymg0wmg0wmAwtPnnzSO7lnBrGcm1jOpabjwej7oHicOiICPEhnA4HAgEAhodEdIB1GpECCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBRkhhBC7RkFGCCHErtldkGk0GvTr1w8cDgeZmZmtHrt48SJGjBgBsViM0NBQfPDBB+wUSQghxGrsLsheeeUVBAUF3XK/UqnExIkTER4ejvT0dHz44Yd466238PXXX7NQJSGEEGvhs11AR/zxxx84dOgQfv75Z/zxxx+tHtu6dSu0Wi2+//57CIVC9OzZE5mZmfjkk0/w1FNPsVQxIYQQS7ObEVllZSUWLFiAH3/8ES4uLrc8npqaipEjR0IoFJrumzRpErKzsyGXy9s8pkajgVKpbHUjhBBiX+wiyBiGwdy5c/H0009j4MCBbT6noqICAQEBre5r+bqioqLN71m9ejWkUqnpFhoaat7CCSGEWByrQbZs2TJwOJw73rKysrBu3TqoVCosX77crK+/fPlyKBQK0624uNisxyeEEGJ5rJ4jW7p0KebOnXvH50RFReHIkSNITU2FSCRq9djAgQPxyCOPYPPmzZDJZKisrGz1eMvXMpmszWOLRKJbjkkIIcS+sBpkfn5+8PPzu+vzPv/8c7z77rumr8vKyjBp0iTs2LEDQ4YMAQAkJibi9ddfh06ng0AgAAAkJSUhLi4OXl5elvkfIIQQwjq76FoMCwtr9bWbmxsAoFu3bggJCQEAPPzww1i1ahXmz5+PV199FZcvX8Znn32GtWvXWr1eQggh1mMXQdYeUqkUhw4dwqJFizBgwAD4+vpixYoV1HpPCCEOjsMwDMN2EbZCqVRCKpVCoVDAw8OD7XIIIcRpdeTz2C7a7wkhhJDboSAjhBBi1yjICCGE2DWHafYwh5bThbRUFSGEsKvlc7g9bRwUZDdRqVQAQEtVEUKIjVCpVJBKpXd8DnUt3sRoNKKsrAzu7u7gcDi3PK5UKhEaGori4mLqarwLeq/ah96n9qH3qX0c6X1iGAYqlQpBQUHgcu98FoxGZDfhcrmmC6zvxMPDw+5/SKyF3qv2ofepfeh9ah9HeZ/uNhJrQc0ehBBC7BoFGSGEELtGQdYBIpEIK1eupBXz24Heq/ah96l96H1qH2d9n6jZgxBCiF2jERkhhBC7RkFGCCHErlGQEUIIsWsUZIQQQuwaBVkHaTQa9OvXDxwOB5mZma0eu3jxIkaMGAGxWIzQ0FB88MEH7BTJkoKCAsyfPx+RkZGQSCTo1q0bVq5cCa1W2+p5zv4+tfjyyy8REREBsViMIUOG4MyZM2yXxKrVq1dj0KBBcHd3h7+/P6ZPn47s7OxWz1Gr1Vi0aBF8fHzg5uaGmTNnorKykqWKbcOaNWvA4XDw0ksvme5ztveJgqyDXnnlFQQFBd1yv1KpxMSJExEeHo709HR8+OGHeOutt/D111+zUCU7srKyYDQa8dVXX+HKlStYu3YtNmzYgNdee830HHqf/rJjxw4sWbIEK1euREZGBvr27YtJkyahqqqK7dJYc/z4cSxatAhpaWlISkqCTqfDxIkT0djYaHrO4sWL8fvvv2Pnzp04fvw4ysrKMGPGDBarZtfZs2fx1VdfoU+fPq3ud7r3iSHttn//fqZ79+7MlStXGADM+fPnTY/997//Zby8vBiNRmO679VXX2Xi4uJYqNR2fPDBB0xkZKTpa3qf/jJ48GBm0aJFpq8NBgMTFBTErF69msWqbEtVVRUDgDl+/DjDMAxTX1/PCAQCZufOnabnXLt2jQHApKamslUma1QqFRMTE8MkJSUxo0aNYl588UWGYZzzfaIRWTtVVlZiwYIF+PHHH+Hi4nLL46mpqRg5ciSEQqHpvkmTJiE7OxtyudyapdoUhUIBb29v09f0PgFarRbp6ekYP3686T4ul4vx48cjNTWVxcpsi0KhAADTz096ejp0Ol2r96179+4ICwtzyvdt0aJFmDp1aqv3A3DO94mCrB0YhsHcuXPx9NNPY+DAgW0+p6KiAgEBAa3ua/m6oqLC4jXaouvXr2PdunVYuHCh6T56n4CamhoYDIY23wdneQ/uxmg04qWXXsKwYcPQq1cvAH/9fAiFQnh6erZ6rjO+b9u3b0dGRgZWr159y2PO+D45dZAtW7YMHA7njresrCysW7cOKpUKy5cvZ7tkVrT3fbpZaWkp7r33XsyaNQsLFixgqXJirxYtWoTLly9j+/btbJdic4qLi/Hiiy9i69atEIvFbJdjE5x6G5elS5di7ty5d3xOVFQUjhw5gtTU1FvWLxs4cCAeeeQRbN68GTKZ7JauoJavZTKZWeu2tva+Ty3KysowZswYDB069JYmDkd+n9rL19cXPB6vzffBWd6DO3nuueewd+9enDhxotW2SjKZDFqtFvX19a1GG872vqWnp6Oqqgrx8fGm+wwGA06cOIEvvvgCBw8edL73ie2TdPagsLCQuXTpkul28OBBBgCza9cupri4mGGY/2ti0Gq1pu9bvny50zUxlJSUMDExMcyDDz7I6PX6Wx6n9+kvgwcPZp577jnT1waDgQkODnbqZg+j0cgsWrSICQoKYnJycm55vKWJYdeuXab7srKyHLqJoS1KpbLV59GlS5eYgQMHMo8++ihz6dIlp3yfKMg6IT8//5auxfr6eiYgIIB57LHHmMuXLzPbt29nXFxcmK+++oq9Qq2spKSEiY6OZsaNG8eUlJQw5eXlplsLep/+sn37dkYkEjGbNm1irl69yjz11FOMp6cnU1FRwXZprHnmmWcYqVTKHDt2rNXPTlNTk+k5Tz/9NBMWFsYcOXKEOXfuHJOYmMgkJiayWLVtuLlrkWGc732iIOuEtoKMYRjmwoULzPDhwxmRSMQEBwcza9asYadAlmzcuJEB0ObtZs7+PrVYt24dExYWxgiFQmbw4MFMWloa2yWx6nY/Oxs3bjQ9p7m5mXn22WcZLy8vxsXFhfnnP//Z6hclZ/X3IHO294m2cSGEEGLXnLprkRBCiP2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIIYTYNQoyQgghdo2CjBBCiF2jICOEEGLXKMgIcRDV1dWQyWR4//33TfedOnUKQqEQycnJLFZGiGXRWouEOJD9+/dj+vTpOHXqFOLi4tCvXz/84x//wCeffMJ2aYRYDAUZIQ5m0aJFOHz4MAYOHIhLly7h7Nmzt2wKS4gjoSAjxME0NzejV69eKC4uRnp6Onr37s12SYRYFJ0jI8TB5OXloaysDEajEQUFBWyXQ4jF0YiMEAei1WoxePBg9OvXD3Fxcfj0009x6dIl+Pv7s10aIRZDQUaIA3n55Zexa9cuXLhwAW5ubhg1ahSkUin27t3LdmmEWAxNLRLiII4dO4ZPP/0UP/74Izw8PMDlcvHjjz8iJSUF69evZ7s8QiyGRmSEEELsGo3ICCGE2DUKMkIIIXaNgowQQohdoyAjhBBi1yjICCGE2DUKMkIIIXaNgowQQohdoyAjhBBi1yjICCGE2DUKMkIIIXaNgowQQohdoyAjhBBi1/4ftxxTLMPJCl0AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import math\n", + "import numpy as np\n", + "import rebound, rebound.data\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "sim = rebound.Simulation()\n", + "rebound.data.add_outer_solar_system(sim) # add some particles for testing\n", + "for i in range(1,sim.N):\n", + " sim.particles[i].m *= 50.\n", + "sim.integrator = \"WHFast\" # This will end badly!\n", + "sim.dt = sim.particles[1].P * 0.002 # Timestep a small fraction of innermost planet's period\n", + "sim.move_to_com()\n", + "E0 = sim.energy() # Calculate initial energy \n", + "rebound.OrbitPlot(sim)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us integrate this system for a few hundred years. An instability will occur. We can then measure the energy error, which is a good estimate as to how accurate the integration was." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Relative energy error with WHFast: 42.011050\n" + ] + } + ], + "source": [ + "sim.integrate(600*2.*math.pi)\n", + "E1 = sim.energy()\n", + "print(\"Relative energy error with WHFast: %f\"%((E0-E1)/E0))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "An energy error that large means we basically go it wrong completely. Let's try this again but use TRACE." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Relative energy error with TRACE: -2.443053e-03\n" + ] + } + ], + "source": [ + "sim = rebound.Simulation()\n", + "rebound.data.add_outer_solar_system(sim) # add some particles for testing\n", + "for i in range(1,sim.N):\n", + " sim.particles[i].m *= 50.\n", + "sim.integrator = \"trace\" \n", + "sim.dt = sim.particles[1].P * 0.002 # Timestep a small fraction of innermost planet's period\n", + "sim.move_to_com()\n", + "E0 = sim.energy() # Calculate initial energy\n", + "\n", + "Nout = 1000\n", + "times = np.linspace(0, 600*2*math.pi, Nout)\n", + "errors = np.zeros(Nout)\n", + "\n", + "for i, t in enumerate(times):\n", + " errors[i] = (sim.energy()-E0)/E0\n", + " sim.integrate(t, exact_finish_time=0)\n", + "#sim.integrate(600*2.*math.pi)\n", + "E1 = sim.energy()\n", + "print(\"Relative energy error with TRACE: %e\"%((E1-E0)/E0))" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABhhElEQVR4nO3deXxU9b0//tfs2XeSEHYEQQiLgIS4L1FApbbWXqrWcrXFpdLl0tqqvWqXb6/etj+vbS8t1l6lrdalraIVxQVFXEAECftO2ElCAtmTmczM5/fHzDk558w5k9mSnMm8no8HD5NZz0lizivvz/vz+ViEEAJEREREScI60AdAREREFA2GFyIiIkoqDC9ERESUVBheiIiIKKkwvBAREVFSYXghIiKipMLwQkREREmF4YWIiIiSin2gDyDR/H4/Tp48iezsbFgsloE+HCIiIoqAEAKtra0oKyuD1Rq+tjLowsvJkycxYsSIgT4MIiIiisGxY8cwfPjwsI8ZdOElOzsbQODkc3JyBvhoiIiIKBItLS0YMWKEfB0PZ9CFF2moKCcnh+GFiIgoyUTS8sGGXSIiIkoqDC9ERESUVBheiIiIKKkwvBAREVFSYXghIiKipMLwQkREREmF4YWIiIiSCsMLERERJRWGFyIiIkoqDC9ERESUVBheiIiIKKkwvBAREVFSGXQbMxIRpYrmzm489+kRnG33oCwvHcPy0jGyMAMTS5N/U1ohBBraPLBbLcjPdA704ZDJMLwQEfWx9/fWY92+0+j2+XHdlDLMHJUPh80S0e654Sz/4CD+sPag6jab1YJ1P7wCw/LS43rtgXCquRMrPj6M2pYuvFp9Ur59wbQy/Par0+P+etHgwfBCRNQHTjR14uY/bsDRMx2q25/dcBQAMHV4Lv55z4Vw2KIfvd96rAn3/WMrahraAQDzJpei3ePFh/sb4PMLbD/ejHd21uLxd/bBarXgR/Mm4ubZI+M/qT6y62QL3t1dh8ff2ad7/7+2nsRH+0/jonFF+PVXpiHNYevnIySzYc8LEVEfWLu3XhVcLjt3CFz2nl+52443Y/yP38Tc/1mHupauqF77l2/twb66NnT7BLLT7Hj0xin46zcq8IVpZQCAu5/djJ/8axdaurxo6ujGX9cfgd8vcLSxAwfqW+H3i8ScZIIsfalaFVwuPXcIMpw23HXZWHxl5nAAwNmObry+7RQ+2t8wUIdJJsLKCxFRHzjSGAgu104pxc9uKEdRlgvNnd04dqYDy94/gDd31AIA9ta14v099fhqFJWRY2c6AQAPXz8JXzx/mNwTMqE0G9ga+vhdp1ow9sE35M9njcrHC3fOgT2Gqk9fkELedVOG4pLxRfjq7JEQQsBisaDD48U1k0vxnee3oLPbh0MNbfD7i2G1cggplZnjJ5eIaAC9/Plx3P3XzfiPF6tx6HRbXK/12eEzuP53H+KP6w4BACrGFKIoywUAyE13oHxYLr5/zQTcNHM4SnPSAAD3v7wdX1z2MTYfORvRezS0uQEAV51XjAJFM+tXZg7HTTOH4/IJQ+TbxhVnyR/bgxf8TUfO4rPDkb1XOEca23HzHzfg4v9+D0tfrEa72wtflFWddrcXHR4fAOCXN02VQ5zU35LhtOPqSSW44+LRAID/emMPJj2yGu/vqY/7+Cl5sfJCRCnN4/XjR//chm5f4KJ78HQbXltyccyv9/dNx7DjRIv8+fkj80IeM644C7/+yjQ8/s4+/HbNfgBA9bEmvLDxKGaOyg/7+sqLvRSKJMU5afj1V6ah0+PDt57bjCnD8/D1ylGoa+lCptOOkpw03PPcZqzdexo1De2oPKcwqnM71dyJ/33vANrdXsyfMhQ7T7Zg/aFGAMDxsyfw8pYTGJLtwurvXoJCzbHp6er2YdvxZgBAhtOGTJfxJWlyWa7ieX78ffMxXDGxOKrjp8GD4YWIUtqp5k45uADA7lMt8PtFzMMSZ9q7AQBfnF6GxZeOVV10taYMU99X3+oO+9q7Trbgnuc2AwAyw1zs0502PHP7bPlzZcgZXZgJ4DQefGU7nHYrbgr2lOipb+nCc58ehRAC88qH4l/bTuK5TwMNx58cbETF2ED4sVoAp92Krm4/Tre68dKm4/jC9LKwM55au7pR9fgHqGsJnHNnty/suc+dXIo/3jYTbW4vlr60FRsOnZGHlij1cNiIiFLaiaZA/8jIggxYLEC3T6Cx3RPz653tCDx3Xnlp2OACAFXnFeNXN03FIwsmAQDqWrrg9fkNH//8xqNyL824kuyYju+8oT3P+9E/t6HTYxwanvrwEH6zZj9++94BLPnb5zja2NOAXN/qxr+2BqYzP3nbLGz/yVxcPK4IAPDfq/fgosfew0ufHcPe2lY8/VENnt1wBC1d3fLz39pZJwcXINDQHI7NasE1k0tx7ZShsFktONPu6TXs0eDFygsRmZLX58d3XtiCfXVtGJGfjmW3zkCGM7G/sn7z7n78z7uBWS6jCjPQ1e1DfasbnxxswBemlcX0V/3ZYPDJz+h9YTWLxYKvzBqBXScDw0x7alsx7sdvYlHlKPz0hvKQx58JBqOZo/Lxh1tnRH1sAPDF84fB5wcefGU7fH6B+tYujCrM1H1sY1tPiDvU0I5DwanZSjarBZPKcuCwWZHpUk9h/uE/t6k+r2/pwtJrJgAAthwN9NzcddlYzBlTiEllkS2sl+awYXRhBg6ebse3n9+CacNz8cD889jAm2JYeSEiU9pb14o3ttfiQH0b3t97GhtrzsDr82PT4TPYcKgRXb0MM/SmqcODZWsPyJ/PGJmPocFhju++UI3HVu+J6vV8wanIUtUmmlVhS3LU/SGvbDkBIUIbX1s6A5WLWytGojjY7Bstl92GWypGYlRhBgDgpU3HsPNks+ox++pa8Ye1B/HylhO9vt4r37pQHh5SHnKBzvkrw48UjIbnpeOKicUoieJ8po8I9AVtrDmDpz6swedH428+puTCygsRmVJrl1f1+b8/85nq8+umDsWyW2KrPgCBC5/H60eG04aV916EcUOyUJKThq3HmgAAu0+1RvV69z73OVbvrJU/j6TyIinMcuG2OaNwutWN1Ttr0dLlRX2rG8XZLlgsFqzecQovfHYMHwbXOMlNd0R1bHrc3YHhqWXvH8Sy9w/i9W9fjMON7ThxthOPvqkObrfNGQWv34+VW04iw2mTA1r5sBxMHZ4nP+67VePx2eEzWHLleNwyeyRqW7rw3IYj8Pj8+Mv6I6hXDBM1tgc+LsjsvbFX68fXnYcLRufjjx8ewqHT7TjVHN06OZT8GF6IyFTqW7vw89d3y/0URlZtO4XysoOYOSofs8cURPUePr/AhkNnAAQaQc8N9o/cUjESOel2LPnbFri7fej2+fHPzcfx4YEGbKw5g+bOblw6vghPfX1WyJDStuNN8seTy3JQGOV+PD//YmCY6PJfvY/DjR2o+K81KM524X9vmYH/emOPasG7vIz4w4tfU9n519aTeDI4vVvr4vFFmDu5FA9fPxl+ITD5kbcAADbN12ByWS4+f+hq+WszpigT/3n9JHx2+Az+sv4INh4+gwW/+wi3VIyUKy+FWdHvW1SQ6cRXZ4/Exwcbceh0e9SL/FHyY3ghIlN5Y9upXoOL5L9X70GG04bqh6+B0x7ZKLjPL3Ddbz/EntpAZUW5DgoQGFYBALfXjyfe3Ydl76v3Dnp3dz021pzB6KJM1VBHmztQKXr92xdjcllOzLNgvnXFODz48nZ4/QL1rW48/VENjp1VbzGQmx7/RoX//eWpeOzNPSjKduLjA42GwQUAMoO9RulOdU+LXp+J3nmXZPd8nbafaMZ/rdoNb3A9mKIYwkvP6waqNpsOn4XbewC56Q6cMyQLTrsVLrsVTrsVLZ3d6Oz2IcNph9fnh91mhdUSmN1ks1iQ7rShrcsLr1/ALwL/Wru8yHTa4bRbceh0G84pzkKWyw6rxYJunx9+EZja7fULdPv8cNmt8Hj98PoFbFaLvIJxh8eHNIcNaQ4rvH6Bs+0eFGa5YLUEfr46PT7kpNsBqL9m2i+hRXWfBR6vH0II2G2Be4QAHDYrMpw2tHt88Pn9yElzwGKxoLnTAyECr+n1CXx8oAHZaQ6cPzIPJ5u7kJ/hQIbTDr/oOReX3YY2txdpDhuEEOjq9iM/0yFXQ60WCwoznRhdpN8r1R8YXojIVLTDRb3p8PjQ3NmNIdmRDT+0dHbLwSXbZUfVeSWq+6Ul/N1ev7wU/XlDczBrVD7+uuEIAGDhHzfAagHWfP9yjCnKhBAC7cFZO4VZzrim7/7brBG4fupQvL2zDt97sRrv7q6Dtv0lEcNGV0wsxhUTi/GvrSfx8YFG+fZ7Lj8HX68cha3HmnD3s58DADJc+nsJ2SNskh2Wn46LxhWiqaMbtc1d8rCT1QIMyYqtdwcASnMDz129s1Y1ZEd97/qpQ/G/cQzbxovhhYhidqKpE89tOIKzHR6MKszEnZeMjXvWR0cMjbgtXZGHl+7gVGSLBdj6yDUhxyuFl92nehaa+93N0zGuOFsOLwDgF4F1V8YUZcLt9csry4ZbaC1S0qqyuekONAebdC87dwjOGZKFzm5fxOcaiXnlpXhkwSS0dHqRlWbHTTOHIzfdoeojyTSY5WWNMKTZrBY89805AAIrEL9WfRJ+ITBzVD5y4xgCu3bKULy/tx5tbh/Otntgt1pgtQYqEx6vH26vDz6/QFleOk42dcJusyLLZUdDmxv5GU64HFa0dHqRn+GA3WaFzRo4J78Q6PYK+ITA6VY38jMc6PYJNHd2IzfdAZfdipYuLxw2C6wWC040daI424UMpw1HznQgL92B/Ewnun1+eH2BagyAYEXHBqvVAosFyHDY0eVV/7xrg6qA+gafL7C2jcthhc8v0NbVUyE52dwFp92KkhwXapu7kGa3IS/TIX+fhOjZimFEQTqECMze8nj9sFktsFstONnUCZ8IfM2On+1ETpoduekO1LW4kZvugNMeeN9E/gzGguGFiKImgr/Uf/DSVnmFVQB47M09OG9oDn5101SUDwu/xomRcOuOGJFm4UTC7Q2EF6fNqhu0XJodi+1WC0YWGEwlDjadtrt7qkVGF/poZbrseHfpZdhf3wqP14/zR+YnpOKi5bBZcftFY0Juz1O8l3YKdIbThg6PD5dNCL82i54LRhfggtHR9SgZKctLl0NRb/pyQbs2dyCUWCwWtLu9yAh+3N+aOjzBYSpbzOfr8fohIOCy2+QZb2ZcCJDhhYiidvezm/HWzjrd+3afasGq7aeiDi+rtp3CY6t3y5sORuNLv/8EVecV4zdfPb/XyodUeTHqkXFpbr9v7gTDxza0SuElELjSHTb5r+xEGJLtGrC/cLPSer6OUh+Q5K3vXYqPDjTgyzOMV+c1m768AGcpfuYSUXmLVZ5ihlus56v8WTdjaJEwvBBR1KSZOgAwoSQbj315ChY+uQGeYDBo6oh+hdqXPz+uG1yqzivBDdPLMGt0Pn762i6U5aVjyvAcfHroDF7ecgKeYCXl3d31+HB/A+aVl4Z9H2krAKfBjsrK8PKdK8fhrsvOMXyt3753AENy0jCxNDBbaSAvXIlWmOnC0GBPiXbNlhEFGbg5il2wiRJt8PyfRkT9psMTGCZ5d+mlGJaXgXSnDbt/Pg9/23gUD63cgX11bTjR1Imy3LSI/3qThnMkdwSHMu65/By5+rD8tpny/V86fzhe+OyY6jn1raFTZutaurDik8PocHtxbmk2GloDwcphFF4Uw0baMPK7m8/HL1btxpyxBVhZHZgR9dDKHfL9WQaNrcnIZrXgg/uugIBIaDWJKBEYXogoKm6vT65eFOekydNnbVYLCoJl681HzuKix97DrFH5eOHOObAbBAUljya8TBuRixumDwv7HGdwiqpEb72PJz84hKc/rgl9P4M9hJSVF214WTCtDAumlWHnyWa8sb025DWUQy2DQaTTz4n6G38yiSgqHe6ehtoMTXNrvmbmyKYjZ7HthHrpeSNuTRAwGtZRWnH7BZg1Kh/XTR0KILBa7LSfvo0P9p2WH3OooQ1AoPH2/JF58u1GQ1vK8GJ08Z5cloutj1yDmkevxX9UnYurJ5XgonGFuPNS4yEmIkqcwfVnAhH1ufbgkJHLbg2pqGQoKhVDsl043eqOeCaQtvJiVBlRuvCcIlx4TxH+uuEIVm07BQBo7uzGmt112HL0LN7cXovjwQXeVtw+GxePL8Lo+1cBCEx11qNsTnXYjIdLpIrTd6vG93qcRJRYDC9EFLF/bT0ph4QMZ2h/x4SSbJTmpGFkQQYsFuB0q1teebY3Hs16F9LeO5G4bspQrNp2Evvr2tDY7kFzZzf+sv6I6jHD89Mjei1lYDHqiyGigcX/M4koIsfOdODbz2+RVzI92xFaUUl32vDBDy/HC3fOQXaw/6MtwhVztZWWaNYQKch04oU7K/GjeRMBhK77MqowA8MiDC/KBmOGFyJzYuWFiCJyoimy9VekYRdp7YvIKy+B8LLy3oswqiAD+VFubAgguE8M0NTZDYslsKLoqu9cjHNLsqMKIlIjsLJHhojMg+GFiCJytj26tVukmTeR7lUkhZdMpy2m4AIAOcFVYU+c7ZSXWR9XnBV1BWXTf1ah3e1FcXbs++4QUd8xXU20qakJs2bNwvTp01FeXo6nnnpqoA+JiAB5Mz3JeM1uzFpZrkCQ6K3ycqKpE29sPyUPQ8UzPVdaPr8+uPJtmsMasjpsJHLSHBiaG9kwExH1P9NVXrKzs7Fu3TpkZGSgvb0d5eXluPHGG1FYWDjQh0aU0qTKy7/NGo4rJxZjxsj8sI+Xel7+76ManDjbiYcWTMKwvNBA8LU/fYqahnb583jCS06aeqp2XnpsFRwiMjfTVV5sNhsyMjIAAG63G0IIeXMoIho4UuWlMMuFeeVDUZwTfkhFWloeAFbvrMWqbSd1H3csuMutJJL1XYyU5qahKKsnsOTFsWMxEZlX1L8l1q1bhwULFqCsrAwWiwUrV64MecyyZcswevRopKWloaKiAhs3bozqPZqamjBt2jQMHz4c9913H4qKiqI9TCJKsDNSeImwH+W6qUPxq5umypWULp2pz16fH17NgivxVF4cNiseu3EqJg3NwcTSbNxxcehuyUSU/KIeNmpvb8e0adNwxx134MYbbwy5/8UXX8TSpUuxfPlyVFRU4IknnsDcuXOxd+9eFBcXAwCmT58Orzd0HPztt99GWVkZ8vLysHXrVtTV1eHGG2/ETTfdhJKSkhhOj4gS5WxwRdr8jMjCi8tuw1dmjcD2E834y/oj8OosOqfdzwiIf0n6qkklqJrE3xdEg1nU4WX+/PmYP3++4f2PP/44Fi9ejNtvvx0AsHz5cqxatQpPP/007r//fgBAdXV1RO9VUlKCadOm4cMPP8RNN92k+xi32w232y1/3tLSEuGZEFE0GtsC4aUgK7o+Ers1EEa6dZa07er2hdwWz7BRJBbOGoEXNx3D1Qw4REkrob8lPB4PNm/ejKqqqp43sFpRVVWF9evXR/QadXV1aG1tBQA0Nzdj3bp1mDBhguHjH330UeTm5sr/RowYEd9JEJEuqfJSEGHlRSKtWKtXeenSqbxEugt1rH56w2Qs/9pM/M/C6X36PkTUdxI626ihoQE+ny9kiKekpAR79uyJ6DWOHDmCO++8U27U/fa3v40pU6YYPv6BBx7A0qVL5c9bWloYYIgS6F9bT2LDoUacag7s2FwQ5Ros9mB4kXaiVtKrvPS1NIcN88pL+/19iShxTDdVevbs2REPKwGAy+WCy+XquwMiSnGPv7NPNZW5MMphI1tw2Mjr16m8BMNLcbYL/3vLjJBdqYmI9CQ0vBQVFcFms6Gurk51e11dHUpL+ZcOUTI6HVzwLctlx12XjkWGM7pfGw6rNGykV3kJBJo0hw2zxxTEeaRElCoS2vPidDoxc+ZMrFmzRr7N7/djzZo1qKysTORbEVE/kVatffabFfj2VeOjfr7dJlVeQsOLO1h5SXOYbskpIjKxqCsvbW1tOHDggPx5TU0NqqurUVBQgJEjR2Lp0qVYtGgRZs2ahdmzZ+OJJ55Ae3u7PPuIiJKLtEhkrG204Rt2pfAS/RL+RJS6og4vmzZtwhVXXCF/LjXLLlq0CCtWrMDChQtx+vRpPPzww6itrcX06dOxevVqrtNClKSkekmsk4DswWEj/anSwWGjGPYfIqLUFXV4ufzyy3tdrn/JkiVYsmRJzAdFROYh/e9uibH2Ig8b6VVegsNGLg4bEVEU+BuDiMISwdpLrJWXnmEj48pLLDs/E1HqYnghorDi3Rc13Aq73T4pvPBXERFFjr8xiCisuHtewjTseoIr7Ma7nxERpRb+xiCisOLteXEEe14+OdiIh1buwM6TzfJ9nmCg6ev9jIhocOFvDCLqRXw9L9JsIwD464YjeOLd/fLnUuXFYe/b/YyIaHBheCGisOTKS8wNu+pfM52env2MeiovbNglosgxvBBRWHLPS8xTpdXP8ykad7tZeSGiGDC8EFFY8gq7MeYLm1X9RL9i+pJUeXGx54WIosDfGEQUVk/lJTbaYSPl1GtpqrT2MURE4fA3BhGFFW/Pi11TefEp0oubU6WJKAb8jUFEYfVsB5KYvhTVsBHDCxHFgL8xiCiseBepU84uAgA/h42IKE78jUFE4QXDhjXG9HLBmAJcPK4IE0uzAQB+PysvRBQf/sYgorAS0bD77Dcr8MC15wFQDxt1Bzdr5Aq7RBQN/sYgorDinSotkfp2fay8EFGc7AN9AERkbvEuUieRhp2EANYfbMRTHx7C4cZ2AOx5IaLoMLwQUVjxTpWWSOHFLwTufnYzmju75ftYeSGiaPA3BhGFJSB6f1AE5GEjIVTBBWDPCxFFh78xiCishFVerD3DRheeUyjf7rRZMXZIZnwvTkQphcNGRBRWzzoviel58QuBNEdgF+kfzZuIr14wAvmZzrhem4hSCysvRBSeVHmJ82WUs428wRlHJTkuBhciihrDCxGFJfW8JKphVwjA5w9MkdbuOE1EFAmGFyIKS+55ibP2IgUVn1/AG1ycjuGFiGLB8EJEYcW7t5FEer5fCHmhOu2O00REkWB4IaKw5BV243ydnobdwHRpALBZ+SuIiKLH3xxEFJa8ykuc6UUaImLlhYjixfBCRGElqufFqhg2Ys8LEcWD4YWIDAnFDtAJ2x7A31N5YXgholgwvBCRIUV2SWjPi5dTpYkoDgwvRGRIuatRIlfYZc8LEcWD4YWIDKmGjeJ8LWlikXKFXVZeiCgWDC9EZEhdeYnvtZQr7Prlygt/BRFR9Pibg4gMqXteEjdsxMoLEcWD4YWIDAkkrmNXHjYSnG1ERPFheCEiQ6rKSwKHjbp9nG1ERLFjeCGiiCRqqjQAediIs42IKBYML0RkSF15iXNXacXzPV5WXogodgwvRGRI2fMSb8ywKH7byJUXG8MLEUWP4YWIDPVFz4uSLd4XJaKUxPBCRIZU67zEWXvRCyocNiKiWDC8EKWQ/31vP374j62qlXPDSeTGjHrP5yJ1RBQL/uYgSiG/fnsfXtp0HFuONUX0+MgiTmR0h43Y80JEMWB4IUpBXd2+iB6XyJ4XvSEiTpUmolgwvBCRsYRuDxB6G3teiCgWDC9EZEg1VTrunhfONiKixGB4ISJDInFbG4WwWAArKy9EFAOGFyIypJoqneAqSUl2WkJfj4hSB8MLUYpQTXuOsI6ifk5ivXjXnAS/IhGlCoYXohThj2Hes7rykrBDwb/NGo5RhZmJe0EiSikML0Qpwh/hwnRKidyYUSnTZU/YaxFR6mF4IUoRMWQX1WyjRMpmeCGiODC8EKWIWCovUnZJ9Izm7DRHYl+QiFIKwwtRiogjuyS8WZfDRkQUD4YXohThj2GTRSFXXhIbX2aOyk/o6xFRauGfP0QpIqaG3WDtJVHRZe0PLkddSxcmlGYn6BWJKBWx8kKUIiKNLmfbPfjNu/tx7EyHovKSmGMYXZSJirGFiXkxIkpZDC9EKUL4I3vcD/6+Ff/z7j7ctPwTRc8Ll/EnIvNgeCFKEZEOG31ysBEAUNfi7llhl9mFiEyE4YUoRcSzSB2zCxGZCcMLUYpQbg8QbY5J9DovRETxYHghShHK1XK1K+cePN2GX721B00dHvVz5MoL0wsRmQenShOlCBGm8nLtbz6E2+vH4YYO9XOkqdLMLkRkIqy8EKUIZc+Ltv/F7Q1MRdp85Kzqdva8EJEZMbwQpYhIel60oUaeKs3SCxGZCMMLUYrw+9WVl5aubuw62aJ6jIB6iEiaKs3sQkRmwp4XohSh7Xm5+vEPUNfixt8WV+g+Bui7jRmJiOLBygtRitDONqprcQMA3t5Z13O7dtiojzZmJCKKB8MLUYpQ9rz4FVsFKANLaCsMh42IyHwYXohShNFsI3Ujr0HlpU+PjIgoOgwvRCnCqMKiHE7yC3VQ4WwjIjIjhheiFGFUYRGsvBBRkmF4IUoRylyiCjLQ/zjwOXteiMh8GF6IUoRRz0u4bQN6Pmd6ISLzYHghShHGgUXofqx8HCsvRGQmDC9EKUI9bGQQZLTPkYaN+vC4iIiixfBClCKMA4vxho2svBCRGTG8EKUIo8ASrudFYmHthYhMxJR7G40ePRo5OTmwWq3Iz8/H+++/P9CHRJT0VA27fuXtPR+HDBux8kJEJmTK8AIAn3zyCbKysgb6MIgGDb9Rzwu0DbuWkPuYXYjITDhsRJQiDPcw0gwbKVfT5caMRGRGUYeXdevWYcGCBSgrK4PFYsHKlStDHrNs2TKMHj0aaWlpqKiowMaNG6N6D4vFgssuuwwXXHABnnvuuWgPkYh0GK6wq3hM6GwjIiLziXrYqL29HdOmTcMdd9yBG2+8MeT+F198EUuXLsXy5ctRUVGBJ554AnPnzsXevXtRXFwMAJg+fTq8Xm/Ic99++22UlZXho48+wrBhw3Dq1ClUVVVhypQpmDp1agynR0QSoVqkzuh27WwjrrBLROYTdXiZP38+5s+fb3j/448/jsWLF+P2228HACxfvhyrVq3C008/jfvvvx8AUF1dHfY9hg0bBgAYOnQorr32Wnz++eeG4cXtdsPtdsuft7S0RHM6RCnDbzCryOh2QLkxY58dFhFR1BLa8+LxeLB582ZUVVX1vIHViqqqKqxfvz6i12hvb0draysAoK2tDe+99x4mT55s+PhHH30Uubm58r8RI0bEdxJEg5RRhSXc0FDPxoxML0RkHgkNLw0NDfD5fCgpKVHdXlJSgtra2oheo66uDhdffDGmTZuGOXPm4Otf/zouuOACw8c/8MADaG5ulv8dO3YsrnMgGqyMd5UOG18AsPJCROZiuqnSY8eOxdatWyN+vMvlgsvl6sMjIhoc/EY9L2Ge01N5ISIyj4RWXoqKimCz2VBXV6e6va6uDqWlpYl8KyKKkjK8+CJMLz09L4wvRGQeCQ0vTqcTM2fOxJo1a+Tb/H4/1qxZg8rKykS+FRFFSZlRvIoldrUzjFTPYeWFiEwo6mGjtrY2HDhwQP68pqYG1dXVKCgowMiRI7F06VIsWrQIs2bNwuzZs/HEE0+gvb1dnn1ERAND2dvi9St7XiJ4DtMLEZlI1OFl06ZNuOKKK+TPly5dCgBYtGgRVqxYgYULF+L06dN4+OGHUVtbi+nTp2P16tUhTbxE1L+U+xn5fPrbAwDqnCJ0biMiGmhRh5fLL7+8l9kJwJIlS7BkyZKYD4qIEk85PNQdceUl8F/2vBCRmXBvI6IUoezR9al6Xoyfw40ZiciMGF6IUoZ+z0sk041YeCEiM2F4IUoRqsqLL8Jho+B/ucIuEZkJwwtRivAbzDYKmSqtyCmClRciMiGGF6IUoay8KNd58YZpetHORCIiMgOGF6IUIQxW2PWFCy+cbUREJsTwQpQilKNDXp9R867mOcH/MroQkZkwvBClCKO9jfxhKy/cVZqIzIfhhShFKDOKcpG6iCovDC9EZCIML0QpQl158Ss+1oQXnSVgOFWaiMyE4YUoRag2ZgzT8+Lx9QQbeYVdZhciMhGGF6IUoWrY1fS8KMNJtzK8cFNpIjIhhheiFOE3CC9evx92q0X3cYJNL0RkQgwvRCki4p4XBU6VJiIzYnghShFexXCQtufFKL9wqjQRmZF9oA+AiPrWWztrkZ/hVDXiejUr7IbsbxTEygsRmRHDC9EgduxMB+7662YAwH1zJ8i3ezXrvBjtLM3tAYjIjDhsRDSI1bZ0yR+7vco+F+UQkh9GpH4YG8MLEZkIwwvRIKasqHR1++SPlT0v3T7jhl2PL/Acl4O/KojIPPgbiWgQUy5M1+b2yh8rh408YSov7u7AfS47f1UQkXnwNxLRIKasqbQrwotyenS4YSNpqMlltyX82IiIYsXwQjSIKXeMbncrho0UPS9hlnmBx8vKCxGZD38jEQ1ibkVVRVV5CdPnonq+NxB4nAwvRGQi/I1ENIh5FDOM2j094aU7XLlFwc3KCxGZEH8jEQ1iyvDSZtDzEo4cXhzseSEi82B4IRrElGu7KIeNusM06SpJ4cdp468KIjIP/kYiGsRUw0aKht3IKy/BdV44bEREJsLfSESDmMfbE1iUPS/eSMOLtM4LF6kjIhPhbySiQUy5AJ1ytd2IKy8+DhsRkfnwNxLRIKYcNlKKOLx0s2GXiMyH4YVoEDMKL3qyXHYUZjpVt7HnhYjMiL+RiAYxdxThxW6zYFxxlu7zuUgdEZkJfyMRDWLhNl3UslutGFWYobqttSvQ5JvptCf0uIiI4sHwQjSI+SPsbQEAp82CUYWZqtvqW7oAAEOyXQk9LiKieDC8EA1ikUcXwG6zoiQnTXVbY7sHAFDM8EJEJsLwQjSIiSjSi91m0W3MtVqAwiyGFyIyD4YXokFMRFF7cVituuElN90Bm9WSyMMiIooLwwvRIBZFy0ug8qKznsvF44ck8IiIiOLH8EI0iEUzbOSwhVZerBbgoevOS/BRERHFh+GFKIn1PpsoimEjnZ6XZ79RgWJNEy8R0UBjeCFKUm6vD1c9/gHu/utmw8dE1bBrtSJNM2zEKdJEZEZceYooSW06fBY1De2oaWiHEAIWS2hTbbyzjcYUZRo8moho4LDyQmRy3Qar5CqDhtE2ANrZRrnpDsP3cdisqobd1799MezcTZqITIi/mYhMbE9tC8ofeQuPv7Mv5D6HIli0dHXrPl/bEqPdeFH9eurKi04hh4jIFBheiEzsF6t2w+3147dr9ofc51Ukk5ZOr+7ztcNG+WHCi10z28gXzTxrIqJ+xPBCZGL2MIvDKYeTjCov2mGjgnCVF6sFLnvPsJEnih2piYj6E8MLkYkZ9Zy0dnWrqjEtnfrhRTtTOttl3KNvt1nhsPWEpWh2pCYi6k+cbURkYk6D8PLVP27AzpMt8uctXQbDRprPXQ7jv1ccNqtqxhIrL0RkVgwvSaCpw4O8DONyP5mf3y+w8fAZnGruhNViwQ3Th8n31bd24ZmPD2PGyHzMHlOgmhFkt+lNfxaq4AIA7W6jnhd1fFEOCwFAdpodrcHg49C81/D8jAjOjIio/zG8mNyzG47gP1fuwM+/WI7b5owa6MOhCHx6qBH1rW5cP3WoXMn45+fHcd8/tsmPyXLZcdV5JQCAX7+1Fy9tOg4gMP35n/dciHHFWfjGnz/DxwcaQ17/VHNXyG1eg+Za7c1OzTouw/MzsPtUIAjZrYH7Xv/2xahv7cK44qxITpeIqN+x58Xk/nPlDgDAQ8H/krkdOt2GhX/cgG8/vwVr952Wb3+1+qTqcau2nZI/VoYRt9ePH/x9K9btOx0SXKqPNcHvF9hX1xryvj6D/hRtpNEOQw3PT5c/zkkP/C1TPiwXV04s0X09IiIzYHghSqBaRRD5/MhZ+WOvXx0utp1oBgCcbnXjSGOH6r6ahnZ8uL8h5LW/uOxjLF93UDe8GFVeQoeNjMPLjJH5uq9BRGQ2HDYyObvVYnhhIvPp8vrkj3cEAwoQumbKgfo2fHb4DG596tOQWT1urx8vbTqm+/q/XL0XX54xPOR2ozVZQiovmvCSnebALRUj0djmxoXnFOq+BhGR2TC8mFyG02Y4k4TMp9PTE0R2KJpq9QLovc99bjgd2Wi5fwD419aTIbcZBtxeel5sFgv+60tTDN+LiMiMOGxkcllh1uUg8+nq7qm8nG5149DpNtz73OfYcrQp5LH1re6Y3kMv8BhXXsLPNgqzBh4RkWnxymhyyo3yjHYOJvPoVIQXAFjwu4/Q7lHflum0hdwmGZaXjoml2Who92D68Fz8ef0R3ccVZTnR0OaRPzecbaTJOdrKi5XphYiSEMOLySn/ovb5he66H9T/1h9sRGtXN66ZXIp2txe/X3sA44uzVZUXALohZV75UJTmuvDh/gbYrRZ8rqjKpDtt+L9/vwAA8PzGoyHPXTCtDJ8caMAvb5qKO1Zskm/3Gs42Ct+wa2UYJqIkxPBicsrw0u0T0FT9aQC4vT7cseIzdHb78JuvTse+ulYse/8gAODm2SMBALNHF6D6WJPuEI/DZsF9cyfivrnAb97drwovSpk6Q4a/u/l83Qqc4bBRbz0vHDgmoiTE8GJyfsXVx+PzIx2Rpxe/X+Dnq3Zh+og81YquFJ+D9e3y8NCv3tqr2ldIqpZMKsvB0mvOxb66Vry7ux7rFGu+KKtn4cKD0T5EekOHhlOlNZ+z8kJEgwHDi8kpL0rR7jXz5o5aPPPxYQBgeEkg5Torx8926j7G5bBizthCzBlbqFrvBehZyRYIDSLKz/QqL0ZirbwwvBBRMmJ4MTllD0V3lLv8Hjvb0fuDKGrKFXFtVgvS7Fb4hEBXd8/3J00xvpfmUFfLlHsIhcsOma7Iq2zaRfB6qNOLdoVdGxt2iSgJccTb5NyKC2K04aVT0SyqXWmVYicFylsrRmLLw1djw4NXYX75UNVj0p09wUM7VGNXBIhwlY9sl8PwPgBYdssM+WOjyov2Zm2lh9mFiJIRw4uJCSFUDZ/Rhhflaq/hFj0zm1PNnfj+S1vx+Nt74Tfh6sLS1zXNYUNOmgPZaQ7kZaiDRpoisLi0lRdFYtCGB2W26K3yct3UofjhvAkAAK8vsu0BtO/HqdJElIwYXkxMe91ef+gMln9wEPUtXWhoc2PnyWb9JwYpKy8dBuuKmM2Zdg8WPrkB//z8OH773gG8uvVE1K/x8ufHcf8/t6HNHboy8d8+PYqpP3kLd/91c8zVKKkalubo+d8nL92pekxWWk+Y0VZebNbIKi+R9Lw4gq8V6fYAoZUXhhciSj7seTExv+biKu0s/ds1++H1Baoyf7h1BuZPGar3dDR19MyCaXd7UZDp1H2cWQgh8L0Xq3H0TE+vzru76vGl80P38jGy7XgTlr60FQCQl+HE/fMnyvf9+ZPDeOS1nQCA1Ttr8fGBRlw8viii112zuw4fHWjAt68cLw8bKftatLOG8jOMw4td1fNiHB60z9Mj9awYb8yo/lxbaLExvBBREmJ4MTGjwoCyivL3zccxf8pQ+P0Cr209CbvNgisnFiPDaVc1+37/pa043NiOhxdMwvVTy/r60GOy5VgT1u07Dafdih9fex4eeW0ndvRSXQICmxw+tHIHLpswBJsOn5Fv31jTCAD49FAjfvHGbmw7rn6tt3bWRhRenvm4Bj/91y4AgAWWnvCiGA7ShpA8VXgJ07CreS8LIgs2EikIRVp50VZamF2IKBkxvJiYdnVUAPjVTVPhFwIHT7fjj+sOYfepwOZ/y9cdxC9X7wUAlOS48H+LLlD1yGwMXtTv+/s2XDJuCHIzwjeDDoS3d9YBAOaXl+LKicV45LWdONXcpVqUTW+Btp/+ayfWH2rE+kONqtvrWtwQQuDmpzbIQ3D/Nms4rp5UisV/2YSPDzT0ekw7Tzbj56/vkj9fu7ce55ZkA1APG2nlZfRUubSPs6uGjXo9hLCkyotRP5RyaGzskEzD5xMRJRP2vJiYXuVlSLYLCy8YidvmjAIAnO0I7G/zWnXPTsN1LW7c8tQGHG4MnSrd2e3Db9bsR11LV0w9H7tOtmDrsaaonyc5eLoNyz84iDW760Luk6omF40rQklOGiyWwNo2je0e+b0v+MW7eOzNPfJzun1+fHrojOp15k0uBQDUtXThUEO7HFxKc9Lwky9MxsxR+QCAQw3taOkKDK39+ZPDuOix9/D9l7bCrWh0fumzY/ALYM7YAvk5p9sCGyq6VJUX9bnkpfeEQ21TrLLyEm/DrN0avvIieej6SVj93UtDKi8ML0SUjBheTKqhzY3/eXdfyO3SEITUv9LV7Ud9Sxf21AYWTlv7g8sxZVguWrq8qGloVz1Xuk49/XENKv5rDe77x7aojumj/Q249rcf4oZlH+PN7aeiPSV8sO805j/xIR57cw++8edNWLlF3YwrLf42bXgenHYrirJcAIBTTYF1Vb72f5+ioc2D5R8clJ9zpLEdHp8fGU4bXrxzDp68bSZ+d8v5sFoCfSD7FQvKvfHdS5DhtKMg04lheekAgJ0nWrD9eDMeeW0nTjR14p+fH8eK4MJ+ALBuf6A6882Lx2JobhoAyOFNNWykGQDKVYQX5XR3QB1YehsaeuVbF4a9X2r+NdyYMRhQCzOdcNqtsGr+j+dGn0SUjBheTKLT48Oy9w/gSGMgcPzbk+vx5AeHQh7nCg5BZDht8mqpUnDJdtkxuihT3l9H68YZw3H7RaPlz/+x+TgOnm6L+Bj/trFnh+NH39xjuBmgHp9f4MevbJeDBgD830c18v2tXd1o6QrMDhqeHwgWJTmB8HK6rQtd3T6cafdA60B94PjHF2ehYmwh5k4uhcNmRX5w2KY2uKDc+OIsVcPyeUMDQz/761vx+7UHVK/5RjCYeX1+uXm4fFiuPFwkBQXldGhtBlCu5aKcsg6og064qdIAcP7IfGSnGY/u9lZ5kYpr0uuGVF4YXogoCTG8mMQv39qDX721F7c/8xk6PF4cOt2u+zhpBorFYkFh8GK8P3gBL8wKfF6a69J9bpbLjkcWTMbun83D+SPzAADbNU2s4Xx+pEn++OiZDvxr20njB2tsrDmD42c7kZvuwGtLLgIA7KltkZtfTzQFltnPy3DIU4SlRdpau7whzbbS8860B4Z9hmSrz1mqikjTpR2a6UDnFGcBAD450IjVO2sBAH/6+iwAwO5TrXB7fTjZ1AWfX8Blt6I42yWHKe17AOrG2+9cOU71OG3lRZkXtBUbPeEeITXsGq2w2xNeLLqvxVEjIkpGDC8m8Y/NxwEEeire3F5r+DjlBVOqLkhDI4XBYZbi7DTd50p/wac7bSgvywUA7K5tiej4Wrq6UdsSqGLcddlYAMDTHx2OuG9m3f7AxoRXnVeMc4ZkIT/DgW6fwP66QPA6GQwv0nAOAGQFj7e1yys3JktOt7rl4wKAnHR1A7IU8n67JlBV0e7pM744UEVZvbMWQgBTh+fiqvOKkZNmh8fnx+GGDrnqMqIgA1arRf76StIMel6067NUnlOo+lyZF+IND71WXoJN39LbhKzzwvRCREmI4cUklHvO/PKtPYaPU679IQ2/SFULqRJTnK1feVFeVM8Jzjw5qtPUq+dgsLpTnO3CXZeeA5fdiu0nmvH50Z5NBz8+0IDFf9kkD9UoSZsTzh5dAIvFgpEFGQACq+kCQEOrJ+TYpbDV5vZijyZkSWFH2tE5J00dXqSwIq1QrN3TZ1yw8iK5aFwRLBYLhuYGwlN9a2AhQKBn+KooJLwoho1UU5zV5z5nbCH++o3ZuvdHskhcdprxzLDeel5Ch43U93OROiJKRgwv/aDD45UvhEZaFavB1rUYP1a5Zoh0gZYu5NKwkdFidFmK8CJVEfT6SPRImxGOKMhAQaYTXwzuUv2Morn11j99ind21eGR13aoniuEwM6TgfAxPThcJR2j9P7SjKJ8xbFLgaS1qxuNberj3BWsxPRWeZFoKy/a8HLJuMB6L9Lw0+lWN5qDwUhqvi3KUn9dlUNRvQ0FXXROke79kWSHP3xtBkYVZmD512aE3Nd75UX9ntrKi3ZxPSKiZMBfXf3gm3/ehMpH1+DYmdAqx0ufHcMFv3gXngj3HlL+tS9dkKWVdKW/0O02q+4UWGXjZ6EmPPSmPjhkJFUhbqsMTNV+e1cdOjzqZfiPnulUfX6quQttbi/sVgvGFgVCQ0Fm4HWk0CJN+S5QrI8iha22Li86gz0uowoDFZstR5sAAM2dgffODQkvxgvDKV9bMiM4fVoKL/Wt7pCqjjYUGk0z1gskquGZKCsvU4fn4YP7rsC88tCVlOUVdnvZ24iVFyIaTBhe+sEnBxvR7RP48yeHQ+774T+3yf0bADAjWJn48ozh8l/VSqrKS/DP5la5KbXn8dqLNQBkOnsu2FKFQwoNvakLHqPUTzO5LAcjCtLh8fqxQbM4nPa4pRlBo4sy5cAlVTGk8CT9t0BR3chW9LxI+zRdNbEEAPDennp0enyKgKEOIy7NwnDahl0A+PkNkwEAiypHyf0r0rBVfUto5SVds8GiLYopz0rKR2qfFu3U5UhnG1nl8MK9jYgo+TG89DFlVeLY2fD9JRlOG/62eA5euqsSv7xpqu5f8MpQIvd1eKW+DpvicaHfWmXPS6EcXrrlnZu7fX7DBty6YOWlOFh5sVgsmDEyUK3YW6uebn261Y17nt2M//f6Lvj8Qp4NNV4xVKMdNjrbrlN5kcKLu6fycsn4IpTlpqHN7cWGmkbDYSNtj4t22AgAbqscjTe+cwkevO48+TYp1DV1eEJeW/saymnGqsJKL4HAonpefOFBClB761ox74l1IZtR9nw39d+Hi9QRUTJieEmQp9YdwrW/+RAba9SrvSp7NXq7UGW57Ehz2DB7TEGYIYnQ8KL3ufbiHbi/57nS8vU+v0BLVzc+3H8akx95C/+5ckfI85TnoWxaHV0YaPo93NCuWpW2tqULb+6oxZ8+qsELnx3FgfrAbChleJH2/pGqG9KqtcqhGWlop10RXjKcNlw2YQgA4OP9DYYNu9rKi97XAwAmleWoqllSwGv3eOVjMwwvysqL4vbe4kC4yku0lKv87qltxWeH1T9/IcNGmp8rFl6IKBmZLrzs3bsX06dPl/+lp6dj5cqVA31YYR0704FfvLEbu0614Bt//gxPf1SDf209CY/Xr+op6a2vJSvMYmR6wjWl6oUfbfCR+mdau7z4xard8Hj9eO7To/jrhiMhz23tCg0JY4oC4aWmsR1n27tDngMAL2w8hn3B6dDjgou8BY49cNH1eAPVHmk2k3L/Hekxbq9fHjbKcNoxdXgeAGBffVsUPS+R/ahnBmdwdXh8aAm+tjQkpQ1ARsNGvQWCaGcbhaPdN+lAnboK1tOwK72f+vlcpI6IkpHpNmacMGECqqurAQBtbW0YPXo0rr766oE9qF5UK/b6ae3y4mfBjfwuGJ2Puy49R77P3Ut4yXZF9+2IZGhESXuhCuw87cHZDo+8ND8APLRyByrHFqpm5LR2qS/kAFCSE+h/aWxzo7FdPUPKYglcMLef6Flc7rzSnvAiHavb68Op5i60e3ywWy0YVagIL8ELc1e3T668pDutGCuFpoY2xdCOpudF87Ww6/QA6clw9lR7pD4SqQIUtvKi34+rK1x4iTZKpGlCmvL7COgtUsd1Xogo+ZkuvCi99tpruOqqq5CZGbobrpnsOBm4QM8alY8urw/u7sCy8p8dPouctKPy49yaZeK1tBfHwIUmcPW58fxhOD84I8bo8U5b+L/+tRdKqQF1y9EmaPs9//n5cfxo3kT5cym8KNcckaodzZ3ekFlLl44PDO18sC+wOF2m04axQ3rCkBQuPjt8Vt5oceyQTFWFJE2n8pLm6HmdY4pZTSE9L5qvjdGuy1pSUOnw9HyvpCpOxMNGBtWML50/DJ8dPoO5wY0jA4+N6LAMpWmaiPfVayovQr1IHWcbEdFgEPWw0bp167BgwQKUlZXBYrHoDuksW7YMo0ePRlpaGioqKrBx48aYDu6ll17CwoULY3puf5I2DpxXXorXv30J3ll6Gb4enEosrSwL9F55CXcheXzhdHknaUm4nhe9vlvty0uL3El7I80eU4CffiEwA+eA5iIoDRspp1tL1Y6Wzu6Q8DK+OAs3TC+TP9f28SiP9bWtgW0GZo0uUL2GVHnp9Pjkr12G046iLKd6Z2YLkOUMX3lRhpFwMlw92wpI7ykdh8ummW2kKrf0Pmz0PwunY919V8jVHenY4xE6bNSqarqWh40s0n+167wwvBBR8ok6vLS3t2PatGlYtmyZ7v0vvvgili5dikceeQSff/45pk2bhrlz56K+vl5+zPTp01FeXh7y7+TJnr1yWlpa8Mknn+Daa6+N4bT6l7QAnbKZdUqwL6Nbsf6Gdo8b7cyeaP8KdmoupsrP9eYMaV9fCi+Hg7tPl+ak9fSxKHak9vkF2oMXf2VfjlR58fj88iq/AHDZuUOw5MpxmFdeivJhObBbLbhTMXwWONbQH73zhuaoPpcqL8rp3OkOW3BfJ+VKvI6Q4Q9tz0tnhOFFmk7e4fHBHRyqkoKQNixaY2jYDW2Yja+BVlt5aff45AUBAeVUaUvwv5rjYXYhoiQU9bDR/PnzMX/+fMP7H3/8cSxevBi33347AGD58uVYtWoVnn76adx///0AIPe0hPPqq6/immuuQVqa/j49ErfbDbe7p9+ipSWyvXoSSVqnRbk54ARFc6rEoxm66NYsLGaNMkqGq7zorfuh/Ss7PRhepKBSkuOSw8vRxg4IIWCxWNDW1TP9Vll5yXLZYbUAftETgL5eOQo/u6Fcfsyr916Mzm5fyKJwLs1FFwBcmkAjVRXUQzjBdWKynfJeS9p+FyD0axNp5SUzWHlpd3sh/e8hBQTta9oNel6iSSDxDttoK0wA8Pq2UygfFti7StrbCAaVFw4bEVEySuhsI4/Hg82bN6OqqqrnDaxWVFVVYf369VG9VqRDRo8++ihyc3PlfyNGjIj6uOOlV3nRW6Jf+kve7fXB7xch1YCoKy9hwku3zhCV9q9safiiVl49N01e58Tj88tNslJTrNNuVVU0LBaL3Gsi9Z9oZ/3YrJaQ4ALoV1604UobcJx2q1y5GKL4WmvfEwhdUK6jO7rKi9vrDwaYnoCgPT51z4t+FaY38VY+9PprDp7uGfIT6uzCnhciGhQSGl4aGhrg8/lQUlKiur2kpAS1tcY7JWs1Nzdj48aNmDt3bq+PfeCBB9Dc3Cz/O3bsWNTHHQ+vz4+zweX5lXvf6FUD3F4/Xt92EhP+czXGPvgGHlttvAEjgF6vgtrwouwD0VZ5gNALnVR5kQzJdiHTaZMvytI6J3ozjSRScGgKPjbSFWL1ZkZpZwSlac9PceVVBkXtGi9AT++KxG+wAm2458k9L/bQKhGgvvAbtL/0KnSF3cifqyXtyH2ksWfIL2S2EXteiGgQMOVso9zcXNTV1UX0WJfLBZdLfxfl/qBc0VQ548VltyHNYUWXos/F7fXjN+/ulz9/fmPPTCQg+r+CtcMsyiEEvdk12gtVhqY6kZPmgMViQW66A2faPWju7MbQ3HT5HPV2N5YqHF3BykakZ6A33GHXjJtpKy92xfkWZfcSXjTB7BdfKg95jB6nzQqLRd3wrG2K7Tleo56XyL+P0W4HEM6YokycaOrEkcYO+P0CVqsFfs1sI23Fi9mFiJJRQisvRUVFsNlsIcGjrq4OpaWlBs9KbtKF3Wm3hiyEpr2otrm98lL5eqK9kIROle65YOsVGkKHjdQXeKmfRZ4C3SFVXkJnGkmkECKHlwjPQS+8aMNVSOXFpl950Rs2Us7oefXei+SF7XpjsVhC1k5RD5X13G416HmJJo8kcthmWF467FYL3F4/6lq7VPdJb5PutGGiYr0drvNCRMkooeHF6XRi5syZWLNmjXyb3+/HmjVrUFlZmci3Mo12d+CirbfAnN5FFQDKctMwfUReyO2J7HnRE7LOi2Z6sVRZyUlXL93fs8aLcWOsNMQSadVBbyhGu6Gj3WZV3aYMh70N0Sk3oYx0dV352BzGFS29zTKBxPW8RFO10bLbLBieHxg6OtwQ2EdLO9sIgOpnjyvsElEyijq8tLW1obq6Wp4xVFNTg+rqahw9GhgCWbp0KZ566in8+c9/xu7du3HPPfegvb1dnn002LS5Axf4TJ3wol04TVKU7cLQ3NBZVNEOIUS7wq7ROi+SkMpLp6by4go9Hzm8RFl50TtWm84quKrgoLh/SBQ9L3o7bIeTptm5W1mdMOwRGcDKizR0d/mEYowOzhY7HOx7kWYbKd9F+TVlwy4RJaOoe142bdqEK664Qv586dKlAIBFixZhxYoVWLhwIU6fPo2HH34YtbW1mD59OlavXh3SxDtYtAUrL3rhRa9SAQRmIklL6yuFNG/28t7hGnb1hEyVduiHF6kxV6q4tISpvCj3H4rkmCW6Dbs6wSDNYZPXmHEoemKk3a0B/ZCYqMpLyDYDVisAnWZo1cdR9LxEdWT6Prjvcuyta8XF44qwZndgyLYhOH1faDc3grp3KNrp+UREZhB1eLn88stDFlfTWrJkCZYsWRLzQSUTaTqt3rCR0UWzINOpWhNGEm37gfb1ox82UvdySBd8KdRIU6X1tgaQ3zN4DF6/evfi3uhVMPRuU1ZblB+PLcrCjTOG4WB9G6omhQZjZVUp0n2NJMrKi3YIybDwEs3mRkbPi1FxThqKg2FYCnLS9Pae7KI//MbKCxElI1PONkom0gJuma7eezgkRVmukJVRgfh7XpTL198wvQyvVp9U3R9u2CjLaZeHR6TbpXVopGEjvV2vQ/ZjiuNiqJ1tpL1N+bHVasHj/zbd8LVU4SXK8oK68mI840kp0hV2tbQ/Isr9pGKhrZrJs40U76Os0HGqNBElIxaN4yTNNtIbNjK60OVnOHWn34bsMNzLdUU7pKEMEo/dOBXP/PsFuHl2z6J9obtK91yYlUNCaU515UUKMVk6Aa23ak809C6kygttNL0ryqqSr5dKoZay8pKfGbronh71bKMoVthVvN5f7piNi8cXRfxcPVJ1rLHdg//7qAaHTgd6X5RHpKy8sPBCRMmIlZc4dXiC4cWpE14MLnSZLltIvwkQ/YUk3GyjdKcNV0wslnd1BsLPNlIOCWU4ArdL4UVaWl/vmLUBKp6LoV44UQbAaHpXMpx2fGFaGdrdXpTpNEeHo9ynSbmHEmA8OyfGUSNV5cWowTsa0syrd3bV4Z1dPUsWKAOVMoBxthERJSOGlzhJjap6lRSjv9JddqvusFE8s42slt6HAIw2ZgTUlZd0Z89uzkDP0vraqdWAzrBRHC2o+pUX/dlGkfjtzefHdBzKhQeVHwNhKi/KqdJRHWZsU6yN6M0IA9QhyRbJ7CkiIhPjsFGc5CXkdcKI0TCHy65feYl2zQ9lcIhk+MaieYjyGNThJVh5CYaXruB/tVOr9d43nj/k9XpT1MNG/fPjKjW7AsDowkzVfUYBKvZF6pQfxx8kjKo3ypdWVlsSucIvEVF/YXiJk0fe/ya6yot2FgsQX8Ou0YVdOTMsfM+LctgocLtUcenoDlQfdIeNEhgodGcbKW4zGoZLNGWLzI/mTVDdF8kwSzTVJ6P9kWJlND0fBpUhVl6IKBkxvMTJ7Q1c4CPZq0ficlgT3vMSyVNDh416LnRZqspLcL8ij7phV7uRIxBacYrnUqgXTuyqYaP+/3Et1qzHY9ywG1sISXR40RuO1L628hyYXYgoGTG8xMkd3Hgx0kXXAGnTRp3worn093YxU/a8RDKfRvt6BZk9S+zXt7jlj6WQIlVcOsM07GpX+Y3nAtzbbCPte/WV5V+bgbwMB/58x+yQ+4x7XmKj2ispAelFL0QD6uNThxemFyJKPmzYjZPc86KzV4/ecvdAIOgY/YUcjWgv5toLlbpnpuc+eZE6qfLSHUXPS5x784TcZo29YTdW88qHYu7kUt1+kERPlY61V8aIceVFP7AwvBBRMmLlJU7hho0cRsNGdv1ho2hFsiOwsiKjd+H96zdm48JzCvGDa3p6O7ThRZoqrXdhTGTDbq+zjfpxLXujAGJUTVMN/0TxPokOEnqz3gKv3fMxZxsRUbJj5SVOPbONomnYtRleZPqS3uFcMn4ILhk/RHWbdC4enx8+v5DPUa/yksgZQL3PNhr4C20kw0Yx97zEeExKehXAwGvrr+3C7EJEyYiVlzhJPS96Fw3lxTZNs9lfIoaNohXpcIY0HOX2+tEVHDIC1A2+Em0lIq7tAXpZpK6/ho3CMar+qBepi2a2kfI14j8/m9WiG/JUvTWKN+VUaSJKRgwvcQo3bGRTXOiUw0QuhzWifpWoLivRrYAfljSDyO31y0NGQGTTwRM928ihmio98D+uxocQ26whSx9UQXqrhg3ApC0iooTir7E4hWvYVVdeFOHFbouoXyUarZqVYCVRbusDoKfy4vH6VTON9I45tPIS/ftJettVOpH7KMUqsspL5GJt9A1H73ue6FlNREQDaeCvBklOWqRO78Jq1IxpNJ3VLJT9O9Jqs3prvACJrrzoTDdXNewO/EW3L3teEnV6fp30ohzKYnghomRn7qtoEnCHWWFXeQlRXpcTGV6G56cn7LUkyiGtpo5geDHo0dH2ocQ6TRjQv3g7B3iROq3xxVm6t6vPO8ael4S07PZeeeEMIyJKdpxtFCe550Vn9pByaX7lX7uRXoQjCQJ/vG0W7np2E75/9YReHxspZbhq6vQA0J9pBKj7eoDoqg7Thueh+liT4rnhh6UcJrjofu/qc9Ht8+P6aWWq2xMy26gPKy9c24WIBhOGlziFm22kdxFJtEllOfjwh1ca3i9i6OS1WCxw2qzw+Pw42xF+2Cik5yWK9/n9rTNw81MbcKSxw/AxA709gFaWy46f3lAecnsiel4S1QelO2ykeOkJpdkJeR8iooEy8FeDJBd22EhxDfH5+z7IJJJ0Ps0dgcqL0bBRyBBEFH/Vl+Wl4x93Xxj2MWZb58VI7CvsJnadF8P3UXw8pigTL945B+8uvawP35GIqO+w8hIHv1/A4zMOLz5FeukOPi4cM1XznXYr4Fb0vPRB5QUAhmS7sPYHlyPTpf+jqGziTeSCeImm7FeJboVd5ceJ+QF48Nrz8P9W7Ua6wyZv7aB96YqxhQl5LyKigcDwEgePIpBod1cG1JWXbl/vlZd+GGWKmBTGmjoD4cW45yX+qdKjizIN71M2BJthkTpDMe5R1Bc9L9+8ZCxumD4M1ceasPgvm0IPkIgoyZn3T9kkIA0ZAfqbJCobdru9vVdetBJxuYk1EElhTKq8GK0IrJ3enKgZM/JxKCpaRntFmUGsDbuJ3phRMiTbpframamqR0QUL1ZeYuD1+fHjV3ZgSLYLQODCoNePoWxz8UQwbGQmUhhr7oyu5yXRF0nl+jlmrryoe1ei6Hnpw/VXlEN6nGFERIMJw0uUjjZ24NXqE3hx0zH5Npfdqtuk6Y+y58VMpKnfrV2BlXuNNvwLWecl0cdhN9dsIyMWw096eV4fZgrVHkZ99zZERP2O4SVKl/7q/ZDbjC7s+RlO+eNIJhuFXMgScMWJtY1Gqry0ewLhxWHXP5j+rLw4TV15UXw8cIehYlNtwDiAB0JElGDm/VPWhIyqJ0Yr5n519ghcP3Uo/r+vTIvo9fuiYTfW15QufF3BdWyMNpLs6yX7lcHQDBszGlHNNooxKST6+68KL6aJVERE8TPv1cCEurp9urfrra4LBC68/3vLDHx55vC+PKw+IYeX4MaMRsEhdG+jxF4k1dsDmPcCHPMidQk/kh62PpjJRERkBgwvUZCqEFpGw0ZKWcF1TIwaX81GDi/B7Q+Mho1CQk0fDhuZe50XxcdRfA36cnY89zAiosHKvFcDEzKqvBgNqSg9v3gOLhpXiL/fXZnow+pFbJdH6cInrU9jdI6J3FVaj6ph18wX4z7YYDFe7HkhosGK4SUK0Q4bKU0ZnovnvjkH5cNy5dv+9s0KTC7LkT/XXmAG8nqjDQpGVY+QFXYTfJVUVV4SuBt3oql7XiJ/Xn6GE1kuOzKdNuRnOBJ6TDZOlSaiQYqzjaJgPGwU20X1wnFFWPWdSzD6/lUAzLXCrvZiZ9RvYuvzqdI9w2ymXqQuxp4Xm9WCzQ9VAUj8VHBWXohosGJ4iYLU/6HljKDnJdlow0rklZfEHkfSLFJn+EnvIumZioUtxoXziIjMzrx/yppQp0c/vGS5+ubik4ghmFirOdrKS8Q9Lwm+Rqq2BzBzeDFhUGDlhYgGK4aXKEg9L8Py0vHkbTPl27MMdkROZpH3vPTf3kamXuelj/YoioeNK+wS0SBl3quBCXUFN1ccWZCBUYUZ8u3ZaYlptDTLRQ8AbJqgYDRko50AlPjKS09Vy0xfHy1rjD0vfclmxoMiIkoAhpcoSJWXNIcVOYrAkmniystIRciKhrbQYjRslOjZRVrKnpdItlgYOPGvsJtoVhMOZRERJYJ5r7omJIWXdKcNOek94cWRoPVHtP0pibgG3nHRGDS0elB1XnFUz9NWXiJdIK4vp0qbeXjO7MNGok+XwyMi6l/mvRqYkLRgm91qRaazZzjDzBWBNIcNDy+YFPXztFkl0mbZRF+3bVYLnrxtJjo9PgzJdiX41RPHYvDxQFJman9ybWpORBQWw0sMLBZ1hcGfoAVazPIXOxDaHBvpGiR9cQ5zJ5cm/kUTTDXbyCTfR6s18T+jRERmwJ6XBHCaeOXXWEU6VRoArps6VP44VXsr+napvtgov4c+M5cHiYiiNPiuuv3oe1XjMaEkG1+rGDXQh5JwIYvUGWzMCABXn1cif2yWqkN/M2PPS5oiVBdmOQfwSIiIEovDRnH4XtW5+F7VuQl7vZCG3YS9cvS0lZdwDbuxLo0/mKj2NhrA41Cy26z48IdXwOsXyHDyf3UiGjz4G410hSxSF+ECcWapOvQ3deXFPF+EEQWxTZUnIjIzDhuZiImueapmTwBIcxr/qKirNCY6iQHCrwARUd9ieDGxgfwLXlt5SXcY799kxn6P/savARFR/2F4IV3aDRfDhhfWGjQ9L/x6EBH1JYYX0qVt2A23zgsbdll5ISLqTwwvpCua5eRVHS8peuVO0dMmIhoQDC+ky+eLIryw8qIeNkrVLwIRUT9heDGxgbwGeqNakZUXbnWAS9EvAhFRP2F4IV3RLCdvZb+HZuhswA6DiCglMLyYyMXjigb6EGTRVF5UmxKmaNWBDbtERP2HK+yawPoHrsSOEy2oOq9YdXt2mh2N7Z4BOSaf3x/xYy2Gn6QSBjgiov7CyosJDM1Nx9WTSkJm6vz+1pk4tyQLT942s9+PKbrKi+LjPjiWZMDKCxFR/2HlxcQmleXg7f+4bEDeW9nz8rMbJod9rFn39elP3CCBiKj/sPJCupSVl69Xjg77WDPuqNzfVH0/qfpFICLqJwwvpCuadV7AIRNNaEvRLwIRUT9heCFdF44rBKCeBm1EPWSSmhdu9rwQEfUf9ryQrgVTy5DusGHK8NxeH2vlkAmHzoiI+hHDC+myWi24ZnJpRI/lbCM2LRMR9ScOG1HcVENFKXrdZoAjIuo/DC8UN+7rw9lGRET9ieGF4sZ9fdi0TETUnxheKH4cMuFsIyKifsTwQnFTzbRJ0Ss3qy1ERP2H4YXiFslaMIMdKy9ERP2H4YXixmZVbd9Pin4RiIj6CcMLxY3ThMG+HyKifsTwQnHjbCNt388AHggRUQpgeKG4qS/WqXnl5lo3RET9h+GFEoBVhxQ9bSKiAcHwQnFjzwubdImI+hPDC8VNvat0al7EU/OsiYgGBsMLxc1i8HEqSdHMRkQ0IBheKG5coI1NukRE/YnhheKmmibMizgREfUxhheKGysvRETUnxheiIiIKKmYMrz8+te/xuTJk1FeXo5nn312oA+HesHKCxER9Sf7QB+A1vbt2/G3v/0NmzdvhhACV1xxBa6//nrk5eUN9KGRAfa8EBFRfzJd5WX37t2orKxEWloa0tPTMW3aNKxevXqgD4vCsCp+ilh5ISKivhZ1eFm3bh0WLFiAsrIyWCwWrFy5MuQxy5Ytw+jRo5GWloaKigps3Lgx4tcvLy/H2rVr0dTUhLNnz2Lt2rU4ceJEtIdJ/YibEhIRUX+Ketiovb0d06ZNwx133IEbb7wx5P4XX3wRS5cuxfLly1FRUYEnnngCc+fOxd69e1FcXAwAmD59Orxeb8hz3377bUyaNAnf+c53cOWVVyI3Nxdz5syBzWaL4dSov3BTQiArzY50hw1evx+FWc6BPhwiokEt6vAyf/58zJ8/3/D+xx9/HIsXL8btt98OAFi+fDlWrVqFp59+Gvfffz8AoLq6Oux73HXXXbjrrrsAAN/85jcxfvx4w8e63W643W7585aWlkhPhRJEtcJuamYX2KwWbHn4agCAw2a60VgiokElob9lPR4PNm/ejKqqqp43sFpRVVWF9evXR/w69fX1AIC9e/di48aNmDt3ruFjH330UeTm5sr/RowYEfsJUEy4MWNAmsOGNAerhEREfS2hs40aGhrg8/lQUlKiur2kpAR79uyJ+HVuuOEGNDc3IzMzE8888wzsduPDfOCBB7B06VL585aWFgaYfpfKkYWIiPqb6aZKA4iqSuNyueByufrwaKg3XOeFiIj6U0KHjYqKimCz2VBXV6e6va6uDqWlpYl8KzIRq4X7ShMRUf9JaHhxOp2YOXMm1qxZI9/m9/uxZs0aVFZWJvKtyETYsEtERP0p6mGjtrY2HDhwQP68pqYG1dXVKCgowMiRI7F06VIsWrQIs2bNwuzZs/HEE0+gvb1dnn1Egw8bdomIqD9FHV42bdqEK664Qv5capZdtGgRVqxYgYULF+L06dN4+OGHUVtbi+nTp2P16tUhTbw0eKgXqWN8ISKivhV1eLn88sshhAj7mCVLlmDJkiUxHxQlF1ZeiIioP3E1LUooFl6IiKivMbxQ3Lg9ABER9SeGF4qbcqo0Ky9ERNTXGF4obgwsRETUnxheKG7q2UYDeCBERJQSGF4oburtAZheiIiobzG8UNy4OQAREfUnhheKHzdmJCKifsTwQnFT9byw9kJERH2M4YXiZmHlhYiI+hHDC8VNuc5LLztHEBERxY3hheKmLLYIML0QEVHfYnihuHGoiIiI+hPDC8WNTbpERNSfGF4ofswuRETUjxheiIiIKKkwvBAREVFSYXihhOJUaSIi6msMLxQ3zjYiIqL+xPBCRERESYXhhYiIiJIKwwvFzWnr+THKdNkH8EiIiCgV8EpDcUtz2PCXO2bD5xfITXcM9OEQEdEgx/BCCXHpuUMG+hCIiChFcNiIiIiIkgrDCxERESUVhhciIiJKKgwvRERElFQYXoiIiCipMLwQERFRUmF4ISIioqTC8EJERERJheGFiIiIkgrDCxERESUVhhciIiJKKgwvRERElFQYXoiIiCipDLpdpYUQAICWlpaEv3Zneyv87g64O9r65PWJiIhSlXRdla7j4VhEJI9KIsePH8eIESMG+jCIiIgoBseOHcPw4cPDPmbQhRe/34+TJ08iOzsbFosloa/d0tKCESNG4NixY8jJyUnoa5tRKp1vKp0rwPMd7Hi+g9dgPlchBFpbW1FWVgarNXxXy6AbNrJarb0mtnjl5OQMuh+acFLpfFPpXAGe72DH8x28Buu55ubmRvQ4NuwSERFRUmF4ISIioqTC8BIFl8uFRx55BC6Xa6APpV+k0vmm0rkCPN/Bjuc7eKXSuYYz6Bp2iYiIaHBj5YWIiIiSCsMLERERJRWGFyIiIkoqDC9ERESUVBheIrRs2TKMHj0aaWlpqKiowMaNGwf6kKL2k5/8BBaLRfVv4sSJ8v1dXV249957UVhYiKysLHz5y19GXV2d6jWOHj2K6667DhkZGSguLsZ9990Hr9fb36eia926dViwYAHKyspgsViwcuVK1f1CCDz88MMYOnQo0tPTUVVVhf3796sec+bMGdx6663IyclBXl4evvGNb6CtrU31mG3btuGSSy5BWloaRowYgV/+8pd9fWq6ejvff//3fw/5fs+bN0/1mGQ530cffRQXXHABsrOzUVxcjC9+8YvYu3ev6jGJ+vldu3YtZsyYAZfLhXHjxmHFihV9fXohIjnfyy+/POT7e/fdd6sekyzn+4c//AFTp06VF16rrKzEm2++Kd8/mL63QO/nO5i+t31GUK9eeOEF4XQ6xdNPPy127twpFi9eLPLy8kRdXd1AH1pUHnnkETF58mRx6tQp+d/p06fl+++++24xYsQIsWbNGrFp0yYxZ84cceGFF8r3e71eUV5eLqqqqsSWLVvEG2+8IYqKisQDDzwwEKcT4o033hA//vGPxcsvvywAiFdeeUV1/2OPPSZyc3PFypUrxdatW8UXvvAFMWbMGNHZ2Sk/Zt68eWLatGliw4YN4sMPPxTjxo0TN998s3x/c3OzKCkpEbfeeqvYsWOHeP7550V6erp48skn++s0Zb2d76JFi8S8efNU3+8zZ86oHpMs5zt37lzxzDPPiB07dojq6mpx7bXXipEjR4q2tjb5MYn4+T106JDIyMgQS5cuFbt27RK/+93vhM1mE6tXrzbd+V522WVi8eLFqu9vc3NzUp7va6+9JlatWiX27dsn9u7dKx588EHhcDjEjh07hBCD63sbyfkOpu9tX2F4icDs2bPFvffeK3/u8/lEWVmZePTRRwfwqKL3yCOPiGnTpune19TUJBwOh/j73/8u37Z7924BQKxfv14IEbhYWq1WUVtbKz/mD3/4g8jJyRFut7tPjz1a2ou53+8XpaWl4le/+pV8W1NTk3C5XOL5558XQgixa9cuAUB89tln8mPefPNNYbFYxIkTJ4QQQvz+978X+fn5qvP90Y9+JCZMmNDHZxSeUXi54YYbDJ+TzOdbX18vAIgPPvhACJG4n98f/vCHYvLkyar3WrhwoZg7d25fn1JY2vMVInCB++53v2v4nGQ+XyGEyM/PF3/6058G/fdWIp2vEIP/e5sIHDbqhcfjwebNm1FVVSXfZrVaUVVVhfXr1w/gkcVm//79KCsrw9ixY3Hrrbfi6NGjAIDNmzeju7tbdZ4TJ07EyJEj5fNcv349pkyZgpKSEvkxc+fORUtLC3bu3Nm/JxKlmpoa1NbWqs4vNzcXFRUVqvPLy8vDrFmz5MdUVVXBarXi008/lR9z6aWXwul0yo+ZO3cu9u7di7Nnz/bT2URu7dq1KC4uxoQJE3DPPfegsbFRvi+Zz7e5uRkAUFBQACBxP7/r169XvYb0mIH+f117vpLnnnsORUVFKC8vxwMPPICOjg75vmQ9X5/PhxdeeAHt7e2orKwc9N9b7flKBuP3NpEG3caMidbQ0ACfz6f6IQGAkpIS7NmzZ4COKjYVFRVYsWIFJkyYgFOnTuGnP/0pLrnkEuzYsQO1tbVwOp3Iy8tTPaekpAS1tbUAgNraWt2vg3SfmUnHp3f8yvMrLi5W3W+321FQUKB6zJgxY0JeQ7ovPz+/T44/FvPmzcONN96IMWPG4ODBg3jwwQcxf/58rF+/HjabLWnP1+/343vf+x4uuugilJeXy8eSiJ9fo8e0tLSgs7MT6enpfXFKYemdLwDccsstGDVqFMrKyrBt2zb86Ec/wt69e/Hyyy8DSL7z3b59OyorK9HV1YWsrCy88sormDRpEqqrqwfl99bofIHB973tCwwvKWT+/Pnyx1OnTkVFRQVGjRqFl156Kel/kCnUV7/6VfnjKVOmYOrUqTjnnHOwdu1aXHXVVQN4ZPG59957sWPHDnz00UcDfSj9wuh877zzTvnjKVOmYOjQobjqqqtw8OBBnHPOOf19mHGbMGECqqur0dzcjH/84x9YtGgRPvjgg4E+rD5jdL6TJk0adN/bvsBho14UFRXBZrOFdLbX1dWhtLR0gI4qMfLy8nDuuefiwIEDKC0thcfjQVNTk+oxyvMsLS3V/TpI95mZdHzhvo+lpaWor69X3e/1enHmzJlB8TUYO3YsioqKcODAAQDJeb5LlizB66+/jvfffx/Dhw+Xb0/Uz6/RY3JycgYk4Budr56KigoAUH1/k+l8nU4nxo0bh5kzZ+LRRx/FtGnT8Jvf/GbQfm+NzldPsn9v+wLDSy+cTidmzpyJNWvWyLf5/X6sWbNGNT6ZjNra2nDw4EEMHToUM2fOhMPhUJ3n3r17cfToUfk8KysrsX37dtUF75133kFOTo5c7jSrMWPGoLS0VHV+LS0t+PTTT1Xn19TUhM2bN8uPee+99+D3++VfHpWVlVi3bh26u7vlx7zzzjuYMGGCqYaM9Bw/fhyNjY0YOnQogOQ6XyEElixZgldeeQXvvfdeyFBWon5+KysrVa8hPaa//1/v7Xz1VFdXA4Dq+5ss56vH7/fD7XYPuu+tEel89Qy2721CDHTHcDJ44YUXhMvlEitWrBC7du0Sd955p8jLy1N1eieD73//+2Lt2rWipqZGfPzxx6KqqkoUFRWJ+vp6IURgOuLIkSPFe++9JzZt2iQqKytFZWWl/Hxpet4111wjqqurxerVq8WQIUNMM1W6tbVVbNmyRWzZskUAEI8//rjYsmWLOHLkiBAiMFU6Ly9PvPrqq2Lbtm3ihhtu0J0qff7554tPP/1UfPTRR2L8+PGqqcNNTU2ipKRE3HbbbWLHjh3ihRdeEBkZGQMyVTrc+ba2toof/OAHYv369aKmpka8++67YsaMGWL8+PGiq6sr6c73nnvuEbm5uWLt2rWq6aMdHR3yYxLx8ytNL73vvvvE7t27xbJlywZkemlv53vgwAHxs5/9TGzatEnU1NSIV199VYwdO1ZceumlSXm+999/v/jggw9ETU2N2LZtm7j//vuFxWIRb7/9thBicH1vezvfwfa97SsMLxH63e9+J0aOHCmcTqeYPXu22LBhw0AfUtQWLlwohg4dKpxOpxg2bJhYuHChOHDggHx/Z2en+Na3viXy8/NFRkaG+NKXviROnTqleo3Dhw+L+fPni/T0dFFUVCS+//3vi+7u7v4+FV3vv/++ABDyb9GiRUKIwHTphx56SJSUlAiXyyWuuuoqsXfvXtVrNDY2iptvvllkZWWJnJwccfvtt4vW1lbVY7Zu3Souvvhi4XK5xLBhw8Rjjz3WX6eoEu58Ozo6xDXXXCOGDBkiHA6HGDVqlFi8eHFI4E6W89U7TwDimWeekR+TqJ/f999/X0yfPl04nU4xduxY1Xv0l97O9+jRo+LSSy8VBQUFwuVyiXHjxon77rtPtRaIEMlzvnfccYcYNWqUcDqdYsiQIeKqq66Sg4sQg+t7K0T48x1s39u+YhFCiP6r8xARERHFhz0vRERElFQYXoiIiCipMLwQERFRUmF4ISIioqTC8EJERERJheGFiIiIkgrDCxERESUVhhciIiJKKgwvRERElFQYXoiIiCipMLwQERFRUmF4ISIioqTy/wOQcJJqKsN7KgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(times, np.abs(errors))\n", + "plt.yscale('log')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As you can see, TRACE is able to integrate this system with much better accuracy. When a close encounter occurs, it automatically (and reversibly!) switches to the BS integrator. When there is no close encounter, you still get all the benefits in terms of speed and accuracy from WHFAST." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are a few options to adjust TRACE. First of all, because it uses BS internally, you may want to set the BS tolerances. This ensures that IAS15 never stalls while it is tries to resolve one very close encounter and can be done with the following command:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# Sets the minimal timestep to a fraction of the global timestep\n", + "sim.ri_ias15.min_dt = 1e-4 * sim.dt " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You also may want to change the critical distance at which MERCURIUS switches over from pure WHFast to IAS15. This is expressed in units of Hill radii. The default is 3 Hill radii, in the following we change it to 5 Hill radii:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "sim.ri_mercurius.r_crit_hill = 5" + ] + } + ], + "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.7.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/rebound/tests/test_trace.py b/rebound/tests/test_trace.py index 5d637c6ec..beb42879e 100644 --- a/rebound/tests/test_trace.py +++ b/rebound/tests/test_trace.py @@ -404,7 +404,7 @@ def jacobi(sim): self.assertLess(dE_trace,1e-5) # reasonable precision for trace self.assertLess(time_trace,time_ias15) # faster than ias15 -''' + def test_pericenter(self): sim = pericenter_sim() @@ -437,7 +437,7 @@ def test_independent_ho(self): ode_ho.y[0] = 1. ode_ho.y[1] = 0. # zero velocity - sim.integrate(20.*math.pi) + sim.integrate(20.*math.pi, exact_finish_time=0) self.assertLess(math.fabs(ode_ho.y[0]-1.),2e-9) # One order of magnitude higher than BS? Probably timestep self.assertLess(math.fabs(ode_ho.y[1]),2e-8) @@ -463,6 +463,6 @@ def test_trace_simulationarchive(self): x0 = sim.particles[1].x self.assertEqual(x0,x1) -''' + if __name__ == "__main__": unittest.main()