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

Last change on this file since 1770 was 1756, checked in by lemeur, 24 years ago

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