Skip to content

GitLab

Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

Commit 1aa27558 authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

new healpix

parent c0558e69
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 937 additions and 1669 deletions
healpix/Healpix_cxx/Healpix_cxx.dox
View file @ 1aa27558
......@@ -125,9 +125,9 @@ Transform 1: Equatorial (2000) -> Galactic (2000)
\verbatim Usage: calc_powspec <almfile1> [<almfile2>] <powspec_file>\endverbatim
\section median_filter
\section median_filter_cxx
This program inputs a HEALPix map, runs a median filter with the desired
radius on it and saves the result to another file.
\verbatim Usage: median_filter <input map> <output map> <radius in arcmin>\endverbatim
\verbatim Usage: median_filter_cxx <input map> <output map> <radius in arcmin>\endverbatim
*/
healpix/Healpix_cxx/alm.h
View file @ 1aa27558
......@@ -27,7 +27,7 @@
/*! \file alm.h
* Class for storing spherical harmonic coefficients.
*
* Copyright (C) 2003, 2004, 2005, 2006, 2007 Max-Planck-Society
* Copyright (C) 2003 - 2009 Max-Planck-Society
* \author Martin Reinecke
*/
......@@ -36,35 +36,77 @@
#include "arr.h"
/*! Class for storing spherical harmonic coefficients. */
template<typename T> class Alm
/*! Base class for calculating the storage layout of spherical harmonic
coefficients. */
class Alm_Base
{
private:
protected:
int lmax, mmax, tval;
arr<T> alm;
public:
/*! Returns the number of coefficients in an Alm object with maximum
quantum numbers \a l and \a m. */
static long Num_Alms (int l, int m)
/*! Returns the total number of coefficients for maximum quantum numbers
\a l and \a m. */
static tsize 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);
}
/*! Constructs an Alm_Base object with given \a lmax and \a mmax. */
Alm_Base (int lmax_=0, int mmax_=0)
: lmax(lmax_), mmax(mmax_), tval(2*lmax+1) {}
/*! Changes the object's maximum quantum numbers to \a lmax and \a mmax. */
void Set (int lmax_, int mmax_)
{
lmax=lmax_;
mmax=mmax_;
tval=2*lmax+1;
}
/*! Returns the maximum \a l */
int Lmax() const { return lmax; }
/*! Returns the maximum \a m */
int Mmax() const { return mmax; }
/*! Returns an array index for a given m, from which the index of a_lm
can be obtained by adding l. */
int index_l0 (int m) const
{ return ((m*(tval-m))>>1); }
/*! Returns the array index of the specified coefficient. */
int index (int l, int m) const
{ return index_l0(m) + l; }
/*! Returns \a true, if both objects have the same \a lmax and \a mmax,
else \a false. */
bool conformable (const Alm_Base &other) const
{ return ((lmax==other.lmax) && (mmax==other.mmax)); }
/*! Swaps the contents of two Alm_Base objects. */
void swap (Alm_Base &other)
{
std::swap(lmax, other.lmax);
std::swap(mmax, other.mmax);
std::swap(tval, other.tval);
}
};
/*! Class for storing spherical harmonic coefficients. */
template<typename T> class Alm: public Alm_Base
{
private:
arr<T> alm;
public:
/*! Constructs an Alm object with given \a lmax and \a mmax. */
Alm (int lmax_=0, int mmax_=0)
: lmax(lmax_), mmax(mmax_), tval(2*lmax+1),
alm (Num_Alms(lmax,mmax))
{}
: Alm_Base(lmax_,mmax_), alm (Num_Alms(lmax,mmax)) {}
/*! Deletes the old coefficients and allocates storage according to
\a lmax and \a mmax. */
void Set (int lmax_, int mmax_)
{
lmax=lmax_;
mmax=mmax_;
tval=2*lmax+1;
Alm_Base::Set(lmax_, mmax_);
alm.alloc(Num_Alms(lmax,mmax));
}
......@@ -73,9 +115,7 @@ template<typename T> class Alm
void Set (arr<T> &data, int lmax_, int mmax_)
{
planck_assert (Num_Alms(lmax_,mmax_)==data.size(),"wrong array size");
lmax=lmax_;
mmax=mmax_;
tval=2*lmax+1;
Alm_Base::Set(lmax_, mmax_);
alm.transfer(data);
}
......@@ -85,11 +125,12 @@ template<typename T> class Alm
/*! Multiplies all coefficients by \a factor. */
template<typename T2> void Scale (const T2 &factor)
{ for (int m=0; m<alm.size(); ++m) alm[m]*=factor; }
{ for (tsize m=0; m<alm.size(); ++m) alm[m]*=factor; }
/*! \a a(l,m) *= \a factor[l] for all \a l,m. */
template<typename T2> void ScaleL (const arr<T2> &factor)
{
planck_assert(factor.size()>lmax, "alm.ScaleL: factor array too short");
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];
......@@ -100,47 +141,35 @@ template<typename T> class Alm
/*! Returns a reference to the specified coefficient. */
T &operator() (int l, int m)
{ return alm[((m*(tval-m))>>1) + l]; }
{ return alm[index(l,m)]; }
/*! Returns a constant reference to the specified coefficient. */
const T &operator() (int l, int m) const
{ return alm[((m*(tval-m))>>1) + l]; }
{ return alm[index(l,m)]; }
/*! Returns a pointer for a given m, from which the address of a_lm
can be obtained by adding l. */
T *mstart (int m)
{ return &alm[(m*(tval-m))>>1]; }
{ return &alm[index_l0(m)]; }
/*! Returns a pointer for a given m, from which the address of a_lm
can be obtained by adding l. */
const T *mstart (int m) const
{ return &alm[(m*(tval-m))>>1]; }
/*! Returns the maximum \a l */
int Lmax() const { return lmax; }
/*! Returns the maximum \a m */
int Mmax() const { return mmax; }
{ return &alm[index_l0(m)]; }
/*! Returns a constant reference to the a_lm data. */
const arr<T> &Alms () const { return alm; }
/*! Swaps the contents of two A_lm objects. */
/*! Swaps the contents of two Alm objects. */
void swap (Alm &other)
{
std::swap(lmax, other.lmax);
std::swap(mmax, other.mmax);
std::swap(tval, other.tval);
Alm_Base::swap(other);
alm.swap(other.alm);
}
/*! Returns \a true, if both objects have the same \a lmax and \a mmax,
else \a false. */
bool conformable (const Alm &other) const
{ return ((lmax==other.lmax) && (mmax==other.mmax)); }
/*! Adds all coefficients from \a other to the own coefficients. */
void Add (const Alm &other)
{
planck_assert (conformable(other), "A_lm are not conformable");
for (int m=0; m<alm.size(); ++m)
for (tsize m=0; m<alm.size(); ++m)
alm[m] += other.alm[m];
}
};
......
healpix/Healpix_cxx/alm2map_cxx.cc
View file @ 1aa27558
/*
* 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) 2003, 2004, 2005 Max-Planck-Society
* Author: Martin Reinecke
*/
#include "xcomplex.h"
#include "cxxutils.h"
#include "paramfile.h"
#include "simparams.h"
#include "healpix_data_io.h"
#include "alm.h"
#include "alm_fitsio.h"
#include "healpix_map.h"
#include "healpix_map_fitsio.h"
#include "alm_healpix_tools.h"
#include "alm_powspec_tools.h"
#include "fitshandle.h"
using namespace std;
template<typename T> void alm2map_cxx (paramfile &params, simparams &par)
{
par.add_comment("-----------------------");
par.add_comment(" ** alm2map_cxx 1.0 **");
par.add_comment("-----------------------");
int nlmax = params.template find<int>("nlmax");
par.add("nlmax","NLMAX",nlmax,"maximum l of the alms");
int nmmax = params.template find<int>("nmmax",nlmax);
par.add("nmmax","NMMAX",nmmax,"maximum m of the alms");
planck_assert(nmmax<=nlmax,"nmmax must not be larger than nlmax");
string infile = params.template find<string>("infile");
par.add_source_file (infile, 2);
string outfile = params.template find<string>("outfile");
int nside = params.template find<int>("nside");
double fwhm_arcmin = params.template find<double>("fwhm_arcmin",0);
arr<double> temp, pol;
get_pixwin (params,par,nlmax,nside,temp,pol);
bool deriv = params.template find<bool>("derivatives",false);
if (deriv)
{
Alm<xcomplex<T> > alm;
read_Alm_from_fits(infile,alm,nlmax,nmmax,2);
if (fwhm_arcmin>0) smooth_with_Gauss (alm, fwhm_arcmin);
Healpix_Map<T> map(nside,RING,SET_NSIDE),
mapdth(nside,RING,SET_NSIDE),
mapdph(nside,RING,SET_NSIDE);
alm.ScaleL(temp);
double offset = alm(0,0).real()/sqrt(fourpi);
alm(0,0) = 0;
alm2map_der1(alm,map,mapdth,mapdph);
for (int m=0; m<map.Npix(); ++m) map[m]+=offset;
fitshandle out;
out.create (outfile);
write_Healpix_map_to_fits (out,map,mapdth,mapdph,FITSUTIL<T>::DTYPE);
par.add_keys(out);
return;
}
bool polarisation = params.template find<bool>("polarisation");
if (!polarisation)
{
Alm<xcomplex<T> > alm;
read_Alm_from_fits(infile,alm,nlmax,nmmax,2);
if (fwhm_arcmin>0) smooth_with_Gauss (alm, fwhm_arcmin);
Healpix_Map<T> map(nside,RING,SET_NSIDE);
alm.ScaleL(temp);
double offset = alm(0,0).real()/sqrt(fourpi);
alm(0,0) = 0;
alm2map(alm,map);
map.add(offset);
fitshandle out;
out.create (outfile);
write_Healpix_map_to_fits (out,map,FITSUTIL<T>::DTYPE);
par.add_keys(out);
}
else
{
Alm<xcomplex<T> > almT, almG, almC;
read_Alm_from_fits(infile,almT,nlmax,nmmax,2);
read_Alm_from_fits(infile,almG,nlmax,nmmax,3);
read_Alm_from_fits(infile,almC,nlmax,nmmax,4);
if (fwhm_arcmin>0) smooth_with_Gauss (almT, almG, almC, fwhm_arcmin);
Healpix_Map<T> mapT(nside,RING,SET_NSIDE), mapQ(nside,RING,SET_NSIDE),
mapU(nside,RING,SET_NSIDE);
almT.ScaleL(temp);
almG.ScaleL(pol); almC.ScaleL(pol);
double offset = almT(0,0).real()/sqrt(fourpi);
almT(0,0) = 0;
alm2map_pol(almT,almG,almC,mapT,mapQ,mapU);
mapT.add(offset);
fitshandle out;
out.create (outfile);
write_Healpix_map_to_fits (out,mapT,mapQ,mapU,FITSUTIL<T>::DTYPE);
par.add_keys(out);
}
}
#include "levels_facilities.h"
#include "error_handling.h"
int main (int argc, const char **argv)
{
PLANCK_DIAGNOSIS_BEGIN
module_startup ("alm2map_cxx", argc, argv, 2, "<parameter file>");
paramfile params (argv[1]);
simparams par;
bool dp = params.find<bool> ("double_precision",false);
dp ? alm2map_cxx<double>(params,par) : alm2map_cxx<float>(params,par);
PLANCK_DIAGNOSIS_END
PLANCK_DIAGNOSIS_BEGIN
alm2map_cxx_module (argc, argv);
PLANCK_DIAGNOSIS_END
}
healpix/Healpix_cxx/alm2map_cxx_module.cc 0 → 100644
View file @ 1aa27558
/*
* 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) 2003-2010 Max-Planck-Society
* Author: Martin Reinecke
*/
#include "xcomplex.h"
#include "cxxutils.h"
#include "paramfile.h"
#include "healpix_data_io.h"
#include "alm.h"
#include "alm_fitsio.h"
#include "healpix_map.h"
#include "healpix_map_fitsio.h"
#include "alm_healpix_tools.h"
#include "alm_powspec_tools.h"
#include "fitshandle.h"
#include "levels_facilities.h"
#include "lsconstants.h"
using namespace std;
namespace {
template<typename T> void alm2map_cxx (paramfile &params)
{
int nlmax = params.template find<int>("nlmax");
int nmmax = params.template find<int>("nmmax",nlmax);
planck_assert(nmmax<=nlmax,"nmmax must not be larger than nlmax");
string infile = params.template find<string>("infile");
string outfile = params.template find<string>("outfile");
int nside = params.template find<int>("nside");
double fwhm = arcmin2rad*params.template find<double>("fwhm_arcmin",0);
arr<double> temp, pol;
get_pixwin (params,nlmax,nside,temp,pol);
bool deriv = params.template find<bool>("derivatives",false);
if (deriv)
{
Alm<xcomplex<T> > alm;
read_Alm_from_fits(infile,alm,nlmax,nmmax,2);
if (fwhm>0) smoothWithGauss (alm, fwhm);
Healpix_Map<T> map(nside,RING,SET_NSIDE),
mapdth(nside,RING,SET_NSIDE),
mapdph(nside,RING,SET_NSIDE);
alm.ScaleL(temp);
double offset = alm(0,0).real()/sqrt(fourpi);
alm(0,0) = 0;
alm2map_der1(alm,map,mapdth,mapdph);
map.Add(T(offset));
write_Healpix_map_to_fits (outfile,map,mapdth,mapdph,planckType<T>());
return;
}
bool polarisation = params.template find<bool>("polarisation");
if (!polarisation)
{
Alm<xcomplex<T> > alm;
read_Alm_from_fits(infile,alm,nlmax,nmmax,2);
if (fwhm>0) smoothWithGauss (alm, fwhm);
Healpix_Map<T> map(nside,RING,SET_NSIDE);
alm.ScaleL(temp);
double offset = alm(0,0).real()/sqrt(fourpi);
alm(0,0) = 0;
alm2map(alm,map);
map.Add(T(offset));
write_Healpix_map_to_fits (outfile,map,planckType<T>());
}
else
{
Alm<xcomplex<T> > almT, almG, almC;
read_Alm_from_fits(infile,almT,almG,almC,nlmax,nmmax,2);
if (fwhm>0) smoothWithGauss (almT, almG, almC, fwhm);
Healpix_Map<T> mapT(nside,RING,SET_NSIDE), mapQ(nside,RING,SET_NSIDE),
mapU(nside,RING,SET_NSIDE);
almT.ScaleL(temp);
almG.ScaleL(pol); almC.ScaleL(pol);
double offset = almT(0,0).real()/sqrt(fourpi);
almT(0,0) = 0;
alm2map_pol(almT,almG,almC,mapT,mapQ,mapU);
mapT.Add(T(offset));
write_Healpix_map_to_fits (outfile,mapT,mapQ,mapU,planckType<T>());
}
}
} // unnamed namespace
int alm2map_cxx_module (int argc, const char **argv)
{
module_startup ("alm2map_cxx", argc, argv, 2, "<parameter file>");
paramfile params (argv[1]);
bool dp = params.find<bool> ("double_precision",false);
dp ? alm2map_cxx<double>(params) : alm2map_cxx<float>(params);
return 0;
}
healpix/Healpix_cxx/alm_fitsio.cc
View file @ 1aa27558
......@@ -25,7 +25,7 @@
*/
/*
* Copyright (C) 2003, 2004, 2005 Max-Planck-Society
* Copyright (C) 2003-2010 Max-Planck-Society
* Author: Martin Reinecke
*/
......@@ -34,6 +34,7 @@
#include "alm.h"
#include "fitshandle.h"
#include "xcomplex.h"
#include "safe_cast.h"
using namespace std;
......@@ -46,26 +47,23 @@ void get_almsize(fitshandle &inp, int &lmax, int &mmax)
return;
}
int n_alms = inp.nelems(1);
int n_alms = safe_cast<int>(inp.nelems(1));
arr<int> index;
const int chunksize=1024*256;
int offset=0;
lmax=-1;
mmax=-1;
while (offset<n_alms)
lmax=mmax=-1;
chunkMaker cm(n_alms,inp.efficientChunkSize(1));
uint64 offset,ppix;
while(cm.getNext(offset,ppix))
{
int ppix=min(chunksize,n_alms-offset);
index.alloc(ppix);
inp.read_column(1,index,offset);
for (int i=0; i<ppix; ++i)
for (tsize i=0; i<ppix; ++i)
{
int l = isqrt(index[i]-1);
int m = index[i] - l*l - l - 1;
if (l>lmax) lmax=l;
if (m>mmax) mmax=m;
}
offset+=chunksize;
}
}
......@@ -95,25 +93,24 @@ void get_almsize_pol(const string &filename, int &lmax, int &mmax)
template<typename T> void read_Alm_from_fits
(fitshandle &inp, Alm<xcomplex<T> >&alms, int lmax, int mmax)
{
int n_alms = inp.nelems(1);
int n_alms = safe_cast<int>(inp.nelems(1));
arr<int> index;
arr<T> re, im;
alms.Set(lmax, mmax);
alms.SetToZero();
int max_index = lmax*lmax + lmax + mmax + 1;
const int chunksize=1024*256;
int offset=0;
while (offset<n_alms)
chunkMaker cm(n_alms,inp.efficientChunkSize(1));
uint64 offset,ppix;
while(cm.getNext(offset,ppix))
{
int ppix=min(chunksize,n_alms-offset);
index.alloc(ppix);
re.alloc(ppix); im.alloc(ppix);
inp.read_column(1,index,offset);
inp.read_column(2,re,offset);
inp.read_column(3,im,offset);
for (int i=0; i<ppix; ++i)
for (tsize i=0; i<ppix; ++i)
{
if (index[i]>max_index) return;
......@@ -124,7 +121,6 @@ template<typename T> void read_Alm_from_fits
if ((l<=lmax) && (m<=mmax))
alms(l,m).Set (re[i], im[i]);
}
offset+=chunksize;
}
}
......@@ -152,10 +148,10 @@ template void read_Alm_from_fits (const string &filename,
template<typename T> void write_Alm_to_fits
(fitshandle &out, const Alm<xcomplex<T> > &alms, int lmax, int mmax,
int datatype)
PDT datatype)
{
vector<fitscolumn> cols;
cols.push_back (fitscolumn("index","l*l+l+m+1",1,TINT32BIT));
cols.push_back (fitscolumn("index","l*l+l+m+1",1,PLANCK_INT32));
cols.push_back (fitscolumn("real","unknown",1,datatype));
cols.push_back (fitscolumn("imag","unknown",1,datatype));
out.insert_bintab(cols);
......@@ -166,14 +162,13 @@ template<typename T> void write_Alm_to_fits
int n_alms = ((mmax+1)*(mmax+2))/2 + (mmax+1)*(lmax-mmax);
int l=0, m=0;
const int chunksize=1024*256;
int offset=0;
while (offset<n_alms)
chunkMaker cm(n_alms,out.efficientChunkSize(1));
uint64 offset,ppix;
while(cm.getNext(offset,ppix))
{
int ppix=min(chunksize,n_alms-offset);
index.alloc(ppix);
re.alloc(ppix); im.alloc(ppix);
for (int i=0; i<ppix; ++i)
for (tsize i=0; i<ppix; ++i)
{
index[i] = l*l + l + m + 1;
if ((l<=lm) && (m<=mm))
......@@ -186,27 +181,25 @@ template<typename T> void write_Alm_to_fits
out.write_column(1,index,offset);
out.write_column(2,re,offset);
out.write_column(3,im,offset);
offset+=chunksize;
}
out.add_key("MAX-LPOL",lmax,"highest l in the table");
out.add_key("MAX-MPOL",mmax,"highest m in the table");
out.set_key("MAX-LPOL",lmax,"highest l in the table");
out.set_key("MAX-MPOL",mmax,"highest m in the table");
}