source: trunk/src/SDWriter.cc@ 422

Last change on this file since 422 was 324, checked in by kil064, 20 years ago

the majority of the frequency information for each row was not being
written out. make some effort to

1) implement this
2) get it right - not really sure about some of the assumption of the

underlying writer

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 KB
RevLine 
[192]1 //#---------------------------------------------------------------------------
[28]2//# SDWriter.cc: ASAP class to write out single dish spectra.
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
[125]5//# ATNF
[28]6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id: SDWriter.cc 324 2005-01-30 07:01:46Z kil064 $
30//#---------------------------------------------------------------------------
31
32#include <string>
33
[81]34#include <casa/aips.h>
35#include <casa/Arrays.h>
36#include <casa/BasicSL/Complex.h>
37#include <casa/Utilities/CountedPtr.h>
[28]38
39#include <atnf/PKSIO/PKSMS2writer.h>
40#include <atnf/PKSIO/PKSSDwriter.h>
41
[324]42#include <tables/Tables/ArrayColumn.h>
43
44
[60]45#include "SDContainer.h"
46#include "SDMemTable.h"
47#include "SDWriter.h"
[192]48#include "SDFITSImageWriter.h"
[199]49#include "SDAsciiWriter.h"
[28]50
[125]51using namespace casa;
[83]52using namespace asap;
[28]53
54//--------------------------------------------------------- SDWriter::SDWriter
55
56// Default constructor.
57
[125]58SDWriter::SDWriter(const std::string &format)
[28]59{
60 cFormat = format;
[192]61//
[28]62 if (cFormat == "MS2") {
63 cWriter = new PKSMS2writer();
[192]64 } else if (cFormat == "SDFITS") {
[28]65 cWriter = new PKSSDwriter();
[192]66 } else if (cFormat == "FITS") {
67 cWriter = 0;
[199]68 } else if (cFormat == "ASCII") {
69 cWriter = 0;
[28]70 }
71}
72
73//-------------------------------------------------------- SDWriter::~SDWriter
74
75// Destructor.
76
77SDWriter::~SDWriter()
78{
[192]79 if (cWriter) {
80 delete cWriter;
81 }
[28]82}
83
84//-------------------------------------------------------- SDWriter::setFormat
85
86// Reset the output format.
87
[192]88Int SDWriter::setFormat(const std::string &format)
[28]89{
90 if (format != cFormat) {
[192]91 if (cWriter) delete cWriter;
92//
[28]93 cFormat = format;
94 if (cFormat == "MS2") {
95 cWriter = new PKSMS2writer();
[192]96 } else if (cFormat == "SDFITS") {
[28]97 cWriter = new PKSSDwriter();
[192]98 } else if (cFormat == "FITS") {
99 cWriter = 0;
[199]100 } else if (cFormat == "ASCII") {
101 cWriter = 0;
[28]102 }
103 }
104 return 0;
105}
106
107//------------------------------------------------------------ SDWriter::write
108
109// Write an SDMemTable to file in the desired format, closing the file when
110// finished.
111
[324]112Int SDWriter::write(const CountedPtr<SDMemTable> in,
[192]113 const std::string &filename)
[28]114{
[192]115
116// Image FITS
117
118 if (cFormat=="FITS") {
119 Bool verbose = True;
120 SDFITSImageWriter iw;
[324]121 if (iw.write(*in, filename, verbose)) {
[192]122 return 0;
123 } else {
124 return 1;
125 }
[199]126 } else if (cFormat=="ASCII") {
127 SDAsciiWriter iw;
[324]128 if (iw.write(*in, filename)) {
[199]129 return 0;
130 } else {
131 return 1;
132 }
133 }
[192]134
135// MS or SDFITS
136
[28]137 // Extract the header from the table.
[324]138 SDHeader hdr = in->getSDHeader();
139 const Int nPol = hdr.npol;
140 const Int nChan = hdr.nchan;
[28]141
[324]142// Get Freq table
143
144 SDFrequencyTable sdft = in->getSDFreqTable();
145 Vector<Double> restFreqs;
146 String restFreqUnit;
147 sdft.restFrequencies(restFreqs, restFreqUnit);
148 Double restFreq = 0.0;
149 if (restFreqs.nelements()>0) {
150 Quantum<Double> rF(restFreqs(0), Unit(restFreqUnit));
151 restFreq = rF.getValue(Unit("Hz"));
152 }
153
154// Table columns
155
156 const Table table = in->table();
157 ROArrayColumn<uInt> freqIDCol(table, "FREQID");
158 Vector<uInt> freqIDs;
159
[28]160 // Create the output file and write static data.
161 Int status;
162 if (status = cWriter->create(filename, hdr.observer, hdr.project,
163 hdr.antennaname, hdr.antennaposition,
164 hdr.obstype, hdr.equinox, hdr.freqref,
165 nChan, nPol, False, False)) {
166 cerr << "Failed to create output file." << endl;
167 return 1;
168 }
169
170 Int scanNo = -1;
171 Int cycleNo;
172 Double mjd0 = 0.0;
173
174 Int count = 0;
[324]175 for (Int iRow = 0; iRow < in->nRow(); iRow++) {
[28]176 // Extract the next integration from the table.
[324]177 SDContainer sd = in->getSDContainer(iRow);
[28]178 if (sd.scanid != scanNo) {
179 scanNo = sd.scanid;
180 mjd0 = sd.timestamp;
181 cycleNo = 1;
182 } else if (fabs(sd.timestamp-mjd0) > sd.interval) {
183 cycleNo++;
184 }
185
[324]186 // Get FreqID vector
187 freqIDCol.get(iRow, freqIDs);
188
[28]189 // Write it out beam by beam.
190 for (Int iBeam = 0; iBeam < hdr.nbeam; iBeam++) {
191
192 // Write it out IF by IF.
193 for (Int iIF = 0; iIF < hdr.nif; iIF++) {
[324]194 uInt freqID = freqIDs(iIF);
195
[28]196 // None of these are stored in SDMemTable by SDReader.
[112]197 //String fieldName = "";
[79]198 //Vector<Double> srcDir(2, 0.0);
[28]199 Vector<Double> srcPM(2, 0.0);
200 Double srcVel = 0.0;
[324]201
202// The writer will assume refPix = nChan/2 + 1. So recompute
203// the frequency at this location.
204
205 Double cdelt = sdft.increment(freqID);
206 Double crval = sdft.referenceValue(freqID);
207 Double crpix = sdft.referencePixel(freqID);
208 Double pixel = nChan/2 + 1;
209 Double refFreqNew = (pixel-crpix)*cdelt + crval;
210//
[112]211 //Vector<Float> tcal(2, 0.0f);
212 //String tcalTime = "";
213 //Float azimuth = 0.0f;
214 //Float elevation = 0.0f;
215 //Float parAngle = 0.0f;
[28]216 Float focusAxi = 0.0f;
217 Float focusTan = 0.0f;
218 Float focusRot = 0.0f;
219 Float temperature = 0.0f;
220 Float pressure = 0.0f;
221 Float humidity = 0.0f;
222 Float windSpeed = 0.0f;
223 Float windAz = 0.0f;
[112]224 //Int refBeam = 0;
[79]225 //Vector<Double> direction(2, 0.0);
[28]226 Vector<Double> scanRate(2, 0.0);
227 Vector<Float> sigma(nPol, 0.0f);
228 Vector<Float> calFctr(nPol, 0.0f);
229 Matrix<Float> baseLin(nPol,2, 0.0f);
230 Matrix<Float> baseSub(nPol,9, 0.0f);
231 Complex xCalFctr;
232 Vector<Complex> xPol;
[324]233//
[28]234 if (status = cWriter->write(sd.scanid, cycleNo, sd.timestamp,
[112]235 sd.interval, sd.fieldname, sd.sourcename,
[94]236 sd.getDirection(iBeam),
[324]237 srcPM, srcVel, iIF+1, refFreqNew,
238 nChan*abs(cdelt), cdelt, restFreq, sd.tcal,
[112]239 sd.tcaltime, sd.azimuth, sd.elevation,
240 sd.parangle,
[28]241 focusAxi, focusTan, focusRot, temperature,
242 pressure, humidity, windSpeed, windAz,
[112]243 sd.refbeam, iBeam+1,
[94]244 sd.getDirection(iBeam),
245 scanRate,
[75]246 sd.getTsys(iBeam, iIF), sigma, calFctr,
247 baseLin, baseSub,
[28]248 sd.getSpectrum(iBeam, iIF),
249 sd.getFlags(iBeam, iIF),
250 xCalFctr, xPol)) {
251 cerr << "Error writing output file." << endl;
252 return 1;
253 }
254
255 count++;
256 }
257 }
258 }
259
260 cout << "SDWriter: wrote " << count << " rows to " << filename << endl;
261 cWriter->close();
262
263 return 0;
264}
Note: See TracBrowser for help on using the repository browser.