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

@@ -1,38 +1,33 @@
-!
-! 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
 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

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

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(:))