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

Last change on this file since 1635 was 1635, checked in by Malte Marquarding, 15 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-2009
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.15 2009-03-24 06:15:33 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,24,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() == 24) {
472      cBaseSubCol->put(irow, pksrec.baseSub);
473
474    } else {
475      Matrix<Float> tmp(24, 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.