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

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

ajustement de type double r_8 etc.
File size: 28.5 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_4 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 for (int_4 ith = 0; ith < map.NbThetaSlices(); ith++)
300 {
301 int_4 nph;
302 r_8 phi0;
303 r_8 theta;
304 map.GetThetaSlice(ith,theta,phi0,pixNumber ,data);
305 for (int i=0;i< nmmax+1;i++)
306 {
307 phase(i)=0;
308 }
309 nph = pixNumber.NElts();
310 double cth = cos(theta);
311
312 //part of the sky out of the symetric cut
313 bool keep_it = (abs(cth) >= cos_theta_cut);
314
315 if (keep_it)
316 {
317 // tableau datain a supprimer
318 // TVector<complex<T> > datain(nph);
319 // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(data(kk),(T)0.);
320
321 // phase = CFromFourierAnalysis(nmmax,datain,phi0);
322 phase = CFromFourierAnalysis(nmmax,data,phi0);
323
324 }
325
326 /*-----------------------------------------------------------------------
327 computes the a_lm by integrating over theta
328 lambda_lm(theta) * phas_m(theta)
329 for each m and l
330 -----------------------------------------------------------------------*/
331 // LambdaBuilder lb(theta,nlmax,nmmax);
332 LambdaLMBuilder lb(theta,nlmax,nmmax);
333 r_8 domega=map.PixSolAngle(map.PixIndexSph(theta,phi0));
334 for (int m = 0; m <= nmmax; m++)
335 {
336 alm(m,m) += (T)lb.lamlm(m,m) * phase(m) * (T)domega; //m,m even
337 for (int l = m+1; l<= nlmax; l++)
338 {
339 alm(l,m) += (T)lb.lamlm(l,m) * phase(m)*(T)domega;
340 }
341 }
342 }
343 return alm;
344}
345template<class T>
346TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<complex<T> >datain, r_8 phi0) const
347{
348 /*=======================================================================
349 integrates (data * phi-dependence-of-Ylm) over phi
350 --> function of m can be computed by FFT
351
352 datain est modifie
353 =======================================================================*/
354 int_4 nph=datain.NElts();
355 if (nph <= 0)
356 {
357 throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
358 }
359 TVector<complex<T> > transformedData(nph);
360 fftIntfPtr_-> FFTForward(datain, transformedData);
361
362 //dataout.ReSize(nmmax+1);
363 TVector< complex<T> > dataout(nmmax+1);
364 // dataout.SetTemp(true);
365
366 int im_max=min(nph,nmmax+1);
367 for (int i=0;i< dataout.NElts();i++) dataout(i)=complex<T>((T)0.,(T)0.);
368 for (int i=0;i<im_max;i++) dataout(i)=transformedData(i);
369
370
371 // for (int i = 0;i <im_max;i++){
372 // dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
373 // }
374 for (int kk=nph; kk<dataout.NElts(); kk++) dataout(kk)=dataout(kk%nph);
375 for (int i = 0;i <dataout.NElts();i++){
376 dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
377 }
378 return dataout;
379}
380
381//&&&&&&&&& nouvelle version
382template<class T>
383TVector< complex<T> > SphericalTransformServer<T>::CFromFourierAnalysis(int_4 nmmax, const TVector<T> datain, r_8 phi0) const
384{
385 //=======================================================================
386 // integrates (data * phi-dependence-of-Ylm) over phi
387 // --> function of m can be computed by FFT
388 // ! with 0<= m <= npoints/2 (: Nyquist)
389 // ! because the data is real the negative m are the conjugate of the
390 // ! positive ones
391
392 // datain est modifie
393 //
394 // =======================================================================
395 int_4 nph=datain.NElts();
396 if (nph <= 0)
397 {
398 throw PException("bizarre : vecteur datain de longueur nulle (CFromFourierAnalysis)");
399 }
400 TVector<complex<T> > transformedData;
401 // a remodifier
402 //FFTPackServer ffts;
403 //ffts.setNormalize(false);
404 //ffts.FFTForward(datain, transformedData);
405
406 fftIntfPtr_-> FFTForward(datain, transformedData);
407 //
408
409 //dataout.ReSize(nmmax+1);
410 TVector< complex<T> > dataout(nmmax+1);
411 // dataout.SetTemp(true);
412
413 // on transfere le resultat de la fft dans dataout.
414 // on s'assure que ca ne depasse pas la taille de dataout
415 int sizeOfTransformToGet = min(transformedData.NElts(),nmmax+1);
416 // int im_max=min(transformedData.NElts()-1,nmmax);
417 for (int i=0;i<sizeOfTransformToGet;i++) dataout(i)=transformedData(i);
418
419
420 // si dataout n'est pas plein, on complete jusqu'a nph valeurs (a moins
421 // que dataout ne soit plein avant d'atteindre nph)
422 if (sizeOfTransformToGet == (transformedData.NElts()))
423 {
424 for (int i=transformedData.NElts(); i<min(nph,dataout.NElts()); i++)
425 {
426
427 // dataout(i) = conj(dataout(2*sizeOfTransformToGet-i-2) );
428 dataout(i) = conj(dataout(nph-i) );
429 }
430 // on conplete, si necessaire, par periodicite
431 for (int kk=nph; kk<dataout.NElts(); kk++)
432 {
433 dataout(kk)=dataout(kk%nph);
434 }
435 }
436 for (int i = 0;i <dataout.NElts();i++){
437 dataout(i)*= (complex<T>)(complex<double>(cos(-i*phi0),sin(-i*phi0)));
438 }
439 return dataout;
440}
441
442template<class T>
443void SphericalTransformServer<T>::GenerateFromAlm(SphericalMap<T>& mapq,
444 SphericalMap<T>& mapu,
445 int_4 pixelSizeIndex,
446 const Alm<T>& alme,
447 const Alm<T>& almb) const
448{
449 /*=======================================================================
450 computes a map form its alm for the HEALPIX pixelisation
451 map(theta,phi) = sum_l_m a_lm Y_lm(theta,phi)
452 = sum_m {e^(i*m*phi) sum_l a_lm*lambda_lm(theta)}
453
454 where Y_lm(theta,phi) = lambda(theta) * e^(i*m*phi)
455
456 * the recurrence of Ylm is the standard one (cf Num Rec)
457 * the sum over m is done by FFT
458
459 =======================================================================*/
460 int_4 nlmax=alme.Lmax();
461 if (nlmax != almb.Lmax())
462 {
463 cout << " SphericalTransformServer: les deux tableaux alm n'ont pas la meme taille" << endl;
464 throw SzMismatchError("SphericalTransformServer: les deux tableaux alm n'ont pas la meme taille");
465 }
466 int_4 nmmax=nlmax;
467 int_4 nsmax=0;
468 mapq.Resize(pixelSizeIndex);
469 mapu.Resize(pixelSizeIndex);
470 char* sphere_type=mapq.TypeOfMap();
471 if (strncmp(sphere_type,mapu.TypeOfMap(),4) != 0)
472 {
473 cout << " SphericalTransformServer: les deux spheres ne sont pas de meme type" << endl;
474 cout << " type 1 " << sphere_type << endl;
475 cout << " type 2 " << mapu.TypeOfMap() << endl;
476 throw SzMismatchError("SphericalTransformServer: les deux spheres ne sont pas de meme type");
477
478 }
479 if (strncmp(sphere_type,"RING",4) == 0)
480 {
481 nsmax=mapq.SizeIndex();
482 }
483 else
484 // pour une sphere Gorski le nombre de pixels est 12*nsmax**2
485 // on calcule une quantite equivalente a nsmax pour la sphere-theta-phi
486 // en vue de l'application du critere Healpix : nlmax<=3*nsmax-1
487 // c'est approximatif ; a raffiner.
488 if (strncmp(sphere_type,"TETAFI",6) == 0)
489 {
490 nsmax=(int_4)sqrt(mapq.NbPixels()/12.);
491 }
492 else
493 {
494 cout << " unknown type of sphere : " << sphere_type << endl;
495 throw IOExc(" unknown type of sphere ");
496 }
497 cout << "GenerateFromAlm: the spheres are of type : " << sphere_type << endl;
498 cout << "GenerateFromAlm: size indices (nside) of spheres= " << nsmax << endl;
499 cout << "GenerateFromAlm: nlmax (from Alm) = " << nlmax << endl;
500 if (nlmax>3*nsmax-1)
501 {
502 cout << "GenerateFromAlm: nlmax should be <= 3*nside-1" << endl;
503 if (strncmp(sphere_type,"TETAFI",6) == 0)
504 {
505 cout << " (for this criterium, nsmax is computed as sqrt(nbPixels/12))" << endl;
506 }
507 }
508 if (alme.Lmax()!=almb.Lmax())
509 {
510 cout << "GenerateFromAlm: arrays Alme and Almb have not the same size ? " << endl;
511 throw SzMismatchError("SphericalTransformServer: arrays Alme and Almb have not the same size ? ");
512 }
513 mapFromWX(nlmax, nmmax, mapq, mapu, alme, almb);
514 // mapFromPM(nlmax, nmmax, mapq, mapu, alme, almb);
515}
516
517
518template<class T>
519void SphericalTransformServer<T>::DecomposeToAlm(const SphericalMap<T>& mapq,
520 const SphericalMap<T>& mapu,
521 Alm<T>& alme,
522 Alm<T>& almb,
523 int_4 nlmax,
524 r_8 cos_theta_cut) const
525{
526 int_4 nmmax = nlmax;
527 // resize et remise a zero
528 alme.ReSizeToLmax(nlmax);
529 almb.ReSizeToLmax(nlmax);
530
531
532 TVector<T> dataq;
533 TVector<T> datau;
534 TVector<int_4> pixNumber;
535
536 char* sphere_type=mapq.TypeOfMap();
537 if (strncmp(sphere_type,mapu.TypeOfMap(),4) != 0)
538 {
539 cout << " SphericalTransformServer: les deux spheres ne sont pas de meme type" << endl;
540 cout << " type 1 " << sphere_type << endl;
541 cout << " type 2 " << mapu.TypeOfMap() << endl;
542 throw SzMismatchError("SphericalTransformServer: les deux spheres ne sont pas de meme type");
543
544 }
545 if (mapq.NbPixels()!=mapu.NbPixels())
546 {
547 cout << " DecomposeToAlm: map Q and map U have not same size ?" << endl;
548 throw SzMismatchError("SphericalTransformServer::DecomposeToAlm: map Q and map U have not same size ");
549 }
550 for (int_4 ith = 0; ith < mapq.NbThetaSlices(); ith++)
551 {
552 int_4 nph;
553 r_8 phi0;
554 r_8 theta;
555 mapq.GetThetaSlice(ith,theta,phi0, pixNumber,dataq);
556 mapu.GetThetaSlice(ith,theta,phi0, pixNumber,datau);
557 if (dataq.NElts() != datau.NElts() )
558 {
559 throw SzMismatchError("the spheres have not the same pixelization");
560 }
561 nph = pixNumber.NElts();
562 r_8 domega=mapq.PixSolAngle(mapq.PixIndexSph(theta,phi0));
563 double cth = cos(theta);
564 //part of the sky out of the symetric cut
565 bool keep_it = (abs(cth) >= cos_theta_cut);
566 if (keep_it)
567 {
568 // almFromPM(nph, nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb);
569 almFromWX(nlmax, nmmax, phi0, domega, theta, dataq, datau, alme, almb);
570 }
571 }
572}
573
574
575template<class T>
576void SphericalTransformServer<T>::almFromWX(int_4 nlmax, int_4 nmmax,
577 r_8 phi0, r_8 domega,
578 r_8 theta,
579 const TVector<T>& dataq,
580 const TVector<T>& datau,
581 Alm<T>& alme,
582 Alm<T>& almb) const
583{
584 TVector< complex<T> > phaseq(nmmax+1);
585 TVector< complex<T> > phaseu(nmmax+1);
586 // TVector<complex<T> > datain(nph);
587 for (int i=0;i< nmmax+1;i++)
588 {
589 phaseq(i)=0;
590 phaseu(i)=0;
591 }
592 // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),0.);
593
594 // phaseq = CFromFourierAnalysis(nmmax,datain,phi0);
595 phaseq = CFromFourierAnalysis(nmmax,dataq,phi0);
596
597 // for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(datau(kk),0.);
598
599 // phaseu= CFromFourierAnalysis(nlmax,nmmax,datain,phi0);
600 phaseu= CFromFourierAnalysis(nmmax,datau,phi0);
601
602 LambdaWXBuilder lwxb(theta,nlmax,nmmax);
603
604 r_8 sqr2inv=1/Rac2;
605 for (int m = 0; m <= nmmax; m++)
606 {
607 r_8 lambda_w=0.;
608 r_8 lambda_x=0.;
609 lwxb.lam_wx(m, m, lambda_w, lambda_x);
610 complex<T> zi_lam_x((T)0., (T)lambda_x);
611 alme(m,m) += ( (T)(lambda_w)*phaseq(m)-zi_lam_x*phaseu(m) )*(T)(domega*sqr2inv);
612 almb(m,m) += ( (T)(lambda_w)*phaseu(m)+zi_lam_x*phaseq(m) )*(T)(domega*sqr2inv);
613
614 for (int l = m+1; l<= nlmax; l++)
615 {
616 lwxb.lam_wx(l, m, lambda_w, lambda_x);
617 zi_lam_x = complex<T>((T)0., (T)lambda_x);
618 alme(l,m) += ( (T)(lambda_w)*phaseq(m)-zi_lam_x*phaseu(m) )*(T)(domega*sqr2inv);
619 almb(l,m) += ( (T)(lambda_w)*phaseu(m)+zi_lam_x*phaseq(m) )*(T)(domega*sqr2inv);
620 }
621 }
622}
623
624
625template<class T>
626void SphericalTransformServer<T>::almFromPM(int_4 nph, int_4 nlmax, int_4 nmmax,
627 r_8 phi0, r_8 domega,
628 r_8 theta,
629 const TVector<T>& dataq,
630 const TVector<T>& datau,
631 Alm<T>& alme,
632 Alm<T>& almb) const
633{
634 TVector< complex<T> > phasep(nmmax+1);
635 TVector< complex<T> > phasem(nmmax+1);
636 TVector<complex<T> > datain(nph);
637 for (int i=0;i< nmmax+1;i++)
638 {
639 phasep(i)=0;
640 phasem(i)=0;
641 }
642
643 for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),datau(kk));
644
645 phasep = CFromFourierAnalysis(nmmax,datain,phi0);
646
647 for(int kk=0; kk<nph; kk++) datain(kk)=complex<T>(dataq(kk),-datau(kk));
648 phasem = CFromFourierAnalysis(nmmax,datain,phi0);
649 LambdaPMBuilder lpmb(theta,nlmax,nmmax);
650
651 for (int m = 0; m <= nmmax; m++)
652 {
653 r_8 lambda_p=0.;
654 r_8 lambda_m=0.;
655 complex<T> im((T)0.,(T)1.);
656 lpmb.lam_pm(m, m, lambda_p, lambda_m);
657
658 alme(m,m) += -( (T)(lambda_p)*phasep(m) + (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
659 almb(m,m) += im*( (T)(lambda_p)*phasep(m) - (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
660 for (int l = m+1; l<= nlmax; l++)
661 {
662 lpmb.lam_pm(l, m, lambda_p, lambda_m);
663 alme(l,m) += -( (T)(lambda_p)*phasep(m) + (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
664 almb(l,m) += im* ( (T)(lambda_p)*phasep(m) - (T)(lambda_m)*phasem(m) )*(T)(domega*0.5);
665 }
666 }
667}
668
669
670template<class T>
671void SphericalTransformServer<T>::mapFromWX(int_4 nlmax, int_4 nmmax,
672 SphericalMap<T>& mapq,
673 SphericalMap<T>& mapu,
674 const Alm<T>& alme,
675 const Alm<T>& almb) const
676{
677 Bm<complex<T> > b_m_theta_q(nmmax);
678 Bm<complex<T> > b_m_theta_u(nmmax);
679
680 for (int_4 ith = 0; ith < mapq.NbThetaSlices();ith++)
681 {
682 int_4 nph;
683 r_8 phi0;
684 r_8 theta;
685 TVector<int_4> pixNumber;
686 TVector<T> datan;
687
688 mapq.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
689 nph = pixNumber.NElts();
690 // -----------------------------------------------------
691 // for each theta, and each m, computes
692 // b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
693 // ------------------------------------------------------
694 LambdaWXBuilder lwxb(theta,nlmax,nmmax);
695 // LambdaPMBuilder lpmb(theta,nlmax,nmmax);
696 r_8 sqr2inv=1/Rac2;
697 for (int m = 0; m <= nmmax; m++)
698 {
699 r_8 lambda_w=0.;
700 r_8 lambda_x=0.;
701 lwxb.lam_wx(m, m, lambda_w, lambda_x);
702 complex<T> zi_lam_x((T)0., (T)lambda_x);
703
704 b_m_theta_q(m) = ( (T)(lambda_w) * alme(m,m) - zi_lam_x * almb(m,m))*(T)sqr2inv ;
705 b_m_theta_u(m) = ( (T)(lambda_w) * almb(m,m) + zi_lam_x * alme(m,m))*(T)sqr2inv;
706
707
708 for (int l = m+1; l<= nlmax; l++)
709 {
710
711 lwxb.lam_wx(l, m, lambda_w, lambda_x);
712 zi_lam_x= complex<T>((T)0., (T)lambda_x);
713
714 b_m_theta_q(m) += ((T)(lambda_w)*alme(l,m)-zi_lam_x *almb(l,m))*(T)sqr2inv;
715 b_m_theta_u(m) += ((T)(lambda_w)*almb(l,m)+zi_lam_x *alme(l,m))*(T)sqr2inv;
716
717 }
718 }
719 // obtains the negative m of b(m,theta) (= complex conjugate)
720 for (int m=1;m<=nmmax;m++)
721 {
722 b_m_theta_q(-m) = conj(b_m_theta_q(m));
723 b_m_theta_u(-m) = conj(b_m_theta_u(m));
724 }
725
726 // TVector<complex<T> > Tempq = fourierSynthesisFromB(b_m_theta_q,nph,phi0);
727 // TVector<complex<T> > Tempu = fourierSynthesisFromB(b_m_theta_u,nph,phi0);
728 TVector<T> Tempq = RfourierSynthesisFromB(b_m_theta_q,nph,phi0);
729 TVector<T> Tempu = RfourierSynthesisFromB(b_m_theta_u,nph,phi0);
730 for (int i=0;i< nph;i++)
731 {
732 // mapq(pixNumber(i))=Tempq(i).real();
733 // mapu(pixNumber(i))=Tempu(i).real();
734 mapq(pixNumber(i))=Tempq(i);
735 mapu(pixNumber(i))=Tempu(i);
736
737 }
738 }
739}
740template<class T>
741void SphericalTransformServer<T>::mapFromPM(int_4 nlmax, int_4 nmmax,
742 SphericalMap<T>& mapq,
743 SphericalMap<T>& mapu,
744 const Alm<T>& alme,
745 const Alm<T>& almb) const
746{
747 Bm<complex<T> > b_m_theta_p(nmmax);
748 Bm<complex<T> > b_m_theta_m(nmmax);
749 for (int_4 ith = 0; ith < mapq.NbThetaSlices();ith++)
750 {
751 int_4 nph;
752 r_8 phi0;
753 r_8 theta;
754 TVector<int_4> pixNumber;
755 TVector<T> datan;
756
757 mapq.GetThetaSlice(ith,theta,phi0, pixNumber,datan);
758 nph = pixNumber.NElts();
759
760 // -----------------------------------------------------
761 // for each theta, and each m, computes
762 // b(m,theta) = sum_over_l>m (lambda_l_m(theta) * a_l_m)
763 //------------------------------------------------------
764
765 LambdaPMBuilder lpmb(theta,nlmax,nmmax);
766 for (int m = 0; m <= nmmax; m++)
767 {
768 r_8 lambda_p=0.;
769 r_8 lambda_m=0.;
770 lpmb.lam_pm(m, m, lambda_p, lambda_m);
771 complex<T> im((T)0.,(T)1.);
772
773 b_m_theta_p(m) = (T)(lambda_p )* (-alme(m,m) - im * almb(m,m));
774 b_m_theta_m(m) = (T)(lambda_m) * (-alme(m,m) + im * almb(m,m));
775
776
777 for (int l = m+1; l<= nlmax; l++)
778 {
779 lpmb.lam_pm(l, m, lambda_p, lambda_m);
780 b_m_theta_p(m) += (T)(lambda_p)*(-alme(l,m)-im *almb(l,m));
781 b_m_theta_m(m) += (T)(lambda_m)*(-alme(l,m)+im *almb(l,m));
782 }
783 }
784
785 // obtains the negative m of b(m,theta) (= complex conjugate)
786 for (int m=1;m<=nmmax;m++)
787 {
788 b_m_theta_p(-m) = conj(b_m_theta_m(m));
789 b_m_theta_m(-m) = conj(b_m_theta_p(m));
790 }
791
792 TVector<complex<T> > Tempp = fourierSynthesisFromB(b_m_theta_p,nph,phi0);
793 TVector<complex<T> > Tempm = fourierSynthesisFromB(b_m_theta_m,nph,phi0);
794
795 for (int i=0;i< nph;i++)
796 {
797 mapq(pixNumber(i))=0.5*(Tempp(i)+Tempm(i)).real();
798 mapu(pixNumber(i))=0.5*(Tempp(i)-Tempm(i)).imag();
799 }
800 }
801}
802
803
804template<class T>
805void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sphq,
806 SphericalMap<T>& sphu,
807 int_4 pixelSizeIndex,
808 const TVector<T>& Cle,
809 const TVector<T>& Clb,
810 const r_8 fwhm) const
811{
812 if (Cle.NElts() != Clb.NElts())
813 {
814 cout << " SphericalTransformServer: les deux tableaux Cl n'ont pas la meme taille" << endl;
815 throw SzMismatchError("SphericalTransformServer::GenerateFromCl : two Cl arrays have not same size");
816 }
817
818 // Alm<T> a2lme,a2lmb;
819 // almFromCl(a2lme, Cle, fwhm);
820 // almFromCl(a2lmb, Clb, fwhm);
821 // Alm<T> a2lme = almFromCl(Cle, fwhm);
822 // Alm<T> a2lmb = almFromCl(Clb, fwhm);
823 Alm<T> a2lme(Cle, fwhm);
824 Alm<T> a2lmb(Clb, fwhm);
825
826 GenerateFromAlm(sphq,sphu,pixelSizeIndex,a2lme,a2lmb);
827}
828template<class T>
829void SphericalTransformServer<T>::GenerateFromCl(SphericalMap<T>& sph,
830 int_4 pixelSizeIndex,
831 const TVector<T>& Cl,
832 const r_8 fwhm) const
833{
834
835 Alm<T> alm(Cl, fwhm);
836 GenerateFromAlm(sph,pixelSizeIndex, alm );
837}
838
839
840
841template <class T>
842TVector<T> SphericalTransformServer<T>::DecomposeToCl(const SphericalMap<T>& sph, int_4 nlmax, r_8 cos_theta_cut) const
843{
844 Alm<T> alm=DecomposeToAlm( sph, nlmax, cos_theta_cut);
845 // power spectrum
846 return alm.powerSpectrum();
847}
848
849#ifdef __CXX_PRAGMA_TEMPLATES__
850#pragma define_template SphericalTransformServer<r_8>
851#pragma define_template SphericalTransformServer<r_4>
852#endif
853#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
854template class SphericalTransformServer<r_8>;
855template class SphericalTransformServer<r_4>;
856#endif
Note: See TracBrowser for help on using the repository browser.