source: trunk/src/SDPol.cc@ 463

Last change on this file since 463 was 459, checked in by kil064, 20 years ago

make templatyed function SDPolUtil::stokesData to handle stokesMask
and stokesTsys calculation. Delete the non-templated versions.

In SDStokesEngone track interface changes in Table system (some bug fixes for me)

Mo
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.8 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDPol.cc: Polarimetric functionality
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31
32#include "SDPol.h"
33#include "SDDefs.h"
34
35#include <casa/Arrays/Array.h>
36#include <casa/Arrays/ArrayMath.h>
37#include <casa/Arrays/ArrayLogical.h>
38#include <casa/Containers/Record.h>
39#include <casa/BasicSL/Constants.h>
40#include <casa/BasicSL/String.h>
41#include <casa/Utilities/Assert.h>
42
43#include <tables/Tables/Table.h>
44#include <tables/Tables/ScalarColumn.h>
45#include <tables/Tables/ArrayColumn.h>
46#include <tables/Tables/ColumnDesc.h>
47#include <tables/Tables/TableRecord.h>
48#include <tables/Tables/DataManError.h>
49
50
51using namespace casa;
52using namespace asap;
53
54
55
56SDStokesEngine::SDStokesEngine (const String& outputColumnName,
57 const String& inputColumnName)
58: BaseMappedArrayEngine<Float,Float> (outputColumnName, inputColumnName)
59{
60 setWritable(False);
61}
62
63
64SDStokesEngine::SDStokesEngine (const Record& spec)
65: BaseMappedArrayEngine<Float,Float> ()
66{
67 setWritable(False);
68 if (spec.isDefined("OUTPUTNAME") && spec.isDefined("INPUTNAME")) {
69 setNames (spec.asString("OUTPUTNAME"), spec.asString("INPUTNAME"));
70 }
71}
72
73SDStokesEngine::SDStokesEngine (const SDStokesEngine& that)
74: BaseMappedArrayEngine<Float,Float> (that)
75{}
76
77SDStokesEngine::~SDStokesEngine()
78{}
79
80
81DataManager* SDStokesEngine::clone() const
82{
83 DataManager* dmPtr = new SDStokesEngine (*this);
84 return dmPtr;
85}
86
87
88String SDStokesEngine::dataManagerType() const
89{
90 return className();
91}
92
93String SDStokesEngine::className()
94{
95 return "SDStokesEngine";
96}
97
98String SDStokesEngine::dataManagerName() const
99{
100 return virtualName();
101}
102
103Record SDStokesEngine::dataManagerSpec() const
104{
105 Record spec;
106 spec.define ("OUTPUTNAME", virtualName());
107 spec.define ("INPUTNAME", storedName());
108 return spec;
109}
110
111DataManager* SDStokesEngine::makeObject (const String&, const Record& spec)
112{
113 DataManager* dmPtr = new SDStokesEngine(spec);
114 return dmPtr;
115}
116
117
118void SDStokesEngine::registerClass()
119{
120 DataManager::registerCtor (className(), makeObject);
121}
122
123
124void SDStokesEngine::create (uInt initialNrrow)
125{
126 BaseMappedArrayEngine<Float,Float>::create (initialNrrow);
127}
128
129void SDStokesEngine::prepare()
130{
131 BaseMappedArrayEngine<Float,Float>::prepare();
132}
133
134Bool SDStokesEngine::canAccessArrayColumnCells (Bool& reask) const
135{
136 reask = False;
137 return True;
138}
139
140
141void SDStokesEngine::getArray (uInt rownr, Array<Float>& output)
142{
143 Array<Float> input;
144 roColumn().get(rownr, input);
145//
146 computeOnGet (output, input);
147}
148
149void SDStokesEngine::putArray (uInt rownr, const Array<Float>& input)
150{
151 throw(AipsError("This Virtual Column is not writable"));
152}
153
154
155IPosition SDStokesEngine::shape (uInt rownr)
156{
157 IPosition inputShape = roColumn().shape (rownr);
158 return findOutputShape(inputShape);
159}
160
161
162
163void SDStokesEngine::computeOnGet(Array<Float>& output,
164 const Array<Float>& input)
165//
166// array of shape (nBeam,nIF,nPol,nChan)
167//
168// We use the scaling convention I=(XX+YY)
169//
170{
171
172// Checks
173
174 const uInt nDim = input.ndim();
175 AlwaysAssert(nDim==4,AipsError);
176 AlwaysAssert(output.ndim()==4,AipsError);
177//
178 const IPosition inputShape = input.shape();
179 const uInt polAxis = asap::PolAxis;
180 const uInt nPol = inputShape(polAxis);
181 AlwaysAssert(nPol==1 || nPol==2 || nPol==4, AipsError);
182
183// The silly Array slice operator does not give me back
184// a const reference so have to caste it away
185
186 Array<Float>& input2 = const_cast<Array<Float>&>(input);
187
188// Slice coordnates
189
190 IPosition start(nDim,0);
191 IPosition end(input.shape()-1);
192
193// Generate Slices
194
195 start(polAxis) = 0;
196 end(polAxis) = 0;
197 Array<Float> C1 = input2(start,end); // Input : C1
198//
199 start(polAxis) = 0;
200 end(polAxis) = 0;
201 Array<Float> I = output(start,end); // Output : I
202//
203 if (nPol==1) {
204 I = Float(2.0)*C1;
205 return;
206 }
207//
208 start(polAxis) = 1;
209 end(polAxis) = 1;
210 Array<Float> C2 = input2(start,end); // Input : C1
211//
212 I = C1 + C2;
213 if (nPol <= 2) return;
214//
215 start(polAxis) = 2;
216 end(polAxis) = 2;
217 Array<Float> C3 = input2(start,end); // Input : C3
218//
219 start(polAxis) = 3;
220 end(polAxis) = 3;
221 Array<Float> C4 = input2(start,end); // Input : C4
222//
223 start(polAxis) = 1;
224 end(polAxis) = 1;
225 Array<Float> Q = output(start,end); // Output : Q
226 Q = C1 - C2;
227//
228 start(polAxis) = 2;
229 end(polAxis) = 2;
230 Array<Float> U = output(start,end); // Output : U
231 U = Float(2.0)*C3;
232//
233 start(polAxis) = 3;
234 end(polAxis) = 3;
235 Array<Float> V = output(start,end); // Output : V
236 V = Float(2.0)*C4;
237}
238
239
240
241IPosition SDStokesEngine::findOutputShape (const IPosition& inputShape) const
242{
243 uInt axis = 2;
244 uInt nPol = inputShape(axis);
245 IPosition outputShape = inputShape;
246 if (nPol==1) {
247 outputShape(axis) = 1; // XX -> I
248 } else if (nPol==2) {
249 outputShape(axis) = 1; // XX YY -> I
250 } else if (nPol==4) {
251 outputShape(axis) = 4; // XX YY R(XY) I(XY) -> I Q U V
252 }
253 return outputShape;
254}
255
256
257
258// SDPolUtil
259
260Array<Float> SDPolUtil::polarizedIntensity (const Array<Float>& Q,
261 const Array<Float>& U)
262{
263 Array<Float> t1 = pow(Q,Double(2.0));
264 Array<Float> t2 = pow(U,Double(2.0));
265 return sqrt(t1+t2);
266}
267
268
269Array<Float> SDPolUtil::positionAngle (const Array<Float>& Q,
270 const Array<Float>& U)
271{
272 return Float(180.0/C::pi/2.0)*atan2(Q,U); // Degrees
273}
274
275
276void SDPolUtil::rotateXYPhase (Array<Float>& C3,
277 Array<Float>& C4,
278 Float phase)
279{
280 Float cosVal = cos(C::pi/180.0*phase);
281 Float sinVal = sin(C::pi/180.0*phase);
282//
283 C3 = C3*cosVal - C4*sinVal;
284 C4 = C3*sinVal + C4*cosVal;
285}
286
287
288
289Array<Float> SDPolUtil::getStokesSlice (Array<Float>& in, const IPosition& start,
290 const IPosition& end, const String& stokes)
291{
292 IPosition s(start);
293 IPosition e(end);
294//
295 if (stokes=="I") {
296 s(asap::PolAxis) = 0;
297 e(asap::PolAxis) = 0;
298 } else if (stokes=="Q") {
299 s(asap::PolAxis) = 1;
300 e(asap::PolAxis) = 1;
301 } else if (stokes=="U") {
302 s(asap::PolAxis) = 2;
303 e(asap::PolAxis) = 2;
304 } else if (stokes=="V") {
305 s(asap::PolAxis) = 3;
306 e(asap::PolAxis) = 3;
307 }
308//
309 return in(s,e);
310}
311
312
313Array<Float> SDPolUtil::circularPolarizationFromStokes (Array<Float>& I,
314 Array<Float>& V,
315 Bool doRR)
316{
317 if (doRR) {
318 return Float(0.5)*(I+V);
319 } else {
320 return Float(0.5)*(I-V);
321 }
322}
323
324Stokes::StokesTypes SDPolUtil::convertStokes(Int val, Bool toStokes, Bool linear)
325{
326 Stokes::StokesTypes stokes = Stokes::Undefined;
327 if (toStokes) {
328 if (val==0) {
329 stokes = Stokes::I;
330 } else if (val==1) {
331 stokes = Stokes::Q;
332 } else if (val==2) {
333 stokes = Stokes::U;
334 } else if (val==3) {
335 stokes = Stokes::V;
336 }
337 } else if (linear) {
338 if (val==0) {
339 stokes = Stokes::XX;
340 } else if (val==1) {
341 stokes = Stokes::YY;
342 } else if (val==2) {
343 stokes = Stokes::XY; // Real(XY)
344 } else if (val==3) {
345 stokes = Stokes::XY; // Imag(XY)
346 }
347 } else {
348 if (val==0) {
349 stokes = Stokes::RR;
350 } else if (val==1) {
351 stokes = Stokes::LL;
352 } else if (val==2) {
353 stokes = Stokes::RL;
354 } else if (val==3) {
355 stokes = Stokes::RL;
356 }
357 }
358//
359 return stokes;
360}
Note: See TracBrowser for help on using the repository browser.