Fixed initial conditions generator for RMVS comparison run (but now t…

…he test particles are not the same

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

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

examples/rmvs_swifter_comparison/9pl_18tp_encounters/init_cond.py

@@ -20,67 +20,27 @@
 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'RHILL_PRESENT  no')
-print(f'GR             no')
-print(f'MU2KG          {MU2KG}')
-print(f'DU2M           {DU2M}')
-print(f'TU2S           {TU2S}')
-sys.stdout = sys.__stdout__
-sim = swiftest.Simulation(param_file=swiftest_input)
-param = sim.param
-
-# Dates to fetch planet ephemerides from JPL Horizons
-tstart = '2021-06-15'
+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,
@@ -106,7 +66,7 @@
 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)
@@ -155,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            {swiftest.J2Sun}")
-print(f"J4            {swiftest.J4Sun}")
-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 @@
-! VERSION     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            2.198e-07
-J4            -4.805e-09
-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,30 +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
-RHILL_PRESENT  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.swifter.in

@@ -2,35 +2,35 @@
 0 0.00029591220819207775568
 0.0 0.0 0.0
 0.0 0.0 0.0
-1 4.9125474498983625e-11 0.0014751244996981092
-1.6306381826061646e-05
-0.3424231950105548 0.04541859803875346 -0.02769885218647143
--0.009135456552973951 0.02914778346055028 0.003219860945895814
-2 7.243452483873647e-10 0.00675910592748785
-4.0453784346544176e-05
--0.7187409324847692 0.008460121607685698 0.04159114061029223
--0.000355507214371255 -0.02031492463311098 -0.0002582876464863964
-3 8.997011382166019e-10 0.010044781409779763
-4.25875607065041e-05
-0.3408913444637776 -0.9577324614296441 4.435598231277718e-05
-0.01592829229600895 0.005704842145070232 -3.222239146362505e-07
-4 9.549535102761465e-11 0.007246745952808949
-2.2657408050928896e-05
--1.518194954784943 0.6839686166809457 0.05157497542651887
--0.005224888295856506 -0.01156432042916445 -0.0001141828894241638
-5 2.825345908631355e-07 0.3552712944567904
-0.0004673261703049093
-4.04554302433205 -2.997516585076586 -0.07806209765007878
-0.004406596292841761 0.006425243736188545 -0.0001252737469736559
-6 8.459715183006416e-08 0.43765252895139073
-0.00038925687730393614
-6.294914485420795 -7.709936089816689 -0.1164782091116569
-0.004015969235853927 0.003521900151979878 -0.000221006899984545
-7 1.2920249163736674e-08 0.46953353749301263
-0.00016953449859497232
-14.85869768503949 13.00480689249104 -0.1442220329069338
--0.002615247360000599 0.002782629142130045 4.407243103045647e-05
-8 1.5243589003230834e-08 0.7812837944287938
-0.00016458790412449367
-29.55697963607957 -4.6325049341728 -0.5858344271811813
-0.0004702098511244648 0.0031273464291932 -7.514820514613305e-05
+1 4.9125474498983623693e-11 0.0014751244996981091688
+1.6306381826061645943e-05
+0.3424231950105547928 0.045418598038753463242 -0.027698852186471431547
+-0.009135456552973951483 0.029147783460550278495 0.003219860945895814102
+2 7.243452483873646905e-10 0.006759105927487850036
+4.0453784346544178454e-05
+-0.7187409324847692238 0.008460121607685697903 0.041591140610292232083
+-0.00035550721437125498313 -0.020314924633110978403 -0.00025828764648639637377
+3 8.9970113821660187435e-10 0.010044781409779763954
+4.25875607065040958e-05
+0.34089134446377761245 -0.95773246142964407746 4.435598231277717939e-05
+0.015928292296008950135 0.005704842145070231768 -3.222239146362504855e-07
+4 9.549535102761465607e-11 0.007246745952808948377
+2.265740805092889601e-05
+-1.5181949547849429294 0.68396861668094566244 0.051574975426518870902
+-0.0052248882958565064094 -0.011564320429164449966 -0.0001141828894241637938
+5 2.825345908631354893e-07 0.35527129445679039879
+0.00046732617030490929307
+4.0455430243320495975 -2.9975165850765859155 -0.07806209765007877943
+0.0044065962928417608604 0.006425243736188544774 -0.00012527374697365590813
+6 8.459715183006415395e-08 0.43765252895139074356
+0.00038925687730393611812
+6.29491448542079457 -7.7099360898166890976 -0.11647820911165690516
+0.004015969235853926976 0.0035219001519798780186 -0.00022100689998454499602
+7 1.2920249163736673626e-08 0.4695335374930126262
+0.00016953449859497231466
+14.8586976850394894 13.004806892491039605 -0.14422203290693380584
+-0.002615247360000599007 0.0027826291421300451516 4.4072431030456472162e-05
+8 1.5243589003230834323e-08 0.78128379442879379807
+0.000164587904124493665
+29.55697963607957135 -4.6325049341728004038 -0.5858344271811812831
+0.00047020985112446482199 0.0031273464291932001093 -7.514820514613305166e-05