Commit 05de536f authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

last healpix

parent a16d3231
/*
* This file is part of Healpix_cxx.
*
* Healpix_cxx 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.
*
* Healpix_cxx 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 Healpix_cxx; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
* For more information about HEALPix, see http://healpix.jpl.nasa.gov
*/
/*
* Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*! \file alm.cc
* Class for storing spherical harmonic coefficients.
*
* Copyright (C) 2003-2011 Max-Planck-Society
* \author Martin Reinecke
*/
#include "alm.h"
using namespace std;
//static
tsize Alm_Base::Num_Alms (int l, int m)
{
planck_assert(m<=l,"mmax must not be larger than lmax");
return ((m+1)*(m+2))/2 + (m+1)*(l-m);
}
void Alm_Base::swap (Alm_Base &other)
{
std::swap(lmax, other.lmax);
std::swap(mmax, other.mmax);
std::swap(tval, other.tval);
}
/*
* This file is part of Healpix_cxx.
*
* Healpix_cxx 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.
*
* Healpix_cxx 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 Healpix_cxx; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
* For more information about HEALPix, see http://healpix.jpl.nasa.gov
*/
/*
* Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*
* Copyright (C) 2011 Max-Planck-Society
* Author: Martin Reinecke
*/
#include "healpix_tables.h"
#include "string_utils.h"
using namespace std;
const nside_dummy SET_NSIDE=nside_dummy();
Healpix_Ordering_Scheme string2HealpixScheme (const string &inp)
{
string tmp=trim(inp);
if (equal_nocase(tmp,"RING")) return RING;
if (equal_nocase(tmp,"NESTED")) return NEST;
planck_fail ("bad Healpix ordering scheme '"+tmp+
"': expected 'RING' or 'NESTED'");
}
const uint16 Healpix_Tables::utab[] = {
#define Z(a) 0x##a##0, 0x##a##1, 0x##a##4, 0x##a##5
#define Y(a) Z(a##0), Z(a##1), Z(a##4), Z(a##5)
#define X(a) Y(a##0), Y(a##1), Y(a##4), Y(a##5)
X(0),X(1),X(4),X(5)
#undef X
#undef Y
#undef Z
};
const uint16 Healpix_Tables::ctab[] = {
#define Z(a) a,a+1,a+256,a+257
#define Y(a) Z(a),Z(a+2),Z(a+512),Z(a+514)
#define X(a) Y(a),Y(a+4),Y(a+1024),Y(a+1028)
X(0),X(8),X(2048),X(2056)
#undef X
#undef Y
#undef Z
};
const int Healpix_Tables::jrll[] = { 2,2,2,2,3,3,3,3,4,4,4,4 },
Healpix_Tables::jpll[] = { 1,3,5,7,0,2,4,6,1,3,5,7 };
const uint8 Healpix_Tables::peano_subpix[2][8][4] =
{ { {0,1,3,2}, {3,0,2,1}, {2,3,1,0}, {1,2,0,3},
{0,3,1,2}, {1,0,2,3}, {2,1,3,0}, {3,2,0,1} },
{ {0,1,3,2}, {1,3,2,0}, {3,2,0,1}, {2,0,1,3},
{0,2,3,1}, {1,0,2,3}, {3,1,0,2}, {2,3,1,0} } };
const uint8 Healpix_Tables::peano_subpath[2][8][4] =
{ { {4,0,6,0}, {7,5,1,1}, {2,4,2,6}, {3,3,7,5},
{0,2,4,4}, {5,1,5,3}, {6,6,0,2}, {1,7,3,7} },
{ {4,0,0,6}, {5,1,1,7}, {6,2,2,4}, {7,3,3,5},
{0,4,4,2}, {1,5,5,3}, {2,6,6,0}, {3,7,7,1} } };
const uint8 Healpix_Tables::peano_face2path[2][12] =
{ { 2,5,2,5,3,6,3,6,2,3,2,3 }, { 2,6,2,3,3,5,2,6,2,3,3,5 } };
const uint8 Healpix_Tables::peano_face2face[2][12] =
{ { 0,5,6,11,10,1,4,7,2,3,8,9 }, { 0,5,8,9,6,1,2,7,10,11,4,3 } };
const int Healpix_Tables::nb_xoffset[] = { -1,-1, 0, 1, 1, 1, 0,-1 },
Healpix_Tables::nb_yoffset[] = { 0, 1, 1, 1, 0,-1,-1,-1 };
const int Healpix_Tables::nb_facearray[][12] =
{ { 8, 9,10,11,-1,-1,-1,-1,10,11, 8, 9 }, // S
{ 5, 6, 7, 4, 8, 9,10,11, 9,10,11, 8 }, // SE
{ -1,-1,-1,-1, 5, 6, 7, 4,-1,-1,-1,-1 }, // E
{ 4, 5, 6, 7,11, 8, 9,10,11, 8, 9,10 }, // SW
{ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11 }, // center
{ 1, 2, 3, 0, 0, 1, 2, 3, 5, 6, 7, 4 }, // NE
{ -1,-1,-1,-1, 7, 4, 5, 6,-1,-1,-1,-1 }, // W
{ 3, 0, 1, 2, 3, 0, 1, 2, 4, 5, 6, 7 }, // NW
{ 2, 3, 0, 1,-1,-1,-1,-1, 0, 1, 2, 3 } }; // N
const int Healpix_Tables::nb_swaparray[][3] =
{ { 0,0,3 }, // S
{ 0,0,6 }, // SE
{ 0,0,0 }, // E
{ 0,0,5 }, // SW
{ 0,0,0 }, // center
{ 5,0,0 }, // NE
{ 0,0,0 }, // W
{ 6,0,0 }, // NW
{ 3,0,0 } }; // N
const int Healpix_Tables::swap_clen[] =
{ 0,7,5,4,12,10,13,18,14,19,18,17,27,21 };
const int Healpix_Tables::swap_cycle[] =
{ 0,1,8,12,16,21,40,
0,1,2,40,114,
0,4,160,263,
0,4,30,49,51,87,526,1027,1105,1387,1807,2637,
0,8,10,18,39,74,146,307,452,4737,
0,1,2,7,9,17,80,410,1526,1921,32859,33566,38931,
0,5,6,10,12,24,27,95,372,494,924,1409,3492,4248,9137,66043,103369,156899,
0,1,2,3,4,45,125,351,697,24337,102940,266194,341855,419857,
0,1,2,3,9,16,1705,2082,2126,8177,12753,15410,52642,80493,83235,88387,99444,
1675361,2495125,
0,2,6,8,9,11,20,50,93,152,183,2137,13671,44794,486954,741908,4803258,
5692573,
0,1,5,6,44,53,470,2847,3433,4906,13654,14710,400447,1797382,2744492,
18775974,23541521,
0,4,9,10,16,33,83,117,318,451,5759,10015,128975,171834,211256,347608,
1278690,2154097,2590798,3427694,5581717,21012301,27023976,72522811,
95032729,139166747,171822389,
0,5,10,267,344,363,2968,3159,9083,18437,76602,147614,1246902,1593138,
2035574,6529391,9511830,11340287,29565945,281666026,677946848 };
/*
* This file is part of Healpix_cxx.
*
* Healpix_cxx 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.
*
* Healpix_cxx 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 Healpix_cxx; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*
* For more information about HEALPix, see http://healpix.jpl.nasa.gov
*/
/*
* Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*! \file healpix_tables.h
* Copyright (C) 2011 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef HEALPIX_TABLES_H
#define HEALPIX_TABLES_H
#include "datatypes.h"
/*! The two possible ordering schemes of a HEALPix map. */
enum Healpix_Ordering_Scheme { RING, /*!< RING scheme */
NEST /*!< NESTED scheme */
};
Healpix_Ordering_Scheme string2HealpixScheme (const std::string &inp);
class nside_dummy {};
extern const nside_dummy SET_NSIDE;
class Healpix_Tables
{
protected:
static const uint16 ctab[], utab[];
static const int jrll[], jpll[];
static const uint8 peano_subpix[2][8][4], peano_subpath[2][8][4],
peano_face2path[2][12], peano_face2face[2][12];
static const int nb_xoffset[], nb_yoffset[],
nb_facearray[][12], nb_swaparray[][3];
static const int swap_clen[], swap_cycle[];
};
#endif
#include "geom_utils.h"
#include "arr.h"
using namespace std;
namespace {
void get_circle (const arr<vec3> &point, tsize q1, tsize q2, vec3 &center,
double &cosrad)
{
center = (point[q1]+point[q2]).Norm();
cosrad = dotprod(point[q1],center);
for (tsize i=0; i<q1; ++i)
if (dotprod(point[i],center)<cosrad) // point outside the current circle
{
center=crossprod(point[q1]-point[i],point[q2]-point[i]).Norm();
cosrad=dotprod(point[i],center);
if (cosrad<0)
{ center.Flip(); cosrad=-cosrad; }
}
}
void get_circle (const arr<vec3> &point, tsize q, vec3 &center,
double &cosrad)
{
center = (point[0]+point[q]).Norm();
cosrad = dotprod(point[0],center);
for (tsize i=1; i<q; ++i)
if (dotprod(point[i],center)<cosrad) // point outside the current circle
get_circle(point,i,q,center,cosrad);
}
} // unnamed namespace
void find_enclosing_circle (const arr<vec3> &point, vec3 &center,
double &cosrad)
{
tsize np=point.size();
planck_assert(np>=3,"too few points");
center = (point[0]+point[1]).Norm();
cosrad = dotprod(point[0],center);
for (tsize i=2; i<np; ++i)
if (dotprod(point[i],center)<cosrad) // point outside the current circle
get_circle(point,i,center,cosrad);
}
/*
* 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 math_utils.h
* Various convenience mathematical functions.
*
* Copyright (C) 2002-2011 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef PLANCK_MATH_UTILS_H
#define PLANCK_MATH_UTILS_H
#include <cmath>
#include "datatypes.h"
/*! \defgroup mathutilsgroup Mathematical helper functions */
/*! \{ */
/*! 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)
{
using namespace std;
return abs(a-b) < (epsilon*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)
{
using namespace std;
return 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)
{
using namespace std;
return I(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+0.5); }
/*! Returns \a v1+v2 if \a v1<0, \a v1-v2 if \a v1>=v2, else \a v1.
\a v1 can be positive or negative; \a v2 must be positive. */
template<typename T> inline T weak_modulo (T v1, T v2)
{ return (v1>=0) ? ((v1<v2) ? v1 : (v1-v2)) : (v1+v2); }
/*! 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. */
inline double fmodulo (double v1, double v2)
{
using namespace std;
return (v1>=0) ? ((v1<v2) ? v1 : fmod(v1,v2)) : (fmod(v1,v2)+v2);
}
/*! 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; }
template <typename I, bool g4> struct isqrt_helper__
{};
template <typename I> struct isqrt_helper__ <I, false>
{
static uint32 isqrt (I arg)
{
using namespace std;
return uint32 (sqrt(arg+0.5));
}
};
template <typename I> struct isqrt_helper__ <I, true>
{
static uint32 isqrt (I arg)
{
using namespace std;
long double arg2 = static_cast<long double>(arg)+0.5;
return uint32 (sqrt(arg2));
}
};
/*! Returns the integer \a n, which fulfills \a n*n<=arg<(n+1)*(n+1). */
template<typename I> inline uint32 isqrt (I arg)
{ return isqrt_helper__<I,(sizeof(I)>4)>::isqrt(arg); }
/*! Returns the largest integer \a n that fulfills \a 2^n<=arg. */
template<typename I> inline unsigned int ilog2 (I arg)
{
unsigned int res=0;
while (arg > 0x0000FFFF) { res+=16; arg>>=16; }
if (arg > 0x000000FF) { res|=8; arg>>=8; }
if (arg > 0x0000000F) { res|=4; arg>>=4; }
if (arg > 0x00000003) { res|=2; arg>>=2; }
if (arg > 0x00000001) { res|=1; }
return res;
}
/*! Returns \a atan2(y,x) if \a x!=0 or \a y!=0; else returns 0. */
inline double safe_atan2 (double y, double x)
{
using namespace std;
return ((x==0.) && (y==0.)) ? 0.0 : atan2(y,x);
}
/*! Helper function for linear interpolation (or extrapolation).
The array must be ordered in ascending order; no two values may be equal. */
template<typename T, typename Iter, typename Comp> inline void interpol_helper
(const Iter &begin, const Iter &end, const T &val, Comp comp, tsize &idx,
T &frac)
{
using namespace std;
planck_assert((end-begin)>1,"sequence too small for interpolation");
idx = lower_bound(begin,end,val,comp)-begin;
if (idx>0) --idx;
idx = min(tsize(end-begin-2),idx);
frac = (val-begin[idx])/(begin[idx+1]-begin[idx]);
}
/*! Helper function for linear interpolation (or extrapolation).
The array must be ordered in ascending order; no two values may be equal. */
template<typename T, typename Iter> inline void interpol_helper
(const Iter &begin, const Iter &end, const T &val, tsize &idx, T &frac)
{ interpol_helper (begin,end,val,std::less<T>(),idx,frac); }
/*! \} */
template<typename T> inline bool multiequal (const T &a, const T &b, const T &c)
{ return (a==b) && (a==c); }
template<typename T> inline bool multiequal (const T &a, const T &b, const T &c,
const T &d)
{ return (a==b) && (a==c) && (a==d); }
template<typename T> inline bool multiequal (const T &a, const T &b, const T &c,
const T &d, const T &e)
{ return (a==b) && (a==c) && (a==d) && (a==e); }
template<typename T> inline bool multiequal (const T &a, const T &b, const T &c,
const T &d, const T &e, const T &f)
{ return (a==b) && (a==c) && (a==d) && (a==e) && (a==f); }
#endif
/*
* 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).
*/
/*
* Class for parsing parameter files
*
* Copyright (C) 2003-2011 Max-Planck-Society
* Authors: Martin Reinecke, Reinhard Hell
*/
#include <iostream>
#include "paramfile.h"
#include "string_utils.h"
using namespace std;
string paramfile::get_valstr(const string &key) const
{
params_type::const_iterator loc=params.find(key);
if (loc!=params.end()) return loc->second;
planck_fail ("Cannot find the key '" + key + "'.");
}
bool paramfile::param_unread (const string &key) const
{ return (read_params.find(key)==read_params.end()); }
paramfile::paramfile (const string &filename, bool verbose_)
: verbose(verbose_)
{ parse_file (filename, params); }
paramfile::paramfile (const params_type &par, bool verbose_)
: params (par), verbose(verbose_)
{}
paramfile::~paramfile()
{
if (verbose)
for (params_type::const_iterator loc=params.begin();
loc!=params.end(); ++loc)
if (param_unread(loc->first))
cout << "Parser warning: unused parameter '"
<< loc->first << "'" << endl;
}
bool paramfile::param_present(const string &key) const
{ return (params.find(key)!=params.end()); }
void paramfile::findhelper (const string &key, const string &value, NDT type,
bool deflt) const
{
string output = value;
if (type==NAT_STRING) output = "'"+output+"'";
if (verbose && param_unread(key))
cout << "Parser: " << key << " = " << output
<< (deflt ? " <default>" : "") << endl;
read_params.insert(key);
}
template<typename T> T paramfile::find (const string &key) const
{
T result = stringToData<T>(get_valstr(key));
findhelper (key, dataToString(result), nativeType<T>(), false);
return result;
}
template unsigned char paramfile::find (const string &key) const;
template signed char paramfile::find (const string &key) const;
template unsigned short paramfile::find (const string &key) const;
template short paramfile::find (const string &key) const;
template unsigned int paramfile::find (const string &key) const;
template int paramfile::find (const string &key) const;
template unsigned long paramfile::find (const string &key) const;
template long paramfile::find (const string &key) const;
template unsigned long long paramfile::find (const string &key) const;
template long long paramfile::find (const string &key) const;
template float paramfile::find (const string &key) const;
template double paramfile::find (const string &key) const;
template long double paramfile::find (const string &key) const;
template bool paramfile::find (const string &key) const;
template string paramfile::find (const string &key) const;
template<typename T> T paramfile::find (const string &key, const T &deflt)
{
if (param_present(key)) return find<T>(key);
string sdeflt=dataToString(deflt);
findhelper (key, sdeflt, nativeType<T>(), true);
params[key]=sdeflt;
return deflt;
}
template unsigned char paramfile::find (const string &key,
const unsigned char &deflt);
template signed char paramfile::find (const string &key,
const signed char &deflt);
template unsigned short paramfile::find (const string &key,
const unsigned short &deflt);
template short paramfile::find (const string &key, const short &deflt);
template unsigned int paramfile::find (const string &key,
const unsigned int &deflt);