From 271fd096f3106c1c29573cbce0f3de11ee9b9ac4 Mon Sep 17 00:00:00 2001 From: Jerry Ma Date: Wed, 24 Jun 2020 22:31:02 -0400 Subject: [PATCH] first tagged version. --- docs/index.rst | 29 +++++++ szpack_wrapper/__init__.py | 4 + szpack_wrapper/extern/SZpack.py | 52 ++++++++++-- szpack_wrapper/sz.py | 139 +++++++++++++++++++++++++++++++- 4 files changed, 218 insertions(+), 6 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index 182d14d..2c6f069 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -39,6 +39,35 @@ The license section from the README file of the ``SZpack.v1.1.1`` tarball: //================================================================================================== +============ +Installation +============ + +The SZpack depends on the +`GNU Scientific Library `_ (GSL), which has +to be installed and made available in the environment. + +The low level python interface comes along with SZPack also depends on +`SWIG `_. Please refer to the documentation +for how to install. + +To install ``szpack_wrapper``, + +.. code:: text + + $ gsl-config --version # check that GSL is installed + 2.6 + $ which swig # check that swig is installed + /usr/local/bin/swig + $ pip install szpack_wrapper + +Behind the scene, the ``SZpack`` v1.1.1 is bundled with this package in the +``szpack_wrapper/extern`` folder, and an accompanying ``setup_package.py`` file +is used to leverage the +`extension-helpers `_ +to build the code as a C extension. + + ============= Reference/API ============= diff --git a/szpack_wrapper/__init__.py b/szpack_wrapper/__init__.py index e22b484..6155bb3 100644 --- a/szpack_wrapper/__init__.py +++ b/szpack_wrapper/__init__.py @@ -6,4 +6,8 @@ from ._astropy_init import * # noqa # ---------------------------------------------------------------------------- +_excluded_from_all = set(globals().keys()) + from .sz import * # noqa + +__all__ = list(set(globals().keys()).difference(_excluded_from_all)) diff --git a/szpack_wrapper/extern/SZpack.py b/szpack_wrapper/extern/SZpack.py index c42b2f4..231caed 100644 --- a/szpack_wrapper/extern/SZpack.py +++ b/szpack_wrapper/extern/SZpack.py @@ -1,10 +1,52 @@ #! /usr/bin/env python -import re -from pathlib import Path +from . import _SZpack -__all__ = ['__version__', ] -_szpack_dir = list(Path(__file__).parent.glob('SZpack.v*'))[0] +_excluded_from_all = set(globals().keys()) -__version__ = re.match(r'SZpack\.(v.+)', _szpack_dir.name).group(1) +__version__ = 'v1.1.1' + + +# the below are copied from the SWIG generated SZpack.py file. +# TODO fix the args, add docs. + + +def compute_5d(*args): + return _SZpack.compute_5d(*args) + + +def compute_3d(*args): + return _SZpack.compute_3d(*args) + + +def compute_asym(*args): + return _SZpack.compute_asym(*args) + + +def compute_CNSN(*args): + return _SZpack.compute_CNSN(*args) + + +def compute_CNSN_opt(*args): + return _SZpack.compute_CNSN_opt(*args) + + +def compute_combo(*args): + return _SZpack.compute_combo(*args) + + +def compute_combo_means(*args): + return _SZpack.compute_combo_means(*args) + + +def compute_combo_means_ex(*args): + return _SZpack.compute_combo_means_ex(*args) + + +def compute_null(tau, TeSZ, betac_para, omega1, sigma, kappa, betac2_perp): + return _SZpack.compute_null( + tau, TeSZ, betac_para, omega1, sigma, kappa, betac2_perp) + + +__all__ = list(set(globals().keys()).difference(_excluded_from_all)) diff --git a/szpack_wrapper/sz.py b/szpack_wrapper/sz.py index 3c70f28..339a7dc 100644 --- a/szpack_wrapper/sz.py +++ b/szpack_wrapper/sz.py @@ -1,5 +1,142 @@ #! /usr/bin/env python +from astropy.cosmology import default_cosmology +from .extern import SZpack +from astropy import constants as const +import numpy as np +import astropy.units as u +from astropy import log + + +__all__ = ['SZ', ] + class SZ(object): - pass + """ This class defines high level interface to call SZpack. + + The SZpack works with the dimensionless (yet cosmology dependent) + wavelength ``x`` defined as + + .. math:: + + x \\equiv \\frac{h\\nu}{k_BT_{CMB0}} + + :meth:`wavelength_to_x` and :meth:`x_to_wavelength` can be used + to convert between ``x`` and `~astropy.units.Quantity` instance. + + Parameters + ---------- + cosmo : `~astropy.cosmology.Cosmology` + + The cosmology to use. `~astropy.cosmology.default_cosmology` is + used if not set. + """ + + def __init__(self, cosmo=None): + if cosmo is None: + cosmo = default_cosmology.get() + self._cosmo = cosmo + + @property + def cosmo(self): + return self._cosmo + + def wavelength_to_x(self, wavelength): + """Return the dimensionless wavelength.""" + frequency = wavelength.to(u.GHz, equivalencies=u.spectral()) + x = (const.h * frequency) / (const.k_B * self.cosmo.Tcmb0) + return x.to_value(u.dimensionless_unscaled) + + def x_to_wavelength(self, x): + """Return the wavelength as `~astropy.units.Quantity`.""" + return (x * (const.k_B * self.cosmo.Tcmb0) / const.h).to( + u.mm, u.spectral()) + + def surface_brightness( + self, runmode='3D', + x=None, + wavelength=None, **kwargs): + """ + Compute the SZ surface_brightness. + + The input can be either given as the dimensionless ``x`` or + as `~astropy.units.Quantity`. + + The output is the SZ surface brightness as `~astropy.units.Quantity`. + + The `runmode` and `kwargs` is passed to SZpack wrapper to execute + the calculation. + + Refer to the README file of SZpack for more information of all + run-modes supported and the associated arguments. + + Parameters + ---------- + x : `~numpy.ndarray` + + Dimensionless Observer frame photon frequency (h nu / k T0) + wavelength : `~astropy.units.Quantity` + + Observer frame photon wavelength/frequency. + + Returns + ------- + `~astropy.units.Quantity` + + The SZ surface brightness. + """ + _compute = getattr(self, f'_compute_{runmode}', None) + if _compute is None: + raise NotImplementedError( + f"compute run mode {runmode} is not implemented.") + if sum([x is None, wavelength is None]) != 1: + raise ValueError( + "specify observing frequency by one of " + "x, frequency, or wavelength.") + # compute x + if x is None: + x = self.wavelength_to_x(wavelength) + result = _compute(x, **kwargs) # in Dn(x) + # dndi value is from SZpack.h line 73 + DnDI = 13.33914078 * self.cosmo.Tcmb0.to_value(u.K) ** 3 + return result * x ** 3 * DnDI * u.MJy / u.sr + + @staticmethod + def _compute_3D( + x, Dtau, Te_keV, betac, muc, betao, muo, eps_Int=1e-4, + check_range=True + ): + if check_range: + if np.min(x) < 0.1: + raise ValueError(f"x < 0.1 at {np.argmin(x)}") + if np.min(x) > 50.0: + raise ValueError(f"x > 50. at {np.argmax(x)}") + result = np.copy(x) + log.debug( + f'szpack.compute_3d(Dtau={Dtau}, Te_keV={Te_keV}, ' + f'betac={betac}, muc={muc}, betao={betao}, muo={muo})') + SZpack.compute_3d(result, Dtau, Te_keV, betac, muc, betao, muo) + return result + + @staticmethod + def _compute_combo_means( + x, tau, TeSZ_keV, + betac_para, omega, sigma, kappa, betac_perp, + check_range=True + ): + if check_range: + if np.min(x) < 0.1: + raise ValueError(f"x < 0.1 at {np.argmin(x)}") + if np.min(x) > 50.0: + raise ValueError(f"x > 50. at {np.argmax(x)}") + if TeSZ_keV > 75.: + raise ValueError(f'TeSZ_keV > 75.') + result = np.copy(x) + log.debug( + f'szpack.compute_combo_means(tau={tau}, TeSZ_keV={TeSZ_keV}, ' + f'betac_para={betac_para}, omega={omega}, ' + f'sigma={sigma}, kappa={kappa}, betac_perp={betac_perp})') + SZpack.compute_combo_means( + result, tau, TeSZ_keV, + betac_para, omega, sigma, kappa, betac_perp) + return result