Implemented Integrated Sachs-Wolfe kernel.

Implemented the derivative of the growth function with respect to redshift. Code NOT tested yet!!!!
NumCosmo · Feb 11, 2022 · d09cf72 · d09cf72
1 parent 408d93b
commit d09cf72
Show file tree

Hide file tree

Showing 11 changed files with 784 additions and 30 deletions.
diff --git a/aminclude_static.am b/aminclude_static.am
@@ -1,6 +1,6 @@
 
 # aminclude_static.am generated automatically by Autoconf
-# from AX_AM_MACROS_STATIC on Fri Sep 10 17:53:39 -03 2021
+# from AX_AM_MACROS_STATIC on Tue Sep 28 14:59:39 -03 2021
 
 
 # Code coverage
@@ -60,7 +60,7 @@ CODE_COVERAGE_GENHTML_OPTIONS_DEFAULT ?=$(if $(CODE_COVERAGE_BRANCH_COVERAGE),--
 CODE_COVERAGE_GENHTML_OPTIONS ?= $(CODE_COVERAGE_GENHTML_OPTIONS_DEFAULT)
 CODE_COVERAGE_IGNORE_PATTERN ?=
 
-GITIGNOREFILES := $(GITIGNOREFILES) $(CODE_COVERAGE_OUTPUT_FILE) $(CODE_COVERAGE_OUTPUT_DIRECTORY)
+GITIGNOREFILES = $(GITIGNOREFILES) $(CODE_COVERAGE_OUTPUT_FILE) $(CODE_COVERAGE_OUTPUT_DIRECTORY)
 code_coverage_v_lcov_cap = $(code_coverage_v_lcov_cap_$(V))
 code_coverage_v_lcov_cap_ = $(code_coverage_v_lcov_cap_$(AM_DEFAULT_VERBOSITY))
 code_coverage_v_lcov_cap_0 = @echo "  LCOV   --capture" $(CODE_COVERAGE_OUTPUT_FILE);
@@ -97,7 +97,7 @@ code-coverage-clean:
 
 code-coverage-dist-clean:
 
-AM_DISTCHECK_CONFIGURE_FLAGS := $(AM_DISTCHECK_CONFIGURE_FLAGS) --disable-code-coverage
+AM_DISTCHECK_CONFIGURE_FLAGS = $(AM_DISTCHECK_CONFIGURE_FLAGS) --disable-code-coverage
  else # ifneq ($(abs_builddir), $(abs_top_builddir))
 check-code-coverage:
 

diff --git a/notebooks/NC_CCL_Colossus_density_profiles.ipynb b/notebooks/NC_CCL_Colossus_density_profiles.ipynb
diff --git a/notebooks/tutorials/DarkEnergyIntro.ipynb b/notebooks/tutorials/DarkEnergyIntro.ipynb
@@ -887,7 +887,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -901,7 +901,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.3"
+   "version": "3.9.7"
   }
  },
  "nbformat": 4,

diff --git a/numcosmo/math/ncm_powspec.c b/numcosmo/math/ncm_powspec.c
@@ -168,6 +168,14 @@ _ncm_powspec_eval (NcmPowspec *powspec, NcmModel *model, const gdouble z, const
   return 0.0;
 }
 
+static gdouble
+_ncm_powspec_deriv_z (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k)
+{
+  g_error ("_ncm_powspec_deriv_z: no default implementation, all children must implement it.");
+
+  return 0.0;
+}
+
 static void _ncm_powspec_eval_vec (NcmPowspec *powspec, NcmModel *model, const gdouble z, NcmVector *k, NcmVector *Pk);
 
 static void
@@ -236,6 +244,7 @@ ncm_powspec_class_init (NcmPowspecClass *klass)
   klass->prepare  = &_ncm_powspec_prepare;
   klass->eval     = &_ncm_powspec_eval;
   klass->eval_vec = &_ncm_powspec_eval_vec;
+  klass->deriv_z  = &_ncm_powspec_deriv_z;
 }
 
 static void

diff --git a/numcosmo/math/ncm_powspec.h b/numcosmo/math/ncm_powspec.h
@@ -51,6 +51,7 @@ struct _NcmPowspecClass
   void (*prepare) (NcmPowspec *powspec, NcmModel *model);
   gdouble (*eval) (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k);
   void (*eval_vec) (NcmPowspec *powspec, NcmModel *model, const gdouble z, NcmVector *k, NcmVector *Pk);
+  gdouble (*deriv_z) (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k);
   void (*get_nknots) (NcmPowspec *powspec, guint *Nz, guint *Nk);
 };
 
@@ -140,6 +141,12 @@ ncm_powspec_eval_vec (NcmPowspec *powspec, NcmModel *model, const gdouble z, Ncm
   return NCM_POWSPEC_GET_CLASS (powspec)->eval_vec (powspec, model, z, k, Pk);
 }
 
