more

.gitignore

+# never store the git version file
+git_version.py
+
+# custom
+setup.cfg
+.idea
+.DS_Store
+*.pyc
+*.html
+.document
+.svn/
+*.csv
+.pytest_cache/
+*.png
+
+# from https://github.com/github/gitignore/blob/master/Python.gitignore
+
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.c
+*.o
+*.so
+*.lo
+*.la
+
+# Distribution / packaging
+.Python
+env/
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*,cover
+.hypothesis/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/build/
+docs/source/mod
+
+# PyBuilder
+target/
+
+# IPython Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# dotenv
+.env
+
+# virtualenv
+venv/
+ENV/
+
+# Spyder project settings
+.spyderproject
+
+# Rope project settings
+.ropeproject

Makefile.am

@@ -14,6 +14,10 @@ libmrutil_la_SOURCES = \
   mr_util/cmplx.h \
   mr_util/aligned_array.h \
   mr_util/simd.h \
+  mr_util/pointing.cc \
+  mr_util/pointing.h \
+  mr_util/vec3.h \
+  mr_util/constants.h \
   mr_util/error_handling.cc \
   mr_util/error_handling.h \
   mr_util/morton_utils.cc \

mr_util/constants.h

+/*
+ *  This file is part of libcxxsupport.
+ *
+ *  libcxxsupport is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  libcxxsupport is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with libcxxsupport; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
+ *  (DLR).
+ */
+
+/*! \file lsconstants.h
+ *  Mathematical, physical and technical constants for LevelS.
+ */
+
+#ifndef MRUTIL_CONSTANTS_H
+#define MRUTIL_CONSTANTS_H
+
+namespace mr {
+
+/*! \defgroup mathconstgroup Mathematical constants */
+/*! \{ */
+
+constexpr double pi=3.141592653589793238462643383279502884197;
+constexpr double twopi=6.283185307179586476925286766559005768394;
+constexpr double inv_twopi=1.0/twopi;
+constexpr double fourpi=12.56637061435917295385057353311801153679;
+constexpr double halfpi=1.570796326794896619231321691639751442099;
+constexpr double inv_halfpi=0.6366197723675813430755350534900574;
+constexpr double inv_sqrt4pi = 0.2820947917738781434740397257803862929220;
+
+constexpr double ln2 = 0.6931471805599453094172321214581766;
+constexpr double inv_ln2 = 1.4426950408889634073599246810018921;
+constexpr double ln10 = 2.3025850929940456840179914546843642;
+
+constexpr double onethird=1.0/3.0;
+constexpr double twothird=2.0/3.0;
+constexpr double fourthird=4.0/3.0;
+
+constexpr double degr2rad=pi/180.0;
+constexpr double arcmin2rad=degr2rad/60;
+constexpr double rad2degr=180.0/pi;
+
+//! Ratio between FWHM and sigma of a Gauss curve (\f$\sqrt{8\ln2}\f$).
+constexpr double sigma2fwhm=2.3548200450309493; // sqrt(8*log(2.))
+constexpr double fwhm2sigma=1/sigma2fwhm;
+
+/*! \} */
+
+/*! \defgroup physconstgroup Physical constants */
+/*! \{ */
+
+constexpr double Jansky2SI=1.0e-26;
+constexpr double SI2Jansky=1.0e+26;
+
+//! Light speed in m/s (CODATA 2006)
+constexpr double speedOfLight=2.99792458e8;
+
+//! Boltzmann's constant in J/K (CODATA 2006)
+constexpr double kBoltzmann=1.3806504e-23;
+
+//! Stefan-Boltzmann constant in W/m^2/K^4 (CODATA 2006)
+constexpr double sigmaStefanBoltzmann=5.6704e-8;
+
+//! Planck's constant in J s (CODATA 2006)
+constexpr double hPlanck=6.62606896e-34;
+
+//! Astronomical unit in m
+constexpr double astronomicalUnit=1.49597870691e11;
+
+//! Solar constant in W/m^2
+constexpr double solarConstant=1368.0;
+
+//! Average CMB temperature in K (Fixsen, ApJ 709 (2009), arXiv:0911.1955)
+constexpr double tcmb = 2.72548;
+
+//! offset (in seconds) between Jan 1, 1958 and Jan 1, 1970
+constexpr double sec_58_70 = 378691200.;
+
+/*! \} */
+
+}
+
+#endif

mr_util/math_utils.h

