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

Last change on this file since 2299 was 2299, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: Yes CAS-1913

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...

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