From 3afb39fb3899a5882a2a1f2db8e94fad0888054b Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 7 Jul 2021 13:51:25 -0400 Subject: [PATCH] Updated swiftest package to include new methods to save swifter files from a swiftest run --- .../helio_swifter_comparison/.idea/.gitignore | 3 + .../helio_swifter_comparison/cb.swiftest.in | Bin 64 -> 87 bytes .../helio_swifter_comparison/init_cond.py | 355 +++--------------- .../helio_swifter_comparison/param.swifter.in | 62 +-- .../param.swiftest.in | 64 ++-- .../helio_swifter_comparison/pl.swifter.in | 56 ++- .../helio_swifter_comparison/pl.swiftest.in | Bin 700 -> 1525 bytes .../helio_swifter_comparison/tp.swifter.in | 24 +- .../helio_swifter_comparison/tp.swiftest.in | Bin 280 -> 565 bytes python/swiftest/swiftest/init_cond.py | 65 +++- python/swiftest/swiftest/io.py | 89 ++++- python/swiftest/swiftest/simulation_class.py | 36 +- 12 files changed, 317 insertions(+), 437 deletions(-) create mode 100644 examples/helio_swifter_comparison/.idea/.gitignore 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 7e13534c0bc881a8cc6e551b3a47dbfdd785d49e..058975b813b2ede68061d83c64275e30fd103ddb 100644 GIT binary patch literal 87 zcmWN`u@L|<2m`^KUd&(t5)}0Px9|=w*~|437p$0B5w!4#V!s5&61QdL>g)+_&4sf2 RI~R7~DCJW7wlL81u^;=Z5pDng literal 64 ncmd;JU|=xH*zksXud@ROkPX6j{SWxW@f6#`_&^~rz{UpvSndZF 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 6ee1ca38380823cc94be008e9a26b0914f75ea10..c1f0be29eec4f9cd2343b9b942ef8f00390f7b88 100644 GIT binary patch literal 1525 zcmXYx$&uqg3^^I*S6)F=es~UOpa+g_*De*baT)*tzr^;R3cXKpPv3{!0D@aRU6CRZjT#JnN z2^F1o(=Dnrs41(-pt7YI=cu<^CL*`43(}F5DJt zaipPSpao$Zw-Cf==la=$MGV%N$Xhg;&rYc@)*DJAi zEda3GAY{?`U{Iuqf=p=(m4e9@8ys&LmU98DG6`G}1aOWqCbrQk7f(HI<{@=WX7=vv|1Q6ju{Guibr1906QUmc~W~EA~=@; z1`syvJU-x>;&KL`hlEjmiznz*Rg*w$SO(kItE*^}pX#xbft?(G;FV+^0oo0G%3dl~ zQasZNnG?(hq8Rj6d$9X3D~m^22vGQ|KKRfkSutl)PUl5e*c;tYl=RrlmDd%LZHy=A zu`tE*_vH(D49Om0tg+N!+v|805yaujrl59~9=(ly?9aL+&apR5IQ!$WWz(|$V7#My zjEuRlD?0+mqLt6EaLIh20WJaDqJaQDc1cx7XpWU(;!s7mFmEH>8`aanEre)~X> PntQBRET&EiN*(_Jhiomc literal 700 zcmd;JU|`?`Vi4c}Vih1}0%B$$W&vVWAZ7z%b|41HgVb^XF-V;U5bt;@Zu%sCm;H}j z2VS4>SZ`lcKjF%y)W`e3*-hm+Ip?DNwaxP<9U{S5LBDY21z|Sw!Ui=f*bV$fm&bFQ|a$wKroGHdC_6K0*$?0jR%c=C* zi_dpUT`{$LKj#na2sR6*{WrFn-uSm+rv00QXUi-c*6d#@CpvLy{FePW1{<-<%tkJ zpv?0oN9hTJ!xs0j`IC%(+F#_f2uRL`cB`-Y*8bKeevfZAH0`g(t=2zx>X7|6 z{>cY>)z|EwVYmC~d%3Oqzk1EsI7xnjeZV>sE34)?`?DS^O1Y(6h5Pq76C>xpn%i~@ z&s?7QAw<$%$|sudewn=eI%yBSh1*}*Jvvb(oaVFG{`qx*tNgO<_AfN=Ztk&-v2VJP zA$n>>g8f}@CgmgZC)>gNyLZlMHZ#eOb^)1YMRQGA?S;%I6f%ae*;_|!x8yza-EP+{ tvw3+7?Ckrdo$OoMWN9B0!S1-oLd`yW#(|LKlT_`au1r>HTk!yF9su;c8`b~- 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 5fe4b5618a2c2358d645755a4759220addef5211..62acd79fc49307ef00eac51b85a99ee95f19be06 100644 GIT binary patch literal 565 zcmXYu%Z&vw2n6qrViI^3&@`X?6Wg}5nj2%Fy2#_^CoUWaA_oAAll0GTp={V@)3olL zitlF^owOb74uNb(^ZepQa;7VbT^n>;d~pXowHs~4g7-{(c4^3pMpTc+d0LlP;^3Ms zN#_JM#?<=#{JO)un1RG*ujafl?j7zX9uI4ygU4=L5~U~tP2UV!2}NrrwRfuMT@$P~ zddwMsuZc}RLm)T`%5YWE3j#9MmVZB7$10?PcC`RDyoJ&b8#ykE2KL6jFp6RDlDllK zOlgE&!O}u%6w0k#HuKeOYg5J>q+JuwOueCyr?i0OIFg35GrFBDNcDn~-YzQZs tDWBB0=G&V&$MxEuTfiXi4|HFol^Ca2t2Wvro`0OseP%e%sIDct^9P%YZ!iD= literal 280 zcmd;JU|?VYVi4c}VgVpt*v81P2#6O0@e&{gi75c_`sW|MZct};_*xaLG|!ybp~pvG z<$=+k{d{L`h(t&*9ssEU;nUN$%=Jh(Wgm6;mVa~;!+|v`d5*qPU_0<4Fo%IT@q;}` z4G7QKbyQ3C_JjRfR`!3F?*43Ve0gitifPaGFLUKVo}uFl9Jvog5r`u=G(%BQ3r z|F_>Fy 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