source: branches/Release-1-fixes/src/SDPol.cc@ 2468

Last change on this file since 2468 was 540, checked in by mar637, 20 years ago

Bug fix: missing break;
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.4 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/Arrays/VectorIter.h>
39#include <casa/Containers/Record.h>
40#include <casa/BasicSL/Constants.h>
41#include <casa/BasicSL/String.h>
42#include <casa/Utilities/Assert.h>
43
44#include <tables/Tables/Table.h>
45#include <tables/Tables/ScalarColumn.h>
46#include <tables/Tables/ArrayColumn.h>
47#include <tables/Tables/ColumnDesc.h>
48#include <tables/Tables/TableRecord.h>
49#include <tables/Tables/DataManError.h>
50
51
52using namespace casa;
53using namespace asap;
54
55
56
57SDStokesEngine::SDStokesEngine (const String& outputColumnName,
58 const String& inputColumnName)
59: BaseMappedArrayEngine<Float,Float> (outputColumnName, inputColumnName)
60{
61 setWritable(False);
62}
63
64
65SDStokesEngine::SDStokesEngine (const Record& spec)
66: BaseMappedArrayEngine<Float,Float> ()
67{
68 setWritable(False);
69 if (spec.isDefined("OUTPUTNAME") && spec.isDefined("INPUTNAME")) {
70 setNames (spec.asString("OUTPUTNAME"), spec.asString("INPUTNAME"));
71 }
72}
73
74SDStokesEngine::SDStokesEngine (const SDStokesEngine& that)
75: BaseMappedArrayEngine<Float,Float> (that)
76{}
77
78SDStokesEngine::~SDStokesEngine()
79{}
80
81
82DataManager* SDStokesEngine::clone() const
83{
84 DataManager* dmPtr = new SDStokesEngine (*this);
85 return dmPtr;
86}
87
88
89String SDStokesEngine::dataManagerType() const
90{
91 return className();
92}
93
94String SDStokesEngine::className()
95{
96 return "SDStokesEngine";
97}
98
99String SDStokesEngine::dataManagerName() const
100{
101 return virtualName();
102}
103
104Record SDStokesEngine::dataManagerSpec() const
105{
106 Record spec;
107 spec.define ("OUTPUTNAME", virtualName());
108 spec.define ("INPUTNAME", storedName());
109 return spec;
110}
111
112DataManager* SDStokesEngine::makeObject (const String&, const Record& spec)
113{
114 DataManager* dmPtr = new SDStokesEngine(spec);
115 return dmPtr;
116}
117
118
119void SDStokesEngine::registerClass()
120{
121 DataManager::registerCtor (className(), makeObject);
122}
123
124
125void SDStokesEngine::create (uInt initialNrrow)
126{
127 BaseMappedArrayEngine<Float,Float>::create (initialNrrow);
128}
129
130void SDStokesEngine::prepare()
131{
132 BaseMappedArrayEngine<Float,Float>::prepare();
133}
134
135Bool SDStokesEngine::canAccessArrayColumnCells (Bool& reask) const
136{
137 reask = False;
138 return True;
139}
140
141
142void SDStokesEngine::getArray (uInt rownr, Array<Float>& output)
143{
144 Array<Float> input;
145 roColumn().get(rownr, input);
146//
147 computeOnGet (output, input);
148}
149
150void SDStokesEngine::putArray (uInt rownr, const Array<Float>& input)
151{
152 throw(AipsError("This Virtual Column is not writable"));
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)/2
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 = 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 = Float(0.5) * (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 = Float(0.5) * (C1 - C2);
228//
229 start(polAxis) = 2;
230 end(polAxis) = 2;
231 Array<Float> U = output(start,end); // Output : U
232 U = C3;
233//
234 start(polAxis) = 3;
235 end(polAxis) = 3;
236 Array<Float> V = output(start,end); // Output : V
237 V = 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(U,Q); // Degrees
274}
275
276
277void SDPolUtil::rotatePhase (Array<Float>& R,
278 Array<Float>& I,
279 Float phase)
280//
281// Apply phase rotation to Z = (R + iI)
282//
283{
284 Float cosVal = cos(C::pi/180.0*phase);
285 Float sinVal = sin(C::pi/180.0*phase);
286//
287 Array<Float> R2 = R*cosVal - I*sinVal;
288 I = R*sinVal + I*cosVal;
289 R = R2;
290}
291
292
293void SDPolUtil::rotateLinPolPhase (Array<Float>& C1,
294 Array<Float>& C2,
295 Array<Float>& C3,
296 Array<Float>& I,
297 Array<Float>& Q,
298 Array<Float>& U,
299 Float phase)
300//
301// Rotate P = Q + iU but do it directly on the linear
302// correlations.
303//
304// We are using I=(XX+YY)/2 convention
305// C1 = XX; C2 = YY, C3 = Real(XY)
306//
307{
308// Rotate Q & U (factor of 2 for polarization)
309
310 rotatePhase(Q, U, 2.0*phase);
311
312// Now recompute C1,C2,C3
313// C4 unchanged
314
315 C1 = I + Q;
316 C2 = I - Q;
317 C3 = U;
318}
319
320
321
322
323
324Array<Float> SDPolUtil::getStokesSlice (Array<Float>& in, const IPosition& start,
325 const IPosition& end, const String& stokes)
326{
327 IPosition s(start);
328 IPosition e(end);
329//
330 if (stokes=="I") {
331 s(asap::PolAxis) = 0;
332 e(asap::PolAxis) = 0;
333 } else if (stokes=="Q") {
334 s(asap::PolAxis) = 1;
335 e(asap::PolAxis) = 1;
336 } else if (stokes=="U") {
337 s(asap::PolAxis) = 2;
338 e(asap::PolAxis) = 2;
339 } else if (stokes=="V") {
340 s(asap::PolAxis) = 3;
341 e(asap::PolAxis) = 3;
342 }
343//
344 return in(s,e);
345}
346
347
348Array<Float> SDPolUtil::circularPolarizationFromStokes (Array<Float>& I,
349 Array<Float>& V,
350 Bool doRR)
351//
352// We use the convention
353// I = (RR+LL)/2
354//
355{
356 if (doRR) {
357 return I + V;
358 } else {
359 return I - V;
360 }
361}
362
363Stokes::StokesTypes SDPolUtil::convertStokes(Int val, Bool toStokes, Bool linear)
364{
365 Stokes::StokesTypes stokes = Stokes::Undefined;
366 if (toStokes) {
367 if (val==0) {
368 stokes = Stokes::I;
369 } else if (val==1) {
370 stokes = Stokes::Q;
371 } else if (val==2) {
372 stokes = Stokes::U;
373 } else if (val==3) {
374 stokes = Stokes::V;
375 }
376 } else if (linear) {
377 if (val==0) {
378 stokes = Stokes::XX;
379 } else if (val==1) {
380 stokes = Stokes::YY;
381 } else if (val==2) {
382 stokes = Stokes::XY; // Real(XY)
383 } else if (val==3) {
384 stokes = Stokes::XY; // Imag(XY)
385 }
386 } else {
387 if (val==0) {
388 stokes = Stokes::RR;
389 } else if (val==1) {
390 stokes = Stokes::LL;
391 } else if (val==2) {
392 stokes = Stokes::RL;
393 } else if (val==3) {
394 stokes = Stokes::RL;
395 }
396 }
397//
398 return stokes;
399}
400
401
402
403String SDPolUtil::polarizationLabel (uInt polIdx, Bool linear, Bool stokes, Bool linPol)
404{
405 Stokes::StokesTypes type = Stokes::Undefined;
406 if (stokes) {
407 switch (polIdx) {
408 case 0:
409 {
410 type = Stokes::I;
411 }
412 break;
413 case 1:
414 {
415 if (linPol) {
416 type = Stokes::Plinear;
417 } else {
418 type = Stokes::Q;
419 }
420 }
421 break;
422 case 2:
423 {
424 if (linPol) {
425 type = Stokes::Pangle;
426 } else {
427 type = Stokes::U;
428 }
429 }
430 break;
431 case 3:
432 {
433 type = Stokes::V;
434 }
435 break;
436 default:
437 {
438 throw(AipsError("Unknown Stokes type"));
439 }
440 }
441 } else {
442 if (linear) {
443 switch (polIdx) {
444 case 0:
445 {
446 type = Stokes::XX;
447 }
448 break;
449 case 1:
450 {
451 type = Stokes::YY;
452 }
453 break;
454 case 2:
455 {
456 type = Stokes::XY; // Really Real(XY)
457 return String("Real(XY)");
458 }
459 break;
460 case 3:
461 {
462 type = Stokes::YX; // Really Imag(XY)
463 return String("Imag(XY)");
464 }
465 break;
466 default:
467 {
468 throw(AipsError("Unknown linear polarization type"));
469 }
470 }
471 } else {
472 switch (polIdx) {
473 case 0:
474 {
475 type = Stokes::RR;
476 }
477 break;
478 case 1:
479 {
480 type = Stokes::LL;
481 }
482 break;
483 case 2:
484 {
485 type = Stokes::RL; // Really Real(RL)
486 return String("Real(RL)");
487 }
488 break;
489 case 3:
490 {
491 type = Stokes::LR; // Really Imag(RL)
492 return String("Imag(RL)");
493 }
494 break;
495 default:
496 {
497 throw(AipsError("Unknown circular polarization type"));
498 }
499 }
500 }
501 }
502//
503 return SDPolUtil::stokesString(type);
504}
505
506
507
508
509// private
510
511String SDPolUtil::stokesString (Stokes::StokesTypes type)
512{
513 return Stokes::name (type);
514}
515
516
517Array<casa::uChar> SDPolUtil::andArrays (const Array<casa::uChar>& in1,
518 const Array<casa::uChar>& in2)
519{
520 Array<uChar> out(in1.shape());
521//
522 Array<uChar>::const_iterator in1Iter;
523 Array<uChar>::const_iterator in2Iter;
524 Array<uChar>::iterator outIter;
525//
526 for (in1Iter=in1.begin(),in2Iter=in2.begin(),outIter=out.begin();
527 in1Iter!=in1.end(); ++in1Iter,++in2Iter,++outIter) {
528 *outIter = *in1Iter & *in2Iter;
529 }
530 return out;
531}
532
533
534Array<Float> SDPolUtil::extractStokesForWriter (Array<Float>& in, const IPosition& start, const IPosition& end)
535//
536// start/end must already have applied the cursor selection of beam and IF
537// Extract specified Stokes for beam/IF and flip nChan and nPol for bloody SDwriter
538//
539{
540 IPosition shapeIn = in.shape();
541 uInt nChan = shapeIn(asap::ChanAxis);
542 uInt nPol = shapeIn(asap::PolAxis);
543//
544 IPosition shapeOut(2,nChan,nPol);
545 Array<Float> out(shapeOut);
546//
547 Array<Float> sliceRef = in(start,end); // Beam and IF now degenerate axes
548 ReadOnlyVectorIterator<Float> itIn(sliceRef, asap::ChanAxis);
549 VectorIterator<Float> itOut(out,0);
550 while (!itIn.pastEnd()) {
551 itOut.vector() = itIn.vector();
552//
553 itIn.next();
554 itOut.next();
555 }
556//
557 return out;
558}
559
560
561
562
563
Note: See TracBrowser for help on using the repository browser.