Checked for 1/0 term when sin_theta = 0

MintonGroup · Nov 16, 2023 · d8e23dc · d8e23dc
1 parent 9675e14
commit d8e23dc
Showing 1 changed file with 8 additions and 3 deletions.
diff --git a/src/swiftest/swiftest_sph.f90 b/src/swiftest/swiftest_sph.f90
@@ -138,11 +138,16 @@ module subroutine swiftest_sph_g_acc_one(GMcb, r_0, phi, theta, rh, c_lm, g_sph,
                 ! Condensed form
 
                 ! fac0 = -(m * cos_tmp * plm / sin_tmp - plm1) / sin_tmp ! dplm
-                fac1 = m * plm / sin_theta
-                fac2 = plm * (l + m + 1) * sin_theta + plm1 * cos_theta
-                fac3 = fac2 - fac1
                 ! fac3 = plm * (l + m + 1) * cos_theta
                 ! fac4 = plm1 * sin_theta
+                if(sin_theta .eq. 0) then
+                    fac1 = 0.0_DP
+                else
+                    fac1 = m * plm / sin_theta
+                end if
+
+                fac2 = plm * (l + m + 1) * sin_theta + plm1 * cos_theta
+                fac3 = fac2 - fac1
                 r_fac = -GMcb * r_0**l / r_mag**(l + 2)
 
                 ! g_sph(:) = 0.0_DP