source: branches/Release2.1/src/STPolLinear.cpp@ 2492

Last change on this file since 2492 was 1189, checked in by mar637, 18 years ago

Fix for ticket #61 - changing of default polarisation. Done ticket #62 as far as possible. Tid data gets automatically set as 'circular'\nAlso add circular pol class which was missing so far.

File size: 4.2 KB
Line 
1//
2// C++ Implementation: STPolLinear
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <Malte.Marquarding@csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <casa/Arrays/ArrayMath.h>
14#include <casa/BasicMath/Math.h>
15#include <casa/Exceptions/Error.h>
16#include <casa/BasicSL/Constants.h>
17#include "STPolLinear.h"
18
19using namespace casa;
20
21namespace asap {
22
23Factory<STPol,STPolLinear> STPolLinear::myFactory;
24
25STPolLinear::~STPolLinear()
26{
27}
28
29
30Vector<Float> asap::STPolLinear::getStokes( uint index )
31{
32 if ( index < 0 || index >4 ) throw(AipsError("Stokes index out of range"));
33 Vector<Float> out;
34 Float phase = getTotalPhase();
35 Vector<Float> q(getSpectrum(0) - getSpectrum(1));
36 if ( nspec() == 4 ) {
37 switch(index) {
38 case 0:
39 out = Vector<Float>(getSpectrum(0) + getSpectrum(1));
40 break;
41 case 1:
42 out = Vector<Float>(q * cos(phase) - Float(2.0)*getSpectrum(2) * sin(phase));
43 break;
44 case 2:
45 out = Vector<Float>(q * sin(phase) + Float(2.0)*getSpectrum(2) * cos(phase));
46 break;
47 case 3:
48 out = getFeedHand() * Float(2.0) * Vector<Float>(getSpectrum(3));
49 break;
50 }
51 } else if (nspec() == 2) {
52 if ( index == 0 )
53 out = Vector<Float>(getSpectrum(0) + getSpectrum(1));
54 }
55 return out;
56}
57
58Vector<Float> asap::STPolLinear::getLinPol( uInt index )
59{
60 if (nspec() != 4) {
61 throw(AipsError("You must have 4 linear polarizations to run this function"));
62 }
63 if ( index < 0 || index >4 ) throw(AipsError("LinPol index out of range"));
64 Vector<Float> out,q,u;
65 if ( nspec() == 4) {
66 switch(index) {
67 case 1:
68 q = getStokes(1);
69 u = getStokes(2);
70 out = Vector<Float>(sqrt(pow(q,Float(2.0))+pow(u, Float(2.0))));
71 break;
72 case 2:
73 q = getStokes(1);
74 u = getStokes(2);
75 out = Vector<Float>(Float(180.0/C::pi/2.0) * atan2(u,q));
76 break;
77 default:
78 out = getStokes(index);
79 }
80 }
81 return out;
82}
83
84Vector<Float> asap::STPolLinear::getLinear(uInt index )
85{
86 return getSpectrum(index);
87}
88
89Vector<Float> asap::STPolLinear::getCircular(uInt index )
90{
91 // convert to stokes I/ V first
92 //
93 // We use the convention
94 // I = (RR+LL) // definition changed
95 if (nspec() != 4) {
96 throw(AipsError("You must have 4 linear polarizations to run this function"));
97 }
98 if ( index == 2 || index ==3 ) throw(AipsError("Re/Imag RL not implemented"));
99 Vector<Float> I,V,out;
100 I = getStokes(0);
101 V = getStokes(3);
102 switch(index) {
103 case 0:
104 out = (I + V)/Float(2.0);
105 break;
106 case 1:
107 out = (I - V)/Float(2.0);
108 break;
109 default:
110 out = Vector<Float>();
111 }
112 return out;
113}
114
115void asap::STPolLinear::rotatePhase( Float phase )
116{
117 if (nspec() != 4) {
118 throw(AipsError("You must have 4 linear polarizations to run this function"));
119 }
120 Float cosVal = cos(C::pi/180.0*phase);
121 Float sinVal = sin(C::pi/180.0*phase);
122 Matrix<Float>& specs = getSpectra();
123 Vector<Float> R2 = specs.column(2)*cosVal - specs.column(3)*sinVal;
124 specs.column(3) = specs.column(2)*sinVal + specs.column(3)*cosVal;
125 specs.column(2) = R2;
126}
127
128void asap::STPolLinear::invertPhase( Float phase )
129{
130 // phase isnt used, just ro keep interface the same for all pol operations
131 if (nspec() != 4) {
132 throw(AipsError("You must have 4 linear polarizations to run this function"));
133 }
134 Matrix<Float>& specs = getSpectra();
135 Vector<Float> I = specs.column(3);
136 I *= Float(-1.0);
137}
138
139void asap::STPolLinear::rotateLinPolPhase( casa::Float phase )
140{
141//
142// Rotate P = Q + iU but do it directly on the linear
143// correlations.
144//
145// We are using I=(XX+YY)/2 convention
146// C1 = XX; C2 = YY, C3 = Real(XY)
147//
148 if (nspec() != 4) {
149 throw(AipsError("You must have 4 linear polarizations to run this function"));
150 }
151 Vector<Float> I,Q,U;
152 I = getStokes(0);
153 Q = getStokes(1);
154 U = getStokes(2);
155 // Rotate Q & U (factor of 2 for polarization)
156 Float cosVal = cos(C::pi/180.0*2.0*phase);
157 Float sinVal = sin(C::pi/180.0*2.0*phase);
158 Vector<Float> Q2 = Q*cosVal - U*sinVal;
159 U = Q*sinVal + U*cosVal;
160 Q = Q2;
161 Matrix<Float>& specs = getSpectra();
162 specs.column(0) = (I+Q)/Float(2.0);
163 specs.column(1) = (I-Q)/Float(2.0);
164 specs.column(2) = U/Float(2.0);
165
166}
167
168}
Note: See TracBrowser for help on using the repository browser.