Line data Source code
1 : #include <synthesis/Utilities/FixVis.h>
2 : #include <msvis/MSVis/MSUVWGenerator.h>
3 : #include <msvis/MSVis/SubMS.h>
4 : #include <casacore/measures/Measures/MeasTable.h>
5 : #include <casacore/measures/Measures/UVWMachine.h>
6 : #include <casacore/casa/Logging/LogIO.h>
7 : #include <casacore/casa/Exceptions/Error.h>
8 : #include <casacore/casa/Quanta/MVAngle.h>
9 : #include <casacore/coordinates/Coordinates/Projection.h>
10 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
11 : #include <casacore/coordinates/Coordinates/ObsInfo.h>
12 : #include <casacore/images/Images/PagedImage.h> // Used to pass coords
13 : #include <casacore/images/Images/ImageInfo.h> // to FTMachine.
14 : #include <casacore/ms/MeasurementSets/MSColumns.h>
15 : #include <casacore/ms/MeasurementSets/MSDopplerUtil.h>
16 : #include <casacore/ms/MSSel/MSSelection.h>
17 : #include <casacore/ms/MSSel/MSSelectionTools.h>
18 : #include <msvis/MSVis/VisibilityIterator.h>
19 : #include <msvis/MSVis/VisBuffer.h>
20 : #include <casacore/casa/BasicSL/String.h> // for parseColumnNames()
21 :
22 : #include <iostream>
23 :
24 : using namespace casacore;
25 : namespace casa {
26 :
27 32 : FixVis::FixVis(MeasurementSet& ms, const String& dataColName) :
28 : FTMachine(),
29 32 : ms_p(ms),
30 32 : msc_p(NULL),
31 32 : nsel_p(0),
32 32 : nAllFields_p(1),
33 32 : npix_p(2),
34 32 : cimageShape_p(4, npix_p, npix_p, 1, 1), // Can we get away with
35 32 : tileShape_p(4, npix_p, npix_p, 1, 1), // (1, 1, 1, 1)? Does it matter?
36 32 : tiledShape_p(cimageShape_p, tileShape_p),
37 32 : antennaSel_p(false),
38 64 : freqFrameValid_p(false)
39 : // obsString_p("")
40 : {
41 32 : logSink() << LogOrigin("FixVis", "") << LogIO::NORMAL3;
42 :
43 32 : antennaId_p.resize();
44 32 : antennaSelStr_p.resize();
45 32 : distances_p.resize();
46 32 : dataCols_p = SubMS::parseColumnNames(dataColName, ms);
47 32 : nDataCols_p = dataCols_p.nelements();
48 :
49 32 : nchan_p = 1; // imageNchan_p;
50 :
51 32 : spectralwindowids_p.resize(ms_p.spectralWindow().nrow());
52 32 : indgen(spectralwindowids_p);
53 :
54 32 : lockCounter_p = 0;
55 32 : }
56 :
57 : // Destructor
58 128 : FixVis::~FixVis()
59 : {
60 32 : if(!ms_p.isNull())
61 32 : ms_p.flush();
62 :
63 32 : delete msc_p;
64 32 : }
65 :
66 : // Interpret field indices (MSSelection)
67 0 : Vector<Int> FixVis::getFieldIdx(const String& fields)
68 : {
69 0 : MSSelection mssel;
70 :
71 0 : mssel.setFieldExpr(fields);
72 0 : return mssel.getFieldList(&ms_p);
73 0 : }
74 :
75 32 : uInt FixVis::setFields(const Vector<Int>& fieldIds)
76 : {
77 32 : logSink() << LogOrigin("FixVis", "setFields");
78 32 : logSink() << LogIO::NORMAL << "Selecting fields ";
79 :
80 32 : nsel_p = fieldIds.nelements();
81 32 : nAllFields_p = ms_p.field().nrow();
82 32 : FieldIds_p.resize(nAllFields_p);
83 :
84 110 : for(Int i = 0; i < static_cast<Int>(nAllFields_p); ++i){
85 78 : FieldIds_p(i) = -1;
86 130 : for(uInt j = 0; j < nsel_p; ++j){
87 90 : if(fieldIds[j] == i){
88 38 : FieldIds_p(i) = i;
89 38 : logSink() << i << " " << LogIO::NORMAL;
90 38 : break;
91 : }
92 : }
93 : }
94 32 : logSink() << LogIO::POST;
95 :
96 32 : return nsel_p;
97 : }
98 :
99 15 : void FixVis::setPhaseDirs(const Vector<MDirection>& phaseDirs)
100 : {
101 15 : phaseDirs_p = phaseDirs;
102 :
103 : // Do a consistency check between fieldIds and FieldIds_p.
104 15 : logSink() << LogOrigin("FixVis", "setPhaseDirs");
105 15 : uInt n2change = phaseDirs.nelements();
106 15 : if(n2change != nsel_p){
107 : logSink() << LogIO::SEVERE
108 : << "Mismatch between the number of new directions and the fields to change"
109 0 : << LogIO::POST;
110 : }
111 15 : }
112 :
113 :
114 1 : void FixVis::convertFieldDirs(const MDirection::Types outType)
115 : {
116 1 : logSink() << LogOrigin("FixVis", "convertFieldDirs");
117 :
118 : // Note that the each direction column in the FIELD table only allows one
119 : // reference frame for the entire column, but polynomials can be assigned on
120 : // a row-wise basis for objects moving in that frame.
121 :
122 : Muvw::Types uvwtype;
123 1 : Muvw::getType(uvwtype, MDirection::showType(outType));
124 :
125 1 : ready_msc_p();
126 :
127 1 : msc_p->uvw().rwKeywordSet().asrwRecord("MEASINFO").define("Ref", MDirection::showType(outType));
128 :
129 1 : MSFieldColumns& msfcs = msc_p->field();
130 :
131 1 : MDirection pd0(msfcs.phaseDirMeas(0));
132 : logSink() << LogIO::DEBUG1
133 1 : << "PHASE_DIR[0] before = " << pd0.tellMe() << " ";
134 : {
135 1 : ostringstream os;
136 1 : os << *(pd0.getData()) << endl;
137 1 : logSink() << os.str() << LogIO::POST;
138 1 : }
139 :
140 : // There is no physical or known programming need to change the delay and
141 : // reference direction frames as well, but for aesthetic reasons we keep them
142 : // all in the same frame if they start in the same frame.
143 : //ArrayMeasColumn<MDirection> msDelayDirCol;
144 : //msDelayDirCol.reference(msfcs.delayDirMeasCol());
145 : //ArrayMeasColumn<MDirection> msReferenceDirCol;
146 : //msReferenceDirCol.reference(msfcs.referenceDirMeasCol());
147 1 : Bool doAll3 = (msfcs.phaseDirMeasCol().getMeasRef().getType() ==
148 2 : msfcs.delayDirMeasCol().getMeasRef().getType() &&
149 1 : msfcs.phaseDirMeasCol().getMeasRef().getType() ==
150 1 : msfcs.referenceDirMeasCol().getMeasRef().getType());
151 :
152 : // Setup conversion machines.
153 : // Set the frame - choose the first antenna. For e.g. VLBI, we
154 : // will need to reset the frame per antenna
155 1 : mLocation_p = msc_p->antenna().positionMeas()(0);
156 :
157 : // Expect problems if a moving object is involved!
158 1 : mFrame_p = MeasFrame(msfcs.timeMeas()(0), mLocation_p);
159 :
160 : //MDirection::Ref startref(msfcs.phaseDirMeasCol().getMeasRef());
161 : // If the either of the start or destination frames refers to a finite
162 : // distance, then the conversion has to be done in two steps:
163 : // MDirection::Convert start2app(msfcs.phaseDirMeasCol()(0), MDirection::APP);
164 : //
165 : // Otherwise the conversion can be done directly.
166 2 : Bool haveMovingFrame = (MDirection::castType(msfcs.phaseDirMeasCol().getMeasRef().getType()) > MDirection::N_Types ||
167 1 : outType > MDirection::N_Types);
168 :
169 : const MDirection::Ref newFrame(haveMovingFrame ? MDirection::APP : outType,
170 1 : mFrame_p);
171 :
172 1 : convertFieldCols(msfcs, newFrame, doAll3);
173 :
174 1 : if(haveMovingFrame){
175 : // Since ArrayMeasCol most likely uses one frame for the whole column, do
176 : // the second half of the conversion with a second pass through the column
177 : // instead of on a row-by-row basis.
178 : logSink() << LogIO::WARN
179 : << "Switching to or from accelerating frames is not well tested."
180 0 : << LogIO::POST;
181 :
182 : // Using msfcs.phaseDirMeasCol()(0)[0] to initialize converter will only
183 : // work if msfcs.phaseDirMeasCol()'s type has been set to APP.
184 0 : const MDirection::Ref newerFrame(outType, mFrame_p);
185 :
186 0 : convertFieldCols(msfcs, newerFrame, doAll3);
187 0 : }
188 :
189 1 : pd0 = msfcs.phaseDirMeas(0);
190 : logSink() << LogIO::DEBUG1
191 1 : << "PHASE_DIR[0] after = " << pd0.tellMe() << " ";
192 : {
193 1 : ostringstream os;
194 1 : os << *(pd0.getData()) << endl;
195 1 : logSink() << os.str() << LogIO::POST;
196 1 : }
197 1 : }
198 :
199 :
200 1 : void FixVis::convertFieldCols(MSFieldColumns& msfcs,
201 : const MDirection::Ref& newFrame,
202 : const Bool doAll3)
203 : {
204 1 : logSink() << LogOrigin("FixVis", "convertFieldCols");
205 : // Unfortunately ArrayMeasColumn::doConvert() is private, which moots the
206 : // point of making a conversion machine here.
207 : // Vector<MDirection> dummyV;
208 : // dummyV.assign(pdc(0));
209 : // MDirection::Convert *converter = new MDirection::Convert(dummyV[0],
210 : // newFrame);
211 : // if(!converter)
212 : // logSink() << "Cannot make direction conversion machine"
213 : // << LogIO::SEVERE;
214 :
215 1 : uInt nrows = msfcs.nrow();
216 :
217 : // Convert each phase tracking center. This will make them numerically
218 : // correct in the new frame, but the column will still be labelled with the
219 : // old frame.
220 :
221 : uInt nOrders;
222 1 : Vector<MDirection> mdarr; // direction for each order
223 1 : Array<Double> darr; // longitude and latitude for each order
224 1 : Vector<Double> dirV;
225 2 : for(uInt i = 0; i < nrows; ++i){
226 1 : nOrders = msfcs.numPoly()(i) + 1;
227 : logSink() << LogIO::DEBUG1 << "numPoly(" << i << ") = " << nOrders - 1
228 1 : << LogIO::POST;
229 :
230 : //pdc.put(i, pdc.doConvert(i, *converter));
231 1 : mdarr = msfcs.phaseDirMeasCol().convert(i, newFrame);
232 1 : darr.resize(IPosition(2, 2, nOrders));
233 2 : for(uInt orderNumber = 0; orderNumber < nOrders; ++orderNumber){
234 1 : dirV = mdarr[orderNumber].getAngle().getValue();
235 1 : darr(IPosition(2, 0, orderNumber)) = dirV[0];
236 1 : darr(IPosition(2, 1, orderNumber)) = dirV[1];
237 : }
238 1 : msfcs.phaseDir().put(i, darr);
239 :
240 : //msfcs.phaseDirMeasCol().put(i, mdarr);
241 :
242 1 : if(doAll3){
243 : //ddc.put(i, ddc.doConvert(i, *converter));
244 1 : mdarr = msfcs.delayDirMeasCol().convert(i, newFrame);
245 2 : for(uInt orderNumber = 0; orderNumber < nOrders; ++orderNumber){
246 1 : dirV = mdarr[orderNumber].getAngle().getValue();
247 1 : darr(IPosition(2, 0, orderNumber)) = dirV[0];
248 1 : darr(IPosition(2, 1, orderNumber)) = dirV[1];
249 : }
250 1 : msfcs.delayDir().put(i, darr);
251 : //rdc.put(i, rdc.doConvert(i, *converter));
252 : //msfcs.referenceDirMeasCol().put(i,
253 : // msfcs.referenceDirMeasCol().convert(i, newFrame));
254 1 : mdarr = msfcs.referenceDirMeasCol().convert(i, newFrame);
255 2 : for(uInt orderNumber = 0; orderNumber < nOrders; ++orderNumber){
256 1 : dirV = mdarr(IPosition(1, orderNumber)).getAngle().getValue();
257 1 : darr(IPosition(2, 0, orderNumber)) = dirV[0];
258 1 : darr(IPosition(2, 1, orderNumber)) = dirV[1];
259 : }
260 1 : msfcs.referenceDir().put(i, darr);
261 : }
262 : }
263 :
264 : // Update the reference frame label(s).
265 1 : msfcs.phaseDirMeasCol().setDescRefCode(newFrame.getType(), false);
266 1 : if(doAll3){
267 1 : msfcs.delayDirMeasCol().setDescRefCode(newFrame.getType(), false);
268 1 : msfcs.referenceDirMeasCol().setDescRefCode(newFrame.getType(), false);
269 : }
270 :
271 : //delete converter;
272 1 : }
273 :
274 :
275 15 : void FixVis::setDistances(const Vector<Double>& distances)
276 : {
277 15 : logSink() << LogOrigin("FixVis", "setDistances");
278 15 : if(distances.nelements() != nsel_p)
279 : logSink() << LogIO::SEVERE
280 : << "Mismatch between the # of specified distances and selected fields."
281 0 : << LogIO::POST;
282 15 : distances_p = distances;
283 15 : }
284 :
285 : // Calculate the (u, v, w)s and store them in ms_p.
286 17 : Bool FixVis::calc_uvw(const String& refcode, const Bool reuse)
287 : {
288 17 : Bool retval = false;
289 :
290 17 : logSink() << LogOrigin("FixVis", "calc_uvw");
291 :
292 17 : if(!ready_msc_p())
293 0 : return false;
294 :
295 17 : if(nsel_p > 0){
296 :
297 : // Get the PHASE_DIR reference frame type for the input ms.
298 17 : MSFieldColumns& msfcs(msc_p->field());
299 17 : MDirection startDir = msfcs.phaseDirMeas(0);
300 17 : MDirection::Types startDirType = MDirection::castType(msfcs.phaseDirMeasCol().getMeasRef().getType());
301 : MDirection::Types wantedDirType;
302 17 : MDirection::getType(wantedDirType, refcode);
303 :
304 17 : if(startDirType != wantedDirType){
305 1 : if(nsel_p < nAllFields_p){
306 : logSink() << LogIO::SEVERE
307 : << "The reference frame must either be changed for all fields or not at all."
308 0 : << LogIO::POST;
309 0 : return false;
310 : }
311 : else
312 1 : convertFieldDirs(wantedDirType);
313 : }
314 16 : else if(reuse){
315 : logSink() << LogIO::NORMAL
316 : << "The UVWs are already in the desired frame - leaving them as is."
317 0 : << LogIO::POST;
318 0 : return true;
319 : }
320 :
321 : try{
322 17 : if(reuse){
323 1 : const MDirection::Ref outref(wantedDirType);
324 :
325 1 : rotateUVW(startDir, outref);
326 1 : }
327 : else{
328 : Muvw::Types uvwtype;
329 : MBaseline::Types bltype;
330 :
331 : try{
332 16 : MBaseline::getType(bltype, refcode);
333 16 : Muvw::getType(uvwtype, refcode);
334 : }
335 0 : catch(AipsError x){
336 : logSink() << LogIO::SEVERE
337 : << "refcode \"" << refcode << "\" is not valid for baselines."
338 0 : << LogIO::POST;
339 0 : return false;
340 0 : }
341 :
342 16 : MSUVWGenerator uvwgen(*msc_p, bltype, uvwtype);
343 16 : retval = uvwgen.make_uvws(FieldIds_p);
344 16 : }
345 :
346 : }
347 0 : catch(AipsError x){
348 0 : logSink() << LogIO::SEVERE << "Error " << x.getMesg() << LogIO::POST;
349 0 : }
350 :
351 17 : }
352 : else{
353 : logSink() << LogIO::SEVERE
354 : << "There is a problem with the selected field IDs and phase tracking centers."
355 0 : << LogIO::POST;
356 : }
357 17 : return retval;
358 : }
359 :
360 : // Convert the UVW column to a new reference frame by rotating the old
361 : // baselines instead of calculating fresh ones.
362 : //
363 : // oldref must be supplied instead of extracted from msc_p->uvw(), because
364 : // the latter might be wrong (use the field direction).
365 1 : void FixVis::rotateUVW(const MDirection &indir, const MDirection::Ref& newref)
366 : {
367 1 : ArrayColumn<Double>& UVWcol = msc_p->uvw();
368 :
369 : // Setup a machine for converting a UVW vector from the old frame to
370 : // uvwtype's frame
371 1 : UVWMachine uvm(newref, indir);
372 1 : RotMatrix rm(uvm.rotationUVW());
373 :
374 1 : uInt nRows = UVWcol.nrow();
375 41419 : for(uInt row = 0; row < nRows; ++row){
376 41418 : UVWcol.put(row, (rm * MVuvw(UVWcol(row))).getVector());
377 : }
378 2 : return;
379 1 : }
380 :
381 : // Don't just calculate the (u, v, w)s, do everything and store them in ms_p.
382 15 : Bool FixVis::fixvis(const String& refcode)
383 : {
384 15 : logSink() << LogOrigin("FixVis", "fixvis");
385 :
386 15 : Bool retval = false;
387 :
388 15 : if(!ready_msc_p())
389 0 : return false;
390 :
391 15 : if(nsel_p > 0){
392 15 : if(phaseDirs_p.nelements() == static_cast<uInt>(nsel_p)){
393 :
394 15 : String telescop = msc_p->observation().telescopeName()(0);
395 15 : MPosition obsPosition;
396 15 : if(!(MeasTable::Observatory(obsPosition, telescop))){
397 : logSink() << LogIO::WARN << "Did not get the position of " << telescop
398 0 : << " from data repository" << LogIO::POST;
399 : logSink() << LogIO::WARN << "Please contact CASA to add it to the repository."
400 0 : << LogIO::POST;
401 0 : logSink() << LogIO::WARN << "Frequency conversion will not work " << LogIO::POST;
402 0 : freqFrameValid_p = false;
403 : }
404 : else{
405 15 : mLocation_p = obsPosition;
406 15 : freqFrameValid_p = true;
407 : }
408 :
409 15 : MSFieldColumns& msfcs = msc_p->field();
410 :
411 15 : mFrame_p = MeasFrame(msfcs.timeMeas()(0), mLocation_p);
412 :
413 15 : msc_p->uvw().rwKeywordSet().asrwRecord("MEASINFO").define("Ref", refcode);
414 :
415 : //**** Adjust the phase tracking centers and distances. ****
416 : // Loop through each selected field.
417 : Int selectedField;
418 15 : Int selFldCounter=0;
419 45 : for(uInt fldCounter = 0; fldCounter < nAllFields_p; ++fldCounter){
420 30 : selectedField = FieldIds_p[fldCounter];
421 30 : if(selectedField<0){
422 15 : continue;
423 : }
424 15 : setImageField(selectedField);
425 15 : if(makeSelection(selectedField)){
426 15 : logSink() << LogIO::NORMAL << "Processing field " << selectedField << LogIO::POST;
427 15 : processSelected(selFldCounter);
428 15 : selFldCounter++;
429 :
430 : // Update FIELD (and/or optional tables SOURCE, OBSERVATION, but not
431 : // POINTING?) to new PTC.
432 :
433 15 : retval = true;
434 : }
435 : else{
436 : logSink() << LogIO::SEVERE
437 : << "Field " << selectedField << " could not be selected for phase tracking center or"
438 : << " distance adjustment."
439 0 : << LogIO::POST;
440 : }
441 : }
442 15 : }
443 0 : else if(phaseDirs_p.nelements() > 0){
444 : logSink() << LogIO::SEVERE
445 : << "There is a problem with the selected field IDs and phase tracking centers.\n"
446 : << "No adjustments of phase tracking centers or distances will be done."
447 0 : << LogIO::POST;
448 : }
449 : }
450 : else{
451 0 : logSink() << LogIO::SEVERE << "No fields are selected." << LogIO::POST;
452 : }
453 15 : return retval;
454 : }
455 :
456 :
457 15 : Bool FixVis::setImageField(const Int fieldid,
458 : const Bool dotrackDir //, const MDirection& trackDir
459 : )
460 : {
461 15 : logSink() << LogOrigin("FixVis", "setImageField()");
462 :
463 : try{
464 :
465 15 : doTrackSource_p = dotrackDir;
466 :
467 15 : fieldid_p = fieldid;
468 :
469 15 : return true;
470 : }
471 : catch(AipsError x){
472 : this->unlock();
473 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
474 : << LogIO::POST;
475 : return false;
476 : }
477 : return true;
478 : }
479 :
480 0 : Bool FixVis::lock()
481 : {
482 0 : Bool ok = true;
483 :
484 0 : if(lockCounter_p == 0)
485 0 : ok = ms_p.lock();
486 0 : ++lockCounter_p;
487 :
488 0 : return ok;
489 : }
490 :
491 0 : void FixVis::unlock()
492 : {
493 0 : if(lockCounter_p == 1)
494 0 : ms_p.unlock();
495 :
496 0 : if(lockCounter_p > 0)
497 0 : --lockCounter_p;
498 0 : }
499 :
500 :
501 15 : Bool FixVis::makeSelection(const Int selectedField)
502 : {
503 15 : logSink() << LogOrigin("FixVis", "makeSelection()");
504 :
505 : //Vis/MSIter will check if SORTED_TABLE exists and resort if necessary.
506 15 : MSSelection thisSelection;
507 15 : if(selectedField >= 0 && nAllFields_p > 1){
508 15 : Vector<Int> wrapper;
509 15 : wrapper.resize(1);
510 15 : wrapper[0] = selectedField;
511 15 : thisSelection.setFieldExpr(MSSelection::indexExprStr(wrapper));
512 15 : }
513 0 : else if(antennaSel_p){
514 0 : if(antennaId_p.nelements() > 0)
515 0 : thisSelection.setAntennaExpr(MSSelection::indexExprStr(antennaId_p));
516 0 : if(antennaSelStr_p[0] != "")
517 0 : thisSelection.setAntennaExpr(MSSelection::nameExprStr(antennaSelStr_p));
518 : }
519 : // if(obsString_p != "")
520 : // thisSelection.setObservationExpr(obsString_p);
521 :
522 15 : TableExprNode exprNode = thisSelection.toTableExprNode(&ms_p);
523 :
524 : // Now remake the selected ms
525 15 : if(!(exprNode.isNull())){
526 15 : mssel_p = MeasurementSet(ms_p(exprNode));
527 : }
528 0 : else if(selectedField < 0 || nsel_p == nAllFields_p){
529 : // Null take all the ms ...setdata() blank means that
530 0 : mssel_p = MeasurementSet(ms_p);
531 : }
532 : else{
533 : logSink() << LogIO::SEVERE
534 : << "Error selecting field " << selectedField << ". "
535 : << "Are you trying to adjust antenna positions at the same time?"
536 0 : << LogIO::POST;
537 0 : return false;
538 : }
539 :
540 : //mssel_p.rename(ms_p.tableName()+"/SELECTED_TABLE", Table::Scratch);
541 15 : if(mssel_p.nrow() == 0)
542 0 : return false;
543 :
544 15 : if(mssel_p.nrow() < ms_p.nrow())
545 : logSink() << LogIO::NORMAL
546 : << mssel_p.nrow() << " rows selected out of " << ms_p.nrow() << "."
547 15 : << LogIO::POST;
548 :
549 15 : delete msc_p;
550 15 : msc_p = NULL;
551 15 : return ready_msc_p();
552 15 : }
553 :
554 48 : Bool FixVis::ready_msc_p()
555 : {
556 48 : Bool retval = false;
557 :
558 48 : if(!msc_p){
559 : try{
560 47 : msc_p = new MSColumns(mssel_p.isNull() ? ms_p : mssel_p);
561 47 : retval = true; // Assume msc_p is OK.
562 : }
563 0 : catch(AipsError& e){
564 : logSink() << LogIO::SEVERE
565 : << "Error getting the columns from the MS."
566 0 : << LogIO::POST;
567 0 : }
568 0 : catch(std::bad_alloc& e){
569 : logSink() << LogIO::SEVERE
570 : << "Error allocating memory for the MS columns."
571 0 : << LogIO::POST;
572 0 : }
573 : // Just let any other exceptions, of the unexpected kind,
574 : // propagate up where they can be seen.
575 : }
576 : else
577 1 : retval = true; // Assume msc_p is OK.
578 :
579 48 : return retval;
580 : }
581 :
582 :
583 15 : void FixVis::processSelected(uInt numInSel)
584 : {
585 15 : logSink() << LogOrigin("FixVis", "processSelected()");
586 :
587 15 : mImage_p = phaseDirs_p[numInSel];
588 :
589 15 : ArrayColumn<Double>& UVWcol = msc_p->uvw();
590 :
591 15 : Block<Int> sort(0);
592 15 : sort.resize(4);
593 15 : sort[0] = MS::ARRAY_ID; // use default sort order
594 15 : sort[1] = MS::FIELD_ID;
595 15 : sort[2] = MS::DATA_DESC_ID;
596 15 : sort[3] = MS::TIME;
597 :
598 15 : VisibilityIterator vi(mssel_p, sort);
599 :
600 : // Loop over all visibilities in the selected field.
601 15 : VisBuffer vb(vi);
602 :
603 15 : vi.origin();
604 :
605 30 : for(vi.originChunks(); vi.moreChunks(); vi.nextChunk()){
606 2715 : for(vi.origin(); vi.more(); ++vi){
607 :
608 2700 : uInt numRows = vb.nRow();
609 :
610 2700 : Matrix<Double> uvw(3, numRows);
611 2700 : uvw=0.0;
612 2700 : Vector<Double> dphase(numRows);
613 2700 : dphase=0.0;
614 :
615 326700 : for(uInt i = 0; i < numRows; ++i){
616 1296000 : for (Int idim=0;idim<3;idim++){
617 972000 : uvw(idim,i)=vb.uvw()(i)(idim);
618 : }
619 : }
620 :
621 : // the following call requires the member variables
622 : // lastFieldId_p
623 : // lastMSId_p
624 : // tangentSpecified_p
625 : // MeasFrame mFrame_p == output ref frame for the UVW coordinates
626 : // MDirection mImage_p == output phase center
627 : // (input phase center is taken from the visbuffer, i.e. from the FIELD table)
628 :
629 2700 : FTMachine::rotateUVW(uvw, dphase, vb);
630 :
631 : // Immediately returns if not needed.
632 2700 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
633 :
634 : // update vis buffer cache
635 9540 : for(uInt datacol = 0; datacol < nDataCols_p; datacol++){
636 6840 : if(dataCols_p[datacol] == MS::DATA){
637 2700 : vb.visCube(); // return value not needed
638 : }
639 4140 : else if(dataCols_p[datacol] == MS::CORRECTED_DATA){
640 2160 : vb.correctedVisCube();
641 : }
642 1980 : else if(dataCols_p[datacol] == MS::MODEL_DATA){
643 1980 : vb.modelVisCube();
644 : }
645 : }
646 :
647 : // apply phase center shift to vis buffer
648 2700 : vb.phaseCenterShift(dphase);
649 :
650 : // write changed UVWs
651 2700 : auto origRows = vb.rowIds();
652 326700 : for(uInt row = 0; row < numRows; row++){
653 324000 : UVWcol.put(origRows(row), uvw.column(row));
654 : }
655 :
656 : // write changed visibilities
657 9540 : for(uInt datacol = 0; datacol < nDataCols_p; datacol++){
658 6840 : if(dataCols_p[datacol] == MS::DATA){
659 2700 : vi.setVis(vb.visCube(),VisibilityIterator::Observed);
660 : }
661 4140 : else if(dataCols_p[datacol] == MS::CORRECTED_DATA){
662 2160 : vi.setVis(vb.correctedVisCube(),VisibilityIterator::Corrected);
663 : }
664 1980 : else if(dataCols_p[datacol] == MS::MODEL_DATA){
665 1980 : vi.setVis(vb.modelVisCube(),VisibilityIterator::Model);
666 : }
667 : }
668 :
669 :
670 2700 : }
671 :
672 : }
673 :
674 15 : }
675 :
676 :
677 5400 : void FixVis::ok() {
678 : // AlwaysAssert(image, AipsError);
679 5400 : }
680 :
681 0 : void FixVis::init()
682 : {
683 0 : logSink() << LogOrigin("FixVis", "init") << LogIO::NORMAL;
684 :
685 : //ok();
686 :
687 : //npol = image->shape()(2);
688 : //nchan = image->shape()(3);
689 0 : }
690 :
691 : // Initialize FTMachine.
692 0 : void FixVis::initializeToVis(ImageInterface<Complex>& iimage,
693 : const VisBuffer&)
694 : {
695 0 : image = &iimage;
696 :
697 0 : ok();
698 0 : init();
699 :
700 : // Initialize the maps for polarization and channel. These maps
701 : // translate visibility indices into image indices
702 : //RR initMaps(vb);
703 0 : }
704 :
705 0 : void FixVis::finalizeToVis()
706 : {
707 0 : }
708 :
709 : // Initialize the FFT to the Sky. Here we have to setup and initialize the
710 : // grid.
711 0 : void FixVis::initializeToSky(ImageInterface<Complex>& iimage,
712 : Matrix<Float>&, //weight,
713 : const VisBuffer&)
714 : {
715 : // image always points to the image
716 0 : image = &iimage;
717 :
718 0 : init();
719 :
720 : // Initialize the maps for polarization and channel. These maps
721 : // translate visibility indices into image indices
722 : //RR initMaps(vb);
723 0 : }
724 :
725 : } // End of casa namespace.
726 :
|
