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

Last change on this file since 1449 was 1399, checked in by Malte Marquarding, 17 years ago

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