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

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

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