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

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

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