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

Last change on this file since 2621 was 2621, checked in by ansari, 21 years ago

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