6817479045

Committed 09 Nov 2023 09:11PM UTC coverage: 53.274% (-0.2%) from 53.454%

Build # 6817479045

Build Type

Pull #107

github

Committed by

web-flow

Commit Message

Merge 034f37fb5 into 7269172d1

Pull Request Pull Request #107: Sample variance

Run Details

8770 of 20813 branches covered (0.0%)

Branch coverage included in aggregate %.

55 of 534 new or added lines in 6 files covered. (10.3%)

49 existing lines in 7 files now uncovered.

44633 of 79430 relevant lines covered (56.19%)

3038914.66 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

58.62

/numcosmo/math/ncm_fftlog_sbessel_j.c

/***************************************************************************
 *            ncm_fftlog_sbessel_j.c
 *
 *  Wed July 19 10:00:26 2017
 *  Copyright  2017  Fernando de Simoni
 *  <fernando.saliby@gmail.com>
 ****************************************************************************/

/***************************************************************************
 *            ncm_fftlog_sbessel_j.c
 *
 *  Sat September 02 18:11:00 2017
 *  Copyright  2017  Sandro Dias Pinto Vitenti
 *  <vitenti@uel.br>
 ****************************************************************************/

/*
 * ncm_fftlog_sbessel_j.c
 *
 * Copyright (C) 2017 - Fernando de Simoni
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 */

/**
 * SECTION:ncm_fftlog_sbessel_j
 * @title: NcmFftlogSBesselJ
 * @short_description: Logarithm fast fourier transform for a kernel given by the spatial correlation function multipoles.
 * @stability: Stable
 * @include: numcosmo/math/ncm_fftlog_sbessel_j.h
 *
 *
 * This object computes the function (see #NcmFftlog)
 * $$Y_n = \int_0^\infty t^{\frac{2\pi i n}{L}} K(t) dt,$$
 * where the kernel are the spherical Bessel function
 * of the first kind multiplied by a power law,
 *
 * \begin{equation}\label{eq:kerneljl}
 * K(t) = t^q j_{\ell}(t).
 * \end{equation}
 *
 * Note that the spherical Bessel function's order, $\ell$ (#NcmFftlogSBesselJ:ell), must be an integer number.
 *
 * The spatial correlation function multipoles, $\xi_{\ell}^{(n)}(r)$, can be defined as
 * (see [Matsubara (2004)][XMatsubara2004a] [[arXiv](https://arxiv.org/abs/astro-ph/0408349)])
 *
 * \begin{equation}\label{eq:xi_multipoles}
 * \xi_{\ell}^{(n)}(r) = \frac{(-1)^{n+\ell}}{r^{2n-\ell}} \int_{0}^{\infty} \frac{\mathrm{d} k}{2\pi^2} \frac{k^2}{k^{2n-\ell}} j_{\ell}(kr) P(k) \,\, .
 * \end{equation}
 * Where, $P(k)$ is the power spectrum (see #NcmPowspec).
 *
 * The multipoles integral can be written in the following format
 *
 * \begin{equation*}
 * \xi_{\ell}^{(n)}(r) = \frac{(-1)^{n+\ell}}{r^2} \int_{0}^{\infty} \mathrm{d}k \, (kr)^{2-2n+\ell} \, j_{\ell}(kr) P(k) \,\, .
 * \end{equation*}
 *
 * The object #NcmFftlogSBesselJ can be used to evaluate the above integral in several ways.
 * For example, the integral can be evaluated by defining the function (see #NcmFftlog for more information)
 * \begin{equation*}
 * F(k) = k^{2-2n+\ell} \, P(k)
 * \end{equation*}
 * and the kernel
 * \begin{equation*}
 * K(t) = j_{\ell}(t) \,\, .
 * \end{equation*}
 * Where, $t=kr$ and $r^{(2-2n+\ell)}$ was taken out of the integral.
 * Comparing this kernel with the one defined in Eq. \eqref{eq:kerneljl}, we have $q=0$.
 *
 * But instead, one might choose another format for the function,
 * \begin{equation*}
 * F(k) = k^{\ell} \, P(k)
 * \end{equation*}
 * and the kernel
 * \begin{equation*}
 * K(t) = t^{2-2n} \, j_{\ell}(t) \,\, ,
 * \end{equation*}
 * which evaluates the same integral, but now with $q=2-2n$, and
 * in this case, the term $r^{\ell}$ was the one taken out of the integral.
 * Therefore, the parameter $q$ is the power of the wavenumber $k$ times the distance $r$, $t=kr$,
 * included to the kernel with the spherical Bessel function.
 * Hereafter, it will be referred to as "spherical Bessel power" (#NcmFftlogSBesselJ:q).
 *
 * In general, $q=0$ is an accurate and fast choice to make, but it is interesting to
 * perform tests to evaluate which kernel format fits best for each type of integral.
 *
 * The #NcmPowspecCorr3d object already evaluates Eq. \eqref{eq:xi_multipoles}
 * for the case of the monopole, $n=\ell=0$, with support for redshift evolution.
 *
 */

