source: trunk/src/SDWriter.cc@ 459

Last change on this file since 459 was 450, checked in by kil064, 20 years ago

make format case insensitive

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.2 KB
Line 
1 //#---------------------------------------------------------------------------
2//# SDWriter.cc: ASAP class to write out single dish spectra.
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
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 450 2005-02-15 09:06:11Z kil064 $
30//#---------------------------------------------------------------------------
31
32#include <string>
33
34#include <casa/aips.h>
35#include <casa/Arrays.h>
36#include <casa/BasicSL/Complex.h>
37#include <casa/Utilities/CountedPtr.h>
38
39#include <atnf/PKSIO/PKSMS2writer.h>
40#include <atnf/PKSIO/PKSSDwriter.h>
41
42#include <tables/Tables/ArrayColumn.h>
43
44
45#include "SDContainer.h"
46#include "SDMemTable.h"
47#include "SDWriter.h"
48#include "SDFITSImageWriter.h"
49#include "SDAsciiWriter.h"
50
51using namespace casa;
52using namespace asap;
53
54//--------------------------------------------------------- SDWriter::SDWriter
55
56// Default constructor.
57
58SDWriter::SDWriter(const std::string &format)
59{
60 cFormat = format;
61//
62 String t(cFormat);
63 t.upcase();
64 if (t== "MS2") {
65 cWriter = new PKSMS2writer();
66 } else if (t== "SDFITS") {
67 cWriter = new PKSSDwriter();
68 } else if (t== "FITS") {
69 cWriter = 0;
70 } else if (t== "ASCII") {
71 cWriter = 0;
72 } else {
73 throw (AipsError("Unrecognized Format"));
74 }
75}
76
77//-------------------------------------------------------- SDWriter::~SDWriter
78
79// Destructor.
80
81SDWriter::~SDWriter()
82{
83 if (cWriter) {
84 delete cWriter;
85 }
86}
87
88//-------------------------------------------------------- SDWriter::setFormat
89
90// Reset the output format.
91
92Int SDWriter::setFormat(const std::string &format)
93{
94 if (format != cFormat) {
95 if (cWriter) delete cWriter;
96 }
97//
98 cFormat = format;
99 String t(cFormat);
100 t.upcase();
101 if (t== "MS2") {
102 cWriter = new PKSMS2writer();
103 } else if (t== "SDFITS") {
104 cWriter = new PKSSDwriter();
105 } else if (t== "FITS") {
106 cWriter = 0;
107 } else if (t== "ASCII") {
108 cWriter = 0;
109 } else {
110 throw (AipsError("Unrecognized Format"));
111 }
112 return 0;
113}
114
115//------------------------------------------------------------ SDWriter::write
116
117// Write an SDMemTable to file in the desired format, closing the file when
118// finished.
119
120Int SDWriter::write(const CountedPtr<SDMemTable> in,
121 const std::string &filename, Bool toStokes)
122{
123
124// Image FITS
125
126 if (cFormat=="FITS") {
127 Bool verbose = True;
128 SDFITSImageWriter iw;
129 if (iw.write(*in, filename, verbose, toStokes)) {
130 return 0;
131 } else {
132 return 1;
133 }
134 } else if (cFormat=="ASCII") {
135 SDAsciiWriter iw;
136 if (iw.write(*in, filename, toStokes)) {
137 return 0;
138 } else {
139 return 1;
140 }
141 }
142
143// MS or SDFITS
144
145 // Extract the header from the table.
146 SDHeader hdr = in->getSDHeader();
147 const Int nPol = hdr.npol;
148 const Int nChan = hdr.nchan;
149 if (toStokes) {
150 cerr << "Stokes conversion not yet available with SDFITS or MS" << endl;
151 }
152
153// Get Freq table
154
155 SDFrequencyTable sdft = in->getSDFreqTable();
156 Vector<Double> restFreqs;
157 String restFreqUnit;
158 sdft.restFrequencies(restFreqs, restFreqUnit);
159 Double restFreq = 0.0;
160 if (restFreqs.nelements()>0) {
161 Quantum<Double> rF(restFreqs(0), Unit(restFreqUnit));
162 restFreq = rF.getValue(Unit("Hz"));
163 }
164
165// Table columns
166
167 const Table table = in->table();
168 ROArrayColumn<uInt> freqIDCol(table, "FREQID");
169 Vector<uInt> freqIDs;
170
171 // Create the output file and write static data.
172 Int status;
173 if (status = cWriter->create(filename, hdr.observer, hdr.project,
174 hdr.antennaname, hdr.antennaposition,
175 hdr.obstype, hdr.equinox, hdr.freqref,
176 nChan, nPol, False, False)) {
177 cerr << "Failed to create output file." << endl;
178 return 1;
179 }
180
181 Int scanNo = -1;
182 Int cycleNo;
183 Double mjd0 = 0.0;
184
185 Int count = 0;
186 for (Int iRow = 0; iRow < in->nRow(); iRow++) {
187 // Extract the next integration from the table.
188 SDContainer sd = in->getSDContainer(iRow);
189 if (sd.scanid != scanNo) {
190 scanNo = sd.scanid;
191 mjd0 = sd.timestamp;
192 cycleNo = 1;
193 } else if (fabs(sd.timestamp-mjd0) > sd.interval) {
194 cycleNo++;
195 }
196
197 // Get FreqID vector
198 freqIDCol.get(iRow, freqIDs);
199
200 // Write it out beam by beam.
201 for (Int iBeam = 0; iBeam < hdr.nbeam; iBeam++) {
202
203 // Write it out IF by IF.
204 for (Int iIF = 0; iIF < hdr.nif; iIF++) {
205 uInt freqID = freqIDs(iIF);
206
207 // None of these are stored in SDMemTable by SDReader.
208 //String fieldName = "";
209 //Vector<Double> srcDir(2, 0.0);
210 Vector<Double> srcPM(2, 0.0);
211 Double srcVel = 0.0;
212
213// The writer will assume refPix = nChan/2 (0-rel). So recompute
214// the frequency at this location.
215
216 Double cdelt = sdft.increment(freqID);
217 Double crval = sdft.referenceValue(freqID);
218 Double crpix = sdft.referencePixel(freqID);
219 Double pixel = nChan/2;
220 Double refFreqNew = (pixel-crpix)*cdelt + crval;
221//
222 //Vector<Float> tcal(2, 0.0f);
223 //String tcalTime = "";
224 //Float azimuth = 0.0f;
225 //Float elevation = 0.0f;
226 //Float parAngle = 0.0f;
227 Float focusAxi = 0.0f;
228 Float focusTan = 0.0f;
229 Float focusRot = 0.0f;
230 Float temperature = 0.0f;
231 Float pressure = 0.0f;
232 Float humidity = 0.0f;
233 Float windSpeed = 0.0f;
234 Float windAz = 0.0f;
235 //Int refBeam = 0;
236 //Vector<Double> direction(2, 0.0);
237 Vector<Double> scanRate(2, 0.0);
238 Vector<Float> sigma(nPol, 0.0f);
239 Vector<Float> calFctr(nPol, 0.0f);
240 Matrix<Float> baseLin(nPol,2, 0.0f);
241 Matrix<Float> baseSub(nPol,9, 0.0f);
242 Complex xCalFctr;
243 Vector<Complex> xPol;
244//
245 if (status = cWriter->write(sd.scanid, cycleNo, sd.timestamp,
246 sd.interval, sd.fieldname, sd.sourcename,
247 sd.getDirection(iBeam),
248 srcPM, srcVel, iIF+1, refFreqNew,
249 nChan*abs(cdelt), cdelt, restFreq, sd.tcal,
250 sd.tcaltime, sd.azimuth, sd.elevation,
251 sd.parangle,
252 focusAxi, focusTan, focusRot, temperature,
253 pressure, humidity, windSpeed, windAz,
254 sd.refbeam, iBeam+1,
255 sd.getDirection(iBeam),
256 scanRate,
257 sd.getTsys(iBeam, iIF), sigma, calFctr,
258 baseLin, baseSub,
259 sd.getSpectrum(iBeam, iIF),
260 sd.getFlags(iBeam, iIF),
261 xCalFctr, xPol)) {
262 cerr << "Error writing output file." << endl;
263 return 1;
264 }
265
266 count++;
267 }
268 }
269 }
270
271 cout << "SDWriter: wrote " << count << " rows to " << filename << endl;
272 cWriter->close();
273
274 return 0;
275}
Note: See TracBrowser for help on using the repository browser.