source: trunk/external/atnf/PKSIO/PKSMS2writer.cc@ 1627

Last change on this file since 1627 was 1452, checked in by Malte Marquarding, 16 years ago

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