source: branches/alma/src/STWriter.cpp@ 1414

Last change on this file since 1414 was 1387, checked in by TakTsutsumi, 17 years ago

merged from NRAO version of ASAP 2.1 with ALMA specific modifications

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