source: branches/alma/external/atnf/PKSIO/PKSMS2writer.cc@ 1500

Last change on this file since 1500 was 1453, checked in by TakTsutsumi, 16 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: many

Test Programs: sd.scantable(), sd.scantable.save()

Put in Release Notes: N/A

Description: copied from current casapy code tree


File size: 37.0 KB
Line 
1//#---------------------------------------------------------------------------
2//# PKSMS2writer.cc: Class to write Parkes multibeam data to a measurementset.
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2000-2006
5//# Associated Universities, Inc. Washington DC, USA.
6//#
7//# This library is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU Library General Public License as published by
9//# the Free Software Foundation; either version 2 of the License, or (at your
10//# option) any later version.
11//#
12//# This library is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14//# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15//# License for more details.
16//#
17//# You should have received a copy of the GNU Library General Public License
18//# along with this library; if not, write to the Free Software Foundation,
19//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning AIPS++ should be addressed as follows:
22//# Internet email: aips2-request@nrao.edu.
23//# Postal address: AIPS++ Project Office
24//# National Radio Astronomy Observatory
25//# 520 Edgemont Road
26//# Charlottesville, VA 22903-2475 USA
27//#
28//# $Id: PKSMS2writer.cc,v 19.11 2006/11/06 22:25:22 mmarquar Exp $
29//#---------------------------------------------------------------------------
30
31#include <atnf/PKSIO/PKSMS2writer.h>
32
33#include <casa/Arrays/ArrayUtil.h>
34#include <casa/Arrays/ArrayMath.h>
35#include <casa/Arrays/ArrayLogical.h>
36#include <casa/Arrays/Cube.h>
37#include <casa/BasicSL/Complex.h>
38#include <casa/BasicSL/Constants.h>
39#include <casa/Quanta/QC.h>
40#include <casa/Logging/LogIO.h>
41#include <measures/Measures/Stokes.h>
42#include <tables/Tables/ArrColDesc.h>
43#include <tables/Tables/IncrementalStMan.h>
44#include <tables/Tables/ScaColDesc.h>
45#include <tables/Tables/SetupNewTab.h>
46#include <tables/Tables/StandardStMan.h>
47#include <tables/Tables/Table.h>
48#include <tables/Tables/TableDesc.h>
49#include <tables/Tables/TiledShapeStMan.h>
50
51//------------------------------------------------- PKSMS2writer::PKSMS2writer
52
53// Default constructor.
54
55PKSMS2writer::PKSMS2writer()
56{
57}
58
59//------------------------------------------------ PKSMS2writer::~PKSMS2writer
60
61// Destructor.
62
63PKSMS2writer::~PKSMS2writer()
64{
65 close();
66}
67
68//------------------------------------------------------- PKSMS2writer::create
69
70// Create the output MS and and write static data.
71
72Int PKSMS2writer::create(
73 const String msName,
74 const String observer,
75 const String project,
76 const String antName,
77 const Vector<Double> antPosition,
78 const String obsMode,
79 const Float equinox,
80 const String dopplerFrame,
81 const Vector<uInt> nChan,
82 const Vector<uInt> nPol,
83 const Vector<Bool> haveXPol,
84 const Bool haveBase,
85 const String fluxUnit)
86{
87 // Open a MS table.
88 TableDesc pksDesc = MS::requiredTableDesc();
89
90 cNChan.assign(nChan);
91 cNPol.assign(nPol);
92 cHaveXPol.assign(haveXPol);
93
94 Int maxNPol = max(cNPol);
95 cGBT = cAPEX = cSMT = cALMA = False;
96
97 // check if it is GBT data
98 if (antName.contains("GBT")) {
99 cGBT = True;
100 }
101 else if (antName.contains("APEX")) {
102 cAPEX = True;
103 }
104 else if (antName.contains("HHT") || antName.contains("SMT")) {
105 cSMT = True;
106 }
107 else if (antName.contains("ALMA")) {
108 cALMA = True;
109 }
110
111
112
113 //cGBT = antName.contains("GBT");
114 //cAPEX = antName.contains("APEX");
115 //cSMT = antName.contains("HHT");
116 //cALMA = antName.contains("ALMA");
117
118 // Add the non-standard CALFCTR column.
119 pksDesc.addColumn(ArrayColumnDesc<Float>("CALFCTR", "Calibration factors",
120 IPosition(1,maxNPol), ColumnDesc::Direct));
121
122 // Add the optional FLOAT_DATA column.
123 MS::addColumnToDesc(pksDesc, MS::FLOAT_DATA, 2);
124 //pksDesc.rwColumnDesc(MS::columnName(MS::FLOAT_DATA)).rwKeywordSet().
125 // define("UNIT", String("Jy"));
126 pksDesc.rwColumnDesc(MS::columnName(MS::FLOAT_DATA)).rwKeywordSet().
127 define("UNIT", fluxUnit);
128 pksDesc.rwColumnDesc(MS::columnName(MS::FLOAT_DATA)).rwKeywordSet().
129 define("MEASURE_TYPE", "");
130
131 if ((cHaveBase = haveBase)) {
132 // Add the non-standard BASELIN and BASESUB columns.
133 pksDesc.addColumn(ArrayColumnDesc<Float>("BASELIN", "Linear baseline fit",
134 IPosition(2,2,maxNPol), ColumnDesc::Direct));
135 pksDesc.addColumn(ArrayColumnDesc<Float>("BASESUB", "Baseline subtracted",
136 IPosition(2,9,maxNPol), ColumnDesc::Direct));
137 }
138
139 // Add the optional DATA column if cross-polarizations are to be recorded.
140 if (ntrue(cHaveXPol)) {
141 // Add the non-standard XCALFCTR column.
142 pksDesc.addColumn(ScalarColumnDesc<Complex>("XCALFCTR",
143 "Cross-polarization calibration factor"));
144
145 MS::addColumnToDesc(pksDesc, MS::DATA, 2);
146 //pksDesc.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet().
147 // define("UNIT", "Jy");
148 pksDesc.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet().
149 define("UNIT", fluxUnit);
150 pksDesc.rwColumnDesc(MS::columnName(MS::DATA)).rwKeywordSet().
151 define("MEASURE_TYPE", "");
152 }
153
154 // Define hypercube for the float data (without coordinates).
155 pksDesc.defineHypercolumn("TiledData", 3,
156 stringToVector(MS::columnName(MS::FLOAT_DATA)));
157
158 SetupNewTable newtab(msName, pksDesc, Table::New);
159
160 // Set Incremental Storage Manager as the default.
161 IncrementalStMan incrStMan("ISMData");
162 newtab.bindAll(incrStMan, True);
163
164 // Use TiledShapeStMan for the FLOAT_DATA hypercube with tile size 16 kiB.
165 TiledShapeStMan tiledStMan("TiledData", IPosition(3,1,128,32));
166 newtab.bindColumn(MS::columnName(MS::FLOAT_DATA), tiledStMan);
167
168 // Use Standard Storage Manager to handle columns that change for each row.
169 StandardStMan stdStMan;
170 newtab.bindColumn(MS::columnName(MS::SCAN_NUMBER), stdStMan);
171 newtab.bindColumn(MS::columnName(MS::TIME), stdStMan);
172 newtab.bindColumn(MS::columnName(MS::SIGMA), stdStMan);
173 if (maxNPol > 2) {
174 newtab.bindColumn(MS::columnName(MS::DATA), stdStMan);
175 }
176
177 // Create the measurementset.
178 cPKSMS = new MeasurementSet(newtab, 0);
179
180 // Create subtables.
181 TableDesc antennaDesc = MSAntenna::requiredTableDesc();
182 SetupNewTable antennaSetup(cPKSMS->antennaTableName(), antennaDesc,
183 Table::New);
184 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::ANTENNA),
185 Table(antennaSetup));
186
187 TableDesc dataDescDesc = MSDataDescription::requiredTableDesc();
188 SetupNewTable dataDescSetup(cPKSMS->dataDescriptionTableName(), dataDescDesc,
189 Table::New);
190 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::DATA_DESCRIPTION),
191 Table(dataDescSetup));
192
193 TableDesc dopplerDesc = MSDoppler::requiredTableDesc();
194 SetupNewTable dopplerSetup(cPKSMS->dopplerTableName(), dopplerDesc,
195 Table::New);
196 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::DOPPLER),
197 Table(dopplerSetup));
198
199 TableDesc feedDesc = MSFeed::requiredTableDesc();
200 MSFeed::addColumnToDesc(feedDesc, MSFeedEnums::FOCUS_LENGTH);
201 SetupNewTable feedSetup(cPKSMS->feedTableName(), feedDesc, Table::New);
202 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::FEED),
203 Table(feedSetup));
204
205 TableDesc fieldDesc = MSField::requiredTableDesc();
206 SetupNewTable fieldSetup(cPKSMS->fieldTableName(), fieldDesc, Table::New);
207 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::FIELD),
208 Table(fieldSetup));
209
210 TableDesc flagCmdDesc = MSFlagCmd::requiredTableDesc();
211 SetupNewTable flagCmdSetup(cPKSMS->flagCmdTableName(), flagCmdDesc,
212 Table::New);
213 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::FLAG_CMD),
214 Table(flagCmdSetup));
215
216 TableDesc historyDesc = MSHistory::requiredTableDesc();
217 SetupNewTable historySetup(cPKSMS->historyTableName(), historyDesc,
218 Table::New);
219 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::HISTORY),
220 Table(historySetup));
221
222 TableDesc observationDesc = MSObservation::requiredTableDesc();
223 SetupNewTable observationSetup(cPKSMS->observationTableName(),
224 observationDesc, Table::New);
225 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::OBSERVATION),
226 Table(observationSetup));
227
228 TableDesc pointingDesc = MSPointing::requiredTableDesc();
229 SetupNewTable pointingSetup(cPKSMS->pointingTableName(), pointingDesc,
230 Table::New);
231 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::POINTING),
232 Table(pointingSetup));
233
234 TableDesc polarizationDesc = MSPolarization::requiredTableDesc();
235 SetupNewTable polarizationSetup(cPKSMS->polarizationTableName(),
236 polarizationDesc, Table::New);
237 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::POLARIZATION),
238 Table(polarizationSetup));
239
240 TableDesc processorDesc = MSProcessor::requiredTableDesc();
241 SetupNewTable processorSetup(cPKSMS->processorTableName(), processorDesc,
242 Table::New);
243 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::PROCESSOR),
244 Table(processorSetup));
245
246 TableDesc sourceDesc = MSSource::requiredTableDesc();
247 MSSource::addColumnToDesc(sourceDesc, MSSourceEnums::TRANSITION, 1);
248 MSSource::addColumnToDesc(sourceDesc, MSSourceEnums::REST_FREQUENCY,
249 1);
250 MSSource::addColumnToDesc(sourceDesc, MSSourceEnums::SYSVEL, 1);
251 SetupNewTable sourceSetup(cPKSMS->sourceTableName(), sourceDesc,
252 Table::New);
253 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::SOURCE),
254 Table(sourceSetup));
255
256 TableDesc spectralWindowDesc = MSSpectralWindow::requiredTableDesc();
257 MSSpectralWindow::addColumnToDesc(spectralWindowDesc,
258 MSSpectralWindowEnums::DOPPLER_ID);
259 SetupNewTable spectralWindowSetup(cPKSMS->spectralWindowTableName(),
260 spectralWindowDesc, Table::New);
261 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::SPECTRAL_WINDOW),
262 Table(spectralWindowSetup));
263
264 TableDesc stateDesc = MSState::requiredTableDesc();
265 SetupNewTable stateSetup(cPKSMS->stateTableName(), stateDesc, Table::New);
266 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::STATE),
267 Table(stateSetup));
268
269 TableDesc sysCalDesc = MSSysCal::requiredTableDesc();
270 MSSysCal::addColumnToDesc(sysCalDesc, MSSysCalEnums::TCAL, 1);
271 MSSysCal::addColumnToDesc(sysCalDesc, MSSysCalEnums::TSYS, 1);
272 SetupNewTable sysCalSetup(cPKSMS->sysCalTableName(), sysCalDesc,
273 Table::New);
274 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::SYSCAL),
275 Table(sysCalSetup));
276
277 TableDesc weatherDesc = MSWeather::requiredTableDesc();
278 MSWeather::addColumnToDesc(weatherDesc, MSWeatherEnums::PRESSURE);
279 MSWeather::addColumnToDesc(weatherDesc, MSWeatherEnums::REL_HUMIDITY);
280 MSWeather::addColumnToDesc(weatherDesc, MSWeatherEnums::TEMPERATURE);
281 SetupNewTable weatherSetup(cPKSMS->weatherTableName(), weatherDesc,
282 Table::New);
283 cPKSMS->rwKeywordSet().defineTable(MS::keywordName(MS::WEATHER),
284 Table(weatherSetup));
285
286 cPKSMS->initRefs();
287
288 // Measurementset subtables.
289 cAntenna = cPKSMS->antenna();
290 cDataDescription = cPKSMS->dataDescription();
291 cDoppler = cPKSMS->doppler();
292 cFeed = cPKSMS->feed();
293 cField = cPKSMS->field();
294 cFlagCmd = cPKSMS->flagCmd();
295 cHistory = cPKSMS->history();
296 cObservation = cPKSMS->observation();
297 cPointing = cPKSMS->pointing();
298 cPolarization = cPKSMS->polarization();
299 cProcessor = cPKSMS->processor();
300 cSource = cPKSMS->source();
301 cSpectralWindow = cPKSMS->spectralWindow();
302 cState = cPKSMS->state();
303 cSysCal = cPKSMS->sysCal();
304 cWeather = cPKSMS->weather();
305
306 // Measurementset table columns;
307 cMSCols = new MSColumns(*cPKSMS);
308 cAntennaCols = new MSAntennaColumns(cAntenna);
309 cDataDescCols = new MSDataDescColumns(cDataDescription);
310 cDopplerCols = new MSDopplerColumns(cDoppler);
311 cFeedCols = new MSFeedColumns(cFeed);
312 cFieldCols = new MSFieldColumns(cField);
313 cFlagCmdCols = new MSFlagCmdColumns(cFlagCmd);
314 cHistoryCols = new MSHistoryColumns(cHistory);
315 cObservationCols = new MSObservationColumns(cObservation);
316 cPointingCols = new MSPointingColumns(cPointing);
317 cPolarizationCols = new MSPolarizationColumns(cPolarization);
318 cProcessorCols = new MSProcessorColumns(cProcessor);
319 cSourceCols = new MSSourceColumns(cSource);
320 cSpWindowCols = new MSSpWindowColumns(cSpectralWindow);
321 cStateCols = new MSStateColumns(cState);
322 cSysCalCols = new MSSysCalColumns(cSysCal);
323 cWeatherCols = new MSWeatherColumns(cWeather);
324
325 cCalFctrCol = new ArrayColumn<Float>(*cPKSMS, "CALFCTR");
326 if (cHaveBase) {
327 cBaseLinCol = new ArrayColumn<Float>(*cPKSMS, "BASELIN");
328 cBaseSubCol = new ArrayColumn<Float>(*cPKSMS, "BASESUB");
329 }
330 if (ntrue(cHaveXPol)) {
331 cXCalFctrCol = new ScalarColumn<Complex>(*cPKSMS, "XCALFCTR");
332 }
333
334
335 // Define Measure references.
336 Vector<String> flagCat(1, "BAD");
337 cMSCols->sigma().rwKeywordSet().define("UNIT", "K");
338 cMSCols->flagCategory().rwKeywordSet().define("CATEGORY", flagCat);
339
340 String dirref;
341 if (equinox == 1950.0f) {
342 dirref = "B1950";
343 } else {
344 dirref = "J2000";
345 }
346
347 cFieldCols->delayDir().rwKeywordSet().asrwRecord("MEASINFO").
348 define("Ref", dirref);
349 cFieldCols->phaseDir().rwKeywordSet().asrwRecord("MEASINFO").
350 define("Ref", dirref);
351 cFieldCols->referenceDir().rwKeywordSet().asrwRecord("MEASINFO").
352 define("Ref", dirref);
353
354 cPointingCols->direction().rwKeywordSet().asrwRecord("MEASINFO").
355 define("Ref", dirref);
356 cPointingCols->target().rwKeywordSet().asrwRecord("MEASINFO").
357 define("Ref", dirref);
358
359 cSourceCols->direction().rwKeywordSet().asrwRecord("MEASINFO").
360 define("Ref", dirref);
361 cSourceCols->restFrequency().rwKeywordSet().asrwRecord("MEASINFO").
362 define("Ref", "REST");
363
364 // Translate Doppler frame name.
365 if (dopplerFrame == "TOPOCENT") {
366 MFrequency::getType(cDopplerFrame, "TOPO");
367 } else if (dopplerFrame == "GEOCENTR") {
368 MFrequency::getType(cDopplerFrame, "GEO");
369 } else if (dopplerFrame == "BARYCENT") {
370 MFrequency::getType(cDopplerFrame, "BARY");
371 } else if (dopplerFrame == "GALACTOC") {
372 MFrequency::getType(cDopplerFrame, "GALACTO");
373 } else if (dopplerFrame == "LOCALGRP") {
374 MFrequency::getType(cDopplerFrame, "LGROUP");
375 } else if (dopplerFrame == "CMBDIPOL") {
376 MFrequency::getType(cDopplerFrame, "CMB");
377 } else if (dopplerFrame == "SOURCE") {
378 MFrequency::getType(cDopplerFrame, "REST");
379 } else if (dopplerFrame == "LSRK") {
380 MFrequency::getType(cDopplerFrame, "LSRK");
381 }
382
383
384 // Store static data.
385 addAntennaEntry(antName, antPosition);
386 addDopplerEntry();
387 addFeedEntry();
388 //addObservationEntry(observer, project);
389 addObservationEntry(observer, project, antName);
390 addProcessorEntry();
391
392 return 0;
393}
394
395//-------------------------------------------------------- PKSMS2writer::write
396
397// Write the next data record.
398
399Int PKSMS2writer::write(
400 const Int scanNo,
401 const Int cycleNo,
402 const Double mjd,
403 const Double interval,
404 const String fieldName,
405 const String srcName,
406 const Vector<Double> srcDir,
407 const Vector<Double> srcPM,
408 const Double srcVel,
409 const String obsMode,
410 const Int IFno,
411 const Double refFreq,
412 const Double bandwidth,
413 const Double freqInc,
414 //const Double restFreq,
415 const Vector<Double> restFreq,
416 const Vector<Float> tcal,
417 const String tcalTime,
418 const Float azimuth,
419 const Float elevation,
420 const Float parAngle,
421 const Float focusAxi,
422 const Float focusTan,
423 const Float focusRot,
424 const Float temperature,
425 const Float pressure,
426 const Float humidity,
427 const Float windSpeed,
428 const Float windAz,
429 const Int refBeam,
430 const Int beamNo,
431 const Vector<Double> direction,
432 const Vector<Double> scanRate,
433 const Vector<Float> tsys,
434 const Vector<Float> sigma,
435 const Vector<Float> calFctr,
436 const Matrix<Float> baseLin,
437 const Matrix<Float> baseSub,
438 const Matrix<Float> &spectra,
439 const Matrix<uChar> &flagged,
440 const Complex xCalFctr,
441 const Vector<Complex> &xPol)
442{
443 // Extend the time range in the OBSERVATION subtable.
444 Vector<Double> timerange(2);
445 cObservationCols->timeRange().get(0, timerange);
446 Double time = mjd*86400.0;
447 if (timerange(0) == 0.0) {
448 timerange(0) = time;
449 }
450 timerange(1) = time;
451 cObservationCols->timeRange().put(0, timerange);
452
453 Int iIF = IFno - 1;
454 Int nChan = cNChan(iIF);
455 Int nPol = cNPol(iIF);
456
457 // IFno is the 1-relative row number in the DATA_DESCRIPTION,
458 // SPECTRAL_WINDOW, and POLARIZATION subtables.
459 if (Int(cDataDescription.nrow()) < IFno) {
460 // Add a new entry to each subtable.
461 addDataDescriptionEntry(IFno);
462 addSpectralWindowEntry(IFno, nChan, refFreq, bandwidth, freqInc);
463 addPolarizationEntry(IFno, nPol);
464 }
465
466 // Find or add the source to the SOURCE subtable.
467 Int srcId = addSourceEntry(srcName, srcDir, srcPM, restFreq, srcVel);
468
469 // Find or add the obsMode to the STATE subtable.
470 Int stateId = addStateEntry(obsMode);
471
472 // FIELD subtable.
473 Int fieldId = addFieldEntry(fieldName, time, direction, scanRate, srcId);
474
475 // POINTING subtable.
476 addPointingEntry(time, interval, fieldName, direction, scanRate);
477
478 // SYSCAL subtable.
479 addSysCalEntry(beamNo, iIF, time, interval, tcal, tsys, nPol);
480
481
482 // Handle weather information.
483 ROScalarColumn<Double> wTime(cWeather, "TIME");
484 Int nWeather = wTime.nrow();
485 if (nWeather == 0 || time > wTime(nWeather-1)) {
486 addWeatherEntry(time, interval, pressure, humidity, temperature);
487 }
488
489
490 // Extend the main table.
491 cPKSMS->addRow();
492 Int irow = cPKSMS->nrow() - 1;
493
494 // Keys.
495 cMSCols->time().put(irow, time);
496 cMSCols->antenna1().put(irow, 0);
497 cMSCols->antenna2().put(irow, 0);
498 cMSCols->feed1().put(irow, beamNo-1);
499 cMSCols->feed2().put(irow, beamNo-1);
500 cMSCols->dataDescId().put(irow, iIF);
501 cMSCols->processorId().put(irow, 0);
502 cMSCols->fieldId().put(irow, fieldId);
503
504 // Non-key attributes.
505 cMSCols->interval().put(irow, interval);
506 cMSCols->exposure().put(irow, interval);
507 cMSCols->timeCentroid().put(irow, time);
508 cMSCols->scanNumber().put(irow, scanNo);
509 cMSCols->arrayId().put(irow, 0);
510 cMSCols->observationId().put(irow, 0);
511 cMSCols->stateId().put(irow, stateId);
512
513 Vector<Double> uvw(3, 0.0);
514 cMSCols->uvw().put(irow, uvw);
515
516 // Baseline fit parameters.
517 if (cHaveBase) {
518 cBaseLinCol->put(irow, baseLin);
519
520 if (baseSub.nrow() == 9) {
521 cBaseSubCol->put(irow, baseSub);
522
523 } else {
524 Matrix<Float> tmp(9, 2, 0.0f);
525 for (Int ipol = 0; ipol < nPol; ipol++) {
526 for (uInt j = 0; j < baseSub.nrow(); j++) {
527 tmp(j,ipol) = baseSub(j,ipol);
528 }
529 }
530 cBaseSubCol->put(irow, tmp);
531 }
532 }
533 // Transpose spectra.
534 Matrix<Float> tmpData(nPol, nChan);
535 Matrix<Bool> tmpFlag(nPol, nChan);
536 for (Int ipol = 0; ipol < nPol; ipol++) {
537 for (Int ichan = 0; ichan < nChan; ichan++) {
538 tmpData(ipol,ichan) = spectra(ichan,ipol);
539 tmpFlag(ipol,ichan) = flagged(ichan,ipol);
540 }
541 }
542 cCalFctrCol->put(irow, calFctr);
543 cMSCols->floatData().put(irow, tmpData);
544 cMSCols->flag().put(irow, tmpFlag);
545
546 // Cross-polarization spectra.
547 if (cHaveXPol(iIF)) {
548 cXCalFctrCol->put(irow, xCalFctr);
549 cMSCols->data().put(irow, xPol);
550 }
551
552 cMSCols->sigma().put(irow, sigma);
553
554 //Vector<Float> weight(1, 1.0f);
555 Vector<Float> weight(nPol, 1.0f);
556 cMSCols->weight().put(irow, weight);
557 //imaging weight
558 //Vector<Float> imagingWeight(nChan);
559 //cMSCols->imagingWeight().put(irow, imagingWeight);
560
561 // Flag information.
562 Cube<Bool> flags(nPol, nChan, 1, False);
563 //cMSCols->flag().put(irow, flags.xyPlane(0));
564 cMSCols->flagCategory().put(irow, flags);
565 cMSCols->flagRow().put(irow, False);
566
567 return 0;
568}
569
570//-------------------------------------------------------- PKSMS2writer::close
571
572// Close the measurementset, flushing all associated tables.
573
574void PKSMS2writer::close()
575{
576 // Delete table column accessors.
577 delete cMSCols; cMSCols=0;
578 delete cAntennaCols; cAntennaCols=0;
579 delete cDataDescCols; cDataDescCols=0;
580 delete cDopplerCols; cDopplerCols=0;
581 delete cFeedCols; cFeedCols=0;
582 delete cFieldCols; cFieldCols=0;
583 delete cFlagCmdCols; cFlagCmdCols=0;
584 delete cHistoryCols; cHistoryCols=0;
585 delete cObservationCols; cObservationCols=0;
586 delete cPointingCols; cPointingCols=0;
587 delete cPolarizationCols; cPolarizationCols=0;
588 delete cProcessorCols; cProcessorCols=0;
589 delete cSourceCols; cSourceCols=0;
590 delete cSpWindowCols; cSpWindowCols=0;
591 delete cStateCols; cStateCols=0;
592 delete cSysCalCols; cSysCalCols=0;
593 delete cWeatherCols; cWeatherCols=0;
594
595 delete cCalFctrCol; cCalFctrCol=0;
596 if (cHaveBase) {
597 delete cBaseLinCol; cBaseLinCol=0;
598 delete cBaseSubCol; cBaseSubCol=0;
599 }
600 if (ntrue(cHaveXPol)) {
601 delete cXCalFctrCol; cXCalFctrCol=0;
602 }
603
604 // Release all subtables.
605 cAntenna = MSAntenna();
606 cDataDescription = MSDataDescription();
607 cDoppler = MSDoppler();
608 cFeed = MSFeed();
609 cField = MSField();
610 cFlagCmd = MSFlagCmd();
611 cHistory = MSHistory();
612 cObservation = MSObservation();
613 cPointing = MSPointing();
614 cPolarization = MSPolarization();
615 cProcessor = MSProcessor();
616 cSource = MSSource();
617 cSpectralWindow = MSSpectralWindow();
618 cState = MSState();
619 cSysCal = MSSysCal();
620 cWeather = MSWeather();
621 // Release the main table.
622 delete cPKSMS; cPKSMS=0;
623}
624
625//---------------------------------------------- PKSMS2writer::addAntennaEntry
626
627// Add an entry to the ANTENNA subtable.
628
629Int PKSMS2writer::addAntennaEntry(
630 const String antName,
631 const Vector<Double> &antPosition)
632{
633 // Extend the ANTENNA subtable.
634 cAntenna.addRow();
635 Int n = cAntenna.nrow() - 1;
636
637 // do specific things for GBT
638 // Data.
639 // plus some more telescopes
640 cAntennaCols->name().put(n, antName);
641 //cAntennaCols->station().put(n, "ATNF_PARKES");
642 if (cGBT) {
643 cAntennaCols->station().put(n, "GREENBANK");
644 cAntennaCols->dishDiameter().put(n, 110.0);
645 }
646 else if (cAPEX) {
647 cAntennaCols->station().put(n, "CHAJNANTOR");
648 cAntennaCols->dishDiameter().put(n, 12.0);
649 }
650 else if (cALMA) {
651 cAntennaCols->station().put(n, "CHAJNANTOR");
652 cAntennaCols->dishDiameter().put(n, 12.0);
653 }
654 else if (cSMT) {
655 cAntennaCols->station().put(n, "MT_GRAHAM");
656 cAntennaCols->dishDiameter().put(n, 10.0);
657 }
658 else {
659 cAntennaCols->station().put(n, "ATNF_PARKES");
660 cAntennaCols->dishDiameter().put(n, 64.0);
661 }
662 cAntennaCols->type().put(n, "GROUND-BASED");
663 cAntennaCols->mount().put(n, "ALT-AZ");
664 cAntennaCols->position().put(n, antPosition);
665 Vector<Double> antOffset(3, 0.0);
666 cAntennaCols->offset().put(n, antOffset);
667 //cAntennaCols->dishDiameter().put(n, 64.0);
668 //if (cGBT) {
669 // cAntennaCols->dishDiameter().put(n, 110.0);
670 //}
671 //else {
672 // cAntennaCols->dishDiameter().put(n, 64.0);
673 //}
674 // Flags.
675 cAntennaCols->flagRow().put(n, False);
676
677 return n;
678}
679
680//-------------------------------------- PKSMS2writer::addDataDescriptionEntry
681
682// Add an entry to the DATA_DESCRIPTION subtable.
683
684Int PKSMS2writer::addDataDescriptionEntry(
685 const Int IFno)
686{
687 // Extend the DATA_DESCRIPTION subtable.
688 while (Int(cDataDescription.nrow()) < IFno) {
689 cDataDescription.addRow();
690 }
691 Int n = IFno - 1;
692
693 // Data.
694 cDataDescCols->spectralWindowId().put(n, n);
695 cDataDescCols->polarizationId().put(n, n);
696
697 // Flags.
698 cDataDescCols->flagRow().put(n, False);
699
700 return n;
701}
702
703//---------------------------------------------- PKSMS2writer::addDopplerEntry
704
705// Add an entry to the DOPPLER subtable.
706
707Int PKSMS2writer::addDopplerEntry()
708{
709 // Extend the DOPPLER subtable.
710 cDoppler.addRow();
711 Int n = cDoppler.nrow() - 1;
712
713 // Keys.
714 cDopplerCols->dopplerId().put(n, n);
715 cDopplerCols->sourceId().put(n, 0);
716
717 // Data.
718 cDopplerCols->transitionId().put(n, 0);
719
720 return n;
721}
722
723//------------------------------------------------- PKSMS2writer::addFeedEntry
724
725// Add an entry to the FEED subtable.
726
727Int PKSMS2writer::addFeedEntry()
728{
729 Int n = cFeed.nrow() - 1;
730 for (Int iBeam = 0; iBeam < 13; iBeam++) {
731 // Extend the FEED subtable.
732 cFeed.addRow();
733 n++;
734
735 // Keys.
736 cFeedCols->antennaId().put(n, cAntenna.nrow()-1);
737 cFeedCols->feedId().put(n, iBeam);
738 cFeedCols->spectralWindowId().put(n, -1);
739 cFeedCols->time().put(n, 0.0);
740 cFeedCols->interval().put(n, -1.0);
741
742 // Data description.
743 cFeedCols->numReceptors().put(n, 2);
744
745 // Data.
746 cFeedCols->beamId().put(n, -1);
747
748 Matrix<Double> beamOffset(2, 2, 0.0);
749 cFeedCols->beamOffset().put(n, beamOffset);
750
751 cFeedCols->focusLength().put(n, 26.0);
752
753 Vector<String> polarizationType(2);
754 polarizationType(0) = "X";
755 polarizationType(1) = "Y";
756 cFeedCols->polarizationType().put(n, polarizationType);
757
758 Matrix<Complex> polResponse(2, 2, Complex(0.0));
759 for (Int i = 0; i < 2; i++) {
760 polResponse(i,i) = Complex(1.0, 0.0);
761 }
762 cFeedCols->polResponse().put(n, polResponse);
763
764 Vector<Double> position(3, 0.0);
765 cFeedCols->position().put(n, position);
766
767 Vector<Double> receptorAngle(2, C::pi_4);
768 receptorAngle(1) += C::pi_2;
769 cFeedCols->receptorAngle().put(n, receptorAngle);
770 }
771
772 return n;
773}
774
775//------------------------------------------------ PKSMS2writer::addFieldEntry
776
777// Add an entry to the FIELD subtable.
778
779Int PKSMS2writer::addFieldEntry(
780 const String fieldName,
781 const Double time,
782 const Vector<Double> direction,
783 const Vector<Double> scanRate,
784 const Int srcId)
785{
786
787 ROScalarColumn<String> fldn(cField, "NAME");
788 ROScalarColumn<Int> sourceid(cField, "SOURCE_ID");
789 Int n;
790 Int nFld = cField.nrow();
791 for (n = 0; n < nFld; n++) {
792 if (fldn(n) == fieldName && sourceid(n) == srcId) {
793 break;
794 }
795 }
796
797 // Extend the FIELD subtable.
798 if (n == nFld) {
799 cField.addRow();
800 //Int n = cField.nrow() - 1;
801
802 // Data.
803 cFieldCols->name().put(n, fieldName);
804 if (cGBT) {
805 cFieldCols->code().put(n, " ");
806 }
807 else {
808 cFieldCols->code().put(n, "DRIFT");
809 }
810 cFieldCols->time().put(n, time);
811
812 //Matrix<Double> track(2, 2);
813 Matrix<Double> track(2, 1);
814 track.column(0) = direction;
815 //track.column(1) = scanRate;
816 cFieldCols->numPoly().put(n, 1);
817 cFieldCols->delayDir().put(n, track);
818 cFieldCols->phaseDir().put(n, track);
819 cFieldCols->referenceDir().put(n, track);
820 cFieldCols->sourceId().put(n, srcId);
821
822 // Flags.
823 cFieldCols->flagRow().put(n, False);
824 }
825
826 return n;
827}
828
829//------------------------------------------ PKSMS2writer::addObservationEntry
830
831// Add an entry to the OBSERVATION subtable.
832
833Int PKSMS2writer::addObservationEntry(
834 const String observer,
835 const String project,
836 const String antName)
837{
838 // Extend the OBSERVATION subtable.
839 cObservation.addRow();
840 Int n = cObservation.nrow() - 1;
841
842 // Data.
843 //cObservationCols->telescopeName().put(n, "Parkes");
844 cObservationCols->telescopeName().put(n, antName);
845 Vector<Double> timerange(2, 0.0);
846 cObservationCols->timeRange().put(n, timerange);
847 cObservationCols->observer().put(n, observer);
848 Vector<String> log(1, "none");
849 cObservationCols->log().put(n, log);
850 //cObservationCols->scheduleType().put(n, "ATNF");
851 cObservationCols->scheduleType().put(n, "");
852 Vector<String> schedule(1, "Not available");
853 cObservationCols->schedule().put(n, schedule);
854 cObservationCols->project().put(n, project);
855 cObservationCols->releaseDate().put(n, 0.0);
856
857 // Flags.
858 cObservationCols->flagRow().put(n, False);
859
860 return n;
861}
862
863//--------------------------------------------- PKSMS2writer::addPointingEntry
864
865// Modified to fill pointing data if the direction is the pointing direction.
866// So the following comment is no longer true.
867
868// Add an entry to the POINTING subtable. This compulsory subtable simply
869// duplicates information in the FIELD subtable.
870
871Int PKSMS2writer::addPointingEntry(
872 const Double time,
873 const Double interval,
874 const String fieldName,
875 const Vector<Double> direction,
876 const Vector<Double> scanRate)
877{
878
879 ROScalarColumn<Double> tms(cPointing, "TIME");
880 Int n;
881 Int ntm = cPointing.nrow();
882 for (n = 0; n < ntm; n++) {
883 if (tms(n) == time) {
884 break;
885 }
886 }
887
888 if (n == ntm) {
889 // Extend the POINTING subtable.
890 cPointing.addRow();
891 //Int n = cPointing.nrow() - 1;
892
893 // Keys.
894 cPointingCols->antennaId().put(n, 0);
895 cPointingCols->time().put(n, time);
896 cPointingCols->interval().put(n, interval);
897
898 // Data.
899 cPointingCols->name().put(n, fieldName);
900 cPointingCols->numPoly().put(n, 1);
901 cPointingCols->timeOrigin().put(n, time);
902
903 //Matrix<Double> track(2, 2);
904 Matrix<Double> track(2, 1);
905 track.column(0) = direction;
906 //track.column(1) = scanRate;
907 cPointingCols->direction().put(n, track);
908 cPointingCols->target().put(n, track);
909 cPointingCols->tracking().put(n, True);
910 }
911 return n;
912}
913
914//----------------------------------------- PKSMS2writer::addPolarizationEntry
915
916// Add an entry to the POLARIZATION subtable.
917
918Int PKSMS2writer::addPolarizationEntry(
919 const Int IFno,
920 const Int nPol)
921{
922 // Extend the POLARIZATION subtable.
923 while (Int(cPolarization.nrow()) < IFno) {
924 cPolarization.addRow();
925 }
926 Int n = IFno - 1;
927
928 // Data description.
929 cPolarizationCols->numCorr().put(n, nPol);
930
931 // Data.
932 Vector<Int> corrType(2);
933 if (nPol == 1) {
934 corrType.resize(1);
935 corrType(0) = Stokes::XX;
936 }
937 else {
938 //Vector<Int> corrType(2);
939 corrType(0) = Stokes::XX;
940 corrType(1) = Stokes::YY;
941 }
942 cPolarizationCols->corrType().put(n, corrType);
943
944 Matrix<Int> corrProduct(2,2,1);
945 if (nPol == 1) {
946 corrProduct.resize(2,1,1);
947 corrProduct(1,0) = 0;
948 }
949 if (nPol == 2) {
950 corrProduct(1,0) = 0;
951 corrProduct(0,1) = 0;
952 }
953 cPolarizationCols->corrProduct().put(n, corrProduct);
954
955 // Flags.
956 cPolarizationCols->flagRow().put(n, False);
957
958 return n;
959}
960
961
962//-------------------------------------------- PKSMS2writer::addProcessorEntry
963
964// Add an entry to the PROCESSOR subtable.
965
966Int PKSMS2writer::addProcessorEntry()
967{
968 // Extend the PROCESSOR subtable.
969 cProcessor.addRow();
970 Int n = cProcessor.nrow() - 1;
971
972 // Data.
973 cProcessorCols->type().put(n, "SPECTROMETER");
974 cProcessorCols->subType().put(n, "MULTIBEAM");
975 cProcessorCols->typeId().put(n, -1);
976 cProcessorCols->modeId().put(n, -1);
977
978 // Flags.
979 cProcessorCols->flagRow().put(n, False);
980
981 return n;
982}
983
984//----------------------------------------------- PKSMS2writer::addSourceEntry
985
986// Add an entry to the SOURCE subtable.
987
988Int PKSMS2writer::addSourceEntry(
989 const String name,
990 const Vector<Double> direction,
991 const Vector<Double> properMotion,
992 //const Double restFreq,
993 const Vector<Double> restFreq,
994 const Double radialVelocity)
995{
996 // Look for an entry in the SOURCE subtable.
997 ROScalarColumn<String> sources(cSource, "NAME");
998 Int n;
999 Int nSrc = sources.nrow();
1000 for (n = 0; n < nSrc; n++) {
1001 if (sources(n) == name) {
1002 break;
1003 }
1004 }
1005
1006 if (n == nSrc) {
1007 // Not found, add a new entry to the SOURCE subtable.
1008 cSource.addRow();
1009
1010 // Keys.
1011 cSourceCols->sourceId().put(n, n);
1012 cSourceCols->time().put(n, 0.0);
1013 cSourceCols->interval().put(n, -1.0);
1014 cSourceCols->spectralWindowId().put(n, -1);
1015
1016 // Data description.
1017 cSourceCols->numLines().put(n, 1);
1018
1019 // Data.
1020 cSourceCols->name().put(n, name);
1021 cSourceCols->calibrationGroup().put(n, 0);
1022 cSourceCols->code().put(n, "");
1023 cSourceCols->direction().put(n, direction);
1024// Vector<Double> position(3, 0.0);
1025// cSourceCols->position().put(n, position);
1026 cSourceCols->properMotion().put(n, properMotion);
1027// Vector<Double> restFrequency(1, restFreq);
1028// cSourceCols->restFrequency().put(n, restFrequency);
1029 cSourceCols->restFrequency().put(n, restFreq);
1030 Vector<Double> sysvel(1, radialVelocity);
1031 cSourceCols->sysvel().put(n, sysvel);
1032 }
1033
1034 return n;
1035}
1036
1037//--------------------------------------- PKSMS2writer::addSpectralWindowEntry
1038
1039// Add an entry to the SPECTRAL_WINDOW subtable.
1040
1041Int PKSMS2writer::addSpectralWindowEntry(
1042 const Int IFno,
1043 const Int nChan,
1044 const Double refFreq,
1045 const Double bandwidth,
1046 const Double freqInc)
1047{
1048 // Extend the SPECTRAL_WINDOW subtable.
1049 while (Int(cSpectralWindow.nrow()) < IFno) {
1050 cSpectralWindow.addRow();
1051 }
1052 Int n = IFno - 1;
1053
1054 // Data description.
1055 cSpWindowCols->numChan().put(n, nChan);
1056
1057 // Data.
1058 //cSpWindowCols->name().put(n, "L-band");
1059 cSpWindowCols->name().put(n, " ");
1060 cSpWindowCols->refFrequency().put(n, refFreq);
1061
1062 // 0-relative reference channel number.
1063 Double refChan = nChan / 2;
1064 Vector<Double> freqs(nChan);
1065 for (Int i = 0; i < nChan; i++) {
1066 freqs(i) = refFreq + (i - refChan)*freqInc;
1067 }
1068 cSpWindowCols->chanFreq().put(n, freqs);
1069
1070 Vector<Double> chanWidths(nChan, freqInc);
1071 cSpWindowCols->chanWidth().put(n, chanWidths);
1072
1073 cSpWindowCols->measFreqRef().put(n, cDopplerFrame);
1074 cSpWindowCols->effectiveBW().put(n, chanWidths);
1075
1076 Vector<Double> resolution(nChan, fabs(freqInc));
1077 cSpWindowCols->resolution().put(n, resolution);
1078
1079 cSpWindowCols->totalBandwidth().put(n, bandwidth);
1080 cSpWindowCols->netSideband().put(n, 0);
1081 cSpWindowCols->ifConvChain().put(n, -1);
1082 cSpWindowCols->freqGroup().put(n, 0);
1083 cSpWindowCols->freqGroupName().put(n, " ");
1084 cSpWindowCols->dopplerId().put(n, 0);
1085
1086 // Flags.
1087 cSpWindowCols->flagRow().put(n, False);
1088
1089 return n;
1090}
1091
1092//------------------------------------------------ PKSMS2writer::addStateEntry
1093
1094// Add an entry to the STATE subtable.
1095
1096Int PKSMS2writer::addStateEntry(
1097 const String obsMode)
1098{
1099 // Look for an entry in the STATE subtable.
1100 for (uInt n = 0; n < cStateCols->nrow(); n++) {
1101 if (cStateCols->obsMode()(n) == obsMode) {
1102 return n;
1103 }
1104 }
1105
1106 // Not found, extend the STATE subtable.
1107 cState.addRow();
1108 uInt n = cStateCols->nrow() - 1;
1109
1110 // Data.
1111 if (obsMode.contains("RF")) {
1112 cStateCols->sig().put(n, False);
1113 cStateCols->ref().put(n, True);
1114 } else if (!obsMode.contains("PA")) {
1115 // Signal and reference are both false for "paddle" data.
1116 cStateCols->sig().put(n, True);
1117 cStateCols->ref().put(n, False);
1118 }
1119
1120 cStateCols->load().put(n, 0.0);
1121 cStateCols->cal().put(n, 0.0);
1122 cStateCols->subScan().put(n, 0);
1123 cStateCols->obsMode().put(n, obsMode);
1124
1125 // Flags.
1126 cStateCols->flagRow().put(n, False);
1127
1128 return n;
1129}
1130
1131//----------------------------------------------- PKSMS2writer::addSysCalEntry
1132
1133// Add an entry to the SYSCAL subtable.
1134
1135Int PKSMS2writer::addSysCalEntry(
1136 const Int beamNo,
1137 const Int spWinId,
1138 const Double time,
1139 const Double interval,
1140 const Vector<Float> tcal,
1141 const Vector<Float> tsys,
1142 const Int nPol)
1143{
1144 LogIO os(LogOrigin("PKSMS2writer", "addSysCalEntry()", WHERE));
1145
1146 // Extend the SYSCAL subtable.
1147 cSysCal.addRow();
1148 Int n = cSysCal.nrow() - 1;
1149
1150 //check fo consistency with n pol
1151 //here assume size of Tcal vector = npol
1152 Vector<Float> inTcal(nPol,0);
1153 Int ndim = tcal.shape()(0);
1154 Vector<Float> tmpTcal = tcal;
1155 if (nPol != ndim) {
1156 os << LogIO::WARN
1157 << "Found "<< ndim <<" Tcal value(s) for the data with "<<nPol<<" polarization(s)"
1158 << "(expecting one Tcal per pol)."<<endl
1159 << "First "<< nPol << " Tcal value(s) will be filled." << LogIO::POST;
1160 tmpTcal.resize(nPol, True);
1161 inTcal = tmpTcal;
1162 }
1163 // Keys.
1164 cSysCalCols->antennaId().put(n, 0);
1165 cSysCalCols->feedId().put(n, beamNo-1);
1166 cSysCalCols->spectralWindowId().put(n, spWinId);
1167 cSysCalCols->time().put(n, time);
1168 cSysCalCols->interval().put(n, interval);
1169
1170 // Data.
1171 //cSysCalCols->tcal().put(n, tcal);
1172 cSysCalCols->tcal().put(n, inTcal);
1173 cSysCalCols->tsys().put(n, tsys);
1174
1175 return n;
1176}
1177
1178//---------------------------------------------- PKSMS2writer::addWeatherEntry
1179
1180// Add an entry to the WEATHER subtable.
1181
1182Int PKSMS2writer::addWeatherEntry(
1183 const Double time,
1184 const Double interval,
1185 const Double pressure,
1186 const Double relHumidity,
1187 const Double temperature)
1188{
1189 // Extend the WEATHER subtable.
1190 cWeather.addRow();
1191 Int n = cWeather.nrow() - 1;
1192
1193 // Keys.
1194 cWeatherCols->antennaId().put(n, 0);
1195 cWeatherCols->time().put(n, time);
1196 cWeatherCols->interval().put(n, interval);
1197
1198 // Data.
1199 cWeatherCols->pressure().put(n, pressure);
1200 cWeatherCols->relHumidity().put(n, relHumidity);
1201 cWeatherCols->temperature().put(n, temperature);
1202
1203 return n;
1204}
Note: See TracBrowser for help on using the repository browser.