source: trunk/src/SDPol.cc@ 458

Last change on this file since 458 was 446, checked in by kil064, 20 years ago

add function convertStokes
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.2 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&, const Record& spec)
109{
110 DataManager* dmPtr = new SDStokesEngine(spec);
111 return dmPtr;
112}
113
114
115void SDStokesEngine::registerClass()
116{
117 DataManager::registerCtor (className(), makeObject);
118}
119
120
121void SDStokesEngine::create (uInt initialNrrow)
122{
123 BaseMappedArrayEngine<Float,Float>::create (initialNrrow);
124}
125
126void SDStokesEngine::prepare()
127{
128 BaseMappedArrayEngine<Float,Float>::prepare();
129}
130
131Bool SDStokesEngine::canAccessArrayColumnCells (Bool& reask) const
132{
133 reask = False;
134 return True;
135}
136
137
138void SDStokesEngine::getArray (uInt rownr, Array<Float>& output)
139{
140 Array<Float> input;
141 roColumn().get(rownr, input);
142//
143 computeOnGet (output, input);
144}
145
146void SDStokesEngine::putArray (uInt rownr, const Array<Float>& input)
147{
148 throw(AipsError("This Virtual Column is not writable"));
149}
150
151
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
324Array<Bool> SDPolUtil::stokesMask (Array<Bool> rawFlags,
325 Bool doLinear)
326//
327// Generate mask for each Stokes parameter from the
328// raw flags. This is a lot of computational work and may
329// not be worth the effort.
330//
331{
332 IPosition shapeIn = rawFlags.shape();
333 uInt nPol = shapeIn(asap::PolAxis);
334 const uInt nDim = shapeIn.nelements();
335 Array<Bool> stokesFlags;
336//
337 IPosition start(nDim,0);
338 IPosition end(shapeIn-1);
339 IPosition shapeOut = shapeIn;
340//
341 if (doLinear) {
342 if (nPol==1) {
343 stokesFlags.resize(shapeOut);
344 stokesFlags = rawFlags;
345 } else if (nPol==2 || nPol==4) {
346
347// Set shape of output array
348
349 if (nPol==2) {
350 shapeOut(asap::PolAxis) = 1;
351 } else {
352 shapeOut(asap::PolAxis) = 4;
353 }
354 stokesFlags.resize(shapeOut);
355
356// Get reference slices and assign/compute
357
358 start(asap::PolAxis) = 0;
359 end(asap::PolAxis) = 0;
360 Array<Bool> M1In = rawFlags(start,end);
361//
362 start(asap::PolAxis) = 1;
363 end(asap::PolAxis) = 1;
364 Array<Bool> M2In = rawFlags(start,end);
365//
366 start(asap::PolAxis) = 0;
367 end(asap::PolAxis) = 0;
368 Array<Bool> M1Out = stokesFlags(start,end);
369 M1Out = M1In && M2In; // I
370//
371 if (nPol==4) {
372 start(asap::PolAxis) = 2;
373 end(asap::PolAxis) = 2;
374 Array<Bool> M3In = rawFlags(start,end);
375//
376 start(asap::PolAxis) = 3;
377 end(asap::PolAxis) = 3;
378 Array<Bool> M4In = rawFlags(start,end);
379//
380 start(asap::PolAxis) = 1;
381 end(asap::PolAxis) = 1;
382 Array<Bool> M2Out = stokesFlags(start,end);
383 M2Out = M1Out; // Q
384//
385 start(asap::PolAxis) = 2;
386 end(asap::PolAxis) = 2;
387 Array<Bool> M3Out = stokesFlags(start,end);
388 M3Out = M3In; // U
389//
390 start(asap::PolAxis) = 3;
391 end(asap::PolAxis) = 3;
392 Array<Bool> M4Out = stokesFlags(start,end);
393 M4Out = M4In; // V
394 }
395 } else {
396 throw(AipsError("Can only handle 1,2 or 4 polarizations"));
397 }
398 } else {
399 throw (AipsError("Only implemented for Linear polarizations"));
400 }
401//
402 return stokesFlags;
403}
404
405Stokes::StokesTypes SDPolUtil::convertStokes(Int val, Bool toStokes, Bool linear)
406{
407 Stokes::StokesTypes stokes = Stokes::Undefined;
408 if (toStokes) {
409 if (val==0) {
410 stokes = Stokes::I;
411 } else if (val==1) {
412 stokes = Stokes::Q;
413 } else if (val==2) {
414 stokes = Stokes::U;
415 } else if (val==3) {
416 stokes = Stokes::V;
417 }
418 } else if (linear) {
419 if (val==0) {
420 stokes = Stokes::XX;
421 } else if (val==1) {
422 stokes = Stokes::YY;
423 } else if (val==2) {
424 stokes = Stokes::XY; // Real(XY)
425 } else if (val==3) {
426 stokes = Stokes::XY; // Imag(XY)
427 }
428 } else {
429 if (val==0) {
430 stokes = Stokes::RR;
431 } else if (val==1) {
432 stokes = Stokes::LL;
433 } else if (val==2) {
434 stokes = Stokes::RL;
435 } else if (val==3) {
436 stokes = Stokes::RL;
437 }
438 }
439//
440 return stokes;
441}
Note: See TracBrowser for help on using the repository browser.