source: trunk/external-alma/asdm2ASAP/ASDMReader.cc @ 3106

Last change on this file since 3106 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No?

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...


Check-in asap modifications from Jim regarding casacore namespace conversion.
File size: 65.3 KB
Line 
1#include <iostream>
2#include <sstream>
3#include "limits.h"
4//#include "float.h"
5
6#include <measures/Measures/MEpoch.h>
7#include <measures/Measures/MPosition.h>
8#include <measures/Measures/MDirection.h>
9#include <measures/Measures/MCDirection.h>
10#include <measures/Measures/MeasFrame.h>
11#include <measures/Measures/MeasConvert.h>
12#include <casa/Logging/LogMessage.h>
13#include <casa/BasicSL/Constants.h>
14
15#include "ASDMReader.h"
16#include <atnf/PKSIO/SrcType.h>
17
18using namespace std ;
19//using namespace casa ;
20using namespace asdm ;
21using namespace sdmbin ;
22
23// sec to day
24double s2d = 1.0 / 86400.0 ;
25
26ASDMReader::ASDMReader()
27  : asdm_(0),
28    sdmBin_(0),
29    vmsData_(0),
30    antennaId_( -1 ),
31    antennaName_( "" ),
32    stationName_( "" ),
33    row_(-1),
34    apc_(AP_CORRECTED),
35    className_("ASDMReader"),
36    dataIndex_( UINT_MAX ),
37    //antennaRow_p( 0 ),
38    //stationRow_p( 0 ),
39    specWinRow_p( 0 ),
40    polarizationRow_p( 0 ),
41    fieldRow_p( 0 )
42{
43  configDescIdList_.resize(0) ;
44  feedIdList_.resize(0) ;
45  fieldIdList_.resize(0) ;
46  mainRow_.resize(0) ;
47  ifno_.clear() ;
48  corrMode_.reset() ;
49  timeSampling_.reset() ;
50  initConvert() ;
51}
52
53ASDMReader::~ASDMReader()
54{
55  close() ;
56  logsink_ = 0 ;
57}
58
59bool ASDMReader::open( const string &filename, const casacore::Record &rec )
60{
61  casacore::String funcName = "open" ;
62
63  // return value
64  bool status = true ;
65
66  // set default
67  timeSampling_.reset() ;
68  corrMode_.reset() ;
69  resolutionType_.reset() ;
70  apc_ = AP_CORRECTED ;
71
72  // parsing ASDM options
73  if ( rec.isDefined( "asdm" ) ) {
74    casacore::Record asdmrec = rec.asRecord( "asdm" ) ;
75
76    // antenna
77    if ( asdmrec.isDefined( "antenna" ) ) {
78      if ( asdmrec.type( asdmrec.fieldNumber( "antenna" ) ) == casacore::TpInt ) {
79        antennaId_ = asdmrec.asInt( "antenna" ) ;
80      }
81      else {
82        antennaName_ = asdmrec.asString( "antenna" ) ;
83      }
84    }
85    else {
86      antennaId_ = 0 ;
87    }
88
89    // ATM phase correction
90    if ( asdmrec.isDefined( "apc" ) ) {
91      if ( asdmrec.asBool( "apc" ) )
92        apc_ = AP_CORRECTED ;
93      else
94        apc_ = AP_UNCORRECTED ;
95    }
96
97    // time sampling
98    String timeSampling = "all" ; // take all time sampling by default
99    if ( asdmrec.isDefined( "sampling" ) ) {
100      timeSampling = asdmrec.asString( "sampling" ) ;
101    }
102    if ( timeSampling == "all" ) {
103      timeSampling_.set( INTEGRATION ) ;
104      timeSampling_.set( SUBINTEGRATION ) ;
105    }
106    else if ( timeSampling == "integration" ) {
107      timeSampling_.set( INTEGRATION ) ;
108    }
109    else if ( timeSampling == "subintegration" ) {
110      timeSampling_.set( SUBINTEGRATION ) ;
111    }
112    else {
113      //throw AipsError( "Unrecognized option for sampling" ) ;
114      logsink_->postLocally( LogMessage( "Unrecognized option for time sampling: "+String::toString(timeSampling), LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
115      status = false ;
116    }     
117
118    // spectral resolution type
119    string resolutionType = "all" ;
120    if ( asdmrec.isDefined( "srt" ) ) {
121      resolutionType = string( asdmrec.asString( "srt" ) ) ;
122    }
123    string resolutionTypes[3] ;
124    Int numType = split( resolutionType, resolutionTypes, 3, "," ) ;
125    for ( Int it = 0 ; it < numType ; it++ ) {
126      if ( resolutionTypes[it] == "all" ) {
127        resolutionType_.reset() ;
128        resolutionType_.set( FULL_RESOLUTION ) ;
129        resolutionType_.set( BASEBAND_WIDE ) ;
130        resolutionType_.set( CHANNEL_AVERAGE ) ;
131        break ;
132      }
133      else if ( resolutionTypes[it] == "fr" ) {
134        resolutionType_.set( FULL_RESOLUTION ) ;
135      }
136      else if ( resolutionTypes[it] == "bw" ) {
137        resolutionType_.set( BASEBAND_WIDE ) ;
138      }
139      else if ( resolutionTypes[it] == "ca" ) {
140        resolutionType_.set( CHANNEL_AVERAGE ) ;
141      }
142    }
143    if ( resolutionType_.size() == 0 ) {
144      logsink_->postLocally( LogMessage( "Unrecognized option for spectral resolution type: "+String::toString(resolutionType), LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
145      status = false ;
146    }
147   
148    // input correlation mode
149    string corrMode = "ao,ca" ;
150    if ( asdmrec.isDefined( "corr" ) ) {
151      corrMode = string( asdmrec.asString( "corr" ) ) ;
152      //logsink_->postLocally( LogMessage("corrMode = "+String(corrMode),LogOrigin(className_,funcName,WHERE)) ) ;
153    }
154    string corrModes[3] ;
155    Int numCorr = split( corrMode, corrModes, 3, "," ) ;
156    for ( Int ic = 0 ; ic < numCorr ; ic++ ) {
157      if ( corrModes[ic] == "ao" ) {
158        corrMode_.set( AUTO_ONLY ) ;
159      }
160      else if ( corrModes[ic] == "ca" ) {
161        corrMode_.set( CROSS_AND_AUTO ) ;
162      }
163    }
164    if ( corrMode_.size() == 0 ) {
165      logsink_->postLocally( LogMessage( "Invalid option for correlation mode: "+String::toString(corrMode), LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
166      status = false ;
167    }
168
169//     logsink_->postLocally( LogMessage( "### asdmrec summary ###", LogOrigin(className_,funcName,WHERE) ) ) ;
170//     ostringstream oss ;
171//     asdmrec.print( oss ) ;
172//     logsink_->postLocally( LogMessage( oss.str(), LogOrigin(className_,funcName,WHERE) ) ) ;
173//     logsink_->postLocally( LogMessage( "#######################", LogOrigin(className_,funcName,WHERE) ) ) ;
174  }
175
176  // create ASDM object
177  asdm_ = new ASDM() ;
178  asdm_->setFromFile( filename ) ;
179
180  AntennaTable &atab = asdm_->getAntenna() ;
181  AntennaRow *antennaRow ;
182  if ( antennaId_ == -1 ) {
183    vector<AntennaRow *> rows = atab.get() ;
184    int idx = -1 ;
185    for ( casacore::uInt irow = 0 ; irow < rows.size() ; irow++ ) {
186      if ( casacore::String(rows[irow]->getName()) == antennaName_ ) {
187        idx = rows[irow]->getAntennaId().getTagValue() ;
188        break ;
189      }
190    }
191    if ( idx == -1 ) {
192      close() ;
193      throw (casacore::AipsError( antennaName_ + " not found." )) ;
194    }
195    else {
196      antennaId_ = idx ;
197    }
198    antennaTag_ = Tag( antennaId_, TagType::Antenna ) ;
199    antennaRow = rows[antennaId_] ;
200  }
201  else {
202    antennaTag_ = Tag( antennaId_, TagType::Antenna ) ;
203    antennaRow = atab.getRowByKey( antennaTag_ ) ;
204    if ( antennaRow == 0 ) {
205      close() ;
206      throw (casacore::AipsError( "AntennaId " + casacore::String::toString(antennaId_) + " is invalid." ) ) ;
207    }
208  }
209 
210
211  // set antenna name
212  if ( antennaName_.size() == 0 ) {
213    antennaName_ = casacore::String( antennaRow->getName() ) ;
214  }
215
216  // get Station row
217  StationRow *stationRow = antennaRow->getStationUsingStationId() ;
218
219  // station name
220  stationName_ = casacore::String( stationRow->getName() ) ;
221
222  // antenna position
223  antennaPosition_.resize( 3 ) ;
224  vector<Length> antpos = stationRow->getPosition() ;
225  for ( casacore::uInt i = 0 ; i < 3 ; i++ )
226    antennaPosition_[i] = Quantity( casacore::Double( antpos[i].get() ), Unit( "m" ) ) ;
227  mp_ = casacore::MPosition( casacore::MVPosition( antennaPosition_ ),
228                         casacore::MPosition::ITRF ) ;
229  mf_.set( mp_ ) ;
230
231  // create SDMBinData object
232  sdmBin_ = new SDMBinData( asdm_, filename ) ;
233
234  // get Main rows
235  //mainRow_ = casacore::Vector<MainRow *>(asdm_->getMain().get()) ;
236
237  // set up IFNO
238  setupIFNO() ;
239
240  // process Station table
241  processStation() ;
242
243  //logsink_->postLocally( LogMessage(  "antennaId_ = "+String::toString(antennaId_), LogOrigin(className_,funcName,WHERE) ) ) ;
244  //logsink_->postLocally( LogMessage(  "antennaName_ = "+antennaName_, LogOrigin(className_,funcName,WHERE) ) ) ;
245
246  return true ;
247}
248
249// void ASDMReader::fill()
250// {
251// }
252
253void ASDMReader::close()
254{
255  clearMainRow() ;
256
257  if ( sdmBin_ )
258    delete sdmBin_ ;
259  sdmBin_ = 0 ;
260
261  if ( asdm_ )
262    delete asdm_ ;
263  asdm_ = 0 ;
264
265  return ;
266}
267
268void ASDMReader::fillHeader( casacore::Int &nchan,
269                             casacore::Int &npol,
270                             casacore::Int &nif,
271                             casacore::Int &nbeam,
272                             casacore::String &observer,
273                             casacore::String &project,
274                             casacore::String &obstype,
275                             casacore::String &antennaname,
276                             casacore::Vector<casacore::Double> &antennaposition,
277                             casacore::Float &equinox,
278                             casacore::String &freqref,
279                             casacore::Double &reffreq,
280                             casacore::Double &bandwidth,
281                             casacore::Double &utc,
282                             casacore::String &fluxunit,
283                             casacore::String &epoch,
284                             casacore::String &poltype )
285{
286  casacore::String funcName = "fillHeader" ;
287
288  ExecBlockTable &ebtab = asdm_->getExecBlock() ;
289  // at the moment take first row of ExecBlock table
290  ExecBlockRow *ebrow = ebtab.get()[0] ;
291  casacore::String telescopeName( ebrow->getTelescopeName() ) ;
292  //casacore::String stationName( stationRow_p->getName() ) ;
293
294  // antennaname
295  // <telescopeName>//<antennaName>@stationName
296  antennaname = telescopeName + "//" + antennaName_ + "@" + stationName_ ;
297  //logsink_->postLocally( LogMessage("antennaName = "+antennaname,LogOrigin(className_,funcName,WHERE)) ) ;
298
299  // antennaposition
300  antennaposition.resize( 3 ) ;
301  for ( casacore::uInt i = 0 ; i < 3 ; i++ )
302    antennaposition[i] = antennaPosition_[i].getValue( Unit("m") ) ;
303
304  // observer
305  observer = ebrow->getObserverName() ;
306
307  // project
308  // T.B.D. (project UID?)
309  project = "T.B.D. (" + ebrow->getProjectUID().toString() + ")" ;
310
311  // utc
312  // start time of the project
313  utc = casacore::Double( ebrow->getStartTime().getMJD() ) ;
314 
315
316  SpectralWindowTable &spwtab = asdm_->getSpectralWindow() ;
317  vector<SpectralWindowRow *> spwrows = spwtab.get() ;
318  int nspwrow = spwrows.size() ;
319
320  // nif
321  //nif = spwrows.size() ;
322  nif = getNumIFs() ;
323
324  // nchan
325  int refidx = -1 ;
326  vector<int> nchans ;
327  for ( int irow = 0 ; irow < nspwrow ; irow++ ) {
328    nchans.push_back( spwrows[irow]->getNumChan() ) ;
329    if ( refidx == -1 && nchans[irow] != 1 && nchans[irow] != 4 )
330      refidx = irow ;
331  }
332  nchan = casacore::Int( *max_element( nchans.begin(), nchans.end() ) ) ;
333
334  //logsink_->postLocally( LogMessage("refidx = "+String::toString(refidx),LogOrigin(className_,funcName,WHERE)) ) ;
335
336  // bandwidth
337  vector<double> bws ;
338  for ( int irow = 0 ; irow < nspwrow ; irow++ ) {
339    if ( nchans[irow] != 4 ) { // exclude WVR data
340      bws.push_back( spwrows[irow]->getTotBandwidth().get() ) ;
341    }
342  }
343  bandwidth = casacore::Double( *max_element( bws.begin(), bws.end() ) ) ;
344
345  // reffreq
346  reffreq = casacore::Double( spwrows[refidx]->getRefFreq().get() ) ;
347
348  // freqref
349  if ( spwrows[refidx]->isMeasFreqRefExists() ) {
350    string mfr = CFrequencyReferenceCode::name( spwrows[refidx]->getMeasFreqRef() ) ;
351//     if (mfr == "TOPO") {
352//       freqref = "TOPOCENT";
353//     } else if (mfr == "GEO") {
354//       freqref = "GEOCENTR";
355//     } else if (mfr == "BARY") {
356//       freqref = "BARYCENT";
357//     } else if (mfr == "GALACTO") {
358//       freqref = "GALACTOC";
359//     } else if (mfr == "LGROUP") {
360//       freqref = "LOCALGRP";
361//     } else if (mfr == "CMB") {
362//       freqref = "CMBDIPOL";
363//     } else if (mfr == "REST") {
364//       freqref = "SOURCE";
365//     }
366    freqref = String( mfr ) ;
367  }
368  else {
369    // frequency reference is TOPOCENT by default
370    //freqref = "TOPOCENT" ;
371    freqref = "TOPO" ;
372  }
373
374
375  PolarizationTable &ptab = asdm_->getPolarization() ;
376  vector<PolarizationRow *> prows = ptab.get() ;
377  vector<int> npols ;
378  refidx = -1 ;
379  for ( unsigned int irow = 0 ; irow < prows.size() ; irow++ ) {
380    npols.push_back( prows[irow]->getNumCorr() ) ;
381    if ( refidx == -1 && npols[irow] != 1 )
382      refidx = irow ;
383  }
384  if ( refidx == -1 )
385    refidx = 0 ;
386
387  // npol
388  npol = casacore::Int( *max_element( npols.begin(), npols.end() ) ) ;
389
390  // poltype
391  vector<StokesParameter> corrType = prows[refidx]->getCorrType() ;
392  if ( corrType[0] == I ||
393       corrType[0] == Q ||
394       corrType[0] == U ||
395       corrType[0] == V )
396    poltype = "stokes" ;
397  else if ( corrType[0] == RR ||
398            corrType[0] == RL ||
399            corrType[0] == LR ||
400            corrType[0] == LL )
401    poltype = "circular" ;
402  else if ( corrType[0] == XX ||
403            corrType[0] == XY ||
404            corrType[0] == YX ||
405            corrType[0] == YY )
406    poltype = "linear" ;
407  else if ( corrType[0] == PLINEAR ||
408            corrType[0] == PANGLE ) {
409    poltype = "linpol" ;
410  }
411  else {
412    poltype = "notype" ;
413  }
414
415
416  FeedTable &ftab = asdm_->getFeed() ;
417  vector<FeedRow *> frows = ftab.get() ;
418  vector<int> nbeams ;
419  for ( unsigned int irow = 0 ; irow < frows.size() ; irow++ ) {
420    if ( frows[irow]->isFeedNumExists() )
421      nbeams.push_back( frows[irow]->getFeedNum() ) ;
422    else
423      nbeams.push_back( 1 ) ;
424  }
425
426  // nbeam
427  nbeam = casacore::Int( *max_element( nbeams.begin(), nbeams.end() ) ) ;
428
429  // fluxunit
430  // tentatively set 'K'? or empty?
431  fluxunit = "K" ;
432
433  // equinox
434  // tentatively set 2000.0
435  equinox = 2000.0 ;
436
437  // epoch
438  // tentatively set "UTC"
439  epoch = "UTC" ;
440
441  // obstype
442  // at the moment take observingMode attribute in SBSummary table
443  SBSummaryRow *sbrow = ebrow->getSBSummaryUsingSBSummaryId() ;
444  vector<string> obsmode = sbrow->getObservingMode() ;
445  obstype = "" ;
446  for ( unsigned int imode = 0 ; imode < obsmode.size() ; imode++ ) {
447    obstype += casacore::String(obsmode[imode]) ;
448    if ( imode != obsmode.size()-1 )
449      obstype += "#" ;
450  }
451}
452
453void ASDMReader::selectConfigDescription()
454{
455  casacore::String funcName = "selectConfigDescription" ;
456
457  vector<ConfigDescriptionRow *> cdrows = asdm_->getConfigDescription().get() ;
458  vector<Tag> cdidTags ;
459  for ( unsigned int irow = 0 ; irow < cdrows.size() ; irow++ ) {
460    //logsink_->postLocally( LogMessage("correlationMode["+String::toString(irow)+"] = "+String::toString(cdrows[irow]->getCorrelationMode()),LogOrigin(className_,funcName,WHERE)) ) ;
461    if ( cdrows[irow]->getCorrelationMode() != CROSS_ONLY ) {
462      cdidTags.push_back( cdrows[irow]->getConfigDescriptionId() ) ;
463    }
464  }
465
466  configDescIdList_.resize( cdidTags.size() ) ;
467  for ( unsigned int i = 0 ; i < cdidTags.size() ; i++ ) {
468    configDescIdList_[i] = casacore::uInt( cdidTags[i].getTagValue() ) ;
469  }
470}
471
472void ASDMReader::selectFeed()
473{
474  feedIdList_.resize(0) ;
475  vector<FeedRow *> frows = asdm_->getFeed().get() ;
476  Tag atag( antennaId_, TagType::Antenna ) ;
477  for ( unsigned int irow = 0 ; irow < frows.size() ; irow++ ) {
478    casacore::uInt feedId = (casacore::uInt)(frows[irow]->getFeedId() ) ;
479    if ( casacore::anyEQ( feedIdList_, feedId ) )
480      continue ;
481    if ( frows[irow]->getAntennaId() == atag ) {
482      unsigned int oldsize = feedIdList_.size() ;
483      feedIdList_.resize( oldsize+1, true ) ;
484      feedIdList_[oldsize] = feedId ;
485    }
486  }
487}
488
489casacore::Vector<casacore::uInt> ASDMReader::getFieldIdList()
490{
491  casacore::String funcName = "getFieldIdList" ;
492
493  vector<FieldRow *> frows = asdm_->getField().get() ;
494  fieldIdList_.resize( frows.size() ) ;
495  for ( unsigned int irow = 0 ; irow < frows.size() ; irow++ ) {
496    //logsink_->postLocally( LogMessage("fieldId["+String::toString(irow)+"]="+String(frows[irow]->getFieldId().toString()),LogOrigin(className_,funcName,WHERE)) ) ;
497    fieldIdList_[irow] = frows[irow]->getFieldId().getTagValue() ;
498  }
499
500  return fieldIdList_ ;
501}
502
503casacore::uInt ASDMReader::getNumMainRow()
504{
505  casacore::uInt nrow = casacore::uInt( mainRow_.size() ) ;
506
507  return nrow ;
508}
509
510void ASDMReader::select()
511{
512  // selection by input CorrelationMode
513  EnumSet<CorrelationMode> esCorrs ;
514  sdmBin_->select( corrMode_ ) ;
515
516  // selection by TimeSampling
517  sdmBin_->select( timeSampling_ ) ;
518
519  // selection by SpectralResolutionType
520  sdmBin_->select( resolutionType_ ) ;
521
522  // selection by AtmPhaseCorrection and output CorrelationMode
523  EnumSet<AtmPhaseCorrection> esApcs ;
524  esApcs.set( apc_ ) ;
525  // always take only autocorrelation data
526  Enum<CorrelationMode> esCorr = AUTO_ONLY ;
527  sdmBin_->selectDataSubset( esCorr, esApcs ) ;
528
529  // select valid configDescriptionId
530  selectConfigDescription() ;
531
532  // select valid feedId
533  selectFeed() ;
534}
535
536casacore::Bool ASDMReader::setMainRow( casacore::uInt irow )
537{
538  casacore::Bool status = true ;
539  row_ = irow ;
540  execBlockTag_ = mainRow_[row_]->getExecBlockId() ;
541
542  unsigned int cdid = mainRow_[row_]->getConfigDescriptionId().getTagValue() ;
543  if ( (int)count(configDescIdList_.begin(),configDescIdList_.end(),cdid) == 0 )
544    status = false ;
545  else {
546    status = (casacore::Bool)(sdmBin_->acceptMainRow( mainRow_[row_] )) ;
547  }
548  return status ;
549}
550
551casacore::Bool ASDMReader::setMainRow( casacore::uInt configDescId, casacore::uInt fieldId )
552{
553  clearMainRow() ;
554
555  Tag configDescTag( (unsigned int)configDescId, TagType::ConfigDescription ) ;
556  Tag fieldTag( (unsigned int)fieldId, TagType::Field ) ;
557  vector<MainRow *> *rows = asdm_->getMain().getByContext( configDescTag, fieldTag );
558  if (rows == 0)
559    return false;
560  mainRow_ = casacore::Vector<MainRow *>( *rows ) ;
561 
562  return true ;
563}
564
565void ASDMReader::clearMainRow()
566{
567  mainRow_.resize(0) ;
568}
569
570void ASDMReader::setupIFNO()
571{
572  casacore::String funcName = "setupIFNO" ;
573
574  vector<SpectralWindowRow *> spwrows = asdm_->getSpectralWindow().get() ;
575  unsigned int nrow = spwrows.size() ;
576  ifno_.clear() ;
577  casacore::uInt idx = 0 ;
578  casacore::uInt wvridx = 0 ;
579  for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
580    casacore::uInt index ;
581    if ( isWVR( spwrows[irow] ) ) {
582      //logsink_->postLocally( LogMessage(spwrows[irow]->getSpectralWindowId().toString()+" is WVR",LogOrigin(className_,funcName,WHERE)) ) ;
583      index = wvridx ;
584    }
585    else {
586      index = ++idx ;
587    }
588    ifno_.insert( pair<Tag,casacore::uInt>(spwrows[irow]->getSpectralWindowId(),index) ) ;
589    //logsink_->postLocally( LogMessage(spwrows[irow]->getSpectralWindowId().toString()+": IFNO="+String::toString(index),LogOrigin(className_,funcName,WHERE)) ) ;
590  }
591}
592
593bool ASDMReader::isWVR( SpectralWindowRow *row )
594{
595  BasebandName bbname = row->getBasebandName() ;
596  int nchan = row->getNumChan() ;
597  if ( bbname == NOBB && nchan == 4 )
598    return true ;
599  else
600    return false ;
601}
602
603casacore::Bool ASDMReader::setData()
604{
605  casacore::String funcName = "setData" ;
606
607  //logsink_->postLocally( LogMessage("try to retrieve binary data",LogOrigin(className_,funcName,WHERE)) ) ;
608 
609//   EnumSet<AtmPhaseCorrection> esApcs ;
610//   esApcs.set( apc_ ) ;
611//   // always take only autocorrelation data
612//   Enum<CorrelationMode> esCorr = AUTO_ONLY ;
613//   vmsData_ = sdmBin_->getDataCols( esCorr, esApcs ) ;
614
615  // 2011/07/06 TN
616  // Workaround to avoid unwanted message from SDMBinData::getDataCols()
617  ostringstream oss ;
618  streambuf *buforg = cout.rdbuf(oss.rdbuf()) ;
619  vmsData_ = sdmBin_->getDataCols() ;
620  cout.rdbuf(buforg) ;
621
622//   logsink_->postLocally( LogMessage("oss.str() = "+oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
623//   cout << "This is test: oss.str()=" << oss.str() << endl ;
624
625
626  //logsink_->postLocally( LogMessage("processorId = "+String::toString(vmsData_->processorId),LogOrigin(className_,funcName,WHERE)) ) ;
627  //logsink_->postLocally( LogMessage("v_time.size() = "+String::toString(vmsData_->v_time.size()),LogOrigin(className_,funcName,WHERE)) ) ;
628  //logsink_->postLocally( LogMessage("   v_time[0] = "+String::toString(vmsData_->v_time[0]),LogOrigin(className_,funcName,WHERE)) ) ;
629  //logsink_->postLocally( LogMessage("v_interval.size() = "+String::toString(vmsData_->v_interval.size()),LogOrigin(className_,funcName,WHERE)) ) ;
630  //logsink_->postLocally( LogMessage("   v_interval[0] = "+String::toString(vmsData_->v_interval[0]),LogOrigin(className_,funcName,WHERE)) ) ;
631  //logsink_->postLocally( LogMessage("v_atmPhaseCorrection.size() = "+String::toString(vmsData_->v_atmPhaseCorrection.size()),LogOrigin(className_,funcName,WHERE)) ) ;
632  //logsink_->postLocally( LogMessage("binNum = "+String::toString(vmsData_->binNum),LogOrigin(className_,funcName,WHERE)) ) ;
633  //logsink_->postLocally( LogMessage("v_projectPath.size() = "+String::toString(vmsData_->v_projectPath.size()),LogOrigin(className_,funcName,WHERE)) ) ;
634  //logsink_->postLocally( LogMessage("v_antennaId1.size() = "+String::toString(vmsData_->v_antennaId1.size()),LogOrigin(className_,funcName,WHERE)) ) ;
635  //logsink_->postLocally( LogMessage("v_antennaId2.size() = "+String::toString(vmsData_->v_antennaId2.size()),LogOrigin(className_,funcName,WHERE)) ) ;
636  //logsink_->postLocally( LogMessage("v_feedId1.size() = "+String::toString(vmsData_->v_feedId1.size()),LogOrigin(className_,funcName,WHERE)) ) ;
637  //logsink_->postLocally( LogMessage("v_feedId2.size() = "+String::toString(vmsData_->v_feedId2.size()),LogOrigin(className_,funcName,WHERE)) ) ;
638  //logsink_->postLocally( LogMessage("v_dataDescId.size() = "+String::toString(vmsData_->v_dataDescId.size()),LogOrigin(className_,funcName,WHERE)) ) ;
639  //logsink_->postLocally( LogMessage("v_timeCentroid.size() = "+String::toString(vmsData_->v_timeCentroid.size()),LogOrigin(className_,funcName,WHERE)) ) ;
640  //logsink_->postLocally( LogMessage("v_exposure.size() = "+String::toString(vmsData_->v_exposure.size()),LogOrigin(className_,funcName,WHERE)) ) ;
641  //logsink_->postLocally( LogMessage("v_numData.size() = "+String::toString(vmsData_->v_numData.size()),LogOrigin(className_,funcName,WHERE)) ) ;
642  //logsink_->postLocally( LogMessage("vv_dataShape.size() = "+String::toString(vmsData_->vv_dataShape.size()),LogOrigin(className_,funcName,WHERE)) ) ;
643  //logsink_->postLocally( LogMessage("v_m_data.size() = "+String::toString(vmsData_->v_m_data.size()),LogOrigin(className_,funcName,WHERE)) ) ;
644  //logsink_->postLocally( LogMessage("v_phaseDir.size() = "+String::toString(vmsData_->v_phaseDir.size()),LogOrigin(className_,funcName,WHERE)) ) ;
645  //logsink_->postLocally( LogMessage("v_stateId.size() = "+String::toString(vmsData_->v_stateId.size()),LogOrigin(className_,funcName,WHERE)) ) ;
646  //logsink_->postLocally( LogMessage("v_msState.size() = "+String::toString(vmsData_->v_msState.size()),LogOrigin(className_,funcName,WHERE)) ) ;
647  //logsink_->postLocally( LogMessage("v_flag.size() = "+String::toString(vmsData_->v_flag.size()),LogOrigin(className_,funcName,WHERE)) ) ;
648
649  dataIdList_.clear() ;
650  unsigned int numTotalData = vmsData_->v_m_data.size() ;
651  for ( unsigned int idata = 0 ; idata < numTotalData ; idata++ ) {
652    if ( vmsData_->v_antennaId1[idata] == (int)antennaId_
653         && vmsData_->v_antennaId2[idata] == (int)antennaId_ )
654      dataIdList_.push_back( idata ) ;
655  }
656  numData_ = dataIdList_.size() ;
657  //logsink_->postLocally( LogMessage("numData_ = "+String::toString(numData_),LogOrigin(className_,funcName,WHERE)) ) ;
658  //logsink_->postLocally( LogMessage("dataSize = "+String::toString(mainRow_[row_]->getDataSize()),LogOrigin(className_,funcName,WHERE)) ) ;
659
660  return true ;
661}
662
663void ASDMReader::prepareData( unsigned int idx )
664{
665  unsigned int i = dataIdList_[idx] ;
666  if ( i != dataIndex_ ) {
667    dataIndex_ = dataIdList_[idx] ;
668    Tag dataDescTag( vmsData_->v_dataDescId[dataIndex_], TagType::DataDescription ) ;
669    DataDescriptionRow *dataDescRow = asdm_->getDataDescription().getRowByKey( dataDescTag ) ;
670    specWinTag_ = dataDescRow->getSpectralWindowId() ;
671    specWinRow_p = dataDescRow->getSpectralWindowUsingSpectralWindowId() ;
672    polarizationRow_p = dataDescRow->getPolarizationUsingPolOrHoloId() ;
673    Tag fieldTag( vmsData_->v_fieldId[dataIndex_], TagType::Field ) ;
674    fieldRow_p = asdm_->getField().getRowByKey( fieldTag ) ;
675    double startSec = vmsData_->v_time[dataIndex_] - 0.5 * vmsData_->v_interval[dataIndex_] ;
676    timeInterval_ = ArrayTimeInterval( startSec*s2d, vmsData_->v_interval[dataIndex_]*s2d ) ;
677  }
678}
679
680casacore::uInt ASDMReader::getIFNo( unsigned int idx )
681{
682  prepareData( idx ) ;
683  return getIFNo() ;
684}
685
686casacore::uInt ASDMReader::getIFNo()
687{
688  map<Tag,casacore::uInt>::iterator iter = ifno_.find( specWinTag_ ) ;
689  if ( iter != ifno_.end() )
690    return iter->second ;
691  else {
692    return 0 ;
693  }
694}
695
696int ASDMReader::getNumPol( unsigned int idx )
697{
698  prepareData( idx ) ;
699  return getNumPol() ;
700}
701
702int ASDMReader::getNumPol()
703{
704  return polarizationRow_p->getNumCorr() ;
705}
706
707void ASDMReader::getFrequency( unsigned int idx,
708                               double &refpix,
709                               double &refval,
710                               double &incr,
711                               string &freqref )
712{
713  prepareData( idx ) ;
714  getFrequency( refpix, refval, incr, freqref ) ;
715}
716
717void ASDMReader::getFrequency( double &refpix,
718                               double &refval,
719                               double &incr,
720                               string &freqref )
721{
722  casacore::String funcName = "getFrequency" ;
723
724  int nchan = specWinRow_p->getNumChan() ;
725  freqref = "TOPO" ;
726  if ( specWinRow_p->isMeasFreqRefExists() )
727    freqref = CFrequencyReferenceCode::toString( specWinRow_p->getMeasFreqRef() ) ;
728  if ( nchan == 1 ) {
729    //logsink_->postLocally( LogMessage("channel averaged data",LogOrigin(className_,funcName,WHERE)) ) ;
730    refpix = 0.0 ;
731    incr = specWinRow_p->getTotBandwidth().get() ;
732    if ( specWinRow_p->isChanFreqStartExists() ) {
733      refval = specWinRow_p->getChanFreqStart().get() ;
734    }
735    else if ( specWinRow_p->isChanFreqArrayExists() ) {
736      refval = specWinRow_p->getChanFreqArray()[0].get() ;
737    }
738    else {
739      throw (casacore::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
740    }     
741  }
742  else if ( nchan % 2 ) {
743    // odd
744    //logsink_->postLocally( LogMessage("odd case",LogOrigin(className_,funcName,WHERE)) ) ;
745    refpix = 0.5 * ( (double)nchan - 1.0 ) ;
746    int ic = ( nchan - 1 ) / 2 ;
747    if ( specWinRow_p->isChanWidthExists() ) {
748      incr = specWinRow_p->getChanWidth().get() ;
749    }
750    else if ( specWinRow_p->isChanWidthArrayExists() ) {
751      incr = specWinRow_p->getChanWidthArray()[0].get() ;
752    }
753    else {
754      throw (casacore::AipsError( "Either chanWidthArray or chanWidth must exist." )) ;
755    }     
756    if ( specWinRow_p->isChanFreqStepExists() ) {
757      if ( specWinRow_p->getChanFreqStep().get() < 0.0 )
758        incr *= -1.0 ;
759    }
760    else if ( specWinRow_p->isChanFreqArrayExists() ) {
761      vector<Frequency> chanFreqArr = specWinRow_p->getChanFreqArray() ;
762      if ( chanFreqArr[0].get() > chanFreqArr[1].get() )
763        incr *= -1.0 ;
764    }
765    else {
766      throw (casacore::AipsError( "Either chanFreqArray or chanFreqStep must exist." )) ;
767    }         
768    if ( specWinRow_p->isChanFreqStartExists() ) {
769      refval = specWinRow_p->getChanFreqStart().get() + refpix * incr ;
770    }
771    else if ( specWinRow_p->isChanFreqArrayExists() ) {
772      refval = specWinRow_p->getChanFreqArray()[ic].get() ;
773    }
774    else {
775      throw (casacore::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
776    }     
777  }
778  else {
779    // even
780    //logsink_->postLocally( LogMessage("even case",LogOrigin(className_,funcName,WHERE)) ) ;
781    refpix = 0.5 * ( (double)nchan - 1.0 ) ;
782    int ic = nchan / 2 ;
783    if ( specWinRow_p->isChanWidthExists() ) {
784      incr = specWinRow_p->getChanWidth().get() ;
785    }
786    else if ( specWinRow_p->isChanWidthArrayExists() ) {
787      incr = specWinRow_p->getChanWidthArray()[0].get() ;
788    }
789    else {
790      throw (casacore::AipsError( "Either chanWidthArray or chanWidth must exist." )) ;
791    }
792    if ( specWinRow_p->isChanFreqStepExists() ) {
793      if ( specWinRow_p->getChanFreqStep().get() < 0.0 )
794        incr *= -1.0 ;
795    }
796    else if ( specWinRow_p->isChanFreqArrayExists() ) {
797      vector<Frequency> chanFreqArr = specWinRow_p->getChanFreqArray() ;
798      if ( chanFreqArr[0].get() > chanFreqArr[1].get() )
799        incr *= -1.0 ;
800    }
801    else {
802      throw (casacore::AipsError( "Either chanFreqArray or chanFreqStep must exist." )) ;
803    }         
804    if ( specWinRow_p->isChanFreqStartExists() ) {
805      refval = specWinRow_p->getChanFreqStart().get() + refpix * incr ;
806    }
807    else if ( specWinRow_p->isChanFreqArrayExists() ) {
808      vector<Frequency> freqs = specWinRow_p->getChanFreqArray() ;
809      refval = 0.5 * ( freqs[ic-1].get() + freqs[ic].get() ) ;     
810    }
811    else {
812      throw (casacore::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
813    }     
814  }
815}
816
817double ASDMReader::getTime( unsigned int idx )
818{
819  prepareData( idx ) ;
820  return getTime() ;
821}
822
823double ASDMReader::getTime()
824{
825  double tsec = vmsData_->v_time[dataIndex_] ;
826  return tsec * s2d ;
827}
828
829double ASDMReader::getInterval( unsigned int idx )
830{
831  prepareData( idx ) ;
832  return getInterval() ;
833}
834
835double ASDMReader::getInterval()
836{
837  return vmsData_->v_interval[dataIndex_] ;
838}
839
840void ASDMReader::getSourceProperty( unsigned int idx,
841                                    string &srcname,
842                                    string &fieldname,
843                                    vector<double> &srcdir,
844                                    vector<double> &srcpm,
845                                    double &sysvel,
846                                    vector<double> &restfreq )
847{
848  prepareData( idx ) ;
849  getSourceProperty( srcname, fieldname, srcdir, srcpm, sysvel, restfreq ) ;
850}
851
852void ASDMReader::getSourceProperty( string &srcname,
853                                    string &fieldname,
854                                    vector<double> &srcdir,
855                                    vector<double> &srcpm,
856                                    double &sysvel,
857                                    vector<double> &restfreq )
858{
859  String funcName = "getSourceProperty" ;
860  ostringstream oss ;
861  oss << fieldRow_p->getFieldName() << "__" << vmsData_->v_fieldId[dataIndex_] ;
862  fieldname = oss.str() ;
863  SourceRow *srow = 0 ;
864  if ( fieldRow_p->isSourceIdExists() ) {
865    int sourceId = fieldRow_p->getSourceId() ;
866    srow = asdm_->getSource().getRowByKey( sourceId, timeInterval_, specWinTag_ ) ;
867  }
868  if ( srow != 0 ) {
869    //logsink_->postLocally( LogMessage("timeInterval_="+String::toString(timeInterval_.getStart().getMJD()),LogOrigin(className_,funcName,WHERE)) ) ;
870    //logsink_->postLocally( LogMessage("specWinTag_="+String::toString(specWinTag_.toString()),LogOrigin(className_,funcName,WHERE)) ) ;
871    //SourceRow *srow = asdm_->getSource().getRowByKey( sourceId, timeInterval_, specWinTag_ ) ;
872    //logsink_->postLocally( LogMessage("sourceId="+String::toString(sourceId),LogOrigin(className_,funcName,WHERE)) ) ;
873    //if ( srow == 0 )
874    //logsink_->postLocally( LogMessage("nullpo",LogOrigin(className_,funcName,WHERE)) ) ;
875
876    // source name
877    srcname = srow->getSourceName() ;
878    //logsink_->postLocally( LogMessage("srcname="+String::toString(srcname),LogOrigin(className_,funcName,WHERE)) ) ;
879
880    // source direction
881    vector<Angle> srcdirA = srow->getDirection() ;
882    srcdir.resize( 2 ) ;
883    srcdir[0] = limitedAngle( srcdirA[0].get() ) ;
884    srcdir[1] = limitedAngle( srcdirA[1].get() ) ;
885    //logsink_->postLocally( LogMessage("srcdir=["+String::toString(srcdir[0])+","+String::toString(srcdir[1])+"]",LogOrigin(className_,funcName,WHERE)) ) ;
886    if ( srow->isDirectionCodeExists() ) {
887      DirectionReferenceCodeMod::DirectionReferenceCode dircode = srow->getDirectionCode() ;
888      //logsink_->postLocally( LogMessage("dircode="+CDirectionReferenceCode::toString(dircode),LogOrigin(className_,funcName,WHERE)) ) ;
889      if ( dircode != DirectionReferenceCodeMod::J2000 ) {
890        // if not J2000, need direction conversion
891        String ref( CDirectionReferenceCode::toString( dircode ) ) ;
892        double mjd = vmsData_->v_time[dataIndex_] * s2d ;
893        srcdir = toJ2000( srcdir, ref, mjd ) ;
894      }
895    }
896
897    // source proper motion
898    srcpm.resize( 2 ) ;
899    vector<AngularRate> srcpmA = srow->getProperMotion() ;
900    srcpm[0] = srcpmA[0].get() ;
901    srcpm[1] = srcpmA[1].get() ;
902    //logsink_->postLocally( LogMessage("srcpm=["+String::toString(srcpm[0])+","+String::toString(srcpm[1])+"]",LogOrigin(className_,funcName,WHERE)) ) ;
903    // systemic velocity
904    if ( srow->isSysVelExists() ) {
905      vector<Speed> sysvelV = srow->getSysVel() ;
906      if ( sysvelV.size() > 0 )
907        sysvel = sysvelV[0].get() ;
908    }
909    else {
910      sysvel = 0.0 ;
911    }
912    //logsink_->postLocally( LogMessage("sysvel="+String::toString(sysvel),LogOrigin(className_,funcName,WHERE)) ) ;
913    // rest frequency
914    if ( srow->isRestFrequencyExists() ) {
915      //logsink_->postLocally( LogMessage("restFrequency exists",LogOrigin(className_,funcName,WHERE)) ) ;
916      vector<Frequency> rf = srow->getRestFrequency() ;
917      restfreq.resize( rf.size() ) ;
918      for ( unsigned int i = 0 ; i < restfreq.size() ; i++ )
919        restfreq[i] = rf[i].get() ;
920    }
921    else {
922      restfreq.resize( 0 ) ;
923    }
924    //logsink_->postLocally( LogMessage("restfreq.size()="+String::toString(restfreq.size()),LogOrigin(className_,funcName,WHERE)) ) ;
925  }
926  else {
927    srcname = fieldRow_p->getFieldName() ;
928    srcdir.resize( 2 ) ;
929    srcdir[0] = 0.0 ;
930    srcdir[1] = 0.0 ;
931    srcpm = srcdir ;
932    sysvel = 0.0 ;
933    restfreq.resize( 0 ) ;
934  }
935}
936
937int ASDMReader::getSrcType( unsigned int scan,
938                            unsigned int subscan )
939{
940  int srctype = SrcType::NOTYPE ;
941  ScanRow *scanrow = asdm_->getScan().getRowByKey( execBlockTag_, (int)scan ) ;
942  ScanIntentMod::ScanIntent scanIntent = scanrow->getScanIntent()[0] ;
943  SubscanRow *subrow = asdm_->getSubscan().getRowByKey( execBlockTag_, (int)scan, (int)subscan ) ;
944  if ( subrow == 0 ) {
945    return srctype ;
946  }
947  SubscanIntentMod::SubscanIntent subIntent = subrow->getSubscanIntent() ;
948  SwitchingModeMod::SwitchingMode swmode = SwitchingModeMod::NO_SWITCHING ;
949  if ( subrow->isSubscanModeExists() )
950    swmode = subrow->getSubscanMode() ;
951  if ( scanIntent == ScanIntentMod::OBSERVE_TARGET ) {
952    // on sky scan
953    if ( swmode == SwitchingModeMod::NO_SWITCHING || swmode == SwitchingModeMod::POSITION_SWITCHING ) {
954      // position switching
955      // tentatively set NO_SWITCHING = POSITION_SWITCHING
956      if ( subIntent == SubscanIntentMod::ON_SOURCE ) {
957        srctype = SrcType::PSON ;
958      }
959      else if ( subIntent == SubscanIntentMod::OFF_SOURCE ) {
960        srctype = SrcType::PSOFF ;
961      }
962    }
963    else if ( swmode == SwitchingModeMod::FREQUENCY_SWITCHING ) {
964      // frequency switching
965      if ( subIntent == SubscanIntentMod::ON_SOURCE ) {
966        srctype = SrcType::FSON ;
967      }
968      else if ( subIntent == SubscanIntentMod::OFF_SOURCE ) {
969        srctype = SrcType::FSOFF ;
970      }
971    }
972    else if ( swmode == SwitchingModeMod::NUTATOR_SWITCHING ) {
973      // nutator switching
974      if ( subIntent == SubscanIntentMod::ON_SOURCE ) {
975        srctype = SrcType::PSON ;
976      }
977      else if ( subIntent == SubscanIntentMod::OFF_SOURCE ) {
978        srctype = SrcType::PSOFF ;
979      }
980    }
981    else {
982      // other switching mode
983      if ( subIntent == SubscanIntentMod::ON_SOURCE ) {
984        srctype = SrcType::SIG ;
985      }
986      else if ( subIntent == SubscanIntentMod::OFF_SOURCE ) {
987        srctype = SrcType::REF ;
988      }
989    }
990  }
991  else if ( scanIntent == ScanIntentMod::CALIBRATE_ATMOSPHERE ) {
992    if ( swmode == SwitchingModeMod::NO_SWITCHING || swmode == SwitchingModeMod::POSITION_SWITCHING ) {
993      // position switching
994      // tentatively set NO_SWITCHING = POSITION_SWITCHING
995      if ( subIntent == SubscanIntentMod::ON_SOURCE ) {
996        srctype = SrcType::PONCAL ;
997      }
998      else if ( subIntent == SubscanIntentMod::OFF_SOURCE ) {
999        srctype = SrcType::POFFCAL ;
1000      }
1001    }
1002    else if ( swmode == SwitchingModeMod::FREQUENCY_SWITCHING ) {
1003      // frequency switching
1004      if ( subIntent == SubscanIntentMod::ON_SOURCE ) {
1005        srctype = SrcType::FONCAL ;
1006      }
1007      else if ( subIntent == SubscanIntentMod::OFF_SOURCE ) {
1008        srctype = SrcType::FOFFCAL ;
1009      }
1010    }
1011    else if ( swmode == SwitchingModeMod::NUTATOR_SWITCHING ) {
1012      // nutator switching
1013      if ( subIntent == SubscanIntentMod::ON_SOURCE ) {
1014        srctype = SrcType::PONCAL ;
1015      }
1016      else if ( subIntent == SubscanIntentMod::OFF_SOURCE ) {
1017        srctype = SrcType::POFFCAL ;
1018      }
1019    }
1020    else {
1021      // other switching mode
1022      if ( subIntent == SubscanIntentMod::ON_SOURCE ) {
1023        srctype = SrcType::CAL ;
1024      }
1025      else if ( subIntent == SubscanIntentMod::OFF_SOURCE ) {
1026        srctype = SrcType::CAL ;
1027      }
1028    }
1029  }
1030  else {
1031    // off sky (e.g. calibrator device) scan
1032    if ( subIntent == SubscanIntentMod::ON_SOURCE ) {
1033      srctype = SrcType::SIG ;
1034    }
1035    else if ( subIntent == SubscanIntentMod::OFF_SOURCE ) {
1036      srctype = SrcType::REF ;
1037    }
1038    else if ( subIntent == SubscanIntentMod::HOT ) {
1039      srctype = SrcType::HOT ;
1040    }
1041    else if ( subIntent == SubscanIntentMod::AMBIENT ) {
1042      srctype = SrcType::SKY ;
1043    }
1044    else {
1045      srctype = SrcType::CAL ;
1046    }
1047  }
1048
1049  return srctype ;
1050}
1051
1052unsigned int ASDMReader::getSubscanNo( unsigned int idx )
1053{
1054  prepareData( idx ) ;
1055  return getSubscanNo() ;
1056}
1057
1058unsigned int ASDMReader::getSubscanNo()
1059{
1060  return mainRow_[row_]->getSubscanNumber() ;
1061}
1062
1063void ASDMReader::getSourceDirection( vector<double> &dir, string &ref )
1064{
1065  dir.resize( 2 ) ;
1066  ref = "J2000" ; // default is J2000
1067  SourceTable &tab = asdm_->getSource() ;
1068  SourceRow *row = tab.get()[0] ;
1069  vector<Angle> dirA = row->getDirection() ;
1070  dir[0] = dirA[0].get() ;
1071  dir[1] = dirA[1].get() ;
1072  if ( row->isDirectionCodeExists() ) {
1073    ref = CDirectionReferenceCode::toString( row->getDirectionCode() ) ;
1074  }
1075}
1076
1077unsigned int ASDMReader::getFlagRow( unsigned int idx )
1078{
1079  prepareData( idx ) ;
1080  return getFlagRow() ;
1081}
1082
1083unsigned int ASDMReader::getFlagRow()
1084{
1085  return vmsData_->v_flag[dataIndex_] ;
1086}
1087
1088vector<unsigned int> ASDMReader::getDataShape( unsigned int idx )
1089{
1090  prepareData( idx ) ;
1091  return getDataShape() ;
1092}
1093
1094vector<unsigned int> ASDMReader::getDataShape()
1095{
1096  return vmsData_->vv_dataShape[dataIndex_] ;
1097}
1098
1099float * ASDMReader::getSpectrum( unsigned int idx )
1100{
1101  prepareData( idx ) ;
1102  return getSpectrum() ;
1103}
1104
1105float * ASDMReader::getSpectrum()
1106{
1107  map<AtmPhaseCorrection, float*> data = vmsData_->v_m_data[dataIndex_] ;
1108  //map<AtmPhaseCorrection, float*>::iterator iter = data.find(AP_UNCORRECTED) ;
1109  map<AtmPhaseCorrection, float*>::iterator iter = data.find(apc_) ;
1110  float *autoCorr = iter->second ;
1111  return autoCorr ;
1112}
1113
1114vector< vector<float> > ASDMReader::getTsys( unsigned int idx )
1115{
1116  prepareData( idx ) ;
1117  return getTsys() ;
1118}
1119
1120vector< vector<float> > ASDMReader::getTsys()
1121{
1122  vector< vector<float> > defaultTsys( 1, vector<float>( 1, 1.0 ) ) ;
1123  SysCalRow *scrow = getSysCalRow() ;
1124  if ( scrow != 0 && scrow->isTsysSpectrumExists() ) {
1125    vector< vector<Temperature> > tsysSpec = scrow->getTsysSpectrum() ;
1126    unsigned int numReceptor = tsysSpec.size() ;
1127    unsigned int numChan = tsysSpec[0].size() ;
1128    vector< vector<float> > tsys( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1129    for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1130      for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1131        tsys[ir][ic] = tsysSpec[ir][ic].get() ;
1132      }
1133    }
1134    return tsys ;
1135  }
1136//   else if ( scrow->isTsysExists() ) {
1137//     vector<Temperature> tsysScalar = scrow->getTsys() ;
1138//     unsigned int numReceptor = tsysScalar.size() ;
1139//     vector< vector<float> > tsys( numReceptor ) ;
1140//     for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1141//       tsys[ir] = vector<float>( 1, tsysScalar[ir].get() ) ;
1142//     return tsys ;
1143//   }
1144  else {
1145    return defaultTsys ;
1146  }
1147}
1148
1149vector< vector<float> > ASDMReader::getTcal( unsigned int idx )
1150{
1151  prepareData( idx ) ;
1152  return getTcal() ;
1153}
1154
1155vector< vector<float> > ASDMReader::getTcal()
1156{
1157  vector< vector<float> > defaultTcal( 1, vector<float>( 1, 1.0 ) ) ;
1158  SysCalRow *scrow = getSysCalRow() ;
1159  if ( scrow != 0 && scrow->isTcalSpectrumExists() ) {
1160    vector< vector<Temperature> > tcalSpec = scrow->getTcalSpectrum() ;
1161    unsigned int numReceptor = tcalSpec.size() ;
1162    unsigned int numChan = tcalSpec[0].size() ;
1163    vector< vector<float> > tcal( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1164    for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1165      for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1166        tcal[ir][ic] = tcalSpec[ir][ic].get() ;
1167      }
1168    }
1169    return tcal ;
1170  }
1171//   else if ( scrow->isTcalExists() ) {
1172//     vector<Temperature> tcalScalar = scrow->getTcal() ;
1173//     unsigned int numReceptor = tcalScalar.size() ;
1174//     vector< vector<float> > tcal( numReceptor, vector<float>( 1, 1.0 ) ) ;
1175//     for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1176//       tcal[ir][0] = tcalScalar[ir][0].get() ;
1177//     return tcal ;
1178//   }
1179  else {
1180    return defaultTcal ;
1181  }
1182}
1183
1184void ASDMReader::getTcalAndTsys( unsigned int idx,
1185                                 vector< vector<float> > &tcal,
1186                                 vector< vector<float> > &tsys )
1187{
1188  prepareData( idx ) ;
1189  getTcalAndTsys( tcal, tsys ) ;
1190}
1191
1192void ASDMReader::getTcalAndTsys( vector< vector<float> > &tcal,
1193                                 vector< vector<float> > &tsys )
1194{
1195  String funcName = "getTcalAndTsys" ;
1196
1197  vector< vector<float> > defaultT( 1, vector<float>( 1, 1.0 ) ) ;
1198  SysCalRow *scrow = getSysCalRow() ;
1199  //logsink_->postLocally( LogMessage("scrow = "+String::toString((long)scrow),LogOrigin(className_,funcName,WHERE)) ) ;
1200  if ( scrow == 0 ) {
1201    tcal = defaultT ;
1202    tsys = defaultT ;
1203  }
1204  else {
1205    if ( scrow->isTsysSpectrumExists() ) {
1206      vector< vector<Temperature> > tsysSpec = scrow->getTsysSpectrum() ;
1207      unsigned int numReceptor = tsysSpec.size() ;
1208      unsigned int numChan = tsysSpec[0].size() ;
1209      //logsink_->postLocally( LogMessage("TSYS: numReceptor = "+String::toString(numReceptor),LogOrigin(className_,funcName,WHERE)) ) ;
1210      //logsink_->postLocally( LogMessage("TSYS:numChan = "+String::toString(numChan),LogOrigin(className_,funcName,WHERE)) ) ;
1211      tsys.resize( numReceptor ) ;
1212      for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1213        tsys[ir].resize( numChan ) ;
1214        for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1215          tsys[ir][ic] = tsysSpec[ir][ic].get() ;
1216        }
1217      }
1218    }
1219    else {
1220      tsys = defaultT ;
1221    }
1222    if ( scrow->isTcalSpectrumExists() ) {
1223      vector< vector<Temperature> > tcalSpec = scrow->getTcalSpectrum() ;
1224      unsigned int numReceptor = tcalSpec.size() ;
1225      unsigned int numChan = tcalSpec[0].size() ;
1226      //logsink_->postLocally( LogMessage("TCAL: numReceptor = "+String::toString(numReceptor),LogOrigin(className_,funcName,WHERE)) ) ;
1227      //logsink_->postLocally( LogMessage("TCAL: numChan = "+String::toString(numChan),LogOrigin(className_,funcName,WHERE)) ) ;
1228      tcal.resize( numReceptor ) ;
1229      for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1230        tcal[ir].resize( numChan ) ;
1231        for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1232          tcal[ir][ic] = tcalSpec[ir][ic].get() ;
1233        }
1234      }
1235    }
1236    else {
1237      tcal = defaultT ;
1238    }
1239  }
1240}
1241
1242vector<float> ASDMReader::getOpacity( unsigned int idx )
1243{
1244  prepareData( idx ) ;
1245  return getOpacity() ;
1246}
1247
1248vector<float> ASDMReader::getOpacity()
1249{
1250  vector<float> tau(0) ;
1251  CalAtmosphereTable &atmtab = asdm_->getCalAtmosphere() ;
1252  unsigned int nrow = atmtab.size() ;
1253  if ( nrow > 0 ) {
1254    //int feedid = vmsData_->v_feedId1[index] ;
1255    //BasebandName bbname = spwrow->getBasebandName() ;
1256    //FeedRow *frow = asdm_->getFeed().getRowByKey( atag, spwtag, tint, feedid ) ;
1257    //int nfeed = frow->getNumReceptor() ;
1258    //vector<ReceiverRow *> rrows = frow->getReceivers() ;
1259    vector<CalAtmosphereRow *> atmrows = atmtab.get() ;
1260    //ReceiverBand rb = rrows[0]->getFrequencyBand() ;
1261    int row0 = -1 ;
1262    //double eps = DBL_MAX ;
1263    double eps = 1.0e10 ;
1264    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1265      CalAtmosphereRow *atmrow = atmrows[irow] ;
1266      if ( casacore::String(atmrow->getAntennaName()) != antennaName_
1267           //|| atmrow->getReceiverBand() != rb
1268           //|| atmrow->getBasebandName() != bbname
1269           || atmrow->getCalDataUsingCalDataId()->getCalType() != CalTypeMod::CAL_ATMOSPHERE )
1270        continue ;
1271      else {
1272        double dt = timeInterval_.getStart().getMJD() - atmrow->getEndValidTime().getMJD() ;
1273        if ( dt >= 0 && dt < eps ) {
1274          eps = dt ;
1275          row0 = (int)irow ;
1276        }
1277      }
1278    }
1279    if ( row0 != -1 ) {
1280      CalAtmosphereRow *atmrow = atmrows[row0] ;
1281      tau = atmrow->getTau() ;
1282    }
1283    else {
1284      tau.resize( 1 ) ;
1285      tau[0] = 0.0 ;
1286    }
1287  }
1288  else {
1289    tau.resize( 1 ) ;
1290    tau[0] = 0.0 ;
1291  }
1292  return tau ;
1293}
1294
1295void ASDMReader::getWeatherInfo( unsigned int idx,
1296                                 float &temperature,
1297                                 float &pressure,
1298                                 float &humidity,
1299                                 float &windspeed,
1300                                 float &windaz )
1301{
1302  prepareData( idx ) ;
1303  getWeatherInfo( temperature, pressure, humidity, windspeed, windaz ) ;
1304}
1305
1306void ASDMReader::getWeatherInfo( float &temperature,
1307                                 float &pressure,
1308                                 float &humidity,
1309                                 float &windspeed,
1310                                 float &windaz )
1311{
1312  casacore::String funcName = "getWeatherInfo" ;
1313
1314  temperature = 0.0 ;
1315  pressure = 0.0 ;
1316  humidity = 0.0 ;
1317  windspeed = 0.0 ;
1318  windaz = 0.0 ;
1319
1320  //logsink_->postLocally( LogMessage("weatherStationId_ = "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1321
1322  WeatherTable &wtab = asdm_->getWeather() ;
1323  if ( wtab.size() == 0 || weatherStationId_ == -1 )
1324    return ;
1325
1326  Tag sttag( (unsigned int)weatherStationId_, TagType::Station ) ;
1327  //WeatherRow *wrow = wtab.getRowByKey( sttag, tint ) ;
1328  vector<WeatherRow *> *wrows = wtab.getByContext( sttag ) ;
1329  WeatherRow *wrow = (*wrows)[0] ;
1330  unsigned int nrow = wrows->size() ;
1331  //logsink_->postLocally( LogMessage("There are "+String::toString(nrow)+" rows for given context: stationId "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1332  ArrayTime startTime = getMidTime( timeInterval_ ) ;
1333  if ( startTime < (*wrows)[0]->getTimeInterval().getStart() ) {
1334    temperature = (*wrows)[0]->getTemperature().get() ;
1335    pressure = (*wrows)[0]->getPressure().get() ;
1336    humidity = (*wrows)[0]->getRelHumidity().get() ;
1337    windspeed = (*wrows)[0]->getWindSpeed().get() ;
1338    windaz = (*wrows)[0]->getWindDirection().get() ;
1339  }
1340  else if ( startTime > getEndTime( (*wrows)[nrow-1]->getTimeInterval() ) ) {
1341    temperature = (*wrows)[nrow-1]->getTemperature().get() ;
1342    pressure = (*wrows)[nrow-1]->getPressure().get() ;
1343    humidity = (*wrows)[nrow-1]->getRelHumidity().get() ;
1344    windspeed = (*wrows)[nrow-1]->getWindSpeed().get() ;
1345    windaz = (*wrows)[nrow-1]->getWindDirection().get() ;
1346  }
1347  else {
1348    for ( unsigned int irow = 1 ; irow < wrows->size() ; irow++ ) {
1349      wrow = (*wrows)[irow-1] ;
1350      ArrayTime tStart = wrow->getTimeInterval().getStart() ;
1351      ArrayTime tEnd = (*wrows)[irow]->getTimeInterval().getStart() ;
1352      if ( startTime >= tStart && startTime <= tEnd ) {
1353        temperature = wrow->getTemperature().get() ;
1354        pressure = wrow->getPressure().get() ;
1355        humidity = wrow->getRelHumidity().get() ;
1356        windspeed = wrow->getWindSpeed().get() ;
1357        windaz = wrow->getWindDirection().get() ;
1358        break ;
1359      }
1360    }
1361  }
1362
1363  // Pa -> hPa
1364  pressure /= 100.0 ;
1365
1366  return ;
1367}
1368
1369void ASDMReader::processStation()
1370{
1371  antennaPad_.resize(0) ;
1372  weatherStation_.resize(0) ;
1373  StationTable &stab = asdm_->getStation() ;
1374  vector<StationRow *> srows = stab.get() ;
1375  for ( unsigned int irow = 0 ; irow < srows.size() ; irow++ ) {
1376    StationTypeMod::StationType stype = srows[irow]->getType() ;
1377    Tag stag = srows[irow]->getStationId() ;
1378    if ( stype == StationTypeMod::ANTENNA_PAD )
1379      antennaPad_.push_back( stag ) ;
1380    else if ( stype == StationTypeMod::WEATHER_STATION )
1381      weatherStation_.push_back( stag ) ;
1382  }
1383
1384   weatherStationId_ = getClosestWeatherStation() ;
1385}
1386
1387int ASDMReader::getClosestWeatherStation()
1388{
1389  if ( weatherStation_.size() == 0 )
1390    return -1 ;
1391
1392  vector<double> apos( 3 ) ;
1393  StationTable &stab = asdm_->getStation() ;
1394  apos[0] = antennaPosition_[0].getValue( Unit("m") ) ;
1395  apos[1] = antennaPosition_[1].getValue( Unit("m") ) ;
1396  apos[2] = antennaPosition_[2].getValue( Unit("m") ) ;
1397 
1398  //double eps = DBL_MAX ;
1399  double eps = 1.0e10 ;
1400  int retval = -1 ;
1401  for ( unsigned int ir = 0 ; ir < weatherStation_.size() ; ir++ ) {
1402    StationRow *srow = stab.getRowByKey( weatherStation_[ir] ) ;
1403    vector<Length> wpos = srow->getPosition() ;
1404    double dist = (apos[0]-wpos[0].get())*(apos[0]-wpos[0].get())
1405      + (apos[1]-wpos[1].get())*(apos[1]-wpos[1].get())
1406      + (apos[2]-wpos[2].get())*(apos[2]-wpos[2].get()) ;
1407    if ( dist < eps ) {
1408      retval = (int)(weatherStation_[ir].getTagValue()) ;
1409    }
1410  }
1411
1412  return retval ;
1413}
1414
1415void ASDMReader::getPointingInfo( unsigned int idx,
1416                                  vector<double> &dir,
1417                                  double &az,
1418                                  double &el,
1419                                  vector<double> &srate )
1420{
1421  prepareData( idx ) ;
1422  getPointingInfo( dir, az, el, srate ) ;
1423}
1424
1425void ASDMReader::getPointingInfo( vector<double> &dir,
1426                                  double &az,
1427                                  double &el,
1428                                  vector<double> &srate )
1429{
1430  String funcName = "getPointingInfo" ;
1431
1432  dir.resize(0) ;
1433  az = -1.0 ;
1434  el = -1.0 ;
1435  srate.resize(0) ;
1436
1437  //Tag atag( antennaId_, TagType::Antenna ) ;
1438  //unsigned int index = dataIdList_[idx] ;
1439  //vector<PointingRow *> *prows = asdm_->getPointing().getByContext( atag ) ;
1440  vector<PointingRow *> *prows = asdm_->getPointing().getByContext( antennaTag_ ) ;
1441
1442  if ( prows == 0 )
1443    return ;
1444
1445  PointingRow *prow ;
1446  PointingRow *qrow ;
1447  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1448  //double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1449  //ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1450
1451  unsigned int nrow = prows->size() ;
1452  //logsink_->postLocally( LogMessage("There are " << nrow << " rows for given context: antennaId "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
1453
1454//   for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1455//     prow = (*prows)[irow] ;
1456//     ArrayTimeInterval ati = prow->getTimeInterval() ;
1457//     ArrayTime pst = ati.getStart() ;
1458//     ArrayTime pet( ati.getStartInMJD()+ati.getDurationInDays() ) ;
1459//     logsink_->postLocally( LogMessage("start: "+pst.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1460//     logsink_->postLocally( LogMessage("end: "+pet.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1461//   }
1462 
1463  srate.resize(2) ;
1464  srate[0] = 0.0 ;
1465  srate[1] = 0.0 ;
1466  az = 0.0 ;
1467  el = 0.0 ;
1468  //double tcen = 0.0 ;
1469  //double tcen = getMidTime( tint ).getMJD() ;
1470  double tcen = getMidTime( timeInterval_ ).getMJD() ;
1471
1472  //
1473  // shape of pointingDirection is (numSample,2) if usePolynomial = False, while it is
1474  // (numTerm,2) if usePolynomial = True.
1475  //
1476  // In the former case, typical sampled time interval is 48msec, which is very small
1477  // compared with typical integration time (~a few sec).
1478  // Scan rate for this case is always [0,0] (or get slope?).
1479  //
1480  // In the later case, scan rate is (pointingDirection[1][0],pointingDirection[1][1])
1481  //
1482  // PointingDirection seems to store AZEL
1483  //
1484  ArrayTimeInterval pTime0 = (*prows)[0]->getTimeInterval() ;
1485  ArrayTimeInterval pTime1 = (*prows)[nrow-1]->getTimeInterval() ;
1486  //if ( tint.getStartInMJD()+tint.getDurationInDays() < pTime0.getStartInMJD() ) {
1487  //if ( getEndTime( tint ) < getStartTime( pTime0 ) ) {
1488  if ( getEndTime( timeInterval_ ) < getStartTime( pTime0 ) ) {
1489    logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint < ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1490    prow = (*prows)[0] ;
1491    vector< vector<double> > dirA = pointingDir( prow ) ;
1492    az = dirA[0][0] ;
1493    el = dirA[0][1] ;
1494    if ( prow->getUsePolynomials() && prow->getNumTerm() > 1 ) {
1495      srate[0] = dirA[1][0] ;
1496      srate[1] = dirA[1][1] ;
1497    }     
1498  }
1499  //else if ( tint.getStartInMJD() > pTime1.getStartInMJD()+pTime1.getDurationInDays() ) {
1500  //else if ( getStartTime( tint ) > getEndTime( pTime1 ) ) {
1501  else if ( getStartTime( timeInterval_ ) > getEndTime( pTime1 ) ) {
1502    logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint > ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1503    prow = (*prows)[nrow-1] ;
1504    int numSample = prow->getNumSample() ;
1505    vector< vector<double> > dirA = pointingDir( prow ) ;
1506    if ( prow->getUsePolynomials() ) {
1507      az = dirA[0][0] ;
1508      el = dirA[0][1] ;
1509      if ( prow->getNumTerm() > 1 ) {
1510        srate[0] = dirA[1][0] ;
1511        srate[1] = dirA[1][1] ;
1512      }
1513    }
1514    else {
1515      az = dirA[numSample-1][0] ;
1516      el = dirA[numSample-1][1] ;
1517    }
1518  }
1519  else {
1520    //ArrayTime startTime = tint.getStart() ;
1521    //ArrayTime endTime = getEndTime( tint ) ;
1522    ArrayTime startTime = timeInterval_.getStart() ;
1523    ArrayTime endTime = getEndTime( timeInterval_ ) ;
1524    int row0 = -1 ;
1525    int row1 = -1 ;
1526    int row2 = -1 ;
1527    //double dt0 = getMidTime( tint ).getMJD() ;
1528    double dt0 = getMidTime( timeInterval_ ).getMJD() ;
1529    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1530      prow = (*prows)[irow] ;
1531      //double dt = getMidTime( tint ).getMJD() - getMidTime( prow->getTimeInterval() ).getMJD() ;
1532      double dt = getMidTime( timeInterval_ ).getMJD() - getMidTime( prow->getTimeInterval() ).getMJD() ;
1533      if ( dt > 0 && dt < dt0 ) {
1534        dt0 = dt ;
1535        row2 = irow ;
1536      }
1537      if ( prow->getTimeInterval().contains( startTime ) )
1538        row0 = irow ;
1539      else if ( prow->getTimeInterval().contains( endTime ) )
1540        row1 = irow ;
1541      if ( row0 != -1 && row1 != -1 )
1542        break ;
1543    }
1544    //logsink_->postLocally( LogMessage("row0 = "+String::toString(row0)+", row1 = "+String::toString(row1)+", row2 = "+String::toString(row2),LogOrigin(className_,funcName,WHERE)) ) ;
1545    if ( row0 == -1 && row1 == -1 ) {
1546      prow = (*prows)[row2] ;
1547      qrow = (*prows)[row2+1] ;
1548      double t0 = getEndTime( prow->getTimeInterval() ).getMJD() ;
1549      double t1 = qrow->getTimeInterval().getStartInMJD() ;
1550      vector< vector<double> > dirA = pointingDir( prow ) ;
1551      vector< vector<double> > dirB = pointingDir( qrow ) ;
1552      double da0 = dirA[0][0] ;
1553      double db0 = dirB[0][0] ;
1554      double da1 = dirA[0][1] ;
1555      double db1 = dirB[0][1] ;
1556      if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1557        double dt = ( tcen - t0 ) / ( t1 - t0 ) ;
1558        az = da0 + ( db0 - da0 ) * dt ;
1559        el = da1 + ( db1 - da1 ) * dt ;
1560        if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1561          double ra0 = dirA[1][0] ;
1562          double rb0 = dirB[1][0] ;
1563          double ra1 = dirA[1][1] ;
1564          double rb1 = dirB[1][1] ;
1565          srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1566          srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1567        }
1568      }
1569      else if ( !(qrow->getUsePolynomials()) ) {
1570        double dt = ( tcen - t0 ) / ( t1 - t0 ) ;
1571        az = da0 + ( db0 - da0 ) * dt ;
1572        el = da1 + ( db1 - da1 ) * dt ;
1573      }
1574      //else {
1575      // nothing to do
1576      //}
1577    }
1578    else {
1579      int ndir = 0 ;
1580      if ( row0 == -1 ) {
1581        row0 = row1 ;
1582      }
1583      else if ( row1 == -1 ) {
1584        row1 = row0 ;
1585      }
1586      prow = (*prows)[row0] ;
1587      qrow = (*prows)[row1] ;
1588      if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1589        //logsink_->postLocally( LogMessage("usePolynomial = True",LogOrigin(className_,funcName,WHERE)) ) ;
1590        if ( row0 == row1 ) {
1591          prow = (*prows)[row0] ;
1592          vector< vector<double> > dirA = pointingDir( prow ) ;
1593          az = dirA[0][0] ;
1594          el = dirA[0][1] ;
1595          if ( prow->getNumTerm() > 1 ) {
1596            srate[0] = dirA[1][0] ;
1597            srate[1] = dirA[1][1] ;
1598          }
1599        }
1600        else {
1601          prow = (*prows)[row0] ;
1602          qrow = (*prows)[row1] ;
1603          vector< vector<double> > dirA = pointingDir( prow ) ;
1604          vector< vector<double> > dirB = pointingDir( qrow ) ;
1605          double t0 = getMidTime( prow->getTimeInterval() ).getMJD() ;
1606          double t1 = getMidTime( qrow->getTimeInterval() ).getMJD() ;
1607          double dt = ( tcen - t0 ) / ( t1 - t0 ) ;
1608          double da0 = dirA[0][0] ;
1609          double db0 = dirB[0][0] ;
1610          double da1 = dirA[0][1] ;
1611          double db1 = dirB[0][1] ;
1612          az = da0 + ( db0 - da0 ) * dt ;
1613          el = da1 + ( db1 - db0 ) * dt ;
1614          if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1615            double ra0 = dirA[1][0] ;
1616            double rb0 = dirB[1][0] ;
1617            double ra1 = dirA[1][1] ;
1618            double rb1 = dirB[1][1] ;
1619            srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1620            srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1621          }
1622        }
1623      }
1624      else if ( prow->getUsePolynomials() == qrow->getUsePolynomials() ) {
1625        //logsink_->postLocally( LogMessage("numSample == numTerm",LogOrigin(className_,funcName,WHERE)) ) ;
1626        tcen = 0.0 ;
1627        for ( int irow = row0 ; irow <= row1 ; irow++ ) {
1628          prow = (*prows)[irow] ;
1629          int numSample = prow->getNumSample() ;
1630          //logsink_->postLocally( LogMessage("numSample = "+String::toString(numSample),LogOrigin(className_,funcName,WHERE)) ) ;
1631          vector< vector<double> > dirA = pointingDir( prow ) ;
1632          if ( prow->isSampledTimeIntervalExists() ) {
1633            //logsink_->postLocally( LogMessage("sampledTimeIntervalExists",LogOrigin(className_,funcName,WHERE)) ) ;
1634            vector<ArrayTimeInterval> stime = prow->getSampledTimeInterval() ;
1635            for ( int isam = 0 ; isam < numSample ; isam++ ) {
1636              //if ( tint.overlaps( stime[isam] ) ) {
1637              //if ( tint.contains( stime[isam] ) ) {
1638              if ( timeInterval_.contains( stime[isam] ) ) {
1639                az += dirA[isam][0] ;
1640                el += dirA[isam][1] ;
1641                tcen += getMidTime( stime[isam] ).getMJD() ;
1642                ndir++ ;
1643              }
1644            }
1645          }
1646          else {
1647            double sampleStart = prow->getTimeInterval().getStartInMJD() ;
1648            double sampleInterval = prow->getTimeInterval().getDurationInDays() / (double)numSample ;
1649            //logsink_->postLocally( LogMessage("sampleStart = "+String::toString(sampleStart),LogOrigin(className_,funcName,WHERE)) )
1650            //logsink_->postLocally( LogMessage("sampleInterval = "+String::toString(sampleInterval),LogOrigin(className_,funcName,WHERE)) ) ;
1651            //logsink_->postLocally( LogMessage("tint = "+tint.toString(),LogOrigin(className_,funcName,WHERE)) ) ;
1652            for ( int isam = 0 ; isam < numSample ; isam++ ) {
1653              ArrayTimeInterval stime( sampleStart+isam*sampleInterval, sampleInterval ) ;
1654              //if ( tint.overlaps( stime ) ) {
1655              //if ( tint.contains( stime ) ) {
1656              if ( timeInterval_.contains( stime ) ) {
1657                az += dirA[isam][0] ;
1658                el += dirA[isam][1] ;
1659                tcen += getMidTime( stime ).getMJD() ;
1660                ndir++ ;
1661              }
1662            }
1663          }
1664        }
1665        if ( ndir > 1 ) {
1666          az /= (double)ndir ;
1667          el /= (double)ndir ;
1668          tcen /= (double)ndir ;
1669        }
1670      }
1671      //else {
1672      // nothing to do
1673      //}
1674    }
1675   
1676    toJ2000( dir, az, el, tcen ) ;
1677  }
1678
1679  return ;
1680}
1681
1682ArrayTime ASDMReader::getMidTime( const ArrayTimeInterval &t )
1683{
1684  return ArrayTime( t.getStartInMJD() + 0.5 * t.getDurationInDays() ) ;
1685}
1686
1687ArrayTime ASDMReader::getEndTime( const ArrayTimeInterval &t )
1688{
1689  return ArrayTime( t.getStartInMJD() + t.getDurationInDays() ) ;
1690}
1691
1692ArrayTime ASDMReader::getStartTime( const ArrayTimeInterval &t )
1693{
1694  return ArrayTime( t.getStartInMJD() ) ;
1695}
1696
1697void ASDMReader::toJ2000( vector<double> &dir,
1698                          double &az,
1699                          double &el,
1700                          double &mjd )
1701{
1702  String funcName = "toJ2000" ;
1703
1704  String ref = "AZELGEO" ;
1705  vector<double> azel( 2 ) ;
1706  azel[0] = az ;
1707  azel[1] = el ;
1708  dir = toJ2000( azel, ref, mjd ) ;
1709}
1710
1711vector<double> ASDMReader::toJ2000( vector<double> &dir,
1712                                    casacore::String &dirref,
1713                                    double &mjd )
1714{
1715  casacore::String funcName = "toJ2000" ;
1716
1717  vector<double> newd( dir ) ;
1718  if ( dirref != "J2000" ) {
1719    me_ = casacore::MEpoch( casacore::Quantity( (casacore::Double)mjd, "d" ), casacore::MEpoch::UTC ) ;
1720    mf_.set( me_ ) ;
1721    casacore::MDirection::Types dirtype ;
1722    casacore::Bool b = casacore::MDirection::getType( dirtype, dirref ) ;
1723    if ( b ) {
1724      casacore::Vector<casacore::Double> cdir = toj2000_( dir ).getAngle( "rad" ).getValue() ;
1725      //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
1726      newd[0] = (double)(cdir[0]) ;
1727      newd[1] = (double)(cdir[1]) ;
1728    }
1729  }
1730  return newd ;
1731}
1732
1733void ASDMReader::setLogger( CountedPtr<LogSinkInterface> &logsink )
1734{
1735  logsink_ = logsink ;
1736}
1737
1738string ASDMReader::getFrame()
1739{
1740  casacore::String funcName = "getFrame" ;
1741 
1742  // default is TOPO
1743  string frame = "TOPO" ;
1744
1745  SpectralWindowTable &spwtab = asdm_->getSpectralWindow() ;
1746  vector<SpectralWindowRow *> rows = spwtab.get() ;
1747  vector<FrequencyReferenceCodeMod::FrequencyReferenceCode> measFreqRef( rows.size() ) ;
1748  int nref = 0 ;
1749  for ( unsigned int irow = 0 ; irow < rows.size() ; irow++ ) {
1750    int nchan = rows[irow]->getNumChan() ;
1751    if ( nchan != 4 ) {
1752      if ( rows[irow]->isMeasFreqRefExists() ) {
1753        measFreqRef[nref] = rows[irow]->getMeasFreqRef() ;
1754        nref++ ;
1755      }
1756    }
1757  }
1758  if ( nref != 0 ) {
1759    frame = CFrequencyReferenceCode::toString( measFreqRef[0] ) ;
1760  }
1761 
1762  //logsink_->postLocally( LogMessage("frame = "+String::toString(frame),LogOrigin(className_,funcName,WHERE)) ) ;
1763
1764  return frame ;
1765}
1766
1767int ASDMReader::getNumIFs()
1768{
1769  casacore::String funcName = "getNumIFs" ;
1770
1771  int nif = 0 ;
1772  vector<SpectralWindowRow *> rows = asdm_->getSpectralWindow().get() ;
1773  unsigned int nrow = rows.size() ;
1774  // check if all rows have freqGroup attribute
1775  bool freqGroupExists = true ;
1776  bool countedWvr = false ;
1777  for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1778    freqGroupExists &= rows[irow]->isFreqGroupExists() ;
1779    if ( rows[irow]->getNumChan() == 4 ) {
1780      if ( !countedWvr ) {
1781        countedWvr = true ;
1782        nif++ ;
1783      }
1784    }
1785    else {
1786      nif++ ;
1787    }
1788  }
1789 
1790  if ( freqGroupExists ) {
1791    vector<int> freqGroup(0) ;
1792    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1793      int fg = rows[irow]->getFreqGroup() ;
1794      if ( (int)count( freqGroup.begin(), freqGroup.end(), fg ) == 0 ) {
1795        freqGroup.push_back( fg ) ;
1796      }
1797    }
1798    nif = freqGroup.size() ;
1799  }
1800
1801  //logsink_->postLocally( LogMessage("nif = "+String::toString(nif),LogOrigin(className_,funcName,WHERE)) ) ;
1802
1803  return nif ;
1804}
1805
1806SysCalRow *ASDMReader::getSysCalRow( unsigned int idx )
1807{
1808  prepareData( idx ) ;
1809  return getSysCalRow() ;
1810}
1811
1812SysCalRow *ASDMReader::getSysCalRow()
1813{
1814  String funcName = "getSysCalRow" ;
1815
1816  SysCalRow *row = 0 ;
1817  int feedid = vmsData_->v_feedId1[dataIndex_] ;
1818  SysCalTable &sctab = asdm_->getSysCal() ;
1819  vector<SysCalRow *> *rows = sctab.getByContext( antennaTag_, specWinTag_, feedid ) ;
1820  //if ( nrow == 0 ) {
1821  if ( rows == 0 ) {
1822    //logsink_->postLocally( LogMessage("no rows in SysCal table",LogOrigin(className_,funcName,WHERE)) ) ;
1823    row = 0 ;
1824  }
1825  else {
1826    unsigned int nrow = rows->size() ;
1827    //logsink_->postLocally( LogMessage("nrow = "+String::toString(nrow),LogOrigin(className_,funcName,WHERE)) ) ;
1828    int scindex = -1 ;
1829    if ( nrow == 1 ) {
1830      scindex = 0 ;
1831    }
1832    else if ( getEndTime( timeInterval_ ) <= getStartTime( (*rows)[0]->getTimeInterval() ) )
1833      scindex = 0 ;
1834    else {
1835      for ( unsigned int irow = 0 ; irow < nrow-1 ; irow++ ) {
1836        ArrayTime t = getMidTime( timeInterval_ ) ;
1837        if ( t > getStartTime( (*rows)[irow]->getTimeInterval() )
1838             && t <= getStartTime( (*rows)[irow+1]->getTimeInterval() ) ) {
1839          scindex = irow ;
1840          break ;
1841        }
1842      }
1843      if ( scindex == -1 )
1844        scindex = nrow-1 ;
1845    }
1846    //logsink_->postLocally( LogMessage("irow = "+String::toString(scindex),LogOrigin(className_,funcName,WHERE)) ) ;
1847    row = (*rows)[scindex] ;
1848  }
1849  return row ;
1850}
1851
1852double ASDMReader::limitedAngle( double angle )
1853{
1854  if ( angle > C::pi )
1855    while ( angle > C::pi ) angle -= C::_2pi ;
1856  else if ( angle < -C::pi )
1857    while ( angle < -C::pi ) angle += C::_2pi ;
1858  return angle ;
1859}
1860
1861vector< vector<double> > ASDMReader::pointingDir( PointingRow *row )
1862{
1863  //String funcName = "pointingDir" ;
1864  vector< vector<Angle> > aTar = row->getTarget() ;
1865  vector< vector<Angle> > aOff = row->getOffset() ;
1866  vector< vector<Angle> > aDir = row->getPointingDirection() ;
1867  vector< vector<Angle> > aEnc = row->getEncoder() ;
1868  unsigned int n = aTar.size() ;
1869  vector< vector<double> > dir( n ) ;
1870  double factor = 1.0 / cos( aTar[0][1].get() ) ;
1871  for ( unsigned int i = 0 ; i < n ; i++ ) {
1872    dir[i].resize( 2 ) ;
1873    /**
1874     * This is approximate way to add offset taking tracking error
1875     * into account
1876     *
1877     * az = dir[0][0] = target[0][0] + offset[0][0] / cos(el)
1878     *                 + encorder[0][0] - direction[0][0]
1879     * el = dir[0][1] = target[0][1] + offset[0][1]
1880     *                 + encorder[0][1] - direction[0][1]
1881     **/
1882    dir[i][0] = aTar[i][0].get() + factor * aOff[i][0].get()
1883               + aEnc[i][0].get() - aDir[i][0].get() ;
1884    dir[i][1] = aTar[i][1].get() + aOff[i][1].get()
1885               + aEnc[i][1].get() - aDir[i][1].get() ;
1886    //logsink_->postLocally( LogMessage("tracking offset: ["+String::toString(aEnc[i][0].get()-aDir[i][0].get())+","+String::toString(aEnc[i][0]-aDir[i][0].get())+"]",LogOrigin(className_,funcName,WHERE)) ) ;
1887  }
1888  return dir ;
1889}
1890
1891void ASDMReader::initConvert()
1892{
1893  toj2000_ = MDirection::Convert(MDirection::AZELGEO,
1894                                 MDirection::Ref(MDirection::J2000, mf_));
1895}
Note: See TracBrowser for help on using the repository browser.