alm_healpix_tools.cc

00001 /*
00002  *  This file is part of Healpix_cxx.
00003  *
00004  *  Healpix_cxx is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  Healpix_cxx is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with Healpix_cxx; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  *
00018  *  For more information about HEALPix, see http://healpix.jpl.nasa.gov
00019  */
00020 
00021 /*
00022  *  Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
00023  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00024  *  (DLR).
00025  */
00026 
00027 /*
00028  *  Copyright (C) 2003, 2004, 2005 Max-Planck-Society
00029  *  Author: Martin Reinecke
00030  */
00031 
00032 #include "alm_healpix_tools.h"
00033 #include "alm_map_tools.h"
00034 #include "alm.h"
00035 #include "healpix_map.h"
00036 #include "xcomplex.h"
00037 
00038 using namespace std;
00039 
00040 namespace {
00041 
00042 void healpix2ringpairs (const Healpix_Base &base,
00043   const arr<double> &weight, vector<ringpair> &pair)
00044   {
00045   pair.clear();
00046   int startpix, ringpix;
00047   double theta, wgt, phi0;
00048   bool shifted;
00049   int nside = base.Nside();
00050   for (int m=0; m<2*nside-1; ++m)
00051     {
00052     base.get_ring_info2 (m+1, startpix, ringpix, theta, shifted);
00053     wgt = weight[m]*fourpi/base.Npix();
00054     phi0 = shifted ? pi/ringpix : 0;
00055     pair.push_back (ringpair(
00056       ringinfo(theta, phi0, wgt, ringpix, startpix),
00057       ringinfo(pi-theta, phi0, wgt, ringpix, base.Npix()-startpix-ringpix)));
00058     }
00059   base.get_ring_info2 (2*nside, startpix, ringpix, theta, shifted);
00060   wgt = weight[2*nside-1]*fourpi/base.Npix();
00061   phi0 = shifted ? pi/ringpix : 0;
00062   pair.push_back (ringpair(ringinfo(theta, phi0, wgt, ringpix, startpix)));
00063   }
00064 
00065 void healpix2ringpairs (const Healpix_Base &base, vector<ringpair> &pair)
00066   {
00067   arr<double> wgt(2*base.Nside());
00068   wgt.fill(1);
00069   healpix2ringpairs(base,wgt,pair);
00070   }
00071 
00072 } // namespace
00073 
00074 
00075 template<typename T> void map2alm (const Healpix_Map<T> &map,
00076   Alm<xcomplex<T> > &alm, const arr<double> &weight, bool add_alm)
00077   {
00078   planck_assert (map.Scheme()==RING, "map2alm: map must be in RING scheme");
00079   planck_assert (weight.size()>=2*map.Nside(),
00080     "map2alm: weight array has too few entries");
00081 
00082   vector<ringpair> pair;
00083   healpix2ringpairs(map,weight,pair);
00084   map2alm(pair,&map[0],alm,add_alm);
00085   }
00086 
00087 template void map2alm (const Healpix_Map<float> &map,
00088   Alm<xcomplex<float> > &alm, const arr<double> &weight,
00089   bool add_alm);
00090 template void map2alm (const Healpix_Map<double> &map,
00091   Alm<xcomplex<double> > &alm, const arr<double> &weight,
00092   bool add_alm);
00093 
00094 template<typename T> void map2alm_iter (const Healpix_Map<T> &map,
00095   Alm<xcomplex<T> > &alm, int num_iter, const arr<double> &weight)
00096   {
00097   map2alm(map,alm,weight);
00098   for (int iter=1; iter<=num_iter; ++iter)
00099     {
00100     Healpix_Map<T> map2(map.Nside(),map.Scheme(),SET_NSIDE);
00101     alm2map(alm,map2);
00102     for (int m=0; m<map.Npix(); ++m)
00103       map2[m] = map[m]-map2[m];
00104     map2alm(map2,alm,weight,true);
00105     }
00106   }
00107 
00108 template void map2alm_iter (const Healpix_Map<float> &map,
00109   Alm<xcomplex<float> > &alm, int num_iter,
00110   const arr<double> &weight);
00111 template void map2alm_iter (const Healpix_Map<double> &map,
00112   Alm<xcomplex<double> > &alm, int num_iter,
00113   const arr<double> &weight);
00114 
00115 template<typename T> void map2alm_iter2 (const Healpix_Map<T> &map,
00116   Alm<xcomplex<T> > &alm, double err_abs, double err_rel)
00117   {
00118   double x_err_abs=1./err_abs, x_err_rel=1./err_rel;
00119   arr<double> wgt(2*map.Nside());
00120   wgt.fill(1);
00121   Healpix_Map<T> map2(map);
00122   alm.SetToZero();
00123   while(true)
00124     {
00125     map2alm(map2,alm,wgt,true);
00126     alm2map(alm,map2);
00127     double errmeasure=0;
00128     for (int m=0; m<map.Npix(); ++m)
00129       {
00130       double err = abs(map[m]-map2[m]);
00131       double rel = (map[m]!=0) ? abs(err/map[m]) : 1e300;
00132       errmeasure = max(errmeasure,min(err*x_err_abs,rel*x_err_rel));
00133       map2[m] = map[m]-map2[m];
00134       }
00135 cout << "map error measure: " << errmeasure << endl;
00136     if (errmeasure<1) break;
00137     }
00138   }
00139 
00140 template void map2alm_iter2 (const Healpix_Map<double> &map,
00141   Alm<xcomplex<double> > &alm, double err_abs, double err_rel);
00142 
00143 
00144 template<typename T> void map2alm_pol
00145   (const Healpix_Map<T> &mapT,
00146    const Healpix_Map<T> &mapQ,
00147    const Healpix_Map<T> &mapU,
00148    Alm<xcomplex<T> > &almT,
00149    Alm<xcomplex<T> > &almG,
00150    Alm<xcomplex<T> > &almC,
00151    const arr<double> &weight,
00152    bool add_alm)
00153   {
00154   planck_assert (mapT.Scheme()==RING,
00155     "map2alm_pol: maps must be in RING scheme");
00156   planck_assert (mapT.conformable(mapQ) && mapT.conformable(mapU),
00157     "map2alm_pol: maps are not conformable");
00158   planck_assert (weight.size()>=2*mapT.Nside(),
00159     "map2alm_pol: at least one weight array has too few entries");
00160 
00161   vector<ringpair> pair;
00162   healpix2ringpairs(mapT,weight,pair);
00163   map2alm_pol(pair,&mapT[0],&mapQ[0],&mapU[0],almT,almG,almC,add_alm);
00164   }
00165 
00166 template void map2alm_pol
00167   (const Healpix_Map<float> &mapT,
00168    const Healpix_Map<float> &mapQ,
00169    const Healpix_Map<float> &mapU,
00170    Alm<xcomplex<float> > &almT,
00171    Alm<xcomplex<float> > &almG,
00172    Alm<xcomplex<float> > &almC,
00173    const arr<double> &weight,
00174    bool add_alm);
00175 template void map2alm_pol
00176   (const Healpix_Map<double> &mapT,
00177    const Healpix_Map<double> &mapQ,
00178    const Healpix_Map<double> &mapU,
00179    Alm<xcomplex<double> > &almT,
00180    Alm<xcomplex<double> > &almG,
00181    Alm<xcomplex<double> > &almC,
00182    const arr<double> &weight,
00183    bool add_alm);
00184 
00185 template<typename T> void map2alm_pol_iter
00186   (const Healpix_Map<T> &mapT,
00187    const Healpix_Map<T> &mapQ,
00188    const Healpix_Map<T> &mapU,
00189    Alm<xcomplex<T> > &almT,
00190    Alm<xcomplex<T> > &almG,
00191    Alm<xcomplex<T> > &almC,
00192    int num_iter,
00193    const arr<double> &weight)
00194   {
00195   map2alm_pol(mapT,mapQ,mapU,almT,almG,almC,weight);
00196   for (int iter=1; iter<=num_iter; ++iter)
00197     {
00198     Healpix_Map<T> mapT2(mapT.Nside(),mapT.Scheme(),SET_NSIDE),
00199                    mapQ2(mapT.Nside(),mapT.Scheme(),SET_NSIDE),
00200                    mapU2(mapT.Nside(),mapT.Scheme(),SET_NSIDE);
00201 
00202     alm2map_pol(almT,almG,almC,mapT2,mapQ2,mapU2);
00203     for (int m=0; m<mapT.Npix(); ++m)
00204       {
00205       mapT2[m] = mapT[m]-mapT2[m];
00206       mapQ2[m] = mapQ[m]-mapQ2[m];
00207       mapU2[m] = mapU[m]-mapU2[m];
00208       }
00209     map2alm_pol(mapT2,mapQ2,mapU2,almT,almG,almC,weight,true);
00210     }
00211   }
00212 
00213 template void map2alm_pol_iter
00214   (const Healpix_Map<float> &mapT,
00215    const Healpix_Map<float> &mapQ,
00216    const Healpix_Map<float> &mapU,
00217    Alm<xcomplex<float> > &almT,
00218    Alm<xcomplex<float> > &almG,
00219    Alm<xcomplex<float> > &almC,
00220    int num_iter,
00221    const arr<double> &weight);
00222 template void map2alm_pol_iter
00223   (const Healpix_Map<double> &mapT,
00224    const Healpix_Map<double> &mapQ,
00225    const Healpix_Map<double> &mapU,
00226    Alm<xcomplex<double> > &almT,
00227    Alm<xcomplex<double> > &almG,
00228    Alm<xcomplex<double> > &almC,
00229    int num_iter,
00230    const arr<double> &weight);
00231 
00232 template<typename T> void map2alm_pol_iter2
00233   (const Healpix_Map<T> &mapT,
00234    const Healpix_Map<T> &mapQ,
00235    const Healpix_Map<T> &mapU,
00236    Alm<xcomplex<T> > &almT,
00237    Alm<xcomplex<T> > &almG,
00238    Alm<xcomplex<T> > &almC,
00239    double err_abs, double err_rel)
00240   {
00241   arr<double> wgt(2*mapT.Nside());
00242   wgt.fill(1);
00243   Healpix_Map<T> mapT2(mapT), mapQ2(mapQ), mapU2(mapU);
00244   almT.SetToZero(); almG.SetToZero(); almC.SetToZero();
00245   while(true)
00246     {
00247     map2alm_pol(mapT2,mapQ2,mapU2,almT,almG,almC,wgt,true);
00248     alm2map_pol(almT,almG,almC,mapT2,mapQ2,mapU2);
00249     double errmeasure=0;
00250     for (int m=0; m<mapT.Npix(); ++m)
00251       {
00252       double err = abs(mapT[m]-mapT2[m]);
00253       double rel = (mapT[m]!=0) ? abs(err/mapT[m]) : 1e300;
00254       errmeasure = max(errmeasure,min(err/err_abs,rel/err_rel));
00255       mapT2[m] = mapT[m]-mapT2[m];
00256       err = abs(mapQ[m]-mapQ2[m]);
00257       rel = (mapQ[m]!=0) ? abs(err/mapQ[m]) : 1e300;
00258       errmeasure = max(errmeasure,min(err/err_abs,rel/err_rel));
00259       mapQ2[m] = mapQ[m]-mapQ2[m];
00260       err = abs(mapU[m]-mapU2[m]);
00261       rel = (mapU[m]!=0) ? abs(err/mapU[m]) : 1e300;
00262       errmeasure = max(errmeasure,min(err/err_abs,rel/err_rel));
00263       mapU2[m] = mapU[m]-mapU2[m];
00264       }
00265 cout << "map error measure: " << errmeasure << endl;
00266     if (errmeasure<1) break;
00267     }
00268   }
00269 
00270 template void map2alm_pol_iter2
00271   (const Healpix_Map<double> &mapT,
00272    const Healpix_Map<double> &mapQ,
00273    const Healpix_Map<double> &mapU,
00274    Alm<xcomplex<double> > &almT,
00275    Alm<xcomplex<double> > &almG,
00276    Alm<xcomplex<double> > &almC,
00277    double err_abs, double err_rel);
00278 
00279 
00280 template<typename T> void alm2map (const Alm<xcomplex<T> > &alm,
00281   Healpix_Map<T> &map)
00282   {
00283   planck_assert (map.Scheme()==RING, "alm2map: map must be in RING scheme");
00284 
00285   vector<ringpair> pair;
00286   healpix2ringpairs(map,pair);
00287   alm2map(alm,pair,&map[0]);
00288   }
00289 
00290 template void alm2map (const Alm<xcomplex<double> > &alm,
00291   Healpix_Map<double> &map);
00292 template void alm2map (const Alm<xcomplex<float> > &alm,
00293   Healpix_Map<float> &map);
00294 
00295 
00296 template<typename T> void alm2map_pol
00297   (const Alm<xcomplex<T> > &almT,
00298    const Alm<xcomplex<T> > &almG,
00299    const Alm<xcomplex<T> > &almC,
00300    Healpix_Map<T> &mapT,
00301    Healpix_Map<T> &mapQ,
00302    Healpix_Map<T> &mapU)
00303   {
00304   planck_assert (mapT.Scheme()==RING,
00305     "alm2map_pol: maps must be in RING scheme");
00306   planck_assert (mapT.conformable(mapQ) && mapT.conformable(mapU),
00307     "alm2map_pol: maps are not conformable");
00308 
00309   vector<ringpair> pair;
00310   healpix2ringpairs(mapT,pair);
00311   alm2map_pol(almT,almG,almC,pair,&mapT[0],&mapQ[0],&mapU[0]);
00312   }
00313 
00314 template void alm2map_pol (const Alm<xcomplex<double> > &almT,
00315                            const Alm<xcomplex<double> > &almG,
00316                            const Alm<xcomplex<double> > &almC,
00317                            Healpix_Map<double> &mapT,
00318                            Healpix_Map<double> &mapQ,
00319                            Healpix_Map<double> &mapU);
00320 
00321 template void alm2map_pol (const Alm<xcomplex<float> > &almT,
00322                            const Alm<xcomplex<float> > &almG,
00323                            const Alm<xcomplex<float> > &almC,
00324                            Healpix_Map<float> &mapT,
00325                            Healpix_Map<float> &mapQ,
00326                            Healpix_Map<float> &mapU);
00327 
00328 
00329 template<typename T> void alm2map_der1
00330   (const Alm<xcomplex<T> > &alm,
00331    Healpix_Map<T> &map,
00332    Healpix_Map<T> &mapdth,
00333    Healpix_Map<T> &mapdph)
00334   {
00335   planck_assert (map.Scheme()==RING,
00336     "alm2map_der1: maps must be in RING scheme");
00337   planck_assert (map.conformable(mapdth) && map.conformable(mapdph),
00338     "alm2map_der1: maps are not conformable");
00339 
00340   vector<ringpair> pair;
00341   healpix2ringpairs(map,pair);
00342   alm2map_der1(alm,pair,&map[0],&mapdth[0],&mapdph[0]);
00343   }
00344 
00345 template void alm2map_der1 (const Alm<xcomplex<double> > &alm,
00346                             Healpix_Map<double> &map,
00347                             Healpix_Map<double> &map_dth,
00348                             Healpix_Map<double> &map_dph);
00349 
00350 template void alm2map_der1 (const Alm<xcomplex<float> > &alm,
00351                             Healpix_Map<float> &map,
00352                             Healpix_Map<float> &map_dth,
00353                             Healpix_Map<float> &map_dph);