- Timestamp:
- Jun 21, 2006, 10:38:26 AM (19 years ago)
- Location:
- trunk/SophyaLib/SkyMap
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/SkyMap/fiospherehealpix.cc
r2869 r2978 81 81 uint_8 itab[3]; 82 82 is.Get(itab, 3); 83 84 // Juin 2006 : NumVer 1->2 avec mise en place de HEALPix NESTED 85 // Si itab[2]==1 --> NESTED 86 bool fgring = true; 87 if (itab[2] == 1) fgring = false; 88 83 89 // Let's Read the SphereCoordSys object -- ATTENTIOn - $CHECK$ 84 90 FIO_SphereCoordSys fio_scs( dobj_->GetCoordSys()); … … 100 106 double Omega; 101 107 is.GetR8(Omega); 102 dobj_->setParameters(nSide,nPix, Omega );108 dobj_->setParameters(nSide,nPix, Omega, fgring); 103 109 104 110 // On lit les DataBlocks; … … 124 130 // On ecrit 3 uint_8 125 131 // 0 : Numero de version, 1 : Size index, 2 reserve a l 132 // Juin 2006 : NumVer 1->2 avec mise en place de HEALPix NESTED 133 // totalement back-compatible , on utilise itab[2] qu'on met a 1 134 // pour NESTED 126 135 uint_8 itab[3]; 127 itab[0] = 1;136 itab[0] = 2; 128 137 itab[1] = dobj_->SizeIndex(); 129 itab[2] = 0;138 itab[2] = (dobj_->IfRING()) ? 0 : 1; 130 139 os.Put(itab, 3); 131 140 // Let's write the SphereCoordSys object … … 164 173 165 174 #ifdef __CXX_PRAGMA_TEMPLATES__ 175 #pragma define_template FIO_SphereHEALPix<uint_2> 166 176 #pragma define_template FIO_SphereHEALPix<int_4> 167 177 #pragma define_template FIO_SphereHEALPix<r_8> … … 172 182 #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES) 173 183 namespace SOPHYA { 184 template class FIO_SphereHEALPix<uint_2>; 174 185 template class FIO_SphereHEALPix<int_4>; 175 186 template class FIO_SphereHEALPix<r_8>; -
trunk/SophyaLib/SkyMap/pixelmap.h
r2973 r2978 161 161 { 162 162 os << "PixelMap<T>::Show() class: " << typeid(*this).name() 163 << " T= " << DataTypeInfo<T>::getTypeName() << endl; 163 << " T= " << DataTypeInfo<T>::getTypeName() 164 << " TypeOfMap= " << TypeOfMap() << endl; 164 165 double solang = 0.; 165 166 if (NbPixels() > 0) solang = PixSolAngle(NbPixels()/2); -
trunk/SophyaLib/SkyMap/skymapinit.cc
r2615 r2978 54 54 55 55 // ---------- Les SphereHEALPix --------- 56 57 PPRegister(FIO_SphereHEALPix<uint_2>); 58 DObjRegister(FIO_SphereHEALPix<uint_2>, SphereHEALPix<uint_2>); 56 59 57 60 PPRegister(FIO_SphereHEALPix<int_4>); -
trunk/SophyaLib/SkyMap/sphereecp.h
r2973 r2978 25 25 //! return the size index value corresponding to resolution res (in radian) 26 26 static inline int_4 ResolToSizeIndex(double res) 27 { return ( M_PI/res); }27 { return (int_4)((M_PI/res)+0.5); } 28 28 29 29 // Constructeur par defaut -
trunk/SophyaLib/SkyMap/spherehealpix.cc
r2973 r2978 62 62 /* --Methode-- */ 63 63 64 template<class T> 65 SphereHEALPix<T>::SphereHEALPix() : pixels_(), sliceBeginIndex_(), 66 sliceLenght_() 67 64 //! Default constructor - optional pixelisation scheme parameter 65 template<class T> 66 SphereHEALPix<T>::SphereHEALPix(bool fgring) : fgring_(fgring), pixels_(), 67 sliceBeginIndex_(), sliceLenght_() 68 68 69 69 { … … 71 71 } 72 72 73 /*! \ fn SOPHYA::SphereHEALPix::SphereHEALPix(int_4 m)74 73 /*! \brief Constructor with specification of nside and optional pixelisation scheme 74 75 75 \param <m> : "nside" of the Healpix algorithm 76 \param <fgring> : if true -> RING pixelisation (default), if not NESTED 76 77 77 78 The total number of pixels will be Npix = 12*nside**2 … … 81 82 82 83 template<class T> 83 SphereHEALPix<T>::SphereHEALPix(int_4 m) 84 85 { 84 SphereHEALPix<T>::SphereHEALPix(int_4 m, bool fgring) 85 86 { 87 fgring_ = fgring; 86 88 InitNul(); 87 89 Pixelize(m); … … 98 100 nPix_ = s.nPix_; 99 101 omeg_ = s.omeg_; 102 fgring_ = s.fgring_; 100 103 if(s.mInfo_) this->mInfo_= new DVList(*s.mInfo_); 101 104 } … … 111 114 nPix_ = s.nPix_; 112 115 omeg_ = s.omeg_; 116 fgring_ = s.fgring_; 113 117 if(s.mInfo_) this->mInfo_= new DVList(*s.mInfo_); 114 118 // CloneOrShare(s); … … 123 127 nPix_ = a.nPix_; 124 128 omeg_ = a.omeg_; 129 fgring_ = a.fgring_; 125 130 pixels_.CloneOrShare(a.pixels_); 126 131 sliceBeginIndex_.CloneOrShare(a.sliceBeginIndex_); … … 137 142 nPix_ = a.nPix_; 138 143 omeg_ = a.omeg_; 144 fgring_ = a.fgring_; 139 145 pixels_.Share(a.pixels_); 140 146 sliceBeginIndex_.Share(a.sliceBeginIndex_); … … 148 154 SphereHEALPix<T>& SphereHEALPix<T>::Set(const SphereHEALPix<T>& a) 149 155 { 150 if (this != &a) 151 { 152 156 if (this != &a) { 153 157 if (a.NbPixels() < 1) 154 158 throw RangeCheckError("SphereHEALPix<T>::Set(a ) - SphereHEALPix a not allocated ! "); … … 156 160 else CopyElt(a); 157 161 158 159 162 if (this->mInfo_) delete this->mInfo_; 160 163 this->mInfo_ = NULL; 161 164 if (a.mInfo_) this->mInfo_ = new DVList(*(a.mInfo_)); 162 165 } 163 166 return(*this); 164 167 } … … 168 171 { 169 172 if (NbPixels() < 1) 170 throw RangeCheckError("SphereHEALPix<T>::CopyElt( const SphereHEALPix<T>&) - Not Allocated SphereHEALPix ! ");173 throw RangeCheckError("SphereHEALPix<T>::CopyElt(a) - Not Allocated SphereHEALPix ! "); 171 174 if (NbPixels() != a.NbPixels()) 172 throw(SzMismatchError("SphereHEALPix<T>::CopyElt( const SphereHEALPix<T>&) SizeMismatch")) ;175 throw(SzMismatchError("SphereHEALPix<T>::CopyElt(a) SizeMismatch")) ; 173 176 nSide_= a.nSide_; 174 177 nPix_ = a.nPix_; 175 178 omeg_ = a.omeg_; 176 179 int k; 177 for (k=0; k< nPix_; k++) pixels_(k) = a.pixels_(k); 180 if (fgring_ != a.fgring_) 181 for (k=0; k< nPix_; k++) pixels_(k) = a.pixels_(k); 182 else { 183 if (fgring_) for (k=0; k< nPix_; k++) 184 pixels_(k) = a.pixels_(ring2nest(nSide_, k)); 185 else for (k=0; k< nPix_; k++) 186 pixels_(k) = a.pixels_(nest2ring(nSide_, k)); 187 } 178 188 for (k=0; k< a.sliceBeginIndex_.Size(); k++) sliceBeginIndex_(k) = a.sliceBeginIndex_(k); 179 189 for (k=0; k< a.sliceLenght_.Size(); k++) sliceLenght_(k) = a.sliceLenght_(k); … … 181 191 } 182 192 183 184 185 //++186 // Titre Destructor187 //--188 //++189 193 template<class T> 190 194 SphereHEALPix<T>::~SphereHEALPix() 191 192 //-- 193 { 194 } 195 196 //++ 197 // Titre Public Methods 198 //-- 195 { 196 } 199 197 200 198 /*! \fn void SOPHYA::SphereHEALPix::Resize(int_4 m) … … 207 205 template<class T> 208 206 void SphereHEALPix<T>::Resize(int_4 m) 209 210 211 207 { 212 208 if ((m <= 0 && nSide_ > 0)) { … … 219 215 } 220 216 217 218 //! return type of the map pixelisation : RING or NESTED 219 template<class T> 220 string SphereHEALPix<T>::TypeOfMap() const 221 { 222 if (fgring_) return string("RING"); 223 else return string("NESTED"); 224 } 225 221 226 template<class T> 222 227 void SphereHEALPix<T>::Pixelize( int_4 m) 223 224 228 // prépare la pixelisation Gorski (m a la même signification 225 229 // que pour le constructeur) 226 //227 //228 230 { 229 231 if (m<=0 || m> 8192) { … … 291 293 { 292 294 uint_4 nbslices = uint_4(4*nSide_-1); 293 if (nSide_<=0) 294 { 295 nbslices = 0; 296 throw PException(" sphere not pixelized, NbSlice=0 "); 297 } 295 if (nSide_<=0) { 296 nbslices = 0; 297 throw PException(" sphere not pixelized, NbSlice=0 "); 298 } 298 299 return nbslices; 299 300 } … … 338 339 template<class T> 339 340 void SphereHEALPix<T>::GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const 340 341 { 342 343 if (index<0 || index >= NbThetaSlices()) 344 { 341 { 342 if (index<0 || index >= NbThetaSlices()) { 345 343 cout << " SphereHEALPix::GetThetaSlice : Pixel index out of range" <<endl; 346 344 throw RangeCheckError(" SphereHEALPix::GetThetaSlice() index out of range"); 347 } 348 345 } 349 346 350 347 int_4 iring= sliceBeginIndex_(index); 351 348 int_4 lring = sliceLenght_(index); 352 349 353 354 350 phi.ReSize(lring); 351 value.ReSize(lring); 355 352 356 353 double TH= 0.; 357 354 double FI= 0.; 358 for(int_4 kk = 0; kk < lring;kk++)359 {355 if (fgring_) { // RING pixelisation scheme 356 for(int_4 kk = 0; kk < lring;kk++) { 360 357 PixThetaPhi(kk+iring,TH,FI); 361 358 phi(kk)= FI; 362 value(kk)= PixVal(kk+iring); 363 } 359 value(kk)= pixels_(kk+iring); 360 } 361 } 362 else { // NESTED pixelisation scheme 363 for(int_4 kk = 0; kk < lring;kk++) { 364 int kkn = ring2nest(nSide_, kk+iring); 365 PixThetaPhi(kkn,TH,FI); 366 phi(kk)= FI; 367 value(kk)= pixels_(kkn); 368 } 369 } 364 370 theta= TH; 365 371 } … … 394 400 value.ReSize(lring); 395 401 396 for(int_4 kk = 0; kk < lring;kk++)397 {402 if (fgring_) { // RING pixelisation scheme 403 for(int_4 kk = 0; kk < lring;kk++) { 398 404 pixelIndices(kk)= kk+iring; 399 value(kk)= PixVal(kk+iring); 400 } 401 PixThetaPhi(iring, theta, phi0); 405 value(kk)= pixels_(kk+iring); 406 } 407 PixThetaPhi(iring, theta, phi0); 408 } 409 else { // NESTED pixelisation scheme 410 for(int_4 kk = 0; kk < lring;kk++) { 411 int_4 kkn = ring2nest(nSide_, kk+iring); 412 pixelIndices(kk)= kkn; 413 value(kk)= pixels_(kkn); 414 } 415 PixThetaPhi(ring2nest(nSide_,iring), theta, phi0); 416 } 402 417 } 403 418 … … 427 442 428 443 429 /*! \fn T& SOPHYA::SphereHEALPix::PixVal(int_4 k) 430 431 \return value of pixel with "RING" index k 432 */ 444 445 //! \return value of the \b k th pixel 433 446 template<class T> 434 447 T& SphereHEALPix<T>::PixVal(int_4 k) … … 436 449 { 437 450 if((k < 0) || (k >= nPix_)) 438 { 439 throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range"); 440 } 451 throw RangeCheckError("SphereHEALPix::PixVal() Pixel index out of range"); 441 452 return pixels_(k); 442 453 } 443 454 444 /*! \fn T const& SOPHYA::SphereHEALPix::PixVal(int_4 k) const 445 446 \return value of pixel with "RING" index k 447 */ 455 //! \return value of the \b k th pixel 448 456 template<class T> 449 457 T const& SphereHEALPix<T>::PixVal(int_4 k) const … … 451 459 { 452 460 if((k < 0) || (k >= nPix_)) 453 { 454 throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range"); 455 } 461 throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range"); 456 462 return *(pixels_.Data()+k); 457 463 } … … 460 466 461 467 \return value of pixel with "NESTED" index k 462 */ 468 463 469 template<class T> 464 470 T& SphereHEALPix<T>::PixValNest(int_4 k) … … 472 478 return pixels_(nest2ring(nSide_,k)); 473 479 } 474 480 */ 475 481 /*! \fn T const& SOPHYA::SphereHEALPix::PixValNest(int_4 k) const 476 482 477 483 \return value of pixel with "NESTED" index k 478 */479 484 480 485 template<class T> … … 489 494 return *(pixels_.Data()+pix); 490 495 } 496 */ 491 497 492 498 /*! \fn bool SOPHYA::SphereHEALPix::ContainsSph(double theta, double phi) const … … 508 514 509 515 { 510 return ang2pix_ring(nSide_,theta,phi); 516 if (fgring_) return ang2pix_ring(nSide_,theta,phi); 517 else return ang2pix_nest(nSide_,theta,phi); 511 518 } 512 519 … … 514 521 515 522 \return "NESTED" index of the pixel corresponding to direction (theta, phi). 516 */ 523 517 524 template<class T> 518 525 int_4 SphereHEALPix<T>::PixIndexSphNest(double theta,double phi) const … … 521 528 return ang2pix_nest(nSide_,theta,phi); 522 529 } 523 524 525 /* --Methode-- */ 526 /*! \fn void SOPHYA::SphereHEALPix::PixThetaPhi(int_4 k,double& theta,double& phi) const 527 \return (theta,phi) coordinates of middle of pixel with "RING" index k 528 */ 530 */ 531 532 //! \return (theta,phi) coordinates of middle of pixel with "RING" index k 529 533 template<class T> 530 534 void SphereHEALPix<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const 531 535 { 532 pix2ang_ring(nSide_,k,theta,phi); 533 } 534 535 /*! \fn T SOPHYA::SphereHEALPix::SetPixels(T v) 536 537 Set all pixels to value v 538 */ 536 if (fgring_) pix2ang_ring(nSide_,k,theta,phi); 537 else pix2ang_nest(nSide_,k,theta,phi); 538 } 539 540 //! Set all pixels to value v 539 541 template <class T> 540 542 T SphereHEALPix<T>::SetPixels(T v) … … 549 551 550 552 \return (theta,phi) coordinates of middle of pixel with "NESTED" index k 551 */552 553 template<class T> 553 554 void SphereHEALPix<T>::PixThetaPhiNest(int_4 k,double& theta,double& phi) const … … 556 557 pix2ang_nest(nSide_,k,theta,phi); 557 558 } 558 559 560 /*! \fn int_4 SOPHYA::SphereHEALPix::NestToRing(int_4 k) const 561 562 translation from NESTED index into RING index 563 */ 559 */ 560 561 //! Conversion from NESTED index into RING index 564 562 template<class T> 565 563 int_4 SphereHEALPix<T>::NestToRing(int_4 k) const 566 567 564 { 568 565 return nest2ring(nSide_,k); 569 566 } 570 567 571 572 /*! \fn int_4 SOPHYA::SphereHEALPix::RingToNest(int_4 k) const 573 574 translation from RING index into NESTED index 575 */ 568 //! Conversion from RING index into NESTED index 576 569 template<class T> 577 570 int_4 SphereHEALPix<T>::RingToNest(int_4 k) const … … 592 585 } 593 586 594 / *! Add a constant value \b x to a SphereHEALPix */587 //! Add a constant value \b x to a SphereHEALPix 595 588 template <class T> 596 589 SphereHEALPix<T>& SphereHEALPix<T>::Add(T a) … … 641 634 { 642 635 if (NbPixels() != a.NbPixels() ) 643 { 644 throw(SzMismatchError("SphereHEALPix<T>::AddElt(const SphereHEALPix<T>&) SizeMismatch")) ; 645 } 636 throw(SzMismatchError("SphereHEALPix<T>::AddElt(a) SizeMismatch")) ; 637 if (fgring_ != a.fgring_) 638 throw(ParmError("SphereHEALPix<T>::AddElt(a) different pixelisation RING<>NESTED")) ; 639 646 640 pixels_ += a.pixels_; 647 641 return (*this); … … 653 647 { 654 648 if (NbPixels() != a.NbPixels() ) 655 { 656 throw(SzMismatchError("SphereHEALPix<T>::SubElt(const SphereHEALPix<T>&) SizeMismatch")) ; 657 } 649 throw(SzMismatchError("SphereHEALPix<T>::SubElt(a) SizeMismatch")) ; 650 if (fgring_ != a.fgring_) 651 throw(ParmError("SphereHEALPix<T>::SubElt(a) different pixelisation RING<>NESTED")) ; 652 658 653 pixels_ -= a.pixels_; 659 654 return (*this); … … 665 660 { 666 661 if (NbPixels() != a.NbPixels() ) 667 { 668 throw(SzMismatchError("SphereHEALPix<T>::MulElt(const SphereHEALPix<T>&) SizeMismatch")) ; 669 } 662 throw(SzMismatchError("SphereHEALPix<T>::MulElt(a) SizeMismatch")) ; 663 if (fgring_ != a.fgring_) 664 throw(ParmError("SphereHEALPix<T>::MulElt(a) different pixelisation RING<>NESTED")) ; 665 670 666 pixels_ *= a.pixels_; 671 667 return (*this); … … 677 673 { 678 674 if (NbPixels() != a.NbPixels() ) 679 {680 throw(SzMismatchError("SphereHEALPix<T>::DivElt(const SphereHEALPix<T>&) SizeMismatch")) ;681 }675 throw(SzMismatchError("SphereHEALPix<T>::DivElt(a) SizeMismatch")) ; 676 if (fgring_ != a.fgring_) 677 throw(ParmError("SphereHEALPix<T>::DivElt(a) different pixelisation RING<>NESTED")) ; 682 678 pixels_ /= a.pixels_; 683 679 return (*this); … … 691 687 { 692 688 Show(os); 689 os << "SphereHEALPix<T>(" << TypeOfMap() << ") NSide= " 690 << nSide_ << " nPix_ = " << nPix_ << " omeg_ = " << omeg_ << endl; 691 693 692 if(this->mInfo_) os << " DVList Info= " << *(this->mInfo_) << endl; 694 // 695 os << " nSide_ = " << nSide_ << endl; 696 os << " nPix_ = " << nPix_ << endl; 697 os << " omeg_ = " << omeg_ << endl; 698 699 os << " content of pixels : "; 693 os << "... Pixel values : "; 700 694 for(int i=0; i < nPix_; i++) 701 695 { -
trunk/SophyaLib/SkyMap/spherehealpix.h
r2973 r2978 52 52 { return HEALPix::ResolToSizeIndex(res); } 53 53 54 SphereHEALPix( );55 SphereHEALPix(int_4 m );54 SphereHEALPix(bool fgring=true); 55 SphereHEALPix(int_4 m, bool fgring=true); 56 56 SphereHEALPix(const SphereHEALPix<T>& s, bool share); 57 57 SphereHEALPix(const SphereHEALPix<T>& s); … … 93 93 virtual T SetPixels(T v); 94 94 95 /*! Pixel Solid angle (steradians)95 /*! \brief Pixel Solid angle (steradians) 96 96 97 97 All the pixels have the same solid angle. The dummy argument is … … 108 108 109 109 virtual void Resize(int_4 m); 110 111 /*! 112 113 \return type of storage of the map : RING or NESTED 114 115 at the moment, always RING 116 */ 117 inline virtual string TypeOfMap() const {return string("RING");}; 118 119 110 virtual string TypeOfMap() const; 111 112 inline bool IfRING() const { return fgring_; } 113 inline bool IfNESTED() const { return ( (fgring_) ? false : true ); } 114 115 /* 120 116 virtual T& PixValNest(int_4 k); 121 117 virtual T const& PixValNest(int_4 k) const; 122 123 118 virtual int_4 PixIndexSphNest(double theta,double phi) const; 124 125 119 virtual void PixThetaPhiNest(int_4 k,double& theta,double& phi) const; 120 */ 126 121 127 122 void Pixelize(int_4); … … 138 133 inline void Print(ostream& os) const { print(os); } 139 134 140 141 142 143 144 145 146 // Operations diverses = , +=, ... 147 135 //--------------------- Operations diverses = , +=, ... 148 136 149 137 SphereHEALPix<T>& Set(const SphereHEALPix<T>& a); … … 153 141 // A += -= *= /= x (ajoute, soustrait, ... x a tous les elements) 154 142 155 143 //! Fill SphereHEALPix with all elements equal to \b x 156 144 virtual SphereHEALPix<T>& SetT(T a); 157 145 inline SphereHEALPix<T>& operator = (T a) {return SetT(a);} … … 202 190 void SetThetaSlices(); 203 191 204 inline void setParameters(int_4 nside, int_4 nbpixels, double solangle) 205 { 206 nSide_= nside; 207 nPix_= nbpixels; 208 omeg_= solangle; 209 } 192 inline void setParameters(int_4 nside, int_4 nbpixels, double solangle, 193 bool fgring) 194 { 195 nSide_= nside; 196 nPix_= nbpixels; 197 omeg_= solangle; 198 fgring_ = fgring; 199 } 210 200 211 201 // ------------- variables internes ----------------------- … … 214 204 int_4 nPix_; 215 205 double omeg_; 206 bool fgring_; // true -> RING pixelisation , false -> NESTED 216 207 217 208 NDataBlock<T> pixels_; 218 NDataBlock<int_4> sliceBeginIndex_; // Rationalisation Mac. D.Y.209 NDataBlock<int_4> sliceBeginIndex_; // Indices in RING scheme 219 210 NDataBlock<int_4> sliceLenght_; 220 211 -
trunk/SophyaLib/SkyMap/spherethetaphi.h
r2973 r2978 29 29 //! return the size index value corresponding to resolution res (in radian) 30 30 static inline int_4 ResolToSizeIndex(double res) 31 { return ( M_PI/2./res); }31 { return (int_4)((M_PI/2./res)+0.5); } 32 32 33 33 SphereThetaPhi();
Note:
See TracChangeset
for help on using the changeset viewer.