Restructuring and removing cruft for clarity and consistency

src/symba/symba_energy_eucl.f90

@@ -21,7 +21,7 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe
 ! internals
    integer(I4B)              :: i, j
    integer(I8B)              :: k
-   real(DP)                  :: rmag, v2, rot2, oblpot, hx, hy, hz
+   real(DP)                  :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz
    real(DP), dimension(npl)  :: irh, kepl, kespinpl, pecb
    real(DP), dimension(npl) :: Lplx, Lply, Lplz
    real(DP), dimension(symba_plA%helio%swiftest%num_plpl_comparisons) :: pepl 
@@ -48,14 +48,18 @@ 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) 
+         hsz = Ip(3,i) * radius(i)**2 * rot(3,i) 
+
          ! Angular momentum from orbit and spin
-         Lplx(i) = mass(i) * (Ip(3,i) * radius(i)**2 * rot(1,i) + hx)
-         Lply(i) = mass(i) * (Ip(3,i) * radius(i)**2 * rot(2,i) + hy)
-         Lplz(i) = mass(i) * (Ip(3,i) * radius(i)**2 * rot(3,i) + hz)
+         Lplx(i) = mass(i) * (hsx + hx)
+         Lply(i) = mass(i) * (hsy + hy)
+         Lplz(i) = mass(i) * (hsz + hz)
 
          ! Kinetic energy from orbit and spin
          kepl(i) = mass(i) * v2
-         kespinpl(i) = mass(i) * Ip(3,i) * radius(i)**2 * rot2
+         kespinpl(i) = mass(i) * (hsx * rot(1, i) + hsy * rot(2, i) + hsz * rot(3, i))
       end do
 
       ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:))

src/symba/symba_frag_pos.f90

@@ -145,12 +145,19 @@ subroutine symba_frag_pos (param, symba_plA, family, x, v, L_spin, Ip, mass, rad
          L_residual(:) = Ltot_before(:) - Ltot_after(:)
          L_spin_frag(:) = L_residual(:) / nfrag
          do i = 1, nfrag
-            rot_frag(:,i) = rot_frag(:,i) + L_spin_frag(:) / (Ip_frag(3, i) * m_frag(i) * rad_frag(i)**2)
+            rot_frag(:,i) = rot_frag(:,i) + L_spin_frag(:) / (Ip_frag(:, i) * m_frag(i) * rad_frag(i)**2)
+         end do
+         L_spin_frag(:) = 0._DP
+         do i = 1, nfrag
+            L_spin_frag(:) = L_spin_frag(:) + Ip_frag(:,i) * rad_frag(i)**2 * m_frag(i) * rot_frag(:,i) 
          end do
       end if
+
+
       call symba_frag_pos_energy_calc(npl, symba_plA, lexclude, ke_after, ke_spin_after, pe_after, Ltot_after,&
          nfrag=nfrag, Ip_frag=Ip_frag, m_frag=m_frag, rad_frag=rad_frag, xb_frag=xb_frag, vb_frag=vb_frag, rot_frag=rot_frag)
          Etot_after = ke_after + ke_spin_after + pe_after
+      L_residual(:) = Ltot_before(:) - Ltot_after(:)
       Lmag_after = norm2(Ltot_after(:))
 
       write(*,        "(' ---------------------------------------------------------------------------')")
@@ -189,13 +196,9 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L
       real(DP), dimension(:,:), intent(out)   :: x_frag, v_frag, rot_frag !! Fragment position, velocities, and spin states
       ! Internals
       real(DP)                                :: mtot, theta, v_frag_norm, r_frag_norm, v_col_norm, r_col_norm
-      real(DP)                                :: b2a, phase_ang
       real(DP), dimension(NDIM)               :: Ltot, delta_r, delta_v
       real(DP), dimension(NDIM)               :: x_col_unit, y_col_unit, z_col_unit
       integer(I4B)                            :: i
-      real(DP), dimension(2,2)                :: orientation
-      real(DP), dimension(2)                  :: frag_vec
-      real(DP), parameter                     :: ecc_ellipse = 0.0_DP ! Eccentricity of the fragment distribution ellipse
 
       ! Now create the fragment distribution
       mtot = sum(mass(:))
@@ -215,17 +218,9 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L
       ! The cross product of the z- by x-axis will give us the y-axis
       call util_crossproduct(y_col_unit, z_col_unit, x_col_unit)
 
-      ! Place the fragments on the collision planea on an ellipse, but with the distance proportional to mass wrt the collisional barycenter
-      ! This gets updated later after the new potential energy is calculated
-      b2a = 1.0_DP / sqrt(1.0_DP - ecc_ellipse**2)
-
       ! The angular spacing of fragments on the ellipse - We will only use half the ellipse
       theta = 2 * PI / nfrag
 
-      ! Adjust the orientation of the ellipse. This mainly affects the disruption case as the first body in the distribution (theta=0) is the largest
-      phase_ang = 0.0_DP
-      orientation = reshape([cos(phase_ang), sin(phase_ang), -sin(phase_ang), cos(phase_ang)], shape(orientation))
-
       ! Re-normalize position and velocity vectors by the fragment number so that for our initial guess we weight each
       ! fragment position by the mass and assume equipartition of energy for the velocity
       r_col_norm = 1.5_DP * 4 * sum(rad_frag(:)) / theta / (nfrag - 1)
@@ -235,14 +230,6 @@ subroutine symba_frag_pos_initialize_fragments(nfrag, xcom, vcom, x, v, L_orb, L
       call random_number(x_frag(:,2:nfrag)) 
 
       x_frag(:, :) = x_frag(:, :) * r_col_norm 
-      ! The rest of the fragments will go in a ring on the x-y plane
-      !do i = 3, nfrag
-      !   frag_vec(:) = [b2a * cos(theta * (i - 1)), sin(theta * (i - 1))]
-      !   frag_vec(:) = matmul(orientation(:,:), frag_vec(:))
-!
-!         r_frag_norm = r_col_norm 
-!         x_frag(:, i) = r_frag_norm * (frag_vec(1) * x_col_unit(:) + frag_vec(2) * y_col_unit(:))
-!      end do
       call symba_frag_pos_com_adjust(xcom, m_frag, x_frag)
 
       call symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, mass, radius, m_frag, rad_frag, Ip_frag, x_frag, v_frag, rot_frag)
@@ -279,7 +266,7 @@ subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, mass, radius
       L_spin_new(:) = L_spin(:,1) + L_spin(:, 2)
       L_spin_frag(:) = L_spin_new(:) / nfrag
       do i = 1, nfrag
-         rot_frag(:,i) = L_spin_frag(:) / (Ip_frag(3, i) * m_frag(i) * rad_frag(i)**2)
+         rot_frag(:,i) = L_spin_frag(:) / (Ip_frag(:, i) * m_frag(i) * rad_frag(i)**2)
       end do
 
       L_orb_frag_mag = norm2(L_orb_old(:)) / nfrag
@@ -293,7 +280,6 @@ subroutine symba_frag_pos_ang_mtm(nfrag, xcom, vcom, L_orb, L_spin, mass, radius
       return 
    end subroutine symba_frag_pos_ang_mtm
 
-
    subroutine symba_frag_pos_kinetic_energy(xcom, vcom, L_orb, m_frag, x_frag, v_frag, ke_target, lmerge)
       !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton
       !!