Changeset 1066 for trunk/src


Ignore:
Timestamp:
07/04/06 11:22:42 (18 years ago)
Author:
mar637
Message:

Enhancement Ticket #50; arbitrary quotients.

Location:
trunk/src
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/STMath.cpp

    r1033 r1066  
    2121#include <casa/Arrays/ArrayLogical.h>
    2222#include <casa/Arrays/ArrayMath.h>
     23#include <casa/Arrays/Slice.h>
     24#include <casa/Arrays/Slicer.h>
    2325#include <casa/Containers/RecordField.h>
    2426#include <tables/Tables/TableRow.h>
     
    6668{
    6769  if ( avmode == "SCAN" && in.size() != 1 )
    68     throw(AipsError("Can't perform 'SCAN' averaging on multiple tables"));
     70    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
     71                    "Use merge first."));
    6972  WeightType wtype = stringToWeight(weight);
    7073
     
    186189}
    187190
     191
    188192CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
    189193                                             bool droprows)
     
    248252}
    249253
    250 CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable >& in,
    251                                           const std::string & mode,
    252                                           bool preserve )
     254CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
     255                                              const std::string & mode,
     256                                              bool preserve )
    253257{
    254258  /// @todo make other modes available
     
    320324  return out;
    321325}
     326
     327
     328CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
     329                                          const CountedPtr< Scantable > & off,
     330                                          bool preserve )
     331{
     332  bool insitu = insitu_;
     333  setInsitu(false);
     334  CountedPtr< Scantable > out = getScantable(on, false);
     335  setInsitu(insitu);
     336  Table& tout = out->table();
     337  const Table& toff = off->table();
     338  TableIterator sit(tout, "SCANNO");
     339  TableIterator s2it(toff, "SCANNO");
     340  while ( !sit.pastEnd() ) {
     341    Table ton = sit.table();
     342    TableRow row(ton);
     343    Table t = s2it.table();
     344    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
     345    ROArrayColumn<Float> outtsysCol(ton, "TSYS");
     346    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
     347    for (uInt i=0; i < ton.nrow(); ++i) {
     348      const TableRecord& rec = row.get(i);
     349      Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
     350                          && t.col("IFNO") == Int(rec.asuInt("IFNO"))
     351                          && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
     352      TableRow offrow(offsel);
     353      const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
     354      RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
     355      RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
     356      RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
     357      Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
     358      Vector<Float> specon, tsyson;
     359      outtsysCol.get(i, tsyson);
     360      outspecCol.get(i, specon);
     361      Vector<uChar> flagon;
     362      outflagCol.get(i, flagon);
     363      MaskedArray<Float> mon = maskedArray(specon, flagon);
     364      MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
     365      MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
     366      if (preserve) {
     367        quot -= tsysoffscalar;
     368      } else {
     369        quot -= tsyson[0];
     370      }
     371      outspecCol.put(i, quot.getArray());
     372      outflagCol.put(i, flagsFromMA(quot));
     373    }
     374    ++sit;
     375    ++s2it;
     376    // take the first off for each on scan which doesn't have a
     377    // matching off scan
     378    // non <= noff:  matching pairs, non > noff matching pairs then first off
     379    if ( s2it.pastEnd() ) s2it.reset();
     380  }
     381  return out;
     382}
     383
    322384
    323385CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
     
    11961258  return out;
    11971259}
     1260
  • trunk/src/STMath.h

    r995 r1066  
    5858             const std::string& weight = "NONE",
    5959             const std::string& avmode = "SCAN");
     60
    6061  casa::CountedPtr< Scantable >
    6162    averagePolarisations( const casa::CountedPtr< Scantable > & in,
     
    6768                  const std::string& mode, bool tsys=false );
    6869
    69   casa::CountedPtr<Scantable> quotient( const casa::CountedPtr<Scantable>& in,
    70                                         const std::string& mode = "NEAREST",
     70  casa::CountedPtr<Scantable> autoQuotient(const casa::CountedPtr<Scantable>& in,
     71                                           const std::string& mode = "NEAREST",
     72                                           bool preserve = true);
     73
     74  casa::CountedPtr<Scantable> quotient( const casa::CountedPtr<Scantable>& on,
     75                                        const casa::CountedPtr<Scantable>& off,
    7176                                        bool preserve = true );
    7277
  • trunk/src/STMathWrapper.h

    r996 r1066  
    5757  { return ScantableWrapper(STMath::unaryOperate(in.getCP(), val, mode, tsys)); }
    5858
    59   ScantableWrapper quotient( const ScantableWrapper& in,
    60                              const std::string& mode = "NEAREST",
     59  ScantableWrapper autoQuotient( const ScantableWrapper& in,
     60                                 const std::string& mode = "NEAREST",
     61                                 bool preserve = true )
     62  { return ScantableWrapper(STMath::autoQuotient(in.getCP(), mode, preserve)); }
     63
     64  ScantableWrapper quotient( const ScantableWrapper& on,
     65                             const ScantableWrapper& off,
    6166                             bool preserve = true )
    62   { return ScantableWrapper(STMath::quotient(in.getCP(), mode, preserve)); }
     67  { return ScantableWrapper( STMath::quotient( on.getCP(), off.getCP(),
     68                                               preserve ) ); }
    6369
    6470  ScantableWrapper
  • trunk/src/python_STMath.cpp

    r993 r1066  
    4747        .def("_averagepol", &STMathWrapper::averagePolarisations)
    4848        .def("_unaryop", &STMathWrapper::unaryOperate)
     49        .def("_auto_quotient", &STMathWrapper::autoQuotient)
    4950        .def("_quotient", &STMathWrapper::quotient)
    5051        .def("_stats", &STMathWrapper::statistic)
Note: See TracChangeset for help on using the changeset viewer.