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

Last change on this file since 3890 was 3836, checked in by ansari, 15 years ago

Declaration extern des classes templ des classes template avec instantiation explicite (avec flag NEED_EXT_DECL_TEMP) pour regler le pb de dynamic_cast sur Mac OS X 10.6, Reza 05/08/2010
File size: 13.5 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 cout << " SphereECP<T>::Resize(int_4 " << m << ") m<0 ou m=NTheta - Ne fait rien " << endl;
316 return;
317 }
318 int mphi = m;
319 if (_pixels.Size() > 0) mphi = m*_pixels.SizeX()/_pixels.SizeY();
320 cout << " SphereECP<T>::Resize(" << m
321 << ") -> _pixels.Resize(" << mphi << "," << m << ")" << endl;
322 sa_size_t sz[5] = {0,0,0,0,0};
323 sz[0] = mphi; sz[1] = m;
324 _pixels.ReSize(2,sz);
325 _outofmapidx = _pixels.Size()+659;
326 _dtheta = (_theta2-_theta1)/m;
327 _dphi = (_phi2-_phi1)/mphi;
328
329 _outofmapntet = M_PI/_dtheta;
330 _outofmapnphi = 2*M_PI/_dphi;
331 return;
332}
333
334template <class T>
335uint_4 SphereECP<T>::NbThetaSlices() const
336{
337 return _pixels.SizeY();
338}
339
340template <class T>
341r_8 SphereECP<T>::ThetaOfSlice(int_4 index) const
342{
343 if( (index < 0) || (index >= _pixels.SizeY()) )
344 throw RangeCheckError("SphereECP::ThetaOfSlice() theta index out of range");
345 return (_theta1 + (index+0.5)*_dtheta);
346
347}
348
349template <class T>
350void SphereECP<T>::GetThetaSlice(int_4 index,r_8& theta,
351 TVector<r_8>& phi, TVector<T>& value) const
352{
353 if( (index < 0) || (index >= _pixels.SizeY()) )
354 throw RangeCheckError("SphereECP::GetThetaSlice() index out of range");
355
356 theta = _theta1 + (index+0.5)*_dtheta;
357 int_4 nphit = 2.*M_PI/_dphi;
358 phi.ReSize(nphit);
359 value.ReSize(nphit);
360 int_4 ioff = _phi1/_dphi;
361 int_4 i;
362 for(i=0; i<ioff; i++) {
363 phi(i) = _phi1 + (i-ioff+0.5)*_dphi;
364 value(i) = _outofmapval;
365 }
366 for(i=ioff; i<ioff+_pixels.SizeX(); i++) {
367 phi(i) = _phi1 + (i-ioff+0.5)*_dphi;
368 value(i) = _pixels(i-ioff,index,0);
369 }
370 for(i=ioff+_pixels.SizeX(); i<nphit; i++) {
371 phi(i) = _phi1 + (i-ioff+0.5)*_dphi;
372 value(i) = _outofmapval;
373 }
374 return;
375}
376
377template <class T>
378void SphereECP<T>::GetThetaSlice(int_4 index, r_8& theta, r_8& phi0,
379 TVector<int_4>& pixidx, TVector<T>& value) const
380{
381 if( (index < 0) || (index >= _pixels.SizeY()) )
382 throw RangeCheckError("SphereECP::GetThetaSlice() index out of range");
383
384 theta = _theta1 + (index+0.5)*_dtheta;
385 int_4 ioff = _phi1/_dphi;
386 phi0 = _phi1 - (ioff-0.5)*_dphi;
387
388 int_4 nphit = 2.*M_PI/_dphi;
389 pixidx.ReSize(nphit);
390 value.ReSize(nphit);
391
392 int_4 i;
393 for(i=0; i<ioff; i++) {
394 pixidx(i) = _outofmapidx+theta/_dtheta*_outofmapnphi+phi0/_dphi+i;
395 value(i) = _outofmapval;
396 }
397 for(i=ioff; i<ioff+_pixels.SizeX(); i++) {
398 pixidx(i) = index*_pixels.SizeX()+(i-ioff);
399 value(i) = _pixels(i-ioff,index,0);
400 }
401 for(i=ioff+_pixels.SizeX(); i<nphit; i++) {
402 pixidx(i) = _outofmapidx+theta/_dtheta*_outofmapnphi+phi0/_dphi+i;
403 value(i) = _outofmapval;
404 }
405 return;
406
407}
408
409template <class T>
410T* SphereECP<T>::GetThetaSliceDataPtr(int_4 index)
411{
412 if( (index < 0) || (index >= _pixels.SizeY()) )
413 throw RangeCheckError("SphereECP::GetThetaSliceDataPtr() index out of range");
414 int_4 ioff = _phi1/_dphi;
415 int_4 nphit = 2.*M_PI/_dphi;
416// Correction bug signale par cmv le 27/08/2008, debordement memoire lors
417// de synthese par SphericalTransformServer ...
418 if ( (ioff != 0) || (nphit!=_pixels.SizeX()) ) return NULL;
419 return _pixels.Data() + index*_pixels.SizeX();
420}
421
422template <class T>
423void SphereECP<T>::Show(ostream& os) const
424{
425 PixelMap<T>::Show(os);
426 if (_partial)
427 os << "SphereECP<T>::Print() - Partial ECP Map NPhi=" << _pixels.SizeX()
428 << " x NTheta= " << _pixels.SizeY() << " SolidAngle=" << PixSolAngle() << "\n"
429 << "ThetaLimits= " << _theta1 << " .. " << _theta2
430 << " PhiLimits= " << _phi1 << " .. " << _phi2 << endl;
431 else
432 os << "SphereECP<T>::Print() - Full ECP Map NPhi=" << _pixels.SizeX()
433 << " x NTheta= " << _pixels.SizeY() << " SolidAngle=" << PixSolAngle() << endl;
434}
435
436template <class T>
437void SphereECP<T>::Print(ostream& os) const
438{
439 Show(os);
440 os << _pixels;
441}
442
443template <class T>
444SphereECP<T>& SphereECP<T>::Set(const SphereECP<T>& a)
445{
446 _pixels.Set(a._pixels);
447 _partial = a._partial;
448 _theta1 = a._theta1;
449 _theta2 = a._theta2;
450 _phi1 = a._phi1;
451 _phi2 = a._phi2;
452 _outofmapidx = a._outofmapidx;
453 _outofmapnphi = a._outofmapnphi;
454 _outofmapntet = a._outofmapntet;
455 _outofmappix = a._outofmappix;
456 _outofmapval = a._outofmapval;
457 _dtheta = a._dtheta;
458 _dphi = a._dphi;
459 return *this ;
460}
461
462template <class T>
463SphereECP<T>& SphereECP<T>::SetCst(T x)
464{
465 _pixels.SetCst(x);
466 return *this ;
467}
468
469template <class T>
470SphereECP<T>& SphereECP<T>::AddCst(T x)
471{
472 _pixels += x;
473 return *this ;
474}
475
476template <class T>
477SphereECP<T>& SphereECP<T>::MulCst(T x)
478{
479 _pixels *= x;
480 return *this ;
481}
482
483
484
485}// Fin du namespace
486
487
488#ifdef __CXX_PRAGMA_TEMPLATES__
489#pragma define_template SphereECP<int_4>
490#pragma define_template SphereECP<r_4>
491#pragma define_template SphereECP<r_8>
492#pragma define_template SphereECP< complex<r_4> >
493#pragma define_template SphereECP< complex<r_8> >
494#endif
495#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
496namespace SOPHYA {
497template class SphereECP<int_4>;
498template class SphereECP<r_4>;
499template class SphereECP<r_8>;
500template class SphereECP< complex<r_4> >;
501template class SphereECP< complex<r_8> >;
502}
503#endif
504
Note: See TracBrowser for help on using the repository browser.