#ifdef HAVE_CONFIG_H
#  include "config.h"
#endif /* HAVE_CONFIG_H */
#include "build_cfg.h"

#include "math/ncm_fftlog_sbessel_j.h"
#include "math/ncm_cfg.h"
#include "math/ncm_c.h"

#ifndef NUMCOSMO_GIR_SCAN
#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_math.h>
#include <complex.h>
#ifdef NUMCOSMO_HAVE_FFTW3
#include <fftw3.h>
#endif /* NUMCOSMO_HAVE_FFTW3 */
#include <math.h>
#endif /* NUMCOSMO_GIR_SCAN */

typedef struct _NcmFftlogSBesselJPrivate
{
  guint ell;
  gdouble q;
} NcmFftlogSBesselJPrivate;

enum
{
  PROP_0,
  PROP_ELL,
  PROP_Q,
  PROP_SIZE,
};

struct _NcmFftlogSBesselJ {
  NcmFftlog parent_instance;
};

G_DEFINE_TYPE_WITH_PRIVATE (NcmFftlogSBesselJ, ncm_fftlog_sbessel_j, NCM_TYPE_FFTLOG);

static void
ncm_fftlog_sbessel_j_init (NcmFftlogSBesselJ *fftlog_jl)
{
  NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);
  
  self->ell = 0;
  self->q   = 0.0;
}

static void
_ncm_fftlog_sbessel_j_set_property (GObject *object, guint prop_id, const GValue *value, GParamSpec *pspec)
{
  NcmFftlogSBesselJ *fftlog_jl = NCM_FFTLOG_SBESSEL_J (object);
  
  g_return_if_fail (NCM_IS_FFTLOG_SBESSEL_J (object));
  
  switch (prop_id)
  {
    case PROP_ELL:
      ncm_fftlog_sbessel_j_set_ell (fftlog_jl, g_value_get_uint (value));
      break;
    case PROP_Q:
      ncm_fftlog_sbessel_j_set_q (fftlog_jl, g_value_get_double (value));
      break;
    default:
      G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
      break;
  }
}

static void
_ncm_fftlog_sbessel_j_get_property (GObject *object, guint prop_id, GValue *value, GParamSpec *pspec)
{
  NcmFftlogSBesselJ *fftlog_jl = NCM_FFTLOG_SBESSEL_J (object);
  
  g_return_if_fail (NCM_IS_FFTLOG_SBESSEL_J (object));
  
  switch (prop_id)
  {
    case PROP_ELL:
      g_value_set_uint (value, ncm_fftlog_sbessel_j_get_ell (fftlog_jl));
      break;
    case PROP_Q:
      g_value_set_double (value, ncm_fftlog_sbessel_j_get_q (fftlog_jl));
      break;
    default:
      G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
      break;
  }
}

static void
_ncm_fftlog_sbessel_j_finalize (GObject *object)
{
  /* Chain up : end */
  G_OBJECT_CLASS (ncm_fftlog_sbessel_j_parent_class)->finalize (object);
}

static void _ncm_fftlog_sbessel_j_get_Ym (NcmFftlog *fftlog, gpointer Ym_0);

static void
ncm_fftlog_sbessel_j_class_init (NcmFftlogSBesselJClass *klass)
{
  GObjectClass *object_class   = G_OBJECT_CLASS (klass);
  NcmFftlogClass *fftlog_class = NCM_FFTLOG_CLASS (klass);
  
  object_class->set_property = &_ncm_fftlog_sbessel_j_set_property;
  object_class->get_property = &_ncm_fftlog_sbessel_j_get_property;
  object_class->finalize     = &_ncm_fftlog_sbessel_j_finalize;
  
  /**
   * NcmFftlogSBesselJ:ell:
   *
   * The spherical Bessel integer order.
   *
   */
  g_object_class_install_property (object_class,
                                   PROP_ELL,
                                   g_param_spec_uint ("ell",
                                                      NULL,
                                                      "Spherical Bessel integer order",
                                                      0, G_MAXUINT32, 0,
                                                      G_PARAM_READWRITE | G_PARAM_CONSTRUCT | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));
  
  /**
   * NcmFftlogSBesselJ:q:
   *
   * The spherical Bessel power, i.e., the power of the variable $t=kr$,
   * included to the kernel $K(t)$ multiplying the spherical Bessel function.
   *
   */
  g_object_class_install_property (object_class,
                                   PROP_Q,
                                   g_param_spec_double ("q",
                                                        NULL,
                                                        "Spherical Bessel power",
                                                        -G_MAXDOUBLE, G_MAXDOUBLE, 0.0,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));
  
  fftlog_class->name   = "sbessel_j";
  fftlog_class->get_Ym = &_ncm_fftlog_sbessel_j_get_Ym;
}