@@ -33,15 +33,152 @@
 #define MRUTIL_MATH_UTILS_H
 
 #include <cmath>
+#include <cstdint>
 
 namespace mr {
 
 /*! Returns \e true if | \a a-b | <= \a epsilon * | \a b |, else \e false. */
 template<typename F> inline bool approx (F a, F b, F epsilon=1e-5)
+  { return std::abs(a-b) <= (epsilon*std::abs(b)); }
+
+/*! Returns \e true if | \a a-b | <= \a epsilon, else \e false. */
+template<typename F> inline bool abs_approx (F a, F b, F epsilon=1e-5)
+  { return std::abs(a-b) <= epsilon; }
+
+/*! Returns the largest integer which is smaller than (or equal to) \a arg. */
+template<typename I, typename F> inline I ifloor (F arg)
+  { return I(std::floor(arg)); }
+
+/*! Returns the integer which is nearest to \a arg. */
+template<typename I, typename F> inline I nearest (F arg)
+  { return ifloor<I>(arg+F(0.5)); }
+
+/*! Returns the remainder of the division \a v1/v2.
+    The result is non-negative.
+    \a v1 can be positive or negative; \a v2 must be positive. */
+template<typename F> inline double fmodulo (F v1, F v2)
   {
-  return std::abs(a-b) <= (epsilon*std::abs(b));
+  using namespace std;
+  if (v1>=0)
+    return (v1<v2) ? v1 : std::fmod(v1,v2);
+  double tmp=std::fmod(v1,v2)+v2;
+  return (tmp==v2) ? F(0) : tmp;
   }
 
+/*! Returns the remainder of the division \a v1/v2.
+    The result is non-negative.
+    \a v1 can be positive or negative; \a v2 must be positive. */
+template<typename I> inline I imodulo (I v1, I v2)
+  { I v=v1%v2; return (v>=0) ? v : v+v2; }
+
+/*! Returns -1 if \a signvalue is negative, else +1. */
+template<typename T> inline T sign (const T& signvalue)
+  { return (signvalue>=0) ? 1 : -1; }
+
+/*! Returns \a val*pow(-1,m) */
+template<typename T, typename I> inline T xpow (I m, T val)
+  { return (m&1) ? -val : val; }
+
+namespace math_utils_detail {
+
+template<typename I, bool g4> struct isqrt_helper__
+  {};
+template<typename I> struct isqrt_helper__ <I, false>
+  {
+  static std::uint32_t isqrt (I arg)
+    {
+    using namespace std;
+    return std::uint32_t (sqrt(arg+0.5));
+    }
+  };
+template<typename I> struct isqrt_helper__ <I, true>
+  {
+  static std::uint32_t isqrt (I arg)
+    {
+    using namespace std;
+    I res = std::sqrt(double(arg)+0.5);
+    if (arg<(std::uint64_t(1)<<50)) return std::uint32_t(res);
+    if (res*res>arg)
+      --res;
+    else if ((res+1)*(res+1)<=arg)
+      ++res;
+    return std::uint32_t(res);
+    }
+  };
+
+/*! Returns the integer \a n, which fulfills \a n*n<=arg<(n+1)*(n+1). */
+template<typename I> inline std::uint32_t isqrt (I arg)
+  { return isqrt_helper__<I,(sizeof(I)>4)>::isqrt(arg); }
+}
+
+using math_utils_detail::isqrt;
+
+/*! Returns the largest integer \a n that fulfills \a 2^n<=arg. */
+template<typename I> inline int ilog2 (I arg)
+  {
+#ifdef __GNUC__
+  if (arg==0) return 0;
+  if (sizeof(I)==sizeof(int))
+    return 8*sizeof(int)-1-__builtin_clz(arg);
+  if (sizeof(I)==sizeof(long))
+    return 8*sizeof(long)-1-__builtin_clzl(arg);
+  if (sizeof(I)==sizeof(long long))
+    return 8*sizeof(long long)-1-__builtin_clzll(arg);
+#endif
+  int res=0;
+  while (arg > 0xFFFF) { res+=16; arg>>=16; }
+  if (arg > 0x00FF) { res|=8; arg>>=8; }
+  if (arg > 0x000F) { res|=4; arg>>=4; }
+  if (arg > 0x0003) { res|=2; arg>>=2; }
+  if (arg > 0x0001) { res|=1; }
+  return res;
+  }
+
+template<typename I> inline int ilog2_nonnull (I arg)
+  {
+#ifdef __GNUC__
+  if (sizeof(I)<=sizeof(int))
+    return 8*sizeof(int)-1-__builtin_clz(arg);
+  if (sizeof(I)==sizeof(long))
+    return 8*sizeof(long)-1-__builtin_clzl(arg);
+  if (sizeof(I)==sizeof(long long))
+    return 8*sizeof(long long)-1-__builtin_clzll(arg);
+#endif
+  return ilog2 (arg);
+  }
+
+template<typename I> inline int trailingZeros(I arg)
+  {
+  if (arg==0) return sizeof(I)<<3;
+#ifdef __GNUC__
+  if (sizeof(I)<=sizeof(int))
+    return __builtin_ctz(arg);
+  if (sizeof(I)==sizeof(long))
+    return __builtin_ctzl(arg);
+  if (sizeof(I)==sizeof(long long))
+    return __builtin_ctzll(arg);
+#endif
+  int res=0;
+  while ((arg&0xFFFF)==0) { res+=16; arg>>=16; }
+  if ((arg&0x00FF)==0) { res|=8; arg>>=8; }
+  if ((arg&0x000F)==0) { res|=4; arg>>=4; }
+  if ((arg&0x0003)==0) { res|=2; arg>>=2; }
+  if ((arg&0x0001)==0) { res|=1; }
+  return res;
+  }
+
+template<typename T1, typename T2>
+inline bool multiequal (const T1 &a, const T2 &b)
+  { return (a==b); }
+
+template<typename T1, typename T2, typename... Args>
+inline bool multiequal (const T1 &a, const T2 &b, Args... args)
+  { return (a==b) && multiequal (a, args...); }
+
+/*! Returns \a atan2(y,x) if \a x!=0 or \a y!=0; else returns 0. */
+inline double safe_atan2 (double y, double x)
+  { return ((x==0.) && (y==0.)) ? 0.0 : std::atan2(y,x); }
+
 }
 
 #endif

