More maintenence on initial conditions generator and fixes to i/o pro…

…blems

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

@@ -1,4 +1,5 @@
 0.0002959122081920778
+0.0002959122081920778
 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
@@ -68,21 +70,39 @@
 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__
-param = swio.read_swiftest_param(swiftest_input)
+sim = swiftest.Simulation(param_file=swiftest_input)
+param = sim.param
 
 # Dates to fetch planet ephemerides from JPL Horizons
 tstart = '2021-06-15'
-ds = swio.solar_system_pl(param, tstart)
+
+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])
@@ -167,8 +187,8 @@
 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"J2            {swiftest.J2Sun}")
+print(f"J4            {swiftest.J4Sun}")
 print(f"CHK_CLOSE     yes")
 print(f"CHK_RMIN      {rmin}")
 print(f"CHK_RMAX      {rmax}")

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

@@ -1,4 +1,4 @@
-! Swifter input file generated using init_cond.py
+! VERSION     Swifter input file generated using init_cond.py
 T0            0 
 TSTOP         365.25
 DT            1.0
@@ -11,8 +11,8 @@ BIN_OUT       bin.swifter.dat
 OUT_TYPE      REAL8
 OUT_FORM      XV
 OUT_STAT      UNKNOWN
-J2            4.7535806948127355e-12
-J4            -2.2473967953572827e-18
+J2            2.198e-07
+J4            -4.805e-09
 CHK_CLOSE     yes
 CHK_RMIN      0.004650467260962157
 CHK_RMAX      1000.0

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

@@ -23,6 +23,7 @@ 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

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

@@ -1,40 +1,36 @@
-10
-0 0.00029591220819207776388
+9
+0 0.00029591220819207775568
 0.0 0.0 0.0
 0.0 0.0 0.0
-1 4.9125474498983625e-11 0.0014751258227142052
+1 4.9125474498983625e-11 0.0014751244996981092
 1.6306381826061646e-05
-0.008059842448018334 -0.4616051037329109 -0.03846017738329229
-0.02248719132054853 0.001934639213990692 -0.001904656977422976
-2 7.243452483873647e-10 0.006759134232034941
+0.3424231950105548 0.04541859803875346 -0.02769885218647143
+-0.009135456552973951 0.02914778346055028 0.003219860945895814
+2 7.243452483873647e-10 0.00675910592748785
 4.0453784346544176e-05
--0.5115875215389065 0.5030818749037324 0.03642547299277956
--0.01425515725454357 -0.01452868630179309 0.0006232072038298823
-3 8.997011382166019e-10 0.010044625087011915
+-0.7187409324847692 0.008460121607685698 0.04159114061029223
+-0.000355507214371255 -0.02031492463311098 -0.0002582876464863964
+3 8.997011382166019e-10 0.010044781409779763
 4.25875607065041e-05
--0.1090020607540907 -1.009893805009766 4.823302918632528e-05
-0.01682491922568941 -0.001910549762056979 3.992660742687128e-08
-4 9.549535102761465e-11 0.0072467897902424765
+0.3408913444637776 -0.9577324614296441 4.435598231277718e-05
+0.01592829229600895 0.005704842145070232 -3.222239146362505e-07
+4 9.549535102761465e-11 0.007246745952808949
 2.2657408050928896e-05
--1.342897929331636 0.9778655112682739 0.05343398538723887
--0.007712315645393206 -0.01011917844182223 -2.287744801261131e-05
-5 2.825345908631355e-07 0.3552720805286442
+-1.518194954784943 0.6839686166809457 0.05157497542651887
+-0.005224888295856506 -0.01156432042916445 -0.0001141828894241638
+5 2.825345908631355e-07 0.3552712944567904
 0.0004673261703049093
-3.923184193414315 -3.168419770483168 -0.0746147877972047
-0.004655552638985802 0.006232623300954468 -0.0001300429201057457
-6 8.459715183006416e-08 0.4376460836930155
+4.04554302433205 -2.997516585076586 -0.07806209765007878
+0.004406596292841761 0.006425243736188545 -0.0001252737469736559
+6 8.459715183006416e-08 0.43765252895139073
 0.00038925687730393614
