source: trunk/src/STWriter.cpp@ 969

Last change on this file since 969 was 901, checked in by mar637, 19 years ago

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