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

Commit

Permalink
Merge branch 'debug'
Browse files Browse the repository at this point in the history
  • Loading branch information
@daminton
daminton committed Aug 19, 2021
2 parents 75119bd + 47d879c commit c6e00da
Show file tree
Hide file tree
Showing 8 changed files with 127 additions and 36 deletions.
88 changes: 60 additions & 28 deletions python/swiftest/swiftest/io.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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")

Expand Down
1 change: 1 addition & 0 deletions python/swiftest/swiftest/simulation_class.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
5 changes: 2 additions & 3 deletions src/io/io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
13 changes: 13 additions & 0 deletions src/modules/swiftest_classes.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
16 changes: 12 additions & 4 deletions src/orbel/orbel.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
32 changes: 32 additions & 0 deletions src/python/orbel.f90
View file
Original file line number Diff line number Diff line change
@@ -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
7 changes: 7 additions & 0 deletions src/symba/symba_discard.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 0 additions & 1 deletion src/symba/symba_io.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit c6e00da

Please sign in to comment.