source: tags/Release2.1.0b/src/STWriter.cpp @ 1241

Last change on this file since 1241 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
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 1241 2006-09-05 10:50:02Z 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"
[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== "FITS") {
[988]68    writer_ = 0;
[450]69  } else if (t== "ASCII") {
[988]70    writer_ = 0;
[450]71  } else {
[988]72    throw (AipsError("Unrecognized export format"));
[28]73  }
74}
75
[822]76STWriter::~STWriter()
[28]77{
[988]78   if (writer_) {
79     delete writer_;
[192]80   }
[28]81}
82
[822]83Int STWriter::setFormat(const std::string &format)
[28]84{
[988]85  if (format != format_) {
86    if (writer_) delete writer_;
[450]87  }
[827]88
[988]89  format_ = format;
90  String t(format_);
[450]91  t.upcase();
92  if (t== "MS2") {
[988]93    writer_ = new PKSMS2writer();
[450]94  } else if (t== "SDFITS") {
[988]95    writer_ = new PKSSDwriter();
[450]96  } else if (t== "FITS") {
[988]97    writer_ = 0;
[450]98  } else if (t== "ASCII") {
[988]99    writer_ = 0;
[450]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
[988]112  if (format_=="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//      }
[988]120  } else if (format_=="ASCII") {
121     STAsciiWriter iw;
[827]122     if (iw.write(*in, filename)) {
[199]123        return 0;
124     } else {
125        return 1;
[988]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();
[1060]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  }
[28]143
[324]144  const Table table = in->table();
[827]145//   ROArrayColumn<uInt> freqIDCol(table, "FREQ_ID");
146//   Vector<uInt> freqIDs;
[324]147
[28]148  // Create the output file and write static data.
149  Int status;
[1060]150  //Bool havexpol = Bool(in->npol() > 2);
151  status = writer_->create(String(filename), hdr.observer, hdr.project,
[28]152                               hdr.antennaname, hdr.antennaposition,
153                               hdr.obstype, hdr.equinox, hdr.freqref,
[1060]154                               nChan, nPol, havexpol, False);
[996]155  if ( status ) {
[717]156    throw(AipsError("Failed to create output file"));
[28]157  }
158
[999]159  Double          srcVel;
[822]160
161  String          fieldName, srcName, tcalTime;
162  Vector<Float>   calFctr, sigma, tcal, tsys;
[999]163  Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2);
[822]164  Matrix<Float>   spectra;
165  Matrix<uChar>   flagtra;
166  Complex         xCalFctr;
[28]167  Int count = 0;
[822]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();
[988]173    TableIterator beamit(stable, "BEAMNO");
[822]174    Int beamno = 1;
175    while (!beamit.pastEnd() ) {
176      Table btable = beamit.table();
[827]177      // position only varies by beam
178      MDirection::ScalarColumn dirCol(btable, "DIRECTION");
179      Vector<Double> direction = dirCol(0).getAngle("rad").getValue();
[988]180      TableIterator cycit(btable, "CYCLENO");
[999]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);
[988]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);
[822]197          // use the first row to fill in all the "metadata"
[827]198          const TableRecord& rec = row.get(0);
[988]199          ROArrayColumn<Float> specCol(itable, "SPECTRA");
[822]200          uInt nchan = specCol(0).nelements();
201          Double cdelt,crval,crpix, restfreq;
[827]202          Float focusAxi, focusTan, focusRot,
203                temperature, pressure, humidity, windSpeed, windAz;
[988]204          Float tmp0,tmp1,tmp2,tmp3,tmp4;
[827]205          Vector<Float> tcalval;
[988]206          String stmp0,stmp1, tcalt;
[822]207          in->frequencies().getEntry(crpix,crval,cdelt, rec.asuInt("FREQ_ID"));
[988]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"));
[827]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"));
[822]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;
[988]222          polConversion(specs, flags, xpol, itable);
223          Vector<Float> tsys = tsysFromTable(itable);
[822]224          // dummy data
[827]225          uInt npol = specs.ncolumn();
[988]226
[822]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);
[996]233          status = writer_->write(scanno, cycno, rec.asDouble("TIME"),
[822]234                                      rec.asDouble("INTERVAL"),
235                                      rec.asString("FIELDNAME"),
236                                      rec.asString("SRCNAME"),
[1072]237                                      srcDir, srcPM, srcVel,hdr.obstype,
[822]238                                      ifno,
239                                      refFreqNew, nchan*abs(cdelt), cdelt,
240                                      restfreq,
[827]241                                      tcal,
242                                      tcalt,
[822]243                                      rec.asFloat("AZIMUTH"),
244                                      rec.asFloat("ELEVATION"),
245                                      rec.asFloat("PARANGLE"),
[988]246                                      focusAxi, focusTan, focusRot,
[827]247                                      temperature,
248                                      pressure, humidity, windSpeed, windAz,
[988]249                                      rec.asInt("REFBEAMNO")+1, beamno,
[827]250                                      direction,
[999]251                                      scanRate,
[827]252                                      tsys,
[822]253                                      sigma, calFctr,// not in scantable
254                                      baseLin, baseSub,// not in scantable
255                                      specs, flags,
256                                      xCalFctr,//
[996]257                                      xpol);
258          if ( status ) {
[988]259            writer_->close();
260            throw(AipsError("STWriter: Failed to export Scantable."));
[822]261          }
[999]262          ++count;
[988]263          ++ifno;
264          ++ifit;
[470]265        }
[988]266        ++cycno;
267        ++cycit;
[28]268      }
[822]269      ++beamno;
270      ++beamit;
[28]271    }
[822]272    ++scanno;
273    ++scanit;
[28]274  }
[717]275  ostringstream oss;
[999]276  oss << "STWriter: wrote " << count << " rows to " << filename;
[717]277  pushLog(String(oss));
[988]278  writer_->close();
[28]279
280  return 0;
281}
[470]282
[827]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) {
[988]289    tmp.resize();
[827]290    tmp = tsysCol(i);
291    out[i] = tmp[0];
292  }
[988]293  return out;
[827]294}
295
296void STWriter::polConversion( Matrix< Float >& specs, Matrix< uChar >& flags,
[822]297                              Vector< Complex > & xpol, const Table & tab )
298{
299  String poltype = tab.keywordSet().asString("POLTYPE");
[1241]300  if ( poltype == "stokes") {
[822]301    String msg = "poltype = " + poltype + " not yet supported in output.";
[988]302    throw(AipsError(msg));
303  }
[827]304  ROArrayColumn<Float> specCol(tab, "SPECTRA");
305  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
[822]306  uInt nchan = specCol(0).nelements();
[988]307  uInt ncol = (tab.nrow()==1 ? 1: 2 );
308  specs.resize(nchan, ncol);
309  flags.resize(nchan, ncol);
[822]310  // the linears
[988]311  for (uInt i=0; i<ncol; ++i) {
[822]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();
[827]318  if ( tab.nrow() == 4 ) {
[822]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) {
[827]324      xpol[k] = Complex(reals[k], imags[k]);
[822]325    }
326  }
327}
[827]328
329
330}
Note: See TracBrowser for help on using the repository browser.