source: trunk/src/SDPol.cc@ 435

Last change on this file since 435 was 429, checked in by kil064, 20 years ago

add functions

getStokesSLice
circularPolarization
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.7 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/Containers/Record.h>
38#include <casa/BasicSL/Constants.h>
39#include <casa/BasicSL/String.h>
40#include <casa/Utilities/Assert.h>
41
42#include <tables/Tables/Table.h>
43#include <tables/Tables/ScalarColumn.h>
44#include <tables/Tables/ArrayColumn.h>
45#include <tables/Tables/ColumnDesc.h>
46#include <tables/Tables/TableRecord.h>
47#include <tables/Tables/DataManError.h>
48
49
50using namespace casa;
51using namespace asap;
52
53
54
55SDStokesEngine::SDStokesEngine (const String& outputColumnName,
56 const String& inputColumnName)
57: BaseMappedArrayEngine<Float,Float> (outputColumnName, inputColumnName)
58{}
59
60
61SDStokesEngine::SDStokesEngine (const Record& spec)
62: BaseMappedArrayEngine<Float,Float> ()
63{
64 if (spec.isDefined("OUTPUTNAME") && spec.isDefined("INPUTNAME")) {
65 setNames (spec.asString("OUTPUTNAME"), spec.asString("INPUTNAME"));
66 }
67}
68
69SDStokesEngine::SDStokesEngine (const SDStokesEngine& that)
70: BaseMappedArrayEngine<Float,Float> (that)
71{}
72
73SDStokesEngine::~SDStokesEngine()
74{}
75
76
77DataManager* SDStokesEngine::clone() const
78{
79 DataManager* dmPtr = new SDStokesEngine (*this);
80 return dmPtr;
81}
82
83
84String SDStokesEngine::dataManagerType() const
85{
86 return className();
87}
88
89String SDStokesEngine::className()
90{
91 return "SDStokesEngine";
92}
93
94String SDStokesEngine::dataManagerName() const
95{
96 return sourceName();
97}
98
99Record SDStokesEngine::dataManagerSpec() const
100{
101 Record spec;
102 spec.define ("OUTPUTNAME", sourceName()); // Ger uses opposite meaning for source/target
103 spec.define ("INPUTNAME", targetName());
104 return spec;
105}
106
107DataManager* SDStokesEngine::makeObject (const String&,
108 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 IPosition inputShape = findInputShape (output.shape());
141 Array<Float> input(inputShape);
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
154void SDStokesEngine::setShape (uInt rownr, const IPosition& outputShape)
155{
156 BaseMappedArrayEngine<Float,Float>::setShape (rownr, findInputShape(outputShape));
157}
158
159void SDStokesEngine::setShapeColumn (const IPosition& outputShape)
160{
161 BaseMappedArrayEngine<Float,Float>::setShapeColumn (findInputShape(outputShape));
162}
163
164
165IPosition SDStokesEngine::shape (uInt rownr)
166{
167 IPosition inputShape = roColumn().shape (rownr);
168 return findOutputShape(inputShape);
169}
170
171
172
173void SDStokesEngine::computeOnGet(Array<Float>& output,
174 const Array<Float>& input)
175//
176// array of shape (nBeam,nIF,nPol,nChan)
177//
178// We use the scaling convention I=(XX+YY)
179//
180{
181
182// Checks
183
184 const uInt nDim = input.ndim();
185 DebugAssert(nDim==4,AipsError);
186 DebugAssert(array.ndim()==4,AipsError);
187 const IPosition inputShape = input.shape();
188 const uInt polAxis = asap::PolAxis;
189 const uInt nPol = inputShape(polAxis);
190 DebugAssert(nPol==1 || nPol==2 || nPol==3, AipsError);
191
192// The silly Array slice operator does not give me back
193// a const reference so have to caste it away
194
195 Array<Float>& input2 = const_cast<Array<Float>&>(input);
196
197// Slice coordnates
198
199 IPosition start(nDim,0);
200 IPosition end(input.shape()-1);
201
202// Generate Slices
203
204 start(polAxis) = 0;
205 end(polAxis) = 0;
206 Array<Float> C1 = input2(start,end); // Input : C1
207//
208 start(polAxis) = 0;
209 end(polAxis) = 0;
210 Array<Float> I = output(start,end); // Output : I
211//
212 if (nPol==1) {
213 I = Float(2.0)*C1;
214 return;
215 }
216//
217 start(polAxis) = 1;
218 end(polAxis) = 1;
219 Array<Float> C2 = input2(start,end); // Input : C1
220//
221 I = C1 + C2;
222 if (nPol <= 2) return;
223//
224 start(polAxis) = 2;
225 end(polAxis) = 2;
226 Array<Float> C3 = input2(start,end); // Input : C3
227//
228 start(polAxis) = 3;
229 end(polAxis) = 3;
230 Array<Float> C4 = input2(start,end); // Input : C4
231//
232 start(polAxis) = 1;
233 end(polAxis) = 1;
234 Array<Float> Q = output(start,end); // Output : Q
235 Q = C1 - C2;
236//
237 start(polAxis) = 2;
238 end(polAxis) = 2;
239 Array<Float> U = output(start,end); // Output : U
240 U = Float(2.0)*C3;
241//
242 start(polAxis) = 3;
243 end(polAxis) = 3;
244 Array<Float> V = output(start,end); // Output : V
245 V = Float(2.0)*C4;
246}
247
248
249
250IPosition SDStokesEngine::findInputShape (const IPosition& outputShape) const
251//
252// Don't know how to handle the degeneracy that both
253// XX -> I
254// XX,YY -> I
255//
256{
257 uInt axis = asap::PolAxis;
258 uInt nPol = outputShape(axis);
259 IPosition inputShape = outputShape;
260 if (nPol==1) {
261 inputShape(axis) = 2; // XX YY -> I
262 } else if (nPol==4) {
263 inputShape(axis) = 4; // XX YY R(XY) I(XY) -> I Q U V
264 }
265 return inputShape;
266}
267
268IPosition SDStokesEngine::findOutputShape (const IPosition& inputShape) const
269{
270 uInt axis = 2;
271 uInt nPol = inputShape(axis);
272 IPosition outputShape = inputShape;
273 if (nPol==1) {
274 outputShape(axis) = 1; // XX -> I
275 } else if (nPol==2) {
276 outputShape(axis) = 1; // XX YY -> I
277 } else if (nPol==4) {
278 outputShape(axis) = 4; // XX YY R(XY) I(XY) -> I Q U V
279 }
280 return outputShape;
281}
282
283
284
285// SDPolUtil
286
287Array<Float> SDPolUtil::polarizedIntensity (const Array<Float>& Q,
288 const Array<Float>& U)
289{
290 Array<Float> t1 = pow(Q,Double(2.0));
291 Array<Float> t2 = pow(U,Double(2.0));
292 return sqrt(t1+t2);
293}
294
295
296Array<Float> SDPolUtil::positionAngle (const Array<Float>& Q,
297 const Array<Float>& U)
298{
299 return Float(180.0/C::pi/2.0)*atan2(Q,U); // Degrees
300}
301
302
303void SDPolUtil::rotateXYPhase (Array<Float>& C3,
304 Array<Float>& C4,
305 Float phase)
306{
307 Float cosVal = cos(C::pi/180.0*phase);
308 Float sinVal = sin(C::pi/180.0*phase);
309//
310 C3 = C3*cosVal - C4*sinVal;
311 C4 = C3*sinVal + C4*cosVal;
312}
313
314
315
316Array<Float> SDPolUtil::getStokesSlice (Array<Float>& in, const IPosition& start,
317 const IPosition& end, const String& stokes)
318{
319 IPosition s(start);
320 IPosition e(end);
321//
322 if (stokes=="I") {
323 s(asap::PolAxis) = 0;
324 e(asap::PolAxis) = 0;
325 } else if (stokes=="Q") {
326 s(asap::PolAxis) = 1;
327 e(asap::PolAxis) = 1;
328 } else if (stokes=="U") {
329 s(asap::PolAxis) = 2;
330 e(asap::PolAxis) = 2;
331 } else if (stokes=="V") {
332 s(asap::PolAxis) = 3;
333 e(asap::PolAxis) = 3;
334 }
335//
336 return in(s,e);
337}
338
339
340Array<Float> SDPolUtil::circularPolarizationFromStokes (Array<Float>& I,
341 Array<Float>& V,
342 Bool doRR)
343{
344 if (doRR) {
345 return Float(0.5)*(I+V);
346 } else {
347 return Float(0.5)*(I-V);
348 }
349}
350
Note: See TracBrowser for help on using the repository browser.