static void
_ncm_fftlog_sbessel_j_get_Ym (NcmFftlog *fftlog, gpointer Ym_0)
{
  NcmFftlogSBesselJ *fftlog_jl          = NCM_FFTLOG_SBESSEL_J (fftlog);
  NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);
  
  const gdouble pi_sqrt  = sqrt (M_PI);
  const gdouble twopi_Lt = 2.0 * M_PI / ncm_fftlog_get_full_length (fftlog);
  const gint Nf          = ncm_fftlog_get_full_size (fftlog);
  
#ifdef NUMCOSMO_HAVE_FFTW3
  fftw_complex *Ym_base = (fftw_complex *) Ym_0;
  gint i;
  
  if (self->q == 0.5)
  {
    for (i = 0; i < Nf; i++)
    {
      const gint phys_i             = ncm_fftlog_get_mode_index (fftlog, i);
      const complex double a        = twopi_Lt * phys_i * I;
      const complex double A        = a + 0.5;
      const complex double xup      = 0.5 * (1.0 + 1.0 * self->ell + A);
      const complex double two_x_m1 = cpow (2.0, A - 1.0);
      complex double U;
      
      gsl_sf_result lngamma_rho_up, lngamma_theta_up;
      
      gsl_sf_lngamma_complex_e (creal (xup), cimag (xup), &lngamma_rho_up, &lngamma_theta_up);
      
      U = cexp (2.0 * I * lngamma_theta_up.val);
      
      Ym_base[i] = pi_sqrt * two_x_m1 * U;
    }
  }
  else
  {
    const gdouble q = self->q;
    
    for (i = 0; i < Nf; i++)
    {
      const gint phys_i             = ncm_fftlog_get_mode_index (fftlog, i);
      const complex double a        = twopi_Lt * phys_i * I;
      const complex double A        = a + q;
      const complex double xup      = 0.5 * (1.0 + 1.0 * self->ell + A);
      const complex double xdw      = 0.5 * (2.0 + 1.0 * self->ell - A);
      const complex double two_x_m1 = cpow (2.0, A - 1.0);
      complex double U;
      
      gsl_sf_result lngamma_rho_up, lngamma_theta_up;
      gsl_sf_result lngamma_rho_dw, lngamma_theta_dw;
      
      gsl_sf_lngamma_complex_e (creal (xup), cimag (xup), &lngamma_rho_up, &lngamma_theta_up);
      gsl_sf_lngamma_complex_e (creal (xdw), cimag (xdw), &lngamma_rho_dw, &lngamma_theta_dw);
      
      U = cexp ((lngamma_rho_up.val - lngamma_rho_dw.val) + I * (lngamma_theta_up.val - lngamma_theta_dw.val));
      
      /*printf ("% 22.15g % 22.15g % 22.15g % 22.15g\n", lngamma_rho_up.val, lngamma_rho_dw.val, lngamma_theta_up.val, lngamma_theta_dw.val);*/
      
      Ym_base[i] = pi_sqrt * two_x_m1 * U;
    }
  }
  
#endif /* NUMCOSMO_HAVE_FFTW3 */
}

/**
 * ncm_fftlog_sbessel_j_new:
 * @ell: Spherical Bessel Integer order
 * @lnr0: output center $\ln(r_0)$
 * @lnk0: input center $\ln(k_0)$
 * @Lk: input/output interval size
 * @N: number of knots
 *
 * Creates a new fftlog Spherical Bessel J object.
 *
 * Returns: (transfer full): a new #NcmFftlogSBesselJ
 */
NcmFftlogSBesselJ *
ncm_fftlog_sbessel_j_new (guint ell, gdouble lnr0, gdouble lnk0, gdouble Lk, guint N)
{
  NcmFftlogSBesselJ *fftlog_jl = g_object_new (NCM_TYPE_FFTLOG_SBESSEL_J,
                                               "ell",  ell,
                                               "lnr0", lnr0,
                                               "lnk0", lnk0,
                                               "Lk",   Lk,
                                               "N",    N,
                                               NULL);
  
  return fftlog_jl;
}

/**
 * ncm_fftlog_sbessel_j_set_ell:
 * @fftlog_jl: a #NcmFftlogSBesselJ
 * @ell: Spherical Bessel integer order $\ell$
 *
 * Sets @ell as the Spherical Bessel integer order $\ell$.
 *
 */
void
ncm_fftlog_sbessel_j_set_ell (NcmFftlogSBesselJ *fftlog_jl, const guint ell)
{
  NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);
  
  if (self->ell != ell)
  {
    NcmFftlog *fftlog = NCM_FFTLOG (fftlog_jl);
    
    self->ell = ell;
    ncm_fftlog_reset (fftlog);
  }
}

/**
 * ncm_fftlog_sbessel_j_get_ell:
 * @fftlog_jl: a #NcmFftlogSBesselJ
 *
 * Returns: the current Spherical Bessel integer order $\ell$.
 */
guint
ncm_fftlog_sbessel_j_get_ell (NcmFftlogSBesselJ *fftlog_jl)
{
  NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);
  
  return self->ell;
}

