Commit 01b54810 authored by Leo Singer's avatar Leo Singer

New upstream version 3.40.0

parent 2743f9e9
......@@ -27,13 +27,14 @@
/*! \file alm.h
* Class for storing spherical harmonic coefficients.
*
* Copyright (C) 2003-2012 Max-Planck-Society
* Copyright (C) 2003-2017 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef PLANCK_ALM_H
#define PLANCK_ALM_H
#include <vector>
#include "arr.h"
/*! Base class for calculating the storage layout of spherical harmonic
......@@ -127,6 +128,15 @@ template<typename T> class Alm: public Alm_Base
for (int l=m; l<=lmax; ++l)
operator()(l,m)*=factor[l];
}
/*! \a a(l,m) *= \a factor[l] for all \a l,m. */
template<typename T2> void ScaleL (const std::vector<T2> &factor)
{
planck_assert(factor.size()>tsize(lmax),
"alm.ScaleL: factor array too short");
for (int m=0; m<=mmax; ++m)
for (int l=m; l<=lmax; ++l)
operator()(l,m)*=factor[l];
}
/*! \a a(l,m) *= \a factor[m] for all \a l,m. */
template<typename T2> void ScaleM (const arr<T2> &factor)
{
......
......@@ -65,7 +65,7 @@ template<typename T> void alm2map_cxx (paramfile &params)
}
arr<double> temp, pol;
get_pixwin (params,nlmax,nside,temp,pol);
get_pixwin (params,nlmax,temp,pol);
bool deriv = params.template find<bool>("derivatives",false);
if (deriv)
......
......@@ -27,7 +27,7 @@
/*! \file alm_fitsio.h
* FITS I/O for spherical harmonic coefficients
*
* Copyright (C) 2003-2010 Max-Planck-Society
* Copyright (C) 2003-2017 Max-Planck-Society
* \author Martin Reinecke
*/
......@@ -61,6 +61,16 @@ void get_almsize_pol(const std::string &filename, int &lmax, int &mmax);
the requested (l,m) range are ignored. */
template<typename T> void read_Alm_from_fits
(fitshandle &inp, Alm<xcomplex<T> > &alms, int lmax, int mmax);
/*! Returns the a_lm of the FITS binary table pointed to by \a inp.
Values not present in the FITS table are set to zero; values outside
the requested (l,m) range are ignored. */
template<typename T> inline Alm<xcomplex<T> > read_Alm_from_fits
(fitshandle &inp, int lmax, int mmax)
{
Alm<xcomplex<T> > res;
read_Alm_from_fits (inp, res, lmax, mmax);
return res;
}
/*! Opens the FITS file \a filename, jumps to the HDU \a hdunum, then reads
the a_lm from the FITS binary table there into \a alms. \a alms is
reallocated with the parameters \a lmax and \a mmax.
......@@ -69,6 +79,17 @@ template<typename T> void read_Alm_from_fits
template<typename T> void read_Alm_from_fits
(const std::string &filename, Alm<xcomplex<T> > &alms,
int lmax, int mmax, int hdunum=2);
/*! Opens the FITS file \a filename, jumps to the HDU \a hdunum, then returns
the a_lm from the FITS binary table there.
Values not present in the FITS table are set to zero; values outside
the requested \a (l,m) range are ignored. */
template<typename T> inline Alm<xcomplex<T> > read_Alm_from_fits
(const std::string &filename, int lmax, int mmax, int hdunum=2)
{
Alm<xcomplex<T> > res;
read_Alm_from_fits (filename, res, lmax, mmax, hdunum);
return res;
}
template<typename T> inline void read_Alm_from_fits
(const std::string &filename, Alm<xcomplex<T> > &almT,
......
......@@ -25,7 +25,7 @@
*/
/*
* Copyright (C) 2003-2015 Max-Planck-Society
* Copyright (C) 2003-2017 Max-Planck-Society
* Author: Martin Reinecke
*/
......@@ -70,6 +70,25 @@ template void map2alm (const Healpix_Map<double> &map,
Alm<xcomplex<double> > &alm, const arr<double> &weight,
bool add_alm);
template<typename T> void alm2map_adjoint (const Healpix_Map<T> &map,
Alm<xcomplex<T> > &alm, bool add_alm)
{
planck_assert (map.Scheme()==RING,
"alm2map_adjoint: map must be in RING scheme");
planck_assert (map.fullyDefined(),"map contains undefined pixels");
checkLmaxNside(alm.Lmax(), map.Nside());
sharp_cxxjob<T> job;
job.set_Healpix_geometry (map.Nside());
job.set_triangular_alm_info (alm.Lmax(), alm.Mmax());
job.alm2map_adjoint(&map[0], &alm(0,0), add_alm);
}
template void alm2map_adjoint (const Healpix_Map<float> &map,
Alm<xcomplex<float> > &alm, bool add_alm);
template void alm2map_adjoint (const Healpix_Map<double> &map,
Alm<xcomplex<double> > &alm, bool add_alm);
template<typename T> void map2alm_iter (const Healpix_Map<T> &map,
Alm<xcomplex<T> > &alm, int num_iter, const arr<double> &weight)
{
......@@ -152,6 +171,37 @@ template void map2alm_spin
Alm<xcomplex<double> > &alm1, Alm<xcomplex<double> > &alm2,
int spin, const arr<double> &weight, bool add_alm);
template<typename T> void alm2map_spin_adjoint
(const Healpix_Map<T> &map1, const Healpix_Map<T> &map2,
Alm<xcomplex<T> > &alm1, Alm<xcomplex<T> > &alm2,
int spin, bool add_alm)
{
planck_assert (map1.Scheme()==RING,
"alm2map_spin_adjoint: maps must be in RING scheme");
planck_assert (map1.conformable(map2),
"alm2map_spin_adjoint: maps are not conformable");
planck_assert (alm1.conformable(alm1),
"alm2map_spin_adjoint: a_lm are not conformable");
planck_assert (map1.fullyDefined()&&map2.fullyDefined(),
"map contains undefined pixels");
checkLmaxNside(alm1.Lmax(), map1.Nside());
sharp_cxxjob<T> job;
job.set_Healpix_geometry (map1.Nside());
job.set_triangular_alm_info (alm1.Lmax(), alm1.Mmax());
job.alm2map_spin_adjoint(&map1[0],&map2[0],&alm1(0,0),&alm2(0,0),spin,
add_alm);
}
template void alm2map_spin_adjoint
(const Healpix_Map<float> &map1, const Healpix_Map<float> &map2,
Alm<xcomplex<float> > &alm1, Alm<xcomplex<float> > &alm2,
int spin, bool add_alm);
template void alm2map_spin_adjoint
(const Healpix_Map<double> &map1, const Healpix_Map<double> &map2,
Alm<xcomplex<double> > &alm1, Alm<xcomplex<double> > &alm2,
int spin, bool add_alm);
template<typename T> void map2alm_spin_iter2
(const Healpix_Map<T> &map1, const Healpix_Map<T> &map2,
Alm<xcomplex<T> > &alm1, Alm<xcomplex<T> > &alm2,
......@@ -235,6 +285,48 @@ template void map2alm_pol
const arr<double> &weight,
bool add_alm);
template<typename T> void alm2map_pol_adjoint
(const Healpix_Map<T> &mapT,
const Healpix_Map<T> &mapQ,
const Healpix_Map<T> &mapU,
Alm<xcomplex<T> > &almT,
Alm<xcomplex<T> > &almG,
Alm<xcomplex<T> > &almC,
bool add_alm)
{
planck_assert (mapT.Scheme()==RING,
"alm2map_pol_adjoint: maps must be in RING scheme");
planck_assert (mapT.conformable(mapQ) && mapT.conformable(mapU),
"alm2map_pol_adjoint: maps are not conformable");
planck_assert (almT.conformable(almG) && almT.conformable(almC),
"alm2map_pol_adjoint: a_lm are not conformable");
planck_assert (mapT.fullyDefined()&&mapQ.fullyDefined()&&mapU.fullyDefined(),
"map contains undefined pixels");
checkLmaxNside(almT.Lmax(), mapT.Nside());
sharp_cxxjob<T> job;
job.set_Healpix_geometry (mapT.Nside());
job.set_triangular_alm_info (almT.Lmax(), almT.Mmax());
job.alm2map_adjoint(&mapT[0], &almT(0,0), add_alm);
job.alm2map_spin_adjoint(&mapQ[0],&mapU[0],&almG(0,0),&almC(0,0),2,add_alm);
}
template void alm2map_pol_adjoint
(const Healpix_Map<float> &mapT,
const Healpix_Map<float> &mapQ,
const Healpix_Map<float> &mapU,
Alm<xcomplex<float> > &almT,
Alm<xcomplex<float> > &almG,
Alm<xcomplex<float> > &almC,
bool add_alm);
template void alm2map_pol_adjoint
(const Healpix_Map<double> &mapT,
const Healpix_Map<double> &mapQ,
const Healpix_Map<double> &mapU,
Alm<xcomplex<double> > &almT,
Alm<xcomplex<double> > &almG,
Alm<xcomplex<double> > &almC,
bool add_alm);
template<typename T> void map2alm_pol_iter
(const Healpix_Map<T> &mapT,
......@@ -332,24 +424,24 @@ template void map2alm_pol_iter2
template<typename T> void alm2map (const Alm<xcomplex<T> > &alm,
Healpix_Map<T> &map)
Healpix_Map<T> &map, bool add_map)
{
planck_assert (map.Scheme()==RING, "alm2map: map must be in RING scheme");
sharp_cxxjob<T> job;
job.set_Healpix_geometry (map.Nside());
job.set_triangular_alm_info (alm.Lmax(), alm.Mmax());
job.alm2map(&alm(0,0), &map[0], false);
job.alm2map(&alm(0,0), &map[0], add_map);
}
template void alm2map (const Alm<xcomplex<double> > &alm,
Healpix_Map<double> &map);
Healpix_Map<double> &map, bool add_map);
template void alm2map (const Alm<xcomplex<float> > &alm,
Healpix_Map<float> &map);
Healpix_Map<float> &map, bool add_map);
template<typename T> void alm2map_spin
(const Alm<xcomplex<T> > &alm1, const Alm<xcomplex<T> > &alm2,
Healpix_Map<T> &map1, Healpix_Map<T> &map2, int spin)
Healpix_Map<T> &map1, Healpix_Map<T> &map2, int spin, bool add_map)
{
planck_assert (map1.Scheme()==RING,
"alm2map_spin: maps must be in RING scheme");
......@@ -361,15 +453,15 @@ template<typename T> void alm2map_spin
sharp_cxxjob<T> job;
job.set_Healpix_geometry (map1.Nside());
job.set_triangular_alm_info (alm1.Lmax(), alm1.Mmax());
job.alm2map_spin(&alm1(0,0),&alm2(0,0),&map1[0],&map2[0],spin,false);
job.alm2map_spin(&alm1(0,0),&alm2(0,0),&map1[0],&map2[0],spin,add_map);
}
template void alm2map_spin
(const Alm<xcomplex<double> > &alm1, const Alm<xcomplex<double> > &alm2,
Healpix_Map<double> &map, Healpix_Map<double> &map2, int spin);
Healpix_Map<double> &map, Healpix_Map<double> &map2, int spin, bool add_map);
template void alm2map_spin
(const Alm<xcomplex<float> > &alm1, const Alm<xcomplex<float> > &alm2,
Healpix_Map<float> &map, Healpix_Map<float> &map2, int spin);
Healpix_Map<float> &map, Healpix_Map<float> &map2, int spin, bool add_map);
template<typename T> void alm2map_pol
......@@ -378,7 +470,8 @@ template<typename T> void alm2map_pol
const Alm<xcomplex<T> > &almC,
Healpix_Map<T> &mapT,
Healpix_Map<T> &mapQ,
Healpix_Map<T> &mapU)
Healpix_Map<T> &mapU,
bool add_map)
{
planck_assert (mapT.Scheme()==RING,
"alm2map_pol: maps must be in RING scheme");
......@@ -390,8 +483,8 @@ template<typename T> void alm2map_pol
sharp_cxxjob<T> job;
job.set_Healpix_geometry (mapT.Nside());
job.set_triangular_alm_info (almT.Lmax(), almT.Mmax());
job.alm2map(&almT(0,0), &mapT[0], false);
job.alm2map_spin(&almG(0,0), &almC(0,0), &mapQ[0], &mapU[0], 2, false);
job.alm2map(&almT(0,0), &mapT[0], add_map);
job.alm2map_spin(&almG(0,0), &almC(0,0), &mapQ[0], &mapU[0], 2, add_map);
}
template void alm2map_pol (const Alm<xcomplex<double> > &almT,
......@@ -399,14 +492,16 @@ template void alm2map_pol (const Alm<xcomplex<double> > &almT,
const Alm<xcomplex<double> > &almC,
Healpix_Map<double> &mapT,
Healpix_Map<double> &mapQ,
Healpix_Map<double> &mapU);
Healpix_Map<double> &mapU,
bool add_map);
template void alm2map_pol (const Alm<xcomplex<float> > &almT,
const Alm<xcomplex<float> > &almG,
const Alm<xcomplex<float> > &almC,
Healpix_Map<float> &mapT,
Healpix_Map<float> &mapQ,
Healpix_Map<float> &mapU);
Healpix_Map<float> &mapU,
bool add_map);
template<typename T> void alm2map_der1
......
......@@ -25,7 +25,7 @@
*/
/*! \file alm_healpix_tools.h
* Copyright (C) 2003-2011 Max-Planck-Society
* Copyright (C) 2003-2017 Max-Planck-Society
* \author Martin Reinecke
*/
......@@ -158,11 +158,23 @@ template<typename T> void map2alm_pol_iter2
determined from this object.
\param map the output map, which must have RING ordering. */
template<typename T> void alm2map (const Alm<xcomplex<T> > &alm,
Healpix_Map<T> &map);
Healpix_Map<T> &map, bool add_map=false);
/*! Adjoint of the alm2map transform.
\param map the input map, which must have RING ordering
\param alm the output a_lms. l_max and m_max of the conversion are
determined from this object. */
template<typename T> void alm2map_adjoint (const Healpix_Map<T> &map,
Alm<xcomplex<T> > &alm, bool add_alm=false);
template<typename T> void alm2map_spin
(const Alm<xcomplex<T> > &alm1, const Alm<xcomplex<T> > &alm2,
Healpix_Map<T> &map1, Healpix_Map<T> &map2, int spin);
Healpix_Map<T> &map1, Healpix_Map<T> &map2, int spin, bool add_map=false);
template<typename T> void alm2map_spin_adjoint
(const Healpix_Map<T> &map1, const Healpix_Map<T> &map2,
Alm<xcomplex<T> > &alm1, Alm<xcomplex<T> > &alm2,
int spin, bool add_alm=false);
/*! Converts a a set of polarised a_lm to a HEALPix map.
\param almT the input temperature a_lms
......@@ -177,7 +189,17 @@ template<typename T> void alm2map_pol
const Alm<xcomplex<T> > &almC,
Healpix_Map<T> &mapT,
Healpix_Map<T> &mapQ,
Healpix_Map<T> &mapU);
Healpix_Map<T> &mapU,
bool add_map=false);
template<typename T> void alm2map_pol_adjoint
(const Healpix_Map<T> &mapT,
const Healpix_Map<T> &mapQ,
const Healpix_Map<T> &mapU,
Alm<xcomplex<T> > &almT,
Alm<xcomplex<T> > &almG,
Alm<xcomplex<T> > &almC,
bool add_alm=false);
/*! Converts a a set of a_lm to a HEALPix map and its first derivatives.
\param alm the input a_lms. l_max and m_max of the conversion are
......
......@@ -25,7 +25,7 @@
*/
/*
* Copyright (C) 2003-2015 Max-Planck-Society
* Copyright (C) 2003-2017 Max-Planck-Society
* Author: Martin Reinecke
*/
......@@ -38,6 +38,7 @@
#include "healpix_map_fitsio.h"
#include "powspec.h"
#include "powspec_fitsio.h"
#include "weight_utils.h"
#include "alm_healpix_tools.h"
#include "alm_powspec_tools.h"
#include "fitshandle.h"
......@@ -61,6 +62,11 @@ template<typename T> void anafast_cxx (paramfile &params)
bool polarisation = params.template find<bool>("polarisation");
int num_iter = params.template find<int>("iter_order",0);
bool remove_mono = params.template find<bool>("remove_monopole",false);
bool do_pwgt = params.param_present("pixelweights");
planck_assert(!(params.param_present("ringweights")&&do_pwgt),
"both pixel and ring weights specified");
planck_assert((num_iter==0)||(!do_pwgt),
"iterative analysis cannot be done in combination with pixel weights");
if (!polarisation)
{
......@@ -79,6 +85,13 @@ template<typename T> void anafast_cxx (paramfile &params)
map.Add(T(-avg));
}
if (do_pwgt)
{
auto pwgt = read_fullweights_from_fits(
params.template find<string>("pixelweights"), map.Nside());
apply_fullweights(map,pwgt);
}
arr<double> weight;
get_ring_weights (params,map.Nside(),weight);
......@@ -115,6 +128,15 @@ template<typename T> void anafast_cxx (paramfile &params)
mapT.Add(T(-avg));
}
if (do_pwgt)
{
auto pwgt = read_fullweights_from_fits(
params.template find<string>("pixelweights"), mapT.Nside());
apply_fullweights(mapT,pwgt);
apply_fullweights(mapQ,pwgt);
apply_fullweights(mapU,pwgt);
}
arr<double> weight;
get_ring_weights (params,mapT.Nside(),weight);
......
#include "levels_facilities.h"
#include "error_handling.h"
int main (int argc, const char **argv)
{
PLANCK_DIAGNOSIS_BEGIN
compute_weights_module (argc, argv);
PLANCK_DIAGNOSIS_END
}
/*
* 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.sourceforge.net
*/
/*
* 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) 2016-2018 Max-Planck-Society
* Author: Martin Reinecke
*/
#include "paramfile.h"
#include "fitshandle.h"
#include "announce.h"
#include "weight_utils.h"
using namespace std;
int compute_weights_module (int argc, const char **argv)
{
module_startup ("compute_weights", argc, argv);
paramfile params (getParamsFromCmdline(argc,argv));
planck_assert (params.param_present("outfile_ring")
|| params.param_present("outfile_full"), "no job requested");
int nside = params.find<int>("nside");
int nlmax = params.find<int>("nlmax");
if (nlmax&1)
{
cout << "Warning: specified nlmax is odd. Reducing by 1" << endl;
--nlmax;
}
double epsilon = params.find<double>("epsilon",1e-7);
int itmax = params.find<int>("max_iter",10000);
if (params.param_present("outfile_ring"))
{
double epsilon_out;
vector<double> wgt=get_ringweights(nside,nlmax,epsilon,itmax,epsilon_out);
fitshandle out;
out.create (params.find<string>("outfile_ring"));
vector<fitscolumn> cols;
int repcount = 1;
cols.push_back (fitscolumn ("TEMPERATURE WEIGHTS","1",repcount,
PLANCK_FLOAT64));
cols.push_back (fitscolumn ("Q-POLARISATION WEIGHTS","1",repcount,
PLANCK_FLOAT64));
cols.push_back (fitscolumn ("U-POLARISATION WEIGHTS","1",repcount,
PLANCK_FLOAT64));
out.insert_bintab(cols);
out.set_key ("CREATOR",string("compute_weights"),
"Software creating the FITS file");
out.set_key ("NSIDE",nside,"Resolution parameter for HEALPIX");
out.set_key ("MAX_LPOL",nlmax,"Maximum multipole l used in map synthesis");
double val=*max_element(wgt.begin(),wgt.end());
out.set_key("MAXVAL1",val,"maximum value of T weights");
out.set_key("MAXVAL2",val,"maximum value of Q weights");
out.set_key("MAXVAL3",val,"maximum value of U weights");
val=*min_element(wgt.begin(),wgt.end());
out.set_key("MINVAL1",val,"minimum value of T weights");
out.set_key("MINVAL2",val,"minimum value of Q weights");
out.set_key("MINVAL3",val,"minimum value of U weights");
out.set_key("EPSILON",epsilon_out,"epsilon reached after minimization");
out.write_column(1,wgt);
out.write_column(2,wgt);
out.write_column(3,wgt);
}
if (params.param_present("outfile_full"))
{
double epsilon_out;
vector<double> wgt=get_fullweights(nside,nlmax,epsilon,itmax,epsilon_out);
fitshandle out;
out.create (params.find<string>("outfile_full"));
vector<fitscolumn> cols;
int repcount = 1;
cols.push_back (fitscolumn ("COMPRESSED PIXEL WEIGHTS","1",repcount,
PLANCK_FLOAT64));
out.insert_bintab(cols);
out.set_key ("CREATOR",string("compute_weights"),
"Software creating the FITS file");
out.set_key ("NSIDE",nside,"Resolution parameter for HEALPix");
out.set_key ("MAX_LPOL",nlmax,"Maximum l multipole");
double val=*max_element(wgt.begin(),wgt.end());
out.set_key("MAXVAL1",val,"maximum value of pixel weights");
val=*min_element(wgt.begin(),wgt.end());
out.set_key("MINVAL1",val,"minimum value of pixel weights");
out.set_key("EPSILON",epsilon_out,"epsilon reached after minimization");
out.write_column(1,wgt);
}
return 0;
}
......@@ -25,7 +25,7 @@
*/
/*
* Copyright (C) 2003-2015 Max-Planck-Society
* Copyright (C) 2003-2016 Max-Planck-Society
* Author: Martin Reinecke
*/
......@@ -587,6 +587,8 @@ template<typename I> void T_Healpix_Base<I>::query_multidisc_general
}
}
#ifndef __BMI2__
template<> inline int T_Healpix_Base<int>::spread_bits (int v) const
{ return utab[v&0xff] | (utab[(v>>8)&0xff]<<16); }
template<> inline int64 T_Healpix_Base<int64>::spread_bits (int v) const
......@@ -608,6 +610,21 @@ template<> inline int T_Healpix_Base<int64>::compress_bits (int64 v) const
| (ctab[(raw>>32)&0xff]<<16) | (ctab[(raw>>40)&0xff]<<20);
}
#else
#include <x86intrin.h>
template<> inline int T_Healpix_Base<int>::spread_bits (int v) const
{ return _pdep_u32(v,0x55555555u); }
template<> inline int64 T_Healpix_Base<int64>::spread_bits (int v) const
{ return _pdep_u64(v,0x5555555555555555ull); }
template<> inline int T_Healpix_Base<int>::compress_bits (int v) const
{ return _pext_u32(v,0x55555555u); }
template<> inline int T_Healpix_Base<int64>::compress_bits (int64 v) const
{ return _pext_u64(v,0x5555555555555555ull); }
#endif
template<typename I> inline void T_Healpix_Base<I>::nest2xyf (I pix, int &ix,
int &iy, int &face_num) const
{
......@@ -621,6 +638,18 @@ template<typename I> inline I T_Healpix_Base<I>::xyf2nest (int ix, int iy,
int face_num) const
{ return (I(face_num)<<(2*order_)) + spread_bits(ix) + (spread_bits(iy)<<1); }
namespace {
// low-level hack to accelerate divisions with a result of [0;3]
template<typename I> inline I special_div(I a, I b)
{
I t=(a>=(b<<1));
a-=t*(b<<1);
return (t<<1)+(a>=b);
}
} // unnamed namespace
template<typename I> void T_Healpix_Base<I>::ring2xyf (I pix, int &ix, int &iy,
int &face_num) const
{
......@@ -633,7 +662,7 @@ template<typename I> void T_Healpix_Base<I>::ring2xyf (I pix, int &ix, int &iy,
iphi = (pix+1) - 2*iring*(iring-1);
kshift = 0;
nr = iring;
face_num=(iphi-1)/nr;
face_num=special_div(iphi-1,nr);
}
else if (pix<(npix_-ncap_)) // Equatorial region
{
......@@ -643,10 +672,10 @@ template<typename I> void T_Healpix_Base<I>::ring2xyf (I pix, int &ix, int &iy,
iphi = ip-tmp*4*nside_ + 1;
kshift = (iring+nside_)&1;
nr = nside_;
I ire = iring-nside_+1,
irm = nl2+2-ire;
I ifm = iphi - ire/2 + nside_ -1,
ifp = iphi - irm/2 + nside_ -1;
I ire = tmp+1,
irm = nl2+1-tmp;
I ifm = iphi - (ire>>1) + nside_ -1,
ifp = iphi - (irm>>1) + nside_ -1;
if (order_>=0)
{ ifm >>= order_; ifp >>= order_; }
else
......@@ -661,11 +690,12 @@ template<typename I> void T_Healpix_Base<I>::ring2xyf (I pix, int &ix, int &iy,
kshift = 0;
nr = iring;
iring = 2*nl2-iring;
face_num = 8 + (iphi-1)/nr;
face_num=special_div(iphi-1,nr)+8;
}
I irt = iring - (jrll[face_num]*nside_) + 1;
I irt = iring - ((2+(face_num>>2))*nside_) + 1;
I ipt = 2*iphi- jpll[face_num]*nr - kshift -1;
// I ipt = 2*iphi- (((face_num&3)<<1)+1-((face_num>>2)&1))*nr - kshift -1;
if (ipt>=nl2) ipt-=8*nside_;
ix = (ipt-irt) >>1;
......@@ -958,7 +988,7 @@ template<typename I> template<typename I2>
tsize nv=vertex.size();
tsize ncirc = inclusive ? nv+1 : nv;
planck_assert(nv>=3,"not enough vertices in polygon");
arr<vec3> vv(nv);
vector<vec3> vv(nv);
for (tsize i=0; i<nv; ++i)
vv[i]=vertex[i].to_vec3();
arr<vec3> normal(ncirc);
......
......@@ -25,7 +25,7 @@
*/
/*! \file healpix_base.h
* Copyright (C) 2003-2014 Max-Planck-Society
* Copyright (C) 2003-2017 Max-Planck-Society
* \author Martin Reinecke
*/
......@@ -228,6 +228,17 @@ template<typename I> class T_Healpix_Base: public Healpix_Tables
whose centers lie inside the disk
\note This method is more efficient in the RING scheme. */
void query_disc (pointing ptg, double radius, rangeset<I> &pixset) const;
/*! Returns the range set of all pixels whose centers lie within the disk
defined by \a dir and \a radius.
\param dir the angular coordinates of the disk center
\param radius the radius (in radians) of the disk
\note This method is more efficient in the RING scheme. */
rangeset<I> query_disc (pointing ptg, double radius) const
{
rangeset<I> res;
query_disc(ptg, radius, res);
return res;
}
/*! Returns the range set of all pixels which overlap with the disk
defined by \a dir and \a radius.
\param dir the angular coordinates of the disk center
......@@ -243,6 +254,24 @@ template<typename I> class T_Healpix_Base: public Healpix_Tables
\note This method is more efficient in the RING scheme. */
void query_disc_inclusive (pointing ptg, double radius, rangeset<I> &pixset,
int fact=1) const;
/*! Returns the range set of all pixels which overlap with the disk
defined by \a dir and \a radius.
\param dir the angular coordinates of the disk center
\param radius the radius (in radians) of the disk
\param fact The overlapping test will be done at the resolution
\a fact*nside. For NESTED ordering, \a fact must be a power of 2,
else it can be any positive integer. A typical choice would be 4.
\note This method may return some pixels which don't overlap with
the disk at all. The higher \a fact is chosen, the fewer false
positives are returned, at the cost of increased run time.
\note This method is more efficient in the RING scheme. */
rangeset<I> query_disc_inclusive (pointing ptg, double radius,
int fact=1) const