-6.185794462795267 -7.804174837804826 -0.110498432926239
-0.004066833203985018 0.003458637040736611 -0.0002219310939327014
-7 1.2920249163736674e-08 0.46946272948265794
+6.294914485420795 -7.709936089816689 -0.1164782091116569
+0.004015969235853927 0.003521900151979878 -0.000221006899984545
+7 1.2920249163736674e-08 0.46953353749301263
 0.00016953449859497232
-14.9290976575471 12.92949673572929 -0.1454099139559955
--0.002599557960646664 0.002795888198858545 4.391864857782088e-05
-8 1.5243589003230834e-08 0.78119478483336
+14.85869768503949 13.00480689249104 -0.1442220329069338
+-0.002615247360000599 0.002782629142130045 4.407243103045647e-05
+8 1.5243589003230834e-08 0.7812837944287938
 0.00016458790412449367
-29.54416169025338 -4.716921603714237 -0.5838030174427992
-0.0004792636209523189 0.00312573757291745 -7.53264045199501e-05
-9 2.1919422829042796e-12 0.05379680851617536
-7.943294877391593e-06
-14.54448346259197 -31.05223519593471 -0.8828000265625595
-0.002923077617691739 0.0006625916902153526 -0.0009142553677224461
+29.55697963607957 -4.6325049341728 -0.5858344271811813
+0.0004702098511244648 0.0031273464291932 -7.514820514613305e-05

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

@@ -1,37 +1,33 @@
-9
+8
 1 4.9125474498983625056e-11
 1.6306381826061645943e-05
-0.008059842448018333591 -0.46160510373291091524 -0.038460177383292291908
-0.02248719132054852949 0.0019346392139906921279 -0.0019046569774229759606
+0.3424231950105547928 0.045418598038753463242 -0.027698852186471431547
+-0.0091354565529739517565 0.029147783460550279367 0.0032198609458958141984
 2 7.243452483873647106e-10
 4.0453784346544178454e-05
--0.51158752153890652004 0.5030818749037323512 0.036425472992779560355
--0.01425515725454356988 -0.014528686301793089943 0.00062320720382988232425
+-0.7187409324847692238 0.008460121607685697903 0.041591140610292232083
+-0.00035550721437125499377 -0.02031492463311097901 -0.0002582876464863963815
 3 8.997011382166018993e-10
 4.25875607065040958e-05
--0.109002060754090704386 -1.0098938050097661101 4.8233029186325282966e-05
-0.016824919225689409317 -0.0019105497620569790936 3.9926607426871282392e-08
+0.34089134446377761245 -0.95773246142964407746 4.435598231277717939e-05
+0.015928292296008950611 0.0057048421450702319388 -3.2222391463625049515e-07
 4 9.549535102761465872e-11
 2.265740805092889601e-05
--1.3428979293316360977 0.97786551126827392366 0.053433985387238869258
--0.007712315645393206092 -0.0101191784418222296525 -2.2877448012611311785e-05
+-1.5181949547849429294 0.68396861668094566244 0.051574975426518870902
+-0.0052248882958565065657 -0.011564320429164450312 -0.000114182889424163797216
 5 2.8253459086313549713e-07
 0.00046732617030490929307
-3.923184193414314791 -3.1684197704831680298 -0.07461478779720470689
-0.0046555526389858022113 0.006232623300954467766 -0.00013004292010574569365
+4.0455430243320495975 -2.9975165850765859155 -0.07806209765007877943
+0.004406596292841760992 0.0064252437361885449665 -0.00012527374697365591187
 6 8.45971518300641563e-08
 0.00038925687730393611812
-6.1857944627952665684 -7.804174837804826126 -0.11049843292623899582
-0.0040668332039850179674 0.0034586370407366113193 -0.00022193109393270141328
+6.29491448542079457 -7.7099360898166890976 -0.11647820911165690516
+0.0040159692358539270963 0.003521900151979878124 -0.00022100689998454500264
 7 1.2920249163736673984e-08
 0.00016953449859497231466
