source: branches/Release12/src/SDPol2.cc@ 2660

Last change on this file since 2660 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: 7.1 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDPol2.cc: Templated 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/BasicSL/String.h>
39#include <casa/Utilities/Assert.h>
40#include <casa/Utilities/DataType.h>
41#include <casa/Utilities/Assert.h>
42
43
44using namespace casa;
45using namespace asap;
46
47
48template <class T>
49Array<T> SDPolUtil::stokesData (Array<T>& rawData, Bool doLinear)
50//
51// Generate data for each Stokes parameter from the
52// raw flags. This is a lot of computational work and may
53// not be worth the effort.
54//
55// Designed for use
56// Bool or uChar for mask
57// Float for TSys
58//
59// Input arrays should be of shape
60// [nBeam,nIF,nPol,nChan]
61//
62{
63 T* t;
64 DataType type = whatType(t);
65 AlwaysAssert(type==TpFloat || type==TpBool || type==TpUChar, AipsError);
66//
67 IPosition shapeIn = rawData.shape();
68 const uInt nDim = shapeIn.nelements();
69 uInt polAxis = 0;
70 AlwaysAssert(nDim==asap::nAxes, AipsError);
71//
72 uInt nPol = shapeIn(asap::PolAxis);
73 Array<T> stokesData;
74//
75 IPosition start(nDim,0);
76 IPosition end(shapeIn-1);
77 IPosition shapeOut = shapeIn;
78//
79 if (doLinear) {
80 if (nPol==1) {
81 stokesData.resize(shapeOut);
82 stokesData = rawData;
83 } else if (nPol==2 || nPol==4) {
84
85// Set shape of output array
86
87 if (nPol==2) {
88 shapeOut(asap::PolAxis) = 1;
89 } else {
90 shapeOut(asap::PolAxis) = 4;
91 }
92 stokesData.resize(shapeOut);
93
94// Get reference slices and assign/compute
95
96 start(asap::PolAxis) = 0;
97 end(asap::PolAxis) = 0;
98 Array<T> M1In = rawData(start,end);
99//
100 start(asap::PolAxis) = 1;
101 end(asap::PolAxis) = 1;
102 Array<T> M2In = rawData(start,end);
103//
104 start(asap::PolAxis) = 0;
105 end(asap::PolAxis) = 0;
106 Array<T> M1Out = stokesData(start,end);
107 M1Out = SDPolUtil::andArrays (M1In, M2In); // I
108//
109 if (nPol==4) {
110 start(asap::PolAxis) = 2;
111 end(asap::PolAxis) = 2;
112 Array<T> M3In = rawData(start,end);
113//
114 start(asap::PolAxis) = 3;
115 end(asap::PolAxis) = 3;
116 Array<T> M4In = rawData(start,end);
117//
118 start(asap::PolAxis) = 1;
119 end(asap::PolAxis) = 1;
120 Array<T> M2Out = stokesData(start,end);
121 M2Out = M1Out; // Q
122//
123 start(asap::PolAxis) = 2;
124 end(asap::PolAxis) = 2;
125 Array<T> M3Out = stokesData(start,end);
126 M3Out = M3In; // U
127//
128 start(asap::PolAxis) = 3;
129 end(asap::PolAxis) = 3;
130 Array<T> M4Out = stokesData(start,end);
131 M4Out = M4In; // V
132 }
133 } else {
134 throw(AipsError("Can only handle 1,2 or 4 polarizations"));
135 }
136 } else {
137 throw (AipsError("Only implemented for Linear polarizations"));
138 }
139//
140 return stokesData;
141}
142
143
144
145template <class T>
146Array<T> SDPolUtil::computeStokesDataForWriter (Array<T>& rawData, Bool doLinear)
147//
148// Generate data for each Stokes parameter from the
149// raw flags. This is a lot of computational work and may
150// not be worth the effort. This function is specifically
151// for the SDWriter
152//
153// Designed for use
154// Bool or uChar for mask
155// Float for TSys
156//
157// Input arrays should be of shape
158// [nChan,nPol] Bool/uChar
159// [nPol] Float
160//
161{
162 T* t;
163 DataType type = whatType(t);
164 AlwaysAssert(type==TpFloat || type==TpBool || type==TpUChar, AipsError);
165//
166 IPosition shapeIn = rawData.shape();
167 const uInt nDim = shapeIn.nelements();
168 uInt polAxis = 0;
169 if (type==TpFloat) {
170 AlwaysAssert(nDim==1,AipsError);
171 } else {
172 AlwaysAssert(nDim==2,AipsError);
173 polAxis = 1;
174 }
175//
176 uInt nPol = shapeIn(polAxis);
177 Array<T> stokesData;
178//
179 IPosition start(nDim,0);
180 IPosition end(shapeIn-1);
181 IPosition shapeOut = shapeIn;
182//
183 if (doLinear) {
184 if (nPol==1) {
185 stokesData.resize(shapeOut);
186 stokesData = rawData;
187 } else if (nPol==2 || nPol==4) {
188
189// Set shape of output array
190
191 if (nPol==2) {
192 shapeOut(polAxis) = 1;
193 } else {
194 shapeOut(polAxis) = 4;
195 }
196 stokesData.resize(shapeOut);
197
198// Get reference slices and assign/compute
199
200 start(polAxis) = 0;
201 end(polAxis) = 0;
202 Array<T> M1In = rawData(start,end);
203//
204 start(polAxis) = 1;
205 end(polAxis) = 1;
206 Array<T> M2In = rawData(start,end);
207//
208 start(polAxis) = 0;
209 end(polAxis) = 0;
210 Array<T> M1Out = stokesData(start,end);
211 M1Out = SDPolUtil::andArrays (M1In, M2In); // I
212//
213 if (nPol==4) {
214 start(polAxis) = 2;
215 end(polAxis) = 2;
216 Array<T> M3In = rawData(start,end);
217//
218 start(polAxis) = 3;
219 end(polAxis) = 3;
220 Array<T> M4In = rawData(start,end);
221//
222 start(polAxis) = 1;
223 end(polAxis) = 1;
224 Array<T> M2Out = stokesData(start,end);
225 M2Out = M1Out; // Q
226//
227 start(polAxis) = 2;
228 end(polAxis) = 2;
229 Array<T> M3Out = stokesData(start,end);
230 M3Out = M3In; // U
231//
232 start(polAxis) = 3;
233 end(polAxis) = 3;
234 Array<T> M4Out = stokesData(start,end);
235 M4Out = M4In; // V
236 }
237 } else {
238 throw(AipsError("Can only handle 1,2 or 4 polarizations"));
239 }
240 } else {
241 throw (AipsError("Only implemented for Linear polarizations"));
242 }
243//
244 return stokesData;
245}
246
247
Note: See TracBrowser for help on using the repository browser.