From 6658513060dad9dee17f797aed7c85198c40d6a1 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 17 Jun 2021 18:14:27 -0400 Subject: [PATCH] Wrote a new swiftestio method that can create a Dataset from JPL/Horizons data. Updated inititial conditions generator to use it in one example. Only the Swiftest version has been completed --- .../1pl_1tp_encounter/.idea/.gitignore | 3 + .../9pl_18tp_encounters/cb.swiftest.in | 8 +- .../9pl_18tp_encounters/cb_test.in | 4 + .../9pl_18tp_encounters/config.swiftest.in | 2 +- .../9pl_18tp_encounters/init_cond.py | 252 ++++--------- .../9pl_18tp_encounters/pl.swiftest.in | 72 ++-- .../9pl_18tp_encounters/pl_test.in | 37 ++ .../swiftest_rmvs_vs_swifter_rmvs.ipynb | 13 +- .../9pl_18tp_encounters/tp.in | 57 ++- .../9pl_18tp_encounters/tp_test.in | 1 + .../mars_ejecta/.idea/.gitignore | 3 + python/swiftestio/swiftestio.py | 340 +++++++++++++++--- 12 files changed, 525 insertions(+), 267 deletions(-) create mode 100644 examples/rmvs_swifter_comparison/1pl_1tp_encounter/.idea/.gitignore create mode 100644 examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb_test.in create mode 100644 examples/rmvs_swifter_comparison/9pl_18tp_encounters/pl_test.in create mode 100644 examples/rmvs_swifter_comparison/9pl_18tp_encounters/tp_test.in create mode 100644 examples/rmvs_swifter_comparison/mars_ejecta/.idea/.gitignore diff --git a/examples/rmvs_swifter_comparison/1pl_1tp_encounter/.idea/.gitignore b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/.idea/.gitignore new file mode 100644 index 000000000..26d33521a --- /dev/null +++ b/examples/rmvs_swifter_comparison/1pl_1tp_encounter/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb.swiftest.in b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb.swiftest.in index aa4a2fa2a..689d47628 100644 --- a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb.swiftest.in +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb.swiftest.in @@ -1,4 +1,4 @@ -0.00029591220819207775568 -0.0046504672609621575315 -4.753580694812735539e-12 --2.2473967953572828617e-18 +0.0002959122081920778 +0.004650467260962157 +4.7535806948127355e-12 +-2.2473967953572827e-18 diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb_test.in b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb_test.in new file mode 100644 index 000000000..689d47628 --- /dev/null +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb_test.in @@ -0,0 +1,4 @@ +0.0002959122081920778 +0.004650467260962157 +4.7535806948127355e-12 +-2.2473967953572827e-18 diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/config.swiftest.in b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/config.swiftest.in index bce02272d..4b7e90915 100644 --- a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/config.swiftest.in +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/config.swiftest.in @@ -1,6 +1,6 @@ ! Swiftest input file generated using init_cond.py T0 0 -TSTOP 80.0 +TSTOP 365.25 DT 1.0 CB_IN cb.swiftest.in PL_IN pl.swiftest.in diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/init_cond.py b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/init_cond.py index a4b662dd4..776684e05 100644 --- a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/init_cond.py +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/init_cond.py @@ -3,11 +3,25 @@ To use the script, modify the variables just after the "if __name__ == '__main__':" line """ import numpy as np -from astroquery.jplhorizons import Horizons -import astropy.constants as const import swiftestio as swio -from scipy.io import FortranFile +import astropy.constants as const import sys +import xarray as xr + +# Swiftest configuration file +# Both codes use the same tp input file +tpin = "tp.in" + +swifter_input = "param.swifter.in" +swifter_pl = "pl.swifter.in" +swifter_bin = "bin.swifter.dat" +swifter_enc = "enc.swifter.dat" + +swiftest_input = "config.swiftest.in" +swiftest_pl = "pl.swiftest.in" +swiftest_cb = "cb.swiftest.in" +swiftest_bin = "bin.swiftest.dat" +swiftest_enc = "enc.swiftest.dat" #Values from JPL Horizons AU2M = np.longdouble(const.au.value) @@ -24,70 +38,6 @@ GU = np.longdouble(GC / (DU2M**3 / (MU2KG * TU2S**2))) GMSun = np.longdouble(GMSunSI / (DU2M**3 / TU2S**2)) -# Solar oblatenes values: From Mecheri et al. (2004), using Corbard (b) 2002 values (Table II) -J2 = np.longdouble(2.198e-7) * (Rsun / DU2M)**2 -J4 = np.longdouble(-4.805e-9) * (Rsun / DU2M)**4 - -npl = 9 -ntp = 2 * npl - -# Planet ids -plid = { - 'mercury' : '1', - 'venus' : '2', - 'earthmoon' : '3', - 'mars' : '4', - 'jupiter' : '5', - 'saturn' : '6', - 'uranus' : '7', - 'neptune' : '8', - 'plutocharon' : '9' -} - -#Planet Msun/M ratio -MSun_over_Mpl = { - 'mercury' : np.longdouble(6023600.0), - 'venus' : np.longdouble(408523.71), - 'earthmoon' : np.longdouble(328900.56), - 'mars' : np.longdouble(3098708.), - 'jupiter' : np.longdouble(1047.3486), - 'saturn' : np.longdouble(3497.898), - 'uranus' : np.longdouble(22902.98), - 'neptune' : np.longdouble(19412.24), - 'plutocharon' : np.longdouble(1.35e8) -} - -#Planet radii in meters -Rpl = { - 'mercury' : np.longdouble(2439.4e3), - 'venus' : np.longdouble(6051.8e3), - 'earthmoon' : np.longdouble(6371.0084e3), # Earth only for radius - 'mars' : np.longdouble(3389.50e3), - 'jupiter' : np.longdouble(69911e3), - 'saturn' : np.longdouble(58232.0e3), - 'uranus' : np.longdouble(25362.e3), - 'neptune' : np.longdouble(24622.e3), - 'plutocharon' : np.longdouble(1188.3e3) -} - -THIRDLONG = np.longdouble(1.0) / np.longdouble(3.0) - -# Both codes use the same tp input file -tpin = "tp.in" - -swifter_input = "param.swifter.in" -swifter_pl = "pl.swifter.in" -swifter_bin = "bin.swifter.dat" -swifter_enc = "enc.swifter.dat" - -swiftest_input = "config.swiftest.in" -swiftest_pl = "pl.swiftest.in" -swiftest_cb = "cb.swiftest.in" -swiftest_bin = "bin.swiftest.dat" -swiftest_enc = "enc.swiftest.dat" - -# Simple initial conditions of a circular planet with one test particle in a close encounter state -# Simulation start, stop, and output cadence times t_0 = 0 # simulation start time deltaT = 1.00 * JD / TU2S # simulation step size end_sim = year / TU2S # simulation end time @@ -102,112 +52,6 @@ tend = '2021-06-16' tstep = '1d' -####################################################### -# Start generating initial conditions for the planets # -####################################################### -pldata = {} -tpdata = {} -plvec = {} -tpvec = {} -Rhill = {} - -for key,val in plid.items(): - pldata[key] = Horizons(id=val, id_type='majorbody',location='@sun', - epochs={'start': tstart, 'stop': tend, - 'step': tstep}) - plvec[key] = np.array([ - pldata[key].vectors()['x'][0], - pldata[key].vectors()['y'][0], - pldata[key].vectors()['z'][0], - pldata[key].vectors()['vx'][0], - pldata[key].vectors()['vy'][0], - pldata[key].vectors()['vz'][0] - ]) - Rhill[key] = pldata[key].elements()['a'][0] * (3 * MSun_over_Mpl[key])**(-1.0 / 3.0) - -# For each planet, we will initialize a pair of test particles. One on its way in, and one on its way out. We will also initialize two additional particles that don't encounter anything -tpnames = range(101, 101 + ntp) -tpxv1 = np.empty((6)) -tpxv2 = np.empty((6)) -for idx,key in enumerate(plid): - rstart = 2 * Rpl[key] / DU2M # Start the test particles at a multiple of the planet radius away - vstart = 1.5 * np.sqrt(2 * GMSun / MSun_over_Mpl[key] / rstart) # Start the test particle velocities at a multiple of the escape speed - xvstart = np.array([rstart / np.sqrt(2.0), rstart / np.sqrt(2.0), 0.0, vstart, 0.0, 0.0]) - # The positions and velocities of each pair of test particles will be in reference to a planet - tpxv1 = plvec[key] + xvstart - tpxv2 = plvec[key] - xvstart - tpvec[tpnames[idx]] = tpxv1 - tpvec[tpnames[idx + npl]] = tpxv2 - -# TP file -tpfile = open(tpin, 'w') -print(ntp, file=tpfile) -for tpid, tp in tpvec.items(): - print(tpid, file=tpfile) - print(f'{tp[0]} {tp[1]} {tp[2]}', file=tpfile) - print(f'{tp[3]} {tp[4]} {tp[5]}', file=tpfile) -tpfile.close() - -# Swifter PL file -plfile = open(swifter_pl, 'w') -print(npl + 1, file=plfile) -print(1,GMSun,file=plfile) -print('0.0 0.0 0.0',file=plfile) -print('0.0 0.0 0.0',file=plfile) -for key, pl in plvec.items(): - print(f'{int(plid[key])+1} {GMSun / MSun_over_Mpl[key]} {Rhill[key]} ! {key}', file=plfile) - print(f'{Rpl[key] / DU2M}', file=plfile) - print(f'{pl[0]} {pl[1]} {pl[2]}', file=plfile) - print(f'{pl[3]} {pl[4]} {pl[5]}', file=plfile) -plfile.close() - -# Swiftest Central body and pl file -cbfile = open(swiftest_cb, 'w') -print(GMSun, file=cbfile) -print(rmin, file=cbfile) -print(J2, file=cbfile) -print(J4, file=cbfile) -cbfile.close() -plfile = open(swiftest_pl, 'w') -print(npl, file=plfile) -for key, pl in plvec.items(): - print(f'{int(plid[key]) + 1} {GMSun / MSun_over_Mpl[key]} ! {key}', file=plfile) - print(f'{Rpl[key] / DU2M}', file=plfile) - print(f'{pl[0]} {pl[1]} {pl[2]}', file=plfile) - print(f'{pl[3]} {pl[4]} {pl[5]}', file=plfile) -plfile.close() - -# Swifter parameter file -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 {tpin}') -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 XV') -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') -sys.stdout = sys.__stdout__ - -# Swiftest configuration file sys.stdout = open(swiftest_input, "w") print(f'! Swiftest input file generated using init_cond.py') print(f'T0 {t_0} ') @@ -238,6 +82,68 @@ print(f'MU2KG {MU2KG}') print(f'DU2M {DU2M}') print(f'TU2S {TU2S}') +sys.stdout = sys.__stdout__ +config = swio.read_swiftest_config(swiftest_input) +ds = swio.solar_system_pl(config, tstart) + +ntp = 18 +dims = ['time', 'id', 'vec'] +tp = [] +t = np.array([0.0]) +clab, plab, tlab = swio.make_swiftest_labels(config) + +# For each planet, we will initialize a pair of test particles. One on its way in, and one on its way out. We will also initialize two additional particles that don't encounter anything +tpnames = np.arange(101, 101 + ntp) +tpxv1 = np.empty((6)) +tpxv2 = np.empty((6)) +cb = ds.sel(id=0) +pl = ds.where(ds.id > 0, drop=True) + +p1 = [] +p2 = [] +p3 = [] +p4 = [] +p5 = [] +p6 = [] + +for i in pl.id: + pli = pl.sel(id=i) + rstart = 2 * np.double(pli['Radius']) # Start the test particles at a multiple of the planet radius away + vstart = 1.5 * np.sqrt(2 * np.double(pli['Mass']) / rstart) # Start the test particle velocities at a multiple of the escape speed + xvstart = np.array([rstart / np.sqrt(2.0), rstart / np.sqrt(2.0), 0.0, vstart, 0.0, 0.0]) + # The positions and velocities of each pair of test particles will be in reference to a planet + plvec = np.array([np.double(pli['px']), + np.double(pli['py']), + np.double(pli['pz']), + np.double(pli['vx']), + np.double(pli['vy']), + np.double(pli['vz'])]) + tpxv1 = plvec + xvstart + tpxv2 = plvec - xvstart + p1.append(tpxv1[0]) + p1.append(tpxv2[0]) + p2.append(tpxv1[1]) + p2.append(tpxv2[1]) + p3.append(tpxv1[2]) + p3.append(tpxv2[2]) + p4.append(tpxv1[3]) + p4.append(tpxv2[3]) + p5.append(tpxv1[4]) + p5.append(tpxv2[4]) + p6.append(tpxv1[5]) + p6.append(tpxv2[5]) + +tvec = np.vstack([p1, p2, p3, p4, p5, p6]) +tpframe = np.expand_dims(tvec.T, axis=0) +tpxr = xr.DataArray(tpframe, dims = dims, coords = {'time' : t, 'id' : tpnames, 'vec' : tlab}) + +tp = [tpxr] +tpda = xr.concat(tp,dim='time') +tpds = tpda.to_dataset(dim = 'vec') + +ds = xr.combine_by_coords([ds, tpds]) +swio.swiftest_xr2_infile(ds, config) + diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in index b93994014..385f84350 100644 --- a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/pl.swiftest.in @@ -1,37 +1,37 @@ 9 -2 4.9125474498983625e-11 ! mercury -1.6306381826061646e-05 -0.008059842448018334 -0.4616051037329109 -0.03846017738329229 -0.02248719132054853 0.001934639213990692 -0.001904656977422976 -3 7.243452483873647e-10 ! venus -4.0453784346544176e-05 --0.5115875215389065 0.5030818749037324 0.03642547299277956 --0.01425515725454357 -0.01452868630179309 0.0006232072038298823 -4 8.997011382166019e-10 ! earthmoon -4.25875607065041e-05 --0.1090020607540907 -1.009893805009766 4.823302918632528e-05 -0.01682491922568941 -0.001910549762056979 3.992660742687128e-08 -5 9.549535102761465e-11 ! mars -2.2657408050928896e-05 --1.342897929331636 0.9778655112682739 0.05343398538723887 --0.007712315645393206 -0.01011917844182223 -2.287744801261131e-05 -6 2.8253459086313547e-07 ! jupiter -0.0004673261703049093 -3.923184193414315 -3.168419770483168 -0.0746147877972047 -0.004655552638985802 0.006232623300954468 -0.0001300429201057457 -7 8.459715183006416e-08 ! saturn -0.00038925687730393614 -6.185794462795267 -7.804174837804826 -0.110498432926239 -0.004066833203985018 0.003458637040736611 -0.0002219310939327014 -8 1.2920249163736674e-08 ! uranus -0.00016953449859497232 -14.9290976575471 12.92949673572929 -0.1454099139559955 --0.002599557960646664 0.002795888198858545 4.391864857782088e-05 -9 1.5243589003230834e-08 ! neptune -0.00016458790412449367 -29.54416169025338 -4.716921603714237 -0.5838030174427992 -0.0004792636209523189 0.00312573757291745 -7.53264045199501e-05 -10 2.1919422829042796e-12 ! plutocharon -7.943294877391593e-06 -14.54448346259197 -31.05223519593471 -0.8828000265625595 -0.002923077617691739 0.0006625916902153526 -0.0009142553677224461 +1 4.9125474498983625056e-11 +2439400.0 +0.008059842448018333591 -0.46160510373291091524 -0.038460177383292291908 +0.02248719132054852949 0.0019346392139906921279 -0.0019046569774229759606 +2 7.243452483873647106e-10 +6051800.0 +-0.51158752153890652004 0.5030818749037323512 0.036425472992779560355 +-0.01425515725454356988 -0.014528686301793089943 0.00062320720382988232425 +3 8.997011382166018993e-10 +6371008.4000000003725 +-0.109002060754090704386 -1.0098938050097661101 4.8233029186325282966e-05 +0.016824919225689409317 -0.0019105497620569790936 3.9926607426871282392e-08 +4 9.549535102761465872e-11 +3389500.0 +-1.3428979293316360977 0.97786551126827392366 0.053433985387238869258 +-0.007712315645393206092 -0.0101191784418222296525 -2.2877448012611311785e-05 +5 2.8253459086313549713e-07 +69911000.0 +3.923184193414314791 -3.1684197704831680298 -0.07461478779720470689 +0.0046555526389858022113 0.006232623300954467766 -0.00013004292010574569365 +6 8.45971518300641563e-08 +58232000.0 +6.1857944627952665684 -7.804174837804826126 -0.11049843292623899582 +0.0040668332039850179674 0.0034586370407366113193 -0.00022193109393270141328 +7 1.2920249163736673984e-08 +25362000.0 +14.929097657547099942 12.9294967357292893695 -0.14540991395599550673 +-0.0025995579606466640267 0.0027958881988585450113 4.391864857782088156e-05 +8 1.5243589003230834746e-08 +24622000.0 +29.544161690253378794 -4.7169216037142369657 -0.58380301744279916587 +0.00047926362095231893815 0.00312573757291745008 -7.532640451995010825e-05 +9 2.1919422829042797324e-12 +1188300.0 +14.544483462591969669 -31.052235195934709822 -0.88280002656255951443 +0.0029230776176917390448 0.0006625916902153525834 -0.0009142553677224461557 diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/pl_test.in b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/pl_test.in new file mode 100644 index 000000000..02247b2d1 --- /dev/null +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/pl_test.in @@ -0,0 +1,37 @@ +9 +1 4.9125474498983625056e-11 +2439400.0 +-0.22865268463418078615 -0.40059887459934212517 -0.011759565382932279237 +0.018760454180373070532 -0.012597043336814979538 -0.002750330975414844079 +2 7.243452483873647106e-10 +6051800.0 +0.055053104402360240888 -0.72501289869517293596 -0.013125811934180060814 +0.020032783376563709821 0.0014585960987979049513 -0.0011360168603540260125 +3 8.997011382166018993e-10 +6371008.4000000003725 +-0.0709344628437889313 -1.013495829860109998 4.7192964279974951866e-05 +0.01688220529198479184 -0.0012655943533160610018 2.2557847569638601276e-08 +4 9.549535102761465872e-11 +3389500.0 +0.82577293992211286966 -1.1262220087443419736 -0.043856929963955730567 +0.011815540916412620865 0.00947446278988147595 -9.133327435156872944e-05 +5 2.8253459086313549713e-07 +69911000.0 +1.7456184157900889176 -4.862064645732346868 -0.018861599055565350658 +0.007018704991297365326 0.0029091005451762219292 -0.00016911683181977029383 +6 8.45971518300641563e-08 +58232000.0 +4.5973377487897115756 -8.8976699986865348535 -0.028288312634969701304 +0.0046547417429797477434 0.0025487687199604531773 -0.000229509604752460194 +7 1.2920249163736673984e-08 +25362000.0 +15.833750897756699416 11.883335170209360143 -0.16096291202065848847 +-0.0023831013046656290505 0.0029651972909487032979 4.1744435577025254117e-05 +8 1.5243589003230834746e-08 +24622000.0 +29.348213377604238872 -5.8472894170428588723 -0.5560318242823257817 +0.0006000134768904796416 0.003101182064622641169 -7.765884757985849474e-05 +9 2.1919422829042797324e-12 +1188300.0 +13.476001769002650121 -31.277013495742810534 -0.5505401447793714098 +0.0029633486183866481637 0.00057521023255640399144 -0.00091622996928541103474 diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb index 28aaaa618..56bf48a02 100644 --- a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/swiftest_rmvs_vs_swifter_rmvs.ipynb @@ -12,6 +12,13 @@ "import matplotlib.pyplot as plt" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": 2, @@ -722,9 +729,9 @@ ], "metadata": { "kernelspec": { - "display_name": "swiftestOOF", + "display_name": "Python 3", "language": "python", - "name": "swiftestoof" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -736,7 +743,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.10" + "version": "3.8.6" } }, "nbformat": 4, diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/tp.in b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/tp.in index a71e3325c..63d1af3d3 100644 --- a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/tp.in +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/tp.in @@ -1,4 +1,55 @@ -1 +18 +101 +3449832.5721127702855 3449832.102447824087 -0.038460177383292291908 +0.02248719805191334506 0.00193463921399069207 -0.0019046569774229759036 +102 +-3449832.5559930852614 -3449833.02565803146 -0.038460177383292291908 +0.022487184589183712574 0.00193463921399069207 -0.0019046569774229759036 +103 +8558537.125181974843 8558538.139851370826 0.036425472992779560355 +-0.0142551408440447086035 -0.014528686301793089508 0.0006232072038298823056 +104 +-8558538.148357016966 -8558537.1336876209825 0.036425472992779560355 +-0.014255173665042430303 -0.014528686301793089508 0.0006232072038298823056 +105 +9009966.37627085112 9009965.47537910752 4.8233029186325282966e-05 +0.016824937050967138374 -0.0019105497620569790364 3.9926607426871281197e-08 +106 +-9009966.594274973497 -9009967.495166717097 4.8233029186325282966e-05 +0.016824901400411679253 -0.0019105497620569790364 3.9926607426871281197e-08 +107 +4793475.5267656762153 4793477.847529117018 0.053433985387238869258 +-0.0077123076835328777112 -0.01011917844182222935 -2.28774480126113111e-05 +108 +-4793478.2125615347177 -4793475.891798093915 0.053433985387238869258 +-0.0077123236072535340108 -0.01011917844182222935 -2.28774480126113111e-05 +109 +98869088.2822496295 98869081.19064567983 -0.07461478779720470689 +0.0046556479963669027827 0.0062326233009544675795 -0.00013004292010574568976 +110 +-98869080.43588125706 -98869087.52748520672 -0.07461478779720470689 +0.004655457281604701361 0.0062326233009544675795 -0.00013004292010574568976 +111 +82352490.3499045223 82352476.359935224056 -0.11049843292623899582 +0.0040668903766289789606 0.0034586370407366112158 -0.00022193109393270140663 +112 +-82352477.97831560671 -82352491.96828490496 -0.11049843292623899582 +0.004066776031341056731 0.0034586370407366112158 -0.00022193109393270140663 +113 +35867299.298004090786 35867297.298403166234 -0.14540991395599550673 +-0.0025995241047005760923 0.0027958881988585449277 4.3918648577820880246e-05 +114 +-35867269.439808771014 -35867271.439409695566 -0.14540991395599550673 +-0.0025995918165927518056 0.0027958881988585449277 4.3918648577820880246e-05 +115 +34820795.876912035048 34820761.615828737617 -0.58380301744279916587 +0.00047930094366581383431 0.0031257375729174499863 -7.532640451995010599e-05 +116 +-34820736.788588650525 -34820771.049671947956 -0.58380301744279916587 +0.0004792262982388240133 0.0031257375729174499863 -7.532640451995010599e-05 117 -29.543928927807162 -4.717154366160452 -0.5838030174427992 --0.013956372192325811 0.00312573757291745 -7.53264045199501e-05 +1680524.5206514112651 1680478.9239327528048 -0.88280002656255951443 +0.0029230796549344268721 0.00066259169021535256356 -0.0009142553677224461283 +118 +-1680495.4316844861023 -1680541.0284031445626 -0.88280002656255951443 +0.0029230755804490510426 0.00066259169021535256356 -0.0009142553677224461283 diff --git a/examples/rmvs_swifter_comparison/9pl_18tp_encounters/tp_test.in b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/tp_test.in new file mode 100644 index 000000000..573541ac9 --- /dev/null +++ b/examples/rmvs_swifter_comparison/9pl_18tp_encounters/tp_test.in @@ -0,0 +1 @@ +0 diff --git a/examples/rmvs_swifter_comparison/mars_ejecta/.idea/.gitignore b/examples/rmvs_swifter_comparison/mars_ejecta/.idea/.gitignore new file mode 100644 index 000000000..26d33521a --- /dev/null +++ b/examples/rmvs_swifter_comparison/mars_ejecta/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/python/swiftestio/swiftestio.py b/python/swiftestio/swiftestio.py index 4c7a51bd0..e38a2cf04 100644 --- a/python/swiftestio/swiftestio.py +++ b/python/swiftestio/swiftestio.py @@ -4,15 +4,20 @@ import xarray as xr from astroquery.jplhorizons import Horizons import astropy.constants as const +import datetime -# Constants, including values from JPL Horizons +# Constants in SI units AU2M = np.longdouble(const.au.value) GMSunSI = np.longdouble(const.GM_sun.value) -Rsun = np.longdouble(const.R_sun.value) +RSun = np.longdouble(const.R_sun.value) GC = np.longdouble(const.G.value) -JD = 86400 -year = np.longdouble(365.25 * JD) +JD2S = 86400 +year = np.longdouble(365.25 * JD2S) einsteinC = np.longdouble(299792458.0) +# Solar oblatenes values: From Mecheri et al. (2004), using Corbard (b) 2002 values (Table II) +J2Sun = np.longdouble(2.198e-7) +J4Sun = np.longdouble(-4.805e-9) + #I/O Routines for reading in Swifter and Swiftest parameter and binary data files def read_swifter_param(inparfile): @@ -127,6 +132,7 @@ def read_swiftest_config(config_file_name): 'DT' : 0.0, 'PL_IN' : "", 'TP_IN' : "", + 'CB_IN' : "", 'IN_TYPE' : "ASCII", 'ISTEP_OUT' : -1, 'BIN_OUT' : "", @@ -297,6 +303,47 @@ def swifter_stream(f, param): yield t, npl, plid, pvec.T, plab, \ ntp, tpid, tvec.T, tlab + +def make_swiftest_labels(config): + tlab = [] + if config['OUT_FORM'] == 'XV': + tlab.append('px') + tlab.append('py') + tlab.append('pz') + tlab.append('vx') + tlab.append('vy') + tlab.append('vz') + elif config['OUT_FORM'] == 'EL': + tlab.append('a') + tlab.append('e') + tlab.append('inc') + tlab.append('capom') + tlab.append('omega') + tlab.append('capm') + plab = tlab.copy() + plab.append('Mass') + plab.append('Radius') + clab = ['Mass', 'Radius', 'J_2', 'J_4'] + if config['ROTATION'] == 'YES': + clab.append('Ip_x') + clab.append('Ip_y') + clab.append('Ip_z') + clab.append('rot_x') + clab.append('rot_y') + clab.append('rot_z') + plab.append('Ip_x') + plab.append('Ip_y') + plab.append('Ip_z') + plab.append('rot_x') + plab.append('rot_y') + plab.append('rot_z') + if config['TIDES'] == 'YES': + clab.append('k2') + clab.append('Q') + plab.append('k2') + plab.append('Q') + return clab, plab, tlab + def swiftest_stream(f, config): """ Reads in a Swifter bin.dat file and returns a single frame of data as a datastream @@ -384,43 +431,7 @@ def swiftest_stream(f, config): t6 = f.read_reals(np.float64) cbid = np.array([0]) - tlab = [] - if config['OUT_FORM'] == 'XV': - tlab.append('px') - tlab.append('py') - tlab.append('pz') - tlab.append('vx') - tlab.append('vy') - tlab.append('vz') - elif config['OUT_FORM'] == 'EL': - tlab.append('a') - tlab.append('e') - tlab.append('inc') - tlab.append('capom') - tlab.append('omega') - tlab.append('capm') - plab = tlab.copy() - plab.append('Mass') - plab.append('Radius') - clab = ['Mass', 'Radius', 'J_2', 'J_4'] - if config['ROTATION'] == 'YES': - clab.append('Ip_x') - clab.append('Ip_y') - clab.append('Ip_z') - clab.append('rot_x') - clab.append('rot_y') - clab.append('rot_z') - plab.append('Ip_x') - plab.append('Ip_y') - plab.append('Ip_z') - plab.append('rot_x') - plab.append('rot_y') - plab.append('rot_z') - if config['TIDES'] == 'YES': - clab.append('k2') - clab.append('Q') - plab.append('k2') - plab.append('Q') + clab, plab, tlab = make_swiftest_labels(config) if npl > 0: pvec = np.vstack([p1,p2,p3,p4,p5,p6,Mpl,Rpl]) @@ -519,18 +530,253 @@ def swiftest2xr(config): ds = xr.combine_by_coords([cbds, plds, tpds]) return ds +def solar_system_pl(config, ephemerides_start_date): + """ + Initializes a Swiftest dataset containing the major planets of the Solar System at a particular data from JPL/Horizons + + Parameters + ---------- + config : dict + Swiftest Configuration parameters. This method uses the unit conversion factors to convert from JPL's AU-day system into the system specified in the config file + ephemerides_start_date : string + Date to use when obtaining the ephemerides in the format YYYY-MM-DD + + Returns + ------- + xarray dataset + """ + # Planet ids + planetid = { + 'mercury': '1', + 'venus': '2', + 'earthmoon': '3', + 'mars': '4', + 'jupiter': '5', + 'saturn': '6', + 'uranus': '7', + 'neptune': '8', + 'plutocharon': '9' + } + + # Planet Msun/M ratio + MSun_over_Mpl = { + 'mercury': np.longdouble(6023600.0), + 'venus': np.longdouble(408523.71), + 'earthmoon': np.longdouble(328900.56), + 'mars': np.longdouble(3098708.), + 'jupiter': np.longdouble(1047.3486), + 'saturn': np.longdouble(3497.898), + 'uranus': np.longdouble(22902.98), + 'neptune': np.longdouble(19412.24), + 'plutocharon': np.longdouble(1.35e8) + } + + # Planet radii in meters + planetradius = { + 'mercury': np.longdouble(2439.4e3), + 'venus': np.longdouble(6051.8e3), + 'earthmoon': np.longdouble(6371.0084e3), # Earth only for radius + 'mars': np.longdouble(3389.50e3), + 'jupiter': np.longdouble(69911e3), + 'saturn': np.longdouble(58232.0e3), + 'uranus': np.longdouble(25362.e3), + 'neptune': np.longdouble(24622.e3), + 'plutocharon': np.longdouble(1188.3e3) + } + + # Unit conversion factors + DCONV = AU2M / config['DU2M'] + VCONV = (AU2M / JD2S) / (config['DU2M'] / config['TU2S']) + THIRDLONG = np.longdouble(1.0) / np.longdouble(3.0) + + # Central body value vectors + GMcb = np.array([GMSunSI * config['TU2S']**2 / config['DU2M']**3]) + Rcb = np.array([RSun / config['DU2M']]) + J2RP2 = np.array([J2Sun * (RSun / config['DU2M'])**2]) + J4RP4 = np.array([J4Sun * (RSun / config['DU2M'])**4]) + cbid = np.array([0]) + cvec = np.vstack([GMcb, Rcb, J2RP2, J4RP4]) + + # Horizons date time internal variables + tstart = datetime.date.fromisoformat(ephemerides_start_date) + tstep = datetime.timedelta(days=1) + tend = tstart + tstep + ephemerides_end_date = tend.isoformat() + ephemerides_step = '1d' + + pldata = {} + p1 = [] + p2 = [] + p3 = [] + p4 = [] + p5 = [] + p6 = [] + Rhill = [] + Rpl = [] + GMpl = [] + + # Fetch solar system ephemerides from Horizons + for key, val in planetid.items(): + pldata[key] = Horizons(id=val, id_type='majorbody', location='@sun', + epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, + 'step': ephemerides_step}) + if config['OUT_FORM'] == 'XV': + p1.append(pldata[key].vectors()['x'][0] * DCONV) + p2.append(pldata[key].vectors()['y'][0] * DCONV) + p3.append(pldata[key].vectors()['z'][0] * DCONV) + p4.append(pldata[key].vectors()['vx'][0] * VCONV) + p5.append(pldata[key].vectors()['vy'][0] * VCONV) + p6.append(pldata[key].vectors()['vz'][0] * VCONV) + elif config['OUT_FORM'] == 'EL': + p1.append(pldata[key].elements()['a'][0] * DCONV) + p2.append(pldata[key].elements()['e'][0]) + p3.append(pldata[key].elements()['inc'][0] * np.pi / 180.0) + p4.append(pldata[key].elements()['Omega'][0] * np.pi / 180.0) + p5.append(pldata[key].elements()['w'][0] * np.pi / 180.0) + p6.append(pldata[key].elements()['M'][0] * np.pi / 180.0) + Rhill.append(pldata[key].elements()['a'][0] * (3 * MSun_over_Mpl[key]) ** (-THIRDLONG)) + Rpl.append(planetradius[key] * DCONV) + GMpl.append(GMcb[0] / MSun_over_Mpl[key]) + # Generate planet value vectors + plid = np.fromiter(planetid.values(), dtype=int) + pvec = np.vstack([p1, p2, p3, p4, p5, p6, GMpl, Rpl]) + + dims = ['time', 'id', 'vec'] + cb = [] + pl = [] + tp = [] + t = np.array([0.0]) + + clab, plab, tlab = make_swiftest_labels(config) + + #Prepare frames by adding an extra axis for the time coordinate + cbframe = np.expand_dims(cvec.T, axis=0) + plframe = np.expand_dims(pvec.T, axis=0) + + #Create xarray DataArrays out of each body type + cbxr = xr.DataArray(cbframe, dims = dims, coords = {'time' : t, 'id' : cbid, 'vec' : clab}) + plxr = xr.DataArray(plframe, dims = dims, coords = {'time' : t, 'id' : plid, 'vec' : plab}) + + cb.append(cbxr) + pl.append(plxr) + + cbda = xr.concat(cb,dim='time') + plda = xr.concat(pl,dim='time') + + cbds = cbda.to_dataset(dim = 'vec') + plds = plda.to_dataset(dim = 'vec') + ds = xr.combine_by_coords([cbds, plds]) + return ds + +def swiftest_xr2_infile(ds, config, framenum=-1): + """ + Writes a set of Swiftest input files from a single frame of a Swiftest xarray dataset + + Parameters + ---------- + ds : xarray dataset + Dataset containing Swiftest n-body data in XV format + framenum : int + Time frame to use to generate the initial conditions. If this argument is not passed, the default is to use the last frame in the dataset. + config : dict + Swiftest Configuration parameters. This method uses the names of the cb, pl, and tp files from the configuration + + Returns + ------- + A set of three input files for a Swiftest run + """ + frame = ds.isel(time=framenum) + cb = frame.where(frame.id == 0, drop=True) + pl = frame.where(frame.id > 0, drop=True) + pl = pl.where(np.invert(np.isnan(pl['Mass'])), drop=True).drop_vars(['J_2', 'J_4']) + tp = frame.where(np.isnan(frame['Mass']), drop=True).drop_vars(['Mass', 'Radius', 'J_2', 'J_4']) + + GMSun = np.double(cb['Mass']) + RSun = np.double(cb['Radius']) + J2 = np.double(cb['J_2']) + J4 = np.double(cb['J_4']) + + if config['IN_TYPE'] == 'ASCII': + # Swiftest Central body file + cbfile = open(config['CB_IN'], 'w') + print(GMSun, file=cbfile) + print(RSun, file=cbfile) + print(J2, file=cbfile) + print(J4, file=cbfile) + cbfile.close() + + plfile = open(config['PL_IN'], 'w') + print(pl.id.count().values, file=plfile) + for i in pl.id: + pli = pl.sel(id=i) + print(i.values, pli['Mass'].values, file=plfile) + print(pli['Radius'].values, file=plfile) + print(pli['px'].values, pli['py'].values, pli['pz'].values, file=plfile) + print(pli['vx'].values, pli['vy'].values, pli['vz'].values, file=plfile) + plfile.close() + + # TP file + tpfile = open(config['TP_IN'], 'w') + print(tp.id.count().values, file=tpfile) + for i in tp.id: + tpi = tp.sel(id=i) + print(i.values, file=tpfile) + print(tpi['px'].values, tpi['py'].values, tpi['pz'].values, file=tpfile) + print(tpi['vx'].values, tpi['vy'].values, tpi['vz'].values, file=tpfile) + tpfile.close() + elif config['IN_TYPE'] == 'REAL8': + # Now make Swiftest files + cbfile = FortranFile(swiftest_cb, 'w') + Msun = np.double(1.0) + cbfile.write_record(np.double(GMSun)) + cbfile.write_record(np.double(rmin)) + cbfile.write_record(np.double(J2)) + cbfile.write_record(np.double(J4)) + cbfile.close() + + plfile = FortranFile(swiftest_pl, 'w') + plfile.write_record(npl) + + plfile.write_record(plid) + plfile.write_record(p_pl[0]) + plfile.write_record(p_pl[1]) + plfile.write_record(p_pl[2]) + plfile.write_record(v_pl[0]) + plfile.write_record(v_pl[1]) + plfile.write_record(v_pl[2]) + plfile.write_record(mass) + plfile.write_record(radius) + plfile.close() + tpfile = FortranFile(swiftest_tp, 'w') + ntp = 1 + tpfile.write_record(ntp) + tpfile.write_record(tpid) + tpfile.write_record(p_tp[0]) + tpfile.write_record(p_tp[1]) + tpfile.write_record(p_tp[2]) + tpfile.write_record(v_tp[0]) + tpfile.write_record(v_tp[1]) + tpfile.write_record(v_tp[2]) + else: + print(f"{config['IN_TYPE']} is an unknown file type") if __name__ == '__main__': - workingdir = '/Users/daminton/git/swiftest/examples/rmvs_swifter_comparison/mars_ejecta/' - inparfile = workingdir + 'param.swifter.in' - param = read_swifter_param(inparfile) - param['BIN_OUT'] = workingdir + param['BIN_OUT'] + workingdir = '/Users/daminton/git/swiftest/examples/rmvs_swifter_comparison/9pl_18tp_encounters/' + #inparfile = workingdir + 'param.swifter.in' + #param = read_swifter_param(inparfile) + #param['BIN_OUT'] = workingdir + param['BIN_OUT'] config_file_name = workingdir + 'config.swiftest.in' config = read_swiftest_config(config_file_name) config['BIN_OUT'] = workingdir + config['BIN_OUT'] - - swiftestdat = swiftest2xr(config) + ds = solar_system_pl(config, '2020-06-17') + ds + + #swiftestdat = swiftest2xr(config) + config['CB_IN'] = workingdir + 'cb_test.in' + config['PL_IN'] = workingdir + 'pl_test.in' + config['TP_IN'] = workingdir + 'tp_test.in' + swiftest_xr2_infile(ds, config) #swifterdat = swifter2xr(param) #print(swiftestdat['a']) #print(swiftestdf.head())