source: trunk/src/STFrequencies.cpp @ 847

//
// C++ Implementation: STFrequencies
//
// Description:
//
//
// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include <casa/iostream.h>
#include <casa/iomanip.h>
#include <casa/Exceptions/Error.h>
#include <casa/Containers/RecordField.h>
#include <casa/Arrays/IPosition.h>

#include <tables/Tables/TableDesc.h>
#include <tables/Tables/SetupNewTab.h>
#include <tables/Tables/ScaColDesc.h>
#include <tables/Tables/TableRecord.h>
#include <tables/Tables/TableParse.h>
#include <tables/Tables/TableRow.h>

#include <coordinates/Coordinates/CoordinateSystem.h>
#include <coordinates/Coordinates/CoordinateUtil.h>

#include "STFrequencies.h"


using namespace casa;

namespace asap {

const String STFrequencies::name_ = "FREQUENCIES";

STFrequencies::STFrequencies(Table::TableType tt) :
  STSubTable( name_, tt )
{
  setup();
}


STFrequencies::~STFrequencies()
{
}

void STFrequencies::setup( )
{
  // add to base class table
  table_.addColumn(ScalarColumnDesc<Double>("REFPIX"));
  table_.addColumn(ScalarColumnDesc<Double>("REFVAL"));
  table_.addColumn(ScalarColumnDesc<Double>("INCREMENT"));

  table_.rwKeywordSet().define("FRAME", String("TOPO"));
  table_.rwKeywordSet().define("BASEFRAME", String("TOPO"));
  table_.rwKeywordSet().define("EQUINOX",String( "J2000"));
  table_.rwKeywordSet().define("UNIT", String("Hz"));
  table_.rwKeywordSet().define("DOPPLER", String("RADIO"));

  // new cached columns
  refpixCol_.attach(table_,"REFPIX");
  refvalCol_.attach(table_,"REFVAL");
  incrCol_.attach(table_,"INCREMENT");
}

uInt STFrequencies::addEntry( Double refpix, Double refval, Double inc )
{
  // test if this already exists
  Table result = table_( near(table_.col("REFVAL"), refval)
                    && near(table_.col("REFPIX"), refpix)
                    && near(table_.col("INCREMENT"), inc) );
  uInt resultid = 0;
  if ( result.nrow() > 0) {
    ROScalarColumn<uInt> c(result, "ID");
    c.get(0, resultid);

  } else {
    uInt rno = table_.nrow();
    table_.addRow();
    // get last assigned freq_id and increment
    if ( rno > 0 ) {
      idCol_.get(rno-1, resultid);
      resultid++;
    }
    refpixCol_.put(rno, refpix);
    refvalCol_.put(rno, refval);
    incrCol_.put(rno, inc);
    idCol_.put(rno, resultid);
  }
  return resultid;
}



void STFrequencies::getEntry( Double& refpix, Double& refval, Double& inc,
                              uInt id )
{
  Table t = table_(table_.col("ID") == Int(id) );
  if (t.nrow() == 0 ) {
    throw(AipsError("STFrequencies::getEntry - freqID out of range"));
  }
  ROTableRow row(t);
  // get first row - there should only be one matching id
  const TableRecord& rec = row.get(0);
  refpix = rec.asDouble("REFPIX");
  refval = rec.asDouble("REFVAL");
  inc = rec.asDouble("INCREMENT");
}

SpectralCoordinate STFrequencies::getSpectralCoordinate( uInt id ) const
{
  Table t = table_(table_.col("ID") == Int(id) );

  if (t.nrow() == 0 ) {
    throw(AipsError("STFrequencies::getSpectralCoordinate - ID out of range"));
  }

  // get the data
  ROTableRow row(t);
  // get first row - there should only be one matching id
  const TableRecord& rec = row.get(0);

  return SpectralCoordinate( getFrame(true), rec.asDouble("REFVAL"),
                             rec.asDouble("INCREMENT"),
                             rec.asDouble("REFPIX"));
}

SpectralCoordinate
  asap::STFrequencies::getSpectralCoordinate( const MDirection& md,
                                              const MPosition& mp,
                                              const MEpoch& me,
                                              Double restfreq, uInt id ) const
{
  SpectralCoordinate spc = getSpectralCoordinate(id);
  spc.setRestFrequency(restfreq, True);
  if ( !spc.setReferenceConversion(getFrame(), me, mp, md) ) {
    throw(AipsError("Couldn't convert frequency frame."));
  }
  String unitstr = getUnitString();
  if ( !unitstr.empty() ) {
    Unit unitu(unitstr);
    if ( unitu == Unit("Hz") ) {
    } else {
      spc.setVelocity(unitstr, getDoppler());
    }
  }
  return spc;
}


void STFrequencies::rescale( Float factor, const std::string& mode )
{
  TableRow row(table_);
  TableRecord& outrec = row.record();
  RecordFieldPtr<Double> rv(outrec, "REFVAL");
  RecordFieldPtr<Double> rp(outrec, "REFPIX");
  RecordFieldPtr<Double> inc(outrec, "INCREMENT");
  for (uInt i=0; i<table_.nrow(); ++i) {

    const TableRecord& rec = row.get(i);

    SpectralCoordinate sc ( getFrame(true), rec.asDouble("REFVAL"),
                            rec.asDouble("INCREMENT"), rec.asDouble("REFPIX") );

    SpectralCoordinate scout;
    if (mode == "BIN") {
      scout = binCsys(sc, Int(factor));
    } else if (mode == "RESAMPLE") {
      scout = resampleCsys(sc, factor);
    }
    *rv = scout.referenceValue()[0];
    *rp = scout.referencePixel()[0];
    *inc = scout.increment()[0];
    row.put(i);
  }
}

SpectralCoordinate STFrequencies::binCsys(const SpectralCoordinate& sc,
                                          Int factor)
{
  CoordinateSystem csys;
  csys.addCoordinate(sc);
  IPosition factors(1, factor);
  CoordinateSystem binnedcs =
    CoordinateUtil::makeBinnedCoordinateSystem(factors, csys, False);
  return binnedcs.spectralCoordinate(0);
}

SpectralCoordinate STFrequencies::resampleCsys(const SpectralCoordinate& sc,
                                               Float width)
{
  Vector<Float> offset(1,0.0);
  Vector<Float> factors(1,1.0/width);
  Vector<Int> newshape;
  CoordinateSystem csys;
  csys.addCoordinate(sc);
  CoordinateSystem csys2 = csys.subImage(offset, factors, newshape);
  return csys2.spectralCoordinate(0);
}


MFrequency::Types STFrequencies::getFrame(bool base) const
{
  // get the ref frame
  String rf = table_.keywordSet().asString("BASEFRAME");

  // Create SpectralCoordinate (units Hz)
  MFrequency::Types mft;
  if (!MFrequency::getType(mft, rf)) {
    ostringstream oss;
    pushLog("WARNING: Frequency type unknown assuming TOPO");
    mft = MFrequency::TOPO;
  }

  return mft;
}

std::string asap::STFrequencies::getFrameString( bool base ) const
{
  if ( base ) return table_.keywordSet().asString("BASEFRAME");
  else return table_.keywordSet().asString("FRAME");
}

std::string asap::STFrequencies::getUnitString( ) const
{
  return table_.keywordSet().asString("UNIT");
}

Unit asap::STFrequencies::getUnit( ) const
{
  return Unit(table_.keywordSet().asString("UNIT"));
}

std::string asap::STFrequencies::getDopplerString( ) const
{
  return table_.keywordSet().asString("DOPPLER");
}

MDoppler::Types asap::STFrequencies::getDoppler( ) const
{
  String dpl = table_.keywordSet().asString("DOPPLER");

  // Create SpectralCoordinate (units Hz)
  MDoppler::Types mdt;
  if (!MDoppler::getType(mdt, dpl)) {
    throw(AipsError("Doppler type unknown"));
  }
  return mdt;
}

std::string asap::STFrequencies::print( int id )
{
  Table t;
  ostringstream oss;
  if ( id < 0 ) t = table_;
  else  t = table_(table_.col("ID") == Int(id) );
  ROTableRow row(t);
  for (uInt i=0; i<t.nrow(); ++i) {
    const TableRecord& rec = row.get(i);
    oss <<  setw(8)
    << "frame" << setw(16) << setprecision(8)
    << rec.asDouble("REFVAL") << setw(10)
    << rec.asDouble("REFPIX") << setw(12)
    << rec.asDouble("INCREMENT") << endl;
  }
  return String(oss);
}

float STFrequencies::getRefFreq( uInt id, uInt channel )
{
  Table t = table_(table_.col("ID") == Int(id) );
  if ( t.nrow() == 0 ) throw(AipsError("Selected Illegal frequency id"));
  ROTableRow row(t);
  const TableRecord& rec = row.get(0);
  return (Double(channel/2) - rec.asDouble("REFPIX"))
          * rec.asDouble("INCREMENT") + rec.asDouble("REFVAL");
}

bool asap::STFrequencies::conformant( const STFrequencies& other ) const
{
  const Record& r = table_.keywordSet();
  const Record& ro = other.table_.keywordSet();
  return ( r.asString("FRAME") == ro.asString("FRAME") &&
           r.asString("EQUINOX") == ro.asString("EQUINOX") &&
           r.asString("UNIT") == ro.asString("UNIT") &&
           r.asString("DOPPLER") == ro.asString("DOPPLER")
          );
}

std::vector< std::string > asap::STFrequencies::getInfo( ) const
{
  const Record& r = table_.keywordSet();
  std::vector<std::string> out;
  out.push_back(r.asString("UNIT"));
  out.push_back(r.asString("FRAME"));
  out.push_back(r.asString("DOPPLER"));
}

void asap::STFrequencies::setInfo( const std::vector< std::string >& theinfo )
{
  if ( theinfo.size() != 3 ) throw(AipsError("setInfo needs three parameters"));
  String un,rfrm,dpl;
  un = theinfo[0];rfrm = theinfo[1];dpl = theinfo[2];
  TableRecord& r = table_.rwKeywordSet();
  setFrame(rfrm);
  MDoppler::Types dtype;
  dpl.upcase();
  if (!MDoppler::getType(dtype, dpl)) {
    throw(AipsError("Doppler type unknown"));
  } else {
    r.define("DOPPLER",dpl);
  }
}
void asap::STFrequencies::setFrame( const std::string & frame )
{
  MFrequency::Types mdr;
  if (!MFrequency::getType(mdr, frame)) {
    Int a,b;const uInt* c;
    const String* valid = MFrequency::allMyTypes(a, b, c);
    String pfix = "Please specify a legal frame type. Types are\n";
    throw(AipsError(pfix+(*valid)));
  } else {
    table_.rwKeywordSet().define("FRAME", frame);
  }
}

} // namespace