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

Last change on this file since 1427 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.