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

Last change on this file since 1720 was 1720, checked in by Malte Marquarding, 14 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.