Changeset 1757 for branches/alma/src


Ignore:
Timestamp:
06/09/10 19:03:06 (15 years ago)
Author:
Kana Sugimoto
Message:

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.


Location:
branches/alma
Files:
6 added
28 edited

Legend:

Unmodified
Added
Removed
  • branches/alma

  • branches/alma/src/Makefile

    r1704 r1757  
    123123             STAsciiWriter.o \
    124124             STFITSImageWriter.o \
     125             STAtmosphere.o \
    125126             Scantable.o \
    126127             Templates.o
     
    136137             python_LineCatalog.o \
    137138             python_SrcType.o \
     139             python_STAtmosphere.o \
     140             python_STCoordinate.o \
    138141             python_asap.o
    139142
     
    172175             STWriter.h \
    173176             STAsciiWriter.h \
    174              STFITSImageWriter.h
     177             STFITSImageWriter.h \
     178             IndexedCompare.h \
     179             STAtmosphere.h \
     180             STCoordinate.h
    175181
    176182STATICCCLIB := libasap.a
  • branches/alma/src/MathUtils.cpp

    r1603 r1757  
    3939#include <scimath/Mathematics/MedianSlider.h>
    4040#include <casa/Exceptions/Error.h>
     41
     42#include <scimath/Fitting/LinearFit.h>
     43#include <scimath/Functionals/Polynomial.h>
     44#include <scimath/Mathematics/AutoDiff.h>
     45
    4146
    4247#include "MathUtils.h"
     
    183188  MedianSlider ms(hwidth);
    184189  Slice sl(0, fwidth-1);
    185   Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl), 
     190  Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
    186191                  const_cast<Vector<Bool>& >(flag)(sl));
    187192  uInt n = in.nelements();
    188193  for (uInt i=hwidth; i<(n-hwidth); ++i) {
    189194    // add data value
    190     out[i] = ms.add(in[i+hwidth], flag[i+hwidth]); 
    191     outflag[i] = (ms.nval() == 0);   
    192   }
    193   // replicate edge values from fisrt value with full width of values
     195    out[i] = ms.add(in[i+hwidth], flag[i+hwidth]);
     196    outflag[i] = (ms.nval() == 0);
     197  }
     198  // replicate edge values from first value with full width of values
    194199  for (uInt i=0;i<hwidth;++i) {
    195200    out[i] = out[hwidth];
    196     outflag[i] = outflag[hwidth];   
     201    outflag[i] = outflag[hwidth];
    197202    out[n-1-i] = out[n-1-hwidth];
    198     outflag[n-1-i] = outflag[n-1-hwidth];   
    199   }
    200 }
     203    outflag[n-1-i] = outflag[n-1-hwidth];
     204  }
     205}
     206
     207void mathutil::polyfit(Vector<Float>& out, Vector<Bool>& outmask,
     208                       const Vector<Float>& in, const Vector<Bool>& mask,
     209                       float width, int order)
     210{
     211  Int hwidth = Int(width+0.5);
     212  Int fwidth = hwidth*2+1;
     213  out.resize(in.nelements());
     214  outmask.resize(mask.nelements());
     215  LinearFit<Float> fitter;
     216  Polynomial<Float> poly(order);
     217  fitter.setFunction(poly);
     218  Vector<Float> sigma(fwidth);
     219  sigma = 1.0;
     220  Vector<Float> parms;
     221  Vector<Float> x(fwidth);
     222  indgen(x);
     223
     224  uInt n = in.nelements();
     225
     226  for (uInt i=hwidth; i<(n-hwidth); ++i) {
     227    // add data value
     228    if (mask[i]) {
     229      Slice sl(i-hwidth, fwidth);
     230      const Vector<Float> &y = const_cast<Vector<Float>& >(in)(sl);
     231      const Vector<Bool> &m = const_cast<Vector<Bool>& >(mask)(sl);
     232      parms = fitter.fit(x, y, sigma, &m);
     233
     234      poly.setCoefficients(parms);
     235      out[i] = poly(x[hwidth]);//cout << in[i] <<"->"<<out[i]<<endl;
     236    } else {
     237      out[i] = in[i];
     238    }
     239    outmask[i] = mask[i];
     240  }
     241  // replicate edge values from first value with full width of values
     242  for (uInt i=0;i<hwidth;++i) {
     243    out[i] = out[hwidth];
     244    outmask[i] = outmask[hwidth];
     245    out[n-1-i] = out[n-1-hwidth];
     246    outmask[n-1-i] = outmask[n-1-hwidth];
     247  }
     248}
  • branches/alma/src/MathUtils.h

    r1603 r1757  
    5151 * @param ignoreOther drop every second channel (NYI)
    5252 */
    53 void hanning(casa::Vector<casa::Float>& out, 
     53void hanning(casa::Vector<casa::Float>& out,
    5454             casa::Vector<casa::Bool>& outmask,
    55              const casa::Vector<casa::Float>& in, 
     55             const casa::Vector<casa::Float>& in,
    5656             const casa::Vector<casa::Bool>& mask,
    5757             casa::Bool relaxed=casa::False,
     
    6060/**
    6161 * Apply a running median to  a masked vector.
    62  * Edge solution:  The first and last hwidth channels will be replicated 
     62 * Edge solution:  The first and last hwidth channels will be replicated
    6363 * from the first/last value from a full window.
    6464 * @param out the smoothed vector
     
    6868 * @param hwidth half-width of the smoothing window
    6969 */
    70  void runningMedian(casa::Vector<casa::Float>& out,
     70 void runningMedian(casa::Vector<casa::Float>& out,
     71                   casa::Vector<casa::Bool>& outflag,
     72                   const casa::Vector<casa::Float>& in,
     73                   const casa::Vector<casa::Bool>& flag,
     74                   float hwidth);
     75
     76 void polyfit(casa::Vector<casa::Float>& out,
    7177                   casa::Vector<casa::Bool>& outmask,
    72                    const casa::Vector<casa::Float>& in, 
     78                   const casa::Vector<casa::Float>& in,
    7379                   const casa::Vector<casa::Bool>& mask,
    74                    float hwidth);
     80                   float hwidth, int order);
    7581
    7682// Generate specified statistic
  • branches/alma/src/RowAccumulator.cpp

    r1603 r1757  
    6767  Float totalweight = weight;
    6868  MaskedArray<Float> data(v,m);
    69   if ( weightType_ == asap::VAR ) {
     69  if ( weightType_ == asap::W_VAR ) {
    7070    if (m.nelements() == userMask_.nelements()) {
    7171      Float fac = 1.0/variance(data(userMask_));
     
    9090  Float w = 1.0;
    9191  tsysSum_ += v[0];
    92   if ( weightType_ == asap::TSYS  || weightType_ == asap::TINTSYS ) {
     92  if ( weightType_ == asap::W_TSYS  || weightType_ == asap::W_TINTSYS ) {
    9393    w /= (v[0]*v[0]);
    9494  }
     
    105105  Float w = 1.0;
    106106  intervalSum_ += inter;
    107   if ( weightType_ == asap::TINT || weightType_ == asap::TINTSYS ) {
     107  if ( weightType_ == asap::W_TINT || weightType_ == asap::W_TINTSYS ) {
    108108    w /= Float(inter);
    109109  }
  • branches/alma/src/RowAccumulator.h

    r1603 r1757  
    3333   * Constructor taking a weight type as defined in @ref STDefs
    3434   */
    35   explicit RowAccumulator(WeightType wt = asap::NONE);
     35  explicit RowAccumulator(WeightType wt = asap::W_NONE);
    3636
    3737 ~RowAccumulator();
  • branches/alma/src/STAsciiWriter.cpp

    r1657 r1757  
    8989
    9090   String rootName(fileName);
    91    if (rootName.length()==0) rootName = String("ascii");
    9291
    9392  Block<String> cols(4);
     
    104103    String dirtype = stable.getDirectionRefString();
    105104    ostringstream onstr;
    106     onstr << "SCAN" << rec.asuInt("SCANNO")
    107     << "_CYCLE" << rec.asuInt("CYCLENO")
    108     << "_BEAM" << rec.asuInt("BEAMNO")
    109     << "_IF" << rec.asuInt("IFNO");
     105
     106    if (rootName.length()==0) {
     107      rootName = String("ascii");
     108    }
     109    if (tab.nrow() > 1) {
     110      if (stable.nscan() > 1)
     111        onstr << "_SCAN" << rec.asuInt("SCANNO");
     112      if (stable.ncycle(rec.asuInt("SCANNO")) > 1)
     113        onstr << "_CYCLE" << rec.asuInt("CYCLENO");
     114      if (stable.nbeam(rec.asuInt("SCANNO")) > 1)
     115        onstr << "_BEAM" << rec.asuInt("BEAMNO");
     116      if (stable.nif(rec.asuInt("SCANNO")) > 1)
     117        onstr << "_IF" << rec.asuInt("IFNO");
     118    }
     119
    110120    String fName = rootName + String(onstr) + String(".txt");
    111121    ofstream of(fName.chars(), ios::trunc);
  • branches/alma/src/STAttr.cpp

    r1387 r1757  
    331331   TidGainElPoly_(2) = -3.219093e-4;
    332332   
     333   // 2009-09-15 - 13mm (22.2GHz) receiver
    333334   ParkesGainElPoly_.resize(3);
    334    ParkesGainElPoly_(0) = 0.296759e-1;
    335    ParkesGainElPoly_(1) = -0.293124e-3;
    336    ParkesGainElPoly_(2) = 0.264295e-6;
     335   ParkesGainElPoly_(0) = -0.194031;
     336   ParkesGainElPoly_(1) = 0.457724e-1;
     337   ParkesGainElPoly_(2) = -0.438659e-3;
    337338}
    338339
  • branches/alma/src/STDefs.h

    r1388 r1757  
    3434
    3535namespace asap {
    36   enum AxisNo { BeamAxis=0,
    37                 IFAxis,
    38                 PolAxis,
    39                 ChanAxis,
    40                 nAxes};
    4136
    4237  enum Instrument {UNKNOWNINST=0,
     
    5348  enum FeedPolType {UNKNOWNFEED, LINEAR, CIRCULAR, N_POL};
    5449
    55   enum WeightType {NONE=0, VAR, TSYS, TINT, TINTSYS, N_WEIGHTTYPES};
    56 
    57   enum TableType {MEMORY=0, PERSISTENT};
    58 
     50  enum WeightType {W_NONE=0, W_VAR, W_TSYS, W_TINT, W_TINTSYS, W_N_WEIGHTTYPES};
    5951
    6052
  • branches/alma/src/STFITSImageWriter.cpp

    r1604 r1757  
    121121  TableIterator iter(tab, cols);
    122122  // Open data file
     123
    123124  while ( !iter.pastEnd() ) {
    124125    Table t = iter.table();
     
    127128    String dirtype = stable.getDirectionRefString();
    128129    ostringstream onstr;
    129     onstr << "SCAN" << rec.asuInt("SCANNO")
    130     << "_CYCLE" << rec.asuInt("CYCLENO")
    131     << "_BEAM" << rec.asuInt("BEAMNO")
    132     << "_IF" << rec.asuInt("IFNO")
    133     << "_POL" << rec.asuInt("POLNO");
     130    if (rootName.length()==0) {
     131      rootName = "fits";
     132    }
     133    if (tab.nrow() > 1) {
     134      if (stable.nscan() > 1)
     135        onstr << "_SCAN" << rec.asuInt("SCANNO");
     136      if (stable.ncycle(rec.asuInt("SCANNO")) > 1)
     137        onstr << "_CYCLE" << rec.asuInt("CYCLENO");
     138      if (stable.nbeam(rec.asuInt("SCANNO")) > 1)
     139        onstr << "_BEAM" << rec.asuInt("BEAMNO");
     140      if (stable.nif(rec.asuInt("SCANNO")) > 1)
     141        onstr << "_IF" << rec.asuInt("IFNO");
     142      if (stable.npol(rec.asuInt("SCANNO")) > 1)
     143        onstr << "_POL" << rec.asuInt("POLNO");
     144    }
    134145    String fileName = rootName + String(onstr) + String(".fits");
    135146    int row0 = t.rowNumbers(tab)[0];
     
    260271    }
    261272    fits_close_file(fptr, &status);
    262  
     273    ostringstream oss;
     274    oss << "Wrote " << fileName;
     275    pushLog(String(oss)); 
    263276    //pushLog(String(oss));
    264277    ++iter;
  • branches/alma/src/STFiller.cpp

    r1684 r1757  
    264264    freqFrame = "REST";
    265265  }
    266   table_->frequencies().setFrame(freqFrame);
    267      
     266  // set both "FRAME" and "BASEFRAME"
     267  table_->frequencies().setFrame(freqFrame, false);
     268  table_->frequencies().setFrame(freqFrame,true);
     269  //table_->focus().setParallactify(true);
    268270}
    269271
     
    329331    if ( status != 0 ) break;
    330332    n += 1;
    331 
    332333    Regex filterrx(".*[SL|PA]$");
    333334    Regex obsrx("^AT.+");
     
    391392    /// @todo this has to change when nchan isn't global anymore
    392393    //id = table_->frequencies().addEntry(Double(header_->nchan/2),
    393     //                                        pksrec.refFreq, pksrec.freqInc);
     394    //                                    pksrec.refFreq, pksrec.freqInc);
    394395    if ( pksrec.nchan == 1 ) {
    395396      id = table_->frequencies().addEntry(Double(0),
     
    416417    RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
    417418    *mweatheridCol = id;
     419
    418420    RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
    419     id = table_->focus().addEntry(pksrec.focusAxi, pksrec.focusTan,
    420                                   pksrec.focusRot);
     421    id = table_->focus().addEntry(pksrec.parAngle, pksrec.focusAxi,
     422                                  pksrec.focusTan, pksrec.focusRot);
    421423    *mfocusidCol = id;
    422424    RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
     
    426428    RecordFieldPtr<Float> elCol(rec, "ELEVATION");
    427429    *elCol = pksrec.elevation;
    428 
    429     RecordFieldPtr<Float> parCol(rec, "PARANGLE");
    430     *parCol = pksrec.parAngle;
    431430
    432431    RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
     
    542541  }
    543542
    544   // set frame keyword of FREQUENCIES table
     543  // set FRAME and BASEFRAME keyword of FREQUENCIES table
    545544  if ( header_->freqref != "TOPO" ) {
    546545    table_->frequencies().setFrame( header_->freqref, false ) ;
     546    table_->frequencies().setFrame( header_->freqref, true ) ;
    547547  }
    548548
  • branches/alma/src/STFocus.cpp

    r957 r1757  
    3434}
    3535
    36 asap::STFocus::STFocus( casa::Table tab ) : STSubTable(tab, name_)
     36STFocus::STFocus( casa::Table tab ) :
     37  STSubTable(tab, name_)
    3738{
     39  parangleCol_.attach(table_,"PARANGLE");
    3840  rotationCol_.attach(table_,"ROTATION");
    3941  axisCol_.attach(table_,"AXIS");
     
    5052}
    5153
    52 STFocus& asap::STFocus::operator =( const STFocus & other )
     54STFocus& STFocus::operator =( const STFocus & other )
    5355{
    5456  if (this != &other) {
    5557    static_cast<STSubTable&>(*this) = other;
     58    parangleCol_.attach(table_,"PARANGLE");
    5659    rotationCol_.attach(table_,"ROTATION");
    5760    axisCol_.attach(table_,"AXIS");
     
    6568  return *this;
    6669}
    67 void asap::STFocus::setup( )
     70void STFocus::setup( )
    6871{
    6972  // add to base class table
     73  table_.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
    7074  table_.addColumn(ScalarColumnDesc<Float>("ROTATION"));
    7175  table_.addColumn(ScalarColumnDesc<Float>("AXIS"));
     
    7680  table_.addColumn(ScalarColumnDesc<Float>("XYPHASE"));
    7781  table_.addColumn(ScalarColumnDesc<Float>("XYPHASEOFFSET"));
     82  table_.rwKeywordSet().define("PARALLACTIFY", False);
    7883
    7984  // new cached columns
     85  parangleCol_.attach(table_,"PARANGLE");
    8086  rotationCol_.attach(table_,"ROTATION");
    8187  axisCol_.attach(table_,"AXIS");
     
    8894}
    8995
    90 uInt STFocus::addEntry( Float fax, Float ftan, Float frot, Float hand,
    91                         Float user, Float mount,
    92                         Float xyphase, Float xyphaseoffset)
     96  uInt STFocus::addEntry( Float pa, Float fax, Float ftan, Float frot, Float hand,
     97                          Float user, Float mount,
     98                          Float xyphase, Float xyphaseoffset)
    9399{
    94   Table result = table_( near(table_.col("ROTATION"), frot)
    95                     && near(table_.col("AXIS"), fax)
    96                     && near(table_.col("TAN"), ftan)
    97                     && near(table_.col("HAND"), hand)
    98                     && near(table_.col("USERPHASE"), user)
    99                     && near(table_.col("MOUNT"), mount)
    100                     && near(table_.col("XYPHASE"), xyphase)
    101                     && near(table_.col("XYPHASEOFFSET"), xyphaseoffset)
    102                     );
     100  Table result = table_(  near(table_.col("PARANGLE"), pa)
     101                          && near(table_.col("ROTATION"), frot)
     102                          && near(table_.col("AXIS"), fax)
     103                          && near(table_.col("TAN"), ftan)
     104                          && near(table_.col("HAND"), hand)
     105                          && near(table_.col("USERPHASE"), user)
     106                          && near(table_.col("MOUNT"), mount)
     107                          && near(table_.col("XYPHASE"), xyphase)
     108                          && near(table_.col("XYPHASEOFFSET"), xyphaseoffset)
     109                          );
    103110  uInt resultid = 0;
    104111  if ( result.nrow() > 0) {
     
    113120      resultid++;
    114121    }
     122    parangleCol_.put(rno, pa);
    115123    rotationCol_.put(rno, frot);
    116124    axisCol_.put(rno, fax);
     
    126134}
    127135
    128 void asap::STFocus::getEntry( Float& rotation, Float& angle, Float& ftan,
    129                               Float& hand, Float& user, Float& mount,
    130                               Float& xyphase, Float& xyphaseoffset,
    131                               uInt id) const
     136  void STFocus::getEntry( Float& pa, Float& rotation, Float& angle, Float& ftan,
     137                                Float& hand, Float& user, Float& mount,
     138                                Float& xyphase, Float& xyphaseoffset,
     139                                uInt id) const
    132140{
    133141  Table t = table_(table_.col("ID") == Int(id) );
     
    138146  // get first row - there should only be one matching id
    139147  const TableRecord& rec = row.get(0);
     148  pa = rec.asFloat("PARANGLE");
    140149  rotation = rec.asFloat("ROTATION");
    141150  angle = rec.asFloat("AXIS");
     
    149158
    150159
    151 casa::Float asap::STFocus::getTotalFeedAngle( casa::uInt id ) const
     160casa::Float STFocus::getTotalAngle( casa::uInt id ) const
    152161{
    153162  Float total = 0.0f;
    154163  Table t = table_(table_.col("ID") == Int(id) );
    155164  if (t.nrow() == 0 ) {
    156     throw(AipsError("STFocus::getEntry - id out of range"));
     165    throw(AipsError("STFocus::getTotalAngle - id out of range"));
     166  }
     167  if (table_.keywordSet().asBool("PARALLACTIFY")) {
     168    return 0.0f;
    157169  }
    158170  ROTableRow row(t);
    159171  // get first row - there should only be one matching id
    160172  const TableRecord& rec = row.get(0);
     173  total += rec.asFloat("PARANGLE"); 
    161174  total += rec.asFloat("ROTATION");
    162175  total += rec.asFloat("USERPHASE");
     
    164177  return total;
    165178}
    166 }
    167179
    168 casa::Float asap::STFocus::getFeedHand( casa::uInt id ) const
     180
     181casa::Float STFocus::getFeedHand( casa::uInt id ) const
    169182{
    170183  Table t = table_(table_.col("ID") == Int(id) );
     
    177190}
    178191
     192void STFocus::setParallactify(bool istrue) {
     193  table_.rwKeywordSet().define("PARALLACTIFY", Bool(istrue));
     194}
     195
     196}
  • branches/alma/src/STFocus.h

    r1353 r1757  
    3737  STFocus& operator=(const STFocus& other);
    3838
    39   casa::uInt addEntry( casa::Float faxis, casa::Float ftan,
     39  casa::uInt addEntry( casa::Float pa, casa::Float faxis, casa::Float ftan,
    4040                       casa::Float frot, casa::Float hand=1.0f,
    4141                       casa::Float mount=0.0f, casa::Float user=0.0f,
    42                        casa::Float xyphase=0.0f, casa::Float xyphaseoffset=0.0f);
     42                       casa::Float xyphase=0.0f,
     43                       casa::Float xyphaseoffset=0.0f);
    4344
    44   void getEntry( casa::Float& fax, casa::Float& ftan,
     45  void getEntry( casa::Float& pa, casa::Float& fax, casa::Float& ftan,
    4546                 casa::Float& frot, casa::Float& hand,
    4647                 casa::Float& mount, casa::Float& user,
     
    4849                 casa::uInt id) const;
    4950
    50   casa::Float getTotalFeedAngle(casa::uInt id) const;
     51  casa::Float getTotalAngle(casa::uInt id) const;
     52
     53  casa::Float getParAngle(casa::uInt id) const {
     54    return parangleCol_(id);
     55  }
    5156  casa::Float getFeedHand(casa::uInt id) const;
     57
     58  void setParallactify(bool istrue=false);
    5259
    5360  const casa::String& name() const { return name_; }
     
    5764  static const casa::String name_;
    5865  casa::ScalarColumn<casa::Float> rotationCol_, axisCol_,
    59                                   tanCol_,handCol_,
    60                                   mountCol_,userCol_,
    61                                   xyphCol_,xyphoffCol_;
     66    tanCol_,handCol_, parangleCol_,
     67    mountCol_,userCol_, xyphCol_,xyphoffCol_,;
    6268};
    6369
  • branches/alma/src/STFrequencies.cpp

    r1603 r1757  
    170170**/
    171171SpectralCoordinate
    172   asap::STFrequencies::getSpectralCoordinate( const MDirection& md,
     172  STFrequencies::getSpectralCoordinate( const MDirection& md,
    173173                                              const MPosition& mp,
    174174                                              const MEpoch& me,
     
    242242{
    243243  Vector<Float> offset(1,0.0);
    244   Vector<Float> factors(1,1.0/width);
     244  Vector<Float> factors(1,width);
    245245  Vector<Int> newshape;
    246246  CoordinateSystem csys;
     
    317317        << rec.asDouble("REFVAL") << setw(7)
    318318        << rec.asDouble("REFPIX")
    319         << setw(12)
     319        << setw(15)
    320320        << rec.asDouble("INCREMENT");
    321321  }
     
    324324    int f = outstr.find_first_not_of(' ');
    325325    int l = outstr.find_last_not_of(' ', outstr.size());
    326     if (f < 0) { 
    327       f = 0; 
    328     }
    329     if ( l < f  || l < f ) { 
     326    if (f < 0) {
     327      f = 0;
     328    }
     329    if ( l < f  || l < f ) {
    330330      l = outstr.size();
    331331    }
     
    439439}
    440440
    441 void STFrequencies::shiftRefPix(int npix, uInt id) 
     441void STFrequencies::shiftRefPix(int npix, uInt id)
    442442{
    443443  Table t = table_(table_.col("ID") == Int(id) );
  • branches/alma/src/STLineFinder.cpp

    r1603 r1757  
    3434#include "STLineFinder.h"
    3535#include "STFitter.h"
     36#include "IndexedCompare.h"
    3637
    3738// STL
     
    110111
    111112protected:
    112    // supplementary function to control running mean calculations.
    113    // It adds a specified channel to the running mean box and
     113   // supplementary function to control running mean/median calculations.
     114   // It adds a specified channel to the running box and
    114115   // removes (ch-maxboxnchan+1)'th channel from there
    115116   // Channels, for which the mask is false or index is beyond the
     
    152153                                           // last point of the detected line
    153154                                           //
     155   bool itsUseMedian;                      // true if median statistics is used
     156                                           // to determine the noise level, otherwise
     157                                           // it is the mean of the lowest 80% of deviations
     158                                           // (default)
     159   int itsNoiseSampleSize;                 // sample size used to estimate the noise statistics
     160                                           // Negative value means the whole spectrum is used (default)
    154161public:
    155162
     
    157164   LFAboveThreshold(std::list<pair<int,int> > &in_lines,
    158165                    int in_min_nchan = 3,
    159                     casa::Float in_threshold = 5) throw();
     166                    casa::Float in_threshold = 5,
     167                    bool use_median = false,
     168                    int noise_sample_size = -1) throw();
    160169   virtual ~LFAboveThreshold() throw();
    161170
     
    200209///////////////////////////////////////////////////////////////////////////////
    201210
     211///////////////////////////////////////////////////////////////////////////////
     212//
     213// LFNoiseEstimator  a helper class designed to estimate off-line variance
     214//                   using statistics depending on the distribution of
     215//                   values (e.g. like a median)
     216//
     217//                   Two statistics are supported: median and an average of
     218//                   80% of smallest values.
     219//
     220
     221struct LFNoiseEstimator {
     222   // construct an object
     223   // size - maximum sample size. After a size number of elements is processed
     224   // any new samples would cause the algorithm to drop the oldest samples in the
     225   // buffer.
     226   explicit LFNoiseEstimator(size_t size);
     227
     228   // add a new sample
     229   // in - the new value
     230   void add(float in);
     231
     232   // median of the distribution
     233   float median() const;
     234
     235   // mean of lowest 80% of the samples
     236   float meanLowest80Percent() const;
     237
     238   // return true if the buffer is full (i.e. statistics are representative)
     239   inline bool filledToCapacity() const { return itsBufferFull;}
     240
     241protected:
     242   // update cache of sorted indices
     243   // (it is assumed that itsSampleNumber points to the newly
     244   // replaced element)
     245   void updateSortedCache() const;
     246
     247   // build sorted cache from the scratch
     248   void buildSortedCache() const;
     249
     250   // number of samples accumulated so far
     251   // (can be less than the buffer size)
     252   size_t numberOfSamples() const;
     253
     254   // this helper method builds the cache if
     255   // necessary using one of the methods
     256   void fillCacheIfNecessary() const;
     257
     258private:
     259   // buffer with samples (unsorted)
     260   std::vector<float> itsVariances;
     261   // current sample number (<=itsVariances.size())
     262   size_t itsSampleNumber;
     263   // true, if the buffer all values in the sample buffer are used
     264   bool itsBufferFull;
     265   // cached indices into vector of samples
     266   mutable std::vector<size_t> itsSortedIndices;
     267   // true if any of the statistics have been obtained at least
     268   // once. This flag allows to implement a more efficient way of
     269   // calculating statistics, if they are needed at once and not
     270   // after each addition of a new element
     271   mutable bool itsStatisticsAccessed;
     272};
     273
     274//
     275///////////////////////////////////////////////////////////////////////////////
     276
     277
    202278} // namespace asap
     279
     280///////////////////////////////////////////////////////////////////////////////
     281//
     282// LFNoiseEstimator  a helper class designed to estimate off-line variance
     283//                   using statistics depending on the distribution of
     284//                   values (e.g. like a median)
     285//
     286//                   Two statistics are supported: median and an average of
     287//                   80% of smallest values.
     288//
     289
     290// construct an object
     291// size - maximum sample size. After a size number of elements is processed
     292// any new samples would cause the algorithm to drop the oldest samples in the
     293// buffer.
     294LFNoiseEstimator::LFNoiseEstimator(size_t size) : itsVariances(size),
     295     itsSampleNumber(0), itsBufferFull(false), itsSortedIndices(size),
     296     itsStatisticsAccessed(false)
     297{
     298   AlwaysAssert(size>0,AipsError);
     299}
     300
     301
     302// add a new sample
     303// in - the new value
     304void LFNoiseEstimator::add(float in)
     305{
     306   if (isnan(in)) {
     307       // normally it shouldn't happen
     308       return;
     309   }
     310   itsVariances[itsSampleNumber] = in;
     311
     312   if (itsStatisticsAccessed) {
     313       // only do element by element addition if on-the-fly
     314       // statistics are needed
     315       updateSortedCache();
     316   }
     317
     318   // advance itsSampleNumber now
     319   ++itsSampleNumber;
     320   if (itsSampleNumber == itsVariances.size()) {
     321       itsSampleNumber = 0;
     322       itsBufferFull = true;
     323   }
     324   AlwaysAssert(itsSampleNumber<itsVariances.size(),AipsError);
     325}
     326
     327// number of samples accumulated so far
     328// (can be less than the buffer size)
     329size_t LFNoiseEstimator::numberOfSamples() const
     330{
     331  // the number of samples accumulated so far may be less than the
     332  // buffer size
     333  const size_t nSamples = itsBufferFull ? itsVariances.size(): itsSampleNumber;
     334  AlwaysAssert( (nSamples > 0) && (nSamples <= itsVariances.size()), AipsError);
     335  return nSamples;
     336}
     337
     338// this helper method builds the cache if
     339// necessary using one of the methods
     340void LFNoiseEstimator::fillCacheIfNecessary() const
     341{
     342  if (!itsStatisticsAccessed) {
     343      if ((itsSampleNumber!=0) || itsBufferFull) {
     344          // build the whole cache efficiently
     345          buildSortedCache();
     346      } else {
     347          updateSortedCache();
     348      }
     349      itsStatisticsAccessed = true;
     350  } // otherwise, it is updated in 'add' using on-the-fly method
     351}
     352
     353// median of the distribution
     354float LFNoiseEstimator::median() const
     355{
     356  fillCacheIfNecessary();
     357  // the number of samples accumulated so far may be less than the
     358  // buffer size
     359  const size_t nSamples = numberOfSamples();
     360  const size_t medSample = nSamples / 2;
     361  AlwaysAssert(medSample < itsSortedIndices.size(), AipsError);
     362  return itsVariances[itsSortedIndices[medSample]];
     363}
     364
     365// mean of lowest 80% of the samples
     366float LFNoiseEstimator::meanLowest80Percent() const
     367{
     368  fillCacheIfNecessary();
     369  // the number of samples accumulated so far may be less than the
     370  // buffer size
     371  const size_t nSamples = numberOfSamples();
     372  float result = 0;
     373  size_t numpt=size_t(0.8*nSamples);
     374  if (!numpt) {
     375      numpt=nSamples; // no much else left,
     376                     // although it is very inaccurate
     377  }
     378  AlwaysAssert( (numpt > 0) && (numpt<itsSortedIndices.size()), AipsError);
     379  for (size_t ch=0; ch<numpt; ++ch) {
     380       result += itsVariances[itsSortedIndices[ch]];
     381  }
     382  result /= float(numpt);
     383  return result;
     384}
     385
     386// update cache of sorted indices
     387// (it is assumed that itsSampleNumber points to the newly
     388// replaced element)
     389void LFNoiseEstimator::updateSortedCache() const
     390{
     391  // the number of samples accumulated so far may be less than the
     392  // buffer size
     393  const size_t nSamples = numberOfSamples();
     394
     395  if (itsBufferFull) {
     396      // first find the index of the element which is being replaced
     397      size_t index = nSamples;
     398      for (size_t i=0; i<nSamples; ++i) {
     399           AlwaysAssert(i < itsSortedIndices.size(), AipsError);
     400           if (itsSortedIndices[i] == itsSampleNumber) {
     401               index = i;
     402               break;
     403           }
     404      }
     405      AlwaysAssert( index < nSamples, AipsError);
     406
     407      const vector<size_t>::iterator indStart = itsSortedIndices.begin();
     408      // merge this element with preceeding block first
     409      if (index != 0) {
     410          // merge indices on the basis of variances
     411          inplace_merge(indStart,indStart+index,indStart+index+1,
     412                        indexedCompare<size_t>(itsVariances.begin()));
     413      }
     414      // merge with the following block
     415      if (index + 1 != nSamples) {
     416          // merge indices on the basis of variances
     417          inplace_merge(indStart,indStart+index+1,indStart+nSamples,
     418                        indexedCompare<size_t>(itsVariances.begin()));
     419      }
     420  } else {
     421      // itsSampleNumber is the index of the new element
     422      AlwaysAssert(itsSampleNumber < itsSortedIndices.size(), AipsError);
     423      itsSortedIndices[itsSampleNumber] = itsSampleNumber;
     424      if (itsSampleNumber >= 1) {
     425          // we have to place this new sample in
     426          const vector<size_t>::iterator indStart = itsSortedIndices.begin();
     427          // merge indices on the basis of variances
     428          inplace_merge(indStart,indStart+itsSampleNumber,indStart+itsSampleNumber+1,
     429                        indexedCompare<size_t>(itsVariances.begin()));
     430      }
     431  }
     432}
     433
     434// build sorted cache from the scratch
     435void LFNoiseEstimator::buildSortedCache() const
     436{
     437  // the number of samples accumulated so far may be less than the
     438  // buffer size
     439  const size_t nSamples = numberOfSamples();
     440  AlwaysAssert(nSamples <= itsSortedIndices.size(), AipsError);
     441  for (size_t i=0; i<nSamples; ++i) {
     442       itsSortedIndices[i]=i;
     443  }
     444
     445  // sort indices, but check the array of variances
     446  const vector<size_t>::iterator indStart = itsSortedIndices.begin();
     447  stable_sort(indStart,indStart+nSamples, indexedCompare<size_t>(itsVariances.begin()));
     448}
     449
     450//
     451///////////////////////////////////////////////////////////////////////////////
    203452
    204453///////////////////////////////////////////////////////////////////////////////
     
    275524}
    276525
    277 // supplementary function to control running mean calculations.
    278 // It adds a specified channel to the running mean box and
     526// supplementary function to control running mean/median calculations.
     527// It adds a specified channel to the running box and
    279528// removes (ch-max_box_nchan+1)'th channel from there
    280529// Channels, for which the mask is false or index is beyond the
     
    339588                (meanch2-square(meanch));
    340589      linmean=coeff*(Float(cur_channel)-meanch)+mean;
    341       linvariance=sqrt(sumf2/Float(box_chan_cntr)-square(mean)-
    342                     square(coeff)*(meanch2-square(meanch)));
     590      linvariance=sumf2/Float(box_chan_cntr)-square(mean)-
     591                    square(coeff)*(meanch2-square(meanch));
     592      if (linvariance<0.) {
     593          // this shouldn't happen normally, but could be due to round-off error
     594          linvariance = 0;
     595      } else {
     596          linvariance = sqrt(linvariance);
     597      }
    343598  }
    344599  need2recalculate=False;
     
    351606///////////////////////////////////////////////////////////////////////////////
    352607//
    353 // LFAboveThreshold - a running mean algorithm for line detection
     608// LFAboveThreshold - a running mean/median algorithm for line detection
    354609//
    355610//
     
    359614LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines,
    360615                                   int in_min_nchan,
    361                                    casa::Float in_threshold) throw() :
     616                                   casa::Float in_threshold,
     617                                   bool use_median,
     618                                   int noise_sample_size) throw() :
    362619             min_nchan(in_min_nchan), threshold(in_threshold),
    363              lines(in_lines), running_box(NULL) {}
     620             lines(in_lines), running_box(NULL), itsUseMedian(use_median),
     621             itsNoiseSampleSize(noise_sample_size) {}
    364622
    365623LFAboveThreshold::~LFAboveThreshold() throw()
     
    474732      running_box=new RunningBox(spectrum,mask,edge,max_box_nchan);
    475733
    476 
    477734      // determine the off-line variance first
    478735      // an assumption made: lines occupy a small part of the spectrum
    479736
    480       std::vector<float> variances(edge.second-edge.first);
    481       DebugAssert(variances.size(),AipsError);
    482 
    483       for (;running_box->haveMore();running_box->next())
    484            variances[running_box->getChannel()-edge.first]=
    485                                 running_box->getLinVariance();
    486 
    487       // in the future we probably should do a proper Chi^2 estimation
    488       // now a simple 80% of smaller values will be used.
    489       // it may degrade the performance of the algorithm for weak lines
    490       // due to a bias of the Chi^2 distribution.
    491       stable_sort(variances.begin(),variances.end());
    492 
    493       Float offline_variance=0;
    494       uInt offline_cnt=uInt(0.8*variances.size());
    495       if (!offline_cnt) offline_cnt=variances.size(); // no much else left,
    496                                     // although it is very inaccurate
    497       for (uInt n=0;n<offline_cnt;++n)
    498            offline_variance+=variances[n];
    499       offline_variance/=Float(offline_cnt);
     737      const size_t noiseSampleSize = itsNoiseSampleSize<0 ? size_t(edge.second-edge.first) :
     738                      std::min(size_t(itsNoiseSampleSize), size_t(edge.second-edge.first));
     739      DebugAssert(noiseSampleSize,AipsError);
     740      const bool globalNoise = (size_t(edge.second - edge.first) == noiseSampleSize);
     741      LFNoiseEstimator ne(noiseSampleSize);
     742
     743      for (;running_box->haveMore();running_box->next()) {
     744           ne.add(running_box->getLinVariance());
     745           if (ne.filledToCapacity()) {
     746               break;
     747           }
     748      }
     749
     750      Float offline_variance = -1; // just a flag that it is unset
     751           
     752      if (globalNoise) {
     753          offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
     754      }
    500755
    501756      // actual search algorithm
     
    510765                                 running_box->next()) {
    511766           const int ch=running_box->getChannel();
    512            if (running_box->getNumberOfBoxPoints()>=minboxnchan)
     767           if (!globalNoise) {
     768               // add a next point for a local noise estimate
     769               ne.add(running_box->getLinVariance());
     770           }   
     771           if (running_box->getNumberOfBoxPoints()>=minboxnchan) {
     772               if (!globalNoise) {
     773                   offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
     774               }
     775               AlwaysAssert(offline_variance>0.,AipsError);
    513776               processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
    514777                  threshold*offline_variance), mask);
    515            else processCurLine(mask); // just finish what was accumulated before
     778           } else processCurLine(mask); // just finish what was accumulated before
    516779
    517780           signs[ch]=getAboveMeanSign();
    518            // os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<<
    519            // threshold*offline_variance<<endl;
    520 
    521            const Float buf=running_box->aboveMean();
    522            if (buf>0) signs[ch]=1;
    523            else if (buf<0) signs[ch]=-1;
    524            else if (buf==0) signs[ch]=0;
    525            //os<<ch<<" "<<spectrum[ch]<<" "<<running_box->getLinMean()<<" "<<
    526            //             threshold*offline_variance<<endl;
     781            //os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<<
     782            //threshold*offline_variance<<endl;
    527783      }
    528784      if (lines.size())
     
    649905//              Setting a very large value doesn't usually provide
    650906//              valid detections.
    651 // in_box_size  the box size for running mean calculation. Default is
     907// in_box_size  the box size for running mean/median calculation. Default is
    652908//              1./5. of the whole spectrum size
     909// in_noise_box the box size for off-line noise estimation (if working with
     910//              local noise. Negative value means use global noise estimate
     911//              Default is -1 (i.e. estimate using the whole spectrum)
     912// in_median    true if median statistics is used as opposed to average of
     913//              the lowest 80% of deviations (default)
    653914void STLineFinder::setOptions(const casa::Float &in_threshold,
    654915                              const casa::Int &in_min_nchan,
    655916                              const casa::Int &in_avg_limit,
    656                               const casa::Float &in_box_size) throw()
     917                              const casa::Float &in_box_size,
     918                              const casa::Float &in_noise_box,
     919                              const casa::Bool &in_median) throw()
    657920{
    658921  threshold=in_threshold;
     
    660923  avg_limit=in_avg_limit;
    661924  box_size=in_box_size;
     925  itsNoiseBox = in_noise_box;
     926  itsUseMedian = in_median;
    662927}
    663928
     
    683948                const casa::uInt &whichRow) throw(casa::AipsError)
    684949{
    685   //const int minboxnchan=4;
    686950  if (scan.null())
    687951      throw AipsError("STLineFinder::findLines - a scan should be set first,"
     
    700964      throw AipsError("STLineFinder::findLines - in_scan and in_mask have different"
    701965            "number of spectral channels.");
     966
     967  // taking flagged channels into account
     968  vector<bool> flaggedChannels = scan->getMask(whichRow);
     969  if (flaggedChannels.size()) {
     970      // there is a mask set for this row
     971      if (flaggedChannels.size() != mask.nelements()) {
     972          throw AipsError("STLineFinder::findLines - internal inconsistency: number of mask elements do not match the number of channels");
     973      }
     974      for (size_t ch = 0; ch<mask.nelements(); ++ch) {
     975           mask[ch] &= flaggedChannels[ch];
     976      }
     977  }
     978
    702979  // number of elements in in_edge
    703980  if (in_edge.size()>2)
     
    7321009      throw AipsError("STLineFinder::findLines - box_size is too small");
    7331010
     1011  // number of elements in the sample for noise estimate
     1012  const int noise_box = itsNoiseBox<0 ? -1 : int(nchan * itsNoiseBox);
     1013
     1014  if ((noise_box!= -1) and (noise_box<2))
     1015      throw AipsError("STLineFinder::findLines - noise_box is supposed to be at least 2 elements");
     1016
    7341017  spectrum.resize();
    7351018  spectrum = Vector<Float>(scan->getSpectrum(whichRow));
     
    7551038     try {
    7561039         // line find algorithm
    757          LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold);
     1040         LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold, itsUseMedian,noise_box);
    7581041         lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
    7591042         signs.resize(lfalg.getSigns().nelements());
  • branches/alma/src/STLineFinder.h

    r1353 r1757  
    150150   // in_box_size  the box size for running mean calculation. Default is
    151151   //              1./5. of the whole spectrum size
     152   // in_noise_box the box size for off-line noise estimation (if working with
     153   //              local noise. Negative value means use global noise estimate
     154   //              Default is -1 (i.e. estimate using the whole spectrum)
     155   // in_median    true if median statistics is used as opposed to average of
     156   //              the lowest 80% of deviations (default)
    152157   void setOptions(const casa::Float &in_threshold=sqrt(3.),
    153158                   const casa::Int &in_min_nchan=3,
    154159                   const casa::Int &in_avg_limit=8,
    155                    const casa::Float &in_box_size=0.2) throw();
     160                   const casa::Float &in_box_size=0.2,
     161                   const casa::Float &in_noise_box=-1.,
     162                   const casa::Bool &in_median = casa::False) throw();
    156163
    157164   // set the scan to work with (in_scan parameter)
     
    241248   // a buffer for the spectrum
    242249   mutable casa::Vector<casa::Float>  spectrum;
     250
     251   // the box size for off-line noise estimation (if working with
     252   // local noise. Negative value means use global noise estimate
     253   // Default is -1 (i.e. estimate using the whole spectrum)
     254   casa::Float itsNoiseBox;
     255
     256   // true if median statistics is used as opposed to average of
     257   // the lowest 80% of deviations (default)
     258   casa::Bool itsUseMedian;
    243259};
    244260
  • branches/alma/src/STMath.cpp

    r1719 r1757  
    5252#include "RowAccumulator.h"
    5353#include "STAttr.h"
     54#include "STSelector.h"
     55
    5456#include "STMath.h"
    55 #include "STSelector.h"
    56 
    5757using namespace casa;
    5858
     
    740740}
    741741
    742 CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left, 
    743                                             const CountedPtr<Scantable>& right, 
     742CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
     743                                            const CountedPtr<Scantable>& right,
    744744                                            const std::string& mode)
    745745{
     
    13791379{
    13801380
    1381  
     1381
    13821382  STSelector sel;
    13831383  CountedPtr< Scantable > ws = getScantable(s, false);
     
    15181518  // for each row
    15191519  // assume that the data order are same between sig and ref
    1520   RowAccumulator acc( asap::TINTSYS ) ;
     1520  RowAccumulator acc( asap::W_TINTSYS ) ;
    15211521  for ( int i = 0 ; i < sig->nrow() ; i++ ) {
    15221522    // get values
     
    19521952  // initialize the lookup table if necessary
    19531953  if ( lookup.empty() ) {
    1954     lookup["NONE"]   = asap::NONE;
    1955     lookup["TINT"] = asap::TINT;
    1956     lookup["TINTSYS"]  = asap::TINTSYS;
    1957     lookup["TSYS"]  = asap::TSYS;
    1958     lookup["VAR"]  = asap::VAR;
     1954    lookup["NONE"]   = asap::W_NONE;
     1955    lookup["TINT"] = asap::W_TINT;
     1956    lookup["TINTSYS"]  = asap::W_TINTSYS;
     1957    lookup["TSYS"]  = asap::W_TSYS;
     1958    lookup["VAR"]  = asap::W_VAR;
    19591959  }
    19601960
     
    20002000    String msg;
    20012001    if ( nc > 0 ) {
    2002       ppoly = new Polynomial<Float>(nc);
     2002      ppoly = new Polynomial<Float>(nc-1);
    20032003      coeff = coeffs;
    20042004      msg = String("user");
     
    20062006      STAttr sdAttr;
    20072007      coeff = sdAttr.gainElevationPoly(inst);
    2008       ppoly = new Polynomial<Float>(3);
     2008      ppoly = new Polynomial<Float>(coeff.nelements()-1);
    20092009      msg = String("built in");
    20102010    }
     
    21442144    if (d < 0) {
    21452145      Instrument inst =
    2146         STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), 
     2146        STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
    21472147                                  True);
    21482148      STAttr sda;
     
    22072207
    22082208CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
    2209                                          float tau )
     2209                                         const std::vector<float>& tau )
    22102210{
    22112211  CountedPtr< Scantable > out = getScantable(in, false);
    22122212
    2213   Table tab = out->table();
    2214   ROScalarColumn<Float> elev(tab, "ELEVATION");
    2215   ArrayColumn<Float> specCol(tab, "SPECTRA");
    2216   ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
    2217   ArrayColumn<Float> tsysCol(tab, "TSYS");
    2218   for ( uInt i=0; i<tab.nrow(); ++i) {
    2219     Float zdist = Float(C::pi_2) - elev(i);
    2220     Float factor = exp(tau/cos(zdist));
    2221     MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
    2222     ma *= factor;
    2223     specCol.put(i, ma.getArray());
    2224     flagCol.put(i, flagsFromMA(ma));
    2225     Vector<Float> tsys;
    2226     tsysCol.get(i, tsys);
    2227     tsys *= factor;
    2228     tsysCol.put(i, tsys);
     2213  Table outtab = out->table();
     2214
     2215  const uInt ntau = uInt(tau.size());
     2216  std::vector<float>::const_iterator tauit = tau.begin();
     2217  AlwaysAssert((ntau == 1 || ntau == in->nif() || ntau == in->nif() * in->npol()),
     2218               AipsError);
     2219  TableIterator iiter(outtab, "IFNO");
     2220  while ( !iiter.pastEnd() ) {
     2221    Table itab = iiter.table();
     2222    TableIterator piter(outtab, "POLNO");
     2223    while ( !piter.pastEnd() ) {
     2224      Table tab = piter.table();
     2225      ROScalarColumn<Float> elev(tab, "ELEVATION");
     2226      ArrayColumn<Float> specCol(tab, "SPECTRA");
     2227      ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
     2228      ArrayColumn<Float> tsysCol(tab, "TSYS");
     2229      for ( uInt i=0; i<tab.nrow(); ++i) {
     2230        Float zdist = Float(C::pi_2) - elev(i);
     2231        Float factor = exp(*tauit/cos(zdist));
     2232        MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
     2233        ma *= factor;
     2234        specCol.put(i, ma.getArray());
     2235        flagCol.put(i, flagsFromMA(ma));
     2236        Vector<Float> tsys;
     2237        tsysCol.get(i, tsys);
     2238        tsys *= factor;
     2239        tsysCol.put(i, tsys);
     2240      }
     2241      if (ntau == in->nif()*in->npol() ) {
     2242        tauit++;
     2243      }
     2244      piter++;
     2245    }
     2246    if (ntau >= in->nif() ) {
     2247      tauit++;
     2248    }
     2249    iiter++;
    22292250  }
    22302251  return out;
     
    22332254CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
    22342255                                             const std::string& kernel,
    2235                                              float width )
     2256                                             float width, int order)
    22362257{
    22372258  CountedPtr< Scantable > out = getScantable(in, false);
     
    22542275      mathutil::runningMedian(specout, maskout, spec , mask, width);
    22552276      convertArray(flag, maskout);
     2277    } else if ( kernel == "poly" ) {
     2278      mathutil::polyfit(specout, maskout, spec, !mask, width, order);
     2279      convertArray(flag, !maskout);
    22562280    }
    22572281    flagCol.put(i, flag);
     
    22622286
    22632287CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
    2264                                         const std::string& kernel, float width )
    2265 {
    2266   if (kernel == "rmedian"  || kernel == "hanning") {
    2267     return smoothOther(in, kernel, width);
     2288                                        const std::string& kernel, float width,
     2289                                        int order)
     2290{
     2291  if (kernel == "rmedian"  || kernel == "hanning" || kernel == "poly") {
     2292    return smoothOther(in, kernel, width, order);
    22682293  }
    22692294  CountedPtr< Scantable > out = getScantable(in, false);
     
    23512376          id = out->molecules().addEntry(rf, name, fname);
    23522377          molidcol.put(k, id);
    2353           Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
    2354           (*it)->focus().getEntry(fax, ftan, frot, fhand,
     2378          Float fpa,frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
     2379          (*it)->focus().getEntry(fpa, fax, ftan, frot, fhand,
    23552380                                  fmount,fuser, fxy, fxyp,
    23562381                                  rec.asuInt("FOCUS_ID"));
    2357           id = out->focus().addEntry(fax, ftan, frot, fhand,
     2382          id = out->focus().addEntry(fpa, fax, ftan, frot, fhand,
    23582383                                     fmount,fuser, fxy, fxyp);
    23592384          focusidcol.put(k, id);
     
    23992424  cols[3] = String("CYCLENO");
    24002425  TableIterator iter(tout, cols);
    2401   CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_, 
     2426  CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
    24022427                                               out->getPolType() );
    24032428  while (!iter.pastEnd()) {
     
    24052430    ArrayColumn<Float> speccol(t, "SPECTRA");
    24062431    ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
    2407     ScalarColumn<Float> parancol(t, "PARANGLE");
    24082432    Matrix<Float> pols(speccol.getColumn());
    24092433    try {
    24102434      stpol->setSpectra(pols);
    2411       Float fang,fhand,parang;
    2412       fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
     2435      Float fang,fhand;
     2436      fang = in->focusTable_.getTotalAngle(focidcol(0));
    24132437      fhand = in->focusTable_.getFeedHand(focidcol(0));
    2414       parang = parancol(0);
    2415       /// @todo re-enable this
    2416       // disable total feed angle to support paralactifying Caswell style
    2417       stpol->setPhaseCorrections(parang, -parang, fhand);
     2438      stpol->setPhaseCorrections(fang, fhand);
    24182439      // use a member function pointer in STPol.  This only works on
    24192440      // the STPol pointer itself, not the Counted Pointer so
     
    26812702      uInt row = tab.rowNumbers()[0];
    26822703      stpol->setSpectra(in->getPolMatrix(row));
    2683       Float fang,fhand,parang;
    2684       fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
     2704      Float fang,fhand;
     2705      fang = in->focusTable_.getTotalAngle(in->mfocusidCol_(row));
    26852706      fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
    2686       parang = in->paraCol_(row);
    2687       /// @todo re-enable this
    2688       // disable total feed angle to support paralactifying Caswell style
    2689       stpol->setPhaseCorrections(parang, -parang, fhand);
     2707      stpol->setPhaseCorrections(fang, fhand);
    26902708      Int npolout = 0;
    26912709      for (uInt i=0; i<tab.nrow(); ++i) {
     
    27372755CountedPtr< Scantable >
    27382756  asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
    2739                           double frequency, double width )
     2757                         double start, double end,
     2758                         const std::string& mode)
    27402759{
    27412760  CountedPtr< Scantable > out = getScantable(in, false);
     
    27552774      Vector<Float> spec = specCol(i);
    27562775      Vector<uChar> flag = flagCol(i);
    2757       Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
    2758       Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
     2776      int fstart = -1;
     2777      int fend = -1;
    27592778      for (unsigned int k=0; k < flag.nelements(); ++k ) {
    27602779        if (flag[k] > 0) {
    2761           spec[k] = 0.0;
     2780          fstart = k;
     2781          while (flag[k] > 0 && k < flag.nelements()) {
     2782            fend = k;
     2783            k++;
     2784          }
    27622785        }
     2786        Float interp = 0.0;
     2787        if (fstart-1 > 0 ) {
     2788          interp = spec[fstart-1];
     2789          if (fend+1 < spec.nelements()) {
     2790            interp = (interp+spec[fend+1])/2.0;
     2791          }
     2792        } else {
     2793          interp = spec[fend+1];
     2794        }
     2795        if (fstart > -1 && fend > -1) {
     2796          for (int j=fstart;j<=fend;++j) {
     2797            spec[j] = interp;
     2798          }
     2799        }
     2800        fstart =-1;
     2801        fend = -1;
    27632802      }
    27642803      Vector<Complex> lags;
    27652804      ffts.fft0(lags, spec);
    2766       Int start =  max(0, lag0);
    2767       Int end =  min(Int(lags.nelements()-1), lag1);
    2768       if (start == end) {
    2769         lags[start] = Complex(0.0);
     2805      Int lag0(start+0.5);
     2806      Int lag1(end+0.5);
     2807      if (mode == "frequency") {
     2808        lag0 = Int(spec.nelements()*abs(inc)/(start)+0.5);
     2809        lag1 = Int(spec.nelements()*abs(inc)/(end)+0.5);
     2810      }
     2811      Int lstart =  max(0, lag0);
     2812      Int lend =  min(Int(lags.nelements()-1), lag1);
     2813      if (lstart == lend) {
     2814        lags[lstart] = Complex(0.0);
    27702815      } else {
    2771         for (int j=start; j <=end ;++j) {
     2816        if (lstart > lend) {
     2817          Int tmp = lend;
     2818          lend = lstart;
     2819          lstart = tmp;
     2820        }
     2821        for (int j=lstart; j <=lend ;++j) {
    27722822          lags[j] = Complex(0.0);
    27732823        }
  • branches/alma/src/STMath.h

    r1680 r1757  
    146146
    147147  casa::CountedPtr<Scantable>
    148     binaryOperate( const casa::CountedPtr<Scantable>& left, 
    149                    const casa::CountedPtr<Scantable>& right, 
     148    binaryOperate( const casa::CountedPtr<Scantable>& left,
     149                   const casa::CountedPtr<Scantable>& right,
    150150                   const std::string& mode);
    151151
     
    290290  casa::CountedPtr<Scantable>
    291291    smooth(const casa::CountedPtr<Scantable>& in, const std::string& kernel,
    292                       float width);
     292                      float width, int order=2);
    293293
    294294  casa::CountedPtr<Scantable>
     
    302302
    303303  casa::CountedPtr<Scantable> opacity(const casa::CountedPtr<Scantable>& in,
    304                                       float tau);
     304                                      const std::vector<float>& tau);
    305305
    306306  casa::CountedPtr<Scantable>
     
    338338   */
    339339  casa::CountedPtr<Scantable>
    340     lagFlag( const casa::CountedPtr<Scantable>& in, double frequency,
    341               double width);
     340    lagFlag( const casa::CountedPtr<Scantable>& in, double start,
     341             double end, const std::string& mode="frequency");
    342342
    343343  // test for average spectra with different channel/resolution
     
    374374                              bool tokelvin, float cfac);
    375375
    376   casa::CountedPtr< Scantable > 
     376  casa::CountedPtr< Scantable >
    377377    smoothOther( const casa::CountedPtr< Scantable >& in,
    378378                 const std::string& kernel,
    379                  float width );
     379                 float width, int order=2 );
    380380
    381381  casa::CountedPtr< Scantable >
  • branches/alma/src/STMathWrapper.h

    r1680 r1757  
    146146
    147147  ScantableWrapper
    148     smooth(const ScantableWrapper& in, const std::string& kernel, float width)
    149   { return ScantableWrapper(STMath::smooth(in.getCP(), kernel, width)); }
     148    smooth(const ScantableWrapper& in, const std::string& kernel, float width,
     149           int order=2)
     150  { return ScantableWrapper(STMath::smooth(in.getCP(), kernel, width, order)); }
    150151
    151152  ScantableWrapper
     
    164165
    165166  ScantableWrapper opacity(const ScantableWrapper& in,
    166                                       float tau)
     167                          const std::vector<float>& tau)
    167168  { return ScantableWrapper(STMath::opacity(in.getCP(), tau)); }
    168169
     
    201202
    202203  ScantableWrapper lagFlag( const ScantableWrapper& in,
    203                             double frequency, double width )
    204   { return ScantableWrapper(STMath::lagFlag(in.getCP(), frequency, width)); }
     204                            double start, double end,
     205                            const std::string& mode="frequency" )
     206  { return ScantableWrapper(STMath::lagFlag(in.getCP(), start, end,
     207                                            mode)); }
    205208
    206209  // test for average spectra with different channel/resolution
  • branches/alma/src/STPol.h

    r1015 r1757  
    3535
    3636  typedef  void (STPol::*polOperation)( casa::Float phase );
    37   STPol(): totalfeed_(0.0),parangle_(0.0),feedhand_(1.0) {}
     37  STPol(): totalangle_(0.0),feedhand_(1.0) {}
    3838  virtual ~STPol() {}
    3939
     
    8383
    8484
    85   void setPhaseCorrections(casa::Float parangle=0.0, casa::Float totalfeed=0.0,
    86                            casa::Float feedhand=1.0)
    87     { totalfeed_=totalfeed;parangle_=parangle;feedhand_=feedhand;}
     85  void setPhaseCorrections(casa::Float totalang=0.0, casa::Float feedhand=1.0)
     86    { totalangle_=totalang;feedhand_=feedhand;}
    8887
    89   casa::Float getTotalPhase() const { return totalfeed_+parangle_; }
     88  casa::Float getTotalPhase() const { return totalangle_; }
    9089  casa::Float getFeedHand() const { return feedhand_; }
    9190
     
    9998  static std::map<std::string, std::map<int, std::string> > labelmap_;
    10099
    101   casa::Float totalfeed_,parangle_,feedhand_;
     100  casa::Float totalangle_,feedhand_;
    102101  std::string mode_;
    103102  casa::Matrix<casa::Float> basespectra_;
  • branches/alma/src/STWriter.cpp

    r1683 r1757  
    4242#include <atnf/PKSIO/PKSMS2writer.h>
    4343#include <atnf/PKSIO/PKSSDwriter.h>
     44#include <atnf/PKSIO/SrcType.h>
    4445
    4546#include <tables/Tables/Table.h>
     
    5152#include "STAsciiWriter.h"
    5253#include "STHeader.h"
     54#include "STMath.h"
     55
    5356
    5457#include "STWriter.h"
     
    104107                    const std::string &filename)
    105108{
     109  // If we write out foreign formats we have to convert the frequency system
     110  // into the output frame, as we do everything related to SPectarlCoordinates
     111  // in asap on-the-fly.
     112
     113  CountedPtr<Scantable> inst = in;
     114  if (in->frequencies().getFrame(true) != in->frequencies().getFrame(false)) {
     115    STMath stm(false);
     116    inst = stm.frequencyAlign(in);
     117  }
    106118
    107119  if (format_=="ASCII") {
    108120    STAsciiWriter iw;
     121    // ASCII calls SpectralCoordinate::toWorld so no freqAlign use 'in'
    109122    if (iw.write(*in, filename)) {
    110123      return 0;
     
    117130      iw.setClass(True);
    118131    }
    119     if (iw.write(*in, filename)) {
     132    if (iw.write(*inst, filename)) {
    120133      return 0;
    121134    }
     
    127140  // this is a little different from what I have done
    128141  // before. Need to check with the Offline User Test data
    129   STHeader hdr = in->getHeader();
     142  STHeader hdr = inst->getHeader();
    130143  //const Int nPol  = hdr.npol;
    131144  //const Int nChan = hdr.nchan;
    132   std::vector<uint> ifs = in->getIFNos();
    133   int nIF = in->nif();//ifs.size();
     145  std::vector<uint> ifs = inst->getIFNos();
     146  int nIF = inst->nif();//ifs.size();
    134147  Vector<uInt> nPol(nIF),nChan(nIF);
    135148  Vector<Bool> havexpol(nIF);
    136149  String fluxUnit = hdr.fluxunit;
    137 
     150  fluxUnit.upcase();
    138151  nPol = 0;nChan = 0; havexpol = False;
    139152  for (uint i=0;i<ifs.size();++i) {
    140     nPol(ifs[i]) = in->npol();
    141     nChan(ifs[i]) = in->nchan(ifs[i]);
     153    nPol(ifs[i]) = inst->npol();
     154    nChan(ifs[i]) = inst->nchan(ifs[i]);
    142155    havexpol(ifs[i]) = nPol(ifs[i]) > 2;
    143156  }
    144 
    145   const Table table = in->table();
     157//   Vector<String> obstypes(2);
     158//   obstypes(0) = "TR";//on
     159//   obstypes(1) = "RF TR";//off
     160  const Table table = inst->table();
     161
    146162
    147163  // Create the output file and write static data.
     
    155171    throw(AipsError("Failed to create output file"));
    156172  }
    157 
    158173
    159174  Int count = 0;
     
    206221          String tcalt;
    207222          Vector<String> stmp0, stmp1;
    208           in->frequencies().getEntry(crpix,crval, pksrec.freqInc,
     223          inst->frequencies().getEntry(crpix,crval, pksrec.freqInc,
    209224                                     rec.asuInt("FREQ_ID"));
    210           in->focus().getEntry(pksrec.focusAxi, pksrec.focusTan,
    211                                pksrec.focusRot, tmp0,tmp1,tmp2,tmp3,tmp4,
    212                                rec.asuInt("FOCUS_ID"));
    213           in->molecules().getEntry(pksrec.restFreq,stmp0,stmp1,
     225          inst->focus().getEntry(pksrec.parAngle, pksrec.focusAxi, pksrec.focusTan,
     226                                 pksrec.focusRot, tmp0,tmp1,tmp2,tmp3,tmp4,
     227                                 rec.asuInt("FOCUS_ID"));
     228          inst->molecules().getEntry(pksrec.restFreq,stmp0,stmp1,
    214229                                   rec.asuInt("MOLECULE_ID"));
    215           in->tcal().getEntry(pksrec.tcalTime, pksrec.tcal,
     230          inst->tcal().getEntry(pksrec.tcalTime, pksrec.tcal,
    216231                              rec.asuInt("TCAL_ID"));
    217           in->weather().getEntry(pksrec.temperature, pksrec.pressure,
     232          inst->weather().getEntry(pksrec.temperature, pksrec.pressure,
    218233                                 pksrec.humidity, pksrec.windSpeed,
    219234                                 pksrec.windAz, rec.asuInt("WEATHER_ID"));
     
    230245          pksrec.fieldName = rec.asString("FIELDNAME");
    231246          pksrec.srcName   = rec.asString("SRCNAME");
    232           pksrec.obsType   = hdr.obstype;
     247          //pksrec.obsType   = obstypes[rec.asInt("SRCTYPE")];
     248          pksrec.obsType = getObsTypes( rec.asInt("SRCTYPE") ) ;
    233249          pksrec.bandwidth = nchan * abs(pksrec.freqInc);
    234250          pksrec.azimuth   = rec.asFloat("AZIMUTH");
    235251          pksrec.elevation = rec.asFloat("ELEVATION");
    236           pksrec.parAngle  = rec.asFloat("PARANGLE");
    237252          pksrec.refBeam   = rec.asInt("REFBEAMNO") + 1;
    238253          pksrec.sigma.resize(npol);
     
    293308{
    294309  String poltype = tab.keywordSet().asString("POLTYPE");
    295   if ( poltype == "stokes") {
     310// Full stokes is not supported. Just allow stokes I
     311  if ( poltype == "stokes" && tab.nrow() != 1) {
    296312    String msg = "poltype = " + poltype + " not yet supported in output.";
    297313    throw(AipsError(msg));
     
    340356}
    341357
    342 }
     358// get obsType string from SRCTYPE value
     359String STWriter::getObsTypes( Int srctype )
     360{
     361  String obsType ;
     362  switch( srctype ) {
     363  case Int(SrcType::PSON):
     364    obsType = "PSON" ;
     365    break ;
     366  case Int(SrcType::PSOFF):
     367    obsType = "PSOFF" ;
     368    break ;
     369  case Int(SrcType::NOD):
     370    obsType = "NOD" ;
     371    break ;
     372  case Int(SrcType::FSON):
     373    obsType = "FSON" ;
     374    break ;
     375  case Int(SrcType::FSOFF):
     376    obsType = "FSOFF" ;
     377    break ;
     378  case Int(SrcType::SKY):
     379    obsType = "SKY" ;
     380    break ;
     381  case Int(SrcType::HOT):
     382    obsType = "HOT" ;
     383    break ;
     384  case Int(SrcType::WARM):
     385    obsType = "WARM" ;
     386    break ;
     387  case Int(SrcType::COLD):
     388    obsType = "COLD" ;
     389    break ;
     390  case Int(SrcType::PONCAL):
     391    obsType = "PSON:CALON" ;
     392    break ;
     393  case Int(SrcType::POFFCAL):
     394    obsType = "PSOFF:CALON" ;
     395    break ;
     396  case Int(SrcType::NODCAL):
     397    obsType = "NOD:CALON" ;
     398    break ;
     399  case Int(SrcType::FONCAL):
     400    obsType = "FSON:CALON" ;
     401    break ;
     402  case Int(SrcType::FOFFCAL):
     403    obsType = "FSOFF:CALOFF" ;
     404    break ;
     405  case Int(SrcType::FSLO):
     406    obsType = "FSLO" ;
     407    break ;
     408  case Int(SrcType::FLOOFF):
     409    obsType = "FS:LOWER:OFF" ;
     410    break ;
     411  case Int(SrcType::FLOSKY):
     412    obsType = "FS:LOWER:SKY" ;
     413    break ;
     414  case Int(SrcType::FLOHOT):
     415    obsType = "FS:LOWER:HOT" ;
     416    break ;
     417  case Int(SrcType::FLOWARM):
     418    obsType = "FS:LOWER:WARM" ;
     419    break ;
     420  case Int(SrcType::FLOCOLD):
     421    obsType = "FS:LOWER:COLD" ;
     422    break ;
     423  case Int(SrcType::FSHI):
     424    obsType = "FSHI" ;
     425    break ;
     426  case Int(SrcType::FHIOFF):
     427    obsType = "FS:HIGHER:OFF" ;
     428    break ;
     429  case Int(SrcType::FHISKY):
     430    obsType = "FS:HIGHER:SKY" ;
     431    break ;
     432  case Int(SrcType::FHIHOT):
     433    obsType = "FS:HIGHER:HOT" ;
     434    break ;
     435  case Int(SrcType::FHIWARM):
     436    obsType = "FS:HIGHER:WARM" ;
     437    break ;
     438  case Int(SrcType::FHICOLD):
     439    obsType = "FS:HIGHER:COLD" ;
     440    break ;
     441  default:
     442    obsType = "NOTYPE" ;
     443  }
     444
     445  return obsType ;
     446}
     447
     448}
  • branches/alma/src/STWriter.h

    r1388 r1757  
    3636#include <casa/aips.h>
    3737#include <casa/Utilities/CountedPtr.h>
     38#include <casa/BasicSL/String.h>
    3839
    3940#include "Logger.h"
     
    8485  void replacePtTab(const casa::Table& tab, const std::string& fname);
    8586
     87  casa::String getObsTypes( casa::Int srctype ) ;
     88
    8689  std::string     format_;
    8790  PKSwriter* writer_;
  • branches/alma/src/Scantable.cpp

    r1707 r1757  
    261261  td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
    262262  td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
    263   td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
    264263  td.addColumn(ScalarColumnDesc<Float>("OPACITY"));
    265264
     
    303302  elCol_.attach(table_, "ELEVATION");
    304303  dirCol_.attach(table_, "DIRECTION");
    305   paraCol_.attach(table_, "PARANGLE");
    306304  fldnCol_.attach(table_, "FIELDNAME");
    307305  rbeamCol_.attach(table_, "REFBEAMNO");
    308306
     307  mweatheridCol_.attach(table_,"WEATHER_ID");
    309308  mfitidCol_.attach(table_,"FIT_ID");
    310309  mfreqidCol_.attach(table_, "FREQ_ID");
     
    501500}
    502501
    503 MPosition Scantable::getAntennaPosition () const
     502MPosition Scantable::getAntennaPosition() const
    504503{
    505504  Vector<Double> antpos;
     
    515514  /// @todo reindex SCANNO, recompute nbeam, nif, npol
    516515  inname = path.expandedName();
    517   table_.deepCopy(inname, Table::New);
     516  // WORKAROUND !!! for Table bug
     517  // Remove when fixed in casacore
     518  if ( table_.tableType() == Table::Memory  && !selector_.empty() ) {
     519    Table tab = table_.copyToMemoryTable(generateName());
     520    tab.deepCopy(inname, Table::New);
     521    tab.markForDelete();
     522
     523  } else {
     524    table_.deepCopy(inname, Table::New);
     525  }
    518526}
    519527
     
    646654}
    647655
    648 std::vector<uint> Scantable::getNumbers(ScalarColumn<uInt>& col)
     656std::vector<uint> Scantable::getNumbers(const ScalarColumn<uInt>& col) const
    649657{
    650658  Vector<uInt> nos(col.getColumn());
     
    840848    specCol_.get(whichrow, arr);
    841849  } else {
    842     CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_, basetype));
     850    CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_,
     851                                               basetype));
    843852    uInt row = uInt(whichrow);
    844853    stpol->setSpectra(getPolMatrix(row));
    845854    Float fang,fhand,parang;
    846     fang = focusTable_.getTotalFeedAngle(mfocusidCol_(row));
     855    fang = focusTable_.getTotalAngle(mfocusidCol_(row));
    847856    fhand = focusTable_.getFeedHand(mfocusidCol_(row));
    848     parang = paraCol_(row);
    849     /// @todo re-enable this
    850     // disable total feed angle to support paralactifying Caswell style
    851     stpol->setPhaseCorrections(parang, -parang, fhand);
     857    stpol->setPhaseCorrections(fang, fhand);
    852858    arr = stpol->getSpectrum(requestedpol, ptype);
    853859  }
     
    921927      << setw(15) << "Polarisations:" << setw(4) << npol()
    922928      << "(" << getPolType() << ")" << endl
    923       << setw(15) << "Channels:"  << setw(4) << nchan() << endl;
    924   oss << endl;
     929      << setw(15) << "Channels:" << nchan() << endl;
    925930  String tmp;
    926931  oss << setw(15) << "Observer:"
     
    967972      << setw(10) << "Time" << setw(18) << "Integration" << endl;
    968973  oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
    969   oss << setw(10) << "" << setw(3) << "IF" << setw(6) << ""
     974  oss << setw(10) << "" << setw(3) << "IF" << setw(3) << ""
    970975      << setw(8) << "Frame" << setw(16)
    971       << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment" <<endl;
     976      << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment"
     977      << setw(7) << "Channels"
     978      << endl;
    972979  oss << asap::SEPERATOR << endl;
    973980  TableIterator iter(table_, "SCANNO");
     
    10041011        ROTableRow irow(isubt);
    10051012        const TableRecord& irec = irow.get(0);
    1006         oss << setw(10) << "";
     1013        oss << setw(9) << "";
    10071014        oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
    1008             << setw(2) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
     1015            << setw(1) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
     1016            << setw(3) << "" << nchan(irec.asuInt("IFNO"))
    10091017            << endl;
    10101018
     
    10401048    Double tm;
    10411049    table_.keywordSet().get("UTC",tm);
    1042     return MEpoch(MVEpoch(tm)); 
     1050    return MEpoch(MVEpoch(tm));
    10431051  }
    10441052}
     
    10471055{
    10481056  return formatDirection(getDirection(uInt(whichrow)));
     1057}
     1058
     1059
     1060SpectralCoordinate Scantable::getSpectralCoordinate(int whichrow) const {
     1061  const MPosition& mp = getAntennaPosition();
     1062  const MDirection& md = getDirection(whichrow);
     1063  const MEpoch& me = timeCol_(whichrow);
     1064  //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
     1065  Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
     1066  return freqTable_.getSpectralCoordinate(md, mp, me, rf,
     1067                                          mfreqidCol_(whichrow));
    10491068}
    10501069
     
    10611080    return stlout;
    10621081  }
    1063 
    1064   const MPosition& mp = getAntennaPosition();
    1065   const MDirection& md = getDirection(whichrow);
    1066   const MEpoch& me = timeCol_(whichrow);
    1067   //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
    1068   Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
    1069   SpectralCoordinate spc =
    1070     freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
     1082  SpectralCoordinate spc = getSpectralCoordinate(whichrow);
    10711083  Vector<Double> pixel(nchan);
    10721084  Vector<Double> world;
     
    12251237           Sort::QuickSort|Sort::NoDuplicates );
    12261238  for (uInt i=0; i<fids.nelements(); ++i) {
    1227     frequencies().shiftRefPix(npix, i);
    1228   }
    1229 }
    1230 
    1231 std::string asap::Scantable::getAntennaName() const
     1239    frequencies().shiftRefPix(npix, fids[i]);
     1240  }
     1241}
     1242
     1243std::string Scantable::getAntennaName() const
    12321244{
    12331245  String out;
     
    12361248}
    12371249
    1238 int asap::Scantable::checkScanInfo(const std::vector<int>& scanlist) const
     1250int Scantable::checkScanInfo(const std::vector<int>& scanlist) const
    12391251{
    12401252  String tbpath;
     
    13211333}
    13221334
    1323 std::vector<double>  asap::Scantable::getDirectionVector(int whichrow) const
     1335std::vector<double> Scantable::getDirectionVector(int whichrow) const
    13241336{
    13251337  Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
     
    16811693}
    16821694
     1695std::vector<float> Scantable::getWeather(int whichrow) const
     1696{
     1697  std::vector<float> out(5);
     1698  //Float temperature, pressure, humidity, windspeed, windaz;
     1699  weatherTable_.getEntry(out[0], out[1], out[2], out[3], out[4],
     1700                         mweatheridCol_(uInt(whichrow)));
     1701
     1702
     1703  return out;
     1704}
     1705
    16831706}
    16841707//namespace asap
  • branches/alma/src/Scantable.h

    r1707 r1757  
    2222#include <casa/Utilities/CountedPtr.h>
    2323
    24 #include <casa/Exceptions/Error.h>
    25 
    26 #include <coordinates/Coordinates/SpectralCoordinate.h>
    27 
    2824#include <tables/Tables/Table.h>
    2925#include <tables/Tables/ArrayColumn.h>
     
    3228#include <measures/TableMeasures/ScalarMeasColumn.h>
    3329
     30#include <coordinates/Coordinates/SpectralCoordinate.h>
     31
    3432#include <casa/Arrays/Vector.h>
    3533#include <casa/Quanta/Quantum.h>
     34
     35#include <casa/Exceptions/Error.h>
    3636
    3737#include "Logger.h"
     
    173173   */
    174174  casa::MDirection getDirection( int whichrow ) const;
    175  
     175
    176176  /**
    177177   * get the direction type as a string, e.g. "J2000"
     
    191191   * @return a string describing the direction reference
    192192   */
    193   std::string getDirectionRefString() const;   
     193  std::string getDirectionRefString() const;
    194194
    195195  /**
     
    298298
    299299  int getBeam(int whichrow) const;
    300   std::vector<uint> getBeamNos() { return getNumbers(beamCol_); }
     300  std::vector<uint> getBeamNos() const { return getNumbers(beamCol_); }
    301301
    302302  int getIF(int whichrow) const;
    303   std::vector<uint> getIFNos() { return getNumbers(ifCol_); }
     303  std::vector<uint> getIFNos() const { return getNumbers(ifCol_); }
    304304
    305305  int getPol(int whichrow) const;
    306   std::vector<uint> getPolNos() { return getNumbers(polCol_); }
    307 
    308   std::vector<uint> getScanNos() { return getNumbers(scanCol_); }
     306  std::vector<uint> getPolNos() const { return getNumbers(polCol_); }
     307
     308  std::vector<uint> getScanNos() const { return getNumbers(scanCol_); }
    309309  int getScan(int whichrow) const { return scanCol_(whichrow); }
    310310
     
    335335    { return azCol_(whichrow); }
    336336  float getParAngle(int whichrow) const
    337     { return paraCol_(whichrow); }
     337    { return focus().getParAngle(mfocusidCol_(whichrow)); }
    338338  int getTcalId(int whichrow) const
    339339    { return mtcalidCol_(whichrow); }
     
    384384
    385385  std::vector<double> getAbcissa(int whichrow) const;
     386
     387  std::vector<float> getWeather(int whichrow) const;
    386388
    387389  std::string getAbcissaLabel(int whichrow) const;
     
    408410
    409411  void shift(int npix);
     412
     413  casa::SpectralCoordinate getSpectralCoordinate(int whichrow) const;
    410414
    411415  void convertDirection(const std::string& newframe);
     
    443447   * against the information found in GBT_GO table for
    444448   * scan number orders to get correct pairs.
    445    * @param[in] scan list 
    446    * @return status 
     449   * @param[in] scan list
     450   * @return status
    447451   */
    448452  int checkScanInfo(const std::vector<int>& scanlist) const;
    449453
    450454  /**
    451    * Get the direction as a vector, for a specific row 
     455   * Get the direction as a vector, for a specific row
    452456   * @param[in] whichrow the row numbyyer
    453    * @return the direction in a vector 
     457   * @return the direction in a vector
    454458   */
    455459  std::vector<double> getDirectionVector(int whichrow) const;
     460
     461  /**
     462   * Set a flag indicating whether the data was parallactified
     463   * @param[in] flag true or false
     464   */
     465  void parallactify(bool flag)
     466    { focus().setParallactify(flag); }
    456467
    457468  /**
     
    523534  int rowToScanIndex(int therow);
    524535
    525   std::vector<uint> getNumbers(casa::ScalarColumn<casa::uInt>& col);
    526 
    527   static const casa::uInt version_ = 2;
     536  std::vector<uint> getNumbers(const casa::ScalarColumn<casa::uInt>& col) const;
     537
     538  static const casa::uInt version_ = 3;
    528539
    529540  STSelector selector_;
     
    549560  casa::ScalarColumn<casa::Float> azCol_;
    550561  casa::ScalarColumn<casa::Float> elCol_;
    551   casa::ScalarColumn<casa::Float> paraCol_;
    552562  casa::ScalarColumn<casa::String> srcnCol_, fldnCol_;
    553563  casa::ScalarColumn<casa::uInt> scanCol_, beamCol_, ifCol_, polCol_, cycleCol_, flagrowCol_;
  • branches/alma/src/ScantableWrapper.h

    r1707 r1757  
    2020#include "STFit.h"
    2121#include "Scantable.h"
     22#include "STCoordinate.h"
    2223
    2324namespace asap {
     
    233234    { return table_->checkScanInfo(scanlist); }
    234235 
    235   std::vector<double>  getDirectionVector(int whichrow) const
     236  std::vector<double> getDirectionVector(int whichrow) const
    236237    { return table_->getDirectionVector(whichrow); }
     238
     239  void parallactify(bool flag)
     240    { table_->parallactify(flag); }
     241
     242  STCoordinate getCoordinate(int row) {
     243    return STCoordinate(table_->getSpectralCoordinate(row));
     244  }
     245
     246  std::vector<float> getWeather(int whichrow) const
     247    { return table_->getWeather(whichrow); }
    237248
    238249  void reshapeSpectrum( int nmin, int nmax )
  • branches/alma/src/python_Scantable.cpp

    r1707 r1757  
    133133    .def("_setsourcetype", &ScantableWrapper::setSourceType)
    134134    .def("_getdirectionvec", &ScantableWrapper::getDirectionVector)
     135    .def("_parallactify", &ScantableWrapper::parallactify)
     136    .def("get_coordinate", &ScantableWrapper::getCoordinate)
     137    .def("_get_weather", &ScantableWrapper::getWeather)
    135138    .def("_reshape", &ScantableWrapper::reshapeSpectrum,
    136139         (boost::python::arg("nmin")=-1,
  • branches/alma/src/python_asap.cpp

    r1693 r1757  
    7575  asap::python::python_LineCatalog();
    7676  asap::python::python_Logger();
     77  asap::python::python_STCoordinate();
     78  asap::python::python_STAtmosphere();
    7779  asap::python::python_SrcType();
    7880
  • branches/alma/src/python_asap.h

    r1693 r1757  
    4747    void python_LineCatalog();
    4848    void python_Logger();
     49    void python_STCoordinate();
     50    void python_STAtmosphere();
    4951    void python_SrcType();
    5052
Note: See TracChangeset for help on using the changeset viewer.