Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STMath.cpp

    r2473 r2345  
    5353#include "STMath.h"
    5454#include "STSelector.h"
    55 #include "Accelerator.h"
    5655
    5756using namespace casa;
    5857using namespace asap;
    5958
    60 
    61 // 2012/02/17 TN
    62 // Since STGrid is implemented, average doesn't consider direction
    63 // when accumulating
    6459// tolerance for direction comparison (rad)
    65 // #define TOL_OTF    1.0e-15
    66 // #define TOL_POINT  2.9088821e-4  // 1 arcmin
     60#define TOL_OTF    1.0e-15
     61#define TOL_POINT  2.9088821e-4  // 1 arcmin
    6762
    6863STMath::STMath(bool insitu) :
    … …  
    8883  WeightType wtype = stringToWeight(weight);
    8984
    90   // 2012/02/17 TN
    91   // Since STGrid is implemented, average doesn't consider direction
    92   // when accumulating
    9385  // check if OTF observation
    94 //   String obstype = in[0]->getHeader().obstype ;
    95 //   Double tol = 0.0 ;
    96 //   if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
    97 //     tol = TOL_OTF ;
    98 //   }
    99 //   else {
    100 //     tol = TOL_POINT ;
    101 //   }
     86  String obstype = in[0]->getHeader().obstype ;
     87  Double tol = 0.0 ;
     88  if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
     89    tol = TOL_OTF ;
     90  }
     91  else {
     92    tol = TOL_POINT ;
     93  }
    10294
    10395  // output
    … …  
    150142  while (!iter.pastEnd()) {
    151143    Table subt = iter.table();
    152     // copy the first row of this selection into the new table
    153     tout.addRow();
    154     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
    155     // re-index to 0
    156     if ( avmode != "SCAN" && avmode != "SOURCE" ) {
    157       scanColOut.put(outrowCount, uInt(0));
    158     }
    159     ++outrowCount;
    160     // 2012/02/17 TN
    161     // Since STGrid is implemented, average doesn't consider direction
    162     // when accumulating
    163 //     MDirection::ScalarColumn dircol ;
    164 //     dircol.attach( subt, "DIRECTION" ) ;
    165 //     Int length = subt.nrow() ;
    166 //     vector< Vector<Double> > dirs ;
    167 //     vector<int> indexes ;
    168 //     for ( Int i = 0 ; i < length ; i++ ) {
    169 //       Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
    170 //       //os << << count++ << ": " ;
    171 //       //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
    172 //       bool adddir = true ;
    173 //       for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
    174 //         //if ( allTrue( t == dirs[j] ) ) {
    175 //         Double dx = t[0] - dirs[j][0] ;
    176 //         Double dy = t[1] - dirs[j][1] ;
    177 //         Double dd = sqrt( dx * dx + dy * dy ) ;
    178 //         //if ( allNearAbs( t, dirs[j], tol ) ) {
    179 //         if ( dd <= tol ) {
    180 //           adddir = false ;
    181 //           break ;
    182 //         }
    183 //       }
    184 //       if ( adddir ) {
    185 //         dirs.push_back( t ) ;
    186 //         indexes.push_back( i ) ;
    187 //       }
     144//     // copy the first row of this selection into the new table
     145//     tout.addRow();
     146//     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     147//     // re-index to 0
     148//     if ( avmode != "SCAN" && avmode != "SOURCE" ) {
     149//       scanColOut.put(outrowCount, uInt(0));
    188150//     }
    189 //     uInt rowNum = dirs.size() ;
    190 //     tout.addRow( rowNum ) ;
    191 //     for ( uInt i = 0 ; i < rowNum ; i++ ) {
    192 //       TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
    193 //       // re-index to 0
    194 //       if ( avmode != "SCAN" && avmode != "SOURCE" ) {
    195 //         scanColOut.put(outrowCount+i, uInt(0));
    196 //       }       
    197 //     }
    198 //     outrowCount += rowNum ;
     151//     ++outrowCount;
     152    MDirection::ScalarColumn dircol ;
     153    dircol.attach( subt, "DIRECTION" ) ;
     154    Int length = subt.nrow() ;
     155    vector< Vector<Double> > dirs ;
     156    vector<int> indexes ;
     157    for ( Int i = 0 ; i < length ; i++ ) {
     158      Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
     159      //os << << count++ << ": " ;
     160      //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
     161      bool adddir = true ;
     162      for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
     163        //if ( allTrue( t == dirs[j] ) ) {
     164        Double dx = t[0] - dirs[j][0] ;
     165        Double dy = t[1] - dirs[j][1] ;
     166        Double dd = sqrt( dx * dx + dy * dy ) ;
     167        //if ( allNearAbs( t, dirs[j], tol ) ) {
     168        if ( dd <= tol ) {
     169          adddir = false ;
     170          break ;
     171        }
     172      }
     173      if ( adddir ) {
     174        dirs.push_back( t ) ;
     175        indexes.push_back( i ) ;
     176      }
     177    }
     178    uInt rowNum = dirs.size() ;
     179    tout.addRow( rowNum ) ;
     180    for ( uInt i = 0 ; i < rowNum ; i++ ) {
     181      TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
     182      // re-index to 0
     183      if ( avmode != "SCAN" && avmode != "SOURCE" ) {
     184        scanColOut.put(outrowCount+i, uInt(0));
     185      }       
     186    }
     187    outrowCount += rowNum ;
    199188    ++iter;
    200189  }
    … …  
    216205      ROScalarColumn<Double> tmp(tin, "TIME");
    217206      Double td;tmp.get(0,td);
    218 
    219 #if 1
    220       static char const*const colNames1[] = { "IFNO", "BEAMNO", "POLNO" };
    221       uInt const values1[] = { rec.asuInt("IFNO"), rec.asuInt("BEAMNO"), rec.asuInt("POLNO") };
    222       SingleTypeEqPredicate<uInt, 3> myPred(tin, colNames1, values1);
    223       CustomTableExprNodeRep myNodeRep(tin, myPred);
    224       myNodeRep.link(); // to avoid automatic delete when myExpr is destructed.
    225       CustomTableExprNode myExpr(myNodeRep);
    226       Table basesubt = tin(myExpr);
    227 #else
    228207      Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
    229208                         && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
    230209                         && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
    231 #endif
    232210      Table subt;
    233211      if ( avmode == "SOURCE") {
    … …  
    240218      }
    241219
    242       // 2012/02/17 TN
    243       // Since STGrid is implemented, average doesn't consider direction
    244       // when accumulating
    245 //       vector<uInt> removeRows ;
    246 //       uInt nrsubt = subt.nrow() ;
    247 //       for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
    248 //         //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
    249 //         Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
    250 //         Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
    251 //         double dx = x0[0] - x1[0];
    252 //         double dy = x0[1] - x1[1];
    253 //         Double dd = sqrt( dx * dx + dy * dy ) ;
    254 //         //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
    255 //         if ( dd > tol ) {
    256 //           removeRows.push_back( irow ) ;
    257 //         }
    258 //       }
    259 //       if ( removeRows.size() != 0 ) {
    260 //         subt.removeRow( removeRows ) ;
    261 //       }
     220      vector<uInt> removeRows ;
     221      uInt nrsubt = subt.nrow() ;
     222      for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
     223        //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
     224        Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
     225        Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
     226        double dx = x0[0] - x1[0];
     227        double dy = x0[1] - x1[1];
     228        Double dd = sqrt( dx * dx + dy * dy ) ;
     229        //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
     230        if ( dd > tol ) {
     231          removeRows.push_back( irow ) ;
     232        }
     233      }
     234      if ( removeRows.size() != 0 ) {
     235        subt.removeRow( removeRows ) ;
     236      }
    262237     
    263 //       if ( nrsubt == removeRows.size() )
    264 //         throw(AipsError("Averaging data is empty.")) ;
     238      if ( nrsubt == removeRows.size() )
     239        throw(AipsError("Averaging data is empty.")) ;
    265240
    266241      specCol.attach(subt,"SPECTRA");
    … …  
    346321{
    347322  (void) mode; // currently unused
    348   // 2012/02/17 TN
    349   // Since STGrid is implemented, average doesn't consider direction
    350   // when accumulating
    351323  // check if OTF observation
    352 //   String obstype = in->getHeader().obstype ;
    353 //   Double tol = 0.0 ;
    354 //   if ( obstype.find( "OTF" ) != String::npos ) {
    355 //     tol = TOL_OTF ;
    356 //   }
    357 //   else {
    358 //     tol = TOL_POINT ;
    359 //   }
     324  String obstype = in->getHeader().obstype ;
     325  Double tol = 0.0 ;
     326  if ( obstype.find( "OTF" ) != String::npos ) {
     327    tol = TOL_OTF ;
     328  }
     329  else {
     330    tol = TOL_POINT ;
     331  }
    360332
    361333  // clone as this is non insitu
    … …  
    390362    flagCol.attach(subt,"FLAGTRA");
    391363    tsysCol.attach(subt,"TSYS");
    392 
    393     tout.addRow();
    394     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
    395     if ( avmode != "SCAN") {
    396       scanColOut.put(outrowCount, uInt(0));
    397     }
    398     Vector<Float> tmp;
    399     specCol.get(0, tmp);
    400     uInt nchan = tmp.nelements();
    401     // have to do channel by channel here as MaskedArrMath
    402     // doesn't have partialMedians
    403     Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
    404     Vector<Float> outspec(nchan);
    405     Vector<uChar> outflag(nchan,0);
    406     Vector<Float> outtsys(1);/// @fixme when tsys is channel based
    407     for (uInt i=0; i<nchan; ++i) {
    408       Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
    409       MaskedArray<Float> ma = maskedArray(specs,flags);
    410       outspec[i] = median(ma);
    411       if ( allEQ(ma.getMask(), False) )
    412         outflag[i] = userflag;// flag data
    413     }
    414     outtsys[0] = median(tsysCol.getColumn());
    415     specColOut.put(outrowCount, outspec);
    416     flagColOut.put(outrowCount, outflag);
    417     tsysColOut.put(outrowCount, outtsys);
    418     Double intsum = sum(intCol.getColumn());
    419     intColOut.put(outrowCount, intsum);
    420     ++outrowCount;
     364//     tout.addRow();
     365//     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     366//     if ( avmode != "SCAN") {
     367//       scanColOut.put(outrowCount, uInt(0));
     368//     }
     369//     Vector<Float> tmp;
     370//     specCol.get(0, tmp);
     371//     uInt nchan = tmp.nelements();
     372//     // have to do channel by channel here as MaskedArrMath
     373//     // doesn't have partialMedians
     374//     Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
     375//     Vector<Float> outspec(nchan);
     376//     Vector<uChar> outflag(nchan,0);
     377//     Vector<Float> outtsys(1);/// @fixme when tsys is channel based
     378//     for (uInt i=0; i<nchan; ++i) {
     379//       Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
     380//       MaskedArray<Float> ma = maskedArray(specs,flags);
     381//       outspec[i] = median(ma);
     382//       if ( allEQ(ma.getMask(), False) )
     383//         outflag[i] = userflag;// flag data
     384//     }
     385//     outtsys[0] = median(tsysCol.getColumn());
     386//     specColOut.put(outrowCount, outspec);
     387//     flagColOut.put(outrowCount, outflag);
     388//     tsysColOut.put(outrowCount, outtsys);
     389//     Double intsum = sum(intCol.getColumn());
     390//     intColOut.put(outrowCount, intsum);
     391//     ++outrowCount;
     392//     ++iter;
     393    MDirection::ScalarColumn dircol ;
     394    dircol.attach( subt, "DIRECTION" ) ;
     395    Int length = subt.nrow() ;
     396    vector< Vector<Double> > dirs ;
     397    vector<int> indexes ;
     398    // Handle MX mode averaging
     399    if (in->nbeam() > 1 ) {     
     400      length = 1;
     401    }
     402    for ( Int i = 0 ; i < length ; i++ ) {
     403      Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
     404      bool adddir = true ;
     405      for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
     406        //if ( allTrue( t == dirs[j] ) ) {
     407        Double dx = t[0] - dirs[j][0] ;
     408        Double dy = t[1] - dirs[j][1] ;
     409        Double dd = sqrt( dx * dx + dy * dy ) ;
     410        //if ( allNearAbs( t, dirs[j], tol ) ) {
     411        if ( dd <= tol ) {
     412          adddir = false ;
     413          break ;
     414        }
     415      }
     416      if ( adddir ) {
     417        dirs.push_back( t ) ;
     418        indexes.push_back( i ) ;
     419      }
     420    }
     421    uInt rowNum = dirs.size() ;
     422    tout.addRow( rowNum );
     423    for ( uInt i = 0 ; i < rowNum ; i++ ) {
     424      TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
     425      // Handle MX mode averaging
     426      if ( avmode != "SCAN") {
     427        scanColOut.put(outrowCount+i, uInt(0));
     428      }
     429    }
     430    MDirection::ScalarColumn dircolOut ;
     431    dircolOut.attach( tout, "DIRECTION" ) ;
     432    for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
     433      Vector<Double> t = \
     434        dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
     435      Vector<Float> tmp;
     436      specCol.get(0, tmp);
     437      uInt nchan = tmp.nelements();
     438      // have to do channel by channel here as MaskedArrMath
     439      // doesn't have partialMedians
     440      Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
     441      // mask spectra for different DIRECTION
     442      for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
     443        Vector<Double> direction = \
     444          dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
     445        //if ( t[0] != direction[0] || t[1] != direction[1] ) {
     446        Double dx = t[0] - direction[0];
     447        Double dy = t[1] - direction[1];
     448        Double dd = sqrt(dx*dx + dy*dy);
     449        //if ( !allNearAbs( t, direction, tol ) ) {
     450        if ( dd > tol &&  in->nbeam() < 2 ) {
     451          flags[jrow] = userflag ;
     452        }
     453      }
     454      Vector<Float> outspec(nchan);
     455      Vector<uChar> outflag(nchan,0);
     456      Vector<Float> outtsys(1);/// @fixme when tsys is channel based
     457      for (uInt i=0; i<nchan; ++i) {
     458        Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
     459        MaskedArray<Float> ma = maskedArray(specs,flags);
     460        outspec[i] = median(ma);
     461        if ( allEQ(ma.getMask(), False) )
     462          outflag[i] = userflag;// flag data
     463      }
     464      outtsys[0] = median(tsysCol.getColumn());
     465      specColOut.put(outrowCount+irow, outspec);
     466      flagColOut.put(outrowCount+irow, outflag);
     467      tsysColOut.put(outrowCount+irow, outtsys);
     468      Vector<Double> integ = intCol.getColumn() ;
     469      MaskedArray<Double> mi = maskedArray( integ, flags ) ;
     470      Double intsum = sum(mi);
     471      intColOut.put(outrowCount+irow, intsum);
     472    }
     473    outrowCount += rowNum ;
    421474    ++iter;
    422 
    423     // 2012/02/17 TN
    424     // Since STGrid is implemented, average doesn't consider direction
    425     // when accumulating
    426 //     MDirection::ScalarColumn dircol ;
    427 //     dircol.attach( subt, "DIRECTION" ) ;
    428 //     Int length = subt.nrow() ;
    429 //     vector< Vector<Double> > dirs ;
    430 //     vector<int> indexes ;
    431 //     // Handle MX mode averaging
    432 //     if (in->nbeam() > 1 ) {     
    433 //       length = 1;
    434 //     }
    435 //     for ( Int i = 0 ; i < length ; i++ ) {
    436 //       Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
    437 //       bool adddir = true ;
    438 //       for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
    439 //         //if ( allTrue( t == dirs[j] ) ) {
    440 //         Double dx = t[0] - dirs[j][0] ;
    441 //         Double dy = t[1] - dirs[j][1] ;
    442 //         Double dd = sqrt( dx * dx + dy * dy ) ;
    443 //         //if ( allNearAbs( t, dirs[j], tol ) ) {
    444 //         if ( dd <= tol ) {
    445 //           adddir = false ;
    446 //           break ;
    447 //         }
    448 //       }
    449 //       if ( adddir ) {
    450 //         dirs.push_back( t ) ;
    451 //         indexes.push_back( i ) ;
    452 //       }
    453 //     }
    454 //     uInt rowNum = dirs.size() ;
    455 //     tout.addRow( rowNum );
    456 //     for ( uInt i = 0 ; i < rowNum ; i++ ) {
    457 //       TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
    458 //       // Handle MX mode averaging
    459 //       if ( avmode != "SCAN") {
    460 //         scanColOut.put(outrowCount+i, uInt(0));
    461 //       }
    462 //     }
    463 //     MDirection::ScalarColumn dircolOut ;
    464 //     dircolOut.attach( tout, "DIRECTION" ) ;
    465 //     for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
    466 //       Vector<Double> t = \
    467 //      dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
    468 //       Vector<Float> tmp;
    469 //       specCol.get(0, tmp);
    470 //       uInt nchan = tmp.nelements();
    471 //       // have to do channel by channel here as MaskedArrMath
    472 //       // doesn't have partialMedians
    473 //       Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
    474 //       // mask spectra for different DIRECTION
    475 //       for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
    476 //         Vector<Double> direction = \
    477 //        dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
    478 //         //if ( t[0] != direction[0] || t[1] != direction[1] ) {
    479 //         Double dx = t[0] - direction[0];
    480 //         Double dy = t[1] - direction[1];
    481 //         Double dd = sqrt(dx*dx + dy*dy);
    482 //         //if ( !allNearAbs( t, direction, tol ) ) {
    483 //         if ( dd > tol &&  in->nbeam() < 2 ) {
    484 //           flags[jrow] = userflag ;
    485 //         }
    486 //       }
    487 //       Vector<Float> outspec(nchan);
    488 //       Vector<uChar> outflag(nchan,0);
    489 //       Vector<Float> outtsys(1);/// @fixme when tsys is channel based
    490 //       for (uInt i=0; i<nchan; ++i) {
    491 //         Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
    492 //         MaskedArray<Float> ma = maskedArray(specs,flags);
    493 //         outspec[i] = median(ma);
    494 //         if ( allEQ(ma.getMask(), False) )
    495 //           outflag[i] = userflag;// flag data
    496 //       }
    497 //       outtsys[0] = median(tsysCol.getColumn());
    498 //       specColOut.put(outrowCount+irow, outspec);
    499 //       flagColOut.put(outrowCount+irow, outflag);
    500 //       tsysColOut.put(outrowCount+irow, outtsys);
    501 //       Vector<Double> integ = intCol.getColumn() ;
    502 //       MaskedArray<Double> mi = maskedArray( integ, flags ) ;
    503 //       Double intsum = sum(mi);
    504 //       intColOut.put(outrowCount+irow, intsum);
    505 //     }
    506 //     outrowCount += rowNum ;
    507 //     ++iter;
    508475  }
    509476  return out;
    … …  
    29872954                    "Use merge first."));
    29882955 
    2989   // 2012/02/17 TN
    2990   // Since STGrid is implemented, average doesn't consider direction
    2991   // when accumulating
    29922956  // check if OTF observation
    2993 //   String obstype = in[0]->getHeader().obstype ;
    2994 //   Double tol = 0.0 ;
    2995 //   if ( obstype.find( "OTF" ) != String::npos ) {
    2996 //     tol = TOL_OTF ;
    2997 //   }
    2998 //   else {
    2999 //     tol = TOL_POINT ;
    3000 //   }
     2957  String obstype = in[0]->getHeader().obstype ;
     2958  Double tol = 0.0 ;
     2959  if ( obstype.find( "OTF" ) != String::npos ) {
     2960    tol = TOL_OTF ;
     2961  }
     2962  else {
     2963    tol = TOL_POINT ;
     2964  }
    30012965
    30022966  CountedPtr<Scantable> out ;     // processed result
    … …  
    35933557    ScalarColumn<uInt> freqidColOut ;
    35943558    freqidColOut.attach( out->table(), "FREQ_ID" ) ;
    3595 //     MDirection::ScalarColumn dirColOut ;
    3596 //     dirColOut.attach( out->table(), "DIRECTION" ) ;
     3559    MDirection::ScalarColumn dirColOut ;
     3560    dirColOut.attach( out->table(), "DIRECTION" ) ;
    35973561    Table &tab = tmpout->table() ;
    35983562    Block<String> cols(1);
    … …  
    36113575      ScalarColumn<uInt> polnos ;
    36123576      polnos.attach( iter.table(), "POLNO" ) ;
    3613 //       MDirection::ScalarColumn dircol ;
    3614 //       dircol.attach( iter.table(), "DIRECTION" ) ;
     3577      MDirection::ScalarColumn dircol ;
     3578      dircol.attach( iter.table(), "DIRECTION" ) ;
    36153579      uInt polno = polnos( 0 ) ;
    36163580      //os << "POLNO iteration: " << polno << LogIO::POST ;
    … …  
    36713635        if ( polout == polno ) {
    36723636          uInt ifout = ifnoColOut( irow ) ;
    3673 //           Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
     3637          Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
    36743638          uInt igrp ;
    36753639          for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
    … …  
    36823646            for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
    36833647              uInt ifno = ifnoCol( jrow ) ;
    3684               // 2012/02/17 TN
    3685               // Since STGrid is implemented, average doesn't consider direction
    3686               // when accumulating
    3687 //               Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
    3688 //               //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction  ) ) {
    3689 //               Double dx = tdir[0] - direction[0] ;
    3690 //               Double dy = tdir[1] - direction[1] ;
    3691 //               Double dd = sqrt( dx * dx + dy * dy ) ;
     3648              Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
     3649              //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction  ) ) {
     3650              Double dx = tdir[0] - direction[0] ;
     3651              Double dy = tdir[1] - direction[1] ;
     3652              Double dd = sqrt( dx * dx + dy * dy ) ;
    36923653              //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
    3693 //               if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
    3694               if ( ifno == freqgrp[igrp][imem] ) {
     3654              if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
    36953655                Vector<Float> spec = specCols( jrow ) ;
    36963656                Vector<uChar> flag = flagCols( jrow ) ;
    … …  
    44634423}
    44644424
    4465 vector<float> STMath::getSpectrumFromTime( double reftime,
    4466                                            CountedPtr<Scantable>& s,
    4467                                            string mode )
    4468 {
    4469   LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
    4470   vector<float> sp ;
    4471 
    4472   if ( s->nrow() == 0 ) {
    4473     os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
    4474     return sp ;
    4475   }
    4476   else if ( s->nrow() == 1 ) {
    4477     //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
    4478     return s->getSpectrum( 0 ) ;
    4479   }
    4480   else {
    4481     ROScalarColumn<Double> timeCol( s->table(), "TIME" ) ;
    4482     Vector<Double> timeVec = timeCol.getColumn() ;
    4483     vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;
    4484     if ( mode == "before" ) {
    4485       int id = -1 ;
    4486       if ( idx[0] != -1 ) {
    4487         id = idx[0] ;
    4488       }
    4489       else if ( idx[1] != -1 ) {
    4490         os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
    4491         id = idx[1] ;
    4492       }
    4493       //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4494       sp = s->getSpectrum( id ) ;
    4495     }
    4496     else if ( mode == "after" ) {
    4497       int id = -1 ;
    4498       if ( idx[1] != -1 ) {
    4499         id = idx[1] ;
    4500       }
    4501       else if ( idx[0] != -1 ) {
    4502         os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
    4503         id = idx[1] ;
    4504       }
    4505       //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4506       sp = s->getSpectrum( id ) ;
    4507     }
    4508     else if ( mode == "nearest" ) {
    4509       int id = -1 ;
    4510       if ( idx[0] == -1 ) {
    4511         id = idx[1] ;
    4512       }
    4513       else if ( idx[1] == -1 ) {
    4514         id = idx[0] ;
    4515       }
    4516       else if ( idx[0] == idx[1] ) {
    4517         id = idx[0] ;
    4518       }
    4519       else {
    4520         //double t0 = getMJD( s->getTime( idx[0] ) ) ;
    4521         //double t1 = getMJD( s->getTime( idx[1] ) ) ;
    4522 //         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
    4523 //         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
    4524         double t0 = timeVec[idx[0]] ;
    4525         double t1 = timeVec[idx[1]] ;
    4526 //         cout << "t0-t0c=" << t0-t0c << endl ;
    4527 //         cout << "t1-t1c=" << t1-t1c << endl ;
    4528 //         double tref = getMJD( reftime ) ;
    4529         double tref = reftime ;
    4530         if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
    4531           id = idx[1] ;
    4532         }
    4533         else {
    4534           id = idx[0] ;
    4535         }
    4536       }
    4537       //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4538       sp = s->getSpectrum( id ) ;     
    4539     }
    4540     else if ( mode == "linear" ) {
    4541       if ( idx[0] == -1 ) {
    4542         // use after
    4543         os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
    4544         int id = idx[1] ;
    4545         //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4546         sp = s->getSpectrum( id ) ;
    4547       }
    4548       else if ( idx[1] == -1 ) {
    4549         // use before
    4550         os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
    4551         int id = idx[0] ;
    4552         //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4553         sp = s->getSpectrum( id ) ;
    4554       }
    4555       else if ( idx[0] == idx[1] ) {
    4556         // use before
    4557         //os << "No need to interporate." << LogIO::POST ;
    4558         int id = idx[0] ;
    4559         //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    4560         sp = s->getSpectrum( id ) ;
    4561       }
    4562       else {
    4563         // do interpolation
    4564         //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
    4565         //double t0 = getMJD( s->getTime( idx[0] ) ) ;
    4566         //double t1 = getMJD( s->getTime( idx[1] ) ) ;
    4567 //         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
    4568 //         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
    4569         double t0 = timeVec[idx[0]] ;
    4570         double t1 = timeVec[idx[1]] ;
    4571 //         cout << "t0-t0c=" << t0-t0c << endl ;
    4572 //         cout << "t1-t1c=" << t1-t1c << endl ;
    4573 //         double tref = getMJD( reftime ) ;
    4574         double tref = reftime ;
    4575         vector<float> sp0 = s->getSpectrum( idx[0] ) ;
    4576         vector<float> sp1 = s->getSpectrum( idx[1] ) ;
    4577         for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {
    4578           float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
    4579           sp.push_back( v ) ;
    4580         }
    4581       }
    4582     }
    4583     else {
    4584       os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
    4585     }
    4586     return sp ;
    4587   }
    4588 }
    4589 
    45904425double STMath::getMJD( string strtime )
    45914426{
    … …  
    46464481  v.push_back( just_before ) ;
    46474482  v.push_back( just_after ) ;
    4648 
    4649   return v ;
    4650 }
    4651 
    4652 vector<int> STMath::getRowIdFromTime( double reftime, Vector<Double> &t )
    4653 {
    4654 //   double reft = reftime ;
    4655   double dtmin = 1.0e100 ;
    4656   double dtmax = -1.0e100 ;
    4657 //   vector<double> dt ;
    4658   int just_before = -1 ;
    4659   int just_after = -1 ;
    4660   //cout << setprecision(24) << reft << endl ;
    4661 //   ROScalarColumn<Double> timeCol( s->table(), "TIME" ) ;
    4662 //   for ( int i = 0 ; i < s->nrow() ; i++ ) {
    4663 //     cout << setprecision(24) << timeCol(i) << endl ;
    4664 //     //dt.push_back( getMJD( s->getTime( i ) ) - reft ) ;
    4665 //     dt.push_back( timeCol(i) - reft ) ;
    4666 //   }
    4667   Vector<Double> dt = t - reftime ;
    4668   for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
    4669     if ( dt[i] > 0.0 ) {
    4670       // after reftime
    4671       if ( dt[i] < dtmin ) {
    4672         just_after = i ;
    4673         dtmin = dt[i] ;
    4674       }
    4675     }
    4676     else if ( dt[i] < 0.0 ) {
    4677       // before reftime
    4678       if ( dt[i] > dtmax ) {
    4679         just_before = i ;
    4680         dtmax = dt[i] ;
    4681       }
    4682     }
    4683     else {
    4684       // just a reftime
    4685       just_before = i ;
    4686       just_after = i ;
    4687       dtmax = 0 ;
    4688       dtmin = 0 ;
    4689       break ;
    4690     }
    4691   }
    4692 
    4693   vector<int> v(2) ;
    4694   v[0] = just_before ;
    4695   v[1] = just_after ;
    46964483
    46974484  return v ;
    … …  
    50194806                                            int index )
    50204807{
    5021 //   string reftime = on->getTime( index ) ;
    5022   ROTableColumn timeCol( on->table(), "TIME" ) ;
    5023   double reftime = timeCol.asdouble(index) ;
     4808  string reftime = on->getTime( index ) ;
    50244809  vector<int> ii( 1, on->getIF( index ) ) ;
    50254810  vector<int> ib( 1, on->getBeam( index ) ) ;
Note: See TracChangeset for help on using the changeset viewer.