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

Last change on this file since 729 was 729, checked in by ansari, 26 years ago

addition des Ylm etc.
File size: 28.8 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
10
11template<class T>
12void SphericalTransformServer<T>::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const
13{
14 /*=======================================================================
15 computes a map form its alm for the HEALPIX pixelisation
16 map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
17 = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
18
19 where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
20
21 * the recurrence of Ylm is the standard one (cf Num Rec)
22 * the sum over m is done by FFT
23
24 =======================================================================*/
25 int_4 nlmax=alm.Lmax();
26 int_4 nmmax=nlmax;
27 int_4 nsmax=0;
28 map.Resize(pixelSizeIndex);
29 char* sphere_type=map.TypeOfMap();
30 if (strncmp(sphere_type,"RING",4) == 0)
31 {
32 nsmax=map.SizeIndex();
33 }
34 else
35 // pour une sphere Gorski le nombre de pixels est 12*nsmax**2
36 // on calcule une quantite equivalente a nsmax pour la sphere-theta-phi
37 // en vue de l'application du critere Healpix : nlmax<=3*nsmax-1
38 // c'est approximatif ; a raffiner.
39 if (strncmp(sphere_type,"TETAFI",6) == 0)
40 {
41 nsmax=(int_4)sqrt(map.NbPixels()/12.);
42 }
43 else
44 {
45 cout << " unknown type of sphere : " << sphere_type << endl;
46 throw IOExc(" unknown type of sphere: " + (string)sphere_type );
47 }
48 cout << "GenerateFromAlm: the sphere is of type : " << sphere_type << endl;
49 cout << "GenerateFromAlm: size index (nside) of the sphere= " << nsmax << endl;
50 cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
51 if (nlmax>3*nsmax-1)
52 {
53 cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
54 if (strncmp(sphere_type,"TETAFI",6) == 0)
55 {
56 cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
57 }
58 }
59 Bm<complex<T> > b_m_theta(nmmax);
60
61 // map.Resize(nsmax);
62
63
64 // pour chaque tranche en theta
65 for (int ith = 0; ith < map.NbThetaSlices();ith++)
66 {
67 int_4 nph;
68 r_8 phi0;
69 r_8 theta;
70 TVector<int_4> pixNumber;
71 TVector<T> datan;
72
73 map.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
74 nph = pixNumber.NElts();
75
76 // -----------------------------------------------------
77 // for each theta, and each m, computes
78 // b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
79 // ------------------------------------------------------
80 LambdaLMBuilder lb(theta,nlmax,nmmax);
81 // somme sur m de 0 a l'infini
82 for (int m = 0; m <= nmmax; m++)
83 {
84 // somme sur l de m a l'infini
85 b_m_theta(m) = (T)( lb.lamlm(m,m) ) * alm(m,m);
86 // if (ith==0 && m==0)
87 // {
88 // cout << " guy: lmm= " << lb.lamlm(m,m) << " alm " << alm(m,m) << "b00= " << b_m_theta(m) << endl;
89 // }
90 for (int l = m+1; l<= nlmax; l++)
91 {
92 b_m_theta(m) += (T)( lb.lamlm(l,m) ) * alm(l,m);
93
94
95 // if (ith==0 && m==0)
96 // {
97 // cout << " guy:l=" << l << " m= " << m << " lmm= " << lb.lamlm(l,m) << " alm " << alm(l,m) << "b00= " << b_m_theta(m) << endl;
98
99 // }
100
101 }
102 }
103
104 // obtains the negative m of b(m,theta) (= complex conjugate)
105
106 for (int m=1;m<=nmmax;m++)
107 {
108 //compiler doesn't have conj()
109 b_m_theta(-m) = conj(b_m_theta(m));
110 }
111 // ---------------------------------------------------------------
112 // sum_m b(m,theta)*exp(i*m*phi) -> f(phi,theta)
113 // ---------------------------------------------------------------*/
114 // TVector<complex<T> > Temp = fourierSynthesisFromB(b_m_theta,nph,phi0);
115 TVector<T> Temp = RfourierSynthesisFromB(b_m_theta,nph,phi0);
116 for (int i=0;i< nph;i++)
117 {
118 // map(pixNumber(i))=Temp(i).real();
119 map(pixNumber(i))=Temp(i);
120 }
121 }
122}
123
124
125
126template<class T>
127TVector< complex<T> > SphericalTransformServer<T>::fourierSynthesisFromB(const Bm<complex<T> >& b_m, int_4 nph, r_8 phi0) const
128{
129 /*=======================================================================
130 dataout(j) = sum_m datain(m) * exp(i*m*phi(j))
131 with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1
132
133 as the set of frequencies {m} is larger than nph,
134 we wrap frequencies within {0..nph-1}
135 ie m = k*nph + m' with m' in {0..nph-1}
136 then
137 noting bw(m') = exp(i*m'*phi0)
138 * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
139 with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
140 dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
141 = Fourier Transform of bw
142 is real
143
144 NB nph is not necessarily a power of 2
145
146=======================================================================*/
147 //**********************************************************************
148 // pour une valeur de phi (indexee par j) la temperature est la transformee
149 // de Fourier de bm (somme sur m de -nmax a +nmmax de bm*exp(i*m*phi)).
150 // on demande nph (nombre de pixels sur la tranche) valeurs de transformees, pour nph valeurs de phi, regulierement reparties sur 2*pi. On a:
151 // DT/T(j) = sum_m b(m) * exp(i*m*phi(j))
152 // sommation de -infini a +infini, en fait limitee a -nmamx, +nmmax
153 // On pose m=k*nph + m', avec m' compris entre 0 et nph-1. Alors :
154 // DT/T(j) = somme_k somme_m' b(k*nph + m')*exp(i*(k*nph + m')*phi(j))
155 // somme_k : de -infini a +infini
156 // somme_m' : de 0 a nph-1
157 // On echange les sommations :
158 // DT/T(j) = somme_k (exp(i*m'*phi(j)) somme_m' b(k*nph + m')*exp(i*(k*nph*phi(j))
159 // mais phi(j) est un multiple entier de 2*pi/nph, la seconde exponentielle
160 // vaut 1.
161 // Il reste a calculer les transformees de Fourier de somme_m' b(k*nph + m')
162 // si phi0 n'est pas nul, il y a juste un decalage a faire.
163 //**********************************************************************
164
165 TVector< complex<T> > bw(nph);
166 TVector< complex<T> > dataout(nph);
167 TVector< complex<T> > data(nph);
168
169
170 for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
171 for (int m=-b_m.Mmax();m<=-1;m++)
172 {
173 int maux=m;
174 while (maux<0) maux+=nph;
175 int iw=maux%nph;
176 double aux=(m-iw)*phi0;
177 bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) ) ;
178 }
179 for (int m=0;m<=b_m.Mmax();m++)
180 {
181 // int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
182 int iw=m%nph;
183 double aux=(m-iw)*phi0;
184 bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
185 }
186
187 // applies the shift in position <-> phase factor in Fourier space
188 for (int mprime=0; mprime < nph; mprime++)
189 {
190 complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
191 data(mprime)=bw(mprime)*
192 (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
193 }
194
195 //sortie.ReSize(nph);
196 TVector< complex<T> > sortie(nph);
197 // sortie.SetTemp(true);
198
199 fftIntfPtr_-> FFTBackward(data, sortie);
200
201 return sortie;
202}
203
204//********************************************
205template<class T>
206TVector<T> SphericalTransformServer<T>::RfourierSynthesisFromB(const Bm<complex<T> >& b_m, int_4 nph, r_8 phi0) const
207{
208 /*=======================================================================
209 dataout(j) = sum_m datain(m) * exp(i*m*phi(j))
210 with phi(j) = j*2pi/nph + kphi0*pi/nph and kphi0 =0 or 1
211
212 as the set of frequencies {m} is larger than nph,
213 we wrap frequencies within {0..nph-1}
214 ie m = k*nph + m' with m' in {0..nph-1}
215 then
216 noting bw(m') = exp(i*m'*phi0)
217 * sum_k (datain(k*nph+m') exp(i*k*pi*kphi0))
218 with bw(nph-m') = CONJ(bw(m')) (if datain(-m) = CONJ(datain(m)))
219 dataout(j) = sum_m' [ bw(m') exp (i*j*m'*2pi/nph) ]
220 = Fourier Transform of bw
221 is real
222
223 NB nph is not necessarily a power of 2
224
225=======================================================================*/
226 //**********************************************************************
227 // pour une valeur de phi (indexee par j) la temperature est la transformee
228 // de Fourier de bm (somme sur m de -nmax a +nmmax de bm*exp(i*m*phi)).
229 // on demande nph (nombre de pixels sur la tranche) valeurs de transformees, pour nph valeurs de phi, regulierement reparties sur 2*pi. On a:
230 // DT/T(j) = sum_m b(m) * exp(i*m*phi(j))
231 // sommation de -infini a +infini, en fait limitee a -nmamx, +nmmax
232 // On pose m=k*nph + m', avec m' compris entre 0 et nph-1. Alors :
233 // DT/T(j) = somme_k somme_m' b(k*nph + m')*exp(i*(k*nph + m')*phi(j))
234 // somme_k : de -infini a +infini
235 // somme_m' : de 0 a nph-1
236 // On echange les sommations :
237 // DT/T(j) = somme_k (exp(i*m'*phi(j)) somme_m' b(k*nph + m')*exp(i*(k*nph*phi(j))
238 // mais phi(j) est un multiple entier de 2*pi/nph, la seconde exponentielle
239 // vaut 1.
240 // Il reste a calculer les transformees de Fourier de somme_m' b(k*nph + m')
241 // si phi0 n'est pas nul, il y a juste un decalage a faire.
242 //**********************************************************************
243
244 TVector< complex<T> > bw(nph);
245 TVector< complex<T> > dataout(nph);
246 TVector< complex<T> > data(nph/2+1);
247
248
249 for (int kk=0; kk<bw.NElts(); kk++) bw(kk)=(T)0.;
250 for (int m=-b_m.Mmax();m<=-1;m++)
251 {
252 int maux=m;
253 while (maux<0) maux+=nph;
254 int iw=maux%nph;
255 double aux=(m-iw)*phi0;
256 bw(iw) += b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) ) ;
257 }
258 for (int m=0;m<=b_m.Mmax();m++)
259 {
260 // int iw=((m % nph) +nph) % nph; //between 0 and nph = m'
261 int iw=m%nph;
262 double aux=(m-iw)*phi0;
263 bw(iw)+=b_m(m) * complex<T>( (T)cos(aux),(T)sin(aux) );
264 }
265
266 // applies the shift in position <-> phase factor in Fourier space
267 for (int mprime=0; mprime <= nph/2; mprime++)
268 {
269 complex<double> aux(cos(mprime*phi0),sin(mprime*phi0));
270 data(mprime)=bw(mprime)*
271 (complex<T>)(complex<double>(cos(mprime*phi0),sin(mprime*phi0)));
272 }
273
274 //sortie.ReSize(nph);
275 TVector<T> sortie;
276 // sortie.SetTemp(true);
277
278 fftIntfPtr_-> FFTBackward(data, sortie);
279
280 return sortie;
281}
282//*******************************************
283
284template<class T>
285 Alm<T> SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& map, int_4 nlmax, r_8 cos_theta_cut) const
286{
287
288 /*-----------------------------------------------------------------------
289 computes the integral in phi : phas_m(theta)
290 for each parallele from north to south pole
291 -----------------------------------------------------------------------*/
292 TVector<T> data;
293 TVector<int_4> pixNumber;
294 int_4 nmmax = nlmax;
295 TVector< complex<T> > phase(nmmax+1);
296 Alm<T> alm;
297 // alm.SetTemp(true);
298 alm.ReSizeToLmax(nlmax);
299 // cout << " nombre de tranches " << map.NbThetaSlices() << endl;
300 for (int ith = 0; ith < map.NbThetaSlices(); ith++)
301 {
302 int_4 nph;
303 r_8 phi0;
304 r_8 theta;
305 map.GetThetaSlice(ith,theta,phi0,pixNumber ,data);
306 for (int i=0;i< nmmax+1;i++)
307 {
308 phase(i)=0;
309 }
310 nph = pixNumber.NElts();
311 double cth = cos(theta);
312
313 //part of the sky out of the symetric cut
314 bool keep_it = (abs(cth) >= cos_theta_cut);
315
316 if (keep_it)
317 {
318 // tableau datain a supprimer
319 // TVector<complex<T> > datain(nph);
320 // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(data(kk),(T)0.);
321
322 // phase = CFromFourierAnalysis(nlmax,nmmax,datain,phi0);
323 phase = CFromFourierAnalysis(nlmax,nmmax,data,phi0);
324
325 }
326
327 /*-----------------------------------------------------------------------
328 computes the a_lm by integrating over theta
329 lambda_lm(theta) * phas_m(theta)
330 for each m and l
331 -----------------------------------------------------------------------*/
332 // LambdaBuilder lb(theta,nlmax,nmmax);
333 LambdaLMBuilder lb(theta,nlmax,nmmax);
334 r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
335 for (int m = 0; m <= nmmax; m++)
336 {
337 alm(m,m) += (T)lb.lamlm(m,m) * phase(m) * (T)domega; //m,m even
338 for (int l = m+1; l<= nlmax; l++)
339 {
340 alm(l,m) += (T)lb.lamlm(l,m) * phase(m)*(T)domega;
341 }
342 }
343 }
344 return alm;
345}
346template<class T>
347TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nlmax,int_4 nmmax, const TVector<complex<T> >datain, r_8 phi0) const
348{
349 /*=======================================================================
350 integrates (data * phi-dependence-of-Ylm) over phi
351 --> function of m can be computed by FFT
352
353 datain est modifie
354 =======================================================================*/
355 int_4 nph=datain.NElts();
356 if (nph <= 0)
357 {
358 throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
359 }
360 TVector<complex<T> > transformedData(nph);
361 fftIntfPtr_-> FFTForward(datain, transformedData);
362
363 //dataout.ReSize(nmmax+1);
364 TVector< complex<T> > dataout(nmmax+1);
365 // dataout.SetTemp(true);
366
367 int im_max=min(nph,nmmax+1);
368 for (int i=0;i< dataout.NElts();i++) dataout(i)=complex<T>((T)0.,(T)0.);
369 for (int i=0;i<im_max;i++) dataout(i)=transformedData(i);
370
371
372 // for (int i = 0;i <im_max;i++){
373 // dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
374 // }
375 for (int kk=nph; kk<dataout.NElts(); kk++) dataout(kk)=dataout(kk%nph);
376 for (int i = 0;i <dataout.NElts();i++){
377 dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
378 }
379 return dataout;
380}
381
382//&&&&&&&&& nouvelle version
383template<class T>
384TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nlmax,int_4 nmmax, const TVector<T> datain, r_8 phi0) const
385{
386 //=======================================================================
387 // integrates (data * phi-dependence-of-Ylm) over phi
388 // --> function of m can be computed by FFT
389 // ! with 0<= m <= npoints/2 (: Nyquist)
390 // ! because the data is real the negative m are the conjugate of the
391 // ! positive ones
392
393 // datain est modifie
394 //
395 // =======================================================================
396 int_4 nph=datain.NElts();
397 if (nph <= 0)
398 {
399 throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
400 }
401 TVector<complex<T> > transformedData;
402 // a remodifier
403 //FFTPackServer ffts;
404 //ffts.setNormalize(false);
405 //ffts.FFTForward(datain, transformedData);
406
407 fftIntfPtr_-> FFTForward(datain, transformedData);
408 //
409
410 //dataout.ReSize(nmmax+1);
411 TVector< complex<T> > dataout(nmmax+1);
412 // dataout.SetTemp(true);
413
414 // on transfere le resultat de la fft dans dataout.
415 // on s'assure que ca ne depasse pas la taille de dataout
416 int sizeOfTransformToGet = min(transformedData.NElts(),nmmax+1);
417 // int im_max=min(transformedData.NElts()-1,nmmax);
418 for (int i=0;i<sizeOfTransformToGet;i++) dataout(i)=transformedData(i);
419
420
421 // si dataout n'est pas plein, on complete jusqu'a nph valeurs (a moins
422 // que dataout ne soit plein avant d'atteindre nph)
423 if (sizeOfTransformToGet == (transformedData.NElts()))
424 {
425 for (int i=transformedData.NElts(); i<min(nph,dataout.NElts()); i++)
426 {
427
428 // dataout(i) = conj(dataout(2*sizeOfTransformToGet-i-2) );
429 dataout(i) = conj(dataout(nph-i) );
430 }
431 // on conplete, si necessaire, par periodicite
432 for (int kk=nph; kk<dataout.NElts(); kk++)
433 {
434 dataout(kk)=dataout(kk%nph);
435 }
436 }
437 for (int i = 0;i <dataout.NElts();i++){
438 dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
439 }
440 return dataout;
441}
442
443template<class T>
444void SphericalTransformServer<T>::GenerateFromAlm(SphericalMap<T>& mapq,
445 SphericalMap<T>& mapu,
446 int_4 pixelSizeIndex,
447 const Alm<T>& alme,
448 const Alm<T>& almb) const
449{
450 /*=======================================================================
451 computes a map form its alm for the HEALPIX pixelisation
452 map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
453 = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
454
455 where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
456
457 * the recurrence of Ylm is the standard one (cf Num Rec)
458 * the sum over m is done by FFT
459
460 =======================================================================*/
461 int_4 nlmax=alme.Lmax();
462 if (nlmax != almb.Lmax())
463 {
464 cout << " SphericalTransformServer: les deux tableaux alm n'ont pas la meme taille" << endl;
465 throw SzMismatchError("SphericalTransformServer: les deux tableaux alm n'ont pas la meme taille");
466 }
467 int_4 nmmax=nlmax;
468 int_4 nsmax=0;
469 mapq.Resize(pixelSizeIndex);
470 mapu.Resize(pixelSizeIndex);
471 char* sphere_type=mapq.TypeOfMap();
472 if (strncmp(sphere_type,mapu.TypeOfMap(),4) != 0)
473 {
474 cout << " SphericalTransformServer: les deux spheres ne sont pas de meme type" << endl;
475 cout << " type 1 " << sphere_type << endl;
476 cout << " type 2 " << mapu.TypeOfMap() << endl;
477 throw SzMismatchError("SphericalTransformServer: les deux spheres ne sont pas de meme type");
478
479 }
480 if (strncmp(sphere_type,"RING",4) == 0)
481 {
482 nsmax=mapq.SizeIndex();
483 }
484 else
485 // pour une sphere Gorski le nombre de pixels est 12*nsmax**2
486 // on calcule une quantite equivalente a nsmax pour la sphere-theta-phi
487 // en vue de l'application du critere Healpix : nlmax<=3*nsmax-1
488 // c'est approximatif ; a raffiner.
489 if (strncmp(sphere_type,"TETAFI",6) == 0)
490 {
491 nsmax=(int_4)sqrt(mapq.NbPixels()/12.);
492 }
493 else
494 {
495 cout << " unknown type of sphere : " << sphere_type << endl;
496 throw IOExc(" unknown type of sphere ");
497 }
498 cout << "GenerateFromAlm: the spheres are of type : " << sphere_type << endl;
499 cout << "GenerateFromAlm: size indices (nside) of spheres= " << nsmax << endl;
500 cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
501 if (nlmax>3*nsmax-1)
502 {
503 cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
504 if (strncmp(sphere_type,"TETAFI",6) == 0)
505 {
506 cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
507 }
508 }
509 if (alme.Lmax()!=almb.Lmax())
510 {
511 cout << "GenerateFromAlm: arrays Alme and Almb have not the same size ? " << endl;
512 throw SzMismatchError("SphericalTransformServer: arrays Alme and Almb have not the same size ? ");
513 }
514 mapFromWX(nlmax, nmmax, mapq, mapu, alme, almb);
515 // mapFromPM(nlmax, nmmax, mapq, mapu, alme, almb);
516}
517
518
519template<class T>
520void SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& mapq,
521 const SphericalMap<T>& mapu,
522 Alm<T>& alme,
523 Alm<T>& almb,
524 int_4 nlmax,
525 r_8 cos_theta_cut) const
526{
527 int_4 nmmax = nlmax;
528 // resize et remise a zero
529 alme.ReSizeToLmax(nlmax);
530 almb.ReSizeToLmax(nlmax);
531
532
533 TVector<T> dataq;
534 TVector<T> datau;
535 TVector<int_4> pixNumber;
536
537 char* sphere_type=mapq.TypeOfMap();
538 if (strncmp(sphere_type,mapu.TypeOfMap(),4) != 0)
539 {
540 cout << " SphericalTransformServer: les deux spheres ne sont pas de meme type" << endl;
541 cout << " type 1 " << sphere_type << endl;
542 cout << " type 2 " << mapu.TypeOfMap() << endl;
543 throw SzMismatchError("SphericalTransformServer: les deux spheres ne sont pas de meme type");
544
545 }
546 if (mapq.NbPixels()!=mapu.NbPixels())
547 {
548 cout << " DecomposeToAlm: map Q and map U have not same size ?" << endl;
549 throw SzMismatchError("SphericalTransformServer::DecomposeToAlm: map Q and map U have not same size ");
550 }
551 for (int ith = 0; ith < mapq.NbThetaSlices(); ith++)
552 {
553 int_4 nph;
554 r_8 phi0;
555 r_8 theta;
556 mapq.GetThetaSlice(ith,theta,phi0, pixNumber,dataq);
557 mapu.GetThetaSlice(ith,theta,phi0, pixNumber,datau);
558 if (dataq.NElts() != datau.NElts() )
559 {
560 throw SzMismatchError("the spheres have not the same pixelization");
561 }
562 nph = pixNumber.NElts();
563 r_8 domega=mapq.PixSolAngle(mapq.PixIndexSph(theta,phi0));
564 double cth = cos(theta);
565 //part of the sky out of the symetric cut
566 bool keep_it = (abs(cth) >= cos_theta_cut);
567 if (keep_it)
568 {
569 // almFromPM(nph, nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb);
570 almFromWX(nph, nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb);
571 }
572 }
573}
574
575
576template<class T>
577void SphericalTransformServer<T>::almFromWX(int_4 nph, int_4 nlmax, int_4 nmmax,
578 r_8 phi0, r_8 domega,
579 r_8 theta,
580 const TVector<T>& dataq,
581 const TVector<T>& datau,
582 Alm<T>& alme,
583 Alm<T>& almb) const
584{
585 TVector< complex<T> > phaseq(nmmax+1);
586 TVector< complex<T> > phaseu(nmmax+1);
587 // TVector<complex<T> > datain(nph);
588 for (int i=0;i< nmmax+1;i++)
589 {
590 phaseq(i)=0;
591 phaseu(i)=0;
592 }
593 // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),0.);
594
595 // phaseq = CFromFourierAnalysis(nlmax,nmmax,datain,phi0);
596 phaseq = CFromFourierAnalysis(nlmax,nmmax,dataq,phi0);
597
598 // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(datau(kk),0.);
599
600 // phaseu= CFromFourierAnalysis(nlmax,nmmax,datain,phi0);
601 phaseu= CFromFourierAnalysis(nlmax,nmmax,datau,phi0);
602
603 LambdaWXBuilder lwxb(theta,nlmax,nmmax);
604
605 r_8 sqr2inv=1/Rac2;
606 for (int m = 0; m <= nmmax; m++)
607 {
608 r_8 lambda_w=0.;
609 r_8 lambda_x=0.;
610 lwxb.lam_wx(m, m, lambda_w, lambda_x);
611 complex<T> zi_lam_x((T)0., (T)lambda_x);
612 alme(m,m) += ( (T)(lambda_w)*phaseq(m)-zi_lam_x*phaseu(m) )*(T)(domega*sqr2inv);
613 almb(m,m) += ( (T)(lambda_w)*phaseu(m)+zi_lam_x*phaseq(m) )*(T)(domega*sqr2inv);
614
615 for (int l = m+1; l<= nlmax; l++)
616 {
617 lwxb.lam_wx(l, m, lambda_w, lambda_x);
618 zi_lam_x = complex<T>((T)0., (T)lambda_x);
619 alme(l,m) += ( (T)(lambda_w)*phaseq(m)-zi_lam_x*phaseu(m) )*(T)(domega*sqr2inv);
620 almb(l,m) += ( (T)(lambda_w)*phaseu(m)+zi_lam_x*phaseq(m) )*(T)(domega*sqr2inv);
621 }
622 }
623}
624
625
626template<class T>
627void SphericalTransformServer<T>::almFromPM(int_4 nph, int_4 nlmax, int_4 nmmax,
628 r_8 phi0, r_8 domega,
629 r_8 theta,
630 const TVector<T>& dataq,
631 const TVector<T>& datau,
632 Alm<T>& alme,
633 Alm<T>& almb) const
634{
635 TVector< complex<T> > phasep(nmmax+1);
636 TVector< complex<T> > phasem(nmmax+1);
637 TVector<complex<T> > datain(nph);
638 for (int i=0;i< nmmax+1;i++)
639 {
640 phasep(i)=0;
641 phasem(i)=0;
642 }
643
644 for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),datau(kk));
645
646 phasep = CFromFourierAnalysis(nlmax,nmmax,datain,phi0);
647
648 for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),-datau(kk));
649 phasem = CFromFourierAnalysis(nlmax,nmmax,datain,phi0);
650 LambdaPMBuilder lpmb(theta,nlmax,nmmax);
651
652 for (int m = 0; m <= nmmax; m++)
653 {
654 r_8 lambda_p=0.;
655 r_8 lambda_m=0.;
656 complex<T> im((T)0.,(T)1.);
657 lpmb.lam_pm(m, m, lambda_p, lambda_m);
658
659 alme(m,m) += -( (T)(lambda_p)*phasep(m) + (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
660 almb(m,m) += im*( (T)(lambda_p)*phasep(m) - (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
661 for (int l = m+1; l<= nlmax; l++)
662 {
663 lpmb.lam_pm(l, m, lambda_p, lambda_m);
664 alme(l,m) += -( (T)(lambda_p)*phasep(m) + (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
665 almb(l,m) += im* ( (T)(lambda_p)*phasep(m) - (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
666 }
667 }
668}
669
670
671template<class T>
672void SphericalTransformServer<T>::mapFromWX(int_4 nlmax, int_4 nmmax,
673 SphericalMap<T>& mapq,
674 SphericalMap<T>& mapu,
675 const Alm<T>& alme,
676 const Alm<T>& almb) const
677{
678 Bm<complex<T> > b_m_theta_q(nmmax);
679 Bm<complex<T> > b_m_theta_u(nmmax);
680
681 for (int ith = 0; ith < mapq.NbThetaSlices();ith++)
682 {
683 int_4 nph;
684 r_8 phi0;
685 r_8 theta;
686 TVector<int_4> pixNumber;
687 TVector<T> datan;
688
689 mapq.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
690 nph = pixNumber.NElts();
691 // -----------------------------------------------------
692 // for each theta, and each m, computes
693 // b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
694 // ------------------------------------------------------
695 LambdaWXBuilder lwxb(theta,nlmax,nmmax);
696 // LambdaPMBuilder lpmb(theta,nlmax,nmmax);
697 r_8 sqr2inv=1/Rac2;
698 for (int m = 0; m <= nmmax; m++)
699 {
700 r_8 lambda_w=0.;
701 r_8 lambda_x=0.;
702 lwxb.lam_wx(m, m, lambda_w, lambda_x);
703 complex<T> zi_lam_x((T)0., (T)lambda_x);
704
705 b_m_theta_q(m) = ( (T)(lambda_w) * alme(m,m) - zi_lam_x * almb(m,m))*(T)sqr2inv ;
706 b_m_theta_u(m) = ( (T)(lambda_w) * almb(m,m) + zi_lam_x * alme(m,m))*(T)sqr2inv;
707
708
709 for (int l = m+1; l<= nlmax; l++)
710 {
711
712 lwxb.lam_wx(l, m, lambda_w, lambda_x);
713 zi_lam_x= complex<T>((T)0., (T)lambda_x);
714
715 b_m_theta_q(m) += ((T)(lambda_w)*alme(l,m)-zi_lam_x *almb(l,m))*(T)sqr2inv;
716 b_m_theta_u(m) += ((T)(lambda_w)*almb(l,m)+zi_lam_x *alme(l,m))*(T)sqr2inv;
717
718 }
719 }
720 // obtains the negative m of b(m,theta) (= complex conjugate)
721 for (int m=1;m<=nmmax;m++)
722 {
723 b_m_theta_q(-m) = conj(b_m_theta_q(m));
724 b_m_theta_u(-m) = conj(b_m_theta_u(m));
725 }
726
727 // TVector<complex<T> > Tempq = fourierSynthesisFromB(b_m_theta_q,nph,phi0);
728 // TVector<complex<T> > Tempu = fourierSynthesisFromB(b_m_theta_u,nph,phi0);
729 TVector<T> Tempq = RfourierSynthesisFromB(b_m_theta_q,nph,phi0);
730 TVector<T> Tempu = RfourierSynthesisFromB(b_m_theta_u,nph,phi0);
731 for (int i=0;i< nph;i++)
732 {
733 // mapq(pixNumber(i))=Tempq(i).real();
734 // mapu(pixNumber(i))=Tempu(i).real();
735 mapq(pixNumber(i))=Tempq(i);
736 mapu(pixNumber(i))=Tempu(i);
737
738 }
739 }
740}
741template<class T>
742void SphericalTransformServer<T>::mapFromPM(int_4 nlmax, int_4 nmmax,
743 SphericalMap<T>& mapq,
744 SphericalMap<T>& mapu,
745 const Alm<T>& alme,
746 const Alm<T>& almb) const
747{
748 Bm<complex<T> > b_m_theta_p(nmmax);
749 Bm<complex<T> > b_m_theta_m(nmmax);
750 for (int ith = 0; ith < mapq.NbThetaSlices();ith++)
751 {
752 int_4 nph;
753 r_8 phi0;
754 r_8 theta;
755 TVector<int_4> pixNumber;
756 TVector<T> datan;
757
758 mapq.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
759 nph = pixNumber.NElts();
760
761 // -----------------------------------------------------
762 // for each theta, and each m, computes
763 // b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
764 //------------------------------------------------------
765
766 LambdaPMBuilder lpmb(theta,nlmax,nmmax);
767 for (int m = 0; m <= nmmax; m++)
768 {
769 r_8 lambda_p=0.;
770 r_8 lambda_m=0.;
771 lpmb.lam_pm(m, m, lambda_p, lambda_m);
772 complex<T> im((T)0.,(T)1.);
773
774 b_m_theta_p(m) = (T)(lambda_p )* (-alme(m,m) - im * almb(m,m));
775 b_m_theta_m(m) = (T)(lambda_m) * (-alme(m,m) + im * almb(m,m));
776
777
778 for (int l = m+1; l<= nlmax; l++)
779 {
780 lpmb.lam_pm(l, m, lambda_p, lambda_m);
781 b_m_theta_p(m) += (T)(lambda_p)*(-alme(l,m)-im *almb(l,m));
782 b_m_theta_m(m) += (T)(lambda_m)*(-alme(l,m)+im *almb(l,m));
783 }
784 }
785
786 // obtains the negative m of b(m,theta) (= complex conjugate)
787 for (int m=1;m<=nmmax;m++)
788 {
789 b_m_theta_p(-m) = conj(b_m_theta_m(m));
790 b_m_theta_m(-m) = conj(b_m_theta_p(m));
791 }
792
793 TVector<complex<T> > Tempp = fourierSynthesisFromB(b_m_theta_p,nph,phi0);
794 TVector<complex<T> > Tempm = fourierSynthesisFromB(b_m_theta_m,nph,phi0);
795
796 for (int i=0;i< nph;i++)
797 {
798 mapq(pixNumber(i))=0.5*(Tempp(i)+Tempm(i)).real();
799 mapu(pixNumber(i))=0.5*(Tempp(i)-Tempm(i)).imag();
800 }
801 }
802}
803
804
805template<class T>
806void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sphq,
807 SphericalMap<T>& sphu,
808 int_4 pixelSizeIndex,
809 const TVector<T>& Cle,
810 const TVector<T>& Clb,
811 const r_8 fwhm) const
812{
813 if (Cle.NElts() != Clb.NElts())
814 {
815 cout << " SphericalTransformServer: les deux tableaux Cl n'ont pas la meme taille" << endl;
816 throw SzMismatchError("SphericalTransformServer::GenerateFromCl : two Cl arrays have not same size");
817 }
818
819 // Alm<T> a2lme,a2lmb;
820 // almFromCl(a2lme, Cle, fwhm);
821 // almFromCl(a2lmb, Clb, fwhm);
822 // Alm<T> a2lme = almFromCl(Cle, fwhm);
823 // Alm<T> a2lmb = almFromCl(Clb, fwhm);
824 Alm<T> a2lme(Cle, fwhm);
825 Alm<T> a2lmb(Clb, fwhm);
826
827 GenerateFromAlm(sphq,sphu,pixelSizeIndex,a2lme,a2lmb);
828}
829template<class T>
830void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sph,
831 int_4 pixelSizeIndex,
832 const TVector<T>& Cl,
833 const r_8 fwhm) const
834{
835
836 // Alm<T> alm;
837 // almFromCl(alm, Cl, fwhm);
838 //Alm<T> alm = almFromCl(Cl, fwhm);
839 Alm<T> alm(Cl, fwhm);
840 GenerateFromAlm(sph,pixelSizeIndex, alm );
841}
842
843
844
845template <class T>
846TVector<T> SphericalTransformServer<T>::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
847{
848 // Alm<T> alm;
849 // DecomposeToAlm( sph, alm,nlmax, cos_theta_cut);
850 Alm<T> alm=DecomposeToAlm( sph, nlmax, cos_theta_cut);
851 // power spectrum
852 return alm.powerSpectrum();
853}
854
855#ifdef __CXX_PRAGMA_TEMPLATES__
856#pragma define_template SphericalTransformServer<r_8>
857#pragma define_template SphericalTransformServer<r_4>
858#endif
859#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
860template class SphericalTransformServer<r_8>;
861template class SphericalTransformServer<r_4>;
862#endif
Note: See TracBrowser for help on using the repository browser.