Commit 2743f9e9 authored by Leo Singer's avatar Leo Singer

Import Upstream version 3.30.0

parents
/*
* 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).
*/
/*! \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.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).
*/
/*! \file alm.h
* Class for storing spherical harmonic coefficients.
*
* Copyright (C) 2003-2012 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef PLANCK_ALM_H
#define PLANCK_ALM_H
#include "arr.h"
/*! Base class for calculating the storage layout of spherical harmonic
coefficients. */
class Alm_Base
{
protected:
int lmax, mmax, tval;
public:
/*! Returns the total number of coefficients for maximum quantum numbers
\a l and \a m. */
static tsize Num_Alms (int l, int 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);
};
/*! 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)
: 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_)
{
Alm_Base::Set(lmax_, mmax_);
alm.alloc(Num_Alms(lmax,mmax));
}
/*! Deallocates the old coefficients and uses the content of \a data
for storage. \a data is deallocated during the call. */
void Set (arr<T> &data, int lmax_, int mmax_)
{
planck_assert (Num_Alms(lmax_,mmax_)==data.size(),"wrong array size");
Alm_Base::Set(lmax_, mmax_);
alm.transfer(data);
}
/*! Sets all coefficients to zero. */
void SetToZero ()
{ alm.fill (0); }
/*! Multiplies all coefficients by \a factor. */
template<typename T2> void Scale (const T2 &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()>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)
{
planck_assert(factor.size()>tsize(mmax),
"alm.ScaleM: factor array too short");
for (int m=0; m<=mmax; ++m)
for (int l=m; l<=lmax; ++l)
operator()(l,m)*=factor[m];
}
/*! Adds \a num to a_00. */
template<typename T2> void Add (const T2 &num)
{ alm[0]+=num; }
/*! Returns a reference to the specified coefficient. */
T &operator() (int l, int m)
{ return alm[index(l,m)]; }
/*! Returns a constant reference to the specified coefficient. */
const T &operator() (int l, int m) const
{ 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[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[index_l0(m)]; }
/*! Returns a constant reference to the a_lm data. */
const arr<T> &Alms () const { return alm; }
/*! Swaps the contents of two Alm objects. */
void swap (Alm &other)
{
Alm_Base::swap(other);
alm.swap(other.alm);
}
/*! 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 (tsize m=0; m<alm.size(); ++m)
alm[m] += other.alm[m];
}
};
#endif
#include "levels_facilities.h"
#include "error_handling.h"
int main (int argc, const char **argv)
{
PLANCK_DIAGNOSIS_BEGIN
alm2map_cxx_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) 2003-2014 Max-Planck-Society
* Author: Martin Reinecke
*/
#include "xcomplex.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"
#include "announce.h"
#include "planck_rng.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);
int cw_lmin=-1, cw_lmax=-1;
if (params.param_present("cw_lmin"))
{
cw_lmin = params.template find<int>("cw_lmin");
cw_lmax = params.template find<int>("cw_lmax");
}
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);
if (cw_lmin>=0) applyCosineWindow(alm, cw_lmin, cw_lmax);
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");
bool do_regnoise = params.param_present("regnoiseT");
if (!polarisation)
{
Alm<xcomplex<T> > alm;
read_Alm_from_fits(infile,alm,nlmax,nmmax,2);
if (fwhm>0) smoothWithGauss (alm, fwhm);
if (cw_lmin>=0) applyCosineWindow(alm, cw_lmin, cw_lmax);
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));
if (do_regnoise)
{
planck_rng rng(params.template find<int>("rand_seed",42));
double rms = params.template find<double>("regnoiseT");
for (int i=0; i<map.Npix(); ++i)
map[i] += rms*rng.rand_gauss();
}
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);
if (cw_lmin>=0) applyCosineWindow(almT, almG, almC, cw_lmin, cw_lmax);
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));
if (do_regnoise)
{
planck_rng rng(params.template find<int>("rand_seed",42));
double rmsT = params.template find<double>("regnoiseT"),
rmsQU = params.template find<double>("regnoiseQU");
for (int i=0; i<mapT.Npix(); ++i)
{
mapT[i] += rmsT *rng.rand_gauss();
mapQ[i] += rmsQU*rng.rand_gauss();
mapU[i] += rmsQU*rng.rand_gauss();
}
}
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);
paramfile params (getParamsFromCmdline(argc,argv));
bool dp = params.find<bool> ("double_precision",false);
dp ? alm2map_cxx<double>(params) : alm2map_cxx<float>(params);
return 0;
}
/*
* 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) 2003-2015 Max-Planck-Society
* Author: Martin Reinecke
*/
#include <string>
#include "alm_fitsio.h"
#include "alm.h"
#include "fitshandle.h"
#include "xcomplex.h"
#include "safe_cast.h"
#include "share_utils.h"
using namespace std;
void get_almsize(fitshandle &inp, int &lmax, int &mmax)
{
if (inp.key_present("MAX-LPOL") && inp.key_present("MAX-MPOL"))
{
inp.get_key ("MAX-LPOL",lmax);
inp.get_key ("MAX-MPOL",mmax);
return;
}
int n_alms = safe_cast<int>(inp.nelems(1));
arr<int> index;
lmax=mmax=-1;
chunkMaker cm(n_alms,inp.efficientChunkSize(1));
uint64 offset,ppix;
while(cm.getNext(offset,ppix))
{
index.alloc(ppix);
inp.read_column(1,index,offset);
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;
}
}
}
void get_almsize(const string &filename, int &lmax, int &mmax, int hdunum)
{
fitshandle inp;
inp.open (filename);
inp.goto_hdu(hdunum);
get_almsize (inp, lmax, mmax);
}
void get_almsize_pol(const string &filename, int &lmax, int &mmax)
{
int tlmax, tmmax;
fitshandle inp;
inp.open (filename);
lmax=mmax=0;
for (int hdu=2; hdu<=4; ++hdu)
{
inp.goto_hdu(hdu);
get_almsize (inp,tlmax,tmmax);
if (tlmax>lmax) lmax=tlmax;
if (tmmax>mmax) mmax=tmmax;
}
}
template<typename T> void read_Alm_from_fits
(fitshandle &inp, Alm<xcomplex<T> >&alms, int lmax, int mmax)
{
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;
chunkMaker cm(n_alms,inp.efficientChunkSize(1));
uint64 offset,ppix;
while(cm.getNext(offset,ppix))
{
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 (tsize i=0; i<ppix; ++i)
{
if (index[i]>max_index) continue;
int l = isqrt(index[i]-1);
int m = index[i] - l*l - l - 1;
planck_assert(m>=0,"negative m encountered");
planck_assert(l>=m, "wrong l,m combination");
if ((l<=lmax) && (m<=mmax))
alms(l,m) = xcomplex<T> (re[i], im[i]);
}
}
}
template void read_Alm_from_fits (fitshandle &inp,
Alm<xcomplex<double> > &alms, int lmax, int mmax);
template void read_Alm_from_fits (fitshandle &inp,
Alm<xcomplex<float> > &alms, int lmax, int mmax);
template<typename T> void read_Alm_from_fits
(const string &filename, Alm<xcomplex<T> >&alms, int lmax, int mmax,
int hdunum)
{
fitshandle inp;
inp.open (filename);
inp.goto_hdu(hdunum);
read_Alm_from_fits(inp,alms,lmax,mmax);
}
template void read_Alm_from_fits (const string &filename,
Alm<xcomplex<double> > &alms, int lmax, int mmax, int hdunum);
template void read_Alm_from_fits (const string &filename,
Alm<xcomplex<float> > &alms, int lmax, int mmax, int hdunum);
template<typename T> void write_Alm_to_fits
(fitshandle &out, const Alm<xcomplex<T> > &alms, int lmax, int mmax,
PDT datatype)
{
vector<fitscolumn> cols;
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);
arr<int> index;
arr<double> re, im;
int lm=alms.Lmax(), mm=alms.Mmax();
int n_alms = ((mmax+1)*(mmax+2))/2 + (mmax+1)*(lmax-mmax);
int l=0, m=0;
chunkMaker cm(n_alms,out.efficientChunkSize(1));
uint64 offset,ppix;
while(cm.getNext(offset,ppix))
{
index.alloc(ppix);
re.alloc(ppix); im.alloc(ppix);
for (tsize i=0; i<ppix; ++i)
{
index[i] = l*l + l + m + 1;
if ((l<=lm) && (m<=mm))
{ re[i] = alms(l,m).real(); im[i] = alms(l,m).imag(); }
else
{ re[i] = 0; im[i] = 0; }
++m;
if ((m>l) || (m>mmax)) { ++l; m=0; }
}
out.write_column(1,index,offset);
out.write_column(2,re,offset);
out.write_column(3,im,offset);
}
out.set_key("MAX-LPOL",lmax,"highest l in the table");
out.set_key("MAX-MPOL",mmax,"highest m in the table");
}
template void write_Alm_to_fits
(fitshandle &out, const Alm<xcomplex<double> > &alms, int lmax,
int mmax, PDT datatype);
template void write_Alm_to_fits
(fitshandle &out, const Alm<xcomplex<float> > &alms, int lmax,
int mmax, PDT datatype);
template<typename T> void write_compressed_Alm_to_fits
(fitshandle &out, const Alm<xcomplex<T> > &alms, int lmax, int mmax,
PDT datatype)
{
vector<fitscolumn> cols;
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);
arr<int> index;
arr<double> re, im;
int n_alms = 0;
for (int m=0; m<=mmax; ++m)
for (int l=m; l<=lmax; ++l)
if (norm(alms(l,m))>0) ++n_alms;
int l=0, m=0;
int real_lmax=0, real_mmax=0;
chunkMaker cm(n_alms,out.efficientChunkSize(1));
uint64 offset,ppix;
while(cm.getNext(offset,ppix))
{
index.alloc(ppix);
re.alloc(ppix); im.alloc(ppix);
for (tsize i=0; i<ppix; ++i)
{
while (norm(alms(l,m))==0)
{
++m;
if ((m>l) || (m>mmax)) { ++l; m=0; }
}
index[i] = l*l + l + m + 1;
re[i] = alms(l,m).real();
im[i] = alms(l,m).imag();
if (l>real_lmax) real_lmax=l;
if (m>real_mmax) real_mmax=m;
++m;
if ((m>l) || (m>mmax)) { ++l; m=0; }
}
out.write_column(1,index,offset);
out.write_column(2,re,offset);
out.write_column(3,im,offset);
}
out.set_key("MAX-LPOL",real_lmax,"highest l in the table");
out.set_key("MAX-MPOL",real_mmax,"highest m in the table");
}
template void write_compressed_Alm_to_fits
(fitshandle &out, const Alm<xcomplex<double> > &alms, int lmax,
int mmax, PDT datatype);
template void write_compressed_Alm_to_fits
(fitshandle &out, const Alm<xcomplex<float> > &alms, int lmax,
int mmax, PDT datatype);
/*
* 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.