source: branches/Release2.1/src/STWriter.cpp@ 3037

Last change on this file since 3037 was 1241, checked in by mar637, 18 years ago

enabled write of circular pols
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.6 KB
Line 
1//#---------------------------------------------------------------------------
2//# STWriter.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: STWriter.cpp 1241 2006-09-05 10:50:02Z 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/Table.h>
45#include <tables/Tables/TableIter.h>
46#include <tables/Tables/TableRow.h>
47#include <tables/Tables/ArrayColumn.h>
48
49//#include "SDFITSImageWriter.h"
50#include "STAsciiWriter.h"
51#include "STHeader.h"
52
53#include "STWriter.h"
54
55using namespace casa;
56namespace asap {
57
58STWriter::STWriter(const std::string &format)
59{
60 format_ = format;
61 String t(format_);
62 t.upcase();
63 if (t== "MS2") {
64 writer_ = new PKSMS2writer();
65 } else if (t== "SDFITS") {
66 writer_ = new PKSSDwriter();
67 } else if (t== "FITS") {
68 writer_ = 0;
69 } else if (t== "ASCII") {
70 writer_ = 0;
71 } else {
72 throw (AipsError("Unrecognized export format"));
73 }
74}
75
76STWriter::~STWriter()
77{
78 if (writer_) {
79 delete writer_;
80 }
81}
82
83Int STWriter::setFormat(const std::string &format)
84{
85 if (format != format_) {
86 if (writer_) delete writer_;
87 }
88
89 format_ = format;
90 String t(format_);
91 t.upcase();
92 if (t== "MS2") {
93 writer_ = new PKSMS2writer();
94 } else if (t== "SDFITS") {
95 writer_ = new PKSSDwriter();
96 } else if (t== "FITS") {
97 writer_ = 0;
98 } else if (t== "ASCII") {
99 writer_ = 0;
100 } else {
101 throw (AipsError("Unrecognized Format"));
102 }
103 return 0;
104}
105
106Int STWriter::write(const CountedPtr<Scantable> in,
107 const std::string &filename)
108{
109
110// Image FITS
111
112 if (format_=="FITS") {
113// Bool verbose = True;
114// SDFITSImageWriter iw;
115// if (iw.write(*in, filename, verbose)) {
116// return 0;
117// } else {
118// return 1;
119// }
120 } else if (format_=="ASCII") {
121 STAsciiWriter iw;
122 if (iw.write(*in, filename)) {
123 return 0;
124 } else {
125 return 1;
126 }
127 }
128
129 // MS or SDFITS
130
131 // Extract the header from the table.
132 STHeader hdr = in->getHeader();
133 //const Int nPol = hdr.npol;
134 //const Int nChan = hdr.nchan;
135 int nIF = in->nif();
136 Vector<uInt> nPol(nIF),nChan(nIF);
137 Vector<Bool> havexpol(nIF);
138 for (int i=0;i<nIF;++i) {
139 nPol(i) = in->npol();
140 nChan(i) = in->nchan(i);
141 havexpol(i) = nPol(i) > 2;
142 }
143
144 const Table table = in->table();
145// ROArrayColumn<uInt> freqIDCol(table, "FREQ_ID");
146// Vector<uInt> freqIDs;
147
148 // Create the output file and write static data.
149 Int status;
150 //Bool havexpol = Bool(in->npol() > 2);
151 status = writer_->create(String(filename), hdr.observer, hdr.project,
152 hdr.antennaname, hdr.antennaposition,
153 hdr.obstype, hdr.equinox, hdr.freqref,
154 nChan, nPol, havexpol, False);
155 if ( status ) {
156 throw(AipsError("Failed to create output file"));
157 }
158
159 Double srcVel;
160
161 String fieldName, srcName, tcalTime;
162 Vector<Float> calFctr, sigma, tcal, tsys;
163 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
164 Matrix<Float> spectra;
165 Matrix<uChar> flagtra;
166 Complex xCalFctr;
167 Int count = 0;
168 Int scanno = 1;
169 // use spearate iterators to ensure renumbering of all numbers
170 TableIterator scanit(table, "SCANNO");
171 while (!scanit.pastEnd() ) {
172 Table stable = scanit.table();
173 TableIterator beamit(stable, "BEAMNO");
174 Int beamno = 1;
175 while (!beamit.pastEnd() ) {
176 Table btable = beamit.table();
177 // position only varies by beam
178 MDirection::ScalarColumn dirCol(btable, "DIRECTION");
179 Vector<Double> direction = dirCol(0).getAngle("rad").getValue();
180 TableIterator cycit(btable, "CYCLENO");
181 ROArrayColumn<Double> srateCol(btable, "SCANRATE");
182 srateCol.get(0, scanRate);
183 ROArrayColumn<Double> spmCol(btable, "SRCPROPERMOTION");
184 spmCol.get(0, srcPM);
185 ROArrayColumn <Double> sdirCol(btable, "SRCDIRECTION");
186 sdirCol.get(0, srcDir);
187 ROScalarColumn<Double> svelCol(btable, "SRCVELOCITY");
188 svelCol.get(0, srcVel);
189 Int cycno = 1;
190 while (!cycit.pastEnd() ) {
191 Table ctable = cycit.table();
192 TableIterator ifit(ctable, "IFNO");
193 Int ifno = 1;
194 while (!ifit.pastEnd() ) {
195 Table itable = ifit.table();
196 TableRow row(itable);
197 // use the first row to fill in all the "metadata"
198 const TableRecord& rec = row.get(0);
199 ROArrayColumn<Float> specCol(itable, "SPECTRA");
200 uInt nchan = specCol(0).nelements();
201 Double cdelt,crval,crpix, restfreq;
202 Float focusAxi, focusTan, focusRot,
203 temperature, pressure, humidity, windSpeed, windAz;
204 Float tmp0,tmp1,tmp2,tmp3,tmp4;
205 Vector<Float> tcalval;
206 String stmp0,stmp1, tcalt;
207 in->frequencies().getEntry(crpix,crval,cdelt, rec.asuInt("FREQ_ID"));
208 in->focus().getEntry(focusAxi, focusTan, focusRot,
209 tmp0,tmp1,tmp2,tmp3,tmp4,
210 rec.asuInt("FOCUS_ID"));
211 in->molecules().getEntry(restfreq,stmp0,stmp1,rec.asuInt("MOLECULE_ID"));
212 in->tcal().getEntry(tcalt,tcalval,rec.asuInt("TCAL_ID"));
213 in->weather().getEntry(temperature, pressure, humidity,
214 windSpeed, windAz,
215 rec.asuInt("WEATHER_ID"));
216 Double pixel = Double(nchan/2);
217 Double refFreqNew = (pixel-crpix)*cdelt + crval;
218 // ok, now we have nrows for the n polarizations in this table
219 Matrix<Float> specs;
220 Matrix<uChar> flags;
221 Vector<Complex> xpol;
222 polConversion(specs, flags, xpol, itable);
223 Vector<Float> tsys = tsysFromTable(itable);
224 // dummy data
225 uInt npol = specs.ncolumn();
226
227 Matrix<Float> baseLin(npol,2, 0.0f);
228 Matrix<Float> baseSub(npol,9, 0.0f);
229 Complex xCalFctr;
230 Vector<Double> scanRate(2, 0.0);
231 Vector<Float> sigma(npol, 0.0f);
232 Vector<Float> calFctr(npol, 0.0f);
233 status = writer_->write(scanno, cycno, rec.asDouble("TIME"),
234 rec.asDouble("INTERVAL"),
235 rec.asString("FIELDNAME"),
236 rec.asString("SRCNAME"),
237 srcDir, srcPM, srcVel,hdr.obstype,
238 ifno,
239 refFreqNew, nchan*abs(cdelt), cdelt,
240 restfreq,
241 tcal,
242 tcalt,
243 rec.asFloat("AZIMUTH"),
244 rec.asFloat("ELEVATION"),
245 rec.asFloat("PARANGLE"),
246 focusAxi, focusTan, focusRot,
247 temperature,
248 pressure, humidity, windSpeed, windAz,
249 rec.asInt("REFBEAMNO")+1, beamno,
250 direction,
251 scanRate,
252 tsys,
253 sigma, calFctr,// not in scantable
254 baseLin, baseSub,// not in scantable
255 specs, flags,
256 xCalFctr,//
257 xpol);
258 if ( status ) {
259 writer_->close();
260 throw(AipsError("STWriter: Failed to export Scantable."));
261 }
262 ++count;
263 ++ifno;
264 ++ifit;
265 }
266 ++cycno;
267 ++cycit;
268 }
269 ++beamno;
270 ++beamit;
271 }
272 ++scanno;
273 ++scanit;
274 }
275 ostringstream oss;
276 oss << "STWriter: wrote " << count << " rows to " << filename;
277 pushLog(String(oss));
278 writer_->close();
279
280 return 0;
281}
282
283Vector<Float> STWriter::tsysFromTable(const Table& tab)
284{
285 ROArrayColumn<Float> tsysCol(tab, "TSYS");
286 Vector<Float> out(tab.nrow());
287 Vector<Float> tmp;
288 for (uInt i=0; i<tab.nrow(); ++i) {
289 tmp.resize();
290 tmp = tsysCol(i);
291 out[i] = tmp[0];
292 }
293 return out;
294}
295
296void STWriter::polConversion( Matrix< Float >& specs, Matrix< uChar >& flags,
297 Vector< Complex > & xpol, const Table & tab )
298{
299 String poltype = tab.keywordSet().asString("POLTYPE");
300 if ( poltype == "stokes") {
301 String msg = "poltype = " + poltype + " not yet supported in output.";
302 throw(AipsError(msg));
303 }
304 ROArrayColumn<Float> specCol(tab, "SPECTRA");
305 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
306 uInt nchan = specCol(0).nelements();
307 uInt ncol = (tab.nrow()==1 ? 1: 2 );
308 specs.resize(nchan, ncol);
309 flags.resize(nchan, ncol);
310 // the linears
311 for (uInt i=0; i<ncol; ++i) {
312 specs.column(i) = specCol(i);
313 flags.column(i) = flagCol(i);
314 }
315 // now the complex if exists
316 Bool hasxpol = False;
317 xpol.resize();
318 if ( tab.nrow() == 4 ) {
319 hasxpol = True;
320 xpol.resize(nchan);
321 Vector<Float> reals, imags;
322 reals = specCol(2); imags = specCol(3);
323 for (uInt k=0; k < nchan; ++k) {
324 xpol[k] = Complex(reals[k], imags[k]);
325 }
326 }
327}
328
329
330}
Note: See TracBrowser for help on using the repository browser.