-14.929097657547099942 12.9294967357292893695 -0.14540991395599550673
--0.0025995579606466640267 0.0027958881988585450113 4.391864857782088156e-05
+14.8586976850394894 13.004806892491039605 -0.14422203290693380584
+-0.0026152473600005990854 0.0027826291421300452348 4.4072431030456473482e-05
 8 1.5243589003230834746e-08
 0.000164587904124493665
-29.544161690253378794 -4.7169216037142369657 -0.58380301744279916587
-0.00047926362095231893815 0.00312573757291745008 -7.532640451995010825e-05
-9 2.1919422829042797324e-12
-7.943294877391593783e-06
-14.544483462591969669 -31.052235195934709822 -0.88280002656255951443
-0.0029230776176917390448 0.0006625916902153525834 -0.0009142553677224461557
+29.55697963607957135 -4.6325049341728004038 -0.5858344271811812831
+0.00047020985112446483607 0.0031273464291932002029 -7.514820514613305391e-05

examples/rmvs_swifter_comparison/9pl_18tp_encounters/tp.in

@@ -1,55 +1,49 @@
-18
+16
 101
-0.0080829031543499848395 -0.46158204302657929174 -0.038460177383292291908
-0.025090740792557679473 0.00193463921399069207 -0.0019046569774229759036
+0.3424462557168864163 0.04544165874508511449 -0.027698852186471431547
+-0.006531907080964799092 0.029147783460550278495 0.003219860945895814102
 102
-0.008036781741686682343 -0.46162816443924253873 -0.038460177383292291908
-0.01988364184853937816 0.00193463921399069207 -0.0019046569774229759036
+0.3424001343042231693 0.045395537332421811993 -0.027698852186471431547
+-0.0117390060249831038736 0.029147783460550278495 0.003219860945895814102
 103
--0.51153031124843428845 0.50313908519420458276 0.036425472992779560355
--0.007907923549198175514 -0.014528686301793089508 0.0006232072038298823056
+-0.7186837221942969922 0.0085173318981578965275 0.041591140610292232083
+0.0059917264909741391188 -0.020314924633110978403 -0.00025828764648639637377
 104
--0.5116447318293787516 0.5030246646132601196 0.036425472992779560355
--0.020602390959888965127 -0.014528686301793089508 0.0006232072038298823056
+-0.7187981427752414554 0.008402911317213499279 0.041591140610292232083
+-0.0067027409197166487598 -0.020314924633110978403 -0.00025828764648639637377
 105
--0.10894183284815117663 -1.0098335771038264852 4.8233029186325282966e-05
-0.023719359459963285097 -0.0019105497620569790364 3.9926607426871281197e-08
+0.34095157236971712633 -0.9576722335237045636 4.435598231277717939e-05
+0.022822732530282826419 0.005704842145070231768 -3.222239146362504855e-07
 106
--0.109062288660030232146 -1.009954032915705735 4.8233029186325282966e-05
-0.00993047899141553253 -0.0019105497620569790364 3.9926607426871281197e-08
+0.34083111655783809857 -0.95779268933558359134 4.435598231277717939e-05
+0.009033852061735073852 0.005704842145070231768 -3.222239146362504855e-07
 107
--1.3428658869178822233 0.97789755368202779806 0.053433985387238869258
--0.0046328365471693128824 -0.01011917844182222935 -2.28774480126113111e-05
+-1.518162912371189055 0.68400065909469953684 0.051574975426518870902
+-0.0021454091976326134308 -0.011564320429164449966 -0.0001141828894241637938
 108
--1.3429299717453899721 0.97783346885452004926 0.053433985387238869258
--0.0107917947436170988396 -0.01011917844182222935 -2.28774480126113111e-05
+-1.5182269971986968038 0.68393657426719178805 0.051574975426518870902
+-0.008304367394080400255 -0.011564320429164449966 -0.0001141828894241637938
 109
-3.9238450924224119731 -3.1677588714750708476 -0.07461478779720470689
-0.041537769840633231855 0.0062326233009544675795 -0.00013004292010574568976
+4.0462039233401467797 -2.9968556860684887333 -0.07806209765007877943
+0.041288813494489193245 0.006425243736188544774 -0.00012527374697365590813
 110
