Added temporary trial code to reduce error

MintonGroup · Nov 2, 2023 · f239c2e · f239c2e
1 parent e32bfa6
commit f239c2e
Showing 1 changed file with 25 additions and 0 deletions.
diff --git a/src/swiftest/swiftest_sph.f90 b/src/swiftest/swiftest_sph.f90
@@ -42,13 +42,27 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
         real(DP)       :: plm, dplm         !! Associated Legendre polynomials at a given l, m
         real(DP)       :: ccss, cssc        !! See definition in source code
         real(DP), dimension(:), allocatable  :: p, p_deriv   !! Associated Lengendre Polynomials at a given cos(theta)
+        real(DP)     :: r2, irh, rinv2, t0, t1, t2, t3, fac0, fac1, fac2, j2rp2, j4rp4, r_fac
 
         g_sph(:) = 0.0_DP
         r_mag = sqrt(dot_product(rh(:), rh(:)))
         l_max = size(c_lm, 2) - 1
         N = (l_max + 1) * (l_max + 2) / 2
         allocate(p(N),p_deriv(N))
 
+        ! to compare w s_obl.f90
+        j2rp2 = -c_lm(1, 3, 1) * r_0**2
+        j4rp4 = -c_lm(1, 4, 1) * r_0**4
+        r2 = dot_product(rh(:), rh(:))
+        irh = 1.0_DP / sqrt(r2)
+        rinv2 = irh**2
+        t0 = -GMcb * rinv2 * rinv2 * irh
+        t1 = 1.5_DP * j2rp2
+        t2 = rh(3) * rh(3) * rinv2
+        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)
+
         if(cos(theta) > epsilon(0.0_DP)) then
             call PLegendreA_d1(p, p_deriv, l_max, cos(theta))      ! Unnormalized Associated Legendre Polynomials and the 1st Derivative
         else
@@ -79,6 +93,17 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
                                                                         * (dplm * cos(theta) + plm * (l + 1))) ! g_y
                 g_sph(3) = g_sph(3) - GMcb * r_0**l / r_mag**(l + 2) * ccss * (-1.0_DP * dplm * sin(theta)**2  &
                                                                         + plm * (l + 1) * cos(theta))          ! g_z
+
+                ! !! Condensed form ???? (to reduce floating point error)
+                ! fac0 = -GMcb / r_mag**2                                                        
+                ! fac1 = cssc * m * plm / sin(theta)
+                ! fac2 = ccss * sin(theta) * (dplm * cos(theta) + plm * (l + 1))
+                ! r_fac = (r_0 / r_mag)**l
+
+                ! g_sph(1) = g_sph(1) + fac0 * r_fac * (fac1 * sin(phi) + fac2 * cos(phi))
+                ! g_sph(2) = g_sph(2) + fac0 * r_fac * (fac1 * cos(phi) + fac2 * sin(phi))
+                ! g_sph(3) = g_sph(3) + fac0 * r_fac * ccss * (-1.0_DP * dplm * sin(theta)**2  &
+                !                                         + plm * (l + 1) * cos(theta))          ! g_z
 
             end do
         end do