From f239c2e05ed81c5a87792ea3cf378c30c7b1cbcf Mon Sep 17 00:00:00 2001
From: anand43 <anand43@purdue.edu>
Date: Thu, 2 Nov 2023 13:54:48 -0400
Subject: [PATCH] Added temporary trial code to reduce error

---
 src/swiftest/swiftest_sph.f90 | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/src/swiftest/swiftest_sph.f90 b/src/swiftest/swiftest_sph.f90
index eb1c2dbc1..5a57e82f5 100644
--- a/src/swiftest/swiftest_sph.f90
+++ b/src/swiftest/swiftest_sph.f90
@@ -42,6 +42,7 @@ 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(:)))
@@ -49,6 +50,19 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
         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