36 #include "sharp_cxx.h" 45 cout <<
"\nWARNING: map analysis requested with lmax>4*nside...\n" 46 "is this really what you want?\n\n";
56 "map2alm: weight array has too few entries");
58 checkLmaxNside(alm.Lmax(), map.
Nside());
61 job.set_weighted_Healpix_geometry (map.
Nside(),&weight[0]);
62 job.set_triangular_alm_info (alm.Lmax(), alm.Mmax());
63 job.map2alm(&map[0], &alm(0,0), add_alm);
74 Alm<xcomplex<T> > &alm,
bool add_alm)
77 "alm2map_adjoint: map must be in RING scheme");
79 checkLmaxNside(alm.Lmax(), map.
Nside());
82 job.set_Healpix_geometry (map.
Nside());
83 job.set_triangular_alm_info (alm.Lmax(), alm.Mmax());
84 job.alm2map_adjoint(&map[0], &alm(0,0), add_alm);
88 Alm<xcomplex<float> > &alm,
bool add_alm);
90 Alm<xcomplex<double> > &alm,
bool add_alm);
96 for (
int iter=1; iter<=num_iter; ++iter)
100 for (
int m=0; m<map.
Npix(); ++m)
101 map2[m] = map[m]-map2[m];
107 Alm<xcomplex<float> > &alm,
int num_iter,
110 Alm<xcomplex<double> > &alm,
int num_iter,
113 template<
typename T>
void map2alm_iter2 (
const Healpix_Map<T> &map,
114 Alm<xcomplex<T> > &alm,
double err_abs,
double err_rel)
116 double x_err_abs=1./err_abs, x_err_rel=1./err_rel;
126 for (
int m=0; m<map.
Npix(); ++m)
128 double err = abs(map[m]-map2[m]);
129 double rel = (map[m]!=0) ? abs(err/map[m]) : 1e300;
130 errmeasure = max(errmeasure,min(err*x_err_abs,rel*x_err_rel));
131 map2[m] = map[m]-map2[m];
134 if (errmeasure<1)
break;
139 Alm<xcomplex<double> > &alm,
double err_abs,
double err_rel);
142 template<
typename T>
void map2alm_spin
144 Alm<xcomplex<T> > &alm1,
Alm<xcomplex<T> > &alm2,
148 "map2alm_spin: maps must be in RING scheme");
150 "map2alm_spin: maps are not conformable");
152 "map2alm_spin: a_lm are not conformable");
154 "map2alm_spin: weight array has too few entries");
156 "map contains undefined pixels");
157 checkLmaxNside(alm1.Lmax(), map1.
Nside());
160 job.set_weighted_Healpix_geometry (map1.
Nside(),&weight[0]);
161 job.set_triangular_alm_info (alm1.Lmax(), alm1.Mmax());
162 job.map2alm_spin(&map1[0],&map2[0],&alm1(0,0),&alm2(0,0),spin,add_alm);
165 template void map2alm_spin
167 Alm<xcomplex<float> > &alm1,
Alm<xcomplex<float> > &alm2,
168 int spin,
const arr<double> &weight,
bool add_alm);
169 template void map2alm_spin
171 Alm<xcomplex<double> > &alm1,
Alm<xcomplex<double> > &alm2,
172 int spin,
const arr<double> &weight,
bool add_alm);
174 template<
typename T>
void alm2map_spin_adjoint
176 Alm<xcomplex<T> > &alm1,
Alm<xcomplex<T> > &alm2,
177 int spin,
bool add_alm)
180 "alm2map_spin_adjoint: maps must be in RING scheme");
182 "alm2map_spin_adjoint: maps are not conformable");
184 "alm2map_spin_adjoint: a_lm are not conformable");
186 "map contains undefined pixels");
187 checkLmaxNside(alm1.Lmax(), map1.
Nside());
190 job.set_Healpix_geometry (map1.
Nside());
191 job.set_triangular_alm_info (alm1.Lmax(), alm1.Mmax());
192 job.alm2map_spin_adjoint(&map1[0],&map2[0],&alm1(0,0),&alm2(0,0),spin,
196 template void alm2map_spin_adjoint
198 Alm<xcomplex<float> > &alm1,
Alm<xcomplex<float> > &alm2,
199 int spin,
bool add_alm);
200 template void alm2map_spin_adjoint
202 Alm<xcomplex<double> > &alm1,
Alm<xcomplex<double> > &alm2,
203 int spin,
bool add_alm);
205 template<
typename T>
void map2alm_spin_iter2
207 Alm<xcomplex<T> > &alm1,
Alm<xcomplex<T> > &alm2,
208 int spin,
double err_abs,
double err_rel)
213 alm1.SetToZero(); alm2.SetToZero();
216 map2alm_spin(map1b,map2b,alm1,alm2,spin,wgt,
true);
217 alm2map_spin(alm1,alm2,map1b,map2b,spin);
219 for (
int m=0; m<map1.
Npix(); ++m)
221 double err = abs(map1[m]-map1b[m]);
222 double rel = (map1[m]!=0) ? abs(err/map1[m]) : 1e300;
223 errmeasure = max(errmeasure,min(err/err_abs,rel/err_rel));
224 map1b[m] = map1[m]-map1b[m];
225 err = abs(map2[m]-map2b[m]);
226 rel = (map2[m]!=0) ? abs(err/map2[m]) : 1e300;
227 errmeasure = max(errmeasure,min(err/err_abs,rel/err_rel));
228 map2b[m] = map2[m]-map2b[m];
231 if (errmeasure<1)
break;
235 template void map2alm_spin_iter2
237 Alm<xcomplex<double> > &alm1,
Alm<xcomplex<double> > &alm2,
238 int spin,
double err_abs,
double err_rel);
244 Alm<xcomplex<T> > &almT,
245 Alm<xcomplex<T> > &almG,
246 Alm<xcomplex<T> > &almC,
251 "map2alm_pol: maps must be in RING scheme");
253 "map2alm_pol: maps are not conformable");
254 planck_assert (almT.conformable(almG) && almT.conformable(almC),
255 "map2alm_pol: a_lm are not conformable");
257 "map2alm_pol: weight array has too few entries");
259 "map contains undefined pixels");
260 checkLmaxNside(almT.Lmax(), mapT.
Nside());
263 job.set_weighted_Healpix_geometry (mapT.
Nside(),&weight[0]);
264 job.set_triangular_alm_info (almT.Lmax(), almT.Mmax());
265 job.map2alm(&mapT[0], &almT(0,0), add_alm);
266 job.map2alm_spin(&mapQ[0],&mapU[0],&almG(0,0),&almC(0,0),2,add_alm);
273 Alm<xcomplex<float> > &almT,
274 Alm<xcomplex<float> > &almG,
275 Alm<xcomplex<float> > &almC,
282 Alm<xcomplex<double> > &almT,
283 Alm<xcomplex<double> > &almG,
284 Alm<xcomplex<double> > &almC,
288 template<
typename T>
void alm2map_pol_adjoint
292 Alm<xcomplex<T> > &almT,
293 Alm<xcomplex<T> > &almG,
294 Alm<xcomplex<T> > &almC,
298 "alm2map_pol_adjoint: maps must be in RING scheme");
300 "alm2map_pol_adjoint: maps are not conformable");
301 planck_assert (almT.conformable(almG) && almT.conformable(almC),
302 "alm2map_pol_adjoint: a_lm are not conformable");
304 "map contains undefined pixels");
305 checkLmaxNside(almT.Lmax(), mapT.
Nside());
308 job.set_Healpix_geometry (mapT.
Nside());
309 job.set_triangular_alm_info (almT.Lmax(), almT.Mmax());
310 job.alm2map_adjoint(&mapT[0], &almT(0,0), add_alm);
311 job.alm2map_spin_adjoint(&mapQ[0],&mapU[0],&almG(0,0),&almC(0,0),2,add_alm);
314 template void alm2map_pol_adjoint
318 Alm<xcomplex<float> > &almT,
319 Alm<xcomplex<float> > &almG,
320 Alm<xcomplex<float> > &almC,
322 template void alm2map_pol_adjoint
326 Alm<xcomplex<double> > &almT,
327 Alm<xcomplex<double> > &almG,
328 Alm<xcomplex<double> > &almC,
335 Alm<xcomplex<T> > &almT,
336 Alm<xcomplex<T> > &almG,
337 Alm<xcomplex<T> > &almC,
342 for (
int iter=1; iter<=num_iter; ++iter)
349 for (
int m=0; m<mapT.
Npix(); ++m)
351 mapT2[m] = mapT[m]-mapT2[m];
352 mapQ2[m] = mapQ[m]-mapQ2[m];
353 mapU2[m] = mapU[m]-mapU2[m];
355 map2alm_pol(mapT2,mapQ2,mapU2,almT,almG,almC,weight,
true);
363 Alm<xcomplex<float> > &almT,
364 Alm<xcomplex<float> > &almG,
365 Alm<xcomplex<float> > &almC,
372 Alm<xcomplex<double> > &almT,
373 Alm<xcomplex<double> > &almG,
374 Alm<xcomplex<double> > &almC,
378 template<
typename T>
void map2alm_pol_iter2
382 Alm<xcomplex<T> > &almT,
383 Alm<xcomplex<T> > &almG,
384 Alm<xcomplex<T> > &almC,
385 double err_abs,
double err_rel)
390 almT.SetToZero(); almG.SetToZero(); almC.SetToZero();
393 map2alm_pol(mapT2,mapQ2,mapU2,almT,almG,almC,wgt,
true);
396 for (
int m=0; m<mapT.
Npix(); ++m)
398 double err = abs(mapT[m]-mapT2[m]);
399 double rel = (mapT[m]!=0) ? abs(err/mapT[m]) : 1e300;
400 errmeasure = max(errmeasure,min(err/err_abs,rel/err_rel));
401 mapT2[m] = mapT[m]-mapT2[m];
402 err = abs(mapQ[m]-mapQ2[m]);
403 rel = (mapQ[m]!=0) ? abs(err/mapQ[m]) : 1e300;
404 errmeasure = max(errmeasure,min(err/err_abs,rel/err_rel));
405 mapQ2[m] = mapQ[m]-mapQ2[m];
406 err = abs(mapU[m]-mapU2[m]);
407 rel = (mapU[m]!=0) ? abs(err/mapU[m]) : 1e300;
408 errmeasure = max(errmeasure,min(err/err_abs,rel/err_rel));
409 mapU2[m] = mapU[m]-mapU2[m];
412 if (errmeasure<1)
break;
416 template void map2alm_pol_iter2
420 Alm<xcomplex<double> > &almT,
421 Alm<xcomplex<double> > &almG,
422 Alm<xcomplex<double> > &almC,
423 double err_abs,
double err_rel);
426 template<
typename T>
void alm2map (
const Alm<xcomplex<T> > &alm,
432 job.set_Healpix_geometry (map.
Nside());
433 job.set_triangular_alm_info (alm.Lmax(), alm.Mmax());
434 job.alm2map(&alm(0,0), &map[0], add_map);
437 template void alm2map (
const Alm<xcomplex<double> > &alm,
439 template void alm2map (
const Alm<xcomplex<float> > &alm,
442 template<
typename T>
void alm2map_spin
443 (
const Alm<xcomplex<T> > &alm1,
const Alm<xcomplex<T> > &alm2,
447 "alm2map_spin: maps must be in RING scheme");
449 "alm2map_spin: maps are not conformable");
451 "alm2map_spin: a_lm are not conformable");
454 job.set_Healpix_geometry (map1.
Nside());
455 job.set_triangular_alm_info (alm1.Lmax(), alm1.Mmax());
456 job.alm2map_spin(&alm1(0,0),&alm2(0,0),&map1[0],&map2[0],spin,add_map);
459 template void alm2map_spin
460 (
const Alm<xcomplex<double> > &alm1,
const Alm<xcomplex<double> > &alm2,
462 template void alm2map_spin
463 (
const Alm<xcomplex<float> > &alm1,
const Alm<xcomplex<float> > &alm2,
468 (
const Alm<xcomplex<T> > &almT,
469 const Alm<xcomplex<T> > &almG,
470 const Alm<xcomplex<T> > &almC,
477 "alm2map_pol: maps must be in RING scheme");
479 "alm2map_pol: maps are not conformable");
480 planck_assert (almT.conformable(almG) && almT.conformable(almC),
481 "alm2map_pol: a_lm are not conformable");
484 job.set_Healpix_geometry (mapT.
Nside());
485 job.set_triangular_alm_info (almT.Lmax(), almT.Mmax());
486 job.alm2map(&almT(0,0), &mapT[0], add_map);
487 job.alm2map_spin(&almG(0,0), &almC(0,0), &mapQ[0], &mapU[0], 2, add_map);
491 const Alm<xcomplex<double> > &almG,
492 const Alm<xcomplex<double> > &almC,
499 const Alm<xcomplex<float> > &almG,
500 const Alm<xcomplex<float> > &almC,
508 (
const Alm<xcomplex<T> > &alm,
514 "alm2map_der1: maps must be in RING scheme");
516 "alm2map_der1: maps are not conformable");
519 job.set_Healpix_geometry (map.
Nside());
520 job.set_triangular_alm_info (alm.Lmax(), alm.Mmax());
521 job.alm2map(&alm(0,0), &map[0],
false);
522 job.alm2map_der1(&alm(0,0), &mapdth[0], &mapdph[0],
false);
bool fullyDefined() const
void alm2map_der1(const Alm< xcomplex< T > > &alm, Healpix_Map< T > &map, Healpix_Map< T > &mapdth, Healpix_Map< T > &mapdph)
void fill(const double &val)
void map2alm_pol(const Healpix_Map< T > &mapT, const Healpix_Map< T > &mapQ, const Healpix_Map< T > &mapU, Alm< xcomplex< T > > &almT, Alm< xcomplex< T > > &almG, Alm< xcomplex< T > > &almC, const arr< double > &weight, bool add_alm)
void map2alm(const Healpix_Map< T > &map, Alm< xcomplex< T > > &alm, const arr< double > &weight, bool add_alm)
void map2alm_pol_iter(const Healpix_Map< T > &mapT, const Healpix_Map< T > &mapQ, const Healpix_Map< T > &mapU, Alm< xcomplex< T > > &almT, Alm< xcomplex< T > > &almG, Alm< xcomplex< T > > &almC, int num_iter, const arr< double > &weight)
void alm2map_pol(const Alm< xcomplex< T > > &almT, const Alm< xcomplex< T > > &almG, const Alm< xcomplex< T > > &almC, Healpix_Map< T > &mapT, Healpix_Map< T > &mapQ, Healpix_Map< T > &mapU, bool add_map)
void alm2map_adjoint(const Healpix_Map< T > &map, Alm< xcomplex< T > > &alm, bool add_alm)
Healpix_Ordering_Scheme Scheme() const
void map2alm_iter(const Healpix_Map< T > &map, Alm< xcomplex< T > > &alm, int num_iter, const arr< double > &weight)
void alm2map(const Alm< xcomplex< T > > &alm, Healpix_Map< T > &map, bool add_map)
#define planck_assert(testval, msg)
bool conformable(const T_Healpix_Base &other) const