Line data Source code
1 : //# StokesVector.h: A fast casacore::RigidVector with casacore::Stokes conversions
2 : //# Copyright (C) 1996,1999,2001,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be adressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 :
29 : #ifndef MSVIS_STOKESVECTOR_H
30 : #define MSVIS_STOKESVECTOR_H
31 :
32 : #include <casacore/casa/aips.h>
33 : #include <casacore/casa/IO/AipsIO.h>
34 : //#include <casacore/tables/Tables/TableRecord.h>
35 : #include <casacore/casa/BasicSL/Complex.h>
36 : #include <casacore/casa/Arrays/Vector.h>
37 : #include <casacore/casa/Arrays/Matrix.h>
38 : #include <casacore/casa/Arrays/MatrixMath.h>
39 : #include <casacore/scimath/Mathematics/RigidVector.h>
40 : #include <casacore/scimath/Mathematics/SquareMatrix.h>
41 : #include <casacore/casa/Arrays/IPosition.h>
42 : #include <casacore/casa/BasicMath/Math.h>
43 : #include <casacore/casa/iostream.h>
44 :
45 : namespace casa { //# NAMESPACE CASA - BEGIN
46 :
47 : //# Forward Declarations
48 : class StokesVector;
49 : // <summary>
50 : // Two specialized 4-vector classes for polarization handling
51 : // </summary>
52 :
53 : // <use visibility=export>
54 :
55 : // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
56 : // </reviewed>
57 : // <prerequisite>
58 : // <li> RigidVector
59 : // <li> Stokes
60 : // </prerequisite>
61 : //
62 : // <etymology>
63 : // StokesVector and CStokesVector (casacore::Complex StokesVector) are two classes
64 : // designed to handle a 4-vector of polarization values (I,Q,U,V or
65 : // e.g., RR,RL,LR,LL).
66 : // </etymology>
67 : //
68 : // <synopsis>
69 : // StokesVectors are RigidVectors of length 4. They have a few special
70 : // member functions to do polarization conversions efficiently.
71 : // CStokesVector also has a real() member function to get at the real
72 : // part of the components.
73 : // </synopsis>
74 : //
75 : // <example>
76 : // <srcblock>
77 : // // Create a real valued I,Q,U,V StokesVector
78 : // StokesVector pixel(4.0,2.0,1.0,0.1);
79 : // // convert to a casacore::Complex valued vector of linear polarizations
80 : // CStokesVector cohVec=applySlin(pixel);
81 : // // convert back to I,Q,U,V
82 : // cohVec.applySlinInv();
83 : // // Write out the real part
84 : // cout<< cohVec.real() <<endl;
85 : // </srcblock>
86 : // </example>
87 : //
88 : // <motivation>
89 : // Full polarization processing of interferometry data uses real and complex
90 : // valued 4-vectors. The StokesVector specialization handles this more
91 : // efficiently than a standard casacore::Vector of size 4.
92 : // </motivation>
93 : //
94 : // <thrown>
95 : // <li> No exceptions
96 : // </thrown>
97 : //
98 : // <todo asof="1996/11/07">
99 : // <li> we may want to add Stokesvector Eigenvalues
100 : // </todo>
101 :
102 :
103 : class CStokesVector:public casacore::RigidVector<casacore::Complex,4> {
104 : public:
105 :
106 : static casacore::String dataTypeId() {return "CStokesVector";};
107 : // CStokesVector(casacore::Int n):RigidVector<casacore::Complex,4>(n) {}
108 : // The casacore::Complex data members are automatically initialized to 0.
109 0 : CStokesVector():RigidVector<casacore::Complex,4>() {}
110 : // Construct from scalar, setting all values to a constant.
111 0 : explicit CStokesVector(const casacore::Complex& c):RigidVector<casacore::Complex,4>(c) {}
112 : // Construct with four values specified.
113 3660260 : CStokesVector(const casacore::Complex& v0, const casacore::Complex& v1,
114 3660260 : const casacore::Complex& v2, const casacore::Complex& v3):
115 3660260 : casacore::RigidVector<casacore::Complex,4>(v0,v1,v2,v3) {}
116 : // Construct from c-array
117 0 : CStokesVector(const casacore::Complex v[4]):RigidVector<casacore::Complex,4>(v) {}
118 : // Construct from casacore::Vector (should have length 4)
119 : // CStokesVector(const casacore::Vector<casacore::Complex> & v):RigidVector<casacore::Complex,4>(v) {}
120 : // Copy constructor with copy semantics.
121 0 : CStokesVector(const CStokesVector& v):RigidVector<casacore::Complex,4>(v){}
122 : // Construct from RigidVector
123 : // CStokesVector(const casacore::RigidVector<casacore::Complex,4>& v):RigidVector<casacore::Complex,4>(v) {}
124 : // Assignment
125 3660000 : CStokesVector& operator=(const CStokesVector& v) {
126 3660000 : casacore::RigidVector<casacore::Complex,4>::operator=(v); return *this;
127 : }
128 : // Assign from a Vector
129 0 : CStokesVector& operator=(const casacore::Vector<casacore::Complex>& v) {
130 0 : casacore::RigidVector<casacore::Complex,4>::operator=(v); return *this;
131 : }
132 : // Assign from a scalar, setting all values to a constant
133 : CStokesVector& operator=(const casacore::Complex& c) {
134 : casacore::RigidVector<casacore::Complex,4>::operator=(c); return *this;
135 : }
136 : // Negation
137 : CStokesVector& operator-() {
138 : casacore::RigidVector<casacore::Complex,4>::operator-(); return *this;
139 : }
140 : // Addition
141 0 : CStokesVector& operator+=(const CStokesVector& v) {
142 0 : casacore::RigidVector<casacore::Complex,4>::operator+=(v); return *this;
143 : }
144 : // Subtraction
145 0 : CStokesVector& operator-=(const CStokesVector& v) {
146 0 : casacore::RigidVector<casacore::Complex,4>::operator-=(v); return *this;
147 : }
148 : CStokesVector& operator*=(const CStokesVector& v) {
149 : casacore::RigidVector<casacore::Complex,4>::operator*=(v); return *this;
150 : }
151 : // casacore::Matrix multiplication - v*=m is equivalent to v=m*v
152 0 : CStokesVector& operator*=(const casacore::SquareMatrix<casacore::Complex,4>& m) {
153 0 : casacore::RigidVector<casacore::Complex,4>::operator*=(m); return *this;
154 : }
155 0 : CStokesVector& operator*=(casacore::Float f) {
156 0 : v_p[0]*=f; v_p[1]*=f; v_p[2]*=f; v_p[3]*=f; return *this;
157 : }
158 : // Equality
159 : casacore::Bool operator==(const CStokesVector& v) const {
160 : return (v_p[0]==v.v_p[0] && v_p[1]==v.v_p[1] &&
161 : v_p[2]==v.v_p[2] && v_p[3]==v.v_p[3]);
162 : }
163 : // Inequality
164 : casacore::Bool operator!=(const CStokesVector& v) const {
165 : return (v_p[0]!=v.v_p[0] || v_p[1]!=v.v_p[1] ||
166 : v_p[2]!=v.v_p[2] || v_p[3]!=v.v_p[3]);
167 : }
168 :
169 : // Apply conversion matrix from casacore::Stokes to linear, in place.
170 : CStokesVector& applySlin() {
171 : casacore::Complex i=v_p[0],q=v_p[1], u=v_p[2],iv=v_p[3]*casacore::Complex(0,1);
172 : v_p[0]=(i+q); v_p[1]=(u+iv);
173 : v_p[2]=(u-iv); v_p[3]=(i-q);
174 : return *this;
175 : }
176 : // Apply conversion matrix from casacore::Stokes to circular, in place
177 : CStokesVector& applyScirc() {
178 : casacore::Complex i=v_p[0],q=v_p[1],iu=v_p[2]*casacore::Complex(0,1),v=v_p[3];
179 : v_p[0]=(i+v); v_p[1]=(q+iu);
180 : v_p[2]=(q-iu); v_p[3]=(i-v);
181 : return *this;
182 : }
183 : // Apply conversion matrix from linear to casacore::Stokes, in place
184 : CStokesVector& applySlinInv() {
185 : using namespace casacore;
186 : casacore::Complex xx=v_p[0],xy=v_p[1],yx=v_p[2],yy=v_p[3];
187 : v_p[0]=(xx+yy)/2; v_p[1]=(xx-yy)/2;
188 : v_p[2]=(xy+yx)/2; v_p[3]=casacore::Complex(0,1)*(yx-xy)/2;
189 : return *this;
190 : }
191 : // Apply conversion matrix from circular to casacore::Stokes, in place
192 : CStokesVector& applyScircInv() {
193 : using namespace casacore;
194 : casacore::Complex rr=v_p[0],rl=v_p[1],lr=v_p[2],ll=v_p[3];
195 : v_p[0]=(rr+ll)/2; v_p[3]=(rr-ll)/2;
196 : v_p[1]=(rl+lr)/2; v_p[2]=casacore::Complex(0,1)*(lr-rl)/2;
197 : return *this;
198 : }
199 : // Return a StokesVector with the real part.
200 : // StokesVector real();
201 : // inner product of two complex vectors
202 : friend casacore::Complex innerProduct(const CStokesVector& l,
203 : const CStokesVector& r) {
204 : return l.v_p[0]*conj(r.v_p[0])+ l.v_p[1]*conj(r.v_p[1])+
205 : l.v_p[2]*conj(r.v_p[2])+ l.v_p[3]*conj(r.v_p[3]);
206 : }
207 : friend double norm(const CStokesVector& l) {
208 : using namespace casacore;
209 : return ::sqrt(square(l.v_p[0].real())+square(l.v_p[0].imag())+
210 : square(l.v_p[1].real())+square(l.v_p[1].imag())+
211 : square(l.v_p[2].real())+square(l.v_p[2].imag())+
212 : square(l.v_p[3].real())+square(l.v_p[3].imag()));
213 : }
214 : // Write out a CStokesVector using the casacore::Vector output method.
215 : friend std::ostream& operator<<(std::ostream& os, const CStokesVector& v) {
216 : os << v.vector();
217 : return os;
218 : }
219 : };
220 :
221 : // Multiplication of CStokesVector by a casacore::Complex SquareMatrix
222 : inline CStokesVector operator*(const casacore::SquareMatrix<casacore::Complex,4>& m,
223 : const CStokesVector& v) {
224 : CStokesVector result(v);
225 : return result*=m;
226 : }
227 :
228 : inline void defaultValue(CStokesVector& v) {
229 : v=casacore::Complex(0.0,0.0);
230 : }
231 :
232 : class StokesVector:public casacore::RigidVector<casacore::Float,4> {
233 :
234 : public:
235 : static casacore::String dataTypeId() {return "StokesVector";};
236 : // StokesVector(casacore::Int n):RigidVector<casacore::Float,4>(n) {}
237 : // Default constructor zeroes vector.
238 0 : StokesVector():RigidVector<casacore::Float,4>() {}
239 : // Construct from scalar, setting all values to a constant.
240 : StokesVector(casacore::Float f):RigidVector<casacore::Float,4>(f) {};
241 : // Construct with four values specified.
242 1210260 : StokesVector(casacore::Float v0, casacore::Float v1, casacore::Float v2, casacore::Float v3): casacore::RigidVector<casacore::Float,4>(v0,v1,v2,v3){}
243 : // Construct from casacore::Vector (should have length 4)
244 : // StokesVector(const casacore::Vector<casacore::Float> & v):RigidVector<casacore::Float,4>(v) {}
245 : // Copy constructor with copy semantics.
246 0 : StokesVector(const StokesVector& v):RigidVector<casacore::Float,4>(v) {}
247 : // Construct from RigidVector
248 : // StokesVector(const casacore::RigidVector<casacore::Float,4>& v):RigidVector<casacore::Float,4>(v) {}
249 : // Assignment
250 1210000 : StokesVector& operator=(const StokesVector& v) {
251 1210000 : casacore::RigidVector<casacore::Float,4>::operator=(v); return *this;
252 : }
253 : // Assign from a Vector
254 : StokesVector& operator=(const casacore::Vector<casacore::Float>& v) {
255 : casacore::RigidVector<casacore::Float,4>::operator=(v); return *this;
256 : }
257 : // Assign from a scalar, setting all values to a constant
258 : StokesVector& operator=(casacore::Float f) {
259 : casacore::RigidVector<casacore::Float,4>::operator=(f); return *this;
260 : }
261 : // Negation
262 : StokesVector& operator-() {
263 : casacore::RigidVector<casacore::Float,4>::operator-(); return *this;
264 : }
265 : // Addition
266 : StokesVector& operator+=(const StokesVector& v) {
267 : casacore::RigidVector<casacore::Float,4>::operator+=(v); return *this;
268 : }
269 : // Subtraction
270 : StokesVector& operator-=(const StokesVector& v) {
271 : casacore::RigidVector<casacore::Float,4>::operator-=(v); return *this;
272 : }
273 0 : StokesVector& operator*=(casacore::Float f) {
274 0 : casacore::RigidVector<casacore::Float,4>::operator*=(f); return *this;
275 : }
276 : StokesVector& operator*=(const StokesVector& v) {
277 : casacore::RigidVector<casacore::Float,4>::operator*=(v); return *this;
278 : }
279 : // casacore::Matrix multiplication - v*=m is equivalent to v=m*v
280 : StokesVector& operator*=(const casacore::SquareMatrix<casacore::Float,4>& m) {
281 : casacore::RigidVector<casacore::Float,4>::operator*=(m); return *this;
282 : }
283 : // Equality
284 : casacore::Bool operator==(const StokesVector& v) const {
285 : return (v_p[0]==v.v_p[0] && v_p[1]==v.v_p[1] &&
286 : v_p[2]==v.v_p[2] && v_p[3]==v.v_p[3]);
287 : }
288 : // Inequality
289 : casacore::Bool operator!=(const StokesVector& v) const {
290 : return (v_p[0]!=v.v_p[0] || v_p[1]!=v.v_p[1] ||
291 : v_p[2]!=v.v_p[2] || v_p[3]!=v.v_p[3]);
292 : }
293 : // Compute the maximum EigenValue
294 : casacore::Float maxEigenValue() const;
295 : // Compute the minimum EigenValue
296 : casacore::Float minEigenValue() const;
297 : // Compute the determinant of the coherence matrix
298 : casacore::Float determinant() const;
299 :
300 : // The innerproduct of 2 StokesVectors
301 0 : friend casacore::Float innerProduct(const StokesVector& l, const StokesVector& r) {
302 0 : return l.v_p[0]*r.v_p[0]+ l.v_p[1]*r.v_p[1]+
303 0 : l.v_p[2]*r.v_p[2]+ l.v_p[3]*r.v_p[3];
304 : }
305 : // Write out a StokesVector using the casacore::Vector output method.
306 0 : friend std::ostream& operator<<(std::ostream& os, const StokesVector& v) {
307 0 : os << v.vector();
308 0 : return os;
309 : }
310 :
311 : };
312 :
313 : inline void defaultValue(StokesVector& v) {
314 : v=0.0f;
315 : }
316 :
317 : // Multiply by a scalar
318 : inline StokesVector operator*(casacore::Float f, const StokesVector& v) {
319 : StokesVector r(v);
320 : return r*=f;
321 : }
322 : // Multiply by a scalar
323 0 : inline StokesVector operator*(const StokesVector& v, casacore::Float f) {
324 0 : StokesVector r(v);
325 0 : return r*=f;
326 : }
327 :
328 : // Multiplication of StokesVector by a SquareMatrix
329 : inline StokesVector operator*(const casacore::SquareMatrix<casacore::Float,4>& m,
330 : const StokesVector& v) {
331 : StokesVector result(v);
332 : return result*=m;
333 : }
334 :
335 : // Apply conversion matrix from casacore::Stokes to linear(returns result)
336 490000 : inline CStokesVector& applySlin(CStokesVector& result,
337 : const StokesVector& v) {
338 490000 : casacore::Complex t=casacore::Complex(0.,v(3));
339 490000 : result(0)=v(0)+v(1);
340 490000 : result(1)=v(2)+t;
341 490000 : result(2)=v(2)-t;
342 490000 : result(3)=v(0)-v(1);
343 490000 : return result;
344 : }
345 : // Apply conversion matrix from casacore::Stokes to linear.
346 : inline CStokesVector applySlin(const StokesVector& v) {
347 : CStokesVector result;
348 : return applySlin(result,v);
349 : }
350 : // Apply conversion matrix from casacore::Stokes to circular.
351 720000 : inline CStokesVector& applyScirc(CStokesVector& result,
352 : const StokesVector& v) {
353 720000 : casacore::Complex t=casacore::Complex(0.,1.0)*v(2);
354 720000 : result(0)=v(0)+v(3);
355 720000 : result(1)=v(1)+t;
356 720000 : result(2)=v(1)-t;
357 720000 : result(3)=v(0)-v(3);
358 720000 : return result;
359 : }
360 : // Apply conversion matrix from casacore::Stokes to circular.
361 : inline CStokesVector applyScirc(const StokesVector& v) {
362 : CStokesVector result;
363 : return applyScirc(result,v);
364 : }
365 :
366 : // Apply conversion matrix from linear to Stokes.
367 1290000 : inline StokesVector& applySlinInv(StokesVector& result, const CStokesVector& v) {
368 : using namespace casacore;
369 1290000 : result(0)=real(v(0)+v(3))/2;
370 1290000 : result(1)=real(v(0)-v(3))/2;
371 1290000 : result(2)=real(v(1)+v(2))/2;
372 1290000 : result(3)=real(casacore::Complex(0.,1.0)*(v(2)-v(1))/2);
373 1290000 : return result;
374 : }
375 :
376 : // Apply conversion matrix from linear to Stokes.
377 : inline StokesVector applySlinInv(const CStokesVector& v) {
378 : StokesVector result;
379 : return applySlinInv(result,v);
380 : }
381 :
382 : // Apply conversion matrix from circular to Stokes.
383 2370000 : inline StokesVector& applyScircInv(StokesVector& result, const CStokesVector& v) {
384 : using namespace casacore;
385 2370000 : result(0)=real(v(0)+v(3))/2;
386 2370000 : result(1)=real(v(1)+v(2))/2;
387 2370000 : result(2)=real(casacore::Complex(0.,1.0)*(v(2)-v(1))/2);
388 2370000 : result(3)=real(v(0)-v(3))/2;
389 2370000 : return result;
390 : }
391 :
392 : // Apply conversion matrix from circular to Stokes.
393 : inline StokesVector applyScircInv(const CStokesVector& v) {
394 : StokesVector result;
395 : return applyScircInv(result,v);
396 : }
397 :
398 : // The following are needed until Image no longer has
399 : // sigma images
400 : //StokesVector& sqrt(const StokesVector& v);
401 :
402 : //CStokesVector& sqrt(const CStokesVector& v);
403 :
404 :
405 : } //# NAMESPACE CASA - END
406 :
407 : #endif
|
