source: branches/Release2.1.2/src/STWriter.cpp@ 2031

Last change on this file since 2031 was 1305, checked in by mar637, 18 years ago

merge from asap2.1.1 release tag

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