Ignore:
Timestamp:
06/09/10 19:03:06 (14 years ago)
Author:
Kana Sugimoto
Message:

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.


Location:
branches/alma
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/alma

  • branches/alma/src/STMath.cpp

    r1719 r1757  
    5252#include "RowAccumulator.h"
    5353#include "STAttr.h"
     54#include "STSelector.h"
     55
    5456#include "STMath.h"
    55 #include "STSelector.h"
    56 
    5757using namespace casa;
    5858
    … …  
    740740}
    741741
    742 CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left, 
    743                                             const CountedPtr<Scantable>& right, 
     742CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
     743                                            const CountedPtr<Scantable>& right,
    744744                                            const std::string& mode)
    745745{
    … …  
    13791379{
    13801380
    1381  
     1381
    13821382  STSelector sel;
    13831383  CountedPtr< Scantable > ws = getScantable(s, false);
    … …  
    15181518  // for each row
    15191519  // assume that the data order are same between sig and ref
    1520   RowAccumulator acc( asap::TINTSYS ) ;
     1520  RowAccumulator acc( asap::W_TINTSYS ) ;
    15211521  for ( int i = 0 ; i < sig->nrow() ; i++ ) {
    15221522    // get values
    … …  
    19521952  // initialize the lookup table if necessary
    19531953  if ( lookup.empty() ) {
    1954     lookup["NONE"]   = asap::NONE;
    1955     lookup["TINT"] = asap::TINT;
    1956     lookup["TINTSYS"]  = asap::TINTSYS;
    1957     lookup["TSYS"]  = asap::TSYS;
    1958     lookup["VAR"]  = asap::VAR;
     1954    lookup["NONE"]   = asap::W_NONE;
     1955    lookup["TINT"] = asap::W_TINT;
     1956    lookup["TINTSYS"]  = asap::W_TINTSYS;
     1957    lookup["TSYS"]  = asap::W_TSYS;
     1958    lookup["VAR"]  = asap::W_VAR;
    19591959  }
    19601960
    … …  
    20002000    String msg;
    20012001    if ( nc > 0 ) {
    2002       ppoly = new Polynomial<Float>(nc);
     2002      ppoly = new Polynomial<Float>(nc-1);
    20032003      coeff = coeffs;
    20042004      msg = String("user");
    … …  
    20062006      STAttr sdAttr;
    20072007      coeff = sdAttr.gainElevationPoly(inst);
    2008       ppoly = new Polynomial<Float>(3);
     2008      ppoly = new Polynomial<Float>(coeff.nelements()-1);
    20092009      msg = String("built in");
    20102010    }
    … …  
    21442144    if (d < 0) {
    21452145      Instrument inst =
    2146         STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), 
     2146        STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
    21472147                                  True);
    21482148      STAttr sda;
    … …  
    22072207
    22082208CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
    2209                                          float tau )
     2209                                         const std::vector<float>& tau )
    22102210{
    22112211  CountedPtr< Scantable > out = getScantable(in, false);
    22122212
    2213   Table tab = out->table();
    2214   ROScalarColumn<Float> elev(tab, "ELEVATION");
    2215   ArrayColumn<Float> specCol(tab, "SPECTRA");
    2216   ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
    2217   ArrayColumn<Float> tsysCol(tab, "TSYS");
    2218   for ( uInt i=0; i<tab.nrow(); ++i) {
    2219     Float zdist = Float(C::pi_2) - elev(i);
    2220     Float factor = exp(tau/cos(zdist));
    2221     MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
    2222     ma *= factor;
    2223     specCol.put(i, ma.getArray());
    2224     flagCol.put(i, flagsFromMA(ma));
    2225     Vector<Float> tsys;
    2226     tsysCol.get(i, tsys);
    2227     tsys *= factor;
    2228     tsysCol.put(i, tsys);
     2213  Table outtab = out->table();
     2214
     2215  const uInt ntau = uInt(tau.size());
     2216  std::vector<float>::const_iterator tauit = tau.begin();
     2217  AlwaysAssert((ntau == 1 || ntau == in->nif() || ntau == in->nif() * in->npol()),
     2218               AipsError);
     2219  TableIterator iiter(outtab, "IFNO");
     2220  while ( !iiter.pastEnd() ) {
     2221    Table itab = iiter.table();
     2222    TableIterator piter(outtab, "POLNO");
     2223    while ( !piter.pastEnd() ) {
     2224      Table tab = piter.table();
     2225      ROScalarColumn<Float> elev(tab, "ELEVATION");
     2226      ArrayColumn<Float> specCol(tab, "SPECTRA");
     2227      ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
     2228      ArrayColumn<Float> tsysCol(tab, "TSYS");
     2229      for ( uInt i=0; i<tab.nrow(); ++i) {
     2230        Float zdist = Float(C::pi_2) - elev(i);
     2231        Float factor = exp(*tauit/cos(zdist));
     2232        MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
     2233        ma *= factor;
     2234        specCol.put(i, ma.getArray());
     2235        flagCol.put(i, flagsFromMA(ma));
     2236        Vector<Float> tsys;
     2237        tsysCol.get(i, tsys);
     2238        tsys *= factor;
     2239        tsysCol.put(i, tsys);
     2240      }
     2241      if (ntau == in->nif()*in->npol() ) {
     2242        tauit++;
     2243      }
     2244      piter++;
     2245    }
     2246    if (ntau >= in->nif() ) {
     2247      tauit++;
     2248    }
     2249    iiter++;
    22292250  }
    22302251  return out;
    … …  
    22332254CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
    22342255                                             const std::string& kernel,
    2235                                              float width )
     2256                                             float width, int order)
    22362257{
    22372258  CountedPtr< Scantable > out = getScantable(in, false);
    … …  
    22542275      mathutil::runningMedian(specout, maskout, spec , mask, width);
    22552276      convertArray(flag, maskout);
     2277    } else if ( kernel == "poly" ) {
     2278      mathutil::polyfit(specout, maskout, spec, !mask, width, order);
     2279      convertArray(flag, !maskout);
    22562280    }
    22572281    flagCol.put(i, flag);
    … …  
    22622286
    22632287CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
    2264                                         const std::string& kernel, float width )
    2265 {
    2266   if (kernel == "rmedian"  || kernel == "hanning") {
    2267     return smoothOther(in, kernel, width);
     2288                                        const std::string& kernel, float width,
     2289                                        int order)
     2290{
     2291  if (kernel == "rmedian"  || kernel == "hanning" || kernel == "poly") {
     2292    return smoothOther(in, kernel, width, order);
    22682293  }
    22692294  CountedPtr< Scantable > out = getScantable(in, false);
    … …  
    23512376          id = out->molecules().addEntry(rf, name, fname);
    23522377          molidcol.put(k, id);
    2353           Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
    2354           (*it)->focus().getEntry(fax, ftan, frot, fhand,
     2378          Float fpa,frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
     2379          (*it)->focus().getEntry(fpa, fax, ftan, frot, fhand,
    23552380                                  fmount,fuser, fxy, fxyp,
    23562381                                  rec.asuInt("FOCUS_ID"));
    2357           id = out->focus().addEntry(fax, ftan, frot, fhand,
     2382          id = out->focus().addEntry(fpa, fax, ftan, frot, fhand,
    23582383                                     fmount,fuser, fxy, fxyp);
    23592384          focusidcol.put(k, id);
    … …  
    23992424  cols[3] = String("CYCLENO");
    24002425  TableIterator iter(tout, cols);
    2401   CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_, 
     2426  CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
    24022427                                               out->getPolType() );
    24032428  while (!iter.pastEnd()) {
    … …  
    24052430    ArrayColumn<Float> speccol(t, "SPECTRA");
    24062431    ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
    2407     ScalarColumn<Float> parancol(t, "PARANGLE");
    24082432    Matrix<Float> pols(speccol.getColumn());
    24092433    try {
    24102434      stpol->setSpectra(pols);
    2411       Float fang,fhand,parang;
    2412       fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
     2435      Float fang,fhand;
     2436      fang = in->focusTable_.getTotalAngle(focidcol(0));
    24132437      fhand = in->focusTable_.getFeedHand(focidcol(0));
    2414       parang = parancol(0);
    2415       /// @todo re-enable this
    2416       // disable total feed angle to support paralactifying Caswell style
    2417       stpol->setPhaseCorrections(parang, -parang, fhand);
     2438      stpol->setPhaseCorrections(fang, fhand);
    24182439      // use a member function pointer in STPol.  This only works on
    24192440      // the STPol pointer itself, not the Counted Pointer so
    … …  
    26812702      uInt row = tab.rowNumbers()[0];
    26822703      stpol->setSpectra(in->getPolMatrix(row));
    2683       Float fang,fhand,parang;
    2684       fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
     2704      Float fang,fhand;
     2705      fang = in->focusTable_.getTotalAngle(in->mfocusidCol_(row));
    26852706      fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
    2686       parang = in->paraCol_(row);
    2687       /// @todo re-enable this
    2688       // disable total feed angle to support paralactifying Caswell style
    2689       stpol->setPhaseCorrections(parang, -parang, fhand);
     2707      stpol->setPhaseCorrections(fang, fhand);
    26902708      Int npolout = 0;
    26912709      for (uInt i=0; i<tab.nrow(); ++i) {
    … …  
    27372755CountedPtr< Scantable >
    27382756  asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
    2739                           double frequency, double width )
     2757                         double start, double end,
     2758                         const std::string& mode)
    27402759{
    27412760  CountedPtr< Scantable > out = getScantable(in, false);
    … …  
    27552774      Vector<Float> spec = specCol(i);
    27562775      Vector<uChar> flag = flagCol(i);
    2757       Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
    2758       Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
     2776      int fstart = -1;
     2777      int fend = -1;
    27592778      for (unsigned int k=0; k < flag.nelements(); ++k ) {
    27602779        if (flag[k] > 0) {
    2761           spec[k] = 0.0;
     2780          fstart = k;
     2781          while (flag[k] > 0 && k < flag.nelements()) {
     2782            fend = k;
     2783            k++;
     2784          }
    27622785        }
     2786        Float interp = 0.0;
     2787        if (fstart-1 > 0 ) {
     2788          interp = spec[fstart-1];
     2789          if (fend+1 < spec.nelements()) {
     2790            interp = (interp+spec[fend+1])/2.0;
     2791          }
     2792        } else {
     2793          interp = spec[fend+1];
     2794        }
     2795        if (fstart > -1 && fend > -1) {
     2796          for (int j=fstart;j<=fend;++j) {
     2797            spec[j] = interp;
     2798          }
     2799        }
     2800        fstart =-1;
     2801        fend = -1;
    27632802      }
    27642803      Vector<Complex> lags;
    27652804      ffts.fft0(lags, spec);
    2766       Int start =  max(0, lag0);
    2767       Int end =  min(Int(lags.nelements()-1), lag1);
    2768       if (start == end) {
    2769         lags[start] = Complex(0.0);
     2805      Int lag0(start+0.5);
     2806      Int lag1(end+0.5);
     2807      if (mode == "frequency") {
     2808        lag0 = Int(spec.nelements()*abs(inc)/(start)+0.5);
     2809        lag1 = Int(spec.nelements()*abs(inc)/(end)+0.5);
     2810      }
     2811      Int lstart =  max(0, lag0);
     2812      Int lend =  min(Int(lags.nelements()-1), lag1);
     2813      if (lstart == lend) {
     2814        lags[lstart] = Complex(0.0);
    27702815      } else {
    2771         for (int j=start; j <=end ;++j) {
     2816        if (lstart > lend) {
     2817          Int tmp = lend;
     2818          lend = lstart;
     2819          lstart = tmp;
     2820        }
     2821        for (int j=lstart; j <=lend ;++j) {
    27722822          lags[j] = Complex(0.0);
    27732823        }
Note: See TracChangeset for help on using the changeset viewer.