source: trunk/src/SDPol.cc@ 480

Last change on this file since 480 was 469, checked in by kil064, 20 years ago

new functions for the SDWriter to handle Stokes output
for SDFITS / MS. This is bit of a mess...
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.1 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)
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
325Stokes::StokesTypes SDPolUtil::convertStokes(Int val, Bool toStokes, Bool linear)
326{
327 Stokes::StokesTypes stokes = Stokes::Undefined;
328 if (toStokes) {
329 if (val==0) {
330 stokes = Stokes::I;
331 } else if (val==1) {
332 stokes = Stokes::Q;
333 } else if (val==2) {
334 stokes = Stokes::U;
335 } else if (val==3) {
336 stokes = Stokes::V;
337 }
338 } else if (linear) {
339 if (val==0) {
340 stokes = Stokes::XX;
341 } else if (val==1) {
342 stokes = Stokes::YY;
343 } else if (val==2) {
344 stokes = Stokes::XY; // Real(XY)
345 } else if (val==3) {
346 stokes = Stokes::XY; // Imag(XY)
347 }
348 } else {
349 if (val==0) {
350 stokes = Stokes::RR;
351 } else if (val==1) {
352 stokes = Stokes::LL;
353 } else if (val==2) {
354 stokes = Stokes::RL;
355 } else if (val==3) {
356 stokes = Stokes::RL;
357 }
358 }
359//
360 return stokes;
361}
362
363Array<casa::uChar> SDPolUtil::andArrays (const Array<casa::uChar>& in1,
364 const Array<casa::uChar>& in2)
365{
366 Array<uChar> out(in1.shape());
367//
368 Array<uChar>::const_iterator in1Iter;
369 Array<uChar>::const_iterator in2Iter;
370 Array<uChar>::iterator outIter;
371//
372 for (in1Iter=in1.begin(),in2Iter=in2.begin(),outIter=out.begin();
373 in1Iter!=in1.end(); ++in1Iter,++in2Iter,++outIter) {
374 *outIter = *in1Iter & *in2Iter;
375 }
376 return out;
377}
378
379
380Array<Float> SDPolUtil::extractStokesForWriter (Array<Float>& in, const IPosition& start, const IPosition& end)
381//
382// start/end must already have applied the cursor selection of beam and IF
383// Extract specified Stokes for beam/IF and flip nChan and nPol for bloody SDwriter
384//
385{
386 IPosition shapeIn = in.shape();
387 uInt nChan = shapeIn(asap::ChanAxis);
388 uInt nPol = shapeIn(asap::PolAxis);
389//
390 IPosition shapeOut(2,nChan,nPol);
391 Array<Float> out(shapeOut);
392//
393 Array<Float> sliceRef = in(start,end); // Beam and IF now degenerate axes
394 ReadOnlyVectorIterator<Float> itIn(sliceRef, asap::ChanAxis);
395 VectorIterator<Float> itOut(out,0);
396 while (!itIn.pastEnd()) {
397 itOut.vector() = itIn.vector();
398//
399 itIn.next();
400 itOut.next();
401 }
402//
403 return out;
404}
405
Note: See TracBrowser for help on using the repository browser.