diff --git a/examples/helio_swifter_comparison/.idea/.gitignore b/examples/helio_swifter_comparison/.idea/.gitignore new file mode 100644 index 000000000..26d33521a --- /dev/null +++ b/examples/helio_swifter_comparison/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/examples/helio_swifter_comparison/cb.swiftest.in b/examples/helio_swifter_comparison/cb.swiftest.in index 7e13534c0..058975b81 100644 Binary files a/examples/helio_swifter_comparison/cb.swiftest.in and b/examples/helio_swifter_comparison/cb.swiftest.in differ diff --git a/examples/helio_swifter_comparison/init_cond.py b/examples/helio_swifter_comparison/init_cond.py index 80b36c1c4..ed46b12cf 100644 --- a/examples/helio_swifter_comparison/init_cond.py +++ b/examples/helio_swifter_comparison/init_cond.py @@ -1,323 +1,54 @@ +import swiftest import numpy as np import sys from astroquery.jplhorizons import Horizons import astropy.constants as const from scipy.io import FortranFile -#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) -MSun_over_Mpl = np.array([6023600.0, - 408523.71, - 328900.56, - 3098708., - 1047.3486, - 3497.898, - 22902.98, - 19412.24, - 1.35e8], dtype=np.longdouble) - -MU2KG = np.longdouble(GMSunSI / GC) #Conversion from mass unit to kg -DU2M = np.longdouble(AU2M) #Conversion from radius unit to centimeters -TU2S = np.longdouble(year) #Conversion from time unit to seconds -GU = np.longdouble(GC / (DU2M**3 / (MU2KG * TU2S**2))) - -GMSun = np.longdouble(GMSunSI / (DU2M**3 / TU2S**2)) - -# Simulation start, stop, and output cadence times -t_0 = 0 # simulation start time -deltaT = 0.25 * JD / TU2S # simulation step size -end_sim = 1 * year / TU2S # simulation end time -t_print = deltaT #year / TU2S #output interval to print results - - -# Solar oblatenes values: From Mecheri et al. (2004), using Corbard (b) 2002 values (Table II) -J2 = 0.0 #np.longdouble(2.198e-7) * (Rsun / DU2M)**2 -J4 = 0.0 #np.longdouble(-4.805e-9) * (Rsun / DU2M)**4 - -tstart = '2021-01-28' -tend = '2021-01-29' -tstep = '1d' -planetid = { - 'mercury' : '1', - 'venus' : '2', - 'earthmoon' : '3', - 'mars' : '4', - 'jupiter' : '5', - 'saturn' : '6', - 'uranus' : '7', - 'neptune' : '8', - 'plutocharon' : '9' -} -npl = 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) -} - -pdata = {} -plvec = {} -Rhill = {} -THIRDLONG = np.longdouble(1.0) / np.longdouble(3.0) - -for key,val in planetid.items(): - pdata[key] = Horizons(id=val, id_type='majorbody',location='@sun', - epochs={'start': tstart, 'stop': tend, - 'step': tstep}) - plvec[key] = np.array([pdata[key].vectors()['x'][0], - pdata[key].vectors()['y'][0], - pdata[key].vectors()['z'][0], - pdata[key].vectors()['vx'][0], - pdata[key].vectors()['vy'][0], - pdata[key].vectors()['vz'][0] - ]) - - Rhill[key] = np.longdouble(pdata[key].elements()['a'][0]) * (3 * MSun_over_Mpl[key])**(-THIRDLONG) - -asteroidid = { - '100001' : 'Ceres', - '100002' : 'Pallas', - '100003' : 'Juno', - '100004' : 'Vesta' +sim = swiftest.Simulation() + +sim.param['MU2KG'] = swiftest.MSun +sim.param['TU2S'] = swiftest.YR2S +sim.param['DU2M'] = swiftest.AU2M +sim.param['T0'] = 0.0 +sim.param['TSTOP'] = 1.0 +sim.param['DT'] = 0.25 * swiftest.JD2S / swiftest.YR2S +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['ISTEP_OUT'] = 1 +sim.param['ISTEP_DUMP'] = 1 + +sim.param['GR'] = 'NO' + +bodyid = { + "Sun": 0, + "Mercury": 1, + "Venus": 2, + "Earth": 3, + "Mars": 4, + "Jupiter": 5, + "Saturn": 6, + "Uranus": 7, + "Neptune": 8, + "Ceres": 101, + "Pallas": 102, + "Juno": 103, + "Vesta": 104 } -ntp = 4 -tdata = {} -tpvec = {} -for key,val in asteroidid.items(): - tdata[key] = Horizons(id=val, id_type='smallbody', location='@sun', - epochs={'start': tstart, 'stop': tend, - 'step': tstep}) - tpvec[key] = np.array([tdata[key].vectors()['x'][0], - tdata[key].vectors()['y'][0], - tdata[key].vectors()['z'][0], - tdata[key].vectors()['vx'][0], - tdata[key].vectors()['vy'][0], - tdata[key].vectors()['vz'][0] - ]) - - -if __name__ == '__main__': - # Convert from AU-day to AU-year just because I find it easier to keep track of the sim progress - for plid in plvec: - plvec[plid][3:] *= year / JD - - for tpid in tpvec: - tpvec[tpid][3:] *= year / JD - - # Names of all output files - swifter_input = "param.swifter.in" - swifter_pl = "pl.swifter.in" - swifter_tp = "tp.swifter.in" - swifter_bin = "bin.swifter.dat" - swifter_enc = "enc.swifter.dat" - - swiftest_input = "param.swiftest.in" - swiftest_pl = "pl.swiftest.in" - swiftest_tp = "tp.swiftest.in" - swiftest_cb = "cb.swiftest.in" - swiftest_bin = "bin.swiftest.dat" - swiftest_enc = "enc.swiftest.dat" - - iout = int(np.ceil(t_print / deltaT)) - rmin = Rsun / DU2M - rmax = np.longdouble(1000.0) - #Make Swifter files +for name, id in bodyid.items(): + sim.add(name, idval=id) - plfile = open(swifter_pl, 'w') - print(npl+1, f'! Planet input file generated using init_cond.py using JPL Horizons data for the major planets (and Pluto) for epoch {tstart}' ,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 i, plid in enumerate(plvec): - print(i + 2,"{:.23g}".format(GMSun * MSun_over_Mpl[plid]**-1),Rhill[plid], file=plfile) - print(Rpl[plid] / DU2M, file=plfile) - print(plvec[plid][0],plvec[plid][1],plvec[plid][2], file=plfile) - print(plvec[plid][3],plvec[plid][4],plvec[plid][5], file=plfile) - plfile.close() - - tpfile = open(swifter_tp, 'w') - print(ntp,file=tpfile) - for tpid, tp in tpvec.items(): - print(tpid, file=tpfile) - print(tp[0],tp[1],tp[2], file=tpfile) - print(tp[3],tp[4],tp[5], file=tpfile) - tpfile.close() - - sys.stdout = open(swifter_input, "w") - print('! Swifter input file generated using init_cond.py') - print('T0 ',t_0) - print('TSTOP ',end_sim) - print('DT ',deltaT) - print('PL_IN ',swifter_pl) - print('TP_IN ',swifter_tp) - print('IN_TYPE ASCII') - print('ISTEP_OUT ',iout) - print('ISTEP_DUMP ',iout) - print('BIN_OUT ',swifter_bin) - print('OUT_TYPE REAL8') - print('OUT_FORM XV') - print('OUT_STAT NEW') - print('J2 ',J2) - print('J4 ',J4) - print('CHK_CLOSE yes') - print('CHK_RMIN ',rmin) - print('CHK_RMAX ',rmax) - print('CHK_EJECT ',rmax) - print('CHK_QMIN ',rmin) - print('CHK_QMIN_COORD HELIO') - print('CHK_QMIN_RANGE ',rmin,rmax) - print('ENC_OUT ',swifter_enc) - print('EXTRA_FORCE no') - print('BIG_DISCARD no') - print('RHILL_PRESENT yes') - - sys.stdout = sys.__stdout__ - #Now make Swiftest files - cbfile = open(swiftest_cb, 'w') - #cbfile = FortranFile(swiftest_cb, 'w') - print(1.0,file=cbfile) - print(rmin,file=cbfile) - print(J2,file=cbfile) - print(J4,file=cbfile) - 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 = open(swiftest_pl, 'w') - plfile = FortranFile(swiftest_pl, 'w') - #print(npl,file=plfile) - plfile.write_record(npl) - - name = np.empty(npl, dtype=np.int32) - px = np.empty(npl, dtype=np.double) - py = np.empty(npl, dtype=np.double) - pz = np.empty(npl, dtype=np.double) - vx = np.empty(npl, dtype=np.double) - vy = np.empty(npl, dtype=np.double) - vz = np.empty(npl, dtype=np.double) - mass = np.empty(npl, dtype=np.double) - Gmass = np.empty(npl, dtype=np.double) - radius = np.empty(npl, dtype=np.double) - for i, plid in enumerate(plvec): - name[i] = i + 2 - px[i] = plvec[plid][0] - py[i] = plvec[plid][1] - pz[i] = plvec[plid][2] - vx[i] = plvec[plid][3] - vy[i] = plvec[plid][4] - vz[i] = plvec[plid][5] - Gmass[i] = GMSun * MSun_over_Mpl[plid]**-1 - radius[i] = Rpl[plid] / DU2M - plfile.write_record(name.T) - plfile.write_record(px.T) - plfile.write_record(py.T) - plfile.write_record(pz.T) - plfile.write_record(vx.T) - plfile.write_record(vy.T) - plfile.write_record(vz.T) - plfile.write_record(Gmass.T) - plfile.write_record(radius.T) - #for i, plid in enumerate(plvec): - # print(i + 2,"{:.23g}".format(np.longdouble(MSun_over_Mpl[plid]**-1)), file=plfile) - # print(Rpl[plid] / DU2M, file=plfile) - # print(plvec[plid][0], plvec[plid][1], plvec[plid][2], file=plfile) - # print(plvec[plid][3], plvec[plid][4], plvec[plid][5], file=plfile) - plfile.close() - #tpfile = open(swiftest_tp, 'w') - tpfile = FortranFile(swiftest_tp, 'w') - #print(ntp,file=tpfile) - tpfile.write_record(ntp) - #for tpid, tp in tpvec.items(): - # print(tpid, file=tpfile) - # print(tp[0],tp[1],tp[2], file=tpfile) - # print(tp[3],tp[4],tp[5], file=tpfile) - - name = np.empty(ntp, dtype=np.int32) - px = np.empty(ntp, dtype=np.double) - py = np.empty(ntp, dtype=np.double) - pz = np.empty(ntp, dtype=np.double) - vx = np.empty(ntp, dtype=np.double) - vy = np.empty(ntp, dtype=np.double) - vz = np.empty(ntp, dtype=np.double) - for i, tpid in enumerate(tpvec): - name[i] = int(tpid) - px[i] = tpvec[tpid][0] - py[i] = tpvec[tpid][1] - pz[i] = tpvec[tpid][2] - vx[i] = tpvec[tpid][3] - vy[i] = tpvec[tpid][4] - vz[i] = tpvec[tpid][5] - tpfile.write_record(name.T) - tpfile.write_record(px.T) - tpfile.write_record(py.T) - tpfile.write_record(pz.T) - tpfile.write_record(vx.T) - tpfile.write_record(vy.T) - tpfile.write_record(vz.T) - - tpfile.close() - - sys.stdout = open(swiftest_input, "w") - print('! Swiftest input file generated using init_cond.py') - print('T0 ',t_0) - print('TSTOP ',end_sim) - print('DT ',deltaT) - print('CB_IN ',swiftest_cb) - print('PL_IN ',swiftest_pl) - print('TP_IN ',swiftest_tp) - print('IN_TYPE REAL8') - print('ISTEP_OUT ',iout) - print('ISTEP_DUMP ',iout) - print('BIN_OUT ',swiftest_bin) - print('OUT_TYPE REAL8') - print('OUT_FORM XV') - print('OUT_STAT REPLACE') - print('CHK_CLOSE yes') - print('CHK_RMIN ',rmin) - print('CHK_RMAX ',rmax) - print('CHK_EJECT ',rmax) - print('CHK_QMIN ',rmin) - print('CHK_QMIN_COORD HELIO') - print('CHK_QMIN_RANGE ',rmin,rmax) - print('ENC_OUT ',swiftest_enc) - print('EXTRA_FORCE no') - print('BIG_DISCARD no') - print('ROTATION no') - print('GR no') - print('MU2KG ',MU2KG) - print('DU2M ',DU2M) - print('TU2S ',TU2S) +sim.param['PL_IN'] = "pl.swiftest.in" +sim.param['TP_IN'] = "tp.swiftest.in" +sim.param['CB_IN'] = "cb.swiftest.in" +sim.save("param.swiftest.in") +sim.param['PL_IN'] = "pl.swifter.in" +sim.param['TP_IN'] = "tp.swifter.in" +sim.save("param.swifter.in", codename="Swifter") - sys.stdout = sys.__stdout__ diff --git a/examples/helio_swifter_comparison/param.swifter.in b/examples/helio_swifter_comparison/param.swifter.in index 5156bfbc0..b51b04137 100644 --- a/examples/helio_swifter_comparison/param.swifter.in +++ b/examples/helio_swifter_comparison/param.swifter.in @@ -1,26 +1,36 @@ -! Swifter input file generated using init_cond.py -T0 0 -TSTOP 1.0 -DT 0.0006844626967830253251 -PL_IN pl.swifter.in -TP_IN tp.swifter.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 0.0 -J4 0.0 -CHK_CLOSE yes -CHK_RMIN 0.0046504672609621575315 -CHK_RMAX 1000.0 -CHK_EJECT 1000.0 -CHK_QMIN 0.0046504672609621575315 -CHK_QMIN_COORD HELIO -CHK_QMIN_RANGE 0.0046504672609621575315 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 1.0 +DT 0.0006844626967830253 +ISTEP_OUT 1 +ISTEP_DUMP 1 +OUT_FORM EL +OUT_TYPE REAL8 +OUT_STAT REPLACE +IN_TYPE ASCII +PL_IN pl.swifter.in +TP_IN tp.swifter.in +CB_IN None +BIN_OUT bin.dat +ENC_OUT enc.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 None +TU2S None +DU2M None +EXTRA_FORCE NO +BIG_DISCARD NO +CHK_CLOSE YES +FRAGMENTATION NO +ROTATION NO +TIDES NO +ENERGY NO +YARKOVSKY NO +YORP NO +J2 4.7535806948127355e-12 +J4 -2.2473967953572827e-18 +RHILL_PRESENT YES diff --git a/examples/helio_swifter_comparison/param.swiftest.in b/examples/helio_swifter_comparison/param.swiftest.in index c32a270f5..43e19d180 100644 --- a/examples/helio_swifter_comparison/param.swiftest.in +++ b/examples/helio_swifter_comparison/param.swiftest.in @@ -1,29 +1,35 @@ -! Swiftest input file generated using init_cond.py -T0 0 -TSTOP 1.0 -DT 0.0006844626967830253251 -CB_IN cb.swiftest.in -PL_IN pl.swiftest.in -TP_IN tp.swiftest.in -IN_TYPE REAL8 -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.0046504672609621575315 -CHK_RMAX 1000.0 -CHK_EJECT 1000.0 -CHK_QMIN 0.0046504672609621575315 -CHK_QMIN_COORD HELIO -CHK_QMIN_RANGE 0.0046504672609621575315 1000.0 -ENC_OUT enc.swiftest.dat -EXTRA_FORCE no -BIG_DISCARD no -ROTATION no -GR no -MU2KG 1.988409870698050917e+30 -DU2M 149597870700.0 -TU2S 31557600.0 +! VERSION Swiftest parameter input +T0 0.0 +TSTOP 1.0 +DT 0.0006844626967830253 +ISTEP_OUT 1 +ISTEP_DUMP 1 +OUT_FORM EL +OUT_TYPE REAL8 +OUT_STAT REPLACE +IN_TYPE ASCII +PL_IN pl.swiftest.in +TP_IN tp.swiftest.in +CB_IN cb.swiftest.in +BIN_OUT bin.dat +ENC_OUT enc.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 31557600.0 +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 diff --git a/examples/helio_swifter_comparison/pl.swifter.in b/examples/helio_swifter_comparison/pl.swifter.in index 3f27a30f5..86255d0e2 100644 --- a/examples/helio_swifter_comparison/pl.swifter.in +++ b/examples/helio_swifter_comparison/pl.swifter.in @@ -1,40 +1,36 @@ -10 ! Planet input file generated using init_cond.py using JPL Horizons data for the major planets (and Pluto) for epoch 2021-01-28 -1 39.476926408897625193 +9 +0 39.47692640889762629 0.0 0.0 0.0 0.0 0.0 0.0 -2 6.553709809565313959502e-06 0.0014751229680878679603 +1 6.553709809565314146e-06 0.0014751254963649629896 1.6306381826061645943e-05 -0.1030256871817049 0.2897796042186122 0.0142290454578165 --11.740042086233977 3.834312447819755 1.3902497138650993 -3 9.6633133995815373809037e-05 0.0067591276497677570562 +0.359124056979876094 -0.1001978128323056938 -0.041130148620746292965 +0.7664364270424182397 10.3592906410849091145 0.7762248217818495593 +2 9.6633133995815381836e-05 0.00675911390647655787 4.0453784346544178454e-05 -0.06110217931597178 -0.7245466902000399 -0.01346904334134495 -7.311995450520371 0.5941125626853079 -0.4137913834645579 -4 0.00012002693582795244940133 0.010044756567554867971 +-0.709853246614207567 0.109615461427968005625 0.042466530791895232277 +-1.166834223638398553 -7.334297883841826485 -0.033323414543104576783 +3 0.000120026935827952456416 0.010044751446422201634 4.25875607065040958e-05 --0.6061796332978544 0.7761214562175627 -3.47496925510033e-05 --5.054824318940399 -3.891667462296022 0.0001972008358555219 -5 1.2739802010675940622175e-05 0.0072464490746761575477 +0.26014404284638581455 -0.9828537227999029069 4.5807148740206238052e-05 +5.9724418390973225248 1.5843954077771575533 -9.4205748659356694786e-05 +4 1.2739802010675941808e-05 0.0072467561525263854473 2.265740805092889601e-05 -0.2751944182282239 1.519376889861282 0.02508924669876031 --4.835983592340031 1.3448550964716055 0.14681413009269897 -6 0.037692251088985675999687 0.35528523591470134027 +-1.4908630412685239808 0.7412277078494349247 0.052104480532706012874 +-2.084278892390818102 -4.1405652065758745757 -0.035644761583621103612 +5 0.03769225108898567778 0.35527141892920680066 0.00046732617030490929307 -3.20013535487531 -3.953498299871184 -0.05517746960976119 -2.111393786386104 1.8660266390974023 -0.05498927104920874 -7 0.011285899820091272946487 0.43763064535027282272 +4.0233930071159198505 -3.029555621945668964 -0.077433472926114965684 +1.626590141045528945 2.3340622087669935288 -0.046085347207395002237 +6 0.01128589982009127331 0.43765136932522133578 0.00038925687730393611812 -5.607382154387239 -8.258649115739544 -0.07958445350697703 -1.5748468612768933 1.1414574628814913 -0.0825033126704916 -8 0.0017236589478267731051497 0.46909694410284686372 +6.274810893232299236 -7.7275164380757708216 -0.115372736553069593635 +1.4703000143673246375 1.2821134193800077011 -0.08078666716402813097 +7 0.001723658947826773068 0.46952007057521305831 0.00016953449859497231466 -15.2822545079529 12.53905417292343 -0.1514152417388489 --0.9198471609855067 1.0454391001033898 0.015745441522152032 -9 0.0020336100526728304385693 0.78071933134573728207 +14.871766666738729157 12.9908875920566391216 -0.14444232402201501175 +-0.9541590491729433116 1.0172543087941671172 0.016087073469786578863 +8 0.0020336100526728302882 0.78126715446178642547 0.000164587904124493665 -29.47483430798945 -5.147687601125057 -0.5733452880201296 -0.1919169794537526 1.1385109490430123 -0.027844342997168703 -10 2.9242167710294538257026e-07 0.053834666988627116535 -7.943294877391593783e-06 -14.14001068477669 -31.141412045097 -0.7565711098949569 -1.0733960936923237 0.23003153432240253 -0.3342452292632091 +29.554624389819270647 -4.648140925388063671 -0.5854586034520335991 +0.1723572655485145611 1.1421549698170996955 -0.027459964210413734165 diff --git a/examples/helio_swifter_comparison/pl.swiftest.in b/examples/helio_swifter_comparison/pl.swiftest.in index 6ee1ca383..c1f0be29e 100644 Binary files a/examples/helio_swifter_comparison/pl.swiftest.in and b/examples/helio_swifter_comparison/pl.swiftest.in differ diff --git a/examples/helio_swifter_comparison/tp.swifter.in b/examples/helio_swifter_comparison/tp.swifter.in index 74f1cc8cb..62acd79fc 100644 --- a/examples/helio_swifter_comparison/tp.swifter.in +++ b/examples/helio_swifter_comparison/tp.swifter.in @@ -1,13 +1,13 @@ 4 -100001 -2.89438049451288 0.2060633440061593 -0.5267473063848516 --0.3678177533773833 3.511709555047983 0.17870206084842993 -100002 -2.402157085583321 -2.063650726188297 1.221462218318817 -1.9929454704322658 1.7670859767767524 -1.3920312816007854 -100003 --1.762420789095617 -2.766072787712991 0.6998178860363623 -2.5146401571141386 -1.536845348435182 0.24681698261461876 -100004 --2.136886752000925 1.023684503224212 0.2293900119298221 --1.391819314123066 -3.7994616469964826 0.2829665835356411 +101 +2.3133253483335658451 1.6360857008750779862 -0.37450983998533471375 +-2.2458876465769251093 2.8378699270656317882 0.50346273267874514076 +102 +3.009555158239280992 -1.1130165423439479788 0.51172110509120705135 +0.70453633545041942506 2.5148434686651768256 -1.80152331908826862 +103 +-0.5218824163555056961 -3.1396467647675119217 0.7342355813480357929 +3.0582031698647593751 -0.12050283730110719834 -0.096945705299882042706 +104 +-2.075368356279947868 -0.76569201199778380573 0.27541025252901979448 +1.7615515330387480359 -3.9484151677488075983 -0.096278788580453326945 diff --git a/examples/helio_swifter_comparison/tp.swiftest.in b/examples/helio_swifter_comparison/tp.swiftest.in index 5fe4b5618..62acd79fc 100644 Binary files a/examples/helio_swifter_comparison/tp.swiftest.in and b/examples/helio_swifter_comparison/tp.swiftest.in differ diff --git a/python/swiftest/swiftest/init_cond.py b/python/swiftest/swiftest/init_cond.py index e9b32bda0..6e5048c0f 100644 --- a/python/swiftest/swiftest/init_cond.py +++ b/python/swiftest/swiftest/init_cond.py @@ -5,7 +5,7 @@ from datetime import date import xarray as xr -def solar_system_horizons(plname, param, ephemerides_start_date, ds): +def solar_system_horizons(plname, idval, param, ephemerides_start_date, ds): """ Initializes a Swiftest dataset containing the major planets of the Solar System at a particular data from JPL/Horizons @@ -36,10 +36,15 @@ def solar_system_horizons(plname, param, ephemerides_start_date, ds): 'Pluto': '9' } - if plname not in planetid: - print(f"{plname} not found or not yet supported") - print("Valid inputs are: ".join(f"{key}" for key, value in planetid.items())) - return + if plname in planetid: + ispl = True + else: + ispl = False + print(f"\nMassive body {plname} not found or not yet supported") + print("This will be created as a massless test particle") + if idval is None: + print("ID value required for this input type") + return # Planet MSun/M ratio MSun_over_Mpl = { @@ -97,6 +102,12 @@ def solar_system_horizons(plname, param, ephemerides_start_date, ds): plab.append('capom') plab.append('omega') plab.append('capm') + tlab.append('a') + tlab.append('e') + tlab.append('inc') + tlab.append('capom') + tlab.append('omega') + tlab.append('capm') elif param['OUT_FORM'] == 'EL': plab.append('px') plab.append('py') @@ -104,20 +115,31 @@ def solar_system_horizons(plname, param, ephemerides_start_date, ds): plab.append('vx') plab.append('vy') plab.append('vz') + tlab.append('px') + tlab.append('py') + tlab.append('pz') + tlab.append('vx') + tlab.append('vy') + tlab.append('vz') plab.append('Rhill') dims = ['time', 'id', 'vec'] t = np.array([0.0]) key = plname - val = planetid[key] + if ispl: + val = planetid[key] + else: + val = -1 if key == "Sun" : # Create central body + print("Creating the Sun as a central body") cb = [] cbframe = np.expand_dims(cvec.T, axis=0) cbxr = xr.DataArray(cbframe, dims=dims, coords={'time': t, 'id': cbid, 'vec': clab}) cbds = cbxr.to_dataset(dim='vec') ds = xr.combine_by_coords([ds, cbds]) else: # Fetch solar system ephemerides from Horizons + print(f"Fetching ephemerides data for {key} from JPL/Horizons") pl = [] p1 = [] p2 = [] @@ -136,9 +158,14 @@ def solar_system_horizons(plname, param, ephemerides_start_date, ds): GMpl = [] pldata = {} - pldata[key] = Horizons(id=val, id_type='majorbody', location='@sun', - epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, - 'step': ephemerides_step}) + if ispl: + pldata[key] = Horizons(id=val, id_type='majorbody', location='@sun', + epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, + 'step': ephemerides_step}) + else: + pldata[key] = Horizons(id=key, id_type='smallbody', location='@sun', + epochs={'start': ephemerides_start_date, 'stop': ephemerides_end_date, + 'step': ephemerides_step}) if param['OUT_FORM'] == 'XV': p1.append(pldata[key].vectors()['x'][0] * DCONV) p2.append(pldata[key].vectors()['y'][0] * DCONV) @@ -165,12 +192,20 @@ def solar_system_horizons(plname, param, ephemerides_start_date, ds): p10.append(pldata[key].vectors()['vx'][0] * VCONV) p11.append(pldata[key].vectors()['vy'][0] * VCONV) p12.append(pldata[key].vectors()['vz'][0] * VCONV) - 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.array([planetid[key]], dtype=int) - pvec = np.vstack([p1, p2, p3, p4, p5, p6, GMpl, Rpl, p7, p8, p9, p10, p11, p12, Rhill]) + if ispl: + 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 + pvec = np.vstack([p1, p2, p3, p4, p5, p6, GMpl, Rpl, p7, p8, p9, p10, p11, p12, Rhill]) + else: + pvec = np.vstack([p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12]) + plab = tlab.copy() + + if idval is None: + plid = np.array([planetid[key]], dtype=int) + else: + plid = np.array([idval], dtype=int) # Prepare frames by adding an extra axis for the time coordinate plframe = np.expand_dims(pvec.T, axis=0) diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 03c6506c4..a9200c996 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -288,7 +288,7 @@ def write_labeled_param(param, param_file_name): ptmp = param.copy() # Print the list of key/value pairs in the preferred order for key in keylist: - val = ptmp.pop(key) + val = ptmp.pop(key, None) print(f"{key:<16} {val}", file=outfile) # Print the remaining key/value pairs in whatever order for key, val in ptmp.items(): @@ -636,7 +636,7 @@ def swiftest_xr2infile(ds, param, framenum=-1): 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. param : dict - Swiftest paramuration parameters. This method uses the names of the cb, pl, and tp files from the paramuration + Swiftest input parameters. This method uses the names of the cb, pl, and tp files from the input Returns ------- @@ -717,6 +717,69 @@ def swiftest_xr2infile(ds, param, framenum=-1): else: print(f"{param['IN_TYPE']} is an unknown file type") + +def swifter_xr2infile(ds, param, framenum=-1): + """ + Writes a set of Swifter input files from a single frame of a Swiftest xarray dataset + + Parameters + ---------- + ds : xarray dataset + Dataset containing Swifter 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. + param : dict + Swifter input parameters. This method uses the names of the pl and tp files from the input + + Returns + ------- + A set of input files for a Swifter 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']) + param['J2'] = np.double(cb['J_2']) + param['J4'] = np.double(cb['J_4']) + param['RHILL_PRESENT'] = "YES" + + if param['IN_TYPE'] == 'ASCII': + # Swiftest Central body file + plfile = open(param['PL_IN'], 'w') + print(pl.id.count().values + 1, file=plfile) + print(cb.id.values[0], cb['Mass'].values[0], 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) + if param['RHILL_PRESENT'] == "YES": + print(i.values, pli['Mass'].values, pli['Rhill'].values, file=plfile) + else: + print(i.values, pli['Mass'].values, file=plfile) + if param['CHK_CLOSE'] == "YES": + 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(param['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() + else: + # Now make Swiftest files + print(f"{param['IN_TYPE']} is an unknown input file type") + + def swift2swifter(swift_param, plname="", tpname="", conversion_questions={}): swifter_param = {} intxt = conversion_questions.get('RHILL', None) @@ -989,6 +1052,7 @@ def swifter2swiftest(swifter_param, plname="", tpname="", cbname="", conversion_ tpnew = open(swiftest_param['TP_IN'], 'w') except IOError: print(f"Cannot write to file {swiftest_param['TP_IN']}") + return swifter_param print(f"Converting TP file: {swifter_param['TP_IN']} -> {swiftest_param['TP_IN']}") try: @@ -1110,6 +1174,7 @@ def swifter2swiftest(swifter_param, plname="", tpname="", cbname="", conversion_ cbnew.close() except IOError: print(f"Cannot write to file {swiftest_param['CB_IN']}") + return swifter_param MTINY = conversion_questions.get('MTINY', None) if not MTINY: @@ -1151,4 +1216,22 @@ def swift2swiftest(swift_param, plname="", tpname="", cbname="", conversion_ques swiftest_param = swifter2swiftest(swifter_param, plname, tpname, cbname, conversion_questions) swiftest_param['! VERSION'] = "Swiftest parameter file converted from Swift" return swiftest_param - \ No newline at end of file + +def swiftest2swifter_param(swiftest_param, J2=0.0, J4=0.0): + swifter_param = swiftest_param.copy() + CBIN = swifter_param.pop("CB_IN", None) + MTINY = swifter_param.pop("MTINY", None) + MU2KG = swifter_param.pop("MU2KG", 1.0) + DU2M = swifter_param.pop("DU2M", 1.0) + TU2S = swifter_param.pop("TU2S", 1.0) + GR = swifter_param.pop("GR", None) + if GR is not None: + if GR == 'YES': + swifter_param['C'] = swiftest.einsteinC * np.longdouble(TU2S) / np.longdouble(DU2M) + swifter_param['J2'] = J2 + swifter_param['J4'] = J4 + swifter_param['RHILL_PRESENT'] = "YES" + swifter_param['CHK_CLOSE'] = "YES" + swifter_param['! VERSION'] = "Swifter parameter file converted from Swiftest" + + return swifter_param \ No newline at end of file diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 41d6c345f..e20ef05f7 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -48,11 +48,12 @@ def __init__(self, codename="Swiftest", param_file=""): 'YORP': "NO", 'MTINY' : "0.0" } + self.codename = codename if param_file != "" : self.read_param(param_file, codename) return - def add(self, plname, date=date.today().isoformat()): + def add(self, plname, date=date.today().isoformat(), idval=None): """ Adds a solar system body to an existing simulation DataSet. @@ -66,7 +67,7 @@ def add(self, plname, date=date.today().isoformat()): ------- self.ds : xarray dataset """ - self.ds = init_cond.solar_system_horizons(plname, self.param, date, self.ds) + self.ds = init_cond.solar_system_horizons(plname, idval, self.param, date, self.ds) return def read_param(self, param_file, codename="Swiftest"): @@ -84,13 +85,15 @@ def read_param(self, param_file, codename="Swiftest"): self.codename = "Unknown" return - def write_param(self, param_file): + def write_param(self, param_file, param=None): + if param is None: + param = self.param # Check to see if the parameter type matches the output type. If not, we need to convert - codename = self.param['! VERSION'].split()[0] + codename = param['! VERSION'].split()[0] if codename == "Swifter" or codename == "Swiftest": - io.write_labeled_param(self.param, param_file) + io.write_labeled_param(param, param_file) elif codename == "Swift": - io.write_swift_param(self.param, param_file) + io.write_swift_param(param, param_file) else: print('Cannot process unknown code type. Call the read_param method with a valid code name. Valid options are "Swiftest", "Swifter", or "Swift".') return @@ -119,7 +122,9 @@ def convert(self, param_file, newcodename="Swiftest", plname="pl.swiftest.in", t self.param = io.swift2swiftest(self.param, plname, tpname, cbname, conversion_questions) else: goodconversion = False - + else: + goodconversion = False + if goodconversion: self.write_param(param_file) else: @@ -163,7 +168,18 @@ def follow(self, codestyle="Swifter"): print('follow.out written') return fol - def save(self, param_file, framenum=-1): - io.swiftest_xr2infile(self.ds, self.param, framenum) - self.write_param(param_file) + def save(self, param_file, framenum=-1, codename="Swiftest"): + if codename == "Swiftest": + io.swiftest_xr2infile(self.ds, self.param, framenum) + self.write_param(param_file) + elif codename == "Swifter": + if self.codename == "Swiftest": + swifter_param = io.swiftest2swifter_param(self.param) + else: + swifter_param = self.param + io.swifter_xr2infile(self.ds, swifter_param, framenum) + self.write_param(param_file, param=swifter_param) + else: + print(f'Saving to {codename} not supported') + return \ No newline at end of file