diff --git a/examples/helio_swifter_comparison/init_cond.py b/examples/helio_swifter_comparison/init_cond.py index 7e45bb4bb..4680d9e0a 100644 --- a/examples/helio_swifter_comparison/init_cond.py +++ b/examples/helio_swifter_comparison/init_cond.py @@ -1,9 +1,5 @@ +#!/usr/bin/env python3 import swiftest -import numpy as np -import sys -from astroquery.jplhorizons import Horizons -import astropy.constants as const -from scipy.io import FortranFile sim = swiftest.Simulation() diff --git a/examples/whm_gr_test/cb.swiftest.in b/examples/whm_gr_test/cb.swiftest.in index 46a8d0257..058975b81 100644 --- a/examples/whm_gr_test/cb.swiftest.in +++ b/examples/whm_gr_test/cb.swiftest.in @@ -1,4 +1,4 @@ -1.0 +39.476926408897626 0.004650467260962157 4.7535806948127355e-12 -2.2473967953572827e-18 diff --git a/examples/whm_gr_test/init_cond.py b/examples/whm_gr_test/init_cond.py old mode 100644 new mode 100755 index 7904eb100..c2580cddb --- a/examples/whm_gr_test/init_cond.py +++ b/examples/whm_gr_test/init_cond.py @@ -1,226 +1,51 @@ -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' +#!/usr/bin/env python3 +import swiftest + +sim = swiftest.Simulation() + +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_FORM'] = "EL" +sim.param['OUT_STAT'] = "UNKNOWN" +sim.param['GR'] = 'YES' + +bodyid = { + "Sun": 0, + "Mercury": 1, + "Venus": 2, + "Earth": 3, + "Mars": 4, + "Jupiter": 5, + "Saturn": 6, + "Uranus": 7, + "Neptune": 8, } -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}') +for name, id in bodyid.items(): + sim.add(name, idval=id) + +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.dat" +sim.param['ENC_OUT'] = "enc.swiftest.dat" +sim.save("param.swiftest.in") +sim.param['PL_IN'] = "pl.swifter.in" +sim.param['TP_IN'] = "tp.swifter.in" +sim.param['BIN_OUT'] = "bin.swifter.dat" +sim.param['ENC_OUT'] = "enc.swifter.dat" +sim.save("param.swifter.in", codename="Swifter") - sys.stdout = sys.__stdout__ diff --git a/examples/whm_gr_test/param.swifter.in b/examples/whm_gr_test/param.swifter.in index 6dbf5ae15..789250f41 100644 --- a/examples/whm_gr_test/param.swifter.in +++ b/examples/whm_gr_test/param.swifter.in @@ -1,27 +1,27 @@ -! Swifter input file generated using init_cond.py -T0 0 -TSTOP 1000.0 -DT 0.0006844626967830253 -PL_IN pl.swifter.in -TP_IN tp.swifter.in -IN_TYPE ASCII -ISTEP_OUT 1461 -ISTEP_DUMP 1461 -BIN_OUT bin.swifter.dat -OUT_TYPE REAL8 -OUT_FORM EL -OUT_STAT UNKNOWN -J2 4.7535806948127355e-12 -J4 -2.2473967953572827e-18 -CHK_CLOSE yes -CHK_RMIN 0.004650467260962157 -CHK_RMAX 1000.0 -CHK_EJECT 1000.0 -CHK_QMIN 0.004650467260962157 -CHK_QMIN_COORD HELIO -CHK_QMIN_RANGE 0.004650467260962157 1000.0 -ENC_OUT enc.swifter.dat -EXTRA_FORCE no -BIG_DISCARD no -RHILL_PRESENT yes -C 63241.07708426628 +! VERSION Swifter parameter file converted from Swiftest +T0 0.0 +TSTOP 1000.0 +DT 0.0006844626967830253 +ISTEP_OUT 1461 +ISTEP_DUMP 1461 +OUT_FORM EL +OUT_TYPE REAL8 +OUT_STAT UNKNOWN +IN_TYPE ASCII +PL_IN pl.swifter.in +TP_IN tp.swifter.in +BIN_OUT bin.swifter.dat +ENC_OUT enc.swifter.dat +CHK_QMIN 0.004650467260962157 +CHK_RMIN 0.004650467260962157 +CHK_RMAX 1000.0 +CHK_EJECT 1000.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 1000.0 +EXTRA_FORCE NO +BIG_DISCARD NO +CHK_CLOSE YES +C 63241.07708426628 +J2 4.7535806948127355e-12 +J4 -2.2473967953572827e-18 +RHILL_PRESENT YES diff --git a/examples/whm_gr_test/param.swiftest.in b/examples/whm_gr_test/param.swiftest.in index 6e7c9ff28..ace6f3cad 100644 --- a/examples/whm_gr_test/param.swiftest.in +++ b/examples/whm_gr_test/param.swiftest.in @@ -1,30 +1,35 @@ -! Swiftest input file generated using init_cond.py -T0 0 -TSTOP 1000.0 -DT 0.0006844626967830253 -CB_IN cb.swiftest.in -PL_IN pl.swiftest.in -TP_IN tp.swiftest.in -IN_TYPE ASCII -ISTEP_OUT 1461 -ISTEP_DUMP 1461 -BIN_OUT bin.swiftest.dat -OUT_TYPE REAL8 -OUT_FORM EL -OUT_STAT UNKNOWN -CHK_CLOSE yes -CHK_RMIN 0.004650467260962157 -CHK_RMAX 1000.0 -CHK_EJECT 1000.0 -CHK_QMIN 0.004650467260962157 -CHK_QMIN_COORD HELIO -CHK_QMIN_RANGE 0.004650467260962157 1000.0 -ENC_OUT enc.swiftest.dat -EXTRA_FORCE no -BIG_DISCARD no -ROTATION no -TIDES no -GR yes -MU2KG 1.988409870698051e+30 -DU2M 149597870700.0 -TU2S 31557600.0 +! VERSION Swiftest parameter input +T0 0.0 +TSTOP 1000.0 +DT 0.0006844626967830253 +ISTEP_OUT 1461 +ISTEP_DUMP 1461 +OUT_FORM EL +OUT_TYPE REAL8 +OUT_STAT UNKNOWN +IN_TYPE ASCII +PL_IN pl.swiftest.in +TP_IN tp.swiftest.in +CB_IN cb.swiftest.in +BIN_OUT bin.swiftest.dat +ENC_OUT enc.swiftest.dat +CHK_QMIN 0.004650467260962157 +CHK_RMIN 0.004650467260962157 +CHK_RMAX 1000.0 +CHK_EJECT 1000.0 +CHK_QMIN_COORD HELIO +CHK_QMIN_RANGE 0.004650467260962157 1000.0 +MU2KG 1.988409870698051e+30 +TU2S 31557600.0 +DU2M 149597870700.0 +EXTRA_FORCE NO +BIG_DISCARD NO +CHK_CLOSE YES +FRAGMENTATION NO +ROTATION NO +TIDES NO +ENERGY NO +GR YES +YARKOVSKY NO +YORP NO +MTINY 0.0 diff --git a/examples/whm_gr_test/pl.swifter.in b/examples/whm_gr_test/pl.swifter.in index 60d090453..fb614bcbd 100644 --- a/examples/whm_gr_test/pl.swifter.in +++ b/examples/whm_gr_test/pl.swifter.in @@ -1,40 +1,36 @@ -10 ! Planet input file generated using init_cond.py using JPL Horizons data for the major planets (and Pluto) for epoch 2021-01-28 -1 39.476926408897626 +9 +0 39.476926408897625196 0.0 0.0 0.0 0.0 0.0 0.0 -2 6.553709809565314e-06 0.001475122968086379 -1.6306381826061646e-05 -0.1030256860922895 0.2897796047098886 0.01422904600374035 --11.74004209950937 3.8343124110162736 1.3902496665973592 -3 9.663313399581537e-05 0.006759127649782299 -4.0453784346544176e-05 -0.06110218027254217 -0.7245466901305982 -0.01346904300924688 -7.311995449678243 0.5941125721336201 -0.4137913843379075 -4 0.00012002693582795245 0.010044756567546644 -4.25875607065041e-05 --0.6061796342297583 0.7761214554702035 -3.4750047790977e-05 --5.054824314301841 -3.891667468503358 0.00019720338148272726 -5 1.2739802010675942e-05 0.0072464490746299084 -2.2657408050928896e-05 -0.2751944175855944 1.51937688993241 0.02508924593104206 --4.835983593209577 1.344855094041679 0.14681413000004515 -6 0.037692251088985676 0.3552852357486061 -0.0004673261703049093 -3.200135438345358 -3.953498213518368 -0.05517737289975112 -2.111393749129838 1.8660266890185446 -0.05498941067210089 -7 0.011285899820091273 0.4376306456694341 -0.00038925687730393614 -5.607382165725712 -8.258649105608766 -0.07958445228024298 -1.5748468603228847 1.1414574661825514 -0.08250331331320372 -8 0.001723658947826773 0.4690969274244374 -0.00016953449859497232 -15.28225422201768 12.53905314208462 -0.1514143582550325 --0.9198472198098231 1.0454390993472462 0.01574538863031621 -9 0.0020336100526728304 0.7807192056765467 -0.00016458790412449367 -29.47483071169769 -5.147686530859088 -0.5733441819169969 -0.19191677740340274 1.1385110364087574 -0.027844325148353527 -10 2.924216771029454e-07 0.0538346817277698 -7.943294877391593e-06 -14.14000920780611 -31.14141812522779 -0.7565722591093476 -1.073396108697069 0.23003123192799815 -0.33424529561177047 +1 6.5537098095653139645e-06 0.0014751244996981091688 +1.6306381826061645943e-05 +0.3424231950105547928 0.045418598038753463242 -0.027698852186471431547 +-3.3367255059737357791 10.64622790896598922 1.1760542104884461008 +2 9.663313399581537916e-05 0.006759105927487850036 +4.0453784346544178454e-05 +-0.7187409324847692238 0.008460121607685697903 0.041591140610292232083 +-0.12984901004910088259 -7.4200262222437848615 -0.09433956287915627552 +3 0.000120026935827952453094 0.010044781409779763954 +4.25875607065040958e-05 +0.34089134446377761245 -0.95773246142964407746 4.435598231277717939e-05 +5.817808761117269037 2.0836935934869021533 -0.00011769228482089048983 +4 1.2739802010675941456e-05 0.007246745952808948377 +2.265740805092889601e-05 +-1.5181949547849429294 0.68396861668094566244 0.051574975426518870902 +-1.908390450061588966 -4.22386803675231535 -0.041705300362175825686 +5 0.037692251088985676735 0.35527129445679039879 +0.00046732617030490929307 +4.0455430243320495975 -2.9975165850765859155 -0.07806209765007877943 +1.6095092959604531543 2.3468202746428659788 -0.045756236082127820444 +6 0.011285899820091272997 0.43765252895139074356 +0.00038925687730393611812 +6.29491448542079457 -7.7099360898166890976 -0.11647820911165690516 +1.466832763395646828 1.2863740305106504463 -0.0807227702193550598 +7 0.0017236589478267730203 0.4695335374930126262 +0.00016953449859497231466 +14.8586976850394894 13.004806892491039605 -0.14422203290693380584 +-0.9552190982402187873 1.0163552941629989916 0.016097455433874226457 +8 0.0020336100526728302319 0.78128379442879379807 +0.000164587904124493665 +29.55697963607957135 -4.6325049341728004038 -0.5858344271811812831 +0.17174414812321077623 1.1422632832628163399 -0.02744788192962509712 diff --git a/examples/whm_gr_test/pl.swiftest.in b/examples/whm_gr_test/pl.swiftest.in index 032262e70..874012fc6 100644 --- a/examples/whm_gr_test/pl.swiftest.in +++ b/examples/whm_gr_test/pl.swiftest.in @@ -1,37 +1,33 @@ -9 -2 1.6601367952719304e-07 -1.6306381826061646e-05 -0.1030256860922895 0.2897796047098886 0.01422904600374035 --11.74004209950937 3.8343124110162736 1.3902496665973592 -3 2.4478383396645447e-06 -4.0453784346544176e-05 -0.06110218027254217 -0.7245466901305982 -0.01346904300924688 -7.311995449678243 0.5941125721336201 -0.4137913843379075 -4 3.0404326462685257e-06 -4.25875607065041e-05 --0.6061796342297583 0.7761214554702035 -3.4750047790977e-05 --5.054824314301841 -3.891667468503358 0.00019720338148272726 -5 3.2271514450538743e-07 -2.2657408050928896e-05 -0.2751944175855944 1.51937688993241 0.02508924593104206 --4.835983593209577 1.344855094041679 0.14681413000004515 -6 0.0009547919384243222 -0.0004673261703049093 -3.200135438345358 -3.953498213518368 -0.05517737289975112 -2.111393749129838 1.8660266890185446 -0.05498941067210089 -7 0.0002858859806661029 -0.00038925687730393614 -5.607382165725712 -8.258649105608766 -0.07958445228024298 -1.5748468603228847 1.1414574661825514 -0.08250331331320372 -8 4.3662440433515637e-05 -0.00016953449859497232 -15.28225422201768 12.53905314208462 -0.1514143582550325 --0.9198472198098231 1.0454390993472462 0.01574538863031621 -9 5.151389020535497e-05 -0.00016458790412449367 -29.47483071169769 -5.147686530859088 -0.5733441819169969 -0.19191677740340274 1.1385110364087574 -0.027844325148353527 -10 7.407407407407407e-09 -7.943294877391593e-06 -14.14000920780611 -31.14141812522779 -0.7565722591093476 -1.073396108697069 0.23003123192799815 -0.33424529561177047 +8 +1 6.5537098095653139645e-06 +1.6306381826061645943e-05 +0.3424231950105547928 0.045418598038753463242 -0.027698852186471431547 +-3.3367255059737357791 10.64622790896598922 1.1760542104884461008 +2 9.663313399581537916e-05 +4.0453784346544178454e-05 +-0.7187409324847692238 0.008460121607685697903 0.041591140610292232083 +-0.12984901004910088259 -7.4200262222437848615 -0.09433956287915627552 +3 0.000120026935827952453094 +4.25875607065040958e-05 +0.34089134446377761245 -0.95773246142964407746 4.435598231277717939e-05 +5.817808761117269037 2.0836935934869021533 -0.00011769228482089048983 +4 1.2739802010675941456e-05 +2.265740805092889601e-05 +-1.5181949547849429294 0.68396861668094566244 0.051574975426518870902 +-1.908390450061588966 -4.22386803675231535 -0.041705300362175825686 +5 0.037692251088985676735 +0.00046732617030490929307 +4.0455430243320495975 -2.9975165850765859155 -0.07806209765007877943 +1.6095092959604531543 2.3468202746428659788 -0.045756236082127820444 +6 0.011285899820091272997 +0.00038925687730393611812 +6.29491448542079457 -7.7099360898166890976 -0.11647820911165690516 +1.466832763395646828 1.2863740305106504463 -0.0807227702193550598 +7 0.0017236589478267730203 +0.00016953449859497231466 +14.8586976850394894 13.004806892491039605 -0.14422203290693380584 +-0.9552190982402187873 1.0163552941629989916 0.016097455433874226457 +8 0.0020336100526728302319 +0.000164587904124493665 +29.55697963607957135 -4.6325049341728004038 -0.5858344271811812831 +0.17174414812321077623 1.1422632832628163399 -0.02744788192962509712 diff --git a/examples/whm_gr_test/swiftest_relativity.ipynb b/examples/whm_gr_test/swiftest_relativity.ipynb index ae586907c..fe0e42cc9 100644 --- a/examples/whm_gr_test/swiftest_relativity.ipynb +++ b/examples/whm_gr_test/swiftest_relativity.ipynb @@ -89,8 +89,8 @@ "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", + "varpiswiftest = swiftestdat['varpi'].sel(id=1) * 180.0 / np.pi\n", + "varpiswifter = swifterdat['varpi'].sel(id=1) * 180.0 / np.pi\n", "tsim = swiftestdat['time']" ] }, @@ -116,16 +116,16 @@ "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" + "Swifter GR : 571.6748624042897\n", + "Swiftest GR : 571.6748624128945\n", + "Obs - Swifter : -0.3538117742853119\n", + "Obs - Swiftest : -0.35381178289026777\n", + "Swiftest - Swifter: 8.604956747149117e-09\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", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA40klEQVR4nO3deZxN5R/A8c/3zj72rUWytSjLGIwt2WXPloqUpYRSaUeL5PdTKhVtZFfJElkiS0mhnyyDQgihhiIiDGOW+/39cS+NMcO9zMyZ5ft+ve7LPc89zznf52K+c57nnOcRVcUYY4y5GJfTARhjjMkeLGEYY4zxiSUMY4wxPrGEYYwxxieWMIwxxvgk0OkAMlLRokW1dOnSTodhjDHZRnR09CFVLZbaZzk6YZQuXZp169Y5HYYxxmQbIrI3rc+sS8oYY4xPLGEYY4zxiSUMY4wxPsnRYxipSUhIICYmhri4OKdDyXVCQ0MpUaIEQUFBTodijLkEuS5hxMTEkC9fPkqXLo2IOB1OrqGqHD58mJiYGMqUKeN0OMaYS5DruqTi4uIoUqSIJYtMJiIUKVLEruyMycZyXcIALFk4xL53Y7K3XJkwjDEmp/rs6/d45ZP7SUpKSPdjW8JwQN68edmzZw9hYWFERkZSvnx5+vTpg9vtZs+ePVSsWPGC9QcPHszw4cPPKStdujSHDh3yK46WLVty9OhRf8M3xmRBcaeO0HtCLYbs+5Bv49Zw7MTRdD9Hrhv0zkquu+46Nm7cSGJiIo0aNWLOnDlUrVo1w8+rqqgqX375ZYafyxiT8WZ89TZz93zMT8EJ1DhemG51BlCoQKqze1wWu8LIAgIDA7nlllvYuXNnuhzvrbfeomLFilSsWJERI0YAsGfPHm6++WYefvhhqlatyu+//372qmT06NFERkYSGRlJmTJlaNiwIQBTp06lUqVKVKxYkf79+589ft68eXn++eepXLkytWrV4sCBAwB89tlnVKxYkcqVK1OvXr10aYsxJm0790bTd1xd/rN/Aj8FJ9A0viQf9PmGelVaZMj5cvUVxstfbOHn/cfS9Zjli+fnpdsr+FXn5MmTLF26lCFDhvhc5+233+aTTz45u71//34AoqOjmThxIqtXr0ZVqVmzJvXr16dQoUJs376diRMn8sEHH5xzrD59+tCnTx8SEhJo1KgRTz75JPv376d///5ER0dTqFAhmjZtypw5c2jXrh2xsbHUqlWLoUOH8uyzzzJ27FheeOEFhgwZwuLFi7nmmmusq8uYDOROSuT7TUt5Zd1T7AuEqKNX8uhtz1DlpqaIK+OuA+wKw0G7du0iMjKSOnXq0KpVK1q08P23gieeeIKNGzeefRUvXhyAlStX0r59e/LkyUPevHnp0KEDK1asAKBUqVLUqlUrzWP269ePRo0acfvtt7N27VoaNGhAsWLFCAwMpEuXLixfvhyA4OBgWrduDUC1atXYs2cPAHXq1KF79+6MHTuWpKSkS/lKjDEXsf6X9XSbUIuHf3yamCChV772vNJ9PlXLN8/QZAG5/ArD3yuB9HZmDCM9qWqan+XJkyfNzyZNmsTevXt57733LnqcoKCgs7fIBgQEkJiYCMDo0aNZvXo1CxYsIDIyko0bN1KkSJFLaYYxJoWkxHg+/3Y0Q/aNhWCoGVuIxmXr07n5fzItBrvCyGHq1avHnDlzOHnyJLGxscyePZu6detesE50dDTDhw/nk08+weX9DaVmzZp89913HDp0iKSkJKZOnUr9+vUveJxdu3ZRs2ZNhgwZQtGiRfn999/TrV3G5Faqyu69K2n4URWG7BtLmNvNXa47+E+XxZmaLCATrzBEpBwwPVlRWWAQsAwYDeQF9gBdVPW8gQURaQ6MBAKAcao6LKNjzgiJiYmEhIRccJ/t27dTokSJs9tvv/02d955p0/Hr1q1Kt27d6dGjRoA9OzZkypVqpztNkrNe++9x99//312sDsqKopx48bx6quv0rBhQ1SVli1b0rZt2wue+5lnnmHHjh2oKo0bN6Zy5co+xWyMSduTY2/n65C9BLqESrEhdCx3Lx0aP+5ILHKhrocMO6lIALAPqAnMBJ5W1e9E5H6gjKq+mMr+vwC3ATHAWqCzqv58ofNERUVpygWUtm7dys0335xubfHXjz/+yIMPPsiaNWsci8FJTn//xmQX7826nzlH13Ag0NP9e7erBs/fOy7DZ0wQkWhVjUrtM6fGMBoDu1R1r/fKY7m3/CtgMfBiiv1rADtV9VcAEZkGtAUumDCymtGjR/POO++cvdXVGGNS2n9wN8PnPMpXIXshUKh2Ih+v3TOXKwul/3MV/nIqYXQCpnrfbwbaAHOBO4FrU9n/GiB5h3gMnquT84hIL6AXQMmSJdMp3PRx5vZVY4xJKSExgfFfPMe8wwv5PURwqfLmTS9Qp+qdhAUHOB0e4MCgt4gE40kQn3mL7gf6ikg0kA+IT61aKmWp9qWp6hhVjVLVqGLFnM/IxhhzMT/+PJv+E9vw/rFF/B4kNDlVjqm3vEGTWp2yTLIAZ64wWgDrVfUAgKpuA5oCiMiNQKtU6sRw7pVHCWB/BsdpjDEZKjHJzY9b59Jj3YtosFAqHu4u2ZfWdbpTKG+o0+Gdx4mE0Zl/u6MQkStU9aCIuIAX8NwxldJa4AYRKYNnsLwTcE9mBGuMMRkh7tQRBn7cheVBv+ESqB97E0+0GUzpEs4+H3YhmZowRCQcz51OvZMVdxaRvt73nwMTvfsWx3P7bEtVTRSRR/AMiAcAE1R1SyaGbowx6WbCFw/z9t8rIARujg+kXpFOPNy1Py5X1l4zJlPHMFT1pKoWUdV/kpWNVNUbva8B6r3PV1X3q2rLZPt96d3nOlUdmplxp7ehQ4dSoUIFIiIiiIyMZPXq1T7VGzRoEF9//TUAK1asoEKFCkRGRrJq1ap0mXn2wIED3HPPPZQtW5Zq1apRu3ZtZs+eDcC3335LgQIFqFKlCjfddBNPP/30ZZ/PmNwm9uRxnh3bkw8OeW4Mbe9uyEfd1vJIhwFZPllALp8axAmrVq1i/vz5rF+/npCQEA4dOkR8fGrj/OdLPjnhlClTePrpp+nRoweTJk1i3bp1tGzZ8gK1z5WYmEhg4L9//apKu3bt6NatG59++ikAe/fuZd68eWf3qVu3LvPnz+fUqVNUqVKF9u3bU6dOHZ/PaUxu5XYrz33Ukm/1N2KDXRRPgGcrvUbjKN//z2YFljAy2R9//EHRokXPPu1dtGhRANasWcOwYcP4/PPPmTt3Lp06deKff/7B7XZTvnx5fv31V7p3707r1q05evQoM2bMYPHixSxZsoTvv/+eU6dOsXLlSgYOHEjr1q159NFH2bRpE4mJiQwePJi2bdsyadIkFixYQFxcHLGxsXzzzTdn4/rmm28IDg4+57bfUqVK8eijj57XhjMLP+3bty+Dvy1jsr9V66cya/1UFgfF4ELoHNCCgfcNy/CJAjNC7k4YCwfAn5vS95hXVYIWac9a0rRpU4YMGcKNN95IkyZNuPvuu6lfvz5Vq1Zlw4YNgKe7qWLFiqxdu5bExERq1jz3kZOePXuycuVKWrduTceOHc9eYZyZOPC5556jUaNGTJgwgaNHj1KjRg2aNGkCeK5wfvrpJwoXLnzOMbds2eLz4k1Hjhxhx44dtuaFMRew7+/DfLroZT5KWAZBcNPpAMZ1WUmBPHmdDu2SZb8Ul83lzZuX6OhoxowZQ7Fixbj77ruZNGkSgYGBXH/99WzdupU1a9bw5JNPsnz5clasWHHRyQNTWrJkCcOGDSMyMpIGDRoQFxfHb7/9BsBtt912XrJITd++falcuTLVq1c/W7ZixQoiIiK46qqraN26NVdddZV/jTcmF3AnJbJ951f0n9GMjxKWEeZ283ihdkzplr2TBeT2K4wLXAlkpICAABo0aECDBg2oVKkSkydPpnv37tStW5eFCxcSFBREkyZN6N69O0lJSeet330xqsqsWbMoV67cOeWrV69Oc4rzChUqMGvWrLPb77//PocOHSIq6t8pZc6MYfzyyy/ceuuttG/fnsjISL9iMyYnU7eb/p80YxEHIQwax5WiZ/P3qFiqtNOhpQu7wshk27dvZ8eOHWe3N27cSKlSpQDP1OQjRoygdu3aFCtWjMOHD7Nt2zYqVLjwfdn58uXj+PHjZ7ebNWvGu+++e3ZNizNdXRfSqFEj4uLiGDVq1NmykydPprrvjTfeyMCBA3nttdcuelxjcouV0bPo8mF1FnGQ4gnKYwVaMaL3/ByTLMASRqY7ceIE3bp1o3z58kRERPDzzz8zePBgwLMGxYEDB86ODURERBAREXHR2SkbNmzIzz//TGRkJNOnT+fFF18kISGBiIgIKlasyIsvppzL8Xwiwpw5c/juu+8oU6YMNWrUoFu3bmkmhT59+rB8+XJ2797t3xdgTA5z9PghBkxqwZM/DWJ3SBytk67nsy7reLBdtlyB4YIcmd48s2TF6c1zO/v+TU7y9aoxvL/pPXaGKBVjg+lZYyiNo5o7HdZlyYrTmxtjTLZ1Ii6Odz7vwtSEXyAE7g2qQa9uoyiUJ9jp0DKUJQxjjPGRut0s/t9ohm1/n8OBLirFhtCtcm+a1XnQ6dAyhSUMY4zxQULCSR6a3JDVQScJd0HngIa0bv0iESVzzzIKljCMMeYipix6k89+m8yuEKVO7JW0ibiXlrd2dzqsTGcJwxhj0rD91xUMWvooPwcnERyk3Ha6LIO6z6RgeM4eq0iLJQxjjEkhKSmB92Y9x9zjX3IqAKr9U5THmg2jarnaTofmKHsOwwFZaXrzo0eP8sEHH6T5uU15bnKbiV+8QKfx1Rl3ahF/BbroUqABEx75NtcnC7CEkemST2/+008/8fXXX3PttddevCKe6c3PTCJ4ZnrzjRs3sn379gxJGGemPK9Xrx6//vor0dHRTJs2jZiYmLP71K1blw0bNrBhwwbmz5/P999/f0lxGOO0TVvnMWLGI7z191y2hSRR63gxplQfTN/272aLtSoyQ6YlDBEpJyIbk72OicjjIhIpIj94y9aJSI006j8hIltEZLOITBWRrLfgrQ9Sm968ePHirFmzhg4dOgAwd+5cwsLCiI+PJy4ujrJlywLQvXt3Zs6cybhx45gxYwZDhgyhc+fODBo0iOnTp5990js2Npb777+f6tWrU6VKFebOnQt4ZqStUaMGkZGRREREsGPHDgYMGMCuXbuIjIzkmWeeOSdWm/Lc5AYJCYnMXzGBB38YyPhT3wHwXKnX+G/XL4kof0e2nIY8o2TaGIaqbgciAUQkAM/a3LOBscDLqrpQRFoCrwMNktcVkWuAx4DyqnpKRGbgWdd70uXE9Nqa19j297bLOcR5bip8E/1r9E/zcyenNx89ejT9+vWjS5cuxMfHk5SUxLBhw9i8eTMbN248L1ab8tzkdIcObePteS8zL2AzuFzUii1MmwoduL1+9lrYKLM4NejdGNilqntFRIH83vICwP406gQCYSKSAIRfYL8s7cz05itWrGDZsmXcfffdDBs2jO7du6c6vXlSUtIlTW8+b968s7PcnpnevHbt2gwdOpSYmBg6dOjADTfc4Ndx+/bty8qVKwkODmbt2rXAv1Oeb9++nQEDBtiU5ybbGP/FQ4z4eyUEeLa7BtelX5f3CA60K4q0OJUwOgFTve8fBxaLyHA8XWS3pNxZVfd5P/8NOAUsUdUlqR1YRHoBvQBKlix5wSAudCWQkZya3vzmm2+mZs2aLFiwgGbNmjFu3Liz3V2psSnPTU508tRRnp58ByvCDhLmdlMlrhiPNXuXCmUrOR1alpfpqVREgoE2wGfeooeAJ1T1WuAJYHwqdQoBbYEyQHEgj4jcm9rxVXWMqkapalSxYlnvCUwnpzf/9ddfKVu2LI899hht2rThp59+Oq9ucjbluclpPvi8Fx0+vZUVYQfJl+Tmxeue4MOHvrVk4SMnrr1aAOtV9YB3uxvwuff9Z0Bqg95NgN2q+peqJnj3P+9KJDtwcnrz6dOnU7FiRSIjI9m2bRtdu3alSJEi1KlTh4oVK5436G1TnpucYtuv0Tw5pi2jjq9iX6BQ59gNLO+2idvr9XQ6tGwl06c3F5FpwGJVnejd3go8pKrfikhj4HVVrZaiTk1gAlAdT5fUJGCdqr57oXPZ9OZZj33/JjMlJLn5/NtxfLhnJH8Furg6QRnWcC43X1uasOAAp8PLkrLM9OYiEg7cBvROVvwgMFJEAoE4vOMPIlIcGKeqLVV1tYjMBNYDicAGYExmxm6MyT7U7WbZ6neZu+lLvgnZT1CA0CmwA/06PknePAWcDi/bytSEoaongSIpylYC1VLZdz/QMtn2S8BLGR2jMSZ7S3Ir0xcP4NWDCyEEbo4LoP0Ng7ijXju7A+oy5cq5pFT1ouMCJv3l5NUdjfMSE+L4J/YQz0+5n1V59pPXrdwR3pLeHZ8nX56CToeXI+S6hBEaGsrhw4cpUqSIJY1MpKocPnyY0NBs+YC+yQYe/bgRKwOOQ16IjAulR7XBNIpq7XRYOUquSxglSpQgJiaGv/76y+lQcp3Q0FBKlCjhdBgmh1m8cgSzNs9iVZjn9vD7Q1vzRLdXHY4qZ8p1CSMoKIgyZco4HYYx5jIdPHKAEXN78YX8CmFQ7XQo73b5yrqfMlCuSxjGmOwtKTGepavfZeKWKWwOSaD0aXi40n9pUbut06HleJYwjDHZRlx8PM9+0oRlAUcgBNq7q/FI59Fckc/GxjKDJQxjTLawY+86nl/Uk62hSVQ8FUybMnfQudlzToeVq1jCMMZkaYcO7eK/8x5gpRwiIFi5y1WDR+8bTcE8IU6HlutYwjDGZFlfrhjFh9s+4PcgpeLx/LSp+BAdG97ndFi5liUMY0yWs/ePnYxd8iRz2Q3BcH/YrTzebZQ9O+UwSxjGmCwjMSGOyQuf472/l5AoQsTxMJ5pNobIGyKdDs1gCcMYk0X8fSKWZ6Y3Yk3gSRDhwTxtaNLkGcpfU9Dp0IyXJQxjjKNOnjjIjG8/5MOD0zgR6OKW2GJ0qd6NetW7OR2aScEShjHGMbv3rOSeZb054XKRH+igNejT+V2uLhTudGgmFZYwjDGZTt1uXph0L/MCNhGMEHWsMHdX60XzW7o4HZq5AEsYxphMNWp2L2b8/T2HvGtT3BNWnad6THQ4KuOLTEsYIlIOmJ6sqCwwCPgWGA2E4llN72FVXZNK/YLAOKAioMD9qroqY6M2xqSXVRsmM3vDNBYGxECgixrHrmB4j/kUyhvmdGjGR5mWMFR1OxAJICIBwD5gNjAWeFlVF4pIS+B1oEEqhxgJLFLVjiISDFgnpzHZgDspiXHzX2LUkTkkBgiiyuvlXuDWqh3JG2KdHNmJU39bjYFdqrpXRBTI7y0vAOxPubOI5AfqAd0BVDUeiM+cUI0xl2rtxo9YsnkR05I2gQitEqpzZ83mVKt0l9OhmUvgd8IQkTxAnKomXcZ5OwFTve8fBxaLyHDABdySyv5lgb+AiSJSGYgG+qlqbCrx9QJ6AZQsWfIyQjTGXKq4hCQWr3iTF37/GIDCiW4eKdOfNvW6EBIY4HB05lJddEV0EXGJyD0iskBEDgLbgD9EZIuIvCEiN/hzQm93UhvgM2/RQ8ATqnot8AQwPpVqgUBVYJSqVgFigQGpHV9Vx6hqlKpGFStWzJ/QjDHppP/Elrzw+8cEqtL4VDleq/Eudzbqaskim/PlCmMZ8DUwENisqm4AESkMNASGichsVf3Ex3O2ANar6gHvdjegn/f9Z3gGtlOKAWJUdbV3eyZpJAxjjHPGzXuETw4t43CIi9Kn4YGbnqddvc5Oh2XSiS8Jo4mqJqQsVNW/gVnALBEJ8uOcnfm3Owo8Yxb18dwt1QjYkcq5/hSR30WknHfwvDHwsx/nNMZkoF/2rGf4gqdZH36QAm6oF1eSwZ0/plihwk6HZtLRRRNGasniUvYBEJFw4Dagd7LiB4GRIhIIxOEdfxCR4sA4VW3p3e9RYIq3S+tXoIcv5zTGZBx3UhKvT+vFwtM/cCyPcF18AIMafkjEDbWcDs1kAJ8HvUXkyVSK/wGiVXWjL8dQ1ZNAkRRlK4Fqqey7H2iZbHsjEOVrvMaYjDV+Xl+i/9zCipDDhAm0dd3B4F4vOx2WyUD+3CUV5X194d1uBawF+ojIZ6r6enoHZ4zJek7GJfDRopd5/5/lEAI3nQ7gv62XcePVRS5e2WRr/iSMIkBVVT0BICIv4Rl8rofnNldLGMbkYAmnY1m/ZRofrh3L2tBYCiW6eeDq7nRq1JuQ0PwXP4DJ9vxJGCU592G5BKCUqp4SkdPpG5YxJqvpP7UVX8lhCIX6p4rTs/FrtrBRLuNPwvgU+EFE5nq3bwemeh/kszuWjMmhFq54i0+3TGVjWBxlTkOrK++gd7fBTodlHOBzwlDV/4jIl8CtgAB9VHWd92Obk9iYHOZUXBwvT2vHAtkHYdAssTjPdZ5O4XwFnQ7NOMSfu6QEuBkooKpDRKSkiNRIbWZZY0z25U5K5OtV7/HRlk/4MfQ0lU+E0bPuBzSIsJsUczt/uqQ+ANx4Hq4bAhzH8+Be9QyIyxjjgCOxsbw0rRXLAj1jFR1d5RjQe7pN6WEA/xJGTVWtKiIbAFT1iPchOmNMDrBj91qeW/Ig20KTqHoiL3dHdKNFnV6I66JTzplcwp+EkeBdx0IBRKQYnisOY0w2djz2CE9ObcEPQbEQCp0DbuWBe9/kygK25Iw5lz8J4x08Cx5dKSJDgY7ACxkSlTEmU8xc+gHjf/2AmGCh+vHCNL6+PV2aP+50WCaL8ucuqSkiEo1n4j+Adqq6NWPCMsZkpF92r2Hk0qdYHnSU8EDlvqDaPN7nQ4IDrfvJpO2iCSONOaQAWohIC1V9K51jMsZkkPjTJ/jvtPuYzU4CApUqx/LTq8Fb3FrJJgs0F+fLFUY+75/l8NwRNc+7fTuwPCOCMsakvy+/n8KnP7/Fj8GeCRseLticu+56hYLhdu+K8Y0v05u/DCAiS/DMJXXcuz2Yf1fNM8ZkUVu2f8nnaz5iVtJmkoKFW/65jhfveJES15w3SbQxF3Q5c0nFA6XTNRpjTLpxu5Wf9/xMz/89Q6zLBSL0LdiFtu36cXWBMKfDM9mQPwnjY2CNiMzGc2tte2ByhkRljLks8aeP8/zHXVgUtBtcLm45fjVtKnekVd1eTodmsjF/7pIaKiILgbreoh6qusHX+iJSDpierKgsMAjP0qyjgVAgEXg4relGvM+BrAP2qWprX89tTG7y+dcv8dK+zyEIiia6aZuvLo/dNwqXS5wOzWRzvtwlJaqqAKq6Hlh/oX3S4l2LO9K7fwCwD89zHWOBl1V1oYi0xLOuRoM0DtMP2ArY5PvGpPBbzFreWvgCS4P341Ll1tgy9L9rMiWL2braJn34ctP1MhF5VERKJi8UkWARaSQik4Fufp63MbBLVffi6d46kwAKAPtTqyAiJfCs8jfOz3MZk+ON/KwXrZbez9Lg/YgqT1/Riff7fmHJwqQrX7qkmgP341n7ogxwFAjDk2yWAG/7uqZ3Mp2Aqd73jwOLRWS495i3pFFnBPAs/97mmyoR6QX0AihZsuSFdjUm29u4ZS5z1s1gFj8BUO9EY17r2p+8+a52ODKTE8lFepLO3VkkCCgKnFLVo5d0Qs+EhfuBCqp6QETeAb5T1VkichfQS1WbpKjTGmipqg+LSAPgaV/GMKKionTdunUX282YbOdY7ElGzxvItNNLSRAhb5Kb0fWmUq5UBUKDbGZZc+lEJFpVU53L3p+7pFDVBOCPy4ynBbBeVQ94t7vhGZsAz3MdqXU51QHaeMc4QoH8IvKJqt57mbEYk+38sT+aN758nq+C9hGmSlupT7tbOlH5+ginQzM5nF8JI5105t/uKPBcbdTHc7dUI2BHygqqOhAYCJDsCsOShclVDh6LY+mqN3j1j+lokFD+dCCP1XqHOhF1L17ZmHSQqQlDRMKB24DeyYofBEaKSCAQh3f8QUSKA+NUtWVmxmhMVnT02J88OaU1P4afJr9bqasR9G33Gtdeca3ToZlcxJ8lWh9Q1fGXczJVPQkUSVG2EjhvjgJV3Q+clyxU9Vs8VyPG5AoDJrZiges3XGFK1bhwelZ/jrpV2zkdlsmF/LnCeFNEuuB5uG4NMFVVt2RMWMaYCV88wvJ9a4kOO0nRRDeNg1vxfK/XELEH8Iwz/EkYh4H/AsF4HsCbISLvqOqHGRGYMbnV4X+OMWb+U3ya+AOEQaXTQYy5dyl5wws5HZrJ5fxJGP+o6jfe94tEZCSwGrCEYUw6SEqMZ/ayQXy2ezE/hyRybbzySOVRtKxhg9oma/B70FtE+uN5FqMAcDzdIzImF0pyK8NmdGJawg4IgY4SQZ87x3JlQVtX22Qdl3KX1Cw8U3u0BV5J33CMyX2Wr53IqPXvsjk0gRviXHS/6THa1H/A6bCMOY8/CaOQiFyrqjuBnSIyFtgALMiY0IzJ2ZKSEnjhk9uZzz5CgpXb9ToevWsyVxcq6HRoxqTKn4SRH/hWRA4BPwMFgaSMCMqYnEzdbhasGMmMbZ+yITSOqOMF6XLrSJpEVnU6NGMuyJ+E0RDYDNTEs763YlcXxvgl5uA+XplzDytC/oZQaJ50Da/2WUBgoM3/ZLI+fxZQ+sn7dpX3ZYzxw+qNn/LSulfYFyJE/VOM3g0HULNCE8TlyyoDxjjPibmkjMlVdh84zJC5LVgXcgqChAfzNOXejq9SOE+w06EZ4xdLGMZkoE8XvcGbf0wmPkSoebIAjcu2pnOzAU6HZcwl8WcuqUeAKap6JAPjMSZH2P1bNI8t6c6eICjsVu7I15z7Ow0jb4j9jmayL3/+9V4FrBWR9cAEYPHF1vE2JrdxJyXyzKSWLAn8g+AApcbREtxZsw/Na7VzOjRjLps/g94viMiLQFOgB/CeiMwAxqvqrowK0JjsYvy8/nx5YBG/BLsBuC9vTR7vcVkTPBuTpfi74p6KyJ/An3hmrS0EzBSRr1T12YwI0Jis7pddSxnz3TAWB/wJwVD96FW89cB8CuYNcTo0Y9KVP2MYj+FZTvUQnmVUn1HVBBFx4Vkl74IJQ0TKAdOTFZUFBuFZ22I0nqVXE4GHVXVNirrXAh/h6RZzA2NUdaSvsRuTEZISk5i5bDRjfxvFgUAhQJW3K75K1fLNKRAe5HR4xqQ7nxKGeCbgrwx0UNW9yT9TVbeItL7YMVR1O55p0RGRAGAfMBsYC7ysqgu9a3a/DjRIUT0ReEpV14tIPiDae1Xzsy/xG5Petu5YyNT/TWA22yBQaBhbgR4NulKlvC0QaXIunxKGtyuqSspkkezzrX6etzGwS1X3iojimXYEPDPg7k/l+H8Af3jfHxeRrcA1eKYoMSbTqCrrN83g/vX/wS1CuNtNryKtueOOIdYFZXI8f8YwVolIdVVdmw7n7QRM9b5/HFgsIsMBF3DLhSqKSGmgCp61OIzJNPGnj/PUpA6sCv4Dt0uo808lHm//Gjdda+tqm9zB37mk+ojIHiAWEDwXHxH+nFBEgoE2wEBv0UPAE6o6S0TuAsYDTdKomxfP9OqPq+qxNPbpBfQCKFmypD+hGZOmkZ/1ZHzsD2ioUCIB2hZsTe+ur9hyqSZXEV8fpRCRUqmVp9VNdYHjtAX6qmpT7/Y/QEFvt5fgWdkvfyr1goD5eJ7/eMuXc0VFRem6dev8Cc+Yc+zcvYJx34zga9d2glWpHleHV3t+QHiwPYBnciYRiVbVqNQ+8+dffbc0yof4GU9n/u2OAs+YRX08d0s1wnPH1Tm8iWQ8sNXXZGHM5UhISOSZia1ZGrIPAiHUDa9FDKZu1J1Oh2aMY/xJGLHJ3ocCrQG/BrtFJBy4DeidrPhBYKSIBAJxeLuTRKQ4ME5VWwJ1gPuATSKy0VvvOVX90p/zG+OLb1aNZO6WuXwT8hcAd7o68vSdfQkPL+pwZMY4y+cuqfMqioQA81S1WfqGlH6sS8r4Y9+R40xb/BKTTn8FwHXxwuiOy7nKVsAzuUh6dUmlFI7n4TtjsjV1u9m99zuGLBpIdHgseZPcdM7fmnvbPU7hAgWdDs+YLMOfJ7034VllDyAAKIb/4xfGZDnPftSMRfInhEON+Dw8dMswoio0cDosY7Icf64wkj/NnQgcUNXEdI7HmEyzbPX7TNowmfVhp7jxlIsGRVvxSNehdqusMWnwZ7Zav26fNSarOnU6gTemd2Vu0ibiw4QG8UX4b5c5FMhX0OnQjMnSfF5MWEQmi0jBZNuFRGRChkRlTAZZtGI490+qwWe6mVKnA3i2xKu8++C3liyM8YE/XVIRqnr0zIaqHhGRKukfkjHp72RcAsNm3MVs3QmhcHfgTTzZ/VPCQ21WWWN85U/CcIlIoTNLtIpIYT/rG5Pp1O1mwfJXGf/LdHaGKJEnwulxyxs0qlLP6dCMyXb8+YH/JvA/EZmJ526pu4ChGRKVMekgKSmB/h+1YLHrAOFBbjoQwQOdPqRkkfNmnjHG+MCf9TCWAevwTN8heNbGsOnFTZY0+5uRzPrlY34MO03N2EJ0jHqR5jVuczosY7I1f9bDmKOq1bA1KEwWtnffVl6e34O1obEQBk2TrmBoz0WEBttYhTGXy58uqR/ScT0MY9JVYkIcExY8y/gjS0kIEar9XYKnOoymUqnSTodmTI6R6ethGJPe5i7/iFlb32VDaBy4XPTOW5eHu76Py2UP4BmTnvxJGC0yLApjLsGBAz/x1fqpvHngCxJDhVrHr+Demt2oX72r06EZkyP5kzB+A7oAZVV1iIiUBK4C7Alwk+l+27+ZuxZ3JtblQoCBVz1Aq7v6UiDMxiqMySj+JIwPADeeu6SGAMfxLJdaPQPiMiZNL05uzxx2gstF/RPlaF65Da1vtasKYzKaPwmjpqpWFZENcPZJ7+AMisuY88xY/AJv7Z9NrMvF1QlK8/xRPNltktNhGZNr+JMwEkQkAO8U5yJSDM8Vh09EpBwwPVlRWWAQnqVZR+NZxS8ReFhV16RSvzkwEs/U6uNUdZgfsZts7K+/tvPK3Ef5OugPcLmocawoT7WfQfkSxZwOzZhcxZ+E8Q4wG7hSRIYCHYEXfa2sqtuBSABv4tnnPd5Y4GVVXSgiLYHXgQbJ63r3fx/P8q4xwFoRmWcPDuZ8oz5/idlHZvFHkCCqvFjiAe7s9oTTYRmTK/kzvfkUEYkGGnuL2qrqtks8b2Ngl6ruFREFzszVUADYn8r+NYCdqvorgIhMA9piDxHmWGt/msGUNR+yNOggBAq1j5bnjR4jKZD/KqdDMybXumjCEJF5KYu8fzYTEVS1zSWctxMw1fv+cWCxiAzHM936Lansfw3we7LtGKBmGvH2AnoBlCxZ8hJCM06KT3Tz7cblDNr0MrFBLq5McPN8lbeocGN9CuQLdTo8Y3I1X64wauP5YT0VWM2/CeOSeAfK2wADvUUPAU+o6iwRuQsYDzRJWS2VQ2kqZajqGGAMQFRUVKr7mKxH3W7W/jiBmRvmsDBgL7hc3J7YnD6tH6HkVaWcDs8Yg28J4yo8YwedgXuABcBUVd1yiedsAaxX1QPe7W5AP+/7z4BxqdSJAa5Ntl2C1LuuTDb16eKXGHZwDgTAdaeF26+6i+6tnyfAntY2Jsu4aMJQ1SRgEbBIRELwJI5vRWSIqr57CefszL/dUeD5wV8fz91SjYAdqdRZC9wgImXwDJZ3wpO8TDb318FtvPR5X1aFHiAQaJHUiMfuHMxVhQo7HZoxJgVfpzcPAVrh+WFfGs8dU5/7ezIRCcdztdI7WfGDwEgRCQTi8I4/iEhxPLfPtlTVRBF5BFiM57baCZdxhWOyAFWlz9hb+V/IMYJClWsShN439uH2+n2dDs0YkwZRvXA3v4hMBioCC4Fpqro5MwJLD1FRUbpu3TqnwzApTFv8NF/tXcGakJPkTXLTLKAtL3UdimfZFWOMk0QkWlWjUvvMlyuM+/DMTnsj8Fiy/9RnZqu15cuMTw4fO8a7sx9lFushBErFw5RO31EgX1GnQzPG+MCXMQxXZgRici53UiLL133A5B8/YV3IKQolunn0uiF0uPV2AgJtdhljsgt/nvQ2xm9xCUm88ukdzOZXCIEm8cV5psPHFC9yhdOhGWP8ZAnDZJhdu1fw8uKn2BB2ipvjAmh1bTfua/64LWxkTDZlCcOku4TEBIZM6cR83U5gCLTR63i2y0cUyFvA6dCMMZfBEoZJV9+vn8H7a19jU2g8EbFhNL3xWbrddqfTYRlj0oElDJMuTp1OYPj0Psx2ryYhVOgg1zP4oVmIy+6ZMCansIRhLos7KZEJC/oy+dBKjga4uPlUEE/Un0Tt8pWdDs0Yk84sYZhLFp+QyPOftGCR608IcHGXqxJd7nifslcUcjo0Y0wGsIRh/JaUGM+S/73PqG0T2B0CNWMLcF+Np6gf1d7p0IwxGcgShvHLb39sofeCu4kJEsKD3HRwXU//B2YQHhLidGjGmAxmCcP4xJ2UyMtTOrI4aQfxgULNv8tzX+P+1I+o6nRoxphMYgnDXNSsb95l1s6JbApJAJeLrsFVeeaJyU6HZYzJZJYwTJq2/fIlM9aMYXbiThJDhBrHi/DafZ9TtICtVWFMbmQJw6Rq2bp5vPDTQI4FuHABr18/gAbVOxEWHOB0aMYYh1jCMOc4fmwfM759hxFHvoQAF/WOV+Lu2h2oV62j06EZYxyWaQlDRMoB05MVlQUGAbWBct6ygsBRVY1Mpf4TQE9AgU1AD1WNy8CQc53layfTb8sbJIrgUqVfkYb06PqOLWxkjAEyMWGo6nYgEkBEAvCszT1bVUec2UdE3gT+SVlXRK4BHgPKq+opEZmBZ13vSRkeeC6QlBjPUxNvZ2nwfhCh7omy3FXvvzSoVMnp0IwxWYhTXVKNgV2quvdMgXh+jb0LaJRGnUAgTEQSgHBgf4ZHmQuMnfcKn/31KX8EC1cnKPdd3Z77uv3H6bCMMVmQUwmjEzA1RVld4ICq7ki5s6ruE5HhwG/AKWCJqi5J7cAi0gvoBVCyZMl0DTonid40k1H/e5PVoScgUKhxpASv9JjFlQXCnQ7NGJNFZXrCEJFgoA0wMMVHnTk/iZypUwhoC5QBjgKfici9qvpJyn1VdQwwBiAqKkrTL/KcQVV57/MhfHpsBidCXYS73Qyv/B+qV2pLaJDdAWWMSZsTVxgtgPWqeuBMgYgEAh2AamnUaQLsVtW/vPt/DtwCnJcwTNpmff0Cy/d8zzdBh8Dlok1iU57rNIA8eYo5HZoxJhtwImGkdiXRBNimqjFp1PkNqCUi4Xi6pBoD6zIuxJwl9nQiX62Zy+B9cyEISsVDv4iXaFz9Dlsu1Rjjs0xNGN4f+LcBvVN8dN6YhogUB8apaktVXS0iM4H1QCKwAW+3k0lbUmI8O35dzBvfvMKasBOEu93cFX47XVr146qiVzsdnjEmmxHVnNvNHxUVpevW5d4LkUEf385s9x4Aap3Kzx0VH6D5rfc7G5QxJksTkWhVjUrtM3vSOwf65n9vMGnLdH4MigMROgc05tkH3yIwwJZLNcZcOksYOUhikptXpjzI5+7VJAULlU+H8J+WH1OmRHmnQzPG5ACWMHKIH9aPZ3L0RFYG/0OJBGhX4j/0btHO6bCMMTmIJYxs7mjsaYbPuJe5rm0QDHVPF2ZE9yUEB9sKeMaY9GUJIxv7+ZclDF02kJ9C44mIDeHOSgNpW6c94rKxCmNM+rOEkQ2dOp3Aq9PuZDa7CAhR7nZF0rfrOArlDXU6NGNMDmYJI5tZ+sNEPvxxJFtDk6h8PJwm1z9C9+b3OR2WMSYXsISRTfx19C9en9mDrwL3kBQqdAy4jhf6zCIg0OZ/MsZkDksYWVxiQhyvz7yHqfE7IAgijofTu/6H1KtUxenQjDG5jCWMLGz3n78z8suuLA04BMADYTW5p/17XJHfxiqMMZnPEkYWdOL4H8z87g3GH1zM0QAXtU4U49Emg4go18Dp0IwxuZgljCzmyD/76DizKQcDXQS5hJ7hNeh59wfkCbXnKowxzrKEkUW4kxJ5/pPbmU8MBLpoGFuBFtW60aJmC6dDM8YYwBJGlvD1qgmM2PI2e4Mg3O2mReANDH54mtNhGWPMOSxhOOj3mNWM/OoFFrv+hCCodrQo/e+czc3FCzodmjHGnMcShkOmLH6fj34fxf4gz4p3/Qq2oMd9rxFgK+AZY7KoTEsYIlIOmJ6sqCwwCKgNlPOWFQSOqmpkKvULAuOAioAC96vqqoyLOGPE7FvNnB/G8uGJ1RAk1DsewdNt+lGmZA2nQzPGmAvKtIShqtuBSAARCQD2AbNVdcSZfUTkTeCfNA4xElikqh1FJBgIz9CA01lcQhKzl89i5G8vE+tyEepWht7Un1pVOpE/NMjp8Iwx5qKc6pJqDOxS1b1nCkREgLuARil3FpH8QD2gO4CqxgPxmRJpOog7dYRXZvRlNpvA5aLJyeo0rnw7TWu3dzo0Y4zxmVMJoxMwNUVZXeCAqu5IZf+ywF/ARBGpDEQD/VQ1NuWOItIL6AVQsmTJdA3aX6rKjKVjeCXmXdwiFEpy07dER+5u+rKjcRljzKXI9IUTvN1JbYDPUnzUmfOTyBmBQFVglKpWAWKBAantqKpjVDVKVaOKFSuWTlH7b9fuFXQaU5X/7nsPtwh1j1/H2KZfWbIwxmRbTlxhtADWq+qBMwUiEgh0AKqlUScGiFHV1d7tmaSRMJymqrw0qRtfazTHQ12UjoeHyvWhZb2+TodmjDGXxYmEkdqVRBNgm6rGpFZBVf8Ukd9FpJx38Lwx8HMGx+m3CV/04qs/17I5OBFRoVlcdd7oNR7P8IwxxmRvmZowRCQcuA3oneKj88Y0RKQ4ME5VW3qLHgWmeLu0fgV6ZHC4PktISGT03MGMiV0FwVA8QZlyxyKKFirhdGjGGJNuMjVhqOpJoEgq5d1TKdsPtEy2vRGIysDw/KZuN6s3jGXUhrGsDzpNqFt5sFg/ujW5m5DQ/E6HZ4wx6cqe9L5EqsqoOQ8x6vj/IAganCpOn6ZDqVA2S+U0Y4xJN5YwLkH0T9MYuepNNoWc4spEpX3hbjzY5WmCAzP9pjNjjMk0ljD89PzE9iyQHSSFCvXii/BY0/cpV6qi02EZY0yGs4Tho1XrJzAheiw/BJ8g4ng4rSLe4J569Z0OyxhjMo0ljIs4ciKOoZ/ex1fBW3EHC7VO5+Odnt8QFmrrahtjchdLGBewesNkRqwdyeaQBCJiw+hUZRitazZEXDZWYYzJfSxhpOLw8ThGzLqXObIdQqA91/FMjxnkCwt2OjRjjHGMJYwUjsTGc8eMahwOdFHxVBC3XXsf97d6wumwjDHGcZYwUsgX7KaqXEWJsGt4rMt4AgNtrQpjjAFLGOcJDArlrfuXOh2GMcZkOTZ6a4wxxieWMIwxxvjEEoYxxhifWMIwxhjjE0sYxhhjfGIJwxhjjE8sYRhjjPGJJQxjjDE+EVV1OoYMIyJ/AXsvsXpR4FA6hpMdWJtzvtzWXrA2+6uUqhZL7YMcnTAuh4isU9Vctd6qtTnny23tBWtzerIuKWOMMT6xhGGMMcYnljDSNsbpABxgbc75clt7wdqcbmwMwxhjjE/sCsMYY4xPLGEYY4zxiSWMFESkuYhsF5GdIjLA6XjSi4hcKyLLRGSriGwRkX7e8sIi8pWI7PD+WShZnYHe72G7iDRzLvpLJyIBIrJBROZ7t3N0ewFEpKCIzBSRbd6/79o5ud0i8oT33/RmEZkqIqE5sb0iMkFEDorI5mRlfrdTRKqJyCbvZ++IiPgchKray/sCAoBdQFkgGPgRKO90XOnUtquBqt73+YBfgPLA68AAb/kA4DXv+/Le9ocAZbzfS4DT7biEdj8JfArM927n6PZ62zIZ6Ol9HwwUzKntBq4BdgNh3u0ZQPec2F6gHlAV2JyszO92AmuA2oAAC4EWvsZgVxjnqgHsVNVfVTUemAa0dTimdKGqf6jqeu/748BWPP/Z2uL5AYP3z3be922Baap6WlV3AzvxfD/ZhoiUAFoB45IV59j2AohIfjw/WMYDqGq8qh4lZ7c7EAgTkUAgHNhPDmyvqi4H/k5R7Fc7ReRqIL+qrlJP9vgoWZ2LsoRxrmuA35Ntx3jLchQRKQ1UAVYDV6rqH+BJKsAV3t1ywncxAngWcCcry8ntBc/V8V/ARG9X3DgRyUMObbeq7gOGA78BfwD/qOoScmh7U+FvO6/xvk9Z7hNLGOdKrS8vR913LCJ5gVnA46p67EK7plKWbb4LEWkNHFTVaF+rpFKWbdqbTCCebotRqloFiMXTVZGWbN1ub599WzzdLsWBPCJy74WqpFKWbdrrh7TaeVntt4Rxrhjg2mTbJfBc3uYIIhKEJ1lMUdXPvcUHvJepeP886C3P7t9FHaCNiOzB07XYSEQ+Iee294wYIEZVV3u3Z+JJIDm13U2A3ar6l6omAJ8Dt5Bz25uSv+2M8b5PWe4TSxjnWgvcICJlRCQY6ATMczimdOG9E2I8sFVV30r20Tygm/d9N2BusvJOIhIiImWAG/AMlmULqjpQVUuoamk8f4/fqOq95ND2nqGqfwK/i0g5b1Fj4Gdybrt/A2qJSLj333hjPONzObW9KfnVTm+31XERqeX9vromq3NxTo/8Z7UX0BLPHUS7gOedjicd23UrnkvPn4CN3ldLoAiwFNjh/bNwsjrPe7+H7fhxJ0VWewEN+PcuqdzQ3khgnffveg5QKCe3G3gZ2AZsBj7Gc2dQjmsvMBXPOE0CniuFBy6lnUCU97vaBbyHd8YPX142NYgxxhifWJeUMcYYn1jCMMYY4xNLGMYYY3xiCcMYY4xPLGEYY4zxiSUMYy5CRIqIyEbv608R2ed9f0JEPsigcz4uIl0vss80EbkhI85vTGrstlpj/CAig4ETqjo8A88RCKzHM7tw4gX2qw/cq6oPZlQsxiRnVxjGXCIRaZBsnY3BIjJZRJaIyB4R6SAir3vXHVjknZblzFoE34lItIgsPjOtQwqNgPWqmigi14nI+mTnvEFEzsyPtQJo4k0wxmQ4SxjGpJ/r8Eyn3hb4BFimqpWAU0Arb9J4F+ioqtWACcDQVI5TB4gGUNVdwD8iEun9rAcwyfuZG8+01ZUzqD3GnMN+MzEm/SxU1QQR2YRnMa5F3vJNQGmgHFAR+Mq7yFkAnqkeUroaz3xIZ4wDeojIk8DdnLt+w0E8s7T6OiuvMZfMEoYx6ec0eH7zF5EE/XeA0I3n/5oAW1S19kWOcwoITbY9C3gJ+AaIVtXDyT4L9e5vTIazLiljMs92oJiI1AbPdPMiUiGV/bYC15/ZUNU4YDEwCpiYYt8bgS0ZE64x57KEYUwmUc+yvx2B10TkRzwzBt+Syq4L8SyzmtwUPLMNLzlTICJXAqfUu+KaMRnNbqs1JgsSkdnAs6q6w7v9NFBAVV9Mts8TwDFVHe9QmCaXsTEMY7KmAXgGv3d4k8d1eG63Te4onvUfjMkUdoVhjDHGJzaGYYwxxieWMIwxxvjEEoYxxhifWMIwxhjjE0sYxhhjfPJ/ZMu6uD4iMSIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 6c81a1a11..d6d147634 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -99,8 +99,8 @@ module swiftest_classes !******************************************************************************************************************************** !> A concrete lass for the central body in a Swiftest simulation type, abstract, public, extends(swiftest_base) :: swiftest_cb - character(len=STRMAX) :: name !! Non-unique name - integer(I4B) :: id !! External identifier (unique) + character(len=STRMAX) :: name !! Non-unique name + integer(I4B) :: id = 0 !! External identifier (unique) real(DP) :: mass = 0.0_DP !! Central body mass (units MU) real(DP) :: Gmass = 0.0_DP !! Central mass gravitational term G * mass (units GU * MU) real(DP) :: radius = 0.0_DP !! Central body radius (units DU) diff --git a/src/whm/whm_getacch.f90 b/src/whm/whm_getacch.f90 index 26a3acb19..c4ee90592 100644 --- a/src/whm/whm_getacch.f90 +++ b/src/whm/whm_getacch.f90 @@ -32,7 +32,7 @@ module subroutine whm_getacch_pl(self, system, param, t, lbeg) call whm_getacch_ah2(cb, pl) call whm_getacch_ah3(pl) - if (param%loblatecb) then + if (param%loblatecb) then cb%aoblbeg = cb%aobl call pl%accel_obl(system) cb%aoblend = cb%aobl @@ -97,13 +97,14 @@ function whm_getacch_ah0(mu, xhp, n) result(ah0) ! Result real(DP), dimension(NDIM) :: ah0 ! Internals - real(DP) :: fac, r2, ir3h + real(DP) :: fac, r2, ir3h, irh integer(I4B) :: i ah0(:) = 0.0_DP do i = 1, n r2 = dot_product(xhp(:, i), xhp(:, i)) - ir3h = 1.0_DP / (r2 * sqrt(r2)) + irh = 1.0_DP / sqrt(r2) + ir3h = irh / r2 fac = mu(i) * ir3h ah0(:) = ah0(:) - fac * xhp(:, i) end do