casa
$Rev:20696$
|
00001 //# StokesVector.h: A fast RigidVector with Stokes conversions 00002 //# Copyright (C) 1996,1999,2001,2003 00003 //# Associated Universities, Inc. Washington DC, USA. 00004 //# 00005 //# This library is free software; you can redistribute it and/or modify it 00006 //# under the terms of the GNU Library General Public License as published by 00007 //# the Free Software Foundation; either version 2 of the License, or (at your 00008 //# option) any later version. 00009 //# 00010 //# This library is distributed in the hope that it will be useful, but WITHOUT 00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public 00013 //# License for more details. 00014 //# 00015 //# You should have received a copy of the GNU Library General Public License 00016 //# along with this library; if not, write to the Free Software Foundation, 00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. 00018 //# 00019 //# Correspondence concerning AIPS++ should be adressed as follows: 00020 //# Internet email: aips2-request@nrao.edu. 00021 //# Postal address: AIPS++ Project Office 00022 //# National Radio Astronomy Observatory 00023 //# 520 Edgemont Road 00024 //# Charlottesville, VA 22903-2475 USA 00025 //# 00026 //# 00027 //# $Id$ 00028 00029 #ifndef MSVIS_STOKESVECTOR_H 00030 #define MSVIS_STOKESVECTOR_H 00031 00032 #include <casa/aips.h> 00033 #include <casa/IO/AipsIO.h> 00034 //#include <tables/Tables/TableRecord.h> 00035 #include <casa/BasicSL/Complex.h> 00036 #include <casa/Arrays/Vector.h> 00037 #include <casa/Arrays/Matrix.h> 00038 #include <casa/Arrays/MatrixMath.h> 00039 #include <scimath/Mathematics/RigidVector.h> 00040 #include <scimath/Mathematics/SquareMatrix.h> 00041 #include <casa/Arrays/IPosition.h> 00042 #include <casa/BasicMath/Math.h> 00043 #include <casa/iostream.h> 00044 00045 namespace casa { //# NAMESPACE CASA - BEGIN 00046 00047 //# Forward Declarations 00048 class StokesVector; 00049 // <summary> 00050 // Two specialized 4-vector classes for polarization handling 00051 // </summary> 00052 00053 // <use visibility=export> 00054 00055 // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos=""> 00056 // </reviewed> 00057 // <prerequisite> 00058 // <li> RigidVector 00059 // <li> Stokes 00060 // </prerequisite> 00061 // 00062 // <etymology> 00063 // StokesVector and CStokesVector (Complex StokesVector) are two classes 00064 // designed to handle a 4-vector of polarization values (I,Q,U,V or 00065 // e.g., RR,RL,LR,LL). 00066 // </etymology> 00067 // 00068 // <synopsis> 00069 // StokesVectors are RigidVectors of length 4. They have a few special 00070 // member functions to do polarization conversions efficiently. 00071 // CStokesVector also has a real() member function to get at the real 00072 // part of the components. 00073 // </synopsis> 00074 // 00075 // <example> 00076 // <srcblock> 00077 // // Create a real valued I,Q,U,V StokesVector 00078 // StokesVector pixel(4.0,2.0,1.0,0.1); 00079 // // convert to a Complex valued vector of linear polarizations 00080 // CStokesVector cohVec=applySlin(pixel); 00081 // // convert back to I,Q,U,V 00082 // cohVec.applySlinInv(); 00083 // // Write out the real part 00084 // cout<< cohVec.real() <<endl; 00085 // </srcblock> 00086 // </example> 00087 // 00088 // <motivation> 00089 // Full polarization processing of interferometry data uses real and complex 00090 // valued 4-vectors. The StokesVector specialization handles this more 00091 // efficiently than a standard Vector of size 4. 00092 // </motivation> 00093 // 00094 // <thrown> 00095 // <li> No exceptions 00096 // </thrown> 00097 // 00098 // <todo asof="1996/11/07"> 00099 // <li> we may want to add Stokesvector Eigenvalues 00100 // </todo> 00101 00102 00103 class CStokesVector:public RigidVector<Complex,4> { 00104 public: 00105 static String dataTypeId() {return "CStokesVector";}; 00106 // CStokesVector(Int n):RigidVector<Complex,4>(n) {} 00107 // The Complex data members are automatically initialized to 0. 00108 CStokesVector():RigidVector<Complex,4>() {} 00109 // Construct from scalar, setting all values to a constant. 00110 CStokesVector(const Complex& c):RigidVector<Complex,4>(c) {} 00111 // Construct with four values specified. 00112 CStokesVector(const Complex& v0, const Complex& v1, 00113 const Complex& v2, const Complex& v3): 00114 RigidVector<Complex,4>(v0,v1,v2,v3) {} 00115 // Construct from c-array 00116 CStokesVector(const Complex v[4]):RigidVector<Complex,4>(v) {} 00117 // Construct from Vector (should have length 4) 00118 // CStokesVector(const Vector<Complex> & v):RigidVector<Complex,4>(v) {} 00119 // Copy constructor with copy semantics. 00120 CStokesVector(const CStokesVector& v):RigidVector<Complex,4>(v){} 00121 // Construct from RigidVector 00122 // CStokesVector(const RigidVector<Complex,4>& v):RigidVector<Complex,4>(v) {} 00123 // Assignment 00124 CStokesVector& operator=(const CStokesVector& v) { 00125 RigidVector<Complex,4>::operator=(v); return *this; 00126 } 00127 // Assign from a Vector 00128 CStokesVector& operator=(const Vector<Complex>& v) { 00129 RigidVector<Complex,4>::operator=(v); return *this; 00130 } 00131 // Assign from a scalar, setting all values to a constant 00132 CStokesVector& operator=(const Complex& c) { 00133 RigidVector<Complex,4>::operator=(c); return *this; 00134 } 00135 // Negation 00136 CStokesVector& operator-() { 00137 RigidVector<Complex,4>::operator-(); return *this; 00138 } 00139 // Addition 00140 CStokesVector& operator+=(const CStokesVector& v) { 00141 RigidVector<Complex,4>::operator+=(v); return *this; 00142 } 00143 // Subtraction 00144 CStokesVector& operator-=(const CStokesVector& v) { 00145 RigidVector<Complex,4>::operator-=(v); return *this; 00146 } 00147 CStokesVector& operator*=(const CStokesVector& v) { 00148 RigidVector<Complex,4>::operator*=(v); return *this; 00149 } 00150 // Matrix multiplication - v*=m is equivalent to v=m*v 00151 CStokesVector& operator*=(const SquareMatrix<Complex,4>& m) { 00152 RigidVector<Complex,4>::operator*=(m); return *this; 00153 } 00154 CStokesVector& operator*=(Float f) { 00155 v_p[0]*=f; v_p[1]*=f; v_p[2]*=f; v_p[3]*=f; return *this; 00156 } 00157 // Equality 00158 Bool operator==(const CStokesVector& v) const { 00159 return (v_p[0]==v.v_p[0] && v_p[1]==v.v_p[1] && 00160 v_p[2]==v.v_p[2] && v_p[3]==v.v_p[3]); 00161 } 00162 // Inequality 00163 Bool operator!=(const CStokesVector& v) const { 00164 return (v_p[0]!=v.v_p[0] || v_p[1]!=v.v_p[1] || 00165 v_p[2]!=v.v_p[2] || v_p[3]!=v.v_p[3]); 00166 } 00167 00168 // Apply conversion matrix from Stokes to linear, in place. 00169 CStokesVector& applySlin() { 00170 Complex i=v_p[0],q=v_p[1], u=v_p[2],iv=v_p[3]*Complex(0,1); 00171 v_p[0]=(i+q); v_p[1]=(u+iv); 00172 v_p[2]=(u-iv); v_p[3]=(i-q); 00173 return *this; 00174 } 00175 // Apply conversion matrix from Stokes to circular, in place 00176 CStokesVector& applyScirc() { 00177 Complex i=v_p[0],q=v_p[1],iu=v_p[2]*Complex(0,1),v=v_p[3]; 00178 v_p[0]=(i+v); v_p[1]=(q+iu); 00179 v_p[2]=(q-iu); v_p[3]=(i-v); 00180 return *this; 00181 } 00182 // Apply conversion matrix from linear to Stokes, in place 00183 CStokesVector& applySlinInv() { 00184 Complex xx=v_p[0],xy=v_p[1],yx=v_p[2],yy=v_p[3]; 00185 v_p[0]=(xx+yy)/2; v_p[1]=(xx-yy)/2; 00186 v_p[2]=(xy+yx)/2; v_p[3]=Complex(0,1)*(yx-xy)/2; 00187 return *this; 00188 } 00189 // Apply conversion matrix from circular to Stokes, in place 00190 CStokesVector& applyScircInv() { 00191 Complex rr=v_p[0],rl=v_p[1],lr=v_p[2],ll=v_p[3]; 00192 v_p[0]=(rr+ll)/2; v_p[3]=(rr-ll)/2; 00193 v_p[1]=(rl+lr)/2; v_p[2]=Complex(0,1)*(lr-rl)/2; 00194 return *this; 00195 } 00196 // Return a StokesVector with the real part. 00197 // StokesVector real(); 00198 // inner product of two complex vectors 00199 friend Complex innerProduct(const CStokesVector& l, 00200 const CStokesVector& r) { 00201 return l.v_p[0]*conj(r.v_p[0])+ l.v_p[1]*conj(r.v_p[1])+ 00202 l.v_p[2]*conj(r.v_p[2])+ l.v_p[3]*conj(r.v_p[3]); 00203 } 00204 friend double norm(const CStokesVector& l) { 00205 return ::sqrt(square(l.v_p[0].real())+square(l.v_p[0].imag())+ 00206 square(l.v_p[1].real())+square(l.v_p[1].imag())+ 00207 square(l.v_p[2].real())+square(l.v_p[2].imag())+ 00208 square(l.v_p[3].real())+square(l.v_p[3].imag())); 00209 } 00210 // Write out a CStokesVector using the Vector output method. 00211 friend ostream& operator<<(ostream& os, const CStokesVector& v) { 00212 os << v.vector(); 00213 return os; 00214 } 00215 }; 00216 00217 // Multiplication of CStokesVector by a Complex SquareMatrix 00218 inline CStokesVector operator*(const SquareMatrix<Complex,4>& m, 00219 const CStokesVector& v) { 00220 CStokesVector result(v); 00221 return result*=m; 00222 } 00223 00224 inline void defaultValue(CStokesVector& v) { 00225 v=Complex(0.0,0.0); 00226 } 00227 00228 class StokesVector:public RigidVector<Float,4> { 00229 00230 public: 00231 static String dataTypeId() {return "StokesVector";}; 00232 // StokesVector(Int n):RigidVector<Float,4>(n) {} 00233 // Default constructor zeroes vector. 00234 StokesVector():RigidVector<Float,4>() {} 00235 // Construct from scalar, setting all values to a constant. 00236 StokesVector(Float f):RigidVector<Float,4>(f) {}; 00237 // Construct with four values specified. 00238 StokesVector(Float v0, Float v1, Float v2, Float v3): RigidVector<Float,4>(v0,v1,v2,v3){} 00239 // Construct from Vector (should have length 4) 00240 // StokesVector(const Vector<Float> & v):RigidVector<Float,4>(v) {} 00241 // Copy constructor with copy semantics. 00242 StokesVector(const StokesVector& v):RigidVector<Float,4>(v) {} 00243 // Construct from RigidVector 00244 // StokesVector(const RigidVector<Float,4>& v):RigidVector<Float,4>(v) {} 00245 // Assignment 00246 StokesVector& operator=(const StokesVector& v) { 00247 RigidVector<Float,4>::operator=(v); return *this; 00248 } 00249 // Assign from a Vector 00250 StokesVector& operator=(const Vector<Float>& v) { 00251 RigidVector<Float,4>::operator=(v); return *this; 00252 } 00253 // Assign from a scalar, setting all values to a constant 00254 StokesVector& operator=(Float f) { 00255 RigidVector<Float,4>::operator=(f); return *this; 00256 } 00257 // Negation 00258 StokesVector& operator-() { 00259 RigidVector<Float,4>::operator-(); return *this; 00260 } 00261 // Addition 00262 StokesVector& operator+=(const StokesVector& v) { 00263 RigidVector<Float,4>::operator+=(v); return *this; 00264 } 00265 // Subtraction 00266 StokesVector& operator-=(const StokesVector& v) { 00267 RigidVector<Float,4>::operator-=(v); return *this; 00268 } 00269 StokesVector& operator*=(Float f) { 00270 RigidVector<Float,4>::operator*=(f); return *this; 00271 } 00272 StokesVector& operator*=(const StokesVector& v) { 00273 RigidVector<Float,4>::operator*=(v); return *this; 00274 } 00275 // Matrix multiplication - v*=m is equivalent to v=m*v 00276 StokesVector& operator*=(const SquareMatrix<Float,4>& m) { 00277 RigidVector<Float,4>::operator*=(m); return *this; 00278 } 00279 // Equality 00280 Bool operator==(const StokesVector& v) const { 00281 return (v_p[0]==v.v_p[0] && v_p[1]==v.v_p[1] && 00282 v_p[2]==v.v_p[2] && v_p[3]==v.v_p[3]); 00283 } 00284 // Inequality 00285 Bool operator!=(const StokesVector& v) const { 00286 return (v_p[0]!=v.v_p[0] || v_p[1]!=v.v_p[1] || 00287 v_p[2]!=v.v_p[2] || v_p[3]!=v.v_p[3]); 00288 } 00289 // Compute the maximum EigenValue 00290 Float maxEigenValue() const; 00291 // Compute the minimum EigenValue 00292 Float minEigenValue() const; 00293 // Compute the determinant of the coherence matrix 00294 Float determinant() const; 00295 00296 // The innerproduct of 2 StokesVectors 00297 friend Float innerProduct(const StokesVector& l, const StokesVector& r) { 00298 return l.v_p[0]*r.v_p[0]+ l.v_p[1]*r.v_p[1]+ 00299 l.v_p[2]*r.v_p[2]+ l.v_p[3]*r.v_p[3]; 00300 } 00301 // Multiplication of StokesVector by a Complex SquareMatrix 00302 friend CStokesVector operator*(const SquareMatrix<Complex,4>& m, 00303 const StokesVector& v) { 00304 // return m * (RigidVector<Float,4>&)v; 00305 return m * v; 00306 } 00307 // Write out a StokesVector using the Vector output method. 00308 friend ostream& operator<<(ostream& os, const StokesVector& v) { 00309 os << v.vector(); 00310 return os; 00311 } 00312 00313 }; 00314 00315 inline void defaultValue(StokesVector& v) { 00316 v=0.0f; 00317 } 00318 00319 // Multiply by a scalar 00320 inline StokesVector operator*(Float f, const StokesVector& v) { 00321 StokesVector r(v); 00322 return r*=f; 00323 } 00324 // Multiply by a scalar 00325 inline StokesVector operator*(const StokesVector& v, Float f) { 00326 StokesVector r(v); 00327 return r*=f; 00328 } 00329 00330 // Multiplication of StokesVector by a SquareMatrix 00331 inline StokesVector operator*(const SquareMatrix<Float,4>& m, 00332 const StokesVector& v) { 00333 StokesVector result(v); 00334 return result*=m; 00335 } 00336 00337 // Apply conversion matrix from Stokes to linear(returns result) 00338 inline CStokesVector& applySlin(CStokesVector& result, 00339 const StokesVector& v) { 00340 Complex t=Complex(0.,v(3)); 00341 result(0)=v(0)+v(1); 00342 result(1)=v(2)+t; 00343 result(2)=v(2)-t; 00344 result(3)=v(0)-v(1); 00345 return result; 00346 } 00347 // Apply conversion matrix from Stokes to linear. 00348 inline CStokesVector applySlin(const StokesVector& v) { 00349 CStokesVector result; 00350 return applySlin(result,v); 00351 } 00352 // Apply conversion matrix from Stokes to circular. 00353 inline CStokesVector& applyScirc(CStokesVector& result, 00354 const StokesVector& v) { 00355 Complex t=Complex(0.,1.0)*v(2); 00356 result(0)=v(0)+v(3); 00357 result(1)=v(1)+t; 00358 result(2)=v(1)-t; 00359 result(3)=v(0)-v(3); 00360 return result; 00361 } 00362 // Apply conversion matrix from Stokes to circular. 00363 inline CStokesVector applyScirc(const StokesVector& v) { 00364 CStokesVector result; 00365 return applyScirc(result,v); 00366 } 00367 00368 // Apply conversion matrix from linear to Stokes. 00369 inline StokesVector& applySlinInv(StokesVector& result, const CStokesVector& v) { 00370 result(0)=real(v(0)+v(3))/2; 00371 result(1)=real(v(0)-v(3))/2; 00372 result(2)=real(v(1)+v(2))/2; 00373 result(3)=real(Complex(0.,1.0)*(v(2)-v(1))/2); 00374 return result; 00375 } 00376 00377 // Apply conversion matrix from linear to Stokes. 00378 inline StokesVector applySlinInv(const CStokesVector& v) { 00379 StokesVector result; 00380 return applySlinInv(result,v); 00381 } 00382 00383 // Apply conversion matrix from circular to Stokes. 00384 inline StokesVector& applyScircInv(StokesVector& result, const CStokesVector& v) { 00385 result(0)=real(v(0)+v(3))/2; 00386 result(1)=real(v(1)+v(2))/2; 00387 result(2)=real(Complex(0.,1.0)*(v(2)-v(1))/2); 00388 result(3)=real(v(0)-v(3))/2; 00389 return result; 00390 } 00391 00392 // Apply conversion matrix from circular to Stokes. 00393 inline StokesVector applyScircInv(const CStokesVector& v) { 00394 StokesVector result; 00395 return applyScircInv(result,v); 00396 } 00397 00398 // The following are needed until Image no longer has 00399 // sigma images 00400 //StokesVector& sqrt(const StokesVector& v); 00401 00402 //CStokesVector& sqrt(const CStokesVector& v); 00403 00404 00405 } //# NAMESPACE CASA - END 00406 00407 #endif