mr_util/pointing.cc

+/*
+ *  This file is part of libcxxsupport.
+ *
+ *  libcxxsupport is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  libcxxsupport is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with libcxxsupport; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
+ *  (DLR).
+ */
+
+/*! \file pointing.cc
+ *  Class representing a direction in 3D space
+ *
+ *  Copyright (C) 2003-2012 Max-Planck-Society
+ *  \author Martin Reinecke
+ */
+
+#include "mr_util/pointing.h"
+#include "mr_util/constants.h"
+#include "mr_util/math_utils.h"
+
+using namespace std;
+
+namespace mr {
+
+vec3 pointing::to_vec3() const
+  {
+  double st=std::sin(theta);
+  return vec3 (st*std::cos(phi), st*std::sin(phi), std::cos(theta));
+  }
+void pointing::from_vec3 (const vec3 &inp)
+  {
+  theta = std::atan2(sqrt(inp.x*inp.x+inp.y*inp.y),inp.z);
+  phi = safe_atan2 (inp.y,inp.x);
+  if (phi<0.) phi += twopi;
+  }
+void pointing::normalize_theta()
+  {
+  theta=fmodulo(theta,twopi);
+  if (theta>pi)
+    {
+    phi+=pi;
+    theta=twopi-theta;
+    }
+  }
+void pointing::normalize()
+  {
+  normalize_theta();
+  phi=fmodulo(phi,twopi);
+  }
+
+ostream &operator<< (ostream &os, const pointing &p)
+  {
+  os << p.theta << ", " << p.phi << std::endl;
+  return os;
+  }
+
+}

mr_util/pointing.h