+NCM_INLINE gdouble 
+ncm_powspec_deriv_z (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k)
+{
+  return NCM_POWSPEC_GET_CLASS (powspec)->deriv_z (powspec, model, z, k);
+}
+
 G_END_DECLS
 
 #endif /* __GTK_DOC_IGNORE__ */

diff --git a/numcosmo/nc_powspec_ml_cbe.c b/numcosmo/nc_powspec_ml_cbe.c
@@ -167,6 +167,7 @@ _nc_powspec_ml_cbe_finalize (GObject *object)
 
 static void _nc_powspec_ml_cbe_prepare (NcmPowspec *powspec, NcmModel *model);
 static gdouble _nc_powspec_ml_cbe_eval (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k);
+static gdouble _nc_powspec_ml_cbe_deriv_z (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k);
 static void _nc_powspec_ml_cbe_get_nknots (NcmPowspec *powspec, guint *Nz, guint *Nk);
 
 static void
@@ -226,6 +227,7 @@ nc_powspec_ml_cbe_class_init (NcPowspecMLCBEClass *klass)
 
   powspec_class->prepare    = &_nc_powspec_ml_cbe_prepare;
   powspec_class->eval       = &_nc_powspec_ml_cbe_eval;
+  powspec_class->deriv_z    = &_nc_powspec_ml_cbe_deriv_z;
   powspec_class->get_nknots = &_nc_powspec_ml_cbe_get_nknots;
 }
 
@@ -278,6 +280,32 @@ _nc_powspec_ml_cbe_eval (NcmPowspec *powspec, NcmModel *model, const gdouble z,
   }
 }
 
+static gdouble
+_nc_powspec_ml_cbe_deriv_z (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k)
+{
+  NcPowspecMLCBE *ps_cbe             = NC_POWSPEC_ML_CBE (powspec);
+  NcPowspecMLCBEPrivate * const self = ps_cbe->priv;
+
+  if (k < self->intern_k_min)
+  {
+    const gdouble lnkmin = log (self->intern_k_min);
+    const gdouble match  = exp (ncm_spline2d_eval (self->lnPk, lnkmin, z)) / ncm_powspec_eval (NCM_POWSPEC (self->eh), model, z, self->intern_k_min);
+
+    return match * ncm_powspec_eval (NCM_POWSPEC (self->eh), model, z, k) * ncm_spline2d_deriv_dzdy (self->lnPk, lnkmin, z);
+  }
+  else if (k > self->intern_k_max)
+  {
+    const gdouble lnkmax = log (self->intern_k_max);
+    const gdouble match  = exp (ncm_spline2d_eval (self->lnPk, lnkmax, z)) / ncm_powspec_eval (NCM_POWSPEC (self->eh), model, z, self->intern_k_max);
+
+    return match * ncm_powspec_eval (NCM_POWSPEC (self->eh), model, z, k) * ncm_spline2d_deriv_dzdy (self->lnPk, lnkmax, z);
+  }
+  else
+  {
+    return exp (ncm_spline2d_eval (self->lnPk, log (k), z)) * ncm_spline2d_deriv_dzdy (self->lnPk, log (k), z);
+  }
+}
+
 static void
 _nc_powspec_ml_cbe_get_nknots (NcmPowspec *powspec, guint *Nz, guint *Nk)
 {

diff --git a/numcosmo/nc_powspec_ml_fix_spline.c b/numcosmo/nc_powspec_ml_fix_spline.c
@@ -131,6 +131,7 @@ static void _nc_powspec_ml_fix_spline_prepare (NcmPowspec *powspec, NcmModel *mo
 static gdouble _nc_powspec_ml_fix_spline_eval (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k);
 static void _nc_powspec_ml_fix_spline_eval_vec (NcmPowspec *powspec, NcmModel *model, const gdouble z, NcmVector *k, NcmVector *Pk);
 static void _nc_powspec_ml_fix_spline_get_nknots (NcmPowspec *powspec, guint *Nz, guint *Nk);
+static gdouble _nc_powspec_ml_fix_spline_eval_deriv_z (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k);
 
 static void
 nc_powspec_ml_fix_spline_class_init (NcPowspecMLFixSplineClass *klass)
@@ -164,6 +165,7 @@ nc_powspec_ml_fix_spline_class_init (NcPowspecMLFixSplineClass *klass)
   powspec_class->prepare    = &_nc_powspec_ml_fix_spline_prepare;
   powspec_class->eval       = &_nc_powspec_ml_fix_spline_eval;
   powspec_class->eval_vec   = &_nc_powspec_ml_fix_spline_eval_vec;
+  powspec_class->deriv_z    = &_nc_powspec_ml_fix_spline_eval_deriv_z;
   powspec_class->get_nknots = &_nc_powspec_ml_fix_spline_get_nknots;
 }
 
@@ -209,6 +211,20 @@ _nc_powspec_ml_fix_spline_eval_vec (NcmPowspec *powspec, NcmModel *model, const
   }
 }
 
