Ignore:
Timestamp:
07/29/10 19:13:46 (14 years ago)
Author:
Kana Sugimoto
Message:

New Development: Yes

JIRA Issue: No (test merging alma branch)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description:


Location:
branches/mergetest
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/mergetest

  • branches/mergetest/external/atnf/PKSIO/PKSMS2reader.cc

    r1720 r1779  
    3434//#---------------------------------------------------------------------------
    3535
    36 #include <atnf/pks/pks_maths.h>
    37 #include <atnf/PKSIO/PKSmsg.h>
    38 #include <atnf/PKSIO/PKSrecord.h>
    39 #include <atnf/PKSIO/PKSMS2reader.h>
    40 
     36// AIPS++ includes.
    4137#include <casa/stdio.h>
    4238#include <casa/Arrays/ArrayMath.h>
    … …  
    4440#include <ms/MeasurementSets/MSColumns.h>
    4541#include <tables/Tables.h>
     42#include <casa/Quanta/MVTime.h>
     43#include <casa/Quanta/MVAngle.h>
     44#include <casa/BasicMath/Math.h>
     45#include <casa/Logging/LogIO.h>
     46#include <casa/Utilities/Sort.h>
     47#include <measures/Measures/MeasConvert.h>
     48#include <measures/Measures/MEpoch.h>
     49#include <measures/Measures/MeasRef.h>
     50
     51
     52// Parkes includes.
     53#include <atnf/pks/pks_maths.h>
     54#include <atnf/PKSIO/PKSrecord.h>
     55#include <atnf/PKSIO/PKSMS2reader.h>
     56
    4657
    4758//------------------------------------------------- PKSMS2reader::PKSMS2reader
    … …  
    5263{
    5364  cMSopen = False;
    54 
    55   // By default, messages are written to stderr.
    56   initMsg();
    5765}
    5866
    … …  
    7078Int PKSMS2reader::open(
    7179        const String msName,
     80        const String antenna,
    7281        Vector<Bool> &beams,
    7382        Vector<Bool> &IFs,
    … …  
    8897
    8998  cPKSMS  = MeasurementSet(msName);
     99
     100  // data selection by antenna
     101  if ( antenna.length() == 0 ) {
     102    cAntId.resize( 1 ) ;
     103    cAntId[0] = 0 ;
     104  }
     105  else {
     106    setupAntennaList( antenna ) ;
     107    if ( cAntId.size() > 1 ) {
     108      LogIO os( LogOrigin( "PKSMS2reader", "open()", WHERE ) ) ;
     109      os << LogIO::WARN << "PKSMS2reader is not ready for multiple antenna selection. Use first antenna id " << cAntId[0] << "."<< LogIO::POST ;
     110      Int tmp = cAntId[0] ;
     111      cAntId.resize( 1 ) ;
     112      cAntId[0] = tmp ;
     113    }
     114    stringstream ss ;
     115    ss << "SELECT FROM $1 WHERE ANTENNA1 == ANTENNA2 && ANTENNA1 IN [" ;
     116    for ( uInt i = 0 ; i < cAntId.size() ; i++ ) {
     117      ss << cAntId[i] ;
     118      if ( i == cAntId.size()-1 ) {
     119        ss << "]" ;
     120      }
     121      else {
     122        ss << "," ;
     123      }
     124    }
     125    string taql = ss.str() ;
     126    //cerr << "taql = " << taql << endl ;
     127    cPKSMS = MeasurementSet( tableCommand( taql, cPKSMS ) ) ;
     128  }
     129
     130  // taql access to the syscal table
     131  cHaveSysCal = False;
     132  if (cHaveSysCal=Table::isReadable(cPKSMS.sysCalTableName())) {
     133    cSysCalTab = Table(cPKSMS.sysCalTableName());
     134  }
     135
     136  // Lock the table for read access.
     137  cPKSMS.lock(False);
     138
    90139  cIdx    = 0;
     140  lastmjd = 0.0;
    91141  cNRow   = cPKSMS.nrow();
    92142  cMSopen = True;
    93 
    94   // Lock the table for read access.
    95   cPKSMS.lock(False);
    96143
    97144  // Main MS table and subtable column access.
    … …  
    107154  ROMSSysCalColumns       sysCalCols(cPKSMS.sysCal());
    108155  ROMSWeatherColumns      weatherCols(cPKSMS.weather());
     156  ROMSAntennaColumns      antennaCols(cPKSMS.antenna());
    109157
    110158  // Column accessors for required columns.
    … …  
    115163  cFieldIdCol.reference(msCols.fieldId());
    116164  cFieldNameCol.reference(fieldCols.name());
     165  cFieldDelayDirCol.reference(fieldCols.delayDir());
    117166
    118167  cSrcIdCol.reference(fieldCols.sourceId());
     168  cSrcId2Col.reference(sourceCols.sourceId());
    119169  cSrcNameCol.reference(sourceCols.name());
    120170  cSrcDirCol.reference(sourceCols.direction());
    … …  
    124174  cStateIdCol.reference(msCols.stateId());
    125175  cObsModeCol.reference(stateCols.obsMode());
     176  cCalCol.reference(stateCols.cal());
     177  cSigStateCol.reference(stateCols.sig());
     178  cRefStateCol.reference(stateCols.ref());
    126179
    127180  cDataDescIdCol.reference(msCols.dataDescId());
     181  cSpWinIdCol.reference(dataDescCols.spectralWindowId());
    128182  cChanFreqCol.reference(spWinCols.chanFreq());
     183  cTotBWCol.reference(spWinCols.totalBandwidth());
    129184
    130185  cWeatherTimeCol.reference(weatherCols.time());
    … …  
    135190  cBeamNoCol.reference(msCols.feed1());
    136191  cPointingCol.reference(pointingCols.direction());
     192  cPointingTimeCol.reference(pointingCols.time());
    137193  cSigmaCol.reference(msCols.sigma());
    138194  cNumReceptorCol.reference(feedCols.numReceptors());
    139195
    140196  // Optional columns.
     197  cHaveTsys = False;
     198  cHaveTcal = False;
    141199  if ((cHaveSrcVel = cPKSMS.source().tableDesc().isColumn("SYSVEL"))) {
    142200    cSrcVelCol.attach(cPKSMS.source(), "SYSVEL");
    143201  }
    144202
    145   if ((cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {
     203  if (cHaveSysCal && (cHaveTsys = cPKSMS.sysCal().tableDesc().isColumn("TSYS"))) {
    146204    cTsysCol.attach(cPKSMS.sysCal(), "TSYS");
     205  }
     206 
     207  if (cHaveSysCal && (cHaveTcal = cPKSMS.sysCal().tableDesc().isColumn("TCAL"))) {
     208    cTcalCol.attach(cPKSMS.sysCal(), "TCAL");
    147209  }
    148210
    … …  
    158220  // Spectral data should always be present.
    159221  haveSpectra = True;
    160   cFloatDataCol.reference(msCols.floatData());
     222  cHaveDataCol = False;
     223  cHaveCorrectedDataCol = False;
     224  ROMSObservationColumns observationCols(cPKSMS.observation());
     225  //String telName = observationCols.telescopeName()(0);
     226  cTelName = observationCols.telescopeName()(0);
     227  //cATF = cTelName.contains("ATF");
     228  //cOSF = cTelName.contains("OSF");
     229  //cALMA = cTelName.contains("ALMA");
     230  cALMA = cTelName.contains("ATF")||cTelName.contains("OSF")||
     231           cTelName.contains("ALMA");
     232
     233  if (cHaveDataCol = cPKSMS.isColumn(MSMainEnums::DATA)) {
     234    if (cALMA) {
     235      //try to read a single baseline interferometeric data
     236      //and treat it as single dish data
     237      //maybe extended for ALMA commissioning later
     238      cDataCol.reference(msCols.data());
     239      if (cHaveCorrectedDataCol = cPKSMS.isColumn(MSMainEnums::CORRECTED_DATA)) {
     240        //cerr<<"Do have CORRECTED_DATA column"<<endl;
     241        cCorrectedDataCol.reference(msCols.correctedData());
     242      }
     243    }
     244  }
     245  else {
     246    cFloatDataCol.reference(msCols.floatData());
     247  }
    161248  cFlagCol.reference(msCols.flag());
    162 
    163   if ((cGetXPol = cPKSMS.isColumn(MSMainEnums::DATA))) {
     249  cFlagRowCol.reference(msCols.flagRow());
     250
     251  if (cGetXPol = (cPKSMS.isColumn(MSMainEnums::DATA) && (!cALMA))) {
    164252    if ((cHaveXCalFctr = cPKSMS.tableDesc().isColumn("XCALFCTR"))) {
    165253      cXCalFctrCol.attach(cPKSMS, "XCALFCTR");
    … …  
    181269
    182270  // Number of IFs.
    183   uInt nIF = dataDescCols.nrow();
     271  //uInt nIF = dataDescCols.nrow();
     272  uInt nIF =spWinCols.nrow();
     273  Vector<Int> spWinIds = cSpWinIdCol.getColumn() ;
    184274  IFs.resize(nIF);
    185275  IFs = True;
     276  for ( Int ispw = 0 ; ispw < nIF ; ispw++ ) {
     277    if ( allNE( ispw, spWinIds ) ) {
     278      IFs(ispw) = False ;
     279    }
     280  }
    186281
    187282  // Number of polarizations and channels in each IF.
    188   ROScalarColumn<Int> spWinIdCol(dataDescCols.spectralWindowId());
    189283  ROScalarColumn<Int> numChanCol(spWinCols.numChan());
    190284
    … …  
    195289  nPol.resize(nIF);
    196290  for (uInt iIF = 0; iIF < nIF; iIF++) {
    197     nChan(iIF) = numChanCol(spWinIdCol(iIF));
    198     nPol(iIF)  = numPolCol(polIdCol(iIF));
     291    if ( IFs(iIF) ) {
     292      nChan(iIF) = numChanCol(cSpWinIdCol(iIF)) ;
     293      nPol(iIF) = numPolCol(polIdCol(iIF)) ;
     294    }
     295    else {
     296      nChan(iIF) = 0 ;
     297      nPol(iIF) = 0 ;
     298    }
    199299  }
    200300
    … …  
    203303  haveXPol = False;
    204304
    205   if (cGetXPol) {
     305  if (cGetXPol && !(cALMA)) {
    206306    for (Int irow = 0; irow < cNRow; irow++) {
    207307      if (cDataCol.isDefined(irow)) {
    … …  
    252352        String &antName,
    253353        Vector<Double> &antPosition,
     354        // before merge...
     355        //String &obsMode,
    254356        String &obsType,
    255357        String &bunit,
    … …  
    258360        Double &mjd,
    259361        Double &refFreq,
    260         Double &bandwidth)
     362        Double &bandwidth) 
    261363{
    262364  if (!cMSopen) {
    … …  
    271373  // Antenna name and ITRF coordinates.
    272374  ROMSAntennaColumns antennaCols(cPKSMS.antenna());
    273   antName = antennaCols.name()(0);
    274   antPosition = antennaCols.position()(0);
     375  //antName = antennaCols.name()(0);
     376  antName = antennaCols.name()(cAntId[0]);
     377  if (cALMA) {
     378     antName = cTelName + "-" + antName;
     379  }
     380  //antPosition = antennaCols.position()(0);
     381  antPosition = antennaCols.position()(cAntId[0]);
    275382
    276383  // Observation type.
    … …  
    282389  }
    283390
    284   // Brightness units.
    285   bunit = cPKSMS.unit(MSMainEnums::FLOAT_DATA);
    286 
     391  bunit = "";
     392  if (cHaveDataCol) {
     393    const TableRecord& keywordSet2
     394       = cDataCol.columnDesc().keywordSet();
     395    if(keywordSet2.isDefined("UNIT")) {
     396      bunit = keywordSet2.asString("UNIT");
     397    }
     398  } else {
     399    const TableRecord& keywordSet
     400       = cFloatDataCol.columnDesc().keywordSet();
     401    if(keywordSet.isDefined("UNIT")) {
     402      bunit = keywordSet.asString("UNIT");
     403    }
     404  }
     405
     406/***
     407  const TableRecord& keywordSet
     408       = cFloatDataCol.columnDesc().keywordSet();
     409  if(keywordSet.isDefined("UNIT")) {
     410    fluxunit = keywordSet.asString("UNIT");
     411  }
     412***/
    287413  // Coordinate equinox.
    288414  ROMSPointingColumns pointingCols(cPKSMS.pointing());
    289415  String dirref = pointingCols.direction().keywordSet().asRecord("MEASINFO").
    290416                    asString("Ref");
     417  cDirRef = dirref;
     418  if (dirref =="AZELGEO" || dirref == "AZEL") {
     419     dirref = "J2000";
     420  }
    291421  sscanf(dirref.chars()+1, "%f", &equinox);
    292422
    … …  
    357487        const Bool getSpectra,
    358488        const Bool getXPol,
     489        const Bool getFeedPos,
     490        const Bool getPointing,
    359491        const Int  coordSys)
    360492{
    … …  
    436568  cGetXPol = cGetXPol && getXPol;
    437569
     570  // Get feed positions?  (Not available.)
     571  cGetFeedPos = False;
     572
     573  // Get Pointing data (for MS)
     574  cGetPointing = getPointing;
     575
    438576  // Coordinate system?  (Only equatorial available.)
    439577  cCoordSys = 0;
    … …  
    490628// Read the next data record.
    491629
     630/**
     631Int PKSMS2reader::read(
     632        Int             &scanNo,
     633        Int             &cycleNo,
     634        Double          &mjd,
     635        Double          &interval,
     636        String          &fieldName,
     637        String          &srcName,
     638        Vector<Double>  &srcDir,
     639        Vector<Double>  &srcPM,
     640        Double          &srcVel,
     641        String          &obsMode,
     642        Int             &IFno,
     643        Double          &refFreq,
     644        Double          &bandwidth,
     645        Double          &freqInc,
     646        Vector<Double>  &restFreq,
     647        Vector<Float>   &tcal,
     648        String          &tcalTime,
     649        Float           &azimuth,
     650        Float           &elevation,
     651        Float           &parAngle,
     652        Float           &focusAxi,
     653        Float           &focusTan,
     654        Float           &focusRot,
     655        Float           &temperature,
     656        Float           &pressure,
     657        Float           &humidity,
     658        Float           &windSpeed,
     659        Float           &windAz,
     660        Int             &refBeam,
     661        Int             &beamNo,
     662        Vector<Double>  &direction,
     663        Vector<Double>  &scanRate,
     664        Vector<Float>   &tsys,
     665        Vector<Float>   &sigma,
     666        Vector<Float>   &calFctr,
     667        Matrix<Float>   &baseLin,
     668        Matrix<Float>   &baseSub,
     669        Matrix<Float>   &spectra,
     670        Matrix<uChar>   &flagged,
     671        uInt            &flagrow,
     672        Complex         &xCalFctr,
     673        Vector<Complex> &xPol)
     674**/
    492675Int PKSMS2reader::read(PKSrecord &pksrec)
    493676{
     677  LogIO os( LogOrigin( "PKSMS2reader", "read()", WHERE ) ) ;
     678
    494679  if (!cMSopen) {
    495680    return 1;
    … …  
    504689  Int ibeam;
    505690  Int iIF;
     691  Int iDataDesc;
     692
    506693  while (True) {
    507694    ibeam = cBeamNoCol(cIdx);
    508     iIF   = cDataDescIdCol(cIdx);
     695    iDataDesc   = cDataDescIdCol(cIdx);
     696    iIF   =cSpWinIdCol(iDataDesc);
    509697    if (cBeams(ibeam) && cIFs(iIF)) {
    510698      break;
    … …  
    516704    }
    517705  }
    518 
    519706  // Renumerate scan no. Here still is 1-based
    520   pksrec.scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
     707  //scanNo = cScanNoCol(cIdx) - cScanNoCol(0) + 1;
     708  //scanNo = cScanNoCol(cIdx);
     709  pksrec.scanNo = cScanNoCol(cIdx);
    521710
    522711  if (pksrec.scanNo != cScanNo) {
    … …  
    542731
    543732  Int srcId = cSrcIdCol(fieldId);
    544   pksrec.srcName = cSrcNameCol(srcId);
     733  //For source with multiple spectral window setting, this is not
     734  // correct. Source name of srcId may not be at 'srcId'th row of SrcNameCol
     735  //srcName = cSrcNameCol(srcId);
     736  for (uInt irow = 0; irow < cSrcId2Col.nrow(); irow++) {
     737    if (cSrcId2Col(irow) == srcId) {
     738      //srcName = cSrcNameCol(irow);
     739      pksrec.srcName = cSrcNameCol(irow);
     740    }
     741  }
     742
    545743  pksrec.srcDir  = cSrcDirCol(srcId);
    546744  pksrec.srcPM   = cSrcPMCol(srcId);
    547745
    548746  // Systemic velocity.
    549   if (!cHaveSrcVel) {
     747  if (!cHaveSrcVel || cALMA) {
    550748    pksrec.srcVel = 0.0f;
    551749  } else {
    … …  
    553751  }
    554752
     753  ROMSAntennaColumns antennaCols(cPKSMS.antenna());
     754  //String telescope = antennaCols.name()(0);
     755  String telescope = antennaCols.name()(cAntId[0]);
     756  Bool cGBT = telescope.contains("GBT");
     757  //Bool cPM = telescope.contains("PM"); // ACA TP antenna
     758  //Bool cDV = telescope.contains("DV"); // VERTEX
     759  //Bool cCM = telescope.contains("CM"); // ACA 7m antenna
     760  //Bool cALMA = cPM || cDV || cCM ;
    555761  // Observation type.
    556   Int stateId = cStateIdCol(cIdx);
    557   pksrec.obsType = cObsModeCol(stateId);
     762  // check if State Table exist
     763  //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
     764  Int stateId = 0;
     765  Int StateNRow = 0;
     766  StateNRow=cObsModeCol.nrow();
     767  if (Table::isReadable(cPKSMS.stateTableName())) {
     768        pksrec.obsType = " ";
     769    if (StateNRow > 0) {
     770      stateId = cStateIdCol(cIdx);
     771      if (stateId == -1) {
     772        //pksrec.obsType = " ";
     773      } else {
     774        pksrec.obsType = cObsModeCol(stateId);
     775        Bool sigState =cSigStateCol(stateId);
     776        Bool refState =cRefStateCol(stateId);
     777        //DEBUG
     778        //cerr <<"stateid="<<stateId<<" obsmode="<<pksrec.obsType<<endl;
     779        if (cGBT) {
     780          // split the obsType string and append a proper label
     781          // (these are GBT specific)
     782          int epos = pksrec.obsType.find_first_of(':');
     783          int nextpos = pksrec.obsType.find_first_of(':',epos+1);
     784          string obsMode1 = pksrec.obsType.substr(0,epos);
     785          string obsMode2 = pksrec.obsType.substr(epos+1,nextpos-epos-1);
     786     
     787          //cerr <<"obsMode2= "<<obsMode2<<endl;
     788          if (!pksrec.srcName.contains("_ps")
     789              &&!pksrec.srcName.contains("_psr")
     790              &&!pksrec.srcName.contains("_nod")
     791              &&!pksrec.srcName.contains("_fs")
     792              &&!pksrec.srcName.contains("_fsr")) {
     793            // if Nod mode observation , append '_nod'
     794            if (obsMode1 == "Nod") {
     795              //pksrec.srcName.append("_nod");
     796              pksrec.srcType = SrcType::NOD ;
     797            } else if (obsMode1 == "OffOn") {
     798            // for GBT position switch observations (OffOn or OnOff)
     799              //if (obsMode2 == "PSWITCHON") pksrec.srcName.append("_ps");
     800              //if (obsMode2 == "PSWITCHOFF") pksrec.srcName.append("_psr");
     801              if (obsMode2 == "PSWITCHON") pksrec.srcType = SrcType::PSON ;
     802              if (obsMode2 == "PSWITCHOFF") pksrec.srcType = SrcType::PSOFF ;
     803            } else {
     804              if (obsMode2 == "FSWITCH") {
     805              // for GBT frequency switch mode
     806                //if (sigState) pksrec.srcName.append("_fs");
     807                //if (refState) pksrec.srcName.append("_fsr");
     808                if (sigState) pksrec.srcType = SrcType::FSON ;
     809                if (refState) pksrec.srcType = SrcType::FSOFF ;
     810              }
     811            }
     812          }
     813        }
     814        else if (cALMA) {
     815          // ALMA tag
     816          // split the obsType string and append a proper label
     817          string substr[1] ;
     818          int numSubstr = split( pksrec.obsType, substr, 1, "," );
     819          String obsType = String( substr[0] );
     820          int epos = obsType.find_first_of('.');
     821          int nextpos = obsType.find_first_of('.',epos+1);
     822          string obsMode1 = obsType.substr(0,epos);
     823          string obsMode2 = obsType.substr(epos+1,nextpos-epos-1);
     824     
     825          //cerr <<"obsMode2= "<<obsMode2<<endl;
     826          // Current OBS_MODE format:
     827          //
     828          //     ON: OBSERVE_TARGET.ON_SOURCE
     829          //    OFF: OBSERVE_TARGET.OFF_SOURCE
     830          //
     831          if (obsMode1 == "OBSERVE_TARGET") {
     832            //if (obsMode2 == "ON_SOURCE") pksrec.srcName.append("_pson");
     833            //if (obsMode2 == "OFF_SOURCE") pksrec.srcName.append("_psoff");
     834            if (obsMode2 == "ON_SOURCE") pksrec.srcType = SrcType::PSON ;
     835            if (obsMode2 == "OFF_SOURCE") pksrec.srcType = SrcType::PSOFF ;
     836          }
     837        }
     838      }
     839    }
     840  }
     841  // CAL state
     842  // this should be apply just for GBT data?
     843  Double Cal;
     844  if (stateId==-1 || StateNRow==0) {
     845    Cal = 0;
     846  } else {
     847    Cal = cCalCol(stateId);
     848  }
     849  if (cGBT) {
     850    if (Cal > 0 && !pksrec.srcName.contains("_calon")) {
     851      //pksrec.srcName.append("_calon");
     852      if ( pksrec.srcType == SrcType::NOD )
     853        pksrec.srcType = SrcType::NODCAL ;
     854      else if ( pksrec.srcType == SrcType::PSON )
     855        pksrec.srcType = SrcType::PONCAL ;
     856      else if ( pksrec.srcType == SrcType::PSOFF )
     857        pksrec.srcType = SrcType::POFFCAL ;
     858      else if ( pksrec.srcType == SrcType::FSON )
     859        pksrec.srcType = SrcType::FONCAL ;
     860      else if ( pksrec.srcType == SrcType::FSOFF )
     861        pksrec.srcType = SrcType::FOFFCAL ;
     862      else
     863        pksrec.srcName.append("_calon");
     864    }
     865  }
    558866
    559867  pksrec.IFno = iIF + 1;
    560868  Int nChan = abs(cEndChan(iIF) - cStartChan(iIF)) + 1;
    561 
     869 
    562870  // Minimal handling on continuum data.
    563871  Vector<Double> chanFreq = cChanFreqCol(iIF);
     872  pksrec.nchan = nChan;
    564873  if (nChan == 1) {
    565     cout << "The input is continuum data. "<< endl;
    566     pksrec.freqInc  = chanFreq(0);
     874    //pksrec.freqInc  = chanFreq(0);
     875    pksrec.freqInc  = cTotBWCol(iIF);
    567876    pksrec.refFreq  = chanFreq(0);
    568     pksrec.restFreq = 0.0f;
     877    pksrec.restFreq.resize(1);
     878    pksrec.restFreq[0] = 0.0f;
    569879  } else {
     880 
    570881    if (cStartChan(iIF) <= cEndChan(iIF)) {
    571882      pksrec.freqInc = chanFreq(1) - chanFreq(0);
    … …  
    575886
    576887    pksrec.refFreq  = chanFreq(cRefChan(iIF)-1);
    577     pksrec.restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
    578   }
    579   pksrec.bandwidth = abs(pksrec.freqInc * nChan);
     888
     889    Bool HaveSrcRestFreq= cSrcRestFrqCol.isDefined(srcId);
     890    if (HaveSrcRestFreq) {
     891      //restFreq = cSrcRestFrqCol(srcId)(IPosition(1,0));
     892      //restFreq = cSrcRestFrqCol(srcId);
     893      pksrec.restFreq = cSrcRestFrqCol(srcId);
     894    } else {
     895      pksrec.restFreq.resize(1);
     896      pksrec.restFreq[0] = 0.0f;
     897    }
     898  }
     899  //pksrec.bandwidth = abs(pksrec.freqInc * nChan);
     900  pksrec.bandwidth = abs(cTotBWCol(0));
    580901
    581902  pksrec.tcal.resize(cNPol(iIF));
    582903  pksrec.tcal      = 0.0f;
    583904  pksrec.tcalTime  = "";
    584   pksrec.azimuth   = 0.0f;
    585   pksrec.elevation = 0.0f;
     905//  pksrec.azimuth   = 0.0f;
     906//  pksrec.elevation = 0.0f;
    586907  pksrec.parAngle  = 0.0f;
    587908
    … …  
    591912
    592913  // Find the appropriate entry in the WEATHER subtable.
    593   Vector<Double> wTimes = cWeatherTimeCol.getColumn();
    594   Int weatherIdx;
    595   for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
    596     if (cWeatherTimeCol(weatherIdx) <= time) {
    597       break;
    598     }
    599   }
    600 
    601   if (weatherIdx < 0) {
     914  //Bool cHaveStateTab=Table::isReadable(cPKSMS.stateTableName());
     915  Bool cHaveWeatherTab = Table::isReadable(cPKSMS.weatherTableName());
     916  Int weatherIdx=-1;
     917  if (cHaveWeatherTab) {
     918    Vector<Double> wTimes = cWeatherTimeCol.getColumn();
     919    for (weatherIdx = wTimes.nelements()-1; weatherIdx >= 0; weatherIdx--) {
     920      if (cWeatherTimeCol(weatherIdx) <= time) {
     921        break;
     922      }
     923    }
     924  }
     925
     926  if (weatherIdx < 0 || !cHaveWeatherTab) {
    602927    // No appropriate WEATHER entry.
    603928    pksrec.temperature = 0.0f;
    … …  
    616941  pksrec.beamNo  = ibeam + 1;
    617942
    618   Matrix<Double> pointingDir = cPointingCol(fieldId);
    619   pksrec.direction = pointingDir.column(0);
     943  //pointing/azel
     944  MVPosition mvpos(antennaCols.position()(cAntId[0]));
     945  MPosition mp(mvpos);
     946  Quantum<Double> qt(time,"s");
     947  MVEpoch mvt(qt);
     948  MEpoch me(mvt);
     949  MeasFrame frame(mp, me);
     950  MDirection md;
    620951  pksrec.pCode = 0;
    621952  pksrec.rateAge = 0.0f;
    622   uInt ncols = pointingDir.ncolumn();
    623   if (ncols == 1) {
    624     pksrec.scanRate = 0.0f;
    625   } else {
    626     pksrec.scanRate(0) = pointingDir.column(1)(0);
    627     pksrec.scanRate(1) = pointingDir.column(1)(1);
    628   }
    629953  pksrec.paRate = 0.0f;
     954  if (cGetPointing) {
     955    //cerr << "get pointing data ...." << endl;
     956    ROScalarColumn<Int> pAntIdCol ;
     957    ROScalarColumn<Double> psTimeCol ;
     958    Table ptTable = cPKSMS.pointing() ;
     959    MSPointing selPtTab( ptTable( ptTable.col("ANTENNA_ID") == cAntId[0] ) ) ;
     960    pAntIdCol.attach( selPtTab, "ANTENNA_ID" ) ;
     961    Vector<Int> antIds = pAntIdCol.getColumn() ;
     962    psTimeCol.attach( selPtTab, "TIME" ) ;
     963    Vector<Double> pTimes = psTimeCol.getColumn();
     964    Bool doInterp = False ;
     965    Int PtIdx=-1;
     966    for (PtIdx = pTimes.nelements()-1; PtIdx >= 0; PtIdx--) {
     967      if ( pTimes[PtIdx] == time ) {
     968        break ;
     969      }
     970      else if ( pTimes[PtIdx] < time ) {
     971        if ( PtIdx != pTimes.nelements()-1 ) {
     972          doInterp = True ;
     973        }
     974        break ;
     975      }
     976    }
     977    if ( PtIdx == -1 ) {
     978      PtIdx = 0 ;
     979    }
     980    //cerr << "got index=" << PtIdx << endl;
     981    Matrix<Double> pointingDir = cPointingCol(PtIdx);
     982    ROMSPointingColumns PtCols( selPtTab ) ;
     983    Vector<Double> pointingDirVec ;
     984    if ( doInterp ) {
     985      Double dt1 = time - pTimes[PtIdx] ;
     986      Double dt2 = pTimes[PtIdx+1] - time ;
     987      Vector<Double> dirVec1 = pointingDir.column(0) ;
     988      Matrix<Double> pointingDir2 = cPointingCol(PtIdx+1) ;
     989      Vector<Double> dirVec2 = pointingDir2.column(0) ;
     990      pointingDirVec = (dt1*dirVec2+dt2*dirVec1)/(dt1+dt2) ;
     991      Vector<MDirection> vmd1(1) ;
     992      Vector<MDirection> vmd2(1) ;
     993      PtCols.directionMeasCol().get(PtIdx,vmd1) ;
     994      Vector<Double> angle1 = vmd1(0).getAngle().getValue("rad") ;
     995      PtCols.directionMeasCol().get(PtIdx+1,vmd2) ;
     996      Vector<Double> angle2 = vmd2(0).getAngle().getValue("rad") ;
     997      Vector<Double> angle = (dt1*angle2+dt2*angle1)/(dt1+dt2) ;
     998      Quantum< Vector<Double> > qangle( angle, "rad" ) ;
     999      String typeStr = vmd1(0).getRefString() ;
     1000      //cerr << "vmd1.getRefString()=" << typeStr << endl ;
     1001      MDirection::Types mdType ;
     1002      MDirection::getType( mdType, typeStr ) ;
     1003      //cerr << "mdType=" << mdType << endl ;
     1004      md = MDirection( qangle, mdType ) ;
     1005      //cerr << "md=" << md.getAngle().getValue("rad") << endl ;
     1006    }
     1007    else {
     1008      pointingDirVec = pointingDir.column(0) ;
     1009      Vector<MDirection> vmd(1);
     1010      PtCols.directionMeasCol().get(PtIdx,vmd);
     1011      md = vmd[0];
     1012    }
     1013    // put J2000 coordinates in "direction"
     1014    if (cDirRef =="J2000") {
     1015      pksrec.direction = pointingDirVec ;
     1016    }
     1017    else {
     1018      pksrec.direction =
     1019        MDirection::Convert(md, MDirection::Ref(MDirection::J2000,
     1020                                                frame)
     1021                            )().getAngle("rad").getValue();
     1022     
     1023    }
     1024    uInt ncols = pointingDir.ncolumn();
     1025    pksrec.scanRate.resize(2);
     1026    if (ncols == 1) {
     1027      pksrec.scanRate = 0.0f;
     1028    } else {
     1029      pksrec.scanRate(0) = pointingDir.column(1)(0);
     1030      pksrec.scanRate(1) = pointingDir.column(1)(1);
     1031    }
     1032  }
     1033  else {
     1034  // Get direction from FIELD table
     1035  // here, assume direction to be the field direction not pointing
     1036    Matrix<Double> delayDir = cFieldDelayDirCol(fieldId);
     1037    pksrec.direction = delayDir.column(0);
     1038    uInt ncols = delayDir.ncolumn();
     1039    pksrec.scanRate.resize(2);
     1040    if (ncols == 1) {
     1041      pksrec.scanRate = 0.0f;
     1042    } else {
     1043      pksrec.scanRate(0)  = delayDir.column(1)(0);
     1044      pksrec.scanRate(1)  = delayDir.column(1)(1);
     1045    }
     1046  }
     1047  // caluculate azimuth and elevation
     1048  // first, get the reference frame
     1049 /**
     1050  MVPosition mvpos(antennaCols.position()(0));
     1051  MPosition mp(mvpos);
     1052  Quantum<Double> qt(time,"s");
     1053  MVEpoch mvt(qt);
     1054  MEpoch me(mvt);
     1055  MeasFrame frame(mp, me);
     1056  **/
     1057  //
     1058  ROMSFieldColumns fldCols(cPKSMS.field());
     1059  Vector<MDirection> vmd(1);
     1060  //MDirection md;
     1061  fldCols.delayDirMeasCol().get(fieldId,vmd);
     1062  md = vmd[0];
     1063  //Vector<Double> dircheck = md.getAngle("rad").getValue();
     1064  //cerr<<"dircheck="<<dircheck<<endl;
     1065
     1066  Vector<Double> azel =
     1067        MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
     1068                                                frame)
     1069                            )().getAngle("rad").getValue();
     1070  //cerr<<"azel="<<azel<<endl;
     1071  pksrec.azimuth = azel[0];
     1072  pksrec.elevation = azel[1];
    6301073
    6311074  // Get Tsys assuming that entries in the SYSCAL table match the main table.
    … …  
    6461089  cSigmaCol.get(cIdx, pksrec.sigma, True);
    6471090
     1091  //get Tcal if available
     1092  if (cHaveTcal) {
     1093    Int nTcalColRow = cTcalCol.nrow();
     1094    uInt nBeam = cBeams.nelements();
     1095    uInt nIF = cIFs.nelements();
     1096    uInt nrws = nBeam * nIF;
     1097    if (nTcalColRow > 0) { 
     1098    // find tcal match with the data with the data time stamp
     1099      Double mjds = pksrec.mjd*(24*3600);
     1100      Double dtcalTime;
     1101      if ( pksrec.mjd > lastmjd || cIdx==0 ) {
     1102        //Table tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds));
     1103        tmptab = cSysCalTab(near(cSysCalTab.col("TIME"),mjds), nrws);
     1104        //DEBUG
     1105        //if (cIdx == 0) {
     1106        //  cerr<<"inital table retrieved"<<endl;
     1107        //}
     1108       
     1109      }
     1110
     1111      if (nBeam == 1) {
     1112        tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF, 1);
     1113      } else {
     1114        tmptab2 = tmptab( tmptab.col("SPECTRAL_WINDOW_ID") == iIF &&
     1115                              tmptab.col("FEED_ID") == ibeam , 1);
     1116      }
     1117      //cerr<<"first subtab rows="<<tmptab.nrow()<<endl;
     1118      int syscalrow = tmptab2.nrow();
     1119      ROArrayColumn<Float> tcalCol(tmptab2, "TCAL");
     1120      ROScalarColumn<Double> tcalTimeCol(tmptab2, "TIME");
     1121      if (syscalrow==0) {
     1122        os << LogIO::NORMAL
     1123           <<"Cannot find any matching Tcal at/near the data timestamp."
     1124           << " Set Tcal=0.0" << LogIO::POST ;
     1125      } else {
     1126        tcalCol.get(0, pksrec.tcal);
     1127        tcalTimeCol.get(0,dtcalTime);
     1128        pksrec.tcalTime = MVTime(dtcalTime/(24*3600)).string(MVTime::YMD);
     1129        //DEBUG
     1130        //cerr<<"cIdx:"<<cIdx<<" tcal="<<tcal<<" tcalTime="<<tcalTime<<endl;
     1131        tmptab.markForDelete();
     1132        tmptab2.markForDelete();
     1133      }
     1134    }
     1135    lastmjd = pksrec.mjd;
     1136  }
     1137
    6481138  // Calibration factors (if available).
    6491139  pksrec.calFctr.resize(cNPol(iIF));
    … …  
    6721162    Matrix<Float> tmpData;
    6731163    Matrix<Bool>  tmpFlag;
    674     cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
     1164    if (cHaveDataCol) {
     1165      Matrix<Complex> tmpCmplxData;
     1166      Matrix<Float> tmpReData;
     1167      Matrix<Float> tmpImData;
     1168      //cerr<<"reading spectra..."<<endl;
     1169      //# TODO - should have a flag to user to select DATA or CORRECTED_DATA
     1170      //# currently just automatically determined, --- read CORRECTED one
     1171      //# if the column exist.
     1172      if (cHaveCorrectedDataCol) {
     1173        cCorrectedDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
     1174      } else {
     1175        cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
     1176      }
     1177      tmpReData = real(tmpCmplxData);
     1178      tmpImData = imag(tmpCmplxData);
     1179      tmpData = sqrt(tmpReData*tmpReData + tmpImData*tmpImData);
     1180    } else {
     1181      cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
     1182    }
    6751183    cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
    6761184
    … …  
    6981206      }
    6991207    }
     1208
     1209    // Row-based flagging info. (True:1, False:0)
     1210    pksrec.flagrow = (cFlagRowCol(cIdx) ? 1 : 0);
    7001211  }
    7011212
    7021213  // Get cross-polarization data.
    7031214  if (cGetXPol) {
     1215    //cerr<<"cGetXPol="<<cGetXPol<<endl;
     1216    //cerr<<"cHaveXCalFctr="<<cHaveXCalFctr<<endl;
     1217
    7041218    if (cHaveXCalFctr) {
    7051219      cXCalFctrCol.get(cIdx, pksrec.xCalFctr);
    … …  
    7081222    }
    7091223
    710     cDataCol.get(cIdx, pksrec.xPol, True);
    711 
    712     if (cEndChan(iIF) < cStartChan(iIF)) {
    713       Complex ctmp;
     1224    if(!cALMA) {
     1225      cDataCol.get(cIdx, pksrec.xPol, True);
     1226
     1227      if (cEndChan(iIF) < cStartChan(iIF)) {
     1228        Complex ctmp;
     1229        Int jchan = nChan - 1;
     1230        for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
     1231          ctmp = pksrec.xPol(ichan);
     1232          pksrec.xPol(ichan) = pksrec.xPol(jchan);
     1233          pksrec.xPol(jchan) = ctmp;
     1234        }
     1235      }
     1236    }
     1237  }
     1238  /**
     1239  cerr<<"scanNo="<<scanNo<<endl;
     1240  cerr<<"cycleNo="<<cycleNo<<endl;
     1241  cerr<<"mjd="<<mjd<<endl;
     1242  cerr<<"interval="<<interval<<endl;
     1243  cerr<<"fieldName="<<fieldName<<endl;
     1244  cerr<<"srcNmae="<<srcName<<endl;
     1245  cerr<<"srcDir="<<srcDir<<endl;
     1246  cerr<<"srcPM="<<srcPM<<endl;
     1247  cerr<<"srcVel="<<srcVel<<endl;
     1248  cerr<<"obsMode="<<obsMode<<endl;
     1249  cerr<<"IFno="<<IFno<<endl;
     1250  cerr<<"refFreq="<<refFreq<<endl;
     1251  cerr<<"tcal="<<tcal<<endl;
     1252  cerr<<"direction="<<direction<<endl;
     1253  cerr<<"scanRate="<<scanRate<<endl;
     1254  cerr<<"tsys="<<tsys<<endl;
     1255  cerr<<"sigma="<<sigma<<endl;
     1256  cerr<<"calFctr="<<calFctr<<endl;
     1257  cerr<<"baseLin="<<baseLin<<endl;
     1258  cerr<<"baseSub="<<baseSub<<endl;
     1259  cerr<<"spectra="<<spectra.shape()<<endl;
     1260  cerr<<"flagged="<<flagged.shape()<<endl;
     1261  cerr<<"xCalFctr="<<xCalFctr<<endl;
     1262  cerr<<"xPol="<<xPol<<endl;
     1263  **/
     1264  cIdx++;
     1265
     1266  return 0;
     1267}
     1268
     1269//--------------------------------------------------------- PKSMS2reader::read
     1270
     1271// Read the next data record, just the basics.
     1272
     1273Int PKSMS2reader::read(
     1274        Int           &IFno,
     1275        Vector<Float> &tsys,
     1276        Vector<Float> &calFctr,
     1277        Matrix<Float> &baseLin,
     1278        Matrix<Float> &baseSub,
     1279        Matrix<Float> &spectra,
     1280        Matrix<uChar> &flagged)
     1281{
     1282  if (!cMSopen) {
     1283    return 1;
     1284  }
     1285
     1286  // Check for EOF.
     1287  if (cIdx >= cNRow) {
     1288    return -1;
     1289  }
     1290
     1291  // Find the next selected beam and IF.
     1292  Int ibeam;
     1293  Int iIF;
     1294  Int iDataDesc;
     1295  while (True) {
     1296    ibeam = cBeamNoCol(cIdx);
     1297    //iIF   = cDataDescIdCol(cIdx);
     1298    iDataDesc   = cDataDescIdCol(cIdx);
     1299    iIF   = cSpWinIdCol(iDataDesc);
     1300    if (cBeams(ibeam) && cIFs(iIF)) {
     1301      break;
     1302    }
     1303
     1304    // Check for EOF.
     1305    if (++cIdx >= cNRow) {
     1306      return -1;
     1307    }
     1308  }
     1309
     1310  IFno = iIF + 1;
     1311  // Get Tsys assuming that entries in the SYSCAL table match the main table.
     1312  cTsysCol.get(cIdx, tsys, True);
     1313
     1314  // Calibration factors (if available).
     1315  if (cHaveCalFctr) {
     1316    cCalFctrCol.get(cIdx, calFctr, True);
     1317  } else {
     1318    calFctr.resize(cNPol(iIF));
     1319    calFctr = 0.0f;
     1320  }
     1321
     1322  // Baseline parameters (if available).
     1323  if (cHaveBaseLin) {
     1324    baseLin.resize(2,cNPol(iIF));
     1325    cBaseLinCol.get(cIdx, baseLin);
     1326
     1327    baseSub.resize(24,cNPol(iIF));
     1328    cBaseSubCol.get(cIdx, baseSub);
     1329
     1330  } else {
     1331    baseLin.resize(0,0);
     1332    baseSub.resize(0,0);
     1333  }
     1334
     1335  if (cGetSpectra) {
     1336    // Get spectral data.
     1337    Matrix<Float> tmpData;
     1338    Matrix<Bool>  tmpFlag;
     1339    if (cHaveDataCol) {
     1340      Matrix<Complex> tmpCmplxData;
     1341      cDataCol.getSlice(cIdx, cDataSel(iIF), tmpCmplxData, True);
     1342      tmpData = real(tmpCmplxData);
     1343    } else {
     1344      cFloatDataCol.getSlice(cIdx, cDataSel(iIF), tmpData, True);
     1345    }
     1346    cFlagCol.getSlice(cIdx, cDataSel(iIF), tmpFlag, True);
     1347
     1348    // Transpose spectra.
     1349    Int nChan = tmpData.ncolumn();
     1350    Int nPol  = tmpData.nrow();
     1351    spectra.resize(nChan, nPol);
     1352    flagged.resize(nChan, nPol);
     1353    if (cEndChan(iIF) >= cStartChan(iIF)) {
     1354      // Simple transposition.
     1355      for (Int ipol = 0; ipol < nPol; ipol++) {
     1356        for (Int ichan = 0; ichan < nChan; ichan++) {
     1357          spectra(ichan,ipol) = tmpData(ipol,ichan);
     1358          flagged(ichan,ipol) = tmpFlag(ipol,ichan);
     1359        }
     1360      }
     1361
     1362    } else {
     1363      // Transpose with inversion.
    7141364      Int jchan = nChan - 1;
    715       for (Int ichan = 0; ichan < nChan/2; ichan++, jchan--) {
    716         ctmp = pksrec.xPol(ichan);
    717         pksrec.xPol(ichan) = pksrec.xPol(jchan);
    718         pksrec.xPol(jchan) = ctmp;
     1365      for (Int ipol = 0; ipol < nPol; ipol++) {
     1366        for (Int ichan = 0; ichan < nChan; ichan++, jchan--) {
     1367          spectra(ichan,ipol) = tmpData(ipol,jchan);
     1368          flagged(ichan,ipol) = tmpFlag(ipol,jchan);
     1369        }
    7191370      }
    7201371    }
    … …  
    7351386  cMSopen = False;
    7361387}
     1388
     1389//-------------------------------------------------------- PKSMS2reader::splitAntenanSelectionString
     1390
     1391// split antenna selection string
     1392// delimiter is ','
     1393
     1394Vector<String> PKSMS2reader::splitAntennaSelectionString( const String s )
     1395{
     1396  Char delim = ',' ;
     1397  Int n = s.freq( delim ) + 1 ;
     1398  Vector<String> antlist ;
     1399  string sl[n] ;
     1400  Int numSubstr = split( s, sl, n, "," );
     1401  antlist.resize( numSubstr ) ;
     1402  for ( Int i = 0 ; i < numSubstr ; i++ ) {
     1403    antlist[i] = String( sl[i] ) ;
     1404    antlist[i].trim() ;
     1405  }
     1406  //cerr << "antlist = " << antlist << endl ;
     1407  return antlist ;
     1408}
     1409
     1410//-------------------------------------------------------- PKSMS2reader::setupAntennaList
     1411
     1412// Fill cAntenna and cAntId
     1413
     1414void PKSMS2reader::setupAntennaList( const String s )
     1415{
     1416  LogIO os( LogOrigin( "PKSMS2reader", "setupAntennaList()", WHERE ) ) ;
     1417  //cerr << "antenna specification: " << s << endl ;
     1418  ROMSAntennaColumns antennaCols(cPKSMS.antenna());
     1419  ROScalarColumn<String> antNames = antennaCols.name();
     1420  Int nrow = antNames.nrow() ;
     1421  Vector<String> antlist = splitAntennaSelectionString( s ) ;
     1422  Int len = antlist.size() ;
     1423  Vector<Int> AntId( len ) ;
     1424  Regex re( "[0-9]+" ) ;
     1425  for ( Int i = 0 ; i < len ; i++ ) {
     1426    if ( antlist[i].matches( re ) ) {
     1427      AntId[i] = atoi( antlist[i].c_str() ) ;
     1428      if ( AntId[i] >= nrow ) {
     1429        os << LogIO::SEVERE << "Antenna index out of range: " << AntId[i] << LogIO::EXCEPTION ;
     1430      }
     1431    }
     1432    else {
     1433      AntId[i] = -1 ;
     1434      for ( uInt j = 0 ; j < antNames.nrow() ; j++ ) {
     1435        if ( antlist[i] == antNames(j) ) {
     1436          AntId[i] = j ;
     1437          break ;
     1438        }
     1439      }
     1440      if ( AntId[i] == -1 ) {
     1441        os << LogIO::SEVERE << "Specified antenna name not found: " << antlist[i] << LogIO::EXCEPTION ;
     1442      }
     1443    }
     1444  }
     1445  //cerr << "AntId = " << AntId << endl ;
     1446  vector<Int> uniqId ;
     1447  uniqId.push_back( AntId(0) ) ;
     1448  for ( uInt i = 1 ; i < AntId.size() ; i++ ) {
     1449    if ( count(uniqId.begin(),uniqId.end(),AntId[i]) == 0 ) {
     1450      uniqId.push_back( AntId[i] ) ;
     1451    }
     1452  }
     1453  Vector<Int> newAntId( uniqId ) ;
     1454  cAntId.assign( newAntId ) ;
     1455  //cerr << "cAntId = " << cAntId << endl ;
     1456}
Note: See TracChangeset for help on using the changeset viewer.