source: Sophya/trunk/SophyaLib/SkyMap/spherehealpix.cc@ 1287

Last change on this file since 1287 was 1217, checked in by ansari, 25 years ago

doc dans les .cc
File size: 14.3 KB
Line 
1#include "machdefs.h"
2#include <math.h>
3#include <complex>
4
5#include "pexceptions.h"
6#include "fiondblock.h"
7#include "spherehealpix.h"
8#include "strutil.h"
9
10extern "C"
11{
12#include <stdio.h>
13#include <stdlib.h>
14#include <unistd.h>
15}
16
17using namespace SOPHYA;
18
19/*!
20 \class SOPHYA::SphereHEALPix
21
22*******************************************************************
23 Pixelisation Gorski
24
25\verbatim
26
27 -----------------------------------------------------------------------
28 version 0.8.2 Aug97 TAC Eric Hivon, Kris Gorski
29 -----------------------------------------------------------------------
30
31 the sphere is split in 12 diamond-faces containing nside**2 pixels each
32 the numbering of the pixels (in the nested scheme) is similar to
33 quad-cube
34 In each face the first pixel is in the lowest corner of the diamond
35 the faces are (x,y) coordinate on each face
36
37 . . . . <--- North Pole
38 / \ / \ / \ / \ ^ ^
39 . 0 . 1 . 2 . 3 . <--- z = 2/3 \ /
40 \ / \ / \ / \ / y \ / x
41 4 . 5 . 6 . 7 . 4 <--- equator \ /
42 / \ / \ / \ / \ \/
43 . 8 . 9 .10 .11 . <--- z = -2/3 (0,0) : lowest corner
44 \ / \ / \ / \ /
45 . . . . <--- South Pole
46\endverbatim
47 phi:0 2Pi
48
49 in the ring scheme pixels are numbered along the parallels
50 the first parallel is the one closest to the north pole and so on
51 on each parallel, pixels are numbered starting from the one closest
52 to phi = 0
53
54 nside MUST be a power of 2 (<= 8192)
55
56*/
57
58/* --Methode-- */
59
60template<class T>
61SphereHEALPix<T>::SphereHEALPix() : pixels_(), sliceBeginIndex_(),
62 sliceLenght_()
63
64
65{
66 InitNul();
67}
68
69/*! \fn SOPHYA::SphereHEALPix::SphereHEALPix(int_4 m)
70
71 \param <m> : "nside" of the Healpix algorithm
72
73 The total number of pixels will be Npix = 12*nside**2
74
75 nside MUST be a power of 2 (<= 8192)
76*/
77
78template<class T>
79SphereHEALPix<T>::SphereHEALPix(int_4 m)
80
81{
82
83 if(m <= 0 || m > 8192)
84 {
85 cout << "SphereHEALPix : m hors bornes [0,8192], m= " << m << endl;
86 throw RangeCheckError("SphereHEALPix<T>::SphereHEALPix() - Out of bound nside (< 8192)!");
87 }
88 // verifier que m est une puissance de deux
89 int x= m;
90 while(x%2 == 0) x/=2;
91 if(x != 1)
92 {
93 cout<<"SphereHEALPix: m doit etre une puissance de deux, m= "<<m<<endl;
94 throw ParmError("SphereHEALPix<T>::SphereHEALPix() - nside != 2^n !");
95 }
96 InitNul();
97 // SetTemp(false);
98 Pixelize(m);
99 SetThetaSlices();
100}
101//++
102template<class T>
103SphereHEALPix<T>::SphereHEALPix(const SphereHEALPix<T>& s, bool share)
104 : pixels_(s.pixels_, share), sliceBeginIndex_(s.sliceBeginIndex_, share),
105 sliceLenght_(s.sliceLenght_, share)
106// copy constructor
107//--
108{
109 nSide_= s.nSide_;
110 nPix_ = s.nPix_;
111 omeg_ = s.omeg_;
112 if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
113}
114//++
115template<class T>
116SphereHEALPix<T>::SphereHEALPix(const SphereHEALPix<T>& s)
117 : pixels_(s.pixels_), sliceBeginIndex_(s.sliceBeginIndex_),
118 sliceLenght_(s.sliceLenght_)
119// copy constructor
120//--
121{
122 nSide_= s.nSide_;
123 nPix_ = s.nPix_;
124 omeg_ = s.omeg_;
125 if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
126 // CloneOrShare(s);
127}
128
129template<class T>
130void SphereHEALPix<T>::CloneOrShare(const SphereHEALPix<T>& a)
131{
132 nSide_= a.nSide_;
133 nPix_ = a.nPix_;
134 omeg_ = a.omeg_;
135 pixels_.CloneOrShare(a.pixels_);
136 sliceBeginIndex_.CloneOrShare(a.sliceBeginIndex_);
137 sliceLenght_.CloneOrShare(a.sliceLenght_);
138
139 // pas forcement a conserver, pas forcement a cet endroit (GLM)
140 // if (a.IsTemp() ) SetTemp(true);
141}
142
143////////////////////////// methodes de copie/share
144template<class T>
145SphereHEALPix<T>& SphereHEALPix<T>::Set(const SphereHEALPix<T>& a)
146{
147 if (this != &a)
148 {
149
150 if (a.NbPixels() < 1)
151 throw RangeCheckError("SphereHEALPix<T>::Set(a ) - Array a not allocated ! ");
152 if (NbPixels() < 1) CloneOrShare(a);
153 else CopyElt(a);
154
155
156 if (mInfo_) delete mInfo_;
157 mInfo_ = NULL;
158 if (a.mInfo_) mInfo_ = new DVList(*(a.mInfo_));
159 }
160 return(*this);
161}
162
163template<class T>
164SphereHEALPix<T>& SphereHEALPix<T>::CopyElt(const SphereHEALPix<T>& a)
165{
166 if (NbPixels() < 1)
167 throw RangeCheckError("SphereHEALPix<T>::CopyElt(const SphereHEALPix<T>& ) - Not Allocated Array ! ");
168 if (NbPixels() != a.NbPixels())
169 throw(SzMismatchError("SphereHEALPix<T>::CopyElt(const SphereHEALPix<T>&) SizeMismatch")) ;
170 nSide_= a.nSide_;
171 nPix_ = a.nPix_;
172 omeg_ = a.omeg_;
173 int k;
174 for (k=0; k< nPix_; k++) pixels_(k) = a.pixels_(k);
175 for (k=0; k< a.sliceBeginIndex_.Size(); k++) sliceBeginIndex_(k) = a.sliceBeginIndex_(k);
176 for (k=0; k< a.sliceLenght_.Size(); k++) sliceLenght_(k) = a.sliceLenght_(k);
177 return(*this);
178}
179//++
180// Titre Destructor
181//--
182//++
183template<class T>
184SphereHEALPix<T>::~SphereHEALPix()
185
186//--
187{
188}
189
190//++
191// Titre Public Methods
192//--
193
194/*! \fn void SOPHYA::SphereHEALPix::Resize(int_4 m)
195 \param <m> "nside" of the Gorski algorithm
196
197 The total number of pixels will be Npix = 12*nside**2
198
199 nside MUST be a power of 2 (<= 8192)
200*/
201template<class T>
202void SphereHEALPix<T>::Resize(int_4 m)
203
204
205{
206 if (m<=0 || m> 8192) {
207 cout << "SphereHEALPix : m hors bornes [0,8192], m= " << m << endl;
208 exit(1);
209 }
210 // verifier que m est une puissance de deux
211 int x= m;
212 while (x%2==0) x/=2;
213 if(x != 1)
214 {
215 cout<<"SphereHEALPix: m doit etre une puissance de deux, m= "<<m<<endl;
216 exit(1);
217 }
218 InitNul();
219 Pixelize(m);
220 SetThetaSlices();
221}
222
223template<class T>
224void SphereHEALPix<T>::Pixelize( int_4 m)
225
226// prépare la pixelisation Gorski (m a la même signification
227// que pour le constructeur)
228//
229//
230{
231 // On memorise les arguments d'appel
232 nSide_= m;
233
234 // Nombre total de pixels sur la sphere entiere
235 nPix_= 12*nSide_*nSide_;
236
237 // pour le moment les tableaux qui suivent seront ranges dans l'ordre
238 // de l'indexation GORSKY "RING"
239 // on pourra ulterieurement changer de strategie et tirer profit
240 // de la dualite d'indexation GORSKY (RING et NEST) : tout dependra
241 // de pourquoi c'est faire
242
243 // Creation et initialisation du vecteur des contenus des pixels
244 pixels_.ReSize(nPix_);
245 pixels_.Reset();
246
247 // solid angle per pixel
248 omeg_= 4.0*Pi/nPix_;
249}
250
251template<class T>
252void SphereHEALPix<T>::InitNul()
253//
254// initialise à zéro les variables de classe
255{
256 nSide_= 0;
257 nPix_ = 0;
258 omeg_ = 0.;
259// pixels_.Reset(); - Il ne faut pas mettre les pixels a zero si share !
260}
261
262/* --Methode-- */
263/* Nombre de pixels du decoupage */
264/*! \fn int_4 SOPHYA::SphereHEALPix::NbPixels() const
265
266 Return number of pixels of the splitting
267*/
268template<class T>
269int_4 SphereHEALPix<T>::NbPixels() const
270{
271 return(nPix_);
272}
273
274
275/*! \fn uint_4 SOPHYA::SphereHEALPix::NbThetaSlices() const
276
277 \return number of slices in theta direction on the sphere
278*/
279template<class T>
280uint_4 SphereHEALPix<T>::NbThetaSlices() const
281{
282 uint_4 nbslices = uint_4(4*nSide_-1);
283 if (nSide_<=0)
284 {
285 nbslices = 0;
286 throw PException(" sphere not pixelized, NbSlice=0 ");
287 }
288 return nbslices;
289}
290
291/*! \fn void SOPHYA::SphereHEALPix::GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const
292
293 For a theta-slice with index 'index', return :
294
295 the corresponding "theta"
296
297 a vector containing the phi's of the pixels of the slice
298
299 a vector containing the corresponding values of pixels
300*/
301template<class T>
302void SphereHEALPix<T>::GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const
303
304{
305
306 if (index<0 || index >= NbThetaSlices())
307 {
308 cout << " SphereHEALPix::GetThetaSlice : Pixel index out of range" <<endl;
309 throw RangeCheckError(" SphereHEALPix::GetThetaSlice : Pixel index out of range");
310 }
311
312
313 int_4 iring= sliceBeginIndex_(index);
314 int_4 lring = sliceLenght_(index);
315
316 phi.ReSize(lring);
317 value.ReSize(lring);
318
319 double TH= 0.;
320 double FI= 0.;
321 for(int_4 kk = 0; kk < lring;kk++)
322 {
323 PixThetaPhi(kk+iring,TH,FI);
324 phi(kk)= FI;
325 value(kk)= PixVal(kk+iring);
326 }
327 theta= TH;
328}
329/*! \fn void SOPHYA::SphereHEALPix::GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const
330
331 For a theta-slice with index 'index', return :
332
333 the corresponding "theta"
334
335 the corresponding "phi" for first pixel of the slice
336
337 a vector containing indices of the pixels of the slice
338
339 (equally distributed in phi)
340
341 a vector containing the corresponding values of pixels
342*/
343
344template<class T>
345void SphereHEALPix<T>::GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const
346
347{
348
349 if (sliceIndex<0 || sliceIndex >= NbThetaSlices())
350 {
351 cout << " SphereHEALPix::GetThetaSlice : Pixel index out of range" <<endl;
352 throw RangeCheckError(" SphereHEALPix::GetThetaSlice : Pixel index out of range");
353 }
354 int_4 iring= sliceBeginIndex_(sliceIndex);
355 int_4 lring = sliceLenght_(sliceIndex);
356 pixelIndices.ReSize(lring);
357 value.ReSize(lring);
358
359 for(int_4 kk = 0; kk < lring;kk++)
360 {
361 pixelIndices(kk)= kk+iring;
362 value(kk)= PixVal(kk+iring);
363 }
364 PixThetaPhi(iring, theta, phi0);
365}
366
367template<class T>
368void SphereHEALPix<T>::SetThetaSlices()
369{
370 sliceBeginIndex_.ReSize(4*nSide_-1);
371 sliceLenght_.ReSize(4*nSide_-1);
372 int sliceIndex;
373 for (sliceIndex=0; sliceIndex< nSide_-1; sliceIndex++)
374 {
375 sliceBeginIndex_(sliceIndex) = 2*sliceIndex*(sliceIndex+1);
376 sliceLenght_(sliceIndex) = 4*(sliceIndex+1);
377 }
378 for (sliceIndex= nSide_-1; sliceIndex< 3*nSide_; sliceIndex++)
379 {
380 sliceBeginIndex_(sliceIndex) = 2*nSide_*(2*sliceIndex-nSide_+1);
381 sliceLenght_(sliceIndex) = 4*nSide_;
382 }
383 for (sliceIndex= 3*nSide_; sliceIndex< 4*nSide_-1; sliceIndex++)
384 {
385 int_4 nc= 4*nSide_-1-sliceIndex;
386 sliceBeginIndex_(sliceIndex) = nPix_-2*nc*(nc+1);
387 sliceLenght_(sliceIndex) = 4*nc;
388 }
389}
390
391
392/*! \fn T& SOPHYA::SphereHEALPix::PixVal(int_4 k)
393
394 \return value of pixel with "RING" index k
395*/
396template<class T>
397T& SphereHEALPix<T>::PixVal(int_4 k)
398
399{
400 if((k < 0) || (k >= nPix_))
401 {
402 throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range");
403 }
404 return pixels_(k);
405}
406
407/*! \fn T const& SOPHYA::SphereHEALPix::PixVal(int_4 k) const
408
409 \return value of pixel with "RING" index k
410*/
411template<class T>
412T const& SphereHEALPix<T>::PixVal(int_4 k) const
413
414{
415 if((k < 0) || (k >= nPix_))
416 {
417 throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range");
418 }
419 return *(pixels_.Data()+k);
420}
421
422/*! \fn T& SOPHYA::SphereHEALPix::PixValNest(int_4 k)
423
424 \return value of pixel with "NESTED" index k
425*/
426template<class T>
427T& SphereHEALPix<T>::PixValNest(int_4 k)
428
429//--
430{
431 if((k < 0) || (k >= nPix_))
432 {
433 throw RangeCheckError("SphereHEALPix::PIxValNest Pixel index out of range");
434 }
435 return pixels_(nest2ring(nSide_,k));
436}
437
438/*! \fn T const& SOPHYA::SphereHEALPix::PixValNest(int_4 k) const
439
440 \return value of pixel with "NESTED" index k
441*/
442
443template<class T>
444T const& SphereHEALPix<T>::PixValNest(int_4 k) const
445
446{
447 if((k < 0) || (k >= nPix_))
448 {
449 throw RangeCheckError("SphereHEALPix::PIxValNest Pixel index out of range");
450 }
451 int_4 pix= nest2ring(nSide_,k);
452 return *(pixels_.Data()+pix);
453}
454
455/*! \fn bool SOPHYA::SphereHEALPix::ContainsSph(double theta, double phi) const
456
457\return true if teta,phi in map
458*/
459template<class T>
460bool SphereHEALPix<T>::ContainsSph(double theta, double phi) const
461{
462return(true);
463}
464
465/*! \fn int_4 SOPHYA::SphereHEALPix::PixIndexSph(double theta,double phi) const
466
467 \return "RING" index of the pixel corresponding to direction (theta, phi).
468 */
469template<class T>
470int_4 SphereHEALPix<T>::PixIndexSph(double theta,double phi) const
471
472{
473 return ang2pix_ring(nSide_,theta,phi);
474}
475
476/*! \fn int_4 SOPHYA::SphereHEALPix::PixIndexSphNest(double theta,double phi) const
477
478 \return "NESTED" index of the pixel corresponding to direction (theta, phi).
479 */
480template<class T>
481int_4 SphereHEALPix<T>::PixIndexSphNest(double theta,double phi) const
482
483{
484 return ang2pix_nest(nSide_,theta,phi);
485}
486
487
488/* --Methode-- */
489/*! \fn void SOPHYA::SphereHEALPix::PixThetaPhi(int_4 k,double& theta,double& phi) const
490 \return (theta,phi) coordinates of middle of pixel with "RING" index k
491 */
492template<class T>
493void SphereHEALPix<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
494{
495 pix2ang_ring(nSide_,k,theta,phi);
496}
497
498/*! \fn T SOPHYA::SphereHEALPix::SetPixels(T v)
499
500Set all pixels to value v
501*/
502template <class T>
503T SphereHEALPix<T>::SetPixels(T v)
504{
505pixels_.Reset(v);
506return(v);
507}
508
509
510
511/*! \fn void SOPHYA::SphereHEALPix::PixThetaPhiNest(int_4 k,double& theta,double& phi) const
512
513 \return (theta,phi) coordinates of middle of pixel with "NESTED" index k
514 */
515template<class T>
516void SphereHEALPix<T>::PixThetaPhiNest(int_4 k,double& theta,double& phi) const
517
518{
519 pix2ang_nest(nSide_,k,theta,phi);
520}
521
522
523/*! \fn int_4 SOPHYA::SphereHEALPix::NestToRing(int_4 k) const
524
525 translation from NESTED index into RING index
526*/
527template<class T>
528int_4 SphereHEALPix<T>::NestToRing(int_4 k) const
529
530{
531 return nest2ring(nSide_,k);
532}
533
534
535/*! \fn int_4 SOPHYA::SphereHEALPix::RingToNest(int_4 k) const
536
537 translation from RING index into NESTED index
538*/
539template<class T>
540int_4 SphereHEALPix<T>::RingToNest(int_4 k) const
541{
542 return ring2nest(nSide_,k);
543}
544
545
546template <class T>
547void SphereHEALPix<T>::print(ostream& os) const
548{
549 if(mInfo_) os << " DVList Info= " << *mInfo_ << endl;
550 //
551 os << " nSide_ = " << nSide_ << endl;
552 os << " nPix_ = " << nPix_ << endl;
553 os << " omeg_ = " << omeg_ << endl;
554
555 os << " content of pixels : ";
556 for(int i=0; i < nPix_; i++)
557 {
558 if(i%5 == 0) os << endl;
559 os << pixels_(i) <<", ";
560 }
561 os << endl;
562
563
564}
565
566
567
568//*******************************************************************
569
570#ifdef __CXX_PRAGMA_TEMPLATES__
571#pragma define_template SphereHEALPix<uint_2>
572#pragma define_template SphereHEALPix<r_8>
573#pragma define_template SphereHEALPix<r_4>
574#pragma define_template SphereHEALPix< complex<r_4> >
575#pragma define_template SphereHEALPix< complex<r_8> >
576#endif
577#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
578template class SphereHEALPix<uint_2>;
579template class SphereHEALPix<r_8>;
580template class SphereHEALPix<r_4>;
581template class SphereHEALPix< complex<r_4> >;
582template class SphereHEALPix< complex<r_8> >;
583#endif
584
Note: See TracBrowser for help on using the repository browser.