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

Last change on this file since 2882 was 2882, checked in by ansari, 20 years ago

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