diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS index 5719030bf..750d8064c 100644 --- a/.github/CODEOWNERS +++ b/.github/CODEOWNERS @@ -1,3 +1,5 @@ skypy @skypyproject/skypy-infrastructure skypy/galaxy @skypyproject/galaxy-developers +skypy/halos @skypyproject/halo-developers +skypy/power_spectrum @skypyproject/power-spectrum-developers CODE_OF_CONDUCT.md @skypyproject/skypy-board diff --git a/README.rst b/README.rst index c02fcd2be..1869bd98b 100644 --- a/README.rst +++ b/README.rst @@ -8,6 +8,8 @@ This package contains methods for modelling the Universe, galaxies and the Milky Way. Also included are methods for generating observed data. * Galaxies_: morphology, luminosity and redshift distributions +* Halo_ and subhalo mass distributions +* `Power Spectra`_ using CAMB and Halofit * Pipelines_ to generate populations of astronomical objects The full list of features can be found in the `SkyPy Documentation`_. @@ -16,6 +18,8 @@ If you use SkyPy for work or research presented in a publication please follow our `Citation Guidelines`_. .. _Galaxies: https://skypy.readthedocs.io/en/latest/galaxies.html +.. _Halo: https://skypy.readthedocs.io/en/latest/halos/index.html +.. _Power Spectra: https://skypy.readthedocs.io/en/latest/power_spectrum/index.html .. _Pipelines: https://skypy.readthedocs.io/en/latest/pipeline/index.html .. _SkyPy Documentation: https://skypy.readthedocs.io/en/latest/ .. _Citation Guidelines: CITATION.rst diff --git a/docs/halos/examples/halos.yml b/docs/halos/examples/halos.yml new file mode 100644 index 000000000..949b58b22 --- /dev/null +++ b/docs/halos/examples/halos.yml @@ -0,0 +1,24 @@ +cosmology: !astropy.cosmology.default_cosmology.get +power_spectrum: !skypy.power_spectrum.eisenstein_hu + wavenumber: !numpy.logspace [-3, 1, 100] + A_s: 2.1982e-09 + n_s: 0.969453 + cosmology: $cosmology +sheth-tormen: !skypy.halos.mass.sheth_tormen + m_min: 1.0E+10 + m_max: 1.0E+13 + resolution: 10000 + wavenumber: !numpy.logspace [-3, 1, 100] + power_spectrum: $power_spectrum + growth_function: 1.0 + cosmology: $cosmology + size: 10000 +press-schechter: !skypy.halos.mass.press_schechter + m_min: 1.0E+10 + m_max: 1.0E+13 + resolution: 10000 + wavenumber: !numpy.logspace [-3, 1, 100] + power_spectrum: $power_spectrum + growth_function: 1.0 + cosmology: $cosmology + size: 10000 diff --git a/docs/halos/index.rst b/docs/halos/index.rst new file mode 100644 index 000000000..d0d82efad --- /dev/null +++ b/docs/halos/index.rst @@ -0,0 +1,113 @@ +********************************* +Dark Matter Halos (`skypy.halos`) +********************************* + +You can reproduce figure 2 in Sheth and Tormen 1999 +and plot the collapse functions for different halo models. For this, you can use +a python script, for example. + +.. plot:: + :include-source: true + :nofigs: + :context: close-figs + + import numpy as np + from astropy.cosmology import Planck15 + from skypy.power_spectrum import eisenstein_hu + from skypy.halos.mass import _sigma_squared + + # Power spectrum and amplitude of perturbations at redshift 0 + growth_0 = 1.0 + A_s, n_s = 2.1982e-09, 0.969453 + k = np.logspace(-3, 1, num=1000, base=10.0) + mass = 10**np.arange(9.0, 15.0, 0.1) + + pk0 = eisenstein_hu(k, A_s, n_s, Planck15, kwmap=0.02, wiggle=True) + sigma = np.sqrt(_sigma_squared(mass, k, pk0, growth_0, Planck15)) + + # Collapse functions + from skypy.halos import mass + + params_ellipsoidal = (0.3, 0.7, 0.3, 1.686) + + ST = mass.sheth_tormen_collapse_function(sigma) + PS = mass.press_schechter_collapse_function(sigma) + EM = mass.ellipsoidal_collapse_function(sigma, params=params_ellipsoidal) + + +.. plot:: + :include-source: false + :context: close-figs + + import matplotlib.pyplot as plt + + delta_c = 1.69 + nu = np.square(delta_c / sigma) + + # plot different collapse functions + plt.loglog(nu, ST, label='Sheth-Tormen') + plt.loglog(nu, PS, label='Press-Schechter') + plt.loglog(nu, EM, label='Ellipsoidal') + + # axes labels and title + plt.xlabel(r'$\nu \equiv (\delta_c / \sigma)^2$') + plt.ylabel(r'$f_c(\nu)$') + plt.title('Collapse function') + + # show plot labels + plt.legend() + plt.show() + + +You can also sample halos using their mass function. For this, you can use a config +file and run the pipeline, for example. + +.. literalinclude:: examples/halos.yml + :language: yaml + +.. plot:: + :include-source: false + :context: close-figs + + from skypy.pipeline import Pipeline + pipeline = Pipeline.read('examples/halos.yml') + pipeline.execute() + + # Draw from different halo mass samplers + halo_massST = pipeline['sheth-tormen'] + halo_massPS = pipeline['press-schechter'] + + plt.hist(np.log10(halo_massST), histtype='step', label='Sheth-Tormen') + plt.hist(np.log10(halo_massPS), histtype='step', label='Press-Schechter') + + # axis label and title + plt.xlabel(r'$log(mass)$') + plt.title('Halo sampler') + + # show plot labels + plt.legend() + plt.show() + + +Halos (`skypy.halos`) +===================== + +.. automodule:: skypy.halos + + +Abundance Matching (`skypy.halos.abundance_matching`) +===================================================== + +.. automodule:: skypy.halos.abundance_matching + + +Mass (`skypy.halos.mass`) +========================= + +.. automodule:: skypy.halos.mass + + +Quenching (`skypy.halos.quenching`) +=================================== + +.. automodule:: skypy.halos.quenching diff --git a/docs/index.rst b/docs/index.rst index 84ee5e262..11458ce4d 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -32,6 +32,8 @@ Packages :maxdepth: 1 galaxies + halos/index + power_spectrum/index utils/index Pipeline diff --git a/docs/power_spectrum/examples/power_spectrum.yml b/docs/power_spectrum/examples/power_spectrum.yml new file mode 100644 index 000000000..c27afe67d --- /dev/null +++ b/docs/power_spectrum/examples/power_spectrum.yml @@ -0,0 +1,14 @@ +cosmology: !astropy.cosmology.default_cosmology.get +wavenumber: !numpy.logspace [-3, 1, 100] +eisenstein_hu_wiggle: !skypy.power_spectrum.eisenstein_hu + wavenumber: $wavenumber + A_s: 2.1982e-09 + n_s: 0.969453 + cosmology: $cosmology + kwmap: 0.02 + wiggle: True +halofit: !skypy.power_spectrum.halofit_smith + wavenumber: $wavenumber + redshift: 0.0 + linear_power_spectrum: $eisenstein_hu_wiggle + cosmology: $cosmology diff --git a/docs/power_spectrum/index.rst b/docs/power_spectrum/index.rst new file mode 100644 index 000000000..a138981cb --- /dev/null +++ b/docs/power_spectrum/index.rst @@ -0,0 +1,38 @@ +*************************************** +Power Spectrum (`skypy.power_spectrum`) +*************************************** +This module contains methods to model the matter power spectrum. + +SkyPy provides wrappers to a number of external codes for calculating the +matter power spectrum, including `~skypy.power_spectrum.camb` and +`~skypy.power_spectrum.classy`. Here we demonstrate calculating the linear +matter power spectrum using `~skypy.power_spectrum.eisenstein_hu` and the +non-linear corrections using `~skypy.power_spectrum.halofit_smith`: + +.. literalinclude:: examples/power_spectrum.yml + :language: yaml + +.. plot:: + :include-source: false + :context: close-figs + + import matplotlib.pyplot as plt + from skypy.pipeline import Pipeline + + pipeline = Pipeline.read('examples/power_spectrum.yml') + pipeline.execute() + + # Eisenstein and Hu power spectrum and Halofit matter power spectra + k = pipeline['wavenumber'] + power_EH_w = pipeline['eisenstein_hu_wiggle'] + hf_Smith = pipeline['halofit'] + + plt.loglog(k, power_EH_w, label='Eisenstein & Hu') + plt.loglog(k, hf_Smith, '--', label='Halofit') + plt.xlabel(r'Wavenumber $(1/Mpc)$') + plt.ylabel(r'Power spectrum $(Mpc^3)$') + plt.legend(frameon=False, loc='lower left'); + + + +.. automodule:: skypy.power_spectrum diff --git a/setup.cfg b/setup.cfg index d2c0d7a6b..44eb05279 100644 --- a/setup.cfg +++ b/setup.cfg @@ -34,9 +34,11 @@ test = speclite>=0.11 all = speclite>=0.11 + colossus docs = sphinx-astropy matplotlib + colossus speclite>=0.11 [options.package_data] diff --git a/skypy/__init__.py b/skypy/__init__.py index 3e41f6c23..2cbe46e3b 100644 --- a/skypy/__init__.py +++ b/skypy/__init__.py @@ -15,4 +15,6 @@ __all__ = [] from . import galaxies # noqa +from . import halos # noqa from . import pipeline # noqa +from . import power_spectrum # noqa diff --git a/skypy/halo/properties.py b/skypy/halo/properties.py new file mode 100644 index 000000000..e4fa6bd98 --- /dev/null +++ b/skypy/halo/properties.py @@ -0,0 +1,71 @@ +"""Halo properties module. + +This module provides methods to add simple properties to halos + +Models +====== + +.. autosummary:: + :nosignatures: + :toctree: ../api/ + + halo_circular_velocity + +""" + +import numpy as np +from astropy.units import Unit + +__all__ = [ + 'halo_circular_velocity', + ] + +def halo_circular_velocity(halo_virial_mass, Delta_v, redshift, cosmology): + """Halo circular velocity. + This function computes the halo circular velocity, setting it + equal to the virial velocity using equation (3) and footnote 2 from [1]_. + + Parameters + ---------- + halo_virial_mass : (nm,) array_like + Array for the virial mass, in units of solar mass. + Delta_v : (nm,) array_like + The mean overdensity of the halo. + redshift : (nm,) array_like + The redshift of each halo. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing methods for the evolution history of + omega_matter and omega_lambda with redshift. + + Returns + -------- + circular_velocity: (nm,) array_like + Halo circular velocity for a given mass array, cosmology and redshift, in + units of km s-1. + + Examples + --------- + >>> import numpy as np + >>> from skypy.halo import properties + + This example will compute the halo circular velocity, for a Planck15 cosmology at redshift 0. + + >>> from astropy.cosmology import Planck15 + >>> cosmology = Planck15 + >>> halo_virial_mass = 10**np.arange(9.0, 12.0, 2) * Unit.Msun + >>> Delta_v = np.arange(201.0, 201.1, 0.1) + >>> redshift = np.arange(0.3, 1, 0.5) + >>> properties.halo_circular_velocity(halo_virial_mass, Delta_v, redshift, cosmology) + + + References + ---------- + .. [1] Barnes and Haehnelt 2010 equation 3 https://arxiv.org/pdf/1403.1873.pdf + """ + + virial_velocity = 96.6 * Unit('km s-1') * \ + np.power(Delta_v * cosmology.Om0 * cosmology.h**2 / 24.4, 1.0/6.0) * \ + np.sqrt((1 + redshift) / 3.3) * \ + np.power(halo_virial_mass / (1.0e11 * Unit('Msun')), 1.0/3.0) + + return virial_velocity diff --git a/skypy/halos/__init__.py b/skypy/halos/__init__.py new file mode 100644 index 000000000..e5ded3234 --- /dev/null +++ b/skypy/halos/__init__.py @@ -0,0 +1,22 @@ +"""Halos module. + +This module contains methods that model the properties of dark matter halo +populations. + +Models +====== +.. autosummary:: + :nosignatures: + :toctree: ../api/ + + colossus_mf +""" + +__all__ = [ + 'colossus_mf', +] + +from . import abundance_matching # noqa F401,F403 +from . import mass # noqa F401,F403 +from . import quenching # noqa F401,F403 +from ._colossus import colossus_mf # noqa F401,F403 diff --git a/skypy/halos/_colossus.py b/skypy/halos/_colossus.py new file mode 100644 index 000000000..2ecd048d7 --- /dev/null +++ b/skypy/halos/_colossus.py @@ -0,0 +1,205 @@ +"""Colossus dark matter halo properties. + +This module contains functions that interfaces with the external code +`Colossus `_. + +""" + +from astropy.cosmology import z_at_value +from astropy import units +import numpy as np +from scipy import integrate +from skypy.galaxies.redshift import redshifts_from_comoving_density + +__all__ = [ + 'colossus_mass_sampler', +] + +try: + import colossus # noqa F401 +except ImportError: + HAS_COLOSSUS = False +else: + HAS_COLOSSUS = True + + +def colossus_mass_sampler(redshift, model, mdef, m_min, m_max, + cosmology, sigma8, ns, size=None, resolution=1000): + """Colossus halo mass sampler. + + This function generate a sample of halos from a mass function which + is available in COLOSSUS [1]_. + + Parameters + ---------- + redshift : float + The redshift values at which to sample the halo mass. + model : string + Mass function model which is available in colossus. + mdef : str + Halo mass definition for spherical overdensities used by colossus. + m_min, m_max : float + Lower and upper bounds for halo mass in units of Solar mass, Msun. + cosmology : astropy.cosmology.Cosmology + Astropy cosmology object + sigma8 : float + Cosmology parameter, amplitude of the (linear) power spectrum on the + scale of 8 Mpc/h. + ns : float + Cosmology parameter, spectral index of scalar perturbation power spectrum. + size : int, optional + Number of halos to sample. If size is None (default), a single value is returned. + resolution : int, optional + Resolution of the inverse transform sampling spline. Default is 1000. + + Returns + ------- + sample : (size,) array_like + Samples drawn from the mass function, in units of solar masses. + + References + ---------- + .. [1] Diemer B., 2018, ApJS, 239, 35 + + """ + from colossus.cosmology.cosmology import fromAstropy + from colossus.lss import mass_function + fromAstropy(cosmology, sigma8=sigma8, ns=ns) + h0 = cosmology.h + m_h0 = np.logspace(np.log10(m_min*h0), np.log10(m_max*h0), resolution) # unit: Msun/h + dndm = mass_function.massFunction(m_h0, redshift, mdef=mdef, model=model, + q_out='dndlnM', q_in='M')/m_h0 + m = m_h0/h0 + CDF = integrate.cumtrapz(dndm, (m), initial=0) + CDF = CDF / CDF[-1] + n_uniform = np.random.uniform(size=size) + masssample = np.interp(n_uniform, CDF, m) + return masssample + + +@units.quantity_input(sky_area=units.sr) +def colossus_mf_redshift(redshift, model, mdef, m_min, m_max, sky_area, + cosmology, sigma8, ns, resolution=1000, noise=True): + r'''Sample redshifts from a COLOSSUS halo mass function. + + Sample the redshifts of dark matter halos following a mass function + implemented in COLOSSUS [1]_ within given mass and redshift ranges and for + a given area of the sky. + + Parameters + ---------- + redshift : array_like + Input redshift grid on which the mass function is evaluated. Halos are + sampled over this redshift range. + model : string + Mass function model which is available in colossus. + mdef : str + Halo mass definition for spherical overdensities used by colossus. + m_min, m_max : float + Lower and upper bounds for the halo mass in units of Solar mass, Msun. + sky_area : `~astropy.units.Quantity` + Sky area over which halos are sampled. Must be in units of solid angle. + cosmology : Cosmology + Cosmology object to convert comoving density. + sigma8 : float + Cosmology parameter, amplitude of the (linear) power spectrum on the + scale of 8 Mpc/h. + ns : float + Cosmology parameter, spectral index of scalar perturbation power + spectrum. + noise : bool, optional + Poisson-sample the number of halos. Default is `True`. + + Returns + ------- + redshifts : array_like + Redshifts of the halo sample. + + References + ---------- + .. [1] Diemer B., 2018, ApJS, 239, 35 + + ''' + from colossus.cosmology.cosmology import fromAstropy + from colossus.lss.mass_function import massFunction + + # Set the cosmology in COLOSSUS + fromAstropy(cosmology, sigma8, ns) + + # Integrate the mass function to get the number density of halos at each redshift + def dndlnM(lnm, z): + return massFunction(np.exp(lnm), z, 'M', 'dndlnM', mdef, model) + lnmmin = np.log(m_min*cosmology.h) + lnmmax = np.log(m_max*cosmology.h) + density = [integrate.quad(dndlnM, lnmmin, lnmmax, args=(z))[0] for z in redshift] + density = np.array(density) * np.power(cosmology.h, 3) + + # Sample halo redshifts and assign to bins + return redshifts_from_comoving_density(redshift, density, sky_area, cosmology, noise) + + +@units.quantity_input(sky_area=units.sr) +def colossus_mf(redshift, model, mdef, m_min, m_max, sky_area, cosmology, + sigma8, ns, resolution=1000, noise=True): + r'''Sample halo redshifts and masses from a COLOSSUS mass function. + + Sample the redshifts and masses of dark matter halos following a mass + function implemented in COLOSSUS [1]_ within given mass and redshift ranges + and for a given area of the sky. + + Parameters + ---------- + redshift : array_like + Defines the edges of a set of redshift bins for which the mass function + is evaluated. Must be a monotonically-increasing one-dimensional array + of values. Halo redshifts are sampled between the minimum and maximum + values in this array. + model : string + Mass function model which is available in colossus. + mdef : str + Halo mass definition for spherical overdensities used by colossus. + m_min, m_max : float + Lower and upper bounds for the halo mass in units of Solar mass, Msun. + sky_area : `~astropy.units.Quantity` + Sky area over which halos are sampled. Must be in units of solid angle. + cosmology : Cosmology + Cosmology object to calculate comoving densities. + sigma8 : float + Cosmology parameter, amplitude of the (linear) power spectrum on the + scale of 8 Mpc/h. + ns : float + Cosmology parameter, spectral index of scalar perturbation power + spectrum. + noise : bool, optional + Poisson-sample the number of halos. Default is `True`. + + Returns + ------- + redshift, mass : tuple of array_like + Redshifts and masses of the halo sample. + + References + ---------- + .. [1] Diemer B., 2018, ApJS, 239, 35 + + ''' + + # Sample halo redshifts and assign to bins + z = colossus_mf_redshift(redshift, model, mdef, m_min, m_max, sky_area, + cosmology, sigma8, ns, resolution, noise) + redshift_bin = np.digitize(z, redshift) + + # Calculate the redshift at the centre of each bin + comoving_distance = cosmology.comoving_distance(redshift) + d_mid = 0.5 * (comoving_distance[:-1] + comoving_distance[1:]) + z_mid = [z_at_value(cosmology.comoving_distance, d) for d in d_mid] + + # Sample halo masses in each redshift bin + m = np.empty_like(z) + for i, zm in enumerate(z_mid): + mask = redshift_bin == i + 1 + size = np.count_nonzero(mask) + m[mask] = colossus_mass_sampler(zm, model, mdef, m_min, m_max, + cosmology, sigma8, ns, size, resolution) + + return z, m diff --git a/skypy/halos/abundance_matching.py b/skypy/halos/abundance_matching.py new file mode 100644 index 000000000..13cee0b1f --- /dev/null +++ b/skypy/halos/abundance_matching.py @@ -0,0 +1,99 @@ +"""Abundance matching module. + +This module provides methods to perform abundance matching between catalogs of +galaxies and dark matter halos. + +Models +====== + +.. autosummary:: + :nosignatures: + :toctree: ../api/ + + vale_ostriker + +""" + +import numpy as np +from skypy.halos.mass import press_schechter, number_subhalos, subhalo_mass_sampler +from skypy.galaxies.luminosity import schechter_lf_magnitude + +__all__ = [ + 'vale_ostriker', +] + + +def vale_ostriker(halo_kwargs, subhalo_kwargs, galaxy_kwargs, cosmology): + """Vale & Ostriker abundance matching. + Generate matched arrays of (sub)halos masses and galaxy absolute magnitudes + following the abundance matching model in [1]_. + + Parameters + ---------- + halo_kwargs : dict + Dictionary of keyword arguments for `~skypy.halos.mass.press_schechter`. + subhalo_kwargs : dict + Dictionary of keyword arguments for `~skypy.halos.mass.number_subhalos` + and `~skypy.halos.mass.subhalo_mass_sampler`. + galaxy_kwargs : dict + Dictionary of keyword arguments for + `~skypy.galaxies.luminosity.schechter_lf_magnitude`. + cosmology : astropy.cosmology.Cosmology + Cosmology argument for `~skypy.halos.mass.press_schechter` and + `~skypy.galaxies.luminosity.schechter_lf_magnitude`. + + Returns + ------- + mass : array_like + Array of (sub)halo masses in units of solar masses. + group : array_like + Array of halo group ID numbers. + parent : array_like + Array of boolean values indicating if the halo is a parent. + magnitude : array_like + Array of galaxy absolute magnitudes. + + References + ---------- + .. [1] Vale A., Ostriker J. P., 2004, MNRAS, 353, 189 + """ + + # Sample halo and subhalo masses + halo_mass = press_schechter(**halo_kwargs, cosmology=cosmology) + halo_mass[::-1].sort() # Sort in-place from high to low for indexing + n_subhalos = number_subhalos(halo_mass, **subhalo_kwargs) + sampler_kwargs = {k: v for k, v in subhalo_kwargs.items() if k not in ['gamma_M', 'noise']} + subhalo_mass = subhalo_mass_sampler(halo_mass, n_subhalos, **sampler_kwargs) + + # Assign subhalos to groups with their parent halos + n_halos = len(halo_mass) + halo_group = np.arange(n_halos) + indexing = np.zeros(n_halos + 1, dtype=int) + indexing[1:] = np.cumsum(n_subhalos) + total_subhalos = indexing[-1] + subhalo_group = np.empty(total_subhalos, dtype=int) + for first, last, id in zip(indexing[:-1], indexing[1:], halo_group): + subhalo_group[first:last] = id + + # Concatenate halos and subhalos + mass = np.concatenate([halo_mass, subhalo_mass]) + group = np.concatenate([halo_group, subhalo_group]) + parent = np.empty(n_halos+total_subhalos, dtype=bool) + parent[:n_halos] = True + parent[n_halos:] = False + + # Sample galaxy magnitudes + magnitude = schechter_lf_magnitude(**galaxy_kwargs, cosmology=cosmology) + n_galaxies = len(magnitude) + + # Sort halos and galaxies by mass and magnitude + n_matches = min(n_halos + total_subhalos, n_galaxies) + sort_mass = np.argsort(mass)[-n_matches:][::-1] + sort_magnitudes = np.argsort(magnitude)[:n_matches] + + mass = mass[sort_mass] + group = group[sort_mass] + parent = parent[sort_mass] + magnitude = magnitude[sort_magnitudes] + + return mass, group, parent, magnitude diff --git a/skypy/halos/mass.py b/skypy/halos/mass.py new file mode 100644 index 000000000..919b97e1f --- /dev/null +++ b/skypy/halos/mass.py @@ -0,0 +1,448 @@ +r'''Halo mass sampler. +This code samples halos from their mass function. + +Models +====== +.. autosummary:: + :nosignatures: + :toctree: ../api/ + + colossus_mass_sampler + ellipsoidal_collapse_function + halo_mass_function + halo_mass_sampler + number_subhalos + press_schechter + press_schechter_collapse_function + press_schechter_mass_function + sheth_tormen + sheth_tormen_collapse_function + sheth_tormen_mass_function + subhalo_mass_sampler +''' + +import numpy as np +from scipy import integrate +from scipy.special import gamma +from skypy.utils.special import gammaincc +from functools import partial +from astropy import units +from skypy.utils.random import schechter +from ._colossus import colossus_mass_sampler # noqa F401,F403 + +__all__ = [ + 'colossus_mass_sampler', + 'ellipsoidal_collapse_function', + 'halo_mass_function', + 'halo_mass_sampler', + 'number_subhalos', + 'press_schechter', + 'press_schechter_collapse_function', + 'press_schechter_mass_function', + 'sheth_tormen', + 'sheth_tormen_collapse_function', + 'sheth_tormen_mass_function', + 'subhalo_mass_sampler', + ] + + +def halo_mass_function(M, wavenumber, power_spectrum, growth_function, + cosmology, collapse_function, params): + r'''Halo mass function. + This function computes the halo mass function, defined + in equation 7.46 in [1]_. + + Parameters + ----------- + M : (nm,) array_like + Array for the halo mass, in units of solar mass. + wavenumber : (nk,) array_like + Array of wavenumbers at which the power spectrum is evaluated, + in units of Mpc-1. + power_spectrum: (nk,) array_like + Linear power spectrum at redshift 0 in Mpc3. + growth_function : float + The growth function evaluated at a given redshift for the given + cosmology. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing methods for the evolution history of + omega_matter and omega_lambda with redshift. + collapse_function: function + Collapse function to choose from a variety of models: + `sheth_tormen_collapse_function`, `press_schechter_collapse_function`. + params: tuple + List of parameters that determines the model used for + the collapse function. + + Returns + -------- + mass_function: (nm,) array_like + Halo mass function for a given mass array, cosmology and redshift, in + units of Mpc-3 Msun-1. + + Examples + --------- + >>> import numpy as np + >>> from skypy.halos import mass + >>> from skypy.power_spectrum import eisenstein_hu + + This example will compute the halo mass function for elliptical and + spherical collapse, for a Planck15 cosmology at redshift 0. + The power spectrum is given by the Eisenstein and Hu fitting formula: + + >>> from astropy.cosmology import Planck15 + >>> cosmo = Planck15 + >>> D0 = 1.0 + >>> k = np.logspace(-3, 1, num=1000, base=10.0) + >>> A_s, n_s = 2.1982e-09, 0.969453 + >>> Pk = eisenstein_hu(k, A_s, n_s, cosmo, kwmap=0.02, wiggle=True) + + The Sheth and Tormen mass function at redshift 0: + + >>> m = 10**np.arange(9.0, 12.0, 2) + >>> mass.sheth_tormen_mass_function(m, k, Pk, D0, cosmo) + array([3.07523240e-11, 6.11387743e-13]) + + And the Press-Schechter mass function at redshift 0: + + >>> mass.press_schechter_mass_function(m, k, Pk, D0, cosmo) + array([3.46908809e-11, 8.09874945e-13]) + + For any other collapse models: + + >>> params_model = (0.3, 0.7, 0.3, 1.686) + >>> mass.halo_mass_function(m, k, Pk, D0, cosmo, + ... mass.ellipsoidal_collapse_function, params=params_model) + array([2.85598921e-11, 5.67987501e-13]) + + References + ---------- + .. [1] Mo, H. and van den Bosch, F. and White, S. (2010), Cambridge + University Press, ISBN: 9780521857932. + ''' + sigma = np.sqrt(_sigma_squared(M, wavenumber, power_spectrum, + growth_function, cosmology)) + f_c = collapse_function(sigma, params) + + dlognu_dlogm = _dlns_dlnM(sigma, M) + rho_bar = (cosmology.critical_density0.to(units.Msun / units.Mpc ** 3)).value + rho_m0 = cosmology.Om0 * rho_bar + + return rho_m0 * f_c * dlognu_dlogm / np.power(M, 2) + + +def halo_mass_sampler(m_min, m_max, resolution, wavenumber, power_spectrum, + growth_function, cosmology, + collapse_function, params, size=None): + r'''Halo mass sampler. + This function samples haloes from their mass function, + see equation 7.46 in [1]_. + + Parameters + ----------- + m_min, m_max : array_like + Lower and upper bounds for the random variable m. + resolution: int + Resolution of the inverse transform sampling spline. + wavenumber : (nk,) array_like + Array of wavenumbers at which the power spectrum is evaluated, + in units of Mpc-1. + power_spectrum: (nk,) array_like + Linear power spectrum at redshift 0 in Mpc3. + growth_function : float + The growth function evaluated at a given redshift for the given + cosmology. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing methods for the evolution history of + omega_matter and omega_lambda with redshift. + collapse_function: function + Collapse function to choose from a variety of models: + `sheth_tormen_collapse_function`, `press_schechter_collapse_function`. + params: tuple + List of parameters that determines the model used for + the collapse function. + size: int, optional + Output shape of samples. Default is None. + + + Returns + -------- + sample: (size,) array_like + Samples drawn from the mass function, in units of solar masses. + + Examples + --------- + >>> import numpy as np + >>> from skypy.halos import mass + >>> from skypy.power_spectrum import eisenstein_hu + + This example will sample from the halo mass function for + a Planck15 cosmology at redshift 0. The power spectrum is given + by the Eisenstein and Hu fitting formula: + + >>> from astropy.cosmology import Planck15 + >>> cosmo = Planck15 + >>> D0 = 1.0 + >>> k = np.logspace(-3, 1, num=100, base=10.0) + >>> A_s, n_s = 2.1982e-09, 0.969453 + >>> Pk = eisenstein_hu(k, A_s, n_s, cosmo, kwmap=0.02, wiggle=True) + + Sampling from the Sheth and Tormen mass function: + + >>> halo_mass = mass.sheth_tormen(1e9, 1e12, 100, k, Pk, D0, cosmo) + + And from the Press-Schechter mass function: + + >>> halo_mass = mass.press_schechter(1e9, 1e12, 100, k, Pk, D0, cosmo) + + For any other collapse models: + + >>> params_model = (0.3, 0.7, 0.3, 1.686) + >>> halo_mass = mass.halo_mass_sampler(1e9, 1e12, 100, k, Pk, D0, cosmo, + ... mass.ellipsoidal_collapse_function, params=params_model) + + References + ---------- + .. [1] Mo, H. and van den Bosch, F. and White, S. (2010), Cambridge + University Press, ISBN: 9780521857932. + ''' + m = np.logspace(np.log10(m_min), np.log10(m_max), resolution) + + massf = halo_mass_function( + m, wavenumber, power_spectrum, growth_function, + cosmology, collapse_function, params=params) + + CDF = integrate.cumtrapz(massf, m, initial=0) + CDF = CDF / CDF[-1] + n_uniform = np.random.uniform(size=size) + return np.interp(n_uniform, CDF, m) + + +def ellipsoidal_collapse_function(sigma, params): + r'''Ellipsoidal collapse function. + This function computes the mass function for ellipsoidal + collapse, see equation 10 in [1]_ or [2]_. + + Parameters + ----------- + sigma: (ns,) array_like + Array of the mass variance at different scales and at a given redshift. + params: float + The :math:`{A,a,p, delta_c}` parameters of the Sheth-Tormen formalism. + + Returns + -------- + f_c: array_like + Array with the values of the collapse function. + + Examples + --------- + >>> import numpy as np + >>> from skypy.halos import mass + >>> from skypy.power_spectrum import eisenstein_hu + >>> from skypy.power_spectrum import growth_function + + This example will compute the mass function for + ellipsoidal collapse and a Planck15 cosmology at redshift 0. + The power spectrum is given by the Eisenstein and Hu fitting formula: + + >>> from astropy.cosmology import Planck15 + >>> cosmo = Planck15 + >>> D0 = 1.0 + >>> k = np.logspace(-3, 1, num=5, base=10.0) + >>> A_s, n_s = 2.1982e-09, 0.969453 + >>> Pk = eisenstein_hu(k, A_s, n_s, cosmo, kwmap=0.02, wiggle=True) + + The Sheth-Tormen collapse function at redshift 0: + + >>> m = 10**np.arange(9.0, 12.0, 2) + >>> sigma = np.sqrt(_sigma_squared(m, k, Pk, D0, cosmo)) + >>> mass.sheth_tormen_collapse_function(sigma) + array([0.17947815, 0.19952375]) + + And the Press-Schechter collapse function at redshift 0: + + >>> mass.press_schechter_collapse_function(sigma) + array([0.17896132, 0.21613726]) + + For any other collapse models: + + >>> params_model = (0.3, 0.7, 0.3, 1.686) + >>> mass.ellipsoidal_collapse_function(sigma, params=params_model) + array([0.16667541, 0.18529452]) + + References + ---------- + .. [1] R. K. Sheth and G. Tormen, Mon. Not. Roy. Astron. Soc. 308, 119 + (1999), astro-ph/9901122. + .. [2] https://www.slac.stanford.edu/econf/C070730/talks/ + Wechsler_080207.pdf + ''' + A, a, p, delta_c = params + + return A * np.sqrt(2.0 * a / np.pi) * (delta_c / sigma) * \ + np.exp(- 0.5 * a * (delta_c / sigma)**2) * \ + (1.0 + np.power(np.power(sigma / delta_c, 2.0) / a, p)) + + +press_schechter_collapse_function = partial(ellipsoidal_collapse_function, + params=(0.5, 1, 0, 1.69)) +sheth_tormen_collapse_function = partial(ellipsoidal_collapse_function, + params=(0.3222, 0.707, 0.3, 1.686)) +sheth_tormen_mass_function = partial( + halo_mass_function, + collapse_function=ellipsoidal_collapse_function, + params=(0.3222, 0.707, 0.3, 1.686)) +press_schechter_mass_function = partial( + halo_mass_function, + collapse_function=ellipsoidal_collapse_function, + params=(0.5, 1, 0, 1.69)) +sheth_tormen = partial(halo_mass_sampler, + collapse_function=ellipsoidal_collapse_function, + params=(0.3222, 0.707, 0.3, 1.686)) +press_schechter = partial(halo_mass_sampler, + collapse_function=ellipsoidal_collapse_function, + params=(0.5, 1, 0, 1.69)) + + +def _sigma_squared(M, k, Pk, growth_function, cosmology): + M = np.atleast_1d(M)[:, np.newaxis] + + # Growth function + Dz2 = np.power(growth_function, 2) + + # Matter mean density today + rho_bar = (cosmology.critical_density0.to(units.Msun / units.Mpc ** 3)).value + rho_m0 = cosmology.Om0 * rho_bar + + R = np.power(3 * M / (4 * np.pi * rho_m0), 1.0 / 3.0) + top_hat = 3. * (np.sin(k * R) - k * R * np.cos(k * R)) / ((k * R)**3.) + integrand = Pk * np.power(top_hat * k, 2) + + return Dz2 * integrate.simps(integrand, k) / (2. * np.pi**2.) + + +def _dlns_dlnM(sigma, M): + ds = np.gradient(sigma, M) + return np.absolute((M / sigma) * ds) + + +def number_subhalos(halo_mass, alpha, beta, gamma_M, x, m_min, noise=True): + r'''Number of subhalos. + This function calculates the number of subhalos above a given initial mass + for a parent halo of given mass according to the model of Vale & Ostriker + 2004 [1]_ equation (7). The number of subhalos returned can optionally be + sampled from a Poisson distribution with that mean. + + Parameters + ----------- + halo_mass : (nm, ) array_like + The mass of the halo parent, in units of solar mass. + alpha, beta : float + Parameters that determines the subhalo Schechter function. Its the amplitude + is defined by equation (2) in [1]. + gamma_M : float + Present day mass fraction in subhalos. + x : float + Parameter that accounts for the added mass of the original, unstripped + subhalos. + m_min : array_like + Original mass of the least massive subhalo, in units of solar mass. + Current stripped mass is given by :math:`x m_{min}`. + noise : bool, optional + Poisson-sample the number of subhalos. Default is `True`. + + Returns + -------- + nsubhalos: array_like + Array of the number of subhalos assigned to parent halos with mass halo_mass. + + Examples + --------- + >>> import numpy as np + >>> from skypy.halos import mass + + This gives the number of subhalos in a parent halo of mass :math:`10^{12} M_\odot` + + >>> halo, min_sh = 1.0e12, 1.0e6 + >>> alpha, beta, gamma_M = 1.9, 1.0, 0.3 + >>> x = 3 + >>> nsh = mass.number_subhalos(halo, alpha, beta, gamma_M, x, min_sh) + + References + ---------- + .. [1] Vale, A. and Ostriker, J.P. (2005), arXiv: astro-ph/0402500. + ''' + # m_min is the minimum original subhalo mass + x_low = m_min / (x * beta * halo_mass) + # Subhalo amplitude from equation [2] + A = gamma_M / (beta * gamma(2.0 - alpha)) + # The mean number of subhalos above a mass threshold + # can be obtained by integrating equation (1) in [1] + n_subhalos = A * gammaincc(1.0 - alpha, x_low) * gamma(1.0 - alpha) + + # Random number of subhalos following a Poisson distribution + # with mean n_subhalos + return np.random.poisson(n_subhalos) if noise else n_subhalos + + +def subhalo_mass_sampler(halo_mass, nsubhalos, alpha, beta, + x, m_min, resolution=100): + r'''Subhalo mass sampler. + This function samples the original, unstriped masses of subhaloes from the + subhalo mass function of their parent halo with a constant mass stripping + factor given by equation (1) in [1]_ and the upper subhalo mass limit from + [2]_. + + Parameters + ----------- + halo_mass : (nm, ) array_like + The mass of the halo parent, in units of solar mass. + nsubhalos: (nm, ) array_like + Array of the number of subhalos assigned to parent halos with mass `halo_mass`. + alpha, beta : float + Parameters that determines the subhalo Schechter function. Its the amplitude + is defined by equation 2 in [1]. + x : float + Parameter that accounts for the added mass of the original, unstripped + subhalos. + m_min : float + Original mass of the least massive subhalo, in units of solar mass. + Current stripped mass is given by :math:`m_{min} / x`. + resolution: int, optional + Resolution of the inverse transform sampling spline. Default is 100. + + Returns + -------- + sample: (nh, ) array_like + List of original masses drawn from the subhalo mass function for each + parent halo, in units of solar mass. The length corresponds to the total + number of subhalos for all parent halos, i.e. `np.sum(nsubhalos)`. + + Examples + --------- + >>> import numpy as np + >>> from skypy.halos import mass + + This example samples 100 subhalos for a parent halo of mass 1.0E12 Msun: + + >>> halo, min_sh = 1.0e12, 1.0e6 + >>> alpha, beta, gamma_M = 1.9, 1.0, 0.3 + >>> x = 3 + >>> nsh = 100 + >>> sh = mass.subhalo_mass_sampler(halo, nsh, alpha, beta, x, min_sh) + + References + ---------- + .. [1] Vale, A. and Ostriker, J.P. (2004), arXiv: astro-ph/0402500. + .. [2] Vale, A. and Ostriker, J.P. (2004), arXiv: astro-ph/0511816. + ''' + halo_mass = np.atleast_1d(halo_mass) + nsubhalos = np.atleast_1d(nsubhalos) + subhalo_list = [] + for M, n in zip(halo_mass, nsubhalos): + x_min = m_min / (x * beta * M) + x_max = 0.5 / (x * beta) + subhalo_mass = schechter(alpha, x_min, x_max, resolution, size=n, scale=x * beta * M) + subhalo_list.append(subhalo_mass) + return np.concatenate(subhalo_list) diff --git a/skypy/halos/quenching.py b/skypy/halos/quenching.py new file mode 100644 index 000000000..54bee05d0 --- /dev/null +++ b/skypy/halos/quenching.py @@ -0,0 +1,115 @@ +"""Galaxy quenching. + +This module implements models for environment and mass +quenching by dark matter halos. + +Models +====== +.. autosummary:: + :nosignatures: + :toctree: ../api/ + + environment_quenched + mass_quenched +""" + +import numpy as np +from scipy import special + +__all__ = [ + 'environment_quenched', + 'mass_quenched', + ] + + +def environment_quenched(nh, probability): + r'''Environment quenching function. + This function implements the model proposed by A.Amara where the + probability of a subhalo being quenched is a fixed + probability. The model is inspired on [1]_ and [2]_. + + Parameters + ---------- + nh: integer + Number of subhalos. + probability: float + Quenching probability. + + Returns + ------- + quenched: (nh,) array_like, boolean + Boolean array indicating which subhalo's host galaxies are + (satellite) environment-quenched. + + Examples + --------- + + This example shows how many subhalos are environtment quenched (True) + and how many survive (False) from a list of 1000 halos: + + >>> import numpy as np + >>> from skypy.halos.quenching import environment_quenched + >>> from collections import Counter + >>> quenched = environment_quenched(1000, 0.5) + >>> Counter(quenched) + Counter({...}) + + References + ---------- + .. [1] Peng et al. 2010, doi 10.1088/0004-637X/721/1/193. + .. [2] Birrer et al. 2014, arXiv 1401.3162. + + ''' + + return np.random.uniform(size=nh) < probability + + +def mass_quenched(halo_mass, offset, width): + r'''Mass quenching function. + This function implements the model proposed by A.Amara where the + probability of a halo being quenched is related to the error function + of the logarithm of the halo's mass standardised by an offset and width + parameter. The model is inspired on [1]_ and [2]_. + + Parameters + ---------- + halo_mass: (nh,) array_like + Array of halo masses in units of solar mass, Msun. + offset: float + Halo mass in Msun at which quenching probability is 50%. + width: float + Width of the error function. + + Returns + ------- + quenched: (nh,) array_like, boolean + Boolean array indicating which halo's host galaxies are mass-quenched. + + Examples + --------- + + This example shows how many halos are mass quenched (True) + and how many survive (False) from a list of 1000 halos: + + >>> import numpy as np + >>> from astropy import units + >>> from skypy.halos.quenching import mass_quenched + >>> from collections import Counter + >>> offset, width = 1.0e12, 0.5 + >>> halo_mass = np.random.lognormal(mean=np.log(offset), sigma=width, + ... size=1000) + >>> quenched = mass_quenched(halo_mass, offset, width) + >>> Counter(quenched) + Counter({...}) + + References + ---------- + .. [1] Peng et al. 2010, doi 10.1088/0004-637X/721/1/193. + .. [2] Birrer et al. 2014, arXiv 1401.3162. + + ''' + + standardised_mass = np.log10(halo_mass / offset) / width + probability = 0.5 * (1.0 + special.erf(standardised_mass/np.sqrt(2))) + nh = len(halo_mass) + return np.random.uniform(size=nh) < probability diff --git a/skypy/halos/tests/__init__.py b/skypy/halos/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/skypy/halos/tests/test_abundance_matching.py b/skypy/halos/tests/test_abundance_matching.py new file mode 100644 index 000000000..4bb321452 --- /dev/null +++ b/skypy/halos/tests/test_abundance_matching.py @@ -0,0 +1,56 @@ +from astropy.cosmology import default_cosmology +import numpy as np +from skypy.halos.abundance_matching import vale_ostriker + + +def test_vale_ostriker(): + """ Test Vale & Ostriker abundance matching algorithm""" + + cosmology = default_cosmology.get() + k = np.logspace(-4, 2, 20) + Pk = np.array([7.0554997e+02, 1.4269495e+03, 2.8806238e+03, 5.7748426e+03, + 1.1311605e+04, 2.0794882e+04, 3.3334904e+04, 4.1292112e+04, + 3.0335636e+04, 1.6623714e+04, 5.9353726e+03, 1.5235534e+03, + 3.3850042e+02, 6.5466128e+01, 1.1470471e+01, 1.8625402e+00, + 2.8532746e-01, 4.1803753e-02, 5.9166346e-03, 8.1485945e-04]) + + halo_kwargs = {'m_min': 1.0E+9, + 'm_max': 1.0E+12, + 'resolution': 1000, + 'size': 100, + 'wavenumber': k, + 'power_spectrum': Pk, + 'growth_function': 0.40368249700456954, } + subhalo_kwargs = {'alpha': -1.91, + 'beta': 0.39, + 'gamma_M': 0.18, + 'x': 3, + 'm_min': 1.0E+9, + 'noise': True, } + galaxy_kwargs = {'redshift': 1, + 'M_star': -21.07994198, + 'alpha': -0.5, + 'm_lim': 35, + 'size': 100, } + + mass, group, parent, mag = vale_ostriker(halo_kwargs, subhalo_kwargs, galaxy_kwargs, cosmology) + + # Check monotonic mass-magnitude relation + np.testing.assert_array_equal(np.argsort(mass)[::-1], np.argsort(mag)) + + # Check group indexing + assert np.min(group) == 0 + assert np.max(group) == np.unique(group).size - 1 + + # Check minimum and maximum mass of parent halos + np.testing.assert_array_less(halo_kwargs['m_min'], mass[parent]) + np.testing.assert_array_less(mass[parent], halo_kwargs['m_max']) + + # Check minimum and maximum mass of child halos + for m, g in zip(mass[parent], group[parent]): + children = np.logical_and(group == g, ~parent) + np.testing.assert_array_less(subhalo_kwargs['m_min'], mass[children]) + np.testing.assert_array_less(mass[children], 0.5*m) + + # Check magnitude limit + np.testing.assert_array_less(mag, galaxy_kwargs['m_lim']) diff --git a/skypy/halos/tests/test_colossus.py b/skypy/halos/tests/test_colossus.py new file mode 100644 index 000000000..f78155251 --- /dev/null +++ b/skypy/halos/tests/test_colossus.py @@ -0,0 +1,107 @@ +import numpy as np +from scipy import integrate +from scipy.stats import kstest +import pytest +from skypy.halos._colossus import HAS_COLOSSUS + + +@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus') +@pytest.mark.flaky +def test_colossus_mass_sampler(): + from astropy.cosmology import WMAP9 + from colossus.cosmology.cosmology import fromAstropy + from colossus.lss import mass_function + from skypy.halos.mass import colossus_mass_sampler + m_min, m_max, size = 1e+12, 1e+16, 100 + sample = colossus_mass_sampler(redshift=0.1, model='despali16', + mdef='500c', m_min=m_min, m_max=m_max, + cosmology=WMAP9, sigma8=0.8200, ns=0.9608, + size=size, resolution=1000) + assert np.all(sample >= m_min) + assert np.all(sample <= m_max) + assert np.shape(sample) == (size,) + fromAstropy(WMAP9, sigma8=0.8200, ns=0.9608) + h0 = WMAP9.h + m_h0 = np.logspace(np.log10(1e+12*h0), np.log10(1e+16*h0), 1000) + dndm = mass_function.massFunction(m_h0, 0.1, mdef='500c', model='despali16', + q_out='dndlnM', q_in='M')/m_h0 + m = m_h0/h0 + CDF = integrate.cumtrapz(dndm, (m), initial=0) + CDF = CDF / CDF[-1] + + D, p = kstest(sample, lambda t: np.interp(t, m, CDF)) + + assert p > 0.01, 'D = {}, p = {}'.format(D, p) + + +@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus') +@pytest.mark.flaky +def test_colossus_mf_redshift(): + + from skypy.halos._colossus import colossus_mf_redshift + from astropy.cosmology import Planck15 + from astropy import units + from scipy import integrate + from colossus.lss.mass_function import massFunction + + # Parameters + redshift = np.linspace(0., 2., 100) + model, mdef = 'despali16', '500c' + m_min, m_max = 1e+12, 1e+16 + sky_area = 1.0 * units.deg**2 + cosmology = Planck15 + sigma8, ns = 0.82, 0.96 + + # Sample redshifts + z_halo = colossus_mf_redshift(redshift, model, mdef, m_min, m_max, sky_area, + cosmology, sigma8, ns, noise=False) + + # Integrate the mass function to get the number density of halos at each redshift + def dndlnM(lnm, z): + return massFunction(np.exp(lnm), z, 'M', 'dndlnM', mdef, model) + lnmmin = np.log(m_min*cosmology.h) + lnmmax = np.log(m_max*cosmology.h) + density = [integrate.quad(dndlnM, lnmmin, lnmmax, args=(z))[0] for z in redshift] + density = np.array(density) * np.power(cosmology.h, 3) + + # Integrate total number of halos for the given area of sky + density *= (sky_area * cosmology.differential_comoving_volume(redshift)).to_value('Mpc3') + n_halo = np.trapz(density, redshift, axis=-1) + + # make sure noise-free sample has right size + assert np.isclose(len(z_halo), n_halo, atol=1.0) + + # Halo redshift CDF + cdf = density # same memory + np.cumsum((density[1:]+density[:-1])/2*np.diff(redshift), out=cdf[1:]) + cdf[0] = 0 + cdf /= cdf[-1] + + # Check distribution of sampled halo redshifts + D, p = kstest(z_halo, lambda z_: np.interp(z_, redshift, cdf)) + assert p > 0.01, 'D = {}, p = {}'.format(D, p) + + +@pytest.mark.skipif(not HAS_COLOSSUS, reason='test requires colossus') +def test_colossus_mf(): + + from skypy.halos import colossus_mf + from astropy.cosmology import Planck15 + from astropy import units + + # redshift and mass distributions are tested separately + # only test that output is consistent here + + z = np.linspace(0., 1., 100) + model, mdef = 'despali16', '500c' + m_min, m_max = 1e+12, 1e+16 + sky_area = 1.0 * units.deg**2 + cosmo = Planck15 + sigma8, ns = 0.82, 0.96 + z_halo, m_halo = colossus_mf(z, model, mdef, m_min, m_max, sky_area, cosmo, sigma8, ns) + + assert len(z_halo) == len(m_halo) + assert np.all(z_halo >= np.min(z)) + assert np.all(z_halo <= np.max(z)) + assert np.all(m_halo >= m_min) + assert np.all(m_halo <= m_max) diff --git a/skypy/halos/tests/test_import.py b/skypy/halos/tests/test_import.py new file mode 100644 index 000000000..89ec19195 --- /dev/null +++ b/skypy/halos/tests/test_import.py @@ -0,0 +1,5 @@ +def test_import(): + import skypy.halos # noqa: F401 + import skypy.halos.abundance_matching # noqa: F401 + import skypy.halos.mass # noqa: F401 + import skypy.halos.quenching # noqa: F401 diff --git a/skypy/halos/tests/test_mass.py b/skypy/halos/tests/test_mass.py new file mode 100644 index 000000000..84b331c2b --- /dev/null +++ b/skypy/halos/tests/test_mass.py @@ -0,0 +1,149 @@ +import numpy as np +from astropy.cosmology import Planck15 +from astropy.units import allclose +from skypy.power_spectrum import eisenstein_hu + + +import skypy.halos.mass as mass + +# Precomputed values for the test, for a Planck15 cosmology at redshift 0 and a +# power spectrum given by the Eisenstein and Hu fitting formula +# Models: Press-Schechter and Sheth-Tormen +mass_array = [1.00000000e+10, 3.16227766e+10, 1.00000000e+11, 3.16227766e+11, + 1.00000000e+12, 3.16227766e+12, 1.00000000e+13, 3.16227766e+13, + 1.00000000e+14, 3.16227766e+14] +ST_fsigma = [0.22851856, 0.24659885, 0.27788, 0.30138876, 0.31338251, + 0.32004239, 0.31068427, 0.2676586, 0.16714694, 0.0957242] +PS_fsigma = [0.27373595, 0.31173523, 0.38096078, 0.4358908, 0.4651783, + 0.48375949, 0.46491057, 0.37505525, 0.18939675, 0.08453155] +E_fsigma = [0.21224081, 0.22905786, 0.25820042, 0.28018667, 0.29148879, + 0.29805329, 0.28973637, 0.25038455, 0.1574158, 0.09077113] +ST_massf = [5.91918389e-12, 2.04006848e-12, 3.90033154e-13, 3.43105199e-14, + 2.47935792e-15, 3.85619829e-16, 2.77447365e-17, 3.53176402e-18, + 2.69047690e-19, 1.13766655e-20] +PS_massf = [7.09042387e-12, 2.57893025e-12, 5.34717624e-13, 4.96224201e-14, + 3.68030594e-15, 5.82882944e-16, 4.15174583e-17, 4.94886632e-18, + 3.04862034e-19, 1.00464376e-20] +E_massf = [5.49755090e-12, 1.89495500e-12, 3.62410841e-13, 3.18968432e-14, + 2.30614353e-15, 3.59125102e-16, 2.58740463e-17, 3.30383237e-18, + 2.53383973e-19, 1.07880011e-20] + +cosmo = Planck15 +k = np.logspace(-3, 1, num=10, base=10.0) +A_s, n_s = 2.1982e-09, 0.969453 +Pk = eisenstein_hu(k, A_s, n_s, cosmo, kwmap=0.02, wiggle=True) + + +def test_halo_mass_function(): + # Test the output and shape is correct given an array of masses + m_array = np.asarray(mass_array) + + # Any particular ellipsoidal collapse model + params_E = (0.3, 0.7, 0.3, 1.686) + fE = mass.ellipsoidal_collapse_function + array_output_E = mass.halo_mass_function(m_array, k, Pk, 1.0, cosmo, + fE, params=params_E) + assert array_output_E.shape == m_array.shape + assert allclose(array_output_E, E_massf) + + # Sheth and Tormen collapse model + array_output_ST = mass.sheth_tormen_mass_function(m_array, k, Pk, + 1.0, cosmo) + assert array_output_ST.shape == m_array.shape + assert allclose(array_output_ST, ST_massf) + + # Press-Schechter collapse model + array_output_PS = mass.press_schechter_mass_function(m_array, k, Pk, + 1.0, cosmo) + assert array_output_PS.shape == m_array.shape + assert allclose(array_output_PS, PS_massf) + + +def test_halo_mass_sampler(): + # Test the output shape is correct given the sample size + n_samples = 1000 + m_min, m_max, resolution = 10**9, 10**12, 100 + # Any particular ellipsoidal collapse model + params_E = (0.3, 0.7, 0.3, 1.686) + fE = mass.ellipsoidal_collapse_function + array_output_E = mass.halo_mass_sampler(m_min, m_max, resolution, k, Pk, + 1.0, cosmo, fE, params=params_E, + size=n_samples) + assert len(array_output_E) == n_samples + + # Sheth and Tormen collapse model + array_output_PS = mass.sheth_tormen(10**9, 10**12, 100, k, + Pk, 1.0, cosmo, size=n_samples) + + assert len(array_output_PS) == n_samples + + # Press-Schechter collapse model + array_output_PS = mass.press_schechter(10**9, 10**12, 100, k, + Pk, 1.0, cosmo, size=n_samples) + + assert len(array_output_PS) == n_samples + + +def test_ellipsoidal_collapse_function(): + # Test any ellipsoidal model against precomputed values + m_array = np.asarray(mass_array) + sigma = np.sqrt(mass._sigma_squared(m_array, k, Pk, 1.0, cosmo)) + params_E = (0.3, 0.7, 0.3, 1.686) + fE = mass.ellipsoidal_collapse_function(sigma, params_E) + assert allclose(fE, E_fsigma) + + # Test the Sheth and Tormen model against precomputed values + m_array = np.asarray(mass_array) + sigma = np.sqrt(mass._sigma_squared(m_array, k, Pk, 1.0, cosmo)) + fst = mass.sheth_tormen_collapse_function(sigma) + assert allclose(fst, ST_fsigma) + + # Test the Press-Schechter model against precomputed values + m_array = np.asarray(mass_array) + sigma = np.sqrt(mass._sigma_squared(m_array, k, Pk, 1.0, cosmo)) + fps = mass.press_schechter_collapse_function(sigma) + assert allclose(fps, PS_fsigma) + + +def test_number_subhalos(): + # Test analytic solution for the mean number of subhalos + halo_mass = 1.0e12 + shm_min = halo_mass / 100 + alpha, beta, gamma_M, x = 0.0, 0.39, 0.18, 3.0 + nsh_output = mass.number_subhalos(halo_mass, alpha, beta, gamma_M, x, shm_min, noise=False) + nsh_mean = (gamma_M / beta) * np.exp(- shm_min / (x * beta * halo_mass)) + + assert round(nsh_output) == round(nsh_mean) + + # Test for array input of halo parents + halo_parents = np.array([1.0e12, 1.0e14]) + shm_min = 1.0e6 + alpha, beta, gamma_M, x = 1.9, 1.0, 0.3, 1.0 + array_nsh = mass.number_subhalos(halo_parents, alpha, beta, gamma_M, x, shm_min) + + assert len(array_nsh) == len(halo_parents) + + +def test_subhalo_mass_sampler(): + # Test the output shape is correct given the sample size + halo_mass, shm_min = 1.0e12, 1.0e6 + alpha, beta, x = 1.9, 1.0, 1.0 + nsh = 20 + array_output = mass.subhalo_mass_sampler(halo_mass, nsh, alpha, beta, x, shm_min, 100) + + assert len(array_output) == np.sum(nsh) + + # For each halo test that each subhalo satisfy shm_min < m < shm_max + shm_max = 0.5 * halo_mass + + assert np.all(array_output > shm_min) and np.all(array_output < shm_max) + + # Repeat the tests for arrays of halos + halo_mass = np.array([1.0e12, 1.0e14]) + nsh = np.array([10, 100]) + shm_max = 0.5 * halo_mass + array_output = mass.subhalo_mass_sampler(halo_mass, nsh, alpha, beta, x, shm_min, 100) + + assert len(array_output) == np.sum(nsh) + assert np.all(array_output[:10] > shm_min) and np.all(array_output[:10] < 0.5 * halo_mass[0]) + assert np.all(array_output[10:] > shm_min) and np.all(array_output[10:] < 0.5 * halo_mass[1]) diff --git a/skypy/halos/tests/test_quenching.py b/skypy/halos/tests/test_quenching.py new file mode 100644 index 000000000..efa0c674d --- /dev/null +++ b/skypy/halos/tests/test_quenching.py @@ -0,0 +1,33 @@ +import numpy as np +import pytest +from scipy import stats +from collections import Counter + +from skypy.halos.quenching import environment_quenched, mass_quenched + + +@pytest.mark.flaky +def test_environment_quenched(): + # Test the quenching process follows a binomial distribution + n, p = 1000, 0.7 + quenched = environment_quenched(n, p) + number_quenched = Counter(quenched)[True] + + p_value = stats.binom_test(number_quenched, n=n, p=p, + alternative='two-sided') + assert p_value > 0.01 + + +@pytest.mark.flaky +def test_mass_quenched(): + # Test the quenching process follows a binomial distribution if the + # logarithmic halo masses are symmetrical about the offset + n = 1000 + offset, width = 1.0e12, 0.5 + halo_mass = 10 ** np.random.uniform(11, 13, n) + quenched = mass_quenched(halo_mass, offset, width) + number_quenched = Counter(quenched)[True] + + p_value = stats.binom_test(number_quenched, n=n, p=0.5, + alternative='two-sided') + assert p_value > 0.01 diff --git a/skypy/power_spectrum/__init__.py b/skypy/power_spectrum/__init__.py new file mode 100644 index 000000000..226934565 --- /dev/null +++ b/skypy/power_spectrum/__init__.py @@ -0,0 +1,50 @@ +""" +This module contains methods that model the matter power spectrum. + +Linear Power Spectrum +===================== + +.. autosummary:: + :nosignatures: + :toctree: ../api/ + + camb + classy + eisenstein_hu + transfer_no_wiggles + transfer_with_wiggles + + +Nonlinear Power Spectrum +======================== + +.. autosummary:: + :nosignatures: + :toctree: ../api/ + + HalofitParameters + halofit + halofit_smith + halofit_takahashi + halofit_bird + + +Growth Functions +================ + +.. autosummary:: + :nosignatures: + :toctree: ../api/ + + growth_factor + growth_function + growth_function_carroll + growth_function_derivative + +""" + +from ._classy import * # noqa F401,F403 +from ._camb import * # noqa F401,F403 +from ._eisenstein_hu import * # noqa F401,F403 +from ._halofit import * # noqa F401,F403 +from ._growth import * # noqa F401,F403 diff --git a/skypy/power_spectrum/_camb.py b/skypy/power_spectrum/_camb.py new file mode 100644 index 000000000..e7fac5ab0 --- /dev/null +++ b/skypy/power_spectrum/_camb.py @@ -0,0 +1,96 @@ +import numpy as np +from astropy import units + + +__all__ = [ + 'camb', +] + + +def camb(wavenumber, redshift, cosmology, A_s, n_s): + r'''CAMB linear matter power spectrum. + Compute the linear matter power spectrum on a two dimensional grid of + redshift and wavenumber using CAMB [1]_. + + Parameters + ---------- + wavenumber : (nk,) array_like + Array of wavenumbers in units of Mpc-1 at which to + evaluate the linear matter power spectrum. + redshift : (nz,) array_like + Array of redshifts at which to evaluate the linear matter power + spectrum. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing omega_matter, omega_baryon, Hubble + parameter and CMB temperature in the present day + A_s : float + Cosmology parameter, amplitude normalisation of curvature perturbation + power spectrum + n_s : float + Cosmology parameter, spectral index of scalar perturbation power + spectrum + + Returns + ------- + power_spectrum : (nz, nk) array_like + Array of values for the linear matter power spectrum in Mpc3 + evaluated at the input wavenumbers for the given primordial power + spectrum parameters, cosmology. For nz redshifts and nk wavenumbers + the returned array will have shape (nz, nk). + + Examples + -------- + >>> import numpy as np + >>> from astropy.cosmology import default_cosmology + >>> cosmology = default_cosmology.get() + >>> redshift = np.array([0, 1]) + >>> wavenumber = np.array([1.e-2, 1.e-1, 1e0]) + >>> A_s = 2.e-9 + >>> n_s = 0.965 + >>> power_spectrum = camb(wavenumber, redshift, cosmology, A_s, n_s) # doctest: +SKIP + + References + ---------- + .. [1] Lewis, A. and Challinor, A. and Lasenby, A. (2000), + doi : 10.1086/309179. + + ''' + + try: + from camb import CAMBparams, model, get_matter_power_interpolator + except ImportError: + raise Exception("camb is required to use skypy.power_spectrum.camb") + + return_shape = (*np.shape(redshift), *np.shape(wavenumber)) + redshift = np.atleast_1d(redshift) + + h2 = cosmology.h * cosmology.h + + pars = CAMBparams() + pars.set_cosmology(H0=cosmology.H0.value, + ombh2=cosmology.Ob0 * h2, + omch2=cosmology.Odm0 * h2, + omk=cosmology.Ok0, + TCMB=cosmology.Tcmb0.value, + mnu=np.sum(cosmology.m_nu.to_value(units.eV)), + standard_neutrino_neff=cosmology.Neff + ) + + # camb interpolator requires redshifts to be in increasing order + redshift_order = np.argsort(redshift) + wavenumber_order = np.argsort(wavenumber) + + pars.InitPower.ns = n_s + pars.InitPower.As = A_s + + pars.NonLinear = model.NonLinear_none + + pk_interp = get_matter_power_interpolator(pars, + nonlinear=False, + hubble_units=False, k_hunit=False, + kmax=np.max(wavenumber), + zmax=np.max(redshift)) + + pzk = pk_interp.P(redshift[redshift_order], wavenumber[wavenumber_order]) + + return pzk[redshift_order].reshape(return_shape) diff --git a/skypy/power_spectrum/_classy.py b/skypy/power_spectrum/_classy.py new file mode 100644 index 000000000..26bfecc1e --- /dev/null +++ b/skypy/power_spectrum/_classy.py @@ -0,0 +1,90 @@ +import numpy as np + +__all__ = [ + 'classy', +] + + +def classy(wavenumber, redshift, cosmology, **kwargs): + """ Return the CLASS computation of the linear matter power spectrum, on a + two dimensional grid of wavenumber and redshift. + + Additional CLASS parameters can be passed via keyword arguments. + + Parameters + ---------- + wavenumber : (nk,) array_like + Array of wavenumbers in units of Mpc-1 at which to + evaluate the linear matter power spectrum. + redshift : (nz,) array_like + Array of redshifts at which to evaluate the linear matter power + spectrum. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing omega_matter, omega_baryon, Hubble + parameter and CMB temperature in the present day + + Returns + ------- + power_spectrum : (nz, nk) array_like + Array of values for the linear matter power spectrum in Mpc3 + evaluated at the input wavenumbers for the given primordial power + spectrum parameters, cosmology. For nz redshifts and nk wavenumbers + the returned array will have shape (nz, nk). + + Examples + -------- + >>> import numpy as np + >>> from astropy.cosmology import default_cosmology + >>> cosmology = default_cosmology.get() + >>> redshift = np.array([0, 1]) + >>> wavenumber = np.array([1.e-2, 1.e-1, 1e0]) + >>> A_s = 2.e-9 + >>> n_s = 0.965 + >>> z_reio = 10. + >>> classy(wavenumber, redshift, cosmology, A_s, n_s, z_reio) # doctest: +SKIP + array([[2.34758952e+04, 8.70837957e+03], + [3.03660813e+03, 1.12836115e+03], + [2.53124880e+01, 9.40802814e+00]]) + + References + ---------- + doi : 10.1088/1475-7516/2011/07/034 + arXiv: 1104.2932, 1104.2933 + + """ + try: + from classy import Class + except ImportError: + raise Exception("classy is required to use skypy.linear.classy") + + h2 = cosmology.h * cosmology.h + + params = { + 'output': 'mPk', + 'P_k_max_1/Mpc': np.max(wavenumber), + 'z_pk': ', '.join(str(z) for z in np.atleast_1d(redshift)), + 'H0': cosmology.H0.value, + 'omega_b': cosmology.Ob0 * h2, + 'omega_cdm': cosmology.Odm0 * h2, + 'T_cmb': cosmology.Tcmb0.value, + 'N_eff': cosmology.Neff, + } + + params.update(kwargs) + + classy_obj = Class() + classy_obj.set(params) + classy_obj.compute() + + z = np.expand_dims(redshift, (-1,)*np.ndim(wavenumber)) + k = np.expand_dims(wavenumber, (0,)*np.ndim(redshift)) + z, k = np.broadcast_arrays(z, k) + pzk = np.empty(z.shape) + + for i in np.ndindex(*pzk.shape): + pzk[i] = classy_obj.pk_lin(k[i], z[i]) + + if pzk.ndim == 0: + pzk = pzk.item() + + return pzk diff --git a/skypy/power_spectrum/_eisenstein_hu.py b/skypy/power_spectrum/_eisenstein_hu.py new file mode 100644 index 000000000..a1c8ceae9 --- /dev/null +++ b/skypy/power_spectrum/_eisenstein_hu.py @@ -0,0 +1,297 @@ +"""Eisenstein and Hu. + +This implements the Eisenstein and Hu fitting +formula for the matter power spectrum. +""" + +from astropy.utils import isiterable +import numpy as np + + +__all__ = [ + 'eisenstein_hu', + 'transfer_with_wiggles', + 'transfer_no_wiggles', +] + + +def transfer_with_wiggles(wavenumber, A_s, n_s, cosmology, kwmap=0.02): + r''' Eisenstein & Hu transfer function with wiggles. + This function returns the Eisenstein & Hu fitting formula for the transfer + function with baryon acoustic oscillation wiggles. This is described in + [1]_ and [2]_. + + Parameters + ---------- + wavenumber : (nk,) array_like + Array of wavenumbers in units of Mpc-1 at which to evaluate + the linear matter power spectrum. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing omega_matter, omega_baryon, Hubble parameter + and CMB temperature at the present day. + A_s, n_s: float + Amplitude and spectral index of primordial scalar fluctuations. + kwmap : float + WMAP normalization for the amplitude of primordial scalar fluctuations, + as described in [3]_, in units of Mpc-1. + Default is 0.02. + + Returns + ------- + transfer : (nk,) array_like + Transfer function evaluated at the given array of wavenumbers for the + input primordial power spectrum parameters A_s and n_s, cosmology and + kwmap normalization. + + Examples + -------- + + This returns the transfer function with wiggles for the Planck 2015 + cosmology at a given array of wavenumbers: + + >>> import numpy as np + >>> from astropy.cosmology import Planck15 + >>> wavenumber = np.logspace(-3, 1, num=5, base=10.0) + >>> A_s, n_s = 2.1982e-09, 0.969453 + >>> transfer_with_wiggles(wavenumber, A_s, n_s, Planck15, kwmap=0.02) + array([9.92144790e-01, 7.78548704e-01, 1.29998169e-01, 4.63863054e-03, + 8.87918075e-05]) + + References + ---------- + .. [1] Eisenstein D. J., Hu W., ApJ, 496, 605 (1998) + .. [2] Eisenstein D. J., Hu W., ApJ, 511, 5 (1999) + .. [3] Komatsu et al., ApJS, 180, 330 (2009) + + ''' + + if isiterable(wavenumber): + wavenumber = np.asarray(wavenumber) + if np.any(wavenumber <= 0): + raise ValueError('Wavenumbers must be positive') + + om0 = cosmology.Om0 + ob0 = cosmology.Ob0 + + if not ob0: + raise ValueError("Ob0 for input cosmology must be non-zero if " + "wiggles = True") + + h0 = cosmology.H0.value / 100 + Tcmb0 = cosmology.Tcmb0.value + + if not Tcmb0: + raise ValueError("Tcmb0 for input cosmology must be non-zero if" + "wiggles = True") + + ak = wavenumber * h0 + om0h2 = om0 * h0**2 + ob0h2 = ob0 * h0**2 + f_baryon = ob0 / om0 + + # redshift and wavenumber equality + k_eq = 7.46e-2 * om0h2 * (Tcmb0 / 2.7)**-2 + z_eq = 2.5e4 * om0h2 * (Tcmb0 / 2.7)**-4 + + # sound horizon and k_silk + z_drag_b1 = 0.313 * om0h2**-0.419 * (1 + 0.607 * om0h2**0.674) + z_drag_b2 = 0.238 * om0h2**0.223 + z_drag = 1291 * om0h2**0.251 / (1 + 0.659 * om0h2**0.828) \ + * (1 + z_drag_b1 * ob0h2**z_drag_b2) + + r_drag = 31.5 * ob0h2 * (Tcmb0 / 2.7)**-4 * (1000. / z_drag) + r_eq = 31.5 * ob0h2 * (Tcmb0 / 2.7)**-4 * (1000. / z_eq) + + sound_horizon = 2 / (3 * k_eq) * np.sqrt(6 / r_eq) * \ + np.log((np.sqrt(1 + r_drag) + np.sqrt(r_drag + r_eq)) / + (1 + np.sqrt(r_eq))) + k_silk = 1.6 * ob0h2**0.52 * om0h2**0.73 * (1 + (10.4 * om0h2)**-0.95) + + # alpha c + alpha_c_a1 = (46.9 * om0h2)**0.670 * (1 + (32.1 * om0h2)**-0.532) + alpha_c_a2 = (12.0 * om0h2)**0.424 * (1 + (45.0 * om0h2)**-0.582) + alpha_c = alpha_c_a1 ** -f_baryon * alpha_c_a2 ** (-f_baryon**3) + + # beta_c + beta_c_b1 = 0.944 / (1 + (458 * om0h2)**-0.708) + beta_c_b2 = (0.395 * om0h2)**-0.0266 + beta_c = 1 / (1 + beta_c_b1 * ((1 - f_baryon)**beta_c_b2 - 1)) + + y = (1.0 + z_eq) / (1 + z_drag) + alpha_b_G = y * (-6 * np.sqrt(1 + y) + (2 + 3 * y) + * np.log((np.sqrt(1 + y) + 1) / (np.sqrt(1 + y) - 1))) + alpha_b = 2.07 * k_eq * sound_horizon * (1 + r_drag)**-0.75 * alpha_b_G + + beta_node = 8.41 * om0h2 ** 0.435 + beta_b = 0.5 + f_baryon + (3 - 2 * f_baryon) * np.sqrt((17.2 * om0h2)**2 + + 1.0) + + q = ak / (13.41 * k_eq) + ks = ak * sound_horizon + + T_c_ln_beta = np.log(np.e + 1.8 * beta_c * q) + T_c_ln_nobeta = np.log(np.e + 1.8 * q) + T_c_C_alpha = 14.2 / alpha_c + 386. / (1 + 69.9 * q ** 1.08) + T_c_C_noalpha = 14.2 + 386. / (1 + 69.9 * q ** 1.08) + + T_c_f = 1 / (1 + (ks / 5.4) ** 4) + + def f(a, b): + return a / (a + b * q**2) + + T_c = T_c_f * f(T_c_ln_beta, T_c_C_noalpha) + \ + (1 - T_c_f) * f(T_c_ln_beta, T_c_C_alpha) + + s_tilde = sound_horizon * (1 + (beta_node / ks)**3)**(-1 / 3) + ks_tilde = ak * s_tilde + + T_b_T0 = f(T_c_ln_nobeta, T_c_C_noalpha) + T_b_1 = T_b_T0 / (1 + (ks / 5.2)**2) + T_b_2 = alpha_b / (1 + (beta_b / ks)**3) * np.exp(-(ak / k_silk)**1.4) + T_b = np.sinc(ks_tilde / np.pi) * (T_b_1 + T_b_2) + + transfer = f_baryon * T_b + (1 - f_baryon) * T_c + + return transfer + + +def transfer_no_wiggles(wavenumber, A_s, n_s, cosmology): + r'''Eisenstein & Hu transfer function without wiggles. + Eisenstein & Hu fitting formula for the transfer function without + baryon acoustic oscillation wiggles. This is described in + [1]_ and [2]_. + + Parameters + ---------- + wavenumber : (nk,) array_like + Array of wavenumbers in units of Mpc-1 at which to evaluate + the linear matter power spectrum. + A_s, n_s: float + Amplitude and spectral index of primordial scalar fluctuations. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing omega_matter, omega_baryon, Hubble parameter + and CMB temperature in the present day. + + Returns + ------- + transfer : (nk, ) array_like + Transfer function evaluated at the given wavenumbers for the input + primordial power spectrum parameters A_s and n_s, cosmology and kwmap + normalization. + + Examples + -------- + + This returns the transfer function without wiggles for the Planck 2015 + cosmology at a given array of wavenumbers: + + >>> import numpy as np + >>> from astropy.cosmology import Planck15 + >>> wavenumber = np.logspace(-3, 1, num=5, base=10.0) + >>> A_s, n_s = 2.1982e-09, 0.969453 + >>> transfer_no_wiggles(wavenumber, A_s, n_s, Planck15) + array([9.91959695e-01, 7.84518347e-01, 1.32327555e-01, 4.60773671e-03, + 8.78447096e-05]) + + References + ---------- + .. [1] Eisenstein D. J., Hu W., ApJ, 496, 605 (1998) + .. [2] Eisenstein D. J., Hu W., ApJ, 511, 5 (1999) + .. [3] Komatsu et al., ApJS, 180, 330 (2009) + + ''' + + if isiterable(wavenumber): + wavenumber = np.asarray(wavenumber) + if np.any(wavenumber <= 0): + raise ValueError('Wavenumbers must be positive') + + om0 = cosmology.Om0 + ob0 = cosmology.Ob0 + h0 = cosmology.H0.value / 100 + Tcmb0 = cosmology.Tcmb0.value + ak = wavenumber * h0 + om0h2 = om0 * h0**2 + f_baryon = ob0 / om0 + + alpha = 1 - 0.328 * np.log(431 * om0h2) * f_baryon + 0.38 * \ + np.log(22.3 * om0h2) * f_baryon**2 + sound = 44.5 * np.log(9.83 / om0h2) / \ + np.sqrt(1 + 10 * (f_baryon * om0h2)**(0.75)) + shape = om0h2 * (alpha + (1 - alpha) / (1 + (0.43 * ak * sound)**4)) + aq = ak * (Tcmb0 / 2.7)**2 / shape + transfer = np.log(2 * np.e + 1.8 * aq) / \ + (np.log(2 * np.e + 1.8 * aq) + + (14.2 + 731 / (1 + 62.5 * aq)) * aq * aq) + + return transfer + + +def eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap=0.02, wiggle=True): + """ Eisenstein & Hu matter power spectrum. + This function returns the Eisenstein and Hu fitting function for the linear + matter power spectrum with (or without) baryon acoustic oscillations, c.f. + [1]_ and [2]_, using + formulation from Komatsu et al (2009) in [3]_. + + Parameters + ---------- + wavenumber : (nk, ) array_like + Array of wavenumbers in units of Mpc-1 at which to evaluate + the linear matter power spectrum. + A_s, n_s: float + Amplitude and spectral index of primordial scalar fluctuations. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing omega_matter, omega_baryon, Hubble parameter + and CMB temperature in the present day. + kwmap : float + WMAP normalization for the amplitude of primordial scalar fluctuations, + as described in [3], in units of Mpc-1. Default is 0.02. + wiggle : bool + Boolean flag to set the use of baryion acoustic oscillations wiggles. + Default is True, for which the power spectrum is computed with the + wiggles. + + Returns + ------- + power_spectrum : array_like + Linear matter power spectrum in units of Mpc3, + evaluated at the given wavenumbers for the input primordial + power spectrum parameters + A_s and n_s, cosmology, and kwmap normalization. + + Examples + -------- + + This example returns the Eisenstein and Hu matter power spectrum with + baryon acoustic oscillations for the Planck 2015 + cosmology at a given array of wavenumbers: + + >>> import numpy as np + >>> from astropy.cosmology import Planck15 + >>> wavenumber = np.logspace(-3, 1, num=5, base=10.0) + >>> A_s, n_s = 2.1982e-09, 0.969453 + >>> eisenstein_hu(wavenumber, A_s, n_s, Planck15, kwmap=0.02, wiggle=True) + array([6.47460158e+03, 3.71610099e+04, 9.65702614e+03, 1.14604456e+02, + 3.91399918e-01]) + + References + ---------- + .. [1] Eisenstein D. J., Hu W., ApJ, 496, 605 (1998) + .. [2] Eisenstein D. J., Hu W., ApJ, 511, 5 (1999) + .. [3] Komatsu et al., ApJS, 180, 330 (2009) + """ + om0 = cosmology.Om0 + h0 = cosmology.H0.value / 100 + if wiggle: + transfer = transfer_with_wiggles(wavenumber, A_s, n_s, cosmology, + kwmap) + else: + transfer = transfer_no_wiggles(wavenumber, A_s, n_s, cosmology) + + # Eq [74] in [3] + power_spectrum = A_s * (2 * wavenumber**2 * 2998**2 / 5 / om0)**2 * \ + transfer**2 * (wavenumber * h0 / kwmap)**(n_s - 1) * 2 * \ + np.pi**2 / wavenumber**3 + + return power_spectrum diff --git a/skypy/power_spectrum/_growth.py b/skypy/power_spectrum/_growth.py new file mode 100644 index 000000000..22ca95ffa --- /dev/null +++ b/skypy/power_spectrum/_growth.py @@ -0,0 +1,232 @@ +"""Growth function. + +This computes the linear growth function in +perturbation theory. +""" + +from astropy.utils import isiterable +import numpy as np +from scipy import integrate + + +__all__ = [ + 'growth_factor', + 'growth_function', + 'growth_function_carroll', + 'growth_function_derivative', +] + + +def growth_function_carroll(redshift, cosmology): + '''Growth function. + + This function returns the growth function as a function of redshift for a + given cosmology as approximated by Carroll, Press & Turner (1992), + equation 29 in [1]_. + + Parameters + ---------- + redshift : (nz,) array_like + Array of redshifts at which to evaluate the growth function. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing methods for the evolution history of + omega_matter and omega_lambda with redshift. + + Returns + ------- + growth : (nz,) array_like + The growth function evaluated at the input redshifts for the given + cosmology. + + Examples + -------- + + This example returns the growth function for a given array of redshifts + and for the Planck 2015 cosmology: + + >>> import numpy as np + >>> from astropy.cosmology import Planck15 + >>> redshift = np.array([0, 1, 2]) + >>> growth_function_carroll(redshift, Planck15) + array([0.781361..., 0.476280..., 0.327549...]) + + References + ---------- + .. [1] Carroll, M. and Press, W. and Turner, E., (1992), + doi : 10.1146/annurev.aa.30.090192.002435 + ''' + + if isiterable(redshift): + redshift = np.asarray(redshift) + if np.any(redshift < 0): + raise ValueError('Redshifts must be non-negative') + + Om = cosmology.Om(redshift) + Ode = cosmology.Ode(redshift) + Dz = 2.5 * Om / (1 + redshift) + return Dz / (np.power(Om, 4.0/7.0) - Ode + (1 + 0.5*Om) * (1.0 + Ode/70.0)) + + +def growth_factor(redshift, cosmology, gamma=6.0/11.0): + r'''Growth factor. + + Function used to calculate :math:`f(z)`, parametrised growth factor as a + function of redshift, as described in [1]_ equation 17. + + Parameters + ---------- + redshift : (nz,) array_like + Array of redshifts at which to evaluate the growth function. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing methods for the evolution history of + omega_matter and omega_lambda with redshift. + gamma : float + Growth index providing an efficient parametrization of the matter + perturbations. + + Returns + ------- + growth_factor : (nz,) array_like + The redshift scaling of the growth factor. + + Examples + -------- + + This example returns the growth factor for a given array of redshifts + and for a given cosmology: + + >>> import numpy as np + >>> from astropy.cosmology import FlatLambdaCDM + >>> cosmology = FlatLambdaCDM(H0=67.04, Om0=0.3183, Ob0=0.047745) + >>> growth_factor(0, cosmology) + 0.5355746155304598 + + References + ---------- + .. [1] E. V. Linder, Phys. Rev. D 72, 043529 (2005) + ''' + z = redshift + + omega_m_z = cosmology.Om(z) + growth_factor = np.power(omega_m_z, gamma) + + return growth_factor + + +def growth_function(redshift, cosmology, gamma=6.0/11.0, z_upper=1100): + r'''Growth function. + + Function used to calculate :math:`D(z)`, growth function at different + redshifts, as described in [1]_ equation 16. + + Parameters + ---------- + redshift : (nz,) array_like + Array of redshifts at which to evaluate the growth function. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing methods for the evolution history of + omega_matter and omega_lambda with redshift. + gamma : float, optional + Growth index providing an efficient parametrization of the matter + perturbations. Default is the 6/11 LCDM value. + z_upper : float, optional + Redshift for the early-time integral cutoff. Default is 1100. + + Returns + ------- + growth_function : (nz,) array_like + The redshift scaling of the growth function. + + Examples + -------- + + This example returns the growth function for a given array of redshifts + and for a given cosmology: + + >>> import numpy as np + >>> from scipy import integrate + >>> from astropy.cosmology import FlatLambdaCDM + >>> cosmology = FlatLambdaCDM(H0=67.04, Om0=0.3183, Ob0=0.047745) + >>> growth_function(0, cosmology) + 0.7909271056297236 + + References + ---------- + .. [1] E. V. Linder, Phys. Rev. D 72, 043529 (2005) + ''' + + def integrand(x): + integrand = (growth_factor(x, cosmology, gamma) - 1) / (1 + x) + return integrand + + z_flat = np.ravel(redshift) + g_flat = np.empty(z_flat.size) + + for i, z in enumerate(z_flat): + integral = integrate.quad(integrand, z, z_upper)[0] + g = np.exp(integral) + g_flat[i] = g / (1 + z) + + if np.isscalar(redshift): + growth_function = g_flat.item() + else: + growth_function = g_flat.reshape(np.shape(redshift)) + + return growth_function + + +def growth_function_derivative(redshift, cosmology, gamma=6.0/11.0): + r'''First derivative of the growth function. + + Function used to calculate D'(z), derivative of the growth function + with respect to redshift as in [1]_ equation 16. + + Parameters + ---------- + redshift : (nz,) array_like + Array of redshifts at which to evaluate the growth function. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing methods for the evolution history of + omega_matter and omega_lambda with redshift. + gamma : float + Growth index providing an efficient parametrization of the matter + perturbations. + + Returns + ------- + growth_function_derivative : (nz,) array_like + The redshift scaling of the derivative of the growth function. + + Notes + ----- + The first derivative of the growth function, :math:`D(z)`, + with respect to redshift reads + + .. math:: D'(z) = - \frac{D(z) f(z)}{1 + z} \;. + + With :math:`f(z)` the growth factor. + + + Examples + -------- + + This example returns the first derivative of the growth function for a + given array of redshifts and cosmology: + + >>> import numpy as np + >>> from scipy import integrate + >>> from astropy.cosmology import FlatLambdaCDM + >>> cosmology = FlatLambdaCDM(H0=67.04, Om0=0.3183, Ob0=0.047745) + >>> growth_function_derivative(0, cosmology) + -0.42360048051025856 + + References + ---------- + .. [1] E. V. Linder, Phys. Rev. D 72, 043529 (2005) + ''' + z = redshift + + growth_function_derivative = - growth_function(z, cosmology, gamma) * \ + growth_factor(z, cosmology, gamma) / (1.0 + z) + + return growth_function_derivative diff --git a/skypy/power_spectrum/_halofit.py b/skypy/power_spectrum/_halofit.py new file mode 100644 index 000000000..3c2d3941a --- /dev/null +++ b/skypy/power_spectrum/_halofit.py @@ -0,0 +1,217 @@ +from astropy.utils import isiterable +from collections import namedtuple +from functools import partial +import numpy as np +from scipy import optimize + + +__all__ = [ + 'HalofitParameters', + 'halofit', + 'halofit_smith', + 'halofit_takahashi', + 'halofit_bird', +] + + +HalofitParameters = namedtuple( + 'HalofitParameters', + ['a', 'b', 'c', 'gamma', 'alpha', 'beta', 'mu', 'nu', 'fa', 'fb', + 'l', 'm', 'p', 'r', 's', 't']) + +_smith_parameters = HalofitParameters( + [0.1670, 0.7940, 1.6762, 1.8369, 1.4861, -0.6206, 0.0], + [0.3084, 0.9466, 0.9463, -0.9400, 0.0], + [0.3214, 0.6669, -0.2807, -0.0793], + [0.2989, 0.8649, 0.1631], + [-0.1452, 0.3700, 1.3884, 0.0], + [0.0, 0.0, 0.3401, 0.9854, 0.8291, 0.0], + [0.1908, -3.5442], + [1.2857, 0.9589], + [-0.0732, -0.1423, 0.0725], + [-0.0307, -0.0585, 0.0743], + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) + +_takahashi_parameters = HalofitParameters( + [0.2250, 0.9903, 2.3706, 2.8553, 1.5222, -0.6038, 0.1749], + [0.5716, 0.5864, -0.5642, -1.5474, 0.2279], + [0.8161, 2.0404, 0.3698, 0.5869], + [-0.0843, 0.1971, 0.8460], + [-0.1959, 1.3373, 6.0835, -5.5274], + [0.3980, 1.2490, 0.3157, -0.7354, 2.0379, -0.1682], + [0.0, -np.inf], + [3.6902, 5.2105], + [-0.0732, -0.1423, 0.0725], + [-0.0307, -0.0585, 0.0743], + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) + +_bird_parameters = HalofitParameters( + [0.1670, 0.7940, 1.6762, 1.8369, 1.4861, -0.6206, 0.0], + [0.3084, 0.9466, 0.9463, -0.9400, 0.0], + [0.3214, 0.6669, -0.2807, -0.0793], + [0.2224, 1.18075, -0.6719], + [-0.1452, 0.3700, 1.3884, 0.0], + [0.0, 0.0, 0.3401, 0.9854, 0.8291, 0.0], + [0.1908, -3.5442], + [1.2857, 0.9589], + [-0.0732, -0.1423, 0.0725], + [-0.0307, -0.0585, 0.0743], + 2.080, 1.2e-3, 26.3, -6.49, 1.44, 12.4) + + +def halofit(wavenumber, redshift, linear_power_spectrum, + cosmology, parameters): + r'''Computation of the non-linear halo power spectrum. + + This function computes the non-linear halo power spectrum, as a function + of redshift and wavenumbers, following [1]_, [2]_ and [3]_. + + Parameters + ---------- + k : (nk,) array_like + Input wavenumbers in units of Mpc-1. + z : (nz,) array_like + Input redshifts + P : (nz, nk) array_like + Linear power spectrum for given wavenumbers + and redshifts Mpc3. + cosmology : astropy.cosmology.Cosmology + Cosmology object providing method for the evolution of + omega_matter with redshift. + parameters : HalofitParameters + namedtuple containing the free parameters of the model. + + Returns + ------- + pknl : (nz, nk) array_like + Non-linear halo power spectrum in units of Mpc3. + + References + ---------- + .. [1] R. E. Smith it et al., VIRGO Consortium, + Mon. Not. Roy. Astron. Soc. 341, 1311 (2003). + .. [2] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya and M. Oguri, + Astrophys. J. 761, 152 (2012). + .. [3] S. Bird, M. Viel and M. G. Haehnelt, + Mon. Not. Roy. Astron. Soc. 420, 2551 (2012). + + Examples + -------- + >>> import numpy as np + >>> from astropy.cosmology import default_cosmology + >>> from skypy.power_spectrum import growth_function, eisenstein_hu, halofit_smith + >>> k = np.logspace(-4, 2, 100, base=10) + >>> z, A_s, n_s = 0, 2.2e-09, 0.97 + >>> cosmology = default_cosmology.get() + >>> dz = growth_function(z, cosmology) + >>> linear_power = eisenstein_hu(k, A_s, n_s, cosmology) * np.square(dz) + >>> nonlinear_power = halofit_smith(k, z, linear_power, cosmology) + ''' + + # Manage shapes of input arrays + return_shape = np.shape(linear_power_spectrum) + redshift = np.atleast_1d(redshift) + if np.ndim(linear_power_spectrum) == 1: + linear_power_spectrum = linear_power_spectrum[np.newaxis, :] + + # Declaration of variables + if isiterable(redshift): + redshift = np.asarray(redshift) + if isiterable(wavenumber): + wavenumber = np.asarray(wavenumber) + if isiterable(linear_power_spectrum): + linear_power_spectrum = np.asarray(linear_power_spectrum) + if np.any(redshift < 0): + raise ValueError('Redshifts must be non-negative') + if np.any(wavenumber <= 0): + raise ValueError('Wavenumbers must be strictly positive') + if np.any(linear_power_spectrum < 0): + raise ValueError('Linear power spectrum must be non-negative') + if not np.all(sorted(wavenumber) == wavenumber): + raise ValueError('Wavenumbers must be provided in ascending order') + + # Redshift-dependent quantities from cosmology + omega_m_z = cosmology.Om(redshift)[:, np.newaxis] + omega_nu_z = cosmology.Onu(redshift)[:, np.newaxis] + omega_de_z = cosmology.Ode(redshift)[:, np.newaxis] + wp1_z = 1.0 + cosmology.w(redshift)[:, np.newaxis] + ode_1pw_z = omega_de_z * wp1_z + + # Linear power spectrum interpolated at each redshift + lnk = np.log(wavenumber) + k2 = np.square(wavenumber) + k3 = np.power(wavenumber, 3) + dl2kz = (linear_power_spectrum * k3) / (2 * np.pi * np.pi) + + # Integrals required to evaluate Smith et al. 2003 equations C5, C7 & C8 + def integral_kn(lnR, n): + R2 = np.exp(2*lnR)[:, np.newaxis] + integrand = dl2kz * np.power(k2, n/2) * np.exp(-k2*R2) + return np.trapz(integrand, lnk, axis=1) + + # Find root at which sigma^2(R) == 1.0 for each redshift + # Smith et al. 2003 equation C5 & C6 + def log_sigma_squared(lnR): + ik0 = integral_kn(lnR, 0) + return np.log(ik0) + guess = np.zeros_like(redshift) + root = optimize.fsolve(log_sigma_squared, guess) + R = np.exp(root)[:, np.newaxis] + ksigma = 1.0 / R + y = wavenumber / ksigma + + # Evaluate integrals at lnR = root for each redshift + ik0 = integral_kn(root, 0)[:, np.newaxis] + ik2 = integral_kn(root, 2)[:, np.newaxis] + ik4 = integral_kn(root, 4)[:, np.newaxis] + + # Effective spectral index neff and curvature C + # Smith et al. 2003 equations C7 & C8 + neff = (2 * R * R * ik2 / ik0) - 3 + c = (4 * R * R / ik0) * (ik2 + R * R * (ik2 * ik2 / ik0 - ik4)) + + # Smith et al. 2003 equations C9-C16 + # With higher order terms from Takahashi et al. 2012 equations A6-A13 + p = parameters + an = np.power(10, np.polyval(p.a[:5], neff) + p.a[5]*c + p.a[6]*ode_1pw_z) + bn = np.power(10, np.polyval(p.b[:3], neff) + p.b[3]*c + p.a[4]*ode_1pw_z) + cn = np.power(10, np.polyval(p.c[:3], neff) + p.c[3]*c) + gamman = np.polyval(p.gamma[:2], neff) + p.gamma[2]*c + alphan = np.abs(np.polyval(p.alpha[:3], neff) + p.alpha[3]*c) + betan = np.polyval(p.beta[:5], neff) + p.beta[5]*c + mun = np.power(10, np.polyval(p.mu, neff)) + nun = np.power(10, np.polyval(p.nu, neff)) + + # Smith et al. 2003 equations C17 & C18 + fa = np.power(omega_m_z, np.asarray(p.fa)[:, np.newaxis, np.newaxis]) + fb = np.power(omega_m_z, np.asarray(p.fb)[:, np.newaxis, np.newaxis]) + f = np.ones((3, np.size(redshift), 1)) + mask = omega_m_z != 1 + fraction = omega_de_z[mask] / (1.0 - omega_m_z[mask]) + f[:, mask] = fraction * fb[:, mask] + (1.0 - fraction) * fa[:, mask] + + # Massive neutrino terms; Bird et al. 2012 equations A6, A9 and A10 + fnu = omega_nu_z / omega_m_z + Qnu = fnu * (p.l - p.t * (omega_m_z - 0.3)) / (1 + p.m * np.power(y, 3)) + dl2kz = dl2kz * (1 + (p.p * fnu * k2) / (1 + 1.5 * k2)) + betan = betan + fnu * (p.r + p.s * np.square(neff)) + + # Two-halo term, Smith et al. 2003 equation C2 + fy = 0.25 * y + 0.125 * np.square(y) + dq2 = dl2kz * (np.power(1+dl2kz, betan) / (1 + alphan*dl2kz)) * np.exp(-fy) + + # One-halo term, Smith et al. 2003 equations C3 and C4 + # With massive neutrino factor Q_nu, Bird et al. 2012 equation A7 + dh2p = an * np.power(y, 3 * f[0])\ + / (1.0 + bn * np.power(y, f[1]) + np.power(cn * f[2] * y, 3 - gamman)) + dh2 = (1 + Qnu) * dh2p / (1.0 + mun / y + nun / (y * y)) + + # Halofit non-linear power spectrum, Smith et al. 2003 equation C1 + pknl = 2 * np.pi * np.pi * (dq2 + dh2) / k3 + + return pknl.reshape(return_shape) + + +halofit_smith = partial(halofit, parameters=_smith_parameters) +halofit_takahashi = partial(halofit, parameters=_takahashi_parameters) +halofit_bird = partial(halofit, parameters=_bird_parameters) diff --git a/skypy/power_spectrum/tests/__init__.py b/skypy/power_spectrum/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/skypy/power_spectrum/tests/data/camb_result.txt b/skypy/power_spectrum/tests/data/camb_result.txt new file mode 100644 index 000000000..62f5af514 --- /dev/null +++ b/skypy/power_spectrum/tests/data/camb_result.txt @@ -0,0 +1,2 @@ +1.938237119497023514e+03 2.033480877732498129e+03 2.133392361215547226e+03 2.238199557662894222e+03 2.348140165780609550e+03 2.463461934271990685e+03 2.584425564318605211e+03 2.711305353383179863e+03 2.844386600502214606e+03 2.983966037629961647e+03 3.130356118426995636e+03 3.283885706961626056e+03 3.444896196361852617e+03 3.613742071863042383e+03 3.790796226269676026e+03 3.976450602317721859e+03 4.171111310496261467e+03 4.375199375636912919e+03 4.589156816082082514e+03 4.813447196225026346e+03 5.048550233739435498e+03 5.294962698503289175e+03 5.553203973099104587e+03 5.823816572887521943e+03 6.107362241921599889e+03 6.404422796309818295e+03 6.715602453952121323e+03 7.041528422230898286e+03 7.382849443325404536e+03 7.740236529852999411e+03 8.114383683174457474e+03 8.506008264109852462e+03 8.915848229233473830e+03 9.344662837033341930e+03 9.793233172329199078e+03 1.026236227395377318e+04 1.075287189380179552e+04 1.126560252952013070e+04 1.180140871696771683e+04 1.236115934974771517e+04 1.294573542229604936e+04 1.355602931667879238e+04 1.419293683224413871e+04 1.485735687757466258e+04 1.555018532401242010e+04 1.627231401032959184e+04 1.702462246242709807e+04 1.780797614018037348e+04 1.862321415258925117e+04 1.947114783561477088e+04 2.035255057221979223e+04 2.126815437614812981e+04 2.221863026244337016e+04 2.320458680502490097e+04 2.422655895182141467e+04 2.528500405861194304e+04 2.638028735243949995e+04 2.751267521303944523e+04 2.868230568908906935e+04 2.988918762532247274e+04 3.113319216862500616e+04 3.241404309916475540e+04 3.373127907583692286e+04 3.508425232988927019e+04 3.647210564252726181e+04 3.789376639598687325e+04 3.934790135253586050e+04 4.083291748009042931e+04 4.234691232737500832e+04 4.388768924672165303e+04 4.545276065604363976e+04 4.703931487822531199e+04 4.864400804544963466e+04 5.026311085654659109e+04 5.189291114155262767e+04 5.352937212295935024e+04 5.516690721634055808e+04 5.679906706545926863e+04 5.841962112008470285e+04 6.002224786865645729e+04 6.160034261237797909e+04 6.314687867814932542e+04 6.465366680613913195e+04 6.611179071123150061e+04 6.751165491678293620e+04 6.884318637208463042e+04 7.009572362174701993e+04 7.125824591953711933e+04 7.231925810725476185e+04 7.326700014909518359e+04 7.408900398171934648e+04 7.477267488089254766e+04 7.530529293028927350e+04 7.567444039118601358e+04 7.586781909309758339e+04 7.587387361757487815e+04 7.568188409962291189e+04 7.528254519445350161e+04 7.466827712552377488e+04 7.383369214652816299e+04 7.277626744111823791e+04 7.149651414886189741e+04 6.999922189002826053e+04 6.829293048596933659e+04 6.639181805428148073e+04 6.431404652589472244e+04 6.208430789542764251e+04 5.973076499842629710e+04 5.728850234947272838e+04 5.479409780120193318e+04 5.228838819919160596e+04 4.981179316813583864e+04 4.740984297836804762e+04 4.512087838092840684e+04 4.296838895033735753e+04 4.097437420007140463e+04 3.917681738631798362e+04 3.759408894518086890e+04 3.618954489654158533e+04 3.492138575760865933e+04 3.375124656008202874e+04 3.263532396794898523e+04 3.151861714397486503e+04 3.032874809207860017e+04 2.898592463560428223e+04 2.743618690896175031e+04 2.567485616834729444e+04 2.375313608352699521e+04 2.176574652607356620e+04 1.982156135466891283e+04 1.802855502204020013e+04 1.648584998787313452e+04 1.526245648914335470e+04 1.435930307142175116e+04 1.369156249758880585e+04 1.311307356588952280e+04 1.246842268501784565e+04 1.165150226515151189e+04 1.065089294878422879e+04 9.552852823646444449e+03 8.512766571801246755e+03 7.689387883508730738e+03 7.126239184267302335e+03 6.728480753163247755e+03 6.346269348664473000e+03 5.847220679792134433e+03 5.242262135427591602e+03 4.663497251668834906e+03 4.225393261710382831e+03 3.923343745987492639e+03 3.644788834620067973e+03 3.305362436940259613e+03 2.950085713003068577e+03 2.664102401738479784e+03 2.447719375757351372e+03 2.232660257730331523e+03 2.002799949971260048e+03 1.804485614656616235e+03 1.642591298910032492e+03 1.483817687190139623e+03 1.332625688197262662e+03 1.204258152905350926e+03 1.086007144020557007e+03 9.756008619701484577e+02 8.786693625790009037e+02 7.895858760686617188e+02 7.090076882760166654e+02 6.366239446159752333e+02 5.708447618392178811e+02 5.118166720036654169e+02 4.585337540134033816e+02 4.104785671321039899e+02 3.672044322748638479e+02 3.282914405032612990e+02 2.933432063424322109e+02 2.619767200135979692e+02 2.338392824752652928e+02 2.086115041644046642e+02 1.860043549136560728e+02 1.657579017463316120e+02 1.476381055732929042e+02 1.314325250447624285e+02 1.169484506462966067e+02 1.040111608198491950e+02 9.246213060839650666e+01 8.215810665546331393e+01 7.296987630955514703e+01 6.478098473720514505e+01 5.748654483751767685e+01 5.099216596582981254e+01 4.521298338871749678e+01 4.007275335221453361e+01 3.550302140904990011e+01 3.144235814230841086e+01 2.783569442230885471e+01 2.463370031501177237e+01 2.179221524325155812e+01 1.927173113130510984e+01 1.703692405555526079e+01 1.505623024674271448e+01 +7.176815422031822891e+02 7.529516239833551481e+02 7.899507365680883595e+02 8.287635290195048583e+02 8.694781463118131342e+02 9.121863388292655372e+02 9.569847783726260104e+02 1.003975317113930373e+03 1.053263888675638327e+03 1.104960659240095310e+03 1.159181568159036487e+03 1.216048588530812822e+03 1.275688509824594576e+03 1.338233150685432975e+03 1.403821020471465772e+03 1.472597559326912005e+03 1.544713646482534841e+03 1.620325881713448553e+03 1.699598604361462549e+03 1.782704099061792704e+03 1.869820676975560445e+03 1.961133014122634904e+03 2.056834207282181069e+03 2.157125958775804520e+03 2.262216969410314505e+03 2.372323267571116048e+03 2.487669282258685598e+03 2.608488032106867649e+03 2.735020223299448389e+03 2.867514521228495141e+03 3.006227754413428556e+03 3.151425116315464948e+03 3.303379930013642024e+03 3.462373712809760946e+03 3.628694742113631492e+03 3.802638372856077694e+03 3.984508060103565185e+03 4.174615057643571163e+03 4.373275147143160211e+03 4.580808957835160072e+03 4.797542778952536537e+03 5.023807975067462394e+03 5.259936571046727295e+03 5.506261468025467366e+03 5.763116365850467446e+03 6.030834678510655067e+03 6.309742189530077667e+03 6.600157574259610556e+03 6.902392996469196987e+03 7.216752310388699698e+03 7.543521735151048233e+03 7.882970020405715331e+03 8.235346525683242362e+03 8.600879529078552878e+03 8.979768955442204970e+03 9.372185265092310146e+03 9.778264010367980518e+03 1.019810382899468459e+04 1.063175766637230481e+04 1.107923155909895104e+04 1.154047842696860062e+04 1.201539563383648601e+04 1.250381453902390967e+04 1.300549876812895855e+04 1.352013238015632669e+04 1.404731861402069262e+04 1.458656471325802158e+04 1.513728201441073725e+04 1.569876755327228602e+04 1.627020940052037986e+04 1.685068659212797138e+04 1.743915755680112125e+04 1.803438477029585920e+04 1.863498893415778366e+04 1.923960028568133202e+04 1.984672907190576007e+04 2.045430128688762125e+04 2.105991930192214568e+04 2.166127417750818859e+04 2.225602627074845077e+04 2.284172872810389163e+04 2.341577736755484148e+04 2.397513807291571356e+04 2.451650730857579401e+04 2.503632836410122763e+04 2.553086602041504375e+04 2.599616279860441500e+04 2.642812695840868764e+04 2.682249649764079004e+04 2.717491227604469168e+04 2.748074806495666053e+04 2.773532783061756709e+04 2.793392501942117815e+04 2.807192402717313962e+04 2.814475803501449991e+04 2.814813560779631734e+04 2.807806949779771094e+04 2.793109452756302926e+04 2.770438411811922197e+04 2.739592469740668457e+04 2.700476824808461242e+04 2.653109209266393736e+04 2.597665664964097596e+04 2.534461430264738738e+04 2.464021541153754515e+04 2.387019045218198153e+04 2.304369193572026415e+04 2.217116235477637383e+04 2.126561766396100211e+04 2.034063058087829268e+04 1.941135496399542535e+04 1.849279348095068053e+04 1.760185306669557031e+04 1.675277454950492756e+04 1.595428466518949790e+04 1.521455840412883845e+04 1.454770890480856178e+04 1.396056825600984848e+04 1.343953801421319076e+04 1.296910354641307276e+04 1.253502425005816076e+04 1.212103445402327998e+04 1.170670807537353539e+04 1.126516528472776008e+04 1.076676305431263972e+04 1.019145379846073592e+04 9.537496008884312687e+03 8.823902572883138419e+03 8.085857061351669472e+03 7.363812492939808180e+03 6.697882775753157148e+03 6.124904150814820241e+03 5.670522442172639785e+03 5.335096438965158086e+03 5.087119124798823577e+03 4.872283976379614614e+03 4.632853052553450652e+03 4.329394818557745566e+03 3.957667300437422455e+03 3.549718024855537806e+03 3.163287086144263867e+03 2.857369671889140591e+03 2.648142452411787417e+03 2.500367499347449211e+03 2.358364267621969702e+03 2.172937343709643301e+03 1.948145014573460458e+03 1.733081348791279197e+03 1.570285534033192562e+03 1.458048195538035770e+03 1.354539518933324189e+03 1.228405921555774739e+03 1.096379053268354937e+03 9.901024298876520788e+02 9.096905982651049953e+02 8.297692101845592560e+02 7.443458318945880592e+02 6.706452813726227760e+02 6.104795108035727935e+02 5.514728263910532178e+02 4.952831359210021560e+02 4.475757716932886296e+02 4.036279116100415649e+02 3.625953271144339283e+02 3.265704165570433020e+02 2.934620405706189672e+02 2.635145609452568465e+02 2.366124947679755337e+02 2.121649202393778069e+02 1.902265416790131098e+02 1.704232825308015151e+02 1.525628513712574090e+02 1.364792977614004030e+02 1.220166020433432124e+02 1.090274684080932133e+02 9.736954473740587446e+01 8.691174008573896970e+01 7.753534586815945318e+01 6.913293889969325789e+01 6.160791462299716414e+01 5.487329482869714781e+01 4.885012677493897115e+01 4.346678988662728926e+01 3.865834761419159094e+01 3.436588108661515406e+01 3.053614710697587498e+01 2.712112112603330871e+01 2.407752039804179134e+01 2.136636186453862152e+01 1.895256383539820533e+01 1.680458494636301126e+01 1.489408785003885249e+01 1.319563013091273262e+01 1.168638024909456341e+01 1.034587054961751917e+01 9.155766446575828610e+00 8.099654638551152530e+00 7.162851014683812245e+00 6.332226587634928627e+00 5.596049909588766802e+00 diff --git a/skypy/power_spectrum/tests/data/class_result.txt b/skypy/power_spectrum/tests/data/class_result.txt new file mode 100644 index 000000000..a7573522c --- /dev/null +++ b/skypy/power_spectrum/tests/data/class_result.txt @@ -0,0 +1,2 @@ +2.203299707906652657e+03 2.311210914981189035e+03 2.424393684325056256e+03 2.543089193030403749e+03 2.667574062997171495e+03 2.798140540231350315e+03 2.935080938986529873e+03 3.078687638310352668e+03 3.229271237537464003e+03 3.387185448005580838e+03 3.552790215760692263e+03 3.726445486848785094e+03 3.908513461867943988e+03 4.099397373680083547e+03 4.299536247108531825e+03 4.509370307330771539e+03 4.729339779524279948e+03 4.959897197280032742e+03 5.201556695420678807e+03 5.454851640441136624e+03 5.720315405246146838e+03 5.998481828323795526e+03 6.289921141018928211e+03 6.595265473515572012e+03 6.915152978353057733e+03 7.250221808070721636e+03 7.601115871995027192e+03 7.968544737566537151e+03 8.353259739478413394e+03 8.756012882388815342e+03 9.177556172713413616e+03 9.618663422506968345e+03 1.008018464209605372e+04 1.056298647603732388e+04 1.106793556888744570e+04 1.159589983541520269e+04 1.214779161177151218e+04 1.272457822412805035e+04 1.332723041282564009e+04 1.395671891820516612e+04 1.461401776069832886e+04 1.530012587190215527e+04 1.601605361603371057e+04 1.676281141722149187e+04 1.754140968777489979e+04 1.835285471050156411e+04 1.919814258993274780e+04 2.007826789029309293e+04 2.099422517580720159e+04 2.194700397868958680e+04 2.293749637501292818e+04 2.396650617334505296e+04 2.503483399362274213e+04 2.614328045578273668e+04 2.729259429677140724e+04 2.848325299910834656e+04 2.971564569218039469e+04 3.099016145726011018e+04 3.230718663595217004e+04 3.366686606780228612e+04 3.506891752657855977e+04 3.651301510933019745e+04 3.799883291310634377e+04 3.952599173980314663e+04 4.109347187926741026e+04 4.269983021441728488e+04 4.434361609991840669e+04 4.602337888221160392e+04 4.773735819384234492e+04 4.948266387547630438e+04 5.125614822379151155e+04 5.305466353546595928e+04 5.487503947816657455e+04 5.671324621678657422e+04 5.856414893918840244e+04 6.042254013598637539e+04 6.228321229779491841e+04 6.414073675795610325e+04 6.598792490647920931e+04 6.781674740951151762e+04 6.961916971888908301e+04 7.138715543467328825e+04 7.311184622833224421e+04 7.478227292396956182e+04 7.638713219119867426e+04 7.791512069963278191e+04 7.935487574958508776e+04 8.069380910843485617e+04 8.191818006509119004e+04 8.301420316735940287e+04 8.396809296304482268e+04 8.476596707271525520e+04 8.539341364253904612e+04 8.583584092846100975e+04 8.607865703851671424e+04 8.610729170393798267e+04 8.590933926982842968e+04 8.547637465997479740e+04 8.480040042909038311e+04 8.387341913188845501e+04 8.268824371897983656e+04 8.124795941468629462e+04 7.956270677115973376e+04 7.764276474784030870e+04 7.549841232508860412e+04 7.314485985675227130e+04 7.061609340540325502e+04 6.795056684212065011e+04 6.518673403798350046e+04 6.236314876564850420e+04 5.952225037343182339e+04 5.671152681037741422e+04 5.397880354273861303e+04 5.137190569198501908e+04 4.893301963605062338e+04 4.668438320178809954e+04 4.464382491179677163e+04 4.282286733123290469e+04 4.120561622712959797e+04 3.976333426253815560e+04 3.844583397269744455e+04 3.718834540955749253e+04 3.591267006813910848e+04 3.453735368925759394e+04 3.299218647913668246e+04 3.122768353566051883e+04 2.923800656535639791e+04 2.706650664939109629e+04 2.480694198903926372e+04 2.258599681029839485e+04 2.054219970322656809e+04 1.879048011091948501e+04 1.739691927815301096e+04 1.635182855588187158e+04 1.557121475777582418e+04 1.490506756614851656e+04 1.417854950618025759e+04 1.325666563576855151e+04 1.211526599709844231e+04 1.086268048974112753e+04 9.686411525071711367e+03 8.750578292490437889e+03 8.103579862276831591e+03 7.651435685782299515e+03 7.216533364341084962e+03 6.646355637802166711e+03 5.956882203927921182e+03 5.299959839175694469e+03 4.802199836566413978e+03 4.457697701354544733e+03 4.139654836024139513e+03 3.754950150243903863e+03 3.349971739425377109e+03 3.025732553133077545e+03 2.779480408073161470e+03 2.534085896673992011e+03 2.273283519140113185e+03 2.048407774800599327e+03 1.863905651898254064e+03 1.683239975066495617e+03 1.512103238244111026e+03 1.365845002323603239e+03 1.231687614022032221e+03 1.106296770328023285e+03 9.962316709697163333e+02 8.951082549061678719e+02 8.036204997247619986e+02 7.214650741281963064e+02 6.468558396140351761e+02 5.798436574496209914e+02 5.192541227897269209e+02 4.647884285137804454e+02 4.157953487813219340e+02 3.717093036554255718e+02 3.320982192367305856e+02 2.965471788441847139e+02 2.646499678715630921e+02 2.360525028365239564e+02 2.104401004508137305e+02 1.875006872124833137e+02 1.669658283411027071e+02 1.486143340446830621e+02 1.322273745219924024e+02 1.175865065563287004e+02 1.045019894751125094e+02 9.283206557884389554e+01 8.243954197113073690e+01 7.318722575558081189e+01 6.494086889808514229e+01 5.759530676072823496e+01 5.106195070551011383e+01 4.525240484419311571e+01 4.007827884660451900e+01 3.546218033477512677e+01 3.136071796576518977e+01 2.773705079657367989e+01 2.455433788420007346e+01 2.177499035005325823e+01 1.934046371409294096e+01 1.716893289416953650e+01 +8.169921322832720989e+02 8.570108171563794031e+02 8.989851780876145995e+02 9.430047744261788694e+02 9.891723745659867291e+02 1.037596569460824412e+03 1.088385975497979416e+03 1.141649233256627895e+03 1.197501786149590316e+03 1.256075174141015395e+03 1.317503265022443884e+03 1.381919926585416988e+03 1.449459869676962626e+03 1.520273148535485007e+03 1.594523201259327834e+03 1.672373914801023602e+03 1.753989176113104577e+03 1.839537471519060546e+03 1.929210306305467157e+03 2.023206372192383924e+03 2.121724363295075818e+03 2.224963146546800999e+03 2.333135097397821482e+03 2.446475567277645951e+03 2.565222143029439394e+03 2.689612411496368622e+03 2.819886062058289099e+03 2.956306764240710436e+03 3.099153442116342831e+03 3.248705264447505215e+03 3.405241400627338408e+03 3.569048846695836346e+03 3.740441948359834441e+03 3.919741021811328665e+03 4.107266383242315896e+03 4.303338818269211515e+03 4.508295529043181887e+03 4.722494040036389379e+03 4.946293137474315699e+03 5.180051607582468932e+03 5.424129666574232033e+03 5.678898390930728965e+03 5.944733841296318133e+03 6.222012104434339562e+03 6.511109261919752498e+03 6.812399582546237980e+03 7.126252867016717573e+03 7.453038239865608375e+03 7.793124825627318387e+03 8.146879881969554845e+03 8.514634510506639344e+03 8.896687065767293461e+03 9.293334719303858947e+03 9.704874642668663000e+03 1.013158474263526841e+04 1.057364220594392282e+04 1.103119141275331094e+04 1.150437672535670390e+04 1.199334149305647043e+04 1.249813976988283321e+04 1.301866770249661022e+04 1.355480528811245858e+04 1.410643252394500450e+04 1.467340939309225178e+04 1.525535534490507234e+04 1.585173082523826088e+04 1.646199345283542425e+04 1.708560084350971374e+04 1.772190026342948477e+04 1.836983643908216254e+04 1.902826233520995811e+04 1.969603091655509343e+04 2.037198615282936225e+04 2.105463036144880607e+04 2.174202663153278263e+04 2.243220915504928780e+04 2.312321212396638657e+04 2.381298995945843126e+04 2.449886227618843623e+04 2.517784544213424670e+04 2.584695394448654770e+04 2.650320161529113102e+04 2.714331149398989510e+04 2.776325981694914299e+04 2.835890459900706628e+04 2.892610385500179473e+04 2.946069047374664297e+04 2.995797867917509575e+04 3.041279494814273494e+04 3.081994682236396693e+04 3.117424184355315447e+04 3.147046641479885147e+04 3.170329146740376382e+04 3.186734870084862632e+04 3.195726978235703791e+04 3.196769420075986272e+04 3.189404443550058204e+04 3.173318283286548467e+04 3.148212642308034265e+04 3.113789223637092437e+04 3.069779935037323594e+04 3.016299548213412709e+04 2.953725798133557328e+04 2.882441578425204716e+04 2.802829783490212867e+04 2.715455851025648735e+04 2.621580233971173220e+04 2.522628768904375102e+04 2.420027292402841704e+04 2.315205385332625883e+04 2.209738259054013542e+04 2.105390345337244435e+04 2.003938726020992908e+04 1.907160470097129655e+04 1.816622544309325531e+04 1.733148622963127127e+04 1.657398087604394095e+04 1.589796680450112399e+04 1.529754335670806176e+04 1.476205596574797892e+04 1.427292007941060729e+04 1.380614629006946598e+04 1.333264860973936447e+04 1.282197696672218080e+04 1.224824856440880831e+04 1.159317897571869798e+04 1.085434650644811700e+04 1.004861605047984267e+04 9.209416660497074190e+03 8.384955211205602609e+03 7.626166632599562945e+03 6.975925342037266091e+03 6.458500695379613717e+03 6.070530004762256794e+03 5.780781214476804962e+03 5.533570065933360638e+03 5.263738286654737749e+03 4.921557490086091093e+03 4.497790187638189309e+03 4.032732114261302286e+03 3.596017012061617152e+03 3.248624535892444783e+03 3.008406897507815302e+03 2.840568351285121935e+03 2.679115183563475057e+03 2.467475917233809923e+03 2.211457587583752684e+03 1.967594266684561944e+03 1.782802333607400442e+03 1.654916989547497224e+03 1.536832153026960896e+03 1.394014403140221702e+03 1.243662299070582549e+03 1.123298141982568723e+03 1.031871056348076308e+03 9.407724101138215929e+02 8.439562202581603287e+02 7.604686472532837342e+02 6.919655962794691959e+02 6.249017486639138497e+02 5.613656417941500649e+02 5.070708876692686999e+02 4.572630267819340020e+02 4.107075477531604975e+02 3.698512755938050987e+02 3.323040146312443994e+02 2.983437452964743670e+02 2.678443607321100330e+02 2.401432118445162303e+02 2.152646812296342205e+02 1.927722853383984329e+02 1.725520566438248693e+02 1.543636601305944112e+02 1.379973311560080163e+02 1.232915168586189623e+02 1.100924612777798615e+02 9.825012395290012535e+01 8.763327515936649093e+01 7.812490159984299964e+01 6.960895200009956341e+01 6.198566342562281761e+01 5.517290733619083909e+01 4.908943606823985562e+01 4.365414556322568274e+01 3.879659393113936972e+01 3.446416357460656599e+01 3.060593258914536818e+01 2.717097907027387649e+01 2.410947413149482088e+01 2.138238350827838019e+01 1.895682530054453707e+01 1.679998914940483346e+01 1.487906675812324941e+01 1.316533028854911436e+01 1.164266696147453040e+01 1.029739432609936856e+01 9.115829931623688154e+00 8.084013755563756121e+00 7.180208793497410369e+00 6.374038211572703894e+00 diff --git a/skypy/power_spectrum/tests/data/linear_power.txt b/skypy/power_spectrum/tests/data/linear_power.txt new file mode 100644 index 000000000..cad67638d --- /dev/null +++ b/skypy/power_spectrum/tests/data/linear_power.txt @@ -0,0 +1,5 @@ +4.307566660102520473e+02 4.931290635794836135e+02 5.645220945652492901e+02 6.462349969961428542e+02 7.397514716658428142e+02 8.467645674632100281e+02 9.692042938494228110e+02 1.109267951517405209e+03 1.269452996079960030e+03 1.452591979990952950e+03 1.661888720275040669e+03 1.900954271762260532e+03 2.173840498851174743e+03 2.485067982837536874e+03 2.839643632832506682e+03 3.243061664426086281e+03 3.701279594371046187e+03 4.220658669616384032e+03 4.807855912881123913e+03 5.469653073200213839e+03 6.212706719389893806e+03 7.043204131930760923e+03 7.966412192744414824e+03 8.986111670954331203e+03 1.010391732880717973e+04 1.131849470281023241e+04 1.262469605332017272e+04 1.401264880793206976e+04 1.546683715782898980e+04 1.696521824105437190e+04 1.847840567402784654e+04 1.996893357193833799e+04 2.139058689781252178e+04 2.268776614249323393e+04 2.379489826507436010e+04 2.463613891687013165e+04 2.512628121139638097e+04 2.517521160885320933e+04 2.470034647761462111e+04 2.365234461194576579e+04 2.205369123692036737e+04 2.003243514681675151e+04 1.781569366515657020e+04 1.566440880465207556e+04 1.378743280726244848e+04 1.228848024387500482e+04 1.112325477066422718e+04 1.002361767305866124e+04 8.573117377860675333e+03 6.742857875273539321e+03 5.232419772215945159e+03 4.467473668696667119e+03 3.728585248113538000e+03 2.751789933317056239e+03 2.248542864176949934e+03 1.727670922198509970e+03 1.339943586669423667e+03 1.012396154052899192e+03 7.787483047651917332e+02 5.860367302088504857e+02 4.385847963767009787e+02 3.269127434685278786e+02 2.423868885696911377e+02 1.788338783323765142e+02 1.313709709107551475e+02 9.607377220583134658e+01 6.996901040698986662e+01 5.075620411608954896e+01 3.668202452203230024e+01 2.641590711328366581e+01 1.895848926214262065e+01 1.356240291332200876e+01 9.672398201459081690e+00 6.877990552690430803e+00 4.877308250036834991e+00 3.449463096551290242e+00 2.433504093772473897e+00 1.712684066745943756e+00 1.202650160190359863e+00 8.426899095778834781e-01 5.892640665071615658e-01 4.112543256277583215e-01 2.864914185112704681e-01 1.992293059008539424e-01 1.383161501478919753e-01 9.587514508345165720e-02 6.635688243758573746e-02 4.586105557697581214e-02 3.165262408188262294e-02 2.181778660666879727e-02 1.502010555344503000e-02 1.032811507310447442e-02 7.093783639760866885e-03 4.867043401835602703e-03 3.335824931142205262e-03 2.284079270565322808e-03 1.562453159326178658e-03 1.067843796784676176e-03 7.291709233128611852e-04 4.974929598788272334e-04 +2.566487270756076100e+02 2.938107670483471452e+02 3.363473821964690273e+02 3.850326703169853886e+02 4.407506337947823454e+02 5.045100065079157048e+02 5.774607055918372680e+02 6.609119027212760784e+02 7.563516045935871261e+02 8.654674716421956191e+02 9.901683670978836744e+02 1.132605790168111525e+03 1.295193878397878279e+03 1.480626035109163695e+03 1.691885422145670645e+03 1.932245542969373219e+03 2.205255940075554918e+03 2.514706702070842766e+03 2.864564143447303650e+03 3.258868891766866000e+03 3.701586991082202985e+03 4.196404879200796131e+03 4.746460612123860301e+03 5.354006806373603467e+03 6.020005551936877055e+03 6.743661763365764273e+03 7.521908370708660186e+03 8.348863205816769550e+03 9.215281809122470804e+03 1.010803084361985384e+04 1.100960163553931125e+04 1.189767167102497115e+04 1.274470561203827128e+04 1.351757676692895075e+04 1.417721612340544198e+04 1.467843493087597199e+04 1.497046615384278084e+04 1.499961932827630153e+04 1.471669037770708383e+04 1.409228136440496201e+04 1.313978918933174464e+04 1.193550648507623737e+04 1.061475181215785960e+04 9.332996787588223015e+03 8.214677470690985501e+03 7.321587943131705288e+03 6.627336041645168734e+03 5.972162289002957550e+03 5.107941062124727978e+03 4.017455856387148287e+03 3.117523140158444221e+03 2.661761316277589231e+03 2.221524896142763282e+03 1.639541391446476155e+03 1.339702224949945048e+03 1.029361999419959375e+03 7.983505375715925538e+02 6.031948075012107893e+02 4.639853005211925847e+02 3.491659971759394807e+02 2.613127981218888181e+02 1.947775765215795332e+02 1.444163058169547185e+02 1.065508461125153019e+02 7.827201554700847907e+01 5.724162453562973241e+01 4.168817077636705193e+01 3.024100659483251263e+01 2.185548278877186235e+01 1.573883695861455578e+01 1.129563979002863761e+01 8.080602619641789275e+00 5.762902543483304107e+00 4.097969130776124125e+00 2.905944475616126788e+00 2.055221387573610503e+00 1.449903802498835770e+00 1.020432695062599260e+00 7.165498693008858222e-01 5.020822883968446559e-01 3.510888734032030500e-01 2.450290558572332245e-01 1.706941846310342792e-01 1.187026267734973689e-01 8.241001630514574461e-02 5.712328069525785673e-02 3.953603218274234776e-02 2.732443271323491477e-02 1.885892040733864322e-02 1.299923507179229221e-02 8.949115068926606284e-03 6.153584600683980124e-03 4.226540608546418970e-03 2.899828585991315062e-03 1.987514738319446500e-03 1.360875138067905521e-03 9.309237583494098028e-04 6.362310157583220705e-04 4.344466471572545506e-04 2.964108160288772185e-04 +1.600492036921232284e+02 1.832238945350652557e+02 2.097502344853492389e+02 2.401109601511138862e+02 2.748573459508031078e+02 3.146184503479429395e+02 3.601113753673050724e+02 4.121525360615732438e+02 4.716698711340245040e+02 5.397158257362485756e+02 6.174807897194829138e+02 7.063064636205739362e+02 8.076983323724997490e+02 9.233360343736912910e+02 1.055079904888654482e+03 1.204971339673925058e+03 1.375223875715696749e+03 1.568201057420852976e+03 1.786376325757837776e+03 2.032269464210761271e+03 2.308353745098840591e+03 2.616928074955806551e+03 2.959949382888062701e+03 3.338822427394555689e+03 3.754147179252581282e+03 4.205427813704868640e+03 4.690751669391216637e+03 5.206450556179043815e+03 5.746759518951603241e+03 6.303488452292952388e+03 6.865718738657385984e+03 7.419529792473534144e+03 7.947750248908690992e+03 8.429721908407413139e+03 8.841080869470737525e+03 9.153646888890509217e+03 9.335761038536400520e+03 9.353941305418487900e+03 9.177503441276223384e+03 8.788114541919834664e+03 8.194129074388030858e+03 7.443124033246685940e+03 6.619485684902320827e+03 5.820167981875827536e+03 5.122770733183479933e+03 4.565829464312396340e+03 4.132885707837310292e+03 3.724311550524327686e+03 3.185372897869555345e+03 2.505333332487551161e+03 1.944124569638468529e+03 1.659906066719841419e+03 1.385369390533277510e+03 1.022437543763776944e+03 8.354542675157618987e+02 6.419224057540799322e+02 4.978609060775501689e+02 3.761594678915037662e+02 2.893467609188279539e+02 2.177440754962799758e+02 1.629577739602347890e+02 1.214656171280236805e+02 9.005972875662283172e+01 6.644637698906717560e+01 4.881136135923902231e+01 3.569655898691486584e+01 2.599723993243574327e+01 1.885865197735943966e+01 1.362933943412608251e+01 9.814926226087749228e+00 7.044095539404557904e+00 5.039160058819049937e+00 3.593814680281249174e+00 2.555542369561032867e+00 1.812181593867290053e+00 1.281660541395452046e+00 9.041772841201642930e-01 6.363539851809317804e-01 4.468490348425400716e-01 3.131044964121753593e-01 2.189432040190619488e-01 1.528030382937760889e-01 1.064469270366772552e-01 7.402437217490953625e-02 5.139186792852414726e-02 3.562275835821427511e-02 2.465514066676026406e-02 1.703984175929268777e-02 1.176064743464902244e-02 8.106477852251981614e-03 5.580774768888429639e-03 3.837448665394871828e-03 2.635721074786377769e-03 1.808367652237913813e-03 1.239437868322936041e-03 8.486579483716337772e-04 5.805351459156417359e-04 3.967612409250406377e-04 2.709261047834744872e-04 1.848453160540263168e-04 +1.068609742737917543e+02 1.223341536763868476e+02 1.400451479557332846e+02 1.603162686452428431e+02 1.835155883131517385e+02 2.100631140493851774e+02 2.404376374958448253e+02 2.751842591961697622e+02 3.149225413325579552e+02 3.603551760251460792e+02 4.122769577267299610e+02 4.715837073676458431e+02 5.392806007499679026e+02 6.164890917237838721e+02 7.044512810572012995e+02 8.045301592205846646e+02 9.182036512112415494e+02 1.047049838345619719e+03 1.192720177210814427e+03 1.356897066168450237e+03 1.541231849202621106e+03 1.747259450488466200e+03 1.976286464161232288e+03 2.229250813424246189e+03 2.506553084224320628e+03 2.807862226388079307e+03 3.131901202281519090e+03 3.476220850257182519e+03 3.836972049505911400e+03 4.208686465140915971e+03 4.584074000855929626e+03 4.953840219051125132e+03 5.306520215601372001e+03 5.628321011345814441e+03 5.902975419749129287e+03 6.111668175411502489e+03 6.233261379321712411e+03 6.245399902892967475e+03 6.127596617240273190e+03 5.867611086561428237e+03 5.471021386015196185e+03 4.969593521773701468e+03 4.419670158689408709e+03 3.885985101036702872e+03 3.420349860548506058e+03 3.048493673626566306e+03 2.759427620466566623e+03 2.486632558033478745e+03 2.126796281638884011e+03 1.672750345607630607e+03 1.298044856385669618e+03 1.108279050446718884e+03 9.249775655632391818e+02 6.826567676704844416e+02 5.578125659385904100e+02 4.285960323739692512e+02 3.324096605855094140e+02 2.511525599247667628e+02 1.931898195146944204e+02 1.453824418560583354e+02 1.088029561482662899e+02 8.109964865577636317e+01 6.013069815879234881e+01 4.436463548843527605e+01 3.259016296333562934e+01 2.383372727615346776e+01 1.735772702095818332e+01 1.259146485769225343e+01 9.099979612773301696e+00 6.553188361765932335e+00 4.703171867486696023e+00 3.364525039704918452e+00 2.399502960558685061e+00 1.706273702770608081e+00 1.209949729298290144e+00 8.557336805327313556e-01 6.036972585204244401e-01 4.248781329162656384e-01 2.983502705107641262e-01 2.090522836994074263e-01 1.461830709080495783e-01 1.020228852589538660e-01 7.107203328216479821e-02 4.942427920999332980e-02 3.431310465721735398e-02 2.378445238503851603e-02 1.646163986904641136e-02 1.137708935667058885e-02 7.852299255263430269e-03 5.412498789353206892e-03 3.726148054776998782e-03 2.562171467585624540e-03 1.759807080998925617e-03 1.207403252908926129e-03 8.275426250517417147e-04 5.666283436351345721e-04 3.876092467917505234e-04 2.649078644645270781e-04 1.808907938652263073e-04 1.234167375270255073e-04 +7.569752992301025074e+01 8.665832705958763427e+01 9.920433394880144817e+01 1.135638676831241298e+02 1.299976612797983933e+02 1.488032368180824108e+02 1.703197578222284960e+02 1.949333593148076602e+02 2.230829230030551003e+02 2.552662176762990782e+02 2.920462550166060964e+02 3.340576112304971161e+02 3.820123266664012931e+02 4.367048100123470817e+02 4.990149330896443303e+02 5.699082028349117763e+02 6.504315428090467321e+02 7.417028246998615941e+02 8.448918973251101079e+02 9.611905278494376716e+02 1.091768485325565280e+03 1.237713070046858093e+03 1.399949839255538336e+03 1.579143193311343111e+03 1.775576896861765590e+03 1.989016442589466578e+03 2.218557210307613332e+03 2.462464277717820096e+03 2.718011027927357645e+03 2.981323834979370531e+03 3.247238584593266296e+03 3.509171339339017777e+03 3.759000659851371438e+03 3.986955959067539425e+03 4.181513985885901093e+03 4.329346496527181444e+03 4.415479954077784896e+03 4.424078567907582510e+03 4.340629789704643372e+03 4.156462813670078503e+03 3.875529002910140207e+03 3.520330568536883675e+03 3.130779186329009008e+03 2.752730596601364141e+03 2.422886724321528618e+03 2.159473490183688682e+03 1.954706629713728262e+03 1.761465715135592745e+03 1.506567072437695515e+03 1.184932761476930864e+03 9.195011558280567669e+02 7.850760032309696044e+02 6.552300072423062147e+02 4.835762676604721264e+02 3.951398879574271064e+02 3.036062623047691318e+02 2.354703426566132691e+02 1.779099296945897493e+02 1.368506346017953490e+02 1.029851338851056397e+02 7.707317928473455027e+01 5.744887806410242348e+01 4.259501987605424489e+01 3.142675186364000695e+01 2.308602230962368296e+01 1.688319141692470282e+01 1.229576156771645579e+01 8.919465635767227951e+00 6.446188458583363534e+00 4.642107892774103739e+00 3.331604410230888380e+00 2.383341875746340488e+00 1.699745378437842680e+00 1.208679834243294504e+00 8.570968631095334800e-01 6.061794432295478652e-01 4.276434086610874097e-01 3.009725991984494309e-01 2.113435581418399156e-01 1.480871909352470261e-01 1.035522786452222482e-01 7.227035372084925080e-02 5.034557660200295742e-02 3.501087164745471381e-02 2.430650931448716939e-02 1.684828636791148651e-02 1.166099677676786170e-02 8.059233671281917921e-03 5.562364201516683745e-03 3.834073120235489008e-03 2.639506197570061670e-03 1.814975510503212758e-03 1.246601484573145241e-03 8.552930055835086158e-04 5.862096340420900687e-04 4.013847551834076881e-04 2.745722913051341484e-04 1.876538290374015004e-04 1.281383252816582969e-04 8.742520125276057475e-05 diff --git a/skypy/power_spectrum/tests/data/truth_bird.txt b/skypy/power_spectrum/tests/data/truth_bird.txt new file mode 100644 index 000000000..f471a795d --- /dev/null +++ b/skypy/power_spectrum/tests/data/truth_bird.txt @@ -0,0 +1,5 @@ +4.307269013457246842e+02 4.930899095542520172e+02 5.644705908182455687e+02 6.461672506949836361e+02 7.396623643672360231e+02 8.466473702082029149e+02 9.690501620345186211e+02 1.109065262440231209e+03 1.269186480915614084e+03 1.452241587671554498e+03 1.661428132569519903e+03 1.900348965943235271e+03 2.173045225806472899e+03 2.484023493511534980e+03 2.838272443883833148e+03 3.241262604912159077e+03 3.698920820032097254e+03 4.217568772419810557e+03 4.803812695735952730e+03 5.464369511237957340e+03 6.205813602138261558e+03 7.034228892064775209e+03 7.954753498732950902e+03 8.971009518700651824e+03 1.008441863920382275e+04 1.129341484491272786e+04 1.259257728963055888e+04 1.397171741819885392e+04 1.541496196479661194e+04 1.689987449346785797e+04 1.839664989941351814e+04 1.986739749739583203e+04 2.126550274199267733e+04 2.253504322580065127e+04 2.361028473725260847e+04 2.441552547880502607e+04 2.486623482817330296e+04 2.487384749385593022e+04 2.435859293537778285e+04 2.327540438318622546e+04 2.165197393134290178e+04 1.962073978624396477e+04 1.740985921237814546e+04 1.527635727521488116e+04 1.342066657546548049e+04 1.193544634075051363e+04 1.076823790719756289e+04 9.663874283269122316e+03 8.260333090039557646e+03 6.567593631768332671e+03 5.196248783916698812e+03 4.465262712763136733e+03 3.771825174973546382e+03 2.940839819393791004e+03 2.477118534366473341e+03 2.028957449640702407e+03 1.690079789383751859e+03 1.410148845149354429e+03 1.199399958480103578e+03 1.025477489343606067e+03 8.854781345671450481e+02 7.708918638170872555e+02 6.749430720288331713e+02 5.928357729227226400e+02 5.210829556534134781e+02 4.572146441464641953e+02 3.995665866901977097e+02 3.470623371602093243e+02 2.990602734437238723e+02 2.552307844480522760e+02 2.154544821452882388e+02 1.797288756738040547e+02 1.480802602322040968e+02 1.204892768208950002e+02 9.684426413622729513e+01 7.692843971844459361e+01 6.043374424126144362e+01 4.698885891392644965e+01 3.619130953744849677e+01 2.763692728960397105e+01 2.094270840202186434e+01 1.576181945925353745e+01 1.179144917372067525e+01 8.775065012182748703e+00 6.500762336011747244e+00 4.797203433029175734e+00 3.528341362750625265e+00 2.587798156287789197e+00 1.893467228177847961e+00 1.382664253624534823e+00 1.007969648484281722e+00 7.337855196135111280e-01 5.335591340181766729e-01 3.875900097483609974e-01 2.813263320321804928e-01 2.040589100085554164e-01 1.479305405369756576e-01 1.071910855231500148e-01 7.764124920212671366e-02 5.621963589264222100e-02 +2.566372050509418727e+02 2.937956056329187504e+02 3.363274322722764964e+02 3.850064202318699813e+02 4.407160951923852963e+02 5.044645643480864692e+02 5.774009212317585025e+02 6.608332555390960579e+02 7.562481532642976845e+02 8.653314105663698683e+02 9.899894463399404003e+02 1.132370558246708697e+03 1.294884696033375121e+03 1.480219794240075544e+03 1.691351887938967593e+03 1.931545217139445867e+03 2.204337321951985814e+03 2.513502800514278988e+03 2.862988069880992498e+03 3.256808345636588001e+03 3.698897439208063133e+03 4.192901237298452543e+03 4.741907270532071379e+03 5.348105915979500423e+03 6.012383616720258942e+03 6.733854792146437831e+03 7.509346144892929260e+03 8.332853648301292196e+03 9.194997024154421524e+03 1.008249718161225610e+04 1.097769718486012607e+04 1.185813599158647048e+04 1.269616845782909695e+04 1.345861519558747386e+04 1.410645287036274385e+04 1.459469472687603047e+04 1.487300817929185723e+04 1.488845542425277199e+04 1.459296348951615801e+04 1.395859700686417273e+04 1.300027678509019097e+04 1.179534941385216553e+04 1.047915037160502106e+04 9.205964332040335648e+03 8.098213693630070338e+03 7.215480057051692711e+03 6.530223487369859868e+03 5.885948990063721794e+03 5.044945988897807183e+03 3.996285669519063049e+03 3.138532618495913539e+03 2.703322226465698350e+03 2.284516461570183765e+03 1.742566985655430926e+03 1.458781625614147742e+03 1.169413988777075701e+03 9.522511019116781199e+02 7.692036685063405912e+02 6.349400277239973320e+02 5.232468096408149449e+02 4.351050906456204075e+02 3.652323721913830354e+02 3.090758249523728409e+02 2.633838347258957242e+02 2.256570314627952598e+02 1.939819941788163931e+02 1.669623193643859906e+02 1.435775361937998298e+02 1.230982672953774539e+02 1.050122907709820765e+02 8.896734536890564016e+01 7.472574349853043429e+01 6.212947771672407526e+01 5.107206029081356036e+01 4.147408580838982317e+01 3.326128359702440918e+01 2.634720144665046604e+01 2.062447142473543238e+01 1.596609704321752687e+01 1.223350221668543902e+01 9.286250008171895942e+00 6.990175125066887318e+00 5.222981950427214137e+00 3.877434660403003441e+00 2.862550559093461455e+00 2.103281537441191062e+00 1.539187485111150000e+00 1.122562906016414619e+00 8.163806260104576884e-01 5.922989298620008336e-01 4.288724108388075051e-01 3.100271948296682289e-01 2.238093360641297547e-01 1.613858595998066692e-01 1.162643948373295461e-01 8.369368341424217317e-02 6.020907527345398247e-02 4.329154661049863101e-02 3.111402825430063107e-02 2.235396226917294241e-02 +1.600445443906299090e+02 1.832177624840515477e+02 2.097421643284513095e+02 2.401003396046942271e+02 2.748433694966484495e+02 3.146000583896338867e+02 3.600871743282839361e+02 4.121206935247508341e+02 4.716279783439352968e+02 5.396607175017364852e+02 6.174083088458933162e+02 7.062111533280722142e+02 8.075730355762785848e+02 9.231713730757762733e+02 1.054863606160904965e+03 1.204687366784980668e+03 1.374851314019202846e+03 1.567712695095607842e+03 1.785736860207064865e+03 2.031433257330329070e+03 2.307262046995361743e+03 2.615505639420132411e+03 2.958100413777570338e+03 3.336425829496752158e+03 3.751051148628974715e+03 4.201443901728283890e+03 4.685648605260679687e+03 5.199948331523130037e+03 5.738524545952745029e+03 6.293131028155137756e+03 6.852794610808294237e+03 7.403548342518211939e+03 7.928192696512293878e+03 8.406073899882243495e+03 8.812884574901814631e+03 9.120577506390171038e+03 9.297729594801439816e+03 9.311216494063448408e+03 9.130825943432868371e+03 8.738749174250617216e+03 8.143783249262372010e+03 7.393698588420975284e+03 6.572736103625030410e+03 5.777456432525072159e+03 5.085049674463397423e+03 4.533870682606760056e+03 4.107751593626251633e+03 3.707847525826492983e+03 3.179709040644539982e+03 2.511587581916398904e+03 1.962256803655259546e+03 1.689280713546306515e+03 1.425586879613462315e+03 1.072238293025411167e+03 8.935138705849365124e+02 7.065885390313513881e+02 5.674122888755709937e+02 4.489611812225294898e+02 3.637364292755124779e+02 2.923767303288539665e+02 2.366147856253086275e+02 1.931091021424126382e+02 1.588462794702028589e+02 1.317047622042278476e+02 1.100269735689040971e+02 9.250234001102158743e+01 7.815578644855651191e+01 6.624888978176419130e+01 5.623262680172326355e+01 4.770252659771870185e+01 4.036484738213990653e+01 3.400791223610865188e+01 2.847978148530639686e+01 2.367044702482261798e+01 1.949821634567302553e+01 1.589972877787024963e+01 1.282292020930636767e+01 1.022194970536965819e+01 8.053200643496714051e+00 6.272260545448747671e+00 4.832578681816528388e+00 3.686255867226801541e+00 2.786253601591068829e+00 2.088676122513974942e+00 1.554291387860955664e+00 1.149219348942160934e+00 8.450248268484529257e-01 6.184327778808454967e-01 4.508121891336519527e-01 3.275378560290994057e-01 2.373187605403186684e-01 1.715583980231222549e-01 1.237863571894008280e-01 8.917794899880705783e-02 6.416302411290354668e-02 4.611616172058426166e-02 3.311644657928903418e-02 2.376422247196396248e-02 1.704307884652621688e-02 1.221694350330704976e-02 +1.068589186525389323e+02 1.223314481262159887e+02 1.400415870596666537e+02 1.603115821133931718e+02 1.835094205317375042e+02 2.100549972085824777e+02 2.404269562985884079e+02 2.751702045410015103e+02 3.149040495775216755e+02 3.603308494307695469e+02 4.122449603904642004e+02 4.715416293390209148e+02 5.392252809910291944e+02 6.164163880928063008e+02 7.043557725019380769e+02 8.044047618659324144e+02 9.180391260833532669e+02 1.046834164280719961e+03 1.192437758434088892e+03 1.356527740760672714e+03 1.540749663260740590e+03 1.746631164826466829e+03 1.975469768178839104e+03 2.228192237462644698e+03 2.505185631597579686e+03 2.806102794957200331e+03 3.129647940231665416e+03 3.473350700936650355e+03 3.833338897162211651e+03 4.204120509989336824e+03 4.578383279968901661e+03 4.946815553053824260e+03 5.297945279809251588e+03 5.617989622032812804e+03 5.890718269356811106e+03 6.097390631855999345e+03 6.216991570065852102e+03 6.227340590201335544e+03 6.108164330571360551e+03 5.847434943998574454e+03 5.450875265940903773e+03 4.950267251363956348e+03 4.401841015991898530e+03 3.870174794662428212e+03 3.407006406911242266e+03 3.038179881879046661e+03 2.753039836680869030e+03 2.485234047650671528e+03 2.130371433070030889e+03 1.679719373012314918e+03 1.308204710525487144e+03 1.123809432986656475e+03 9.455095617259846676e+02 7.048852048445424998e+02 5.840099203465522351e+02 4.569945990888368783e+02 3.628005742294832885e+02 2.826056365581287650e+02 2.255030917469377698e+02 1.777726156525733643e+02 1.407917593998938344e+02 1.122804390005792214e+02 9.015202284365858532e+01 7.294670747710922853e+01 5.952052646115501489e+01 4.895912976027293695e+01 4.057879347212503518e+01 3.385500865088642541e+01 2.839192878973010892e+01 2.389233961435146014e+01 2.013609333344681573e+01 1.696159941662111237e+01 1.425177420199413625e+01 1.192262152944657494e+01 9.914077758655523098e+00 8.182627753324782205e+00 6.695508158014927602e+00 5.426450687485676383e+00 4.352848395873063936e+00 3.454122160078960224e+00 2.710889819763047104e+00 2.104521166381269293e+00 1.616873882018566411e+00 1.230285159011821117e+00 9.279174278483200977e-01 6.943053719891815190e-01 5.158032599984636590e-01 3.807723487546819108e-01 2.795448463487534285e-01 2.042629483474907026e-01 1.486628057895998112e-01 1.078392854794387457e-01 7.801202190343449472e-02 5.630738044088808708e-02 4.056606942237907937e-02 2.918090436874877514e-02 2.096484709447135608e-02 1.504663881063740376e-02 1.079002170753740920e-02 7.732250577483334129e-03 +7.569654554520784018e+01 8.665703142299147999e+01 9.920262866952434422e+01 1.135616233049270392e+02 1.299947074784229244e+02 1.487993495300450491e+02 1.703146423412462696e+02 1.949266281074130234e+02 2.230740666185803320e+02 2.552545666433662177e+02 2.920309300027826112e+02 3.340374579648135409e+02 3.819858311552784471e+02 4.366699883892551952e+02 4.989691890066966948e+02 5.698481436281476817e+02 6.503527437803829798e+02 7.415995288474113067e+02 8.447566358997962652e+02 9.610136455321628546e+02 1.091537553103515620e+03 1.237412169765083718e+03 1.399558707394340217e+03 1.578636219252471847e+03 1.774921984707324555e+03 1.988173768034433124e+03 2.217477938413597258e+03 2.461089357820826990e+03 2.716270253373651485e+03 2.979135453639252319e+03 3.244509876566112325e+03 3.505800743260301260e+03 3.754882221113598916e+03 3.981987058556846478e+03 4.175607392252891259e+03 4.322447545207881376e+03 4.407588637689520510e+03 4.415273802255197552e+03 4.331088535249802590e+03 4.146461355776842538e+03 3.865413324015332819e+03 3.510456265449885450e+03 3.121447252850642144e+03 2.744159312352199777e+03 2.415241394620433312e+03 2.152948451643580938e+03 1.949631091664078895e+03 1.758267355012300413e+03 1.505284372814509652e+03 1.184971322640226845e+03 9.207559942835440552e+02 7.883868321718200605e+02 6.605016427043894964e+02 4.895449236646347799e+02 4.027971073758608895e+02 3.123285711312957460e+02 2.452652382576368382e+02 1.884701795340540968e+02 1.481485424898141616e+02 1.147177676775025645e+02 8.903896084745807116e+01 6.945575653171330544e+01 5.446072647654650467e+01 4.299294637500825900e+01 3.422172045997492518e+01 2.748441671661727526e+01 2.228480538445809245e+01 1.824176197503834373e+01 1.506672133534975977e+01 1.254189181682893484e+01 1.050475166946161210e+01 8.834914366482303549e+00 7.444345116358795700e+00 6.269572851132196867e+00 5.265582628650760988e+00 4.400928027054657576e+00 3.653793195195512045e+00 3.008861155689652289e+00 2.454866823349546756e+00 1.982791706105456520e+00 1.584657531211403514e+00 1.252853042915384130e+00 9.798754621146523736e-01 7.583249627345184374e-01 5.810007025587953233e-01 4.410184655114193220e-01 3.319453540496589694e-01 2.479537214738308237e-01 1.839516195976669821e-01 1.356379967262777986e-01 9.947598587554040694e-02 7.261564961220443037e-02 5.279876712406497696e-02 3.826322277083853141e-02 2.765365904604897423e-02 1.994110947691092942e-02 1.435320562548800596e-02 1.031566935668742739e-02 7.404817951845461127e-03 5.310034910426884858e-03 diff --git a/skypy/power_spectrum/tests/data/truth_smith.txt b/skypy/power_spectrum/tests/data/truth_smith.txt new file mode 100644 index 000000000..7eb1ada2f --- /dev/null +++ b/skypy/power_spectrum/tests/data/truth_smith.txt @@ -0,0 +1,5 @@ +4.307268924531355196e+02 4.930898975252965215e+02 5.644705745323444717e+02 6.461672286240507219e+02 7.396623344240108509e+02 8.466473295364911564e+02 9.690501067183487294e+02 1.109065187099530249e+03 1.269186378142639114e+03 1.452241447243328366e+03 1.661427940344499348e+03 1.900348702312104933e+03 2.173044863512668144e+03 2.484022994580630893e+03 2.838271755294170362e+03 3.241261652491735731e+03 3.698919499867241484e+03 4.217566938798318006e+03 4.803810144315382786e+03 5.464365955738852790e+03 6.205808642458262511e+03 7.034221971411966479e+03 7.954743847440484387e+03 8.970996083303058185e+03 1.008439999801260819e+04 1.129338911721163640e+04 1.259254205705160712e+04 1.397166969722000249e+04 1.541489830298934248e+04 1.689979131369342213e+04 1.839654428798639492e+04 1.986726871873374694e+04 2.126535483812315579e+04 2.253488905730105034e+04 2.361015159935536576e+04 2.441546172541224223e+04 2.486631541117378219e+04 2.487417383802523545e+04 2.435928017612599069e+04 2.327655061839192058e+04 2.165361805442990226e+04 1.962283556616375063e+04 1.741230059613273625e+04 1.527905863525586574e+04 1.342366198447084207e+04 1.193896959892624000e+04 1.077269347091320378e+04 9.669380723427724661e+03 8.265749528105183344e+03 6.570931177289904554e+03 5.197453799547990457e+03 4.466210302796373981e+03 3.772152735520938222e+03 2.938908440860989231e+03 2.474669128531984825e+03 2.025642891775987891e+03 1.686382479470304588e+03 1.406286320224776546e+03 1.195782657322279420e+03 1.022270971493611114e+03 8.828291567488212195e+02 7.688408083713314909e+02 6.734167349371214186e+02 5.916586916884594984e+02 5.199976980896733494e+02 4.559135133833511873e+02 3.977322815435029497e+02 3.444118677656650789e+02 2.953861350555572471e+02 2.504326543522320776e+02 2.095555439216587104e+02 1.728730507604634283e+02 1.405102565461263850e+02 1.125098265927821473e+02 8.877878963069832707e+01 6.907873568814957821e+01 5.305007979037982579e+01 4.025348430063843352e+01 3.021446068580046074e+01 2.246248982072604150e+01 1.656015298985475681e+01 1.212104196316129467e+01 8.817609251253863079e+00 6.381396580455243495e+00 4.598359600549411397e+00 3.301649435334715221e+00 2.363586489249736466e+00 1.687927183687062849e+00 1.203009738511493065e+00 8.560059531933452703e-01 6.082859937534358385e-01 4.317874232116069577e-01 3.062326423160630084e-01 2.170322408451179630e-01 1.537258462344259224e-01 1.088346057088414032e-01 7.702349878025624430e-02 5.449377824405470666e-02 3.854458849543985349e-02 2.725795436662901972e-02 +2.566372070647824444e+02 2.937956082549575285e+02 3.363274356739386803e+02 3.850064246260228629e+02 4.407161008392331496e+02 5.044645715590253303e+02 5.774009303684969154e+02 6.608332670034329794e+02 7.562481674709129038e+02 8.653314278862086439e+02 9.899894669943378176e+02 1.132370582120972813e+03 1.294884722361138529e+03 1.480219821083515853e+03 1.691351911350992168e+03 1.931545229836732233e+03 2.204337311245250021e+03 2.513502745125577349e+03 2.862987935117602774e+03 3.256808076178578176e+03 3.698896948499309474e+03 4.192900392220923095e+03 4.741905869734218868e+03 5.348103660024029523e+03 6.012380068318095255e+03 6.733849325518218393e+03 7.509337885350612851e+03 8.332841407389925735e+03 9.194979243404755834e+03 1.008247191063525679e+04 1.097766213846122446e+04 1.185808875572474244e+04 1.269610694168860209e+04 1.345853844428342018e+04 1.410636233094699310e+04 1.459459594429773097e+04 1.487291260787819738e+04 1.488838151935294627e+04 1.459293565766150823e+04 1.395864108315190788e+04 1.300041277017523316e+04 1.179558622742311491e+04 1.047948904778120777e+04 9.206410480489948895e+03 8.098792878366159130e+03 7.216250232059411246e+03 6.531271458645319399e+03 5.887313135566280835e+03 5.046433552770399729e+03 3.997547216703016602e+03 3.139561900391413019e+03 2.704380832301299506e+03 2.285511303808846606e+03 1.743271670977121630e+03 1.459394622274818175e+03 1.169929166911697166e+03 9.527472197269701155e+02 7.697754348222712224e+02 6.356608109143477350e+02 5.242121203221487349e+02 4.363778885378634982e+02 3.668417746739735321e+02 3.110104485809366111e+02 2.655837846866302812e+02 2.280131925622416986e+02 1.963426923006732352e+02 1.691463525221310249e+02 1.453936467006655846e+02 1.243682539025089540e+02 1.055945966000173399e+02 8.877730754612444741e+01 7.374868082728224294e+01 6.042350662820458496e+01 4.875937518261158488e+01 3.872216598544438426e+01 3.025816793225992996e+01 2.327636888915649038e+01 1.764434703487938094e+01 1.319722319072422501e+01 9.754045565665814266e+00 7.134592581918813181e+00 5.172192912078578786e+00 3.721322948950293696e+00 2.660571015263709427e+00 1.892255194133234397e+00 1.340042007465594542e+00 9.456596298715985727e-01 6.654509581220939829e-01 4.671959996503354251e-01 3.273996768213129349e-01 2.290928403394646151e-01 1.601131243797833392e-01 1.117963823270008705e-01 7.800043332819341091e-02 5.438778394891596868e-02 3.790479668540509134e-02 2.640694362177795768e-02 1.839108656400300496e-02 1.280528281560912307e-02 8.914259017347641451e-03 +1.600445454886789491e+02 1.832177639030060448e+02 2.097421661538284638e+02 2.401003419400704502e+02 2.748433724645088887e+02 3.146000621299231170e+02 3.600871789925964208e+02 4.121206992627674595e+02 4.716279852766159024e+02 5.396607256725123989e+02 6.174083181360579147e+02 7.062111633148324472e+02 8.075730453001588103e+02 9.231713806663123023e+02 1.054863608243006183e+03 1.204687357548991486e+03 1.374851283780088579e+03 1.567712628083412710e+03 1.785736731212292170e+03 2.031433026662544535e+03 2.307261652987914658e+03 2.615504987415869891e+03 2.958099360151757537e+03 3.336424158818431351e+03 3.751048541559980549e+03 4.201439890784090494e+03 4.685642515255811304e+03 5.199939201651379335e+03 5.738511031411903787e+03 6.293111281599033646e+03 6.852766148448641616e+03 7.403507907234705272e+03 7.928136141482349558e+03 8.405996130777855797e+03 8.812779604219087560e+03 9.120438683067710372e+03 9.297550061366824593e+03 9.310989869241813722e+03 9.130547112630281845e+03 8.738414923747541252e+03 8.143392422562437787e+03 7.393251582426176356e+03 6.572233696912924643e+03 5.776897900315601873e+03 5.084430169154232317e+03 4.533178560788861432e+03 4.106969054172840515e+03 3.706961866080707296e+03 3.178737492817685506e+03 2.510596357717497085e+03 1.961277083060902214e+03 1.688187881683597425e+03 1.424385605999264271e+03 1.071146764640740230e+03 8.923375464071239094e+02 7.054619066893744730e+02 5.663252623211387800e+02 4.479988244488304190e+02 3.628761365821521849e+02 2.917039754964362714e+02 2.361709650221392849e+02 1.929258272383900135e+02 1.589460242856839614e+02 1.320907989079396145e+02 1.106787982124160834e+02 9.337498505863452181e+01 7.918060640330755007e+01 6.733876416561122369e+01 5.728985466910523172e+01 4.862868149996033651e+01 4.107133534394322538e+01 3.442608262145343900e+01 2.856916059923916862e+01 2.342358364769205181e+01 1.894083763650656138e+01 1.508572283763530741e+01 1.182496820578950825e+01 9.119917918658297040e+00 6.923010817939688089e+00 5.177532304073872105e+00 3.820131567781743787e+00 2.785263161702906221e+00 2.010058906742148377e+00 1.438114244220764171e+00 1.021516989930042918e+00 7.213014528012988702e-01 5.068520600950924271e-01 3.547633630385326198e-01 2.475264780853942415e-01 1.722672032410733634e-01 1.196474241551674922e-01 8.296654149057455796e-02 5.745731268654837021e-02 3.975074996509905118e-02 2.747860719089641537e-02 1.898306731580501686e-02 1.310745045369341409e-02 9.046823515999731502e-03 6.242181881462068559e-03 4.305939266274396576e-03 +1.068589188878538749e+02 1.223314484124126409e+02 1.400415874001882628e+02 1.603115825061008195e+02 1.835094209636677078e+02 2.100549976472497633e+02 2.404269566778199874e+02 2.751702047381854754e+02 3.149040493774637639e+02 3.603308484687854047e+02 4.122449580634386166e+02 4.715416246672880334e+02 5.392252724050856614e+02 6.164163731056534061e+02 7.043557472081972719e+02 8.044047201798651940e+02 9.180390585854337360e+02 1.046834056480774507e+03 1.192437588170634854e+03 1.356527474331028998e+03 1.540749249702308362e+03 1.746630527531652660e+03 1.975468792708408728e+03 2.228190754010133787e+03 2.505183390041592247e+03 2.806099429823723312e+03 3.129642922206254298e+03 3.473343270895439218e+03 3.833327978188506677e+03 4.204104592793131815e+03 4.578360278012604795e+03 4.946782626300369884e+03 5.297898631614168153e+03 5.617924284254152553e+03 5.890627909040604209e+03 6.097267441977914132e+03 6.216826347229060957e+03 6.227123142119775366e+03 6.107884350493817692e+03 5.847083336827551648e+03 5.450445496112046385e+03 4.949755702844156986e+03 4.401245211635934538e+03 3.869489260814816589e+03 3.406217264347990749e+03 3.037259836327594257e+03 2.751949166685053115e+03 2.483945655575433193e+03 2.128930625051486004e+03 1.678252127553695345e+03 1.306741056368453656e+03 1.122176256201643810e+03 9.437320696721025115e+02 7.032337126683634096e+02 5.822596297134774659e+02 4.552918270922349961e+02 3.611283546374926345e+02 2.810376970202713665e+02 2.240052880611867749e+02 1.764130995204149031e+02 1.395980508657663677e+02 1.112753819475748429e+02 8.935720349086699343e+01 7.237526338827024119e+01 5.917406451819297075e+01 4.882703002251533775e+01 4.063717196442842550e+01 3.406832230293183983e+01 2.871574675831018553e+01 2.427702740314572338e+01 2.053077514630150802e+01 1.731819652545134858e+01 1.452889839044627784e+01 1.208912100010555157e+01 9.951830780880669636e+00 8.088012260198567560e+00 6.479003691266574272e+00 5.110166876342798759e+00 3.966349446692774894e+00 3.029463047482390259e+00 2.278010983876620710e+00 1.687939408797406227e+00 1.234079979208950961e+00 8.916904573179076365e-01 6.378379042000109633e-01 4.524207196726385272e-01 3.186771611514723213e-01 2.231989529683754969e-01 1.556102158772473543e-01 1.080890842578590694e-01 7.485934696730291160e-02 5.172407378999929056e-02 3.567244573924009676e-02 2.456604917794153950e-02 1.689796459832801356e-02 1.161279382077663215e-02 7.974924767301149256e-03 5.473569032761559268e-03 3.755097378574730367e-03 2.575251966350863668e-03 +7.569654554222167064e+01 8.665703140076294630e+01 9.920262861245024055e+01 1.135616231875471698e+02 1.299947072597673525e+02 1.487993491448931138e+02 1.703146416866574100e+02 1.949266270220325907e+02 2.230740648509327855e+02 2.552545638033934665e+02 2.920309254879463197e+02 3.340374508476671167e+02 3.819858200128091994e+02 4.366699710444571565e+02 4.989691621382403923e+02 5.698481021824939035e+02 6.503526800891851281e+02 7.415994313072965269e+02 8.447564870058025690e+02 9.610134189625468935e+02 1.091537209422325077e+03 1.237411650127502071e+03 1.399557924393697476e+03 1.578635043721100374e+03 1.774920226860842604e+03 1.988171150820222920e+03 2.217474060291238402e+03 2.461083641493130017e+03 2.716261876279684657e+03 2.979123255485607842e+03 3.244492239146854445e+03 3.505775438935838338e+03 3.754846229744550783e+03 3.981936360456334569e+03 4.175536760149389011e+03 4.322350385863453084e+03 4.407456977112778986e+03 4.415098560019019715e+03 4.330860250959804944e+03 4.146171405217351094e+03 3.865055252057541111e+03 3.510026244985114317e+03 3.120942499177726404e+03 2.743574247525914870e+03 2.414562673628659468e+03 2.152150430956434775e+03 1.948676758735185331e+03 1.757131469621372162e+03 1.504008597322270816e+03 1.183667913587031308e+03 9.194472054789081312e+02 7.869241121442331632e+02 6.589099940404860263e+02 4.880480025843961585e+02 4.012149126720584604e+02 3.107803446117438853e+02 2.437404219017894889e+02 1.870263876458986942e+02 1.467616270460216015e+02 1.134377959578459638e+02 8.788478427270899829e+01 6.844225891513558224e+01 5.360182254377985345e+01 4.229714418071976922e+01 3.369061053189702903e+01 2.711270039021696476e+01 2.205990053825261654e+01 1.814497057429624149e+01 1.507487471569938897e+01 1.262925862195758420e+01 1.064476314913351018e+01 9.001623920025169312e+00 7.613624062187170694e+00 6.420184437960062773e+00 5.380239579486063128e+00 4.467322605073089647e+00 3.665431125293586945e+00 2.965418340465647340e+00 2.361806499356021050e+00 1.850151783932875826e+00 1.425187108895539900e+00 1.079909946057611148e+00 8.055682003949958814e-01 5.922611037834141223e-01 4.297672118214306280e-01 3.082911623010851088e-01 2.189939025734843892e-01 1.542969824575534943e-01 1.079887483285981026e-01 7.517055853999345116e-02 5.209894317475911862e-02 3.598356435433126138e-02 2.478475128904975705e-02 1.703415184231726540e-02 1.168724536889349623e-02 8.007897514648740223e-03 5.481069532485032501e-03 3.748450071894476442e-03 2.561862103790463890e-03 1.750002219125025267e-03 diff --git a/skypy/power_spectrum/tests/data/truth_takahashi.txt b/skypy/power_spectrum/tests/data/truth_takahashi.txt new file mode 100644 index 000000000..7ccdde9ab --- /dev/null +++ b/skypy/power_spectrum/tests/data/truth_takahashi.txt @@ -0,0 +1,5 @@ +4.307265640670998437e+02 4.930894566634436842e+02 5.644699826749822478e+02 6.461664340624427041e+02 7.396612677463605223e+02 8.466458975703293390e+02 9.690481844025078999e+02 1.109062606579969952e+03 1.269182914147416795e+03 1.452236797473573006e+03 1.661421699195238716e+03 1.900340325677399733e+03 2.173033621694780322e+03 2.484007909367616776e+03 2.838251516047000223e+03 3.241234504552677208e+03 3.698883096762272544e+03 4.217518147641132600e+03 4.803744791802210784e+03 5.464278500806233751e+03 6.205691762442397703e+03 7.034066054108157914e+03 7.954536396847890501e+03 8.970721080979541512e+03 1.008403732998496889e+04 1.129291430269161174e+04 1.259192673690367155e+04 1.397088370933270562e+04 1.541391475313140108e+04 1.689859690651369237e+04 1.839515833130610190e+04 1.986577553125675331e+04 2.126395531349316661e+04 2.253397670315176583e+04 2.361041404202426929e+04 2.441799611533794450e+04 2.487271952444912313e+04 2.488653151181447174e+04 2.437991796078946936e+04 2.330747132843887084e+04 2.169573485915516358e+04 1.967554028770895457e+04 1.747410093751401655e+04 1.534950196352905368e+04 1.350601701191330540e+04 1.204305813980304629e+04 1.091575058790722323e+04 9.865142370597644913e+03 8.492297563328149408e+03 6.770065873078516233e+03 5.364504445805525393e+03 4.672297731354738971e+03 4.005591448272348316e+03 3.111675403075856593e+03 2.660440224330383899e+03 2.187252619884640353e+03 1.832569945036744002e+03 1.528529210865050118e+03 1.304979309951210780e+03 1.114433278860039763e+03 9.605831502320733080e+02 8.349058236137399263e+02 7.299355661379674984e+02 6.404680351628857125e+02 5.625711504996520489e+02 4.933366579403630681e+02 4.307905543571280873e+02 3.736536959105059736e+02 3.211923679700412322e+02 2.730795860858679021e+02 2.292736588730092819e+02 1.898969589949338683e+02 1.551079350224469238e+02 1.249802049542186353e+02 9.942371664277064269e+01 7.817107322232011768e+01 6.081627046638685385e+01 4.687345309139724492e+01 3.583166773540314409e+01 2.719624674727603164e+01 2.051551794188647548e+01 1.539478157869204722e+01 1.150067444575264375e+01 8.559129309399482466e+00 6.349646757002925135e+00 4.697898079749368350e+00 3.468001037648732687e+00 2.555261087211646132e+00 1.879774486682399326e+00 1.381030541355414698e+00 1.013496806721725596e+00 7.430913291225851314e-01 5.444135032868325785e-01 3.986004224003147911e-01 2.916855683919730713e-01 2.133530995401369534e-01 1.559990349234907092e-01 1.140277578489452825e-01 8.332729811337266645e-02 6.087944257146948174e-02 +2.566374549793677602e+02 2.937959381880738192e+02 3.363278747588973374e+02 3.850070089716037955e+02 4.407168784981828367e+02 5.044656064779655367e+02 5.774023076434120867e+02 6.608350998717271523e+02 7.562506066128376006e+02 8.653346737965007378e+02 9.899937864338784266e+02 1.132376329987141844e+03 1.294892370748510075e+03 1.480229997873758066e+03 1.691365451408511490e+03 1.931563242835537039e+03 2.204361271353182019e+03 2.513534609252618793e+03 2.863030298186096843e+03 3.256864373798410270e+03 3.698971719254930576e+03 4.192999612882504152e+03 4.742037375859893473e+03 5.348277659340700666e+03 6.012609740972447980e+03 6.734151477224823793e+03 7.509733566524773778e+03 8.333356319195196193e+03 9.195643602341226142e+03 1.008331921959524152e+04 1.097872600987671285e+04 1.185939658908412821e+04 1.269766901151056118e+04 1.346033115230911244e+04 1.410830573953881685e+04 1.459652949186358273e+04 1.487458068934260518e+04 1.488945259512788107e+04 1.459308833111326385e+04 1.395773696546067367e+04 1.299870772973169369e+04 1.179382075265963249e+04 1.047867168683310047e+04 9.207292240966866302e+03 8.101260743527091108e+03 7.218916151502538014e+03 6.531817883374426856e+03 5.884961397291072899e+03 5.045774792120902930e+03 4.005368435461730314e+03 3.154382033091373614e+03 2.719093536384417803e+03 2.301963040196355450e+03 1.764399557270857940e+03 1.482414913255337069e+03 1.194975088935287204e+03 9.793745005835763777e+02 7.970833505947541653e+02 6.635677861633905650e+02 5.517910706982075908e+02 4.631059936848866982e+02 3.922993796007798437e+02 3.348336471002563712e+02 2.875129578959206356e+02 2.478972936834242660e+02 2.141287508339894146e+02 1.848789412958264791e+02 1.592007831461585283e+02 1.364423121111404953e+02 1.161666881512263814e+02 9.808803991082135099e+01 8.201923349032446708e+01 6.783307957493508411e+01 5.543410218709757231e+01 4.473734432675032480e+01 3.565102318682105675e+01 2.806332105640140284e+01 2.183681350676216226e+01 1.681260303272659229e+01 1.282119760432201971e+01 9.694607637571394321e+00 7.275926284221493567e+00 5.425334510803968335e+00 4.022865551535909745e+00 2.968711141551550803e+00 2.181910008671726242e+00 1.598145100810399732e+00 1.167194355345001489e+00 8.503956334427152930e-01 6.183333577850390617e-01 4.488429095138153913e-01 3.253562259588215766e-01 2.355700215106088236e-01 1.703974818020464466e-01 1.231573421744894536e-01 8.895514436898650790e-02 6.421625290177684908e-02 4.633647823629277740e-02 3.342252399546118830e-02 2.410026124380630283e-02 +1.600448095248837319e+02 1.832181140972656692e+02 2.097426306184939904e+02 2.401009579590414091e+02 2.748441894846869218e+02 3.146011457278915486e+02 3.600886161341001639e+02 4.121226052757016305e+02 4.716305130946259396e+02 5.396640780864727276e+02 6.174127640282739549e+02 7.062170591633214372e+02 8.075808636499572231e+02 9.231817477986957101e+02 1.054877353964017402e+03 1.204705580900979839e+03 1.374875439586798393e+03 1.567744640587978438e+03 1.785779142634125719e+03 2.031489190166658545e+03 2.307335981025174988e+03 2.615603266999501557e+03 2.958229145099687685e+03 3.336595242958679137e+03 3.751273502290708166e+03 4.201734662523469524e+03 4.686026896938557911e+03 5.200437106962725011e+03 5.739150136650830973e+03 6.293921496163569827e+03 6.853776021736013718e+03 7.404737728671934747e+03 7.929586221012993519e+03 8.407628902935895894e+03 8.814495487246986158e+03 9.122049291480454485e+03 9.298760039602373581e+03 9.311410718664987144e+03 9.129777116453391500e+03 8.736208110460160242e+03 8.139874065983102810e+03 7.389024186377510887e+03 6.568176692569712031e+03 5.773630852942241290e+03 5.081629919950226395e+03 4.528974775694330674e+03 4.097988791913282512e+03 3.691026471389237486e+03 3.161081404978128376e+03 2.502323773652543423e+03 1.961915175291856940e+03 1.683822662900574414e+03 1.416840122726524442e+03 1.073751009408775417e+03 8.927389887259163288e+02 7.088238750199743663e+02 5.710298145089752779e+02 4.549951789413874508e+02 3.705551773578180814e+02 3.005005603741109894e+02 2.456604652292004687e+02 2.026827863801484284e+02 1.686268673530603905e+02 1.414019129497345375e+02 1.193910594128699501e+02 1.013364925748184078e+02 8.631154429392528016e+01 7.362613556984116769e+01 6.277633101503271718e+01 5.339811725725967761e+01 4.523195745204996854e+01 3.809331872034308475e+01 3.184995996128871809e+01 2.640434535818594597e+01 2.168067072424146247e+01 1.761580319753583623e+01 1.415335943983063416e+01 1.123992655623993464e+01 8.822500250636153396e+00 6.846834696742503112e+00 5.257073823785801459e+00 3.996951618779873261e+00 3.012031247656790534e+00 2.251961027604939325e+00 1.672049916213899090e+00 1.234011436985801069e+00 9.060148386309134771e-01 6.622582007600448240e-01 4.822677625709745897e-01 3.500857586925881226e-01 2.534590599481115758e-01 1.830957832490689163e-01 1.320223301792243231e-01 9.504978707413792860e-02 6.834455013297406623e-02 4.909091244956988870e-02 3.523078280797605627e-02 2.526582312348636325e-02 1.810878192078373139e-02 1.297283334788827160e-02 +1.068590681396608062e+02 1.223316460795109890e+02 1.400418491866120689e+02 1.603119292089488965e+02 1.835098801244325557e+02 2.100556057382765687e+02 2.404277619971090587e+02 2.751712712416614863e+02 3.149054617516602548e+02 3.603327188442603983e+02 4.122474348974250233e+02 4.715449044924264399e+02 5.392296153613460774e+02 6.164221234591439043e+02 7.043633604304251321e+02 8.044147986324705926e+02 9.180523984450588841e+02 1.046851709254309981e+03 1.192460940976926395e+03 1.356558353994014624e+03 1.540790056441992192e+03 1.746684404288879477e+03 1.975539835186776145e+03 2.228284263066102540e+03 2.505306159859875152e+03 2.806260048712626030e+03 3.129852032203881663e+03 3.473613679713103011e+03 3.833674432349369454e+03 4.204542894450379208e+03 4.578905257056580922e+03 4.947444283459657527e+03 5.298675620532538233e+03 5.618793985877475279e+03 5.891533057808501326e+03 6.098101275102925683e+03 6.217422465271039982e+03 6.227262960123927769e+03 6.107341069922797033e+03 5.845715592169390220e+03 5.448318296947064482e+03 4.947203649151039826e+03 4.398768771789787934e+03 3.867463760089305651e+03 3.404513382538053065e+03 3.034890456357917628e+03 2.747053491836332796e+03 2.475282946129894754e+03 2.119565129707965298e+03 1.674903607979238132e+03 1.309327874928307665e+03 1.122151540948170350e+03 9.419469237429468649e+02 7.082216812562309087e+02 5.856906842304076690e+02 4.604666849814812508e+02 3.667762738000701574e+02 2.878046976019009549e+02 2.306611573423447794e+02 1.833661610952427452e+02 1.466232123520373989e+02 1.181545763329396266e+02 9.593729183286072271e+01 7.852766078851192333e+01 6.479902201675587037e+01 5.386322253489907297e+01 4.505697602785391354e+01 3.787693955337223173e+01 3.194694029745793173e+01 2.698746083390275885e+01 2.279283364593328542e+01 1.921228478206160162e+01 1.613560476168447977e+01 1.348198408632537415e+01 1.119148718313602586e+01 9.218543801688994677e+00 7.527066551513696524e+00 6.086945017049108841e+00 4.871691301408283437e+00 3.857009779398825167e+00 3.020003398631996738e+00 2.338743865728087457e+00 1.792069864967342774e+00 1.359639806457449751e+00 1.022265793735524397e+00 7.623959016271072420e-01 5.645245500138303951e-01 4.153990885317878634e-01 3.040175786648222633e-01 2.214729288904968285e-01 1.607069213399734409e-01 1.162271691412471442e-01 8.382471578077654861e-02 6.031528806613471738e-02 4.331551504592954377e-02 3.105730987187258971e-02 2.223867755630027504e-02 1.590667844556230281e-02 1.136734901306653972e-02 8.117409412399672480e-03 +7.569661560088660224e+01 8.665712412121320085e+01 9.920275132448260536e+01 1.135617855914026677e+02 1.299949221925355118e+02 1.487996335945859983e+02 1.703150181339045730e+02 1.949271252150973908e+02 2.230747241537352750e+02 2.552554363010440284e+02 2.920320800933119472e+02 3.340389787300423450e+02 3.819878417713579211e+02 4.366726461815823086e+02 4.989727015596757838e+02 5.698527846710610447e+02 6.503588739773073257e+02 7.416076229167205156e+02 8.447673178312820710e+02 9.610277340316514483e+02 1.091556119711737210e+03 1.237436612233084588e+03 1.399590840587212597e+03 1.578678384615464438e+03 1.774977176334740534e+03 1.988245767281655617e+03 2.217571437719482674e+03 2.461210035050604347e+03 2.716424726947648196e+03 2.979331002578181142e+03 3.244753747706088689e+03 3.506098798232423633e+03 3.755236626800945032e+03 3.982392767788485344e+03 4.176047444709057345e+03 4.322888002467423576e+03 4.407975590894411653e+03 4.415537449346890753e+03 4.331161608411184716e+03 4.146316883102971588e+03 3.865112942192857190e+03 3.510176054140917131e+03 3.121438484103347946e+03 2.744635357564802234e+03 2.416227704959256698e+03 2.154138909412273733e+03 1.950371969309190263e+03 1.758130989010129269e+03 1.505512761503495085e+03 1.188238273359719187e+03 9.270201331461443033e+02 7.941443444678848209e+02 6.659561129139799505e+02 4.980278665293651557e+02 4.107220714589107047e+02 3.208873586949387686e+02 2.537687582431594251e+02 1.970726139392105836e+02 1.562331264281119161e+02 1.224233980805752253e+02 9.626528239465909564e+01 7.613177073125643801e+01 6.056137224592987423e+01 4.851087143171739768e+01 3.916134230936003746e+01 3.186156926690134128e+01 2.612198311342227086e+01 2.156755717551892459e+01 1.791514443325376504e+01 1.495185484128012021e+01 1.251903517860096748e+01 1.049878310205313703e+01 8.803820804141428980e+00 7.369623183239951913e+00 6.148461532948975794e+00 5.104872907153407802e+00 4.212206016181656132e+00 3.450004256264182168e+00 2.802025017713299704e+00 2.254783410194939997e+00 1.796540437689756686e+00 1.416669669894781380e+00 1.105320997942948580e+00 8.532792911082427434e-01 6.519219579007310417e-01 4.932240533609070776e-01 3.698071903194494192e-01 2.750290614150741053e-01 2.030753181929662032e-01 1.490062378801289000e-01 1.087404284827924450e-01 7.898722057050032441e-02 5.714879560942737885e-02 4.121080934123119016e-02 2.963488674184570904e-02 2.126105368981827046e-02 1.522393805505509526e-02 1.088365043061141277e-02 7.770497516601956867e-03 5.541801156290020701e-03 diff --git a/skypy/power_spectrum/tests/test_camb.py b/skypy/power_spectrum/tests/test_camb.py new file mode 100644 index 000000000..f26e5d3fb --- /dev/null +++ b/skypy/power_spectrum/tests/test_camb.py @@ -0,0 +1,38 @@ +import numpy as np +from astropy.cosmology import Planck15 +from astropy.units import allclose +from astropy.utils.data import get_pkg_data_filename +import pytest + +# load the external camb result to test against +camb_result_filename = get_pkg_data_filename('data/camb_result.txt') +test_pzk = np.loadtxt(camb_result_filename) + +# try to import the requirement, if it doesn't exist, skip test +try: + __import__('camb') +except ImportError: + CAMB_NOT_FOUND = True +else: + CAMB_NOT_FOUND = False + + +@pytest.mark.skipif(CAMB_NOT_FOUND, reason='CAMB not found') +def test_camb(): + ''' + Test a default astropy cosmology + ''' + from skypy.power_spectrum import camb + + # test shape and compare with the mocked power spectrum + redshift = [0.0, 1.0] + wavenumber = np.logspace(-4.0, np.log10(2.0), 200) + pzk = camb(wavenumber, redshift, Planck15, 2.e-9, 0.965) + assert pzk.shape == (len(redshift), len(wavenumber)) + assert allclose(pzk, test_pzk, rtol=1.e-4) + + # also check redshifts are ordered correctly + redshift = [1.0, 0.0] + pzk = camb(wavenumber, redshift, Planck15, 2.e-9, 0.965) + assert pzk.shape == (len(redshift), len(wavenumber)) + assert allclose(pzk, test_pzk[np.argsort(redshift), :], rtol=1.e-4) diff --git a/skypy/power_spectrum/tests/test_class.py b/skypy/power_spectrum/tests/test_class.py new file mode 100644 index 000000000..f83e16477 --- /dev/null +++ b/skypy/power_spectrum/tests/test_class.py @@ -0,0 +1,46 @@ +import numpy as np +from astropy.cosmology import Planck15 +from astropy.units import allclose +from astropy import units as u +from astropy.utils.data import get_pkg_data_filename +import pytest + +# load the external class result to test against +class_result_filename = get_pkg_data_filename('data/class_result.txt') +test_pzk = np.loadtxt(class_result_filename) + +# try to import the requirement, if it doesn't exist, skip test +try: + __import__('classy') +except ImportError: + CLASS_NOT_FOUND = True +else: + CLASS_NOT_FOUND = False + + +@pytest.mark.skipif(CLASS_NOT_FOUND, reason='classy not found') +def test_classy(): + ''' + Test a default astropy cosmology + ''' + from skypy.power_spectrum import classy + + Pl15massless = Planck15.clone(name='Planck 15 massless neutrino', m_nu=[0., 0., 0.]*u.eV) + + redshift = [0.0, 1.0] + wavenumber = np.logspace(-4.0, np.log10(2.0), 200) + pzk = classy(wavenumber, redshift, Pl15massless) + assert pzk.shape == (len(redshift), len(wavenumber)) + assert allclose(pzk, test_pzk, rtol=1.e-4) + + # also check redshifts are ordered correctly + redshift = [1.0, 0.0] + pzk = classy(wavenumber, redshift, Pl15massless) + assert pzk.shape == (len(redshift), len(wavenumber)) + assert allclose(pzk, test_pzk[np.argsort(redshift)], rtol=1.e-4) + + # also check scalar arguments are treated correctly + redshift = 1.0 + wavenumber = 1.e-1 + pzk = classy(wavenumber, redshift, Pl15massless) + assert np.isscalar(pzk) diff --git a/skypy/power_spectrum/tests/test_eisenstein_hu.py b/skypy/power_spectrum/tests/test_eisenstein_hu.py new file mode 100644 index 000000000..fa670c876 --- /dev/null +++ b/skypy/power_spectrum/tests/test_eisenstein_hu.py @@ -0,0 +1,76 @@ +import numpy as np +import pytest +from astropy.cosmology import Planck15, FlatLambdaCDM +from skypy.power_spectrum import eisenstein_hu + + +def test_eisenstein_hu(): + """ Test Eisenstein & Hu Linear matter power spectrum with + and without wiggles using Planck15 cosmology""" + cosmology = Planck15 + A_s = 2.1982e-09 + n_s = 0.969453 + kwmap = 0.02 + + # Test that a scalar input gives a scalar output + scalar_input = 1 + scalar_output_w = eisenstein_hu(scalar_input, A_s, n_s, cosmology, kwmap, + wiggle=True) + scalar_output_nw = eisenstein_hu(scalar_input, A_s, n_s, cosmology, kwmap, + wiggle=False) + assert np.isscalar(scalar_output_w) + assert np.isscalar(scalar_output_nw) + + # Test that an array input gives an array output + array_shape = (10,) + array_input = np.random.uniform(size=array_shape) + array_output_w = eisenstein_hu(array_input, A_s, n_s, cosmology, kwmap, + wiggle=True) + array_output_nw = eisenstein_hu(array_input, A_s, n_s, cosmology, kwmap, + wiggle=False) + assert array_output_w.shape == array_shape + assert array_output_nw.shape == array_shape + + # Test pk against precomputed values for Planck15 cosmology + wavenumber = np.logspace(-3, 1, num=5, base=10.0) + pk_eisensteinhu_w = eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap, + wiggle=True) + pk_eisensteinhu_nw = eisenstein_hu(wavenumber, A_s, n_s, cosmology, kwmap, + wiggle=False) + pk_cosmosis_w = np.array([6.47460158e+03, 3.71610099e+04, 9.65702614e+03, + 1.14604456e+02, 3.91399918e-01]) + pk_cosmosis_nw = np.array([6.47218600e+03, 3.77330704e+04, 1.00062077e+04, + 1.13082980e+02, 3.83094714e-01]) + + assert np.allclose(pk_eisensteinhu_w, pk_cosmosis_w) + assert np.allclose(pk_eisensteinhu_nw, pk_cosmosis_nw) + + # Test for failure when wavenumber <= 0 + negative_wavenumber_scalar = 0 + with pytest.raises(ValueError): + eisenstein_hu(negative_wavenumber_scalar, A_s, n_s, cosmology, kwmap, + wiggle=True) + with pytest.raises(ValueError): + eisenstein_hu(negative_wavenumber_scalar, A_s, n_s, cosmology, kwmap, + wiggle=False) + negative_wavenumber_array = [0, 1, -2, 3] + with pytest.raises(ValueError): + eisenstein_hu(negative_wavenumber_array, A_s, n_s, cosmology, kwmap, + wiggle=True) + with pytest.raises(ValueError): + eisenstein_hu(negative_wavenumber_array, A_s, n_s, cosmology, kwmap, + wiggle=False) + + # Test for failure when cosmology has Ob0 = 0 and wiggle = True + zero_ob0_cosmology = FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=2.725) + wavenumber = np.logspace(-3, 1, num=5, base=10.0) + with pytest.raises(ValueError): + eisenstein_hu(wavenumber, A_s, n_s, zero_ob0_cosmology, kwmap, + wiggle=True) + + # Test for failure when cosmology has Tcmb = 0 and wiggle = True + zero_Tcmb0_cosmology = FlatLambdaCDM(H0=70, Om0=0.3, Ob0=0.05) + wavenumber = np.logspace(-3, 1, num=5, base=10.0) + with pytest.raises(ValueError): + eisenstein_hu(wavenumber, A_s, n_s, zero_Tcmb0_cosmology, kwmap, + wiggle=True) diff --git a/skypy/power_spectrum/tests/test_growth.py b/skypy/power_spectrum/tests/test_growth.py new file mode 100644 index 000000000..efc65fc5f --- /dev/null +++ b/skypy/power_spectrum/tests/test_growth.py @@ -0,0 +1,97 @@ +import numpy as np +from astropy.cosmology import FlatLambdaCDM, Planck15 +from astropy.units import isclose, allclose +import pytest + +from skypy.power_spectrum import (growth_function_carroll, growth_factor, + growth_function, growth_function_derivative) + + +def test_carroll(): + """ Test a FlatLambdaCDM cosmology with omega_matter = 1""" + cosmology = FlatLambdaCDM(Om0=1.0, H0=70.0) + + # Test that a scalar input gives a scalar output + scalar_input = 1 + scalar_output = growth_function_carroll(scalar_input, cosmology) + assert np.isscalar(scalar_output) + + # Test that an array input gives an array output + array_shape = (10,) + array_input = np.random.uniform(size=array_shape) + array_output = growth_function_carroll(array_input, cosmology) + assert array_output.shape == array_shape + + # Test against theory for omega_matter = 1.0 + redshift = np.linspace(0.0, 10.0, 100) + Dz_carroll = growth_function_carroll(redshift, cosmology) + Dz_theory = 1.0 / (1.0 + redshift) + assert allclose(Dz_carroll, Dz_theory) + + # Test against precomputed values for Planck15 + redshift = np.linspace(0, 5, 4) + Dz_carroll = growth_function_carroll(redshift, Planck15) + Dz_truth = np.array([0.78136173, 0.36635322, 0.22889793, 0.16577711]) + assert allclose(Dz_carroll, Dz_truth) + + # Test for failure when redshift < 0 + negative_redshift_scalar = -1 + with pytest.raises(ValueError): + growth_function_carroll(negative_redshift_scalar, cosmology) + negative_redshift_array = [0, 1, -2, 3] + with pytest.raises(ValueError): + growth_function_carroll(negative_redshift_array, cosmology) + + +def test_growth(): + """ + Test a FlatLambdaCDM cosmology with omega_matter = 1.0 + and Planck15 cosmology + """ + + redshift = np.linspace(0., 10., 101) + cosmology_flat = FlatLambdaCDM(H0=70.0, Om0=1.0) + + fz = growth_factor(redshift, cosmology_flat) + Dz = growth_function(redshift, cosmology_flat) + Dzprime = growth_function_derivative(redshift, cosmology_flat) + + # Test growth factor + assert redshift.shape == fz.shape,\ + "Length of redshift array and growth rate array do not match" + assert isclose(fz[0], 1.0),\ + "Growth factor at redshift 0 is not close to 1.0" + + # Test growth function + assert redshift.shape == Dz.shape,\ + "Length of redshift array and growth function array do not match" + assert allclose(Dz, 1. / (1. + redshift)),\ + "Growth function is not close to the scale factor" + + # make sure that growth_function with scalar returns scalar + Dz2 = growth_function(redshift[2], cosmology_flat) + assert np.isscalar(Dz2), 'growth function with scalar did not produce scalar' + assert Dz2 == Dz[2], 'growth function with scalar produced inconsistent result' + + # Test growth function derivative + assert redshift.shape == Dzprime.shape,\ + "Length of redshift array and growth function array do not match" + assert isclose(Dzprime[0], -1.0),\ + "Derivative of growth function at redshift 0 is not close to -1.0" + + # Test against precomputed values using Planck15 cosmology + zvec = np.linspace(0.0, 1.0, 4) + fz_planck15 = growth_factor(zvec, Planck15) + Dz_planck15 = growth_function(zvec, Planck15) + Dzprime_planck15 = growth_function_derivative(zvec, Planck15) + + precomputed_fz_planck15 = np.array([0.5255848, 0.69412802, 0.80439553, + 0.87179376]) + precomputed_Dz_planck15 = np.array([0.66328939, 0.55638978, 0.4704842, + 0.40368459]) + precomputed_Dzprime_planck15 = np.array([-0.34861482, -0.2896543, + -0.22707323, -0.17596485]) + + assert allclose(fz_planck15, precomputed_fz_planck15) + assert allclose(Dz_planck15, precomputed_Dz_planck15) + assert allclose(Dzprime_planck15, precomputed_Dzprime_planck15) diff --git a/skypy/power_spectrum/tests/test_halofit.py b/skypy/power_spectrum/tests/test_halofit.py new file mode 100644 index 000000000..5dfc2dbff --- /dev/null +++ b/skypy/power_spectrum/tests/test_halofit.py @@ -0,0 +1,95 @@ +import numpy as np +from astropy.cosmology import Planck15 +from astropy.units import allclose +from astropy.utils.data import get_pkg_data_filename +import pytest +from skypy.power_spectrum import halofit_smith +from skypy.power_spectrum import halofit_takahashi +from skypy.power_spectrum import halofit_bird + + +# Power spectrum data for tests +linear_power_filename = get_pkg_data_filename('data/linear_power.txt') +truth_smith_filename = get_pkg_data_filename('data/truth_smith.txt') +truth_takahashi_filename = get_pkg_data_filename('data/truth_takahashi.txt') +truth_bird_filename = get_pkg_data_filename('data/truth_bird.txt') +linear_power = np.loadtxt(linear_power_filename) +truth_smith = np.loadtxt(truth_smith_filename) +truth_takahashi = np.loadtxt(truth_takahashi_filename) +truth_bird = np.loadtxt(truth_bird_filename) + + +def test_halofit(): + """Test Smith, Takahashi and Bird Halofit models with Planck15 cosmology""" + + # Wavenumbers and redshifts for tests + k = np.logspace(-4, 2, 100, base=10) + z = np.linspace(0, 2, 5) + + # Non-linear power spectra from Smith, Takahashi and Bird models + nl_power_smith = halofit_smith(k, z, linear_power, Planck15) + nl_power_takahashi = halofit_takahashi(k, z, linear_power, Planck15) + nl_power_bird = halofit_bird(k, z, linear_power, Planck15) + + # Test shape of outputs + assert np.shape(nl_power_smith) == np.shape(linear_power) + assert np.shape(nl_power_takahashi) == np.shape(linear_power) + assert np.shape(nl_power_bird) == np.shape(linear_power) + + # Test outputs against precomputed values + assert allclose(nl_power_smith, truth_smith) + assert allclose(nl_power_takahashi, truth_takahashi) + assert allclose(nl_power_bird, truth_bird) + + # Test when redshift is a scalar + z_scalar = z[0] + power_1d = linear_power[0, :] + truth_scalar_redshift = truth_smith[0, :] + smith_scalar_redshift = halofit_smith(k, z_scalar, power_1d, Planck15) + assert allclose(smith_scalar_redshift, truth_scalar_redshift) + assert np.shape(smith_scalar_redshift) == np.shape(power_1d) + + # Test for failure when wavenumber is a scalar + k_scalar = k[0] + power_1d = linear_power[:, 0] + with pytest.raises(TypeError): + halofit_smith(k_scalar, z, power_1d, Planck15) + + # Test for failure when wavenumber array is the wrong size + k_wrong_size = np.logspace(-4.0, 2.0, 7) + with pytest.raises(ValueError): + halofit_smith(k_wrong_size, z, linear_power, Planck15) + + # Test for failure when redshift array is the wrong size + z_wrong_size = np.linspace(0.0, 2.0, 3) + with pytest.raises(ValueError): + halofit_smith(k, z_wrong_size, linear_power, Planck15) + + # Test for failure when wavenumber is negative + k_negative = np.copy(k) + k_negative[0] = -1.0 + with pytest.raises(ValueError): + halofit_smith(k_negative, z, linear_power, Planck15) + + # Test for failure when wavenumber is zero + k_zero = np.copy(k) + k_zero[0] = 0.0 + with pytest.raises(ValueError): + halofit_smith(k_zero, z, linear_power, Planck15) + + # Test for failure when redshift is negative + z_negative = np.copy(z) + z_negative[0] = -1.0 + with pytest.raises(ValueError): + halofit_smith(k, z_negative, linear_power, Planck15) + + # Test for failure when linear power spectrum is negative + power_negative = np.copy(linear_power) + power_negative[0, 0] = -1.0 + with pytest.raises(ValueError): + halofit_smith(k, z, power_negative, Planck15) + + # Test for failure when wavenumber array is not in asscending order + k_wrong_order = k[::-1] + with pytest.raises(ValueError): + halofit_smith(k_wrong_order, z, linear_power, Planck15) diff --git a/skypy/power_spectrum/tests/test_import.py b/skypy/power_spectrum/tests/test_import.py new file mode 100644 index 000000000..9f3e2b752 --- /dev/null +++ b/skypy/power_spectrum/tests/test_import.py @@ -0,0 +1,2 @@ +def test_import(): + import skypy.power_spectrum # noqa: F401