00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #ifndef MSVIS_STOKESVECTOR_H
00030 #define MSVIS_STOKESVECTOR_H
00031
00032 #include <casa/aips.h>
00033 #include <casa/IO/AipsIO.h>
00034
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 {
00046
00047
00048 class StokesVector;
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 class CStokesVector:public RigidVector<Complex,4> {
00104 public:
00105 static String dataTypeId() {return "CStokesVector";};
00106
00107
00108 CStokesVector():RigidVector<Complex,4>() {}
00109
00110 CStokesVector(const Complex& c):RigidVector<Complex,4>(c) {}
00111
00112 CStokesVector(const Complex& v0, const Complex& v1,
00113 const Complex& v2, const Complex& v3):
00114 RigidVector<Complex,4>(v0,v1,v2,v3) {}
00115
00116 CStokesVector(const Complex v[4]):RigidVector<Complex,4>(v) {}
00117
00118
00119
00120 CStokesVector(const CStokesVector& v):RigidVector<Complex,4>(v){}
00121
00122
00123
00124 CStokesVector& operator=(const CStokesVector& v) {
00125 RigidVector<Complex,4>::operator=(v); return *this;
00126 }
00127
00128 CStokesVector& operator=(const Vector<Complex>& v) {
00129 RigidVector<Complex,4>::operator=(v); return *this;
00130 }
00131
00132 CStokesVector& operator=(const Complex& c) {
00133 RigidVector<Complex,4>::operator=(c); return *this;
00134 }
00135
00136 CStokesVector& operator-() {
00137 RigidVector<Complex,4>::operator-(); return *this;
00138 }
00139
00140 CStokesVector& operator+=(const CStokesVector& v) {
00141 RigidVector<Complex,4>::operator+=(v); return *this;
00142 }
00143
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
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
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
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
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
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
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
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
00197
00198
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
00211 friend ostream& operator<<(ostream& os, const CStokesVector& v) {
00212 os << v.vector();
00213 return os;
00214 }
00215 };
00216
00217
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
00233
00234 StokesVector():RigidVector<Float,4>() {}
00235
00236 StokesVector(Float f):RigidVector<Float,4>(f) {};
00237
00238 StokesVector(Float v0, Float v1, Float v2, Float v3): RigidVector<Float,4>(v0,v1,v2,v3){}
00239
00240
00241
00242 StokesVector(const StokesVector& v):RigidVector<Float,4>(v) {}
00243
00244
00245
00246 StokesVector& operator=(const StokesVector& v) {
00247 RigidVector<Float,4>::operator=(v); return *this;
00248 }
00249
00250 StokesVector& operator=(const Vector<Float>& v) {
00251 RigidVector<Float,4>::operator=(v); return *this;
00252 }
00253
00254 StokesVector& operator=(Float f) {
00255 RigidVector<Float,4>::operator=(f); return *this;
00256 }
00257
00258 StokesVector& operator-() {
00259 RigidVector<Float,4>::operator-(); return *this;
00260 }
00261
00262 StokesVector& operator+=(const StokesVector& v) {
00263 RigidVector<Float,4>::operator+=(v); return *this;
00264 }
00265
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
00276 StokesVector& operator*=(const SquareMatrix<Float,4>& m) {
00277 RigidVector<Float,4>::operator*=(m); return *this;
00278 }
00279
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
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
00290 Float maxEigenValue() const;
00291
00292 Float minEigenValue() const;
00293
00294 Float determinant() const;
00295
00296
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
00302 friend CStokesVector operator*(const SquareMatrix<Complex,4>& m,
00303 const StokesVector& v) {
00304
00305 return m * v;
00306 }
00307
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
00320 inline StokesVector operator*(Float f, const StokesVector& v) {
00321 StokesVector r(v);
00322 return r*=f;
00323 }
00324
00325 inline StokesVector operator*(const StokesVector& v, Float f) {
00326 StokesVector r(v);
00327 return r*=f;
00328 }
00329
00330
00331 inline StokesVector operator*(const SquareMatrix<Float,4>& m,
00332 const StokesVector& v) {
00333 StokesVector result(v);
00334 return result*=m;
00335 }
00336
00337
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
00348 inline CStokesVector applySlin(const StokesVector& v) {
00349 CStokesVector result;
00350 return applySlin(result,v);
00351 }
00352
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
00363 inline CStokesVector applyScirc(const StokesVector& v) {
00364 CStokesVector result;
00365 return applyScirc(result,v);
00366 }
00367
00368
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
00378 inline StokesVector applySlinInv(const CStokesVector& v) {
00379 StokesVector result;
00380 return applySlinInv(result,v);
00381 }
00382
00383
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
00393 inline StokesVector applyScircInv(const CStokesVector& v) {
00394 StokesVector result;
00395 return applyScircInv(result,v);
00396 }
00397
00398
00399
00400
00401
00402
00403
00404
00405 }
00406
00407 #endif