| 1 | #ifndef SPHEREHEALPIX_SEEN
|
|---|
| 2 | #define SPHEREHEALPIX_SEEN
|
|---|
| 3 |
|
|---|
| 4 | #include "sphericalmap.h"
|
|---|
| 5 | #include "tvector.h"
|
|---|
| 6 | #include "ndatablock.h"
|
|---|
| 7 |
|
|---|
| 8 | #include "anydataobj.h"
|
|---|
| 9 | #include "ppersist.h"
|
|---|
| 10 |
|
|---|
| 11 | #include "HEALPixUtils.h"
|
|---|
| 12 |
|
|---|
| 13 | namespace SOPHYA {
|
|---|
| 14 |
|
|---|
| 15 |
|
|---|
| 16 | // ***************** CLASSE SphereHEALPix *****************************
|
|---|
| 17 |
|
|---|
| 18 | //! class SphereHEALPix
|
|---|
| 19 | /*!
|
|---|
| 20 | Pixelisation Gorski
|
|---|
| 21 |
|
|---|
| 22 |
|
|---|
| 23 | -----------------------------------------------------------------------
|
|---|
| 24 | version 0.8.2 Aug97 TAC Eric Hivon, Kris Gorski
|
|---|
| 25 | -----------------------------------------------------------------------
|
|---|
| 26 |
|
|---|
| 27 | the sphere is split in 12 diamond-faces containing nside**2 pixels each
|
|---|
| 28 |
|
|---|
| 29 | the numbering of the pixels (in the nested scheme) is similar to
|
|---|
| 30 | quad-cube
|
|---|
| 31 | In each face the first pixel is in the lowest corner of the diamond
|
|---|
| 32 |
|
|---|
| 33 | the faces are (x,y) coordinate on each face
|
|---|
| 34 | \verbatim
|
|---|
| 35 | . . . . <--- North Pole
|
|---|
| 36 | / \ / \ / \ / \ ^ ^
|
|---|
| 37 | . 0 . 1 . 2 . 3 . <--- z = 2/3 \ /
|
|---|
| 38 | \ / \ / \ / \ / y \ / x
|
|---|
| 39 | 4 . 5 . 6 . 7 . 4 <--- equator \ /
|
|---|
| 40 | / \ / \ / \ / \ \/
|
|---|
| 41 | . 8 . 9 .10 .11 . <--- z = -2/3 (0,0) : lowest corner
|
|---|
| 42 | \ / \ / \ / \ /
|
|---|
| 43 | . . . . <--- South Pole
|
|---|
| 44 | \endverbatim
|
|---|
| 45 | phi:0 2Pi
|
|---|
| 46 |
|
|---|
| 47 | in the ring scheme pixels are numbered along the parallels
|
|---|
| 48 | the first parallel is the one closest to the north pole and so on
|
|---|
| 49 | on each parallel, pixels are numbered starting from the one closest
|
|---|
| 50 | to phi = 0
|
|---|
| 51 |
|
|---|
| 52 | nside MUST be a power of 2 (<= 8192)
|
|---|
| 53 |
|
|---|
| 54 | */
|
|---|
| 55 |
|
|---|
| 56 |
|
|---|
| 57 | template<class T>
|
|---|
| 58 | class FIO_SphereHEALPix;
|
|---|
| 59 |
|
|---|
| 60 | template<class T>
|
|---|
| 61 | class FITS_SphereHEALPix;
|
|---|
| 62 |
|
|---|
| 63 | template<class T>
|
|---|
| 64 | class SphereHEALPix : public SphericalMap<T>
|
|---|
| 65 | {
|
|---|
| 66 | public :
|
|---|
| 67 | // Static Methods to ease the use of HEALPix index <> angle conversion methods
|
|---|
| 68 |
|
|---|
| 69 | static inline int_4 nest2ring(int_4 nside,int_4 ipnest)
|
|---|
| 70 | { return HEALPix::nest2ring(nside, ipnest); }
|
|---|
| 71 | static inline int_4 ring2nest(int_4 nside,int_4 ipring)
|
|---|
| 72 | { return HEALPix::ring2nest(nside, ipring); }
|
|---|
| 73 | static inline int_4 ang2pix_ring(int_4 nside,double theta,double phi)
|
|---|
| 74 | { return HEALPix::ang2pix_ring(nside, theta, phi); }
|
|---|
| 75 | static inline int_4 ang2pix_nest(int_4 nside,double theta,double phi)
|
|---|
| 76 | { return HEALPix::ang2pix_nest(nside, theta, phi); }
|
|---|
| 77 | static inline void pix2ang_ring(int_4 nside,int_4 ipix,double& theta,double& phi)
|
|---|
| 78 | { HEALPix::pix2ang_ring(nside, ipix, theta, phi); }
|
|---|
| 79 | static inline void pix2ang_nest(int_4 nside,int_4 ipix,double& theta,double& phi)
|
|---|
| 80 | { HEALPix::pix2ang_nest(nside, ipix, theta, phi); }
|
|---|
| 81 |
|
|---|
| 82 | SphereHEALPix();
|
|---|
| 83 | /*!
|
|---|
| 84 | m is the "nside" of the Gorski algorithm
|
|---|
| 85 |
|
|---|
| 86 | The total number of pixels will be Npix = 12*nside**2
|
|---|
| 87 |
|
|---|
| 88 | nside MUST be a power of 2 (<= 8192)
|
|---|
| 89 | */
|
|---|
| 90 | SphereHEALPix(int_4 m);
|
|---|
| 91 | SphereHEALPix(const SphereHEALPix<T>& s, bool share);
|
|---|
| 92 | SphereHEALPix(const SphereHEALPix<T>& s);
|
|---|
| 93 | //! Destructor
|
|---|
| 94 | virtual ~SphereHEALPix();
|
|---|
| 95 |
|
|---|
| 96 | // Temporaire?
|
|---|
| 97 | inline virtual bool IsTemp(void) const {
|
|---|
| 98 |
|
|---|
| 99 | if (sliceBeginIndex_.IsTemp() != pixels_.IsTemp() || sliceLenght_.IsTemp() != pixels_.IsTemp() )
|
|---|
| 100 | throw PException(" l'etat 'temporaire' de la spherehealpix est incoherent");
|
|---|
| 101 | return pixels_.IsTemp();
|
|---|
| 102 | }
|
|---|
| 103 | /*! Setting blockdata to temporary (see ndatablock documentation) */
|
|---|
| 104 | inline virtual void SetTemp(bool temp=false) const
|
|---|
| 105 | {
|
|---|
| 106 | pixels_.SetTemp(temp);
|
|---|
| 107 | sliceBeginIndex_.SetTemp(temp);
|
|---|
| 108 | sliceLenght_.SetTemp(temp);
|
|---|
| 109 | };
|
|---|
| 110 | // ------------------ Definition of PixelMap abstract methods
|
|---|
| 111 |
|
|---|
| 112 | /* Nombre de pixels du decoupage */
|
|---|
| 113 | /*! Return number of pixels of the splitting */
|
|---|
| 114 | virtual int_4 NbPixels() const;
|
|---|
| 115 |
|
|---|
| 116 | /* Valeur du contenu du pixel d'indice "RING" k */
|
|---|
| 117 | /*! Return value of pixel with "RING" index k */
|
|---|
| 118 | virtual T& PixVal(int_4 k);
|
|---|
| 119 | virtual T const& PixVal(int_4 k) const;
|
|---|
| 120 |
|
|---|
| 121 | /* Nombre de tranches en theta */
|
|---|
| 122 | /*! Return number of slices in theta direction on the sphere */
|
|---|
| 123 | uint_4 NbThetaSlices() const;
|
|---|
| 124 | /*! For a theta-slice with index 'index', return :
|
|---|
| 125 |
|
|---|
| 126 | the corresponding "theta"
|
|---|
| 127 |
|
|---|
| 128 | a vector containing the phi's of the pixels of the slice
|
|---|
| 129 |
|
|---|
| 130 | a vector containing the corresponding values of pixels
|
|---|
| 131 | */
|
|---|
| 132 | virtual void GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const;
|
|---|
| 133 | /*! For a theta-slice with index 'index', return :
|
|---|
| 134 |
|
|---|
| 135 | the corresponding "theta"
|
|---|
| 136 |
|
|---|
| 137 | the corresponding "phi" for first pixel of the slice
|
|---|
| 138 |
|
|---|
| 139 | a vector containing indices of the pixels of the slice
|
|---|
| 140 |
|
|---|
| 141 | (equally distributed in phi)
|
|---|
| 142 |
|
|---|
| 143 | a vector containing the corresponding values of pixels
|
|---|
| 144 | */
|
|---|
| 145 | virtual void GetThetaSlice(int_4 sliceIndex,r_8& theta, r_8& phi0, TVector<int_4>& pixelIndices,TVector<T>& value) const ;
|
|---|
| 146 |
|
|---|
| 147 | /* Return true if teta,phi in map */
|
|---|
| 148 | virtual bool ContainsSph(double theta, double phi) const;
|
|---|
| 149 | /* Indice "RING" du pixel vers lequel pointe une direction definie par
|
|---|
| 150 | ses coordonnees spheriques */
|
|---|
| 151 | /*! Return "RING" index of the pixel corresponding to direction (theta, phi).
|
|---|
| 152 | */
|
|---|
| 153 | virtual int_4 PixIndexSph(double theta,double phi) const;
|
|---|
| 154 |
|
|---|
| 155 | /* Coordonnees spheriques du milieu du pixel d'indice "RING" k */
|
|---|
| 156 | virtual void PixThetaPhi(int_4 k,double& theta,double& phi) const;
|
|---|
| 157 |
|
|---|
| 158 | /*! Set all pixels to value v */
|
|---|
| 159 | virtual T SetPixels(T v);
|
|---|
| 160 |
|
|---|
| 161 | /* Pixel Solid angle (steradians) */
|
|---|
| 162 | /*! Pixel Solid angle (steradians)
|
|---|
| 163 |
|
|---|
| 164 | All the pixels have the same solid angle. The dummy argument is
|
|---|
| 165 | for compatibility with eventual pixelizations which would not
|
|---|
| 166 | fulfil this requirement.
|
|---|
| 167 | */
|
|---|
| 168 | virtual double PixSolAngle(int_4 dummy=0) const;
|
|---|
| 169 |
|
|---|
| 170 | /* Acces to the DataBlock */
|
|---|
| 171 | inline NDataBlock<T>& DataBlock() {return pixels_;}
|
|---|
| 172 | inline const NDataBlock<T>& DataBlock() const {return pixels_;}
|
|---|
| 173 |
|
|---|
| 174 | // --------------- Specific methods
|
|---|
| 175 |
|
|---|
| 176 | /*!
|
|---|
| 177 | m is the "nside" of the Gorski algorithm
|
|---|
| 178 |
|
|---|
| 179 | The total number of pixels will be Npix = 12*nside**2
|
|---|
| 180 |
|
|---|
| 181 | nside MUST be a power of 2 (<= 8192)
|
|---|
| 182 | */
|
|---|
| 183 | virtual void Resize(int_4 m);
|
|---|
| 184 |
|
|---|
| 185 | // pour l'instant le tableau est ordonne selon RING, uniquement
|
|---|
| 186 | inline virtual char* TypeOfMap() const {return "RING";};
|
|---|
| 187 |
|
|---|
| 188 |
|
|---|
| 189 | /* Valeur du contenu du pixel d'indice "NEST" k */
|
|---|
| 190 | /*! Return value of pixel with "NESTED" index k */
|
|---|
| 191 | virtual T& PixValNest(int_4 k);
|
|---|
| 192 | /*! Return value of pixel with "NESTED" index k */
|
|---|
| 193 | virtual T const& PixValNest(int_4 k) const;
|
|---|
| 194 |
|
|---|
| 195 | /* Indice "NEST" du pixel vers lequel pointe une direction definie par
|
|---|
| 196 | ses coordonnees spheriques */
|
|---|
| 197 | /*! Return "NESTED" index of the pixel corresponding to direction (theta, phi).
|
|---|
| 198 | */
|
|---|
| 199 | virtual int_4 PixIndexSphNest(double theta,double phi) const;
|
|---|
| 200 |
|
|---|
| 201 | /* Coordonnees spheriques du milieu du pixel d'indice "NEST" k */
|
|---|
| 202 | /*! Return (theta,phi) coordinates of middle of pixel with "NESTED" index k
|
|---|
| 203 | */
|
|---|
| 204 | virtual void PixThetaPhiNest(int_4 k,double& theta,double& phi) const;
|
|---|
| 205 |
|
|---|
| 206 | /* algorithme de pixelisation */
|
|---|
| 207 | void Pixelize(int_4);
|
|---|
| 208 |
|
|---|
| 209 | /* convertit index nested en ring */
|
|---|
| 210 | /*! translation from NESTED index into RING index */
|
|---|
| 211 | int_4 NestToRing(int_4) const;
|
|---|
| 212 |
|
|---|
| 213 | /* convertit index ring en nested" */
|
|---|
| 214 | /*! translation from RING index into NESTED index */
|
|---|
| 215 | int_4 RingToNest(int_4) const;
|
|---|
| 216 |
|
|---|
| 217 |
|
|---|
| 218 | /* retourne la valeur du parametre Gorski */
|
|---|
| 219 | inline virtual int_4 SizeIndex() const {return(nSide_);}
|
|---|
| 220 |
|
|---|
| 221 | /* impression */
|
|---|
| 222 | void print(ostream& os) const;
|
|---|
| 223 |
|
|---|
| 224 |
|
|---|
| 225 |
|
|---|
| 226 | inline SphereHEALPix<T>& operator = (const SphereHEALPix<T>& a)
|
|---|
| 227 | {return Set(a);}
|
|---|
| 228 |
|
|---|
| 229 | void CloneOrShare(const SphereHEALPix<T>& a);
|
|---|
| 230 | SphereHEALPix<T>& Set(const SphereHEALPix<T>& a);
|
|---|
| 231 | SphereHEALPix<T>& CopyElt(const SphereHEALPix<T>& a);
|
|---|
| 232 |
|
|---|
| 233 |
|
|---|
| 234 | // friend declaration for classes which handle persistence and FITS IO
|
|---|
| 235 | friend class FIO_SphereHEALPix<T>;
|
|---|
| 236 | friend class FITS_SphereHEALPix<T>;
|
|---|
| 237 |
|
|---|
| 238 | protected :
|
|---|
| 239 |
|
|---|
| 240 | // ------------- méthodes internes ----------------------
|
|---|
| 241 | void InitNul();
|
|---|
| 242 | void SetThetaSlices();
|
|---|
| 243 |
|
|---|
| 244 | inline void setParameters(int_4 nside, int_4 nbpixels, double solangle)
|
|---|
| 245 | {
|
|---|
| 246 | nSide_= nside;
|
|---|
| 247 | nPix_= nbpixels;
|
|---|
| 248 | omeg_= solangle;
|
|---|
| 249 | }
|
|---|
| 250 |
|
|---|
| 251 | // ------------- variables internes -----------------------
|
|---|
| 252 |
|
|---|
| 253 | int_4 nSide_;
|
|---|
| 254 | int_4 nPix_;
|
|---|
| 255 | double omeg_;
|
|---|
| 256 |
|
|---|
| 257 | NDataBlock<T> pixels_;
|
|---|
| 258 | NDataBlock<int_4> sliceBeginIndex_; // Rationalisation Mac. D.Y.
|
|---|
| 259 | NDataBlock<int_4> sliceLenght_;
|
|---|
| 260 |
|
|---|
| 261 | };
|
|---|
| 262 |
|
|---|
| 263 |
|
|---|
| 264 |
|
|---|
| 265 | } // Fin du namespace
|
|---|
| 266 |
|
|---|
| 267 | #endif
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
