source: Sophya/trunk/SophyaLib/Samba/sphericaltransformserver.cc@ 1715

Last change on this file since 1715 was 1715, checked in by ansari, 24 years ago

compil avec SGI-CC (redeclaration k for(int k ...)) - Reza 22/10/2001
File size: 42.1 KB
Line 
1#include "machdefs.h"
2#include <iostream.h>
3#include <math.h>
4#include <complex>
5#include "sphericaltransformserver.h"
6#include "tvector.h"
7#include "nbrandom.h"
8#include "nbmath.h"
9#include "timing.h"
10//#include "spherehealpix.h"
11
12
13/*! \class SOPHYA::SphericalTransformServer
14
15 Class for performing analysis and synthesis of sky maps using spin-0 or spin-2 spherical harmonics.
16
17Maps must be SOPHYA SphericalMaps (SphereGorski or SphereThetaPhi).
18
19Temperature and polarization (Stokes parameters) can be developped on spherical harmonics :
20\f[
21\frac{\Delta T}{T}(\hat{n})=\sum_{lm}a_{lm}^TY_l^m(\hat{n})
22\f]
23\f[
24Q(\hat{n})=\frac{1}{\sqrt{2}}\sum_{lm}N_l\left(a_{lm}^EW_{lm}(\hat{n})+a_{lm}^BX_{lm}(\hat{n})\right)
25\f]
26\f[
27U(\hat{n})=-\frac{1}{\sqrt{2}}\sum_{lm}N_l\left(a_{lm}^EX_{lm}(\hat{n})-a_{lm}^BW_{lm}(\hat{n})\right)
28\f]
29\f[
30\left(Q \pm iU\right)(\hat{n})=\sum_{lm}a_{\pm 2lm}\, _{\pm 2}Y_l^m(\hat{n})
31\f]
32
33\f[
34Y_l^m(\hat{n})=\lambda_l^m(\theta)e^{im\phi}
35\f]
36\f[
37_{\pm}Y_l^m(\hat{n})=_{\pm}\lambda_l^m(\theta)e^{im\phi}
38\f]
39\f[
40W_{lm}(\hat{n})=\frac{1}{N_l}\,_{w}\lambda_l^m(\theta)e^{im\phi}
41\f]
42\f[
43X_{lm}(\hat{n})=\frac{-i}{N_l}\,_{x}\lambda_l^m(\theta)e^{im\phi}
44\f]
45
46(see LambdaLMBuilder, LambdaPMBuilder, LambdaWXBuilder classes)
47
48power spectra :
49
50\f[
51C_l^T=\frac{1}{2l+1}\sum_{m=0}^{+ \infty }\left|a_{lm}^T\right|^2=\langle\left|a_{lm}^T\right|^2\rangle
52\f]
53\f[
54C_l^E=\frac{1}{2l+1}\sum_{m=0}^{+\infty}\left|a_{lm}^E\right|^2=\langle\left|a_{lm}^E\right|^2\rangle
55\f]
56\f[
57C_l^B=\frac{1}{2l+1}\sum_{m=0}^{+\infty}\left|a_{lm}^B\right|^2=\langle\left|a_{lm}^B\right|^2\rangle
58\f]
59
60\arg
61\b Synthesis : Get temperature and polarization maps from \f$a_{lm}\f$ coefficients or from power spectra, (methods GenerateFrom...).
62
63\b Temperature:
64\f[
65\frac{\Delta T}{T}(\hat{n})=\sum_{lm}a_{lm}^TY_l^m(\hat{n}) = \sum_{-\infty}^{+\infty}b_m(\theta)e^{im\phi}
66\f]
67
68with
69\f[
70b_m(\theta)=\sum_{l=\left|m\right|}^{+\infty}a_{lm}^T\lambda_l^m(\theta)
71\f]
72
73\b Polarisation
74\f[
75Q \pm iU = \sum_{-\infty}^{+\infty}b_m^{\pm}(\theta)e^{im\phi}
76\f]
77
78where :
79\f[
80b_m^{\pm}(\theta) = \sum_{l=\left|m\right|}^{+\infty}a_{\pm 2lm}\,_{\pm}\lambda_l^m(\theta)
81\f]
82
83or :
84\f[
85Q = \sum_{-\infty}^{+\infty}b_m^{Q}(\theta)e^{im\phi}
86\f]
87\f[
88U = \sum_{-\infty}^{+\infty}b_m^{U}(\theta)e^{im\phi}
89\f]
90
91where:
92\f[
93b_m^{Q}(\theta) = \frac{1}{\sqrt{2}}\sum_{l=\left|m\right|}^{+\infty}\left(a_{lm}^E\,_{w}\lambda_l^m(\theta)-ia_{lm}^B\,_{x}\lambda_l^m(\theta)\right)
94\f]
95\f[
96b_m^{U}(\theta) = \frac{1}{\sqrt{2}}\sum_{l=\left|m\right|}^{+\infty}\left(ia_{lm}^E\,_{x}\lambda_l^m(\theta)+a_{lm}^B\,_{w}\lambda_l^m(\theta)\right)
97\f]
98
99Since the pixelization provides "slices" with constant \f$\theta\f$ and \f$\phi\f$ equally distributed on \f$2\pi\f$ \f$\frac{\Delta T}{T}\f$, \f$Q\f$,\f$U\f$ can be computed by FFT.
100
101
102\arg
103\b Analysis : Get \f$a_{lm}\f$ coefficients or power spectra from temperature and polarization maps (methods DecomposeTo...).
104
105\b Temperature:
106\f[
107a_{lm}^T=\int\frac{\Delta T}{T}(\hat{n})Y_l^{m*}(\hat{n})d\hat{n}
108\f]
109
110approximated as :
111\f[
112a_{lm}^T=\sum_{\theta_k}\omega_kC_m(\theta_k)\lambda_l^m(\theta_k)
113\f]
114where :
115\f[
116C_m (\theta _k)=\sum_{\phi _{k\prime}}\frac{\Delta T}{T}(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}
117\f]
118Since the pixelization provides "slices" with constant \f$\theta\f$ and \f$\phi\f$ equally distributed on \f$2\pi\f$ (\f$\omega_k\f$ is the solid angle of each pixel of the slice \f$\theta_k\f$) \f$C_m\f$ can be computed by FFT.
119
120\b polarisation:
121
122\f[
123a_{\pm 2lm}=\sum_{\theta_k}\omega_kC_m^{\pm}(\theta_k)\,_{\pm}\lambda_l^m(\theta_k)
124\f]
125where :
126\f[
127C_m^{\pm} (\theta _k)=\sum_{\phi _{k\prime}}\left(Q \pm iU\right)(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}
128\f]
129or :
130
131\f[
132a_{lm}^E=\frac{1}{\sqrt{2}}\sum_{\theta_k}\omega_k\left(C_m^{Q}(\theta_k)\,_{w}\lambda_l^m(\theta_k)-iC_m^{U}(\theta_k)\,_{x}\lambda_l^m(\theta_k)\right)
133\f]
134\f[
135a_{lm}^B=\frac{1}{\sqrt{2}}\sum_{\theta_k}\omega_k\left(iC_m^{Q}(\theta_k)\,_{x}\lambda_l^m(\theta_k)+C_m^{U}(\theta_k)\,_{w}\lambda_l^m(\theta_k)\right)
136\f]
137
138where :
139\f[
140C_m^{Q} (\theta _k)=\sum_{\phi _{k\prime}}Q(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}
141\f]
142\f[
143C_m^{U} (\theta _k)=\sum_{\phi _{k\prime}}U(\theta _k,\phi_{k\prime})e^{-im\phi _{k\prime}}
144\f]
145
146 */
147
148 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const
149
150 synthesis of a temperature map from Alm coefficients
151*/
152template<class T>
153void SphericalTransformServer<T>::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const
154{
155 /*=======================================================================
156 computes a map form its alm for the HEALPIX pixelisation
157 map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
158 = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
159
160 where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
161
162 * the recurrence of Ylm is the standard one (cf Num Rec)
163 * the sum over m is done by FFT
164
165 =======================================================================*/
166 int_4 nlmax=alm.Lmax();
167 int_4 nmmax=nlmax;
168 int_4 nsmax=0;
169 map.Resize(pixelSizeIndex);
170 char* sphere_type=map.TypeOfMap();
171 if (strncmp(sphere_type,"RING",4) == 0)
172 {
173 nsmax=map.SizeIndex();
174 }
175 else
176 // pour une sphere Gorski le nombre de pixels est 12*nsmax**2
177 // on calcule une quantite equivalente a nsmax pour la sphere-theta-phi
178 // en vue de l'application du critere Healpix : nlmax<=3*nsmax-1
179 // c'est approximatif ; a raffiner.
180 if (strncmp(sphere_type,"TETAFI",6) == 0)
181 {
182 nsmax=(int_4)sqrt(map.NbPixels()/12.);
183 }
184 else
185 {
186 cout << " unknown type of sphere : " << sphere_type << endl;
187 throw IOExc(" unknown type of sphere: " + (string)sphere_type );
188 }
189 // cout << "GenerateFromAlm: the sphere is of type : " << sphere_type << endl;
190 // cout << "GenerateFromAlm: size index (nside) of the sphere= " << nsmax << endl;
191 // cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
192 if (nlmax>3*nsmax-1)
193 {
194 // cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
195 if (strncmp(sphere_type,"TETAFI",6) == 0)
196 {
197 cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
198 cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
199 }
200 }
201 Bm<complex<T> > b_m_theta(nmmax);
202
203 // map.Resize(nsmax);
204
205
206 // pour chaque tranche en theta
207 for (int_4 ith = 0; ith < map.NbThetaSlices();ith++)
208 {
209 int_4 nph;
210 r_8 phi0;
211 r_8 theta;
212 TVector<int_4> pixNumber;
213 TVector<T> datan;
214
215 map.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
216 nph = pixNumber.NElts();
217
218 // -----------------------------------------------------
219 // for each theta, and each m, computes
220 // b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
221 // ------------------------------------------------------
222 LambdaLMBuilder lb(theta,nlmax,nmmax);
223 // somme sur m de 0 a l'infini
224 int m;
225 for (m = 0; m <= nmmax; m++)
226 {
227 b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
228 for (int l = m+1; l<= nlmax; l++)
229 {
230 b_m_theta(m) += (T)( lb.lamlm(l,m) ) * alm(l,m);
231 }
232 }
233
234 // obtains the negative m of b(m,theta) (= complex conjugate)
235
236 for (m=1;m<=nmmax;m++)
237 {
238 b_m_theta(-m) = conj(b_m_theta(m));
239 }
240 // ---------------------------------------------------------------
241 // sum_m b(m,theta)*exp(i*m*phi) -> f(phi,theta)
242 // ---------------------------------------------------------------*/
243 TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
244 for (int i=0;i< nph;i++)
245 {
246 map(pixNumber(i))=Temp(i);
247 }
248 }
249}
250
251
252
253 /*! \fn TVector< complex<T> > SOPHYA::SphericalTransformServer::fourierSynthesisFromB(const Bm<complex<T> >& b_m, int_4 nph, r_8 phi0) const
254
255\return a vector with nph elements which are sums :\f$\sum_{m=-mmax}^{mmax}b_m(\theta)e^{im\varphi}\f$ for nph values of \f$\varphi\f$ regularly distributed in \f$[0,\pi]\f$ ( calculated by FFT)
256
257 The object b_m (\f$b_m\f$) of the class Bm is a special vector which index goes from -mmax to mmax.
258 */
259template<class T>
260TVector< complex<T> > SphericalTransformServer<T>::fourierSynthesisFromB(const Bm<complex<T> >& b_m, int_4 nph, r_8 phi0) const
261{
262 /*=======================================================================
263 dataout(j) = sum_m datain(m) * exp(i*m*phi(j))
264 with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1
265
266 as the set of frequencies {m} is larger than nph,
267 we wrap frequencies within {0..nph-1}
268 ie m = k*nph + m' with m' in {0..nph-1}
269 then
270 noting bw(m') = exp(i*m'*phi0)
271 * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
272 with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
273 dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
274 = Fourier Transform of bw
275 is real
276
277 NB nph is not necessarily a power of 2
278
279=======================================================================*/
280 //**********************************************************************
281 // pour une valeur de phi (indexee par j) la temperature est la transformee
282 // de Fourier de bm (somme sur m de -nmax a +nmmax de bm*exp(i*m*phi)).
283 // on demande nph (nombre de pixels sur la tranche) valeurs de transformees, pour nph valeurs de phi, regulierement reparties sur 2*pi. On a:
284 // DT/T(j) = sum_m b(m) * exp(i*m*phi(j))
285 // sommation de -infini a +infini, en fait limitee a -nmamx, +nmmax
286 // On pose m=k*nph + m', avec m' compris entre 0 et nph-1. Alors :
287 // DT/T(j) = somme_k somme_m' b(k*nph + m')*exp(i*(k*nph + m')*phi(j))
288 // somme_k : de -infini a +infini
289 // somme_m' : de 0 a nph-1
290 // On echange les sommations :
291 // DT/T(j) = somme_k (exp(i*m'*phi(j)) somme_m' b(k*nph + m')*exp(i*(k*nph*phi(j))
292 // mais phi(j) est un multiple entier de 2*pi/nph, la seconde exponentielle
293 // vaut 1.
294 // Il reste a calculer les transformees de Fourier de somme_m' b(k*nph + m')
295 // si phi0 n'est pas nul, il y a juste un decalage a faire.
296 //**********************************************************************
297
298 TVector< complex<T> > bw(nph);
299 TVector< complex<T> > dataout(nph);
300 TVector< complex<T> > data(nph);
301
302
303 for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
304 int m;
305 for (m=-b_m.Mmax();m<=-1;m++)
306 {
307 int maux=m;
308 while (maux<0) maux+=nph;
309 int iw=maux%nph;
310 double aux=(m-iw)*phi0;
311 bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) ) ;
312 }
313 for (m=0;m<=b_m.Mmax();m++)
314 {
315 // int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
316 int iw=m%nph;
317 double aux=(m-iw)*phi0;
318 bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
319 }
320
321 // applies the shift in position <-> phase factor in Fourier space
322 for (int mprime=0; mprime < nph; mprime++)
323 {
324 complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
325 data(mprime)=bw(mprime)*
326 (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
327 }
328
329 //sortie.ReSize(nph);
330 TVector< complex<T> > sortie(nph);
331
332 fftIntfPtr_-> FFTBackward(data, sortie);
333
334 return sortie;
335}
336
337//********************************************
338/*! \fn TVector<T> SOPHYA::SphericalTransformServer::RfourierSynthesisFromB(const Bm<complex<T> >& b_m, int_4 nph, r_8 phi0) const
339
340same as fourierSynthesisFromB, but return a real vector, taking into account the fact that b(-m) is conjugate of b(m) */
341template<class T>
342TVector<T> SphericalTransformServer<T>::RfourierSynthesisFromB(const Bm<complex<T> >& b_m, int_4 nph, r_8 phi0) const
343{
344 /*=======================================================================
345 dataout(j) = sum_m datain(m) * exp(i*m*phi(j))
346 with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1
347
348 as the set of frequencies {m} is larger than nph,
349 we wrap frequencies within {0..nph-1}
350 ie m = k*nph + m' with m' in {0..nph-1}
351 then
352 noting bw(m') = exp(i*m'*phi0)
353 * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
354 with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
355 dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
356 = Fourier Transform of bw
357 is real
358
359 NB nph is not necessarily a power of 2
360
361=======================================================================*/
362 //**********************************************************************
363 // pour une valeur de phi (indexee par j) la temperature est la transformee
364 // de Fourier de bm (somme sur m de -nmax a +nmmax de bm*exp(i*m*phi)).
365 // on demande nph (nombre de pixels sur la tranche) valeurs de transformees, pour nph valeurs de phi, regulierement reparties sur 2*pi. On a:
366 // DT/T(j) = sum_m b(m) * exp(i*m*phi(j))
367 // sommation de -infini a +infini, en fait limitee a -nmamx, +nmmax
368 // On pose m=k*nph + m', avec m' compris entre 0 et nph-1. Alors :
369 // DT/T(j) = somme_k somme_m' b(k*nph + m')*exp(i*(k*nph + m')*phi(j))
370 // somme_k : de -infini a +infini
371 // somme_m' : de 0 a nph-1
372 // On echange les sommations :
373 // DT/T(j) = somme_k (exp(i*m'*phi(j)) somme_m' b(k*nph + m')*exp(i*(k*nph*phi(j))
374 // mais phi(j) est un multiple entier de 2*pi/nph, la seconde exponentielle
375 // vaut 1.
376 // Il reste a calculer les transformees de Fourier de somme_m' b(k*nph + m')
377 // si phi0 n'est pas nul, il y a juste un decalage a faire.
378 //**********************************************************************
379
380 TVector< complex<T> > bw(nph);
381 TVector< complex<T> > dataout(nph);
382 TVector< complex<T> > data(nph/2+1);
383
384
385 for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
386 int m;
387 for (m=-b_m.Mmax();m<=-1;m++)
388 {
389 int maux=m;
390 while (maux<0) maux+=nph;
391 int iw=maux%nph;
392 double aux=(m-iw)*phi0;
393 bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) ) ;
394 }
395 for (m=0;m<=b_m.Mmax();m++)
396 {
397 // int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
398 int iw=m%nph;
399 double aux=(m-iw)*phi0;
400 bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
401 }
402
403 // applies the shift in position <-> phase factor in Fourier space
404 for (int mprime=0; mprime <= nph/2; mprime++)
405 {
406 complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
407 data(mprime)=bw(mprime)*
408 (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
409 }
410
411 TVector<T> sortie;
412
413 fftIntfPtr_-> FFTBackward(data, sortie);
414
415 return sortie;
416}
417//*******************************************
418
419 /*! \fn Alm<T> SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const
420
421\return the Alm coefficients from analysis of a temperature map. THE MAP CAN BE MODIFIED (if iterationOrder >0)
422
423 \param<nlmax> : maximum value of the l index
424
425 \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
426
427\param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. (Default : 0 -> standard analysis)
428 */
429template<class T>
430void SphericalTransformServer<T>::DecomposeToAlm(SphericalMap<T>& map, Alm<T>& alm, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
431{
432 int_4 nmmax = nlmax;
433 // PrtTim("appel carteVersAlm");
434 carteVersAlm(map, nlmax, cos_theta_cut, alm);
435 // PrtTim("retour carteVersAlm");
436 if (iterationOrder > 0)
437 {
438 TVector<int_4> fact(iterationOrder+2);
439 fact(0) = 1;
440 int k;
441 for (k=1; k <= iterationOrder+1; k++)
442 {
443 fact(k) = fact(k-1)*k;
444 }
445 Alm<T> alm2(alm);
446 T Tzero = (T)0.;
447 complex<T> complexZero = complex<T>(Tzero, Tzero);
448 alm = complexZero;
449 int signe = 1;
450 int nbIteration = iterationOrder+1;
451 for (k=1; k <= nbIteration; k++)
452 {
453 T facMult = (T)(0.5*signe*fact(iterationOrder)*(2*nbIteration-k)/(fact(k)*fact(nbIteration-k)));
454 for (int m = 0; m <= nmmax; m++)
455 {
456 for (int l = m; l<= nlmax; l++)
457 {
458 alm(l,m) += facMult*alm2(l,m);
459 }
460 }
461 if (k == nbIteration) break;
462 signe = -signe;
463 for (int k=0; k< map.NbPixels(); k++) map(k) = (T)0.;
464 // synthetize a map from the estimated alm
465 // PrtTim("appel GenerateFromAlm");
466 GenerateFromAlm( map, map.SizeIndex(), alm2);
467 // PrtTim("retour GenerateFromAlm");
468 alm2 = complexZero;
469 // analyse the new map
470 // PrtTim("appel carteVersAlm");
471 carteVersAlm(map, nlmax, cos_theta_cut, alm2);
472 // PrtTim("retour carteVersAlm");
473 }
474 }
475}
476
477template<class T>
478 void SphericalTransformServer<T>::carteVersAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut, Alm<T>& alm) const
479{
480
481 /*-----------------------------------------------------------------------
482 computes the integral in phi : phas_m(theta)
483 for each parallele from north to south pole
484 -----------------------------------------------------------------------*/
485 TVector<T> data;
486 TVector<int_4> pixNumber;
487 int_4 nmmax = nlmax;
488 TVector< complex<T> > phase(nmmax+1);
489
490 alm.ReSizeToLmax(nlmax);
491 for (int_4 ith = 0; ith < map.NbThetaSlices(); ith++)
492 {
493 r_8 phi0;
494 r_8 theta;
495 // PrtTim("debut 1ere tranche ");
496 map.GetThetaSlice(ith,theta,phi0,pixNumber ,data);
497 phase = complex<T>((T)0.,(T)0.);
498 double cth = cos(theta);
499
500 //part of the sky out of the symetric cut
501 bool keep_it = (fabs(cth) >= cos_theta_cut);
502
503 // PrtTim("fin 1ere tranche ");
504
505 if (keep_it)
506 {
507 // phase = CFromFourierAnalysis(nmmax,data,phi0);
508 // PrtTim("avant Fourier ");
509 CFromFourierAnalysis(nmmax,data,phase, phi0);
510 // PrtTim("apres Fourier ");
511
512 }
513
514// ---------------------------------------------------------------------
515// computes the a_lm by integrating over theta
516// lambda_lm(theta) * phas_m(theta)
517// for each m and l
518// -----------------------------------------------------------------------
519 // PrtTim("avant instanciation LM ");
520 LambdaLMBuilder lb(theta,nlmax,nmmax);
521 // PrtTim("apres instanciation LM ");
522 r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
523
524 // PrtTim("avant mise a jour Alm ");
525 complex<T> fi;
526 T facteur;
527 int index;
528 for (int m = 0; m <= nmmax; m++)
529 {
530 fi = phase(m);
531 for (int l = m; l<= nlmax; l++)
532 {
533 index = alm.indexOfElement(l,m);
534 // facteur = (T)(lb.lamlm(l,m) * domega);
535 facteur = (T)(lb.lamlm(index) * domega);
536 // alm(l,m) += facteur * fi ;
537 alm(index) += facteur * fi ;
538 }
539 }
540
541
542
543 //
544 //
545 // PrtTim("apres mise a jour Alm ");
546 }
547}
548 /*! \fn TVector< complex<T> > SOPHYA::SphericalTransformServer::CFromFourierAnalysis(int_4 nmmax, const TVector<complex<T> >datain, r_8 phi0) const
549
550\return a vector with mmax elements which are sums :
551\f$\sum_{k=0}^{nphi}datain(\theta,\varphi_k)e^{im\varphi_k}\f$ for (mmax+1) values of \f$m\f$ from 0 to mmax.
552 */
553template<class T>
554TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<complex<T> >datain, r_8 phi0) const
555{
556 /*=======================================================================
557 integrates (data * phi-dependence-of-Ylm) over phi
558 --> function of m can be computed by FFT
559
560 datain est modifie
561 =======================================================================*/
562 int_4 nph=datain.NElts();
563 if (nph <= 0)
564 {
565 throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
566 }
567 TVector<complex<T> > transformedData(nph);
568 fftIntfPtr_-> FFTForward(datain, transformedData);
569
570 TVector< complex<T> > dataout(nmmax+1);
571
572 int im_max=min(nph,nmmax+1);
573 int i;
574 dataout = complex<T>((T)0.,(T)0.);
575 // for (i=0;i< dataout.NElts();i++) dataout(i)=complex<T>((T)0.,(T)0.);
576 for (i=0;i<im_max;i++) dataout(i)=transformedData(i);
577
578
579 for (int kk=nph; kk<dataout.NElts(); kk++) dataout(kk)=dataout(kk%nph);
580 for (i = 0;i <dataout.NElts();i++){
581 dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
582 }
583 return dataout;
584}
585
586//&&&&&&&&& nouvelle version
587/* \fn TVector< complex<T> > SOPHYA::SphericalTransformServer::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, r_8 phi0) const
588
589same as previous one, but with a "datain" which is real (not complex) */
590template<class T>
591void SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, TVector< complex<T> >& dataout, r_8 phi0) const
592{
593 //=======================================================================
594 // integrates (data * phi-dependence-of-Ylm) over phi
595 // --> function of m can be computed by FFT
596 // ! with 0<= m <= npoints/2 (: Nyquist)
597 // ! because the data is real the negative m are the conjugate of the
598 // ! positive ones
599
600 // datain est modifie
601 //
602 // =======================================================================
603 int_4 nph=datain.NElts();
604 if (nph <= 0)
605 {
606 throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
607 }
608 if (nph%2 != 0 )
609 {
610 throw PException("SphericalTransformServer<T>::CFromFourierAnalysis : longueur de datain impair ?");
611 }
612 TVector<complex<T> > transformedData;
613
614 // la taille du vecteur complexe retourne est nph/2+1 (si la taille
615 // du vecteur reel entre est nph)
616 fftIntfPtr_-> FFTForward(datain, transformedData);
617
618 // TVector< complex<T> > dataout(nmmax+1);
619 dataout.ReSize(nmmax+1);
620
621 // on transfere le resultat de la fft dans dataout.
622
623 int maxFreqAccessiblesParFFT = min(nph/2,nmmax);
624 int i;
625 for (i=0;i<=maxFreqAccessiblesParFFT;i++) dataout(i)=transformedData(i);
626
627
628 // si dataout n'est pas plein, on complete jusqu'a nph+1 valeurs (a moins
629 // que dataout ne soit plein avant d'atteindre nph)
630 if (maxFreqAccessiblesParFFT != nmmax )
631 {
632 int maxMfft = min(nph,nmmax);
633 for (i=maxFreqAccessiblesParFFT+1; i<=maxMfft; i++)
634 {
635 dataout(i) = conj(dataout(nph-i) );
636 }
637 // on conplete, si necessaire, par periodicite
638 if ( maxMfft != nmmax )
639 {
640 for (int kk=nph+1; kk <= nmmax; kk++)
641 {
642 dataout(kk)=dataout(kk%nph);
643 }
644 }
645 }
646 for (i = 0;i <dataout.NElts();i++)
647 {
648 dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
649 }
650 // return dataout;
651}
652
653 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromAlm(SphericalMap<T>& mapq,
654 SphericalMap<T>& mapu,
655 int_4 pixelSizeIndex,
656 const Alm<T>& alme,
657 const Alm<T>& almb) const
658
659synthesis of a polarization map from Alm coefficients. The spheres mapq and mapu contain respectively the Stokes parameters. */
660template<class T>
661void SphericalTransformServer<T>::GenerateFromAlm(SphericalMap<T>& mapq,
662 SphericalMap<T>& mapu,
663 int_4 pixelSizeIndex,
664 const Alm<T>& alme,
665 const Alm<T>& almb) const
666{
667 /*=======================================================================
668 computes a map form its alm for the HEALPIX pixelisation
669 map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
670 = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
671
672 where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
673
674 * the recurrence of Ylm is the standard one (cf Num Rec)
675 * the sum over m is done by FFT
676
677 =======================================================================*/
678 int_4 nlmax=alme.Lmax();
679 if (nlmax != almb.Lmax())
680 {
681 cout << " SphericalTransformServer: les deux tableaux alm n'ont pas la meme taille" << endl;
682 throw SzMismatchError("SphericalTransformServer: les deux tableaux alm n'ont pas la meme taille");
683 }
684 int_4 nmmax=nlmax;
685 int_4 nsmax=0;
686 mapq.Resize(pixelSizeIndex);
687 mapu.Resize(pixelSizeIndex);
688 char* sphere_type=mapq.TypeOfMap();
689 if (strncmp(sphere_type,mapu.TypeOfMap(),4) != 0)
690 {
691 cout << " SphericalTransformServer: les deux spheres ne sont pas de meme type" << endl;
692 cout << " type 1 " << sphere_type << endl;
693 cout << " type 2 " << mapu.TypeOfMap() << endl;
694 throw SzMismatchError("SphericalTransformServer: les deux spheres ne sont pas de meme type");
695
696 }
697 if (strncmp(sphere_type,"RING",4) == 0)
698 {
699 nsmax=mapq.SizeIndex();
700 }
701 else
702 // pour une sphere Gorski le nombre de pixels est 12*nsmax**2
703 // on calcule une quantite equivalente a nsmax pour la sphere-theta-phi
704 // en vue de l'application du critere Healpix : nlmax<=3*nsmax-1
705 // c'est approximatif ; a raffiner.
706 if (strncmp(sphere_type,"TETAFI",6) == 0)
707 {
708 nsmax=(int_4)sqrt(mapq.NbPixels()/12.);
709 }
710 else
711 {
712 cout << " unknown type of sphere : " << sphere_type << endl;
713 throw IOExc(" unknown type of sphere ");
714 }
715 cout << "GenerateFromAlm: the spheres are of type : " << sphere_type << endl;
716 cout << "GenerateFromAlm: size indices (nside) of spheres= " << nsmax << endl;
717 cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
718 if (nlmax>3*nsmax-1)
719 {
720 cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
721 if (strncmp(sphere_type,"TETAFI",6) == 0)
722 {
723 cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
724 }
725 }
726 if (alme.Lmax()!=almb.Lmax())
727 {
728 cout << "GenerateFromAlm: arrays Alme and Almb have not the same size ? " << endl;
729 throw SzMismatchError("SphericalTransformServer: arrays Alme and Almb have not the same size ? ");
730 }
731 mapFromWX(nlmax, nmmax, mapq, mapu, alme, almb);
732 // mapFromPM(nlmax, nmmax, mapq, mapu, alme, almb);
733}
734
735
736 /*! \fn void SOPHYA::SphericalTransformServer::DecomposeToAlm(const SphericalMap<T>& mapq,
737 const SphericalMap<T>& mapu,
738 Alm<T>& alme,
739 Alm<T>& almb,
740 int_4 nlmax,
741 r_8 cos_theta_cut) const
742
743analysis of a polarization map into Alm coefficients.
744
745 The spheres \c mapq and \c mapu contain respectively the Stokes parameters.
746
747 \c a2lme and \c a2lmb will receive respectively electric and magnetic Alm's
748 nlmax : maximum value of the l index
749
750 \c cos_theta_cut : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
751 */
752template<class T>
753void SphericalTransformServer<T>::DecomposeToAlm(SphericalMap<T>& mapq,
754 SphericalMap<T>& mapu,
755 Alm<T>& alme,
756 Alm<T>& almb,
757 int_4 nlmax,
758 r_8 cos_theta_cut,
759 int iterationOrder) const
760{
761 int_4 nmmax = nlmax;
762 carteVersAlm(mapq, mapu, alme, almb, nlmax, cos_theta_cut);
763 if (iterationOrder > 0)
764 {
765 TVector<int_4> fact(iterationOrder+2);
766 fact(0) = 1;
767 int k;
768 for (k=1; k <= iterationOrder+1; k++)
769 {
770 fact(k) = fact(k-1)*k;
771 }
772 Alm<T> alme2(alme);
773 Alm<T> almb2(almb);
774 T Tzero = (T)0.;
775 complex<T> complexZero = complex<T>(Tzero, Tzero);
776 alme = complexZero;
777 almb = complexZero;
778 int signe = 1;
779 int nbIteration = iterationOrder+1;
780 for (k=1; k <= nbIteration; k++)
781 {
782 T facMult = (T)(0.5*signe*fact(iterationOrder)*(2*nbIteration-k)/(fact(k)*fact(nbIteration-k)));
783 for (int m = 0; m <= nmmax; m++)
784 {
785 for (int l = m; l<= nlmax; l++)
786 {
787 alme(l,m) += facMult*alme2(l,m);
788 almb(l,m) += facMult*almb2(l,m);
789 }
790 }
791 if (k == nbIteration) break;
792 signe = -signe;
793 for (int k=0; k< mapq.NbPixels(); k++)
794 {
795 mapq(k) = (T)0.;
796 mapu(k) = (T)0.;
797 }
798 // synthetize a map from the estimated alm
799 GenerateFromAlm(mapq,mapu,mapq.SizeIndex(),alme2,almb2);
800 alme2 = complexZero;
801 almb2 = complexZero;
802 // analyse the new map
803 carteVersAlm(mapq, mapu, alme2, almb2, nlmax, cos_theta_cut);
804 }
805 }
806}
807
808template<class T>
809void SphericalTransformServer<T>::carteVersAlm(const SphericalMap<T>& mapq,
810 const SphericalMap<T>& mapu,
811 Alm<T>& alme,
812 Alm<T>& almb,
813 int_4 nlmax,
814 r_8 cos_theta_cut) const
815{
816 int_4 nmmax = nlmax;
817 // resize et remise a zero
818 alme.ReSizeToLmax(nlmax);
819 almb.ReSizeToLmax(nlmax);
820
821
822 TVector<T> dataq;
823 TVector<T> datau;
824 TVector<int_4> pixNumber;
825
826 char* sphere_type=mapq.TypeOfMap();
827 if (strncmp(sphere_type,mapu.TypeOfMap(),4) != 0)
828 {
829 cout << " SphericalTransformServer: les deux spheres ne sont pas de meme type" << endl;
830 cout << " type 1 " << sphere_type << endl;
831 cout << " type 2 " << mapu.TypeOfMap() << endl;
832 throw SzMismatchError("SphericalTransformServer: les deux spheres ne sont pas de meme type");
833
834 }
835 if (mapq.NbPixels()!=mapu.NbPixels())
836 {
837 cout << " DecomposeToAlm: map Q and map U have not same size ?" << endl;
838 throw SzMismatchError("SphericalTransformServer::DecomposeToAlm: map Q and map U have not same size ");
839 }
840 for (int_4 ith = 0; ith < mapq.NbThetaSlices(); ith++)
841 {
842 r_8 phi0;
843 r_8 theta;
844 mapq.GetThetaSlice(ith,theta,phi0, pixNumber,dataq);
845 mapu.GetThetaSlice(ith,theta,phi0, pixNumber,datau);
846 if (dataq.NElts() != datau.NElts() )
847 {
848 throw SzMismatchError("the spheres have not the same pixelization");
849 }
850 r_8 domega=mapq.PixSolAngle(mapq.PixIndexSph(theta,phi0));
851 double cth = cos(theta);
852 //part of the sky out of the symetric cut
853 bool keep_it = (fabs(cth) >= cos_theta_cut);
854 if (keep_it)
855 {
856 // almFromPM(pixNumber.NElts(), nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb);
857 almFromWX(nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb);
858 }
859 }
860}
861
862
863 /*! \fn void SOPHYA::SphericalTransformServer::almFromWX(int_4 nlmax, int_4 nmmax,
864 r_8 phi0, r_8 domega,
865 r_8 theta,
866 const TVector<T>& dataq,
867 const TVector<T>& datau,
868 Alm<T>& alme,
869 Alm<T>& almb) const
870
871Compute polarized Alm's as :
872\f[
873a_{lm}^E=\frac{1}{\sqrt{2}}\sum_{slices}{\omega_{pix}\left(\,_{w}\lambda_l^m\tilde{Q}-i\,_{x}\lambda_l^m\tilde{U}\right)}
874\f]
875\f[
876a_{lm}^B=\frac{1}{\sqrt{2}}\sum_{slices}{\omega_{pix}\left(i\,_{x}\lambda_l^m\tilde{Q}+\,_{w}\lambda_l^m\tilde{U}\right)}
877\f]
878
879where \f$\tilde{Q}\f$ and \f$\tilde{U}\f$ are C-coefficients computed by FFT (method CFromFourierAnalysis, called by present method) from the Stokes parameters.
880
881\f$\omega_{pix}\f$ are solid angle of each pixel.
882
883dataq, datau : Stokes parameters.
884
885 */
886template<class T>
887void SphericalTransformServer<T>::almFromWX(int_4 nlmax, int_4 nmmax,
888 r_8 phi0, r_8 domega,
889 r_8 theta,
890 const TVector<T>& dataq,
891 const TVector<T>& datau,
892 Alm<T>& alme,
893 Alm<T>& almb) const
894{
895 TVector< complex<T> > phaseq(nmmax+1);
896 TVector< complex<T> > phaseu(nmmax+1);
897 // TVector<complex<T> > datain(nph);
898 for (int i=0;i< nmmax+1;i++)
899 {
900 phaseq(i)=0;
901 phaseu(i)=0;
902 }
903 // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),0.);
904
905 // phaseq = CFromFourierAnalysis(nmmax,dataq,phi0);
906 CFromFourierAnalysis(nmmax,dataq,phaseq, phi0);
907
908 // phaseu= CFromFourierAnalysis(nmmax,datau,phi0);
909 CFromFourierAnalysis(nmmax,datau,phaseu, phi0);
910
911 LambdaWXBuilder lwxb(theta,nlmax,nmmax);
912
913 r_8 sqr2inv=1/Rac2;
914 for (int m = 0; m <= nmmax; m++)
915 {
916 r_8 lambda_w=0.;
917 r_8 lambda_x=0.;
918 lwxb.lam_wx(m, m, lambda_w, lambda_x);
919 complex<T> zi_lam_x((T)0., (T)lambda_x);
920 alme(m,m) += ( (T)(lambda_w)*phaseq(m)-zi_lam_x*phaseu(m) )*(T)(domega*sqr2inv);
921 almb(m,m) += ( (T)(lambda_w)*phaseu(m)+zi_lam_x*phaseq(m) )*(T)(domega*sqr2inv);
922
923 for (int l = m+1; l<= nlmax; l++)
924 {
925 lwxb.lam_wx(l, m, lambda_w, lambda_x);
926 zi_lam_x = complex<T>((T)0., (T)lambda_x);
927 alme(l,m) += ( (T)(lambda_w)*phaseq(m)-zi_lam_x*phaseu(m) )*(T)(domega*sqr2inv);
928 almb(l,m) += ( (T)(lambda_w)*phaseu(m)+zi_lam_x*phaseq(m) )*(T)(domega*sqr2inv);
929 }
930 }
931}
932
933
934 /*! \fn void SOPHYA::SphericalTransformServer::almFromPM(int_4 nph, int_4 nlmax,
935 int_4 nmmax,
936 r_8 phi0, r_8 domega,
937 r_8 theta,
938 const TVector<T>& dataq,
939 const TVector<T>& datau,
940 Alm<T>& alme,
941 Alm<T>& almb) const
942
943Compute polarized Alm's as :
944\f[
945a_{lm}^E=-\frac{1}{2}\sum_{slices}{\omega_{pix}\left(\,_{+}\lambda_l^m\tilde{P^+}+\,_{-}\lambda_l^m\tilde{P^-}\right)}
946\f]
947\f[
948a_{lm}^B=\frac{i}{2}\sum_{slices}{\omega_{pix}\left(\,_{+}\lambda_l^m\tilde{P^+}-\,_{-}\lambda_l^m\tilde{P^-}\right)}
949\f]
950
951where \f$\tilde{P^{\pm}}=\tilde{Q}\pm\tilde{U}\f$ computed by FFT (method CFromFourierAnalysis, called by present method) from the Stokes parameters,\f$Q\f$ and \f$U\f$ .
952
953\f$\omega_{pix}\f$ are solid angle of each pixel.
954
955dataq, datau : Stokes parameters.
956
957 */
958template<class T>
959void SphericalTransformServer<T>::almFromPM(int_4 nph, int_4 nlmax,
960 int_4 nmmax,
961 r_8 phi0, r_8 domega,
962 r_8 theta,
963 const TVector<T>& dataq,
964 const TVector<T>& datau,
965 Alm<T>& alme,
966 Alm<T>& almb) const
967{
968 TVector< complex<T> > phasep(nmmax+1);
969 TVector< complex<T> > phasem(nmmax+1);
970 TVector<complex<T> > datain(nph);
971 for (int i=0;i< nmmax+1;i++)
972 {
973 phasep(i)=0;
974 phasem(i)=0;
975 }
976 int kk;
977 for(kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),datau(kk));
978
979 phasep = CFromFourierAnalysis(nmmax,datain,phi0);
980
981 for(kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),-datau(kk));
982 phasem = CFromFourierAnalysis(nmmax,datain,phi0);
983 LambdaPMBuilder lpmb(theta,nlmax,nmmax);
984
985 for (int m = 0; m <= nmmax; m++)
986 {
987 r_8 lambda_p=0.;
988 r_8 lambda_m=0.;
989 complex<T> im((T)0.,(T)1.);
990 lpmb.lam_pm(m, m, lambda_p, lambda_m);
991
992 alme(m,m) += -( (T)(lambda_p)*phasep(m) + (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
993 almb(m,m) += im*( (T)(lambda_p)*phasep(m) - (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
994 for (int l = m+1; l<= nlmax; l++)
995 {
996 lpmb.lam_pm(l, m, lambda_p, lambda_m);
997 alme(l,m) += -( (T)(lambda_p)*phasep(m) + (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
998 almb(l,m) += im* ( (T)(lambda_p)*phasep(m) - (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
999 }
1000 }
1001}
1002
1003
1004/*! \fn void SOPHYA::SphericalTransformServer::mapFromWX(int_4 nlmax, int_4 nmmax,
1005 SphericalMap<T>& mapq,
1006 SphericalMap<T>& mapu,
1007 const Alm<T>& alme,
1008 const Alm<T>& almb) const
1009
1010synthesis of Stokes parameters following formulae :
1011
1012\f[
1013Q=\sum_{m=-mmax}^{mmax}b_m^qe^{im\varphi}
1014\f]
1015\f[
1016U=\sum_{m=-mmax}^{mmax}b_m^ue^{im\varphi}
1017\f]
1018
1019computed by FFT (method fourierSynthesisFromB called by the present one)
1020
1021with :
1022
1023\f[
1024b_m^q=-\frac{1}{\sqrt{2}}\sum_{l=|m|}^{lmax}{\left(\,_{w}\lambda_l^ma_{lm}^E-i\,_{x}\lambda_l^ma_{lm}^B\right) }
1025\f]
1026\f[
1027b_m^u=\frac{1}{\sqrt{2}}\sum_{l=|m|}^{lmax}{\left(i\,_{x}\lambda_l^ma_{lm}^E+\,_{w}\lambda_l^ma_{lm}^B\right) }
1028\f]
1029 */
1030template<class T>
1031void SphericalTransformServer<T>::mapFromWX(int_4 nlmax, int_4 nmmax,
1032 SphericalMap<T>& mapq,
1033 SphericalMap<T>& mapu,
1034 const Alm<T>& alme,
1035 const Alm<T>& almb) const
1036{
1037 Bm<complex<T> > b_m_theta_q(nmmax);
1038 Bm<complex<T> > b_m_theta_u(nmmax);
1039
1040 for (int_4 ith = 0; ith < mapq.NbThetaSlices();ith++)
1041 {
1042 int_4 nph;
1043 r_8 phi0;
1044 r_8 theta;
1045 TVector<int_4> pixNumber;
1046 TVector<T> datan;
1047
1048 mapq.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
1049 nph = pixNumber.NElts();
1050 // -----------------------------------------------------
1051 // for each theta, and each m, computes
1052 // b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
1053 // ------------------------------------------------------
1054 LambdaWXBuilder lwxb(theta,nlmax,nmmax);
1055 // LambdaPMBuilder lpmb(theta,nlmax,nmmax);
1056 r_8 sqr2inv=1/Rac2;
1057 int m;
1058 for (m = 0; m <= nmmax; m++)
1059 {
1060 r_8 lambda_w=0.;
1061 r_8 lambda_x=0.;
1062 lwxb.lam_wx(m, m, lambda_w, lambda_x);
1063 complex<T> zi_lam_x((T)0., (T)lambda_x);
1064
1065 b_m_theta_q(m) = ( (T)(lambda_w) * alme(m,m) - zi_lam_x * almb(m,m))*(T)sqr2inv ;
1066 b_m_theta_u(m) = ( (T)(lambda_w) * almb(m,m) + zi_lam_x * alme(m,m))*(T)sqr2inv;
1067
1068
1069 for (int l = m+1; l<= nlmax; l++)
1070 {
1071
1072 lwxb.lam_wx(l, m, lambda_w, lambda_x);
1073 zi_lam_x= complex<T>((T)0., (T)lambda_x);
1074
1075 b_m_theta_q(m) += ((T)(lambda_w)*alme(l,m)-zi_lam_x *almb(l,m))*(T)sqr2inv;
1076 b_m_theta_u(m) += ((T)(lambda_w)*almb(l,m)+zi_lam_x *alme(l,m))*(T)sqr2inv;
1077
1078 }
1079 }
1080 // obtains the negative m of b(m,theta) (= complex conjugate)
1081 for (m=1;m<=nmmax;m++)
1082 {
1083 b_m_theta_q(-m) = conj(b_m_theta_q(m));
1084 b_m_theta_u(-m) = conj(b_m_theta_u(m));
1085 }
1086
1087 // TVector<complex<T> > Tempq = fourierSynthesisFromB(b_m_theta_q,nph,phi0);
1088 // TVector<complex<T> > Tempu = fourierSynthesisFromB(b_m_theta_u,nph,phi0);
1089 TVector<T> Tempq = RfourierSynthesisFromB(b_m_theta_q,nph,phi0);
1090 TVector<T> Tempu = RfourierSynthesisFromB(b_m_theta_u,nph,phi0);
1091 for (int i=0;i< nph;i++)
1092 {
1093 // mapq(pixNumber(i))=Tempq(i).real();
1094 // mapu(pixNumber(i))=Tempu(i).real();
1095 mapq(pixNumber(i))=Tempq(i);
1096 mapu(pixNumber(i))=Tempu(i);
1097
1098 }
1099 }
1100}
1101/*! \fn void SOPHYA::SphericalTransformServer::mapFromPM(int_4 nlmax, int_4 nmmax,
1102 SphericalMap<T>& mapq,
1103 SphericalMap<T>& mapu,
1104 const Alm<T>& alme,
1105 const Alm<T>& almb) const
1106
1107synthesis of polarizations following formulae :
1108
1109\f[
1110P^+ = \sum_{m=-mmax}^{mmax} {b_m^+e^{im\varphi} }
1111\f]
1112\f[
1113P^- = \sum_{m=-mmax}^{mmax} {b_m^-e^{im\varphi} }
1114\f]
1115
1116computed by FFT (method fourierSynthesisFromB called by the present one)
1117
1118with :
1119
1120\f[
1121b_m^+=-\sum_{l=|m|}^{lmax}{\,_{+}\lambda_l^m \left( a_{lm}^E+ia_{lm}^B \right) }
1122\f]
1123\f[
1124b_m^-=-\sum_{l=|m|}^{lmax}{\,_{+}\lambda_l^m \left( a_{lm}^E-ia_{lm}^B \right) }
1125\f]
1126 */
1127template<class T>
1128void SphericalTransformServer<T>::mapFromPM(int_4 nlmax, int_4 nmmax,
1129 SphericalMap<T>& mapq,
1130 SphericalMap<T>& mapu,
1131 const Alm<T>& alme,
1132 const Alm<T>& almb) const
1133{
1134 Bm<complex<T> > b_m_theta_p(nmmax);
1135 Bm<complex<T> > b_m_theta_m(nmmax);
1136 for (int_4 ith = 0; ith < mapq.NbThetaSlices();ith++)
1137 {
1138 int_4 nph;
1139 r_8 phi0;
1140 r_8 theta;
1141 TVector<int_4> pixNumber;
1142 TVector<T> datan;
1143
1144 mapq.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
1145 nph = pixNumber.NElts();
1146
1147 // -----------------------------------------------------
1148 // for each theta, and each m, computes
1149 // b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
1150 //------------------------------------------------------
1151
1152 LambdaPMBuilder lpmb(theta,nlmax,nmmax);
1153 int m;
1154 for (m = 0; m <= nmmax; m++)
1155 {
1156 r_8 lambda_p=0.;
1157 r_8 lambda_m=0.;
1158 lpmb.lam_pm(m, m, lambda_p, lambda_m);
1159 complex<T> im((T)0.,(T)1.);
1160
1161 b_m_theta_p(m) = (T)(lambda_p )* (-alme(m,m) - im * almb(m,m));
1162 b_m_theta_m(m) = (T)(lambda_m) * (-alme(m,m) + im * almb(m,m));
1163
1164
1165 for (int l = m+1; l<= nlmax; l++)
1166 {
1167 lpmb.lam_pm(l, m, lambda_p, lambda_m);
1168 b_m_theta_p(m) += (T)(lambda_p)*(-alme(l,m)-im *almb(l,m));
1169 b_m_theta_m(m) += (T)(lambda_m)*(-alme(l,m)+im *almb(l,m));
1170 }
1171 }
1172
1173 // obtains the negative m of b(m,theta) (= complex conjugate)
1174 for (m=1;m<=nmmax;m++)
1175 {
1176 b_m_theta_p(-m) = conj(b_m_theta_m(m));
1177 b_m_theta_m(-m) = conj(b_m_theta_p(m));
1178 }
1179
1180 TVector<complex<T> > Tempp = fourierSynthesisFromB(b_m_theta_p,nph,phi0);
1181 TVector<complex<T> > Tempm = fourierSynthesisFromB(b_m_theta_m,nph,phi0);
1182
1183 for (int i=0;i< nph;i++)
1184 {
1185 mapq(pixNumber(i))=0.5*(Tempp(i)+Tempm(i)).real();
1186 mapu(pixNumber(i))=0.5*(Tempp(i)-Tempm(i)).imag();
1187 }
1188 }
1189}
1190
1191
1192 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromCl(SphericalMap<T>& sphq,
1193 SphericalMap<T>& sphu,
1194 int_4 pixelSizeIndex,
1195 const TVector<T>& Cle,
1196 const TVector<T>& Clb,
1197 const r_8 fwhm) const
1198
1199synthesis of a polarization map from power spectra electric-Cl and magnetic-Cl (Alm's are generated randomly, following a gaussian distribution).
1200 \param fwhm FWHM in arcmin for random generation of Alm's (eg. 5)
1201*/
1202template<class T>
1203void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sphq,
1204 SphericalMap<T>& sphu,
1205 int_4 pixelSizeIndex,
1206 const TVector<T>& Cle,
1207 const TVector<T>& Clb,
1208 const r_8 fwhm) const
1209{
1210 if (Cle.NElts() != Clb.NElts())
1211 {
1212 cout << " SphericalTransformServer: les deux tableaux Cl n'ont pas la meme taille" << endl;
1213 throw SzMismatchError("SphericalTransformServer::GenerateFromCl : two Cl arrays have not same size");
1214 }
1215
1216 // Alm<T> a2lme,a2lmb;
1217 // almFromCl(a2lme, Cle, fwhm);
1218 // almFromCl(a2lmb, Clb, fwhm);
1219 // Alm<T> a2lme = almFromCl(Cle, fwhm);
1220 // Alm<T> a2lmb = almFromCl(Clb, fwhm);
1221 Alm<T> a2lme(Cle, fwhm);
1222 Alm<T> a2lmb(Clb, fwhm);
1223
1224 GenerateFromAlm(sphq,sphu,pixelSizeIndex,a2lme,a2lmb);
1225}
1226 /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromCl(SphericalMap<T>& sph,
1227 int_4 pixelSizeIndex,
1228 const TVector<T>& Cl,
1229 const r_8 fwhm) const
1230
1231synthesis of a temperature map from power spectrum Cl (Alm's are generated randomly, following a gaussian distribution). */
1232template<class T>
1233void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sph,
1234 int_4 pixelSizeIndex,
1235 const TVector<T>& Cl,
1236 const r_8 fwhm) const
1237{
1238
1239 Alm<T> alm(Cl, fwhm);
1240 GenerateFromAlm(sph,pixelSizeIndex, alm );
1241}
1242
1243
1244
1245/*! \fn TVector<T> SOPHYA::SphericalTransformServer::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
1246
1247\return power spectrum from analysis of a temperature map. THE MAP CAN BE MODIFIED (if iterationOrder >0)
1248
1249 \param<nlmax> : maximum value of the l index
1250
1251 \param<cos_theta_cut> : cosinus of the symmetric cut EULER angle theta : cos_theta_cut=0 means no cut ; cos_theta_cut=1 all the sphere is cut.
1252
1253\param<iterationOrder> : 1,2,3,4.... order of an iterative analysis. (Default : 0 -> standard analysis)
1254
1255 */
1256template <class T>
1257TVector<T> SphericalTransformServer<T>::DecomposeToCl(SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut, int iterationOrder) const
1258{
1259 Alm<T> alm;
1260 DecomposeToAlm( sph, alm, nlmax, cos_theta_cut, iterationOrder);
1261 // cout << " SphericalTransformServer : impression alm " << endl;
1262 // alm.Print();
1263 // power spectrum
1264 return alm.powerSpectrum();
1265}
1266
1267#ifdef __CXX_PRAGMA_TEMPLATES__
1268#pragma define_template SphericalTransformServer<r_8>
1269#pragma define_template SphericalTransformServer<r_4>
1270#endif
1271#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
1272template class SphericalTransformServer<r_8>;
1273template class SphericalTransformServer<r_4>;
1274#endif
Note: See TracBrowser for help on using the repository browser.