1 | #include "spherethetaphi.h"
|
---|
2 | #include "smathconst.h"
|
---|
3 | #include <complex>
|
---|
4 | #include "piocmplx.h"
|
---|
5 | #include "fiondblock.h"
|
---|
6 | #include <iostream.h>
|
---|
7 |
|
---|
8 |
|
---|
9 | /*!
|
---|
10 | \class SOPHYA::SphereThetaPhi
|
---|
11 | \ingroup SkyMap
|
---|
12 | Class implementing spherical maps, with equal latitude (iso-theta) rings
|
---|
13 | pixelisation scheme - with template data types (double, float, complex, ...)
|
---|
14 |
|
---|
15 | sphere splitted with respect to theta, phi : each hemisphere is
|
---|
16 | splitted into (m-1) parallels (equator does not enter into account).
|
---|
17 | This operation defines m slices, each of which is splitted into
|
---|
18 | equidistant meridians. This splitting is realized in such a way that
|
---|
19 | all pixels have the same area and are as square as possible.
|
---|
20 |
|
---|
21 | One begins with the hemisphere with positive z, starting from the pole
|
---|
22 | toward the equator. The first pixel is the polar cap ; it is circular
|
---|
23 | and centered on theta=0.
|
---|
24 | */
|
---|
25 |
|
---|
26 | //***************************************************************
|
---|
27 | //++
|
---|
28 | // Class SphereThetaPhi
|
---|
29 | //
|
---|
30 | //
|
---|
31 | // include spherethetaphi.h
|
---|
32 | //--
|
---|
33 | //++
|
---|
34 | //
|
---|
35 | // Links Parents
|
---|
36 | //
|
---|
37 | // SphericalMap
|
---|
38 | //--
|
---|
39 |
|
---|
40 | /* --Methode-- */
|
---|
41 | //++
|
---|
42 | // Titre Constructors
|
---|
43 | //--
|
---|
44 | //++
|
---|
45 |
|
---|
46 | template <class T>
|
---|
47 | SphereThetaPhi<T>::SphereThetaPhi()
|
---|
48 | : NPhi_(), TNphi_(), Theta_(), pixels_()
|
---|
49 |
|
---|
50 | //--
|
---|
51 | {
|
---|
52 | InitNul();
|
---|
53 | }
|
---|
54 |
|
---|
55 |
|
---|
56 | /* --Methode-- */
|
---|
57 |
|
---|
58 | //++
|
---|
59 | template <class T>
|
---|
60 | SphereThetaPhi<T>::SphereThetaPhi(int_4 m)
|
---|
61 |
|
---|
62 | // m is the number of slices in theta on an hemisphere (the polar cap
|
---|
63 | // forms the first slice).
|
---|
64 | //--
|
---|
65 | {
|
---|
66 | InitNul();
|
---|
67 | Pixelize(m);
|
---|
68 | }
|
---|
69 |
|
---|
70 | template <class T>
|
---|
71 | SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s, bool share)
|
---|
72 | : NPhi_(s.NPhi_, share), TNphi_(s.TNphi_, share), Theta_(s.Theta_, share),
|
---|
73 | pixels_(s.pixels_ , share)
|
---|
74 | {
|
---|
75 |
|
---|
76 | NTheta_= s.NTheta_;
|
---|
77 | NPix_ = s.NPix_;
|
---|
78 | Omega_ = s.Omega_;
|
---|
79 | if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
|
---|
80 | }
|
---|
81 |
|
---|
82 | template <class T>
|
---|
83 | SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s)
|
---|
84 | : NPhi_(s.NPhi_), TNphi_(s.TNphi_), Theta_(s.Theta_), pixels_(s.pixels_)
|
---|
85 | {
|
---|
86 |
|
---|
87 | NTheta_= s.NTheta_;
|
---|
88 | NPix_ = s.NPix_;
|
---|
89 | Omega_ = s.Omega_;
|
---|
90 | if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
|
---|
91 | }
|
---|
92 |
|
---|
93 | //++
|
---|
94 | // Titre Destructor
|
---|
95 | //--
|
---|
96 | //++
|
---|
97 | template <class T>
|
---|
98 | SphereThetaPhi<T>::~SphereThetaPhi()
|
---|
99 |
|
---|
100 | //--
|
---|
101 | {}
|
---|
102 |
|
---|
103 | //++
|
---|
104 | // Titre Public Méthods
|
---|
105 | //--
|
---|
106 | template <class T>
|
---|
107 | void SphereThetaPhi<T>::InitNul()
|
---|
108 | //
|
---|
109 | {
|
---|
110 | NTheta_= 0;
|
---|
111 | NPix_ = 0;
|
---|
112 | // pixels_.Reset(); Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$
|
---|
113 | }
|
---|
114 |
|
---|
115 |
|
---|
116 | //++
|
---|
117 | template <class T>
|
---|
118 | void SphereThetaPhi<T>::Resize(int_4 m)
|
---|
119 | // re-pixelize the sphere
|
---|
120 | //--
|
---|
121 | {
|
---|
122 | InitNul();
|
---|
123 | Pixelize(m);
|
---|
124 | }
|
---|
125 |
|
---|
126 | template<class T>
|
---|
127 | void SphereThetaPhi<T>::CloneOrShare(const SphereThetaPhi<T>& a)
|
---|
128 | {
|
---|
129 |
|
---|
130 | NTheta_= a.NTheta_;
|
---|
131 | NPix_ = a.NPix_;
|
---|
132 | Omega_ = a.Omega_;
|
---|
133 | NPhi_.CloneOrShare(a.NPhi_);
|
---|
134 | TNphi_.CloneOrShare(a.TNphi_);
|
---|
135 | Theta_.CloneOrShare(a.Theta_);
|
---|
136 | pixels_.CloneOrShare(a.pixels_);
|
---|
137 | if (mInfo_) {delete mInfo_; mInfo_ = NULL;}
|
---|
138 | if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
|
---|
139 | }
|
---|
140 | template<class T>
|
---|
141 | void SphereThetaPhi<T>::Share(const SphereThetaPhi<T>& a)
|
---|
142 | {
|
---|
143 |
|
---|
144 | NTheta_= a.NTheta_;
|
---|
145 | NPix_ = a.NPix_;
|
---|
146 | Omega_ = a.Omega_;
|
---|
147 | NPhi_.Share(a.NPhi_);
|
---|
148 | TNphi_.Share(a.TNphi_);
|
---|
149 | Theta_.Share(a.Theta_);
|
---|
150 | pixels_.Share(a.pixels_);
|
---|
151 | if (mInfo_) {delete mInfo_; mInfo_ = NULL;}
|
---|
152 | if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
|
---|
153 | }
|
---|
154 |
|
---|
155 | ////////////////////////// methodes de copie/share
|
---|
156 | template<class T>
|
---|
157 | SphereThetaPhi<T>& SphereThetaPhi<T>::Set(const SphereThetaPhi<T>& a)
|
---|
158 | {
|
---|
159 | if (this != &a)
|
---|
160 | {
|
---|
161 |
|
---|
162 |
|
---|
163 | if (a.NbPixels() < 1)
|
---|
164 | throw RangeCheckError("SphereThetaPhi<T>::Set(a ) - Array a not allocated ! ");
|
---|
165 | if (NbPixels() < 1) CloneOrShare(a);
|
---|
166 | else CopyElt(a);
|
---|
167 | if (mInfo_) delete mInfo_;
|
---|
168 | mInfo_ = NULL;
|
---|
169 | if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
|
---|
170 | }
|
---|
171 | return(*this);
|
---|
172 | }
|
---|
173 |
|
---|
174 | template<class T>
|
---|
175 | SphereThetaPhi<T>& SphereThetaPhi<T>::CopyElt(const SphereThetaPhi<T>& a)
|
---|
176 | {
|
---|
177 | if (NbPixels() < 1)
|
---|
178 | throw RangeCheckError("SphereThetaPhi<T>::CopyElt(const SphereThetaPhi<T>& ) - Not Allocated Array ! ");
|
---|
179 | if (NbPixels() != a.NbPixels())
|
---|
180 | throw(SzMismatchError("SphereThetaPhi<T>::CopyElt(const SphereThetaPhi<T>&) SizeMismatch")) ;
|
---|
181 |
|
---|
182 | NTheta_= a.NTheta_;
|
---|
183 | NPix_ = a.NPix_;
|
---|
184 | Omega_ = a.Omega_;
|
---|
185 | int k;
|
---|
186 | for (k=0; k< NPix_; k++) pixels_(k) = a.pixels_(k);
|
---|
187 | for (k=0; k< a.NPhi_.Size(); k++) NPhi_(k) = a.NPhi_(k);
|
---|
188 | for (k=0; k< a.TNphi_.Size(); k++) TNphi_(k) = a.TNphi_(k);
|
---|
189 | for (k=0; k< a.Theta_.Size(); k++) Theta_(k) = a.Theta_(k);
|
---|
190 | return(*this);
|
---|
191 |
|
---|
192 | }
|
---|
193 |
|
---|
194 | /* --Methode-- */
|
---|
195 | //++
|
---|
196 | template <class T>
|
---|
197 | int_4 SphereThetaPhi<T>::NbPixels() const
|
---|
198 |
|
---|
199 | // Return total number of pixels
|
---|
200 | //--
|
---|
201 | {
|
---|
202 | return(NPix_);
|
---|
203 | }
|
---|
204 |
|
---|
205 | /* --Methode-- */
|
---|
206 | //++
|
---|
207 | template <class T>
|
---|
208 | T& SphereThetaPhi<T>::PixVal(int_4 k)
|
---|
209 |
|
---|
210 | // Return value of pixel with index k
|
---|
211 | //--
|
---|
212 | {
|
---|
213 | if((k < 0) || (k >= NPix_))
|
---|
214 | {
|
---|
215 | throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
|
---|
216 | }
|
---|
217 | return pixels_(k);
|
---|
218 | }
|
---|
219 |
|
---|
220 | //++
|
---|
221 | template <class T>
|
---|
222 | T const& SphereThetaPhi<T>::PixVal(int_4 k) const
|
---|
223 |
|
---|
224 | // Return value of pixel with index k
|
---|
225 | //--
|
---|
226 | {
|
---|
227 | if((k < 0) || (k >= NPix_))
|
---|
228 | {
|
---|
229 | throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
|
---|
230 | }
|
---|
231 | return *(pixels_.Data()+k);
|
---|
232 | }
|
---|
233 |
|
---|
234 | /* --Methode-- */
|
---|
235 | //++
|
---|
236 | template <class T>
|
---|
237 | bool SphereThetaPhi<T>::ContainsSph(double /*theta*/, double /*phi*/) const
|
---|
238 | //--
|
---|
239 | {
|
---|
240 | return(true);
|
---|
241 | }
|
---|
242 |
|
---|
243 | /* --Methode-- */
|
---|
244 | //++
|
---|
245 | template <class T>
|
---|
246 | int_4 SphereThetaPhi<T>::PixIndexSph(double theta, double phi) const
|
---|
247 |
|
---|
248 | // Return index of the pixel corresponding to
|
---|
249 | // direction (theta, phi).
|
---|
250 | //--
|
---|
251 | {
|
---|
252 | double dphi;
|
---|
253 | int i,k;
|
---|
254 | bool fgzn = false;
|
---|
255 |
|
---|
256 | if((theta > Pi) || (theta < 0.)) return(-1);
|
---|
257 | if((phi < 0.) || (phi > DeuxPi)) return(-1);
|
---|
258 | if(theta > Pi*0.5) {fgzn = true; theta = Pi-theta;}
|
---|
259 |
|
---|
260 | // La bande d'indice kt est limitée par les valeurs de theta contenues dans
|
---|
261 | // Theta_[kt] et Theta_[kt+1]
|
---|
262 | for( i=1; i< NTheta_; i++ )
|
---|
263 | if( theta < Theta_(i) ) break;
|
---|
264 |
|
---|
265 | dphi= DeuxPi/(double)NPhi_(i-1);
|
---|
266 |
|
---|
267 | if (fgzn) k= NPix_-TNphi_(i)+(int_4)(phi/dphi);
|
---|
268 | else k= TNphi_(i-1)+(int_4)(phi/dphi);
|
---|
269 | return(k);
|
---|
270 | }
|
---|
271 |
|
---|
272 | /* --Methode-- */
|
---|
273 | //++
|
---|
274 | template <class T>
|
---|
275 | void SphereThetaPhi<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
|
---|
276 |
|
---|
277 | // Return (theta,phi) coordinates of middle of pixel with index k
|
---|
278 | //--
|
---|
279 | {
|
---|
280 | int i;
|
---|
281 | bool fgzn = false;
|
---|
282 |
|
---|
283 | if((k < 0) || (k >= NPix_)) {theta = -99999.; phi = -99999.; return; }
|
---|
284 | if( k >= NPix_/2) {fgzn = true; k = NPix_-1-k;}
|
---|
285 |
|
---|
286 | // recupère l'indice i de la tranche qui contient le pixel k
|
---|
287 | for( i=0; i< NTheta_; i++ )
|
---|
288 | if( k < TNphi_(i+1) ) break;
|
---|
289 |
|
---|
290 | // angle theta
|
---|
291 | theta= 0.5*(Theta_(i)+Theta_(i+1));
|
---|
292 | if (fgzn) theta= Pi-theta;
|
---|
293 |
|
---|
294 | // angle phi
|
---|
295 | k -= TNphi_(i);
|
---|
296 | phi= DeuxPi/(double)NPhi_(i)*(double)(k+.5);
|
---|
297 | if (fgzn) phi= DeuxPi-phi;
|
---|
298 | }
|
---|
299 |
|
---|
300 | template <class T>
|
---|
301 | T SphereThetaPhi<T>::SetPixels(T v)
|
---|
302 | {
|
---|
303 | pixels_.Reset(v);
|
---|
304 | return(v);
|
---|
305 | }
|
---|
306 |
|
---|
307 | //++
|
---|
308 | template <class T>
|
---|
309 | double SphereThetaPhi<T>::PixSolAngle(int_4 /*dummy*/) const
|
---|
310 |
|
---|
311 | // Pixel Solid angle (steradians)
|
---|
312 | // All the pixels have the same solid angle. The dummy argument is
|
---|
313 | // for compatibility with eventual pixelizations which would not
|
---|
314 | // fulfil this requirement.
|
---|
315 | //--
|
---|
316 | {
|
---|
317 | return Omega_;
|
---|
318 | }
|
---|
319 |
|
---|
320 | /* --Methode-- */
|
---|
321 | //++
|
---|
322 | template <class T>
|
---|
323 | void SphereThetaPhi<T>::Limits(int_4 k,double& tetMin,double& tetMax,double& phiMin,double& phiMax)
|
---|
324 |
|
---|
325 | // Return values of theta,phi which limit the pixel with index k
|
---|
326 | //--
|
---|
327 | {
|
---|
328 | int j;
|
---|
329 | double dphi;
|
---|
330 | bool fgzn= false;
|
---|
331 |
|
---|
332 | if((k < 0) || (k >= NPix_)) {
|
---|
333 | tetMin= -99999.;
|
---|
334 | phiMin= -99999.;
|
---|
335 | tetMax= -99999.;
|
---|
336 | phiMax= -99999.;
|
---|
337 | return;
|
---|
338 | }
|
---|
339 |
|
---|
340 | // si on se trouve dans l'hémisphère Sud
|
---|
341 | if(k >= NPix_/2) {
|
---|
342 | fgzn= true;
|
---|
343 | k= NPix_-1-k;
|
---|
344 | }
|
---|
345 |
|
---|
346 | // on recupere l'indice i de la tranche qui contient le pixel k
|
---|
347 | int i;
|
---|
348 | for( i=0; i< NTheta_; i++ )
|
---|
349 | if(k < TNphi_(i+1)) break;
|
---|
350 |
|
---|
351 | // valeurs limites de theta dans l'hemisphere Nord
|
---|
352 | tetMin= Theta_(i);
|
---|
353 | tetMax= Theta_(i+1);
|
---|
354 | // valeurs limites de theta dans l'hemisphere Sud
|
---|
355 | if (fgzn) {
|
---|
356 | tetMin= Pi - Theta_(i+1);
|
---|
357 | tetMax= Pi - Theta_(i);
|
---|
358 | }
|
---|
359 |
|
---|
360 | // pixel correspondant dans l'hemisphere Nord
|
---|
361 | if (fgzn) k= TNphi_(i+1)-k+TNphi_(i)-1;
|
---|
362 |
|
---|
363 | // indice j de discretisation ( phi= j*dphi )
|
---|
364 | j= k-TNphi_(i);
|
---|
365 | dphi= DeuxPi/(double)NPhi_(i);
|
---|
366 |
|
---|
367 | // valeurs limites de phi
|
---|
368 | phiMin= dphi*(double)(j);
|
---|
369 | phiMax= dphi*(double)(j+1);
|
---|
370 | return;
|
---|
371 | }
|
---|
372 |
|
---|
373 | /* --Methode-- */
|
---|
374 | //++
|
---|
375 | template <class T>
|
---|
376 | uint_4 SphereThetaPhi<T>::NbThetaSlices() const
|
---|
377 |
|
---|
378 | // Return number of theta-slices on the sphere
|
---|
379 | //--
|
---|
380 | {
|
---|
381 | if (NTheta_<=0)
|
---|
382 | {
|
---|
383 | throw PException(" sphere not pixelized, NbSlice=0 ");
|
---|
384 | }
|
---|
385 | return( 2*NTheta_);
|
---|
386 | }
|
---|
387 |
|
---|
388 | /* --Methode-- */
|
---|
389 | //++
|
---|
390 | template <class T>
|
---|
391 | int_4 SphereThetaPhi<T>::NPhi(int_4 kt) const
|
---|
392 |
|
---|
393 | // Return number of pixels in phi-direction of the kt-th slice
|
---|
394 | //--
|
---|
395 | {
|
---|
396 | int nbpix;
|
---|
397 | // verification
|
---|
398 | if((kt < 0) || (kt >= 2*NTheta_)) return(-1);
|
---|
399 |
|
---|
400 | // si on se trouve dans l'hemisphere Sud
|
---|
401 | if(kt >= NTheta_) {
|
---|
402 | kt= 2*NTheta_-1-kt;
|
---|
403 | }
|
---|
404 |
|
---|
405 | // nombre de pixels
|
---|
406 | nbpix= NPhi_(kt);
|
---|
407 | return(nbpix);
|
---|
408 | }
|
---|
409 |
|
---|
410 |
|
---|
411 | /* --Methode-- */
|
---|
412 | //++
|
---|
413 | template <class T>
|
---|
414 | void SphereThetaPhi<T>::Theta(int_4 kt,double& tetMin,double& tetMax)
|
---|
415 |
|
---|
416 | // Return theta values which limit the slice kt
|
---|
417 | //--
|
---|
418 | {
|
---|
419 | bool fgzn= false;
|
---|
420 | // verification
|
---|
421 | if( (kt< 0) || (kt>= 2*NTheta_) ) {
|
---|
422 | tetMin= -99999.;
|
---|
423 | tetMax= -99999.;
|
---|
424 | return;
|
---|
425 | }
|
---|
426 |
|
---|
427 | // si on se trouve dans l'hemisphere Sud
|
---|
428 | if( kt >= NTheta_ ) {
|
---|
429 | fgzn= true;
|
---|
430 | kt= 2*NTheta_-1-kt;
|
---|
431 | }
|
---|
432 |
|
---|
433 | // valeurs limites de theta dans l'hemisphere Nord
|
---|
434 | tetMin= Theta_(kt);
|
---|
435 | tetMax= Theta_(kt+1);
|
---|
436 | // valeurs limites de theta dans l'hemisphere Sud
|
---|
437 | if (fgzn) {
|
---|
438 | tetMin= Pi - Theta_(kt+1);
|
---|
439 | tetMax= Pi - Theta_(kt);
|
---|
440 | }
|
---|
441 | }
|
---|
442 |
|
---|
443 | /* --Methode-- */
|
---|
444 | //++
|
---|
445 | template <class T>
|
---|
446 | void SphereThetaPhi<T>::Phi(int_4 kt,int_4 jp,double& phiMin,double& phiMax)
|
---|
447 |
|
---|
448 | // Return values of phi which limit the jp-th pixel of the kt-th slice
|
---|
449 | //--
|
---|
450 | {
|
---|
451 | // verification
|
---|
452 | if((kt < 0) || (kt >= 2*NTheta_)) {
|
---|
453 | phiMin= -99999.;
|
---|
454 | phiMax= -99999.;
|
---|
455 | return;
|
---|
456 | }
|
---|
457 |
|
---|
458 | // si on se trouve dans l'hemisphere Sud
|
---|
459 | if(kt >= NTheta_) kt= 2*NTheta_-1-kt;
|
---|
460 |
|
---|
461 | // verifie si la tranche kt contient au moins jp pixels
|
---|
462 | if( (jp< 0) || (jp >= NPhi_(kt)) ) {
|
---|
463 | phiMin= -88888.;
|
---|
464 | phiMax= -88888.;
|
---|
465 | return;
|
---|
466 | }
|
---|
467 |
|
---|
468 | double dphi= DeuxPi/(double)NPhi_(kt);
|
---|
469 | phiMin= dphi*(double)(jp);
|
---|
470 | phiMax= dphi*(double)(jp+1);
|
---|
471 | return;
|
---|
472 | }
|
---|
473 |
|
---|
474 | /* --Methode-- */
|
---|
475 | //++
|
---|
476 | template <class T>
|
---|
477 | int_4 SphereThetaPhi<T>::Index(int_4 kt,int_4 jp) const
|
---|
478 |
|
---|
479 | // Return pixel index with sequence index jp in the slice kt
|
---|
480 | //--
|
---|
481 | {
|
---|
482 | int k;
|
---|
483 | bool fgzn= false;
|
---|
484 |
|
---|
485 | // si on se trouve dans l'hemisphere Sud
|
---|
486 | if(kt >= NTheta_) {
|
---|
487 | fgzn= true;
|
---|
488 | kt= 2*NTheta_-1-kt;
|
---|
489 | }
|
---|
490 |
|
---|
491 | // si la tranche kt contient au moins i pixels
|
---|
492 | if( (jp>=0) && (jp<NPhi_(kt)) )
|
---|
493 | {
|
---|
494 | // dans l'hemisphere Sud
|
---|
495 | if (fgzn) k= NPix_-TNphi_(kt+1)+jp;
|
---|
496 | // dans l'hemisphere Nord
|
---|
497 | else k= TNphi_(kt)+jp;
|
---|
498 | }
|
---|
499 | else
|
---|
500 | {
|
---|
501 | k= 9999;
|
---|
502 | printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
|
---|
503 | }
|
---|
504 | return(k);
|
---|
505 | }
|
---|
506 |
|
---|
507 | /* --Methode-- */
|
---|
508 | //++
|
---|
509 | template <class T>
|
---|
510 | void SphereThetaPhi<T>::ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp)
|
---|
511 |
|
---|
512 | // Return indices kt (theta) and jp (phi) of pixel with index k
|
---|
513 | //--
|
---|
514 | {
|
---|
515 | bool fgzn= false;
|
---|
516 | // si on se trouve dans l'hemisphere Sud
|
---|
517 | if(k >= NPix_/2)
|
---|
518 | {
|
---|
519 | fgzn= true;
|
---|
520 | k= NPix_-1-k;
|
---|
521 | }
|
---|
522 |
|
---|
523 | // on recupere l'indice kt de la tranche qui contient le pixel k
|
---|
524 | int i;
|
---|
525 | for(i = 0; i < NTheta_; i++)
|
---|
526 | if(k < TNphi_(i+1)) break;
|
---|
527 |
|
---|
528 | // indice kt de tranche
|
---|
529 | if (fgzn) kt= 2*NTheta_-1-i;
|
---|
530 | else kt= i;
|
---|
531 |
|
---|
532 | // indice jp de pixel
|
---|
533 | if (fgzn) jp= TNphi_(i+1)-k-1;
|
---|
534 | else jp= k-TNphi_(i);
|
---|
535 | }
|
---|
536 | //++
|
---|
537 | template <class T>
|
---|
538 | void SphereThetaPhi<T>::Pixelize(int_4 m)
|
---|
539 |
|
---|
540 | // achieve the splitting into pixels (m has the same signification
|
---|
541 | // as for the constructor)
|
---|
542 | //
|
---|
543 | // Each theta-slice of the north hemisphere will be spitted starting f
|
---|
544 | // from phi=0 ...
|
---|
545 | //
|
---|
546 | // South hemisphere is scanned in the same direction according to phi
|
---|
547 | // and from equator to the pole (the pixel following the last one of
|
---|
548 | // the slice closest to the equator with z>0, is the pixel with lowest
|
---|
549 | // phi of the slice closest of the equator with z<0).
|
---|
550 | //--
|
---|
551 | {
|
---|
552 | int ntotpix,j;
|
---|
553 |
|
---|
554 | // Decodage et controle des arguments d'appel :
|
---|
555 | // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
|
---|
556 | if (m < 2) m = 2;
|
---|
557 | if (m > 16384) m = 16384;
|
---|
558 |
|
---|
559 | // On memorise les arguments d'appel
|
---|
560 | NTheta_ = m;
|
---|
561 |
|
---|
562 | // On commence par decouper l'hemisphere z>0.
|
---|
563 | // Creation des vecteurs contenant :
|
---|
564 | // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
|
---|
565 | // Theta_= new double[m+1];
|
---|
566 | Theta_.ReSize(m+1);
|
---|
567 |
|
---|
568 | // Le nombre de pixels en phi de chacune des bandes en theta
|
---|
569 | // NPhi_ = new int_4[m];
|
---|
570 | // une taille de m suffit, mais je mets m+1 pour que les 3 tableaux aient
|
---|
571 | // la meme taille pour une manipulation plus faciles par la librairie
|
---|
572 | // cfitsio -- GLM (13-04-00)
|
---|
573 | NPhi_.ReSize(m+1);
|
---|
574 |
|
---|
575 | // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
|
---|
576 | // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
|
---|
577 | // TNphi_= new int_4[m+1];
|
---|
578 | TNphi_.ReSize(m+1);
|
---|
579 |
|
---|
580 | // Calcul du nombre total de pixels dans chaque bande pour optimiser
|
---|
581 | // le rapport largeur/hauteur des pixels
|
---|
582 |
|
---|
583 | //calotte polaire
|
---|
584 | TNphi_(0)= 0;
|
---|
585 | NPhi_(0) = 1;
|
---|
586 |
|
---|
587 | //bandes jusqu'a l'equateur
|
---|
588 | for(j = 1; j < m; j++)
|
---|
589 | {
|
---|
590 | TNphi_(j)= TNphi_(j-1)+NPhi_(j-1);
|
---|
591 | NPhi_(j) = (int_4)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
|
---|
592 | }
|
---|
593 |
|
---|
594 | // Nombre total de pixels sur l'hemisphere
|
---|
595 | ntotpix = TNphi_(m-1)+NPhi_(m-1);
|
---|
596 | TNphi_(m)= ntotpix;
|
---|
597 | // et sur la sphere entiere
|
---|
598 | NPix_= 2*ntotpix;
|
---|
599 |
|
---|
600 | // Creation et initialisation du vecteur des contenus des pixels
|
---|
601 | pixels_.ReSize(NPix_);
|
---|
602 | pixels_.Reset();
|
---|
603 |
|
---|
604 | // Determination des limites des bandes en theta :
|
---|
605 | // omeg est l'angle solide couvert par chaque pixel,
|
---|
606 | // une bande donnee kt couvre un angle solide NPhi_[kt]*omeg
|
---|
607 | // egal a 2* Pi*(cos Theta_[kt]-cos Theta_[kt+1]). De meme, l'angle solide
|
---|
608 | //de la
|
---|
609 | // calotte allant du pole a la limite haute de la bande kt vaut
|
---|
610 | // 2* Pi*(1.-cos Theta_[kt+1])= TNphi_[kt]*omeg...
|
---|
611 |
|
---|
612 | double omeg2pi= 1./(double)ntotpix;
|
---|
613 | Omega_ = omeg2pi*DeuxPi;
|
---|
614 |
|
---|
615 | for(j=0; j <= m; j++)
|
---|
616 | {
|
---|
617 | Theta_(j)= acos(1.-(double)TNphi_(j)*omeg2pi);
|
---|
618 | }
|
---|
619 | }
|
---|
620 |
|
---|
621 | //++
|
---|
622 | template <class T>
|
---|
623 | void SphereThetaPhi<T>::GetThetaSlice(int_4 index,r_8& theta, TVector<r_8>& phi, TVector<T>& value) const
|
---|
624 |
|
---|
625 | // For a theta-slice with index 'index', return :
|
---|
626 | // the corresponding "theta"
|
---|
627 | // a vector containing the phi's of the pixels of the slice
|
---|
628 | // a vector containing the corresponding values of pixels
|
---|
629 | //--
|
---|
630 |
|
---|
631 | {
|
---|
632 |
|
---|
633 | if(index < 0 || index >= NbThetaSlices())
|
---|
634 | {
|
---|
635 | throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
|
---|
636 | }
|
---|
637 |
|
---|
638 | int iring= Index(index,0);
|
---|
639 | int lring = NPhi(index);
|
---|
640 |
|
---|
641 | phi.ReSize(lring);
|
---|
642 | value.ReSize(lring);
|
---|
643 | double Te= 0.;
|
---|
644 | double Fi= 0.;
|
---|
645 | for(int kk = 0; kk < lring; kk++)
|
---|
646 | {
|
---|
647 | PixThetaPhi(kk+iring,Te,Fi);
|
---|
648 | phi(kk)= Fi;
|
---|
649 | value(kk)= PixVal(kk+iring);
|
---|
650 | }
|
---|
651 | theta= Te;
|
---|
652 | }
|
---|
653 |
|
---|
654 | //++
|
---|
655 | template <class T>
|
---|
656 | void SphereThetaPhi<T>::GetThetaSlice(int_4 index,r_8& theta, r_8& phi0,TVector<int_4>& pixelIndices, TVector<T>& value) const
|
---|
657 |
|
---|
658 | // For a theta-slice with index 'index', return :
|
---|
659 | // the corresponding "theta"
|
---|
660 | // the corresponding "phi" for first pixel of the slice
|
---|
661 | // a vector containing the indices of the pixels of the slice
|
---|
662 | // (equally distributed in phi)
|
---|
663 | // a vector containing the corresponding values of pixels
|
---|
664 | //--
|
---|
665 |
|
---|
666 | {
|
---|
667 |
|
---|
668 | if(index < 0 || index >= NbThetaSlices())
|
---|
669 | {
|
---|
670 | throw RangeCheckError("SphereThetaPhi::PIxVal Pixel index out of range");
|
---|
671 |
|
---|
672 | }
|
---|
673 |
|
---|
674 | int iring= Index(index,0);
|
---|
675 | int lring = NPhi(index);
|
---|
676 |
|
---|
677 | pixelIndices.ReSize(lring);
|
---|
678 | value.ReSize(lring);
|
---|
679 | double Te= 0.;
|
---|
680 | double Fi= 0.;
|
---|
681 | for(int kk = 0; kk < lring; kk++)
|
---|
682 | {
|
---|
683 | pixelIndices(kk)=kk+iring ;
|
---|
684 | value(kk)= PixVal(kk+iring);
|
---|
685 | }
|
---|
686 | PixThetaPhi(iring,theta,phi0);
|
---|
687 | }
|
---|
688 |
|
---|
689 |
|
---|
690 |
|
---|
691 |
|
---|
692 | template <class T>
|
---|
693 | void SphereThetaPhi<T>::print(ostream& os) const
|
---|
694 | {
|
---|
695 | if(mInfo_) os << " DVList Info= " << *mInfo_ << endl;
|
---|
696 | //
|
---|
697 | os << " NTheta_= " << NTheta_ << endl;
|
---|
698 | os << " NPix_ = " << NPix_ << endl;
|
---|
699 | os << " Omega_ = " << Omega_ << endl;
|
---|
700 |
|
---|
701 | os << " contenu de NPhi_ : ";
|
---|
702 | int i;
|
---|
703 | for(i=0; i < NTheta_; i++)
|
---|
704 | {
|
---|
705 | if(i%5 == 0) os << endl;
|
---|
706 | os << NPhi_(i) <<", ";
|
---|
707 | }
|
---|
708 | os << endl;
|
---|
709 |
|
---|
710 | os << " contenu de Theta_ : ";
|
---|
711 | for(i=0; i < NTheta_+1; i++)
|
---|
712 | {
|
---|
713 | if(i%5 == 0) os << endl;
|
---|
714 | os << Theta_(i) <<", ";
|
---|
715 | }
|
---|
716 | os << endl;
|
---|
717 |
|
---|
718 | os << " contenu de TNphi_ : ";
|
---|
719 | for(i=0; i < NTheta_+1; i++)
|
---|
720 | {
|
---|
721 | if(i%5 == 0) os << endl;
|
---|
722 | os << TNphi_(i) <<", ";
|
---|
723 | }
|
---|
724 | os << endl;
|
---|
725 |
|
---|
726 | os << " contenu de pixels : ";
|
---|
727 | for(i=0; i < NPix_; i++)
|
---|
728 | {
|
---|
729 | if(i%5 == 0) os << endl;
|
---|
730 | os << pixels_(i) <<", ";
|
---|
731 | }
|
---|
732 | os << endl;
|
---|
733 | }
|
---|
734 |
|
---|
735 | // ...... Operations de calcul ......
|
---|
736 |
|
---|
737 |
|
---|
738 | //! Fill a SphereThetaPhi with a constant value \b a
|
---|
739 | template <class T>
|
---|
740 | SphereThetaPhi<T>& SphereThetaPhi<T>::SetT(T a)
|
---|
741 | {
|
---|
742 | if (NbPixels() < 1)
|
---|
743 | throw RangeCheckError("SphereThetaPhi<T>::SetT(T ) - SphereThetaPhi not dimensionned ! ");
|
---|
744 | pixels_ = a;
|
---|
745 | return (*this);
|
---|
746 | }
|
---|
747 |
|
---|
748 | /*! Add a constant value \b x to a SphereThetaPhi */
|
---|
749 | template <class T>
|
---|
750 | SphereThetaPhi<T>& SphereThetaPhi<T>::Add(T a)
|
---|
751 | {
|
---|
752 | if (NbPixels()< 1)
|
---|
753 | throw RangeCheckError("SphereThetaPhi<T>::Add(T ) - SphereThetaPhi not dimensionned ! ");
|
---|
754 | pixels_ += a;
|
---|
755 | return (*this);
|
---|
756 | }
|
---|
757 |
|
---|
758 | /*! Substract a constant value \b a to a SphereThetaPhi */
|
---|
759 | template <class T>
|
---|
760 | SphereThetaPhi<T>& SphereThetaPhi<T>::Sub(T a,bool fginv)
|
---|
761 | {
|
---|
762 | if (NbPixels()< 1)
|
---|
763 | throw RangeCheckError("SphereThetaPhi<T>::Sub(T ) - SphereThetaPhi not dimensionned ! ");
|
---|
764 | pixels_.Sub(a,fginv);
|
---|
765 | return (*this);
|
---|
766 | }
|
---|
767 |
|
---|
768 | /*! multiply a SphereThetaPhi by a constant value \b a */
|
---|
769 | template <class T>
|
---|
770 | SphereThetaPhi<T>& SphereThetaPhi<T>::Mul(T a)
|
---|
771 | {
|
---|
772 | if (NbPixels()< 1)
|
---|
773 | throw RangeCheckError("SphereThetaPhi<T>::Mul(T ) - SphereThetaPhi not dimensionned ! ");
|
---|
774 | pixels_ *= a;
|
---|
775 | return (*this);
|
---|
776 | }
|
---|
777 |
|
---|
778 | /*! divide a SphereThetaPhi by a constant value \b a */
|
---|
779 | template <class T>
|
---|
780 | SphereThetaPhi<T>& SphereThetaPhi<T>::Div(T a)
|
---|
781 | {
|
---|
782 | if (NbPixels()< 1)
|
---|
783 | throw RangeCheckError("SphereThetaPhi<T>::Div(T ) - SphereThetaPhi not dimensionned ! ");
|
---|
784 | pixels_ /= a;
|
---|
785 | return (*this);
|
---|
786 | }
|
---|
787 |
|
---|
788 | // >>>> Operations avec 2nd membre de type SphereThetaPhi
|
---|
789 | //! Add two SphereThetaPhi
|
---|
790 |
|
---|
791 | template <class T>
|
---|
792 | SphereThetaPhi<T>& SphereThetaPhi<T>::AddElt(const SphereThetaPhi<T>& a)
|
---|
793 | {
|
---|
794 | if (NbPixels()!= a.NbPixels())
|
---|
795 | {
|
---|
796 | throw(SzMismatchError("SphereThetaPhi<T>::AddElt(const SphereThetaPhi<T>&) SizeMismatch")) ;
|
---|
797 | }
|
---|
798 | pixels_ += a.pixels_;
|
---|
799 | return (*this);
|
---|
800 | }
|
---|
801 |
|
---|
802 | //! Substract two SphereThetaPhi
|
---|
803 | template <class T>
|
---|
804 | SphereThetaPhi<T>& SphereThetaPhi<T>::SubElt(const SphereThetaPhi<T>& a)
|
---|
805 | {
|
---|
806 | if (NbPixels()!= a.NbPixels())
|
---|
807 | {
|
---|
808 | throw(SzMismatchError("SphereThetaPhi<T>::SubElt(const SphereThetaPhi<T>&) SizeMismatch")) ;
|
---|
809 | }
|
---|
810 | pixels_ -= a.pixels_;
|
---|
811 | return (*this);
|
---|
812 | }
|
---|
813 |
|
---|
814 | //! Multiply two SphereThetaPhi (elements by elements)
|
---|
815 | template <class T>
|
---|
816 | SphereThetaPhi<T>& SphereThetaPhi<T>::MulElt(const SphereThetaPhi<T>& a)
|
---|
817 | {
|
---|
818 | if (NbPixels()!= a.NbPixels())
|
---|
819 | {
|
---|
820 | throw(SzMismatchError("SphereThetaPhi<T>::MulElt(const SphereThetaPhi<T>&) SizeMismatch")) ;
|
---|
821 | }
|
---|
822 | pixels_ *= a.pixels_;
|
---|
823 | return (*this);
|
---|
824 | }
|
---|
825 |
|
---|
826 | //! Divide two SphereThetaPhi (elements by elements) - No protection for divide by 0
|
---|
827 | template <class T>
|
---|
828 | SphereThetaPhi<T>& SphereThetaPhi<T>::DivElt(const SphereThetaPhi<T>& a)
|
---|
829 | {
|
---|
830 | if (NbPixels()!= a.NbPixels())
|
---|
831 | {
|
---|
832 | throw(SzMismatchError("SphereThetaPhi<T>::DivElt(const SphereThetaPhi<T>&) SizeMismatch")) ;
|
---|
833 | }
|
---|
834 | pixels_ /= a.pixels_;
|
---|
835 | return (*this);
|
---|
836 | }
|
---|
837 |
|
---|
838 |
|
---|
839 |
|
---|
840 | #ifdef __CXX_PRAGMA_TEMPLATES__
|
---|
841 | #pragma define_template SphereThetaPhi<int_4>
|
---|
842 | #pragma define_template SphereThetaPhi<r_8>
|
---|
843 | #pragma define_template SphereThetaPhi<r_4>
|
---|
844 | #pragma define_template SphereThetaPhi< complex<r_4> >
|
---|
845 | #pragma define_template SphereThetaPhi< complex<r_8> >
|
---|
846 | #endif
|
---|
847 | #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
|
---|
848 | template class SphereThetaPhi<int_4>;
|
---|
849 | template class SphereThetaPhi<r_8>;
|
---|
850 | template class SphereThetaPhi<r_4>;
|
---|
851 | template class SphereThetaPhi< complex<r_4> >;
|
---|
852 | template class SphereThetaPhi< complex<r_8> >;
|
---|
853 | #endif
|
---|