+/*
+ *  This file is part of libcxxsupport.
+ *
+ *  libcxxsupport is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  libcxxsupport is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with libcxxsupport; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
+ *  (DLR).
+ */
+
+/*! \file pointing.h
+ *  Class representing a direction in 3D space
+ *
+ *  Copyright (C) 2003-2020 Max-Planck-Society
+ *  \author Martin Reinecke
+ */
+
+#ifndef MRUTIL_POINTING_H
+#define MRUTIL_POINTING_H
+
+#include <cmath>
+#include "mr_util/vec3.h"
+
+namespace mr {
+
+/*! \defgroup pointinggroup Pointings */
+/*! \{ */
+
+/*! Class representing a direction in 3D space or a location on the
+    unit sphere. All angles in radians. */
+class pointing
+  {
+  public:
+    /*! Colatitude of the pointing (i.e. the North pole is at \a theta=0). */
+    double theta;
+    /*! Longitude of the pointing. */
+    double phi;
+
+    /*! Default constructor. \a theta and \a phi are not initialized. */
+    pointing() {}
+    /*! Creates a pointing with \a Theta and \a Phi. */
+    pointing (double Theta, double Phi) : theta(Theta), phi(Phi) {}
+
+// FIXME: should become "explicit" some time
+    /*! Creates a pointing from the vector \a inp. \a inp need not be
+        normalized. */
+    pointing (const vec3 &inp)
+      { from_vec3(inp); }
+// FIXME: should be removed some time
+    /*! Returns a normalized vector pointing in the same direction. */
+    operator vec3() const
+      { return to_vec3(); }
+    /*! Returns a normalized vector pointing in the same direction. */
+    vec3 to_vec3() const;
+    /*! Converts \a inp to \a ptg. \a inp need not be normalized. */
+    void from_vec3 (const vec3 &inp);
+    /*! Changes the angles so that \a 0<=theta<=pi. */
+    void normalize_theta();
+    /*! Changes the angles so that \a 0<=theta<=pi and \a 0<=phi<2*pi. */
+    void normalize();
+  };
+
+/*! Writes \a p to \a os.
+    \relates pointing */
+std::ostream &operator<< (std::ostream &os, const pointing &p);
+
+/*! \} */
+
+}
+
+#endif

mr_util/vec3.h

+/*
+ *  This file is part of libcxxsupport.
+ *
+ *  libcxxsupport is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  libcxxsupport is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with libcxxsupport; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+/*
+ *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
+ *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
+ *  (DLR).
+ */
+
+/*! \file vec3.h
+ *  Class representing 3D cartesian vectors
+ *
+ *  Copyright (C) 2003, 2006 Max-Planck-Society
+ *  \author Martin Reinecke
+ */
+
+#ifndef MRUTIL_VEC3_H
+#define MRUTIL_VEC3_H
+
+#include <cmath>
+#include <iostream>
+
+namespace mr {
+
+/*! \defgroup vec3group 3D vectors */
+/*! \{ */
+
+/*! Class representing a 3D cartesian vector. */
+template<typename T>class vec3_t
+  {
+  public:
+    T x, /*!< x-coordinate */
+      y, /*!< y-coordinate */
+      z; /*!< z-coordinate */
+
+    /*! Default constructor. Does not initialize \a x, \a y, and \a z. */
+    vec3_t () {}
+    /*! Creates a vector with the coordinates \a xc, \a yc, and \a zc. */
+    vec3_t (T xc, T yc, T zc)
+      : x(xc), y(yc), z(zc) {}
+    template<typename T2> explicit vec3_t (const vec3_t<T2> &orig)
+      : x(orig.x), y(orig.y), z(orig.z) {}
+
+    /*! Sets the vector components to \a xc, \a yc, and \a zc. */
+    void Set (T xc, T yc, T zc)
+      { x=xc; y=yc; z=zc; }
+    /*! Creates a unit vector from a z coordinate and an azimuthal angle. */
+    void set_z_phi (T z_, T phi)
+      {
+      T sintheta = std::sqrt((T(1)-z_)*(T(1)+z_));
+      x = sintheta*std::cos(phi);
+      y = sintheta*std::sin(phi);
+      z = z_;
+      }
+
+    /*! Normalizes the vector to length 1. */
+    void Normalize ()
+      {
+      T l = T(1)/std::sqrt (x*x + y*y + z*z);
+      x*=l; y*=l; z*=l;
+      }
+
+    vec3_t Norm() const
+      {
+      vec3_t res(*this);
+      res.Normalize();
+      return res;
+      }
+
+    /*! Returns the length of the vector. */
+    T Length () const
+      { return sqrt (x*x + y*y + z*z); }
+
+    /*! Returns the squared length of the vector. */
+    T SquaredLength () const
+      { return (x*x + y*y + z*z); }
+    /*! Returns the vector with the signs of all coordinates flipped. */
+    const vec3_t operator- () const
+      { return vec3_t (-x, -y, -z); }
+    /*! Flips the signs of all coordinates. */
+    void Flip ()
+      { x=-x; y=-y; z=-z; }
+    /*! Returns (\a *this + \a vec). */
+    const vec3_t operator+ (const vec3_t &vec) const
+      { return vec3_t (x+vec.x, y+vec.y, z+vec.z); }