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

Last change on this file since 2501 was 2322, checked in by cmv, 23 years ago
  • passage xxstream.h en xxstream
  • compile avec gcc_3.2, gcc_2.96 et cxx En 3.2 le seek from ::end semble marcher (voir Eval/COS/pbseekios.cc)

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