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

Last change on this file since 2273 was 2273, 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...

Support latest ASDM definition in asdm2ASAP.
Now asdm2ASAP is linked to libalma_v3.so to be able to read latest ASDM.

The importasdm task invokes asdm2ASAP when singledish is True.

Older version of asdm2ASAP is renamec as oldasdm2ASAP.
The oldasdm2ASAP and related classes are built by linking to libalma.so.
It is installed but not used by any tasks.


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 = subrow->getSubscanMode() ;
894  if ( scanIntent == OBSERVE_TARGET ) {
895    // on sky scan
896    if ( swmode == NO_SWITCHING || swmode == POSITION_SWITCHING ) {
897      // position switching
898      // tentatively set NO_SWITCHING = POSITION_SWITCHING
899      if ( subIntent == ON_SOURCE ) {
900        srctype = SrcType::PSON ;
901      }
902      else if ( subIntent == OFF_SOURCE ) {
903        srctype = SrcType::PSOFF ;
904      }
905    }
906    else if ( swmode == FREQUENCY_SWITCHING ) {
907      // frequency switching
908      if ( subIntent == ON_SOURCE ) {
909        srctype = SrcType::FSON ;
910      }
911      else if ( subIntent == OFF_SOURCE ) {
912        srctype = SrcType::FSOFF ;
913      }
914    }
915    else if ( swmode == NUTATOR_SWITCHING ) {
916      // nutator switching
917      if ( subIntent == ON_SOURCE ) {
918        srctype = SrcType::PSON ;
919      }
920      else if ( subIntent == OFF_SOURCE ) {
921        srctype = SrcType::PSOFF ;
922      }
923    }
924    else {
925      // other switching mode
926      if ( subIntent == ON_SOURCE ) {
927        srctype = SrcType::SIG ;
928      }
929      else if ( subIntent == OFF_SOURCE ) {
930        srctype = SrcType::REF ;
931      }
932    }
933  }
934  else if ( scanIntent == CALIBRATE_ATMOSPHERE ) {
935    if ( swmode == NO_SWITCHING || swmode == POSITION_SWITCHING ) {
936      // position switching
937      // tentatively set NO_SWITCHING = POSITION_SWITCHING
938      if ( subIntent == ON_SOURCE ) {
939        srctype = SrcType::PONCAL ;
940      }
941      else if ( subIntent == OFF_SOURCE ) {
942        srctype = SrcType::POFFCAL ;
943      }
944    }
945    else if ( swmode == FREQUENCY_SWITCHING ) {
946      // frequency switching
947      if ( subIntent == ON_SOURCE ) {
948        srctype = SrcType::FONCAL ;
949      }
950      else if ( subIntent == OFF_SOURCE ) {
951        srctype = SrcType::FOFFCAL ;
952      }
953    }
954    else if ( swmode == NUTATOR_SWITCHING ) {
955      // nutator switching
956      if ( subIntent == ON_SOURCE ) {
957        srctype = SrcType::PONCAL ;
958      }
959      else if ( subIntent == OFF_SOURCE ) {
960        srctype = SrcType::POFFCAL ;
961      }
962    }
963    else {
964      // other switching mode
965      if ( subIntent == ON_SOURCE ) {
966        srctype = SrcType::CAL ;
967      }
968      else if ( subIntent == OFF_SOURCE ) {
969        srctype = SrcType::CAL ;
970      }
971    }
972  }
973  else {
974    // off sky (e.g. calibrator device) scan
975    if ( subIntent == ON_SOURCE ) {
976      srctype = SrcType::SIG ;
977    }
978    else if ( subIntent == OFF_SOURCE ) {
979      srctype = SrcType::REF ;
980    }
981    else if ( subIntent == HOT ) {
982      srctype = SrcType::HOT ;
983    }
984    else if ( subIntent == AMBIENT ) {
985      srctype = SrcType::SKY ;
986    }
987    else {
988      srctype = SrcType::CAL ;
989    }
990  }
991
992  return srctype ;
993}
994
995unsigned int ASDMReader::getSubscanNo( unsigned int idx )
996{
997  //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)) ) ;
998  return vmsData_->v_msState[dataIdList_[idx]].subscanNum ;
999}
1000
1001vector<double> ASDMReader::getSourceDirection( unsigned int idx )
1002{
1003  vector<double> dir( 2, 0.0 ) ;
1004  unsigned int index = dataIdList_[idx] ;
1005  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1006  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1007  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1008  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1009  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1010  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1011  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1012  string srcname ;
1013  if ( frow->isSourceIdExists() ) {
1014    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1015    int sid = frow->getSourceId() ;
1016    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1017    vector<Angle> srcdir = srow->getDirection() ;
1018    dir[0] = limitedAngle( srcdir[0].get() ) ;
1019    dir[1] = limitedAngle( srcdir[1].get() ) ;
1020    if ( srow->isDirectionCodeExists() ) {
1021      DirectionReferenceCode dircode = srow->getDirectionCode() ;
1022      //logsink_->postLocally( LogMessage("dircode="+CDirectionReferenceCode::toString(dircode),LogOrigin(className_,funcName,WHERE)) ) ;
1023      if ( dircode != J2000 ) {
1024        // if not J2000, need direction conversion
1025        string ref = CDirectionReferenceCode::toString( dircode ) ;
1026        double mjd = vmsData_->v_time[index] * s2d ;
1027        Tag atag( antennaId_, TagType::Antenna ) ;
1028        AntennaRow *arow = asdm_->getAntenna().getRowByKey( atag ) ;
1029        StationRow *srow = arow->getStationUsingStationId() ;
1030        vector<Length> antposL = srow->getPosition() ;
1031        casa::Vector<casa::Double> antpos( 3 ) ;
1032        for ( int i = 0 ; i < 3 ; i++ )
1033          antpos[i] = antposL[i].get() ;
1034        dir = toJ2000( dir, ref, mjd, antpos ) ;
1035      }
1036    }
1037  }
1038  return dir ;
1039}
1040
1041void ASDMReader::getSourceDirection( unsigned int idx,
1042                                     vector<double> &dir,
1043                                     string &ref )
1044{
1045  dir.resize( 2 ) ;
1046  dir[0] = 0.0 ;
1047  dir[1] = 0.0 ;
1048  ref = "J2000" ;
1049  unsigned int index = dataIdList_[idx] ;
1050  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1051  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1052  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1053  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1054  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1055  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1056  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1057  string srcname ;
1058  if ( frow->isSourceIdExists() ) {
1059    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1060    int sid = frow->getSourceId() ;
1061    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1062    vector<Angle> srcdir = srow->getDirection() ;
1063    if ( srow->isDirectionCodeExists() ) {
1064      ref = CDirectionReferenceCode::toString( srow->getDirectionCode() ) ;
1065    }
1066    dir[0] = limitedAngle( srcdir[0].get() ) ;
1067    dir[1] = limitedAngle( srcdir[1].get() ) ;
1068  }
1069}
1070
1071void ASDMReader::getSourceDirection( vector<double> &dir, string &ref )
1072{
1073  dir.resize( 2 ) ;
1074  ref = "J2000" ; // default is J2000
1075  SourceTable &tab = asdm_->getSource() ;
1076  SourceRow *row = tab.get()[0] ;
1077  vector<Angle> dirA = row->getDirection() ;
1078  dir[0] = dirA[0].get() ;
1079  dir[1] = dirA[1].get() ;
1080  if ( row->isDirectionCodeExists() ) {
1081    ref = CDirectionReferenceCode::toString( row->getDirectionCode() ) ;
1082  }
1083}
1084
1085vector<double> ASDMReader::getSourceProperMotion( unsigned int idx )
1086{
1087  vector<double> pm( 2, 0.0 ) ;
1088  unsigned int index = dataIdList_[idx] ;
1089  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1090  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1091  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1092  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1093  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1094  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1095  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1096  string srcname ;
1097  if ( frow->isSourceIdExists() ) {
1098    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1099    int sid = frow->getSourceId() ;
1100    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1101    vector<AngularRate> srcpm = srow->getProperMotion() ;
1102    pm[0] = srcpm[0].get() ;
1103    pm[1] = srcpm[1].get() ;
1104  }
1105  return pm ;
1106}
1107
1108double ASDMReader::getSysVel( unsigned int idx )
1109{
1110  double sysvel = 0.0 ;
1111  unsigned int index = dataIdList_[idx] ;
1112  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1113  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1114  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1115  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1116  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1117  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1118  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1119  string srcname ;
1120  if ( frow->isSourceIdExists() ) {
1121    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1122    int sid = frow->getSourceId() ;
1123    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1124    if ( srow->isSysVelExists() ) {
1125      vector<Speed> sysvelV = srow->getSysVel() ;
1126      if ( sysvelV.size() > 0 )
1127        sysvel = sysvelV[0].get() ;
1128    }
1129  }
1130  return sysvel ;
1131}
1132
1133unsigned int ASDMReader::getFlagRow( unsigned int idx )
1134{
1135  return vmsData_->v_flag[dataIdList_[idx]] ;
1136}
1137
1138vector<unsigned int> ASDMReader::getDataShape( unsigned int idx )
1139{
1140  return vmsData_->vv_dataShape[dataIdList_[idx]] ;
1141}
1142
1143float * ASDMReader::getSpectrum( unsigned int idx )
1144{
1145  map<AtmPhaseCorrection, float*> data = vmsData_->v_m_data[dataIdList_[idx]] ;
1146  //map<AtmPhaseCorrection, float*>::iterator iter = data.find(AP_UNCORRECTED) ;
1147  map<AtmPhaseCorrection, float*>::iterator iter = data.find(apc_) ;
1148  float *autoCorr = iter->second ;
1149  return autoCorr ;
1150}
1151
1152// bool * ASDMReader::getFlagChannel( unsigned int idx )
1153// {
1154//   return 0 ;
1155// }
1156
1157vector< vector<float> > ASDMReader::getTsys( unsigned int idx )
1158{
1159  vector< vector<float> > defaultTsys( 1, vector<float>( 1, 1.0 ) ) ;
1160  SysCalRow *scrow = getSysCalRow( idx ) ;
1161  if ( scrow != 0 && scrow->isTsysSpectrumExists() ) {
1162    vector< vector<Temperature> > tsysSpec = scrow->getTsysSpectrum() ;
1163    unsigned int numReceptor = tsysSpec.size() ;
1164    unsigned int numChan = tsysSpec[0].size() ;
1165    vector< vector<float> > tsys( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1166    for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1167      for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1168        tsys[ir][ic] = tsysSpec[ir][ic].get() ;
1169      }
1170    }
1171    return tsys ;
1172  }
1173//   else if ( scrow->isTsysExists() ) {
1174//     vector<Temperature> tsysScalar = scrow->getTsys() ;
1175//     unsigned int numReceptor = tsysScalar.size() ;
1176//     vector< vector<float> > tsys( numReceptor ) ;
1177//     for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1178//       tsys[ir] = vector<float>( 1, tsysScalar[ir].get() ) ;
1179//     return tsys ;
1180//   }
1181  else {
1182    return defaultTsys ;
1183  }
1184}
1185
1186vector< vector<float> > ASDMReader::getTcal( unsigned int idx )
1187{
1188  vector< vector<float> > defaultTcal( 1, vector<float>( 1, 1.0 ) ) ;
1189  SysCalRow *scrow = getSysCalRow( idx ) ;
1190  if ( scrow != 0 && scrow->isTcalSpectrumExists() ) {
1191    vector< vector<Temperature> > tcalSpec = scrow->getTcalSpectrum() ;
1192    unsigned int numReceptor = tcalSpec.size() ;
1193    unsigned int numChan = tcalSpec[0].size() ;
1194    vector< vector<float> > tcal( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1195    for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1196      for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1197        tcal[ir][ic] = tcalSpec[ir][ic].get() ;
1198      }
1199    }
1200    return tcal ;
1201  }
1202//   else if ( scrow->isTcalExists() ) {
1203//     vector<Temperature> tcalScalar = scrow->getTcal() ;
1204//     unsigned int numReceptor = tcalScalar.size() ;
1205//     vector< vector<float> > tcal( numReceptor, vector<float>( 1, 1.0 ) ) ;
1206//     for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1207//       tcal[ir][0] = tcalScalar[ir][0].get() ;
1208//     return tcal ;
1209//   }
1210  else {
1211    return defaultTcal ;
1212  }
1213}
1214
1215void ASDMReader::getTcalAndTsys( unsigned int idx,
1216                                 vector< vector<float> > &tcal,
1217                                 vector< vector<float> > &tsys )
1218{
1219  String funcName = "getTcalAndTsys" ;
1220
1221  vector< vector<float> > defaultT( 1, vector<float>( 1, 1.0 ) ) ;
1222  SysCalRow *scrow = getSysCalRow( idx ) ;
1223  if ( scrow == 0 ) {
1224    tcal = defaultT ;
1225    tsys = defaultT ;
1226  }
1227  else {
1228    if ( scrow->isTsysSpectrumExists() ) {
1229      vector< vector<Temperature> > tsysSpec = scrow->getTsysSpectrum() ;
1230      unsigned int numReceptor = tsysSpec.size() ;
1231      unsigned int numChan = tsysSpec[0].size() ;
1232      tsys.resize( numReceptor ) ;
1233      for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1234        tsys[ir].resize( numChan ) ;
1235        for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1236          tsys[ir][ic] = tsysSpec[ir][ic].get() ;
1237        }
1238      }
1239    }
1240    else {
1241      tsys = defaultT ;
1242    }
1243    if ( scrow->isTcalSpectrumExists() ) {
1244      vector< vector<Temperature> > tcalSpec = scrow->getTcalSpectrum() ;
1245      unsigned int numReceptor = tcalSpec.size() ;
1246      unsigned int numChan = tcalSpec[0].size() ;
1247      tcal.resize( numReceptor ) ;
1248      for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1249        tcal[ir].resize( numReceptor ) ;
1250        for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1251          tcal[ir][ic] = tcalSpec[ir][ic].get() ;
1252        }
1253      }
1254    }
1255    else {
1256      tcal = defaultT ;
1257    }
1258  }
1259}
1260
1261vector<float> ASDMReader::getOpacity( unsigned int idx )
1262{
1263  vector<float> tau(0) ;
1264  CalAtmosphereTable &atmtab = asdm_->getCalAtmosphere() ;
1265  unsigned int nrow = atmtab.size() ;
1266  if ( nrow > 0 ) {
1267    unsigned int index = dataIdList_[idx] ;
1268    //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1269    double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1270    ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1271    //int feedid = vmsData_->v_feedId1[index] ;
1272    //Tag atag( antennaId_, TagType::Antenna ) ;
1273    //Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1274    //DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey(ddtag) ;
1275    //Tag spwtag = ddrow->getSpectralWindowId() ;
1276    //SpectralWindowRow *spwrow = ddrow->getSpectralWindowUsingSpectralWindowId() ;
1277    //BasebandName bbname = spwrow->getBasebandName() ;
1278    //FeedRow *frow = asdm_->getFeed().getRowByKey( atag, spwtag, tint, feedid ) ;
1279    //int nfeed = frow->getNumReceptor() ;
1280    //vector<ReceiverRow *> rrows = frow->getReceivers() ;
1281    vector<CalAtmosphereRow *> atmrows = atmtab.get() ;
1282    //ReceiverBand rb = rrows[0]->getFrequencyBand() ;
1283    int row0 = -1 ;
1284    double eps = tint.getStart().getMJD() ;
1285    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1286      CalAtmosphereRow *atmrow = atmrows[irow] ;
1287      if ( casa::String(atmrow->getAntennaName()) != antennaName_
1288           //|| atmrow->getReceiverBand() != rb
1289           //|| atmrow->getBasebandName() != bbname
1290           || atmrow->getCalDataUsingCalDataId()->getCalType() != CAL_ATMOSPHERE )
1291        continue ;
1292      else {
1293        double dt = tint.getStart().getMJD() - atmrow->getEndValidTime().getMJD() ;
1294        if ( dt >= 0 && dt < eps ) {
1295          eps = dt ;
1296          row0 = (int)irow ;
1297        }
1298      }
1299    }
1300    if ( row0 != -1 ) {
1301      CalAtmosphereRow *atmrow = atmrows[row0] ;
1302      tau = atmrow->getTau() ;
1303    }
1304  }
1305  else {
1306    tau.resize( 1 ) ;
1307    tau[0] = 0.0 ;
1308  }
1309  return tau ;
1310}
1311
1312void ASDMReader::getWeatherInfo( unsigned int idx,
1313                                 float &temperature,
1314                                 float &pressure,
1315                                 float &humidity,
1316                                 float &windspeed,
1317                                 float &windaz )
1318{
1319  casa::String funcName = "getWeatherInfo" ;
1320
1321  temperature = 0.0 ;
1322  pressure = 0.0 ;
1323  humidity = 0.0 ;
1324  windspeed = 0.0 ;
1325  windaz = 0.0 ;
1326
1327  //logsink_->postLocally( LogMessage("weatherStationId_ = "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1328
1329  WeatherTable &wtab = asdm_->getWeather() ;
1330  if ( wtab.size() == 0 || weatherStationId_ == -1 )
1331    return ;
1332
1333  unsigned int index = dataIdList_[idx] ;
1334  //Tag anttag( antennaId_, TagType::Antenna ) ;
1335  //Tag sttag = (asdm_->getAntenna().getRowByKey( anttag ))->getStationId() ;
1336  Tag sttag( (unsigned int)weatherStationId_, TagType::Station ) ;
1337  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1338  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1339  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1340  //WeatherRow *wrow = wtab.getRowByKey( sttag, tint ) ;
1341  vector<WeatherRow *> *wrows = wtab.getByContext( sttag ) ;
1342  WeatherRow *wrow = (*wrows)[0] ;
1343  unsigned int nrow = wrows->size() ;
1344  //logsink_->postLocally( LogMessage("There are "+String::toString(nrow)+" rows for given context: stationId "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1345  ArrayTime startTime = getMidTime( tint ) ;
1346  if ( startTime < (*wrows)[0]->getTimeInterval().getStart() ) {
1347    temperature = (*wrows)[0]->getTemperature().get() ;
1348    pressure = (*wrows)[0]->getPressure().get() ;
1349    humidity = (*wrows)[0]->getRelHumidity().get() ;
1350    windspeed = (*wrows)[0]->getWindSpeed().get() ;
1351    windaz = (*wrows)[0]->getWindDirection().get() ;
1352  }
1353  else if ( startTime > getEndTime( (*wrows)[nrow-1]->getTimeInterval() ) ) {
1354    temperature = (*wrows)[nrow-1]->getTemperature().get() ;
1355    pressure = (*wrows)[nrow-1]->getPressure().get() ;
1356    humidity = (*wrows)[nrow-1]->getRelHumidity().get() ;
1357    windspeed = (*wrows)[nrow-1]->getWindSpeed().get() ;
1358    windaz = (*wrows)[nrow-1]->getWindDirection().get() ;
1359  }
1360  else {
1361    for ( unsigned int irow = 1 ; irow < wrows->size() ; irow++ ) {
1362      wrow = (*wrows)[irow-1] ;
1363      ArrayTime tStart = wrow->getTimeInterval().getStart() ;
1364      ArrayTime tEnd = (*wrows)[irow]->getTimeInterval().getStart() ;
1365      if ( startTime >= tStart && startTime <= tEnd ) {
1366        temperature = wrow->getTemperature().get() ;
1367        pressure = wrow->getPressure().get() ;
1368        humidity = wrow->getRelHumidity().get() ;
1369        windspeed = wrow->getWindSpeed().get() ;
1370        windaz = wrow->getWindDirection().get() ;
1371        break ;
1372      }
1373    }
1374  }
1375
1376  // Pa -> hPa
1377  pressure /= 100.0 ;
1378
1379  return ;
1380}
1381
1382void ASDMReader::processStation()
1383{
1384  antennaPad_.resize(0) ;
1385  weatherStation_.resize(0) ;
1386  StationTable &stab = asdm_->getStation() ;
1387  vector<StationRow *> srows = stab.get() ;
1388  for ( unsigned int irow = 0 ; irow < srows.size() ; irow++ ) {
1389    StationType stype = srows[irow]->getType() ;
1390    Tag stag = srows[irow]->getStationId() ;
1391    if ( stype == ANTENNA_PAD )
1392      antennaPad_.push_back( stag ) ;
1393    else if ( stype == WEATHER_STATION )
1394      weatherStation_.push_back( stag ) ;
1395  }
1396
1397   weatherStationId_ = getClosestWeatherStation() ;
1398}
1399
1400int ASDMReader::getClosestWeatherStation()
1401{
1402  if ( weatherStation_.size() == 0 )
1403    return -1 ;
1404
1405  Tag atag( antennaId_, TagType::Antenna ) ;
1406  Tag stag = (asdm_->getAntenna().getRowByKey( atag ))->getStationId() ;
1407  vector<double> apos( 3 ) ;
1408  StationTable &stab = asdm_->getStation() ;
1409  StationRow *srow = stab.getRowByKey( stag ) ;
1410  vector<Length> pos = srow->getPosition() ;
1411  apos[0] = pos[0].get() ;
1412  apos[1] = pos[1].get() ;
1413  apos[2] = pos[2].get() ;
1414 
1415  double eps = 1.0e20 ;
1416  int retval = -1 ;
1417  for ( unsigned int ir = 0 ; ir < weatherStation_.size() ; ir++ ) {
1418    srow = stab.getRowByKey( weatherStation_[ir] ) ;
1419    vector<Length> wpos = srow->getPosition() ;
1420    double dist = (apos[0]-wpos[0].get())*(apos[0]-wpos[0].get())
1421      + (apos[1]-wpos[1].get())*(apos[1]-wpos[1].get())
1422      + (apos[2]-wpos[2].get())*(apos[2]-wpos[2].get()) ;
1423    if ( dist < eps ) {
1424      retval = (int)(weatherStation_[ir].getTagValue()) ;
1425    }
1426  }
1427
1428  return retval ;
1429}
1430
1431void ASDMReader::getPointingInfo( unsigned int idx,
1432                                  vector<double> &dir,
1433                                  double &az,
1434                                  double &el,
1435                                  vector<double> &srate )
1436{
1437  String funcName = "getPointingInfo" ;
1438
1439  dir.resize(0) ;
1440  az = -1.0 ;
1441  el = -1.0 ;
1442  srate.resize(0) ;
1443
1444  Tag atag( antennaId_, TagType::Antenna ) ;
1445  unsigned int index = dataIdList_[idx] ;
1446  vector<PointingRow *> *prows = asdm_->getPointing().getByContext( atag ) ;
1447
1448  if ( prows == 0 )
1449    return ;
1450
1451  PointingRow *prow ;
1452  PointingRow *qrow ;
1453  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1454  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1455  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1456
1457  unsigned int nrow = prows->size() ;
1458  //logsink_->postLocally( LogMessage("There are " << nrow << " rows for given context: antennaId "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
1459
1460//   for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1461//     prow = (*prows)[irow] ;
1462//     ArrayTimeInterval ati = prow->getTimeInterval() ;
1463//     ArrayTime pst = ati.getStart() ;
1464//     ArrayTime pet( ati.getStartInMJD()+ati.getDurationInDays() ) ;
1465//     logsink_->postLocally( LogMessage("start: "+pst.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1466//     logsink_->postLocally( LogMessage("end: "+pet.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1467//   }
1468 
1469  srate.resize(2) ;
1470  srate[0] = 0.0 ;
1471  srate[1] = 0.0 ;
1472  az = 0.0 ;
1473  el = 0.0 ;
1474  //double tcen = 0.0 ;
1475  double tcen = getMidTime( tint ).getMJD() ;
1476
1477  //
1478  // shape of pointingDirection is (numSample,2) if usePolynomial = False, while it is
1479  // (numTerm,2) if usePolynomial = True.
1480  //
1481  // In the former case, typical sampled time interval is 48msec, which is very small
1482  // compared with typical integration time (~a few sec).
1483  // Scan rate for this case is always [0,0] (or get slope?).
1484  //
1485  // In the later case, scan rate is (pointingDirection[1][0],pointingDirection[1][1])
1486  //
1487  // PointingDirection seems to store AZEL
1488  //
1489  ArrayTimeInterval pTime0 = (*prows)[0]->getTimeInterval() ;
1490  ArrayTimeInterval pTime1 = (*prows)[nrow-1]->getTimeInterval() ;
1491  //if ( tint.getStartInMJD()+tint.getDurationInDays() < pTime0.getStartInMJD() ) {
1492  if ( getEndTime( tint ) < getStartTime( pTime0 ) ) {
1493    logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint < ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1494    prow = (*prows)[0] ;
1495    vector< vector<double> > dirA = pointingDir( prow ) ;
1496    az = dirA[0][0] ;
1497    el = dirA[0][1] ;
1498    if ( prow->getUsePolynomials() && prow->getNumTerm() > 1 ) {
1499      srate[0] = dirA[1][0] ;
1500      srate[1] = dirA[1][1] ;
1501    }     
1502  }
1503  //else if ( tint.getStartInMJD() > pTime1.getStartInMJD()+pTime1.getDurationInDays() ) {
1504  else if ( getStartTime( tint ) > getEndTime( pTime1 ) ) {
1505    logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint > ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1506    prow = (*prows)[nrow-1] ;
1507    int numSample = prow->getNumSample() ;
1508    vector< vector<double> > dirA = pointingDir( prow ) ;
1509    if ( prow->getUsePolynomials() ) {
1510      az = dirA[0][0] ;
1511      el = dirA[0][1] ;
1512      if ( prow->getNumTerm() > 1 ) {
1513        srate[0] = dirA[1][0] ;
1514        srate[1] = dirA[1][1] ;
1515      }
1516    }
1517    else {
1518      az = dirA[numSample-1][0] ;
1519      el = dirA[numSample-1][1] ;
1520    }
1521  }
1522  else {
1523    ArrayTime startTime = tint.getStart() ;
1524    ArrayTime endTime = getEndTime( tint ) ;
1525    int row0 = -1 ;
1526    int row1 = -1 ;
1527    int row2 = -1 ;
1528    double dt0 = getMidTime( tint ).getMJD() ;
1529    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1530      prow = (*prows)[irow] ;
1531      double dt = getMidTime( tint ).getMJD() - getMidTime( prow->getTimeInterval() ).getMJD() ;
1532      if ( dt > 0 && dt < dt0 ) {
1533        dt0 = dt ;
1534        row2 = irow ;
1535      }
1536      if ( prow->getTimeInterval().contains( startTime ) )
1537        row0 = irow ;
1538      else if ( prow->getTimeInterval().contains( endTime ) )
1539        row1 = irow ;
1540      if ( row0 != -1 && row1 != -1 )
1541        break ;
1542    }
1543    //logsink_->postLocally( LogMessage("row0 = "+String::toString(row0)+", row1 = "+String::toString(row1)+", row2 = "+String::toString(row2),LogOrigin(className_,funcName,WHERE)) ) ;
1544    if ( row0 == -1 && row1 == -1 ) {
1545      prow = (*prows)[row2] ;
1546      qrow = (*prows)[row2+1] ;
1547      double t0 = getEndTime( prow->getTimeInterval() ).getMJD() ;
1548      double t1 = qrow->getTimeInterval().getStartInMJD() ;
1549      vector< vector<double> > dirA = pointingDir( prow ) ;
1550      vector< vector<double> > dirB = pointingDir( qrow ) ;
1551      double da0 = dirA[0][0] ;
1552      double db0 = dirB[0][0] ;
1553      double da1 = dirA[0][1] ;
1554      double db1 = dirB[0][1] ;
1555      if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1556        double dt = ( tcen - t0 ) / ( t1 - t0 ) ;
1557        az = da0 + ( db0 - da0 ) * dt ;
1558        el = da1 + ( db1 - da1 ) * dt ;
1559        if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1560          double ra0 = dirA[1][0] ;
1561          double rb0 = dirB[1][0] ;
1562          double ra1 = dirA[1][1] ;
1563          double rb1 = dirB[1][1] ;
1564          srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1565          srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1566        }
1567      }
1568      else if ( !(qrow->getUsePolynomials()) ) {
1569        double dt = ( tcen - t0 ) / ( t1 - t0 ) ;
1570        az = da0 + ( db0 - da0 ) * dt ;
1571        el = da1 + ( db1 - da1 ) * dt ;
1572      }
1573      //else {
1574      // nothing to do
1575      //}
1576    }
1577    else {
1578      int ndir = 0 ;
1579      if ( row0 == -1 ) {
1580        row0 = row1 ;
1581      }
1582      else if ( row1 == -1 ) {
1583        row1 = row0 ;
1584      }
1585      prow = (*prows)[row0] ;
1586      qrow = (*prows)[row1] ;
1587      if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1588        //logsink_->postLocally( LogMessage("usePolynomial = True",LogOrigin(className_,funcName,WHERE)) ) ;
1589        if ( row0 == row1 ) {
1590          prow = (*prows)[row0] ;
1591          vector< vector<double> > dirA = pointingDir( prow ) ;
1592          az = dirA[0][0] ;
1593          el = dirA[0][1] ;
1594          if ( prow->getNumTerm() > 1 ) {
1595            srate[0] = dirA[1][0] ;
1596            srate[1] = dirA[1][1] ;
1597          }
1598        }
1599        else {
1600          prow = (*prows)[row0] ;
1601          qrow = (*prows)[row1] ;
1602          vector< vector<double> > dirA = pointingDir( prow ) ;
1603          vector< vector<double> > dirB = pointingDir( qrow ) ;
1604          double t0 = getMidTime( prow->getTimeInterval() ).getMJD() ;
1605          double t1 = getMidTime( qrow->getTimeInterval() ).getMJD() ;
1606          double dt = ( tcen - t0 ) / ( t1 - t0 ) ;
1607          double da0 = dirA[0][0] ;
1608          double db0 = dirB[0][0] ;
1609          double da1 = dirA[0][1] ;
1610          double db1 = dirB[0][1] ;
1611          az = da0 + ( db0 - da0 ) * dt ;
1612          el = da1 + ( db1 - db0 ) * dt ;
1613          if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1614            double ra0 = dirA[1][0] ;
1615            double rb0 = dirB[1][0] ;
1616            double ra1 = dirA[1][1] ;
1617            double rb1 = dirB[1][1] ;
1618            srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1619            srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1620          }
1621        }
1622      }
1623      else if ( prow->getUsePolynomials() == qrow->getUsePolynomials() ) {
1624        //logsink_->postLocally( LogMessage("numSample == numTerm",LogOrigin(className_,funcName,WHERE)) ) ;
1625        tcen = 0.0 ;
1626        for ( int irow = row0 ; irow <= row1 ; irow++ ) {
1627          prow = (*prows)[irow] ;
1628          int numSample = prow->getNumSample() ;
1629          //logsink_->postLocally( LogMessage("numSample = "+String::toString(numSample),LogOrigin(className_,funcName,WHERE)) ) ;
1630          vector< vector<double> > dirA = pointingDir( prow ) ;
1631          if ( prow->isSampledTimeIntervalExists() ) {
1632            //logsink_->postLocally( LogMessage("sampledTimeIntervalExists",LogOrigin(className_,funcName,WHERE)) ) ;
1633            vector<ArrayTimeInterval> stime = prow->getSampledTimeInterval() ;
1634            for ( int isam = 0 ; isam < numSample ; isam++ ) {
1635              //if ( tint.overlaps( stime[isam] ) ) {
1636              if ( tint.contains( stime[isam] ) ) {
1637                az += dirA[isam][0] ;
1638                el += dirA[isam][1] ;
1639                tcen += getMidTime( stime[isam] ).getMJD() ;
1640                ndir++ ;
1641              }
1642            }
1643          }
1644          else {
1645            double sampleStart = prow->getTimeInterval().getStartInMJD() ;
1646            double sampleInterval = prow->getTimeInterval().getDurationInDays() / (double)numSample ;
1647            //logsink_->postLocally( LogMessage("sampleStart = "+String::toString(sampleStart),LogOrigin(className_,funcName,WHERE)) )
1648            //logsink_->postLocally( LogMessage("sampleInterval = "+String::toString(sampleInterval),LogOrigin(className_,funcName,WHERE)) ) ;
1649            //logsink_->postLocally( LogMessage("tint = "+tint.toString(),LogOrigin(className_,funcName,WHERE)) ) ;
1650            for ( int isam = 0 ; isam < numSample ; isam++ ) {
1651              ArrayTimeInterval stime( sampleStart+isam*sampleInterval, sampleInterval ) ;
1652              //if ( tint.overlaps( stime ) ) {
1653              if ( tint.contains( stime ) ) {
1654                az += dirA[isam][0] ;
1655                el += dirA[isam][1] ;
1656                tcen += getMidTime( stime ).getMJD() ;
1657                ndir++ ;
1658              }
1659            }
1660          }
1661        }
1662        if ( ndir > 1 ) {
1663          az /= (double)ndir ;
1664          el /= (double)ndir ;
1665          tcen /= (double)ndir ;
1666        }
1667      }
1668      //else {
1669      // nothing to do
1670      //}
1671    }
1672   
1673    AntennaRow *arow = asdm_->getAntenna().getRowByKey( Tag( antennaId_, TagType::Antenna ) ) ;
1674    StationRow *srow = arow->getStationUsingStationId() ;
1675    vector<Length> antposL = srow->getPosition() ;
1676    casa::Vector<casa::Double> antpos( 3 ) ;
1677    for ( int i = 0 ; i < 3 ; i++ )
1678      antpos[i] = antposL[i].get() ;
1679    //logsink_->postLocally( LogMessage("tcen = "+String::toString(tcen),LogOrigin(className_,funcName,WHERE)) ) ;
1680    //logsink_->postLocally( LogMessage("antpos = "+String::toString(antpos),LogOrigin(className_,funcName,WHERE)) ) ;
1681    toJ2000( dir, az, el, tcen, antpos ) ;
1682
1683  }
1684
1685  return ;
1686}
1687
1688ArrayTime ASDMReader::getMidTime( const ArrayTimeInterval &t )
1689{
1690  return ArrayTime( t.getStartInMJD() + 0.5 * t.getDurationInDays() ) ;
1691}
1692
1693ArrayTime ASDMReader::getEndTime( const ArrayTimeInterval &t )
1694{
1695  return ArrayTime( t.getStartInMJD() + t.getDurationInDays() ) ;
1696}
1697
1698ArrayTime ASDMReader::getStartTime( const ArrayTimeInterval &t )
1699{
1700  return ArrayTime( t.getStartInMJD() ) ;
1701}
1702
1703void ASDMReader::toJ2000( vector<double> &dir,
1704                          double az,
1705                          double el,
1706                          double mjd,
1707                          casa::Vector<casa::Double> antpos )
1708{
1709  String funcName = "toJ2000" ;
1710
1711  vector<double> azel( 2 ) ;
1712  azel[0] = az ;
1713  azel[1] = el ;
1714  dir = toJ2000( azel, "AZELGEO", mjd, antpos ) ;
1715}
1716
1717vector<double> ASDMReader::toJ2000( vector<double> dir,
1718                                    casa::String dirref,
1719                                    double mjd,
1720                                    casa::Vector<casa::Double> antpos )
1721{
1722  casa::String funcName = "toJ2000" ;
1723
1724  vector<double> newd( dir ) ;
1725  if ( dirref != "J2000" ) {
1726    casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
1727    casa::Vector<casa::Quantity> qantpos( 3 ) ;
1728    qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
1729    qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
1730    qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
1731    casa::MPosition mp( casa::MVPosition( qantpos ),
1732                        casa::MPosition::ITRF ) ;
1733    //ostringstream oss ;
1734    //mp.print( oss ) ;
1735    //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
1736   
1737    casa::MeasFrame mf( me, mp ) ;
1738    casa::MDirection::Types dirtype ;
1739    casa::Bool b = casa::MDirection::getType( dirtype, dirref ) ;
1740    if ( b ) {
1741      casa::MDirection::Convert toj2000( dirtype,
1742                                         casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
1743      casa::Vector<casa::Double> cdir = toj2000( dir ).getAngle( "rad" ).getValue() ;
1744      //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
1745      newd[0] = (double)(cdir[0]) ;
1746      newd[1] = (double)(cdir[1]) ;
1747    }
1748  }
1749  return newd ;
1750}
1751
1752void ASDMReader::setLogger( CountedPtr<LogSinkInterface> &logsink )
1753{
1754  logsink_ = logsink ;
1755}
1756
1757string ASDMReader::getFrame()
1758{
1759  casa::String funcName = "getFrame" ;
1760 
1761  // default is TOPO
1762  string frame = "TOPO" ;
1763
1764  SpectralWindowTable &spwtab = asdm_->getSpectralWindow() ;
1765  vector<SpectralWindowRow *> rows = spwtab.get() ;
1766  vector<FrequencyReferenceCode> measFreqRef( rows.size() ) ;
1767  int nref = 0 ;
1768  for ( unsigned int irow = 0 ; irow < rows.size() ; irow++ ) {
1769    int nchan = rows[irow]->getNumChan() ;
1770    if ( nchan != 4 ) {
1771      if ( rows[irow]->isMeasFreqRefExists() ) {
1772        measFreqRef[nref] = rows[irow]->getMeasFreqRef() ;
1773        nref++ ;
1774      }
1775    }
1776  }
1777  if ( nref != 0 ) {
1778    frame = CFrequencyReferenceCode::toString( measFreqRef[0] ) ;
1779  }
1780 
1781  //logsink_->postLocally( LogMessage("frame = "+String::toString(frame),LogOrigin(className_,funcName,WHERE)) ) ;
1782
1783  return frame ;
1784}
1785
1786int ASDMReader::getNumIFs()
1787{
1788  casa::String funcName = "getNumIFs" ;
1789
1790  int nif = 0 ;
1791  vector<SpectralWindowRow *> rows = asdm_->getSpectralWindow().get() ;
1792  unsigned int nrow = rows.size() ;
1793  // check if all rows have freqGroup attribute
1794  bool freqGroupExists = true ;
1795  bool countedWvr = false ;
1796  for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1797    freqGroupExists &= rows[irow]->isFreqGroupExists() ;
1798    if ( rows[irow]->getNumChan() == 4 ) {
1799      if ( !countedWvr ) {
1800        countedWvr = true ;
1801        nif++ ;
1802      }
1803    }
1804    else {
1805      nif++ ;
1806    }
1807  }
1808 
1809  if ( freqGroupExists ) {
1810    vector<int> freqGroup(0) ;
1811    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1812      int fg = rows[irow]->getFreqGroup() ;
1813      if ( (int)count( freqGroup.begin(), freqGroup.end(), fg ) == 0 ) {
1814        freqGroup.push_back( fg ) ;
1815      }
1816    }
1817    nif = freqGroup.size() ;
1818  }
1819
1820  //logsink_->postLocally( LogMessage("nif = "+String::toString(nif),LogOrigin(className_,funcName,WHERE)) ) ;
1821
1822  return nif ;
1823}
1824
1825SysCalRow *ASDMReader::getSysCalRow( unsigned int idx )
1826{
1827  String funcName = "getSysCalRow" ;
1828
1829  SysCalRow *row = 0 ;
1830  unsigned int index = dataIdList_[idx] ;
1831  Tag anttag( antennaId_, TagType::Antenna ) ;
1832  int feedid = vmsData_->v_feedId1[index] ;
1833  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1834  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1835  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1836  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1837  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1838  //int nchan = asdm_->getSpectralWindow().getRowByKey(spwtag)->getNumChan() ;
1839  vector< vector<float> > defaultTsys( 1, vector<float>( 1, 1.0 ) ) ;
1840  SysCalTable &sctab = asdm_->getSysCal() ;
1841  //vector<SysCalRow *> rows = sctab.get() ;
1842  vector<SysCalRow *> *rows = sctab.getByContext( anttag, spwtag, feedid ) ;
1843  //if ( nrow == 0 ) {
1844  if ( rows == 0 ) {
1845    //logsink_->postLocally( LogMessage("no rows in SysCal table",LogOrigin(className_,funcName,WHERE)) ) ;
1846    row = 0 ;
1847  }
1848  else {
1849    unsigned int nrow = rows->size() ;
1850    //logsink_->postLocally( LogMessage("nrow = "+String::toString(nrow),LogOrigin(className_,funcName,WHERE)) ) ;
1851    int scindex = -1 ;
1852    if ( nrow == 1 ) {
1853      scindex = 0 ;
1854    }
1855    else if ( getEndTime( tint ) <= getStartTime( (*rows)[0]->getTimeInterval() ) )
1856      scindex = 0 ;
1857    else {
1858      for ( unsigned int irow = 0 ; irow < nrow-1 ; irow++ ) {
1859        ArrayTime t = getMidTime( tint ) ;
1860        if ( t > getStartTime( (*rows)[irow]->getTimeInterval() )
1861             && t <= getStartTime( (*rows)[irow+1]->getTimeInterval() ) ) {
1862          scindex = irow ;
1863          break ;
1864        }
1865      }
1866      if ( scindex == -1 )
1867        scindex = nrow-1 ;
1868    }
1869    //logsink_->postLocally( LogMessage("irow = "+String::toString(scindex),LogOrigin(className_,funcName,WHERE)) ) ;
1870    row = (*rows)[scindex] ;
1871  }
1872  return row ;
1873}
1874
1875double ASDMReader::limitedAngle( double angle )
1876{
1877  if ( angle > C::pi )
1878    while ( angle > C::pi ) angle -= C::_2pi ;
1879  else if ( angle < -C::pi )
1880    while ( angle < -C::pi ) angle += C::_2pi ;
1881  return angle ;
1882}
1883
1884vector< vector<double> > ASDMReader::pointingDir( PointingRow *row )
1885{
1886  //String funcName = "pointingDir" ;
1887  vector< vector<Angle> > aTar = row->getTarget() ;
1888  vector< vector<Angle> > aOff = row->getOffset() ;
1889  vector< vector<Angle> > aDir = row->getPointingDirection() ;
1890  vector< vector<Angle> > aEnc = row->getEncoder() ;
1891  unsigned int n = aTar.size() ;
1892  vector< vector<double> > dir( n ) ;
1893  double factor = 1.0 / cos( aTar[0][1].get() ) ;
1894  for ( unsigned int i = 0 ; i < n ; i++ ) {
1895    dir[i].resize( 2 ) ;
1896    /**
1897     * This is approximate way to add offset taking tracking error
1898     * into account
1899     *
1900     * az = dir[0][0] = target[0][0] + offset[0][0] / cos(el)
1901     *                 + encorder[0][0] - direction[0][0]
1902     * el = dir[0][1] = target[0][1] + offset[0][1]
1903     *                 + encorder[0][1] - direction[0][1]
1904     **/
1905    dir[i][0] = aTar[i][0].get() + factor * aOff[i][0].get()
1906               + aEnc[i][0].get() - aDir[i][0].get() ;
1907    dir[i][1] = aTar[i][1].get() + aOff[i][1].get()
1908               + aEnc[i][1].get() - aDir[i][1].get() ;
1909    //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)) ) ;
1910  }
1911  return dir ;
1912}
Note: See TracBrowser for help on using the repository browser.