From c044922f13e3bc9cea5cf673102aed1ecf693fef Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 21 May 2021 13:37:35 -0400 Subject: [PATCH] Fixed some issues related to spin kinetic energy and created a new test case for troubleshooting energy gain when colliding with the central body (sun) --- examples/symba_energy_momentum/param.sun.in | 25 +++++++++------------ examples/symba_energy_momentum/sun.in | 24 ++++++++++---------- src/symba/symba_energy_eucl.f90 | 7 +++--- 3 files changed, 26 insertions(+), 30 deletions(-) diff --git a/examples/symba_energy_momentum/param.sun.in b/examples/symba_energy_momentum/param.sun.in index 65b97fd77..d1408253b 100644 --- a/examples/symba_energy_momentum/param.sun.in +++ b/examples/symba_energy_momentum/param.sun.in @@ -1,14 +1,12 @@ -! -! Parameter file for Mars in SI units -! T0 0.0e0 -TSTOP 0.01 ! simulation length in seconds = 100 years -DT 0.01 ! stepsize in seconds +TSTOP 0.00381 ! simulation length in seconds = 100 years +DT 0.00001 ! stepsize in seconds PL_IN sun.in TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ! output cadence every year BIN_OUT bin.sun.dat +PARTICLE_FILE particle.sun.dat OUT_TYPE REAL8 ! double precision real output OUT_FORM XV ! osculating element output OUT_STAT REPLACE @@ -16,23 +14,20 @@ ISTEP_DUMP 1 ! system dump cadence J2 0.0 ! no J2 term J4 0.0 ! no J4 term CHK_CLOSE yes ! check for planetary close encounters -CHK_RMIN 3389.5e3 -CHK_RMAX 3389.5e6 +CHK_RMIN 0.005 +CHK_RMAX 1e6 CHK_EJECT -1.0 ! ignore this check CHK_QMIN -1.0 ! ignore this check !CHK_QMIN_COORD HELIO ! commented out here !CHK_QMIN_RANGE 1.0 1000.0 ! commented out here -ENC_OUT enc.merger.dat +ENC_OUT enc.sun.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT yes ! Hill's sphere radii in input file -MTINY 1.0e-6 +MTINY 1.0e-16 FRAGMENTATION yes -MU2KG 1.0 -TU2S 1.0 -DU2M 1.0 -!MU2KG 1.98908e30 -!TU2S 3.1556925e7 -!DU2M 1.49598e11 +MU2KG 1.98908e30 +TU2S 3.1556925e7 +DU2M 1.49598e11 ENERGY yes ROTATION yes diff --git a/examples/symba_energy_momentum/sun.in b/examples/symba_energy_momentum/sun.in index 206e45e58..8359affdb 100644 --- a/examples/symba_energy_momentum/sun.in +++ b/examples/symba_energy_momentum/sun.in @@ -1,18 +1,18 @@ -3 ! Mars System in SI units -1 4.28388662e+13 -.0 .0 .0 ! x y z -.0 .0 .0 !vx vy vz +2 +1 39.47841760435743 +0.0 0.0 0.0 +0.0 0.0 0.0 0.4 0.4 0.4 !Ip 0.0 0.0 0.0 !rot -2 3.87032350e+03 2.97272135e+03 ! particle number mass Rhill -1.62942641e+03 !particle radius in AU --3.50100037e+05 -8.90344807e+05 -1.16096189e+03 ! x y z -1.96968977e+01 -7.69668748e+01 1.90000055e-01 ! vx vy vz +2 1e-07 0.0009 +7e-06 +0.1 0.0 0.0 +-10.00 2.00 0.0 0.4 0.4 0.4 !Ip 0.0 0.0 0.0 !rot -3 9.35875233e+04 8.60656833e+03 ! particle number mass Rhill -6.94344459e+03 !particle radius in AU --3.50432670e+06 -8.90649894e+06 -6.13283693e+03 ! x y z -1.96781506e+03 -7.75898328e+02 -1.55165450e-01 ! vx vy vz +3 1e-08 0.0004 +3.25e-06 +1.0 4.20E-05 0.0 +0.00 -6.28 0.0 0.4 0.4 0.4 !Ip 0.0 0.0 0.0 !rot diff --git a/src/symba/symba_energy_eucl.f90 b/src/symba/symba_energy_eucl.f90 index 68d6eb6b4..3dc691295 100644 --- a/src/symba/symba_energy_eucl.f90 +++ b/src/symba/symba_energy_eucl.f90 @@ -48,8 +48,9 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe hy = xb(3,i) * vb(1,i) - xb(1,i) * vb(3,i) hz = xb(1,i) * vb(2,i) - xb(2,i) * vb(1,i) - hsx = Ip(1,i) * radius(i)**2 * rot(1,i) - hsy = Ip(2,i) * radius(i)**2 * rot(2,i) + ! For simplicity, we always assume that the rotation pole is the 3rd principal axis + hsx = Ip(3,i) * radius(i)**2 * rot(1,i) + hsy = Ip(3,i) * radius(i)**2 * rot(2,i) hsz = Ip(3,i) * radius(i)**2 * rot(3,i) ! Angular momentum from orbit and spin @@ -59,7 +60,7 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe ! Kinetic energy from orbit and spin kepl(i) = mass(i) * v2 - kespinpl(i) = mass(i) * (hsx * rot(1, i) + hsy * rot(2, i) + hsz * rot(3, i)) + kespinpl(i) = mass(i) * Ip(3, i) * radius(i)**2 * rot2 end do ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:))