-3.9225232944062176088 -3.169080669491265212 -0.07461478779720470689
--0.03222666456266162771 0.0062326233009544675795 -0.00013004292010574568976
+4.0448821253239524154 -2.9981774840846830976 -0.07806209765007877943
+-0.03247562090880566632 0.006425243736188544774 -0.00012527374697365590813
 111
-6.1863449551503970625 -7.8036243454496956318 -0.11049843292623899582
-0.026180002087518607773 0.0034586370407366112158 -0.00022193109393270140663
+6.295464977775925064 -7.7093855974615586035 -0.11647820911165690516
+0.026129138119387516903 0.0035219001519798780186 -0.00022100689998454499602
 112
-6.1852439704401360743 -7.80472533015995662 -0.11049843292623899582
--0.018046335679548570347 0.0034586370407366112158 -0.00022193109393270140663
+6.2943639930656640757 -7.7104865821718195917 -0.11647820911165690516
+-0.018097199647679661216 0.0035219001519798780186 -0.00022100689998454499602
 113
-14.92933741553430238 12.9297364937164918075 -0.14540991395599550673
-0.010495205578015229936 0.0027958881988585449277 4.3918648577820880246e-05
+14.858937443026691838 13.005046650478242043 -0.14422203290693380584
+0.0104795161786612953114 0.0027826291421300451516 4.4072431030456472162e-05
 114
-14.928857899559897504 12.9292569777420869315 -0.14540991395599550673
--0.015694321499308556966 0.0027958881988585449277 4.3918648577820880246e-05
+14.858457927052286962 13.004567134503837167 -0.14422203290693380584
+-0.01571001089866249159 0.0027826291421300451516 4.4072431030456472162e-05
 115
-29.544394452699595632 -4.7166888412680219034 -0.58380301744279916587
-0.014914899434230450767 0.0031257375729174499863 -7.532640451995010599e-05
+29.557212398525788188 -4.6322721717265853414 -0.5858344271811812831
+0.01490584566440259634 0.0031273464291932001093 -7.514820514613305166e-05
 116
-29.543928927807161955 -4.717154366160452028 -0.58380301744279916587
--0.013956372192325811402 0.0031257375729174499863 -7.532640451995010599e-05
-117
-14.544494696107316045 -31.052223962419365222 -0.88280002656255951443
-0.0037110399802695879026 0.00066259169021535256356 -0.0009142553677224461283
-118
-14.544472229076623293 -31.052246429450054421 -0.88280002656255951443
-0.002135115255113890012 0.00066259169021535256356 -0.0009142553677224461283
+29.556746873633354511 -4.632737696619015466 -0.5858344271811812831
+-0.013965425962153665829 0.0031273464291932001093 -7.514820514613305166e-05

python/swiftest/swiftest/io.py

@@ -659,6 +659,7 @@ def swiftest_xr2infile(ds, param, framenum=-1):
         # Swiftest Central body file
         cbfile = open(param['CB_IN'], 'w')
         print(GMSun, file=cbfile)
+        print(GMSun, file=cbfile)
         print(RSun, file=cbfile)
         print(J2, file=cbfile)
         print(J4, file=cbfile)
@@ -687,6 +688,8 @@ def swiftest_xr2infile(ds, param, framenum=-1):
         # Now make Swiftest files
         cbfile = FortranFile(param['CB_IN'], 'w')
         MSun = np.double(1.0)
+        cbid = 0
+        cbfile.write_record(cbid)
         cbfile.write_record(np.double(GMSun))
         cbfile.write_record(np.double(rmin))
         cbfile.write_record(np.double(J2))

src/io/io.f90

@@ -105,6 +105,9 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
             case ("BIG_DISCARD")
                call util_toupper(param_value)
                if (param_value == "YES" .or. param_value == 'T' ) self%lbig_discard = .true.
