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

Last change on this file since 2618 was 2615, checked in by cmv, 21 years ago

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