Line data Source code
1 : //# VLAFiller.cc:
2 : //# Copyright (C) 1999,2000,2001,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: casa-feedback@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: VLAFiller.cc,v 19.20.4.1 2005/12/26 21:09:43 wyoung Exp $
27 :
28 : #include <nrao/VLA/VLAFiller.h>
29 : #include <nrao/VLA/VLAADA.h>
30 : #include <nrao/VLA/VLACDA.h>
31 : #include <nrao/VLA/VLARCA.h>
32 : #include <nrao/VLA/VLASDA.h>
33 : #include <casacore/casa/Arrays/Array.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 : #include <casacore/casa/Arrays/ArrayMath.h>
36 : #include <casacore/casa/Arrays/Cube.h>
37 : #include <casacore/casa/Arrays/Matrix.h>
38 : #include <casacore/casa/Arrays/Slice.h>
39 : #include <casacore/casa/Arrays/Vector.h>
40 : #include <casacore/casa/Containers/RecordFieldId.h>
41 : #include <casacore/casa/Exceptions/Error.h>
42 : #include <casacore/casa/Logging/LogFilter.h>
43 : #include <casacore/casa/Logging/LogMessage.h>
44 : #include <casacore/casa/Logging/LogOrigin.h>
45 : #include <casacore/casa/Logging/LogSink.h>
46 : #include <casacore/casa/BasicSL/Complex.h>
47 : #include <casacore/ms/MeasurementSets/MSAntenna.h>
48 : #include <casacore/ms/MeasurementSets/MSAntennaColumns.h>
49 : #include <casacore/ms/MeasurementSets/MSDataDescColumns.h>
50 : #include <casacore/ms/MeasurementSets/MSDataDescription.h>
51 : #include <casacore/ms/MeasurementSets/MSFeed.h>
52 : #include <casacore/ms/MeasurementSets/MSFeedColumns.h>
53 : #include <casacore/ms/MeasurementSets/MSField.h>
54 : #include <casacore/ms/MeasurementSets/MSFieldColumns.h>
55 : #include <casacore/ms/MSSel/MSFieldIndex.h>
56 : #include <casacore/ms/MeasurementSets/MSObsColumns.h>
57 : #include <casacore/ms/MeasurementSets/MSObservation.h>
58 : #include <casacore/ms/MeasurementSets/MSPolColumns.h>
59 : #include <casacore/ms/MeasurementSets/MSPolarization.h>
60 : #include <casacore/ms/MeasurementSets/MSSpWindowColumns.h>
61 : #include <casacore/ms/MeasurementSets/MSSpectralWindow.h>
62 : #include <casacore/ms/MeasurementSets/MSTileLayout.h>
63 : #include <casacore/measures/Measures/MCBaseline.h>
64 : #include <casacore/measures/Measures/MDirection.h>
65 : #include <casacore/measures/Measures/MDoppler.h>
66 : #include <casacore/measures/Measures/MEpoch.h>
67 : #include <casacore/measures/Measures/MPosition.h>
68 : #include <casacore/measures/Measures/MeasRef.h>
69 : #include <casacore/measures/Measures/MeasTable.h>
70 : #include <casacore/casa/OS/File.h>
71 : #include <casacore/casa/OS/Path.h>
72 : #include <casacore/casa/Quanta/MVAngle.h>
73 : #include <casacore/casa/Quanta/MVBaseline.h>
74 : #include <casacore/casa/Quanta/MVDirection.h>
75 : #include <casacore/casa/Quanta/MVDoppler.h>
76 : #include <casacore/casa/Quanta/MVEpoch.h>
77 : #include <casacore/casa/Quanta/MVFrequency.h>
78 : #include <casacore/casa/Quanta/MVPosition.h>
79 : #include <casacore/casa/Quanta/MVTime.h>
80 : #include <casacore/casa/Quanta/MVuvw.h>
81 : #include <casacore/casa/Quanta/Quantum.h>
82 : #include <casacore/casa/Quanta/Unit.h>
83 : #include <casacore/measures/TableMeasures/ArrayMeasColumn.h>
84 : #include <casacore/measures/TableMeasures/ScalarMeasColumn.h>
85 : #include <casacore/tables/Tables/ArrayColumn.h>
86 : #include <casacore/tables/DataMan/IncrementalStMan.h>
87 : #include <casacore/tables/DataMan/TiledColumnStMan.h>
88 : #include <casacore/tables/Tables/RefRows.h>
89 : #include <casacore/tables/Tables/ScaColDesc.h>
90 : #include <casacore/tables/Tables/ScalarColumn.h>
91 : #include <casacore/tables/Tables/SetupNewTab.h>
92 : #include <casacore/tables/DataMan/StandardStMan.h>
93 : #include <casacore/tables/Tables/Table.h>
94 : #include <casacore/tables/Tables/TableDesc.h>
95 : #include <casacore/tables/Tables/TableInfo.h>
96 : #include <casacore/tables/Tables/TableRecord.h>
97 : #include <casacore/tables/Tables/TableUtil.h>
98 : #include <casacore/tables/DataMan/TiledDataStMan.h>
99 : #include <casacore/tables/DataMan/TiledShapeStMan.h>
100 : #include <casacore/casa/Utilities/Assert.h>
101 : #include <casacore/casa/Utilities/DataType.h>
102 : #include <casacore/casa/BasicSL/String.h>
103 : #include <sstream>
104 : #include <iomanip>
105 :
106 : #include <casacore/tables/TaQL/ExprNode.h>
107 :
108 : const uInt maxCDA(VLAEnum::CDA3+1);
109 : const uInt maxIF(VLAEnum::IFD+1);
110 : const uInt nCat = 6; // Number of Flag categories
111 : const uInt maxSubarrays = 5; // The maximum number of subarrays;
112 : const String sigmaCol = "sigmaHyperColumn";
113 : const String dataCol = "dataHyperColumn";
114 : const String flagCol = "flagHyperColumn";
115 : //New addition
116 : const String modDataCol = "modDataHyperColumn";
117 : const String corrDataCol = "corrDataHyperColumn";
118 : const String chanFlagCol = "flagChanHyperColumn";
119 : //====
120 : const RecordFieldId sigmaTileId("SIGMA_HYPERCUBE_ID");
121 : const RecordFieldId dataTileId("DATA_HYPERCUBE_ID");
122 : const RecordFieldId flagTileId("FLAG_CATEGORY_HYPERCUBE_ID");
123 : //=====new addition
124 : const RecordFieldId modDataTileId("MODEL_DATA_HYPERCUBE_ID");
125 : const RecordFieldId corrDataTileId("CORR_DATA_HYPERCUBE_ID");
126 : const RecordFieldId chanFlagTileId("FLAG_HYPERCUBE_ID");
127 : //=====
128 : const Quantum<Double> dirTol(10.0, "mas"); // Tolerance on matching fields
129 : const Quantum<Double> posTol(1, "m"); // Tolerance on matching antennas
130 : const Double ns2m = QC::c( ).getValue("m/ns");
131 :
132 : // The following struct is just ways of bundling up a bunch of arguments to
133 : // functions to avoid having too many in the function interface.
134 : struct IterationStatus {
135 : uInt nVLARecords;
136 : uInt nRows;
137 : uInt nAnt;
138 : Block<Block<uInt> > lastAnt;
139 : uInt nSpw;
140 : Block<Block<Int> > lastSpw;
141 : uInt nPol;
142 : Block<Block<Int> > lastPol;
143 : uInt nFld;
144 : Block<Int> lastFld;
145 : String lastProject;
146 : String lastObsMode;
147 : };
148 :
149 16 : VLAFiller::VLAFiller(MeasurementSet& output, VLALogicalRecord& input, Double freqTolerance, Bool autocorr, const String& antnamescheme, const Bool& applyTsys):
150 : MSColumns(output),
151 16 : itsRecord(input),
152 16 : itsInputFilters(),
153 16 : itsMS(output),
154 16 : itsFrame(),
155 : //theirFrames(0, MeasFrame()),
156 16 : itsMSDirType(MDirection::castType(field().referenceDirMeasCol().getMeasRef().getType())),
157 16 : itsDirType(MDirection::N_Types),
158 16 : itsDirCtr(),
159 16 : itsAzElCtr(),
160 16 : itsUvwCtr(),
161 16 : itsFreqCtr(),
162 16 : itsBlCtr(),
163 16 : itsFldId(maxSubarrays, -1),
164 16 : itsAntId(),
165 16 : itsSpId(maxCDA, -1),
166 16 : itsPolId(0),
167 16 : itsDataId(maxSubarrays),
168 16 : itsNewScan(true),
169 16 : itsScan(maxSubarrays, 0),
170 16 : itsProject(),
171 16 : itsLog(),
172 16 : itsDataAcc(),
173 16 : itsTileId(),
174 16 : itsSigmaAcc(),
175 16 : itsFlagAcc(),
176 16 : itsModDataAcc(),
177 16 : itsCorrDataAcc(),
178 16 : itsChanFlagAcc(),
179 16 : itsDataShapes(0),
180 16 : itsFreqTolerance(freqTolerance),
181 16 : itsApplyTsys(applyTsys),
182 16 : itsEVLAisOn(false),
183 16 : itsInitEpoch(false),
184 16 : itsRevBeenWarned(false),
185 16 : itsNoPolInfoWarned(false),
186 48 : itsZeroEpochWarned(false)
187 : {
188 16 : String antscheme=antnamescheme;
189 16 : antscheme.downcase();
190 16 : itsNewAntName=false;
191 16 : if(antscheme.contains("new"))
192 10 : itsNewAntName=true;
193 16 : itsKeepAutoCorr=autocorr;
194 16 : checkStop = false;
195 16 : fillStarted = false;
196 16 : AlwaysAssert(itsMS.tableInfo().subType() == "VLA", AipsError);
197 : /*
198 : itsDataAcc = TiledDataStManAccessor(itsMS, dataCol);
199 : itsModDataAcc=TiledDataStManAccessor(itsMS,modDataCol);
200 : itsCorrDataAcc=TiledDataStManAccessor(itsMS,corrDataCol);
201 : itsChanFlagAcc=TiledDataStManAccessor(itsMS,chanFlagCol);
202 : itsSigmaAcc = TiledDataStManAccessor(itsMS, sigmaCol);
203 : itsFlagAcc = TiledDataStManAccessor(itsMS, flagCol);
204 : const uInt nShapes = itsDataAcc.nhypercubes();
205 : itsDataShapes.resize(nShapes);
206 : for (uInt s = 0; s < nShapes; s++) {
207 : itsDataShapes[s] = itsDataAcc.getHypercubeShape(s).getFirst(2);
208 : }
209 : */
210 : // This ms is starting ...make sure its not J2000 when data is B1950
211 : // avoid unecessary conversions
212 16 : if(nrow()==0){
213 :
214 9 : itsInitEpoch=false;
215 :
216 : }
217 : else{
218 7 : itsInitEpoch=true;
219 : }
220 :
221 : // Deduce what the last scan was from the LAST_SCAN keyword in the SCAN
222 : // column. The alternative would be to read all the data in this column!
223 : // Note that the scan number is now per subarray.
224 :
225 16 : if (nrow() != 0) {
226 7 : const RecordFieldId key("LAST_SCAN");
227 7 : const TableRecord& keywords = scanNumber().keywordSet();
228 7 : DebugAssert(keywords.isDefined(key.fieldName()), AipsError);
229 7 : if (keywords.dataType(key) == TpInt) { // Old ms (v16 or earlier)
230 0 : itsScan[0] = keywords.asInt(key);
231 7 : } else if (keywords.dataType(key) == TpArrayInt) {
232 7 : DebugAssert(keywords.shape(key) == IPosition(1, maxSubarrays),AipsError);
233 7 : const Vector<Int> lastscans(keywords.asArrayInt(key));
234 7 : itsScan = makeBlock(lastscans);
235 7 : }
236 7 : }
237 :
238 : // Check if the frame has been initialised. If not then initialise it with a
239 : // position, epoch & direction measure. The direction and time values (but
240 : // not types) will be updated in the fillOne function. The position should
241 : // never be changed.
242 16 : itsFrame.set(MDirection(MVDirection(), itsMSDirType));
243 16 : itsFrame.set(MEpoch(MVEpoch(), timeMeas().getMeasRef()));
244 16 : MPosition vlaCentre;
245 16 : AlwaysAssert(MeasTable::Observatory(vlaCentre, "VLA"), AipsError);
246 16 : itsFrame.set(vlaCentre);
247 :
248 : // For the AzEl converter, we only need AzEl
249 16 : itsAzElCtr.setOut(MeasRef<MDirection>(MDirection::AZEL, itsFrame));
250 :
251 : // For the direction and uvw converter the output type is fixed.
252 16 : itsDirCtr.setOut(MeasRef<MDirection>(itsMSDirType, itsFrame));
253 16 : itsUvwCtr.setOut(MeasRef<Muvw>(Muvw::fromDirType(itsMSDirType), itsFrame));
254 : // For the frequency converter the input type is fixed.
255 16 : itsFreqCtr.setModel(MFrequency(MVFrequency(),
256 32 : MFrequency::Ref(MFrequency::TOPO, itsFrame)));
257 : // For the baseline converter both the input and output frames are known.
258 16 : itsBlCtr.set(MBaseline(MVBaseline(),
259 32 : MBaseline::Ref(MBaseline::HADEC, itsFrame)),
260 : MBaseline::ITRF);
261 16 : }
262 :
263 16 : VLAFiller::~VLAFiller()
264 : {
265 16 : }
266 :
267 16 : void VLAFiller::setFilter(const VLAFilterSet& filter){
268 16 : itsInputFilters = filter;
269 16 : }
270 :
271 0 : void VLAFiller::setStopParams(String &pCode, String &sTime){
272 0 : if(!pCode.empty() || !sTime.empty()){
273 0 : projectCode = upcase(pCode);
274 0 : if(!sTime.empty()){
275 0 : Quantum<Double> t;
276 0 : MVTime::read(t, sTime);
277 0 : stopTime = MVEpoch(t);
278 0 : }
279 0 : fillStarted = false;
280 0 : checkStop = true;
281 : }
282 0 : return;
283 : }
284 : // Well this needs some more work as right now I can only figure
285 : // out how to have it stop filling after a time. Looks like the
286 : // project code can end up being something like system, especially
287 : // if there are subarrays involved.
288 :
289 20569 : Bool VLAFiller::stopFilling(VLALogicalRecord &record)
290 : {
291 20569 : Bool rstat(false);
292 20569 : if(checkStop){
293 : /*
294 : if(fillStarted){
295 : //std::cerr << projectCode << std::endl;
296 : //std::cerr << upcase(record.SDA().obsId()) << std::endl;
297 : String itsCode(upcase(record.SDA().obsId()));
298 : if(!projectCode.empty() && !projectCode.matches(itsCode)){
299 : rstat = true;
300 : }
301 : }
302 : */
303 0 : const Double recordTime = record.RCA().obsDay() +
304 0 : Quantum<Double>(record.SDA().obsTime(), "s").getValue("d");
305 : //std::cerr << stopTime.get() << std::endl;
306 : //std::cerr << recordTime << std::endl;
307 0 : if(recordTime > stopTime.get())
308 0 : rstat = true;
309 : else
310 0 : rstat = false;
311 : }
312 20569 : return rstat;
313 : }
314 :
315 16 : void VLAFiller::fill(Int verbose){
316 16 : itsLog.origin(LogOrigin("VLAFiller", "fill"));
317 16 : IterationStatus counts;
318 16 : counts.nVLARecords = 0;
319 16 : counts.nRows = nrow();
320 16 : counts.nAnt = antenna().nrow();
321 16 : counts.lastAnt.resize(maxSubarrays);
322 16 : counts.nSpw = spectralWindow().nrow();
323 16 : counts.lastSpw = Block<Block<Int> >(maxSubarrays, Block<Int>(2, -1));
324 16 : counts.nPol = polarization().nrow();
325 16 : counts.lastPol = Block<Block<Int> >(maxSubarrays, Block<Int>(2, -1));;
326 16 : counts.nFld = field().nrow();
327 16 : counts.lastFld.resize(maxSubarrays);
328 16 : counts.lastFld = -1;
329 16 : counts.lastProject = "";
330 16 : counts.lastObsMode = "";
331 16 : const uInt initialRow = nrow();
332 16 : itsRevBeenWarned = false;
333 : #if defined(AIPS_DEBUG)
334 16 : const LogFilter saved(LogMessage::DEBUGGING);
335 16 : if (verbose > 0) LogSink::globalSink().filter(saved);
336 : #endif
337 16 : AipsError error;
338 : try {
339 :
340 16 : String notstr("NOT ");
341 16 : if (itsApplyTsys)
342 13 : notstr="";
343 :
344 16 : itsLog << LogIO::NORMAL
345 32 : << "Data and weights will "+notstr+"be scaled by Nominal Sensitivity."
346 32 : << LogIO::POST;
347 :
348 :
349 20585 : while (fillOne()) {
350 20569 : counts.nVLARecords++;
351 20569 : if (nrow() != counts.nRows) {
352 4570 : if (verbose > 0 &&
353 0 : counts.nVLARecords%verbose == 0) {
354 0 : logCurrentRecord(counts);
355 : }
356 4570 : logChanges(counts);
357 : }
358 20569 : counts.nRows = nrow();
359 : }
360 16 : } catch (AipsError x) {
361 0 : itsLog << LogIO::SEVERE
362 : << "An error occurred. The error message is:" << endl
363 : << "'" << x.getMesg() << "'" << endl
364 : << "Perhaps you ran out of disk space or "
365 : << "are using a flaky tape (drive)." << endl
366 : << "Trying to write a valid measurement set with" << endl
367 0 : << "as much data as possible." << LogIO::POST;
368 0 : error = x;
369 0 : } catch (...){
370 0 : itsLog << LogIO::SEVERE << "Something really bad happened!" << LogIO::POST;
371 0 : }
372 : // Now fixup the observation subtable (only if new data has been added).
373 16 : if ((nrow() - initialRow) > 0) {
374 11 : MSObservation& msObs = itsMS.observation();
375 11 : const uInt newRow = msObs.nrow();
376 11 : msObs.addRow();
377 11 : MSObservationColumns& obsCols = observation();
378 11 : obsCols.telescopeName().put(newRow, "VLA");
379 11 : const String unavailable = "unavailable";
380 11 : const String unknown = "unknown";
381 11 : const Vector<String> unavailableVec(1, unavailable);
382 11 : obsCols.observer().put(newRow, unavailable);
383 11 : obsCols.log().put(newRow, unavailableVec);
384 11 : obsCols.scheduleType().put(newRow, unknown);
385 11 : obsCols.schedule().put(newRow, unavailableVec);
386 11 : obsCols.project().put(newRow, itsProject);
387 11 : obsCols.flagRow().put(newRow, false);
388 11 : Vector<MEpoch> obsTimeRange(2);
389 11 : obsTimeRange(0) = timeMeas()(initialRow);
390 11 : const uInt finalRow = nrow() - 1;
391 11 : obsTimeRange(1) = timeMeas()(finalRow);
392 11 : obsCols.timeRangeMeas().put(newRow, obsTimeRange);
393 11 : MVEpoch releaseDate = obsTimeRange(1).getValue();
394 11 : releaseDate += MVEpoch(Quantum<Double>(1.5, "a"));
395 11 : obsTimeRange(1).set(releaseDate);
396 11 : obsCols.releaseDateMeas().put(newRow, obsTimeRange(1));
397 11 : ostringstream oos;
398 :
399 : // Tidy up FIELD name duplicates
400 11 : fixFieldDuplicates(itsMS.field());
401 :
402 11 : } else {
403 5 : if (nrow() == 0) {
404 0 : itsLog << "No data in the measurement set\n";
405 : } else {
406 5 : itsLog << "No data appended to the measurement set\n";
407 : }
408 5 : if (counts.nVLARecords == 0) {
409 0 : itsLog << LogIO::SEVERE << "Your input may not be in VLA archive format"
410 0 : << LogIO::POST;
411 : } else {
412 5 : itsLog << "There may be a problem with your data selection criteria"
413 5 : << LogIO::WARN << LogIO::POST;
414 : }
415 : }
416 :
417 :
418 16 : if (verbose >= 0) {
419 16 : summarise();
420 : }
421 :
422 32 : scanNumber().rwKeywordSet().define(RecordFieldId("LAST_SCAN"),
423 32 : Vector<Int>(itsScan.begin(),itsScan.end()));
424 :
425 : #if defined(AIPS_DEBUG)
426 16 : LogSink::globalSink().filter(saved);
427 : #endif
428 16 : if (error.getMesg().length() > 0) throw(error);
429 16 : }
430 :
431 20585 : Bool VLAFiller::fillOne() {
432 20585 : if (!itsRecord.read()) return false;
433 20569 : if (stopFilling(itsRecord)) return false;
434 20569 : if (!itsInputFilters.passThru(itsRecord)){
435 : // OK here we mark a new scan if we skip anything.
436 15999 : itsNewScan=true;
437 15999 : return true;
438 : }
439 :
440 4570 : fillStarted = true;
441 4570 : const VLARCA& rca = itsRecord.RCA();
442 4570 : const VLASDA& sda = itsRecord.SDA();
443 : //For new ms and first go...make sure to init this to B1950 if data is so
444 : // as default blank ms is in J2000 direction
445 4570 : if(!itsInitEpoch){
446 9 : itsMSDirType=validEpoch(sda.epoch());
447 9 : itsFrame.set(MDirection(MVDirection(), itsMSDirType));
448 9 : setDirectionRef(itsMSDirType);
449 9 : setUVWRef(Muvw::fromDirType(itsMSDirType));
450 : // For the direction and uvw converter the output type is fixed.
451 9 : itsDirCtr.setOut(MeasRef<MDirection>(itsMSDirType, itsFrame));
452 9 : itsUvwCtr.setOut(MeasRef<Muvw>(Muvw::fromDirType(itsMSDirType), itsFrame));
453 9 : itsMS.initRefs();
454 :
455 9 : itsInitEpoch=true;
456 9 : itsMS.flush();
457 : }
458 : // Check if the revision number is supported.
459 4570 : if (rca.revision() < 23 && !itsRevBeenWarned) {
460 0 : itsRevBeenWarned = true;
461 0 : itsLog << LogIO::WARN
462 : << "This function has not been tested on VLA archive data "
463 : << "with revisions less " << endl
464 : << "than 23 & the data in this record has a revision level of "
465 : << rca.revision() << endl
466 : << "It is very likely that the correlation data will be scaled "
467 : << "incorrectly"
468 0 : << LogIO::POST;
469 : }
470 4570 : const uInt subArray = sda.subArray() - 1;
471 : { // Keep track of which projects have been copied.
472 4570 : const String thisProject = sda.obsId();
473 4570 : if (!itsProject.contains(thisProject)) {
474 11 : if (itsProject.length() != 0) itsProject += String(" ");
475 11 : itsProject += thisProject;
476 : }
477 4570 : }
478 : { // set the observing time. Do this now as it may be needed to convert the
479 : // field directions from the observed direction type to the direction type
480 : // used in the MS.
481 9140 : const MVEpoch obsTime(rca.obsDay(),
482 9140 : Quantum<Double>(sda.obsTime(), "s").getValue("d"));
483 : // cerr << MVTime(obsTime.getTime()).string(MVTime::YMD) << endl;
484 4570 : itsFrame.resetEpoch(obsTime);
485 : //NEED to use the exact date the EVLA antenna got in
486 : // If after 2005 EVLA can be in
487 4570 : if(obsTime.getDay() > 53371.0)
488 0 : itsEVLAisOn=true;
489 : else
490 4570 : itsEVLAisOn=false;
491 4570 : }
492 :
493 : { // Workout the field ID.
494 4570 : const MVDirection fieldDirVal(sda.sourceDir());
495 4570 : const MDirection::Types fieldDirRef(validEpoch(sda.epoch()));
496 : // Need to convert the direction to the same type as the MS. Otherwise the
497 : // match will fail.
498 :
499 4570 : MDirection fieldDir(fieldDirVal, fieldDirRef);
500 4570 : if (fieldDirRef != itsMSDirType) { // need to do a conversion
501 0 : if (fieldDirRef == itsDirType) { // no need to setup the converter
502 0 : fieldDir = itsDirCtr(fieldDirVal);
503 : } else {
504 0 : itsDirCtr.setModel(fieldDir);
505 0 : itsDirType = fieldDirRef;
506 0 : fieldDir = itsDirCtr();
507 : // at the same time the UVW converter should be initialised
508 0 : itsUvwCtr.setModel(Muvw(MVuvw(), Muvw::fromDirType(fieldDirRef)));
509 : }
510 : }
511 :
512 : // Determine if field already exists in FIELD subtable
513 4570 : Bool fldMatch=false;
514 :
515 : // First match on name (maybe multiple name matches with diff directions?):
516 : // (this is a clumsy way--and inefficient for large FIELD tables--to
517 : // include name matching, should have andName option in matchDirection)
518 4570 : MSField& msFld=itsMS.field();
519 :
520 : //Damn Damnation...as there is no MSColumns::attach ...need to use
521 : // a refreshed mscolumn ...especially if the epoch of the direction
522 : // is resetted above..for now create a redundant msc....
523 : //to do a matchdirection...need to enhance mscolumns to have an attach
524 4570 : MSColumns msc(itsMS);
525 4570 : MSFieldIndex MSFldIdx(msFld);
526 4570 : Vector<Int> fldNameMatches = MSFldIdx.matchFieldName(sda.sourceName());
527 : Int nfNM;
528 4570 : fldNameMatches.shape(nfNM);
529 : // found at least one name match, verify/discern from direction matching
530 4570 : Int ifNM=0;
531 4570 : if (nfNM > 0) {
532 9080 : while (ifNM < nfNM && !fldMatch) {
533 4540 : fldMatch = (fldNameMatches(ifNM)==msc.field().matchDirection(fieldDir, fieldDir, fieldDir, dirTol,
534 4540 : fldNameMatches(ifNM)));
535 4540 : if (!fldMatch) ifNM++;
536 : }
537 : }
538 :
539 : Int thisfldId;
540 4570 : if (fldMatch) {
541 : // found match:
542 4540 : thisfldId=fldNameMatches(ifNM);
543 : } else {
544 : // found no match, adding a row:
545 30 : addSource(fieldDir);
546 30 : thisfldId=addField(fieldDir);
547 : }
548 :
549 4570 : if (thisfldId != itsFldId[subArray]) {
550 112 : itsFrame.resetDirection(fieldDir.getValue());
551 112 : itsNewScan = true;
552 112 : itsFldId[subArray] = thisfldId;
553 : }
554 4570 : }
555 :
556 4570 : const uInt nAnt = rca.nAntennas();
557 : // Cache the uvw coordinates of each antenna. For holography mode,
558 : // these are the az, el offsets.
559 4570 : Block<Double> antUvw(3*nAnt);
560 : {
561 4570 : uInt elem = 0;
562 4570 : Vector<Double> convertedUvw(3);
563 : Double u, v, w;
564 4570 : const Bool doConversion = (itsMSDirType == validEpoch(sda.epoch())) ? false : true;
565 126638 : for (uInt a = 0; a < nAnt; a++) {
566 122068 : const VLAADA& ada = itsRecord.ADA(a);
567 122068 : u = ada.u();
568 122068 : v = ada.v();
569 122068 : w = ada.w();
570 122068 : if (doConversion) {
571 0 : convertedUvw = itsUvwCtr(MVuvw(u, v, w)).getValue().getValue();
572 0 : u = convertedUvw(0);
573 0 : v = convertedUvw(1);
574 0 : w = convertedUvw(2);
575 : }
576 122068 : antUvw[elem] = u; elem++;
577 122068 : antUvw[elem] = v; elem++;
578 122068 : antUvw[elem] = w; elem++;
579 : }
580 4570 : }
581 :
582 4570 : Block<Bool> shadowed(nAnt, false);
583 4570 : Bool badData = false;
584 : { // find out if any antennae are shadowed
585 4570 : uInt a1Idx = 0;
586 126638 : for (uInt a1 = 0; a1 < nAnt; a1++) {
587 122068 : uInt a2Idx = (a1+1)*3;
588 1691766 : for (uInt a2 = a1+1; a2 < nAnt; a2++) {
589 1569698 : Double u = antUvw[a1Idx] - antUvw[a2Idx]; a2Idx++;
590 1569698 : Double v = antUvw[a1Idx+1] - antUvw[a2Idx]; a2Idx++;
591 1569698 : if (u*u + v*v < 625) {
592 278 : badData = true;
593 278 : Double w = antUvw[a1Idx+2] - antUvw[a2Idx];
594 278 : if (w > 0) {
595 0 : shadowed[a2] = true;
596 : } else {
597 278 : shadowed[a1] = true;
598 : }
599 : }
600 1569698 : a2Idx++;
601 : }
602 122068 : a1Idx += 3;
603 : }
604 : }
605 :
606 : // Workout the antenna ID
607 4570 : if (itsAntId.nelements() != nAnt) {
608 11 : itsAntId.resize(nAnt);
609 11 : itsAntId = -1;
610 : }
611 : {
612 4570 : DebugAssert(itsFrame.position() != 0, AipsError);
613 4570 : DebugAssert(MPosition::castType
614 : (itsFrame.position()->getRefPtr()->getType()) ==
615 : MPosition::ITRF, AipsError);
616 : const MVPosition vlaCentrePos =
617 4570 : dynamic_cast<const MPosition*>(itsFrame.position())->getValue();
618 :
619 :
620 4570 : Vector<Int> antOrder(29);
621 4570 : antOrder=-1;
622 4570 : Vector<MPosition> thisPos(nAnt);
623 126638 : for (uInt a = 0; a < nAnt; a++) {
624 122068 : const VLAADA& ada = itsRecord.ADA(a);
625 : // Need to do the conversion from bx, by, bz (which is the HADEC frame)
626 : // to the ITRF frame prior to adding the reference position.
627 : // However, bx,by,bz differ from HADEC by handedness, thus
628 : // negate the y-component so ant positions come out right-side-out:
629 : // (perhaps HADEC is not the right thing to use)
630 122068 : const MVBaseline thisBl(ada.bx(), -ada.by(), ada.bz());
631 122068 : MVPosition thisAnt = itsBlCtr(thisBl).getValue();
632 122068 : thisAnt += vlaCentrePos;
633 : // const MPosition thisPos(thisAnt, MPosition::ITRF);
634 122068 : thisPos(a) = MPosition(thisAnt, MPosition::ITRF);
635 122068 : String leAntName;
636 :
637 122068 : antOrder(ada.antId()-1)=a;
638 :
639 122068 : if(!itsEVLAisOn){
640 : // ignore the frontend temperature naming
641 122068 : leAntName=ada.antName(false);
642 122068 : if(itsNewAntName){
643 120124 : leAntName=String("VA")+leAntName;
644 : }
645 : }
646 : else{
647 0 : leAntName=ada.antName(itsNewAntName);
648 : }
649 244136 : itsAntId[a] = antenna().matchAntenna(leAntName, thisPos(a), posTol,
650 122068 : itsAntId[a]);
651 :
652 : // if (itsAntId[a]<0)
653 : // cout << a << " " << ada.antId() << " " << leAntName << endl;
654 :
655 122068 : }
656 :
657 : /*
658 : cout << nAnt << " " << antOrder.nelements() << " " << ntrue(antOrder>-1) << " "
659 : << itsAntId.nelements() << " " << max(antOrder) << " "
660 : << min(Vector<Int>(itsAntId,nAnt))
661 : << endl;
662 : cout << "antOrder = " << antOrder << endl;
663 : */
664 : // If there are antennas to add
665 4570 : if (min(Vector<Int>(itsAntId.begin(),nAnt,int(0)))<0) {
666 : // cout << "itsAntId 0 = " << Vector<Int>(itsAntId,nAnt) << endl;
667 :
668 300 : for (uInt ai=0; ai < antOrder.nelements(); ++ai) {
669 290 : if (antOrder(ai)>-1) {
670 268 : Int ao(antOrder(ai));
671 268 : if (itsAntId[ao] < 0) {
672 261 : itsAntId[ao] = addAntenna(thisPos(ao), ao);
673 261 : addFeed(itsAntId[ao]);
674 261 : itsNewScan = true;
675 : }
676 : }
677 : }
678 : }
679 : // cout << "itsAntId 1 = " << Vector<Int>(itsAntId,nAnt) << endl;
680 :
681 :
682 4570 : }
683 :
684 :
685 : // For holography data, add the pointing data which is to be
686 : // found in the u,v parts of the ADA
687 : // Is this Holography data? If so, the UVWs are actually the
688 : // pointing offsets - U = Az, V = El
689 4570 : Bool isHolo=(itsRecord.SDA().obsMode()=="H ");
690 4570 : if(isHolo) {
691 :
692 : // We store AzEl in the pointing table. We only need to do this
693 : // when the table is empty
694 0 : MSPointingColumns& pointingCol = pointing();
695 0 : if(pointingCol.nrow()==0) {
696 0 : pointingCol.directionMeasCol().setDescRefCode(MDirection::AZEL);
697 : }
698 :
699 0 : const MVDirection fieldDirVal(itsRecord.SDA().sourceDir());
700 0 : const MDirection::Types fieldDirRef(validEpoch(itsRecord.SDA().epoch()));
701 0 : MDirection fieldDir(fieldDirVal, fieldDirRef);
702 :
703 : // For the actual direction, we put (Az,El) = (0,0). For the
704 : // target, we put the actual Ra, Dec. The uv data (in ns!) is
705 : // really the pointing offset in radians.
706 0 : for (uInt a = 0; a < nAnt; a++) {
707 0 : if(itsAntId[a]>-1) {
708 0 : const VLAADA& ada = itsRecord.ADA(a);
709 0 : MDirection thisDir(MVDirection(0.0, 0.0), MDirection::AZEL);
710 0 : thisDir.shift(-ada.u()/ns2m, ada.v()/ns2m, true);
711 0 : addPointing(thisDir, fieldDir, itsAntId[a]);
712 0 : }
713 : }
714 0 : }
715 :
716 : // Now create a bunch of blocks that are necessary for reindexing the data
717 : // from different correlator blocks into the appropriate rows of the MS.
718 4570 : Block<Block<VLAEnum::CDA> > CDAId; // Maps the local spectral ID to CDA's
719 : // ie., CDAId[0:nSpID][0:nCDA] = whichCDA
720 4570 : Block<Block<uInt> > polId(maxCDA); // contains the polarisation index for
721 : // each CDA
722 : // ie., polId[0:4][0:nPols] = polIdx;
723 :
724 4570 : Block<Block<uInt> > polTypes(maxCDA);
725 4570 : Block<Block<uInt> > polNewOrder;
726 4570 : Vector<Bool> rotStokesOrder;
727 :
728 :
729 22850 : for (uInt c = 0; c < maxCDA; c++) {
730 18280 : const VLACDA& cda = itsRecord.CDA(c);
731 18280 : if (!cda.isValid()) {
732 10718 : itsSpId[c] = -1;
733 : } else {
734 7562 : const VLAEnum::CDA thisCDA = VLAEnum::CDA(c);
735 : // can not deal with npol = 0, may be arising in poorly understood old correlator modes, needs investigating
736 7562 : if (sda.npol(thisCDA) == 0) {
737 : // warn once and consider as if this is an invalid CDA
738 0 : if (!itsNoPolInfoWarned) {
739 0 : itsNoPolInfoWarned = true;
740 0 : itsLog << LogIO::SEVERE
741 : << "Unable to determine polarization information for some or all correlator modes." << endl
742 : << "That data can not be filled and the resulting visibility file may be empty."
743 0 : << LogIO::POST;
744 : }
745 0 : itsSpId[c] = -1;
746 : } else {
747 : // Firstly, determine the spectral characteristics of the
748 : // data in the current CDA
749 7562 : const uInt nChan = sda.nChannels(thisCDA);
750 7562 : const Unit hz("Hz");
751 7562 : const Double chanWidth = sda.channelWidth(thisCDA);
752 7562 : const Quantum<Double> bandwidth(nChan*chanWidth, hz);
753 : // The tolerance is set at 1/4 of a channel. It is not set smaller
754 : // because when Doppler tracking is used the total bandwidth, when
755 : // converted to the rest type, will be slightly different from the
756 : // topocentric value.
757 : // above is the original comments.
758 : // We reset the default tolerance for frequency to be the value of the
759 : // channel width and also give user the ability to pass in a tolerance
760 : // for frequency into vlafillerfromdisk(). For dataset G192 we need a
761 : // tolerance of 6 times of the channe width. For dataset NGC7538, one
762 : // has to give a tolerance as larger as 60 times its channel width ( 60000000Hz ).
763 :
764 7562 : if( itsFreqTolerance == 0.0 ) itsFreqTolerance = chanWidth;
765 7562 : const Quantum<Double> tolerance( itsFreqTolerance, hz);
766 : // Determine the reference frequency.
767 7562 : MFrequency refFreq;
768 : {
769 7562 : if (sda.dopplerTracking(thisCDA)) {
770 2824 : const MDoppler dop(Quantity(sda.radialVelocity(thisCDA), "m/s"),
771 2824 : sda.dopplerDefn(thisCDA));
772 : refFreq =
773 4236 : MFrequency::fromDoppler(dop,
774 2824 : MVFrequency(sda.restFrequency(thisCDA)),
775 1412 : sda.restFrame(thisCDA));
776 1412 : } else {
777 12300 : refFreq = MFrequency(MVFrequency(sda.obsFrequency(thisCDA)),
778 6150 : MFrequency::TOPO);
779 : }
780 : }
781 : // The local spectral Id is the value that is either zero or one and
782 : // depends on which IF the data in the CDA came from. Be aware that data
783 : // from IF B may have a local spectral Id value of either zero or one,
784 : // depending on whether IF A is also being used.
785 : uInt localSpId;
786 : // See if there is a matching row.
787 : {
788 7562 : const uInt nSpId = CDAId.nelements();
789 7562 : const uInt ifChain = sda.electronicPath(thisCDA);
790 : // set MeasFrame to MeasRef of MFrequency, which is need when converting MFrequency
791 : // a different frame.
792 : // refFreq.getRef().set( itsFrame );
793 : // no, ScalarMeasColumn<M>put() will not accept this! so instead, we do
794 : // Find the channel frequencies and pass the first one to matchSpw().
795 7562 : Vector<Double> chanFreqs(nChan);
796 7562 : indgen(chanFreqs, sda.edgeFrequency( thisCDA )+0.5*chanWidth, chanWidth);
797 7562 : const MFrequency::Types itsFreqType = MFrequency::castType(refFreq.getRef().getType());
798 7562 : if (itsFreqType != MFrequency::TOPO) {
799 : // have to convert the channel frequencies from topocentric to the specifed
800 : // frequency type.
801 1412 : MFrequency::Convert freqCnvtr;
802 1412 : freqCnvtr.setModel( MFrequency(MVFrequency(), MFrequency::Ref( MFrequency::TOPO, itsFrame )) );
803 1412 : freqCnvtr.setOut( itsFreqType );
804 1412 : Double freqInHzCnvtrd = freqCnvtr(chanFreqs(0)).getValue().getValue();
805 1412 : chanFreqs( 0 ) = freqInHzCnvtrd;
806 1412 : }
807 :
808 7562 : MFrequency chanFreq1 = MFrequency( MVFrequency( chanFreqs( 0 ) ), itsFreqType );
809 : // now call the matchSpw() method:
810 15124 : itsSpId[c] = spectralWindow().matchSpw(refFreq, chanFreq1, itsFrame, doppler(), source(), nChan, bandwidth,
811 7562 : ifChain, tolerance, itsSpId[c]);
812 :
813 : // for testing frequency handling
814 : /*
815 : cout.precision(12);
816 : cout << "Field = " << sda.sourceName()
817 : << " " << Int(thisCDA)
818 : << " lo=" << sda.edgeFrequency(thisCDA)
819 : << " (" << sda.obsFrequency(thisCDA)<<")"
820 : << " frame="<< sda.restFrame(thisCDA)
821 : << " v="<< sda.radialVelocity(thisCDA)
822 : << " rest="<< sda.restFrequency(thisCDA)
823 : << " freq1="<<chanFreqs(0)
824 : << " new="<<itsSpId[c]
825 : << endl;
826 : */
827 :
828 7562 : if (itsSpId[c] < 0) {
829 : // add an entry to Dopper subtable before addSpectralWindow! Also make sure addSouce is called before this!
830 31 : addDoppler( thisCDA );
831 31 : itsSpId[c] = addSpectralWindow(thisCDA, refFreq, nChan,
832 : bandwidth.getValue(hz), ifChain);
833 31 : localSpId = nSpId;
834 : } else {
835 7531 : localSpId = 0;
836 7531 : while (localSpId < nSpId &&
837 13489 : CDAId[localSpId].nelements() > 0 &&
838 2979 : itsSpId[CDAId[localSpId][0]] != itsSpId[c]) {
839 2979 : localSpId++;
840 : }
841 : }
842 7562 : if (localSpId == nSpId) {
843 7562 : CDAId.resize(nSpId + 1);
844 : }
845 7562 : }
846 : // Put this CDA into its spot the indexing blocks.
847 7562 : const uInt nCDA = CDAId[localSpId].nelements();
848 7562 : CDAId[localSpId].resize(nCDA + 1);
849 7562 : CDAId[localSpId][nCDA] = thisCDA;
850 7562 : uInt polIdx = 0;
851 7562 : if (nCDA != 0) { // Here is a tricky section for you.
852 : // The debugging statements should help
853 0 : const Block<uInt>& prevPolId = polId[CDAId[localSpId][nCDA-1]];
854 0 : polIdx = prevPolId[prevPolId.nelements()-1] + 1;
855 :
856 : #if defined(AIPS_DEBUG)
857 0 : itsLog << LogIO::DEBUGGING;
858 0 : itsLog << "CDA's containing this spectral ID: [";
859 0 : for (uInt c = 0; c < CDAId[localSpId].nelements(); c++) {
860 0 : itsLog << static_cast<Int>(CDAId[localSpId][c])
861 0 : << ((c+1 < CDAId[localSpId].nelements()) ? ", " : "]\n");
862 : }
863 0 : itsLog << "The previous CDA containing this spectral ID: "
864 0 : << static_cast<Int>(CDAId[localSpId][nCDA-1]) << endl;
865 0 : itsLog << "The polarization map of this CDA: [" ;
866 : {
867 0 : const uInt w = CDAId[localSpId][nCDA-1];
868 0 : for (uInt c = 0; c < polId[w].nelements(); c++) {
869 0 : itsLog << polId[w][c]
870 0 : << ((c+1 < polId[w].nelements()) ? ", " : "]\n");
871 : }
872 : }
873 0 : itsLog << "The last element of the polarization map: "
874 0 : << prevPolId.nelements()-1 << endl;
875 0 : itsLog << "The next polarization starts at index: "
876 0 : << polIdx << endl;
877 0 : itsLog << LogIO::POST << LogIO::NORMAL;
878 : #endif
879 : }
880 7562 : const uInt nPols = sda.npol(thisCDA);
881 7562 : polId[c].resize(nPols);
882 33076 : for (uInt p = 0; p < nPols; p++) {
883 25514 : polId[c][p] = p + polIdx;
884 : }
885 7562 : }
886 : }
887 : }
888 4570 : const uInt nSpId = CDAId.nelements();
889 4570 : if (nSpId == 0) {
890 : // This can occur if there is a single antenna subarray doing VLBI. This
891 : // antenna may not be connected to the VLA correlator and hence the
892 : // auto-correlation cannot be calculated. In this case all the CDA's are
893 : // invalid and there is no data to add to the main table of the MS.
894 0 : DebugAssert(nAnt == 1, AipsError);
895 0 : return true;
896 : }
897 :
898 : // Check if the transfer switch is only set on some antennas. If so warn
899 : // the user that the polarization may be misslabeled
900 : {
901 4570 : const Stokes::StokesTypes ant0Pol = itsRecord.ADA(0).ifPol(VLAEnum::IFA);
902 122068 : for (uInt a = 1; a < nAnt; a++) {
903 117498 : if (itsRecord.ADA(a).ifPol(VLAEnum::IFA) != ant0Pol) {
904 : // only warn if there's been no warning on this antenna yet - only ants with warnings are ever here
905 0 : if (itsTransferWarned.count(itsRecord.ADA(a).antName(itsNewAntName)) == 0) {
906 0 : itsLog << LogIO::WARN
907 : << "The IF transfer switch for antenna "
908 0 : << itsRecord.ADA(a).antName(itsNewAntName)
909 : << " is different from the setting for antenna "
910 0 : << itsRecord.ADA(0).antName(itsNewAntName) << "." << endl
911 : << "Correlations involving this antenna may have "
912 : << "incorrect polarization labelling."
913 0 : << LogIO::POST;
914 0 : itsTransferWarned[itsRecord.ADA(a).antName(itsNewAntName)] = true;
915 : }
916 : }
917 : }
918 : }
919 :
920 : // Now sort out the polarisation subtable
921 4570 : if (nSpId != itsPolId.nelements()) {
922 11 : itsPolId.resize(nSpId, true);
923 11 : itsPolId = -1;
924 : }
925 4570 : polTypes.resize(nSpId);
926 4570 : polNewOrder.resize(nSpId);
927 4570 : rotStokesOrder.resize(nSpId);
928 4570 : rotStokesOrder.set(false);
929 12132 : for (uInt s = 0; s < nSpId; s++) {
930 7562 : const Block<VLAEnum::CDA>& usedCDAs = CDAId[s];
931 7562 : const uInt nCda = usedCDAs.nelements();
932 7562 : uInt nPol = 0;
933 15124 : for (uInt i = 0; i < nCda; i++) {
934 7562 : nPol += polId[usedCDAs[i]].nelements();
935 : }
936 :
937 7562 : Vector<Stokes::StokesTypes> allPols(nPol);
938 7562 : uInt p = 0;
939 15124 : for (uInt i = 0; i < nCda; i++) {
940 7562 : Vector<Stokes::StokesTypes> pol(itsRecord.polarisations(usedCDAs[i]));
941 33076 : for (uInt j = 0; j < pol.nelements(); j++, p++) {
942 25514 : allPols(p) = pol(j);
943 : }
944 7562 : }
945 7562 : polTypes[s].resize(allPols.nelements());
946 7562 : polNewOrder[s].resize(allPols.nelements());
947 :
948 7562 : Bool standard=true;
949 33076 : for (uInt i=0; i < allPols.nelements(); ++i){
950 25514 : polTypes[s][i]=static_cast<uInt> (allPols[i]);
951 25514 : polNewOrder[s][i]=i;
952 25514 : if( (allPols[i] > Stokes::LL) || (allPols[i] < Stokes::RR)){
953 0 : standard=false;
954 : }
955 : }
956 : //Now if the 4-stokes are not in RR RL LR LL order....make sure it is
957 7562 : if((allPols.nelements() == 4) && standard ){
958 5984 : if(allPols[0] != Stokes::RR){
959 0 : rotStokesOrder[s]=true;
960 0 : polNewOrder[s][0]=polIndexer(allPols[0]);
961 : }
962 5984 : if(allPols[1] != Stokes::RL){
963 5984 : rotStokesOrder[s]=true;
964 5984 : polNewOrder[s][1]=polIndexer(allPols[1]);
965 : }
966 5984 : if(allPols[2] != Stokes::LR){
967 5984 : rotStokesOrder[s]=true;
968 5984 : polNewOrder[s][2]=polIndexer(allPols[2]);
969 : }
970 5984 : if(allPols[3] != Stokes::LL){
971 5984 : rotStokesOrder[s]=true;
972 5984 : polNewOrder[s][3]=polIndexer(allPols[3]);
973 : }
974 :
975 5984 : allPols[0]=Stokes::RR;
976 5984 : allPols[1]=Stokes::RL;
977 5984 : allPols[2]=Stokes::LR;
978 5984 : allPols[3]=Stokes::LL;
979 : }
980 7562 : itsPolId[s] = polarization().match(allPols, itsPolId[s]);
981 7562 : if (itsPolId[s] < 0) {
982 9 : itsPolId[s] = addPolarization(allPols);
983 : }
984 7562 : }
985 :
986 : // Keep these values handy, as they are needed in lots of places
987 4570 : Block<uInt> nChannels(nSpId);
988 4570 : Block<uInt> nProducts(nSpId);
989 12132 : for (uInt s = 0; s < nSpId; s++) {
990 7562 : nProducts[s] = polarization().numCorr()(itsPolId[s]);
991 7562 : nChannels[s] = spectralWindow().numChan()(itsSpId[CDAId[s][0]]);
992 : }
993 :
994 4570 : Block<Int>& thisDataId = itsDataId[subArray];
995 : // Now sort out the data description subtable
996 4570 : if (nSpId != thisDataId.nelements()) {
997 11 : thisDataId.resize(nSpId, true);
998 11 : thisDataId = -1;
999 : }
1000 12132 : for (uInt s = 0; s < nSpId; s++) {
1001 7562 : const uInt spwId = itsSpId[CDAId[s][0]];
1002 7562 : Int newId = dataDescription().match(spwId, itsPolId[s], thisDataId[s]);
1003 7562 : if (newId < 0) {
1004 31 : newId = addDataDescription(spwId, static_cast<uInt>(itsPolId[s]));
1005 : // this is a good place to add a hypercube as a new row in the Data
1006 : // Description subtable may specify a different data shape
1007 : //addHypercubes(nProducts[s], nChannels[s]);
1008 : }
1009 7562 : if (newId != thisDataId[s]) {
1010 155 : thisDataId[s] = newId;
1011 155 : itsNewScan = true;
1012 : }
1013 : }
1014 :
1015 :
1016 :
1017 : # if defined(AIPS_DEBUG)
1018 4570 : itsLog << LogIO::DEBUGGING;
1019 4570 : itsLog << CDAId.nelements() << " spectral ID's in this record "
1020 9140 : << "(correlator mode '" << VLAEnum::name(sda. correlatorMode())
1021 9140 : << "')" << endl;
1022 12132 : for (uInt s = 0; s < CDAId.nelements(); s++) {
1023 7562 : itsLog << " Id: " << itsSpId[CDAId[s][0]] << " is in CDA's [";
1024 15124 : for (uInt i = 0; i < CDAId[s].nelements(); i++) {
1025 7562 : itsLog << static_cast<Int>(CDAId[s][i])+1
1026 7562 : << ((i+1 < CDAId[s].nelements()) ? ", " : "]\n");
1027 : }
1028 15124 : for (uInt cda = 0; cda < CDAId[s].nelements(); cda++) {
1029 7562 : itsLog << " CDA: " << static_cast<Int>(CDAId[s][cda]) + 1
1030 7562 : << " contains " << polId[CDAId[s][cda]].nelements()
1031 22686 : << " polarizations [";
1032 33076 : for (uInt i = 0; i < polId[CDAId[s][cda]].nelements(); i++) {
1033 25514 : itsLog << polId[CDAId[s][cda]][i]
1034 25514 : << ((i+1 < polId[CDAId[s][cda]].nelements()) ? ", " : "] (");
1035 : }
1036 7562 : Vector<Stokes::StokesTypes> pol(itsRecord.polarisations(CDAId[s][cda]));
1037 33076 : for (uInt i = 0; i < pol.nelements(); i++) {
1038 51028 : itsLog << Stokes::name(pol(i))
1039 25514 : << ((i+1 < pol.nelements()) ? ", " : ")\n");
1040 : }
1041 7562 : }
1042 : }
1043 4570 : itsLog << LogIO::POST << LogIO::NORMAL;
1044 : #endif
1045 : // decide if this is a new scan
1046 4570 : if (itsNewScan) {
1047 161 : itsNewScan = false;
1048 161 : Int nextScan = itsScan[0];
1049 805 : for (uInt i = 1; i < maxSubarrays; i++) {
1050 644 : nextScan = max(nextScan, itsScan[i]);
1051 : }
1052 161 : itsScan[subArray] = nextScan + 1;
1053 : // flush any data to disk. This gives the user the opportunity to examine
1054 : // the MS while it is being filled. Doing it more frequently than once per
1055 : // scan starts to eat into the performance when filling from disk.
1056 161 : itsMS.flush();
1057 : }
1058 : // add empty rows
1059 4570 : const uInt nCorr = (nAnt*(nAnt-1))/2;
1060 4570 : uInt row = itsMS.nrow();
1061 4570 : uInt rowsPerSpId = nCorr+nAnt;
1062 4570 : if(!itsKeepAutoCorr)
1063 4442 : rowsPerSpId = nCorr;
1064 4570 : const uInt newRows = rowsPerSpId*nSpId;
1065 4570 : itsMS.addRow(newRows);
1066 : //extendHypercubes(nProducts, nChannels, rowsPerSpId);
1067 : // Some variables needed in assorted places from now on.
1068 4570 : Vector<Int> vecInt(newRows);
1069 4570 : const Double intTime = sda.intTime();
1070 : { // fill the columns where all the rows are identical and independent of the
1071 : // data description id
1072 4570 : const RefRows rows(row, row+newRows-1);
1073 4570 : vecInt = itsScan[subArray];
1074 4570 : scanNumber().putColumnCells(rows, vecInt);
1075 4570 : vecInt = itsFldId[subArray];
1076 4570 : fieldId().putColumnCells(rows, vecInt);
1077 4570 : vecInt = subArray;
1078 4570 : arrayId().putColumnCells(rows, vecInt);
1079 4570 : vecInt = itsMS.observation().nrow();
1080 4570 : observationId().putColumnCells(rows, vecInt);
1081 4570 : vecInt = 0;
1082 4570 : feed1().putColumnCells(rows, vecInt);
1083 4570 : feed2().putColumnCells(rows, vecInt);
1084 4570 : vecInt = -1;
1085 4570 : stateId().putColumnCells(rows, vecInt);
1086 4570 : processorId().putColumnCells(rows, vecInt);
1087 4570 : Vector<Double> vecDbl(newRows, intTime);
1088 4570 : exposure().putColumnCells(rows, vecDbl);
1089 4570 : interval().putColumnCells(rows, vecDbl);
1090 4570 : const MEpoch* mep = dynamic_cast<const MEpoch*>(itsFrame.epoch());
1091 4570 : AlwaysAssert( mep != 0, AipsError);
1092 4570 : vecDbl = mep->getValue().getTime("s").getValue();
1093 4570 : time().putColumnCells(rows, vecDbl);
1094 4570 : timeCentroid().putColumnCells(rows, vecDbl);
1095 4570 : }
1096 :
1097 : // Construct a bunch of variables that will be used inside the data writing
1098 : // loop
1099 4570 : Vector<Double> blUvw(3);
1100 4570 : Matrix<Complex> cData(nProducts[0], nChannels[0]);
1101 4570 : Matrix<Complex> modData(nProducts[0], nChannels[0]);
1102 4570 : Matrix<Complex> onePol;
1103 4570 : Matrix<Bool> flags(nProducts[0], nChannels[0]);
1104 4570 : Vector<Float> weights(nProducts[0]), sigmas(nProducts[0]);
1105 4570 : Cube<Bool> flagLevels(nProducts[0], nChannels[0], nCat);
1106 : VLAEnum::CDA cda;
1107 4570 : vecInt.resize(rowsPerSpId);
1108 :
1109 12132 : for (uInt s = 0; s < nSpId; s++) {
1110 7562 : const Block<VLAEnum::CDA>& usedCDAs = CDAId[s];
1111 7562 : const uInt nCDA = usedCDAs.nelements();
1112 7562 : DebugAssert(nCDA > 0, AipsError);
1113 7562 : const uInt nChan = nChannels[s];
1114 7562 : const uInt nPol = nProducts[s];
1115 7562 : const Double channelWidth = itsRecord.SDA().channelWidth(usedCDAs[0]);
1116 :
1117 : // fill the columns where all the rows are identical and dependent on the
1118 : // data description id
1119 : {
1120 7562 : const RefRows rows(row, row+rowsPerSpId-1);
1121 7562 : vecInt = thisDataId[s];
1122 7562 : dataDescId().putColumnCells(rows, vecInt);
1123 7562 : }
1124 :
1125 : // cache the online IF flags and nominal sensitivity of each antenna. It
1126 : // simplifies having to do it for each baseline later on.
1127 7562 : Block<Matrix<VLAEnum::IF> > whichIF(nCDA);
1128 7562 : Cube<Bool> antFlagLevels(maxIF, nAnt, 4, false);
1129 7562 : Matrix<Bool> antFlag(maxIF, nAnt, false);
1130 7562 : Matrix<Float> sens(maxIF, nAnt,0.333);
1131 7562 : Bool isScaledByNS(false);
1132 : { // First work ouk out which IF's are used by this spectral id.
1133 7562 : Block<Bool> usedIFs(maxIF, false);
1134 15124 : for (uInt c = 0; c < nCDA; c++) {
1135 : const Matrix<VLAEnum::IF>& curIF =
1136 7562 : whichIF[c] = sda.ifUsage(usedCDAs[c]);
1137 7562 : const uInt nCorr = curIF.ncolumn();
1138 33076 : for (uInt p = 0; p < nCorr; p++) {
1139 25514 : usedIFs[curIF(0, p)] = usedIFs[curIF(1, p)] = true;
1140 : }
1141 : }
1142 :
1143 : // For each antenna find the IF flags and sensitivity
1144 7562 : Int nAntIF(0);
1145 7562 : Int nAppAntIF(0);
1146 210414 : for (uInt a = 0; a < nAnt; a++) { // set the flag to true if data is Bad
1147 202852 : const VLAADA& ada = itsRecord.ADA(a);
1148 1014260 : for (uInt i = 0; i < maxIF; i++) {
1149 811408 : if (usedIFs[i]) {
1150 364420 : const VLAEnum::IF thisIF = static_cast<VLAEnum::IF>(i);
1151 364420 : const Float ns = ada.nominalSensitivity(thisIF);
1152 364420 : sens(i, a) = (ns > 1.0e-10) ? ns : 0.333;
1153 :
1154 : // count Ant/IF combos that have nom sens applied to amplitudes
1155 364420 : ++nAntIF;
1156 364420 : if (ada.nomSensApplied(thisIF,rca.revision())) ++nAppAntIF;
1157 364420 : const uInt status = ada.ifStatus(thisIF);
1158 364420 : if (status > 0) badData = true;
1159 364420 : if ((status&0x01) != 0) antFlagLevels(i, a, 0) = true;
1160 364420 : if ((status&0x02) != 0) antFlagLevels(i, a, 1) = true;
1161 364420 : if ((status&0x04) != 0) antFlag(i, a)=antFlagLevels(i, a, 2)=true;
1162 364420 : if ((status&0x08) != 0) antFlag(i, a)=antFlagLevels(i, a, 3)=true;
1163 364420 : if(!isHolo) {
1164 364420 : antFlag(i, a) |= shadowed[a];
1165 : }
1166 : }
1167 : }
1168 : }
1169 : // determine global state of nom sens application
1170 7562 : if (nAppAntIF==0) {
1171 0 : isScaledByNS=false;
1172 : // cout << "****DATA has NOT been scaled by NOMINAL SENSITIVITY*****************" << endl;
1173 : }
1174 : else {
1175 : // one or more ant/if combos indicate that NOM SENS has been applied
1176 : // in this case it is true for all even if not indicated for all.
1177 7562 : isScaledByNS=true;
1178 : // cout << "****DATA has been scaled by NOMINAL SENSITIVITY*****************" << endl;
1179 : }
1180 7562 : }
1181 :
1182 7562 : cData.resize(nPol, nChan);
1183 7562 : modData.resize(nPol, nChan);
1184 7562 : if(nPol==4){
1185 5984 : modData.row(0).set(1);
1186 5984 : modData.row(3).set(1);
1187 5984 : modData.row(1).set(0);
1188 5984 : modData.row(2).set(0);
1189 : }
1190 : else{
1191 1578 : modData.set(1);
1192 : }
1193 7562 : weights.resize(nPol);
1194 7562 : sigmas.resize(nPol);
1195 7562 : flags.resize(nPol, nChan);
1196 7562 : flagLevels.resize(nPol, nChan, nCat);
1197 7562 : if (!badData) {
1198 2795 : flags = false;
1199 2795 : flagLevels = false;
1200 : }
1201 7562 : if (nChan != 1 && nPol != 1 ) onePol.resize(1, nChan);
1202 7562 : const Slice allChan(0, nChan);
1203 :
1204 : // Fill in the correlations
1205 7562 : uInt b = 0;
1206 210414 : for (uInt a1 = 0; a1 < nAnt; a1++) {
1207 3025594 : for (uInt a2 = a1; a2 < nAnt; a2++) {
1208 2822742 : const Bool crossCorr = (a1 == a2) ? false : true;
1209 2822742 : if(crossCorr || (!crossCorr && itsKeepAutoCorr)){
1210 5246692 : for (uInt c = 0; c < nCDA; c++) {
1211 2623346 : cda = usedCDAs[c];
1212 2623346 : if (nChan == 1 || nPol == 1) {
1213 2623346 : if (crossCorr) {
1214 2619890 : itsRecord.CDA(cda).crossCorr(b).data(cData);
1215 : } else {
1216 3456 : itsRecord.CDA(cda).autoCorr(a1).data(cData);
1217 : }
1218 2623346 : if(nPol==4 && rotStokesOrder[s]){
1219 2100384 : Vector<Complex> olddata(4);
1220 2100384 : olddata=cData.column(0);
1221 10501920 : for (uInt kk=0; kk < 4; ++kk){
1222 8401536 : cData.column(0)(polNewOrder[s][kk])=olddata[kk];
1223 : }
1224 2100384 : }
1225 2623346 : } else {
1226 0 : DebugAssert(polId[cda].nelements() == 1, AipsError);
1227 0 : if (crossCorr) {
1228 0 : itsRecord.CDA(cda).crossCorr(b).data(onePol);
1229 : } else {
1230 0 : itsRecord.CDA(cda).autoCorr(a1).data(onePol);
1231 : }
1232 0 : const Slice curPol(polNewOrder[s][polId[cda][0]], 1);
1233 0 : cData(curPol, allChan) = onePol;
1234 : }
1235 2623346 : const Matrix<VLAEnum::IF>& curIF = whichIF[c];
1236 2623346 : if (nChan == 1) { // Continuum
1237 2100384 : DebugAssert(curIF.ncolumn() == 4, AipsError);
1238 2100384 : DebugAssert(polId[cda][0] == 0, AipsError);
1239 : uInt p;
1240 10501920 : for (uInt ip = 0; ip < 4; ip++) {
1241 8401536 : const VLAEnum::IF if0 = curIF(0, ip);
1242 8401536 : const VLAEnum::IF if1 = curIF(1, ip);
1243 :
1244 : // re-ordered output poln (ip-->p)!
1245 8401536 : p=polNewOrder[s][ip];
1246 :
1247 8401536 : const Double w = intTime * .12/10000.;
1248 : // The fudge factor of .12/10000 is to make the VLA filler
1249 : // consistent with AIPS. It is discussed in the .help file.
1250 8401536 : weights(p) = w * channelWidth;
1251 8401536 : sigmas(p) = 1.0/ sqrt(w * channelWidth);
1252 :
1253 : // If requested, apply Tsys scaling to data & weights
1254 8401536 : if (itsApplyTsys) {
1255 : // sens already guaranteed > 0.0
1256 8300448 : Float blsens = sens(if0, a1) * sens(if1, a2);
1257 :
1258 : // always apply to weights & sigma
1259 8300448 : weights(p)/=blsens;
1260 8300448 : sigmas(p)*=sqrt(blsens);
1261 :
1262 : // only apply to data if necessary (post-ModComp)
1263 8300448 : if (!isScaledByNS)
1264 0 : cData(p,0)*=sqrt(blsens);
1265 :
1266 : }
1267 : else
1268 : // Raw CCs requested
1269 101088 : if (isScaledByNS) {
1270 : // UN-correct data which was scaled on-line (e.g. pre-EVLA)
1271 101088 : Float blsens = sens(if0, a1) * sens(if1, a2);
1272 : // Only if correction is sane
1273 101088 : if (blsens>1.0e-10)
1274 101088 : cData(p,0)/=sqrt(blsens);
1275 : }
1276 :
1277 :
1278 8401536 : if (badData) {
1279 5872932 : flags(p, 0) = antFlag(if0, a1) || antFlag(if1, a2);
1280 29364660 : for (uInt l = 0; l < 4; l++) {
1281 23491728 : flagLevels(p, 0, l) =
1282 23491728 : antFlagLevels(if0, a1, l) || antFlagLevels(if1, a2, l);
1283 : }
1284 : // Don't flag holography data for apparent shadowing
1285 : // since we don't actually know if the data is
1286 : // shadowed
1287 5872932 : if(isHolo) {
1288 0 : flagLevels(p, 0, 4) = false;
1289 : }
1290 : else {
1291 5872932 : flagLevels(p, 0, 4) = shadowed[a1] || shadowed[a2];
1292 : }
1293 5872932 : flagLevels(p, 0, 5) = false;
1294 : }
1295 : }
1296 : } else {// spectral line
1297 522962 : DebugAssert(curIF.ncolumn() == 1, AipsError);
1298 522962 : const VLAEnum::IF if0 = curIF(0, 0);
1299 522962 : const VLAEnum::IF if1 = curIF(1, 0);
1300 :
1301 : // re-ordered output polarization!
1302 : // const uInt startPol = polId[cda][0];
1303 522962 : const uInt startPol =polNewOrder[s][polId[cda][0]];
1304 :
1305 522962 : const Double w = intTime * .12/10000.;
1306 : // The fudge factor of .12/10000 is to make the VLA filler
1307 : // consitant with AIPS. It is discussed in the .help file.
1308 522962 : weights(startPol) = w * channelWidth;
1309 522962 : sigmas(startPol) = 1.0/sqrt(w * channelWidth);
1310 :
1311 : // If requested, apply Tsys scaling to data & weights
1312 522962 : if (itsApplyTsys) {
1313 522962 : const Float blsens = sens(if0, a1) * sens(if1, a2);
1314 :
1315 : // always apply to weights & sigma
1316 522962 : weights(startPol)/=blsens;
1317 522962 : sigmas(startPol)*=sqrt(blsens);
1318 :
1319 : // only apply to data if necessary (post-ModComp)
1320 522962 : if (!isScaledByNS) {
1321 0 : Array<Complex> thisdat(cData(startPol,allChan));
1322 0 : thisdat*=sqrt(blsens);
1323 0 : }
1324 : }
1325 : else
1326 : // Raw CCs requested
1327 0 : if (isScaledByNS) {
1328 : // UN-correct data which was scaled on-line
1329 0 : const Float blsens = sens(if0, a1) * sens(if1, a2);
1330 0 : Array<Complex> thisdat(cData(startPol,allChan));
1331 0 : thisdat/=sqrt(blsens);
1332 0 : }
1333 :
1334 522962 : if (badData) {
1335 194382 : const Slice curPol(startPol, 1);
1336 194382 : flags(curPol, allChan) = antFlag(if0, a1) || antFlag(if1, a2);
1337 971910 : for (uInt l = 0; l < 4; l++) {
1338 1555056 : flagLevels(curPol, allChan, l) =
1339 2332584 : antFlagLevels(if0, a1, l) || antFlagLevels(if1, a2, l);
1340 : }
1341 194382 : if(isHolo) {
1342 0 : flagLevels(curPol, allChan, 4) = false;
1343 : }
1344 : else {
1345 194382 : flagLevels(curPol, allChan, 4) = shadowed[a1] || shadowed[a2];
1346 : }
1347 194382 : flagLevels(curPol, allChan, 5) = false;
1348 : }
1349 : }
1350 : }
1351 :
1352 : // Some aips++ tools, in particular calibrator, require that the index
1353 : // in ANTENNA1 be less than or equal to the index in the ANTENNA2
1354 : // column. To accommodate this they are swapped here.
1355 2623346 : uInt ant1 = a1;
1356 2623346 : uInt ant2 = a2;
1357 2623346 : if (itsAntId[a1] > itsAntId[a2]) {
1358 1090400 : if(nPol < 4){
1359 224546 : cData=conj(cData);
1360 : }
1361 : else{
1362 865854 : cData.row(0)=conj(cData.row(0));
1363 865854 : cData.row(3)=conj(cData.row(3));
1364 : //R_iL_j has to be moved to L_jR_i
1365 1731708 : Vector<Complex> tmprldata=conj(cData.row(1));
1366 865854 : cData.row(1)=conj(cData.row(2));
1367 865854 : cData.row(2)=tmprldata;
1368 865854 : Float tmpwt=weights(1);
1369 865854 : weights(1)=weights(2);
1370 865854 : weights(2)=tmpwt;
1371 865854 : tmpwt=sigmas(1);
1372 865854 : sigmas(1)=sigmas(2);
1373 865854 : sigmas(2)=tmpwt;
1374 865854 : Vector<Bool> tmpflg=flags.row(1);
1375 865854 : flags.row(1)=flags.row(2);
1376 865854 : flags.row(2)=tmpflg;
1377 865854 : }
1378 1090400 : ant1 = a2;
1379 1090400 : ant2 = a1;
1380 : }
1381 :
1382 2623346 : weight().put(row, weights);
1383 2623346 : sigma().put(row, sigmas);
1384 2623346 : flag().put(row, flags);
1385 2623346 : flagCategory().put(row, flagLevels);
1386 2623346 : if (badData) {
1387 1662615 : flagRow().put(row, allEQ(flags, true));
1388 : } else {
1389 960731 : flagRow().put(row, false);
1390 : }
1391 :
1392 :
1393 2623346 : data().put(row, cData);
1394 2623346 : correctedData().put(row,cData);
1395 2623346 : modelData().put(row, modData);
1396 2623346 : uInt a1Index = ant1*3, a2Index = ant2*3;
1397 2623346 : if(isHolo) {
1398 0 : blUvw=0.0;
1399 : }
1400 : else {
1401 2623346 : blUvw(0) = antUvw[a1Index] - antUvw[a2Index]; a1Index++; a2Index++;
1402 2623346 : blUvw(1) = antUvw[a1Index] - antUvw[a2Index]; a1Index++; a2Index++;
1403 2623346 : blUvw(2) = antUvw[a1Index] - antUvw[a2Index];
1404 : }
1405 2623346 : uvw().put(row, blUvw);
1406 2623346 : antenna1().put(row, itsAntId[ant1]);
1407 2623346 : antenna2().put(row, itsAntId[ant2]); row++;
1408 2623346 : if (crossCorr) b++;
1409 : }
1410 : }
1411 : }
1412 7562 : }
1413 4570 : return true;
1414 4570 : }
1415 :
1416 16 : MeasurementSet VLAFiller::
1417 : getMS(const Path& tableName, const Bool overwrite) {
1418 16 : if (overwrite == false && File(tableName).exists()) {
1419 7 : return openMS(tableName);
1420 : } else {
1421 9 : return emptyMS(tableName, overwrite);
1422 : }
1423 : }
1424 :
1425 9 : MeasurementSet VLAFiller::
1426 : emptyMS(const Path& tableName, const Bool overwrite) {
1427 9 : AlwaysAssert(tableName.isValid(), AipsError);
1428 9 : AlwaysAssert(File(tableName.dirName()).isWritable(), AipsError);
1429 :
1430 : // Add all the required columns
1431 9 : TableDesc msDesc = MeasurementSet::requiredTableDesc();
1432 : // Add the data column (as it is an an optional one)
1433 9 : MeasurementSet::addColumnToDesc(msDesc, MeasurementSet::DATA, 2);
1434 : //Scratch columns
1435 9 : MeasurementSet::addColumnToDesc(msDesc, MeasurementSet::MODEL_DATA, 2);
1436 9 : MeasurementSet::addColumnToDesc(msDesc, MeasurementSet::CORRECTED_DATA, 2);
1437 :
1438 : // Add the tiled id column indices
1439 : /*
1440 : msDesc.addColumn(ScalarColumnDesc<Int>(dataTileId.fieldName(),
1441 : "Index for Data tiling"));
1442 : msDesc.addColumn(ScalarColumnDesc<Int>(sigmaTileId.fieldName(),
1443 : "Index for Sigma tiling"));
1444 : msDesc.addColumn(ScalarColumnDesc<Int>(flagTileId.fieldName(),
1445 : "Index for Flag Category tiling"));
1446 :
1447 : msDesc.addColumn(ScalarColumnDesc<Int>(modDataTileId.fieldName(),
1448 : "Index for Model Data tiling"));
1449 : msDesc.addColumn(ScalarColumnDesc<Int>(corrDataTileId.fieldName(),
1450 : "Index for Corrected Data tiling"));
1451 : msDesc.addColumn(ScalarColumnDesc<Int>(chanFlagTileId.fieldName(),
1452 : "Index for Flag tiling"));
1453 : */
1454 : // setup hypercolumns for the data/flag/flag_catagory/sigma & weight columns.
1455 : {
1456 9 : Vector<String> dataCols(1);
1457 9 : dataCols(0) = MeasurementSet::columnName(MeasurementSet::DATA);
1458 9 : const Vector<String> coordCols(0);
1459 9 : const Vector<String> idCols(1, dataTileId.fieldName());
1460 : // msDesc.defineHypercolumn(dataCol, 3, dataCols, coordCols, idCols);
1461 9 : msDesc.defineHypercolumn(dataCol, 3, dataCols);
1462 9 : }
1463 : {
1464 9 : Vector<String> dataCols(1);
1465 9 : dataCols(0) = MeasurementSet::columnName(MeasurementSet::MODEL_DATA);
1466 9 : const Vector<String> coordCols(0);
1467 9 : const Vector<String> idCols(1, modDataTileId.fieldName());
1468 : // msDesc.defineHypercolumn(modDataCol, 3, dataCols, coordCols, idCols);
1469 9 : msDesc.defineHypercolumn(modDataCol, 3, dataCols);
1470 9 : }
1471 : {
1472 9 : Vector<String> dataCols(1);
1473 9 : dataCols(0) = MeasurementSet::columnName(MeasurementSet::CORRECTED_DATA);
1474 9 : const Vector<String> coordCols(0);
1475 9 : const Vector<String> idCols(1, corrDataTileId.fieldName());
1476 : //msDesc.defineHypercolumn(corrDataCol, 3, dataCols, coordCols, idCols);
1477 9 : msDesc.defineHypercolumn(corrDataCol, 3, dataCols);
1478 9 : }
1479 : {
1480 9 : Vector<String> dataCols(1);
1481 9 : dataCols(0) = MeasurementSet::columnName(MeasurementSet::FLAG);
1482 9 : const Vector<String> coordCols(0);
1483 9 : const Vector<String> idCols(1, chanFlagTileId.fieldName());
1484 : //msDesc.defineHypercolumn(chanFlagCol, 3, dataCols, coordCols, idCols);
1485 9 : msDesc.defineHypercolumn(chanFlagCol, 3, dataCols);
1486 9 : }
1487 : {
1488 9 : Vector<String> dataCols(1);
1489 9 : dataCols(0) = MeasurementSet::columnName(MeasurementSet::SIGMA);
1490 9 : const Vector<String> coordCols(0);
1491 9 : const Vector<String> idCols(1, sigmaTileId.fieldName());
1492 : //msDesc.defineHypercolumn(sigmaCol, 2, dataCols, coordCols, idCols);
1493 9 : msDesc.defineHypercolumn(sigmaCol, 2, dataCols);
1494 9 : }
1495 : //sigma and weight unbound as of moving to tiledshapestman
1496 : {
1497 9 : Vector<String> dataCols(1);
1498 9 : dataCols(0) = MeasurementSet::columnName(MeasurementSet::WEIGHT);
1499 9 : const Vector<String> coordCols(0);
1500 9 : const Vector<String> idCols(1, sigmaTileId.fieldName());
1501 9 : msDesc.defineHypercolumn("TiledWgtCol", 2, dataCols);
1502 9 : }
1503 :
1504 : {
1505 : const Vector<String> dataCols
1506 9 : (1, MeasurementSet::columnName(MeasurementSet::FLAG_CATEGORY));
1507 9 : const Vector<String> coordCols(0);
1508 9 : const Vector<String> idCols(1, flagTileId.fieldName());
1509 : // msDesc.defineHypercolumn(flagCol, 4, dataCols, coordCols, idCols);
1510 9 : msDesc.defineHypercolumn(flagCol, 4, dataCols);
1511 9 : }
1512 :
1513 9 : Table::TableOption option = Table::NewNoReplace;
1514 9 : if (overwrite) option = Table::New;
1515 9 : SetupNewTable newMS(tableName.originalName(), msDesc, option);
1516 :
1517 : // setup storage managers. Use the incremental storage manager for
1518 : // columns where the data is likely to be the same for more than four
1519 : // rows at a time.
1520 : {
1521 9 : IncrementalStMan incrMan("Incremental data manager");
1522 9 : newMS.bindColumn(MeasurementSet::
1523 : columnName(MeasurementSet::ARRAY_ID), incrMan);
1524 9 : newMS.bindColumn(MeasurementSet::
1525 : columnName(MeasurementSet::EXPOSURE), incrMan);
1526 9 : newMS.bindColumn(MeasurementSet::
1527 : columnName(MeasurementSet::FEED1), incrMan);
1528 9 : newMS.bindColumn(MeasurementSet::
1529 : columnName(MeasurementSet::FEED2), incrMan);
1530 9 : newMS.bindColumn(MeasurementSet::
1531 : columnName(MeasurementSet::FIELD_ID), incrMan);
1532 9 : newMS.bindColumn(MeasurementSet::
1533 : columnName(MeasurementSet::FLAG_ROW), incrMan);
1534 9 : newMS.bindColumn(MeasurementSet::
1535 : columnName(MeasurementSet::INTERVAL), incrMan);
1536 9 : newMS.bindColumn(MeasurementSet::
1537 : columnName(MeasurementSet::OBSERVATION_ID), incrMan);
1538 9 : newMS.bindColumn(MeasurementSet::
1539 : columnName(MeasurementSet::PROCESSOR_ID), incrMan);
1540 9 : newMS.bindColumn(MeasurementSet::
1541 : columnName(MeasurementSet::SCAN_NUMBER), incrMan);
1542 9 : newMS.bindColumn(MeasurementSet::
1543 : columnName(MeasurementSet::STATE_ID), incrMan);
1544 9 : newMS.bindColumn(MeasurementSet::
1545 : columnName(MeasurementSet::TIME), incrMan);
1546 9 : newMS.bindColumn(MeasurementSet::
1547 : columnName(MeasurementSet::TIME_CENTROID), incrMan);
1548 9 : }
1549 : // Use a 1 MB tile size
1550 : // The tile length should not be shorter than (npol, nchan)
1551 : // in the first 2 dimensions, otherwise performance suffers
1552 : // due to bookkeeping overhead [yes, that implies that using
1553 : // tiling is pretty pointless for such small data... But then
1554 : // you will suffer from the horrors of BucketCache.]. Therefore
1555 : // set the length to the maximum. Unfortunately npol, nchan is
1556 : // unknown here, so set the lengths to (4,128); the consequence
1557 : // of hardcoding this is that the real tile size will be less
1558 : // than 1MB, in fact 1MB/(4*128) = 2KB for npol=nchan columns.
1559 9 : IPosition tileShape(3, 4, 128, 256);
1560 :
1561 : {
1562 : //TiledDataStMan dataMan(dataCol);
1563 9 : TiledShapeStMan dataMan(dataCol, tileShape);
1564 9 : newMS.bindColumn(MeasurementSet::
1565 : columnName(MeasurementSet::DATA), dataMan);
1566 : // newMS.bindColumn(dataTileId.fieldName(), dataMan);
1567 9 : }
1568 : {
1569 : //TiledDataStMan dataMan(modDataCol);
1570 9 : TiledShapeStMan dataMan(modDataCol, tileShape);
1571 9 : newMS.bindColumn(MeasurementSet::
1572 : columnName(MeasurementSet::MODEL_DATA), dataMan);
1573 : // newMS.bindColumn(modDataTileId.fieldName(), dataMan);
1574 9 : }
1575 : {
1576 : // TiledDataStMan dataMan(corrDataCol);
1577 9 : TiledShapeStMan dataMan(corrDataCol, tileShape);
1578 9 : newMS.bindColumn(MeasurementSet::
1579 : columnName(MeasurementSet::CORRECTED_DATA), dataMan);
1580 : // newMS.bindColumn(corrDataTileId.fieldName(), dataMan);
1581 9 : }
1582 : {
1583 : //TiledDataStMan dataMan(chanFlagCol);
1584 : TiledShapeStMan dataMan(chanFlagCol,
1585 0 : IPosition(3,
1586 9 : tileShape(0),
1587 9 : tileShape(1),
1588 9 : tileShape(2)*64)
1589 9 : );
1590 9 : newMS.bindColumn(MeasurementSet::
1591 : columnName(MeasurementSet::FLAG), dataMan);
1592 : // newMS.bindColumn(chanFlagTileId.fieldName(), dataMan);
1593 9 : }
1594 : {
1595 : //TiledDataStMan dataMan(sigmaCol);
1596 9 : TiledShapeStMan dataMan(sigmaCol, IPosition(2,tileShape(0), tileShape(1) * tileShape(2)));
1597 9 : newMS.bindColumn(MeasurementSet::
1598 : columnName(MeasurementSet::SIGMA), dataMan);
1599 :
1600 : //Hmmm before weight and sigma were bound by the same stman..
1601 18 : TiledShapeStMan dataMan2("TiledWgtCol", IPosition(2,tileShape(0), tileShape(1) * tileShape(2)));
1602 :
1603 9 : newMS.bindColumn(MeasurementSet::
1604 : columnName(MeasurementSet::WEIGHT), dataMan2);
1605 : // newMS.bindColumn(sigmaTileId.fieldName(), dataMan);
1606 9 : }
1607 :
1608 : {
1609 : //TiledDataStMan dataMan(flagCol);
1610 9 : TiledShapeStMan dataMan(flagCol, IPosition(4,tileShape(0),tileShape(1), 1,tileShape(2)));
1611 9 : newMS.bindColumn(MeasurementSet::
1612 : columnName(MeasurementSet::FLAG_CATEGORY), dataMan);
1613 : // newMS.bindColumn(flagTileId.fieldName(), dataMan);
1614 9 : }
1615 :
1616 : {
1617 18 : TiledColumnStMan tiledStUVW("TiledUVW",IPosition(2, 3, tileShape(0) * tileShape(1) * tileShape(2) * 2 / 3));
1618 9 : newMS.bindColumn(MS::columnName(MS::UVW),tiledStUVW);
1619 9 : }
1620 : // The standard storage manager is the default manager but by default it only
1621 : // creates a bucket for every 32 rows. Thats too small for the columns in the
1622 : // main table of a measurement set. So I'll explicitly bind these columns
1623 : // here with a bucket size of (351+27)*128 bytes
1624 : {
1625 9 : StandardStMan stMan("Standard data manager", 32768);
1626 9 : newMS.bindColumn(MeasurementSet::
1627 : columnName(MeasurementSet::ANTENNA1), stMan);
1628 9 : newMS.bindColumn(MeasurementSet::
1629 : columnName(MeasurementSet::ANTENNA2), stMan);
1630 9 : newMS.bindColumn(MeasurementSet::
1631 : columnName(MeasurementSet::DATA_DESC_ID), stMan);
1632 9 : }
1633 :
1634 : // Finally create the MeasurementSet.
1635 9 : MeasurementSet ms(newMS);
1636 :
1637 : { // Set the TableInfo
1638 9 : TableInfo& info(ms.tableInfo());
1639 9 : info.setType(TableInfo::type(TableInfo::MEASUREMENTSET));
1640 9 : info.setSubType(String("VLA"));
1641 : info.readmeAddLine
1642 9 : ("This is a MeasurementSet Table holding measurements from the VLA");
1643 9 : info.readmeAddLine("radio synthesis array (operated by NRAO)");
1644 : }
1645 : {//Create the SOURCE subtable
1646 :
1647 9 : TableDesc sourceTD=MSSource::requiredTableDesc();
1648 9 : MSSource::addColumnToDesc(sourceTD, MSSource::REST_FREQUENCY);
1649 9 : MSSource::addColumnToDesc(sourceTD, MSSource::SYSVEL);
1650 9 : MSSource::addColumnToDesc(sourceTD, MSSource::TRANSITION);
1651 9 : SetupNewTable sourceSetup(ms.sourceTableName(),sourceTD,option);
1652 18 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::SOURCE),
1653 18 : Table(sourceSetup,0));
1654 9 : }
1655 : // create the required subtables.
1656 9 : ms.createDefaultSubtables(option);
1657 : //
1658 : { // add optional column to SPECTRAL_WINDOW. Added by GYL
1659 18 : ms.spectralWindow().addColumn(
1660 18 : ScalarColumnDesc<Int>(
1661 : MSSpectralWindow::columnName(MSSpectralWindow::DOPPLER_ID),
1662 : MSSpectralWindow::columnStandardComment(
1663 : MSSpectralWindow::DOPPLER_ID)));
1664 :
1665 : }
1666 : { // Create the DOPPLER subtable. Added by GYL
1667 9 : TableDesc dopplerTD=MSDoppler::requiredTableDesc();
1668 : //MSDoppler::addColumnToDesc(dopplerTD, MSDoppler::VELDEF);
1669 9 : SetupNewTable dopplerSetup(ms.dopplerTableName(),dopplerTD,option);
1670 18 : ms.rwKeywordSet().defineTable(MS::keywordName(MS::DOPPLER),
1671 18 : Table(dopplerSetup,0));
1672 9 : }
1673 :
1674 : // Adjust the Measure references to ones used by the VLA.
1675 : {
1676 9 : MSColumns msc(ms);
1677 9 : msc.setEpochRef(MEpoch::IAT);
1678 9 : msc.setDirectionRef(MDirection::J2000);
1679 9 : msc.uvwMeas().setDescRefCode(Muvw::J2000);
1680 9 : msc.antenna().setPositionRef(MPosition::ITRF);
1681 9 : msc.antenna().setOffsetRef(MPosition::ITRF);
1682 : { // Put the right values into the CATEGORY keyword
1683 9 : Vector<String> categories(nCat);
1684 9 : categories(0) = "ONLINE_1";
1685 9 : categories(1) = "ONLINE_2";
1686 9 : categories(2) = "ONLINE_4";
1687 9 : categories(3) = "ONLINE_8";
1688 9 : categories(4) = "SHADOW";
1689 9 : categories(5) = "FLAG_CMD";
1690 9 : msc.setFlagCategories(categories);
1691 9 : }
1692 9 : }
1693 18 : return ms;
1694 9 : }
1695 :
1696 7 : MeasurementSet VLAFiller::
1697 : openMS(const Path& tableName, const Bool readonly) {
1698 7 : const String& msName = tableName.absoluteName();
1699 7 : if (!Table::isReadable(msName)) {
1700 0 : throw(AipsError(String("VLAFiller::openMS(...) - cannot read ") + msName +
1701 0 : String(" because it does not exist or is not a table.")));
1702 : }
1703 7 : if (!readonly && !Table::isWritable(msName)) {
1704 0 : throw(AipsError(String("VLAFiller::openMS(...) - cannot write to ") +
1705 0 : msName + "."));
1706 : }
1707 : {
1708 7 : const TableInfo info = TableUtil::tableInfo(msName);
1709 7 : if (info.type() != TableInfo::type(TableInfo::MEASUREMENTSET)) {
1710 0 : throw(AipsError(String("VLAFiller::openMS(...) - the table ") +
1711 0 : msName + String(" is not a measurement set.")));
1712 : }
1713 7 : if (info.subType() != "VLA") {
1714 0 : throw(AipsError(String("VLAFiller::openMS(...) - the table ") +
1715 0 : msName + String(" is not a VLA measurement set.")));
1716 : }
1717 : {
1718 7 : const Table t(msName);
1719 7 : const TableRecord& keys = t.keywordSet();
1720 7 : const String versionString("MS_VERSION");
1721 7 : const RecordFieldId versionKey(versionString);
1722 7 : if (!keys.isDefined(versionString) ||
1723 14 : keys.dataType(versionKey) != TpFloat ||
1724 7 : !near(keys.asFloat(versionKey), 2.0f)) {
1725 0 : throw(AipsError(String("VLAFiller::openMS(...) - the table ") +
1726 0 : msName +
1727 0 : String(" is not a version 2 measurement set.")));
1728 : }
1729 7 : }
1730 7 : }
1731 7 : const Table::TableOption openOption = readonly ? Table::Old : Table::Update;
1732 7 : return MeasurementSet(msName, openOption);
1733 : }
1734 :
1735 0 : void VLAFiller::logCurrentRecord(IterationStatus& counts) {
1736 0 : if (itsRecord.isValid()) {
1737 0 : const uInt curRow = nrow();
1738 0 : const MEpoch obsTime = timeMeas()(curRow-1);
1739 0 : itsLog << "Record " << counts.nVLARecords
1740 : << " was observed at "
1741 0 : << MVTime(obsTime.getValue().getTime()).string(MVTime::YMD)
1742 0 : << " and has " << curRow - counts.nRows
1743 0 : << " rows of data" << LogIO::POST;
1744 0 : }
1745 0 : }
1746 :
1747 4570 : void VLAFiller::logChanges(IterationStatus& counts) {
1748 : {
1749 4570 : const String& curProject = itsRecord.SDA().obsId();
1750 4570 : if (counts.lastProject != curProject) {
1751 11 : counts.lastProject = curProject;
1752 : // log output clean-up (CAS-208)
1753 : //itsLog << "Project changed to " << curProject << LogIO::POST;
1754 11 : itsLog << "Project " << curProject;
1755 :
1756 11 : const uInt lastRow = nrow() - 1;
1757 11 : const MEpoch obsTime = timeMeas()(lastRow);
1758 11 : itsLog << " starting at "
1759 22 : << MVTime(obsTime.getValue().getTime()).string(MVTime::YMD)
1760 11 : << " (" << obsTime.getRefString() << ")"
1761 22 : << LogIO::POST;
1762 11 : }
1763 4570 : }
1764 : {
1765 4570 : const String& curObsMode = itsRecord.SDA().obsMode();
1766 4570 : if (counts.lastObsMode != curObsMode) {
1767 11 : counts.lastObsMode = curObsMode;
1768 : //itsLog << "ObsMode changed to " << itsRecord.SDA().obsModeFull()
1769 22 : itsLog << "ObsMode: " << itsRecord.SDA().obsModeFull()
1770 22 : << LogIO::POST;
1771 : }
1772 4570 : }
1773 4570 : const Int subArray = arrayId()(nrow() - 1);
1774 : {
1775 4570 : const MSAntennaColumns& ant = antenna();
1776 4570 : const uInt nAnt = itsRecord.RCA().nAntennas();
1777 4570 : Bool changed = (nAnt != counts.lastAnt[subArray].nelements());
1778 : {
1779 4570 : Int a = nAnt;
1780 126343 : while (!changed && a > 0) {
1781 121773 : a--;
1782 121773 : changed = (itsAntId[a] != static_cast<Int>(counts.lastAnt[subArray][a]));
1783 : }
1784 : }
1785 4570 : if (changed) {
1786 : // some log output clean-ups (CAS-208)
1787 : //itsLog << "Array configuration";
1788 : //itsLog << " for sub-array " << subArray + 1;
1789 11 : itsLog << "Sub-array " << subArray + 1
1790 11 : << LogIO::POST;
1791 : //if (counts.lastAnt[subArray].nelements() != 0) itsLog << " changed";
1792 11 : const uInt lastRow = nrow() - 1;
1793 11 : const MEpoch obsTime = timeMeas()(lastRow);
1794 : //itsLog << " at "
1795 : // << MVTime(obsTime.getValue().getTime()).string(MVTime::YMD)
1796 : // << " (" << obsTime.getRefString() << ")"
1797 : // << LogIO::POST;
1798 11 : counts.lastAnt[subArray].resize(nAnt, true);
1799 306 : for (uInt i = 0; i < nAnt; i++) {
1800 295 : const uInt a = static_cast<uInt>(itsAntId[i]);
1801 295 : itsLog << "Station: " << ant.station()(a)
1802 590 : << " Antenna: " << ant.name()(a);
1803 : //if (counts.lastAnt[subArray][i] != a) {
1804 : // itsLog << " \tNEW";
1805 : //}
1806 295 : counts.lastAnt[subArray][i] = a;
1807 295 : itsLog << LogIO::POST;
1808 : }
1809 11 : }
1810 4570 : itsLog << LogIO::POST;
1811 : }
1812 : {
1813 4570 : const MSDataDescColumns& dd = dataDescription();
1814 4570 : const MSSpWindowColumns& spw = spectralWindow();
1815 4570 : const MSPolarizationColumns& pol = polarization();
1816 4570 : const Block<Int>& thisDataId = itsDataId[subArray];
1817 4570 : Block<Int>& lastSpw = counts.lastSpw[subArray];
1818 4570 : Block<Int>& lastPol = counts.lastPol[subArray];
1819 13710 : for (uInt d = 0; d < lastSpw.nelements(); d++) {
1820 9140 : if (d < thisDataId.nelements()) {
1821 7562 : const Int curDD = thisDataId[d];
1822 7562 : const Int curSpw = dd.spectralWindowId()(curDD);
1823 7562 : bool logPostPending = false;
1824 7562 : if (lastSpw[d] != curSpw) {
1825 155 : logPostPending = true;
1826 155 : lastSpw[d] = curSpw;
1827 : // also reset lastPol so the pol is logged for this spw
1828 155 : lastPol[d] = -1;
1829 : //itsLog << "Spectral window for IF#"
1830 155 : itsLog << "Spectral window "
1831 155 : << spw.ifConvChain()(curSpw) + 1
1832 : << ": "
1833 : // << " (on sub-array " << subArray + 1 << ")"
1834 : // << " changed to "
1835 310 : << spw.name()(curSpw);
1836 155 : const MEpoch obsTime = timeMeas()(nrow()-1);
1837 : //itsLog << " at "
1838 : // << MVTime(obsTime.getValue().getTime()).string(MVTime::YMD)
1839 : // << " (" << obsTime.getRefString() << ")";
1840 : //if (counts.nSpw != spw.nrow() &&
1841 : // static_cast<uInt>(curSpw) >= counts.nSpw) {
1842 : // counts.nSpw = curSpw + 1;
1843 : // itsLog << " NEW";
1844 : //}
1845 : // itsLog << LogIO::POST;
1846 155 : }
1847 7562 : const Int curPol = dd.polarizationId()(curDD);
1848 7562 : if (lastPol[d] != curPol) {
1849 155 : logPostPending = true;
1850 155 : lastPol[d] = curPol;
1851 : //itsLog << "Polarization setup for IF#"
1852 : // << spw.ifConvChain()(curSpw) + 1
1853 : // << " (on sub-array " << subArray + 1 << ")"
1854 : // << " changed to ";
1855 155 : const Vector<Int> corr = pol.corrType()(curPol);
1856 155 : const Int nCorr = corr.nelements();
1857 155 : itsLog << " ";
1858 605 : for (Int c = 0; c < nCorr - 1; c++) {
1859 450 : itsLog << Stokes::name(Stokes::type(corr(c))) << ", ";
1860 : }
1861 155 : itsLog << Stokes::name(Stokes::type(corr(nCorr-1)));
1862 155 : const MEpoch obsTime = timeMeas()(nrow()-1);
1863 : //itsLog << " at "
1864 : // << MVTime(obsTime.getValue().getTime()).string(MVTime::YMD)
1865 : // << " (" << obsTime.getRefString() << ")";
1866 : //if (counts.nPol != pol.nrow() &&
1867 : // static_cast<uInt>(curPol) >= counts.nPol) {
1868 : // counts.nPol = curPol + 1;
1869 : // itsLog << " NEW";
1870 : //}
1871 : //itsLog << LogIO::POST;
1872 155 : }
1873 7562 : if (logPostPending) {
1874 155 : itsLog << LogIO::POST;
1875 155 : logPostPending = false;
1876 : }
1877 : } else {
1878 1578 : lastSpw[d] = -1;
1879 1578 : lastPol[d] = -1;
1880 : }
1881 : }
1882 : }
1883 : {
1884 4570 : const MSFieldColumns& fld = field();
1885 4570 : const Int thisFld = itsFldId[subArray];
1886 4570 : if (counts.lastFld[subArray] != thisFld) {
1887 112 : counts.lastFld[subArray] = thisFld;
1888 : //itsLog << "Field changed to ";
1889 : {
1890 112 : const String& fldName = fld.name()(thisFld);
1891 112 : if (fldName.length() > 0) {
1892 112 : itsLog << fldName;
1893 : } else {
1894 0 : Array<MDirection> amd;
1895 0 : const Unit rad("rad");
1896 0 : fld.referenceDirMeasCol().get(thisFld, amd, true);
1897 0 : const MDirection& md = amd(IPosition(amd.ndim(), 0));
1898 0 : const MVDirection& mdv = md.getValue();
1899 0 : const MVTime ra(mdv.getLong(rad));
1900 0 : const MVAngle dec(mdv.getLat(rad));
1901 0 : itsLog << "(" << ra.string(MVTime::TIME, 6) << ", "
1902 0 : << dec.string(MVTime::ANGLE, 6) << ")";
1903 0 : itsLog << " " << md.getRefString();
1904 0 : }
1905 112 : }
1906 112 : const uInt lastRow = nrow() - 1;
1907 112 : const MEpoch obsTime = timeMeas()(lastRow);
1908 : //itsLog << " (on sub-array " << subArray+1 << ") at "
1909 112 : itsLog << " "
1910 : // << MVTime(obsTime.getValue().getTime()).string(MVTime::YMD)
1911 112 : << MVTime(obsTime.getValue().getTime()).string(MVTime::YMD);
1912 : // << " (" << obsTime.getRefString() << ")";
1913 : //if (counts.nFld != fld.nrow() &&
1914 : // static_cast<uInt>(thisFld) >= counts.nFld) {
1915 : // counts.nFld = fld.nrow();
1916 : // itsLog << " NEW";
1917 : //}
1918 112 : itsLog << LogIO::POST;
1919 112 : }
1920 : }
1921 4570 : }
1922 :
1923 :
1924 16 : void VLAFiller::summarise() {
1925 16 : itsLog << LogIO::NORMAL;
1926 16 : itsLog << "Finished filling the measurement set." << endl;
1927 16 : itsLog << "The measurement set contains " << nrow() << " rows." << endl;
1928 16 : itsLog << "The antenna sub-table contains "
1929 16 : << antenna().nrow() << " entries" << endl;
1930 16 : itsLog << "The field sub-table contains "
1931 16 : << field().nrow() << " entries" << endl;
1932 16 : itsLog << "The spectral window sub-table contains "
1933 16 : << spectralWindow().nrow() << " entries" << endl;
1934 16 : itsLog << "The polarization sub-table contains "
1935 16 : << polarization().nrow() << " entries" << endl;
1936 16 : itsLog << "The pointing sub-table contains "
1937 16 : << pointing().nrow() << " entries" << endl;
1938 16 : itsLog << LogIO::POST;
1939 16 : }
1940 :
1941 261 : uInt VLAFiller::addAntenna(const MPosition& antPos, uInt whichAnt) {
1942 :
1943 261 : MSAntennaColumns& ant = antenna();
1944 261 : const uInt newRow = ant.nrow();
1945 261 : itsMS.antenna().addRow(1);
1946 :
1947 261 : String leAntName;
1948 261 : if(!itsEVLAisOn){
1949 : // ignore the frontend temperature naming
1950 261 : leAntName=itsRecord.ADA(whichAnt).antName(false);
1951 261 : if(itsNewAntName){
1952 207 : leAntName=String("VA")+leAntName;
1953 : }
1954 : }
1955 : else{
1956 0 : leAntName=itsRecord.ADA(whichAnt).antName(itsNewAntName);
1957 : }
1958 :
1959 261 : ant.name().put(newRow, leAntName);
1960 :
1961 261 : ant.station().put(newRow, itsRecord.ADA(whichAnt).padName());
1962 261 : ant.type().put(newRow, "GROUND-BASED");
1963 261 : ant.mount().put(newRow, "ALT-AZ");
1964 261 : ant.positionMeas().put(newRow, antPos);
1965 261 : ant.offset().put(newRow, Vector<Double>(3, 0.0));
1966 261 : ant.dishDiameter().put(newRow, 25.0);
1967 261 : ant.flagRow().put(newRow, false);
1968 261 : return newRow;
1969 261 : }
1970 :
1971 :
1972 0 : uInt VLAFiller::addPointing(const MDirection& antDir,
1973 : const MDirection& fieldDir,
1974 : uInt whichAnt) {
1975 0 : MSPointingColumns& pointingCol = pointing();
1976 0 : const uInt newRow = pointingCol.nrow();
1977 0 : itsMS.pointing().addRow(1);
1978 :
1979 : // Should the antennaId just be whichAnt?
1980 : // (I.e., the MS ANT Id and not the VLA ant name integer?)
1981 0 : pointingCol.antennaId().put(newRow, itsRecord.ADA(whichAnt).antId());
1982 0 : const MEpoch* mepPtr = dynamic_cast<const MEpoch*>(itsFrame.epoch());
1983 0 : pointingCol.timeMeas().put(newRow, *mepPtr);
1984 0 : pointingCol.numPoly().put(newRow, 0);
1985 0 : pointingCol.timeOriginMeas().put(newRow, itsFrame.epoch());
1986 0 : pointingCol.interval().put(newRow, itsRecord.SDA().intTime());
1987 0 : Array<MDirection> bDir(IPosition(1,1));
1988 0 : bDir(IPosition(1,0))=antDir;
1989 0 : pointingCol.directionMeasCol().put(newRow, bDir);
1990 0 : bDir(IPosition(1,0))=fieldDir;
1991 0 : pointingCol.targetMeasCol().put(newRow, bDir);
1992 0 : pointingCol.tracking().put(newRow, true);
1993 0 : return newRow;
1994 0 : }
1995 :
1996 261 : void VLAFiller::addFeed(uInt whichAnt) {
1997 261 : MSFeedColumns& fd = feed();
1998 261 : const uInt newRow = fd.nrow();
1999 261 : itsMS.feed().addRow(1);
2000 :
2001 261 : fd.antennaId().put(newRow, whichAnt);
2002 261 : fd.feedId().put(newRow, 0);
2003 261 : fd.spectralWindowId().put(newRow, -1);
2004 261 : fd.time().put(newRow, 0.0);
2005 261 : fd.interval().put(newRow, 0.0);
2006 261 : fd.numReceptors().put(newRow, 2);
2007 261 : fd.beamId().put(newRow, -1);
2008 :
2009 261 : fd.beamOffset().put(newRow, Matrix<Double>(2,2, 0.0) );
2010 : {
2011 261 : Vector<String> pt(2);
2012 261 : pt(0) = "R"; pt(1) = "L";
2013 261 : fd.polarizationType().put(newRow, pt);
2014 261 : }
2015 : {
2016 261 : Matrix<Complex> resp(2, 2, Complex(0.0f, 0.0f));
2017 261 : resp.diagonal() = Complex(1.0f, 0.0f);
2018 261 : fd.polResponse().put(newRow, resp);
2019 261 : }
2020 261 : fd.position().put(newRow, Vector<Double>(3, 0.0));
2021 261 : fd.receptorAngle().put(newRow, Vector<Double>(2, 0.0));
2022 261 : }
2023 :
2024 30 : uInt VLAFiller::addField(const MDirection& dir) {
2025 30 : MSFieldColumns& fld = field();
2026 30 : const uInt newRow = fld.nrow();
2027 30 : itsMS.field().addRow(1);
2028 30 : fld.name().put(newRow, itsRecord.SDA().sourceName());
2029 30 : fld.code().put(newRow, itsRecord.SDA().calCode());
2030 30 : const MEpoch* mepPtr = dynamic_cast<const MEpoch*>(itsFrame.epoch());
2031 30 : fld.timeMeas().put(newRow, *mepPtr);
2032 30 : fld.numPoly().put(newRow, 0);
2033 : {
2034 30 : Vector<MDirection> aDir(1, dir);
2035 30 : fld.delayDirMeasCol().put(newRow, aDir);
2036 30 : fld.phaseDirMeasCol().put(newRow, aDir);
2037 30 : fld.referenceDirMeasCol().put(newRow, aDir);
2038 30 : }
2039 30 : fld.sourceId().put(newRow, newRow);
2040 30 : if (!fld.ephemerisId().isNull()) {
2041 0 : fld.ephemerisId().put(newRow, -1);
2042 : }
2043 30 : fld.flagRow().put(newRow, false);
2044 30 : return newRow;
2045 : }
2046 :
2047 31 : uInt VLAFiller::addDoppler( const VLAEnum::CDA cda ) {
2048 31 : const VLASDA& sda = itsRecord.SDA();
2049 31 : MSDopplerColumns& dopc = doppler();
2050 31 : const uInt newRow = dopc.nrow();
2051 31 : itsMS.doppler().addRow(1);
2052 31 : dopc.dopplerId().put(newRow, newRow );
2053 : // find the source_id
2054 31 : MSSourceColumns& srcc = source();
2055 31 : const uInt source_id = srcc.nrow() - 1;
2056 31 : dopc.sourceId().put(newRow, source_id );
2057 : // transition column in SOURCE subtable is not filled. So here we fill TRANSITION_ID with 0.
2058 31 : dopc.transitionId().put(newRow, 0);
2059 31 : if (sda.dopplerTracking( cda )) {
2060 6 : const MDoppler dop(Quantity(1.0*sda.radialVelocity( cda ), "m/s"), sda.dopplerDefn( cda ));
2061 3 : dopc.velDefMeas().put(newRow, dop );
2062 3 : }else {
2063 28 : dopc.velDefMeas().put(newRow, MDoppler(Quantity(0.0,"m/s"),MDoppler::RADIO));
2064 : }
2065 :
2066 31 : return newRow;
2067 : }
2068 :
2069 31 : uInt VLAFiller::addSpectralWindow(const VLAEnum::CDA cda,
2070 : const MFrequency& refFreq,
2071 : const uInt nChan,
2072 : const Double /*bandwidth*/,
2073 : const uInt ifChain) {
2074 31 : MSSpWindowColumns& spw = spectralWindow();
2075 31 : const uInt newRow = spw.nrow();
2076 31 : itsMS.spectralWindow().addRow(1);
2077 : //cout.precision(8);
2078 31 : spw.refFrequencyMeas().put(newRow, refFreq);
2079 31 : spw.numChan().put(newRow, nChan);
2080 31 : spw.ifConvChain().put(newRow, ifChain);
2081 : // write doppler_id to spectral_window
2082 31 : MSDopplerColumns& dopc = doppler();
2083 31 : const uInt doppler_id = dopc.nrow() - 1;
2084 31 : spw.dopplerId().put( newRow, doppler_id );
2085 :
2086 31 : const VLASDA& sda = itsRecord.SDA();
2087 31 : const Double chanWidth = sda.channelWidth(cda);
2088 : const MFrequency::Types itsFreqType =
2089 31 : MFrequency::castType(refFreq.getRef().getType());
2090 :
2091 31 : Vector<Double> chanFreq(nChan);
2092 31 : indgen(chanFreq, sda.edgeFrequency(cda)+0.5*chanWidth, chanWidth);
2093 31 : Vector<Double> chanWidths(nChan, chanWidth);
2094 31 : if (itsFreqType != MFrequency::TOPO) {
2095 : // have to convert the channel frequencies from topocentric to the specifed
2096 : // frequency type.
2097 3 : itsFreqCtr.setOut(itsFreqType);
2098 3 : Vector<Double> chanValInHz(1);
2099 3 : MVFrequency chanVal;
2100 384 : for (uInt c = 0; c < nChan; c++) {
2101 381 : chanValInHz = chanFreq(c);
2102 381 : chanVal.putVector(chanValInHz);
2103 381 : chanFreq(c) = itsFreqCtr(chanVal).getValue().getValue();
2104 : }
2105 : // To calculate the channel widths I just need to convert the topocentric
2106 : // channel width!
2107 3 : chanValInHz = chanWidth;
2108 3 : chanVal.putVector(chanValInHz);
2109 3 : chanWidths = itsFreqCtr(chanVal).getValue().getValue();
2110 3 : }
2111 31 : spw.chanFreq().put(newRow, chanFreq);
2112 31 : spw.chanWidth().put(newRow, chanWidths);
2113 31 : spw.effectiveBW().put(newRow, chanWidths);
2114 :
2115 31 : spw.flagRow().put(newRow, false);
2116 : {
2117 31 : Quantum<Double> qChanWidth(chanWidth, "Hz");
2118 31 : if (chanWidth < 1E6) {
2119 5 : qChanWidth.convert("kHz");
2120 : } else {
2121 26 : qChanWidth.convert("MHz");
2122 : }
2123 31 : Quantum<Double> qRefFreq(refFreq.get("GHz"));
2124 31 : if (qRefFreq.getValue() < 1) {
2125 0 : qRefFreq.convert("MHz");
2126 : }
2127 31 : ostringstream str;
2128 31 : str << nChan
2129 31 : << "*" << std::setprecision(3) << qChanWidth
2130 31 : << " channels @ " << qRefFreq
2131 31 : << " (" << refFreq.getRefString() << ")";
2132 31 : spw.name().put(newRow, String(str));
2133 31 : }
2134 31 : if (sda.smoothed() || nChan == 1) {
2135 : // the effective resolution is just the channel width
2136 31 : spw.resolution().put(newRow, chanWidths);
2137 : } else {
2138 : // 1.21 is the FWHM of the implicitly convolved sync function
2139 : // (see TMS 2nd ed p. 286)
2140 0 : spw.resolution().put(newRow,
2141 0 : static_cast<Array<Double> >(chanWidths) *
2142 0 : static_cast<Double>(1.21));
2143 : // The static cast is for a possible SGI compiler bug. The compiler seems
2144 : // unable to figure out that chanWidths is an Array
2145 :
2146 : }
2147 :
2148 31 : spw.totalBandwidth().put(newRow, sum(chanWidths));
2149 31 : spw.netSideband().put(newRow, 1);
2150 : //
2151 : static uInt curMSRows = 0;
2152 31 : Int freqGroupId = 1;
2153 31 : if (nrow() == curMSRows) {
2154 18 : if (newRow > 0) {
2155 13 : freqGroupId = spw.freqGroup()(newRow-1);
2156 : }
2157 : } else {
2158 13 : freqGroupId = max(spw.freqGroup().getColumn()) + 1;
2159 : }
2160 31 : spw.freqGroup().put(newRow, freqGroupId);
2161 31 : spw.freqGroupName().put(newRow, "Group " + String::toString(freqGroupId));
2162 31 : curMSRows = nrow();
2163 31 : return newRow;
2164 :
2165 31 : }
2166 :
2167 9 : uInt VLAFiller::addPolarization(const Vector<Stokes::StokesTypes>& polTypes) {
2168 9 : MSPolarizationColumns& pol = polarization();
2169 9 : const uInt newRow = pol.nrow();
2170 9 : itsMS.polarization().addRow(1);
2171 9 : const uInt nCorr = polTypes.nelements();
2172 9 : pol.numCorr().put(newRow, nCorr);
2173 9 : Vector<Int> polInt(nCorr);
2174 9 : Matrix<Int> polProd(2, nCorr);
2175 33 : for (uInt p = 0; p < nCorr; p++) {
2176 24 : polInt(p) = polTypes(p);
2177 24 : switch (polTypes(p)) {
2178 9 : case Stokes::RR:
2179 9 : polProd(0, p) = 0;
2180 9 : polProd(1, p) = 0;
2181 9 : break;
2182 5 : case Stokes::RL:
2183 5 : polProd(0, p) = 0;
2184 5 : polProd(1, p) = 1;
2185 5 : break;
2186 5 : case Stokes::LR:
2187 5 : polProd(0, p) = 1;
2188 5 : polProd(1, p) = 0;
2189 5 : break;
2190 5 : case Stokes::LL:
2191 5 : polProd(0, p) = 1;
2192 5 : polProd(1, p) = 1;
2193 5 : break;
2194 0 : default:
2195 0 : throw(AipsError("VLAFiller::addPolarization - Bad polarization value"));
2196 : }
2197 : }
2198 9 : pol.corrType().put(newRow, polInt);
2199 9 : pol.corrProduct().put(newRow, polProd);
2200 9 : pol.flagRow().put(newRow, false);
2201 9 : return newRow;
2202 9 : }
2203 :
2204 31 : uInt VLAFiller::addDataDescription(uInt spwId, uInt polId) {
2205 31 : MSDataDescColumns& dd = dataDescription();
2206 31 : const uInt newRow = dd.nrow();
2207 31 : itsMS.dataDescription().addRow(1);
2208 31 : dd.spectralWindowId().put(newRow, spwId);
2209 31 : dd.polarizationId().put(newRow, polId);
2210 31 : if (!dd.lagId().isNull()) {
2211 0 : dd.lagId().put(newRow, -1);
2212 : }
2213 31 : return newRow;
2214 : }
2215 :
2216 :
2217 30 : uInt VLAFiller::addSource(const MDirection& dir ){
2218 :
2219 : // TBD: this should be revised to handle multiple restfreq/cda
2220 : // (requires careful coordination with addSpw, addFld, etc.)
2221 :
2222 30 : const VLASDA& sda = itsRecord.SDA();
2223 :
2224 30 : MSSourceColumns& src = source();
2225 30 : const uInt newRow = src.nrow();
2226 :
2227 30 : if (newRow==0) {
2228 : // Set frame info...
2229 9 : src.setFrequencyRef(MFrequency::REST);
2230 :
2231 : // This assumes the MS will have one frame
2232 9 : for (uInt c = 0; c < maxCDA; c++) {
2233 9 : const VLACDA& cda = itsRecord.CDA(c);
2234 :
2235 9 : if (cda.isValid()) {
2236 9 : const VLAEnum::CDA thisCDA = VLAEnum::CDA(c);
2237 : // We subtract 1 because RV frames are one less than Freq frames...
2238 9 : src.setRadialVelocityRef(MRadialVelocity::Types(sda.restFrame(thisCDA)-1));
2239 9 : break;
2240 : }
2241 : }
2242 : }
2243 :
2244 30 : itsMS.source().addRow(1);
2245 30 : src.name().put(newRow, sda.sourceName());
2246 30 : src.sourceId().put( newRow, newRow ); // added by GYL
2247 30 : src.spectralWindowId().put(newRow, -1);
2248 30 : const MEpoch* mepPtr = dynamic_cast<const MEpoch*>(itsFrame.epoch());
2249 30 : src.timeMeas().put(newRow, *mepPtr);
2250 30 : src.code().put(newRow, sda.calCode());
2251 30 : src.directionMeas().put(newRow, dir);
2252 30 : MPosition obsPos;
2253 30 : MeasTable::Observatory(obsPos, "VLA");
2254 30 : MeasFrame frame(*mepPtr, obsPos, dir);
2255 30 : Vector<Double> restFreq;
2256 30 : Vector<Double> sysvel;
2257 30 : uInt validRest=0;
2258 150 : for (uInt c = 0; c < maxCDA; c++) {
2259 120 : const VLACDA& cda = itsRecord.CDA(c);
2260 :
2261 120 : if (cda.isValid()) {
2262 46 : const VLAEnum::CDA thisCDA = VLAEnum::CDA(c);
2263 : // check if it has a valid frame...some old data come without one !
2264 : // Who said bad documentation is better than none ...aargh !
2265 46 : if(sda.restFrame(thisCDA) < MFrequency::N_Types){
2266 46 : ++validRest;
2267 46 : restFreq.resize(validRest, true);
2268 46 : restFreq(validRest-1)=sda.restFrequency(thisCDA);
2269 :
2270 46 : sysvel.resize(validRest,true);
2271 46 : sysvel(validRest-1) = sda.radialVelocity(thisCDA);
2272 :
2273 : }
2274 : }
2275 : }
2276 :
2277 30 : Vector<String> transition(restFreq.size(), "Unknown");
2278 30 : src.numLines().put(newRow, validRest);
2279 30 : src.restFrequency().put(newRow, restFreq);
2280 30 : src.sysvel().put(newRow, sysvel);
2281 30 : src.transition().put(newRow, transition);
2282 :
2283 30 : return newRow;
2284 :
2285 30 : }
2286 :
2287 0 : void VLAFiller::extendHypercubes(const Block<uInt>& nPol,
2288 : const Block<uInt>& nChan, uInt rows) {
2289 0 : const uInt nSpId = nChan.nelements();
2290 0 : DebugAssert(nPol.nelements() == nSpId, AipsError);
2291 0 : for (uInt s = 0; s < nSpId; s++) {
2292 : {
2293 0 : itsTileId.define(sigmaTileId, static_cast<Int>(nPol[s]));
2294 0 : itsSigmaAcc.extendHypercube(rows, itsTileId);
2295 0 : itsTileId.removeField(sigmaTileId);
2296 : }
2297 : {
2298 0 : itsTileId.define(dataTileId, static_cast<Int>(10*nChan[s] + nPol[s]));
2299 0 : itsDataAcc.extendHypercube(rows, itsTileId);
2300 0 : itsTileId.removeField(dataTileId);
2301 : }
2302 : {
2303 0 : itsTileId.define(modDataTileId, static_cast<Int>(10*nChan[s] + nPol[s]));
2304 0 : itsModDataAcc.extendHypercube(rows, itsTileId);
2305 0 : itsTileId.removeField(modDataTileId);
2306 : }
2307 : {
2308 0 : itsTileId.define(corrDataTileId, static_cast<Int>(10*nChan[s] + nPol[s])); itsCorrDataAcc.extendHypercube(rows, itsTileId);
2309 0 : itsTileId.removeField(corrDataTileId);
2310 : }
2311 :
2312 : {
2313 0 : itsTileId.define(chanFlagTileId, static_cast<Int>(10*nChan[s] + nPol[s])); itsChanFlagAcc.extendHypercube(rows, itsTileId);
2314 0 : itsTileId.removeField(chanFlagTileId);
2315 : }
2316 : {
2317 0 : itsTileId.define(flagTileId, static_cast<Int>(10*nChan[s] + nPol[s]));
2318 0 : itsFlagAcc.extendHypercube(rows, itsTileId);
2319 0 : itsTileId.removeField(flagTileId);
2320 : }
2321 : }
2322 0 : }
2323 :
2324 0 : void VLAFiller::addHypercubes(uInt nPol, uInt nChan) {
2325 0 : DebugAssert(nChan > 0 && nChan <= 4096, AipsError);
2326 0 : DebugAssert(nPol > 0 && nPol <= 4, AipsError);
2327 0 : Bool addDataCube = true;
2328 0 : Bool addSigmaCube = true;
2329 0 : uInt s = itsDataShapes.nelements();
2330 0 : while (addDataCube && s > 0) {
2331 0 : s--;
2332 0 : const IPosition& curShape = itsDataShapes[s];
2333 0 : if (curShape(0) == static_cast<Int>(nPol)) {
2334 0 : addSigmaCube = false;
2335 : }
2336 0 : if (curShape(0) == static_cast<Int>(nPol) &&
2337 0 : curShape(1) == static_cast<Int>(nChan)) {
2338 0 : addDataCube = false;
2339 : }
2340 :
2341 0 : DebugAssert(addDataCube || !addSigmaCube, AipsError);
2342 : }
2343 :
2344 0 : if (addDataCube) {
2345 0 : if (addSigmaCube) {
2346 0 : itsTileId.define(sigmaTileId, static_cast<Int>(nPol));
2347 :
2348 : // 1 MB tile size
2349 0 : uInt rowTiles = 131072 / nPol;
2350 0 : if (rowTiles < 1) rowTiles = 1;
2351 0 : itsSigmaAcc.addHypercube(IPosition(2, nPol, 0),
2352 0 : IPosition(2, nPol, rowTiles),
2353 0 : itsTileId);
2354 : }
2355 :
2356 0 : itsTileId.define(dataTileId, static_cast<Int>(10*nChan + nPol));
2357 : // 1 MB tile size
2358 0 : uInt rowTiles = 131072 / nChan / nPol;
2359 0 : if (rowTiles < 1) rowTiles = 1;
2360 0 : itsDataAcc.addHypercube(IPosition(3, nPol, nChan, 0),
2361 0 : IPosition(3, nPol, nChan, rowTiles),
2362 0 : itsTileId);
2363 :
2364 0 : itsTileId.removeField(dataTileId);
2365 :
2366 0 : itsTileId.define(modDataTileId, static_cast<Int>(10*nChan + nPol));
2367 0 : itsModDataAcc.addHypercube(IPosition(3, nPol, nChan, 0),
2368 0 : IPosition(3, nPol, nChan, rowTiles),
2369 0 : itsTileId);
2370 0 : itsTileId.removeField(modDataTileId);
2371 :
2372 0 : itsTileId.define(corrDataTileId, static_cast<Int>(10*nChan + nPol));
2373 0 : itsCorrDataAcc.addHypercube(IPosition(3, nPol, nChan, 0),
2374 0 : IPosition(3, nPol, nChan, rowTiles),
2375 0 : itsTileId);
2376 0 : itsTileId.removeField(corrDataTileId);
2377 0 : itsTileId.define(chanFlagTileId, static_cast<Int>(10*nChan + nPol));
2378 0 : itsChanFlagAcc.addHypercube(IPosition(3, nPol, nChan, 0),
2379 0 : IPosition(3, nPol, nChan, rowTiles),
2380 0 : itsTileId);
2381 0 : itsTileId.removeField(chanFlagTileId);
2382 :
2383 0 : itsTileId.define(flagTileId, static_cast<Int>(10*nChan + nPol));
2384 0 : rowTiles = 131072 / (nPol * nChan * nCat);
2385 0 : if (rowTiles < 1) rowTiles = 1;
2386 0 : itsFlagAcc.addHypercube(IPosition(4, nPol, nChan, nCat, 0),
2387 0 : IPosition(4, nPol, nChan, nCat, rowTiles),
2388 0 : itsTileId);
2389 0 : itsTileId.removeField(flagTileId);
2390 0 : const uInt nCubes = itsDataShapes.nelements();
2391 0 : itsDataShapes.resize(nCubes + 1);
2392 0 : itsDataShapes[nCubes] = IPosition(2, nPol, nChan);
2393 : }
2394 0 : }
2395 :
2396 17952 : Int VLAFiller::polIndexer(Stokes::StokesTypes& stokes){
2397 17952 : if (stokes == Stokes::RR)
2398 0 : return 0;
2399 17952 : else if (stokes == Stokes::RL)
2400 5984 : return 1;
2401 11968 : else if (stokes == Stokes::LR)
2402 5984 : return 2;
2403 5984 : else if (stokes == Stokes::LL)
2404 5984 : return 3;
2405 : else
2406 0 : return -1;
2407 : }
2408 :
2409 11 : void VLAFiller::fixFieldDuplicates(MSField& msFld) {
2410 :
2411 11 : MSFieldIndex MSFldIdx(msFld);
2412 11 : MSFieldColumns fldcol(msFld);
2413 11 : Vector<String> name(fldcol.name().getColumn());
2414 11 : Int nFld=fldcol.nrow();
2415 :
2416 43 : for (Int ifld=0;ifld<nFld;++ifld) {
2417 32 : String thisname=name(ifld);
2418 32 : Vector<Int> nameMatches = MSFldIdx.matchFieldName(name(ifld));
2419 32 : Int nMatch=nameMatches.nelements();
2420 32 : if (nMatch>1) {
2421 0 : Int suffix(0);
2422 : {
2423 0 : ostringstream newname;
2424 0 : newname << thisname << "_" << suffix;
2425 0 : name(nameMatches(0))=String(newname);
2426 0 : fldcol.name().put(nameMatches(0),name(nameMatches(0)));
2427 0 : }
2428 0 : for (Int imatch=1;imatch<nMatch;++imatch) {
2429 0 : suffix++;
2430 0 : ostringstream newname;
2431 0 : newname << thisname << "_" << suffix;
2432 0 : name(nameMatches(imatch))=String(newname);
2433 0 : fldcol.name().put(nameMatches(imatch),name(nameMatches(imatch)));
2434 0 : }
2435 : }
2436 32 : }
2437 :
2438 11 : }
2439 :
2440 9149 : casacore::MDirection::Types VLAFiller::validEpoch(casacore::MDirection::Types mdType)
2441 : {
2442 9149 : if (mdType == MDirection::N_Types) {
2443 : // epoch in data is 0, warn and assume B1950_VLA
2444 0 : mdType = MDirection::B1950_VLA;
2445 0 : if (!itsZeroEpochWarned) {
2446 0 : itsLog << LogIO::WARN
2447 : << "epoch is 0, assuming B1950_VLA"
2448 0 : << LogIO::POST;
2449 0 : itsZeroEpochWarned = true;
2450 : }
2451 : }
2452 9149 : return mdType;
2453 : }
2454 :
2455 :
2456 : // Local Variables:
2457 : // compile-command: "gmake VLAFiller; cd ../../apps/vlafiller; gmake OPTLIB=1"
2458 : // End:
|
