Skip to content
This repository was archived by the owner on Aug 28, 2024. It is now read-only.

Commit

Permalink
Fixed bug that was causing test particle encounters to be computed in…
Browse files Browse the repository at this point in the history
…correctly in SyMBA. Fixed another initial conditions file
  • Loading branch information
daminton committed Sep 28, 2021
1 parent b3a897d commit c20c368
Show file tree
Hide file tree
Showing 11 changed files with 103 additions and 105 deletions.

Large diffs are not rendered by default.

Binary file not shown.
108 changes: 48 additions & 60 deletions examples/symba_swifter_comparison/1pl_1tp_encounter/init_cond.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
swiftest_pl = "pl.swiftest.in"
swiftest_tp = "tp.swiftest.in"
swiftest_cb = "cb.swiftest.in"
swiftest_bin = "bin.swiftest.dat"
swiftest_bin = "bin.swiftest.nc"
swiftest_enc = "enc.swiftest.dat"
swiftest_dis = "discard.swiftest.dat"

Expand All @@ -43,11 +43,16 @@
rmax = 1000.0

npl = 1
plid = 2
tpid = 100
ntp = 1
plid = 1
tpid = 2

cbname = "Sun"
plname = "Planet"
tpname = "TestParticle"

radius = np.double(4.25875607065041e-05)
mass = np.double(0.00012002693582795244940133)
Gmass = np.double(0.00012002693582795244940133)
apl = np.longdouble(1.0)
atp = np.longdouble(1.01)
vpl = np.longdouble(2 * np.pi)
Expand All @@ -62,23 +67,21 @@
rhill = np.double(apl * 0.0100447248332378922085)

#Make Swifter files
plfile = open(swifter_pl, 'w')
print(npl+1, f'! Planet input file generated using init_cond.py',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)
print(plid,"{:.23g}".format(mass),rhill, file=plfile)
print(radius, file=plfile)
print(*p_pl, file=plfile)
print(*v_pl, file=plfile)
plfile.close()

tpfile = open(swifter_tp, 'w')
print(1,file=tpfile)
print(tpid, file=tpfile)
print(*p_tp, file=tpfile)
print(*v_tp, file=tpfile)
tpfile.close()
with open(swifter_pl, 'w') as plfile:
print(npl+1, f'! Planet input file generated using init_cond.py',file=plfile)
print(0,GMSun,file=plfile)
print('0.0 0.0 0.0',file=plfile)
print('0.0 0.0 0.0',file=plfile)
print(plid,"{:.23g}".format(Gmass),rhill, file=plfile)
print(radius, file=plfile)
print(*p_pl, file=plfile)
print(*v_pl, file=plfile)

with open(swifter_tp, 'w') as tpfile:
print(1,file=tpfile)
print(tpid, file=tpfile)
print(*p_tp, file=tpfile)
print(*v_tp, file=tpfile)

sys.stdout = open(swifter_input, "w")
print(f'! Swifter input file generated using init_cond.py')
Expand Down Expand Up @@ -110,41 +113,25 @@
sys.stdout = sys.__stdout__

#Now make Swiftest files
cbfile = FortranFile(swiftest_cb, 'w')
Msun = np.double(1.0)
cbfile.write_record(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 = FortranFile(swiftest_pl, 'w')
plfile.write_record(npl)

plfile.write_record(plid)
plfile.write_record(p_pl[0])
plfile.write_record(p_pl[1])
plfile.write_record(p_pl[2])
plfile.write_record(v_pl[0])
plfile.write_record(v_pl[1])
plfile.write_record(v_pl[2])
plfile.write_record(mass)
plfile.write_record(rhill)
plfile.write_record(radius)
plfile.close()
tpfile = FortranFile(swiftest_tp, 'w')
ntp = 1
tpfile.write_record(ntp)
tpfile.write_record(tpid)
tpfile.write_record(p_tp[0])
tpfile.write_record(p_tp[1])
tpfile.write_record(p_tp[2])
tpfile.write_record(v_tp[0])
tpfile.write_record(v_tp[1])
tpfile.write_record(v_tp[2])

tpfile.close()
with open(swiftest_cb, 'w') as cbfile:
print(cbname,file=cbfile)
print(GMSun, file=cbfile)
print(rmin, file=cbfile)
print(J2, file=cbfile)
print(J4, file=cbfile)

with open(swiftest_pl, 'w') as plfile:
print(npl, f'! Planet input file generated using init_cond.py',file=plfile)
print(plname,"{:.23g}".format(Gmass),rhill, file=plfile)
print(radius, file=plfile)
print(*p_pl, file=plfile)
print(*v_pl, file=plfile)

with open(swiftest_tp, 'w') as tpfile:
print(ntp,file=tpfile)
print(tpname, file=tpfile)
print(*p_tp, file=tpfile)
print(*v_tp, file=tpfile)

sys.stdout = open(swiftest_input, "w")
print(f'! Swiftest input file generated using init_cond.py')
Expand All @@ -154,12 +141,13 @@
print(f'CB_IN {swiftest_cb}')
print(f'PL_IN {swiftest_pl}')
print(f'TP_IN {swiftest_tp}')
print(f'IN_TYPE REAL8')
print(f'IN_TYPE ASCII')
print(f'IN_FORM XV')
print(f'ISTEP_OUT {iout:d}')
print(f'ISTEP_DUMP {iout:d}')
print(f'ISTEP_DUMP {10*iout:d}')
print(f'BIN_OUT {swiftest_bin}')
print(f'OUT_TYPE REAL8')
print(f'OUT_FORM XV')
print(f'OUT_TYPE NETCDF_DOUBLE')
print(f'OUT_FORM XVEL')
print(f'OUT_STAT REPLACE')
print(f'CHK_CLOSE yes')
print(f'CHK_RMIN {rmin}')
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,13 @@ DT 0.0006844626967830253
CB_IN cb.swiftest.in
PL_IN pl.swiftest.in
TP_IN tp.swiftest.in
IN_TYPE REAL8
IN_TYPE ASCII
IN_FORM XV
ISTEP_OUT 1
ISTEP_DUMP 1
BIN_OUT bin.swiftest.dat
OUT_TYPE REAL8
OUT_FORM XV
ISTEP_DUMP 10
BIN_OUT bin.swiftest.nc
OUT_TYPE NETCDF_DOUBLE
OUT_FORM XVEL
OUT_STAT REPLACE
CHK_CLOSE yes
CHK_RMIN 0.004650467260962157
Expand All @@ -30,3 +31,4 @@ DU2M 149597870700.0
TU2S 31557600.0
RHILL_PRESENT yes
GMTINY 1e-12
INTERACTION_LOOPS TRIANGULAR
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
2 ! Planet input file generated using init_cond.py
1 39.476926408897625196
0 39.476926408897625196
0.0 0.0 0.0
0.0 0.0 0.0
2 0.00012002693582795244940133 0.010044724833237892
1 0.00012002693582795244940133 0.010044724833237892
4.25875607065041e-05
1.0 0.0 0.0
0.0 6.283185307179586 0.0
Binary file not shown.
Loading

0 comments on commit c20c368

Please sign in to comment.