Changed to 4pi normalisation. Scaled C_lm in python by 4pi for correc…

…t C00 terms.UNSURE, STILL TESTING
MintonGroup · Nov 10, 2023 · 36bc720 · 36bc720
1 parent 6b2c496
commit 36bc720
Show file tree

Hide file tree

Showing 2 changed files with 5 additions and 3 deletions.
diff --git a/src/swiftest/swiftest_sph.f90 b/src/swiftest/swiftest_sph.f90
@@ -62,9 +62,11 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
         t3 = 1.875_DP * j4rp4 * rinv2
         fac1 = t0 * (t1 - t3 - (5 * t1 - (14.0_DP - 21.0_DP * t2) * t3) * t2)
         fac2 = 2 * t0 * (t1 - (2.0_DP - (14.0_DP * t2 / 3.0_DP)) * t3)
+        fac0 = 4 * PI
 
         if(cos(theta) > epsilon(0.0_DP)) then
-            call PlmBar_d1(p, p_deriv, l_max, cos(theta))      ! Unnormalized Associated Legendre Polynomials and the 1st Derivative
+            ! 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) 
         end if
@@ -75,7 +77,7 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
                 ! Associated Legendre Polynomials 
                 lmindex = PlmIndex(l, m)  
                 plm = p(lmindex)                ! p_l,m
-                dplm = p_deriv(lmindex)         ! d(p_l,m)
+                ! dplm = p_deriv(lmindex)         ! d(p_l,m)
 
                 ! C_lm and S_lm with Cos and Sin of m * phi
                 ccss = c_lm(m+1, l+1, 1) * cos(m * phi) & 

diff --git a/swiftest/sph_harmonics.py b/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') # export as array with 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
 
         # Return reference radius EQUALS the radius of the Central Body
         print(f'Ensure that the Central Body radius equals the reference radius.')