/**
 * ncm_fftlog_sbessel_j_set_q:
 * @fftlog_jl: a #NcmFftlogSBesselJ
 * @q: Spherical Bessel power factor $q$
 *
 * Sets @q as the Spherical Bessel power.
 *
 */
void
ncm_fftlog_sbessel_j_set_q (NcmFftlogSBesselJ *fftlog_jl, const gdouble q)
{
  NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);
  
  if (self->q != q)
  {
    NcmFftlog *fftlog = NCM_FFTLOG (fftlog_jl);
    
    self->q = q;
    ncm_fftlog_reset (fftlog);
  }
}

/**
 * ncm_fftlog_sbessel_j_get_q:
 * @fftlog_jl: a #NcmFftlogSBesselJ
 *
 * Returns: the current Spherical Bessel power $q$.
 */
gdouble
ncm_fftlog_sbessel_j_get_q (NcmFftlogSBesselJ *fftlog_jl)
{
  NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);
  
  return self->q;
}

/**
 * ncm_fftlog_sbessel_j_set_best_lnr0:
 * @fftlog_jl: a #NcmFftlogSBesselJ
 *
 * Sets the value of $\ln(r_0)$ which gives the best results for
 * the transformation based on the current value of $\ln(k_0)$,
 * this is based in the rule of thumb $\mathrm{max}_{x^*}(j_l)$
 * where $ x^* \approx l + 1$.
 *
 */
void
ncm_fftlog_sbessel_j_set_best_lnr0 (NcmFftlogSBesselJ *fftlog_jl)
{
  NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);
  NcmFftlog *fftlog                     = NCM_FFTLOG (fftlog_jl);
  
  gint signp = 0;
  
  const gdouble lnk0      = ncm_fftlog_get_lnk0 (fftlog);
  const gdouble Lk        = ncm_fftlog_get_length (fftlog);
  const gdouble ell       = self->ell;
  const gdouble lnc0      = (ell == 0) ? 0.0 : ((ell - 1.0) * Lk + 2.0 * (ell + 1.0) * M_LN2 - ncm_c_lnpi () + 2.0 * lgamma_r (1.5 + ell, &signp)) / (2.0 * (1.0 + ell));
  const gdouble best_lnr0 = -lnk0 + lnc0;
  
  ncm_fftlog_set_lnr0 (fftlog, best_lnr0);
}

/**
 * ncm_fftlog_sbessel_j_set_best_lnk0:
 * @fftlog_jl: a #NcmFftlogSBesselJ
 *
 * Sets the value of $\ln(k_0)$ which gives the best results for
 * the transformation based on the current value of $\ln(r_0)$,
 * this is based in the rule of thumb $\mathrm{max}_{x^*}(j_l)$
 * where $ x^* \approx l + 1$.
 *
 */
void
ncm_fftlog_sbessel_j_set_best_lnk0 (NcmFftlogSBesselJ *fftlog_jl)
{
  NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);
  NcmFftlog *fftlog                     = NCM_FFTLOG (fftlog_jl);
  
  gint signp = 0;
  
  const gdouble lnr0      = ncm_fftlog_get_lnr0 (fftlog);
  const gdouble Lk        = ncm_fftlog_get_length (fftlog);
  const gdouble ell       = self->ell;
  const gdouble lnc0      = (ell == 0) ? 0.0 : ((ell - 1.0) * Lk + 2.0 * (ell + 1.0) * M_LN2 - ncm_c_lnpi () + 2.0 * lgamma_r (1.5 + ell, &signp)) / (2.0 * (1.0 + ell));
  const gdouble best_lnk0 = -lnr0 + lnc0;
  
  ncm_fftlog_set_lnk0 (fftlog, best_lnk0);
}


