update fitting nnotebook to account for corrected monopole term

LSSTDESC · Aug 16, 2023 · 0ee57e9 · 0ee57e9
1 parent e517e13
commit 0ee57e9
Show file tree

Hide file tree

Showing 2 changed files with 198 additions and 149 deletions.
diff --git a/examples/triaxiality/fitting_elliptical_halo.ipynb b/examples/triaxiality/fitting_elliptical_halo.ipynb
diff --git a/examples/triaxiality/model_functions.py b/examples/triaxiality/model_functions.py
@@ -4,7 +4,12 @@
 from scipy.interpolate import InterpolatedUnivariateSpline
 
 
+
+## CALCULATE MONOPOLE COMPONENT OF GAMMA_T
 def gamma_tangential_monopole(r, mdelta, cdelta, z_cl, cosmo, hpmd='nfw'):
+    ## Calculate the monopole component of gamma_T given by
+    ## 2 * r^-2 * Integrate[ r'*Sigma(r', mdelta, cdelta, z_cl, **kwargs), 0, r] - Sigma(r, mdelta, cdelta, z_cl, **kwargs)
+
     kappa_0 = clmm.compute_surface_density(r, mdelta, cdelta, z_cl, cosmo, delta_mdef=200, 
                                      halo_profile_model=hpmd, massdef='mean', alpha_ein=None, 
                                      verbose=False, validate_input=True)
@@ -21,6 +26,32 @@ def gamma_tangential_monopole(r, mdelta, cdelta, z_cl, cosmo, hpmd='nfw'):
 
 
 
+## CALCULATE MONOPOLE COMPONENT OF GAMMA_T WITH EPSILON^2 CORRECTIONS
+def gamma_tangential_monopole_e2corrected(r, ell, mdelta, cdelta, z_cl, cosmo, hpmd='nfw'):
+    ## Calculate the monopole component of gamma_T given by
+    ## 2 * r^-2 * Integrate[ r'*Sigma0(r', mdelta, cdelta, z_cl, **kwargs), 0, r] - Sigma0(r, mdelta, cdelta, z_cl, **kwargs)
+    ## where Sigma0 = Sigma*(1 + e^2/8 * (eta_0 + eta_0^2/2 + r*(deta_0/dr)/2))
+
+    kappa = clmm.compute_surface_density(r, mdelta, cdelta, z_cl, cosmo, delta_mdef=200,
+                                     halo_profile_model=hpmd, massdef='mean', alpha_ein=None, 
+                                     verbose=False, validate_input=True)
+
+    eta_0  = r * np.gradient(np.log(kappa), r)
+    Deta_0 = r * np.gradient(eta_0, r)
+
+    kappa_0 = kappa * (1 + ell**2/8 * (eta_0 + eta_0**2/2 + Deta_0/2))
+
+    f = lambda r_i, kapp_0_i: r_i * kappa_0_i
+
+    integrate = np.vectorize(pyfunc = scipy.integrate.quad_vec)
+
+    integral, err = integrate(f,0,r)
+
+
+    return (2/r**2)*integral - kappa_0
+
+
+
 def g_tangential_quadrupole(r, mdelta, cdelta, z_cl, cosmo, hpmd='nfw'):
     kappa_0 = clmm.compute_surface_density(r, mdelta, cdelta, z_cl, cosmo, delta_mdef=200, 
                                      halo_profile_model=hpmd, massdef='mean', alpha_ein=None, 
@@ -93,6 +124,10 @@ def _delta_sigma_const(ell_, r, mdelta, cdelta, z_cl, cosmo, hpmd='nfw', sample_
     return (ell_/4.0)*(2*I_2 - sigma_0*eta_0), eta_0_interpolation_func
 
 
+def _sigma(r, mdelta, cdelta, z_cl, cosmo, hpmd='nfw', sample_N=1000, delta_mdef=200):
+    return clmm.compute_surface_density(r, mdelta, cdelta, z_cl, cosmo, delta_mdef=delta_mdef, 
+                                     halo_profile_model=hpmd, massdef='mean', alpha_ein=None, 
+                                     verbose=False, validate_input=True)
 
 def _delta_sigma_excess(ell_, r, mdelta, cdelta, z_cl, cosmo, hpmd='nfw', sample_N=1000, delta_mdef=200):
     return clmm.compute_excess_surface_density(r, mdelta, cdelta, z_cl, cosmo, delta_mdef=delta_mdef,