Line data Source code
1 : //# VisBuffer.cc: buffer for iterating through MS in large blocks
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <msvis/MSVis/VisibilityIterator.h>
29 : #include <msvis/MSVis/VisBuffer.h>
30 : #include <msvis/MSVis/VisBufferAsyncWrapper.h>
31 : #include <stdcasa/UtilJ.h>
32 : #include <casacore/casa/Arrays/ArrayMath.h>
33 : #include <casacore/casa/Arrays/ArrayLogical.h>
34 : #include <casacore/casa/Arrays/MaskedArray.h>
35 : #include <casacore/casa/Arrays/MaskArrMath.h>
36 : #include <casacore/casa/Arrays/ArrayUtil.h>
37 : #include <casacore/casa/OS/Path.h>
38 : #include <casacore/casa/Utilities/Assert.h>
39 : #include <casacore/casa/Utilities/GenSort.h>
40 : #include <casacore/casa/OS/Timer.h>
41 : #include <casacore/ms/MeasurementSets/MSColumns.h>
42 : #include <casacore/ms/MeasurementSets/MSIter.h>
43 :
44 : #define CheckVisIter() checkVisIter (__func__, __FILE__, __LINE__)
45 : #define CheckVisIter1(s) checkVisIter (__func__, __FILE__, __LINE__,s)
46 : #define CheckVisIterBase() checkVisIterBase (__func__, __FILE__, __LINE__)
47 :
48 :
49 : // For debugging; remove/comment-out when working
50 : //#include "VLAT.h"
51 : //#define Log(level, ...) \
52 : // {if (VlaData::loggingInitialized_p && level <= VlaData::logLevel_p) \
53 : // Logger::log (__VA_ARGS__);};
54 :
55 : using namespace casacore;
56 :
57 : namespace casa { //# NAMESPACE CASA - BEGIN
58 :
59 407 : VisBuffer::VisBuffer()
60 : : corrSorted_p(false),
61 : lastPointTableRow_p(0),
62 : This(this),
63 : twoWayConnection_p(false),
64 407 : visIter_p(static_cast<ROVisibilityIterator *>(0))
65 : {
66 407 : validate();
67 407 : oldMSId_p = -1;
68 :
69 814 : std::unique_ptr<casa::VisModelDataI> check(VisModelDataI::create());
70 407 : if (check)
71 407 : visModelData_p = check.release();
72 407 : }
73 :
74 827 : VisBuffer::VisBuffer(ROVisibilityIterator & iter)
75 : : This(this),
76 827 : visIter_p(&iter)
77 : {
78 827 : iter.attachVisBuffer(*this);
79 827 : twoWayConnection_p = true;
80 827 : oldMSId_p = -1;
81 827 : corrSorted_p = false;
82 :
83 1654 : std::unique_ptr<casa::VisModelDataI> check(VisModelDataI::create());
84 827 : if (check)
85 827 : visModelData_p = check.release();
86 827 : }
87 :
88 5898 : VisBuffer::VisBuffer(const VisBuffer & vb)
89 : : This(this),
90 5898 : visIter_p(static_cast<ROVisibilityIterator *>(0))
91 : {
92 5898 : corrSorted_p = false;
93 5898 : operator=(vb);
94 5898 : }
95 :
96 5898 : VisBuffer & VisBuffer::operator=(const VisBuffer & other)
97 : {
98 5898 : if (this != &other) {
99 5898 : assign(other);
100 5898 : oldMSId_p = -1;
101 : }
102 5898 : return *this;
103 : }
104 :
105 : VisBuffer &
106 6162 : VisBuffer::assign(const VisBuffer & other, Bool copy)
107 : {
108 6162 : if (not other.visModelData_p.null()) visModelData_p = other.visModelData_p->clone();
109 :
110 6162 : if (other.corrSorted_p) {
111 0 : throw(AipsError("Cannot assign a VisBuffer that has had correlations sorted!"));
112 : }
113 :
114 6162 : if (this != &other) {
115 12324 : if (visIter_p != static_cast<ROVisibilityIterator *>(0) &&
116 6162 : twoWayConnection_p &&
117 0 : visIter_p != other.getVisibilityIterator()) {
118 :
119 : // If this VB is attached to it's visibility iterator and the
120 : // assignment will result in association with a different VI then
121 : // detach it.
122 :
123 0 : visIter_p->detachVisBuffer(*this);
124 : }
125 :
126 6162 : visIter_p = other.getVisibilityIterator ();
127 6162 : other.copyMsInfo(oldMSId_p, msOK_p, newMS_p);
128 :
129 6162 : twoWayConnection_p = false;
130 :
131 6162 : if (visIter_p == static_cast<ROVisibilityIterator *>(0)) {
132 0 : validate();
133 0 : copyCache (other, true); // force copying
134 6162 : } else if (copy) {
135 5898 : copyCache (other, false); // copy only if there's something there
136 : } else {
137 264 : invalidate();
138 : }
139 :
140 : }
141 6162 : return *this;
142 : }
143 :
144 :
145 250 : void VisBuffer::copyCoordInfo(const VisBuffer& other, Bool force)
146 : {
147 : // Just do the nominally non-row-dep values
148 250 : cacheCopyNormal(arrayIdOK_p, other.arrayIdOK(), arrayId_p, other, &VisBuffer::arrayId, force);
149 250 : cacheCopyNormal(dataDescriptionIdOK_p, other.dataDescriptionIdOK(), dataDescriptionId_p, other, &VisBuffer::dataDescriptionId, force);
150 250 : cacheCopyNormal(fieldIdOK_p, other.fieldIdOK(), fieldId_p, other, &VisBuffer::fieldId, force);
151 250 : cacheCopyNormal(spectralWindowOK_p, other.spectralWindowOK(), spectralWindow_p, other,
152 : &VisBuffer::spectralWindow, force);
153 250 : cacheCopyNormal(nCorrOK_p, other.nCorrOK(), nCorr_p, other, &VisBuffer::nCorr, force);
154 250 : cacheCopyNormal(nChannelOK_p, other.nChannelOK(), nChannel_p, other, &VisBuffer::nChannel, force);
155 250 : cacheCopyArray(frequencyOK_p, other.frequencyOK(), frequency_p, other, &VisBuffer::frequency, force);
156 250 : }
157 :
158 : void
159 5898 : VisBuffer::copyCache (const VisBuffer & other, Bool force)
160 : {
161 : // Copies cache status from the other VisBuffer and if the status is true
162 : // then the cached values are copied over from the other VisBuffer as well.
163 :
164 : // Keep in order so that finding omitted ones will be easier
165 : // in the future
166 :
167 5898 : cacheCopyArray (antenna1OK_p, other.antenna1OK (), antenna1_p, other, & VisBuffer::antenna1, force);
168 5898 : cacheCopyArray (antenna2OK_p, other.antenna2OK (), antenna2_p, other, & VisBuffer::antenna2, force);
169 5898 : cacheCopyNormal (arrayIdOK_p, other.arrayIdOK (), arrayId_p, other, & VisBuffer::arrayId, force);
170 : ////cacheCopyArray (chanAveBoundsOK_p, other.chanAveBoundsOK (), chanAveBounds_p, other, & VisBuffer::chanAveBounds, force);
171 5898 : cacheCopyArray (channelOK_p, other.channelOK (), channel_p, other, & VisBuffer::channel, force);
172 5898 : cacheCopyArray (cjonesOK_p, other.cjonesOK (), cjones_p, other, & VisBuffer::CJones, force);
173 5898 : cacheCopyArray (correctedVisCubeOK_p, other.correctedVisCubeOK (),
174 5898 : correctedVisCube_p, other, & VisBuffer::correctedVisCube, force);
175 5898 : cacheCopyArray (correctedVisibilityOK_p, other.correctedVisibilityOK (),
176 5898 : correctedVisibility_p, other, & VisBuffer::correctedVisibility, force);
177 5898 : cacheCopyArray (corrTypeOK_p, other.corrTypeOK (), corrType_p, other, & VisBuffer::corrType, force);
178 5898 : cacheCopyNormal (dataDescriptionIdOK_p, other.dataDescriptionIdOK(), dataDescriptionId_p, other, & VisBuffer::dataDescriptionId, force);
179 5898 : cacheCopyArray (direction1OK_p, other.direction1OK (), direction1_p, other, & VisBuffer::direction1, force);
180 5898 : cacheCopyArray (direction2OK_p, other.direction2OK (), direction2_p, other, & VisBuffer::direction2, force);
181 5898 : cacheCopyArray (exposureOK_p, other.exposureOK (), exposure_p, other, & VisBuffer::exposure, force);
182 5898 : cacheCopyArray (feed1OK_p, other.feed1OK (), feed1_p, other, & VisBuffer::feed1, force);
183 5898 : cacheCopyArray (feed1_paOK_p, other.feed1_paOK (), feed1_pa_p, other, & VisBuffer::feed1_pa, force);
184 5898 : cacheCopyArray (feed2OK_p, other.feed2OK (), feed2_p, other, & VisBuffer::feed2, force);
185 5898 : cacheCopyArray (feed2_paOK_p, other.feed2_paOK (), feed2_pa_p, other, & VisBuffer::feed2_pa, force);
186 5898 : cacheCopyNormal (fieldIdOK_p, other.fieldIdOK (), fieldId_p, other, & VisBuffer::fieldId, force);
187 5898 : cacheCopyArray (flagOK_p, other.flagOK (), flag_p, other, & VisBuffer::flag, force);
188 5898 : cacheCopyArray (flagCategoryOK_p, other.flagCategoryOK (), flagCategory_p, other, & VisBuffer::flagCategory, force);
189 5898 : cacheCopyArray (flagCubeOK_p, other.flagCubeOK (), flagCube_p, other, & VisBuffer::flagCube, force);
190 5898 : cacheCopyArray (flagRowOK_p, other.flagRowOK (), flagRow_p, other, & VisBuffer::flagRow, force);
191 5898 : cacheCopyArray (floatDataCubeOK_p, other.floatDataCubeOK (), floatDataCube_p, other, & VisBuffer::floatDataCube, force);
192 5898 : cacheCopyArray (frequencyOK_p, other.frequencyOK (), frequency_p, other, & VisBuffer::frequency, force);
193 5898 : cacheCopyArray (imagingWeightOK_p, other.imagingWeightOK (), imagingWeight_p, other, & VisBuffer::imagingWeight, force);
194 : //cacheCopyArray (lsrFrequencyOK_p, other.lsrFrequencyOK (), lsrFrequency_p, other, & VisBuffer::lsrFrequency, force);
195 5898 : cacheCopyArray (modelVisCubeOK_p, other.modelVisCubeOK (), modelVisCube_p, other, & VisBuffer::modelVisCube, force);
196 5898 : cacheCopyArray (modelVisibilityOK_p, other.modelVisibilityOK (),
197 5898 : modelVisibility_p, other, & VisBuffer::modelVisibility, force);
198 5898 : cacheCopyNormal (nChannelOK_p, other.nChannelOK (), nChannel_p, other, & VisBuffer::nChannel, force);
199 5898 : cacheCopyNormal (nCorrOK_p, other.nCorrOK (), nCorr_p, other, & VisBuffer::nCorr, force);
200 : //cacheCopyNormal (nCatOK_p, other.nCatOK (), nCat_p, other, & VisBuffer::nCat, force);
201 5898 : cacheCopyNormal (nRowOK_p, other.nRowOK (), nRow_p, other, & VisBuffer::nRow, force);
202 5898 : cacheCopyArray (observationIdOK_p, other.observationIdOK (), observationId_p, other, & VisBuffer::observationId, force);
203 5898 : cacheCopyNormal (phaseCenterOK_p, other.phaseCenterOK (), phaseCenter_p, other, & VisBuffer::phaseCenter, force);
204 5898 : cacheCopyNormal (polFrameOK_p, other.polFrameOK (), polFrame_p, other, & VisBuffer::polFrame, force);
205 5898 : cacheCopyArray (processorIdOK_p, other.processorIdOK (), processorId_p, other, & VisBuffer::processorId, force);
206 5898 : cacheCopyArray (rowIdsOK_p, other.rowIdsOK (), rowIds_p, other, & VisBuffer::rowIds, force);
207 5898 : cacheCopyArray (scanOK_p, other.scanOK (), scan_p, other, & VisBuffer::scan, force);
208 5898 : cacheCopyArray (sigmaOK_p, other.sigmaOK (), sigma_p, other, & VisBuffer::sigma, force);
209 5898 : cacheCopyArray (sigmaMatOK_p, other.sigmaMatOK (), sigmaMat_p, other, & VisBuffer::sigmaMat, force);
210 5898 : cacheCopyNormal (spectralWindowOK_p, other.spectralWindowOK (), spectralWindow_p, other, & VisBuffer::spectralWindow, force);
211 5898 : cacheCopyArray (stateIdOK_p, other.stateIdOK (), stateId_p, other, & VisBuffer::stateId, force);
212 5898 : cacheCopyArray (timeOK_p, other.timeOK (), time_p, other, & VisBuffer::time, force);
213 5898 : cacheCopyArray (timeCentroidOK_p, other.timeCentroidOK (), timeCentroid_p, other, & VisBuffer::timeCentroid, force);
214 5898 : cacheCopyArray (timeIntervalOK_p, other.timeIntervalOK (), timeInterval_p, other, & VisBuffer::timeInterval, force);
215 5898 : cacheCopyArray (uvwOK_p, other.uvwOK (), uvw_p, other, & VisBuffer::uvw, force);
216 5898 : cacheCopyArray (uvwMatOK_p, other.uvwMatOK (), uvwMat_p, other, & VisBuffer::uvwMat, force);
217 5898 : cacheCopyArray (visCubeOK_p, other.visCubeOK (), visCube_p, other, & VisBuffer::visCube, force);
218 5898 : cacheCopyArray (visibilityOK_p, other.visibilityOK (), visibility_p, other, & VisBuffer::visibility, force);
219 5898 : cacheCopyArray (weightOK_p, other.weightOK (), weight_p, other, & VisBuffer::weight, force);
220 : ////cacheCopyArray (weightCubeOK_p, other.weightCubeOK (), weightCube_p, other, & VisBuffer::weightCube, force);
221 5898 : cacheCopyArray (weightMatOK_p, other.weightMatOK (), weightMat_p, other, & VisBuffer::weightMat, force);
222 5898 : cacheCopyArray (weightSpectrumOK_p, other.weightSpectrumOK (),
223 5898 : weightSpectrum_p, other, & VisBuffer::weightSpectrum, force);
224 :
225 5898 : }
226 :
227 7300 : VisBuffer::~VisBuffer()
228 : {
229 6724 : if (visIter_p != static_cast<ROVisibilityIterator *>(0) && twoWayConnection_p) {
230 827 : visIter_p->detachVisBuffer(*this);
231 : }
232 7300 : }
233 :
234 : VisBuffer &
235 0 : VisBuffer::operator-=(const VisBuffer & vb)
236 : {
237 : // check the shapes
238 0 : AlwaysAssert(nRow_p == vb.nRow(), AipsError);
239 0 : AlwaysAssert(nChannel_p == vb.nChannel(), AipsError);
240 0 : AlwaysAssert(nCorr_p == vb.nCorr(), AipsError);
241 : // make sure flag and flagRow are current
242 0 : flag();
243 0 : flagRow();
244 : // flagCategory?
245 :
246 : // do the subtraction, or'ing the flags
247 0 : Int nRows = nRow ();
248 0 : Int nChannels = nChannel ();
249 0 : for (Int row = 0; row < nRows; row++) {
250 0 : if (vb.flagRow()(row)) {
251 0 : flagRow_p(row) = true;
252 : }
253 0 : if (!flagRow_p(row)) {
254 0 : for (Int chn = 0; chn < nChannels; chn++) {
255 0 : if (vb.flag()(chn, row)) {
256 0 : flag_p(chn, row) = true;
257 : }
258 0 : if (!flag_p(chn, row)) {
259 0 : visibility_p(chn, row) -= vb.visibility()(chn, row);
260 : }
261 : }
262 : }
263 : }
264 0 : return *this;
265 : }
266 :
267 : void
268 0 : VisBuffer::attachToVisIter(ROVisibilityIterator & iter)
269 : {
270 0 : if (visIter_p != static_cast<ROVisibilityIterator *>(0) && twoWayConnection_p) {
271 0 : visIter_p->detachVisBuffer(*this);
272 : }
273 0 : visIter_p = &iter;
274 0 : iter.attachVisBuffer(*this);
275 0 : twoWayConnection_p = true;
276 0 : }
277 :
278 : void
279 0 : VisBuffer::detachFromVisIter ()
280 : {
281 0 : if (visIter_p != NULL) {
282 0 : visIter_p->detachVisBuffer(* this);
283 :
284 0 : visIter_p = NULL;
285 : }
286 0 : }
287 :
288 44033 : void VisBuffer::invalidate()
289 : {
290 :
291 44033 : setAllCacheStatuses (false);
292 44033 : lastPointTableRow_p = 0;
293 44033 : }
294 :
295 407 : void VisBuffer::validate()
296 : {
297 407 : setAllCacheStatuses (true);
298 407 : }
299 :
300 : Int
301 0 : VisBuffer::getOldMsId () const
302 : {
303 0 : return oldMSId_p;
304 : }
305 :
306 228 : String VisBuffer::msName(Bool stripPath) const{
307 456 : String name="";
308 228 : if(visIter_p != NULL){
309 228 : name=visIter_p->ms().antenna().tableName();
310 228 : name.erase(name.length()-8);
311 228 : if(stripPath){
312 0 : Path path(name);
313 0 : return path.baseName();
314 : }
315 :
316 : }
317 :
318 228 : return name;
319 : }
320 :
321 : ROVisibilityIterator *
322 6425 : VisBuffer::getVisibilityIterator () const
323 : {
324 6425 : return visIter_p;
325 : }
326 :
327 : Matrix<Float> &
328 0 : VisBuffer::imagingWeight ()
329 : {
330 0 : static_cast<const VisBuffer *> (this) -> imagingWeight ();
331 :
332 0 : return imagingWeight_p;
333 : }
334 :
335 : const Matrix<Float> &
336 263 : VisBuffer::imagingWeight () const
337 : {
338 263 : const VisImagingWeight & weightGenerator = getVisibilityIterator()->getImagingWeightGenerator ();
339 :
340 263 : return imagingWeight (weightGenerator);
341 : }
342 :
343 : const Matrix<Float> &
344 263 : VisBuffer::imagingWeight (const VisImagingWeight & weightGenerator) const
345 : {
346 263 : if (imagingWeightOK_p){
347 0 : return imagingWeight_p;
348 : }
349 :
350 263 : if (weightGenerator.getType () == "none") {
351 0 : throw (AipsError ("Programmer Error... imaging weights not set"));
352 : }
353 :
354 : // Extract data weights correctly [nchan,nrow]
355 : // NB: VB::weight() yields corr-,chan-indep row weights
356 526 : Matrix<Float> wtm;
357 263 : if (weightSpectrum().nelements()==0)
358 263 : wtm.reference(weight().reform(IPosition(2,1,nRow()))); // add degenerate chan axis
359 : else
360 0 : weightGenerator.unPolChanWeight(wtm,weightSpectrum());
361 :
362 526 : Matrix<Bool> flagmat = flag ();
363 263 : imagingWeight_p.resize (flagmat.shape ());
364 :
365 526 : Vector<Double> fvec;
366 526 : Matrix<Double> uvwmat;
367 :
368 263 : String type = weightGenerator.getType();
369 263 : if (weightGenerator.doFilter() || type == "uniform" || type == "radial"){
370 0 : fvec = frequency ();
371 0 : uvwmat = uvwMat ();
372 : }
373 :
374 263 : if (weightGenerator.getType () == "uniform") {
375 :
376 0 : weightGenerator.weightUniform (imagingWeight_p, flagmat, uvwmat, fvec, wtm, msId (), fieldId ());
377 :
378 263 : } else if (weightGenerator.getType () == "radial") {
379 :
380 0 : weightGenerator.weightRadial (imagingWeight_p, flagmat, uvwmat, fvec, wtm);
381 :
382 : } else {
383 :
384 263 : weightGenerator.weightNatural (imagingWeight_p, flagmat, wtm);
385 : }
386 :
387 263 : if (weightGenerator.doFilter ()) {
388 :
389 0 : weightGenerator.filter (imagingWeight_p, flagmat, uvwmat, fvec, wtm);
390 : }
391 :
392 263 : This->imagingWeightOK_p = true;
393 :
394 263 : return imagingWeight_p;
395 : }
396 :
397 :
398 : //const Matrix<Float> &
399 : //VisBuffer::imagingWeight () const
400 : //{
401 : // return imagingWeight ();
402 : //}
403 :
404 : Bool
405 0 : VisBuffer::newArrayId () const
406 : {
407 0 : CheckVisIter ();
408 0 : return visIter_p->newArrayId ();
409 : }
410 :
411 : Bool
412 0 : VisBuffer::newFieldId () const
413 : {
414 0 : CheckVisIter ();
415 0 : return visIter_p->newFieldId ();
416 : }
417 :
418 : Bool
419 0 : VisBuffer::newSpectralWindow () const
420 : {
421 0 : CheckVisIter ();
422 0 : return visIter_p->newSpectralWindow ();
423 : }
424 :
425 :
426 :
427 :
428 : void
429 44440 : VisBuffer::setAllCacheStatuses (bool status)
430 : {
431 44440 : antenna1OK_p = status;
432 44440 : antenna2OK_p = status;
433 44440 : arrayIdOK_p = status;
434 44440 : channelOK_p = status;
435 44440 : cjonesOK_p = status;
436 44440 : correctedVisCubeOK_p = status;
437 44440 : correctedVisibilityOK_p = status;
438 44440 : corrTypeOK_p = status;
439 44440 : dataDescriptionIdOK_p = status;
440 44440 : direction1OK_p = status;
441 44440 : firstDirection1OK_p=status;
442 44440 : direction2OK_p = status;
443 44440 : exposureOK_p = status;
444 44440 : feed1_paOK_p = status;
445 44440 : feed1OK_p = status;
446 44440 : feed2_paOK_p = status;
447 44440 : feed2OK_p = status;
448 44440 : fieldIdOK_p = status;
449 44440 : flagCubeOK_p = status;
450 44440 : flagOK_p = status;
451 44440 : flagRowOK_p = status;
452 44440 : flagCategoryOK_p = status;
453 44440 : floatDataCubeOK_p = status;
454 44440 : frequencyOK_p = status;
455 44440 : imagingWeightOK_p = status;
456 : ///////////lsrFrequencyOK_p = status;
457 44440 : modelVisCubeOK_p = status;
458 44440 : modelVisibilityOK_p = status;
459 44440 : msOK_p = status;
460 44440 : nChannelOK_p = status;
461 44440 : nCorrOK_p = status;
462 : // nCatOK_p = status;
463 44440 : nRowOK_p = status;
464 44440 : observationIdOK_p = status;
465 44440 : phaseCenterOK_p = status;
466 44440 : polFrameOK_p = status;
467 44440 : processorIdOK_p = status;
468 44440 : rowIdsOK_p = status;
469 44440 : scanOK_p = status;
470 44440 : sigmaMatOK_p = status;
471 44440 : sigmaOK_p = status;
472 44440 : spectralWindowOK_p = status;
473 44440 : stateIdOK_p = status;
474 44440 : timeCentroidOK_p = status;
475 44440 : timeIntervalOK_p = status;
476 44440 : timeOK_p = status;
477 44440 : uvwMatOK_p = status;
478 44440 : uvwOK_p = status;
479 44440 : visCubeOK_p = status;
480 44440 : visibilityOK_p = status;
481 44440 : weightMatOK_p = status;
482 44440 : weightOK_p = status;
483 44440 : weightSpectrumOK_p = status;
484 :
485 44440 : }
486 :
487 2616 : Cube<Complex>& VisBuffer::dataCube(const MS::PredefinedColumns whichcol)
488 : {
489 2616 : switch(whichcol){
490 2582 : case MS::DATA:
491 2582 : return This->visCube();
492 34 : case MS::MODEL_DATA:
493 34 : return This->modelVisCube();
494 0 : case MS::CORRECTED_DATA:
495 0 : return This->correctedVisCube();
496 0 : default:
497 0 : throw(AipsError(MS::columnName(whichcol) + " is not supported as a data Cube."));
498 : }
499 : }
500 :
501 0 : const Cube<Complex>& VisBuffer::dataCube(const MS::PredefinedColumns whichcol) const
502 : {
503 0 : switch(whichcol){
504 0 : case MS::DATA:
505 0 : return This->visCube();
506 0 : case MS::MODEL_DATA:
507 0 : return This->modelVisCube();
508 0 : case MS::CORRECTED_DATA:
509 0 : return This->correctedVisCube();
510 0 : default:
511 0 : throw(AipsError(MS::columnName(whichcol) + " is not supported as a data Cube."));
512 : }
513 : }
514 :
515 0 : void VisBuffer::freqAverage()
516 : {
517 0 : Matrix<CStokesVector> newVisibility(1, nRow());
518 0 : Matrix<Bool> newFlag(1, nRow());
519 0 : newFlag = true;
520 : Double newFrequency;
521 0 : newFrequency = 0;
522 0 : Int nfreq = 0;
523 0 : Int nChan = nChannel();
524 0 : for (Int row = 0; row < nRow(); row++) {
525 0 : if (!flagRow()(row)) {
526 0 : Int n = 0;
527 0 : for (Int chn = 0; chn < nChan; chn++) {
528 0 : if (!flag()(chn, row)) {
529 0 : newVisibility(0, row) += visibility()(chn, row);
530 0 : newFlag(0, row) = false;
531 0 : newFrequency += frequency()(chn);
532 0 : n++;
533 0 : nfreq++;
534 : }
535 : }
536 0 : if (n == 0) {
537 0 : flagRow()(row) = true;
538 : }
539 0 : if (n > 0) {
540 0 : newVisibility(0, row) *= 1.0f / n;
541 : }
542 : }
543 : }
544 : // Average frequency for this buffer (should really be row based)
545 0 : if (nfreq > 0) {
546 0 : newFrequency /= Double(nfreq);
547 : }
548 0 : nChannel_p = 1;
549 0 : flag_p.reference(newFlag);
550 0 : visibility_p.reference(newVisibility);
551 0 : frequency_p.resize(1);
552 0 : frequency_p(0) = newFrequency;
553 0 : }
554 :
555 4035 : void VisBuffer::freqAveCubes()
556 : {
557 : // TBD: Use weightSpec, if present, and update weightMat accordingly
558 : // TBD: Provide partial decimation option
559 :
560 : // Ensure visCube filled, at least
561 4035 : visCube();
562 :
563 : // Freq-averaged shape
564 8070 : IPosition csh = visCube().shape();
565 4035 : csh(1) = 1; // One channel in output
566 :
567 8070 : Cube<Complex> newVisCube(csh);
568 4035 : newVisCube = Complex(0.0);
569 4035 : Matrix<Bool> newFlag(1, nRow());
570 4035 : newFlag = true;
571 : Double newFrequency;
572 4035 : newFrequency = 0;
573 4035 : Int nfreq = 0;
574 4035 : Int nChan = nChannel();
575 4035 : Int nCor = nCorr();
576 179520 : for (Int row = 0; row < nRow(); row++) {
577 175485 : if (!flagRow()(row)) {
578 175485 : Int n = 0;
579 1837580 : for (Int chn = 0; chn < nChan; chn++) {
580 1662095 : if (!flag()(chn, row)) {
581 1643961 : newFlag(0, row) = false;
582 1643961 : newFrequency += frequency()(chn);
583 7434062 : for (Int cor = 0; cor < nCor; cor++) {
584 5790101 : newVisCube(cor, 0, row) += visCube()(cor, chn, row);
585 : }
586 1643961 : n++;
587 1643961 : nfreq++;
588 :
589 1643961 : if (row == -1)
590 0 : cout << "V: "
591 0 : << chn << " " << n << " "
592 0 : << visCube()(0, chn, row) << " "
593 0 : << newVisCube(0, 0, row) << " "
594 0 : << endl;
595 :
596 :
597 : }
598 : }
599 175485 : if (n == 0) {
600 0 : flagRow()(row) = true;
601 : }
602 175485 : if (n > 0) {
603 350970 : Matrix<Complex> nVC;
604 175485 : nVC.reference(newVisCube.xyPlane(row));
605 175485 : nVC *= Complex(1.0f / n);
606 :
607 175485 : if (row == -1) {
608 0 : cout << "V:-----> " << n << " " << newVisCube(0, 0, row) << endl;
609 : }
610 :
611 :
612 : }
613 : }
614 : }
615 4035 : visCube_p.reference(newVisCube);
616 :
617 : // Now do model, if present
618 4035 : if (modelVisCubeOK_p) {
619 :
620 8070 : Cube<Complex> newModelVisCube(csh);
621 4035 : newModelVisCube = Complex(0.0);
622 179520 : for (Int row = 0; row < nRow(); row++) {
623 175485 : if (!flagRow()(row)) {
624 175485 : Int n = 0;
625 1837580 : for (Int chn = 0; chn < nChan; chn++) {
626 1662095 : if (!flag()(chn, row)) {
627 :
628 1643961 : n++;
629 7434062 : for (Int cor = 0; cor < nCor; cor++) {
630 5790101 : newModelVisCube(cor, 0, row) += modelVisCube()(cor, chn, row);
631 : }
632 :
633 1643961 : if (row == -1)
634 0 : cout << "M: "
635 0 : << chn << " " << n << " "
636 0 : << modelVisCube()(0, chn, row) << " "
637 0 : << newModelVisCube(0, 0, row) << " "
638 0 : << endl;
639 :
640 : }
641 : }
642 175485 : if (n == 0) {
643 0 : flagRow()(row) = true;
644 : }
645 175485 : if (n > 0) {
646 350970 : Matrix<Complex> nMVC;
647 175485 : nMVC.reference(newModelVisCube.xyPlane(row));
648 175485 : nMVC *= Complex(1.0f / n);
649 :
650 175485 : if (row == -1) {
651 0 : cout << "M:-----> " << n << " " << newModelVisCube(0, 0, row) << endl;
652 : }
653 : }
654 : }
655 : }
656 4035 : modelVisCube_p.reference(newModelVisCube);
657 : }
658 :
659 : // Use averaged flags
660 4035 : flag_p.reference(newFlag);
661 :
662 : // Average frequency for this buffer
663 : // (Strictly, this should really be row based, but doing this
664 : // average here suggests frequency precision isn't so important)
665 4035 : if (nfreq > 0) {
666 4035 : newFrequency /= Double(nfreq);
667 : }
668 4035 : nChannel_p = 1;
669 4035 : frequency_p.resize(1);
670 4035 : frequency_p(0) = newFrequency;
671 :
672 4035 : }
673 :
674 0 : void VisBuffer::formStokes()
675 : {
676 :
677 : // We must form the weights and flags correctly
678 0 : formStokesWeightandFlag();
679 :
680 : // Now do whatever data is present
681 0 : if (visCubeOK_p) {
682 0 : formStokes(visCube_p);
683 : }
684 :
685 0 : if (modelVisCubeOK_p) {
686 0 : formStokes(modelVisCube_p);
687 : }
688 :
689 0 : if (correctedVisCubeOK_p) {
690 0 : formStokes(correctedVisCube_p);
691 : }
692 :
693 0 : if (floatDataCubeOK_p) {
694 0 : formStokes(floatDataCube_p);
695 : }
696 0 : }
697 :
698 : void
699 1841 : VisBuffer::lsrFrequency(const Int & spw, Vector<Double>& freq, Bool & convert,
700 : const Bool ignoreconv) const
701 : {
702 1841 : CheckVisIter ();
703 1841 : visIter_p->lsrFrequency(spw, freq, convert, ignoreconv);
704 1841 : }
705 :
706 :
707 0 : void VisBuffer::formStokesWeightandFlag()
708 : {
709 :
710 : // Ensure corrType, weightMat and flagCube are filled
711 0 : corrType();
712 0 : weightMat();
713 0 : flagCube();
714 :
715 0 : switch (nCorr()) {
716 0 : case 4: {
717 :
718 0 : Slice all = Slice();
719 0 : Slice pp(0, 1, 1), pq(1, 1, 1), qp(2, 1, 1), qq(3, 1, 1);
720 0 : Slice a(0, 1, 1), b(1, 1, 1), c(2, 1, 1), d(3, 1, 1);
721 :
722 : // Sort for linears
723 0 : if (polFrame() == MSIter::Linear) {
724 0 : d = Slice(1, 1, 1); // Q
725 0 : b = Slice(2, 1, 1); // U
726 0 : c = Slice(3, 1, 1); // V
727 : }
728 :
729 0 : Matrix<Float> newWtMat(weightMat_p.shape());
730 0 : newWtMat(a, all) = newWtMat(d, all) = (weightMat_p(pp, all) + weightMat_p(qq, all));
731 0 : newWtMat(b, all) = newWtMat(c, all) = (weightMat_p(pq, all) + weightMat_p(qp, all));
732 0 : weightMat_p.reference(newWtMat);
733 :
734 0 : Cube<Bool> newFlagCube(flagCube_p.shape());
735 0 : newFlagCube(a, all, all) = newFlagCube(d, all, all) = (flagCube_p(pp, all, all) | flagCube_p(qq, all, all));
736 0 : newFlagCube(b, all, all) = newFlagCube(c, all, all) = (flagCube_p(pq, all, all) | flagCube_p(qp, all, all));
737 0 : flagCube_p.reference(newFlagCube);
738 :
739 0 : corrType_p(0) = Stokes::I;
740 0 : corrType_p(1) = Stokes::Q;
741 0 : corrType_p(2) = Stokes::U;
742 0 : corrType_p(3) = Stokes::V;
743 :
744 0 : break;
745 : }
746 0 : case 2: {
747 : // parallel hands only
748 0 : Slice all = Slice();
749 0 : Slice pp(0, 1, 1), qq(1, 1, 1);
750 0 : Slice a(0, 1, 1), d(1, 1, 1);
751 :
752 0 : Matrix<Float> newWtMat(weightMat_p.shape());
753 0 : newWtMat(a, all) = newWtMat(d, all) = weightMat_p(pp, all) + weightMat_p(qq, all);
754 0 : weightMat_p.reference(newWtMat);
755 :
756 0 : Cube<Bool> newFlagCube(flagCube_p.shape());
757 0 : newFlagCube(a, all, all) = newFlagCube(d, all, all) = flagCube_p(pp, all, all) | flagCube_p(qq, all, all);
758 0 : flagCube_p.reference(newFlagCube);
759 :
760 0 : corrType_p(0) = Stokes::I;
761 0 : corrType_p(1) = ((polFrame() == MSIter::Circular) ? Stokes::V : Stokes::Q);
762 :
763 0 : break;
764 : }
765 0 : case 1: {
766 :
767 : // Just need to re-label as I
768 0 : corrType_p(0) = Stokes::I;
769 :
770 : }
771 0 : default: {
772 0 : cout << "Insufficient correlations to form Stokes" << endl;
773 0 : break;
774 : }
775 : }
776 :
777 0 : }
778 :
779 :
780 :
781 0 : void VisBuffer::formStokes(Cube<Complex>& vis)
782 : {
783 :
784 0 : Cube<Complex> newvis(vis.shape());
785 0 : newvis.set(0.0);
786 0 : Slice all = Slice();
787 :
788 0 : switch (nCorr()) {
789 0 : case 4: {
790 :
791 0 : Slice pp(0, 1, 1), pq(1, 1, 1), qp(2, 1, 1), qq(3, 1, 1);
792 0 : Slice a(0, 1, 1), b(1, 1, 1), c(2, 1, 1), d(3, 1, 1);
793 :
794 0 : if (polFrame() == MSIter::Linear) {
795 0 : d = Slice(1, 1, 1); // Q
796 0 : b = Slice(2, 1, 1); // U
797 0 : c = Slice(3, 1, 1); // V
798 : }
799 :
800 0 : newvis(a, all, all) = (vis(pp, all, all) + vis(qq, all, all)); // I / I
801 0 : newvis(d, all, all) = (vis(pp, all, all) - vis(qq, all, all)); // V / Q
802 :
803 0 : newvis(b, all, all) = (vis(pq, all, all) + vis(qp, all, all)); // Q / U
804 0 : newvis(c, all, all) = (vis(pq, all, all) - vis(qp, all, all)) / Complex(0.0, 1.0); // U / V
805 0 : newvis /= Complex(2.0);
806 :
807 0 : vis.reference(newvis);
808 :
809 0 : break;
810 : }
811 0 : case 2: {
812 : // parallel hands only
813 0 : Slice pp(0, 1, 1), qq(1, 1, 1);
814 0 : Slice a(0, 1, 1), d(1, 1, 1);
815 :
816 0 : newvis(a, all, all) = (vis(pp, all, all) + vis(qq, all, all)); // I / I
817 0 : newvis(d, all, all) = (vis(pp, all, all) - vis(qq, all, all)); // V / Q
818 0 : newvis /= Complex(2.0);
819 :
820 0 : vis.reference(newvis);
821 :
822 0 : break;
823 : }
824 0 : case 1: {
825 : // need do nothing for single correlation case
826 0 : break;
827 : }
828 0 : default: {
829 0 : cout << "Insufficient correlations to form Stokes" << endl;
830 0 : break;
831 : }
832 : }
833 0 : }
834 :
835 0 : void VisBuffer::formStokes(Cube<Float>& fcube)
836 : {
837 0 : Cube<Float> newfcube(fcube.shape());
838 0 : newfcube.set(0.0);
839 0 : Slice all = Slice();
840 :
841 0 : switch (nCorr()) {
842 0 : case 4: {
843 : throw(AipsError(
844 0 : "Forming all 4 Stokes parameters out of FLOAT_DATA is not supported."));
845 :
846 : Slice pp(0, 1, 1), pq(1, 1, 1), qp(2, 1, 1), qq(3, 1, 1);
847 : Slice a(0, 1, 1), b(1, 1, 1), c(2, 1, 1), d(3, 1, 1);
848 :
849 : if (polFrame() == MSIter::Linear) {
850 : d = Slice(1, 1, 1); // Q
851 : b = Slice(2, 1, 1); // U
852 : c = Slice(3, 1, 1); // V
853 : }
854 :
855 : newfcube(a, all, all) = (fcube(pp, all, all) + fcube(qq, all, all)); // I / I
856 : newfcube(d, all, all) = (fcube(pp, all, all) - fcube(qq, all, all)); // V / Q
857 :
858 : newfcube(b, all, all) = (fcube(pq, all, all) + fcube(qp, all, all)); // Q / U
859 :
860 : //// U / V
861 : // This clearly isn't going to work. AFAICT it is impossible to
862 : // simultaneously measure all 4 polarizations with the same feed, so
863 : // FLOAT_DATA shouldn't come with 4 polarizations.
864 : //newfcube(c,all,all)=(fcube(pq,all,all)-fcube(qp,all,all))/Complex(0.0,1.0);
865 : newfcube(c, all, all) = (fcube(pq, all, all) - fcube(qp, all, all));
866 :
867 : // The cast is necessary to stop the compiler from promoting it to a Complex.
868 : newfcube *= static_cast<Float>(0.5);
869 :
870 : fcube.reference(newfcube);
871 :
872 : break;
873 : }
874 0 : case 2: {
875 : // parallel hands only
876 0 : Slice pp(0, 1, 1), qq(1, 1, 1);
877 0 : Slice a(0, 1, 1), d(1, 1, 1);
878 :
879 0 : newfcube(a, all, all) = (fcube(pp, all, all) + fcube(qq, all, all)); // I / I
880 0 : newfcube(d, all, all) = (fcube(pp, all, all) - fcube(qq, all, all)); // V / Q
881 :
882 : // The cast is necessary to stop the compiler from promoting it to a Complex.
883 0 : newfcube *= static_cast<Float>(0.5);
884 :
885 0 : fcube.reference(newfcube);
886 :
887 0 : break;
888 : }
889 0 : case 1: {
890 : // need do nothing for single correlation case
891 0 : break;
892 : }
893 0 : default: {
894 0 : cout << "Insufficient correlations to form Stokes" << endl;
895 0 : break;
896 : }
897 : }
898 0 : }
899 :
900 12 : void VisBuffer::channelAve(const Matrix<Int>& chanavebounds,Bool calmode)
901 : {
902 :
903 : // TBD:
904 : // a. nChanAve examination
905 : // b. calmode=T DONE
906 : // c. doWtSp clauses and rowWtFac (not needed?)
907 : // d. divide-by-zero safety DONE
908 : // e. no-averaging case sigma<->weight
909 :
910 :
911 : // Only do something if there is something to do
912 12 : if ( chanavebounds.nelements()>0 ) {
913 :
914 : // refer to the supplied chanavebounds
915 12 : chanAveBounds_p.reference(chanavebounds);
916 :
917 : // Examine chanAveBounds_p for consistency, detect nChanAve
918 : // TBD: improve actual examination...
919 12 : Int nChanAve = abs(chanAveBounds_p(0,1)-chanAveBounds_p(0,0))+1;
920 12 : nChanAve = max(1,nChanAve); // ensure >0
921 :
922 12 : Int nChanOut(chanAveBounds_p.nrow());
923 24 : Vector<Int> chans(channel()); // Ensure channels pre-filled
924 :
925 : // Row weight and sigma currently refer to the number of channels in the
926 : // MS, regardless of any selection.
927 : //Int nAllChan0 = visIter_p->msColumns().spectralWindow().numChan()(spectralWindow());
928 :
929 12 : const Bool doWtSp(visIter_p->existsWeightSpectrum());
930 :
931 : // Apply averaging to whatever data is present
932 12 : if (visCubeOK_p) chanAveVisCube(visCube(),nChanOut);
933 12 : if (modelVisCubeOK_p) chanAveVisCube(modelVisCube(),nChanOut);
934 12 : if (correctedVisCubeOK_p) chanAveVisCube(correctedVisCube(),nChanOut);
935 12 : if (floatDataCubeOK_p) chanAveVisCube(floatDataCube(), nChanOut);
936 :
937 12 : uInt nCor = nCorr();
938 12 : uInt nrows = nRow();
939 24 : Matrix<Float> rowWtFac(nCor, nrows);
940 :
941 12 : uInt nch = flagCube().shape()(1); // # of selected channels
942 : //Double selChanFac;
943 :
944 12 : if(doWtSp){ // Get the total weight spectrum
945 0 : const Cube<Float>& wtsp(weightSpectrum()); // while it is unaveraged.
946 :
947 0 : for(uInt row = 0; row < nrows; ++row){
948 0 : for(uInt icor = 0; icor < nCor; ++icor){
949 0 : rowWtFac(icor, row) = 0.0;
950 0 : for(uInt ichan = 0; ichan < nch; ++ichan)
951 : // Presumably the input row weight was set without taking flagging
952 : // into account.
953 0 : rowWtFac(icor, row) += wtsp(icor, ichan, row);
954 : }
955 : }
956 : }
957 : else
958 12 : rowWtFac = 1.0;
959 :
960 : // The false makes it leave weightSpectrum() averaged if it is present.
961 12 : if(flagCubeOK_p) chanAveFlagCube(flagCube(), nChanOut, false);
962 :
963 12 : if(flagCategoryOK_p)
964 0 : chanAveFlagCategory(flagCategory(), nChanOut);
965 :
966 : // Collapse the frequency values themselves, and count the number of
967 : // selected channels.
968 : // TBD: move this up to bounds calculation loop?
969 24 : Vector<Double> newFreq(nChanOut,0.0);
970 24 : Vector<Int> newChan(nChanOut,0);
971 12 : frequency(); // Ensure frequencies pre-filled
972 12 : Int nChan0(chans.nelements());
973 12 : Int ichan=0;
974 12 : Int totn = 0;
975 60 : for(Int ochan = 0; ochan < nChanOut; ++ochan){
976 48 : Int n = 0;
977 :
978 96 : while(chans[ichan] >= chanAveBounds_p(ochan, 0) &&
979 96 : chans[ichan] <= chanAveBounds_p(ochan, 1) &&
980 : ichan < nChan0){
981 48 : ++n;
982 48 : newFreq[ochan] += (frequency()[ichan] - newFreq[ochan]) / n;
983 48 : newChan[ochan] += chans[ichan];
984 48 : ichan++;
985 : }
986 48 : if (n>0) {
987 48 : newChan[ochan] /= n;
988 48 : totn += n;
989 : }
990 : }
991 :
992 : // Install the new values
993 12 : frequency().reference(newFreq);
994 12 : channel().reference(newChan);
995 12 : nChannel()=nChanOut;
996 :
997 12 : if(doWtSp){
998 : // chanAccCube(weightSpectrum(), nChanOut); already done.
999 0 : const Cube<Float>& wtsp(weightSpectrum());
1000 :
1001 0 : for(uInt row = 0; row < nrows; ++row){
1002 0 : for(uInt icor = 0; icor < nCor; ++icor){
1003 0 : Float totwtsp = rowWtFac(icor, row);
1004 :
1005 0 : rowWtFac(icor, row) = 0.0;
1006 0 : for(Int ochan = 0; ochan < nChanOut; ++ochan){
1007 0 : Float oswt = wtsp(icor, ochan, row); // output spectral
1008 : // weight
1009 0 : if(oswt > 0.0)
1010 0 : rowWtFac(icor, row) += oswt;
1011 : else
1012 0 : flagCube()(icor, ochan, row) = true;
1013 : }
1014 0 : if(totwtsp > 0.0)
1015 0 : rowWtFac(icor, row) /= totwtsp;
1016 : }
1017 : }
1018 : }
1019 : // This is slightly fudgy because it ignores the unselected part of
1020 : // weightSpectrum. Unfortunately now that selection is applied to
1021 : // weightSpectrum, we can't get at the unselected channel weights.
1022 : //Float selChanFac = static_cast<Float>(totn) / nAllChan0;
1023 :
1024 6456 : for(uInt row = 0; row < nrows; ++row){
1025 32220 : for(uInt icor = 0; icor < nCor; ++icor){
1026 : // Float rwfcr = rowWtFac(icor, row);
1027 :
1028 25776 : if (calmode) {
1029 :
1030 : // Downstream processing in calibration will use weightMat only
1031 :
1032 : // Just magnify by channal averaging factor
1033 0 : weightMat()(icor, row) *= nChanAve;
1034 :
1035 : /*
1036 : if(totn < nAllChan0)
1037 : weightMat()(icor, row) *= selChanFac * rwfcr;
1038 : if(rwfcr > 0.0) // Unlike WEIGHT, SIGMA is for a single chan.
1039 : sigmaMat()(icor, row) /= sqrt(nch * rwfcr / nChanOut);
1040 : */
1041 : }
1042 : else {
1043 :
1044 :
1045 : // selectively update sigma and weight info according to which
1046 : // columns were processed
1047 :
1048 25776 : if (visCubeOK_p || floatDataCubeOK_p) {
1049 : // update sigmaMat by inverse sqrt of number of channels averaged
1050 25776 : sigmaMat()(icor,row) /= sqrt(Double(nChanAve)); // nChanAve>0 already ensured
1051 :
1052 25776 : if (!correctedVisCubeOK_p) {
1053 : // calc weightMat from sigmaMat
1054 0 : Float &s= sigmaMat()(icor,row);
1055 0 : weightMat()(icor,row)= (s>0.0 ? 1./square(s) : FLT_EPSILON);
1056 : }
1057 : }
1058 :
1059 25776 : if (correctedVisCubeOK_p || modelVisCubeOK_p) {
1060 : // update weightMat
1061 25776 : weightMat()(icor,row)*=Double(nChanAve);
1062 25776 : if (!visCubeOK_p) {
1063 : // calc sigmaMat from weightMat
1064 0 : Float &w = weightMat()(icor,row);
1065 0 : sigmaMat()(icor,row)= ( w>0.0 ? 1./sqrt(w) : FLT_MAX );
1066 : }
1067 : }
1068 :
1069 : } // calmode
1070 :
1071 : }
1072 : }
1073 : } // chanavebounds
1074 :
1075 : // do we need to do something here that is datacolumn-specific to get sigma/weight reconciled?
1076 :
1077 12 : }
1078 :
1079 12 : void VisBuffer::chanAveFlagCube(Cube<Bool>& flagcube, Int nChanOut,
1080 : const Bool restoreWeightSpectrum)
1081 : {
1082 12 : IPosition csh(flagcube.shape());
1083 12 : Int nChan0 = csh(1);
1084 12 : csh(1) = nChanOut;
1085 :
1086 12 : if(nChan0 < nChanOut)
1087 : // It's possible that data has already been averaged. I could try
1088 : // refilling data if I knew which column to use, but I don't.
1089 : // Chuck it to the caller.
1090 0 : throw(AipsError("Can't average " + String(nChan0) + " channels to " +
1091 0 : String(nChanOut) + " channels!"));
1092 :
1093 12 : Vector<Int>& chans(channel());
1094 12 : if(nChan0 == nChanOut && chans.nelements() > 0 && chans[0] == 0)
1095 0 : return; // No-op.
1096 :
1097 24 : Cube<Bool> newFlag(csh);
1098 12 : newFlag = true;
1099 :
1100 12 : const Bool doWtSp(visIter_p->existsWeightSpectrum());
1101 :
1102 : // weightSpectrum could be either averaged or unaveraged. Choose averaged.
1103 12 : if(doWtSp && weightSpectrum().shape()[1] > nChanOut)
1104 0 : chanAccCube(weightSpectrum(), nChanOut);
1105 :
1106 12 : Int nCor = nCorr();
1107 : Int ichan;
1108 6456 : for (Int row=0; row<nRow(); row++) {
1109 6444 : if (!flagRow()(row)) {
1110 6444 : ichan = 0;
1111 32220 : for (Int ochan = 0; ochan < nChanOut; ++ochan) {
1112 51552 : while (chans[ichan] >= chanAveBounds_p(ochan, 0) &&
1113 51552 : chans[ichan] <= chanAveBounds_p(ochan, 1) &&
1114 : ichan < nChan0) {
1115 128880 : for(Int icor = 0; icor < nCor; ++icor){
1116 103104 : if(!flagcube(icor, ichan, row))
1117 100692 : newFlag(icor, ochan, row) = false;
1118 : }
1119 25776 : ++ichan;
1120 : }
1121 : }
1122 : }
1123 : }
1124 : // Use averaged flags
1125 12 : flagcube.reference(newFlag);
1126 :
1127 12 : if(doWtSp && restoreWeightSpectrum)
1128 0 : fillWeightSpectrum();
1129 : }
1130 :
1131 0 : void VisBuffer::chanAveFlagCategory(Array<Bool>& flagcat, const Int nChanOut)
1132 : {
1133 0 : if(flagcat.nelements() == 0) // Averaging a degenerate Array is fast,
1134 0 : return; // and avoids the conformance check below.
1135 :
1136 0 : IPosition csh(flagcat.shape());
1137 0 : Int nChan0 = csh(1);
1138 0 : csh(1) = nChanOut;
1139 :
1140 0 : if(nChan0 < nChanOut)
1141 : // It's possible that data has already been averaged. I could try
1142 : // refilling data if I knew which column to use, but I don't.
1143 : // Chuck it to the caller.
1144 0 : throw(AipsError("Can't average " + String(nChan0) + " channels to " +
1145 0 : String(nChanOut) + " channels!"));
1146 :
1147 0 : Vector<Int>& chans(channel());
1148 0 : if(nChan0 == nChanOut && chans.nelements() > 0 && chans[0] == 0)
1149 0 : return; // No-op.
1150 :
1151 0 : Array<Bool> newFC(csh);
1152 0 : newFC = true;
1153 :
1154 0 : Int nCor = nCorr();
1155 0 : Int nCat = csh(2);
1156 : Int ichan;
1157 0 : for(Int row = 0; row < nRow(); ++row){
1158 0 : ichan = 0;
1159 0 : for(Int ochan = 0; ochan < nChanOut; ++ochan){
1160 0 : while(chans[ichan] >= chanAveBounds_p(ochan, 0) &&
1161 0 : chans[ichan] <= chanAveBounds_p(ochan, 1) &&
1162 : ichan < nChan0) {
1163 0 : for(Int icor = 0; icor < nCor; ++icor){
1164 0 : for(Int icat = 0; icat < nCat; ++icat){
1165 0 : if(!flagcat(IPosition(4, icor, ichan, icat, row)))
1166 0 : newFC(IPosition(4, icor, ochan, icat, row)) = false;
1167 : }
1168 : }
1169 0 : ++ichan;
1170 : }
1171 : }
1172 : }
1173 : // Use averaged flag Categories
1174 0 : flagcat.reference(newFC);
1175 : }
1176 :
1177 : // Sort correlations: (PP,QQ,PQ,QP) -> (PP,PQ,QP,QQ)
1178 8827 : void VisBuffer::sortCorr()
1179 : {
1180 :
1181 : // This method is for temporarily sorting the correlations
1182 : // into canonical order if the MS data is out-of-order
1183 : // NB: Always sorts the weightMat()
1184 : // NB: Only works on the visCube-style data
1185 : // NB: It only sorts the data columns which are already present
1186 : // (so make sure the ones you need are already read!)
1187 : // NB: It is the user's responsibility to run unSortCorr
1188 : // after using the sorted data to put it back in order
1189 : // NB: corrType_p is NOT changed to match the sorted
1190 : // correlations (it is expected that this sort is
1191 : // temporary, and that we will run unSortCorr
1192 : // NB: This method does nothing if no sort required
1193 :
1194 : // If nominal order is non-canonical (only for nCorr=4)
1195 : // and data not yet sorted
1196 8827 : if (nonCanonCorr() && !corrSorted_p) {
1197 :
1198 :
1199 : // First, do weightMat
1200 : {
1201 0 : weightMat(); // (ensures it is filled)
1202 :
1203 0 : Vector<Float> wtmp(nRow());
1204 0 : Vector<Float> w1, w2, w3;
1205 0 : IPosition wblc(2, 0, 0), wtrc(2, 0, nRow() - 1), vec(1, nRow());
1206 :
1207 0 : wblc(0) = wtrc(0) = 1;
1208 0 : w1.reference(weightMat_p(wblc, wtrc).reform(vec));
1209 0 : wblc(0) = wtrc(0) = 2;
1210 0 : w2.reference(weightMat_p(wblc, wtrc).reform(vec));
1211 0 : wblc(0) = wtrc(0) = 3;
1212 0 : w3.reference(weightMat_p(wblc, wtrc).reform(vec));
1213 0 : wtmp = w1;
1214 0 : w1 = w2;
1215 0 : w2 = w3;
1216 0 : w3 = wtmp;
1217 : }
1218 :
1219 : // Now do data:
1220 :
1221 : // Work space, references, coords
1222 0 : Matrix<Complex> tmp(nChannel(), nRow());
1223 0 : Matrix<Complex> p1, p2, p3;
1224 0 : IPosition blc(3, 0, 0, 0), trc(3, 0, nChannel() - 1, nRow() - 1), mat(2, nChannel(), nRow());
1225 :
1226 : // Do visCube if present
1227 0 : if (visCubeOK_p && visCube_p.nelements() > 0) {
1228 0 : blc(0) = trc(0) = 1;
1229 0 : p1.reference(visCube_p(blc, trc).reform(mat));
1230 0 : blc(0) = trc(0) = 2;
1231 0 : p2.reference(visCube_p(blc, trc).reform(mat));
1232 0 : blc(0) = trc(0) = 3;
1233 0 : p3.reference(visCube_p(blc, trc).reform(mat));
1234 0 : tmp = p1;
1235 0 : p1 = p2;
1236 0 : p2 = p3;
1237 0 : p3 = tmp;
1238 : }
1239 : // Do modelVisCube if present
1240 0 : if (modelVisCubeOK_p && modelVisCube_p.nelements() > 0) {
1241 0 : blc(0) = trc(0) = 1;
1242 0 : p1.reference(modelVisCube_p(blc, trc).reform(mat));
1243 0 : blc(0) = trc(0) = 2;
1244 0 : p2.reference(modelVisCube_p(blc, trc).reform(mat));
1245 0 : blc(0) = trc(0) = 3;
1246 0 : p3.reference(modelVisCube_p(blc, trc).reform(mat));
1247 0 : tmp = p1;
1248 0 : p1 = p2;
1249 0 : p2 = p3;
1250 0 : p3 = tmp;
1251 : }
1252 : // Do correctedVisCube if present
1253 0 : if (correctedVisCubeOK_p && correctedVisCube_p.nelements() > 0) {
1254 0 : blc(0) = trc(0) = 1;
1255 0 : p1.reference(correctedVisCube_p(blc, trc).reform(mat));
1256 0 : blc(0) = trc(0) = 2;
1257 0 : p2.reference(correctedVisCube_p(blc, trc).reform(mat));
1258 0 : blc(0) = trc(0) = 3;
1259 0 : p3.reference(correctedVisCube_p(blc, trc).reform(mat));
1260 0 : tmp = p1;
1261 0 : p1 = p2;
1262 0 : p2 = p3;
1263 0 : p3 = tmp;
1264 : }
1265 : // Do floatDataCube if present
1266 0 : if (floatDataCubeOK_p && floatDataCube_p.nelements() > 0) {
1267 0 : Matrix<Float> tmp(nChannel(), nRow());
1268 0 : Matrix<Float> p1, p2, p3;
1269 :
1270 0 : blc(0) = trc(0) = 1;
1271 0 : p1.reference(floatDataCube_p(blc, trc).reform(mat));
1272 0 : blc(0) = trc(0) = 2;
1273 0 : p2.reference(floatDataCube_p(blc, trc).reform(mat));
1274 0 : blc(0) = trc(0) = 3;
1275 0 : p3.reference(floatDataCube_p(blc, trc).reform(mat));
1276 0 : tmp = p1;
1277 0 : p1 = p2;
1278 0 : p2 = p3;
1279 0 : p3 = tmp;
1280 : }
1281 :
1282 : // Data is now sorted into canonical order
1283 0 : corrSorted_p = true;
1284 : }
1285 :
1286 8827 : }
1287 :
1288 : // Unsort correlations: (PP,PQ,QP,QQ) -> (PP,QQ,PQ,QP)
1289 4668 : void VisBuffer::unSortCorr()
1290 : {
1291 : // This method is for restoring the non-canonical correlation
1292 : // sort order so that the data matches the order indicated
1293 : // by corrType()
1294 : // NB: Always unsorts the weightMat()
1295 : // NB: Only works on the visCube-style data
1296 : // NB: It only unsorts the data columns which are already present
1297 : // (so make sure sortCorr sorted the ones you needed!)
1298 : // NB: This method is a no-op if no sort required, or if
1299 : // sortCorr hadn't been run since the last unSortCorr
1300 :
1301 : // If nominal order is non-canonical (only for nCorr=4)
1302 : // and if data has been sorted
1303 4668 : if (nonCanonCorr() && corrSorted_p) {
1304 :
1305 : // First, do weights
1306 : {
1307 0 : Vector<Float> wtmp(nRow());
1308 0 : Vector<Float> w1, w2, w3;
1309 0 : IPosition wblc(1, 0, 0), wtrc(3, 0, nRow() - 1), vec(1, nRow());
1310 :
1311 0 : wblc(0) = wtrc(0) = 1;
1312 0 : w1.reference(weightMat_p(wblc, wtrc).reform(vec));
1313 0 : wblc(0) = wtrc(0) = 2;
1314 0 : w2.reference(weightMat_p(wblc, wtrc).reform(vec));
1315 0 : wblc(0) = wtrc(0) = 3;
1316 0 : w3.reference(weightMat_p(wblc, wtrc).reform(vec));
1317 0 : wtmp = w3;
1318 0 : w3 = w2;
1319 0 : w2 = w1;
1320 0 : w1 = wtmp;
1321 : }
1322 : // Now do data:
1323 :
1324 : // Work space, references, coords
1325 0 : Matrix<Complex> tmp(nChannel(), nRow());
1326 0 : Matrix<Complex> p1, p2, p3;
1327 0 : IPosition blc(3, 0, 0, 0), trc(3, 0, nChannel() - 1, nRow() - 1), mat(2, nChannel(), nRow());
1328 :
1329 : // Do visCube if present
1330 0 : if (visCubeOK_p && visCube_p.nelements() > 0) {
1331 0 : blc(0) = trc(0) = 1;
1332 0 : p1.reference(visCube_p(blc, trc).reform(mat));
1333 0 : blc(0) = trc(0) = 2;
1334 0 : p2.reference(visCube_p(blc, trc).reform(mat));
1335 0 : blc(0) = trc(0) = 3;
1336 0 : p3.reference(visCube_p(blc, trc).reform(mat));
1337 0 : tmp = p3;
1338 0 : p3 = p2;
1339 0 : p2 = p1;
1340 0 : p1 = tmp;
1341 : }
1342 : // Do modelVisCube if present
1343 0 : if (modelVisCubeOK_p && modelVisCube_p.nelements() > 0) {
1344 0 : blc(0) = trc(0) = 1;
1345 0 : p1.reference(modelVisCube_p(blc, trc).reform(mat));
1346 0 : blc(0) = trc(0) = 2;
1347 0 : p2.reference(modelVisCube_p(blc, trc).reform(mat));
1348 0 : blc(0) = trc(0) = 3;
1349 0 : p3.reference(modelVisCube_p(blc, trc).reform(mat));
1350 0 : tmp = p3;
1351 0 : p3 = p2;
1352 0 : p2 = p1;
1353 0 : p1 = tmp;
1354 : }
1355 : // Do correctedVisCube if present
1356 0 : if (correctedVisCubeOK_p && correctedVisCube_p.nelements() > 0) {
1357 0 : blc(0) = trc(0) = 1;
1358 0 : p1.reference(correctedVisCube_p(blc, trc).reform(mat));
1359 0 : blc(0) = trc(0) = 2;
1360 0 : p2.reference(correctedVisCube_p(blc, trc).reform(mat));
1361 0 : blc(0) = trc(0) = 3;
1362 0 : p3.reference(correctedVisCube_p(blc, trc).reform(mat));
1363 0 : tmp = p3;
1364 0 : p3 = p2;
1365 0 : p2 = p1;
1366 0 : p1 = tmp;
1367 : }
1368 : // Do floatDataCube if present
1369 0 : if (floatDataCubeOK_p && floatDataCube_p.nelements() > 0) {
1370 0 : Matrix<Float> tmp(nChannel(), nRow());
1371 0 : Matrix<Float> p1, p2, p3;
1372 :
1373 0 : blc(0) = trc(0) = 1;
1374 0 : p1.reference(floatDataCube_p(blc, trc).reform(mat));
1375 0 : blc(0) = trc(0) = 2;
1376 0 : p2.reference(floatDataCube_p(blc, trc).reform(mat));
1377 0 : blc(0) = trc(0) = 3;
1378 0 : p3.reference(floatDataCube_p(blc, trc).reform(mat));
1379 0 : tmp = p3;
1380 0 : p3 = p2;
1381 0 : p2 = p1;
1382 0 : p1 = tmp;
1383 : }
1384 :
1385 : // Data is now back to corrType order
1386 0 : corrSorted_p = false;
1387 : }
1388 :
1389 4668 : }
1390 :
1391 13495 : Bool VisBuffer::nonCanonCorr()
1392 : {
1393 13495 : Vector<Int>& corrs(corrType());
1394 : // Only a meaningful question is all 4 corrs present
1395 13495 : if (corrs.nelements() == 4)
1396 : // (assumes corrs(0) is RR or XX)
1397 : {
1398 3828 : return (corrs(1) == Stokes::LL || corrs(1) == Stokes::YY);
1399 : } else
1400 : // Assumed OK (fewer than 4 elements, or in canonical order already)
1401 : {
1402 9667 : return false;
1403 : }
1404 : }
1405 :
1406 : // Fill weight matrix from sigma matrix
1407 4159 : void VisBuffer::resetWeightMat()
1408 : {
1409 :
1410 : // fill sigmaMat_p, size weightMat_p storage
1411 4159 : sigmaMat();
1412 4159 : IPosition ip(sigmaMat_p.shape());
1413 4159 : weightMat_p.resize(ip);
1414 :
1415 4159 : Int nPol(ip(0));
1416 4159 : Int nRow(ip(1));
1417 :
1418 : // Weight is inverse square of sigma (or zero[?])
1419 4159 : Float * w = weightMat_p.data();
1420 4159 : Float * s = sigmaMat_p.data();
1421 181874 : for (Int irow = 0; irow < nRow; ++irow)
1422 873570 : for (Int ipol = 0; ipol < nPol; ++ipol, ++w, ++s)
1423 695855 : if (*s > 0.0f) {
1424 695855 : *w = 1.0f / square(*s);
1425 : } else {
1426 0 : *w = 0.0f;
1427 : }
1428 :
1429 : // As of 2014, we define wt = 1/sigma**2 (indep of nchan)
1430 : // Scale by (unselected!) # of channels
1431 : // (to stay aligned with original nominal weights)
1432 : // Int nchan = msColumns().spectralWindow().numChan()(spectralWindow());
1433 : // weightMat_p *= Float(nchan);
1434 :
1435 : // weightMat_p now OK
1436 4159 : weightMatOK_p = true;
1437 :
1438 4159 : }
1439 :
1440 :
1441 : // Rotate visibility phase for phase center offsets
1442 2700 : void VisBuffer::phaseCenterShift(const Vector<Double>& phase)
1443 : {
1444 :
1445 2700 : AlwaysAssert(static_cast<Int>(phase.nelements()) == nRow(), AipsError);
1446 :
1447 : // phase is in metres
1448 : // phase*(-2*pi*f/c) gives phase for the channel of the given baseline in radian
1449 : // sign convention will _correct_ data
1450 :
1451 5400 : Vector<Double> freq(frequency());
1452 : Double ph, udx;
1453 2700 : Complex cph;
1454 :
1455 326700 : for (Int irow = 0; irow < nRow(); ++irow){
1456 :
1457 324000 : udx = phase(irow) * -2.0 * C::pi/C::c; // in radian/Hz
1458 :
1459 648000 : for (Int ichn = 0; ichn < nChannel(); ++ichn) {
1460 : // Calculate the Complex factor for this row and channel
1461 324000 : ph = udx * freq(ichn);
1462 :
1463 324000 : if(ph!=0.){
1464 324000 : cph = Complex(cos(ph), sin(ph));
1465 : // Shift each correlation:
1466 972000 : for (Int icor = 0; icor < nCorr(); ++icor) {
1467 648000 : if (visCubeOK_p) {
1468 648000 : visCube_p(icor, ichn, irow) *= cph;
1469 : }
1470 648000 : if (modelVisCubeOK_p) {
1471 475200 : modelVisCube_p(icor, ichn, irow) *= cph;
1472 : }
1473 648000 : if (correctedVisCubeOK_p) {
1474 518400 : correctedVisCube_p(icor, ichn, irow) *= cph;
1475 : }
1476 : // Of course floatDataCube does not have a phase to rotate.
1477 : }
1478 : }
1479 :
1480 : }
1481 : }
1482 :
1483 2700 : }
1484 :
1485 : // Rotate visibility phase for phase center offsets
1486 0 : void VisBuffer::phaseCenterShift(Double dx, Double dy)
1487 : {
1488 :
1489 : // no-op if no net shift
1490 0 : if (!(abs(dx) > 0 || abs(dy) > 0)) {
1491 0 : return;
1492 : }
1493 :
1494 : // Offsets in radians (input is arcsec)
1495 0 : dx *= (C::pi / 180.0 / 3600.0);
1496 0 : dy *= (C::pi / 180.0 / 3600.0);
1497 :
1498 : // Extra path as fraction of U and V
1499 0 : Vector<Double> udx;
1500 0 : udx = uvwMat().row(0);
1501 0 : Vector<Double> vdy;
1502 0 : vdy = uvwMat().row(1);
1503 0 : udx *= dx; // in m
1504 0 : vdy *= dy;
1505 :
1506 : // Combine axes
1507 0 : udx += vdy;
1508 :
1509 0 : phaseCenterShift(udx);
1510 :
1511 : }
1512 :
1513 :
1514 :
1515 : // Divide visCube by modelVisCube
1516 3828 : void VisBuffer::normalize(const Bool & /* phaseOnly */)
1517 : {
1518 :
1519 : // NB: phase-only now handled by SolvableVisCal
1520 : // (we will remove phaseOnly parameter later)
1521 :
1522 :
1523 : // NB: Handles pol-dep weights in chan-indep way
1524 : // TBD: optimizations?
1525 : // TBD: Handle channel-dep weights?
1526 :
1527 : // Only if all relevant columns are present
1528 3828 : if (visCubeOK_p && modelVisCubeOK_p && weightMatOK_p) {
1529 :
1530 : // cout << "Normalizing!----------------------------" << endl;
1531 :
1532 3828 : Int nCor = nCorr();
1533 :
1534 : // Amplitude data
1535 3828 : Float amp(1.0);
1536 7656 : Vector<Float> ampCorr(nCor);
1537 7656 : Vector<Int> n(nCor);
1538 :
1539 3828 : Bool * flR = flagRow().data();
1540 3828 : Bool * fl = flag().data();
1541 :
1542 176088 : for (Int irow = 0; irow < nRow(); irow++, flR++) {
1543 172260 : if (!*flR) {
1544 172260 : ampCorr = 0.0f;
1545 172260 : n = 0;
1546 1550340 : for (Int ich = 0; ich < nChannel(); ich++, fl++) {
1547 1378080 : if (!*fl) {
1548 6890400 : for (Int icorr = 0; icorr < nCor; icorr++) {
1549 5512320 : amp = abs(modelVisCube_p(icorr, ich, irow));
1550 5512320 : if (amp > 0.0f) {
1551 5512320 : visCube_p(icorr, ich, irow) = Complex( DComplex(visCube_p(icorr, ich, irow)) /
1552 5512320 : DComplex(modelVisCube_p(icorr, ich, irow)) );
1553 :
1554 2756160 : modelVisCube_p(icorr, ich, irow) = Complex(1.0);
1555 2756160 : ampCorr(icorr) += amp;
1556 2756160 : n(icorr)++;
1557 : } else
1558 : // zero data if model is zero
1559 : {
1560 2756160 : visCube_p(icorr, ich, irow) = 0.0;
1561 : }
1562 : }
1563 : }
1564 : }
1565 :
1566 861300 : for (Int icorr = 0; icorr < nCor; icorr++) {
1567 689040 : if (n(icorr) > 0) {
1568 344520 : weightMat_p(icorr, irow) *= square(ampCorr(icorr) / Float(n(icorr)));
1569 : } else {
1570 344520 : weightMat_p(icorr, irow) = 0.0f;
1571 : }
1572 : }
1573 :
1574 : } else {
1575 : // Zero weight on this flagged row
1576 0 : weightMat_p.column(irow) = 0.0f;
1577 :
1578 : // Advance fl over this row
1579 0 : fl += nChannel();
1580 : }
1581 3828 : }
1582 : } else {
1583 0 : throw(AipsError("Failed to normalize data by model!"));
1584 : }
1585 3828 : }
1586 :
1587 0 : Vector<Int> VisBuffer::vecIntRange(const MSCalEnums::colDef & calEnum) const
1588 : {
1589 : // Return a column range for a generic integer column as
1590 : // identified by the enum specification in class MSCalEnums
1591 :
1592 : // Prepare the flag column masking
1593 0 : LogicalArray mask(!flagRow());
1594 : MaskedArray<Int>* maskArray;
1595 :
1596 : // A dummy vector for columns not yet supported (returns a value of [-1]);
1597 0 : Vector<Int> nullIndex(antenna1().shape(), -1);
1598 :
1599 0 : switch (calEnum) {
1600 : // ANTENNA1
1601 0 : case MSC::ANTENNA1: {
1602 0 : maskArray = new MaskedArray<Int>(antenna1(), mask);
1603 0 : break;
1604 : };
1605 : // ANTENNA2
1606 0 : case MSC::ANTENNA2: {
1607 0 : maskArray = new MaskedArray<Int>(antenna2(), mask);
1608 0 : break;
1609 : };
1610 : // FEED1
1611 0 : case MSC::FEED1: {
1612 0 : maskArray = new MaskedArray<Int>(feed1(), mask);
1613 0 : break;
1614 : };
1615 : // FIELD_ID
1616 0 : case MSC::FIELD_ID: {
1617 0 : Vector<Int> fieldIdVec(antenna1().shape(), fieldId());
1618 0 : maskArray = new MaskedArray<Int>(fieldIdVec, mask);
1619 0 : break;
1620 : };
1621 : // ARRAY_ID
1622 0 : case MSC::ARRAY_ID: {
1623 0 : maskArray = new MaskedArray<Int>(nullIndex, mask);
1624 0 : break;
1625 : };
1626 : // OBSERVATION_ID
1627 0 : case MSC::OBSERVATION_ID: {
1628 0 : maskArray = new MaskedArray<Int>(nullIndex, mask);
1629 0 : break;
1630 : };
1631 : // SCAN_NUMBER
1632 0 : case MSC::SCAN_NUMBER: {
1633 0 : maskArray = new MaskedArray<Int>(nullIndex, mask);
1634 0 : break;
1635 : };
1636 : // PROCESSOR_ID
1637 0 : case MSC::PROCESSOR_ID: {
1638 0 : maskArray = new MaskedArray<Int>(nullIndex, mask);
1639 0 : break;
1640 : };
1641 : // PHASE_ID
1642 0 : case MSC::PHASE_ID: {
1643 0 : maskArray = new MaskedArray<Int>(nullIndex, mask);
1644 0 : break;
1645 : };
1646 : // STATE_ID
1647 0 : case MSC::STATE_ID: {
1648 0 : maskArray = new MaskedArray<Int>(nullIndex, mask);
1649 0 : break;
1650 : };
1651 : // PULSAR_BIN
1652 0 : case MSC::PULSAR_BIN: {
1653 0 : maskArray = new MaskedArray<Int>(nullIndex, mask);
1654 0 : break;
1655 : };
1656 : // PULSAR_GATE_ID
1657 0 : case MSC::PULSAR_GATE_ID: {
1658 0 : maskArray = new MaskedArray<Int>(nullIndex, mask);
1659 0 : break;
1660 : };
1661 : // FREQ_GROUP
1662 0 : case MSC::FREQ_GROUP: {
1663 0 : maskArray = new MaskedArray<Int>(nullIndex, mask);
1664 0 : break;
1665 : };
1666 : // CALIBRATION_GROUP
1667 0 : case MSC::CALIBRATION_GROUP: {
1668 0 : maskArray = new MaskedArray<Int>(nullIndex, mask);
1669 0 : break;
1670 : };
1671 :
1672 0 : default: {
1673 0 : throw(AipsError("Request for non-existent uv-data column"));
1674 : };
1675 : };
1676 :
1677 : // Return only unique indices
1678 0 : Vector<Int> retval = unique(maskArray->getCompressedArray());
1679 0 : if (maskArray) {
1680 0 : delete(maskArray);
1681 : }
1682 :
1683 0 : return retval;
1684 : };
1685 :
1686 0 : Vector<Int> VisBuffer::antIdRange() const
1687 : {
1688 : // Return a column range for ANTENNA_ID, including the
1689 : // union of the ANTENNA1 and ANTENNA2 columns indices
1690 :
1691 0 : Vector<Int> ant1 = vecIntRange(MSC::ANTENNA1);
1692 0 : Vector<Int> ant2 = vecIntRange(MSC::ANTENNA2);
1693 0 : Vector<Int> ant12 = concatenateArray(ant1, ant2);
1694 :
1695 : // Return only unique index values
1696 0 : return unique(ant12);
1697 : };
1698 :
1699 0 : Bool VisBuffer::timeRange(MEpoch & rTime, MVEpoch & rTimeEP,
1700 : MVEpoch & rInterval) const
1701 : {
1702 : // Return the time range of data in the vis buffer
1703 : // (makes simplistic assumptions in the absence of
1704 : // interval information for now)
1705 :
1706 : // Initialization
1707 0 : Bool retval = false;
1708 :
1709 0 : if (nRow() > 0) {
1710 0 : retval = true;
1711 0 : LogicalArray mask(!flagRow());
1712 0 : MaskedArray<Double> maskTime(time(), mask);
1713 0 : Double minTime = min(maskTime);
1714 0 : Double maxTime = max(maskTime);
1715 : // Mean time
1716 0 : rTime = MEpoch(Quantity((minTime + maxTime) / 2, "s"));
1717 : // Extra precision time is always null for now
1718 0 : rTimeEP = MVEpoch(Quantity(0, "s"));
1719 : // Interval
1720 0 : rInterval = MVEpoch(Quantity(maxTime - minTime, "s"));
1721 : };
1722 0 : return retval;
1723 : };
1724 :
1725 :
1726 2704 : Vector<casacore::rownr_t>& VisBuffer::rowIds()
1727 : {
1728 2704 : if (!rowIdsOK_p) {
1729 :
1730 2704 : rowIdsOK_p = true;
1731 2704 : visIter_p->rowIds(rowIds_p);
1732 : }
1733 2704 : return rowIds_p;
1734 : }
1735 :
1736 :
1737 : void
1738 14 : VisBuffer::updateCoordInfo(const VisBuffer * vb, const Bool dirDependent )
1739 : {
1740 14 : updateCoord (vb, vb->antenna1OK (), & VisBuffer::antenna1, antenna1_p, antenna1OK_p);
1741 14 : updateCoord (vb, vb->antenna2OK (), & VisBuffer::antenna2, antenna2_p, antenna2OK_p);
1742 14 : updateCoordS (vb, vb->arrayIdOK (), & VisBuffer::arrayId, arrayId_p, arrayIdOK_p);
1743 14 : updateCoordS (vb, vb->dataDescriptionIdOK(), & VisBuffer::dataDescriptionId, dataDescriptionId_p, dataDescriptionIdOK_p);
1744 14 : updateCoordS (vb, vb->fieldIdOK (), & VisBuffer::fieldId, fieldId_p, fieldIdOK_p);
1745 14 : updateCoordS (vb, vb->spectralWindowOK (), & VisBuffer::spectralWindow, spectralWindow_p, spectralWindowOK_p);
1746 14 : updateCoord (vb, vb->timeOK (), & VisBuffer::time, time_p, timeOK_p);
1747 14 : updateCoord (vb, vb->frequencyOK (), & VisBuffer::frequency, frequency_p, frequencyOK_p);
1748 14 : updateCoordS (vb, vb->nRowOK (), & VisBuffer::nRow, nRow_p, nRowOK_p);
1749 :
1750 14 : vb->copyMsInfo (oldMSId_p, msOK_p, newMS_p);
1751 :
1752 14 : updateCoord (vb, vb->feed1OK (), & VisBuffer::feed1, feed1_p, feed1OK_p);
1753 14 : updateCoord (vb, vb->feed2OK (), & VisBuffer::feed2, feed2_p, feed2OK_p);
1754 :
1755 14 : if(dirDependent){
1756 0 : updateCoord (vb, vb->feed1_paOK (), & VisBuffer::feed1_pa, feed1_pa_p, feed1_paOK_p);
1757 0 : updateCoord (vb, vb->feed2_paOK (), & VisBuffer::feed2_pa, feed2_pa_p, feed2_paOK_p);
1758 0 : updateCoord (vb, vb->direction1OK (), & VisBuffer::direction1, direction1_p, direction1OK_p);
1759 0 : updateCoord (vb, vb->direction2OK (), & VisBuffer::direction2, direction2_p, direction2OK_p);
1760 : }
1761 :
1762 14 : }
1763 :
1764 4532 : void VisBuffer::setVisCube(Complex c)
1765 : {
1766 4532 : visCube_p.resize(visIter_p->visibilityShape());
1767 4532 : visCube_p.set(c);
1768 4532 : visCubeOK_p = true;
1769 4532 : }
1770 12800 : void VisBuffer::setModelVisCube(Complex c)
1771 : {
1772 12800 : modelVisCube_p.resize(visIter_p->visibilityShape());
1773 12800 : modelVisCube_p.set(c);
1774 12800 : modelVisCubeOK_p = true;
1775 12800 : }
1776 0 : void VisBuffer::setCorrectedVisCube(Complex c)
1777 : {
1778 0 : correctedVisCube_p.resize(visIter_p->visibilityShape());
1779 0 : correctedVisCube_p.set(c);
1780 0 : correctedVisCubeOK_p = true;
1781 0 : }
1782 0 : void VisBuffer::setVisCube(const Cube<Complex>& vis)
1783 : {
1784 0 : visCube_p.resize(vis.shape());
1785 0 : visCube_p = vis;
1786 0 : visCubeOK_p = true;
1787 0 : }
1788 5598 : void VisBuffer::setModelVisCube(const Cube<Complex>& vis)
1789 : {
1790 5598 : modelVisCube_p.resize(vis.shape());
1791 5598 : modelVisCube_p = vis;
1792 5598 : modelVisCubeOK_p = true;
1793 5598 : }
1794 :
1795 3828 : void VisBuffer::setModelVisCube(const Vector<Float>& stokes)
1796 : {
1797 :
1798 : /*
1799 : cout << "Specified Stokes Parameters: " << stokes << endl;
1800 :
1801 : cout << "polFrame() = " << polFrame()
1802 : << " " << MSIter::Circular
1803 : << " " << MSIter::Linear
1804 : << endl;
1805 : */
1806 :
1807 : // Stokes parameters, nominally unpolarized, unit I
1808 3828 : Float I(1.0), Q(0.0), U(0.0), V(0.0);
1809 :
1810 : // Only fill as many as are specified, up to 4 (unspecified will be assumed zero)
1811 19140 : for (uInt i = 0; i < stokes.nelements(); ++i)
1812 15312 : switch (i) {
1813 3828 : case 0: {
1814 3828 : I = stokes(i);
1815 3828 : break;
1816 : }
1817 3828 : case 1: {
1818 3828 : Q = stokes(i);
1819 3828 : break;
1820 : }
1821 3828 : case 2: {
1822 3828 : U = stokes(i);
1823 3828 : break;
1824 : }
1825 3828 : case 3: {
1826 3828 : V = stokes(i);
1827 3828 : break;
1828 : }
1829 0 : default: {
1830 0 : break;
1831 : }
1832 : }
1833 :
1834 : // Convert to correlations, according to basis
1835 7656 : Vector<Complex> stkvis(4, Complex(0.0)); // initially all zero
1836 3828 : switch (polFrame()) {
1837 0 : case MSIter::Circular: {
1838 0 : stkvis(0) = Complex(I + V);
1839 0 : stkvis(1) = Complex(Q, U);
1840 0 : stkvis(2) = Complex(Q, -U);
1841 0 : stkvis(3) = Complex(I - V);
1842 0 : break;
1843 : }
1844 3828 : case MSIter::Linear: {
1845 3828 : stkvis(0) = Complex(I + Q);
1846 3828 : stkvis(1) = Complex(U, V);
1847 3828 : stkvis(2) = Complex(U, -V);
1848 3828 : stkvis(3) = Complex(I - Q);
1849 3828 : break;
1850 : }
1851 0 : default:
1852 0 : throw(AipsError("Model-setting only works for CIRCULAR and LINEAR bases, for now."));
1853 : break;
1854 : }
1855 :
1856 : // A map onto the actual correlations in the VisBuffer
1857 3828 : Vector<Int> corrmap;
1858 3828 : corrmap = corrType();
1859 3828 : corrmap -= corrmap(0);
1860 : // This MUST yield indices < 4
1861 3828 : if (max(corrmap) > 3) {
1862 0 : throw(AipsError("HELP! The correlations in the data are not normal!"));
1863 : }
1864 :
1865 :
1866 : // Set the modelVisCube accordingly
1867 3828 : modelVisCube_p.resize(visIter_p->visibilityShape());
1868 3828 : modelVisCube_p.set(0.0);
1869 19140 : for (Int icorr = 0; icorr < nCorr(); ++icorr)
1870 15312 : if (abs(stkvis(corrmap(icorr))) > 0.0) {
1871 7656 : modelVisCube_p(Slice(icorr, 1, 1), Slice(), Slice()).set(stkvis(corrmap(icorr)));
1872 : }
1873 3828 : modelVisCubeOK_p = true;
1874 :
1875 : // Lookup flux density calibrator scaling, and apply it per channel...
1876 : // TBD
1877 :
1878 :
1879 3828 : }
1880 :
1881 : Int
1882 0 : VisBuffer::nRowChunk() const
1883 : {
1884 0 : CheckVisIter ();
1885 0 : return visIter_p->nRowChunk ();
1886 : }
1887 :
1888 :
1889 14 : Int VisBuffer::numberCoh () const
1890 : {
1891 14 : CheckVisIter ();
1892 :
1893 14 : return visIter_p -> numberCoh ();
1894 : }
1895 :
1896 :
1897 : void
1898 344524 : VisBuffer::checkVisIter (const char * func, const char * file, int line, const char * extra) const
1899 : {
1900 344524 : checkVisIterBase (func, file, line, extra);
1901 344524 : }
1902 :
1903 : void
1904 345254 : VisBuffer::checkVisIterBase (const char * func, const char * file, int line, const char * extra) const
1905 : {
1906 345254 : if (visIter_p == NULL) {
1907 0 : throw AipsError (String ("No VisibilityIterator is available to fill this field in (") +
1908 0 : func + extra + ")", file, line);
1909 : }
1910 345254 : }
1911 :
1912 :
1913 0 : void VisBuffer::setCorrectedVisCube(const Cube<Complex>& vis)
1914 : {
1915 0 : correctedVisCube_p.resize(vis.shape());
1916 0 : correctedVisCube_p = vis;
1917 0 : correctedVisCubeOK_p = true;
1918 0 : }
1919 :
1920 0 : void VisBuffer::setFloatDataCube(const Cube<Float>& fcube)
1921 : {
1922 0 : floatDataCube_p.resize(fcube.shape());
1923 0 : floatDataCube_p = fcube;
1924 0 : floatDataCubeOK_p = true;
1925 0 : }
1926 :
1927 0 : void VisBuffer::refModelVis(const Matrix<CStokesVector>& mvis)
1928 : {
1929 0 : modelVisibility_p.resize();
1930 0 : modelVisibility_p.reference(mvis);
1931 0 : modelVisibilityOK_p = true;
1932 0 : }
1933 :
1934 0 : void VisBuffer::removeScratchCols()
1935 : {
1936 : // removes scratch data from the vb
1937 0 : modelVisibility_p.resize();
1938 0 : modelVisibilityOK_p = false;
1939 0 : correctedVisibility_p.resize();
1940 0 : correctedVisibilityOK_p = false;
1941 0 : }
1942 :
1943 11437 : Int& VisBuffer::fillnCorr()
1944 : {
1945 11437 : CheckVisIter ();
1946 11437 : nCorrOK_p = true;
1947 11437 : nCorr_p = corrType().nelements();
1948 11437 : return nCorr_p;
1949 : }
1950 :
1951 5121 : Vector<Int>& VisBuffer::fillObservationId()
1952 : {
1953 5121 : CheckVisIter();
1954 5121 : observationIdOK_p = true;
1955 5121 : return visIter_p->observationId(observationId_p);
1956 : }
1957 :
1958 0 : Vector<Int>& VisBuffer::fillProcessorId()
1959 : {
1960 0 : CheckVisIter();
1961 0 : processorIdOK_p = true;
1962 0 : return visIter_p->processorId(processorId_p);
1963 : }
1964 :
1965 0 : Vector<Int>& VisBuffer::fillStateId()
1966 : {
1967 0 : CheckVisIter();
1968 0 : stateIdOK_p = true;
1969 0 : return visIter_p->stateId(stateId_p);
1970 : }
1971 :
1972 0 : Array<Bool>& VisBuffer::fillFlagCategory()
1973 : {
1974 0 : CheckVisIter();
1975 0 : flagCategoryOK_p = true;
1976 0 : return visIter_p->flagCategory(flagCategory_p);
1977 : }
1978 :
1979 18933 : Int& VisBuffer::fillnChannel()
1980 : {
1981 18933 : CheckVisIter ();
1982 18933 : nChannelOK_p = true;
1983 : // nChannel_p=visIter_p->channelGroupSize();
1984 18933 : nChannel_p = channel().nelements();
1985 18933 : return nChannel_p;
1986 : }
1987 :
1988 18933 : Vector<Int>& VisBuffer::fillChannel()
1989 : {
1990 18933 : CheckVisIter ();
1991 18933 : channelOK_p = true;
1992 18933 : return visIter_p->channel(channel_p);
1993 : }
1994 :
1995 24671 : Int& VisBuffer::fillnRow()
1996 : {
1997 24671 : CheckVisIter ();
1998 24671 : nRowOK_p = true;
1999 24671 : nRow_p = visIter_p->nRow();
2000 24671 : return nRow_p;
2001 : }
2002 :
2003 16843 : Vector<Int>& VisBuffer::fillAnt1()
2004 : {
2005 16843 : CheckVisIter ();
2006 16843 : antenna1OK_p = true;
2007 16843 : return visIter_p->antenna1(antenna1_p);
2008 : }
2009 :
2010 14340 : Vector<Int>& VisBuffer::fillAnt2()
2011 : {
2012 14340 : CheckVisIter ();
2013 14340 : antenna2OK_p = true;
2014 14340 : return visIter_p->antenna2(antenna2_p);
2015 : }
2016 :
2017 379 : Vector<Int>& VisBuffer::fillFeed1()
2018 : {
2019 379 : CheckVisIter ();
2020 379 : feed1OK_p = true;
2021 379 : return visIter_p->feed1(feed1_p);
2022 : }
2023 :
2024 379 : Vector<Int>& VisBuffer::fillFeed2()
2025 : {
2026 379 : CheckVisIter ();
2027 379 : feed2OK_p = true;
2028 379 : return visIter_p->feed2(feed2_p);
2029 : }
2030 :
2031 0 : Vector<SquareMatrix<Complex, 2> >& VisBuffer::fillCjones()
2032 : {
2033 0 : CheckVisIter ();
2034 0 : cjonesOK_p = true;
2035 0 : return visIter_p->CJones(cjones_p);
2036 : }
2037 :
2038 18752 : Vector<Int>& VisBuffer::fillCorrType()
2039 : {
2040 18752 : CheckVisIter ();
2041 18752 : corrTypeOK_p = true;
2042 18752 : return visIter_p->corrType(corrType_p);
2043 : }
2044 :
2045 : // calling fillFeed1_pa or fillFeed2_pa will fill antenna, feed
2046 : // and time caches automatically
2047 365 : Vector<Float>& VisBuffer::fillFeed1_pa()
2048 : {
2049 365 : CheckVisIterBase ();
2050 :
2051 : // fill feed, antenna and time caches, if not filled before
2052 365 : feed1();
2053 365 : antenna1();
2054 365 : time();
2055 365 : feed1_paOK_p = true;
2056 365 : feed1_pa_p.resize(antenna1_p.nelements()); // could also use nRow()
2057 :
2058 : // now actual calculations
2059 292704 : for (uInt row = 0; row < feed1_pa_p.nelements(); ++row) {
2060 584678 : const Vector<Float>& ant_pa = feed_pa(time_p(row)); // caching inside
2061 : // ROVisibilityIterator, if the time doesn't change. Otherwise
2062 : // we should probably fill both buffers for feed1 and feed2
2063 : // simultaneously to speed up things.
2064 :
2065 292339 : DebugAssert((uInt(antenna1_p(row)) < ant_pa.nelements()), AipsError);
2066 292339 : DebugAssert(antenna1_p(row) >= 0, AipsError);
2067 292339 : feed1_pa_p(row) = ant_pa(antenna1_p(row));
2068 : // currently feed_pa returns only the first feed position angle
2069 : // we need to add an offset if this row corresponds to a
2070 : // different feed
2071 292339 : if (feed1_p(row)) // an if-statement to avoid unnecessary operations
2072 : // in the single feed case, everything would
2073 : // work without it.
2074 0 : feed1_pa_p(row) += visIter_p->receptorAngles()(0,
2075 0 : antenna1_p(row), feed1_p(row)) -
2076 0 : visIter_p->receptorAngles()(0, antenna1_p(row), 0);
2077 : }
2078 365 : return feed1_pa_p;
2079 : }
2080 :
2081 365 : Vector<Float>& VisBuffer::fillFeed2_pa()
2082 : {
2083 365 : CheckVisIterBase ();
2084 :
2085 : // fill feed, antenna and time caches, if not filled before
2086 365 : feed2();
2087 365 : antenna2();
2088 365 : time();
2089 365 : feed2_paOK_p = true;
2090 365 : feed2_pa_p.resize(antenna2_p.nelements()); // could also use nRow()
2091 :
2092 : // now actual calculations
2093 292704 : for (uInt row = 0; row < feed2_pa_p.nelements(); ++row) {
2094 584678 : const Vector<Float>& ant_pa = feed_pa(time_p(row)); // caching inside
2095 : // ROVisibilityIterator, if the time doesn't change. Otherwise
2096 : // we should probably fill both buffers for feed1 and feed2
2097 : // simultaneously to speed up things.
2098 292339 : DebugAssert((uInt(antenna2_p(row)) < ant_pa.nelements()), AipsError);
2099 292339 : DebugAssert(antenna2_p(row) >= 0, AipsError);
2100 292339 : feed2_pa_p(row) = ant_pa(antenna2_p(row));
2101 : // currently feed_pa returns only the first feed position angle
2102 : // we need to add an offset if this row correspods to a
2103 : // different feed
2104 292339 : if (feed2_p(row)) // an if-statement to avoid unnecessary operations
2105 : // in the single feed case, everything would
2106 : // work without it.
2107 0 : feed2_pa_p(row) += visIter_p->receptorAngles()(0,
2108 0 : antenna2_p(row), feed2_p(row)) -
2109 0 : visIter_p->receptorAngles()(0, antenna2_p(row), 0);
2110 : }
2111 365 : return feed2_pa_p;
2112 : }
2113 :
2114 0 : Vector<MDirection>& VisBuffer::fillDirection1()
2115 : {
2116 : //Timer tim;
2117 : //tim.mark();
2118 0 : CheckVisIterBase ();
2119 : // fill feed1_pa cache, antenna, feed and time will be filled automatically
2120 : //feed1_pa();
2121 : // commented above as it calculates the Par-angle for all antennas when
2122 : //it may not be used at all...a 300% speedup on some large pointing table
2123 0 : antenna1();
2124 0 : feed1();
2125 0 : time();
2126 0 : direction1OK_p = true;
2127 0 : firstDirection1OK_p=true;
2128 0 : direction1_p.resize(antenna1_p.nelements()); // could also use nRow()
2129 0 : const MSPointingColumns & mspc = msColumns().pointing();
2130 0 : lastPointTableRow_p = mspc.pointingIndex(antenna1()(0),
2131 0 : time()(0), lastPointTableRow_p);
2132 0 : if (visIter_p->allBeamOffsetsZero() && lastPointTableRow_p < 0) {
2133 : // if no true pointing information is found
2134 : // just return the phase center from the field table
2135 0 : direction1_p.set(phaseCenter());
2136 0 : lastPointTableRow_p = 0;
2137 0 : return direction1_p;
2138 : }
2139 0 : for (uInt row = 0; row < antenna1_p.nelements(); ++row) {
2140 0 : DebugAssert(antenna1_p(row) >= 0, AipsError);
2141 0 : DebugAssert(feed1_p(row) >= 0, AipsError);
2142 0 : Int pointIndex1 = mspc.pointingIndex(antenna1()(row), time()(row), lastPointTableRow_p);
2143 :
2144 : //cout << "POINTINDEX " << pointIndex1 << endl;
2145 : // if no true pointing information is found
2146 : // use the phase center from the field table
2147 0 : if (pointIndex1 >= 0) {
2148 0 : lastPointTableRow_p = pointIndex1;
2149 0 : direction1_p(row) = mspc.directionMeas(pointIndex1, timeInterval()(row));
2150 : } else {
2151 0 : direction1_p(row) = phaseCenter();
2152 : }
2153 0 : if (!visIter_p->allBeamOffsetsZero()) {
2154 : //Now we can calculate the Par-angles
2155 0 : feed1_pa();
2156 : RigidVector<Double, 2> beamOffset =
2157 0 : visIter_p->getBeamOffsets()(0, antenna1_p(row), feed1_p(row));
2158 0 : if (visIter_p->antennaMounts()(antenna1_p(row)) == "ALT-AZ" ||
2159 0 : visIter_p->antennaMounts()(antenna1_p(row)) == "alt-az") {
2160 0 : SquareMatrix<Double, 2> xform(SquareMatrix<Double, 2>::General);
2161 : // SquareMatrix' default constructor is a bit strange.
2162 : // We will probably need to change it in the future
2163 0 : Double cpa = cos(feed1_pa_p(row));
2164 0 : Double spa = sin(feed1_pa_p(row));
2165 0 : xform(0, 0) = cpa;
2166 0 : xform(1, 1) = cpa;
2167 0 : xform(0, 1) = -spa;
2168 0 : xform(1, 0) = spa;
2169 0 : beamOffset *= xform; // parallactic angle rotation
2170 : }
2171 : // x direction is flipped to convert az-el type frame to ra-dec
2172 0 : direction1_p(row).shift(-beamOffset(0), beamOffset(1), true);
2173 : }
2174 : }
2175 : //tim.show("fill dir1");
2176 : //cerr << "allbeamOff " << visIter_p->allBeamOffsetsZero() << endl;
2177 0 : firstDirection1_p=direction1_p[0];
2178 0 : return direction1_p;
2179 : }
2180 :
2181 0 : MDirection& VisBuffer::fillFirstDirection1()
2182 : {
2183 : //Timer tim;
2184 : //tim.mark();
2185 0 : CheckVisIterBase ();
2186 : // fill feed1_pa cache, antenna, feed and time will be filled automatically
2187 0 : feed1();
2188 0 : time();
2189 0 : antenna1();
2190 :
2191 : //feed1_pa();
2192 0 : firstDirection1OK_p=true;
2193 0 : const MSPointingColumns & mspc = msColumns().pointing();
2194 0 : lastPointTableRow_p = mspc.pointingIndex(antenna1()(0),
2195 0 : time()(0), lastPointTableRow_p);
2196 0 : if (visIter_p->allBeamOffsetsZero() && lastPointTableRow_p < 0) {
2197 : // if no true pointing information is found
2198 : // just return the phase center from the field table
2199 0 : firstDirection1_p=phaseCenter();
2200 0 : lastPointTableRow_p = 0;
2201 0 : return firstDirection1_p;
2202 : }
2203 :
2204 0 : Int pointIndex1=lastPointTableRow_p;
2205 :
2206 : //cout << "POINTINDEX " << pointIndex1 << endl;
2207 : // if no true pointing information is found
2208 : // use the phase center from the field table
2209 0 : if (pointIndex1 >= 0) {
2210 :
2211 0 : firstDirection1_p = mspc.directionMeas(pointIndex1, timeInterval()(0));
2212 : } else {
2213 0 : firstDirection1_p = phaseCenter();
2214 : }
2215 0 : if (!visIter_p->allBeamOffsetsZero()) {
2216 0 : feed1_pa();
2217 : RigidVector<Double, 2> beamOffset =
2218 0 : visIter_p->getBeamOffsets()(0, antenna1_p(0), feed1_p(0));
2219 0 : if (visIter_p->antennaMounts()(antenna1_p(0)) == "ALT-AZ" ||
2220 0 : visIter_p->antennaMounts()(antenna1_p(0)) == "alt-az") {
2221 0 : SquareMatrix<Double, 2> xform(SquareMatrix<Double, 2>::General);
2222 : // SquareMatrix' default constructor is a bit strange.
2223 : // We will probably need to change it in the future
2224 0 : Double cpa = cos(feed1_pa_p(0));
2225 0 : Double spa = sin(feed1_pa_p(0));
2226 0 : xform(0, 0) = cpa;
2227 0 : xform(1, 1) = cpa;
2228 0 : xform(0, 1) = -spa;
2229 0 : xform(1, 0) = spa;
2230 0 : beamOffset *= xform; // parallactic angle rotation
2231 : }
2232 : // x direction is flipped to convert az-el type frame to ra-dec
2233 0 : firstDirection1_p.shift(-beamOffset(0), beamOffset(1), true);
2234 : }
2235 :
2236 0 : return firstDirection1_p;
2237 : }
2238 0 : Vector<MDirection>& VisBuffer::fillDirection2()
2239 : {
2240 0 : CheckVisIterBase ();
2241 : // fill feed2_pa cache, antenna, feed and time will be filled automatically
2242 0 : feed2_pa();
2243 0 : direction2OK_p = true;
2244 0 : direction2_p.resize(antenna2_p.nelements()); // could also use nRow()
2245 0 : const MSPointingColumns & mspc = msColumns().pointing();
2246 0 : lastPointTableRow_p = mspc.pointingIndex(antenna2()(0), time()(0), lastPointTableRow_p);
2247 0 : if (visIter_p->allBeamOffsetsZero() && lastPointTableRow_p < 0) {
2248 : // if no true pointing information is found
2249 : // just return the phase center from the field table
2250 0 : direction2_p.set(phaseCenter());
2251 0 : lastPointTableRow_p = 0;
2252 0 : return direction2_p;
2253 : }
2254 0 : for (uInt row = 0; row < antenna2_p.nelements(); ++row) {
2255 0 : DebugAssert(antenna2_p(row) >= 0, AipsError);
2256 0 : DebugAssert(feed2_p(row) >= 0, AipsError);
2257 0 : Int pointIndex2 = mspc.pointingIndex(antenna2()(row), time()(row), lastPointTableRow_p);
2258 : // if no true pointing information is found
2259 : // use the phase center from the field table
2260 0 : if (pointIndex2 >= 0) {
2261 0 : lastPointTableRow_p = pointIndex2;
2262 0 : direction2_p(row) = mspc.directionMeas(pointIndex2, timeInterval()(row));
2263 : } else {
2264 0 : direction2_p(row) = phaseCenter();
2265 : }
2266 0 : if (!visIter_p->allBeamOffsetsZero()) {
2267 : RigidVector<Double, 2> beamOffset =
2268 0 : visIter_p->getBeamOffsets()(0, antenna2_p(row), feed2_p(row));
2269 0 : if (visIter_p->antennaMounts()(antenna2_p(row)) == "ALT-AZ" ||
2270 0 : visIter_p->antennaMounts()(antenna2_p(row)) == "alt-az") {
2271 0 : SquareMatrix<Double, 2> xform(SquareMatrix<Double, 2>::General);
2272 : // SquareMatrix' default constructor is a bit strange.
2273 : // We will probably need to change it in the future
2274 0 : Double cpa = cos(feed2_pa_p(row));
2275 0 : Double spa = sin(feed2_pa_p(row));
2276 0 : xform(0, 0) = cpa;
2277 0 : xform(1, 1) = cpa;
2278 0 : xform(0, 1) = -spa;
2279 0 : xform(1, 0) = spa;
2280 0 : beamOffset *= xform; // parallactic angle rotation
2281 : }
2282 : // x direction is flipped to convert az-el type frame to ra-dec
2283 0 : direction2_p(row).shift(-beamOffset(0), beamOffset(1), true);
2284 : }
2285 : }
2286 0 : return direction2_p;
2287 : }
2288 :
2289 19340 : Int& VisBuffer::fillFieldId()
2290 : {
2291 19340 : CheckVisIter ();
2292 19340 : fieldIdOK_p = true;
2293 19340 : fieldId_p = visIter_p->fieldId();
2294 19340 : return fieldId_p;
2295 : }
2296 :
2297 629 : Int& VisBuffer::fillArrayId()
2298 : {
2299 629 : CheckVisIter ();
2300 629 : arrayIdOK_p = true;
2301 629 : arrayId_p = visIter_p->arrayId();
2302 629 : return arrayId_p;
2303 : }
2304 :
2305 264 : Int& VisBuffer::fillDataDescriptionId ()
2306 : {
2307 264 : CheckVisIter ();
2308 264 : dataDescriptionIdOK_p = true;
2309 264 : dataDescriptionId_p = visIter_p->dataDescriptionId();
2310 264 : return dataDescriptionId_p;
2311 : }
2312 :
2313 10959 : Matrix<Bool>& VisBuffer::fillFlag()
2314 : {
2315 10959 : CheckVisIter ();
2316 10959 : flagOK_p = true;
2317 10959 : return visIter_p->flag(flag_p);
2318 : }
2319 :
2320 5393 : Cube<Bool>& VisBuffer::fillFlagCube()
2321 : {
2322 5393 : CheckVisIter ();
2323 5393 : flagCubeOK_p = true;
2324 5393 : return visIter_p->flag(flagCube_p);
2325 : }
2326 :
2327 15969 : Vector<Bool>& VisBuffer::fillFlagRow()
2328 : {
2329 15969 : CheckVisIter ();
2330 15969 : flagRowOK_p = true;
2331 15969 : return visIter_p->flagRow(flagRow_p);
2332 : }
2333 :
2334 4985 : Vector<Int>& VisBuffer::fillScan()
2335 : {
2336 4985 : CheckVisIter ();
2337 4985 : scanOK_p = true;
2338 4985 : return visIter_p->scan(scan_p);
2339 : }
2340 :
2341 24188 : Vector<Double>& VisBuffer::fillFreq()
2342 : {
2343 24188 : CheckVisIter ();
2344 24188 : frequencyOK_p = true;
2345 24188 : return visIter_p->frequency(frequency_p);
2346 : }
2347 :
2348 : //Vector<Double>& VisBuffer::fillLSRFreq()
2349 : //{
2350 : // CheckVisIter ();
2351 : // lsrFrequencyOK_p = true;
2352 : // return visIter_p->lsrFrequency(lsrFrequency_p);
2353 : //}
2354 :
2355 338 : MDirection& VisBuffer::fillPhaseCenter()
2356 : {
2357 338 : CheckVisIter ();
2358 338 : phaseCenterOK_p = true;
2359 338 : return phaseCenter_p = visIter_p->phaseCenter();
2360 : }
2361 5821 : const MDirection VisBuffer::phaseCenter(const Double time) const
2362 : {
2363 5821 : CheckVisIter ();
2364 :
2365 5821 : return This->phaseCenter(This->fieldId(), time);
2366 : }
2367 5821 : const MDirection VisBuffer::phaseCenter(const Int field, const Double time) const
2368 : {
2369 5821 : CheckVisIter ();
2370 :
2371 5821 : return visIter_p->phaseCenter(field, time);
2372 : }
2373 :
2374 12245 : Int& VisBuffer::fillPolFrame()
2375 : {
2376 12245 : CheckVisIter ();
2377 12245 : polFrameOK_p = true;
2378 12245 : polFrame_p = visIter_p->polFrame();
2379 12245 : return polFrame_p;
2380 : }
2381 :
2382 0 : Vector<Float>& VisBuffer::fillSigma()
2383 : {
2384 0 : CheckVisIter ();
2385 0 : sigmaOK_p = true;
2386 0 : return visIter_p->sigma(sigma_p);
2387 : }
2388 :
2389 4171 : Matrix<Float>& VisBuffer::fillSigmaMat()
2390 : {
2391 4171 : CheckVisIter ();
2392 4171 : sigmaMatOK_p = true;
2393 4171 : return visIter_p->sigmaMat(sigmaMat_p);
2394 : }
2395 :
2396 15640 : Int& VisBuffer::fillSpW()
2397 : {
2398 15640 : CheckVisIter ();
2399 15640 : spectralWindowOK_p = true;
2400 15640 : spectralWindow_p = visIter_p->spectralWindow();
2401 15640 : return spectralWindow_p;
2402 : }
2403 :
2404 25514 : Vector<Double>& VisBuffer::fillTime()
2405 : {
2406 25514 : CheckVisIter ();
2407 25514 : timeOK_p = true;
2408 25514 : return visIter_p->time(time_p);
2409 : }
2410 :
2411 0 : Vector<Double>& VisBuffer::fillTimeCentroid()
2412 : {
2413 0 : CheckVisIter ();
2414 0 : timeCentroidOK_p = true;
2415 0 : return visIter_p->timeCentroid(timeCentroid_p);
2416 : }
2417 :
2418 2499 : Vector<Double>& VisBuffer::fillTimeInterval()
2419 : {
2420 2499 : CheckVisIter ();
2421 2499 : timeIntervalOK_p = true;
2422 2499 : return visIter_p->timeInterval(timeInterval_p);
2423 : }
2424 :
2425 0 : Vector<Double>& VisBuffer::fillExposure()
2426 : {
2427 0 : CheckVisIter ();
2428 0 : exposureOK_p = true;
2429 0 : return visIter_p->exposure(exposure_p);
2430 : }
2431 :
2432 :
2433 14923 : Vector<RigidVector<Double, 3> >& VisBuffer::filluvw()
2434 : {
2435 14923 : CheckVisIter ();
2436 14923 : uvwOK_p = true;
2437 14923 : return visIter_p->uvw(uvw_p);
2438 : }
2439 :
2440 0 : Matrix<Double>& VisBuffer::filluvwMat()
2441 : {
2442 0 : CheckVisIter ();
2443 0 : uvwMatOK_p = true;
2444 0 : return visIter_p->uvwMat(uvwMat_p);
2445 : }
2446 :
2447 0 : Matrix<CStokesVector>& VisBuffer::fillVis(VisibilityIterator::DataColumn whichOne)
2448 : {
2449 0 : switch (whichOne) {
2450 0 : case VisibilityIterator::Model:
2451 0 : CheckVisIter1 (" (Model)");
2452 0 : modelVisibilityOK_p = true;
2453 0 : return visIter_p->visibility(modelVisibility_p, whichOne);
2454 : break;
2455 0 : case VisibilityIterator::Corrected:
2456 0 : CheckVisIter1 (" (Corrected)");
2457 0 : correctedVisibilityOK_p = true;
2458 0 : return visIter_p->visibility(correctedVisibility_p, whichOne);
2459 : break;
2460 0 : case VisibilityIterator::Observed:
2461 : default:
2462 0 : CheckVisIter1 (" (Observed)");
2463 0 : visibilityOK_p = true;
2464 0 : return visIter_p->visibility(visibility_p, whichOne);
2465 : break;
2466 : }
2467 : }
2468 :
2469 24036 : Cube<Complex>& VisBuffer::fillVisCube(VisibilityIterator::DataColumn whichOne)
2470 : {
2471 24036 : switch (whichOne) {
2472 2719 : case VisibilityIterator::Model:
2473 : {
2474 2719 : CheckVisIter1 (" (Model)");
2475 2719 : modelVisCubeOK_p = true;
2476 5438 : String modelkey; //=String("definedmodel_field_")+String::toString(fieldId());
2477 : Int snum;
2478 2719 : Bool hasmodkey = False;
2479 2719 : if (not visModelData_p.null()) hasmodkey = visModelData_p->isModelDefinedI(fieldId(), visIter_p->ms(), modelkey, snum);
2480 2719 : if( hasmodkey || !(visIter_p->ms().tableDesc().isColumn("MODEL_DATA")))
2481 : {
2482 :
2483 : //cerr << "HASMOD " << visModelData_p.hasModel(msId(), fieldId(), spectralWindow()) << endl;
2484 0 : visModelData_p->init(*this);
2485 : ////This bit can be removed when we do not support old model vesions saved
2486 0 : if(!(visModelData_p->isVersion2())){
2487 0 : if(visModelData_p->hasModel(msId(), fieldId(), spectralWindow()) ==-1){
2488 0 : if(hasmodkey){
2489 : //String whichrec=visIter_p->ms().keywordSet().asString(modelkey);
2490 0 : TableRecord modrec;
2491 0 : if(visModelData_p->getModelRecordI(modelkey, modrec, visIter_p->ms())){
2492 0 : visModelData_p->addModel(modrec, Vector<Int>(1, msId()), *this);
2493 : }
2494 : }
2495 : }
2496 : }
2497 : ////////////////////////////
2498 0 : visModelData_p->getModelVis(*this);
2499 : }
2500 : else
2501 : {
2502 2719 : visIter_p->visibility(modelVisCube_p, whichOne);
2503 : }
2504 : }
2505 2719 : return modelVisCube_p;
2506 : break;
2507 4678 : case VisibilityIterator::Corrected:
2508 4678 : CheckVisIter1 (" (Corrected)");
2509 4678 : correctedVisCubeOK_p = true;
2510 4678 : return visIter_p->visibility(correctedVisCube_p, whichOne);
2511 : break;
2512 16639 : case VisibilityIterator::Observed:
2513 : default:
2514 16639 : CheckVisIter1 (" (Observed)");
2515 16639 : visCubeOK_p = true;
2516 16639 : return visIter_p->visibility(visCube_p, whichOne);
2517 : break;
2518 : }
2519 : }
2520 :
2521 0 : Cube<Float>& VisBuffer::fillFloatDataCube()
2522 : {
2523 0 : CheckVisIter ();
2524 0 : floatDataCubeOK_p = true;
2525 0 : return visIter_p->floatData(floatDataCube_p);
2526 : }
2527 :
2528 6488 : Vector<Float>& VisBuffer::fillWeight()
2529 : {
2530 6488 : CheckVisIter ();
2531 6488 : weightOK_p = true;
2532 6488 : return visIter_p->weight(weight_p);
2533 : }
2534 :
2535 8953 : Matrix<Float>& VisBuffer::fillWeightMat()
2536 : {
2537 8953 : CheckVisIter ();
2538 8953 : weightMatOK_p = true;
2539 8953 : return visIter_p->weightMat(weightMat_p);
2540 : }
2541 :
2542 4705 : Cube<Float>& VisBuffer::fillWeightSpectrum()
2543 : {
2544 4705 : CheckVisIter ();
2545 4705 : weightSpectrumOK_p = true;
2546 4705 : return visIter_p->weightSpectrum(weightSpectrum_p);
2547 : }
2548 :
2549 : //Matrix<Float>& VisBuffer::fillImagingWeight()
2550 : //{
2551 : // CheckVisIter ();
2552 : // imagingWeightOK_p = true;
2553 : // return visIter_p->imagingWeight(imagingWeight_p);
2554 : //}
2555 :
2556 584678 : Vector<Float> VisBuffer::feed_pa(Double time) const
2557 : {
2558 584678 : return visIter_p->feed_pa(time);
2559 : }
2560 :
2561 0 : Float VisBuffer::parang0(Double time) const
2562 : {
2563 0 : return visIter_p->parang0(time);
2564 : }
2565 :
2566 0 : Vector<Float> VisBuffer::parang(Double time) const
2567 : {
2568 0 : return visIter_p->parang(time);
2569 : }
2570 :
2571 596 : Int VisBuffer::numberAnt () const
2572 : {
2573 596 : return msColumns().antenna().nrow(); /* for single (sub)array only.*/
2574 : }
2575 :
2576 0 : MDirection VisBuffer::azel0(Double time) const
2577 : {
2578 0 : return visIter_p->azel0(time);
2579 : }
2580 :
2581 0 : Vector<Double>& VisBuffer::azel0Vec(Double time, Vector<Double>& azelVec) const
2582 : {
2583 0 : MDirection azelMeas = This->azel0(time);
2584 0 : azelVec.resize(2);
2585 0 : azelVec = azelMeas.getAngle("deg").getValue();
2586 0 : return azelVec;
2587 : }
2588 :
2589 0 : Vector<MDirection> VisBuffer::azel(Double time) const
2590 : {
2591 0 : return visIter_p->azel(time);
2592 : }
2593 :
2594 0 : Matrix<Double>& VisBuffer::azelMat(Double time, Matrix<Double>& azelMat) const
2595 : {
2596 0 : Vector<MDirection> azelMeas = This->azel(time);
2597 0 : azelMat.resize(2, azelMeas.nelements());
2598 0 : for (uInt iant = 0; iant < azelMeas.nelements(); ++iant) {
2599 0 : azelMat.column(iant) = (azelMeas(iant).getAngle("deg").getValue());
2600 : }
2601 0 : return azelMat;
2602 :
2603 : }
2604 :
2605 0 : Double VisBuffer::hourang(Double time) const
2606 : {
2607 0 : return visIter_p->hourang(time);
2608 : }
2609 :
2610 0 : Vector<Int> VisBuffer::unique(const Vector<Int>& indices) const
2611 : {
2612 : // Filter integer index arrays for unique values
2613 : //
2614 0 : uInt n = indices.nelements();
2615 0 : Vector<Int> uniqIndices(n);
2616 0 : if (n > 0) {
2617 : // Sort temporary array in place
2618 0 : Vector<Int> sortedIndices = indices.copy();
2619 0 : GenSort<Int>::sort(sortedIndices);
2620 :
2621 : // Extract unique elements
2622 0 : uniqIndices(0) = sortedIndices(0);
2623 0 : uInt nUniq = 1;
2624 0 : for (uInt i = 1; i < n; i++) {
2625 0 : if (sortedIndices(i) != uniqIndices(nUniq - 1)) {
2626 0 : uniqIndices(nUniq++) = sortedIndices(i);
2627 : };
2628 : };
2629 0 : uniqIndices.resize(nUniq, true);
2630 : };
2631 0 : return uniqIndices;
2632 : }
2633 :
2634 : void
2635 6176 : VisBuffer::copyMsInfo (Int & msID, Bool & msOk, Bool & newMs) const
2636 : {
2637 6176 : msID = msId();
2638 6176 : msOk = msOK_p;
2639 6176 : newMs = newMS_p;
2640 6176 : }
2641 :
2642 :
2643 210909 : Bool VisBuffer::checkMSId()
2644 : {
2645 : //if this is not a new iteration then don't even check;
2646 : //Let the state be
2647 210909 : if (msOK_p) {
2648 193650 : return false;
2649 : }
2650 :
2651 17259 : if (visIter_p != static_cast<ROVisibilityIterator *> (0)) {
2652 :
2653 17259 : if (oldMSId_p != visIter_p->msId()) {
2654 :
2655 460 : oldMSId_p = visIter_p->msId();
2656 460 : newMS_p = true;
2657 : } else {
2658 16799 : newMS_p = false;
2659 : }
2660 :
2661 17259 : msOK_p = true;
2662 17259 : return newMS_p;
2663 :
2664 : } else {
2665 0 : return newMS_p;
2666 : }
2667 :
2668 : return false;
2669 : }
2670 :
2671 : VisBuffer *
2672 0 : VisBuffer::clone () const
2673 : {
2674 0 : return new VisBuffer (* this);
2675 : }
2676 :
2677 : void
2678 0 : VisBuffer::dirtyComponentsAdd (const VbDirtyComponents & dirtyComponents)
2679 : {
2680 0 : dirtyComponents_p = dirtyComponents_p + dirtyComponents;
2681 0 : }
2682 :
2683 : void
2684 0 : VisBuffer::dirtyComponentsAdd (VisBufferComponents::EnumType component)
2685 : {
2686 0 : dirtyComponents_p = dirtyComponents_p + VbDirtyComponents::singleton (component);
2687 0 : }
2688 :
2689 :
2690 : void
2691 0 : VisBuffer::dirtyComponentsClear ()
2692 : {
2693 0 : dirtyComponents_p = VbDirtyComponents::none();
2694 0 : }
2695 :
2696 : VbDirtyComponents
2697 0 : VisBuffer::dirtyComponentsGet () const
2698 : {
2699 0 : return dirtyComponents_p;
2700 : }
2701 :
2702 :
2703 : void
2704 0 : VisBuffer::dirtyComponentsSet (const VbDirtyComponents & dirtyComponents)
2705 : {
2706 0 : dirtyComponents_p = dirtyComponents;
2707 0 : }
2708 :
2709 : void
2710 0 : VisBuffer::dirtyComponentsSet (VisBufferComponents::EnumType component)
2711 : {
2712 0 : dirtyComponents_p = VbDirtyComponents::singleton (component);
2713 0 : }
2714 :
2715 0 : Bool VisBuffer::fetch(const asyncio::PrefetchColumns *pfc)
2716 : {
2717 0 : Bool success = true;
2718 :
2719 0 : for(asyncio::PrefetchColumns::const_iterator c = pfc->begin();
2720 0 : c != pfc->end(); ++c){
2721 0 : switch(*c){
2722 0 : case VisBufferComponents::Ant1:
2723 0 : This->antenna1();
2724 0 : break;
2725 0 : case VisBufferComponents::Ant2:
2726 0 : This->antenna2();
2727 0 : break;
2728 0 : case VisBufferComponents::ArrayId:
2729 0 : This->arrayId();
2730 0 : break;
2731 0 : case VisBufferComponents::Channel:
2732 0 : This->channel();
2733 0 : break;
2734 0 : case VisBufferComponents::CorrType:
2735 0 : This->corrType();
2736 0 : break;
2737 0 : case VisBufferComponents::Corrected:
2738 0 : This->correctedVisibility();
2739 0 : break;
2740 0 : case VisBufferComponents::CorrectedCube:
2741 0 : This->correctedVisCube();
2742 0 : break;
2743 0 : case VisBufferComponents::DataDescriptionId:
2744 0 : This->dataDescriptionId();
2745 0 : break;
2746 0 : case VisBufferComponents::Direction1:
2747 0 : This->direction1();
2748 0 : break;
2749 0 : case VisBufferComponents::Direction2:
2750 0 : This->direction2();
2751 0 : break;
2752 0 : case VisBufferComponents::Exposure:
2753 0 : This->exposure();
2754 0 : break;
2755 0 : case VisBufferComponents::Feed1:
2756 0 : This->feed1();
2757 0 : break;
2758 0 : case VisBufferComponents::Feed1_pa:
2759 0 : This->feed1_pa();
2760 0 : break;
2761 0 : case VisBufferComponents::Feed2:
2762 0 : This->feed2();
2763 0 : break;
2764 0 : case VisBufferComponents::Feed2_pa:
2765 0 : This->feed2_pa();
2766 0 : break;
2767 0 : case VisBufferComponents::FieldId:
2768 0 : This->fieldId();
2769 0 : break;
2770 0 : case VisBufferComponents::Flag:
2771 0 : This->flag();
2772 0 : break;
2773 0 : case VisBufferComponents::FlagCube:
2774 0 : This->flagCube();
2775 0 : break;
2776 0 : case VisBufferComponents::FlagCategory:
2777 0 : This->flagCategory();
2778 0 : break;
2779 0 : case VisBufferComponents::FlagRow:
2780 0 : This->flagRow();
2781 0 : break;
2782 0 : case VisBufferComponents::Freq:
2783 0 : This->frequency();
2784 0 : break;
2785 0 : case VisBufferComponents::ImagingWeight:
2786 : // This->imagingweight(); // do not fill this one
2787 0 : break;
2788 0 : case VisBufferComponents::Model:
2789 0 : This->modelVisibility();
2790 0 : break;
2791 0 : case VisBufferComponents::ModelCube:
2792 0 : This->modelVisCube();
2793 0 : break;
2794 0 : case VisBufferComponents::NChannel:
2795 0 : This->nChannel();
2796 0 : break;
2797 0 : case VisBufferComponents::NCorr:
2798 0 : This->nCorr();
2799 0 : break;
2800 0 : case VisBufferComponents::NRow:
2801 0 : This->nRow();
2802 0 : break;
2803 0 : case VisBufferComponents::ObservationId:
2804 0 : This->observationId();
2805 0 : break;
2806 0 : case VisBufferComponents::Observed:
2807 0 : This->visibility();
2808 0 : break;
2809 0 : case VisBufferComponents::ObservedCube:
2810 0 : This->visCube();
2811 0 : break;
2812 0 : case VisBufferComponents::PhaseCenter:
2813 0 : This->phaseCenter();
2814 0 : break;
2815 0 : case VisBufferComponents::PolFrame:
2816 0 : This->polFrame();
2817 0 : break;
2818 0 : case VisBufferComponents::ProcessorId:
2819 0 : This->processorId();
2820 0 : break;
2821 0 : case VisBufferComponents::Scan:
2822 0 : This->scan();
2823 0 : break;
2824 0 : case VisBufferComponents::Sigma:
2825 0 : This->sigma();
2826 0 : break;
2827 0 : case VisBufferComponents::SigmaMat:
2828 0 : This->sigmaMat();
2829 0 : break;
2830 0 : case VisBufferComponents::SpW:
2831 0 : This->spectralWindow();
2832 0 : break;
2833 0 : case VisBufferComponents::StateId:
2834 0 : This->stateId();
2835 0 : break;
2836 0 : case VisBufferComponents::Time:
2837 0 : This->time();
2838 0 : break;
2839 0 : case VisBufferComponents::TimeCentroid:
2840 0 : This->timeCentroid();
2841 0 : break;
2842 0 : case VisBufferComponents::TimeInterval:
2843 0 : This->timeInterval();
2844 0 : break;
2845 0 : case VisBufferComponents::Uvw:
2846 0 : This->uvw();
2847 0 : break;
2848 0 : case VisBufferComponents::UvwMat:
2849 0 : This->uvwMat();
2850 0 : break;
2851 0 : case VisBufferComponents::Weight:
2852 0 : This->weight();
2853 0 : break;
2854 0 : case VisBufferComponents::WeightMat:
2855 0 : This->weightMat();
2856 0 : break;
2857 0 : case VisBufferComponents::WeightSpectrum:
2858 0 : This->weightSpectrum();
2859 0 : break;
2860 0 : default:
2861 0 : throw(AipsError("Unrecognized column type"));
2862 : }
2863 : }
2864 0 : return success;
2865 : }
2866 :
2867 419 : VisBufferAutoPtr::VisBufferAutoPtr ()
2868 : {
2869 419 : visBuffer_p = NULL;
2870 419 : }
2871 :
2872 0 : VisBufferAutoPtr::VisBufferAutoPtr (VisBufferAutoPtr & other)
2873 : {
2874 : // Take ownership of the other's object
2875 :
2876 0 : visBuffer_p = other.visBuffer_p;
2877 0 : other.visBuffer_p = NULL;
2878 0 : }
2879 :
2880 0 : VisBufferAutoPtr::VisBufferAutoPtr (VisBuffer & vb)
2881 : {
2882 0 : constructVb (& vb);
2883 0 : }
2884 :
2885 0 : VisBufferAutoPtr::VisBufferAutoPtr (VisBuffer * vb)
2886 : {
2887 0 : constructVb (vb);
2888 0 : }
2889 :
2890 419 : VisBufferAutoPtr::VisBufferAutoPtr (ROVisibilityIterator & rovi)
2891 : {
2892 419 : construct (& rovi, true);
2893 419 : }
2894 :
2895 :
2896 0 : VisBufferAutoPtr::VisBufferAutoPtr (ROVisibilityIterator * rovi)
2897 : {
2898 0 : construct (rovi, true);
2899 0 : }
2900 :
2901 1676 : VisBufferAutoPtr::~VisBufferAutoPtr ()
2902 : {
2903 838 : delete visBuffer_p;
2904 838 : }
2905 :
2906 : VisBufferAutoPtr &
2907 157 : VisBufferAutoPtr::operator= (VisBufferAutoPtr & other)
2908 : {
2909 157 : if (this != & other){
2910 :
2911 157 : delete visBuffer_p; // release any currently referenced object
2912 :
2913 : // Take ownership of the other's object
2914 :
2915 157 : visBuffer_p = other.visBuffer_p;
2916 157 : other.visBuffer_p = NULL;
2917 : }
2918 :
2919 157 : return * this;
2920 : }
2921 :
2922 : VisBuffer &
2923 6224 : VisBufferAutoPtr::operator* () const
2924 : {
2925 6224 : assert (visBuffer_p != NULL);
2926 :
2927 6224 : return * visBuffer_p;
2928 : }
2929 :
2930 : VisBuffer *
2931 18162 : VisBufferAutoPtr::operator-> () const
2932 : {
2933 18162 : assert (visBuffer_p != NULL);
2934 :
2935 18162 : return visBuffer_p;
2936 : }
2937 :
2938 : void
2939 576 : VisBufferAutoPtr::construct (ROVisibilityIterator * rovi, Bool attachVi)
2940 : {
2941 576 : if (rovi->isAsynchronous ()){
2942 :
2943 : // Create an asynchronous VisBuffer
2944 :
2945 : VisBufferAsyncWrapper * vba;
2946 :
2947 0 : if (attachVi){
2948 0 : vba = new VisBufferAsyncWrapper (* rovi);
2949 : }
2950 : else{
2951 0 : vba = new VisBufferAsyncWrapper ();
2952 : }
2953 :
2954 0 : visBuffer_p = vba;
2955 : }
2956 : else{
2957 :
2958 : // This is a synchronous VI so just create a synchronous VisBuffer.
2959 :
2960 576 : if (attachVi){
2961 419 : visBuffer_p = new VisBuffer (* rovi);
2962 : }
2963 : else{
2964 157 : visBuffer_p = new VisBuffer ();
2965 : }
2966 : }
2967 576 : }
2968 :
2969 : void
2970 0 : VisBufferAutoPtr::constructVb (VisBuffer * vb)
2971 : {
2972 0 : VisBufferAsync * vba = dynamic_cast<VisBufferAsync *> (vb);
2973 :
2974 0 : if (vba != NULL){
2975 :
2976 : // Create an asynchronous VisBuffer
2977 :
2978 0 : VisBufferAsyncWrapper * vbaNew = new VisBufferAsyncWrapper (* vba);
2979 :
2980 0 : visBuffer_p = vbaNew;
2981 : }
2982 : else{
2983 :
2984 : // This is a synchronous VI so just create a synchronous VisBuffer.
2985 :
2986 0 : visBuffer_p = new VisBuffer (* vb);
2987 : }
2988 0 : }
2989 :
2990 : VisBuffer *
2991 157 : VisBufferAutoPtr::get () const
2992 : {
2993 157 : return visBuffer_p;
2994 : }
2995 :
2996 :
2997 :
2998 :
2999 : VisBuffer *
3000 0 : VisBufferAutoPtr::release ()
3001 : {
3002 0 : VisBuffer * result = visBuffer_p;
3003 0 : visBuffer_p = NULL;
3004 :
3005 0 : return result;
3006 : }
3007 :
3008 :
3009 : void
3010 0 : VisBufferAutoPtr::set (VisBuffer & vb)
3011 : {
3012 0 : delete visBuffer_p;
3013 0 : visBuffer_p = & vb;
3014 0 : }
3015 :
3016 : void
3017 0 : VisBufferAutoPtr::set (VisBuffer * vb)
3018 : {
3019 0 : delete visBuffer_p;
3020 0 : visBuffer_p = vb;
3021 0 : }
3022 :
3023 : void
3024 157 : VisBufferAutoPtr::set (ROVisibilityIterator * rovi, Bool attachIt)
3025 : {
3026 157 : delete visBuffer_p;
3027 157 : construct (rovi, attachIt);
3028 157 : }
3029 :
3030 : void
3031 0 : VisBufferAutoPtr::set (ROVisibilityIterator & rovi, Bool attachIt)
3032 : {
3033 0 : set (& rovi, attachIt);
3034 0 : }
3035 :
3036 : const VbDirtyComponents VbDirtyComponents::all_p = initializeAll ();
3037 :
3038 : VbDirtyComponents
3039 0 : VbDirtyComponents::operator+ (const VbDirtyComponents & other) const
3040 : {
3041 0 : VbDirtyComponents result = * this;
3042 :
3043 0 : result.set_p.insert (other.begin(), other.end());
3044 :
3045 0 : return result;
3046 : }
3047 :
3048 :
3049 :
3050 : VbDirtyComponents::const_iterator
3051 0 : VbDirtyComponents::begin () const
3052 : {
3053 0 : return set_p.begin();
3054 : }
3055 :
3056 : Bool
3057 0 : VbDirtyComponents::contains (VisBufferComponents::EnumType component) const
3058 : {
3059 0 : return utilj::containsKey (component, set_p);
3060 : }
3061 :
3062 : VbDirtyComponents::const_iterator
3063 0 : VbDirtyComponents::end () const
3064 : {
3065 0 : return set_p.end();
3066 : }
3067 :
3068 : VbDirtyComponents
3069 0 : VbDirtyComponents::all ()
3070 : {
3071 0 : return all_p;
3072 : }
3073 :
3074 : VbDirtyComponents
3075 0 : VbDirtyComponents::exceptThese (VisBufferComponents::EnumType component, ...)
3076 : {
3077 : va_list vaList;
3078 :
3079 0 : va_start (vaList, component);
3080 :
3081 0 : VisBufferComponents::EnumType c = component;
3082 0 : VbDirtyComponents dirtyComponents = all();
3083 :
3084 0 : while (c != VisBufferComponents::Unknown){
3085 :
3086 0 : ThrowIf (! all().contains (c), "Not a writable VB component: " + String::toString (c));
3087 :
3088 0 : dirtyComponents.set_p.erase (c);
3089 0 : c = (VisBufferComponents::EnumType) va_arg (vaList, Int);
3090 : }
3091 :
3092 0 : va_end (vaList);
3093 :
3094 0 : return dirtyComponents;
3095 :
3096 : }
3097 :
3098 : VbDirtyComponents
3099 1 : VbDirtyComponents::initializeAll ()
3100 : {
3101 :
3102 1 : VbDirtyComponents all;
3103 :
3104 : VisBufferComponents::EnumType
3105 1 : writableComponents [] = {VisBufferComponents::Corrected,
3106 : VisBufferComponents::CorrectedCube,
3107 : VisBufferComponents::Flag,
3108 : VisBufferComponents::FlagCube,
3109 : VisBufferComponents::FlagRow,
3110 : VisBufferComponents::Model,
3111 : VisBufferComponents::ModelCube,
3112 : VisBufferComponents::Observed,
3113 : VisBufferComponents::ObservedCube,
3114 : VisBufferComponents::Sigma,
3115 : VisBufferComponents::SigmaMat,
3116 : VisBufferComponents::Weight,
3117 : VisBufferComponents::WeightMat,
3118 : VisBufferComponents::Unknown};
3119 :
3120 1 : for (Int i = 0; ; i++){
3121 :
3122 14 : if (writableComponents [i] == VisBufferComponents::Unknown){
3123 1 : break;
3124 : }
3125 :
3126 13 : all.set_p.insert (writableComponents [i]);
3127 : }
3128 :
3129 2 : return all;
3130 : }
3131 :
3132 : VbDirtyComponents
3133 0 : VbDirtyComponents::none ()
3134 : {
3135 0 : return VbDirtyComponents ();
3136 : }
3137 :
3138 : VbDirtyComponents
3139 0 : VbDirtyComponents::singleton (VisBufferComponents::EnumType component)
3140 : {
3141 0 : ThrowIf (! all().contains (component), "Not a writable VB component.");
3142 0 : VbDirtyComponents result;
3143 0 : result.set_p.insert (component);
3144 :
3145 0 : return result;
3146 : }
3147 :
3148 : VbDirtyComponents
3149 0 : VbDirtyComponents::these (VisBufferComponents::EnumType component, ...)
3150 : {
3151 : va_list vaList;
3152 :
3153 0 : va_start (vaList, component);
3154 :
3155 0 : VisBufferComponents::EnumType c = component;
3156 0 : VbDirtyComponents dirtyComponents;
3157 :
3158 0 : while (c != VisBufferComponents::Unknown){
3159 :
3160 0 : ThrowIf (! all().contains (c), "Not a writable VB component: " + String::toString (c));
3161 :
3162 0 : dirtyComponents.set_p.insert (c);
3163 0 : c = (VisBufferComponents::EnumType ) va_arg (vaList, Int);
3164 : }
3165 :
3166 0 : va_end (vaList);
3167 :
3168 0 : return dirtyComponents;
3169 : }
3170 :
3171 :
3172 :
3173 :
3174 : } //# NAMESPACE CASA - END
3175 :
|
