source: trunk/src/SDPol.cc@ 441

Last change on this file since 441 was 440, checked in by kil064, 20 years ago

add function SDPolUtil::stokesMask
fix error in nPol check in SDStokesENgine
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.3 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
61
62SDStokesEngine::SDStokesEngine (const Record& spec)
63: BaseMappedArrayEngine<Float,Float> ()
64{
65 if (spec.isDefined("OUTPUTNAME") && spec.isDefined("INPUTNAME")) {
66 setNames (spec.asString("OUTPUTNAME"), spec.asString("INPUTNAME"));
67 }
68}
69
70SDStokesEngine::SDStokesEngine (const SDStokesEngine& that)
71: BaseMappedArrayEngine<Float,Float> (that)
72{}
73
74SDStokesEngine::~SDStokesEngine()
75{}
76
77
78DataManager* SDStokesEngine::clone() const
79{
80 DataManager* dmPtr = new SDStokesEngine (*this);
81 return dmPtr;
82}
83
84
85String SDStokesEngine::dataManagerType() const
86{
87 return className();
88}
89
90String SDStokesEngine::className()
91{
92 return "SDStokesEngine";
93}
94
95String SDStokesEngine::dataManagerName() const
96{
97 return sourceName();
98}
99
100Record SDStokesEngine::dataManagerSpec() const
101{
102 Record spec;
103 spec.define ("OUTPUTNAME", sourceName()); // Ger uses opposite meaning for source/target
104 spec.define ("INPUTNAME", targetName());
105 return spec;
106}
107
108DataManager* SDStokesEngine::makeObject (const String&,
109 const Record& spec)
110{
111 DataManager* dmPtr = new SDStokesEngine(spec);
112 return dmPtr;
113}
114
115
116void SDStokesEngine::registerClass()
117{
118 DataManager::registerCtor (className(), makeObject);
119}
120
121
122void SDStokesEngine::create (uInt initialNrrow)
123{
124 BaseMappedArrayEngine<Float,Float>::create (initialNrrow);
125}
126
127void SDStokesEngine::prepare()
128{
129 BaseMappedArrayEngine<Float,Float>::prepare();
130}
131
132Bool SDStokesEngine::canAccessArrayColumnCells (Bool& reask) const
133{
134 reask = False;
135 return True;
136}
137
138
139void SDStokesEngine::getArray (uInt rownr, Array<Float>& output)
140{
141 Array<Float> input;
142 roColumn().get(rownr, input);
143//
144 computeOnGet (output, input);
145}
146
147void SDStokesEngine::putArray (uInt rownr, const Array<Float>& input)
148{
149 throw(AipsError("This Virtual Column is not writable"));
150}
151
152
153
154
155
156IPosition SDStokesEngine::shape (uInt rownr)
157{
158 IPosition inputShape = roColumn().shape (rownr);
159 return findOutputShape(inputShape);
160}
161
162
163
164void SDStokesEngine::computeOnGet(Array<Float>& output,
165 const Array<Float>& input)
166//
167// array of shape (nBeam,nIF,nPol,nChan)
168//
169// We use the scaling convention I=(XX+YY)
170//
171{
172
173// Checks
174
175 const uInt nDim = input.ndim();
176 AlwaysAssert(nDim==4,AipsError);
177 AlwaysAssert(output.ndim()==4,AipsError);
178//
179 const IPosition inputShape = input.shape();
180 const uInt polAxis = asap::PolAxis;
181 const uInt nPol = inputShape(polAxis);
182 AlwaysAssert(nPol==1 || nPol==2 || nPol==4, AipsError);
183
184// The silly Array slice operator does not give me back
185// a const reference so have to caste it away
186
187 Array<Float>& input2 = const_cast<Array<Float>&>(input);
188
189// Slice coordnates
190
191 IPosition start(nDim,0);
192 IPosition end(input.shape()-1);
193
194// Generate Slices
195
196 start(polAxis) = 0;
197 end(polAxis) = 0;
198 Array<Float> C1 = input2(start,end); // Input : C1
199//
200 start(polAxis) = 0;
201 end(polAxis) = 0;
202 Array<Float> I = output(start,end); // Output : I
203//
204 if (nPol==1) {
205 I = Float(2.0)*C1;
206 return;
207 }
208//
209 start(polAxis) = 1;
210 end(polAxis) = 1;
211 Array<Float> C2 = input2(start,end); // Input : C1
212//
213 I = C1 + C2;
214 if (nPol <= 2) return;
215//
216 start(polAxis) = 2;
217 end(polAxis) = 2;
218 Array<Float> C3 = input2(start,end); // Input : C3
219//
220 start(polAxis) = 3;
221 end(polAxis) = 3;
222 Array<Float> C4 = input2(start,end); // Input : C4
223//
224 start(polAxis) = 1;
225 end(polAxis) = 1;
226 Array<Float> Q = output(start,end); // Output : Q
227 Q = C1 - C2;
228//
229 start(polAxis) = 2;
230 end(polAxis) = 2;
231 Array<Float> U = output(start,end); // Output : U
232 U = Float(2.0)*C3;
233//
234 start(polAxis) = 3;
235 end(polAxis) = 3;
236 Array<Float> V = output(start,end); // Output : V
237 V = Float(2.0)*C4;
238}
239
240
241
242IPosition SDStokesEngine::findOutputShape (const IPosition& inputShape) const
243{
244 uInt axis = 2;
245 uInt nPol = inputShape(axis);
246 IPosition outputShape = inputShape;
247 if (nPol==1) {
248 outputShape(axis) = 1; // XX -> I
249 } else if (nPol==2) {
250 outputShape(axis) = 1; // XX YY -> I
251 } else if (nPol==4) {
252 outputShape(axis) = 4; // XX YY R(XY) I(XY) -> I Q U V
253 }
254 return outputShape;
255}
256
257
258
259// SDPolUtil
260
261Array<Float> SDPolUtil::polarizedIntensity (const Array<Float>& Q,
262 const Array<Float>& U)
263{
264 Array<Float> t1 = pow(Q,Double(2.0));
265 Array<Float> t2 = pow(U,Double(2.0));
266 return sqrt(t1+t2);
267}
268
269
270Array<Float> SDPolUtil::positionAngle (const Array<Float>& Q,
271 const Array<Float>& U)
272{
273 return Float(180.0/C::pi/2.0)*atan2(Q,U); // Degrees
274}
275
276
277void SDPolUtil::rotateXYPhase (Array<Float>& C3,
278 Array<Float>& C4,
279 Float phase)
280{
281 Float cosVal = cos(C::pi/180.0*phase);
282 Float sinVal = sin(C::pi/180.0*phase);
283//
284 C3 = C3*cosVal - C4*sinVal;
285 C4 = C3*sinVal + C4*cosVal;
286}
287
288
289
290Array<Float> SDPolUtil::getStokesSlice (Array<Float>& in, const IPosition& start,
291 const IPosition& end, const String& stokes)
292{
293 IPosition s(start);
294 IPosition e(end);
295//
296 if (stokes=="I") {
297 s(asap::PolAxis) = 0;
298 e(asap::PolAxis) = 0;
299 } else if (stokes=="Q") {
300 s(asap::PolAxis) = 1;
301 e(asap::PolAxis) = 1;
302 } else if (stokes=="U") {
303 s(asap::PolAxis) = 2;
304 e(asap::PolAxis) = 2;
305 } else if (stokes=="V") {
306 s(asap::PolAxis) = 3;
307 e(asap::PolAxis) = 3;
308 }
309//
310 return in(s,e);
311}
312
313
314Array<Float> SDPolUtil::circularPolarizationFromStokes (Array<Float>& I,
315 Array<Float>& V,
316 Bool doRR)
317{
318 if (doRR) {
319 return Float(0.5)*(I+V);
320 } else {
321 return Float(0.5)*(I-V);
322 }
323}
324
325Array<Bool> SDPolUtil::stokesMask (Array<Bool> rawFlags,
326 Bool doLinear)
327//
328// Generate mask for each Stokes parameter from the
329// raw flags. This is a lot of computational work and may
330// not be worth the effort.
331//
332{
333 IPosition shapeIn = rawFlags.shape();
334 uInt nPol = shapeIn(asap::PolAxis);
335 const uInt nDim = shapeIn.nelements();
336 Array<Bool> stokesFlags;
337//
338 IPosition start(nDim,0);
339 IPosition end(shapeIn-1);
340 IPosition shapeOut = shapeIn;
341//
342 if (doLinear) {
343 if (nPol==1) {
344 stokesFlags.resize(shapeOut);
345 stokesFlags = rawFlags;
346 } else if (nPol==2 || nPol==4) {
347
348// Set shape of output array
349
350 if (nPol==2) {
351 shapeOut(asap::PolAxis) = 1;
352 } else {
353 shapeOut(asap::PolAxis) = 4;
354 }
355 stokesFlags.resize(shapeOut);
356
357// Get reference slices and assign/compute
358
359 start(asap::PolAxis) = 0;
360 end(asap::PolAxis) = 0;
361 Array<Bool> M1In = rawFlags(start,end);
362//
363 start(asap::PolAxis) = 1;
364 end(asap::PolAxis) = 1;
365 Array<Bool> M2In = rawFlags(start,end);
366//
367 start(asap::PolAxis) = 0;
368 end(asap::PolAxis) = 0;
369 Array<Bool> M1Out = stokesFlags(start,end);
370 M1Out = M1In && M2In; // I
371//
372 if (nPol==4) {
373 start(asap::PolAxis) = 2;
374 end(asap::PolAxis) = 2;
375 Array<Bool> M3In = rawFlags(start,end);
376//
377 start(asap::PolAxis) = 3;
378 end(asap::PolAxis) = 3;
379 Array<Bool> M4In = rawFlags(start,end);
380//
381 start(asap::PolAxis) = 1;
382 end(asap::PolAxis) = 1;
383 Array<Bool> M2Out = stokesFlags(start,end);
384 M2Out = M1Out; // Q
385//
386 start(asap::PolAxis) = 2;
387 end(asap::PolAxis) = 2;
388 Array<Bool> M3Out = stokesFlags(start,end);
389 M3Out = M3In; // U
390//
391 start(asap::PolAxis) = 3;
392 end(asap::PolAxis) = 3;
393 Array<Bool> M4Out = stokesFlags(start,end);
394 M4Out = M4In; // V
395 }
396 } else {
397 throw(AipsError("Can only handle 1,2 or 4 polarizations"));
398 }
399 } else {
400 throw (AipsError("Only implemented for Linear polarizations"));
401 }
402//
403 return stokesFlags;
404}
Note: See TracBrowser for help on using the repository browser.