diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 4e46ad808..e8601f133 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -762,8 +762,14 @@ def swiftest_xr2infile(ds, param, framenum=-1): else: print(i.values, pli['GMass'].values, file=plfile) 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) + if param['IN_FORM'] == 'XV': + print(pli['px'].values, pli['py'].values, pli['pz'].values, file=plfile) + print(pli['vx'].values, pli['vy'].values, pli['vz'].values, file=plfile) + elif param['IN_FORM'] == 'EL': + print(pli['a'].values, pli['e'].values, pli['inc'].values, file=plfile) + print(pli['capom'].values, pli['omega'].values, pli['capm'].values, file=plfile) + else: + print(f"{param['IN_FORM']} is not a valid input format type.") plfile.close() # TP file @@ -772,8 +778,14 @@ def swiftest_xr2infile(ds, param, framenum=-1): 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) + if param['IN_FORM'] == 'XV': + print(tpi['px'].values, tpi['py'].values, tpi['pz'].values, file=tpfile) + print(tpi['vx'].values, tpi['vy'].values, tpi['vz'].values, file=tpfile) + elif param['IN_FORM'] == 'EL': + print(tpi['a'].values, tpi['e'].values, tpi['inc'].values, file=tpfile) + print(tpi['capom'].values, tpi['omega'].values, tpi['capm'].values, file=tpfile) + else: + print(f"{param['IN_FORM']} is not a valid input format type.") tpfile.close() elif param['IN_TYPE'] == 'REAL8': # Now make Swiftest files @@ -788,23 +800,33 @@ def swiftest_xr2infile(ds, param, framenum=-1): plfile = FortranFile(param['PL_IN'], 'w') npl = pl.id.count().values plid = pl.id.values - px = pl['px'].values - py = pl['py'].values - pz = pl['pz'].values - vx = pl['vx'].values - vy = pl['vy'].values - vz = pl['vz'].values + if param['IN_FORM'] == 'XV': + v1 = pl['px'].values + v2 = pl['py'].values + v3 = pl['pz'].values + v4 = pl['vx'].values + v5 = pl['vy'].values + v6 = pl['vz'].values + elif param['IN_FORM'] == 'EL': + v1 = pl['a'].values + v2 = pl['e'].values + v3 = pl['inc'].values + v4 = pl['capom'].values + v5 = pl['omega'].values + v6 = pl['capm'].values + else: + print(f"{param['IN_FORM']} is not a valid input format type.") Gmass = pl['GMass'].values radius = pl['Radius'].values plfile.write_record(npl) plfile.write_record(plid) - plfile.write_record(px) - plfile.write_record(py) - plfile.write_record(pz) - plfile.write_record(vx) - plfile.write_record(vy) - plfile.write_record(vz) + plfile.write_record(v1) + plfile.write_record(v2) + plfile.write_record(v3) + plfile.write_record(v4) + plfile.write_record(v5) + plfile.write_record(v6) plfile.write_record(Gmass) if param['RHILL_PRESENT'] == 'YES': rhill = pl['Rhill'].values @@ -814,20 +836,30 @@ def swiftest_xr2infile(ds, param, framenum=-1): tpfile = FortranFile(param['TP_IN'], 'w') ntp = tp.id.count().values tpid = tp.id.values - px = tp['px'].values - py = tp['py'].values - pz = tp['pz'].values - vx = tp['vx'].values - vy = tp['vy'].values - vz = tp['vz'].values + if param['IN_FORM'] == 'XV': + v1 = tp['px'].values + v2 = tp['py'].values + v3 = tp['pz'].values + v4 = tp['vx'].values + v5 = tp['vy'].values + v6 = tp['vz'].values + elif param['IN_FORM'] == 'EL': + v1 = tp['a'].values + v2 = tp['e'].values + v3 = tp['inc'].values + v4 = tp['capom'].values + v5 = tp['omega'].values + v6 = tp['capm'].values + else: + print(f"{param['IN_FORM']} is not a valid input format type.") tpfile.write_record(ntp) tpfile.write_record(tpid) - tpfile.write_record(px) - tpfile.write_record(py) - tpfile.write_record(pz) - tpfile.write_record(vx) - tpfile.write_record(vy) - tpfile.write_record(vz) + tpfile.write_record(v1) + tpfile.write_record(v2) + tpfile.write_record(v3) + tpfile.write_record(v4) + tpfile.write_record(v5) + tpfile.write_record(v6) else: print(f"{param['IN_TYPE']} is an unknown file type") diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 6b737d076..c602cbe8c 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -20,6 +20,7 @@ def __init__(self, codename="Swiftest", param_file=""): 'TP_IN': "tp.in", 'CB_IN': "cb.in", 'IN_TYPE': "ASCII", + 'IN_FORM': "XV", 'ISTEP_OUT': "1", 'ISTEP_DUMP': "1", 'BIN_OUT': "bin.dat", diff --git a/src/io/io.f90 b/src/io/io.f90 index 670d98db8..16af0ca10 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -13,7 +13,6 @@ module subroutine io_conservation_report(self, param, lterminal) logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen ! Internals real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now - real(DP), dimension(NDIM), save :: Ltot_last, Lorbit_last, Lspin_last real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now real(DP) :: Eorbit_error, Etotal_error, Ecoll_error real(DP) :: GMtot_now @@ -68,7 +67,7 @@ module subroutine io_conservation_report(self, param, lterminal) Ecoll_error = param%Ecollisions / abs(param%Eorbit_orig) Etotal_error = (Eorbit_now - param%Ecollisions - param%Eorbit_orig - param%Euntracked) / abs(param%Eorbit_orig) Merror = (GMtot_now - param%GMtot_orig) / param%GMtot_orig - write(*, egytermfmt) Lerror, Ecoll_error, Etotal_error, Merror + write(*, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror end if end associate @@ -800,7 +799,7 @@ module subroutine io_param_writer(self, unit, iotype, v_list, iostat, iomsg) write(unit, '("LSPIN_ORIG ",3(1X,ES25.17))') param%Lspin_orig(:) write(unit, '("LESCAPE ",3(1X,ES25.17))') param%Lescape(:) - write(param_name, Afmt) "GMescape"; write(param_value, Rfmt) param%GMescape; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) + write(param_name, Afmt) "GMESCAPE"; write(param_value, Rfmt) param%GMescape; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "ECOLLISIONS"; write(param_value, Rfmt) param%Ecollisions; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) write(param_name, Afmt) "EUNTRACKED"; write(param_value, Rfmt) param%Euntracked; write(unit, Afmt, err = 667, iomsg = iomsg) adjustl(param_name), adjustl(param_value) end if diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index f08e8eb05..1ee3bdf80 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -775,6 +775,19 @@ module pure subroutine orbel_xv2aqt(mu, x, v, a, q, capm, tperi) real(DP), intent(out) :: tperi !! time of pericenter passage end subroutine orbel_xv2aqt + module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) + implicit none + real(DP), intent(in) :: mu !! Gravitational constant + real(DP), dimension(:), intent(in) :: x !! Position vector + real(DP), dimension(:), intent(in) :: v !! Velocity vector + real(DP), intent(out) :: a !! semimajor axis + real(DP), intent(out) :: e !! eccentricity + real(DP), intent(out) :: inc !! inclination + real(DP), intent(out) :: capom !! longitude of ascending node + real(DP), intent(out) :: omega !! argument of periapsis + real(DP), intent(out) :: capm !! mean anomaly + end subroutine orbel_xv2el + module subroutine orbel_xv2el_vec(self, cb) implicit none class(swiftest_body), intent(inout) :: self !! Swiftest body object diff --git a/src/orbel/orbel.f90 b/src/orbel/orbel.f90 index f1ab88825..c07200b14 100644 --- a/src/orbel/orbel.f90 +++ b/src/orbel/orbel.f90 @@ -890,7 +890,7 @@ module subroutine orbel_xv2el_vec(self, cb) end do end subroutine orbel_xv2el_vec - pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) + module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) !! author: David A. Minton !! !! Compute osculating orbital elements from relative Cartesian position and velocity @@ -906,9 +906,17 @@ pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2el.f90 !! Adapted from Martin Duncan's Swift routine orbel_xv2el.f implicit none - real(DP), intent(in) :: mu - real(DP), dimension(:), intent(in) :: x, v - real(DP), intent(out) :: a, e, inc, capom, omega, capm + ! Arguments + real(DP), intent(in) :: mu !! Gravitational constant + real(DP), dimension(:), intent(in) :: x !! Position vector + real(DP), dimension(:), intent(in) :: v !! Velocity vector + real(DP), intent(out) :: a !! semimajor axis + real(DP), intent(out) :: e !! eccentricity + real(DP), intent(out) :: inc !! inclination + real(DP), intent(out) :: capom !! longitude of ascending node + real(DP), intent(out) :: omega !! argument of periapsis + real(DP), intent(out) :: capm !! mean anomaly + ! Internals integer(I4B) :: iorbit_type real(DP) :: r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, cape, tmpf, capf real(DP), dimension(NDIM) :: hvec diff --git a/src/python/orbel.f90 b/src/python/orbel.f90 new file mode 100644 index 000000000..4dd26a68d --- /dev/null +++ b/src/python/orbel.f90 @@ -0,0 +1,32 @@ +module orbel + use swiftest + private + public :: xv2el +contains + pure elemental subroutine xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm) + use module_interfaces, only : orbel_xv2el + implicit none + ! Arguments + real*8, intent(in) :: mu, px, py, pz, vx, vy, vz + real*8, intent(out) :: a, e, inc, capom, omega, capm + !$f2py intent(in) mu + !$f2py intent(in) px + !$f2py intent(in) py + !$f2py intent(in) pz + !$f2py intent(in) vx + !$f2py intent(in) vy + !$f2py intent(in) vz + !$f2py intent(out) a + !$f2py intent(out) e + !$f2py intent(out) inc + !$f2py intent(out) capom + !$f2py intent(out) omega + !$f2py intent(out) capm + ! Internals + real*8, dimension(3) :: x, v + x = [px, py, pz] + v = [vx, vy, vz] + call orbel_xv2el(x, v, mu, a, e, inc, capom, omega, capm) + return + end subroutine xv2el +end module orbel \ No newline at end of file diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 193461bc7..acd233cba 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -179,14 +179,21 @@ subroutine symba_discard_nonplpl(pl, system, param) class(symba_pl), intent(inout) :: pl !! SyMBA test particle object class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + ! Internals + logical, dimension(pl%nbody) :: ldiscard ! First check for collisions with the central body associate(npl => pl%nbody, cb => system%cb) if (npl == 0) return + ldiscard(1:npl) = pl%ldiscard(1:npl) ! Don't include any bodies that were previously flagged for discard in here if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) then call symba_discard_cb_pl(pl, system, param) end if if (param%qmin >= 0.0_DP .and. npl > 0) call symba_discard_peri_pl(pl, system, param) + if (any(.not.ldiscard(1:npl) .and. pl%ldiscard(1:npl))) then + ldiscard(1:npl) = .not.ldiscard(1:npl) .and. pl%ldiscard(1:npl) + call system%pl_discards%append(pl, ldiscard) + end if end associate return diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 1c0981232..634eb1083 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -298,7 +298,6 @@ module subroutine symba_io_read_particle(system, param) end if end if if (.not.lmatch) then - write(*,*) 'Particle id ',id,' not found. Skipping' read(LUN, err = 667, iomsg = errmsg) tmpinfo end if end do