Improve NFW and Hernquist _fx_projected and _fx_cumul2d (#1096)

LSSTDESC · Jul 10, 2023 · cad503b · cad503b
1 parent b4f6b70
commit cad503b
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 10 deletions.
diff --git a/pyccl/halos/profiles/hernquist.py b/pyccl/halos/profiles/hernquist.py
@@ -106,13 +106,15 @@ def _fx_projected(self, x):
 
         def f1(xx):
             x2m1 = xx * xx - 1
+            sqx2m1 = np.sqrt(-x2m1)
             return (-3 / 2 / x2m1**2
-                    + (x2m1+3) * np.arccosh(1 / xx) / 2 / np.fabs(x2m1)**2.5)
+                    + (x2m1+3) * np.arcsinh(sqx2m1 / xx) / 2 / (-x2m1)**2.5)
 
         def f2(xx):
             x2m1 = xx * xx - 1
+            sqx2m1 = np.sqrt(x2m1)
             return (-3 / 2 / x2m1**2
-                    + (x2m1+3) * np.arccos(1 / xx) / 2 / np.fabs(x2m1)**2.5)
+                    + (x2m1+3) * np.arcsin(sqx2m1 / xx) / 2 / x2m1**2.5)
 
         xf = x.flatten()
         return np.piecewise(xf,
@@ -143,13 +145,15 @@ def _fx_cumul2d(self, x):
 
         def f1(xx):
             x2m1 = xx * xx - 1
+            sqx2m1 = np.sqrt(-x2m1)
             return (1 + 1 / x2m1
-                    + (x2m1 + 1) * np.arccosh(1 / xx) / np.fabs(x2m1)**1.5)
+                    + (x2m1 + 1) * np.arcsinh(sqx2m1 / xx) / (-x2m1)**1.5)
 
         def f2(xx):
             x2m1 = xx * xx - 1
+            sqx2m1 = np.sqrt(x2m1)
             return (1 + 1 / x2m1
-                    - (x2m1 + 1) * np.arccos(1 / xx) / np.fabs(x2m1)**1.5)
+                    - (x2m1 + 1) * np.arcsin(sqx2m1 / xx) / x2m1**1.5)
 
         xf = x.flatten()
         f = np.piecewise(xf,

diff --git a/pyccl/halos/profiles/nfw.py b/pyccl/halos/profiles/nfw.py
@@ -109,11 +109,13 @@ def _fx_projected(self, x):
 
         def f1(xx):
             x2m1 = xx * xx - 1
-            return 1 / x2m1 + np.arccosh(1 / xx) / np.fabs(x2m1)**1.5
+            sqx2m1 = np.sqrt(-x2m1)
+            return 1 / x2m1 + np.arcsinh(sqx2m1 / xx) / (-x2m1)**1.5
 
         def f2(xx):
             x2m1 = xx * xx - 1
-            return 1 / x2m1 - np.arccos(1 / xx) / np.fabs(x2m1)**1.5
+            sqx2m1 = np.sqrt(x2m1)
+            return 1 / x2m1 - np.arcsin(sqx2m1 / xx) / x2m1**1.5
 
         xf = x.flatten()
         return np.piecewise(xf,
@@ -143,12 +145,12 @@ def _projected_analytic(self, cosmo, r, M, a):
     def _fx_cumul2d(self, x):
 
         def f1(xx):
-            sqx2m1 = np.sqrt(np.fabs(xx * xx - 1))
-            return np.log(0.5 * xx) + np.arccosh(1 / xx) / sqx2m1
+            sqx2m1 = np.sqrt(-(xx * xx - 1))
+            return np.log(0.5 * xx) + np.arcsinh(sqx2m1 / xx) / sqx2m1
 
         def f2(xx):
-            sqx2m1 = np.sqrt(np.fabs(xx * xx - 1))
-            return np.log(0.5 * xx) + np.arccos(1 / xx) / sqx2m1
+            sqx2m1 = np.sqrt(xx * xx - 1)
+            return np.log(0.5 * xx) + np.arcsin(sqx2m1 / xx) / sqx2m1
 
         xf = x.flatten()
         f = np.piecewise(xf,