casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
StokesVector.h
Go to the documentation of this file.
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