diff --git a/examples/swifter_conversion/cb.swiftest.in b/examples/swifter_conversion/cb.swiftest.in new file mode 100644 index 000000000..c33aa2b5c --- /dev/null +++ b/examples/swifter_conversion/cb.swiftest.in @@ -0,0 +1,4 @@ +0.00029591220819207774 +0.004650467260962157 +4.7535806948127355e-12 +-2.2473967953572827e-18 diff --git a/examples/swifter_conversion/param.in b/examples/swifter_conversion/param.in new file mode 100644 index 000000000..98cb8a5ce --- /dev/null +++ b/examples/swifter_conversion/param.in @@ -0,0 +1,26 @@ +! Swifter input file generated using init_cond.py +T0 0 +TSTOP 80.0 +DT 1.0 +PL_IN pl.in +TP_IN tp.in +IN_TYPE ASCII +ISTEP_OUT 1 +ISTEP_DUMP 1 +BIN_OUT bin.dat +OUT_TYPE REAL8 +OUT_FORM XV +OUT_STAT NEW +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.dat +EXTRA_FORCE no +BIG_DISCARD no +RHILL_PRESENT yes diff --git a/examples/swifter_conversion/pl.in b/examples/swifter_conversion/pl.in new file mode 100644 index 000000000..b98db89c9 --- /dev/null +++ b/examples/swifter_conversion/pl.in @@ -0,0 +1,40 @@ +10 +1 0.00029591220819207775568 +0.0 0.0 0.0 +0.0 0.0 0.0 +2 4.9125474498983625e-11 0.0014751258227142052 ! mercury +1.6306381826061646e-05 +0.008059842448018334 -0.4616051037329109 -0.03846017738329229 +0.02248719132054853 0.001934639213990692 -0.001904656977422976 +3 7.243452483873647e-10 0.006759134232034942 ! venus +4.0453784346544176e-05 +-0.5115875215389065 0.5030818749037324 0.03642547299277956 +-0.01425515725454357 -0.01452868630179309 0.0006232072038298823 +4 8.997011382166019e-10 0.010044625087011913 ! earthmoon +4.25875607065041e-05 +-0.1090020607540907 -1.009893805009766 4.823302918632528e-05 +0.01682491922568941 -0.001910549762056979 3.992660742687128e-08 +5 9.549535102761465e-11 0.007246789790242477 ! mars +2.2657408050928896e-05 +-1.342897929331636 0.9778655112682739 0.05343398538723887 +-0.007712315645393206 -0.01011917844182223 -2.287744801261131e-05 +6 2.8253459086313547e-07 0.3552720805286442 ! jupiter +0.0004673261703049093 +3.923184193414315 -3.168419770483168 -0.0746147877972047 +0.004655552638985802 0.006232623300954468 -0.0001300429201057457 +7 8.459715183006416e-08 0.4376460836930155 ! saturn +0.00038925687730393614 +6.185794462795267 -7.804174837804826 -0.110498432926239 +0.004066833203985018 0.003458637040736611 -0.0002219310939327014 +8 1.2920249163736674e-08 0.46946272948265794 ! uranus +0.00016953449859497232 +14.9290976575471 12.92949673572929 -0.1454099139559955 +-0.002599557960646664 0.002795888198858545 4.391864857782088e-05 +9 1.5243589003230834e-08 0.7811947848333599 ! neptune +0.00016458790412449367 +29.54416169025338 -4.716921603714237 -0.5838030174427992 +0.0004792636209523189 0.00312573757291745 -7.53264045199501e-05 +10 2.1919422829042796e-12 0.05379680851617536 ! plutocharon +7.943294877391593e-06 +14.54448346259197 -31.05223519593471 -0.8828000265625595 +0.002923077617691739 0.0006625916902153526 -0.0009142553677224461 diff --git a/examples/swifter_conversion/pl.swiftest.in b/examples/swifter_conversion/pl.swiftest.in new file mode 100644 index 000000000..6e41f04b7 --- /dev/null +++ b/examples/swifter_conversion/pl.swiftest.in @@ -0,0 +1,37 @@ +9 +2 4.9125474498983625e-11 +1.6306381826061646e-05 +0.008059842448018334 -0.4616051037329109 -0.03846017738329229 +0.02248719132054853 0.001934639213990692 -0.001904656977422976 +3 7.243452483873647e-10 +4.0453784346544176e-05 +-0.5115875215389065 0.5030818749037324 0.03642547299277956 +-0.01425515725454357 -0.01452868630179309 0.0006232072038298823 +4 8.997011382166019e-10 +4.25875607065041e-05 +-0.1090020607540907 -1.009893805009766 4.823302918632528e-05 +0.01682491922568941 -0.001910549762056979 3.992660742687128e-08 +5 9.549535102761465e-11 +2.2657408050928896e-05 +-1.342897929331636 0.9778655112682739 0.05343398538723887 +-0.007712315645393206 -0.01011917844182223 -2.287744801261131e-05 +6 2.8253459086313547e-07 +0.0004673261703049093 +3.923184193414315 -3.168419770483168 -0.0746147877972047 +0.004655552638985802 0.006232623300954468 -0.0001300429201057457 +7 8.459715183006416e-08 +0.00038925687730393614 +6.185794462795267 -7.804174837804826 -0.110498432926239 +0.004066833203985018 0.003458637040736611 -0.0002219310939327014 +8 1.2920249163736674e-08 +0.00016953449859497232 +14.9290976575471 12.92949673572929 -0.1454099139559955 +-0.002599557960646664 0.002795888198858545 4.391864857782088e-05 +9 1.5243589003230834e-08 +0.00016458790412449367 +29.54416169025338 -4.716921603714237 -0.5838030174427992 +0.0004792636209523189 0.00312573757291745 -7.53264045199501e-05 +10 2.1919422829042796e-12 +7.943294877391593e-06 +14.54448346259197 -31.05223519593471 -0.8828000265625595 +0.002923077617691739 0.0006625916902153526 -0.0009142553677224461 diff --git a/examples/swifter_conversion/swifter2swiftest.py b/examples/swifter_conversion/swifter2swiftest.py new file mode 100755 index 000000000..04c70fed5 --- /dev/null +++ b/examples/swifter_conversion/swifter2swiftest.py @@ -0,0 +1,17 @@ +import sys +import argparse +import swiftestio as swio +""" + Converts initial conditions files from Swifter to Swiftest +""" + +if __name__ == '__main__': + ap = argparse.ArgumentParser() + ap.add_argument("-i", "--input_swifter_param", required=True, help="Input Swifter parameter file to convert") + ap.add_argument("-o", "--output_swiftest_param", required=True, help="Converted Swiftest parameter file") + args = vars(ap.parse_args()) + inparam = args['input_swifter_param'] + outparam = args['output_swiftest_param'] + print(f"Swifter parameter is {inparam}") + print(f"Swiftest parameter file is {outparam}") + swio.swifter2swiftest(inparam,outparam) diff --git a/examples/swifter_conversion/tp.in b/examples/swifter_conversion/tp.in new file mode 100644 index 000000000..9c026369e --- /dev/null +++ b/examples/swifter_conversion/tp.in @@ -0,0 +1,4 @@ +1 +100 +1.01 0.0 0.0 +0.0 6.252003053624663 0.0 diff --git a/examples/swifter_conversion/tp.swiftest.in b/examples/swifter_conversion/tp.swiftest.in new file mode 100644 index 000000000..9c026369e --- /dev/null +++ b/examples/swifter_conversion/tp.swiftest.in @@ -0,0 +1,4 @@ +1 +100 +1.01 0.0 0.0 +0.0 6.252003053624663 0.0 diff --git a/python/swifter2swiftest.py b/python/swifter2swiftest.py deleted file mode 100644 index 8ba09e542..000000000 --- a/python/swifter2swiftest.py +++ /dev/null @@ -1,155 +0,0 @@ -import numpy as np -import sys -import swiftestio as swio - -""" - Converts initial conditions files from Swifter to Swiftest -""" -G = 6.6743E-11 # Universal gravitational constant in SI units - -def read_param(param): - # Read and parse old Swifter parameter file - swifterfile = param['swifterfile'] - swiftestfile = param['swiftestfile'] - - # Read param.in file - print('Reading Swifter file ' + param['swifterfile']) - fold = open(swifterfile, 'r') - swifterlines = fold.readlines() - fold.close() - for line in swifterlines: - fields = line.split() - if len(fields) > 0: - if ('PL_IN' == fields[0].upper()): param['PL_OLD'] = fields[1] - if ('TP_IN' == fields[0].upper()): param['TP_OLD'] = fields[1] - if ('CHK_CLOSE' == fields[0].upper()): param['CHK_CLOSE'] = fields[1].upper() - if ('RHILL_PRESENT'== fields[0].upper()): param['RHILL_PRESENT'] = fields[1].upper() - if ('J2' == fields[0].upper()): param['J2'] = float(fields[1]) - if ('J4' == fields[0].upper()): param['J4'] = float(fields[1]) - - print('Reading Swiftest file ' + param['swiftestfile']) - fnew = open(swiftestfile, 'r') - swiftestlines = fnew.readlines() - for line in swiftestlines: - fields = line.split() - if len(fields) > 0: - if ('MU2KG' == fields[0].upper()): param['MU2KG'] = float(fields[1]) - if ('DU2M' == fields[0].upper()): param['DU2M'] = float(fields[1]) - if ('TU2S' == fields[0].upper()): param['TU2S'] = float(fields[1]) - if ('CB_IN' == fields[0].upper()): param['CB_NEW'] = fields[1] - if ('PL_IN' == fields[0].upper()): param['PL_NEW'] = fields[1] - if ('TP_IN' == fields[0].upper()): param['TP_NEW'] = fields[1] - if ('CHK_CLOSE' == fields[0].upper()): param['CHK_CLOSE_NEW'] = fields[1].upper() - -if __name__ == '__main__': - - param = {'swifterfile' : "param.in", - 'swiftestfile': "param.in"} - - read_param(param) - MU2KG = param['MU2KG'] - DU2M = param['DU2M'] - TU2S = param['TU2S'] - - for key, val in param.items(): - print(key, val) - - GU = G / (DU2M ** 3 / (MU2KG * TU2S ** 2)) - - cbrad = 6.95700e8 / DU2M - - plnew = open(param['PL_NEW'], 'w') - print(f'Writing out new PL file: {param["PL_NEW"]}') - - with open(param['PL_OLD'], 'r') as plold: - line = plold.readline() - line = line.split("!")[0] # Ignore comments - i_list = [i for i in line.split(" ") if i.strip()] - npl = int(i_list[0]) - print(npl) - print(npl - 1, file=plnew) - line = plold.readline() - #line = line[0].split("!")[0] # Ignore comments - i_list = [i for i in line.split(" ") if i.strip()] - GMsun = float(i_list[1]) # Store central body GM for later - line = plold.readline() # Ignore the two zero vector lines - line = plold.readline() - for n in range(1, npl): # Loop over planets - line = plold.readline() - i_list = [i for i in line.split(" ") if i.strip()] - name = int(i_list[0]) - GMpl = float(i_list[1]) - #print(f'{name} {GMpl} -> {GMpl / GU}') - print(name, GMpl,file=plnew) - if param['CHK_CLOSE'] == 'YES': - line = plold.readline() - i_list = [i for i in line.split(" ") if i.strip()] - plrad = float(i_list[0]) - print(plrad) - print(plrad,file=plnew) - line = plold.readline() - i_list = [i for i in line.split(" ") if i.strip()] - xh = float(i_list[0]) - yh = float(i_list[1]) - zh = float(i_list[2]) - print(xh, yh, zh) - print(xh, yh, zh, file=plnew) - line = plold.readline() - i_list = [i for i in line.split(" ") if i.strip()] - vx = float(i_list[0]) - vy = float(i_list[1]) - vz = float(i_list[2]) - print(vx, vy, vz) - print(vx, vy, vz, file=plnew) - - plold.close() - plnew.close() - tpnew = open(param['TP_NEW'], 'w') - - print(f'Writing out new TP file: {param["TP_NEW"]}') - with open(param['TP_OLD'], 'r') as tpold: - line = tpold.readline() - line = line.split("!")[0] # Ignore comments - i_list = [i for i in line.split(" ") if i.strip()] - ntp = int(i_list[0]) - print(ntp) - print(ntp, file=tpnew) - for n in range(0, ntp): # Loop over test particles - line = tpold.readline() - i_list = [i for i in line.split(" ") if i.strip()] - name = int(i_list[0]) - print(name) - print(name, file=tpnew) - line = tpold.readline() - i_list = [i for i in line.split(" ") if i.strip()] - xh = float(i_list[0]) - yh = float(i_list[1]) - zh = float(i_list[2]) - print(xh, yh, zh) - print(xh, yh, zh, file=tpnew) - line = tpold.readline() - i_list = [i for i in line.split(" ") if i.strip()] - vx = float(i_list[0]) - vy = float(i_list[1]) - vz = float(i_list[2]) - print(vx, vy, vz) - print(vx, vy, vz, file=tpnew) - - tpold.close() - tpnew.close() - - print(f'Writing out new CB file: {param["CB_NEW"]}') - # Write out new central body file - cbnew = open(param['CB_NEW'], 'w') - - print(f'GMsun = {GMsun} -> {GMsun / GU}') - print(GMsun, file=cbnew) - print(cbrad) - print(cbrad, file=cbnew) - print(param['J2']) - print(param['J2'], file=cbnew) - print(param['J4']) - print(param['J4'], file=cbnew) - - cbnew.close() - diff --git a/python/swiftestio/swiftestio.py b/python/swiftestio/swiftestio.py index 513d63665..71a945a1d 100644 --- a/python/swiftestio/swiftestio.py +++ b/python/swiftestio/swiftestio.py @@ -1,5 +1,4 @@ import numpy as np -import pandas as pd from scipy.io import FortranFile import xarray as xr from astroquery.jplhorizons import Horizons @@ -7,12 +6,13 @@ import datetime # Constants in SI units +GC = np.longdouble(const.G.value) AU2M = np.longdouble(const.au.value) GMSunSI = np.longdouble(const.GM_sun.value) +MSun = np.longdouble(const.M_sun.value) RSun = np.longdouble(const.R_sun.value) -GC = np.longdouble(const.G.value) JD2S = 86400 -year = np.longdouble(365.25 * JD2S) +YR2S = 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) @@ -156,7 +156,6 @@ def read_swiftest_param(param_file_name): 'TU2S' : -1.0, 'DU2M' : -1.0, 'GU' : -1.0, - 'INV_C2' : -1.0, 'EXTRA_FORCE' : 'NO', 'BIG_DISCARD' : 'NO', 'CHK_CLOSE' : 'NO', @@ -204,7 +203,6 @@ def read_swiftest_param(param_file_name): param['DU2M'] = float(param['DU2M']) param['MU2KG'] = float(param['MU2KG']) param['TU2S'] = float(param['TU2S']) - param['INV_C2'] = float(param['INV_C2']) param['EXTRA_FORCE'] = param['EXTRA_FORCE'].upper() param['BIG_DISCARD'] = param['BIG_DISCARD'].upper() param['CHK_CLOSE'] = param['CHK_CLOSE'].upper() @@ -558,7 +556,7 @@ def solar_system_pl(param, ephemerides_start_date): 'plutocharon': '9' } - # Planet Msun/M ratio + # Planet MSun/M ratio MSun_over_Mpl = { 'mercury': np.longdouble(6023600.0), 'venus': np.longdouble(408523.71), @@ -727,7 +725,7 @@ def swiftest_xr2_infile(ds, param, framenum=-1): elif param['IN_TYPE'] == 'REAL8': # Now make Swiftest files cbfile = FortranFile(swiftest_cb, 'w') - Msun = np.double(1.0) + MSun = np.double(1.0) cbfile.write_record(np.double(GMSun)) cbfile.write_record(np.double(rmin)) cbfile.write_record(np.double(J2)) @@ -760,6 +758,158 @@ def swiftest_xr2_infile(ds, param, framenum=-1): else: print(f"{param['IN_TYPE']} is an unknown file type") +def swifter2swiftest(inparam,outparam): + print(f"Swifter parameter is {inparam}") + print(f"Swiftest parameter file is {outparam}") + swifter_param = read_swifter_param(inparam) + swiftest_param = swifter_param.copy() + print("Select the unit system to use:") + print("1) MSun-AU-year") + print("2) MSun-AU-day") + print("3) SI: kg-m-s") + print("4) CGS: g-cm-s") + print("5) Set units manually") + inval = input("> ") + try: + unit_type = int(inval) + except ValueError: + goodval = False + else: + goodval = (unit_type > 0 and unit_type < 6) + if not goodval: + print(f"{inval} is not a valid menu option") + sys.exit(-1) + if unit_type == 1: + print("Unit system is MSun-AU-year") + swiftest_param['MU2KG'] = MSun + swiftest_param['DU2M'] = AU2M + swiftest_param['TU2S'] = YR2S + elif unit_type == 2: # MSun-AU-day + print("Unit system is MSun-AU-day") + swiftest_param['MU2KG'] = MSun + swiftest_param['DU2M'] = AU2M + swiftest_param['TU2S'] = JD2S + elif unit_type == 3: # SI: kg-m-s + print("Unit system is SI: kg-m-s") + swiftest_param['MU2KG'] = 1.0 + swiftest_param['DU2M'] = 1.0 + swiftest_param['TU2S'] = 1.0 + elif unit_type == 4: # CGS: g-cm-s + print("Unit system is CGS: g-cm-s") + swiftest_param['MU2KG'] = 1e-3 + swiftest_param['DU2M'] = 1.0e-2 + swiftest_param['TU2S'] = 1.0 + elif unit_type == 5: + print("User-defined units.") + print("Define each unit (mass, distance, and time) by its corresponding SI value.") + swiftest_param['MU2KG'] = input("Mass value in kilograms: ") + swiftest_param['DU2M'] = input("Distance value in meters: ") + swiftest_param['TU2S'] = input("Time unit in seconds: ") + + print("Set central body radius:") + print(f"1) Use Swifter parameter value of CHK_RMIN = {swifter_param['CHK_RMIN']}") + print(f"2) Set value manually") + inval = input("> ") + try: + cbrad_type = int(inval) + except ValueError: + goodval = False + else: + goodval = (cbrad_type > 0 and cbrad_type < 3) + if not goodval: + print(f"{inval} is not a valid menu option") + sys.exit(-1) + if cbrad_type == 1: + cbrad = swifter_param['CHK_RMIN'] + elif cbrad_type == 2: + cbrad = input("Enter radius of central body in simulation Distance Units: ") + + swiftest_param['PL_IN'] = 'pl.swiftest.in' + swiftest_param['TP_IN'] = 'tp.swiftest.in' + swiftest_param['CB_IN'] = 'cb.swiftest.in' + + plnew = open(swiftest_param['PL_IN'], 'w') + + print(f'Writing out new PL file: {swiftest_param["PL_IN"]}') + with open(swifter_param['PL_IN'], 'r') as plold: + line = plold.readline() + line = line.split("!")[0] # Ignore comments + i_list = [i for i in line.split(" ") if i.strip()] + npl = int(i_list[0]) + print(npl - 1, file=plnew) + line = plold.readline() + i_list = [i for i in line.split(" ") if i.strip()] + GMcb = float(i_list[1]) # Store central body GM for later + line = plold.readline() # Ignore the two zero vector lines + line = plold.readline() + for n in range(1, npl): # Loop over planets + line = plold.readline() + i_list = [i for i in line.split(" ") if i.strip()] + name = int(i_list[0]) + GMpl = float(i_list[1]) + print(name, GMpl,file=plnew) + if swifter_param['CHK_CLOSE'] == 'YES': + line = plold.readline() + i_list = [i for i in line.split(" ") if i.strip()] + plrad = float(i_list[0]) + print(plrad,file=plnew) + line = plold.readline() + i_list = [i for i in line.split(" ") if i.strip()] + xh = float(i_list[0]) + yh = float(i_list[1]) + zh = float(i_list[2]) + print(xh, yh, zh, file=plnew) + line = plold.readline() + i_list = [i for i in line.split(" ") if i.strip()] + vx = float(i_list[0]) + vy = float(i_list[1]) + vz = float(i_list[2]) + print(vx, vy, vz, file=plnew) + + plold.close() + plnew.close() + tpnew = open(swiftest_param['TP_IN'], 'w') + + print(f'Writing out new TP file: {swiftest_param["TP_IN"]}') + with open(swifter_param['TP_IN'], 'r') as tpold: + line = tpold.readline() + line = line.split("!")[0] # Ignore comments + i_list = [i for i in line.split(" ") if i.strip()] + ntp = int(i_list[0]) + print(ntp, file=tpnew) + for n in range(0, ntp): # Loop over test particles + line = tpold.readline() + i_list = [i for i in line.split(" ") if i.strip()] + name = int(i_list[0]) + print(name, file=tpnew) + line = tpold.readline() + i_list = [i for i in line.split(" ") if i.strip()] + xh = float(i_list[0]) + yh = float(i_list[1]) + zh = float(i_list[2]) + print(xh, yh, zh, file=tpnew) + line = tpold.readline() + i_list = [i for i in line.split(" ") if i.strip()] + vx = float(i_list[0]) + vy = float(i_list[1]) + vz = float(i_list[2]) + print(vx, vy, vz, file=tpnew) + + tpold.close() + tpnew.close() + + print(f'Writing out new CB file: {swiftest_param["CB_IN"]}') + # Write out new central body file + cbnew = open(swiftest_param['CB_IN'], 'w') + + print(GMcb, file=cbnew) + print(cbrad, file=cbnew) + print(swifter_param['J2'], file=cbnew) + print(swifter_param['J4'], file=cbnew) + + cbnew.close() + return + if __name__ == '__main__': workingdir = '/Users/daminton/git/swiftest/examples/rmvs_swifter_comparison/9pl_18tp_encounters/'