1	/***************************************************************************
2	* ncm_fftlog_sbessel_j.c
3	*
4	* Wed July 19 10:00:26 2017
5	* Copyright 2017 Fernando de Simoni
6	* <fernando.saliby@gmail.com>
7	****************************************************************************/
8
9	/***************************************************************************
10	* ncm_fftlog_sbessel_j.c
11	*
12	* Sat September 02 18:11:00 2017
13	* Copyright 2017 Sandro Dias Pinto Vitenti
14	* <vitenti@uel.br>
15	****************************************************************************/
16
17	/*
18	* ncm_fftlog_sbessel_j.c
19	*
20	* Copyright (C) 2017 - Fernando de Simoni
21	*
22	* This program is free software; you can redistribute it and/or modify
23	* it under the terms of the GNU General Public License as published by
24	* the Free Software Foundation; either version 2 of the License, or
25	* (at your option) any later version.
26	*
27	* This program is distributed in the hope that it will be useful,
28	* but WITHOUT ANY WARRANTY; without even the implied warranty of
29	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30	* GNU General Public License for more details.
31	*
32	* You should have received a copy of the GNU General Public License
33	* along with this program. If not, see <http://www.gnu.org/licenses/>.
34	*/
35
36	/**
37	* SECTION:ncm_fftlog_sbessel_j
38	* @title: NcmFftlogSBesselJ
39	* @short_description: Logarithm fast fourier transform for a kernel given by the spatial correlation function multipoles.
40	* @stability: Stable
41	* @include: numcosmo/math/ncm_fftlog_sbessel_j.h
42	*
43	*
44	* This object computes the function (see #NcmFftlog)
45	* $$Y_n = \int_0^\infty t^{\frac{2\pi i n}{L}} K(t) dt,$$
46	* where the kernel are the spherical Bessel function
47	* of the first kind multiplied by a power law,
48	*
49	* \begin{equation}\label{eq:kerneljl}
50	* K(t) = t^q j_{\ell}(t).
51	* \end{equation}
52	*
53	* Note that the spherical Bessel function's order, $\ell$ (#NcmFftlogSBesselJ:ell), must be an integer number.
54	*
55	* The spatial correlation function multipoles, $\xi_{\ell}^{(n)}(r)$, can be defined as
56	* (see [Matsubara (2004)][XMatsubara2004a] [[arXiv](https://arxiv.org/abs/astro-ph/0408349)])
57	*
58	* \begin{equation}\label{eq:xi_multipoles}
59	* \xi_{\ell}^{(n)}(r) = \frac{(-1)^{n+\ell}}{r^{2n-\ell}} \int_{0}^{\infty} \frac{\mathrm{d} k}{2\pi^2} \frac{k^2}{k^{2n-\ell}} j_{\ell}(kr) P(k) \,\, .
60	* \end{equation}
61	* Where, $P(k)$ is the power spectrum (see #NcmPowspec).
62	*
63	* The multipoles integral can be written in the following format
64	*
65	* \begin{equation*}
66	* \xi_{\ell}^{(n)}(r) = \frac{(-1)^{n+\ell}}{r^2} \int_{0}^{\infty} \mathrm{d}k \, (kr)^{2-2n+\ell} \, j_{\ell}(kr) P(k) \,\, .
67	* \end{equation*}
68	*
69	* The object #NcmFftlogSBesselJ can be used to evaluate the above integral in several ways.
70	* For example, the integral can be evaluated by defining the function (see #NcmFftlog for more information)
71	* \begin{equation*}
72	* F(k) = k^{2-2n+\ell} \, P(k)
73	* \end{equation*}
74	* and the kernel
75	* \begin{equation*}
76	* K(t) = j_{\ell}(t) \,\, .
77	* \end{equation*}
78	* Where, $t=kr$ and $r^{(2-2n+\ell)}$ was taken out of the integral.
79	* Comparing this kernel with the one defined in Eq. \eqref{eq:kerneljl}, we have $q=0$.
80	*
81	* But instead, one might choose another format for the function,
82	* \begin{equation*}
83	* F(k) = k^{\ell} \, P(k)
84	* \end{equation*}
85	* and the kernel
86	* \begin{equation*}
87	* K(t) = t^{2-2n} \, j_{\ell}(t) \,\, ,
88	* \end{equation*}
89	* which evaluates the same integral, but now with $q=2-2n$, and
90	* in this case, the term $r^{\ell}$ was the one taken out of the integral.
91	* Therefore, the parameter $q$ is the power of the wavenumber $k$ times the distance $r$, $t=kr$,
92	* included to the kernel with the spherical Bessel function.
93	* Hereafter, it will be referred to as "spherical Bessel power" (#NcmFftlogSBesselJ:q).
94	*
95	* In general, $q=0$ is an accurate and fast choice to make, but it is interesting to
96	* perform tests to evaluate which kernel format fits best for each type of integral.
97	*
98	* The #NcmPowspecCorr3d object already evaluates Eq. \eqref{eq:xi_multipoles}
99	* for the case of the monopole, $n=\ell=0$, with support for redshift evolution.
100	*
101	*/
102
103	#ifdef HAVE_CONFIG_H
104	# include "config.h"
105	#endif /* HAVE_CONFIG_H */
106	#include "build_cfg.h"
107
108	#include "math/ncm_fftlog_sbessel_j.h"
109	#include "math/ncm_cfg.h"
110	#include "math/ncm_c.h"
111
112	#ifndef NUMCOSMO_GIR_SCAN
113	#include <gsl/gsl_sf_result.h>
114	#include <gsl/gsl_sf_gamma.h>
115	#include <gsl/gsl_sf_trig.h>
116	#include <gsl/gsl_math.h>
117	#include <complex.h>
118	#ifdef NUMCOSMO_HAVE_FFTW3
119	#include <fftw3.h>
120	#endif /* NUMCOSMO_HAVE_FFTW3 */
121	#include <math.h>
122	#endif /* NUMCOSMO_GIR_SCAN */
123
124	typedef struct _NcmFftlogSBesselJPrivate
125	{
126	guint ell;
127	gdouble q;
128	} NcmFftlogSBesselJPrivate;
129
130	enum
131	{
132	PROP_0,
133	PROP_ELL,
134	PROP_Q,
135	PROP_SIZE,
136	};
137
138	struct _NcmFftlogSBesselJ {
139	NcmFftlog parent_instance;
140	};
141
142	G_DEFINE_TYPE_WITH_PRIVATE (NcmFftlogSBesselJ, ncm_fftlog_sbessel_j, NCM_TYPE_FFTLOG);	110!
143
144	static void
145	ncm_fftlog_sbessel_j_init (NcmFftlogSBesselJ *fftlog_jl)	8✔
146	{
147	NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);	8✔
148
149	self->ell = 0;	8✔
150	self->q = 0.0;	8✔
151	}	8✔
152
153	static void
154	_ncm_fftlog_sbessel_j_set_property (GObject object, guint prop_id, const GValue value, GParamSpec *pspec)	16✔
155	{
156	NcmFftlogSBesselJ *fftlog_jl = NCM_FFTLOG_SBESSEL_J (object);	16✔
157
158	g_return_if_fail (NCM_IS_FFTLOG_SBESSEL_J (object));	16✔
159
160	switch (prop_id)	16!
161	{
162	case PROP_ELL:	8✔
163	ncm_fftlog_sbessel_j_set_ell (fftlog_jl, g_value_get_uint (value));	8✔
164	break;	8✔
165	case PROP_Q:	8✔
166	ncm_fftlog_sbessel_j_set_q (fftlog_jl, g_value_get_double (value));	8✔
167	break;	8✔
168	default:	×
169	G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
170	break;	×
171	}
172	}
173
174	static void
175	_ncm_fftlog_sbessel_j_get_property (GObject object, guint prop_id, GValue value, GParamSpec *pspec)	×
176	{
177	NcmFftlogSBesselJ *fftlog_jl = NCM_FFTLOG_SBESSEL_J (object);	×
178
179	g_return_if_fail (NCM_IS_FFTLOG_SBESSEL_J (object));	×
180
181	switch (prop_id)	×
182	{
183	case PROP_ELL:	×
184	g_value_set_uint (value, ncm_fftlog_sbessel_j_get_ell (fftlog_jl));	×
185	break;	×
186	case PROP_Q:	×
187	g_value_set_double (value, ncm_fftlog_sbessel_j_get_q (fftlog_jl));	×
188	break;	×
189	default:	×
190	G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
191	break;	×
192	}
193	}
194
195	static void
196	_ncm_fftlog_sbessel_j_finalize (GObject *object)	8✔
197	{
198	/* Chain up : end */
199	G_OBJECT_CLASS (ncm_fftlog_sbessel_j_parent_class)->finalize (object);	8✔
200	}	8✔
201
202	static void _ncm_fftlog_sbessel_j_get_Ym (NcmFftlog *fftlog, gpointer Ym_0);
203
204	static void
205	ncm_fftlog_sbessel_j_class_init (NcmFftlogSBesselJClass *klass)	110✔
206	{
207	GObjectClass *object_class = G_OBJECT_CLASS (klass);	110✔
208	NcmFftlogClass *fftlog_class = NCM_FFTLOG_CLASS (klass);	110✔
209
210	object_class->set_property = &_ncm_fftlog_sbessel_j_set_property;	110✔
211	object_class->get_property = &_ncm_fftlog_sbessel_j_get_property;	110✔
212	object_class->finalize = &_ncm_fftlog_sbessel_j_finalize;	110✔
213
214	/**
215	* NcmFftlogSBesselJ:ell:
216	*
217	* The spherical Bessel integer order.
218	*
219	*/
220	g_object_class_install_property (object_class,	110✔
221	PROP_ELL,
222	g_param_spec_uint ("ell",
223	NULL,
224	"Spherical Bessel integer order",
225	0, G_MAXUINT32, 0,
226	G_PARAM_READWRITE \| G_PARAM_CONSTRUCT \| G_PARAM_STATIC_NAME \| G_PARAM_STATIC_BLURB));
227
228	/**
229	* NcmFftlogSBesselJ:q:
230	*
231	* The spherical Bessel power, i.e., the power of the variable $t=kr$,
232	* included to the kernel $K(t)$ multiplying the spherical Bessel function.
233	*
234	*/
235	g_object_class_install_property (object_class,	110✔
236	PROP_Q,
237	g_param_spec_double ("q",
238	NULL,
239	"Spherical Bessel power",
240	-G_MAXDOUBLE, G_MAXDOUBLE, 0.0,
241	G_PARAM_READWRITE \| G_PARAM_CONSTRUCT \| G_PARAM_STATIC_NAME \| G_PARAM_STATIC_BLURB));
242
243	fftlog_class->name = "sbessel_j";	110✔
244	fftlog_class->get_Ym = &_ncm_fftlog_sbessel_j_get_Ym;	110✔
245	}	110✔
246
247	static void
248	_ncm_fftlog_sbessel_j_get_Ym (NcmFftlog *fftlog, gpointer Ym_0)	58✔
249	{
250	NcmFftlogSBesselJ *fftlog_jl = NCM_FFTLOG_SBESSEL_J (fftlog);	58✔
251	NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);	58✔
252
253	const gdouble pi_sqrt = sqrt (M_PI);	58✔
254	const gdouble twopi_Lt = 2.0 * M_PI / ncm_fftlog_get_full_length (fftlog);	58✔
255	const gint Nf = ncm_fftlog_get_full_size (fftlog);	58✔
256
257	#ifdef NUMCOSMO_HAVE_FFTW3
258	fftw_complex Ym_base = (fftw_complex ) Ym_0;	58✔
259	gint i;
260
261	if (self->q == 0.5)	58!
262	{
263	for (i = 0; i < Nf; i++)	×
264	{
265	const gint phys_i = ncm_fftlog_get_mode_index (fftlog, i);	×
266	const complex double a = twopi_Lt * phys_i * I;	×
267	const complex double A = a + 0.5;	×
268	const complex double xup = 0.5 * (1.0 + 1.0 * self->ell + A);	×
269	const complex double two_x_m1 = cpow (2.0, A - 1.0);	×
270	complex double U;
271
272	gsl_sf_result lngamma_rho_up, lngamma_theta_up;
273
274	gsl_sf_lngamma_complex_e (creal (xup), cimag (xup), &lngamma_rho_up, &lngamma_theta_up);	×
275
276	U = cexp (2.0 * I * lngamma_theta_up.val);	×
277
278	Ym_base[i] = pi_sqrt * two_x_m1 * U;	×
279	}
280	}
281	else
282	{
283	const gdouble q = self->q;	58✔
284
285	for (i = 0; i < Nf; i++)	47,738✔
286	{
287	const gint phys_i = ncm_fftlog_get_mode_index (fftlog, i);	47,680✔
288	const complex double a = twopi_Lt * phys_i * I;	47,680✔
289	const complex double A = a + q;	47,680✔
290	const complex double xup = 0.5 * (1.0 + 1.0 * self->ell + A);	47,680✔
291	const complex double xdw = 0.5 * (2.0 + 1.0 * self->ell - A);	47,680✔
292	const complex double two_x_m1 = cpow (2.0, A - 1.0);	47,680✔
293	complex double U;
294
295	gsl_sf_result lngamma_rho_up, lngamma_theta_up;
296	gsl_sf_result lngamma_rho_dw, lngamma_theta_dw;
297
298	gsl_sf_lngamma_complex_e (creal (xup), cimag (xup), &lngamma_rho_up, &lngamma_theta_up);	47,680✔
299	gsl_sf_lngamma_complex_e (creal (xdw), cimag (xdw), &lngamma_rho_dw, &lngamma_theta_dw);	47,680✔
300
301	U = cexp ((lngamma_rho_up.val - lngamma_rho_dw.val) + I * (lngamma_theta_up.val - lngamma_theta_dw.val));	47,680✔
302
303	/printf ("% 22.15g % 22.15g % 22.15g % 22.15g\n", lngamma_rho_up.val, lngamma_rho_dw.val, lngamma_theta_up.val, lngamma_theta_dw.val);/
304
305	Ym_base[i] = pi_sqrt * two_x_m1 * U;	47,680✔
306	}
307	}
308
309	#endif /* NUMCOSMO_HAVE_FFTW3 */
310	}	58✔
311
312	/**
313	* ncm_fftlog_sbessel_j_new:
314	* @ell: Spherical Bessel Integer order
315	* @lnr0: output center $\ln(r_0)$
316	* @lnk0: input center $\ln(k_0)$
317	* @Lk: input/output interval size
318	* @N: number of knots
319	*
320	* Creates a new fftlog Spherical Bessel J object.
321	*
322	* Returns: (transfer full): a new #NcmFftlogSBesselJ
323	*/
324	NcmFftlogSBesselJ *
325	ncm_fftlog_sbessel_j_new (guint ell, gdouble lnr0, gdouble lnk0, gdouble Lk, guint N)	8✔
326	{
327	NcmFftlogSBesselJ *fftlog_jl = g_object_new (NCM_TYPE_FFTLOG_SBESSEL_J,	8✔
328	"ell", ell,
329	"lnr0", lnr0,
330	"lnk0", lnk0,
331	"Lk", Lk,
332	"N", N,
333	NULL);
334
335	return fftlog_jl;	8✔
336	}
337
338	/**
339	* ncm_fftlog_sbessel_j_set_ell:
340	* @fftlog_jl: a #NcmFftlogSBesselJ
341	* @ell: Spherical Bessel integer order $\ell$
342	*
343	* Sets @ell as the Spherical Bessel integer order $\ell$.
344	*
345	*/
346	void
347	ncm_fftlog_sbessel_j_set_ell (NcmFftlogSBesselJ *fftlog_jl, const guint ell)	8✔
348	{
349	NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);	8✔
350
351	if (self->ell != ell)	8!
352	{
UNCOV 353	NcmFftlog *fftlog = NCM_FFTLOG (fftlog_jl);	×
354
UNCOV 355	self->ell = ell;	×
UNCOV 356	ncm_fftlog_reset (fftlog);	×
357	}
358	}	8✔
359
360	/**
361	* ncm_fftlog_sbessel_j_get_ell:
362	* @fftlog_jl: a #NcmFftlogSBesselJ
363	*
364	* Returns: the current Spherical Bessel integer order $\ell$.
365	*/
366	guint
367	ncm_fftlog_sbessel_j_get_ell (NcmFftlogSBesselJ *fftlog_jl)	×
368	{
369	NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);	×
370
371	return self->ell;	×
372	}
373
374	/**
375	* ncm_fftlog_sbessel_j_set_q:
376	* @fftlog_jl: a #NcmFftlogSBesselJ
377	* @q: Spherical Bessel power factor $q$
378	*
379	* Sets @q as the Spherical Bessel power.
380	*
381	*/
382	void
383	ncm_fftlog_sbessel_j_set_q (NcmFftlogSBesselJ *fftlog_jl, const gdouble q)	8✔
384	{
385	NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);	8✔
386
387	if (self->q != q)	8!
388	{
389	NcmFftlog *fftlog = NCM_FFTLOG (fftlog_jl);	×
390
391	self->q = q;	×
392	ncm_fftlog_reset (fftlog);	×
393	}
394	}	8✔
395
396	/**
397	* ncm_fftlog_sbessel_j_get_q:
398	* @fftlog_jl: a #NcmFftlogSBesselJ
399	*
400	* Returns: the current Spherical Bessel power $q$.
401	*/
402	gdouble
403	ncm_fftlog_sbessel_j_get_q (NcmFftlogSBesselJ *fftlog_jl)	×
404	{
405	NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);	×
406
407	return self->q;	×
408	}
409
410	/**
411	* ncm_fftlog_sbessel_j_set_best_lnr0:
412	* @fftlog_jl: a #NcmFftlogSBesselJ
413	*
414	* Sets the value of $\ln(r_0)$ which gives the best results for
415	* the transformation based on the current value of $\ln(k_0)$,
416	* this is based in the rule of thumb $\mathrm{max}_{x^*}(j_l)$
417	* where $ x^* \approx l + 1$.
418	*
419	*/
420	void
421	ncm_fftlog_sbessel_j_set_best_lnr0 (NcmFftlogSBesselJ *fftlog_jl)	2✔
422	{
423	NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);	2✔
424	NcmFftlog *fftlog = NCM_FFTLOG (fftlog_jl);	2✔
425
426	gint signp = 0;	2✔
427
428	const gdouble lnk0 = ncm_fftlog_get_lnk0 (fftlog);	2✔
429	const gdouble Lk = ncm_fftlog_get_length (fftlog);	2✔
430	const gdouble ell = self->ell;	2✔
431	const gdouble lnc0 = (ell == 0) ? 0.0 : ((ell - 1.0) * Lk + 2.0 * (ell + 1.0) * M_LN2 - ncm_c_lnpi () + 2.0 * lgamma_r (1.5 + ell, &signp)) / (2.0 * (1.0 + ell));	2!
432	const gdouble best_lnr0 = -lnk0 + lnc0;	2✔
433
434	ncm_fftlog_set_lnr0 (fftlog, best_lnr0);	2✔
435	}	2✔
436
437	/**
438	* ncm_fftlog_sbessel_j_set_best_lnk0:
439	* @fftlog_jl: a #NcmFftlogSBesselJ
440	*
441	* Sets the value of $\ln(k_0)$ which gives the best results for
442	* the transformation based on the current value of $\ln(r_0)$,
443	* this is based in the rule of thumb $\mathrm{max}_{x^*}(j_l)$
444	* where $ x^* \approx l + 1$.
445	*
446	*/
447	void
448	ncm_fftlog_sbessel_j_set_best_lnk0 (NcmFftlogSBesselJ *fftlog_jl)	×
449	{
450	NcmFftlogSBesselJPrivate * const self = ncm_fftlog_sbessel_j_get_instance_private (fftlog_jl);	×
451	NcmFftlog *fftlog = NCM_FFTLOG (fftlog_jl);	×
452
453	gint signp = 0;	×
454
455	const gdouble lnr0 = ncm_fftlog_get_lnr0 (fftlog);	×
456	const gdouble Lk = ncm_fftlog_get_length (fftlog);	×
457	const gdouble ell = self->ell;	×
458	const gdouble lnc0 = (ell == 0) ? 0.0 : ((ell - 1.0) * Lk + 2.0 * (ell + 1.0) * M_LN2 - ncm_c_lnpi () + 2.0 * lgamma_r (1.5 + ell, &signp)) / (2.0 * (1.0 + ell));	×
459	const gdouble best_lnk0 = -lnr0 + lnc0;	×
460
461	ncm_fftlog_set_lnk0 (fftlog, best_lnk0);	×
462	}	×
463

NumCosmo / NumCosmo / 6817479045

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous