Fixed initial conditions on 9pl_18tp RMVS test case to track down bug…

… introduced on OOPSymba branch

examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb.in

@@ -0,0 +1,4 @@
+0.00029591220819207774
+0.004650467260962157
+4.7535806948127355e-12
+-2.2473967953572827e-18

examples/rmvs_swifter_comparison/9pl_18tp_encounters/cb.swiftest.in

@@ -1,4 +1,4 @@
-0.0002959122081920778
+0.00029591220819207774
 0.004650467260962157
 4.7535806948127355e-12
 -2.2473967953572827e-18

examples/rmvs_swifter_comparison/9pl_18tp_encounters/init_cond.py

@@ -1,4 +1,6 @@
+#!/usr/bin/env python3
 import numpy as np
+import swiftest
 import swiftest.io as swio
 import astropy.constants as const
 import sys
@@ -18,75 +20,53 @@
 swiftest_bin   = "bin.swiftest.dat"
 swiftest_enc   = "enc.swiftest.dat"
 
-#Values from JPL Horizons
-AU2M = np.longdouble(const.au.value)
-GMSunSI = np.longdouble(const.GM_sun.value)
-Rsun = np.longdouble(const.R_sun.value)
-GC = np.longdouble(const.G.value)
-JD = 86400
-year = np.longdouble(365.25 * JD)
-c = np.longdouble(299792458.0)
-
-MU2KG    = np.longdouble(GMSunSI / GC) #Conversion from mass unit to kg
-DU2M     = np.longdouble(AU2M)         #Conversion from radius unit to centimeters
-TU2S     = np.longdouble(JD)         #Conversion from time unit to seconds
-GU       = np.longdouble(GC / (DU2M**3 / (MU2KG * TU2S**2)))
-GMSun = np.longdouble(GMSunSI / (DU2M**3 / TU2S**2))
-
-t_0	  = 0 # simulation start time
-deltaT	= 1.00 * JD / TU2S   # simulation step size
-end_sim = year / TU2S # simulation end time
-t_print = deltaT  #output interval to print results
-
-iout = int(np.ceil(t_print / deltaT))
-rmin = Rsun / DU2M
-rmax = 1000.0
-
-sys.stdout = open(swiftest_input, "w")
-print(f'! VERSION      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          {tpin}')
-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       XV')
-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             no')
-print(f'MU2KG          {MU2KG}')
-print(f'DU2M           {DU2M}')
-print(f'TU2S           {TU2S}')
-sys.stdout = sys.__stdout__
-param = swio.read_swiftest_param(swiftest_input)
-
-# Dates to fetch planet ephemerides from JPL Horizons
-tstart = '2021-06-15'
-ds = swio.solar_system_pl(param, tstart)
+sim = swiftest.Simulation()
+
+sim.param['T0'] = 0.0
+sim.param['DT'] = 1.0 
+sim.param['TSTOP'] = 365.25
+sim.param['ISTEP_OUT']  = 11
+sim.param['ISTEP_DUMP'] = 1
+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'] = "XV"
+sim.param['OUT_STAT'] = "UNKNOWN"
+sim.param['GR'] = 'NO'
+sim.param['CHK_CLOSE'] = 'YES'
+
+sim.param['MU2KG'] = swiftest.MSun
+sim.param['TU2S'] = swiftest.JD2S
+sim.param['DU2M'] = swiftest.AU2M
+
+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)
+
+ds = sim.ds
 cb = ds.sel(id=0)
 pl = ds.where(ds.id > 0, drop=True)
 npl = pl.id.size
 
-ntp = 18
+ntp = 16
 dims = ['time', 'id', 'vec']
 tp = []
 t = np.array([0.0])
-clab, plab, tlab = swio.make_swiftest_labels(param)
+clab, plab, tlab = swio.make_swiftest_labels(sim.param)
 
 # 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)
@@ -135,50 +115,18 @@
 tpda = xr.concat(tp,dim='time')
 tpds = tpda.to_dataset(dim = 'vec')
 
