source: branches/polybatch/external/atnf/PKSIO/PKSMS2writer.cc@ 2745

Last change on this file since 2745 was 1720, checked in by Malte Marquarding, 15 years ago

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