Had to renormalise plm1. Acceleration matches swiftest_obl.f90

src/swiftest/swiftest_sph.f90

@@ -68,7 +68,7 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
             ! call PlmBar_d1(p, p_deriv, l_max, cos(theta))      ! Associated Legendre Polynomials and the 1st Derivative
             call PlmBar(p, l_max, cos(theta))
         else
-            call PlmBar_d1(p, p_deriv, l_max, 0.0_DP) 
+            call PlmBar(p, l_max, 0.0_DP) 
         end if
 
         do l = 1, l_max ! skipping the l = 0 term; It is the spherical body term
@@ -96,32 +96,38 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
                 ! g_sph(3) = g_sph(3) - GMcb * r_0**l / r_mag**(l + 2) * ccss * (-1 * dplm * sin(theta)**2  &
                 !                                                         + plm * (l + 1) * cos(theta))          ! g_z
 
-                ! cos_tmp = cos(theta)
-                ! sin_tmp = sin(theta)
-                ! sin2_tmp = sin(2 * theta)
-                ! sin_phi = sin(phi)
-                ! cos_phi = cos(phi)
+                cos_tmp = cos(theta)
+                sin_tmp = sin(theta)
+                sin2_tmp = sin(2 * theta)
+                sin_phi = sin(phi)
+                cos_phi = cos(phi)
                 ! g_sph(:) = 0.0_DP 
-
-                ! !! Condensed form to reduce floating point error
-                ! ! fac0 = -GMcb / r_mag**2                                                        
-                ! fac1 = cssc * m * plm / sin(theta)
-                ! fac2 = ccss * sin(theta) 
-                ! fac3 = dplm * cos(theta) + plm * (l + 1)
-                ! r_fac = -GMcb * r_0**l / r_mag**(l + 2)
-
-                ! g_sph(1) = g_sph(1) + r_fac * (fac1 * sin(phi) + fac2 * fac3 * cos(phi))
-                ! g_sph(2) = g_sph(2) + r_fac * (fac1 * cos(phi) + fac2 * fac3 * sin(phi))
-                ! g_sph(3) = g_sph(3) + r_fac * ccss * (fac3 * cos(theta) - dplm)
-                                ! (-1 * dplm * sin(theta)**2  + plm * (l + 1) * cos(theta))          ! g_z
-
-                ! Condensed form with alternative form for g_sph
+
+                ! Alternative form for g_sph
+
                 if ((m+1) .le. l) then
                     lmindex = PlmIndex(l, m+1) 
-                    plm1 = p(lmindex)
+                    plm1 = p(lmindex) 
+                    if(m .eq. 0) then
+                        plm1 = plm1 * sqrt(((l + m + 1) * (l - m)) / 2.0) ! renormalize plm1 to the norm of plm
+                    else 
+                        plm1 = plm1 * sqrt((l + m + 1) * (l - m) * 1.0)       ! renormalize plm1 to the norm of plm
+                    end if
                 else
                     plm1 = 0.0_DP  
                 end if 
+
+                ! !! Alternative form of dplm
+
+                ! g_sph(1) = g_sph(1) - GMcb * r_0**l / r_mag**(l + 2) * (cssc * m / sin(theta) * plm * sin(phi) &
+                !                                                        + ccss * cos(phi) * (plm * ((l + m + 1) * sin(theta) - m / sin(theta)) & 
+                !                                                                              + plm1 * cos(theta)))
+                ! g_sph(2) = g_sph(2) - GMcb * r_0**l / r_mag**(l + 2) * (-cssc * m / sin(theta) * plm * cos(phi) &
+                !                                                         + ccss * sin(phi) * (plm * ((l + m + 1) * sin(theta) - m / sin(theta)) & 
+                !                                                                              + plm1 * cos(theta))) 
+                ! g_sph(3) = g_sph(3) - GMcb * r_0**l / r_mag**(l + 2) * ccss * (plm * (l + m +1) * cos(theta) - plm1 * sin(theta))          
+
+                ! Condensed form
 
                 ! fac0 = -(m * cos_tmp * plm / sin_tmp - plm1) / sin_tmp ! dplm
                 fac1 = m / sin(theta) * plm
@@ -135,6 +141,9 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
                 g_sph(1) = g_sph(1) + r_fac * (cssc * fac1 * sin(phi) + ccss * (fac2 - fac1) * cos(phi))
                 g_sph(2) = g_sph(2) + r_fac * (-cssc * fac1 * cos(phi) + ccss * (fac2 - fac1) * sin(phi))
                 g_sph(3) = g_sph(3) + r_fac * ccss * (plm * (l + m + 1) * cos(theta) - plm1 * sin(theta))
+
+                fac0 = (.mag. g_sph(:))
+                ! fac3 = 3 * c_lm(m+1, l+1, 1) / 2 * r_0**2 / r_mag**4 * (3*(cos(theta))**2 - 1) * GMcb ! g_sph for J2
 
             end do
         end do

swiftest/sph_harmonics.py

@@ -72,7 +72,7 @@ def clm_from_ellipsoid(mass, density, a, b = None, c = None, lmax = 6, ref_radiu
 
         # get gravity coefficients
         clm_class = pysh.SHGravCoeffs.from_shape(shape_SH, rho = density, gm = Gmass) # 4pi normalization
-        clm = clm_class.to_array(normalization = '4pi') * 4 * np.pi # export as array with 4pi normalization and scaling by 4*pi to match normalisation
+        clm = clm_class.to_array(normalization = '4pi') # export as array with 4pi normalization and not scaling by 4*pi to match normalisation
 
         # Return reference radius EQUALS the radius of the Central Body
         print(f'Ensure that the Central Body radius equals the reference radius.')