-ds = xr.combine_by_coords([ds, tpds])
-swio.swiftest_xr2infile(ds, param)
-
-# Swifter PL file
-plfile = open(swifter_pl, 'w')
-print(npl + 1, file=plfile)
-print(0,GMSun,file=plfile)
-print('0.0 0.0 0.0',file=plfile)
-print('0.0 0.0 0.0',file=plfile)
-for i in pl.id:
-    pli = pl.sel(id=i)
-    print(f"{int(i)} {pli['Mass'].values[0]} {pli['Rhill'].values[0]}", file=plfile)
-    print(f"{pli['Radius'].values[0]}", file=plfile)
-    print(f"{pli['px'].values[0]} {pli['py'].values[0]} {pli['pz'].values[0]}", file=plfile)
-    print(f"{pli['vx'].values[0]} {pli['vy'].values[0]} {pli['vz'].values[0]}", file=plfile)
-plfile.close()
-
-# Swifter parameter file
-sys.stdout = open(swifter_input, "w")
-print(f"! VERSION     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      UNKNOWN")
-print(f"J2            {param['J2']}")
-print(f"J4            {param['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__
-
+sim.ds = xr.combine_by_coords([sim.ds, tpds])
+swio.swiftest_xr2infile(sim.ds, sim.param)
+
+sim.param['PL_IN'] = swiftest_pl
+sim.param['TP_IN'] = tpin
+sim.param['CB_IN'] = swiftest_cb
+sim.param['BIN_OUT'] = swiftest_bin
+sim.param['ENC_OUT'] = swiftest_enc
+sim.save(swiftest_input)
+
+sim.param['PL_IN'] = swifter_pl
+sim.param['TP_IN'] = tpin
+sim.param['BIN_OUT'] = swifter_bin
+sim.param['ENC_OUT'] = swifter_enc
+sim.save(swifter_input, codename="Swifter")

examples/rmvs_swifter_comparison/9pl_18tp_encounters/param.swifter.in

@@ -1,26 +1,26 @@
-! Swifter input file generated using init_cond.py
-T0            0 
-TSTOP         365.25
-DT            1.0
-PL_IN         pl.swifter.in
-TP_IN         tp.in
-IN_TYPE       ASCII
-ISTEP_OUT     1
-ISTEP_DUMP    1
-BIN_OUT       bin.swifter.dat
-OUT_TYPE      REAL8
-OUT_FORM      XV
-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
+! VERSION        Swifter parameter file converted from Swiftest
+T0               0.0
+TSTOP            365.25
+DT               1.0
+ISTEP_OUT        11
+ISTEP_DUMP       1
+OUT_FORM         XV
+OUT_TYPE         REAL8
+OUT_STAT         UNKNOWN
+IN_TYPE          ASCII
+PL_IN            pl.swifter.in
+TP_IN            tp.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
+J2               4.7535806948127355e-12
+J4               -2.2473967953572827e-18
+RHILL_PRESENT    YES

examples/rmvs_swifter_comparison/9pl_18tp_encounters/param.swiftest.in

@@ -1,29 +1,35 @@
-! VERSION      Swiftest input file generated using init_cond.py
-T0             0 
-TSTOP          365.25
-DT             1.0
-CB_IN          cb.swiftest.in
-PL_IN          pl.swiftest.in
-TP_IN          tp.in
-IN_TYPE        ASCII
-ISTEP_OUT      1
-ISTEP_DUMP     1
-BIN_OUT        bin.swiftest.dat
-OUT_TYPE       REAL8
-OUT_FORM       XV
-OUT_STAT       REPLACE
-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
-GR             no
-MU2KG          1.988409870698051e+30
-DU2M           149597870700.0
-TU2S           86400.0
+! VERSION        Swiftest parameter input
+T0               0.0
+TSTOP            365.25
+DT               1.0
+ISTEP_OUT        11
+ISTEP_DUMP       1
+OUT_FORM         XV
+OUT_TYPE         REAL8
+OUT_STAT         UNKNOWN
+IN_TYPE          ASCII
+PL_IN            pl.swiftest.in
+TP_IN            tp.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             86400
+DU2M             149597870700.0
+EXTRA_FORCE      NO
+BIG_DISCARD      NO
+CHK_CLOSE        YES
+FRAGMENTATION    NO
+ROTATION         NO
+TIDES            NO
+ENERGY           NO
+GR               NO
+YARKOVSKY        NO
+YORP             NO
+MTINY            0.0

examples/rmvs_swifter_comparison/9pl_18tp_encounters/pl.in

@@ -0,0 +1,33 @@
+8
+1 4.9125474498983623693e-11
+1.6306381826061645943e-05
+0.33206272695596028566 0.07436707001147663254 -0.02438290851908785084
+-0.0115920916602103591525 0.028710618792657981169 0.0034094833969203438596
+2 7.243452483873646905e-10
+4.0453784346544178454e-05
+-0.7188115337296047125 -0.0118554711069603201795 0.041316403191083782287
+0.00021427347881133320621 -0.020313576971905909774 -0.00029114855617710840843
+3 8.9970113821660187435e-10
+4.25875607065040958e-05
+0.35677088372527121507 -0.95189300879814897627 4.4027442504036787155e-05
+0.015830039028334789986 0.0059737936889703449964 -3.3484113013969089573e-07
+4 9.549535102761465607e-11
+2.265740805092889601e-05
+-1.5233712071242269115 0.6723825347339112968 0.051459143378398922164
+-0.0051275613251079554117 -0.011607719813367209372 -0.000117479966462153095864
+5 2.825345908631354893e-07
+0.00046732617030490929307
+4.049944927347420176 -2.9910878677758190314 -0.078187280837353656526
+0.0043972077687938898594 0.006432188574295680597 -0.00012509257442073270106
+6 8.459715183006415395e-08
+0.00038925687730393611812
+6.298929503477405767 -7.706413024510769816 -0.11669919842191249504
+0.0040140666547768266703 0.0035242303011843410798 -0.00022097170940726839814
+7 1.2920249163736673626e-08
+0.00016953449859497231466
+14.856082147529010129 13.007589275314199284 -0.14417795763685259391
+-0.0026158276515510360365 0.0027821364817078499815 4.40781085949555924e-05
+8 1.5243589003230834323e-08
+0.000164587904124493665
+29.55744967800954015 -4.629377558152945049 -0.58590957207831262377
+0.00046987400245862169295 0.0031274056019462009859 -7.51415892482447254e-05