Changeset 822 for trunk


Ignore:
Timestamp:
02/17/06 10:05:16 (19 years ago)
Author:
mar637
Message:

changes to use the new scantable. this revision is not working.

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STWriter.cpp

    r820 r822  
    11//#---------------------------------------------------------------------------
    2 //# SDWriter.cc: ASAP class to write out single dish spectra.
     2//# STWriter.cc: ASAP class to write out single dish spectra.
    33//#---------------------------------------------------------------------------
    44//# Copyright (C) 2004
    … …  
    4949#include "SDContainer.h"
    5050#include "SDMemTable.h"
    51 #include "SDWriter.h"
     51#include "STWriter.h"
    5252#include "SDFITSImageWriter.h"
    5353#include "SDAsciiWriter.h"
    … …  
    5757using namespace asap;
    5858
    59 //--------------------------------------------------------- SDWriter::SDWriter
     59//--------------------------------------------------------- STWriter::STWriter
    6060
    6161// Default constructor.
    6262
    63 SDWriter::SDWriter(const std::string &format)
     63STWriter::STWriter(const std::string &format)
    6464{
    6565  cFormat = format;
    … …  
    8080}
    8181
    82 //-------------------------------------------------------- SDWriter::~SDWriter
     82//-------------------------------------------------------- STWriter::~STWriter
    8383
    8484// Destructor.
    8585
    86 SDWriter::~SDWriter()
     86STWriter::~STWriter()
    8787{
    8888   if (cWriter) {
    … …  
    9191}
    9292
    93 //-------------------------------------------------------- SDWriter::setFormat
     93//-------------------------------------------------------- STWriter::setFormat
    9494
    9595// Reset the output format.
    9696
    97 Int SDWriter::setFormat(const std::string &format)
     97Int STWriter::setFormat(const std::string &format)
    9898{
    9999  if (format != cFormat) {
    … …  
    118118}
    119119
    120 //------------------------------------------------------------ SDWriter::write
     120//------------------------------------------------------------ STWriter::write
    121121
    122122// Write an SDMemTable to file in the desired format, closing the file when
    123123// finished.
    124124
    125 Int SDWriter::write(const CountedPtr<SDMemTable> in,
     125Int STWriter::write(const CountedPtr<SDMemTable> in,
    126126                    const std::string &filename, Bool toStokes)
    127127{
    … …  
    180180  }
    181181
    182   Int scanNo = -1;
    183   Int cycleNo;
    184   Double mjd0 = 0.0;
    185 //
     182  Vector<Double>  srcPM(2, 0.0);
     183  Double          srcVel = 0.0;
     184
    186185  Array<Float> spectra, tSys, stokes;
    187186  Array<uChar> flags;
    188187  Bool doLinear = True;
    189 //
     188
     189  String          fieldName, srcName, tcalTime;
     190  Vector<Float>   calFctr, sigma, tcal, tsys;
     191  Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2,0.0);
     192  Matrix<Float>   spectra;
     193  Matrix<uChar>   flagtra;
     194  Complex         xCalFctr;
    190195  Int count = 0;
    191   for (Int iRow = 0; iRow < in->nRow(); iRow++) {
    192     // Extract the next integration from the table.
    193     SDContainer sd = in->getSDContainer(iRow);
    194     if (sd.scanid != scanNo) {
    195       scanNo = sd.scanid;
    196       mjd0 = sd.timestamp;
    197       cycleNo = 1;
    198     } else if (fabs(sd.timestamp-mjd0) > sd.interval) {
    199       cycleNo++;
     196  Int scanno = 1;
     197  // use spearate iterators to ensure renumbering of all numbers
     198  TableIterator scanit(table, "SCANNO");
     199  while (!scanit.pastEnd() ) {
     200    Table stable = scanit.table();
     201    TableIterator beamit(table, "BEAMNO");
     202    Int beamno = 1;
     203    while (!beamit.pastEnd() ) {
     204      Table btable = beamit.table();
     205      TableIterator ifit(btable, "IFNO");
     206      Int ifno = 1;
     207      while (!ifit.pastEnd() ) {
     208        Table itable = ifit.table();
     209        TableIterator cycit(itable, "CYCLENO");
     210        Int cycno = 1;
     211        while (!cycit.pastEnd() ) {
     212          Table ctable = cycit.table();
     213          TableRow row(ctable);
     214          // use the first row to fill in all the "metadata"
     215          const TableRecord& rec = row.record(0);
     216          ROArrayColumn<Float> specCol(ctable, "SPECTRA");
     217          uInt nchan = specCol(0).nelements();
     218          Double cdelt,crval,crpix, restfreq;
     219          String tmp,tmp2;
     220          in->frequencies().getEntry(crpix,crval,cdelt, rec.asuInt("FREQ_ID"));
     221          in->molecules().getEntry(restfreq,tmp,tmp2,rec.asuInt("RESTFREQ_ID"));
     222          Double pixel = Double(nchan/2);
     223          Double refFreqNew = (pixel-crpix)*cdelt + crval;
     224          // ok, now we have nrows for the n polarizations in this table
     225          Matrix<Float> specs;
     226          Matrix<uChar> flags;
     227          Vector<Complex> xpol;
     228          polConversion(specs, flags, xpol, ctable);
     229          // dummy data
     230          uInt nPol = specs.ncolumns();
     231          Matrix<Float>   baseLin(npol,2, 0.0f);
     232          Matrix<Float>   baseSub(npol,9, 0.0f);
     233          Complex         xCalFctr;
     234          Vector<Double>  scanRate(2, 0.0);
     235          Vector<Float>   sigma(npol, 0.0f);
     236          Vector<Float>   calFctr(npol, 0.0f);
     237
     238
     239          if (status = cWriter->write(scano, cycno, rec.asDouble("TIME"),
     240                                      rec.asDouble("INTERVAL"),
     241                                      rec.asString("FIELDNAME"),
     242                                      rec.asString("SRCNAME"),
     243                                      direction ,//
     244                                      srcPM, srcVel, // not in scantable yet
     245                                      ifno,
     246                                      refFreqNew, nchan*abs(cdelt), cdelt,
     247                                      restfreq,
     248                                      tcal,//
     249                                      tcaltime,//
     250                                      rec.asFloat("AZIMUTH"),
     251                                      rec.asFloat("ELEVATION"),
     252                                      rec.asFloat("PARANGLE"),
     253                                      focusAxi, focusTan, focusRot,//
     254                                      temperature,//
     255                                      pressure, humidity, windSpeed, windAz,//
     256                                      rec.asInt("REFBEAM"), beamno,
     257                                      direction,//
     258                                      scanRate,// not in scantable
     259                                      tSys, //
     260                                      sigma, calFctr,// not in scantable
     261                                      baseLin, baseSub,// not in scantable
     262                                      specs, flags,
     263                                      xCalFctr,//
     264                                      xpol)
     265                                      ) {
     266            cerr << "Error writing output file." << endl;
     267            return 1;
     268          }
     269
     270          ++cycno;
     271          ++cycit;
     272        }
     273        ++ifno;
     274        ++ifit;
     275      }
     276      ++beamno;
     277      ++beamit;
    200278    }
    201 
    202     // Get FreqID vector
    203     freqIDCol.get(iRow, freqIDs);
    204 
    205     // Get Stokes array (all beams/IFs/Stokes).  Its not in SDContainer
    206     if (toStokes) {
    207        stokes = in->getStokesSpectrum(iRow,-1,-1);
    208     }
    209 
    210     IPosition start(asap::nAxes,0);
    211     IPosition end(stokes.shape()-1);
    212 
    213     // Write it out beam by beam.
    214     for (Int iBeam = 0; iBeam < hdr.nbeam; iBeam++) {
    215        start(asap::BeamAxis) = iBeam;
    216        end(asap::BeamAxis) = iBeam;
    217 
    218       // Write it out IF by IF.
    219       for (Int iIF = 0; iIF < hdr.nif; iIF++) {
    220         start(asap::IFAxis) = iIF;
    221         end(asap::IFAxis) = iIF;
    222         uInt freqID = freqIDs(iIF);
    223 
    224         // None of these are stored in SDMemTable by SDReader.
    225         //String          fieldName = "";
    226         //Vector<Double>  srcDir(2, 0.0);
    227         Vector<Double>  srcPM(2, 0.0);
    228         Double          srcVel = 0.0;
    229 
    230 // The writer will assume refPix = nChan/2 (0-rel).  So recompute
    231 // the frequency at this location.
    232 
    233         Double          cdelt = sdft.increment(freqID);
    234         Double          crval = sdft.referenceValue(freqID);
    235         Double          crpix = sdft.referencePixel(freqID);
    236         Double          pixel = nChan/2;
    237         Double          refFreqNew = (pixel-crpix)*cdelt + crval;
    238 //
    239         //Vector<Float>   tcal(2, 0.0f);
    240         //String          tcalTime = "";
    241         //Float           azimuth = 0.0f;
    242         //Float           elevation = 0.0f;
    243         //Float           parAngle = 0.0f;
    244         Float           focusAxi = 0.0f;
    245         Float           focusTan = 0.0f;
    246         Float           focusRot = 0.0f;
    247         Float           temperature = 0.0f;
    248         Float           pressure = 0.0f;
    249         Float           humidity = 0.0f;
    250         Float           windSpeed = 0.0f;
    251         Float           windAz = 0.0f;
    252         //Int             refBeam = 0;
    253         //Vector<Double>  direction(2, 0.0);
    254         Vector<Double>  scanRate(2, 0.0);
    255         Vector<Float>   sigma(nPol, 0.0f);
    256         Vector<Float>   calFctr(nPol, 0.0f);
    257         Matrix<Float>   baseLin(nPol,2, 0.0f);
    258         Matrix<Float>   baseSub(nPol,9, 0.0f);
    259         Complex         xCalFctr;
    260         Vector<Complex> xPol;
    261 
    262 // Get data. 
    263 
    264         if (toStokes) {
    265 
    266 // Big pain in the ass because
    267 //  1) Stokes vector may have less polarizations than input
    268 //  2) Stokes column virtual and not in SDContainer
    269 //  3) Tsys and FLags already selected for IF and Beam
    270 //  3) Have to flip order
    271 
    272            Array<Float> t = sd.getTsys(iBeam,iIF);
    273            tSys = SDPolUtil::computeStokesTSysForWriter (t, doLinear);
    274            Array<uChar> t2 = sd.getFlags(iBeam,iIF);
    275            flags = SDPolUtil::computeStokesFlagsForWriter (t2, doLinear);
    276            spectra = SDPolUtil::extractStokesForWriter (stokes,start,end);
    277         } else {
    278            tSys = sd.getTsys(iBeam,iIF);
    279            flags = sd.getFlags(iBeam,iIF);
    280            spectra = sd.getSpectrum(iBeam,iIF);       
    281         }
    282 //
    283         if (status = cWriter->write(sd.scanid+1, cycleNo, sd.timestamp,
    284                                     sd.interval, sd.fieldname, sd.sourcename,
    285                                     sd.getDirection(iBeam),
    286                                     srcPM, srcVel, iIF+1, refFreqNew,
    287                                     nChan*abs(cdelt), cdelt, restFreq, sd.tcal,
    288                                     sd.tcaltime, sd.azimuth, sd.elevation,
    289                                     sd.parangle,
    290                                     focusAxi, focusTan, focusRot, temperature,
    291                                     pressure, humidity, windSpeed, windAz,
    292                                     sd.refbeam, iBeam+1,
    293                                     sd.getDirection(iBeam),
    294                                     scanRate,
    295                                     tSys, sigma, calFctr,
    296                                     baseLin, baseSub,
    297                                     spectra, flags,
    298                                     xCalFctr, xPol)) {
    299           cerr << "Error writing output file." << endl;
    300           return 1;
    301         }
    302 
    303         count++;
    304       }
    305     }
     279    ++scanno;
     280    ++scanit;
    306281  }
    307282  ostringstream oss;
    308   oss << "SDWriter: wrote " << count << " rows to " << filename << endl;
     283  oss << "STWriter: wrote " << count << " rows to " << filename << endl;
    309284  pushLog(String(oss));
    310285  cWriter->close();
    … …  
    313288}
    314289
    315 
     290void STWriter::polConversion( Matrix< Float >& specs, Matrix< Float >& flags,
     291                              Vector< Complex > & xpol, const Table & tab )
     292{
     293  TableRow row(tab);
     294  String poltype = tab.keywordSet().asString("POLTYPE");
     295  if ( poltype != "linear")
     296    String msg = "poltype = " + poltype + " not yet supported in output.";
     297    throw(AipsError("msg"));
     298  // use the first row to fill in all the "metadata"
     299  const TableRecord& rec = row.record(0);
     300  ROArrayColumn<Float> specCol(ctable, "SPECTRA");
     301  ROArrayColumn<uChar> flagCol(ctable, "FLAGTRA");
     302  uInt nchan = specCol(0).nelements();
     303  uInt ncols = ( ctable.nrow()==1 ? 1: 2 );
     304  specs.resize(nchan, ncols);
     305  flags.resize(nchan, ncols);
     306  // the linears
     307  for (uInt i=0; i<ncols; ++i) {
     308    specs.column(i) = specCol(i);
     309    flags.column(i) = flagCol(i);
     310  }
     311  // now the complex if exists
     312  Vector<Complex> xpol;
     313  Bool hasxpol = False;
     314  xpol.resize();
     315  if ( ctable.nrow() == 4 ) {
     316    hasxpol = True;
     317    xpol.resize(nchan);
     318    Vector<Float> reals, imags;
     319    reals = specCol(2); imags = specCol(3);
     320    for (uInt k=0; k < nchan; ++k) {
     321  `   xpol[k] = Complex(reals[k], imags[k]);
     322    }
     323  }
     324}

  • trunk/src/STWriter.h

    r820 r822  
    11//#---------------------------------------------------------------------------
    2 //# SDWriter.h: ASAP class to write out single dish spectra.
     2//# STWriter.h: ASAP class to write out single dish spectra.
    33//#---------------------------------------------------------------------------
    44//# Copyright (C) 2004
    … …  
    2929//# $Id$
    3030//#---------------------------------------------------------------------------
    31 #ifndef SDWRITER_H
    32 #define SDWRITER_H
     31#ifndef STWRITER_H
     32#define STWRITER_H
    3333
    3434#include <string>
    … …  
    4646namespace asap {
    4747
    48 class SDWriter : public SDLog {
     48class STWriter : public SDLog {
    4949public:
    50   SDWriter(const string &format = "SDFITS");
    51   ~SDWriter();
     50  STWriter(const string &format = "SDFITS");
     51  ~STWriter();
    5252
    5353// Format can be "SDFITS", "FITS", "MS2" or "ASCII"
    … …  
    5858
    5959private:
     60  void polConversion( casa::Matrix<casa::Float>& spec,
     61                      casa::Matrix<casa::Float>& flag,
     62                      casa::Vector<casa::Complex>& xpol,
     63                      const Table& tab);
    6064  std::string     cFormat;
    6165  PKSwriter *cWriter;
Note: See TracChangeset for help on using the changeset viewer.