Get CMB Source function through python

Is there a way to obtain the CMB source function (often denoted $\Delta_l(\eta_0,k)$ or $\Theta_l(\eta_0 , k)$ ) through the python wrapper? I am referring to the following function, which is usually used to obtain the CMB powerspectrum in harmonic space $C_l$:

$$ Cl=\frac{1}{2\pi^2}\int\frac{\mathrm{d}k}{k}P(k)\Theta^2\text{l}(\eta_0,k) $$

I would like to have a function in python that returns $\Theta_l(\eta_0,k)$ as an array over all multipoles $l$ for the argument $k$, something like

def theta(k):
    return cosmo.get_transfer_functions_at_q(k)

I think this is already implemented in CLASS, namely as transfer_functions_at_q in transfer.c, even though the function is never called. Is there a way to make it callable through classy?

Thank you very much in advance for your input. :)

I tried implementing the function I described above, but it does not return the values I was expecting:

huge values for some high l (where I would have expected a strong suppression through the spherical Bessel function),
NaN as max and min values, when iterating through $l$ and $k$ values,
CLASS error when trying to obtain values for wavenumbers q>0.5:

Error in Class: transfer_functions_at_q(L:73) :error in array_interpolate_two( ptr->q, 1, 0, ptr->transfer[index_md] +((index_ic * ptr->tt_size[index_md] + index_tt) * ptr->l_size[index_md] + index_l) * ptr->q_size, 1, ptr->q_size, q, transfer_function, 1, ptr->error_message);
=>array_interpolate_two(L:2783) : x=1.000000e+00 > x_max=5.107672e-01

Question

My question is: is transfer_functions_at_q actually the function I am looking for, or is there another way to obtain the $\Theta_l$ I am interested in?

`classy` Implementation

I made the function accessible through python as follows: In classy.pyx I added the following lines:

def transfer_at_q(self, index_l, q, index_md=0, index_ic=0, index_tt=0):
        """
        transfer_at_q(self, index_l, q, index_md=0, index_ic=0, index_tt=0)

        Transfer function \f$ \Delta_l^{X} (q) \f$ at a given wavenumber q.

         For a given mode (scalar, vector, tensor), initial condition, type
         (temperature, polarization, lensing, etc) and multipole, computes
         the transfer function for the given value of q by interpolating
         between pre-computed values of q.

         Wavenumbers are called q in this module and k in the perturbation
         module. In flat universes k=q. In non-flat universes q and k differ
         through \f$ q2 = k2 + K(1+m)\f$, where m=0,1,2 for scalar, vector,
         tensor.

        Parameters
        ----------
        index_l:    int
                    Define the multipole l
        q:          float
                    The requested wavenumber
        index_md:   int
                    Index of the mode (scalar=0, vector=1, tensor=2)
        index_ic:   int
                    Index of the initial condition (0 for adiabatic)
        index_tt:   int
                    Index of the type (TT=0, EE=1, TE=2, PP=3, TP=4, BB=5)?

        Returns
        -------
        delta_l_q: float
                    The transfer function at the requested wavenumber q

        """

        cdef double delta_l_q

        if transfer_functions_at_q(&self.tr, index_md, index_ic, index_tt, index_l, q, &delta_l_q) == _FAILURE_:
            raise CosmoSevereError(self.tr.error_message)

        return delta_l_q

In cclassy.pxd I added the following:

double transfer_functions_at_q(
                            void * tr,
                            int index_md,
                            int index_ic,
                            int index_tt,
                            int index_l,
                            double q,
                            double * transfer_function
                            )

lesgourg / class_public