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