source: branches/Release12/src/SDWriter.cc@ 1090

Last change on this file since 1090 was 717, checked in by mar637, 20 years ago

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