Changeset 699


Ignore:
Timestamp:
10/20/05 16:35:24 (19 years ago)
Author:
mar637
Message:

whitespace tidy

Location:
branches/Release-2-fixes/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/Release-2-fixes/src/SDMath.cc

    r653 r699  
    110110
    111111
    112 
    113112SDMemTable* SDMath::frequencyAlignment(const SDMemTable& in,
    114113                                       const String& refTime,
     
    116115                                       Bool perFreqID) const
    117116{
    118 // Get frame info from Table
    119 
     117  // Get frame info from Table
    120118   std::vector<std::string> info = in.getCoordInfo();
    121119
    122 // Parse frequency system
    123 
     120   // Parse frequency system
    124121   String systemStr(info[1]);
    125122   String baseSystemStr(info[3]);
    126123   if (baseSystemStr==systemStr) {
    127       throw(AipsError("You have not set a frequency frame different from the initial - use function set_freqframe"));
     124     throw(AipsError("You have not set a frequency frame different from the initial - use function set_freqframe"));
    128125   }
    129 //
     126
    130127   MFrequency::Types freqSystem;
    131128   MFrequency::getType(freqSystem, systemStr);
    132129
    133 // Do it
    134 
    135130   return frequencyAlign(in, freqSystem, refTime, method, perFreqID);
    136131}
     
    138133
    139134
    140 CountedPtr<SDMemTable> SDMath::average(const std::vector<CountedPtr<SDMemTable> >& in,
    141                                        const Vector<Bool>& mask, Bool scanAv,
    142                                        const String& weightStr, Bool alignFreq) const
    143 //
     135CountedPtr<SDMemTable>
     136SDMath::average(const std::vector<CountedPtr<SDMemTable> >& in,
     137                const Vector<Bool>& mask, Bool scanAv,
     138                const String& weightStr, Bool alignFreq) const
    144139// Weighted averaging of spectra from one or more Tables.
    145 //
    146 {
    147 
    148 // Convert weight type
    149  
     140{
     141  // Convert weight type 
    150142  WeightType wtType = NONE;
    151143  convertWeightString(wtType, weightStr, True);
    152144
    153 // Create output Table by cloning from the first table
    154 
     145  // Create output Table by cloning from the first table
    155146  SDMemTable* pTabOut = new SDMemTable(*in[0],True);
    156147  if (in.size() > 1) {
     
    159150    }
    160151  }
    161 // Setup
    162 
     152  // Setup
    163153  IPosition shp = in[0]->rowAsMaskedArray(0).shape();      // Must not change
    164154  Array<Float> arr(shp);
     
    166156  const Bool useMask = (mask.nelements() == shp(asap::ChanAxis));
    167157
    168 // Columns from Tables
    169 
     158  // Columns from Tables
    170159  ROArrayColumn<Float> tSysCol;
    171160  ROScalarColumn<Double> mjdCol;
     
    175164  ROScalarColumn<Int> scanIDCol;
    176165
    177 // Create accumulation MaskedArray. We accumulate for each channel,if,pol,beam
    178 // Note that the mask of the accumulation array will ALWAYS remain ALL True.
    179 // The MA is only used so that when data which is masked Bad is added to it,
    180 // that data does not contribute.
     166  // Create accumulation MaskedArray. We accumulate for each
     167  // channel,if,pol,beam Note that the mask of the accumulation array
     168  // will ALWAYS remain ALL True.  The MA is only used so that when
     169  // data which is masked Bad is added to it, that data does not
     170  // contribute.
    181171
    182172  Array<Float> zero(shp);
     
    186176  MaskedArray<Float> sum(zero,good);
    187177
    188 // Counter arrays
    189 
     178  // Counter arrays
    190179  Array<Float> nPts(shp);             // Number of points
    191180  nPts = 0.0;
     
    193182  nInc = 1.0;
    194183
    195 // Create accumulation Array for variance. We accumulate for
    196 // each if,pol,beam, but average over channel.  So we need
    197 // a shape with one less axis dropping channels.
    198 
     184  // Create accumulation Array for variance. We accumulate for each
     185  // if,pol,beam, but average over channel.  So we need a shape with
     186  // one less axis dropping channels.
    199187  const uInt nAxesSub = shp.nelements() - 1;
    200188  IPosition shp2(nAxesSub);
     
    209197  IPosition pos2(nAxesSub,0);                        // For indexing
    210198
    211 // Time-related accumulators
    212 
     199  // Time-related accumulators
    213200  Double time;
    214201  Double timeSum = 0.0;
     
    216203  Double interval = 0.0;
    217204
    218 // To get the right shape for the Tsys accumulator we need to
    219 // access a column from the first table.  The shape of this
    220 // array must not change.  Note however that since the TSysSqSum
    221 // array is used in a normalization process, and that I ignore the
    222 // channel axis replication of values for now, it loses a dimension
     205  // To get the right shape for the Tsys accumulator we need to access
     206  // a column from the first table.  The shape of this array must not
     207  // change.  Note however that since the TSysSqSum array is used in a
     208  // normalization process, and that I ignore the channel axis
     209  // replication of values for now, it loses a dimension
    223210
    224211  Array<Float> tSysSum, tSysSqSum;
     
    227214    tSysCol.attach(tabIn,"TSYS");
    228215    tSysSum.resize(tSysCol.shape(0));
    229 //
    230216    tSysSqSum.resize(shp2);
    231217  }
    232   tSysSum =0.0;
     218  tSysSum = 0.0;
    233219  tSysSqSum = 0.0;
    234220  Array<Float> tSys;
    235221
    236 // Scan and row tracking
    237 
     222  // Scan and row tracking
    238223  Int oldScanID = 0;
    239224  Int outScanID = 0;
     
    243228  Int tableStart = 0;
    244229
    245 // Source and FreqID
    246 
     230  // Source and FreqID
    247231  String sourceName, oldSourceName, sourceNameStart;
    248232  Vector<uInt> freqID, freqIDStart, oldFreqID;
    249233
    250 // Loop over tables
    251 
     234  // Loop over tables
    252235  Float fac = 1.0;
    253236  const uInt nTables = in.size();
    254237  for (uInt iTab=0; iTab<nTables; iTab++) {
    255238
    256 // Should check that the frequency tables don't change if doing FreqAlignment
    257 
    258 // Attach columns to Table
    259 
     239    // Should check that the frequency tables don't change if doing
     240    // FreqAlignment
     241   
     242    // Attach columns to Table
    260243     const Table& tabIn = in[iTab]->table();
    261244     tSysCol.attach(tabIn, "TSYS");
     
    266249     scanIDCol.attach(tabIn, "SCANID");
    267250
    268 // Loop over rows in Table
    269 
     251     // Loop over rows in Table
    270252     const uInt nRows = in[iTab]->nRow();
    271253     for (uInt iRow=0; iRow<nRows; iRow++) {
    272 
    273 // Check conformance
    274 
     254       // Check conformance
    275255        IPosition shp2 = in[iTab]->rowAsMaskedArray(iRow).shape();
    276256        if (!shp.isEqual(shp2)) {
     257          delete pTabOut;
    277258           throw (AipsError("Shapes for all rows must be the same"));
    278259        }
    279260
    280 // If we are not doing scan averages, make checks for source and
    281 // frequency setup and warn if averaging across them
    282 
     261        // If we are not doing scan averages, make checks for source
     262        // and frequency setup and warn if averaging across them
    283263        scanIDCol.getScalar(iRow, scanID);
    284264
    285 // Get quantities from columns
    286 
     265        // Get quantities from columns
    287266        srcNameCol.getScalar(iRow, sourceName);
    288267        mjdCol.get(iRow, time);
     
    291270        fqIDCol.get(iRow, freqID);
    292271
    293 // Initialize first source and freqID
    294 
     272        // Initialize first source and freqID
    295273        if (iRow==0 && iTab==0) {
    296274          sourceNameStart = sourceName;
     
    298276        }
    299277
    300 // If we are doing scan averages, see if we are at the end of an
    301 // accumulation period (scan).  We must check soutce names too,
    302 // since we might have two tables with one scan each but different
    303 // source names; we shouldn't average different sources together
    304 
     278        // If we are doing scan averages, see if we are at the end of
     279        // an accumulation period (scan).  We must check soutce names
     280        // too, since we might have two tables with one scan each but
     281        // different source names; we shouldn't average different
     282        // sources together
    305283        if (scanAv && ( (scanID != oldScanID)  ||
    306284                        (iRow==0 && iTab>0 && sourceName!=oldSourceName))) {
    307285
    308 // Normalize data in 'sum' accumulation array according to weighting scheme
    309 
    310            normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, asap::ChanAxis, nAxesSub);
    311 
    312 // Get ScanContainer for the first row of this averaged Scan
    313 
     286          // Normalize data in 'sum' accumulation array according to
     287          // weighting scheme
     288           normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType,
     289                     asap::ChanAxis, nAxesSub);
     290
     291           // Get ScanContainer for the first row of this averaged Scan
    314292           SDContainer scOut = in[iTab]->getSDContainer(rowStart);
    315293
    316 // Fill scan container. The source and freqID come from the
    317 // first row of the first table that went into this average (
    318 // should be the same for all rows in the scan average)
     294           // Fill scan container. The source and freqID come from the
     295           // first row of the first table that went into this average
     296           // ( should be the same for all rows in the scan average)
    319297
    320298           Float nR(nAccum);
    321299           fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
    322                     timeSum/nR, intSum, sourceNameStart, freqIDStart);
    323 
    324 // Write container out to Table
    325 
     300                   timeSum/nR, intSum, sourceNameStart, freqIDStart);
     301           
     302           // Write container out to Table
    326303           pTabOut->putSDContainer(scOut);
    327 
    328 // Reset accumulators
    329 
     304           
     305           // Reset accumulators           
    330306           sum = 0.0;
    331307           sumSq = 0.0;
    332308           nAccum = 0;
    333 //
     309
    334310           tSysSum =0.0;
    335311           tSysSqSum =0.0;
     
    338314           nPts = 0.0;
    339315
    340 // Increment
    341 
     316           // Increment
    342317           rowStart = iRow;              // First row for next accumulation
    343318           tableStart = iTab;            // First table for next accumulation
    344            sourceNameStart = sourceName; // First source name for next accumulation
     319           sourceNameStart = sourceName; // First source name for next
     320                                         // accumulation
    345321           freqIDStart = freqID;         // First FreqID for next accumulation
    346 //
     322
    347323           oldScanID = scanID;
    348            outScanID += 1;               // Scan ID for next accumulation period
     324           outScanID += 1;               // Scan ID for next
     325                                         // accumulation period
    349326        }
    350327
    351 // Accumulate
    352 
    353         accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum, tSysSqSum,
    354                    tSys, nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis,
     328        // Accumulate
     329        accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts,
     330                   tSysSum, tSysSqSum, tSys,
     331                  nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis,
    355332                   nAxesSub, useMask, wtType);
    356 //
    357        oldSourceName = sourceName;
    358        oldFreqID = freqID;
     333        oldSourceName = sourceName;
     334        oldFreqID = freqID;
    359335     }
    360336  }
    361337
    362 // OK at this point we have accumulation data which is either
    363 //   - accumulated from all tables into one row
    364 // or
    365 //   - accumulated from the last scan average
    366 //
    367 // Normalize data in 'sum' accumulation array according to weighting scheme
    368 
    369   normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, asap::ChanAxis, nAxesSub);
    370 
    371 // Create and fill container.  The container we clone will be from
    372 // the last Table and the first row that went into the current
    373 // accumulation.  It probably doesn't matter that much really...
    374 
     338  // OK at this point we have accumulation data which is either
     339  //   - accumulated from all tables into one row
     340  // or
     341  //   - accumulated from the last scan average
     342  //
     343  // Normalize data in 'sum' accumulation array according to weighting
     344  // scheme
     345
     346  normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType,
     347            asap::ChanAxis, nAxesSub);
     348
     349  // Create and fill container.  The container we clone will be from
     350  // the last Table and the first row that went into the current
     351  // accumulation.  It probably doesn't matter that much really...
    375352  Float nR(nAccum);
    376353  SDContainer scOut = in[tableStart]->getSDContainer(rowStart);
    377354  fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
    378            timeSum/nR, intSum, sourceNameStart, freqIDStart);
     355          timeSum/nR, intSum, sourceNameStart, freqIDStart);
    379356  pTabOut->putSDContainer(scOut);
    380357  pTabOut->resetCursor();
    381 //
     358
    382359  return CountedPtr<SDMemTable>(pTabOut);
    383360}
     
    385362
    386363
    387 CountedPtr<SDMemTable> SDMath::binaryOperate(const CountedPtr<SDMemTable>&
    388                                              left,
    389                                              const CountedPtr<SDMemTable>&
    390                                              right,
    391                                              const String& op, Bool preserve,
    392                                              Bool doTSys) const
    393 {
    394 
    395 // Check operator
    396 
     364CountedPtr<SDMemTable>
     365SDMath::binaryOperate(const CountedPtr<SDMemTable>& left,
     366                      const CountedPtr<SDMemTable>& right,
     367                      const String& op, Bool preserve, Bool doTSys) const
     368{
     369
     370  // Check operator
    397371  String op2(op);
    398372  op2.upcase();
     
    413387  }
    414388
    415 // Check rows
    416 
     389  // Check rows
    417390  const uInt nRowLeft = left->nRow();
    418391  const uInt nRowRight = right->nRow();
    419   Bool ok = (nRowRight==1&&nRowLeft>0) ||
    420             (nRowRight>=1&&nRowLeft==nRowRight);
     392  Bool ok = (nRowRight==1 && nRowLeft>0) ||
     393            (nRowRight>=1 && nRowLeft==nRowRight);
    421394  if (!ok) {
    422395     throw (AipsError("The right Scan Table can have one row or the same number of rows as the left Scan Table"));
    423396  }
    424397
    425 // Input Tables
    426 
     398  // Input Tables
    427399  const Table& tLeft = left->table();
    428400  const Table& tRight = right->table();
    429401
    430 // TSys columns
    431 
     402  // TSys columns
    432403  ROArrayColumn<Float> tSysLeftCol, tSysRightCol;
    433404  if (doTSys) {
    434      tSysLeftCol.attach(tLeft, "TSYS");
    435      tSysRightCol.attach(tRight, "TSYS");
    436   }
    437 
    438 // First row for right
    439 
     405    tSysLeftCol.attach(tLeft, "TSYS");
     406    tSysRightCol.attach(tRight, "TSYS");
     407  }
     408
     409  // First row for right
    440410  Array<Float> tSysLeftArr, tSysRightArr;
    441411  if (doTSys) tSysRightCol.get(0, tSysRightArr);
    442   MaskedArray<Float>* pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(0));
     412  MaskedArray<Float>* pMRight =
     413    new MaskedArray<Float>(right->rowAsMaskedArray(0));
     414
    443415  IPosition shpRight = pMRight->shape();
    444416
    445 // Output Table cloned from left
    446 
     417  // Output Table cloned from left
    447418  SDMemTable* pTabOut = new SDMemTable(*left, True);
    448419  pTabOut->appendToHistoryTable(right->getHistoryTable());
    449 // Loop over rows
    450 
     420
     421  // Loop over rows
    451422  for (uInt i=0; i<nRowLeft; i++) {
    452 
    453 // Get data
    454 
    455      MaskedArray<Float> mLeft(left->rowAsMaskedArray(i));
    456      IPosition shpLeft = mLeft.shape();
    457      if (doTSys) tSysLeftCol.get(i, tSysLeftArr);
    458 //
    459      if (nRowRight>1) {
    460         delete pMRight;
    461         pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(i));
    462         shpRight = pMRight->shape();
    463         if (doTSys) tSysRightCol.get(i, tSysRightArr);
    464      }
    465 //
    466      if (!shpRight.isEqual(shpLeft)) {
    467         throw(AipsError("left and right scan tables are not conformant"));
    468      }
    469      if (doTSys) {
    470         if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) {
    471            throw(AipsError("left and right Tsys data are not conformant"));
    472         }
    473         if (!shpRight.isEqual(tSysRightArr.shape())) {
    474            throw(AipsError("left and right scan tables are not conformant"));
    475         }
    476      }
    477 
    478 // Make container
    479 
     423   
     424    // Get data
     425    MaskedArray<Float> mLeft(left->rowAsMaskedArray(i));
     426    IPosition shpLeft = mLeft.shape();
     427    if (doTSys) tSysLeftCol.get(i, tSysLeftArr);
     428   
     429    if (nRowRight>1) {
     430      delete pMRight;
     431      pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(i));
     432      shpRight = pMRight->shape();
     433      if (doTSys) tSysRightCol.get(i, tSysRightArr);
     434    }
     435
     436    if (!shpRight.isEqual(shpLeft)) {
     437      delete pTabOut;
     438      delete pMRight;
     439      throw(AipsError("left and right scan tables are not conformant"));
     440    }
     441    if (doTSys) {
     442      if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) {
     443        delete pTabOut;
     444        delete pMRight;
     445        throw(AipsError("left and right Tsys data are not conformant"));
     446      }
     447      if (!shpRight.isEqual(tSysRightArr.shape())) {
     448        delete pTabOut;
     449        delete pMRight;
     450        throw(AipsError("left and right scan tables are not conformant"));
     451      }
     452    }
     453
     454    // Make container
    480455     SDContainer sc = left->getSDContainer(i);
    481456
    482 // Operate on data and TSys
    483 
     457     // Operate on data and TSys
    484458     if (what==0) {                               
    485459        MaskedArray<Float> tmp = mLeft + *pMRight;
     
    511485     }
    512486
    513 // Put new row in output Table
    514 
     487     // Put new row in output Table
    515488     pTabOut->putSDContainer(sc);
    516489  }
    517490  if (pMRight) delete pMRight;
    518491  pTabOut->resetCursor();
    519 
     492 
    520493  return CountedPtr<SDMemTable>(pTabOut);
    521494}
    522 
    523495
    524496
     
    828800  shapeOut(asap::PolAxis) = 1;                          // Average all polarizations
    829801  if (shapeIn(asap::PolAxis)==1) {
    830      throw(AipsError("The input has only one polarisation"));
     802    delete  pTabOut;
     803    throw(AipsError("The input has only one polarisation"));
    831804  }
    832805//
     
    876849      ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
    877850      ReadOnlyArrayIterator<Float>* pItTsysPlane = 0;
    878       if (wtType==TSYS) pItTsysPlane = new ReadOnlyArrayIterator<Float>(tSys, axes);
     851      if (wtType==TSYS)
     852        pItTsysPlane = new ReadOnlyArrayIterator<Float>(tSys, axes);
    879853
    880854// Accumulations
     
    900874           ReadOnlyVectorIterator<Float>* pItTsysVec = 0;
    901875           if (wtType==TSYS) {
    902               pItTsysVec = new ReadOnlyVectorIterator<Float>(pItTsysPlane->array(), 1);
     876              pItTsysVec =
     877                new ReadOnlyVectorIterator<Float>(pItTsysPlane->array(), 1);
    903878           }             
    904879//
     
    908883
    909884              if (useMask) {
    910                  const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
     885                 const MaskedArray<Float> spec(itDataVec.vector(),
     886                                               mask&&itMaskVec.vector());
    911887                 if (wtType==VAR) {
    912888                    fac = 1.0 / variance(spec);
     
    916892                 }                   
    917893              } else {
    918                  const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
     894                 const MaskedArray<Float> spec(itDataVec.vector(),
     895                                               itMaskVec.vector());
    919896                 if (wtType==VAR) {
    920897                    fac = 1.0 / variance(spec);
     
    927904// Normalize spectrum (without OTF mask) and accumulate
    928905
    929               const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
     906              const MaskedArray<Float> spec(fac*itDataVec.vector(),
     907                                            itMaskVec.vector());
    930908              vecSum += spec;
    931909              norm += fac;
     
    1001979{
    1002980
    1003 // Number of channels
    1004 
     981  // Number of channels
    1005982   const uInt nChan = in.nChan();
    1006983
    1007 // Generate Kernel
    1008 
     984   // Generate Kernel
    1009985   VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernelType);
    1010986   Vector<Float> kernel = VectorKernel::make(type, width, nChan, True, False);
    1011987
    1012 // Generate Convolver
    1013 
     988   // Generate Convolver
    1014989   IPosition shape(1,nChan);
    1015990   Convolver<Float> conv(kernel, shape);
    1016991
    1017 // New Table
    1018 
     992   // New Table
    1019993   SDMemTable* pTabOut = new SDMemTable(in,True);
    1020994
    1021 // Output Vectors
    1022 
     995   // Output Vectors
    1023996   Vector<Float> valuesOut(nChan);
    1024997   Vector<Bool> maskOut(nChan);
    1025998
    1026 // Get data slice bounds
    1027 
     999   // Get data slice bounds
    10281000   IPosition start, end;
    10291001   setCursorSlice (start, end, doAll, in);
    10301002
    1031 // Loop over rows in Table
    1032 
     1003   // Loop over rows in Table
    10331004   for (uInt ri=0; ri < in.nRow(); ++ri) {
    10341005
    1035 // Get slice of data
    1036 
     1006     // Get slice of data
    10371007      MaskedArray<Float> dataIn = in.rowAsMaskedArray(ri);
    10381008
    1039 // Deconstruct and get slices which reference these arrays
    1040 
     1009      // Deconstruct and get slices which reference these arrays
    10411010      Array<Float> valuesIn = dataIn.getArray();
    10421011      Array<Bool> maskIn = dataIn.getMask();
    1043 //
     1012
    10441013      Array<Float> valuesIn2 = valuesIn(start,end);       // ref to valuesIn
    10451014      Array<Bool> maskIn2 = maskIn(start,end);
    10461015
    1047 // Iterate through by spectra
    1048 
     1016      // Iterate through by spectra
    10491017      VectorIterator<Float> itValues(valuesIn2, asap::ChanAxis);
    10501018      VectorIterator<Bool> itMask(maskIn2, asap::ChanAxis);
    10511019      while (!itValues.pastEnd()) {
    1052        
    1053 // Smooth
    1054 
    1055          if (kernelType==VectorKernel::HANNING) {
    1056             mathutil::hanning(valuesOut, maskOut, itValues.vector(), itMask.vector());
    1057             itMask.vector() = maskOut;
    1058          } else {
    1059             mathutil::replaceMaskByZero(itValues.vector(), itMask.vector());
    1060             conv.linearConv(valuesOut, itValues.vector());
    1061          }
    1062 //   
    1063          itValues.vector() = valuesOut;
    1064 //
    1065          itValues.next();
    1066          itMask.next();
     1020       
     1021        // Smooth
     1022        if (kernelType==VectorKernel::HANNING) {
     1023          mathutil::hanning(valuesOut, maskOut, itValues.vector(),
     1024                            itMask.vector());
     1025          itMask.vector() = maskOut;
     1026        } else {
     1027          mathutil::replaceMaskByZero(itValues.vector(), itMask.vector());
     1028          conv.linearConv(valuesOut, itValues.vector());
     1029        }
     1030
     1031        itValues.vector() = valuesOut;
     1032        itValues.next();
     1033        itMask.next();
    10671034      }
    1068 
    1069 // Create and put back
    1070 
     1035     
     1036      // Create and put back
    10711037      SDContainer sc = in.getSDContainer(ri);
    10721038      putDataInSDC(sc, valuesIn, maskIn);
    1073 //
     1039
    10741040      pTabOut->putSDContainer(sc);
    10751041   }
    1076 //
     1042
    10771043  return pTabOut;
    10781044}
     
    10911057  SDMemTable* pTabOut = new SDMemTable(in, True);
    10921058
    1093 // Find out how to convert values into Jy and K (e.g. units might be mJy or mK)
    1094 // Also automatically find out what we are converting to according to the
    1095 // flux unit
    1096 
     1059  // Find out how to convert values into Jy and K (e.g. units might be
     1060  // mJy or mK) Also automatically find out what we are converting to
     1061  // according to the flux unit
    10971062  Unit fluxUnit(sh.fluxunit);
    10981063  Unit K(String("K"));
    10991064  Unit JY(String("Jy"));
    1100 //
     1065
    11011066  Bool toKelvin = True;
    11021067  Double cFac = 1.0;   
    11031068  if (fluxUnit==JY) {
    1104      cout << "Converting to K" << endl;
    1105 //
    1106      Quantum<Double> t(1.0,fluxUnit);
    1107      Quantum<Double> t2 = t.get(JY);
    1108      cFac = (t2 / t).getValue();               // value to Jy
    1109 //
    1110      toKelvin = True;
    1111      sh.fluxunit = "K";
     1069    cout << "Converting to K" << endl;
     1070   
     1071    Quantum<Double> t(1.0,fluxUnit);
     1072    Quantum<Double> t2 = t.get(JY);
     1073    cFac = (t2 / t).getValue();               // value to Jy
     1074   
     1075    toKelvin = True;
     1076    sh.fluxunit = "K";
    11121077  } else if (fluxUnit==K) {
    1113      cout << "Converting to Jy" << endl;
    1114 //
    1115      Quantum<Double> t(1.0,fluxUnit);
    1116      Quantum<Double> t2 = t.get(K);
    1117      cFac = (t2 / t).getValue();              // value to K
    1118 //
    1119      toKelvin = False;
    1120      sh.fluxunit = "Jy";
     1078    cout << "Converting to Jy" << endl;
     1079   
     1080    Quantum<Double> t(1.0,fluxUnit);
     1081    Quantum<Double> t2 = t.get(K);
     1082    cFac = (t2 / t).getValue();              // value to K
     1083   
     1084    toKelvin = False;
     1085    sh.fluxunit = "Jy";
    11211086  } else {
    1122      throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
     1087    throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
    11231088  }
    11241089  pTabOut->putSDHeader(sh);
    1125 
    1126 // Make sure input values are converted to either Jy or K first...
    1127 
     1090 
     1091  // Make sure input values are converted to either Jy or K first...
    11281092  Float factor = cFac;
    11291093
    1130 // Select method
    1131 
     1094  // Select method
    11321095  if (JyPerK>0.0) {
    1133      factor *= JyPerK;
    1134      if (toKelvin) factor = 1.0 / JyPerK;
    1135 //
    1136      cout << "Jy/K = " << JyPerK << endl;
    1137      Vector<Float> factors(in.nRow(), factor);
    1138      scaleByVector(pTabOut, in, doAll, factors, False);
     1096    factor *= JyPerK;
     1097    if (toKelvin) factor = 1.0 / JyPerK;
     1098    cout << "Jy/K = " << JyPerK << endl;
     1099    Vector<Float> factors(in.nRow(), factor);
     1100    scaleByVector(pTabOut, in, doAll, factors, False);
    11391101  } else if (etaAp>0.0) {
    1140      Bool throwIt = True;
    1141      Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
    1142      SDAttr sda;
    1143      if (D < 0) D = sda.diameter(inst);
    1144      Float JyPerK = SDAttr::findJyPerK (etaAp,D);
    1145      cout << "Jy/K = " << JyPerK << endl;
    1146      factor *= JyPerK;
    1147      if (toKelvin) {
    1148         factor = 1.0 / factor;
    1149      }
    1150 //
    1151      Vector<Float> factors(in.nRow(), factor);
    1152      scaleByVector(pTabOut, in, doAll, factors, False);
     1102    Bool throwIt = True;
     1103    Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
     1104    SDAttr sda;
     1105    if (D < 0) D = sda.diameter(inst);
     1106    Float JyPerK = SDAttr::findJyPerK (etaAp,D);
     1107    cout << "Jy/K = " << JyPerK << endl;
     1108    factor *= JyPerK;
     1109    if (toKelvin) {
     1110      factor = 1.0 / factor;
     1111    }
     1112
     1113    Vector<Float> factors(in.nRow(), factor);
     1114    scaleByVector(pTabOut, in, doAll, factors, False);
    11531115  } else {
    1154 
    1155 // OK now we must deal with automatic look up of values.
    1156 // We must also deal with the fact that the factors need
    1157 // to be computed per IF and may be different and may
    1158 // change per integration.
    1159 
    1160      cout << "Looking up conversion factors" << endl;
    1161      convertBrightnessUnits (pTabOut, in, toKelvin, cFac, doAll);
    1162   }
    1163 //
     1116   
     1117    // OK now we must deal with automatic look up of values.
     1118    // We must also deal with the fact that the factors need
     1119    // to be computed per IF and may be different and may
     1120    // change per integration.
     1121   
     1122    cout << "Looking up conversion factors" << endl;
     1123    convertBrightnessUnits (pTabOut, in, toKelvin, cFac, doAll);
     1124  }
     1125
    11641126  return pTabOut;
    11651127}
    1166 
    1167 
    1168 
    11691128
    11701129
     
    11751134{
    11761135
    1177 // Get header and clone output table
    1178 
     1136  // Get header and clone output table
    11791137  SDHeader sh = in.getSDHeader();
    11801138  SDMemTable* pTabOut = new SDMemTable(in, True);
    11811139
    1182 // Get elevation data from SDMemTable and convert to degrees
    1183 
     1140  // Get elevation data from SDMemTable and convert to degrees
    11841141  const Table& tab = in.table();
    11851142  ROScalarColumn<Float> elev(tab, "ELEVATION");
    11861143  Vector<Float> x = elev.getColumn();
    11871144  x *= Float(180 / C::pi);                        // Degrees
    1188 //
     1145 
    11891146  const uInt nC = coeffs.nelements();
    11901147  if (fileName.length()>0 && nC>0) {
    1191      throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
    1192   }
    1193 
    1194 // Correct
    1195 
     1148    throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
     1149  }
     1150 
     1151  // Correct
    11961152  if (nC>0 || fileName.length()==0) {
    1197 
    1198 // Find instrument
    1199 
     1153    // Find instrument
    12001154     Bool throwIt = True;
    12011155     Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
    12021156     
    1203 // Set polynomial
    1204 
     1157     // Set polynomial
    12051158     Polynomial<Float>* pPoly = 0;
    12061159     Vector<Float> coeff;
    12071160     String msg;
    12081161     if (nC>0) {
    1209         pPoly = new Polynomial<Float>(nC);
    1210         coeff = coeffs;
    1211         msg = String("user");
     1162       pPoly = new Polynomial<Float>(nC);
     1163       coeff = coeffs;
     1164       msg = String("user");
    12121165     } else {
    1213         SDAttr sdAttr;
    1214         coeff = sdAttr.gainElevationPoly(inst);
    1215         pPoly = new Polynomial<Float>(3);
    1216         msg = String("built in");
     1166       SDAttr sdAttr;
     1167       coeff = sdAttr.gainElevationPoly(inst);
     1168       pPoly = new Polynomial<Float>(3);
     1169       msg = String("built in");
    12171170     }
    1218 //
     1171     
    12191172     if (coeff.nelements()>0) {
    1220         pPoly->setCoefficients(coeff);
     1173       pPoly->setCoefficients(coeff);
    12211174     } else {
    1222         throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
     1175       delete pPoly;
     1176       throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
    12231177     }
    1224 //
    12251178     cout << "Making polynomial correction with " << msg << " coefficients" << endl;
    12261179     const uInt nRow = in.nRow();
    12271180     Vector<Float> factor(nRow);
    12281181     for (uInt i=0; i<nRow; i++) {
    1229         factor[i] = 1.0 / (*pPoly)(x[i]);
     1182       factor[i] = 1.0 / (*pPoly)(x[i]);
    12301183     }
    12311184     delete pPoly;
    1232 //
    12331185     scaleByVector (pTabOut, in, doAll, factor, True);
     1186
    12341187  } else {
    1235 
    1236 // Indicate which columns to read from ascii file
    1237 
    1238      String col0("ELEVATION");
    1239      String col1("FACTOR");
    1240 
    1241 // Read and correct
    1242 
    1243      cout << "Making correction from ascii Table" << endl;
    1244      scaleFromAsciiTable (pTabOut, in, fileName, col0, col1,
    1245                           methodStr, doAll, x, True);
    1246    }
    1247 //
    1248    return pTabOut;
    1249 }
    1250 
     1188   
     1189    // Indicate which columns to read from ascii file
     1190    String col0("ELEVATION");
     1191    String col1("FACTOR");
     1192   
     1193    // Read and correct
     1194   
     1195    cout << "Making correction from ascii Table" << endl;
     1196    scaleFromAsciiTable (pTabOut, in, fileName, col0, col1,
     1197                         methodStr, doAll, x, True);
     1198  }
     1199
     1200  return pTabOut;
     1201}
    12511202 
    12521203
     
    12541205{
    12551206
    1256 // Get header and clone output table
     1207  // Get header and clone output table
    12571208
    12581209  SDHeader sh = in.getSDHeader();
     
    12781229
    12791230  scaleByVector (pTabOut, in, doAll, factor, True);
    1280 //
     1231
    12811232  return pTabOut;
    12821233}
     
    12891240//
    12901241{
    1291    if (in.nPol() != 4) {
    1292       throw(AipsError("You must have 4 polarizations to run this function"));
    1293    }
    1294 //
     1242  if (in.nPol() != 4) {
     1243    throw(AipsError("You must have 4 polarizations to run this function"));
     1244  }
     1245
    12951246   SDHeader sh = in.getSDHeader();
    12961247   Instrument inst = SDAttr::convertInstrument (sh.antennaname, False);
     
    13241275      IPosition shape = data.shape();
    13251276 
    1326 // Get polarization slice references
    1327  
     1277      // Get polarization slice references
    13281278      Array<Float> C3 = data(start3,end3);
    13291279      Array<Float> C4 = data(start4,end4);
    13301280   
    1331 // Rotate
    1332  
     1281      // Rotate
    13331282      SDPolUtil::rotatePhase(C3, C4, value);
    13341283   
    1335 // Put
    1336    
     1284      // Put
    13371285      specCol.put(i,data);
    13381286   }
    13391287}     
    1340 
    13411288
    13421289
     
    15231470
    15241471
    1525 SDMemTable* SDMath::frequencyAlign (const SDMemTable& in,
     1472SDMemTable* SDMath::frequencyAlign(const SDMemTable& in,
    15261473                                   MFrequency::Types freqSystem,
    15271474                                   const String& refTime,
     
    15511498   Matrix<uInt> ddIdx;
    15521499   SDDataDesc dDesc;
    1553    generateDataDescTable (ddIdx, dDesc, nIF, in, tabIn, srcCol, fqIDCol, perFreqID);
     1500   generateDataDescTable(ddIdx, dDesc, nIF, in, tabIn, srcCol,
     1501                          fqIDCol, perFreqID);
    15541502
    15551503// Get reference Epoch to time of first row or given String
    1556 
     1504   
    15571505   Unit DAY(String("d"));
    15581506   MEpoch::Ref epochRef(in.getTimeReference());
    15591507   MEpoch refEpoch;
    15601508   if (refTime.length()>0) {
    1561       refEpoch = epochFromString(refTime, in.getTimeReference());
     1509     refEpoch = epochFromString(refTime, in.getTimeReference());
    15621510   } else {
    1563       refEpoch = in.getEpoch(0);
     1511     refEpoch = in.getEpoch(0);
    15641512   }
    15651513   cout << "Aligning at reference Epoch " << formatEpoch(refEpoch)
    15661514        << " in frame " << MFrequency::showType(freqSystem) << endl;
    15671515   
    1568 // Get Reference Position
    1569 
     1516   // Get Reference Position
     1517   
    15701518   MPosition refPos = in.getAntennaPosition();
    1571 
    1572 // Create FrequencyAligner Block. One FA for each possible
    1573 // source/freqID (perFreqID=True) or source/IF (perFreqID=False) combination
    1574 
     1519   
     1520   // Create FrequencyAligner Block. One FA for each possible
     1521   // source/freqID (perFreqID=True) or source/IF (perFreqID=False)
     1522   // combination
     1523   
    15751524   PtrBlock<FrequencyAligner<Float>* > a(dDesc.length());
    15761525   generateFrequencyAligners (a, dDesc, in, nChan, freqSystem, refPos,
    15771526                              refEpoch, perFreqID);
    1578 
    1579 // Generate and fill output Frequency Table.  WHen perFreqID=True, there is one output FreqID
    1580 // for each entry in the SDDataDesc table.  However, in perFreqID=False mode, there may be
    1581 // some degeneracy, so we need a little translation map
    1582 
     1527   
     1528   // Generate and fill output Frequency Table.  WHen perFreqID=True,
     1529   // there is one output FreqID for each entry in the SDDataDesc
     1530   // table.  However, in perFreqID=False mode, there may be some
     1531   // degeneracy, so we need a little translation map
     1532   
    15831533   SDFrequencyTable freqTabOut = in.getSDFreqTable();
    15841534   freqTabOut.setLength(0);
     
    15861536   units = String("Hz");
    15871537   Bool linear=True;
    1588 //
     1538   //
    15891539   Vector<uInt> ddFQTrans(dDesc.length(),0);
    15901540   for (uInt i=0; i<dDesc.length(); i++) {
    15911541
    1592 // Get Aligned SC in Hz
    1593 
    1594       SpectralCoordinate sC = a[i]->alignedSpectralCoordinate(linear);
    1595       sC.setWorldAxisUnits(units);
    1596 
    1597 // Add FreqID
    1598 
    1599       uInt idx = freqTabOut.addFrequency(sC.referencePixel()[0],
    1600                                          sC.referenceValue()[0],
    1601                                          sC.increment()[0]);
    1602       ddFQTrans(i) = idx;                                       // output FreqID = ddFQTrans(ddIdx)
     1542     // Get Aligned SC in Hz
     1543     
     1544     SpectralCoordinate sC = a[i]->alignedSpectralCoordinate(linear);
     1545     sC.setWorldAxisUnits(units);
     1546     
     1547     // Add FreqID
     1548     
     1549     uInt idx = freqTabOut.addFrequency(sC.referencePixel()[0],
     1550                                        sC.referenceValue()[0],
     1551                                        sC.increment()[0]);
     1552     // output FreqID = ddFQTrans(ddIdx)
     1553     ddFQTrans(i) = idx;
    16031554   }
    1604 
    1605 // Interpolation method
    1606 
     1555   
     1556   // Interpolation method
     1557   
    16071558   InterpolateArray1D<Double,Float>::InterpolationMethod interp;
    16081559   convertInterpString(interp, methodStr);
    1609 
    1610 // New output Table
    1611 
     1560   
     1561   // New output Table
     1562   
    16121563   cout << "Create output table" << endl;
    16131564   SDMemTable* pTabOut = new SDMemTable(in,True);
    16141565   pTabOut->putSDFreqTable(freqTabOut);
    1615 
    1616 // Loop over rows in Table
    1617 
     1566   
     1567   // Loop over rows in Table
     1568   
    16181569   Bool extrapolate=False;
    16191570   const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis);
     
    16261577   uInt ifIdx, faIdx;
    16271578   Vector<Double> xIn;
    1628 //
     1579   //
    16291580   for (uInt iRow=0; iRow<nRows; ++iRow) {
    1630       if (iRow%10==0) {
    1631          cout << "Processing row " << iRow << endl;
    1632       }
    1633 
    1634 // Get EPoch
    1635 
     1581     if (iRow%10==0) {
     1582       cout << "Processing row " << iRow << endl;
     1583     }
     1584     
     1585     // Get EPoch
     1586     
    16361587     Quantum<Double> tQ2(times[iRow],DAY);
    16371588     MVEpoch mv2(tQ2);
    16381589     MEpoch epoch(mv2, epochRef);
    1639 
    1640 // Get copy of data
    1641    
     1590     
     1591     // Get copy of data
     1592     
    16421593     const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow));
    16431594     Array<Float> values = mArrIn.getArray();
    16441595     Array<Bool> mask = mArrIn.getMask();
    1645 
    1646 // For each row, the Frequency abcissa will be the same regardless
    1647 // of polarization.  For all other axes (IF and BEAM) the abcissa
    1648 // will change.  So we iterate through the data by pol-chan planes
    1649 // to mimimize the work.  Probably won't work for multiple beams
    1650 // at this point.
    1651 
     1596     
     1597     // For each row, the Frequency abcissa will be the same
     1598     // regardless of polarization.  For all other axes (IF and BEAM)
     1599     // the abcissa will change.  So we iterate through the data by
     1600     // pol-chan planes to mimimize the work.  Probably won't work for
     1601     // multiple beams at this point.
     1602     
    16521603     ArrayIterator<Float> itValuesPlane(values, polChanAxes);
    16531604     ArrayIterator<Bool> itMaskPlane(mask, polChanAxes);
    16541605     while (!itValuesPlane.pastEnd()) {
    16551606
    1656 // Find the IF index and then the FA PtrBlock index
    1657 
    1658         const IPosition& pos = itValuesPlane.pos();
    1659         ifIdx = pos(asap::IFAxis);
    1660         faIdx = ddIdx(iRow,ifIdx);
    1661 
    1662 // Generate abcissa for perIF.  Could cache this in a Matrix
    1663 // on a per scan basis.   Pretty expensive doing it for every row.
    1664 
    1665         if (!perFreqID) {
    1666           xIn.resize(nChan);
    1667           uInt fqID = dDesc.secID(ddIdx(iRow,ifIdx));
    1668           SpectralCoordinate sC = in.getSpectralCoordinate(fqID);
     1607       // Find the IF index and then the FA PtrBlock index
     1608       
     1609       const IPosition& pos = itValuesPlane.pos();
     1610       ifIdx = pos(asap::IFAxis);
     1611       faIdx = ddIdx(iRow,ifIdx);
     1612       
     1613       // Generate abcissa for perIF.  Could cache this in a Matrix on
     1614       // a per scan basis.  Pretty expensive doing it for every row.
     1615       
     1616       if (!perFreqID) {
     1617        xIn.resize(nChan);
     1618        uInt fqID = dDesc.secID(ddIdx(iRow,ifIdx));
     1619        SpectralCoordinate sC = in.getSpectralCoordinate(fqID);
    16691620           Double w;
    16701621           for (uInt i=0; i<nChan; i++) {
    1671               sC.toWorld(w,Double(i));
     1622             sC.toWorld(w,Double(i));
    16721623              xIn[i] = w;
    16731624           }
    1674         }
    1675 //
    1676         VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1);
    1677         VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
    1678 
    1679 // Iterate through the plane by vector and align
    1680 
     1625       }
     1626       
     1627       VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1);
     1628       VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
     1629
     1630       // Iterate through the plane by vector and align
     1631       
    16811632        first = True;
    16821633        useCachedAbcissa=False;
    16831634        while (!itValuesVec.pastEnd()) {     
    1684            if (perFreqID) {
    1685               ok = a[faIdx]->align (yOut, maskOut, itValuesVec.vector(),
    1686                                     itMaskVec.vector(), epoch, useCachedAbcissa,
    1687                                     interp, extrapolate);
    1688            } else {
    1689               ok = a[faIdx]->align (yOut, maskOut, xIn, itValuesVec.vector(),
    1690                                     itMaskVec.vector(), epoch, useCachedAbcissa,
    1691                                     interp, extrapolate);
    1692            }
    1693 //
    1694            itValuesVec.vector() = yOut;
    1695            itMaskVec.vector() = maskOut;
    1696 //
    1697            itValuesVec.next();
    1698            itMaskVec.next();
    1699 //
    1700            if (first) {
    1701               useCachedAbcissa = True;
    1702               first = False;
    1703            }
     1635          if (perFreqID) {
     1636            ok = a[faIdx]->align (yOut, maskOut, itValuesVec.vector(),
     1637                                  itMaskVec.vector(), epoch, useCachedAbcissa,
     1638                                  interp, extrapolate);
     1639          } else {
     1640            ok = a[faIdx]->align (yOut, maskOut, xIn, itValuesVec.vector(),
     1641                                  itMaskVec.vector(), epoch, useCachedAbcissa,
     1642                                  interp, extrapolate);
     1643          }
     1644          //
     1645          itValuesVec.vector() = yOut;
     1646          itMaskVec.vector() = maskOut;
     1647          //
     1648          itValuesVec.next();
     1649          itMaskVec.next();
     1650          //
     1651          if (first) {
     1652            useCachedAbcissa = True;
     1653            first = False;
     1654          }
    17041655        }
    1705 //
    1706        itValuesPlane.next();
    1707        itMaskPlane.next();
     1656        //
     1657        itValuesPlane.next();
     1658        itMaskPlane.next();
    17081659     }
    1709 
    1710 // Create SDContainer and put back
    1711 
    1712     SDContainer sc = in.getSDContainer(iRow);
    1713     putDataInSDC(sc, values, mask);
    1714 
    1715 // Set output FreqIDs
    1716 
    1717     for (uInt i=0; i<nIF; i++) {
     1660     
     1661     // Create SDContainer and put back
     1662     
     1663     SDContainer sc = in.getSDContainer(iRow);
     1664     putDataInSDC(sc, values, mask);
     1665     
     1666     // Set output FreqIDs
     1667     
     1668     for (uInt i=0; i<nIF; i++) {
    17181669       uInt idx = ddIdx(iRow,i);               // Index into SDDataDesc table
    17191670       freqID(i) = ddFQTrans(idx);             // FreqID in output FQ table
    1720     }
    1721     sc.putFreqMap(freqID);
    1722 //
    1723     pTabOut->putSDContainer(sc);
     1671     }
     1672     sc.putFreqMap(freqID);
     1673     //
     1674     pTabOut->putSDContainer(sc);
    17241675   }
    1725 
    1726 // Now we must set the base and extra frames to the
    1727 // input frame
    1728 
     1676   
     1677   // Now we must set the base and extra frames to the input frame
    17291678   std::vector<string> info = pTabOut->getCoordInfo();
    17301679   info[1] = MFrequency::showType(freqSystem);   // Conversion frame
     
    17321681   pTabOut->setCoordInfo(info);
    17331682
    1734 // Clean up PointerBlock
    1735 
     1683   // Clean up PointerBlock   
    17361684   for (uInt i=0; i<a.nelements(); i++) delete a[i];
    1737 //
     1685
    17381686   return pTabOut;
    17391687}
    17401688
     1689
     1690SDMemTable* SDMath::frequencySwitch(const SDMemTable& in) const
     1691{
     1692  if (in.nIF() != 2) {
     1693    throw(AipsError("nIF != 2 "));
     1694  }
     1695  Bool clear = True;
     1696  SDMemTable* pTabOut = new SDMemTable(in, clear);
     1697  const Table& tabIn = in.table();
     1698
     1699  // Shape of input and output data
     1700  const IPosition& shapeIn = in.rowAsMaskedArray(0).shape();
     1701
     1702  const uInt nRows = in.nRow();
     1703  AlwaysAssert(asap::nAxes==4,AipsError);
     1704
     1705  ROArrayColumn<Float> tSysCol;
     1706  Array<Float> tsys;
     1707  tSysCol.attach(tabIn,"TSYS");
     1708
     1709  for (uInt iRow=0; iRow<nRows; iRow++) {
     1710    // Get data for this row
     1711    MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
     1712    tSysCol.get(iRow, tsys);
     1713
     1714    // whole Array for IF 0
     1715    IPosition start(asap::nAxes,0);
     1716    IPosition end = shapeIn-1;
     1717    end(asap::IFAxis) = 0;
     1718
     1719    MaskedArray<Float> on = marr(start,end);
     1720    Array<Float> ton = tsys(start,end);
     1721    // Make a copy as "src" is a refrence which is manipulated.
     1722    // oncopy is needed for the inverse quotient
     1723    MaskedArray<Float> oncopy = on.copy();
     1724
     1725    // whole Array for IF 1
     1726    start(asap::IFAxis) = 1;
     1727    end(asap::IFAxis) = 1;
     1728
     1729    MaskedArray<Float> off = marr(start,end);
     1730    Array<Float> toff = tsys(start,end);
     1731
     1732    on /= off; on -= 1.0f;
     1733    on *= ton;
     1734    off /= oncopy; off -= 1.0f;
     1735    off *= toff;
     1736
     1737    SDContainer sc = in.getSDContainer(iRow);
     1738    putDataInSDC(sc, marr.getArray(), marr.getMask());
     1739    pTabOut->putSDContainer(sc);
     1740  }
     1741  return pTabOut;
     1742}
    17411743
    17421744void SDMath::fillSDC(SDContainer& sc,
     
    21632165
    21642166
    2165 MEpoch SDMath::epochFromString (const String& str, MEpoch::Types timeRef) const
     2167MEpoch SDMath::epochFromString(const String& str, MEpoch::Types timeRef) const
    21662168{
    21672169   Quantum<Double> qt;
     
    21842186
    21852187
    2186 void SDMath::generateFrequencyAligners (PtrBlock<FrequencyAligner<Float>* >& a,
    2187                                         const SDDataDesc& dDesc,
    2188                                         const SDMemTable& in, uInt nChan,
    2189                                         MFrequency::Types system,
    2190                                         const MPosition& refPos,
    2191                                         const MEpoch& refEpoch,
    2192                                         Bool perFreqID) const
     2188void SDMath::generateFrequencyAligners(PtrBlock<FrequencyAligner<Float>* >& a,
     2189                                       const SDDataDesc& dDesc,
     2190                                       const SDMemTable& in, uInt nChan,
     2191                                       MFrequency::Types system,
     2192                                       const MPosition& refPos,
     2193                                       const MEpoch& refEpoch,
     2194                                       Bool perFreqID) const
    21932195{
    21942196   for (uInt i=0; i<dDesc.length(); i++) {
     
    22142216}
    22152217
    2216 Vector<uInt> SDMath::getRowRange (const SDMemTable& in) const
     2218Vector<uInt> SDMath::getRowRange(const SDMemTable& in) const
    22172219{
    22182220   Vector<uInt> range(2);
     
    22232225   
    22242226
    2225 Bool SDMath::rowInRange (uInt i, const Vector<uInt>& range) const
     2227Bool SDMath::rowInRange(uInt i, const Vector<uInt>& range) const
    22262228{
    22272229   return (i>=range[0] && i<=range[1]);
    22282230}
     2231
  • branches/Release-2-fixes/src/SDMath.h

    r653 r699  
    6363
    6464// Copy Constructor (copy semantics)
    65    SDMath (const SDMath& other);
     65   SDMath(const SDMath& other);
    6666
    6767// Assignment  (copy semantics)
     
    7272
    7373// Binary Table operators. op=ADD, SUB, MUL, DIV, QUOTIENT
    74    casa::CountedPtr<SDMemTable> binaryOperate (const casa::CountedPtr<SDMemTable>& left,
     74   casa::CountedPtr<SDMemTable> binaryOperate(const casa::CountedPtr<SDMemTable>& left,
    7575                                               const casa::CountedPtr<SDMemTable>& right,
    7676                                               const casa::String& op, casa::Bool preserve,
     
    9999
    100100// Smooth
    101    SDMemTable* smooth (const SDMemTable& in, const casa::String& kernel,
     101   SDMemTable* smooth(const SDMemTable& in, const casa::String& kernel,
    102102                       casa::Float width, casa::Bool doAll) const;
    103103
    104104// Flux conversion between Jansky and Kelvin
    105    SDMemTable* convertFlux (const SDMemTable& in, casa::Float D, casa::Float etaAp,
     105   SDMemTable* convertFlux(const SDMemTable& in, casa::Float D, casa::Float etaAp,
    106106                            casa::Float JyPerK, casa::Bool doAll) const;
    107107
    108108// Gain-elevation correction
    109    SDMemTable* gainElevation (const SDMemTable& in, const casa::Vector<casa::Float>& coeffs,
     109   SDMemTable* gainElevation(const SDMemTable& in, const casa::Vector<casa::Float>& coeffs,
    110110                              const casa::String& fileName,
    111111                              const casa::String& method, casa::Bool doAll) const;
    112112
    113113// Frequency Alignment
    114    SDMemTable* frequencyAlignment (const SDMemTable& in, const casa::String& refTime,
     114   SDMemTable* frequencyAlignment(const SDMemTable& in, const casa::String& refTime,
    115115                                   const casa::String& method, casa::Bool perFreqID) const;
    116116
    117117// Opacity correction
    118    SDMemTable* opacity (const SDMemTable& in, casa::Float tau, casa::Bool doAll) const;
     118   SDMemTable* opacity(const SDMemTable& in, casa::Float tau, casa::Bool doAll) const;
    119119
    120120// Simple unary mathematical operations.  what=0 (mul) or 1 (add)
     
    126126                          const casa::String& wtStr) const;
    127127
     128  SDMemTable* frequencySwitch(const SDMemTable& in) const;
     129
     130
    128131// Rotate XY phase. Value in degrees.
    129    void rotateXYPhase (SDMemTable& in, casa::Float value, casa::Bool doAll);
     132   void rotateXYPhase(SDMemTable& in, casa::Float value, casa::Bool doAll);
    130133
    131134// Rotate Q & U by operating on the raw correlations.Value in degrees.
    132    void rotateLinPolPhase (SDMemTable& in, casa::Float value, casa::Bool doAll);
     135   void rotateLinPolPhase(SDMemTable& in, casa::Float value, casa::Bool doAll);
    133136
    134137
     
    158161
    159162// Work out conversion factor for converting Jy<->K per IF per row and apply
    160    void convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in,
     163   void convertBrightnessUnits(SDMemTable* pTabOut, const SDMemTable& in,
    161164                                casa::Bool toKelvin, casa::Float sFac, casa::Bool doAll) const;
    162165
    163166// Convert weight string to enum value
    164167
    165    void convertWeightString (WeightType& wt, const casa::String& weightStr,
     168   void convertWeightString(WeightType& wt, const casa::String& weightStr,
    166169                             casa::Bool listType) const;
    167170
     
    185188
    186189// Scale data and optionally TSys by values in a Vector
    187    void scaleByVector (SDMemTable* pTabOut, const SDMemTable& in,
     190   void scaleByVector(SDMemTable* pTabOut, const SDMemTable& in,
    188191                       casa::Bool doAll, const casa::Vector<casa::Float>& factor,
    189192                       casa::Bool doTSys) const;
    190193
    191194// Convert time String to Epoch
    192    casa::MEpoch epochFromString (const casa::String& str, casa::MEpoch::Types timeRef) const;
     195   casa::MEpoch epochFromString(const casa::String& str, casa::MEpoch::Types timeRef) const;
    193196
    194197// Function to fill Scan Container when averaging in time
    195198
    196   void fillSDC (SDContainer& sc, const casa::Array<casa::Bool>& mask,
     199  void fillSDC(SDContainer& sc, const casa::Array<casa::Bool>& mask,
    197200                const casa::Array<casa::Float>& data,
    198201                const casa::Array<casa::Float>& tSys,
     
    205208
    206209// Align in Frequency
    207    SDMemTable* frequencyAlign (const SDMemTable& in,
     210   SDMemTable* frequencyAlign(const SDMemTable& in,
    208211                              casa::MFrequency::Types system,
    209212                              const casa::String& timeRef,
     
    212215
    213216// Generate frequency aligners
    214    void generateFrequencyAligners (casa::PtrBlock<casa::FrequencyAligner<casa::Float>* >& a,
     217   void generateFrequencyAligners(casa::PtrBlock<casa::FrequencyAligner<casa::Float>* >& a,
    215218                                   const SDDataDesc& dDesc,
    216219                                   const SDMemTable& in, casa::uInt nChan,
     
    220223                                   casa::Bool perFreqID) const;
    221224
    222 // Generate data description table (combines source and freqID)
    223    void generateDataDescTable (casa::Matrix<casa::uInt>& ddIdx,
     225// Generate data description table(combines source and freqID)
     226   void generateDataDescTable(casa::Matrix<casa::uInt>& ddIdx,
    224227                               SDDataDesc& dDesc,
    225228                               casa::uInt nIF,
     
    231234
    232235// Get row range from SDMemTable state
    233    casa::Vector<casa::uInt> getRowRange (const SDMemTable& in) const;
     236   casa::Vector<casa::uInt> getRowRange(const SDMemTable& in) const;
    234237
    235238// Is row in the row range ?
    236    casa::Bool rowInRange (casa::uInt i, const casa::Vector<casa::uInt>& range) const;
     239   casa::Bool rowInRange(casa::uInt i, const casa::Vector<casa::uInt>& range) const;
    237240
    238241// Set slice to cursor or all axes
    239     void setCursorSlice (casa::IPosition& start, casa::IPosition& end,
     242    void setCursorSlice(casa::IPosition& start, casa::IPosition& end,
    240243                         casa::Bool doAll, const SDMemTable& in) const;
    241244
    242245// Function to normalize data when averaging in time
    243246
    244   void normalize (casa::MaskedArray<casa::Float>& data,
     247  void normalize(casa::MaskedArray<casa::Float>& data,
    245248                  const casa::Array<casa::Float>& sumSq,
    246249                  const casa::Array<casa::Float>& tSysSumSq,
     
    250253
    251254// Put the data and mask into the SDContainer
    252    void putDataInSDC (SDContainer& sc, const casa::Array<casa::Float>& data,
     255   void putDataInSDC(SDContainer& sc, const casa::Array<casa::Float>& data,
    253256                      const casa::Array<casa::Bool>& mask) const;
    254257
    255258// Read ascii file into a Table
    256259
    257    casa::Table readAsciiFile (const casa::String& fileName) const;
     260   casa::Table readAsciiFile(const casa::String& fileName) const;
    258261};
    259262
Note: See TracChangeset for help on using the changeset viewer.