Changeset 39 for trunk/src


Ignore:
Timestamp:
07/07/04 18:51:17 (20 years ago)
Author:
mmarquar
Message:

Added handling of frequencies, SpectralCoordinate and an abscissa retrieval function,

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/SDMemTable.cc

    r22 r39  
    4848#include <aips/Tables/ArrayColumn.h>
    4949#include <aips/Tables/TableRecord.h>
    50 
     50#include <aips/Measures/MFrequency.h>
    5151
    5252#include "SDMemTable.h"
     
    114114  td.addColumn(ScalarColumnDesc<Int>("SCANID")); 
    115115  td.addColumn(ScalarColumnDesc<Double>("INTERVAL")); 
     116  td.addColumn(ArrayColumnDesc<uInt>("FREQID"));
    116117  // Now create a new table from the description.
    117118
     
    163164  uInt n = spec.shape(0)(3);
    164165  chanMask_.resize(n,true);
    165   for (it = whichChans.begin(); it != whichChans.end(); it++) {
     166  for (it = whichChans.begin(); it != whichChans.end(); ++it) {
    166167    if (*it < n)
    167168      chanMask_[*it] = false;
     
    216217  return spectrum;
    217218}
     219
     220std::vector<double> SDMemTable::getAbscissa(Int whichRow,
     221                                            const std::string& whichUnit,
     222                                            double restfreq) {
     223  std::vector<double> absc(nChan());
     224  Vector<Double> absc1(nChan());
     225  indgen(absc1);
     226  ROArrayColumn<uInt> fid(table_, "FREQID");
     227  Vector<uInt> v;
     228  fid.get(whichRow, v);
     229  uInt specidx = v(IFSel_);
     230  cerr << "specidx = " << specidx << endl;
     231  Unit u;
     232  if (whichUnit == "") {
     233    // get unit from table
     234  } else {
     235    u = String(whichUnit);
     236  }
     237  SpectralCoordinate spc = getCoordinate(specidx);
     238  cerr << "debug"  << endl;
     239  if ( u == Unit("km/s") ) {
     240    cerr << "vel ??? " << restfreq << endl;
     241    if (Double(restfreq) >  Double(0.000001)) {
     242      cerr << "converting to velocities"<< endl;
     243      spc.setRestFrequency(Double(restfreq));
     244      spc.setVelocity(u.getName());
     245      Vector<Double> wrld;
     246      spc.pixelToVelocity(wrld,absc1);
     247      std::vector<double>::iterator it;
     248      uInt i = 0;
     249      for (it = absc.begin(); it != absc.end(); ++it) {
     250        (*it) = wrld[i];
     251        i++;
     252      }
     253    }
     254  } else if (u == Unit("Hz")) {
     255    Vector<String> wau(1); wau = u.getName();
     256    spc.setWorldAxisUnits(wau);
     257    cerr << " converting in frequency" << endl;
     258    std::vector<double>::iterator it;
     259    Double tmp;
     260    uInt i = 0;
     261    for (it = absc.begin(); it != absc.end(); ++it) {
     262     
     263      spc.toWorld(tmp,absc1[i]);
     264      (*it) = tmp;
     265      i++;
     266    }
     267    cerr << "converted all pic to world" << endl;
     268  }
     269  cerr << "exiting getAbscissa" << endl;
     270  return absc;
     271}
     272
    218273void SDMemTable::getSpectrum(Vector<Float>& spectrum, Int whichRow=0) {
    219274  ROArrayColumn<Float> spec(table_, "SPECTRA");
     
    324379}
    325380
     381SpectralCoordinate SDMemTable::getCoordinate(uInt whichIdx)  const {
     382
     383  Table t = table_.keywordSet().asTable("FREQUENCIES");
     384  if (whichIdx > t.nrow() ) {
     385    cerr << "SDMemTable::getCoordinate - whichIdx out of range" << endl;
     386    return;
     387  }
     388  Double rp,rv,inc;
     389  String rf;
     390  ROScalarColumn<Double> rpc(t, "REFPIX");
     391  ROScalarColumn<Double> rvc(t, "REFVAL");
     392  ROScalarColumn<Double> incc(t, "INCREMENT");
     393  t.keywordSet().get("REFFRAME",rf);
     394 
     395  MFrequency::Types mft;
     396  if (!MFrequency::getType(mft, rf)) {
     397    cerr << "Frequency type unknown assuming TOPO" << endl;
     398    mft = MFrequency::TOPO;   
     399  }   
     400  rpc.get(whichIdx, rp);
     401  rvc.get(whichIdx, rv);
     402  incc.get(whichIdx, inc);
     403  cerr << "creating speccord from " << whichIdx << ": "
     404       << rp <<", " << rv << ", " << inc << ", " << mft <<endl;
     405  SpectralCoordinate spec(mft,rv,inc,rp);
     406  cerr << "debugit" << endl;
     407  return spec;
     408}
     409
     410bool SDMemTable::putSDFreqTable(const SDFrequencyTable& sdft) {
     411  TableDesc td("", "1", TableDesc::Scratch);
     412  td.addColumn(ScalarColumnDesc<Double>("REFPIX"));
     413  td.addColumn(ScalarColumnDesc<Double>("REFVAL"));
     414  td.addColumn(ScalarColumnDesc<Double>("INCREMENT"));
     415  SetupNewTable aNewTab("freqs", td, Table::New);
     416  Table aTable (aNewTab, Table::Memory, sdft.length());
     417  ScalarColumn<Double> sc0(aTable, "REFPIX");
     418  ScalarColumn<Double> sc1(aTable, "REFVAL");
     419  ScalarColumn<Double> sc2(aTable, "INCREMENT");
     420  for (uInt i=0; i < sdft.length(); ++i) {
     421    sc0.put(i,sdft.referencePixel(i));
     422    sc1.put(i,sdft.referenceValue(i));
     423    sc2.put(i,sdft.increment(i));
     424  }
     425  aTable.rwKeywordSet().define("REFFRAME", sdft.refFrame());
     426  aTable.rwKeywordSet().define("EQUINOX", sdft.equinox());
     427  aTable.rwKeywordSet().define("Unit", String("kms-1"));
     428  table_.rwKeywordSet().defineTable ("FREQUENCIES", aTable);
     429  cerr << "debug - putSDFreqTable" << endl;
     430  return True;
     431}
     432
     433SDFrequencyTable SDMemTable::getSDFreqTable() const  {
     434  SDFrequencyTable sdft;
     435 
     436  return sdft;
     437}
     438
    326439bool SDMemTable::putSDContainer(const SDContainer& sdc) {
    327440  ScalarColumn<Double> mjd(table_, "TIME");
     
    332445  ScalarColumn<Int> scan(table_, "SCANID");
    333446  ScalarColumn<Double> integr(table_, "INTERVAL");
     447  ArrayColumn<uInt> freqid(table_, "FREQID");
    334448
    335449  uInt rno = table_.nrow();
     
    343457  scan.put(rno, sdc.scanid);
    344458  integr.put(rno, sdc.interval);
     459  freqid.put(rno, sdc.getFreqMap());
    345460 
    346461  return true;
     
    355470  ROScalarColumn<Int> scan(table_, "SCANID");
    356471  ROScalarColumn<Double> integr(table_, "INTERVAL");
     472  ROArrayColumn<uInt> freqid(table_, "FREQID");
    357473
    358474  SDContainer sdc(nBeam(),nIF(),nPol(),nChan());
     
    364480  Array<Float> tsys;
    365481  Array<uChar> flagtrum;
     482  Vector<uInt> fmap;
    366483  spec.get(whichRow, spectrum);
    367484  sdc.putSpectrum(spectrum);
     
    370487  ts.get(whichRow, tsys);
    371488  sdc.putTsys(tsys);
     489  freqid.get(whichRow, fmap);
     490  sdc.putFreqMap(fmap);
    372491  return sdc;
    373492}
  • trunk/src/SDMemTable.h

    r22 r39  
    4141#include <aips/Arrays/MaskedArray.h>
    4242
     43#include <trial/Coordinates/SpectralCoordinate.h>
     44
    4345namespace atnf_sd {
    4446
     
    6668
    6769  // put data from meta conatiner into the table
    68   virtual bool putSDContainer(const SDContainer& sdc);
    69   virtual bool putSDHeader(const SDHeader& sdh);
    70   virtual bool putSDFreqTable(const SDFrequencyTable& sdft) {;}
     70  bool putSDContainer(const SDContainer& sdc);
     71  bool putSDHeader(const SDHeader& sdh);
     72  bool putSDFreqTable(const SDFrequencyTable& sdft);
    7173
    7274  //get the dat wrapped up in a meta container
    73   virtual SDContainer getSDContainer(uInt whichRow=0) const;
    74   virtual SDHeader getSDHeader() const;
    75  
     75  SDContainer getSDContainer(uInt whichRow=0) const;
     76  SDHeader getSDHeader() const;
     77  SDFrequencyTable getSDFreqTable() const;
    7678  // get spectrum,mask and tsys for the given row, at the selected
    7779  // cursor - all as stl vectors
     
    127129                                      Bool useSelection = False);
    128130
     131  SpectralCoordinate getCoordinate(uInt whichIdx) const;
     132  std::vector<double> getAbscissa(int whichRow,
     133                                  const std::string& whichUnit="GHz",
     134                                  double restfreq=0.0);
    129135private:
    130136  // set up table structure
Note: See TracChangeset for help on using the changeset viewer.