00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "alm_powspec_tools.h"
00033 #include "alm.h"
00034 #include "planck_rng.h"
00035 #include "powspec.h"
00036 #include "xcomplex.h"
00037 #include "rotmatrix.h"
00038 #include "openmp_support.h"
00039
00040 using namespace std;
00041 template<typename T> void create_alm
00042 (const PowSpec &powspec, Alm<xcomplex<T> > &alm, planck_rng &rng)
00043 {
00044 int lmax = alm.Lmax();
00045 int mmax = alm.Mmax();
00046 const double hsqrt2 = 1/sqrt(2.);
00047
00048 for (int l=0; l<=lmax; ++l)
00049 {
00050 double rms_tt = sqrt(powspec.tt(l));
00051 double zeta1_r = rng.rand_gauss();
00052 alm(l,0) = zeta1_r * rms_tt;
00053 for (int m=1; m<=min(l,mmax); ++m)
00054 {
00055 zeta1_r = rng.rand_gauss()*hsqrt2;
00056 double zeta1_i = rng.rand_gauss()*hsqrt2;
00057 alm(l,m).Set (zeta1_r*rms_tt, zeta1_i*rms_tt);
00058 }
00059 }
00060 }
00061
00062 template void create_alm (const PowSpec &powspec,
00063 Alm<xcomplex<float> > &alm, planck_rng &rng);
00064 template void create_alm (const PowSpec &powspec,
00065 Alm<xcomplex<double> > &alm, planck_rng &rng);
00066
00067
00068 template<typename T> void create_alm_pol
00069 (const PowSpec &powspec,
00070 Alm<xcomplex<T> > &almT,
00071 Alm<xcomplex<T> > &almG,
00072 Alm<xcomplex<T> > &almC,
00073 planck_rng &rng)
00074 {
00075 int lmax = almT.Lmax();
00076 int mmax = almT.Mmax();
00077 const double hsqrt2 = 1/sqrt(2.);
00078
00079 for (int l=0; l<=lmax; ++l)
00080 {
00081 double rms_tt=0, rms_g1=0;
00082 if (powspec.tt(l) != 0)
00083 {
00084 rms_tt = sqrt(powspec.tt(l));
00085 rms_g1 = powspec.tg(l)/rms_tt;
00086 }
00087
00088 double zeta1_r = rng.rand_gauss();
00089 almT(l,0) = zeta1_r * rms_tt;
00090 almG(l,0) = zeta1_r * rms_g1;
00091 for (int m=1; m<=min(l,mmax); ++m)
00092 {
00093 zeta1_r = rng.rand_gauss()*hsqrt2;
00094 double zeta1_i = rng.rand_gauss()*hsqrt2;
00095 almT(l,m).Set (zeta1_r*rms_tt, zeta1_i*rms_tt);
00096 almG(l,m).Set (zeta1_r*rms_g1, zeta1_i*rms_g1);
00097 }
00098 }
00099
00100 for (int l=0; l<=lmax; ++l)
00101 {
00102 double rms_g2 = 0;
00103 double rms_cc = 0;
00104 if (powspec.tt(l) != 0)
00105 {
00106 rms_g2 = powspec.gg(l) - (powspec.tg(l)/powspec.tt(l))*powspec.tg(l);
00107 if (rms_g2 <= 0)
00108 {
00109 planck_assert (abs(rms_g2) <= 1e-8*abs(powspec.gg(l)),
00110 "Inconsistent TT, GG and TG spectra at l="+dataToString(l));
00111 rms_g2 = 0;
00112 }
00113 rms_g2 = sqrt(rms_g2);
00114 rms_cc = sqrt(powspec.cc(l));
00115 }
00116 almG(l,0) += rng.rand_gauss()*rms_g2;
00117 almC(l,0) = rng.rand_gauss()*rms_cc;
00118
00119 for (int m=1; m<=min(l,mmax); ++m)
00120 {
00121 double zeta2_r = rng.rand_gauss()*hsqrt2;
00122 double zeta2_i = rng.rand_gauss()*hsqrt2;
00123 double zeta3_r = rng.rand_gauss()*hsqrt2;
00124 double zeta3_i = rng.rand_gauss()*hsqrt2;
00125
00126 almG(l,m) += xcomplex<T> (zeta2_r*rms_g2,zeta2_i*rms_g2);
00127 almC(l,m).Set (zeta3_r*rms_cc,zeta3_i*rms_cc);
00128 }
00129 }
00130 }
00131
00132 template void create_alm_pol
00133 (const PowSpec &powspec,
00134 Alm<xcomplex<float> > &almT,
00135 Alm<xcomplex<float> > &almG,
00136 Alm<xcomplex<float> > &almC,
00137 planck_rng &rng);
00138 template void create_alm_pol
00139 (const PowSpec &powspec,
00140 Alm<xcomplex<double> > &almT,
00141 Alm<xcomplex<double> > &almG,
00142 Alm<xcomplex<double> > &almC,
00143 planck_rng &rng);
00144
00145
00146 template<typename T> void extract_powspec
00147 (const Alm<xcomplex<T> > &alm, PowSpec &powspec)
00148 {
00149 arr<double> tt(alm.Lmax()+1);
00150 for (int l=0; l<=alm.Lmax(); ++l)
00151 {
00152 tt[l] = norm(alm(l,0));
00153 int limit = min(l,alm.Mmax());
00154 for (int m=1; m<=limit; ++m)
00155 tt[l] += 2*norm(alm(l,m));
00156 tt[l] /= (2*l+1);
00157 }
00158 powspec.Set(tt);
00159 }
00160
00161 template void extract_powspec
00162 (const Alm<xcomplex<float> > &alm, PowSpec &powspec);
00163 template void extract_powspec
00164 (const Alm<xcomplex<double> > &alm, PowSpec &powspec);
00165
00166
00167 template<typename T> void extract_crosspowspec
00168 (const Alm<xcomplex<T> > &alm1,
00169 const Alm<xcomplex<T> > &alm2,PowSpec &powspec)
00170 {
00171 planck_assert (alm1.conformable(alm2),
00172 "extract_crosspowspec: a_lms are not conformable");
00173 arr<double> tt(alm1.Lmax()+1);
00174 for (int l=0; l<=alm1.Lmax(); ++l)
00175 {
00176 tt[l] = alm1(l,0).re*alm2(l,0).re;
00177 int limit = min(l,alm1.Mmax());
00178 for (int m=1; m<=limit; ++m)
00179 tt[l] += 2 * (alm1(l,m).re*alm2(l,m).re + alm1(l,m).im*alm2(l,m).im);
00180 tt[l] /= (2*l+1);
00181 }
00182 powspec.Set(tt);
00183 }
00184
00185 template void extract_crosspowspec
00186 (const Alm<xcomplex<float> > &alm1,
00187 const Alm<xcomplex<float> > &alm2, PowSpec &powspec);
00188 template void extract_crosspowspec
00189 (const Alm<xcomplex<double> > &alm1,
00190 const Alm<xcomplex<double> > &alm2, PowSpec &powspec);
00191
00192
00193 template<typename T> void extract_powspec
00194 (const Alm<xcomplex<T> > &almT,
00195 const Alm<xcomplex<T> > &almG,
00196 const Alm<xcomplex<T> > &almC,
00197 PowSpec &powspec)
00198 {
00199 planck_assert (almT.conformable(almG) && almT.conformable(almC),
00200 "extract_powspec: a_lms are not conformable");
00201 int lmax = almT.Lmax();
00202 arr<double> tt(lmax+1), gg(lmax+1), cc(lmax+1), tg(lmax+1);
00203 for (int l=0; l<=lmax; ++l)
00204 {
00205 tt[l] = norm(almT(l,0));
00206 gg[l] = norm(almG(l,0));
00207 cc[l] = norm(almC(l,0));
00208 tg[l] = (almT(l,0)*conj(almG(l,0))).re;
00209 int limit = min(l,almT.Mmax());
00210 for (int m=1; m<=limit; ++m)
00211 {
00212 tt[l] += 2*norm(almT(l,m));
00213 gg[l] += 2*norm(almG(l,m));
00214 cc[l] += 2*norm(almC(l,m));
00215 tg[l] += 2*(almT(l,m)*conj(almG(l,m))).re;
00216 }
00217 tt[l] /= (2*l+1);
00218 gg[l] /= (2*l+1);
00219 cc[l] /= (2*l+1);
00220 tg[l] /= (2*l+1);
00221 }
00222 powspec.Set(tt,gg,cc,tg);
00223 }
00224
00225 template void extract_powspec
00226 (const Alm<xcomplex<float> > &almT,
00227 const Alm<xcomplex<float> > &almG,
00228 const Alm<xcomplex<float> > &almC,
00229 PowSpec &powspec);
00230 template void extract_powspec
00231 (const Alm<xcomplex<double> > &almT,
00232 const Alm<xcomplex<double> > &almG,
00233 const Alm<xcomplex<double> > &almC,
00234 PowSpec &powspec);
00235
00236
00237 template<typename T> void smooth_with_Gauss
00238 (Alm<xcomplex<T> > &alm, double fwhm_arcmin)
00239 {
00240 int fct = (fwhm_arcmin>=0) ? 1 : -1;
00241 double sigma = fwhm_arcmin/60*degr2rad*fwhm2sigma;
00242 arr<double> gb(alm.Lmax()+1);
00243 for (int l=0; l<=alm.Lmax(); ++l)
00244 gb[l] = exp(-.5*fct*l*(l+1)*sigma*sigma);
00245 alm.ScaleL(gb);
00246 }
00247
00248 template void smooth_with_Gauss
00249 (Alm<xcomplex<float> > &alm, double fwhm_arcmin);
00250 template void smooth_with_Gauss
00251 (Alm<xcomplex<double> > &alm, double fwhm_arcmin);
00252
00253
00254 template<typename T> void smooth_with_Gauss
00255 (Alm<xcomplex<T> > &almT,
00256 Alm<xcomplex<T> > &almG,
00257 Alm<xcomplex<T> > &almC,
00258 double fwhm_arcmin)
00259 {
00260 int fct = (fwhm_arcmin>=0) ? 1 : -1;
00261 double sigma = fwhm_arcmin/60*degr2rad*fwhm2sigma;
00262 double fact_pol = exp(2*fct*sigma*sigma);
00263 arr<double> gb(almT.Lmax()+1);
00264 for (int l=0; l<=almT.Lmax(); ++l)
00265 gb[l] = exp(-.5*fct*l*(l+1)*sigma*sigma);
00266 almT.ScaleL(gb);
00267 for (int l=0; l<=almT.Lmax(); ++l)
00268 gb[l] *= fact_pol;
00269 almG.ScaleL(gb); almC.ScaleL(gb);
00270 }
00271
00272 template void smooth_with_Gauss
00273 (Alm<xcomplex<float> > &almT,
00274 Alm<xcomplex<float> > &almG,
00275 Alm<xcomplex<float> > &almC,
00276 double fwhm_arcmin);
00277 template void smooth_with_Gauss
00278 (Alm<xcomplex<double> > &almT,
00279 Alm<xcomplex<double> > &almG,
00280 Alm<xcomplex<double> > &almC,
00281 double fwhm_arcmin);
00282
00283 namespace {
00284
00285 #if 0
00286
00287 class wigner_d
00288 {
00289 private:
00290 double p,q;
00291 arr<double> sqt;
00292 arr2<double> d;
00293 int n;
00294
00295 void do_line0 (double *l1, int j)
00296 {
00297 double xj = 1./j;
00298 l1[j] = -p*l1[j-1];
00299 for (int i=j-1; i>=1; --i)
00300 l1[i] = xj*sqt[j]*(q*sqt[j-i]*l1[i] - p*sqt[i]*l1[i-1]);
00301 l1[0] = q*l1[0];
00302 }
00303 void do_line (const double *l1, double *l2, int j, int k)
00304 {
00305 double xj = 1./j;
00306 double t1 = xj*sqt[j-k]*q, t2 = xj*sqt[j-k]*p;
00307 double t3 = xj*sqt[k]*p, t4 = xj*sqt[k]*q;
00308 l2[j] = sqt[j] * (t4*l1[j-1]-t2*l2[j-1]);
00309 for (int i=j-1; i>=1; --i)
00310 l2[i] = t1*sqt[j-i]*l2[i] - t2*sqt[i]*l2[i-1]
00311 +t3*sqt[j-i]*l1[i] + t4*sqt[i]*l1[i-1];
00312 l2[0] = sqt[j] * (t3*l1[0]+t1*l2[0]);
00313 }
00314
00315 public:
00316 wigner_d(int lmax, double ang)
00317 : p(sin(ang/2)), q(cos(ang/2)), sqt(2*lmax+1), d(lmax+1,2*lmax+1), n(-1)
00318 { for (int m=0; m<sqt.size(); ++m) sqt[m] = sqrt(double(m)); }
00319
00320 const arr2<double> &recurse ()
00321 {
00322 ++n;
00323 if (n==0)
00324 d[0][0] = 1;
00325 else if (n==1)
00326 {
00327 d[0][0] = q*q; d[0][1] = -p*q*sqt[2]; d[0][2] = p*p;
00328 d[1][0] = -d[0][1]; d[1][1] = q*q-p*p; d[1][2] = d[0][1];
00329 }
00330 else
00331 {
00332
00333 int sign = (n&1)? -1 : 1;
00334 for (int i=0; i<=2*n-2; ++i)
00335 {
00336 d[n][i] = sign*d[n-2][2*n-2-i];
00337 sign=-sign;
00338 }
00339 do_line (d[n-1],d[n],2*n-1,n);
00340 for (int k=n; k>=2; --k)
00341 {
00342 do_line (d[k-2],d[k-1],2*n-1,k-1);
00343 do_line (d[k-1],d[k],2*n,k);
00344 }
00345 do_line0 (d[0],2*n-1);
00346 do_line (d[0],d[1],2*n,1);
00347 do_line0 (d[0],2*n);
00348 }
00349 return d;
00350 }
00351 };
00352
00353 #else
00354
00355 class wigner_d
00356 {
00357 private:
00358 double p,q;
00359 arr<double> sqt;
00360 arr2<double> d, dd;
00361 int n;
00362
00363 public:
00364 wigner_d(int lmax, double ang)
00365 : p(sin(ang/2)), q(cos(ang/2)), sqt(2*lmax+1), d(lmax+1,2*lmax+1),
00366 dd(lmax+1,2*lmax+1), n(-1)
00367 { for (int m=0; m<sqt.size(); ++m) sqt[m] = sqrt(double(m)); }
00368
00369 const arr2<double> &recurse ()
00370 {
00371 ++n;
00372 if (n==0)
00373 d[0][0] = 1;
00374 else if (n==1)
00375 {
00376 d[0][0] = q*q; d[0][1] = -p*q*sqt[2]; d[0][2] = p*p;
00377 d[1][0] = -d[0][1]; d[1][1] = q*q-p*p; d[1][2] = d[0][1];
00378 }
00379 else
00380 {
00381
00382 int sign = (n&1)? -1 : 1;
00383 for (int i=0; i<=2*n-2; ++i)
00384 {
00385 d[n][i] = sign*d[n-2][2*n-2-i];
00386 sign=-sign;
00387 }
00388 for (int j=2*n-1; j<=2*n; ++j)
00389 {
00390 double xj = 1./j;
00391 dd[0][0] = q*d[0][0];
00392 for (int i=1;i<j; ++i)
00393 dd[0][i] = xj*sqt[j]*(q*sqt[j-i]*d[0][i] - p*sqt[i]*d[0][i-1]);
00394 dd[0][j] = -p*d[0][j-1];
00395 #pragma omp parallel
00396 {
00397 int k;
00398 #pragma omp for schedule(static)
00399 for (k=1; k<=n; ++k)
00400 {
00401 double t1 = xj*sqt[j-k]*q, t2 = xj*sqt[j-k]*p;
00402 double t3 = xj*sqt[k ]*p, t4 = xj*sqt[k ]*q;
00403 dd[k][0] = xj*sqt[j]*(q*sqt[j-k]*d[k][0] + p*sqt[k]*d[k-1][0]);
00404 for (int i=1; i<j; ++i)
00405 dd[k][i] = t1*sqt[j-i]*d[k ][i] - t2*sqt[i]*d[k ][i-1]
00406 + t3*sqt[j-i]*d[k-1][i] + t4*sqt[i]*d[k-1][i-1];
00407 dd[k][j] = -t2*sqt[j]*d[k][j-1] + t4*sqt[j]*d[k-1][j-1];
00408 }
00409 }
00410 dd.swap(d);
00411 }
00412 }
00413 return d;
00414 }
00415 };
00416
00417 #endif
00418
00419 }
00420
00421 template<typename T> void rotate_alm (Alm<xcomplex<T> > &alm,
00422 double psi, double theta, double phi)
00423 {
00424 planck_assert (alm.Lmax()==alm.Mmax(),
00425 "rotate_alm: lmax must be equal to mmax");
00426 int lmax=alm.Lmax();
00427 arr<xcomplex<double> > exppsi(lmax+1), expphi(lmax+1);
00428 for (int m=0; m<=lmax; ++m)
00429 {
00430 exppsi[m].Set (cos(psi*m),-sin(psi*m));
00431 expphi[m].Set (cos(phi*m),-sin(phi*m));
00432 }
00433
00434 wigner_d rec(lmax,theta);
00435
00436 arr<xcomplex<double> > almtmp(lmax+1);
00437
00438 for (int l=0; l<=lmax; ++l)
00439 {
00440 announce_progress (pow(double(l),3),pow(l-1.,3),pow(double(lmax),3));
00441 const arr2<double> &d(rec.recurse());
00442
00443 for (int m=0; m<=l; ++m)
00444 almtmp[m] = alm(l,0)*d[l][l+m];
00445
00446 #pragma omp parallel
00447 {
00448 int lo,hi;
00449 openmp_calc_share(0,l+1,lo,hi);
00450
00451 bool flip = true;
00452 for (int mm=1; mm<=l; ++mm)
00453 {
00454 xcomplex<double> t1 = alm(l,mm)*exppsi[mm];
00455 bool flip2 = ((mm+lo)&1) ? true : false;
00456 for (int m=lo; m<hi; ++m)
00457 {
00458 double d1 = flip2 ? -d[l-mm][l-m] : d[l-mm][l-m];
00459 double d2 = flip ? -d[l-mm][l+m] : d[l-mm][l+m];
00460 double f1 = d1+d2, f2 = d1-d2;
00461 almtmp[m].re += t1.re*f1; almtmp[m].im += t1.im*f2;
00462 flip2 = !flip2;
00463 }
00464 flip = !flip;
00465 }
00466 }
00467
00468 for (int m=0; m<=l; ++m)
00469 alm(l,m) = almtmp[m]*expphi[m];
00470 }
00471 end_announce_progress();
00472 }
00473
00474 template void rotate_alm (Alm<xcomplex<float> > &alm,
00475 double psi, double theta, double phi);
00476 template void rotate_alm (Alm<xcomplex<double> > &alm,
00477 double psi, double theta, double phi);
00478
00479 template<typename T> void rotate_alm (Alm<xcomplex<T> > &almT,
00480 Alm<xcomplex<T> > &almG, Alm<xcomplex<T> > &almC,
00481 double psi, double theta, double phi)
00482 {
00483 planck_assert (almT.Lmax()==almT.Mmax(),
00484 "rotate_alm: lmax must be equal to mmax");
00485 planck_assert (almG.conformable(almT) && almC.conformable(almT),
00486 "rotate_alm: a_lm are not conformable");
00487 int lmax=almT.Lmax();
00488 arr<xcomplex<double> > exppsi(lmax+1), expphi(lmax+1);
00489 for (int m=0; m<=lmax; ++m)
00490 {
00491 exppsi[m].Set (cos(psi*m),-sin(psi*m));
00492 expphi[m].Set (cos(phi*m),-sin(phi*m));
00493 }
00494
00495 wigner_d rec(lmax,theta);
00496
00497 arr<xcomplex<double> > almtmpT(lmax+1), almtmpG(lmax+1), almtmpC(lmax+1);
00498
00499 for (int l=0; l<=lmax; ++l)
00500 {
00501 announce_progress (pow(double(l),3),pow(l-1.,3),pow(double(lmax),3));
00502 const arr2<double> &d(rec.recurse());
00503
00504 for (int m=0; m<=l; ++m)
00505 {
00506 almtmpT[m] = almT(l,0)*d[l][m+l];
00507 almtmpG[m] = almG(l,0)*d[l][m+l];
00508 almtmpC[m] = almC(l,0)*d[l][m+l];
00509 }
00510
00511 #pragma omp parallel
00512 {
00513 int lo,hi;
00514 openmp_calc_share(0,l+1,lo,hi);
00515
00516 bool flip = true;
00517 for (int mm=1; mm<=l; ++mm)
00518 {
00519 xcomplex<double> t1T = almT(l,mm)*exppsi[mm];
00520 xcomplex<double> t1G = almG(l,mm)*exppsi[mm];
00521 xcomplex<double> t1C = almC(l,mm)*exppsi[mm];
00522 bool flip2 = ((mm+lo)&1) ? true : false;
00523 for (int m=lo; m<hi; ++m)
00524 {
00525 double d1 = flip2 ? -d[l-mm][l-m] : d[l-mm][l-m];
00526 double d2 = flip ? -d[l-mm][l+m] : d[l-mm][l+m];
00527 double f1 = d1+d2, f2 = d1-d2;
00528 almtmpT[m].re += t1T.re*f1; almtmpT[m].im += t1T.im*f2;
00529 almtmpG[m].re += t1G.re*f1; almtmpG[m].im += t1G.im*f2;
00530 almtmpC[m].re += t1C.re*f1; almtmpC[m].im += t1C.im*f2;
00531 flip2 = !flip2;
00532 }
00533 flip = !flip;
00534 }
00535 }
00536
00537 for (int m=0; m<=l; ++m)
00538 {
00539 almT(l,m) = almtmpT[m]*expphi[m];
00540 almG(l,m) = almtmpG[m]*expphi[m];
00541 almC(l,m) = almtmpC[m]*expphi[m];
00542 }
00543 }
00544 end_announce_progress();
00545 }
00546
00547 template void rotate_alm (Alm<xcomplex<float> > &almT,
00548 Alm<xcomplex<float> > &almG, Alm<xcomplex<float> > &almC,
00549 double psi, double theta, double phi);
00550 template void rotate_alm (Alm<xcomplex<double> > &almT,
00551 Alm<xcomplex<double> > &almG, Alm<xcomplex<double> > &almC,
00552 double psi, double theta, double phi);
00553
00554
00555 template<typename T> void rotate_alm (Alm<xcomplex<T> > &alm,
00556 const rotmatrix &mat)
00557 {
00558 double a1, a2, a3;
00559 mat.Extract_CPAC_Euler_Angles (a1, a2, a3);
00560 rotate_alm (alm, a3, a2, a1);
00561 }
00562
00563 template void rotate_alm (Alm<xcomplex<float> > &alm, const rotmatrix &mat);
00564 template void rotate_alm (Alm<xcomplex<double> > &alm, const rotmatrix &mat);
00565
00566 template<typename T> void rotate_alm (Alm<xcomplex<T> > &almT,
00567 Alm<xcomplex<T> > &almG, Alm<xcomplex<T> > &almC,
00568 const rotmatrix &mat)
00569 {
00570 double a1, a2, a3;
00571 mat.Extract_CPAC_Euler_Angles (a1, a2, a3);
00572 rotate_alm (almT, almG, almC, a3, a2, a1);
00573 }
00574
00575 template void rotate_alm (Alm<xcomplex<float> > &almT,
00576 Alm<xcomplex<float> > &almG, Alm<xcomplex<float> > &almC,
00577 const rotmatrix &mat);
00578 template void rotate_alm (Alm<xcomplex<double> > &almT,
00579 Alm<xcomplex<double> > &almG, Alm<xcomplex<double> > &almC,
00580 const rotmatrix &mat);