+            case ("RHILL_PRESENT")
+               call util_toupper(param_value)
+               if (param_value == "YES" .or. param_value == 'T' ) self%lrhill_present = .true.
             case ("MU2KG")
                read(param_value, *) self%MU2KG
             case ("TU2S")
@@ -213,6 +216,7 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg)
       write(*,*) "ENC_OUT        = ",trim(adjustl(self%encounter_file))
       write(*,*) "EXTRA_FORCE    = ",self%lextra_force
       write(*,*) "BIG_DISCARD    = ",self%lbig_discard
+      write(*,*) "RHILL_PRESENT  = ",self%lrhill_present
       write(*,*) "ENERGY         = ",self%lenergy
       write(*,*) "MU2KG          = ",self%MU2KG       
       write(*,*) "TU2S           = ",self%TU2S        
@@ -315,6 +319,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg)
          write(param_name, Afmt) "MU2KG"; write(param_value, Rfmt) param%MU2KG; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "TU2S"; write(param_value, Rfmt) param%TU2S ; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "DU2M"; write(param_value, Rfmt) param%DU2M; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
+         write(param_name, Afmt) "RHILL_PRESENT"; write(param_value, Lfmt) param%lrhill_present; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "EXTRA_FORCE"; write(param_value, Lfmt) param%lextra_force; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "BIG_DISCARD"; write(param_value, Lfmt) param%lbig_discard; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
          write(param_name, Afmt) "CHK_CLOSE"; write(param_value, Lfmt) param%lclose; write(unit, Afmt) adjustl(param_name), adjustl(param_value)
@@ -592,10 +597,15 @@ module subroutine io_read_body_in(self, param)
                select type(self)
                class is (swiftest_pl)
                   read(iu, *, iostat = ierr) self%id(i), val 
-                  self%mass(i) = real(val / param%GU, kind=DP)
                   self%Gmass(i) = real(val, kind=DP)
+                  self%mass(i) = real(val / param%GU, kind=DP)
+
                   if (param%lclose) then
-                     read(iu, *, iostat = ierr) self%radius(i)
+                     if (param%lrhill_present) then
+                        read(iu, *, iostat = ierr) self%radius(i), self%rhill(i)
+                     else
+                        read(iu, *, iostat = ierr) self%radius(i)
+                     end if
                      if (ierr /= 0 ) exit
                   else
                      self%radius(i) = 0.0_DP
@@ -623,6 +633,7 @@ module subroutine io_read_body_in(self, param)
          ierr = -1
       end select
       close(iu)
+
       if (ierr /= 0 ) then
          write(*,*) 'Error reading in initial conditions from ',trim(adjustl(infile))
          call util_exit(FAILURE)
@@ -932,6 +943,7 @@ module subroutine io_read_initialize_system(self, param)
 
       call self%cb%initialize(param)
       call self%pl%initialize(param)
+      if (.not.param%lrhill_present) call self%pl%set_rhill(self%cb)
       call self%tp%initialize(param)
       call self%set_msys()
       call self%pl%set_mu(self%cb) 

src/modules/swiftest_classes.f90

@@ -61,6 +61,7 @@ module swiftest_classes
       real(DP)             :: inv_c2         = -1.0_DP            !! Inverse speed of light squared in the system units
 
       !Logical flags to turn on or off various features of the code
+      logical :: lrhill_present = .false. !! Hill radii are given as an input rather than calculated by the code (can be used to inflate close encounter regions manually)
       logical :: lextra_force   = .false. !! User defined force function turned on
       logical :: lbig_discard   = .false. !! Save big bodies on every discard
       logical :: lclose         = .false. !! Turn on close encounters

src/rmvs/rmvs_step.f90

@@ -29,7 +29,6 @@ module subroutine rmvs_step_system(self, param, t, dt)
                associate(system => self, ntp => tp%nbody, npl => pl%nbody)
                   allocate(xbeg, source=pl%xh)
                   allocate(vbeg, source=pl%vh)
-                  call pl%set_rhill(cb)
                   call pl%set_beg_end(xbeg = xbeg, vbeg = vbeg)
                   ! ****** Check for close encounters ***** !
                   system%rts = RHSCALE