From 36bc720f51f261cc93eb834b740a2b3a351cac11 Mon Sep 17 00:00:00 2001 From: anand43 Date: Fri, 10 Nov 2023 16:21:22 -0500 Subject: [PATCH] Changed to 4pi normalisation. Scaled C_lm in python by 4pi for correct C00 terms.UNSURE, STILL TESTING --- src/swiftest/swiftest_sph.f90 | 6 ++++-- swiftest/sph_harmonics.py | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/swiftest/swiftest_sph.f90 b/src/swiftest/swiftest_sph.f90 index cda7d617f..aba2148a7 100644 --- 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 index 79f0abf03..832abc8d3 100644 --- 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.')