source: trunk/src/SDWriter.cc@ 480

Last change on this file since 480 was 470, checked in by kil064, 20 years ago

handle STokes output now in SDFITS/MS formats. Getting messy.
  • 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 470 2005-02-18 00:11:53Z kil064 $
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 cerr << "Failed to create output file." << endl;
180 return 1;
181 }
182
183 Int scanNo = -1;
184 Int cycleNo;
185 Double mjd0 = 0.0;
186//
187 Array<Float> spectra, tSys, stokes;
188 Array<uChar> flags;
189 Bool doLinear = True;
190//
191 Int count = 0;
192 for (Int iRow = 0; iRow < in->nRow(); iRow++) {
193 // Extract the next integration from the table.
194 SDContainer sd = in->getSDContainer(iRow);
195 if (sd.scanid != scanNo) {
196 scanNo = sd.scanid;
197 mjd0 = sd.timestamp;
198 cycleNo = 1;
199 } else if (fabs(sd.timestamp-mjd0) > sd.interval) {
200 cycleNo++;
201 }
202
203 // Get FreqID vector
204 freqIDCol.get(iRow, freqIDs);
205
206 // Get Stokes array (all beams/IFs/Stokes). Its not in SDContainer
207 if (toStokes) {
208 stokes = in->getStokesSpectrum(iRow,-1,-1);
209 }
210
211 IPosition start(asap::nAxes,0);
212 IPosition end(stokes.shape()-1);
213
214 // Write it out beam by beam.
215 for (Int iBeam = 0; iBeam < hdr.nbeam; iBeam++) {
216 start(asap::BeamAxis) = iBeam;
217 end(asap::BeamAxis) = iBeam;
218
219 // Write it out IF by IF.
220 for (Int iIF = 0; iIF < hdr.nif; iIF++) {
221 start(asap::IFAxis) = iIF;
222 end(asap::IFAxis) = iIF;
223 uInt freqID = freqIDs(iIF);
224
225 // None of these are stored in SDMemTable by SDReader.
226 //String fieldName = "";
227 //Vector<Double> srcDir(2, 0.0);
228 Vector<Double> srcPM(2, 0.0);
229 Double srcVel = 0.0;
230
231// The writer will assume refPix = nChan/2 (0-rel). So recompute
232// the frequency at this location.
233
234 Double cdelt = sdft.increment(freqID);
235 Double crval = sdft.referenceValue(freqID);
236 Double crpix = sdft.referencePixel(freqID);
237 Double pixel = nChan/2;
238 Double refFreqNew = (pixel-crpix)*cdelt + crval;
239//
240 //Vector<Float> tcal(2, 0.0f);
241 //String tcalTime = "";
242 //Float azimuth = 0.0f;
243 //Float elevation = 0.0f;
244 //Float parAngle = 0.0f;
245 Float focusAxi = 0.0f;
246 Float focusTan = 0.0f;
247 Float focusRot = 0.0f;
248 Float temperature = 0.0f;
249 Float pressure = 0.0f;
250 Float humidity = 0.0f;
251 Float windSpeed = 0.0f;
252 Float windAz = 0.0f;
253 //Int refBeam = 0;
254 //Vector<Double> direction(2, 0.0);
255 Vector<Double> scanRate(2, 0.0);
256 Vector<Float> sigma(nPol, 0.0f);
257 Vector<Float> calFctr(nPol, 0.0f);
258 Matrix<Float> baseLin(nPol,2, 0.0f);
259 Matrix<Float> baseSub(nPol,9, 0.0f);
260 Complex xCalFctr;
261 Vector<Complex> xPol;
262
263// Get data.
264
265 if (toStokes) {
266
267// Big pain in the ass because
268// 1) Stokes vector may have less polarizations than input
269// 2) Stokes column virtual and not in SDContainer
270// 3) Tsys and FLags already selected for IF and Beam
271// 3) Have to flip order
272
273 Array<Float> t = sd.getTsys(iBeam,iIF);
274 tSys = SDPolUtil::computeStokesTSysForWriter (t, doLinear);
275 Array<uChar> t2 = sd.getFlags(iBeam,iIF);
276 flags = SDPolUtil::computeStokesFlagsForWriter (t2, doLinear);
277 spectra = SDPolUtil::extractStokesForWriter (stokes,start,end);
278 } else {
279 tSys = sd.getTsys(iBeam,iIF);
280 flags = sd.getFlags(iBeam,iIF);
281 spectra = sd.getSpectrum(iBeam,iIF);
282 }
283//
284 if (status = cWriter->write(sd.scanid, cycleNo, sd.timestamp,
285 sd.interval, sd.fieldname, sd.sourcename,
286 sd.getDirection(iBeam),
287 srcPM, srcVel, iIF+1, refFreqNew,
288 nChan*abs(cdelt), cdelt, restFreq, sd.tcal,
289 sd.tcaltime, sd.azimuth, sd.elevation,
290 sd.parangle,
291 focusAxi, focusTan, focusRot, temperature,
292 pressure, humidity, windSpeed, windAz,
293 sd.refbeam, iBeam+1,
294 sd.getDirection(iBeam),
295 scanRate,
296 tSys, sigma, calFctr,
297 baseLin, baseSub,
298 spectra, flags,
299 xCalFctr, xPol)) {
300 cerr << "Error writing output file." << endl;
301 return 1;
302 }
303
304 count++;
305 }
306 }
307 }
308
309 cout << "SDWriter: wrote " << count << " rows to " << filename << endl;
310 cWriter->close();
311
312 return 0;
313}
314
315
Note: See TracBrowser for help on using the repository browser.