hpxtest.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) 2004, 2005, 2006 Max-Planck-Society
00029  *  Author: Martin Reinecke
00030  */
00031 
00032 /*
00033 
00034 Candidates for testing the validity of the Healpix routines:
00035 
00036 - done: ang2pix(pix2ang(i)) = i for all Healpix_Bases
00037 - done: pix2ang(ang2pix(ptg)) dot ptg > 1-delta for all Healpix_Bases
00038 - done: ring2nest(nest2ring(i)) = i for all hierarchical Healpix_Bases
00039 - done: downgrade(upgrade(map)) = map for all maps
00040 - done: map and downgraded map should have same average
00041 - done: alm2map(map2alm(map)) approx map (same for pol)
00042 - partly done: neighbor tests
00043 - powspec -> alm -> powspec (should produce similar powspecs, also for pol)
00044 - done: two swap_schemes() should produce original map
00045 - done: query_disc tests (dot products etc.)
00046 - a_lms: test Set(), Scale(), Add(), alm(l,m) = alm.mstart(m)[l], etc.
00047 
00048 */
00049 
00050 #include <iostream>
00051 #include "healpix_base.h"
00052 #include "healpix_base2.h"
00053 #include "healpix_map.h"
00054 #include "arr.h"
00055 #include "planck_rng.h"
00056 #include "lsconstants.h"
00057 #include "alm.h"
00058 #include "alm_healpix_tools.h"
00059 #include "alm_powspec_tools.h"
00060 #include "geom_utils.h"
00061 
00062 using namespace std;
00063 
00064 const int nsamples = 1000000;
00065 
00066 planck_rng rng;
00067 
00068 void random_dir (pointing &ptg)
00069   {
00070   ptg.theta = acos(rng.rand_uni()*2-1);
00071   ptg.phi = rng.rand_uni()*twopi;
00072   }
00073 
00074 void check_ringnestring()
00075   {
00076   cout << "testing ring2nest(nest2ring(m))==m" << endl;
00077   for (int order=0; order<=13; ++order)
00078     {
00079     cout << "order = " << order << endl;
00080     Healpix_Base base (order,RING);
00081     for (int m=0; m<nsamples; ++m)
00082       {
00083       int pix = int(rng.rand_uni()*base.Npix());
00084       if (base.ring2nest(base.nest2ring(pix))!=pix)
00085         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00086       }
00087     }
00088   }
00089 void check_ringnestring2()
00090   {
00091   cout << "testing ring2nest(nest2ring(m))==m" << endl;
00092   for (int order=0; order<=29; ++order)
00093     {
00094     cout << "order = " << order << endl;
00095     Healpix_Base2 base (order,RING);
00096     for (int m=0; m<nsamples; ++m)
00097       {
00098       int64 pix = int64(rng.rand_uni()*base.Npix());
00099       if (base.ring2nest(base.nest2ring(pix))!=pix)
00100         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00101       }
00102     }
00103   }
00104 
00105 void check_nestpeanonest()
00106   {
00107   cout << "testing peano2nest(nest2peano(m))==m" << endl;
00108   for (int order=0; order<=13; ++order)
00109     {
00110     cout << "order = " << order << endl;
00111     Healpix_Base base (order,NEST);
00112     for (int m=0; m<nsamples; ++m)
00113       {
00114       int pix = int(rng.rand_uni()*base.Npix());
00115       if (base.peano2nest(base.nest2peano(pix))!=pix)
00116         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00117       }
00118     }
00119   }
00120 void check_nestpeanonest2()
00121   {
00122   cout << "testing peano2nest(nest2peano(m))==m" << endl;
00123   for (int order=0; order<=29; ++order)
00124     {
00125     cout << "order = " << order << endl;
00126     Healpix_Base2 base (order,NEST);
00127     for (int m=0; m<nsamples; ++m)
00128       {
00129       int64 pix = int64(rng.rand_uni()*base.Npix());
00130       if (base.peano2nest(base.nest2peano(pix))!=pix)
00131         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00132       }
00133     }
00134   }
00135 
00136 void check_pixangpix()
00137   {
00138   cout << "testing ang2pix(pix2ang(m))==m" << endl;
00139   for (int order=0; order<=13; ++order)
00140     {
00141     cout << "order = " << order << endl;
00142     Healpix_Base base1 (order,RING);
00143     Healpix_Base base2 (order,NEST);
00144     for (int m=0; m<nsamples; ++m)
00145       {
00146       int pix = int(rng.rand_uni()*base1.Npix());
00147       if (base1.ang2pix(base1.pix2ang(pix))!=pix)
00148         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00149       if (base2.ang2pix(base2.pix2ang(pix))!=pix)
00150         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00151       }
00152     }
00153   for (int nside=3; nside<(1<<13); nside+=nside/2+1)
00154     {
00155     cout << "nside = " << nside << endl;
00156     Healpix_Base base (nside,RING,SET_NSIDE);
00157     for (int m=0; m<nsamples; ++m)
00158       {
00159       int pix = int(rng.rand_uni()*base.Npix());
00160       if (base.ang2pix(base.pix2ang(pix))!=pix)
00161         cout << "  PROBLEM: nside = " << nside << ", pixel = " << pix << endl;
00162       }
00163     }
00164   }
00165 void check_pixangpix2()
00166   {
00167   cout << "testing ang2pix(pix2ang(m))==m" << endl;
00168   for (int order=0; order<=29; ++order)
00169     {
00170     cout << "order = " << order << endl;
00171     Healpix_Base2 base1 (order,RING);
00172     Healpix_Base2 base2 (order,NEST);
00173     for (int m=0; m<nsamples; ++m)
00174       {
00175       int64 pix = int64(rng.rand_uni()*base1.Npix());
00176       if (base1.ang2pix(base1.pix2ang(pix))!=pix)
00177         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00178       if (base2.ang2pix(base2.pix2ang(pix))!=pix)
00179         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00180       }
00181     }
00182   for (int nside=3; nside<(1<<29); nside+=nside/2+1)
00183     {
00184     cout << "nside = " << nside << endl;
00185     Healpix_Base2 base (nside,RING,SET_NSIDE);
00186     for (int m=0; m<nsamples; ++m)
00187       {
00188       int64 pix = int64(rng.rand_uni()*base.Npix());
00189       if (base.ang2pix(base.pix2ang(pix))!=pix)
00190         cout << "  PROBLEM: nside = " << nside << ", pixel = " << pix << endl;
00191       }
00192     }
00193   }
00194 
00195 void check_pixvecpix()
00196   {
00197   cout << "testing vec2pix(pix2vec(m))==m" << endl;
00198   for (int order=0; order<=13; ++order)
00199     {
00200     cout << "order = " << order << endl;
00201     Healpix_Base base1 (order,RING);
00202     Healpix_Base base2 (order,NEST);
00203     for (int m=0; m<nsamples; ++m)
00204       {
00205       int pix = int(rng.rand_uni()*base1.Npix());
00206       if (base1.vec2pix(base1.pix2vec(pix))!=pix)
00207         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00208       if (base2.vec2pix(base2.pix2vec(pix))!=pix)
00209         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00210       }
00211     }
00212   for (int nside=3; nside<(1<<13); nside+=nside/2+1)
00213     {
00214     cout << "nside = " << nside << endl;
00215     Healpix_Base base (nside,RING,SET_NSIDE);
00216     for (int m=0; m<nsamples; ++m)
00217       {
00218       int pix = int(rng.rand_uni()*base.Npix());
00219       if (base.vec2pix(base.pix2vec(pix))!=pix)
00220         cout << "  PROBLEM: nside = " << nside << ", pixel = " << pix << endl;
00221       }
00222     }
00223   }
00224 void check_pixvecpix2()
00225   {
00226   cout << "testing vec2pix(pix2vec(m))==m" << endl;
00227   for (int order=0; order<=29; ++order)
00228     {
00229     cout << "order = " << order << endl;
00230     Healpix_Base2 base1 (order,RING);
00231     Healpix_Base2 base2 (order,NEST);
00232     for (int m=0; m<nsamples; ++m)
00233       {
00234       int64 pix = int64(rng.rand_uni()*base1.Npix());
00235       if (base1.vec2pix(base1.pix2vec(pix))!=pix)
00236         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00237       if (base2.vec2pix(base2.pix2vec(pix))!=pix)
00238         cout << "  PROBLEM: order = " << order << ", pixel = " << pix << endl;
00239       }
00240     }
00241   for (int nside=3; nside<(1<<29); nside+=nside/2+1)
00242     {
00243     cout << "nside = " << nside << endl;
00244     Healpix_Base2 base (nside,RING,SET_NSIDE);
00245     for (int m=0; m<nsamples; ++m)
00246       {
00247       int64 pix = int64(rng.rand_uni()*base.Npix());
00248       if (base.vec2pix(base.pix2vec(pix))!=pix)
00249         cout << "  PROBLEM: nside = " << nside << ", pixel = " << pix << endl;
00250       }
00251     }
00252   }
00253 
00254 void check_angpixang()
00255   {
00256   cout << "testing pix2ang(ang2pix(ptg)) approx ptg" << endl;
00257   for (int order=0; order<=13; ++order)
00258     {
00259     cout << "order = " << order << endl;
00260     Healpix_Base base1 (order,NEST);
00261     Healpix_Base base2 (order,RING);
00262     double maxang = base1.max_pixrad();
00263     for (int m=0; m<nsamples; ++m)
00264       {
00265       pointing ptg;
00266       random_dir (ptg);
00267       if (v_angle(base1.pix2ang(base1.ang2pix(ptg)),ptg)>maxang)
00268         cout << "  PROBLEM: order = " << order << ", ptg = " << ptg << endl;
00269       if (v_angle(base2.pix2ang(base2.ang2pix(ptg)),ptg)>maxang)
00270         cout << "  PROBLEM: order = " << order << ", ptg = " << ptg << endl;
00271       }
00272     }
00273   for (int nside=3; nside<(1<<13); nside+=nside/2+1)
00274     {
00275     cout << "nside = " << nside << endl;
00276     Healpix_Base base (nside,RING,SET_NSIDE);
00277     double maxang = base.max_pixrad();
00278     for (int m=0; m<nsamples; ++m)
00279       {
00280       pointing ptg;
00281       random_dir (ptg);
00282       if (v_angle(base.pix2ang(base.ang2pix(ptg)),ptg)>maxang)
00283         cout << "  PROBLEM: nside = " << nside << ", ptg = " << ptg << endl;
00284       }
00285     }
00286   }
00287 void check_angpixang2()
00288   {
00289   cout << "testing pix2ang(ang2pix(ptg)) approx ptg" << endl;
00290   for (int order=0; order<=29; ++order)
00291     {
00292     cout << "order = " << order << endl;
00293     Healpix_Base2 base1 (order,NEST);
00294     Healpix_Base2 base2 (order,RING);
00295     double maxang = base1.max_pixrad();
00296     for (int m=0; m<nsamples; ++m)
00297       {
00298       pointing ptg;
00299       random_dir (ptg);
00300       if (v_angle(base1.pix2ang(base1.ang2pix(ptg)),ptg)>maxang)
00301         cout << "  PROBLEM: order = " << order << ", ptg = " << ptg << endl;
00302       if (v_angle(base2.pix2ang(base2.ang2pix(ptg)),ptg)>maxang)
00303         cout << "  PROBLEM: order = " << order << ", ptg = " << ptg << endl;
00304       }
00305     }
00306   for (int nside=3; nside<(1<<29); nside+=nside/2+1)
00307     {
00308     cout << "nside = " << nside << endl;
00309     Healpix_Base2 base (nside,RING,SET_NSIDE);
00310     double maxang = base.max_pixrad();
00311     for (int m=0; m<nsamples; ++m)
00312       {
00313       pointing ptg;
00314       random_dir (ptg);
00315       if (v_angle(base.pix2ang(base.ang2pix(ptg)),ptg)>maxang)
00316         cout << "  PROBLEM: nside = " << nside << ", ptg = " << ptg << endl;
00317       }
00318     }
00319   }
00320 
00321 void check_neighbors()
00322   {
00323   cout << "testing neighbor function" << endl;
00324   for (int order=0; order<=13; ++order)
00325     {
00326     cout << "order = " << order << endl;
00327     Healpix_Base base (order,NEST), base2(order,RING);
00328     double maxang = 2.01*base.max_pixrad();
00329     for (int m=0; m<nsamples/10; ++m)
00330       {
00331       int pix = int(rng.rand_uni()*base.Npix());
00332       fix_arr<int,8> nb;
00333       vec3 pixpt = base.pix2ang(pix);
00334       base.neighbors(pix,nb);
00335       sort(&nb[0],&nb[0]+8);
00336       int check=0;
00337       for (int n=0; n<8; ++n)
00338         {
00339         if (nb[n]<0)
00340           {
00341           ++check;
00342           }
00343         else
00344           {
00345           if (v_angle(base.pix2ang(nb[n]),pixpt)>maxang)
00346             cout << " PROBLEM: order = " << order << ", pix = " << pix << endl;
00347           if ((n>0) && (nb[n]==nb[n-1]))
00348             cout << " PROBLEM: order = " << order << ", pix = " << pix << endl;
00349           }
00350         }
00351       planck_assert((check<=1)||((order==0)&&(check<=2)),"too few neighbors");
00352       pixpt = base2.pix2ang(pix);
00353       base2.neighbors(pix,nb);
00354       for (int n=0; n<8; ++n)
00355         if ((nb[n]>=0) && (v_angle(base2.pix2ang(nb[n]),pixpt)>maxang))
00356           cout << "  PROBLEM2: order = " << order << ", pix = " << pix << endl;
00357       }
00358     }
00359   for (int nside=3; nside<(1<<13); nside+=nside/2+1)
00360     {
00361     cout << "nside = " << nside << endl;
00362     Healpix_Base base (nside,RING,SET_NSIDE);
00363     double maxang = 2.01*base.max_pixrad();
00364     for (int m=0; m<nsamples/10; ++m)
00365       {
00366       int pix = int(rng.rand_uni()*base.Npix());
00367       fix_arr<int,8> nb;
00368       vec3 pixpt = base.pix2ang(pix);
00369       base.neighbors(pix,nb);
00370       for (int n=0; n<8; ++n)
00371         if ((nb[n]>=0) && (v_angle(base.pix2ang(nb[n]),pixpt)>maxang))
00372           cout << "  PROBLEM: nside = " << nside << ", pix = " << pix << endl;
00373       }
00374     }
00375   }
00376 void check_neighbors2()
00377   {
00378   cout << "testing neighbor function" << endl;
00379   for (int order=0; order<=29; ++order)
00380     {
00381     cout << "order = " << order << endl;
00382     Healpix_Base2 base (order,NEST), base2(order,RING);
00383     double maxang = 2.01*base.max_pixrad();
00384     for (int m=0; m<nsamples/10; ++m)
00385       {
00386       int64 pix = int64(rng.rand_uni()*base.Npix());
00387       fix_arr<int64,8> nb;
00388       vec3 pixpt = base.pix2ang(pix);
00389       base.neighbors(pix,nb);
00390       sort(&nb[0],&nb[0]+8);
00391       int check=0;
00392       for (int n=0; n<8; ++n)
00393         {
00394         if (nb[n]<0)
00395           {
00396           ++check;
00397           }
00398         else
00399           {
00400           if (v_angle(base.pix2ang(nb[n]),pixpt)>maxang)
00401             cout << " PROBLEM: order = " << order << ", pix = " << pix << endl;
00402           if ((n>0) && (nb[n]==nb[n-1]))
00403             cout << " PROBLEM: order = " << order << ", pix = " << pix << endl;
00404           }
00405         }
00406       planck_assert((check<=1)||((order==0)&&(check<=2)),"too few neighbors");
00407       pixpt = base2.pix2ang(pix);
00408       base2.neighbors(pix,nb);
00409       for (int n=0; n<8; ++n)
00410         if ((nb[n]>=0) && (v_angle(base2.pix2ang(nb[n]),pixpt)>maxang))
00411           cout << "  PROBLEM2: order = " << order << ", pix = " << pix << endl;
00412       }
00413     }
00414   for (int nside=3; nside<(1<<29); nside+=nside/2+1)
00415     {
00416     cout << "nside = " << nside << endl;
00417     Healpix_Base2 base (nside,RING,SET_NSIDE);
00418     double maxang = 2.01*base.max_pixrad();
00419     for (int m=0; m<nsamples/10; ++m)
00420       {
00421       int64 pix = int64(rng.rand_uni()*base.Npix());
00422       fix_arr<int64,8> nb;
00423       vec3 pixpt = base.pix2ang(pix);
00424       base.neighbors(pix,nb);
00425       for (int n=0; n<8; ++n)
00426         if ((nb[n]>=0) && (v_angle(base.pix2ang(nb[n]),pixpt)>maxang))
00427           cout << "  PROBLEM: nside = " << nside << ", pix = " << pix << endl;
00428       }
00429     }
00430   }
00431 
00432 void check_swap_scheme()
00433   {
00434   cout << "testing whether double swap_scheme() returns the original map"
00435        << endl << "(for orders 0 to 10)." << endl;
00436   for (int order=0; order<=10; ++order)
00437     {
00438     cout << "order = " << order << endl;
00439     Healpix_Map<unsigned char> map(order,NEST);
00440     for (int m=0; m<map.Npix(); ++m) map[m]=m&0xFF;
00441     map.swap_scheme();
00442     map.swap_scheme();
00443     for (int m=0; m<map.Npix(); ++m)
00444       if (map[m]!=(m&0xFF))
00445         cout << "  PROBLEM: order = " << order << ", pix = " << m << endl;
00446     }
00447   }
00448 
00449 void check_query_disc()
00450   {
00451   cout << "testing whether all pixels found by query_disc() really" << endl
00452        << "lie inside the disk (and vice versa)" << endl;
00453   for (int order=0; order<=5; ++order)
00454     {
00455     cout << "order = " << order << endl;
00456     Healpix_Map <bool> map (order,RING);
00457     map.fill(false);
00458     vector<int> list;
00459     for (int m=0; m<100000; ++m)
00460       {
00461       pointing ptg;
00462       random_dir (ptg);
00463       double rad = pi/1 * rng.rand_uni();
00464       map.query_disc(ptg,rad,list);
00465       vec3 vptg=ptg;
00466       double cosrad=cos(rad);
00467       for (unsigned int i=0; i<list.size(); ++i)
00468         map[list[i]] = true;
00469       for (int i=0; i<map.Npix(); ++i)
00470         {
00471         bool inside = dotprod(map.pix2ang(i),vptg)>cosrad;
00472         if (inside^map[i])
00473           cout << "  PROBLEM: order = " << order << ", ptg = " << ptg << endl;
00474         }
00475       for (unsigned int i=0; i<list.size(); ++i)
00476         map[list[i]] = false;
00477       }
00478     }
00479   }
00480 
00481 void helper_oop (int order)
00482   {
00483   Healpix_Map<double> map (order,RING), map2 (order,NEST), map3 (order,RING);
00484   for (int m=0; m<map.Npix(); ++m) map[m] = rng.rand_uni()+0.01;
00485   map2.Import(map);
00486   map3.Import(map2);
00487   for (int m=0; m<map.Npix(); ++m)
00488     if (!approx(map[m],map3[m],1e-12))
00489       cout << "PROBLEM: order = " << order << endl;
00490   }
00491 void helper_udgrade (int order, Healpix_Ordering_Scheme s1,
00492   Healpix_Ordering_Scheme s2)
00493   {
00494   Healpix_Map<double> map (order,s1), map2 (order+2,s2), map3 (order, s1);
00495   for (int m=0; m<map.Npix(); ++m) map[m] = rng.rand_uni()+0.01;
00496   map2.Import(map);
00497   map3.Import(map2);
00498   for (int m=0; m<map.Npix(); ++m)
00499     if (!approx(map[m],map3[m],1e-12))
00500       cout << "PROBLEM: order = " << order << endl;
00501   }
00502 void helper_udgrade2 (int nside)
00503   {
00504   Healpix_Map<double> map (nside,RING,SET_NSIDE), map2 (nside*3,RING,SET_NSIDE),
00505     map3 (nside, RING,SET_NSIDE);
00506   for (int m=0; m<map.Npix(); ++m) map[m] = rng.rand_uni()+0.01;
00507   map2.Import(map);
00508   map3.Import(map2);
00509   for (int m=0; m<map.Npix(); ++m)
00510     if (!approx(map[m],map3[m],1e-12))
00511       cout << "PROBLEM: nside = " << nside << endl;
00512   }
00513 
00514 void check_import()
00515   {
00516   cout << "testing out-of-place swapping" << endl;
00517   for (int order=0; order<=7; ++order)
00518     {
00519     cout << "order = " << order << endl;
00520     helper_oop(order);
00521     }
00522   cout << "testing downgrade(upgrade(map)) == map" << endl;
00523   for (int order=0; order<=7; ++order)
00524     {
00525     cout << "order = " << order << endl;
00526     cout << "RING, RING" << endl;
00527     helper_udgrade(order,RING,RING);
00528     cout << "RING, NEST" << endl;
00529     helper_udgrade(order,RING,NEST);
00530     cout << "NEST, NEST" << endl;
00531     helper_udgrade(order,NEST,NEST);
00532     cout << "NEST, RING" << endl;
00533     helper_udgrade(order,NEST,RING);
00534     }
00535   for (int nside=3; nside<500; nside+=nside/2+1)
00536     {
00537     cout << "nside = " << nside << endl;
00538     helper_udgrade2(nside);
00539     }
00540   }
00541 
00542 void check_average()
00543   {
00544   cout << "testing whether average(map) == average(downgraded map)" << endl;
00545   for (int order=1; order<=10; ++order)
00546     {
00547     cout << "order = " << order << endl;
00548     Healpix_Map<double> map (order,RING), map2(1,RING);
00549     for (int m=0; m<map.Npix(); ++m)
00550       map[m] = rng.rand_uni()+0.01;
00551     map2.Import(map);
00552     double avg=map.average(), avg2=map2.average();
00553     if (!approx(avg,avg2,1e-10))
00554       cout << "PROBLEM: order = " << order << " " << avg/avg2-1 << endl;
00555     }
00556   for (int nside=3; nside<1000; nside += nside/2+1)
00557     {
00558     cout << "nside = " << nside << endl;
00559     Healpix_Map<double> map (nside,RING,SET_NSIDE), map2(1,RING,SET_NSIDE);
00560     for (int m=0; m<map.Npix(); ++m)
00561       map[m] = rng.rand_uni()+0.01;
00562     map2.Import(map);
00563     double avg=map.average(), avg2=map2.average();
00564     if (!approx(avg,avg2,1e-10))
00565       cout << "PROBLEM: nside = " << nside << " " << avg/avg2-1 << endl;
00566     }
00567   }
00568 
00569 void random_alm (Alm<xcomplex<double> >&almT, Alm<xcomplex<double> >&almG,
00570   Alm<xcomplex<double> >&almC, int lmax, int mmax)
00571   {
00572   almT.Set(lmax,mmax); almG.Set(lmax,mmax); almC.Set(lmax,mmax);
00573 
00574   for (int l=0; l<=lmax; ++l)
00575     {
00576     almT(l,0).re=rng.rand_gauss(); almT(l,0).im=0.;
00577     almG(l,0).re=rng.rand_gauss(); almG(l,0).im=0.;
00578     almC(l,0).re=rng.rand_gauss(); almC(l,0).im=0.;
00579     }
00580 
00581   for (int m=1; m<=mmax; ++m)
00582     for (int l=m; l<=lmax; ++l)
00583       {
00584       almT(l,m).re=rng.rand_gauss(); almT(l,m).im=rng.rand_gauss();
00585       almG(l,m).re=rng.rand_gauss(); almG(l,m).im=rng.rand_gauss();
00586       almC(l,m).re=rng.rand_gauss(); almC(l,m).im=rng.rand_gauss();
00587       }
00588   almG(0,0)=almC(0,0)=almG(1,0)=almC(1,0)=almG(1,1)=almC(1,1)=0;
00589   }
00590 
00591 void random_alm (Alm<xcomplex<double> >&alm, int lmax, int mmax)
00592   {
00593   alm.Set(lmax,mmax);
00594 
00595   for (int l=0; l<=lmax; ++l)
00596     { alm(l,0).re=rng.rand_gauss(); alm(l,0).im=0.; }
00597 
00598   for (int m=1; m<=mmax; ++m)
00599     for (int l=m; l<=lmax; ++l)
00600       { alm(l,m).re=rng.rand_gauss(); alm(l,m).im=rng.rand_gauss(); }
00601   }
00602 
00603 void check_alm (const Alm<xcomplex<double> >&oalm,
00604   const Alm<xcomplex<double> >&alm, double epsilon)
00605   {
00606   for (int m=0; m<=alm.Mmax(); ++m)
00607     for (int l=m; l<=alm.Lmax(); ++l)
00608       {
00609       if (!abs_approx(oalm(l,m).re,alm(l,m).re,epsilon))
00610         cout << "Problemr " << l << " " << m << endl;
00611       if (!abs_approx(oalm(l,m).im,alm(l,m).im,epsilon))
00612         cout << "Problemi " << l << " " << m <<  endl;
00613       }
00614   }
00615 
00616 void check_alm2map2alm (int lmax, int mmax, int nside)
00617   {
00618   cout << "testing whether a_lm->map->a_lm returns original a_lm" << endl;
00619   cout << "lmax=" << lmax <<", mmax=" << mmax << ", nside=" << nside << endl;
00620   const double epsilon = 1e-8;
00621   Alm<xcomplex<double> > oalmT(lmax,mmax),almT(lmax,mmax),
00622     oalmG(lmax,mmax),almG(lmax,mmax),oalmC(lmax,mmax),almC(lmax,mmax);
00623   Healpix_Map<double> mapT(nside,RING,SET_NSIDE), mapQ(nside,RING,SET_NSIDE),
00624     mapU(nside,RING,SET_NSIDE);
00625 
00626   random_alm(oalmT,oalmG,oalmC,lmax,mmax);
00627   alm2map(oalmT,mapT);
00628   map2alm_iter2(mapT,almT,1e-12,1e-12);
00629   check_alm (oalmT, almT, epsilon);
00630 
00631   alm2map_pol(oalmT,oalmG,oalmC,mapT,mapQ,mapU);
00632   map2alm_pol_iter2(mapT,mapQ,mapU,almT,almG,almC,1e-12,1e-12);
00633   check_alm (oalmT, almT, epsilon);
00634   check_alm (oalmG, almG, epsilon);
00635   check_alm (oalmC, almC, epsilon);
00636   }
00637 
00638 void check_smooth_alm ()
00639   {
00640   cout << "testing whether unsmooth(smooth(a_lm)) returns a_lm" << endl;
00641   const double epsilon = 1e-14;
00642   const double fwhm_arcmin = 100;
00643   const int lmax=300, mmax=300;
00644   Alm<xcomplex<double> > oalmT(lmax,mmax),almT(lmax,mmax),
00645     oalmG(lmax,mmax),almG(lmax,mmax),oalmC(lmax,mmax),almC(lmax,mmax);
00646 
00647   random_alm(oalmT,oalmG,oalmC,lmax,mmax);
00648 
00649   almT=oalmT; almG=oalmG; almC=oalmC;
00650   smooth_with_Gauss (almT, fwhm_arcmin);
00651   smooth_with_Gauss (almT, -fwhm_arcmin);
00652   check_alm (oalmT, almT, epsilon);
00653   almT=oalmT;
00654   smooth_with_Gauss (almT, almG, almC, fwhm_arcmin);
00655   smooth_with_Gauss (almT, almG, almC, -fwhm_arcmin);
00656   check_alm (oalmT, almT, epsilon);
00657   check_alm (oalmG, almG, epsilon);
00658   check_alm (oalmC, almC, epsilon);
00659   }
00660 
00661 void check_rot_alm ()
00662   {
00663   cout << "testing whether rot^-1(rot(a_lm)) returns a_lm" << endl;
00664   const double epsilon = 2e-13;
00665   const int lmax=300;
00666   Alm<xcomplex<double> > oalm(lmax,lmax),alm(lmax,lmax);
00667 
00668   random_alm(oalm,lmax,lmax);
00669 
00670   alm=oalm;
00671   rotate_alm (alm,3,4,5);
00672   rotate_alm (alm,-5,-4,-3);
00673   check_alm (oalm, alm, epsilon);
00674   }
00675 
00676 int main()
00677   {
00678   check_smooth_alm();
00679   check_rot_alm();
00680   check_alm2map2alm(620,620,256);
00681   check_alm2map2alm(620,2,256);
00682   check_average();
00683   check_import();
00684   check_neighbors();
00685   check_neighbors2();
00686   check_pixangpix();
00687   check_pixangpix2();
00688   check_pixvecpix();
00689   check_pixvecpix2();
00690   check_angpixang();
00691   check_angpixang2();
00692   check_ringnestring();
00693   check_ringnestring2();
00694   check_nestpeanonest();
00695   check_nestpeanonest2();
00696   check_query_disc();
00697   check_swap_scheme();
00698   }