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

Last change on this file since 2471 was 2471, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: 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...

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