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

Last change on this file since 1374 was 1325, checked in by mar637, 18 years ago

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