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

Last change on this file since 858 was 833, checked in by ansari, 25 years ago

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