Changeset 1611 for branches/alma/src


Ignore:
Timestamp:
08/03/09 11:41:38 (15 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-1423

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Averaging methods are changed to refer DIRECTION all the time.
This change affects the behavior of sdaverage when timeaverage is
set to True.
Tolerance for direction is set to 1.0e-15 rad for OTF data and
is set to 1.0e-5 rad (around 2 arcsec) for non-OTF data.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/alma/src/STMath.cpp

    r1609 r1611  
    5555
    5656// tolerance for direction comparison (rad)
    57 Double tol = 1.0e-15 ;
    58 //Double tol = 1.0 ;
     57#define TOL_OTF    1.0e-15
     58#define TOL_POINT  1.0e-5  // 2 arcsec
    5959
    6060STMath::STMath(bool insitu) :
    … …  
    8181  // check if OTF observation
    8282  String obstype = in[0]->getHeader().obstype ;
    83   bool otfscan = false ;
     83  //bool otfscan = false ;
     84  Double tol = 0.0 ;
    8485  if ( obstype.find( "OTF" ) != String::npos ) {
    8586    //cout << "OTF scan" << endl ;
    86     otfscan = true ;
     87    //otfscan = true ;
     88    tol = TOL_OTF ;
     89  }
     90  else {
     91    tol = TOL_POINT ;
    8792  }
    8893
    … …  
    144149//     }
    145150//     ++outrowCount;
    146     if ( otfscan ) {
    147       MDirection::ScalarColumn dircol ;
    148       dircol.attach( subt, "DIRECTION" ) ;
    149       Int length = subt.nrow() ;
    150       vector< Vector<Double> > dirs ;
    151       vector<int> indexes ;
    152       for ( Int i = 0 ; i < length ; i++ ) {
    153         Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
    154         //cout << setw( 3 ) << count++ << ": " ;
    155         //cout << "[" << t[0] << "," << t[1] << "]" <<  endl ;
    156         bool adddir = true ;
    157         for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
    158           //if ( allTrue( t == dirs[j] ) ) {
    159           if ( allNearAbs( t, dirs[j], tol ) ) {
    160             adddir = false ;
    161             break ;
    162           }
     151    MDirection::ScalarColumn dircol ;
     152    dircol.attach( subt, "DIRECTION" ) ;
     153    Int length = subt.nrow() ;
     154    vector< Vector<Double> > dirs ;
     155    vector<int> indexes ;
     156    for ( Int i = 0 ; i < length ; i++ ) {
     157      Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
     158      //cout << setw( 3 ) << count++ << ": " ;
     159      //cout << "[" << t[0] << "," << t[1] << "]" <<  endl ;
     160      bool adddir = true ;
     161      for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
     162        //if ( allTrue( t == dirs[j] ) ) {
     163        Double dx = t[0] - dirs[j][0] ;
     164        Double dy = t[1] - dirs[j][1] ;
     165        Double dd = sqrt( dx * dx + dy * dy ) ;
     166        //if ( allNearAbs( t, dirs[j], tol ) ) {
     167        if ( dd <= tol ) {
     168          adddir = false ;
     169          break ;
    163170        }
    164         if ( adddir ) {
    165           dirs.push_back( t ) ;
    166           indexes.push_back( i ) ;
    167         }
    168       }
    169       uInt rowNum = dirs.size() ;
    170       //cout << "dirs.size()=" << dirs.size() << endl ;
    171       tout.addRow( rowNum ) ;
    172       for ( uInt i = 0 ; i < rowNum ; i++ ) {
    173         TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
    174         // re-index to 0
    175         if ( avmode != "SCAN" && avmode != "SOURCE" ) {
    176           scanColOut.put(outrowCount+i, uInt(0));
    177         }       
    178       }
    179       outrowCount += rowNum ;
    180     }
    181     else {
    182       // copy the first row of this selection into the new table
    183       tout.addRow();
    184       TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     171      }
     172      if ( adddir ) {
     173        dirs.push_back( t ) ;
     174        indexes.push_back( i ) ;
     175      }
     176    }
     177    uInt rowNum = dirs.size() ;
     178    //cout << "dirs.size()=" << dirs.size() << endl ;
     179    tout.addRow( rowNum ) ;
     180    for ( uInt i = 0 ; i < rowNum ; i++ ) {
     181      TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
    185182      // re-index to 0
    186183      if ( avmode != "SCAN" && avmode != "SOURCE" ) {
    187         scanColOut.put(outrowCount, uInt(0));
    188       }
    189       ++outrowCount;
    190     }
     184        scanColOut.put(outrowCount+i, uInt(0));
     185      }       
     186    }
     187    outrowCount += rowNum ;
    191188    ++iter;
    192189  }
    … …  
    221218        subt = basesubt;
    222219      }
    223       // for OTF
    224       if ( otfscan ) {
    225         vector<uInt> removeRows ;
    226         uInt nrsubt = subt.nrow() ;
    227         for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
    228           //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
    229           if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
    230             removeRows.push_back( irow ) ;
    231           }
     220
     221      vector<uInt> removeRows ;
     222      uInt nrsubt = subt.nrow() ;
     223      for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
     224        //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
     225        Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
     226        Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
     227        double dx = x0[0] - x1[0] ;
     228        double dy = x0[0] - x1[0] ;
     229        Double dd = sqrt( dx * dx + dy * dy ) ;
     230        //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
     231        if ( dd > tol ) {
     232          removeRows.push_back( irow ) ;
    232233        }
    233         //cout << "removeRows.size()=" << removeRows.size() << endl ;
    234         if ( removeRows.size() != 0 ) {
    235           //cout << "[" ;
    236           //for ( uInt irow=0 ; irow<removeRows.size()-1 ; irow++ )
    237           //cout << removeRows[irow] << "," ;
    238           //cout << removeRows[removeRows.size()-1] << "]" << endl ;
    239           subt.removeRow( removeRows ) ;
    240         }
    241        
    242         if ( nrsubt == removeRows.size() )
    243           throw(AipsError("Averaging data is empty.")) ;
    244       }
     234      }
     235      //cout << "removeRows.size()=" << removeRows.size() << endl ;
     236      if ( removeRows.size() != 0 ) {
     237        //cout << "[" ;
     238        //for ( uInt irow=0 ; irow<removeRows.size()-1 ; irow++ )
     239        //cout << removeRows[irow] << "," ;
     240        //cout << removeRows[removeRows.size()-1] << "]" << endl ;
     241        subt.removeRow( removeRows ) ;
     242      }
     243     
     244      if ( nrsubt == removeRows.size() )
     245        throw(AipsError("Averaging data is empty.")) ;
     246
    245247      specCol.attach(subt,"SPECTRA");
    246248      flagCol.attach(subt,"FLAGTRA");
    … …  
    312314  // check if OTF observation
    313315  String obstype = in->getHeader().obstype ;
    314   bool otfscan = false ;
     316  //bool otfscan = false ;
     317  Double tol = 0.0 ;
    315318  if ( obstype.find( "OTF" ) != String::npos ) {
    316319    //cout << "OTF scan" << endl ;
    317     otfscan = true ;
     320    //otfscan = true ;
     321    tol = TOL_OTF ;
     322  }
     323  else {
     324    tol = TOL_POINT ;
    318325  }
    319326
    … …  
    378385//     ++outrowCount;
    379386//     ++iter;
    380     if ( otfscan ) {
    381       MDirection::ScalarColumn dircol ;
    382       dircol.attach( subt, "DIRECTION" ) ;
    383       Int length = subt.nrow() ;
    384       vector< Vector<Double> > dirs ;
    385       vector<int> indexes ;
    386       for ( Int i = 0 ; i < length ; i++ ) {
    387         Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
    388         bool adddir = true ;
    389         for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
    390           //if ( allTrue( t == dirs[j] ) ) {
    391           if ( allNearAbs( t, dirs[j], tol ) ) {
    392             adddir = false ;
    393             break ;
    394           }
     387    MDirection::ScalarColumn dircol ;
     388    dircol.attach( subt, "DIRECTION" ) ;
     389    Int length = subt.nrow() ;
     390    vector< Vector<Double> > dirs ;
     391    vector<int> indexes ;
     392    for ( Int i = 0 ; i < length ; i++ ) {
     393      Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
     394      bool adddir = true ;
     395      for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
     396        //if ( allTrue( t == dirs[j] ) ) {
     397        Double dx = t[0] - dirs[j][0] ;
     398        Double dy = t[1] - dirs[j][1] ;
     399        Double dd = sqrt( dx * dx + dy * dy ) ;
     400        //if ( allNearAbs( t, dirs[j], tol ) ) {
     401        if ( dd <= tol ) {
     402          adddir = false ;
     403          break ;
    395404        }
    396         if ( adddir ) {
    397           dirs.push_back( t ) ;
    398           indexes.push_back( i ) ;
    399         }
    400       }
    401       uInt rowNum = dirs.size() ;
    402       tout.addRow( rowNum );
    403       for ( uInt i = 0 ; i < rowNum ; i++ ) {
    404         TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
    405         if ( avmode != "SCAN") {
    406           //scanColOut.put(outrowCount+i, uInt(0));
    407         }
    408       }
    409       MDirection::ScalarColumn dircolOut ;
    410       dircolOut.attach( tout, "DIRECTION" ) ;
    411       for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
    412         Vector<Double> t = dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
    413         Vector<Float> tmp;
    414         specCol.get(0, tmp);
    415         uInt nchan = tmp.nelements();
    416         // have to do channel by channel here as MaskedArrMath
    417         // doesn't have partialMedians
    418         Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
    419         // mask spectra for different DIRECTION
    420         //cout << "irow=" << outrowCount+irow << ": flagged [" ;
    421         for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
    422           Vector<Double> direction = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
    423           //if ( t[0] != direction[0] || t[1] != direction[1] ) {
    424           if ( !allNearAbs( t, direction, tol ) ) {
    425             //cout << jrow << " " ;
    426             flags[jrow] = userflag ;
    427           }
    428         }
    429         //cout << "]" << endl ;
    430         //cout << "flags=" << flags << endl ;
    431         Vector<Float> outspec(nchan);
    432         Vector<uChar> outflag(nchan,0);
    433         Vector<Float> outtsys(1);/// @fixme when tsys is channel based
    434         for (uInt i=0; i<nchan; ++i) {
    435           Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
    436           MaskedArray<Float> ma = maskedArray(specs,flags);
    437           outspec[i] = median(ma);
    438           if ( allEQ(ma.getMask(), False) )
    439             outflag[i] = userflag;// flag data
    440         }
    441         outtsys[0] = median(tsysCol.getColumn());
    442         specColOut.put(outrowCount+irow, outspec);
    443         flagColOut.put(outrowCount+irow, outflag);
    444         tsysColOut.put(outrowCount+irow, outtsys);
    445         Vector<Double> integ = intCol.getColumn() ;
    446         MaskedArray<Double> mi = maskedArray( integ, flags ) ;
    447         Double intsum = sum(mi);
    448         intColOut.put(outrowCount+irow, intsum);
    449       }
    450       outrowCount += rowNum ;
    451     }
    452     else {
    453       tout.addRow();
    454       TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     405      }
     406      if ( adddir ) {
     407        dirs.push_back( t ) ;
     408        indexes.push_back( i ) ;
     409      }
     410    }
     411    uInt rowNum = dirs.size() ;
     412    tout.addRow( rowNum );
     413    for ( uInt i = 0 ; i < rowNum ; i++ ) {
     414      TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
    455415      if ( avmode != "SCAN") {
    456         scanColOut.put(outrowCount, uInt(0));
    457       }
     416        //scanColOut.put(outrowCount+i, uInt(0));
     417      }
     418    }
     419    MDirection::ScalarColumn dircolOut ;
     420    dircolOut.attach( tout, "DIRECTION" ) ;
     421    for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
     422      Vector<Double> t = dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
    458423      Vector<Float> tmp;
    459424      specCol.get(0, tmp);
    … …  
    462427      // doesn't have partialMedians
    463428      Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
     429      // mask spectra for different DIRECTION
     430      //cout << "irow=" << outrowCount+irow << ": flagged [" ;
     431      for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
     432        Vector<Double> direction = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
     433        //if ( t[0] != direction[0] || t[1] != direction[1] ) {
     434        Double dx = t[0] - direction[0] ;
     435        Double dy = t[1] - direction[1] ;
     436        Double dd = sqrt( dx * dx + dy * dy ) ;
     437        //if ( !allNearAbs( t, direction, tol ) ) {
     438        if ( dd > tol ) {
     439          //cout << jrow << " " ;
     440          flags[jrow] = userflag ;
     441        }
     442      }
     443      //cout << "]" << endl ;
     444      //cout << "flags=" << flags << endl ;
    464445      Vector<Float> outspec(nchan);
    465446      Vector<uChar> outflag(nchan,0);
    … …  
    473454      }
    474455      outtsys[0] = median(tsysCol.getColumn());
    475       specColOut.put(outrowCount, outspec);
    476       flagColOut.put(outrowCount, outflag);
    477       tsysColOut.put(outrowCount, outtsys);
    478       Double intsum = sum(intCol.getColumn());
    479       intColOut.put(outrowCount, intsum);
    480       ++outrowCount;
    481     }
     456      specColOut.put(outrowCount+irow, outspec);
     457      flagColOut.put(outrowCount+irow, outflag);
     458      tsysColOut.put(outrowCount+irow, outtsys);
     459      Vector<Double> integ = intCol.getColumn() ;
     460      MaskedArray<Double> mi = maskedArray( integ, flags ) ;
     461      Double intsum = sum(mi);
     462      intColOut.put(outrowCount+irow, intsum);
     463    }
     464    outrowCount += rowNum ;
    482465    ++iter;
    483466  }
    … …  
    24072390  // check if OTF observation
    24082391  String obstype = in[0]->getHeader().obstype ;
    2409   bool otfscan = false ;
     2392  //bool otfscan = false ;
     2393  Double tol = 0.0 ;
    24102394  if ( obstype.find( "OTF" ) != String::npos ) {
    24112395    //cout << "OTF scan" << endl ;
    2412     otfscan = true ;
     2396    //otfscan = true ;
     2397    tol = TOL_OTF ;
     2398  }
     2399  else {
     2400    tol = TOL_POINT ;
    24132401  }
    24142402
    … …  
    30913079              Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
    30923080              //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction  ) ) {
    3093               if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
     3081              Double dx = tdir[0] - direction[0] ;
     3082              Double dy = tdir[1] - direction[1] ;
     3083              Double dd = sqrt( dx * dx + dy * dy ) ;
     3084              //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
     3085              if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
    30943086                Vector<Float> spec = specCols( jrow ) ;
    30953087                Vector<uChar> flag = flagCols( jrow ) ;
Note: See TracChangeset for help on using the changeset viewer.