+static gdouble
+_nc_powspec_ml_fix_spline_eval_deriv_z (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k)
+{
+  NcHICosmo *cosmo            = NC_HICOSMO (model);
+  NcHIPrim *prim              = NC_HIPRIM (ncm_model_peek_submodel_by_mid (model, nc_hiprim_id ()));
+  NcPowspecMLFixSpline *ps_fs = NC_POWSPEC_ML_FIX_SPLINE (powspec);
+  const gdouble growth        = nc_growth_func_eval (ps_fs->gf, cosmo, z);
+  const gdouble deriv_growth  = nc_growth_func_eval_deriv (ps_fs->gf, cosmo, z) * (-1.0 / gsl_pow_2(1.0 + z)); /* dD(z)/dz */
+
+  NCM_UNUSED (prim);
+
+  return ncm_spline_eval (ps_fs->Pk, k) * 2.0 * growth * deriv_growth;
+}
+
 static void
 _nc_powspec_ml_fix_spline_get_nknots (NcmPowspec *powspec, guint *Nz, guint *Nk)
 {

diff --git a/numcosmo/nc_powspec_ml_transfer.c b/numcosmo/nc_powspec_ml_transfer.c
@@ -142,6 +142,7 @@ _nc_powspec_ml_transfer_finalize (GObject *object)
 static void _nc_powspec_ml_transfer_prepare (NcmPowspec *powspec, NcmModel *model);
 static gdouble _nc_powspec_ml_transfer_eval (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k);
 static void _nc_powspec_ml_transfer_eval_vec (NcmPowspec *powspec, NcmModel *model, const gdouble z, NcmVector *k, NcmVector *Pk);
+static gdouble _nc_powspec_ml_transfer_deriv_z (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k);
 static void _nc_powspec_ml_transfer_get_nknots (NcmPowspec *powspec, guint *Nz, guint *Nk);
 
 static void
@@ -189,6 +190,7 @@ nc_powspec_ml_transfer_class_init (NcPowspecMLTransferClass *klass)
   powspec_class->prepare    = &_nc_powspec_ml_transfer_prepare;
   powspec_class->eval       = &_nc_powspec_ml_transfer_eval;
   powspec_class->eval_vec   = &_nc_powspec_ml_transfer_eval_vec;
+  powspec_class->deriv_z    = &_nc_powspec_ml_transfer_deriv_z;
   powspec_class->get_nknots = &_nc_powspec_ml_transfer_get_nknots;
 }
 
@@ -252,6 +254,22 @@ _nc_powspec_ml_transfer_eval_vec (NcmPowspec *powspec, NcmModel *model, const gd
   }
 }
 
+static gdouble
+_nc_powspec_ml_transfer_deriv_z (NcmPowspec *powspec, NcmModel *model, const gdouble z, const gdouble k)
+{
+  NcHICosmo *cosmo            = NC_HICOSMO (model);
+  NcHIPrim *prim              = NC_HIPRIM (ncm_model_peek_submodel_by_mid (model, nc_hiprim_id ()));
+  NcPowspecMLTransfer *ps_mlt = NC_POWSPEC_ML_TRANSFER (powspec);
+  const gdouble kh            = k / nc_hicosmo_h (cosmo);
+  const gdouble growth        = nc_growth_func_eval (ps_mlt->gf, cosmo, z);
+  const gdouble deriv_growth  = nc_growth_func_eval_deriv (ps_mlt->gf, cosmo, z) * (-1.0 / gsl_pow_2(1.0 + z));
+  const gdouble tf            = nc_transfer_func_eval (ps_mlt->tf, cosmo, kh);
+  const gdouble tf2           = tf * tf;
+  const gdouble Delta_zeta_k  = nc_hiprim_SA_powspec_k (prim, k);
+
+  return k * Delta_zeta_k * ps_mlt->Pm_k2Pzeta * tf2 * 2.0 * growth * deriv_growth;
+}
+
 static void
 _nc_powspec_ml_transfer_get_nknots (NcmPowspec *powspec, guint *Nz, guint *Nk)
 {

diff --git a/numcosmo/xcor/nc_xcor.c b/numcosmo/xcor/nc_xcor.c
@@ -30,7 +30,7 @@
  *
  * The angular power spectrum between observables $A$ and $B$ with kernels $W^A$ and $W^B$ is
  * \begin{equation}
- * C_{\ell}^{AB} &= \int dz W^A(z) \int dz^\prime W^B (z^\prime) \times \int dk \frac{2}{\pi} k^2 P(k, z, z^\prime) j_{\ell}(k\chi(z)) j_{\ell} (k\chi(z^\prime)).
+ * C_{\ell}^{AB} = \int dz W^A(z) \int dz^\prime W^B (z^\prime) \times \int dk \frac{2}{\pi} k^2 P(k, z, z^\prime) j_{\ell}(k\chi(z)) j_{\ell} (k\chi(z^\prime)).
  * \end{equation}
  * In the Limber approximation, it reduces to
  * \begin{equation}