source: Sophya/trunk/SophyaLib/SkyMap/sphereecp.cc@ 4066

Last change on this file since 4066 was 4066, checked in by ansari, 13 years ago

Ajout gestion de niveau d'impression globale pour les cartes du module SkyMap ds pixelmap.h et prise en compte du niveau pour les prints ds Resize des SphericalMaps (fait en janv 2012 ?), Reza 27/04/2012
File size: 13.6 KB
Line 
1#include "sopnamsp.h"
2#include <math.h>
3
4#define SPHEREECP_CC_BFILE // avoid extern template declarations
5#include "sphereecp.h"
6
7/*!
8 \class SOPHYA::SphereECP
9 \ingroup SkyMap
10 \brief Spherical maps in Equatorial Cylindrical Projection
11
12 Class implementing spherical maps, with a pixelisation
13 corresponding to the an Equatorial Cylindrical Projection (ECP).
14 Pixel size varie from equator to poles. The sphere is divided
15 into a number of slices along the parallels. Each slice is subdivided
16 into the same number of pixels.
17 This class provides the possibility to have partial coverage.
18 Can be used for the following types: <br>
19 <tt> int_4 , r_4 , r_8 , complex<r_4> , complex<r_8> </tt>
20*/
21
22namespace SOPHYA {
23
24//! Default constructor
25template <class T>
26SphereECP<T>::SphereECP()
27 : _pixels()
28{
29 _partial = false;
30 _theta1 = 0.;
31 _theta2 = M_PI;
32 _phi1 = 0.;
33 _phi2 = 2*M_PI;
34 _outofmapidx = -99;
35 _outofmapnphi = 0;
36 _outofmapntet = 0;
37 _outofmappix = 0;
38 _outofmapval = 0;
39 _dtheta = 0.;
40 _dphi = 0.;
41}
42
43/*!
44 \brief Constructor with the pixelisation parameter m
45
46 Complete coverage.
47 The sphere is divided into \b m slices in theta (0..pi). each slice is divided
48 into \b 2m pixels (along phi)
49*/
50template <class T>
51SphereECP<T>::SphereECP(int m)
52 : _pixels(2*m,m,0,0,0,true)
53{
54 _partial = false;
55 _theta1 = 0.;
56 _theta2 = M_PI;
57 _phi1 = 0.;
58 _phi2 = 2*M_PI;
59 _outofmapidx = -99;
60 _outofmapnphi = 0;
61 _outofmapntet = 0;
62 _outofmappix = 0;
63 _outofmapval = 0;
64 _dtheta = _dphi = M_PI/m;
65}
66
67/*!
68 Complete coverage.
69 \b ntet slices in theta (0..pi) and \b nphi pixels (along phi) for
70 each slice in theta
71*/
72template <class T>
73SphereECP<T>::SphereECP(int ntet, int nphi)
74 : _pixels(nphi,ntet,0,0,0,true)
75{
76 _partial = false;
77 _theta1 = 0.;
78 _theta2 = M_PI;
79 _phi1 = 0.;
80 _phi2 = 2*M_PI;
81 _outofmapidx = -99;
82 _outofmapnphi = 0;
83 _outofmapntet = 0;
84 _outofmappix = 0;
85 _outofmapval = 0;
86 _dtheta = M_PI/ntet;
87 _dphi = 2*M_PI/nphi;
88}
89
90/*!
91 Partial coverage :
92 \b ntet slices in theta , in the range (tet1 .. tet2) and \b nphi pixels (along phi) for
93 each slice in theta in the range (phi1 .. phi2)
94*/
95template <class T>
96SphereECP<T>::SphereECP(r_8 tet1, r_8 tet2, int ntet, r_8 phi1, r_8 phi2, int nphi)
97 : _pixels(nphi,ntet,0,0,0,true)
98{
99 if( (tet1 > M_PI) || (tet1 < 0.) || (tet2 > M_PI) || (tet2 < 0.) ||
100 (phi1 < 0.) || (phi1 > 2*M_PI) || (phi1 < 0.) || (phi1 > 2*M_PI) ||
101 (tet2 <= tet1) || (phi2 <= phi1) )
102 throw ParmError("SphereECP::SphereECP() Bad range for theta/phi limits");
103
104 _partial = true;
105 _theta1 = tet1;
106 _theta2 = tet2;
107 _phi1 = phi1;
108 _phi2 = phi2;
109 _outofmapidx = _pixels.Size()+659;
110 // Reza, 23 sep 2004
111 // pour des cartes partielles, et des pixels sur la sphere, mais
112 // en dehors de la couverture, on code le numero de tranche en theta,
113 // compte a partir de theta=0,phi=0 + _outofmapidx
114 // theta -> idx = _outofmapidx + theta/_dtheta*_outofmapnphi + _outofmapntet
115 _outofmapntet = 0;
116 _outofmappix = 0;
117 _outofmapval = 0;
118 _dtheta = (_theta2-_theta1)/ntet;
119 _dphi = (_phi2-_phi1)/nphi;
120
121 _outofmapntet = M_PI/_dtheta;
122 _outofmapnphi = 2*M_PI/_dphi;
123}
124
125//! Copy constructor - shares the pixel data if \b share=true
126template <class T>
127SphereECP<T>::SphereECP(const SphereECP<T>& a, bool share)
128 : _pixels(a._pixels, share)
129{
130 _partial = a._partial;
131 _theta1 = a._theta1;
132 _theta2 = a._theta2;
133 _phi1 = a._phi1;
134 _phi2 = a._phi2;
135 _outofmapidx = a._outofmapidx;
136 _outofmapnphi = a._outofmapnphi;
137 _outofmapntet = a._outofmapntet;
138 _outofmappix = a._outofmappix;
139 _outofmapval = a._outofmapval;
140 _dtheta = a._dtheta;
141 _dphi = a._dphi;
142}
143
144//! Copy constructor - shares the pixel data
145template <class T>
146SphereECP<T>::SphereECP(const SphereECP<T>& a)
147 : _pixels(a._pixels)
148{
149 _partial = a._partial;
150 _theta1 = a._theta1;
151 _theta2 = a._theta2;
152 _phi1 = a._phi1;
153 _phi2 = a._phi2;
154 _outofmapidx = a._outofmapidx;
155 _outofmapnphi = a._outofmapnphi;
156 _outofmapntet = a._outofmapntet;
157 _outofmappix = a._outofmappix;
158 _outofmapval = a._outofmapval;
159 _dtheta = a._dtheta;
160 _dphi = a._dphi;
161}
162
163template <class T>
164SphereECP<T>::~SphereECP()
165{
166}
167
168/*!
169 \brief Extract a partial map from a full ECP skymap
170 Theta and Phi limits are adjusted automatically to be correspond to the source map pixel boundaries
171*/
172template <class T>
173SphereECP<T> SphereECP<T>::ExtractPartial(r_8 tet1, r_8 tet2, r_8 phi1, r_8 phi2)
174{
175 if (IsPartial())
176 throw ParmError("SphereECP::ExtractPartial() source map NOT a full sky map");
177 if( (tet1 > M_PI) || (tet1 < 0.) || (tet2 > M_PI) || (tet2 < 0.) ||
178 (phi1 < 0.) || (phi1 > 2*M_PI) || (phi1 < 0.) || (phi1 > 2*M_PI) ||
179 (tet2 <= tet1) || (phi2 <= phi1) )
180 throw ParmError("SphereECP::ExtractPartial() Bad range for theta/phi limits");
181
182// correction/calcul des limites en theta et nombre de rangees en theta
183 sa_size_t ntet = (sa_size_t)((tet2-tet1)/_dtheta + 0.5);
184 sa_size_t offt = (sa_size_t)(tet1/_dtheta + 1e-3);
185 tet1 = _dtheta*offt;
186 tet2 = tet1+ntet*_dtheta;
187// correction/calcul des limites en phi et nombre de colonnes en phi
188 sa_size_t nphi = (sa_size_t)((phi2-phi1)/_dphi + 0.5);
189 sa_size_t offp = (sa_size_t)(phi1/_dphi + 1e-3);
190 phi1 = _dphi*offp;
191 phi2 = phi1+nphi*_dphi;
192
193// On cree une carte partielle avec les bons parametres
194 SphereECP<T> pmap(tet1, tet2, ntet, phi1, phi2, nphi);
195// On recopie les valeurs de pixels de la zone correspondante
196 pmap.GetPixelArray() = _pixels.SubArray(Range(offp, 0, nphi, 0), Range(offt, 0, ntet, 0),
197 Range::first(), Range::first(), Range::first());
198 return pmap;
199}
200
201template <class T>
202string SphereECP<T>::TypeOfMap() const
203{
204 return string("ECP");
205}
206
207template <class T>
208void SphereECP<T>::SetTemp(bool temp) const
209{
210 _pixels.SetTemp(temp);
211}
212
213template <class T>
214int_4 SphereECP<T>::NbPixels() const
215{
216 return _pixels.Size();
217}
218
219template <class T>
220T& SphereECP<T>::PixVal(int_4 k)
221{
222 // On pourrait etre plus strict en demandant aussi k<_outofmapidx+_outofmapnphi*_outofmapntet (Reza 23/09/04)
223 if ((_partial) && (k>=_outofmapidx)) return _outofmappix;
224 if((k < 0) || (k >= _pixels.Size()) )
225 throw RangeCheckError("SphereECP::PIxVal() Pixel index out of range");
226 return _pixels.DataBlock()(k);
227}
228
229template <class T>
230T const& SphereECP<T>::PixVal(int_4 k) const
231{
232 // On pourrait etre plus strict en demandant aussi k<_outofmapidx+_outofmapnphi*_outofmapntet (Reza 23/09/04)
233 if ((_partial) && (k>=_outofmapidx)) return _outofmapval;
234 if((k < 0) || (k >= _pixels.Size()) )
235 throw RangeCheckError("SphereECP::PIxVal() Pixel index out of range");
236 // return _pixels.DataBlock()(k);
237 return _pixels.Data()[k];
238}
239
240template <class T>
241bool SphereECP<T>::ContainsSph(double theta, double phi) const
242{
243 if (_partial) {
244 if ( (theta >= _theta1) && (theta <= _theta2) &&
245 (phi >= _phi1) && (phi <= _phi2) ) return true;
246 else return false;
247 }
248 else return false;
249}
250
251template <class T>
252int_4 SphereECP<T>::PixIndexSph(double theta, double phi) const
253{
254 if((theta > M_PI) || (theta < 0.)) return(-1);
255 if((phi < 0.) || (phi > 2*M_PI)) return(-1);
256 if (_partial) {
257 if ( (theta < _theta1) || (theta > _theta2) ||
258 (phi < _phi1) || (phi > _phi2) ) return _outofmapidx+theta/_dtheta*_outofmapnphi+phi/_dphi;
259 }
260 int_4 it = (theta-_theta1)/_dtheta;
261 int_4 jp = (phi-_phi1)/_dphi;
262 return (it*_pixels.SizeX()+jp);
263}
264
265template <class T>
266void SphereECP<T>::PixThetaPhi(int_4 k, double& theta, double& phi) const
267{
268 if ((_partial) && (k>=_outofmapidx)) {
269 theta = (((k-_outofmapidx) / _outofmapnphi)+0.5)*_dtheta;
270 phi = (((k-_outofmapidx) % _outofmapnphi)+0.5)*_dphi;;
271 return;
272 }
273 if((k < 0) || (k >= _pixels.Size()) )
274 throw RangeCheckError("SphereECP::PixThetaPhi() Pixel index out of range");
275
276 int_4 it = k / _pixels.SizeX();
277 int_4 jp = k % _pixels.SizeX();
278 theta = _theta1+(it+0.5)*_dtheta;
279 phi = _phi1+(jp+0.5)*_dphi;
280 return;
281}
282
283template <class T>
284T SphereECP<T>::SetPixels(T v)
285{
286 _pixels = v;
287 return v;
288}
289
290template <class T>
291double SphereECP<T>::PixSolAngle(int_4 k) const
292{
293 if ((_partial) && (k>=_outofmapidx)) {
294 double theta = (((k-_outofmapidx) / _outofmapnphi)+0.5)*_dtheta;
295 return (_dtheta*_dphi*sin(theta));
296 }
297 if((k < 0) || (k >= _pixels.Size()) )
298 throw RangeCheckError("SphereECP::PixSolAngle(int_4 k) Pixel index out of range");
299
300 int_4 it = k / _pixels.SizeX();
301 double theta = _theta1+it*_dtheta;
302 return (_dtheta*_dphi*sin(theta));
303}
304
305template <class T>
306int_4 SphereECP<T>::SizeIndex() const
307{
308 return _pixels.SizeY();
309}
310
311template <class T>
312void SphereECP<T>::Resize(int_4 m)
313{
314 if ( (m <= 0) || ( m == _pixels.SizeY()) ) {
315 if (PixelMap_GetGlobalPrintLevel()>0)
316 cout << " SphereECP<T>::Resize(int_4 " << m << ") m<0 ou m=NTheta - Ne fait rien " << endl;
317 return;
318 }
319 int mphi = m;
320 if (_pixels.Size() > 0) mphi = m*_pixels.SizeX()/_pixels.SizeY();
321 if (PixelMap_GetGlobalPrintLevel()>0)
322 cout << " SphereECP<T>::Resize(" << m
323 << ") -> _pixels.Resize(" << mphi << "," << m << ")" << endl;
324 sa_size_t sz[5] = {0,0,0,0,0};
325 sz[0] = mphi; sz[1] = m;
326 _pixels.ReSize(2,sz);
327 _outofmapidx = _pixels.Size()+659;
328 _dtheta = (_theta2-_theta1)/m;
329 _dphi = (_phi2-_phi1)/mphi;
330
331 _outofmapntet = M_PI/_dtheta;
332 _outofmapnphi = 2*M_PI/_dphi;
333 return;
334}
335
336template <class T>
337uint_4 SphereECP<T>::NbThetaSlices() const
338{
339 return _pixels.SizeY();
340}
341
342template <class T>
343r_8 SphereECP<T>::ThetaOfSlice(int_4 index) const
344{
345 if( (index < 0) || (index >= _pixels.SizeY()) )
346 throw RangeCheckError("SphereECP::ThetaOfSlice() theta index out of range");
347 return (_theta1 + (index+0.5)*_dtheta);
348
349}
350
351template <class T>
352void SphereECP<T>::GetThetaSlice(int_4 index,r_8& theta,
353 TVector<r_8>& phi, TVector<T>& value) const
354{
355 if( (index < 0) || (index >= _pixels.SizeY()) )
356 throw RangeCheckError("SphereECP::GetThetaSlice() index out of range");
357
358 theta = _theta1 + (index+0.5)*_dtheta;
359 int_4 nphit = 2.*M_PI/_dphi;
360 phi.ReSize(nphit);
361 value.ReSize(nphit);
362 int_4 ioff = _phi1/_dphi;
363 int_4 i;
364 for(i=0; i<ioff; i++) {
365 phi(i) = _phi1 + (i-ioff+0.5)*_dphi;
366 value(i) = _outofmapval;
367 }
368 for(i=ioff; i<ioff+_pixels.SizeX(); i++) {
369 phi(i) = _phi1 + (i-ioff+0.5)*_dphi;
370 value(i) = _pixels(i-ioff,index,0);
371 }
372 for(i=ioff+_pixels.SizeX(); i<nphit; i++) {
373 phi(i) = _phi1 + (i-ioff+0.5)*_dphi;
374 value(i) = _outofmapval;
375 }
376 return;
377}
378
379template <class T>
380void SphereECP<T>::GetThetaSlice(int_4 index, r_8& theta, r_8& phi0,
381 TVector<int_4>& pixidx, TVector<T>& value) const
382{
383 if( (index < 0) || (index >= _pixels.SizeY()) )
384 throw RangeCheckError("SphereECP::GetThetaSlice() index out of range");
385
386 theta = _theta1 + (index+0.5)*_dtheta;
387 int_4 ioff = _phi1/_dphi;
388 phi0 = _phi1 - (ioff-0.5)*_dphi;
389
390 int_4 nphit = 2.*M_PI/_dphi;
391 pixidx.ReSize(nphit);
392 value.ReSize(nphit);
393
394 int_4 i;
395 for(i=0; i<ioff; i++) {
396 pixidx(i) = _outofmapidx+theta/_dtheta*_outofmapnphi+phi0/_dphi+i;
397 value(i) = _outofmapval;
398 }
399 for(i=ioff; i<ioff+_pixels.SizeX(); i++) {
400 pixidx(i) = index*_pixels.SizeX()+(i-ioff);
401 value(i) = _pixels(i-ioff,index,0);
402 }
403 for(i=ioff+_pixels.SizeX(); i<nphit; i++) {
404 pixidx(i) = _outofmapidx+theta/_dtheta*_outofmapnphi+phi0/_dphi+i;
405 value(i) = _outofmapval;
406 }
407 return;
408
409}
410
411template <class T>
412T* SphereECP<T>::GetThetaSliceDataPtr(int_4 index)
413{
414 if( (index < 0) || (index >= _pixels.SizeY()) )
415 throw RangeCheckError("SphereECP::GetThetaSliceDataPtr() index out of range");
416 int_4 ioff = _phi1/_dphi;
417 int_4 nphit = 2.*M_PI/_dphi;
418// Correction bug signale par cmv le 27/08/2008, debordement memoire lors
419// de synthese par SphericalTransformServer ...
420 if ( (ioff != 0) || (nphit!=_pixels.SizeX()) ) return NULL;
421 return _pixels.Data() + index*_pixels.SizeX();
422}
423
424template <class T>
425void SphereECP<T>::Show(ostream& os) const
426{
427 PixelMap<T>::Show(os);
428 if (_partial)
429 os << "SphereECP<T>::Print() - Partial ECP Map NPhi=" << _pixels.SizeX()
430 << " x NTheta= " << _pixels.SizeY() << " SolidAngle=" << PixSolAngle() << "\n"
431 << "ThetaLimits= " << _theta1 << " .. " << _theta2
432 << " PhiLimits= " << _phi1 << " .. " << _phi2 << endl;
433 else
434 os << "SphereECP<T>::Print() - Full ECP Map NPhi=" << _pixels.SizeX()
435 << " x NTheta= " << _pixels.SizeY() << " SolidAngle=" << PixSolAngle() << endl;
436}
437
438template <class T>
439void SphereECP<T>::Print(ostream& os) const
440{
441 Show(os);
442 os << _pixels;
443}
444
445template <class T>
446SphereECP<T>& SphereECP<T>::Set(const SphereECP<T>& a)
447{
448 _pixels.Set(a._pixels);
449 _partial = a._partial;
450 _theta1 = a._theta1;
451 _theta2 = a._theta2;
452 _phi1 = a._phi1;
453 _phi2 = a._phi2;
454 _outofmapidx = a._outofmapidx;
455 _outofmapnphi = a._outofmapnphi;
456 _outofmapntet = a._outofmapntet;
457 _outofmappix = a._outofmappix;
458 _outofmapval = a._outofmapval;
459 _dtheta = a._dtheta;
460 _dphi = a._dphi;
461 return *this ;
462}
463
464template <class T>
465SphereECP<T>& SphereECP<T>::SetCst(T x)
466{
467 _pixels.SetCst(x);
468 return *this ;
469}
470
471template <class T>
472SphereECP<T>& SphereECP<T>::AddCst(T x)
473{
474 _pixels += x;
475 return *this ;
476}
477
478template <class T>
479SphereECP<T>& SphereECP<T>::MulCst(T x)
480{
481 _pixels *= x;
482 return *this ;
483}
484
485
486
487}// Fin du namespace
488
489
490#ifdef __CXX_PRAGMA_TEMPLATES__
491#pragma define_template SphereECP<int_4>
492#pragma define_template SphereECP<r_4>
493#pragma define_template SphereECP<r_8>
494#pragma define_template SphereECP< complex<r_4> >
495#pragma define_template SphereECP< complex<r_8> >
496#endif
497#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
498namespace SOPHYA {
499template class SphereECP<int_4>;
500template class SphereECP<r_4>;
501template class SphereECP<r_8>;
502template class SphereECP< complex<r_4> >;
503template class SphereECP< complex<r_8> >;
504}
505#endif
506
Note: See TracBrowser for help on using the repository browser.