From 168d2de5a206d4ac7576f57efa8074070d1a2606 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 19:47:51 -0500 Subject: [PATCH 1/2] Fixed some inconsistencies with the new API and the id vs name dimensioning --- examples/helio_gr_test/init_cond.py | 48 ++++ examples/whm_gr_test/init_cond.py | 226 ++++++++++++++++++ .../whm_gr_test/swiftest_relativity.ipynb | 193 +++++++++++++++ python/swiftest/swiftest/io.py | 24 +- python/swiftest/swiftest/simulation_class.py | 77 ++++-- 5 files changed, 539 insertions(+), 29 deletions(-) create mode 100755 examples/helio_gr_test/init_cond.py create mode 100644 examples/whm_gr_test/init_cond.py create mode 100644 examples/whm_gr_test/swiftest_relativity.ipynb diff --git a/examples/helio_gr_test/init_cond.py b/examples/helio_gr_test/init_cond.py new file mode 100755 index 000000000..9e89dc358 --- /dev/null +++ b/examples/helio_gr_test/init_cond.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python3 +import swiftest + +sim = swiftest.Simulation() +sim.param['PL_IN'] = "pl.swiftest.in" +sim.param['TP_IN'] = "tp.swiftest.in" +sim.param['CB_IN'] = "cb.swiftest.in" +sim.param['BIN_OUT'] = "bin.swiftest.nc" + +sim.param['MU2KG'] = swiftest.MSun +sim.param['TU2S'] = swiftest.YR2S +sim.param['DU2M'] = swiftest.AU2M +sim.param['T0'] = 0.0 +sim.param['DT'] = 0.25 * swiftest.JD2S / swiftest.YR2S +sim.param['TSTOP'] = 1000.0 +sim.param['ISTEP_OUT'] = 1461 +sim.param['ISTEP_DUMP'] = 1461 +sim.param['CHK_QMIN_COORD'] = "HELIO" +sim.param['CHK_QMIN'] = swiftest.RSun / swiftest.AU2M +sim.param['CHK_QMIN_RANGE'] = f"{swiftest.RSun / swiftest.AU2M} 1000.0" +sim.param['CHK_RMIN'] = swiftest.RSun / swiftest.AU2M +sim.param['CHK_RMAX'] = 1000.0 +sim.param['CHK_EJECT'] = 1000.0 +sim.param['OUT_STAT'] = "UNKNOWN" +sim.param['IN_FORM'] = "EL" +sim.param['OUT_FORM'] = "XVEL" +sim.param['OUT_TYPE'] = "NETCDF_DOUBLE" +sim.param['RHILL_PRESENT'] = "YES" +sim.param['GR'] = 'YES' + +bodyid = { + "Sun": 0, + "Mercury": 1, + "Venus": 2, + "Earth": 3, + "Mars": 4, + "Jupiter": 5, + "Saturn": 6, + "Uranus": 7, + "Neptune": 8, +} + +for name, id in bodyid.items(): + sim.add(name, idval=id, date="2027-04-30") + +sim.save("param.swiftest.in") + + diff --git a/examples/whm_gr_test/init_cond.py b/examples/whm_gr_test/init_cond.py new file mode 100644 index 000000000..7904eb100 --- /dev/null +++ b/examples/whm_gr_test/init_cond.py @@ -0,0 +1,226 @@ +import numpy as np +import sys +from astroquery.jplhorizons import Horizons +import astropy.constants as const + +#Values from JPL Horizons +AU2M = const.au.value +GMSunSI = const.GM_sun.value +Rsun = const.R_sun.value +GC = const.G.value +JD = 86400 +year = 365.25 * JD +c = 299792458.0 +MSun_over_Mpl = [6023600.0, + 408523.71, + 328900.56, + 3098708., + 1047.3486, + 3497.898, + 22902.98, + 19412.24, + 1.35e8] + +MU2KG = GMSunSI / GC #Conversion from mass unit to kg +DU2M = AU2M #Conversion from radius unit to centimeters +TU2S = year #Conversion from time unit to seconds +GU = GC / (DU2M**3 / (MU2KG * TU2S**2)) + +GMSun = GMSunSI / (DU2M**3 / TU2S**2) + +t_print = 10.e0 * year / TU2S #output interval to print results +deltaT = 0.25 * JD / TU2S #timestep simulation +end_sim = 1.0e3 * year / TU2S + t_print #end time + +# Solar oblatenes values: From Mecheri et al. (2004), using Corbard (b) 2002 values (Table II) +J2 = 2.198e-7 * (Rsun / DU2M)**2 +J4 = -4.805e-9 * (Rsun / DU2M)**4 + +tstart = '2021-01-28' +tend = '2021-01-29' +tstep = '1d' +planetid = { + 'mercury' : '1', + 'venus' : '2', + 'earthmoon' : '3', + 'mars' : '4', + 'jupiter' : '5', + 'saturn' : '6', + 'uranus' : '7', + 'neptune' : '8', + 'plutocharon' : '9' +} +npl = 9 + +#Planet Msun/M ratio +MSun_over_Mpl = { + 'mercury' : 6023600.0, + 'venus' : 408523.71, + 'earthmoon' : 328900.56, + 'mars' : 3098708., + 'jupiter' : 1047.3486, + 'saturn' : 3497.898, + 'uranus' : 22902.98, + 'neptune' : 19412.24, + 'plutocharon' : 1.35e8 +} + +#Planet radii in meters +Rpl = { + 'mercury' : 2439.4e3, + 'venus' : 6051.8e3, + 'earthmoon' : 6371.0084e3, # Earth only for radius + 'mars' : 3389.50e3, + 'jupiter' : 69911e3, + 'saturn' : 58232.0e3, + 'uranus' : 25362.e3, + 'neptune' : 24622.e3, + 'plutocharon' : 1188.3e3 +} + +pdata = {} +plvec = {} +Rhill = {} + +for key,val in planetid.items(): + pdata[key] = Horizons(id=val, id_type='majorbody',location='@sun', + epochs={'start': tstart, 'stop': tend, + 'step': tstep}) + plvec[key] = np.array([pdata[key].vectors()['x'][0], + pdata[key].vectors()['y'][0], + pdata[key].vectors()['z'][0], + pdata[key].vectors()['vx'][0], + pdata[key].vectors()['vy'][0], + pdata[key].vectors()['vz'][0] + ]) + Rhill[key] = pdata[key].elements()['a'][0] * (3 * MSun_over_Mpl[key])**(-1.0 / 3.0) + + +if __name__ == '__main__': + # Convert from AU-day to AU-year just because I find it easier to keep track of the sim progress + for plid in plvec: + plvec[plid][3:] *= year / JD + + # Names of all output files + swifter_input = "param.swifter.in" + swifter_pl = "pl.swifter.in" + swifter_tp = "tp.swifter.in" + swifter_bin = "bin.swifter.dat" + swifter_enc = "enc.swifter.dat" + + swiftest_input = "param.swiftest.in" + swiftest_pl = "pl.swiftest.in" + swiftest_tp = "tp.swiftest.in" + swiftest_cb = "cb.swiftest.in" + swiftest_bin = "bin.swiftest.dat" + swiftest_enc = "enc.swiftest.dat" + + # Simulation start, stop, and output cadence times + t_0 = 0 # simulation start time + end_sim = 1000.0e0 * year / TU2S # simulation end time + deltaT = 0.25 * JD / TU2S # simulation step size + t_print = 1.0 * year / TU2S #output interval to print results + + iout = int(np.ceil(t_print / deltaT)) + rmin = Rsun / DU2M + rmax = 1000.0 + + #Make Swifter files + plfile = open(swifter_pl, 'w') + print(f'{npl+1} ! Planet input file generated using init_cond.py using JPL Horizons data for the major planets (and Pluto) for epoch {tstart}' ,file=plfile) + print(f'1 {GMSun}',file=plfile) + print(f'0.0 0.0 0.0',file=plfile) + print(f'0.0 0.0 0.0',file=plfile) + for i, plid in enumerate(plvec): + print(f'{i + 2} {GMSun / MSun_over_Mpl[plid]} {Rhill[plid]}', file=plfile) + print(f'{Rpl[plid] / DU2M}', file=plfile) + print(f'{plvec[plid][0]} {plvec[plid][1]} {plvec[plid][2]}', file=plfile) + print(f'{plvec[plid][3]} {plvec[plid][4]} {plvec[plid][5]}', file=plfile) + plfile.close() + + tpfile = open(swifter_tp, 'w') + print(0,file=tpfile) + tpfile.close() + + sys.stdout = open(swifter_input, "w") + print(f'! Swifter input file generated using init_cond.py') + print(f'T0 {t_0} ') + print(f'TSTOP {end_sim}') + print(f'DT {deltaT}') + print(f'PL_IN {swifter_pl}') + print(f'TP_IN {swifter_tp}') + print(f'IN_TYPE ASCII') + print(f'ISTEP_OUT {iout:d}') + print(f'ISTEP_DUMP {iout:d}') + print(f'BIN_OUT {swifter_bin}') + print(f'OUT_TYPE REAL8') + print(f'OUT_FORM EL') + print(f'OUT_STAT NEW') + print(f'J2 {J2}') + print(f'J4 {J4}') + print(f'CHK_CLOSE yes') + print(f'CHK_RMIN {rmin}') + print(f'CHK_RMAX {rmax}') + print(f'CHK_EJECT {rmax}') + print(f'CHK_QMIN {rmin}') + print(f'CHK_QMIN_COORD HELIO') + print(f'CHK_QMIN_RANGE {rmin} {rmax}') + print(f'ENC_OUT {swifter_enc}') + print(f'EXTRA_FORCE no') + print(f'BIG_DISCARD no') + print(f'RHILL_PRESENT yes') + print(f'C {c / (DU2M / TU2S)}') + + #Now make Swiftest files + cbfile = open(swiftest_cb, 'w') + print(f'{1.0}',file=cbfile) + print(f'{rmin}',file=cbfile) + print(f'{J2}',file=cbfile) + print(f'{J4}',file=cbfile) + + plfile = open(swiftest_pl, 'w') + print(npl,file=plfile) + + for i, plid in enumerate(plvec): + print(f'{i + 2} {1.0 / MSun_over_Mpl[plid]}', file=plfile) + print(f'{Rpl[plid] / DU2M}', file=plfile) + print(f'{plvec[plid][0]} {plvec[plid][1]} {plvec[plid][2]}', file=plfile) + print(f'{plvec[plid][3]} {plvec[plid][4]} {plvec[plid][5]}', file=plfile) + plfile.close() + tpfile = open(swiftest_tp, 'w') + print(0,file=tpfile) + tpfile.close() + + sys.stdout = open(swiftest_input, "w") + print(f'! Swiftest input file generated using init_cond.py') + print(f'T0 {t_0} ') + print(f'TSTOP {end_sim}') + print(f'DT {deltaT}') + print(f'CB_IN {swiftest_cb}') + print(f'PL_IN {swiftest_pl}') + print(f'TP_IN {swiftest_tp}') + print(f'IN_TYPE ASCII') + print(f'ISTEP_OUT {iout:d}') + print(f'ISTEP_DUMP {iout:d}') + print(f'BIN_OUT {swiftest_bin}') + print(f'OUT_TYPE REAL8') + print(f'OUT_FORM EL') + print(f'OUT_STAT REPLACE') + print(f'CHK_CLOSE yes') + print(f'CHK_RMIN {rmin}') + print(f'CHK_RMAX {rmax}') + print(f'CHK_EJECT {rmax}') + print(f'CHK_QMIN {rmin}') + print(f'CHK_QMIN_COORD HELIO') + print(f'CHK_QMIN_RANGE {rmin} {rmax}') + print(f'ENC_OUT {swiftest_enc}') + print(f'EXTRA_FORCE no') + print(f'BIG_DISCARD no') + print(f'ROTATION no') + print(f'GR yes') + print(f'MU2KG {MU2KG}') + print(f'DU2M {DU2M}') + print(f'TU2S {TU2S}') + + + sys.stdout = sys.__stdout__ diff --git a/examples/whm_gr_test/swiftest_relativity.ipynb b/examples/whm_gr_test/swiftest_relativity.ipynb new file mode 100644 index 000000000..ae586907c --- /dev/null +++ b/examples/whm_gr_test/swiftest_relativity.ipynb @@ -0,0 +1,193 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import swiftest\n", + "from astroquery.jplhorizons import Horizons" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swifter file param.swifter.in\n", + "Reading in time 1.000e+03\n", + "Creating Dataset\n", + "Successfully converted 1001 output frames.\n", + "Swifter simulation data stored as xarray DataSet .ds\n" + ] + } + ], + "source": [ + "swiftersim = swiftest.Simulation(param_file=\"param.swifter.in\", codename=\"Swifter\")\n", + "swiftersim.bin2xr()\n", + "swifterdat = swiftersim.ds" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file param.swiftest.in\n", + "Reading in time 1.000e+03\n", + "Creating Dataset\n", + "Successfully converted 1001 output frames.\n", + "Swiftest simulation data stored as xarray DataSet .ds\n" + ] + } + ], + "source": [ + "swiftestsim = swiftest.Simulation(param_file=\"param.swiftest.in\")\n", + "swiftestsim.bin2xr()\n", + "swiftestdat = swiftestsim.ds" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "swifterdat['varpi'] = swifterdat['omega'] + swifterdat['capom']\n", + "swiftestdat['varpi'] = swiftestdat['omega'] + swiftestdat['capom']" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "obj = Horizons(id='1', id_type='majorbody',location='@sun',\n", + " epochs={'start':'2021-01-28', 'stop':'3021-02-05',\n", + " 'step':'1y'})\n", + "el = obj.elements()\n", + "t = (el['datetime_jd']-el['datetime_jd'][0]) / 365.25\n", + "varpi_obs = el['w'] + el['Omega']" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "varpiswiftest = swiftestdat['varpi'].sel(id=2) * 180.0 / np.pi\n", + "varpiswifter = swifterdat['varpi'].sel(id=2) * 180.0 / np.pi\n", + "tsim = swiftestdat['time']" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "dvarpi_swiftest = np.diff(varpiswiftest) * 3600 * 100 \n", + "dvarpi_swifter = np.diff(varpiswifter) * 3600 * 100 \n", + "dvarpi_obs = np.diff(varpi_obs) / np.diff(t) * 3600 * 100 " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mean precession rate for Mercury long. peri. (arcsec/100 y)\n", + "JPL Horizons : 571.3210506300043\n", + "Swifter GR : 571.1981012667947\n", + "Swiftest GR : 1.5844780122245083\n", + "Obs - Swifter : 0.12294936320971743\n", + "Obs - Swiftest : 569.7365726177798\n", + "Swiftest - Swifter: -569.61362325457\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArP0lEQVR4nO3deXyV5bnv/8+VmTlzyACEeRQjRHCoioh2a1Grp93V7T6V9rSW396n7a4vh7pb3W5ObbU/W92ttW51O/xOu6mtY4/WrW1tHfqzIiCiDDEEApkISSAMCSHTdf5YixgwQBYkayV5vu/XK6+s536m616Bda3nvp/nvs3dERGR4IqLdQAiIhJbSgQiIgGnRCAiEnBKBCIiAadEICIScAmxDuBkZGZmemFhYazDEBEZVNasWVPv7llHlw/KRFBYWMjq1atjHYaIyKBiZtt7KlfTkIhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMApEYiIBJwSgYhIwEXtOQIzmw481a1oEnAH8CfgIWAkUA5c5+77ohWXiMhA5e7sbqijpnI7+yvWE1/zHvmX3UJ+wfg+PU/UEoG7lwBFAGYWD1QBzwFPAze5++tm9mXgZuD2aMUlIhJrB/btprZ8M3u2vUfb7h3E793B6KZtjGrfTb7XkhHers3j2fTRuYM3ERzlIqDM3beHrxTeCJf/HngFJQIRGUIOtTSzb0899RWl7K/aSEfDVuIO7iZnzxpGdu4nk0ZGhrftdKPe0qhPyqd2xAyqcz7PiKwJjBl/GlnjZzB3ZFqfxxerRHANsDL8+kPgCuAF4PPAuJ52MLMbgBsAxo/v22woInKqOtrb2bOrksa6Khp2bCSu6l06Ww+Svncjee2VZFkLhwf5afc4mi2F8pSZ1CbPoTR9Mgmp+WTPPJdRqVlkZ+eRHcXYLdpTVZpZElANzHb3WjObAfwEyAB+C3zD3TOOd4zi4mLXWEMiEnXu7GmopXTT+yTu205HXSnJez5iZHMlaW07SeVA16YtnshBS2FnUiH7RxbC2NOIH57O2Jlnk5FbSMrwkcc+Tz8xszXuXnx0eSyuCC4F1rp7LYC7bwYuATCzacBnYhCTiAgAba2HqK/eRmPNVg7sLKN9TwUJe7eRdqCM7I4a0mhmQXjbDjeqLIe6pALq06bTObaIhFFZZBXOInfyaaQlD6fvG3L6XiwSwbV83CyEmWW7+y4ziwO+S+gOIhGRfuOdndTXVrK7agt7y9fR2bCVpH3lpB2sYHzHdnLNye22/S7S2ZkykY1pp0NaIakFM0nInEz2+KmMHzWawd5YHdVEYGbDgYuBr3UrvtbM/jH8+lng8WjGJCJDk3d2UvbhX2naXcPBqg0k7NtByp6PSG7fz6SObWSZd7XZt3oCtXHZ7EnOZ1XGhcSlT2RYViFp+VPIzp9IdsrwqLbZR1vU+wj6gvoIROSw8k2raawqpbm2lNb9DWTuXkdqWy05nbtItI6u7do8nq2JU7CEZPakzYUxBSSmjyNz4lzGTZ6DxQ3952sHUh+BiEivdR5qpn7HJnbt2klT9WasvoQR+8rIailndOdeCq3tiO3L4idTN3wSlamXEJ8xmZH5M0nLm0R23kSmx8fHqBYDmxKBiMRcc/MBqkrWsr++ggM1pdBQxujWXWQe2kFeRw3Z5l1NM82eTGX8OCqGz6RtRD6WOYXUSfPILZzF6NGjmZw0IqZ1GYyUCEQkKvY3NlC7o4R9FRtg/05aa0sYdaCckW31jPU6plr7x9sygrq4LHamTGbrmL/BsqeTmpFL5oQZZBdMYZq+2fcpJQIR6TON9bU072ugeuNfONRYQ8fuHeTs/5DMthrG+D6mWGfXtnsZQU3iBPaMnknD6CXYuAWMyBpP/sRZjErNYpQZk2JYlyBRIhCRiFRv20zzvgYaKzbQtrOEzsYKRrdUkd1eTQ67SQXywtu2kER5wkS2pJ1H5/Asksedzpj8mSRnjCMvZyxj4od+B+1goEQgIkfwzk721Newq3wDB3aVc6huKykNmxjRUkta+y7yaOjattONhrh0DiaMoTq1mO3p0/Bh6WRMXUD2hOmMTs1ihlkMayO9oUQgEkCdHR3sqtlOw47NHNq1hc6GrSTsLSftwBZGde4lnX2kd9u+2nLYk5jNjjHz2ZY9h6T08aSOm0lqVgFZOQUAg/6hqiBTIhAZojra29m5o5TdFZto3lmK795G8v7tpLZUkttRw1hrY2x42zaPp8ay2D2skNrk0+jMmklS5kTSx81gTGYeedn5Xc09MvQoEYgMYt7Zya7qbVRvfJvWvTvprN9CclMVww41UNC2lXwOkh/etsUT2RmfS2NKAbtGn4ulT2L42KmMzJ1GzrgpjE9J1rf6gFIiEBngvLOTloNNVJW8y77ydbTWlTG6fh1pbbWkeyM51kZOeNsWT6Q+LpP9CelsSl8C+cWMzJtO1vgZZOaOpzBOt13KJykRiMSYu7OvsZ62g03UV3xE45a/wr5qUvaXM7y1gcK2rQyzdqaEt2/1BHYkTKByzBlsH5YFo3NJn3EeozJyySmYQoHusZcIKRGIRElHexu1O0rZU13KnqotJNRtYPiBcsa01jKhsxKAzPC2rR5PTVwuBxIzWJvz3/CUVBLz55I5eT75hVOZkqD/utJ39K9JpA913Xq57UP2V20OTVyybyvpB3eQ21FDnrV3dbo2ezLVCQU0JudRnX05ccPTGJY9iYwJs8ifNJsJMa2JBIkSgUiEDjbtp7GmjJ1bP6Szci1JByoY1lRJ2qFqRngT6dbWdetlqydQHZ9LQ8oEqscsIj5rKsNzJpOVP5mc8dOZEoARL2XgUyIQ6YF3dtBQs4Oq0rUc3LGOlF3vMfrQTsa01pLG/q6JSzrc2EkGdQm51I0+C0seBWmFDMudQeaEWeSMm0phYiKFsa6QyHEoEUhgtbe1sWPrJurK3oN91cRXryax7QBjDlUxtmMnmdbW1WZfRQ47Ewuoy5iJjc6lM3USo7LGMb5oEdkpw8jXUAkyiCkRyJDW2LiHnVvX07j5zdCH/cE6Upu2kdNeRaK3M8lauwY2a2AM++LS2JMyjppRnwrNUpUzhbwZC8jLzidfzTgyRCkRyKC3t6GWmqrt7N+xHna+D/t2En+wgbFt28mjntTwdoc8kd2Wyt7ETDalX4wnDicpdxYZE4tITs8jO28iGbr1UgIoaonAzKYDT3UrmgTcAfyZ0IT1KUA78A/uvipacckg4M7e+ip2bV1P465K2nd9ROLebYxoriKzvYYsGhkT3rTV46mPy+BQwmh2pZ5BddpUEnNnkjVtIbnjp5IbF3fEpOQiEsVE4O4lQBGAmcUDVcBzwCPAv7r7y2Z2GfBDYFG04pKBoaOjk/rGRho2/JlDtSX47nIS91cw6mAlOR07GcOhrg/7Djd2xmWzN2ksO0adw9aM6QzPHMeYcbPInlxEXsqwmNZFZLCJVdPQRUCZu283MwdGh8vHANUxikn6WUd7Owf276Hmo7Xs31lGR80H0NFG+p73KWwrI8c6uoZKaPJkauLGUp+cR1X6WVhaIcNyppKRN4Gc8VPJH5HWNYaOiJyaWCWCa4CV4df/BLxiZvcCccA5Pe1gZjcANwCMH6+hsQaqto5OKnfuoqZ8E4l1m+ioeZ/4lt2Maqogv62cMXbwiGacTuLYmjyDtXnXkDwyncTc2aRNO5fMnDymJKoLSyQazN2je0KzJELf+me7e62Z/QR43d2fMbO/BW5w9yXHO0ZxcbGvXr06GuFKD9raO6ip2k59xWaaakrpbNhG4r7tjGmpYmxHDRm2r2vbQ57IXhtFXVI+B0ZOpDN9KsmZE0gtmMGEGfOJi4vDdDeOSFSY2Rp3Lz66PBZfuS4F1rp7bXj5euCb4de/AR6NQUxylNZDLdRWlLK7ooTm2lK8YRspB3aQ2lLF2M5axtuhriGLO9yoj8+iMTmfylEXUpU+kRFjpxKXNZXx0+eTnZBAdkxrIyLHE4tEcC0fNwtB6OrgAkJ3Dy0GSmMQUyC1t7WyfdNqDjXvY9/WNdjuMlKaKhjZWseE9u2Ms07Ghbc96EnUJoxl3/ACNo4+h7j0SYzMnUrmuBmk5U8mJyG5q31fRAaXqCYCMxsOXAx8rVvxV4F/M7MEoIVwP4D0jfb2DrZ88Ff2NdTQuXsrvqeCEXs2knGogqzOBiZbe9e2zZ5MTXwezYlpvDv2AuIzpzAqbyrZ42eSnjOOQjXhiAxJUU0E7t4MZBxV9hYwP5pxDDVtrYeo3LKexsoSDtV+BE31jNizifxDW0jwdmZY88fbejzbEwrZNWwyFaOXkJB/OgnDx5A/8yzSswuYrOGNRQJH/+sHiX2NDTRUb+VAQzVNFR8QV7eJUQe2kdNWwShvYqJ1dG17eOKSsjFn0ZkwnKRx80kfN50RuVMZNiqdKaNSY1cRERlwlAgGkKb9jdRs3cC+6o9ord9GfMNHDDtYQ/ahHWSzu+thC4A9jKYmcTxb0s6nY1gGibmzGZZeQMGshSSnDGPKsBExq4eIDC5KBFHk7uzZ3UBDxSYat66lfU8F8furGNVUTlpbHWOp65qOEEKDoNUl5LF9TDFlGTNIyphA8ugscqfNJyOngLSY1UREhhIlgj7mnZ007KqkqbGe+q3v01qzgYSD9WTv+4DU1p2kc6Br0pJONxoslfrEPKpHzmRb+lUk585idN40ssZPJyMt88gOFRGRfqBEcBK8s5Nd1dvYu6uC/TvLaCt/h8TmWka21DC6fTe51JEJTADaPY4DDKcyYRyVaUtIyphAQuZEsqacSfa4KWQlp5AV6wqJSKApERzD4blnq0tWc3B3FS27q0mo30Tewc2M6DxADo1d9823egK1cdnsS8ykZuQstucWE5cyivTJxWQXziQ1NbNrKGQRkYEm8ImgsX4nO7d+wP7qj2iv30LS3nJGN+8gp6OadJq7mnEA9tgYKobNoC5hFFty55GUmkfahDlk5E1iXKoacURkcApUIlj3h5W0lPyR5KZqhrXtIaetkjT2dX1b73CjNi6b+uQCNqXNxdMnMTx3JmnjppOelU/amHR10IrIkBOoRNDy0Z84bdf/oS4+i6aEdErTL6AzYyrDcqeTPm4mOeOnk5ecQl6sAxURiaJAJYL5X/kpCQkPaagEEZFuApUIEpOSYx2CiMiAo6/GIiIBp0QgIhJwSgQiIgGnRCAiEnBKBCIiAadEICIScEoEIiIBF7XnCMxsOvBUt6JJwB3A2cD0cFkq0OjuRdGKS0Qk6KKWCNy9BCgCMLN4oAp4zt3vP7yNmf0I2ButmEREJHZPFl8ElLn79sMFZmbA3wKLYxSTiEggxaqP4Bpg5VFl5wG17l4ag3hERAIr6onAzJKAK4DfHLXqWj6ZHLrvd4OZrTaz1XV1df0ZoohIoMTiiuBSYK271x4uMLME4GqO7Ew+grs/7O7F7l6claXJHUVE+kosEkFP3/yXAJvdvTIG8YiIBFrEicDMRoTv+omYmQ0HLgaePWpVT30GIiISBSe8a8jM4gh9UF8HnAkcApLNrA74HfBwbzt43b0Z+MTkvu6+LIKYRUSkD/XmiuBPwGTgNmCsu49z92xCd/n8FbjbzP6+H2MUEZF+1JvnCJa4e9vRhe6+G3gGeMbMEvs8MhERiYoTXhH0lAROZhsRERmYev1ksZnd2EPxXmCNu6/rs4hERCSqIrlrqBhYDuSHf24AFgGPmNktfR+aiIhEQyRjDWUA89z9AICZ/QvwNHA+sAb4Yd+HJyIi/S2SK4LxQGu35TZggrsfJHRLqYiIDEKRXBH8J/BXM3shvHw5sNLMRgAb+zwyERGJil4nAnf/X2b2O+BTgAHL3X11ePV1/RGciIj0v143DYXnC5gJjAlPJrPLzBb0V2AiIhIdkfQRPEhoWslrw8v7gZ/1eUQiIhJVkfQRLHT3eWb2HoC77wnPLSAiIoNYJFcEbeFRRx3AzLKAzn6JSkREoiaSK4KfAM8BOWZ2F/A54Lv9EpWIDGptbW1UVlbS0tIS61ACKSUlhYKCAhITezcMXCR3Df3SzNYQmnge4LPuvukkYhSRIa6yspJRo0ZRWFhI6D4TiRZ3p6GhgcrKSiZOnNirfXozH0FPYwwBXGpml7r7jyMJUkSGvpaWFiWBGDEzMjIyiGRu995cEYwK/55OaGKa34aXLwfeiChCEQkMJYHYifS9P2EicPd/DR/4VUJjDe0PL98J/CbyEEVEZCA5lbGGWoHC3u5sZtPNbF23n31m9k/hdV83sxIz22BmGrxORPrEyJEjKS8vZ9iwYRQVFTFr1iyWL19OZ2cn5eXlzJkz57j733nnndx7771HlBUWFlJfXx9RHJdddhmNjY2Rhh81kdw19L+BVWb2HKFbSK8Cnuztzu5eAhQBhG9DrQKeM7MLgSuBue5+yMyyI4hJROSEJk+ezLp162hvb2fx4sU8//zzzJs3r9/P6+64O7/73e/6/VynotdXBO5+F/AlYA/QCHzJ3X9wkue9CChz9+3A/wPc7e6HwufZdZLHFBE5roSEBM455xy2bNnSJ8f78Y9/zJw5c5gzZw73338/AOXl5cycOZN/+Id/YN68eVRUVHRdRTz00EMUFRVRVFTExIkTufDCCwFYuXIlp512GnPmzOHWW2/tOv7IkSP5zne+w+mnn85ZZ51FbW0tAL/5zW+YM2cOp59+Oueff/4p16M3dw2ZuzuAu68F1h5vm166BlgZfj0NOC/8bEILcJO7vxvBsURkAPvX/7OBjdX7+vSYs/JG8y+Xz454v+bmZv74xz+yYsWKXu9z33338Ytf/KJrubq6GoA1a9bw+OOP88477+DuLFy4kAsuuIC0tDRKSkp4/PHHefDBB4841vLly1m+fDltbW0sXryYG2+8kerqam699VbWrFlDWloal1xyCc8//zyf/exnaWpq4qyzzuKuu+7illtu4ZFHHuG73/0uK1as4JVXXiE/P79Pmpx6c0Xwp3Ab/vjuhWaWZGaLzexJ4PrenjA8LMUVfNzRnACkAWcBNwO/th66vM3sBjNbbWarI7ktSkSkrKyMoqIizj33XD7zmc9w6aWX9nrfb33rW6xbt67rJy8vD4C33nqLq666ihEjRjBy5Eiuvvpq3nzzTQAmTJjAWWeddcxjfvOb32Tx4sVcfvnlvPvuuyxatIisrCwSEhK47rrreOON0A2ZSUlJLF26FID58+dTXl4OwLnnnsuyZct45JFH6OjoOJm35Ai96SP4G+DLhOYemEioWWgYoSTyKnBfhHMWXwqsdffa8HIl8Gz4imKVmXUCmcARn/bu/jDwMEBxcXEkVx8iEkMn8829rx3uI+hLx2sEGTFixDHXPfHEE2zfvp0HHnjghMdJTEzsuhU0Pj6e9vZ2AB566CHeeecdXnrpJYqKili3bh0ZGRknUw2gF1cE7t7i7g+6+7nABELt+2e4+wR3/+pJTFx/LR83CwE8DywGMLNpQBIQWZe8iEiUnX/++Tz//PM0NzfT1NTEc889x3nnnXfcfdasWcO9997LL37xC+LiQh+/Cxcu5PXXX6e+vp6Ojg5WrlzJBRdccNzjlJWVsXDhQlasWEFmZiYVFRWnVJdI7hrC3duAmpM9mZkNBy4Gvtat+DHgMTP7kNAtqddH2N8gIvIJ7e3tJCcnH3ebkpISCgoKupbvu+8+Pv/5z/fq+PPmzWPZsmUsWBCaluUrX/kKZ5xxRlfzTU8eeOABdu/e3dVJXFxczKOPPsoPfvADLrzwQtydyy67jCuvvPK457755pspLS3F3bnooos4/fTTexXzsdhg/MwtLi721atXn3hDEYmJTZs2MXPmzJjG8P777/PVr36VVatWxTSOWOnpb2Bma9y9+OhtI3mgTERkUHjooYe49tpr+d73vhfrUAaFSKaq/B/9GYiISF9Zvnw5Gzdu5JJLLol1KINCJH0EPzKz64B2YBWw0t039E9YIiISLZE0DTUA3wN+DBwgdL//146/i4iIDHSRXBHsdffXwq//y8z+DXgH+Pe+D0tERKIlottHAczsVkIPfI0B9vd5RCIiElUnc9fQM8AWoAD4ft+GIyLSN+666y5mz57N3LlzKSoq4p133unVfnfccQd/+MMfAHjzzTeZPXs2RUVFvP32230yimhtbS1/93d/x6RJk5g/fz5nn302zz33HAB//vOfGTNmDGeccQYzZszgpptuOuXz9UYkiSDNzMa5+xZ3/3dgKUoEIjIAvf3227z44ousXbuW9evX84c//IFx48b1at8VK1awZMkSAH75y19y0003sW7dOkpKSiJOBIeHhDjM3fnsZz/L+eefz9atW1mzZg2/+tWvqKys7NrmvPPO47333uO9997jxRdf5C9/+UtE5zwZkTQNjQb+bGb1wEYgFTj10Y5ERPpYTU0NmZmZXU8WZ2ZmArBq1Sruvvtunn32WV544QWuueYa9u7dS2dnJ7NmzWLr1q0sW7aMpUuX0tjYyK9//WteeeUVXn31Vf7yl79w8OBB3nrrLW677TaWLl3K17/+dT744APa29u58847ufLKK3niiSd46aWXaGlpoampiddee60rrtdee42kpCSWL1/eVTZhwgS+/vWvf6IOhyfTqaqq6ud3K7JEcCHwIbCQ0PzFDrzUH0GJyBDy8rdh5wd9e8yxp8Gldx9z9SWXXMKKFSuYNm0aS5Ys4Qtf+AIXXHAB8+bN47333gNCzT5z5szh3Xffpb29nYULFx5xjK985Su89dZbLF26lM997nM88cQTrF69umuwuH/+539m8eLFPPbYYzQ2NrJgwYKuK4m3336b9evXk56efsQxN2zY0OsJcfbs2UNpaWmfzDdwIpFMTLPe3Tvd/W13f8Ldn3R3DQ4nIgPOyJEjWbNmDQ8//DBZWVl84Qtf4IknniAhIYEpU6awadMmVq1axY033sgbb7zBm2++ecIB44726quvcvfdd1NUVMSiRYtoaWlhx44dAFx88cWfSAI9+cd//EdOP/10zjzzzK6yN998k7lz5zJ27FiWLl3K2LFjI6v8SYj4riERkYgc55t7f4qPj2fRokUsWrSI0047jSeffJJly5Zx3nnn8fLLL5OYmMiSJUtYtmwZHR0dn5ib+ETcnWeeeYbp06cfUf7OO+8ccxjq2bNn88wzz3Qt/+xnP6O+vp7i4o+H/znvvPN48cUX+eijj/jUpz7FVVddRVFRUUSxRUpjDYnIkFNSUkJpaWnX8rp165gwYQIQGj76/vvv5+yzzyYrK4uGhgY2b97M7NnHnzdh1KhR7N//8R3zn/70p/npT3/aNZ/A4San41m8eDEtLS38/Oc/7yprbm7ucdtp06Zx2223cc8995zwuKcqkrGG/qeZpfVnMCIifeHAgQNcf/31zJo1i7lz57Jx40buvPNOIDT+f21tbVfb+9y5c5k7dy49TIx4hAsvvJCNGzdSVFTEU089xe23305bWxtz585lzpw53H777SeMy8x4/vnnef3115k4cSILFizg+uuvP+aH/fLly3njjTfYtm1bZG9AhHo9DLWZfY/QXMNrCc0h8Eqs5g3QMNQiA9tAGIY66PplGGp3/y4wFfgPYBlQambfN7PJpxauiIjEUkR9BOErgJ3hn3ZCk84/bWY/7IfYREQkCnp915CZfQO4ntB8wo8CN7t7m5nFAaXALf0TooiI9KdeJQIL9aKcDlzt7tu7r3P3TjNb2h/BiYhI/+tVInB3N7Mzjk4C3dZvOtExzGw68FS3oknAHYSGqvgqUBcu/2d3P/WRnUREpFci6SN428zOPPFmPXP3EncvcvciYD7QDDwXXn3f4XVKAiIi0RVJIrgQ+KuZlZnZejP7wMzWn+R5LwLKjnWFISJyqgbSMNSNjY08+OCDx1wf66GpI0kElxJqzlkMXE5oGOrLT/K81wAruy3/z3ByeexYD62Z2Q1mttrMVtfV1fW0iYgIMHCGoT7seIlgIAxNHUkiuP4YPxExsyTgCuA34aKfA5OBIqAG+FFP+7n7w+5e7O7FWVlZkZ5WRAKkp2Go8/LyWLVqFVdffTUAL7zwAsOGDaO1tZWWlhYmTZoEwLJly3j66ad59NFH+fWvf82KFSu49tprueOOO3jqqae6nixuamriy1/+MmeeeSZnnHEGL7zwAhAaYXTBggUUFRUxd+5cSktL+fa3v01ZWRlFRUXcfPPNR8Q6EIamjmTQuaZur1MIXRGcsJO4B5cCa929FuDwbwAzewR48SSOKSID1D2r7mHz7s19eswZ6TO4dcGtx1wfy2GoH3roIb75zW9y3XXX0draSkdHB3fffTcffvgh69at+0SsA2Fo6kieLP5Rt5+7gEVA/kmc81q6NQuZWW63dVcRmvNAROSkxXIY6rPPPpvvf//73HPPPWzfvp1hw4ZFdNxYDE19KsNQDyfUZ9BrZjYcuBj4WrfiH5pZEaGJbsqPWicig9zxvrn3p1gNQz1z5kwWLlzISy+9xKc//WkeffTRrmanngyEoakjGX30g3CH7noz2wCUAP8WycncvdndM9x9b7ey/+7up7n7XHe/wt1rIjmmiMjRYjkM9datW5k0aRLf+MY3uOKKK1i/fv0n9u1uIAxNHUln8eG7hC4HLgHy3P2BPo1GRKQPxHIY6qeeeoo5c+ZQVFTE5s2b+eIXv0hGRgbnnnsuc+bM+URn8UAYmrrXw1APJBqGWmRg0zDUsdcvw1Cb2ZNmltptOc3MHjuVQEVEJPYiaRqa6+6NhxfcfQ9wRp9HJCIiURVJIojr/tSvmaVzancdicgQNhibnYeKSN/7SD7IfwT8/2b2NKFbPf8WuCuis4lIIKSkpNDQ0EBGRsYJO2Glb7k7DQ0NpKSk9HqfSOYj+BOwmtBYQ0ZoboKNJxOoiAxtBQUFVFZWonHBYiMlJYWCgoJebx/JfATPu/t8QB/+InJciYmJTJw4MdZhSC9F0kfw11OZj0BERAamSPoILgSWm1k5oQHojNDFwtz+CExERKIjkkRwab9FISIiMRNJ09AO4Dzg+vDMYg7k9EtUIiISNZEkggeBswkNIw2wH/hZn0ckIiJRFUnT0EJ3n2dm70HoyeLwbGMiIjKIRXJF0GZm8YSahDCzLKCzX6ISEZGoiSQR/AR4Dsgxs7uAt4Af9EtUIiISNb1uGnL3X5rZGuCicNGV7t63E5GKiEjUnTARmNlvjy4K//60meHuV/R9WCIiEi29uSI4G6ggNOH8O3ycCCJiZtOBp7oVTQLucPf7w+tvAv5fIMvd60/mHCIiErneJIKxhCacvxb4O+AlYKW7b4jkRO5eAhQBhDudqwj1OWBm48Ln2BHJMUVE5NSdsLPY3Tvc/b/c/XrgLGAL8Gcz+/opnPcioCz8YBrAfcAthO9IEhGR6OntMNTJwGcIXRUUErqD6NlTOO81hJqaMLMrgCp3f/9445ab2Q3ADQDjx48/hVOLiEh3J5y83syeBOYALwO/cvcPT+mEoYfQqoHZhJ5O/hNwibvvDQ9oV3yiPgJNXi8iErljTV7fmyuC/05otNFpwDe6fWs/PPro6AhjuRRY6+61ZnYaMBE4fDVQAKw1swXuvjPC44qIyEk4YSJw90geOuuNawk3C7n7B0D24RW9vSIQEZG+09cf8sdlZsMJ3R10Kv0LIiLShyIZdO6UuXszkHGc9YXRi0ZERCDKVwQiIjLwKBGIiAScEoGISMApEYiIBJwSgYhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMApEYiIBJwSgYhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMApEYiIBJwSgYhIwCkRiIgEXNSmqjSz6cBT3YomAXcQmrrySqAT2AUsc/fqaMUlIhJ0UUsE7l4CFAGYWTxQBTwH7HH328Pl3yCUHJZHKy4RkaCL6uT13VwElLn79qPKRwAeg3hERAIrVongGmDl4QUzuwv4IrAXuLCnHczsBuAGgPHjx0chRBGRYDD36H4BN7MkoBqY7e61R627DUhx93853jGKi4t99erV/RiliMjQY2Zr3L346PJY3DV0KbD26CQQ9p/Af4tyPCIigRaLRHAtRzYLTe227gpgc9QjEhEJsKj2EZjZcOBi4Gvdiu8O31raCWxHdwyJiERVVBOBuzcTem6ge5magkREYkhPFouIBJwSgYhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMApEYiIBJwSgYhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMApEYiIBJwSgYhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMBFbarK8LzET3UrmgTcAeQDlwOtQBnwJXdvjFZcIiJBF7UrAncvcfcidy8C5gPNwHPA74E57j4X+Ai4LVoxiYhI7JqGLgLK3H27u7/q7u3h8r8CBTGKSUQkkGKVCK4BVvZQ/mXg5Z52MLMbzGy1ma2uq6vr1+BERIIk6onAzJKAK4DfHFX+HaAd+GVP+7n7w+5e7O7FWVlZ/R+oiEhARK2zuJtLgbXuXnu4wMyuB5YCF7m7xyAmEZHAikUiuJZuzUJm9jfArcAF7t4cg3hERAItqk1DZjYcuBh4tlvxA8Ao4Pdmts7MHopmTCIiQRfVK4LwN/6Mo8qmROv8D73/EL/b9ju6tz453rV8+LXTbf1Ry0c77rpjtHIdb5/jr4r8eMdraTtuHH18vD5/n6JlQIQQ+yAGRAwDoNV4ILwP9194P+fkndOnx4xF01DMZA/PZmrqVMwMwwBCv42Pl8PrDMPMuvY9vL4n3bf7xLpj7Hcy+xzPyR7vZOvV18c7mfcpWk7m7zEU6W8xMGIYO2Jsnx8zUIng6qlXc/XUq2MdhojIgKKxhkREAk6JQEQk4JQIREQCTolARCTglAhERAJOiUBEJOCUCEREAk6JQEQk4GwgPLYdKTOrA7af5O6ZQH0fhjMYqM7BoDoHw6nUeYK7f2Ic/0GZCE6Fma129+JYxxFNqnMwqM7B0B91VtOQiEjAKRGIiARcEBPBw7EOIAZU52BQnYOhz+scuD4CERE5UhCvCEREpBslAhGRgAtUIjCzvzGzEjPbYmbfjnU8fcHMxpnZn8xsk5ltMLNvhsvTzez3ZlYa/p3WbZ/bwu9BiZl9OnbRnxozizez98zsxfDykK6zmaWa2dNmtjn89z47AHX+Vvjf9YdmttLMUoZanc3sMTPbZWYfdiuLuI5mNt/MPgiv+4lFMqWcuwfiB4gHyoBJQBLwPjAr1nH1Qb1ygXnh16OAj4BZwA+Bb4fLvw3cE349K1z3ZGBi+D2Jj3U9TrLuNwL/CbwYXh7SdQaeBL4Sfp0EpA7lOgP5wDZgWHj518CyoVZn4HxgHvBht7KI6wisAs4GDHgZuLS3MQTpimABsMXdt7p7K/Ar4MoYx3TK3L3G3deGX+8HNhH6D3QloQ8Owr8/G359JfArdz/k7tuALYTem0HFzAqAzwCPdisesnU2s9GEPjD+A8DdW929kSFc57AEYJiZJQDDgWqGWJ3d/Q1g91HFEdXRzHKB0e7+toeywv/XbZ8TClIiyAcqui1XhsuGDDMrBM4A3gFy3L0GQskCyA5vNlTeh/uBW4DObmVDuc6TgDrg8XBz2KNmNoIhXGd3rwLuBXYANcBed3+VIVznbiKtY3749dHlvRKkRNBTe9mQuXfWzEYCzwD/5O77jrdpD2WD6n0ws6XALndf09tdeigbVHUm9M14HvBzdz8DaCLUZHAsg77O4XbxKwk1geQBI8zs74+3Sw9lg6rOvXCsOp5S3YOUCCqBcd2WCwhdZg56ZpZIKAn80t2fDRfXhi8XCf/eFS4fCu/DucAVZlZOqIlvsZn9gqFd50qg0t3fCS8/TSgxDOU6LwG2uXudu7cBzwLnMLTrfFikdawMvz66vFeClAjeBaaa2UQzSwKuAX4b45hOWfjOgP8ANrn7j7ut+i1wffj19cAL3cqvMbNkM5sITCXUyTRouPtt7l7g7oWE/o6vufvfM7TrvBOoMLPp4aKLgI0M4ToTahI6y8yGh/+dX0SoD2wo1/mwiOoYbj7ab2Znhd+rL3bb58Ri3WMe5d75ywjdVVMGfCfW8fRRnT5F6BJwPbAu/HMZkAH8ESgN/07vts93wu9BCRHcWTAQf4BFfHzX0JCuM1AErA7/rZ8H0gJQ538FNgMfAv+b0N0yQ6rOwEpCfSBthL7Z/4+TqSNQHH6fyoAHCI8c0ZsfDTEhIhJwQWoaEhGRHigRiIgEnBKBiEjAKRGIiAScEoGISMApEUigmVmGma0L/+w0s6rw6wNm9mA/nfOfzOyLJ9jmV2Y2tT/OL3I03T4qEmZmdwIH3P3efjxHArCW0Iix7cfZ7gLg7939q/0Vi8hhuiIQ6YGZLeo2z8GdZvakmb1qZuVmdrWZ/TA89vt/hYf4ODwe/OtmtsbMXjk8RMBRFgNr3b3dzCab2dpu55xqZofHT3oTWBJOHCL9SolApHcmExr2+krgF8Cf3P004CDwmXAy+CnwOXefDzwG3NXDcc4F1gC4exmw18yKwuu+BDwRXtdJaIjh0/upPiJd9G1DpHdedvc2M/uA0CRH/xUu/wAoBKYDc4DfhyeGiic0bMDRcgmNl3PYo8CXzOxG4AscOX7+LkKjbvZ2lFWRk6JEINI7hyD0Td3M2vzjzrVOQv+PDNjg7mef4DgHgZRuy88A/wK8Bqxx94Zu61LC24v0KzUNifSNEiDLzM6G0NDgZja7h+02AVMOL7h7C/AK8HPg8aO2nQZs6J9wRT6mRCDSBzw0/enngHvM7H1Co8Ce08OmLxOacrK7XxIaQfbVwwVmlgMc9PAsVSL9SbePikSZmT0H3OLupeHlm4Ax7n57t22+Bexz9/+IUZgSIOojEIm+bxPqNC4NJ4XJhG4r7a6R0Pj7Iv1OVwQiIgGnPgIRkYBTIhARCTglAhGRgFMiEBEJOCUCEZGA+7+0PkMTQghxhAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "ax.plot(t, varpi_obs, label=\"JPL Horizons\")\n", + "ax.plot(tsim, varpiswifter, label=\"Swifter GR\")\n", + "ax.plot(tsim, varpiswiftest, label=\"Swiftest GR\")\n", + "ax.set_xlabel('Time (y)')\n", + "ax.set_ylabel('Mercury $\\\\varpi$ (deg)')\n", + "ax.legend()\n", + "print('Mean precession rate for Mercury long. peri. (arcsec/100 y)')\n", + "print(f'JPL Horizons : {np.mean(dvarpi_obs)}')\n", + "print(f'Swifter GR : {np.mean(dvarpi_swifter)}')\n", + "print(f'Swiftest GR : {np.mean(dvarpi_swiftest)}')\n", + "print(f'Obs - Swifter : {np.mean(dvarpi_obs - dvarpi_swifter)}')\n", + "print(f'Obs - Swiftest : {np.mean(dvarpi_obs - dvarpi_swiftest)}')\n", + "print(f'Swiftest - Swifter: {np.mean(dvarpi_swiftest - dvarpi_swifter)}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "swiftestOOF", + "language": "python", + "name": "swiftestoof" + }, + "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.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 2b8c47cb3..4fd797b76 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -1048,13 +1048,22 @@ def select_active_from_frame(ds, param, framenum=-1): frame = ds.isel(time=[framenum]) iframe = frame.isel(time=0) + if "name" in ds.dims: + count_dim = "name" + elif "id" in ds.dims: + count_dim = "id" + # Select only the active particles at this time step # Remove the inactive particles if param['OUT_FORM'] == 'XV' or param['OUT_FORM'] == 'XVEL': - iactive = iframe['id'].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['xhx'])), drop=True).id + iactive = iframe[count_dim].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['xhx'])), drop=True)[count_dim] else: - iactive = iframe['id'].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['a'])), drop = True).id - frame = frame.sel(id=iactive.values) + iactive = iframe[count_dim].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['a'])), drop = True)[count_dim] + if count_dim == "id": + frame = frame.sel(id=iactive.values) + elif count_dim == "name": + frame = frame.sel(name=iactive.values) + return frame @@ -1079,6 +1088,10 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram param_tmp['OUT_FORM'] = param['IN_FORM'] frame = select_active_from_frame(ds, param_tmp, framenum) + if "name" in frame.dims: + frame = frame.swap_dims({"name" : "id"}) + frame = frame.reset_coords("name") + if in_type == "NETCDF_DOUBLE" or in_type == "NETCDF_FLOAT": # Convert strings back to byte form and save the NetCDF file # Note: xarray will call the character array dimension string32. The Fortran code @@ -1186,8 +1199,13 @@ def swifter_xr2infile(ds, param, framenum=-1): ------- A set of input files for a Swifter run """ + frame = ds.isel(time=framenum) + if "name" in frame.dims: + frame = frame.swap_dims({"name" : "id"}) + frame = frame.reset_coords("name") + cb = frame.where(frame.id == 0, drop=True) pl = frame.where(frame.id > 0, drop=True) pl = pl.where(np.invert(np.isnan(pl['Gmass'])), drop=True).drop_vars(['j2rp2', 'j4rp4']) diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 7c17c3d32..c6fe5e8da 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -68,6 +68,8 @@ def __init__(self,read_param: bool = True, **kwargs: Any): integrator : {"symba","rmvs","whm","helio"}, default "symba" Name of the n-body integrator that will be used when executing a run. Parameter input file equivalent: None + read_param : bool, default False + Read the parameter file given by `param_file`. param_file : str, path-like, or file-lke, default "param.in" Name of the parameter input file that will be passed to the integrator. Parameter input file equivalent: None @@ -291,7 +293,7 @@ def __init__(self,read_param: bool = True, **kwargs: Any): # Set the location of the parameter input file param_file = kwargs.pop("param_file",self.param_file) - read_param = kwargs.pop("read_param",True) + read_param = kwargs.pop("read_param",False) self.set_parameter(verbose=False,param_file=param_file) #----------------------------------------------------------------- @@ -301,9 +303,8 @@ def __init__(self,read_param: bool = True, **kwargs: Any): # If the user asks to read in an old parameter file, override any default parameters with values from the file # If the file doesn't exist, flag it for now so we know to create it if read_param: - if os.path.exists(self.param_file): - self.read_param(self.param_file, codename=self.codename, verbose=self.verbose) - param_file_found = True + #good_param is self.read_param() + if self.read_param(): # We will add the parameter file to the kwarg list. This will keep the set_parameter method from # overriding everything with defaults when there are no arguments passed to Simulation() kwargs['param_file'] = self.param_file @@ -743,7 +744,7 @@ def get_parameter(self, **kwargs): return param_dict def set_integrator(self, - codename: Literal["swiftest", "swifter", "swift"] | None = None, + codename: Literal["Swiftest", "Swifter", "Swift"] | None = None, integrator: Literal["symba","rmvs","whm","helio"] | None = None, mtiny: float | None = None, gmtiny: float | None = None, @@ -2316,58 +2317,82 @@ def _combine_and_fix_dsnew(self,dsnew): Updated Dataset with ntp, npl values and types fixed. """ + if "id" not in self.data.dims: + if len(np.unique(dsnew['name'])) == len(dsnew['name']): + dsnew = dsnew.swap_dims({"id" : "name"}) + dsnew = dsnew.reset_coords("id") + else: + warnings.warn("Non-unique names detected for bodies. The Dataset will be dimensioned by integer id instead of name.") + warnings.warn("Consider using unique names instead.") + + if self.param['OUT_TYPE'] == "NETCDF_DOUBLE": + dsnew = io.fix_types(dsnew, ftype=np.float64) + elif self.param['OUT_TYPE'] == "NETCDF_FLOAT": + dsnew = io.fix_types(dsnew, ftype=np.float32) self.data = xr.combine_by_coords([self.data, dsnew]) def get_nvals(ds): + if "name" in ds.dims: + count_dim = "name" + elif "id" in ds.dims: + count_dim = "id" if "Gmass" in ds: - ds['ntp'] = ds['id'].where(np.isnan(ds['Gmass'])).count(dim="id") - ds['npl'] = ds['id'].where(~(np.isnan(ds['Gmass']))).count(dim="id") - 1 + ds['ntp'] = ds[count_dim].where(np.isnan(ds['Gmass'])).count(dim=count_dim) + ds['npl'] = ds[count_dim].where(~(np.isnan(ds['Gmass']))).count(dim=count_dim) - 1 else: - ds['ntp'] = ds['id'].count(dim="id") + ds['ntp'] = ds[count_dim].count(dim=count_dim) ds['npl'] = xr.full_like(ds['ntp'],0) return ds dsnew = get_nvals(dsnew) self.data = get_nvals(self.data) - if self.param['OUT_TYPE'] == "NETCDF_DOUBLE": - dsnew = io.fix_types(dsnew, ftype=np.float64) - self.data = io.fix_types(self.data, ftype=np.float64) - elif self.param['OUT_TYPE'] == "NETCDF_FLOAT": - dsnew = io.fix_types(dsnew, ftype=np.float32) - self.data = io.fix_types(self.data, ftype=np.float32) - return dsnew - def read_param(self, param_file, codename="Swiftest", verbose=True): + def read_param(self, + param_file : os.PathLike | str | None = None, + codename: Literal["Swiftest", "Swifter", "Swift"] | None = None, + verbose: bool | None = None): """ Reads in an input parameter file and stores the values in the param dictionary. Parameters ---------- - param_file : string + param_file : str or path-like, default is the value of the Simulation object's internal `param_file`. File name of the input parameter file - codename : string + codename : {"Swiftest", "Swifter", "Swift"}, default is the value of the Simulation object's internal`codename` Type of parameter file, either "Swift", "Swifter", or "Swiftest" + verbose : bool, default is the value of the Simulation object's internal `verbose` + If set to True, then more information is printed by Simulation methods as they are executed. Setting to + False suppresses most messages other than errors. Returns ------- - + True if the parameter file exists and is read correctly. False otherwise. """ + if param_file is None: + param_file = self.param_file + + if coename is None: + codename = self.codename + + if verbose is None: + verbose = self.verbose + + if not os.path.exists(param_file): + return False + if codename == "Swiftest": - param_old = self.param.copy() - self.param = io.read_swiftest_param(param_file, param_old, verbose=verbose) - self.codename = "Swiftest" + self.param = io.read_swiftest_param(param_file, param, verbose=verbose) elif codename == "Swifter": self.param = io.read_swifter_param(param_file, verbose=verbose) - self.codename = "Swifter" elif codename == "Swift": self.param = io.read_swift_param(param_file, verbose=verbose) - self.codename = "Swift" else: warnings.warn(f'{codename} is not a recognized code name. Valid options are "Swiftest", "Swifter", or "Swift".') - self.codename = "Unknown" - return + return False + + return True def write_param(self, codename: Literal["Swiftest", "Swifter", "Swift"] | None = None, From be6893487d307a92c876af8e3ed2ffd85d6289a9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 20:26:51 -0500 Subject: [PATCH 2/2] Added the GR/no GR comparison test Notebook to the whm_gr_test example --- .../whm_gr_test/swiftest_relativity.ipynb | 997 ++++++++++++++++-- 1 file changed, 927 insertions(+), 70 deletions(-) diff --git a/examples/whm_gr_test/swiftest_relativity.ipynb b/examples/whm_gr_test/swiftest_relativity.ipynb index ae586907c..0cd73753f 100644 --- a/examples/whm_gr_test/swiftest_relativity.ipynb +++ b/examples/whm_gr_test/swiftest_relativity.ipynb @@ -6,11 +6,11 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "import pandas as pd\n", "import swiftest\n", - "from astroquery.jplhorizons import Horizons" + "from astroquery.jplhorizons import Horizons\n", + "import datetime\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" ] }, { @@ -22,18 +22,442 @@ "name": "stdout", "output_type": "stream", "text": [ - "Reading Swifter file param.swifter.in\n", - "Reading in time 1.000e+03\n", - "Creating Dataset\n", - "Successfully converted 1001 output frames.\n", - "Swifter simulation data stored as xarray DataSet .ds\n" + "Creating the Sun as a central body\n", + "Fetching ephemerides data for Mercury from JPL/Horizons\n", + "Fetching ephemerides data for Venus from JPL/Horizons\n", + "Fetching ephemerides data for Earth from JPL/Horizons\n", + "Fetching ephemerides data for Mars from JPL/Horizons\n", + "Fetching ephemerides data for Jupiter from JPL/Horizons\n", + "Fetching ephemerides data for Saturn from JPL/Horizons\n", + "Fetching ephemerides data for Uranus from JPL/Horizons\n", + "Fetching ephemerides data for Neptune from JPL/Horizons\n", + "Writing initial conditions to file init_cond.nc\n", + "Writing parameter inputs to file /home/daminton/git_debug/swiftest/examples/whm_gr_test/param.gr.in\n" ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 9, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Uranus' 'Neptune'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables: (12/14)\n",
+       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
+       "    id             (name) int64 0 1 2 3 4 5 6 7 8\n",
+       "    Gmass          (time, name) float64 39.48 6.554e-06 ... 0.001724 0.002034\n",
+       "    radius         (time, name) float64 0.00465 1.631e-05 ... 0.0001646\n",
+       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
+       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
+       "    ...             ...\n",
+       "    xhz            (time, name) float64 nan -0.002529 ... -0.0407 -0.7287\n",
+       "    vhx            (time, name) float64 nan -9.241 3.454 ... -1.315 -0.08755\n",
+       "    vhy            (time, name) float64 nan 7.755 6.477 ... 1.932 0.5372 1.149\n",
+       "    vhz            (time, name) float64 nan 1.481 -0.1103 ... 0.01905 -0.02168\n",
+       "    ntp            (time) int64 0\n",
+       "    npl            (time) int64 8
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 9, time: 1)\n", + "Coordinates:\n", + " * name (name) \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 9, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Uranus' 'Neptune'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables: (12/14)\n",
+       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
+       "    id             (name) int64 0 1 2 3 4 5 6 7 8\n",
+       "    Gmass          (time, name) float64 39.48 6.554e-06 ... 0.001724 0.002034\n",
+       "    radius         (time, name) float64 0.00465 1.631e-05 ... 0.0001646\n",
+       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
+       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
+       "    ...             ...\n",
+       "    xhz            (time, name) float64 nan -0.002529 ... -0.0407 -0.7287\n",
+       "    vhx            (time, name) float64 nan -9.241 3.454 ... -1.315 -0.08755\n",
+       "    vhy            (time, name) float64 nan 7.755 6.477 ... 1.932 0.5372 1.149\n",
+       "    vhz            (time, name) float64 nan 1.481 -0.1103 ... 0.01905 -0.02168\n",
+       "    ntp            (time) int64 0\n",
+       "    npl            (time) int64 8
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 9, time: 1)\n", + "Coordinates:\n", + " * name (name) " - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "\n", "ax.plot(t, varpi_obs, label=\"JPL Horizons\")\n", - "ax.plot(tsim, varpiswifter, label=\"Swifter GR\")\n", - "ax.plot(tsim, varpiswiftest, label=\"Swiftest GR\")\n", + "ax.plot(tsim, varpisim_gr, label=\"Swiftest WHM GR\")\n", + "ax.plot(tsim, varpisim_nogr, label=\"Swiftest WHM No GR\")\n", "ax.set_xlabel('Time (y)')\n", "ax.set_ylabel('Mercury $\\\\varpi$ (deg)')\n", "ax.legend()\n", "print('Mean precession rate for Mercury long. peri. (arcsec/100 y)')\n", - "print(f'JPL Horizons : {np.mean(dvarpi_obs)}')\n", - "print(f'Swifter GR : {np.mean(dvarpi_swifter)}')\n", - "print(f'Swiftest GR : {np.mean(dvarpi_swiftest)}')\n", - "print(f'Obs - Swifter : {np.mean(dvarpi_obs - dvarpi_swifter)}')\n", - "print(f'Obs - Swiftest : {np.mean(dvarpi_obs - dvarpi_swiftest)}')\n", - "print(f'Swiftest - Swifter: {np.mean(dvarpi_swiftest - dvarpi_swifter)}')" + "print(f'JPL Horizons : {np.mean(dvarpi_obs)}')\n", + "print(f'Swiftest No GR : {np.mean(dvarpi_nogr)}')\n", + "print(f'Swiftest GR : {np.mean(dvarpi_gr)}')\n", + "print(f'Obs - Swiftest GR : {np.mean(dvarpi_obs - dvarpi_gr)}')\n", + "print(f'Obs - Swiftest No GR : {np.mean(dvarpi_obs - dvarpi_nogr)}')" ] }, { @@ -171,9 +1028,9 @@ ], "metadata": { "kernelspec": { - "display_name": "swiftestOOF", + "display_name": "swiftest", "language": "python", - "name": "swiftestoof" + "name": "swiftest" }, "language_info": { "codemirror_mode": { @@ -185,7 +1042,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.10" + "version": "3.8.5" } }, "nbformat": 4,