Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/hpc33/src/STMath.cpp

    r2345 r2473  
    5353#include "STMath.h"
    5454#include "STSelector.h"
     55#include "Accelerator.h"
    5556
    5657using namespace casa;
    5758using namespace asap;
    5859
     60
     61// 2012/02/17 TN
     62// Since STGrid is implemented, average doesn't consider direction
     63// when accumulating
    5964// tolerance for direction comparison (rad)
    60 #define TOL_OTF    1.0e-15
    61 #define TOL_POINT  2.9088821e-4  // 1 arcmin
     65// #define TOL_OTF    1.0e-15
     66// #define TOL_POINT  2.9088821e-4  // 1 arcmin
    6267
    6368STMath::STMath(bool insitu) :
    … …  
    8388  WeightType wtype = stringToWeight(weight);
    8489
     90  // 2012/02/17 TN
     91  // Since STGrid is implemented, average doesn't consider direction
     92  // when accumulating
    8593  // check if OTF observation
    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   }
     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//   }
    94102
    95103  // output
    … …  
    142150  while (!iter.pastEnd()) {
    143151    Table subt = iter.table();
    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));
     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//       }
    150188//     }
    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 ;
     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 ;
    188199    ++iter;
    189200  }
    … …  
    205216      ROScalarColumn<Double> tmp(tin, "TIME");
    206217      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
    207228      Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
    208229                         && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
    209230                         && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
     231#endif
    210232      Table subt;
    211233      if ( avmode == "SOURCE") {
    … …  
    218240      }
    219241
    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       }
     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//       }
    237262     
    238       if ( nrsubt == removeRows.size() )
    239         throw(AipsError("Averaging data is empty.")) ;
     263//       if ( nrsubt == removeRows.size() )
     264//         throw(AipsError("Averaging data is empty.")) ;
    240265
    241266      specCol.attach(subt,"SPECTRA");
    … …  
    321346{
    322347  (void) mode; // currently unused
     348  // 2012/02/17 TN
     349  // Since STGrid is implemented, average doesn't consider direction
     350  // when accumulating
    323351  // check if OTF observation
    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   }
     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//   }
    332360
    333361  // clone as this is non insitu
    … …  
    362390    flagCol.attach(subt,"FLAGTRA");
    363391    tsysCol.attach(subt,"TSYS");
    364 //     tout.addRow();
    365 //     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
    366 //     if ( avmode != "SCAN") {
    367 //       scanColOut.put(outrowCount, uInt(0));
     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;
     421    ++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;
    368434//     }
    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
     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//       }
    384453//     }
    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;
     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 ;
    392507//     ++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 ;
    474     ++iter;
    475508  }
    476509  return out;
    … …  
    29542987                    "Use merge first."));
    29552988 
     2989  // 2012/02/17 TN
     2990  // Since STGrid is implemented, average doesn't consider direction
     2991  // when accumulating
    29562992  // check if OTF observation
    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   }
     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//   }
    29653001
    29663002  CountedPtr<Scantable> out ;     // processed result
    … …  
    35573593    ScalarColumn<uInt> freqidColOut ;
    35583594    freqidColOut.attach( out->table(), "FREQ_ID" ) ;
    3559     MDirection::ScalarColumn dirColOut ;
    3560     dirColOut.attach( out->table(), "DIRECTION" ) ;
     3595//     MDirection::ScalarColumn dirColOut ;
     3596//     dirColOut.attach( out->table(), "DIRECTION" ) ;
    35613597    Table &tab = tmpout->table() ;
    35623598    Block<String> cols(1);
    … …  
    35753611      ScalarColumn<uInt> polnos ;
    35763612      polnos.attach( iter.table(), "POLNO" ) ;
    3577       MDirection::ScalarColumn dircol ;
    3578       dircol.attach( iter.table(), "DIRECTION" ) ;
     3613//       MDirection::ScalarColumn dircol ;
     3614//       dircol.attach( iter.table(), "DIRECTION" ) ;
    35793615      uInt polno = polnos( 0 ) ;
    35803616      //os << "POLNO iteration: " << polno << LogIO::POST ;
    … …  
    36353671        if ( polout == polno ) {
    36363672          uInt ifout = ifnoColOut( irow ) ;
    3637           Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
     3673//           Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
    36383674          uInt igrp ;
    36393675          for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
    … …  
    36463682            for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
    36473683              uInt ifno = ifnoCol( jrow ) ;
    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 ) ;
     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 ) ;
    36533692              //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
    3654               if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
     3693//               if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
     3694              if ( ifno == freqgrp[igrp][imem] ) {
    36553695                Vector<Float> spec = specCols( jrow ) ;
    36563696                Vector<uChar> flag = flagCols( jrow ) ;
    … …  
    44234463}
    44244464
     4465vector<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
    44254590double STMath::getMJD( string strtime )
    44264591{
    … …  
    44814646  v.push_back( just_before ) ;
    44824647  v.push_back( just_after ) ;
     4648
     4649  return v ;
     4650}
     4651
     4652vector<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 ;
    44834696
    44844697  return v ;
    … …  
    48065019                                            int index )
    48075020{
    4808   string reftime = on->getTime( index ) ;
     5021//   string reftime = on->getTime( index ) ;
     5022  ROTableColumn timeCol( on->table(), "TIME" ) ;
     5023  double reftime = timeCol.asdouble(index) ;
    48095024  vector<int> ii( 1, on->getIF( index ) ) ;
    48105025  vector<int> ib( 1, on->getBeam( index ) ) ;
Note: See TracChangeset for help on using the changeset viewer.