Commit ea0ecc6b authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

update last healpix

parent ad938765
......@@ -27,7 +27,7 @@
/*! \file alm.h
* Class for storing spherical harmonic coefficients.
*
* Copyright (C) 2003 - 2009 Max-Planck-Society
* Copyright (C) 2003-2011 Max-Planck-Society
* \author Martin Reinecke
*/
......@@ -46,11 +46,8 @@ class Alm_Base
public:
/*! 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);
}
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) {}
......@@ -83,12 +80,7 @@ class Alm_Base
{ 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);
}
void swap (Alm_Base &other);
};
/*! Class for storing spherical harmonic coefficients. */
......
#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
}
Parameters read by alm2map_cxx:
nlmax (integer):
maximum order of l
nmmax (integer):
maximum order of m (must not be larger than nlmax, default=nlmax)
infile (string):
input file containing the a_lm
outfile (string):
output file name for the Healpix map(s)
nside (integer):
nside parameter for the output map(s)
polarisation (bool):
if false, only the intensity map is generated,
if true, maps for I, Q and U are generated
fwhm_arcmin (double, default=0):
FWHM in arc minutes of a Gaussian beam, which is used to smooth
the a_lm
pixel_window (bool, default=false):
if true, the a_lm are multiplied by the appropriate pixel window function
if (pixel_window)
healpix_data (string):
directory containing the Healpix data files
endif
double_precision (bool, default=false):
if false, a_lm and maps are read/written in single precision,
otherwise in double precision.
/*
* 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;
}
/*
* 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 <string>
#include "alm_fitsio.h"
#include "alm.h"
#include "fitshandle.h"
#include "xcomplex.h"
#include "safe_cast.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) return;
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).Set (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).re; im[i] = alms(l,m).im; }
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 (alms(l,m).norm()>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 (alms(l,m).norm()==0)
{
++m;
if ((m>l) || (m>mmax)) { ++l; m=0; }
}
index[i] = l*l + l + m + 1;
re[i] = alms(l,m).re;
im[i] = alms(l,m).im;
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.
*
* 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_fitsio.h
* FITS I/O for spherical harmonic coefficients
*
* Copyright (C) 2003-2010 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef PLANCK_ALM_FITSIO_H
#define PLANCK_ALM_FITSIO_H
#include <string>
#include "xcomplex.h"
#include "datatypes.h"
#include "fitshandle.h"
template<typename T> class Alm;
/*! \defgroup alm_fitsio_group FITS-based I/O of a_lm */
/*! \{ */
/*! Returns the maximum \a l and \a m multipole moments found in the FITS HDU
pointed to be \a inp in \a lmax and \a mmax. */
void get_almsize(fitshandle &inp, int &lmax, int &mmax);
/*! Returns the maximum \a l and \a m multipole moments found in the HDU
\a hdunum of file \a filename in \a lmax and \a mmax. */
void get_almsize(const std::string &filename, int &lmax, int &mmax,
int hdunum=2);
/*! Returns the maximum \a l and \a m multipole moments found in the HDUs
2, 3 and 4 of file \a filename in \a lmax and \a mmax. */
void get_almsize_pol(const std::string &filename, int &lmax, int &mmax);
/*! Reads the a_lm of the FITS binary table pointed to by \a inp into
\a alms. \a alms is reallocated with the parameters \a lmax and \a mmax.
Values not present in the FITS table are set to zero; values outside
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);
/*! 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.
Values not present in the FITS table are set to zero; values outside
the requested \a (l,m) range are ignored. */
template<typename T> void read_Alm_from_fits
(const std::string &filename, Alm<xcomplex<T> > &alms,
int lmax, int mmax, int hdunum=2);
template<typename T> inline void read_Alm_from_fits
(const std::string &filename, Alm<xcomplex<T> > &almT,
Alm<xcomplex<T> > &almG, Alm<xcomplex<T> > &almC,
int lmax, int mmax, int firsthdu=2)
{
read_Alm_from_fits (filename, almT, lmax, mmax, firsthdu);
read_Alm_from_fits (filename, almG, lmax, mmax, firsthdu+1);
read_Alm_from_fits (filename, almC, lmax, mmax, firsthdu+2);
}
/*! Inserts a new binary table into \a out, which contains three columns
of type PLANCK_INT32, \a datatype and \a datatype, respectively.
The data in \a alms is written into this table; values outside
the requested (\a lmax, \a mmax) range are omitted. */
template<typename T> void write_Alm_to_fits
(fitshandle &out, const Alm<xcomplex<T> > &alms,
int lmax, int mmax, PDT datatype);
template<typename T> inline void write_Alm_to_fits
(const std::string &outfile, const Alm<xcomplex<T> > &alms,
int lmax, int mmax, PDT datatype)
{
fitshandle out;
out.create(outfile);
write_Alm_to_fits (out, alms, lmax, mmax, datatype);
}
template<typename T> inline void write_Alm_to_fits
(const std::string &outfile, const Alm<xcomplex<T> > &almT,
const Alm<xcomplex<T> > &almG, const Alm<xcomplex<T> > &almC,
int lmax, int mmax, PDT datatype)
{
fitshandle out;
out.create(outfile);
write_Alm_to_fits (out, almT, lmax, mmax, datatype);
write_Alm_to_fits (out, almG, lmax, mmax, datatype);
write_Alm_to_fits (out, almC, lmax, mmax, datatype);
}
/*! Inserts a new binary table into \a out, which contains three columns
of type PLANCK_INT32, \a datatype and \a datatype, respectively.
The data in \a alms is written into this table; values outside
the requested (\a lmax, \a mmax) range are omitted. Values with an absolute
magnitude of zero are not written. */
template<typename T> void write_compressed_Alm_to_fits
(fitshandle &out, const Alm<xcomplex<T> > &alms,
int lmax, int mmax, PDT datatype);
/*! \} */
#endif
......@@ -25,7 +25,7 @@
*/
/*
* Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008 Max-Planck-Society
* Copyright (C) 2003-2011 Max-Planck-Society
* Author: Martin Reinecke
*/
......@@ -99,7 +99,6 @@ template<typename T> void map2alm_iter2 (const Healpix_Map<T> &map,