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

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

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