first tagged version.

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 <https://www.gnu.org/software/gsl/>`_ (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 <http://www.swig.org/index.php>`_. 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 <https://extension-helpers.readthedocs.io/en/latest/>`_
+to build the code as a C extension.
+
+
 =============
 Reference/API
 =============

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))

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))

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