Line data Source code
1 :
2 : //# RFDataMapper.cc: this defines RFDataMapper
3 : //# Copyright (C) 2000,2001
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: casa-feedback@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : #include <flagging/Flagging/RFDataMapper.h>
29 : #include <casacore/casa/Arrays/Slice.h>
30 : #include <casacore/casa/Arrays/Cube.h>
31 : #include <casacore/casa/Arrays/ArrayMath.h>
32 : #include <casacore/casa/Exceptions/Error.h>
33 : #include <casacore/casa/BasicMath/Math.h>
34 : #include <casacore/casa/BasicSL/Constants.h>
35 : #include <msvis/MSVis/VisBuffer.h>
36 : #include <msvis/MSVis/VisibilityIterator.h>
37 : #include <casacore/measures/Measures/MDirection.h>
38 :
39 : using namespace casacore;
40 : namespace casa { //# NAMESPACE CASA - BEGIN
41 :
42 : /*
43 : // Define functions for mapping a VisBuffer to obsevred, model and corrected
44 : // visibility cubes
45 : #define CubeMapper(name,method) \
46 : static Cube<Complex> * CubeMap##name ( VisBuffer &vb ) \
47 : { return &vb.method(); }
48 : CubeMapper(Obs,visCube);
49 : CubeMapper(Model,modelVisCube);
50 : CubeMapper(Corrected,correctedVisCube);
51 : #undef CubeMapper
52 :
53 : // Define functions for mapping a VisBuffer to residual
54 : // visibility cube
55 :
56 : #define CubeMapper1(name,method1,method2) \
57 : static Cube<Complex> * CubeMap##name ( VisBuffer &vb ) \
58 : { return &(vb.method1() - vb.method2()); }
59 : CubeMapper(ResCorrected,correctedVisCube,modelVisCube);
60 : CubeMapper(ResObs,visCube,modelVisCube);
61 : #undef CubeMapper1
62 : */
63 :
64 : /// before assigning pviscube..
65 : /* vb.visCube() = vb.correctedVisCube() - vb.modelVisCube() */
66 0 : static Cube<Complex> * CubeMapObs ( VisBuffer &vb ){ return &vb.visCube(); }
67 0 : static Cube<Complex> * CubeMapModel ( VisBuffer &vb ){ return &vb.modelVisCube(); }
68 0 : static Cube<Complex> * CubeMapCorrected ( VisBuffer &vb ){ return &vb.correctedVisCube(); }
69 0 : static Cube<Complex> * CubeMapResCorrected ( VisBuffer &vb )
70 0 : { vb.visCube() = vb.correctedVisCube() - vb.modelVisCube();
71 0 : return &vb.visCube(); }
72 : //{ return &(vb.correctedVisCube() - vb.modelVisCube()); }
73 0 : static Cube<Complex> * CubeMapResObs ( VisBuffer &vb )
74 0 : { vb.visCube() = vb.visCube() - vb.modelVisCube();
75 0 : return &vb.visCube(); }
76 : //{ return &(vb.visCube() - vb.modelVisCube()); }
77 :
78 :
79 : // a dummyRowMapper for uninitialized row mappings
80 0 : Float RFDataMapper::dummyRowMapper (uInt)
81 : {
82 0 : throw( AipsError("DataMapper: call to unititialized RowMapper") );
83 : }
84 :
85 : // define a bunch of row mapper for various UVW stuff
86 : #define UVW ((*puvw)(ir))
87 : #define UVRowMapper(name,expr) \
88 : Float RFDataMapper::name##_RowMapper (uInt ir) \
89 : { return expr; }
90 0 : UVRowMapper(U,UVW(0));
91 0 : UVRowMapper(V,UVW(1));
92 0 : UVRowMapper(W,UVW(2));
93 0 : UVRowMapper(AbsU,abs(UVW(0)));
94 0 : UVRowMapper(AbsV,abs(UVW(1)));
95 0 : UVRowMapper(AbsW,abs(UVW(2)));
96 0 : UVRowMapper(UVD,sqrt(UVW(0)*UVW(0)+UVW(1)*UVW(1)));
97 0 : UVRowMapper(UVA,atan2(UVW(0),UVW(1))/C::pi*180);
98 0 : UVRowMapper(HA,sin_dec!=0 ? atan2(UVW(1)/sin_dec,UVW(0))/C::pi*180 : 0 );
99 :
100 : // these arrays define a mapping between column names and cube mappers
101 : const String
102 : COL_ID[] = { "OBS", "DATA", "MODEL",
103 : "CORR", "CORRECTED",
104 : "RES",
105 : "RES_CORR", "RES_CORRECTED",
106 : "RES_OBS", "RES_DATA"};
107 : const RFDataMapper::CubeMapperFunc
108 : COL_MAP[] = { &CubeMapObs, &CubeMapObs, &CubeMapModel,
109 : &CubeMapCorrected, &CubeMapCorrected,
110 : &CubeMapResCorrected,
111 : &CubeMapResCorrected, &CubeMapResCorrected,
112 : &CubeMapResObs, &CubeMapResObs};
113 :
114 : // -----------------------------------------------------------------------
115 : // RFDataMapper::getCubeMapper
116 : // Maps column name to a cube mapper function
117 : // -----------------------------------------------------------------------
118 0 : RFDataMapper::CubeMapperFunc RFDataMapper::getCubeMapper( const String &column,Bool throw_excp )
119 : {
120 : // setup cube mapper
121 0 : if( !column.length() )
122 0 : return COL_MAP[0];
123 0 : String col( upcase(column) );
124 0 : for( uInt i=0; i<sizeof(COL_ID)/sizeof(COL_ID[0]); i++ ) {
125 0 : if( col.matches(COL_ID[i]) )
126 0 : return COL_MAP[i];
127 : }
128 0 : if( throw_excp )
129 0 : throw( AipsError("DataMapper: unknown column "+column) );
130 0 : return NULL;
131 0 : }
132 :
133 :
134 : // -----------------------------------------------------------------------
135 : // Constructor 1
136 : // For visibility expressions
137 : // -----------------------------------------------------------------------
138 0 : RFDataMapper::RFDataMapper ( const String &colmn,DDMapper *map )
139 0 : : desc(""),
140 0 : ddm(map),
141 : // rowmapper(NULL),
142 0 : cubemap(getCubeMapper(colmn,true)),
143 0 : mytype(MAPCORR)
144 : {
145 0 : sin_dec = -10; // need not be computed by default, so set to <-1
146 0 : full_cycle = cycle_base = 0; // by default, mapped values are not cyclic
147 0 : rowmapper = &RFDataMapper::dummyRowMapper;
148 0 : }
149 :
150 : // -----------------------------------------------------------------------
151 : // Constructor 2
152 : // For visibility expressions, with explicit column specification
153 : // -----------------------------------------------------------------------
154 0 : RFDataMapper::RFDataMapper ( const Vector<String> &expr0,const String &defcol )
155 0 : : expr_desc(""),
156 0 : ddm(NULL),
157 0 : cubemap(NULL),
158 : // rowmapper(dummyRowMapper),
159 0 : mytype(MAPCORR)
160 : {
161 0 : sin_dec = -10; // need not be computed by default, so set to <-1
162 0 : full_cycle = cycle_base = 0; // by default, mapped values are not cyclic
163 0 : rowmapper = &RFDataMapper::dummyRowMapper;
164 0 : Vector<String> expr( splitExpression(expr0) );
165 : // first, check for per-row expressions (i.e., UVW, etc.)
166 : // at this point, assume we have a per-row expression (i.e. UVW, etc.)
167 0 : Bool absof=false;
168 0 : String el = upcase(expr(0));
169 0 : if( el == "ABS" )
170 : {
171 0 : absof=true;
172 0 : if( expr.nelements()<2 )
173 0 : throw(AipsError("DataMapper: illegal expression "+expr(0)));
174 0 : el = upcase( expr(1) );
175 : }
176 0 : else if( el.matches("ABS",0) )
177 : {
178 0 : absof = true;
179 0 : el = el.after(3);
180 : }
181 :
182 0 : if( el == "U" )
183 0 : rowmapper = absof ? &RFDataMapper::AbsU_RowMapper : &RFDataMapper::U_RowMapper;
184 0 : else if( el == "V" )
185 0 : rowmapper = absof ? &RFDataMapper::AbsV_RowMapper : &RFDataMapper::V_RowMapper;
186 0 : else if( el == "W" )
187 0 : rowmapper = absof ? &RFDataMapper::AbsW_RowMapper : &RFDataMapper::W_RowMapper;
188 0 : else if( el == "UVD" || el == "UVDIST" )
189 0 : rowmapper = &RFDataMapper::UVD_RowMapper;
190 0 : else if( el == "UVA" || el == "UVANGLE" )
191 : {
192 0 : rowmapper = &RFDataMapper::UVA_RowMapper;
193 0 : full_cycle = 360; // mapping into angles
194 0 : cycle_base = -180;
195 : }
196 0 : else if( el == "HA" || el == "HANGLE" )
197 : {
198 0 : sin_dec = 0; // will need to be computed
199 0 : rowmapper = &RFDataMapper::HA_RowMapper;
200 0 : full_cycle = 360; // mapping into angles
201 0 : cycle_base = -180;
202 : }
203 : else
204 0 : rowmapper = NULL;
205 : // have we set up a valid row mapper? Return then
206 0 : if( rowmapper )
207 : {
208 0 : desc = absof ? "|"+el+"|" : el;
209 0 : expr_desc = desc;
210 0 : mytype = MAPROW;
211 0 : return;
212 : }
213 : // at this point, it must be a valid correlation expression
214 0 : String column(defcol);
215 :
216 : // see if expression starts with a non-empty column specification, if so,
217 : // remember the column, and shift it out of the expression vector
218 0 : CubeMapperFunc cm = getCubeMapper(expr(0));
219 :
220 0 : if( cm && expr(0).length() )
221 : {
222 0 : column = expr(0);
223 0 : expr = expr(Slice(1,expr.nelements()-1));
224 : }
225 :
226 : // check if it parses to a valid DDMapper expression
227 0 : ddm = DDFunc::getMapper(expr_desc,expr);
228 : // valid expression? Set ourselves up as a correlation mapper then
229 :
230 0 : if( ddm )
231 : {
232 0 : if( !cm ) // set column from defcol if not set above
233 0 : cm = getCubeMapper(defcol,true);
234 0 : cubemap = cm;
235 0 : desc = (column.length() ? "("+upcase(column)+")" : String("") )+expr_desc;
236 0 : mytype = MAPCORR;
237 :
238 0 : return;
239 : }
240 :
241 : // invalid expression, so throw an exception
242 0 : String s;
243 0 : for( uInt i=0; i<expr.nelements(); i++ )
244 : {
245 0 : if( i )
246 0 : s+=" ";
247 0 : s+=expr(i);
248 : }
249 0 : throw(AipsError("DataMapper: illegal expression "+s));
250 0 : }
251 :
252 0 : RFDataMapper::~RFDataMapper ()
253 : {
254 0 : if( ddm )
255 0 : delete ddm;
256 0 : };
257 :
258 : // computes a correlations mask
259 0 : RFlagWord RFDataMapper::corrMask ( const VisibilityIterator &vi )
260 : {
261 0 : Vector<Int> corr;
262 0 : vi.corrType(corr);
263 : // for a visibilities mapper, ask the DDMapper
264 0 : if( mytype == MAPCORR )
265 : {
266 0 : if( !ddm->reset( corr ) )
267 0 : return 0;
268 0 : return (RFlagWord) ddm->corrMask();
269 : }
270 : // for a row mapper, flag all correlations
271 0 : else if( mytype == MAPROW )
272 : {
273 0 : return (1<<corr.nelements())-1;
274 : }
275 : else // paranoid case
276 0 : throw( AipsError( "DataMapper: unknown mytype") );
277 0 : }
278 :
279 : // sets up for a visbuffer
280 0 : void RFDataMapper::setVisBuffer ( VisBuffer &vb )
281 : {
282 0 : if( mytype == MAPCORR )
283 0 : pviscube = (*cubemap)(vb);
284 0 : else if( mytype == MAPROW )
285 0 : puvw = &vb.uvw();
286 : else
287 0 : throw( AipsError( "DataMapper: unknown mytype") );
288 : // extract sine of declination, if needed
289 0 : if( sin_dec>=-1 )
290 : {
291 0 : sin_dec = sin( MDirection::Convert( vb.phaseCenter(),
292 0 : MDirection::Ref(MDirection::J2000)
293 0 : )().getAngle().getBaseValue()(1) );
294 : }
295 0 : }
296 :
297 :
298 :
299 : } //# NAMESPACE CASA - END
300 :
|
