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

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

Number of IFs is not a number of rows of SpectralWindow? table,
but a number of frequency groups.
When freqGroup is not defined, number of IFs is set to number
of non-WVR spectral windows plus 1 (WVR spectral window).


File size: 63.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
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  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  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
345  }
346  else {
347    // frequency reference is TOPOCENT by default
348    freqref = "TOPOCENT" ;
349  }
350
351
352  PolarizationTable &ptab = asdm_->getPolarization() ;
353  vector<PolarizationRow *> prows = ptab.get() ;
354  vector<int> npols ;
355  refidx = -1 ;
356  for ( unsigned int irow = 0 ; irow < prows.size() ; irow++ ) {
357    npols.push_back( prows[irow]->getNumCorr() ) ;
358    if ( refidx == -1 && npols[irow] != 1 )
359      refidx = irow ;
360  }
361  if ( refidx == -1 )
362    refidx = 0 ;
363
364  // npol
365  npol = casa::Int( *max_element( npols.begin(), npols.end() ) ) ;
366
367  // poltype
368  vector<StokesParameter> corrType = prows[refidx]->getCorrType() ;
369  if ( corrType[0] == I ||
370       corrType[0] == Q ||
371       corrType[0] == U ||
372       corrType[0] == V )
373    poltype = "stokes" ;
374  else if ( corrType[0] == RR ||
375            corrType[0] == RL ||
376            corrType[0] == LR ||
377            corrType[0] == LL )
378    poltype = "circular" ;
379  else if ( corrType[0] == XX ||
380            corrType[0] == XY ||
381            corrType[0] == YX ||
382            corrType[0] == YY )
383    poltype = "linear" ;
384  else if ( corrType[0] == PLINEAR ||
385            corrType[0] == PANGLE ) {
386    poltype = "linpol" ;
387  }
388  else {
389    poltype = "notype" ;
390  }
391
392
393  FeedTable &ftab = asdm_->getFeed() ;
394  vector<FeedRow *> frows = ftab.get() ;
395  vector<int> nbeams ;
396  for ( unsigned int irow = 0 ; irow < frows.size() ; irow++ ) {
397    if ( frows[irow]->isFeedNumExists() )
398      nbeams.push_back( frows[irow]->getFeedNum() ) ;
399    else
400      nbeams.push_back( 1 ) ;
401  }
402
403  // nbeam
404  nbeam = casa::Int( *max_element( nbeams.begin(), nbeams.end() ) ) ;
405
406  // fluxunit
407  // tentatively set 'K'? or empty?
408  fluxunit = "K" ;
409
410  // equinox
411  // tentatively set 2000.0
412  equinox = 2000.0 ;
413
414  // epoch
415  // tentatively set "UTC"
416  epoch = "UTC" ;
417
418  // obstype
419  // at the moment take observingMode attribute in SBSummary table
420  SBSummaryRow *sbrow = ebrow->getSBSummaryUsingSBSummaryId() ;
421  vector<string> obsmode = sbrow->getObservingMode() ;
422  obstype = "" ;
423  for ( unsigned int imode = 0 ; imode < obsmode.size() ; imode++ ) {
424    obstype += casa::String(obsmode[imode]) ;
425    if ( imode != obsmode.size()-1 )
426      obstype += "#" ;
427  }
428}
429
430void ASDMReader::selectConfigDescription()
431{
432  String funcName = "selectConfigDescription" ;
433
434  vector<ConfigDescriptionRow *> cdrows = asdm_->getConfigDescription().get() ;
435  vector<Tag> cdidTags ;
436  for ( unsigned int irow = 0 ; irow < cdrows.size() ; irow++ ) {
437    //logsink_->postLocally( LogMessage("correlationMode["+String::toString(irow)+"] = "+String::toString(cdrows[irow]->getCorrelationMode()),LogOrigin(className_,funcName,WHERE)) ) ;
438    if ( cdrows[irow]->getCorrelationMode() != CROSS_ONLY ) {
439      cdidTags.push_back( cdrows[irow]->getConfigDescriptionId() ) ;
440    }
441  }
442
443  configDescIdList_.resize( cdidTags.size() ) ;
444  for ( unsigned int i = 0 ; i < cdidTags.size() ; i++ ) {
445    configDescIdList_[i] = casa::uInt( cdidTags[i].getTagValue() ) ;
446  }
447}
448
449void ASDMReader::selectFeed()
450{
451  feedIdList_.resize(0) ;
452  vector<FeedRow *> frows = asdm_->getFeed().get() ;
453  Tag atag( antennaId_, TagType::Antenna ) ;
454  for ( unsigned int irow = 0 ; irow < frows.size() ; irow++ ) {
455    casa::uInt feedId = (casa::uInt)(frows[irow]->getFeedId() ) ;
456    if ( casa::anyEQ( feedIdList_, feedId ) )
457      continue ;
458    if ( frows[irow]->getAntennaId() == atag ) {
459      unsigned int oldsize = feedIdList_.size() ;
460      feedIdList_.resize( oldsize+1, true ) ;
461      feedIdList_[oldsize] = feedId ;
462    }
463  }
464}
465
466casa::Vector<casa::uInt> ASDMReader::getFieldIdList()
467{
468  String funcName = "getFieldIdList" ;
469
470  vector<FieldRow *> frows = asdm_->getField().get() ;
471  fieldIdList_.resize( frows.size() ) ;
472  for ( unsigned int irow = 0 ; irow < frows.size() ; irow++ ) {
473    //logsink_->postLocally( LogMessage("fieldId["+String::toString(irow)+"]="+String(frows[irow]->getFieldId().toString()),LogOrigin(className_,funcName,WHERE)) ) ;
474    fieldIdList_[irow] = frows[irow]->getFieldId().getTagValue() ;
475  }
476
477  return fieldIdList_ ;
478}
479
480casa::uInt ASDMReader::getNumMainRow()
481{
482  casa::uInt nrow = casa::uInt( mainRow_.size() ) ;
483
484  return nrow ;
485}
486
487void ASDMReader::select()
488{
489  // selection by input CorrelationMode
490  EnumSet<CorrelationMode> esCorrs ;
491  sdmBin_->select( corrMode_ ) ;
492
493  // selection by TimeSampling
494  sdmBin_->select( timeSampling_ ) ;
495
496  // selection by SpectralResolutionType
497  sdmBin_->select( resolutionType_ ) ;
498
499  // selection by AtmPhaseCorrection and output CorrelationMode
500  EnumSet<AtmPhaseCorrection> esApcs ;
501  esApcs.set( apc_ ) ;
502  // always take only autocorrelation data
503  Enum<CorrelationMode> esCorr = AUTO_ONLY ;
504  sdmBin_->selectDataSubset( esCorr, esApcs ) ;
505
506  // select valid configDescriptionId
507  selectConfigDescription() ;
508
509  // select valid feedId
510  selectFeed() ;
511}
512
513casa::Bool ASDMReader::setMainRow( casa::uInt irow )
514{
515  casa::Bool status = true ;
516  row_ = irow ;
517
518  unsigned int cdid = mainRow_[row_]->getConfigDescriptionId().getTagValue() ;
519  if ( (int)count(configDescIdList_.begin(),configDescIdList_.end(),cdid) == 0 )
520    status = false ;
521  else {
522    status = (casa::Bool)(sdmBin_->acceptMainRow( mainRow_[row_] )) ;
523  }
524  return status ;
525}
526
527casa::Bool ASDMReader::setMainRow( casa::uInt configDescId, casa::uInt fieldId )
528{
529  clearMainRow() ;
530
531  Tag configDescTag( (unsigned int)configDescId, TagType::ConfigDescription ) ;
532  Tag fieldTag( (unsigned int)fieldId, TagType::Field ) ;
533  mainRow_ = casa::Vector<MainRow *>( *(asdm_->getMain().getByContext( configDescTag, fieldTag ) ) ) ;
534 
535  return true ;
536}
537
538void ASDMReader::clearMainRow()
539{
540  mainRow_.resize(0) ;
541}
542
543void ASDMReader::setupIFNO()
544{
545  String funcName = "setupIFNO" ;
546
547  vector<SpectralWindowRow *> spwrows = asdm_->getSpectralWindow().get() ;
548  unsigned int nrow = spwrows.size() ;
549  ifno_.clear() ;
550  casa::uInt idx = 0 ;
551  casa::uInt wvridx = 0 ;
552  for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
553    casa::uInt index ;
554    if ( isWVR( spwrows[irow] ) ) {
555      //logsink_->postLocally( LogMessage(spwrows[irow]->getSpectralWindowId().toString()+" is WVR",LogOrigin(className_,funcName,WHERE)) ) ;
556      index = wvridx ;
557    }
558    else {
559      index = ++idx ;
560    }
561    ifno_.insert( pair<Tag,casa::uInt>(spwrows[irow]->getSpectralWindowId(),index) ) ;
562    //logsink_->postLocally( LogMessage(spwrows[irow]->getSpectralWindowId().toString()+": IFNO="+String::toString(index),LogOrigin(className_,funcName,WHERE)) ) ;
563  }
564}
565
566bool ASDMReader::isWVR( SpectralWindowRow *row )
567{
568  BasebandName bbname = row->getBasebandName() ;
569  int nchan = row->getNumChan() ;
570  if ( bbname == NOBB && nchan == 4 )
571    return true ;
572  else
573    return false ;
574}
575
576// casa::Vector<casa::uInt> ASDMReader::getDataDescIdList( casa::uInt cdid )
577// {
578//   Tag cdTag( (unsigned int)cdid, TagType::ConfigDescription ) ;
579//   ConfigDescriptionRow *cdrow = asdm_->getConfigDescription().getRowByKey( cdTag ) ;
580//   vector<Tag> ddTags = cdrow->getDataDescriptionId() ;
581//   casa::Vector<casa::uInt> ddidList( ddTags.size() ) ;
582//   for ( unsigned int idd = 0 ; idd < ddTags.size() ; idd++ ) {
583//     ddidList[idd] = ddTags[idd].getTagValue() ;
584//   }
585//   return ddidList ;
586// }
587
588// casa::Vector<casa::uInt> ASDMReader::getSwitchCycleIdList( casa::uInt cdid )
589// {
590//   Tag cdTag( (unsigned int)cdid, TagType::ConfigDescription ) ;
591//   ConfigDescriptionRow *cdrow = asdm_->getConfigDescription().getRowByKey( cdTag ) ;
592//   vector<Tag> scTags = cdrow->getSwitchCycleId() ;
593//   casa::Vector<casa::uInt> scidList( scTags.size() ) ;
594//   for ( unsigned int idd = 0 ; idd < scTags.size() ; idd++ ) {
595//     scidList[idd] = scTags[idd].getTagValue() ;
596//   }
597//   return scidList ;
598// }
599
600// casa::Vector<casa::uInt> ASDMReader::getFeedIdList( casa::uInt cdid )
601// {
602//   String funcName = "getFeedIdList" ;
603//   
604//   Tag cdTag( (unsigned int)cdid, TagType::ConfigDescription ) ;
605//   ConfigDescriptionRow *cdrow = asdm_->getConfigDescription().getRowByKey( cdTag ) ;
606//   casa::Vector<casa::uInt> feedIdList ;
607//   vector<int> feedIds = cdrow->getFeedId() ;
608//   for ( unsigned int ife = 0 ; ife < feedIds.size() ; ife++ ) {
609//     logsink_->postLocally( LogMessage("feedIds["+String::toString(ife)+"]="+String::toString(feedIds[ife]),LogOrigin(className_,funcName,WHERE)) ) ;
610//     if ( casa::anyEQ( feedIdList, casa::uInt( feedIds[ife] ) ) )
611//       continue ;
612//     if ( casa::anyEQ( feedIdList_, casa::uInt( feedIds[ife] ) ) ) {
613//       casa::uInt oldsize = feedIdList.size() ;
614//       feedIdList.resize( oldsize+1, true ) ;
615//       feedIdList[oldsize] = casa::uInt( feedIds[ife] ) ;
616//     }
617//   }
618//   logsink_->postLocally( LogMessage("feedIdList.size() = "+String::toString(feedIdList.size()),LogOrigin(className_,funcName,WHERE)) ) ;
619//   return feedIdList ;
620// }
621
622casa::Bool ASDMReader::setData()
623{
624  String funcName = "setData" ;
625
626  //logsink_->postLocally( LogMessage("try to retrieve binary data",LogOrigin(className_,funcName,WHERE)) ) ;
627 
628//   EnumSet<AtmPhaseCorrection> esApcs ;
629//   esApcs.set( apc_ ) ;
630//   // always take only autocorrelation data
631//   Enum<CorrelationMode> esCorr = AUTO_ONLY ;
632//   vmsData_ = sdmBin_->getDataCols( esCorr, esApcs ) ;
633
634  // 2011/07/06 TN
635  // Workaround to avoid unwanted message from SDMBinData::getDataCols()
636  ostringstream oss ;
637  streambuf *buforg = cout.rdbuf(oss.rdbuf()) ;
638  vmsData_ = sdmBin_->getDataCols() ;
639  cout.rdbuf(buforg) ;
640
641//   logsink_->postLocally( LogMessage("oss.str() = "+oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
642//   cout << "This is test: oss.str()=" << oss.str() << endl ;
643
644
645  //logsink_->postLocally( LogMessage("processorId = "+String::toString(vmsData_->processorId),LogOrigin(className_,funcName,WHERE)) ) ;
646  //logsink_->postLocally( LogMessage("v_time.size() = "+String::toString(vmsData_->v_time.size()),LogOrigin(className_,funcName,WHERE)) ) ;
647  //logsink_->postLocally( LogMessage("   v_time[0] = "+String::toString(vmsData_->v_time[0]),LogOrigin(className_,funcName,WHERE)) ) ;
648  //logsink_->postLocally( LogMessage("v_interval.size() = "+String::toString(vmsData_->v_interval.size()),LogOrigin(className_,funcName,WHERE)) ) ;
649  //logsink_->postLocally( LogMessage("   v_interval[0] = "+String::toString(vmsData_->v_interval[0]),LogOrigin(className_,funcName,WHERE)) ) ;
650  //logsink_->postLocally( LogMessage("v_atmPhaseCorrection.size() = "+String::toString(vmsData_->v_atmPhaseCorrection.size()),LogOrigin(className_,funcName,WHERE)) ) ;
651  //logsink_->postLocally( LogMessage("binNum = "+String::toString(vmsData_->binNum),LogOrigin(className_,funcName,WHERE)) ) ;
652  //logsink_->postLocally( LogMessage("v_projectPath.size() = "+String::toString(vmsData_->v_projectPath.size()),LogOrigin(className_,funcName,WHERE)) ) ;
653  //logsink_->postLocally( LogMessage("v_antennaId1.size() = "+String::toString(vmsData_->v_antennaId1.size()),LogOrigin(className_,funcName,WHERE)) ) ;
654  //logsink_->postLocally( LogMessage("v_antennaId2.size() = "+String::toString(vmsData_->v_antennaId2.size()),LogOrigin(className_,funcName,WHERE)) ) ;
655  //logsink_->postLocally( LogMessage("v_feedId1.size() = "+String::toString(vmsData_->v_feedId1.size()),LogOrigin(className_,funcName,WHERE)) ) ;
656  //logsink_->postLocally( LogMessage("v_feedId2.size() = "+String::toString(vmsData_->v_feedId2.size()),LogOrigin(className_,funcName,WHERE)) ) ;
657  //logsink_->postLocally( LogMessage("v_dataDescId.size() = "+String::toString(vmsData_->v_dataDescId.size()),LogOrigin(className_,funcName,WHERE)) ) ;
658  //logsink_->postLocally( LogMessage("v_timeCentroid.size() = "+String::toString(vmsData_->v_timeCentroid.size()),LogOrigin(className_,funcName,WHERE)) ) ;
659  //logsink_->postLocally( LogMessage("v_exposure.size() = "+String::toString(vmsData_->v_exposure.size()),LogOrigin(className_,funcName,WHERE)) ) ;
660  //logsink_->postLocally( LogMessage("v_numData.size() = "+String::toString(vmsData_->v_numData.size()),LogOrigin(className_,funcName,WHERE)) ) ;
661  //logsink_->postLocally( LogMessage("vv_dataShape.size() = "+String::toString(vmsData_->vv_dataShape.size()),LogOrigin(className_,funcName,WHERE)) ) ;
662  //logsink_->postLocally( LogMessage("v_m_data.size() = "+String::toString(vmsData_->v_m_data.size()),LogOrigin(className_,funcName,WHERE)) ) ;
663  //logsink_->postLocally( LogMessage("v_phaseDir.size() = "+String::toString(vmsData_->v_phaseDir.size()),LogOrigin(className_,funcName,WHERE)) ) ;
664  //logsink_->postLocally( LogMessage("v_stateId.size() = "+String::toString(vmsData_->v_stateId.size()),LogOrigin(className_,funcName,WHERE)) ) ;
665  //logsink_->postLocally( LogMessage("v_msState.size() = "+String::toString(vmsData_->v_msState.size()),LogOrigin(className_,funcName,WHERE)) ) ;
666  //logsink_->postLocally( LogMessage("v_flag.size() = "+String::toString(vmsData_->v_flag.size()),LogOrigin(className_,funcName,WHERE)) ) ;
667
668  dataIdList_.clear() ;
669  unsigned int numTotalData = vmsData_->v_m_data.size() ;
670  for ( unsigned int idata = 0 ; idata < numTotalData ; idata++ ) {
671    if ( vmsData_->v_antennaId1[idata] == (int)antennaId_
672         && vmsData_->v_antennaId2[idata] == (int)antennaId_ )
673      dataIdList_.push_back( idata ) ;
674  }
675  numData_ = dataIdList_.size() ;
676  //logsink_->postLocally( LogMessage("numData_ = "+String::toString(numData_),LogOrigin(className_,funcName,WHERE)) ) ;
677  //logsink_->postLocally( LogMessage("dataSize = "+String::toString(mainRow_[row_]->getDataSize()),LogOrigin(className_,funcName,WHERE)) ) ;
678
679  return true ;
680}
681
682casa::uInt ASDMReader::getIFNo( unsigned int idx )
683{
684  Tag ddTag( vmsData_->v_dataDescId[dataIdList_[idx]], TagType::DataDescription ) ;
685  DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey( ddTag ) ;
686  Tag spwid = ddrow->getSpectralWindowId() ;
687  map<Tag,casa::uInt>::iterator iter = ifno_.find( spwid ) ;
688  if ( iter != ifno_.end() )
689    return iter->second ;
690  else {
691    return 0 ;
692  }
693}
694
695int ASDMReader::getNumPol( unsigned int idx )
696{
697  Tag ddTag( vmsData_->v_dataDescId[dataIdList_[idx]], TagType::DataDescription ) ;
698  DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey( ddTag ) ;
699  PolarizationRow *polrow = ddrow->getPolarizationUsingPolOrHoloId() ;
700  return polrow->getNumCorr() ;
701}
702
703void ASDMReader::getFrequency( unsigned int idx,
704                               double &refpix,
705                               double &refval,
706                               double &incr )
707{
708  String funcName = "getFrequency" ;
709
710  Tag ddTag( vmsData_->v_dataDescId[dataIdList_[idx]], TagType::DataDescription ) ;
711  DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey( ddTag ) ;
712  //Tag spwid = ddrow->getSpectralWindowId() ;
713  SpectralWindowRow *spwrow = ddrow->getSpectralWindowUsingSpectralWindowId() ;
714  int nchan = spwrow->getNumChan() ;
715  if ( nchan == 1 ) {
716    //logsink_->postLocally( LogMessage("channel averaged data",LogOrigin(className_,funcName,WHERE)) ) ;
717    refpix = 0.0 ;
718    incr = spwrow->getTotBandwidth().get() ;
719    if ( spwrow->isChanFreqStartExists() ) {
720      refval = spwrow->getChanFreqStart().get() ;
721    }
722    else if ( spwrow->isChanFreqArrayExists() ) {
723      refval = spwrow->getChanFreqArray()[0].get() ;
724    }
725    else {
726      throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
727    }     
728  }
729  else if ( nchan % 2 ) {
730    // odd
731    //logsink_->postLocally( LogMessage("odd case",LogOrigin(className_,funcName,WHERE)) ) ;
732    refpix = 0.5 * ( (double)nchan - 1.0 ) ;
733    int ic = ( nchan - 1 ) / 2 ;
734    if ( spwrow->isChanWidthExists() ) {
735      incr = spwrow->getChanWidth().get() ;
736    }
737    else if ( spwrow->isChanWidthArrayExists() ) {
738      incr = spwrow->getChanWidthArray()[0].get() ;
739    }
740    else {
741      throw (casa::AipsError( "Either chanWidthArray or chanWidth must exist." )) ;
742    }     
743    if ( spwrow->isChanFreqStepExists() ) {
744      if ( spwrow->getChanFreqStep().get() < 0.0 )
745        incr *= -1.0 ;
746    }
747    else if ( spwrow->isChanFreqArrayExists() ) {
748      vector<Frequency> chanFreqArr = spwrow->getChanFreqArray() ;
749      if ( chanFreqArr[0].get() > chanFreqArr[1].get() )
750        incr *= -1.0 ;
751    }
752    else {
753      throw (casa::AipsError( "Either chanFreqArray or chanFreqStep must exist." )) ;
754    }         
755    if ( spwrow->isChanFreqStartExists() ) {
756      refval = spwrow->getChanFreqStart().get() + refpix * incr ;
757    }
758    else if ( spwrow->isChanFreqArrayExists() ) {
759      refval = spwrow->getChanFreqArray()[ic].get() ;
760    }
761    else {
762      throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
763    }     
764  }
765  else {
766    // even
767    //logsink_->postLocally( LogMessage("even case",LogOrigin(className_,funcName,WHERE)) ) ;
768    refpix = 0.5 * ( (double)nchan - 1.0 ) ;
769    int ic = nchan / 2 ;
770    if ( spwrow->isChanWidthExists() ) {
771      incr = spwrow->getChanWidth().get() ;
772    }
773    else if ( spwrow->isChanWidthArrayExists() ) {
774      incr = spwrow->getChanWidthArray()[0].get() ;
775    }
776    else {
777      throw (casa::AipsError( "Either chanWidthArray or chanWidth must exist." )) ;
778    }
779    if ( spwrow->isChanFreqStepExists() ) {
780      if ( spwrow->getChanFreqStep().get() < 0.0 )
781        incr *= -1.0 ;
782    }
783    else if ( spwrow->isChanFreqArrayExists() ) {
784      vector<Frequency> chanFreqArr = spwrow->getChanFreqArray() ;
785      if ( chanFreqArr[0].get() > chanFreqArr[1].get() )
786        incr *= -1.0 ;
787    }
788    else {
789      throw (casa::AipsError( "Either chanFreqArray or chanFreqStep must exist." )) ;
790    }         
791    if ( spwrow->isChanFreqStartExists() ) {
792      refval = spwrow->getChanFreqStart().get() + refpix * incr ;
793    }
794    else if ( spwrow->isChanFreqArrayExists() ) {
795      vector<Frequency> freqs = spwrow->getChanFreqArray() ;
796      refval = 0.5 * ( freqs[ic-1].get() + freqs[ic].get() ) ;     
797    }
798    else {
799      throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
800    }     
801  }
802}
803
804vector<double> ASDMReader::getRestFrequency( unsigned int idx )
805{
806  vector<double> rf( 0 ) ;
807  unsigned int index = dataIdList_[idx] ;
808  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
809  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
810  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
811  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
812  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
813  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
814  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
815  if ( frow->isSourceIdExists() ) {
816    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
817    int sid = frow->getSourceId() ;
818    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
819    if ( srow->isRestFrequencyExists() ) {
820      //logsink_->postLocally( LogMessage("restFrequency exists",LogOrigin(className_,funcName,WHERE)) ) ;
821      vector<Frequency> restfreq = srow->getRestFrequency() ;
822      rf.resize( restfreq.size() ) ;
823      for ( unsigned int i = 0 ; i < restfreq.size() ; i++ )
824        rf[i] = restfreq[i].get() ;
825    }
826  }
827  return rf ;
828}
829
830double ASDMReader::getTime( unsigned int idx )
831{
832  double tsec = vmsData_->v_time[dataIdList_[idx]] ;
833  return tsec * s2d ;
834}
835
836double ASDMReader::getInterval( unsigned int idx )
837{
838  return vmsData_->v_interval[dataIdList_[idx]] ;
839}
840
841string ASDMReader::getSourceName( unsigned int idx )
842{
843  unsigned int index = dataIdList_[idx] ;
844  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
845  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
846  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
847  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
848  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
849  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
850  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
851  string srcname ;
852  if ( frow->isSourceIdExists() ) {
853    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
854    int sid = frow->getSourceId() ;
855    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
856    srcname = srow->getSourceName() ;
857  }
858  else {
859    srcname = frow->getFieldName() ;
860  }
861  return srcname ;
862}
863
864string ASDMReader::getFieldName( unsigned int idx )
865{
866  int fid = vmsData_->v_fieldId[dataIdList_[idx]] ;
867  Tag ftag( fid, TagType::Field ) ;
868  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
869  ostringstream oss ;
870  oss << frow->getFieldName() << "__" << fid ;
871  return oss.str() ;
872}
873
874int ASDMReader::getSrcType( unsigned int scan,
875                            unsigned int subscan )
876{
877  int srctype = SrcType::NOTYPE ;
878  Tag ebtag = mainRow_[row_]->getExecBlockId() ;
879  ScanRow *scanrow = asdm_->getScan().getRowByKey( ebtag, (int)scan ) ;
880  ScanIntent scanIntent = scanrow->getScanIntent()[0] ;
881  SubscanRow *subrow = asdm_->getSubscan().getRowByKey( ebtag, (int)scan, (int)subscan ) ;
882  SubscanIntent subIntent = subrow->getSubscanIntent() ;
883  SwitchingMode swmode = subrow->getSubscanMode() ;
884  if ( scanIntent == OBSERVE_TARGET ) {
885    // on sky scan
886    if ( swmode == NO_SWITCHING || swmode == POSITION_SWITCHING ) {
887      // position switching
888      // tentatively set NO_SWITCHING = POSITION_SWITCHING
889      if ( subIntent == ON_SOURCE ) {
890        srctype = SrcType::PSON ;
891      }
892      else if ( subIntent == OFF_SOURCE ) {
893        srctype = SrcType::PSOFF ;
894      }
895    }
896    else if ( swmode == FREQUENCY_SWITCHING ) {
897      // frequency switching
898      if ( subIntent == ON_SOURCE ) {
899        srctype = SrcType::FSON ;
900      }
901      else if ( subIntent == OFF_SOURCE ) {
902        srctype = SrcType::FSOFF ;
903      }
904    }
905    else if ( swmode == NUTATOR_SWITCHING ) {
906      // nutator switching
907      if ( subIntent == ON_SOURCE ) {
908        srctype = SrcType::PSON ;
909      }
910      else if ( subIntent == OFF_SOURCE ) {
911        srctype = SrcType::PSOFF ;
912      }
913    }
914    else {
915      // other switching mode
916      if ( subIntent == ON_SOURCE ) {
917        srctype = SrcType::SIG ;
918      }
919      else if ( subIntent == OFF_SOURCE ) {
920        srctype = SrcType::REF ;
921      }
922    }
923  }
924  else {
925    // off sky (e.g. calibrator device) scan
926    if ( subIntent == ON_SOURCE ) {
927      srctype = SrcType::SIG ;
928    }
929    else if ( subIntent == OFF_SOURCE ) {
930      srctype = SrcType::REF ;
931    }
932    else if ( subIntent == HOT ) {
933      srctype = SrcType::HOT ;
934    }
935    else if ( subIntent == AMBIENT ) {
936      srctype = SrcType::SKY ;
937    }
938    else {
939      srctype = SrcType::CAL ;
940    }
941  }
942
943  return srctype ;
944}
945
946unsigned int ASDMReader::getSubscanNo( unsigned int idx )
947{
948  //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)) ) ;
949  return vmsData_->v_msState[dataIdList_[idx]].subscanNum ;
950}
951
952vector<double> ASDMReader::getSourceDirection( unsigned int idx )
953{
954  vector<double> dir( 2, 0.0 ) ;
955  unsigned int index = dataIdList_[idx] ;
956  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
957  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
958  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
959  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
960  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
961  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
962  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
963  string srcname ;
964  if ( frow->isSourceIdExists() ) {
965    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
966    int sid = frow->getSourceId() ;
967    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
968    vector<Angle> srcdir = srow->getDirection() ;
969    if ( srow->isDirectionCodeExists() ) {
970      DirectionReferenceCode dircode = srow->getDirectionCode() ;
971      if ( dircode != J2000 ) {
972        // if not J2000, need direction conversion
973      }
974    }
975    dir[0] = srcdir[0].get() ;
976    dir[1] = srcdir[1].get() ;
977  }
978  return dir ;
979}
980
981vector<double> ASDMReader::getSourceProperMotion( unsigned int idx )
982{
983  vector<double> pm( 2, 0.0 ) ;
984  unsigned int index = dataIdList_[idx] ;
985  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
986  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
987  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
988  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
989  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
990  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
991  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
992  string srcname ;
993  if ( frow->isSourceIdExists() ) {
994    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
995    int sid = frow->getSourceId() ;
996    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
997    vector<AngularRate> srcpm = srow->getProperMotion() ;
998    pm[0] = srcpm[0].get() ;
999    pm[1] = srcpm[1].get() ;
1000  }
1001  return pm ;
1002}
1003
1004double ASDMReader::getSysVel( unsigned int idx )
1005{
1006  double sysvel = 0.0 ;
1007  unsigned int index = dataIdList_[idx] ;
1008  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1009  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1010  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1011  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1012  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1013  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1014  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1015  string srcname ;
1016  if ( frow->isSourceIdExists() ) {
1017    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1018    int sid = frow->getSourceId() ;
1019    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1020    if ( srow->isSysVelExists() ) {
1021      vector<Speed> sysvelV = srow->getSysVel() ;
1022      if ( sysvelV.size() > 0 )
1023        sysvel = sysvelV[0].get() ;
1024    }
1025  }
1026  return sysvel ;
1027}
1028
1029unsigned int ASDMReader::getFlagRow( unsigned int idx )
1030{
1031  return vmsData_->v_flag[dataIdList_[idx]] ;
1032}
1033
1034vector<unsigned int> ASDMReader::getDataShape( unsigned int idx )
1035{
1036  return vmsData_->vv_dataShape[dataIdList_[idx]] ;
1037}
1038
1039float * ASDMReader::getSpectrum( unsigned int idx )
1040{
1041  map<AtmPhaseCorrection, float*> data = vmsData_->v_m_data[dataIdList_[idx]] ;
1042  //map<AtmPhaseCorrection, float*>::iterator iter = data.find(AP_UNCORRECTED) ;
1043  map<AtmPhaseCorrection, float*>::iterator iter = data.find(apc_) ;
1044  float *autoCorr = iter->second ;
1045  return autoCorr ;
1046}
1047
1048// bool * ASDMReader::getFlagChannel( unsigned int idx )
1049// {
1050//   return 0 ;
1051// }
1052
1053vector< vector<float> > ASDMReader::getTsys( unsigned int idx )
1054{
1055  unsigned int index = dataIdList_[idx] ;
1056  Tag anttag( antennaId_, TagType::Antenna ) ;
1057  int feedid = vmsData_->v_feedId1[index] ;
1058  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1059  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1060  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1061  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1062  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1063  //int nchan = asdm_->getSpectralWindow().getRowByKey(spwtag)->getNumChan() ;
1064  vector< vector<float> > defaultTsys( 1, vector<float>( 1, 1.0 ) ) ;
1065  SysCalTable &sctab = asdm_->getSysCal() ;
1066  if ( sctab.size() == 0 )
1067    return defaultTsys ;
1068  SysCalRow *scrow = sctab.getRowByKey( anttag, spwtag, tint, feedid ) ;
1069  if ( scrow->isTsysSpectrumExists() ) {
1070    vector< vector<Temperature> > tsysSpec = scrow->getTsysSpectrum() ;
1071    unsigned int numReceptor = tsysSpec.size() ;
1072    unsigned int numChan = tsysSpec[0].size() ;
1073    vector< vector<float> > tsys( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1074    for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1075      for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1076        tsys[ir][ic] = tsysSpec[ir][ic].get() ;
1077      }
1078    }
1079    return tsys ;
1080  }
1081//   else if ( scrow->isTsysExists() ) {
1082//     vector<Temperature> tsysScalar = scrow->getTsys() ;
1083//     unsigned int numReceptor = tsysScalar.size() ;
1084//     vector< vector<float> > tsys( numReceptor ) ;
1085//     for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1086//       tsys[ir] = vector<float>( 1, tsysScalar[ir].get() ) ;
1087//     return tsys ;
1088//   }
1089  else {
1090    return defaultTsys ;
1091  }
1092}
1093
1094vector< vector<float> > ASDMReader::getTcal( unsigned int idx )
1095{
1096  unsigned int index = dataIdList_[idx] ;
1097  Tag anttag( antennaId_, TagType::Antenna ) ;
1098  int feedid = vmsData_->v_feedId1[index] ;
1099  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1100  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1101  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1102  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1103  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1104  //int nchan = asdm_->getSpectralWindow().getRowByKey(spwtag)->getNumChan() ;
1105  vector< vector<float> > defaultTcal( 1, vector<float>( 1, 1.0 ) ) ;
1106  SysCalTable &sctab = asdm_->getSysCal() ;
1107  if ( sctab.size() == 0 )
1108    return defaultTcal ;
1109  SysCalRow *scrow = sctab.getRowByKey( anttag, spwtag, tint, feedid ) ;
1110  if ( scrow->isTcalSpectrumExists() ) {
1111    vector< vector<Temperature> > tcalSpec = scrow->getTcalSpectrum() ;
1112    unsigned int numReceptor = tcalSpec.size() ;
1113    unsigned int numChan = tcalSpec[0].size() ;
1114    vector< vector<float> > tcal( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1115    for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1116      for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1117        tcal[ir][ic] = tcalSpec[ir][ic].get() ;
1118      }
1119    }
1120    return tcal ;
1121  }
1122//   else if ( scrow->isTcalExists() ) {
1123//     vector<Temperature> tcalScalar = scrow->getTcal() ;
1124//     unsigned int numReceptor = tcalScalar.size() ;
1125//     vector< vector<float> > tcal( numReceptor, vector<float>( 1, 1.0 ) ) ;
1126//     for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1127//       tcal[ir][0] = tcalScalar[ir][0].get() ;
1128//     return tcal ;
1129//   }
1130  else {
1131    return defaultTcal ;
1132  }
1133}
1134
1135vector<float> ASDMReader::getOpacity( unsigned int idx )
1136{
1137  vector<float> tau(0) ;
1138  CalAtmosphereTable &atmtab = asdm_->getCalAtmosphere() ;
1139  unsigned int nrow = atmtab.size() ;
1140  if ( nrow > 0 ) {
1141    unsigned int index = dataIdList_[idx] ;
1142    //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1143    double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1144    ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1145    //int feedid = vmsData_->v_feedId1[index] ;
1146    //Tag atag( antennaId_, TagType::Antenna ) ;
1147    //Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1148    //DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey(ddtag) ;
1149    //Tag spwtag = ddrow->getSpectralWindowId() ;
1150    //SpectralWindowRow *spwrow = ddrow->getSpectralWindowUsingSpectralWindowId() ;
1151    //BasebandName bbname = spwrow->getBasebandName() ;
1152    //FeedRow *frow = asdm_->getFeed().getRowByKey( atag, spwtag, tint, feedid ) ;
1153    //int nfeed = frow->getNumReceptor() ;
1154    //vector<ReceiverRow *> rrows = frow->getReceivers() ;
1155    vector<CalAtmosphereRow *> atmrows = atmtab.get() ;
1156    //ReceiverBand rb = rrows[0]->getFrequencyBand() ;
1157    int row0 = -1 ;
1158    double eps = tint.getStart().getMJD() ;
1159    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1160      CalAtmosphereRow *atmrow = atmrows[irow] ;
1161      if ( casa::String(atmrow->getAntennaName()) != antennaName_
1162           //|| atmrow->getReceiverBand() != rb
1163           //|| atmrow->getBasebandName() != bbname
1164           || atmrow->getCalDataUsingCalDataId()->getCalType() != CAL_ATMOSPHERE )
1165        continue ;
1166      else {
1167        double dt = tint.getStart().getMJD() - atmrow->getEndValidTime().getMJD() ;
1168        if ( dt >= 0 && dt < eps ) {
1169          eps = dt ;
1170          row0 = (int)irow ;
1171        }
1172      }
1173    }
1174    if ( row0 != -1 ) {
1175      CalAtmosphereRow *atmrow = atmrows[row0] ;
1176      tau = atmrow->getTau() ;
1177    }
1178  }
1179  return tau ;
1180}
1181
1182void ASDMReader::getWeatherInfo( unsigned int idx,
1183                                 float &temperature,
1184                                 float &pressure,
1185                                 float &humidity,
1186                                 float &windspeed,
1187                                 float &windaz )
1188{
1189  temperature = 0.0 ;
1190  pressure = 0.0 ;
1191  humidity = 0.0 ;
1192  windspeed = 0.0 ;
1193  windaz = 0.0 ;
1194
1195  //logsink_->postLocally( LogMessage("weatherStationId_ = "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1196
1197  WeatherTable &wtab = asdm_->getWeather() ;
1198  if ( wtab.size() == 0 || weatherStationId_ == -1 )
1199    return ;
1200
1201  unsigned int index = dataIdList_[idx] ;
1202  //Tag anttag( antennaId_, TagType::Antenna ) ;
1203  //Tag sttag = (asdm_->getAntenna().getRowByKey( anttag ))->getStationId() ;
1204  Tag sttag( (unsigned int)weatherStationId_, TagType::Station ) ;
1205  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1206  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1207  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1208  //WeatherRow *wrow = wtab.getRowByKey( sttag, tint ) ;
1209  vector<WeatherRow *> *wrows = wtab.getByContext( sttag ) ;
1210  WeatherRow *wrow = (*wrows)[0] ;
1211  unsigned int nrow = wrows->size() ;
1212  //logsink_->postLocally( LogMessage("There are "+String::toString(nrow)+" rows for given context: stationId "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1213  ArrayTime startTime = tint.getStart() ;
1214  if ( startTime < (*wrows)[0]->getTimeInterval().getStart() ) {
1215    temperature = (*wrows)[0]->getTemperature().get() ;
1216    pressure = (*wrows)[0]->getPressure().get() ;
1217    humidity = (*wrows)[0]->getRelHumidity().get() ;
1218    windspeed = (*wrows)[0]->getWindSpeed().get() ;
1219    windaz = (*wrows)[0]->getWindDirection().get() ;
1220  }
1221  else if ( startTime > (*wrows)[nrow-1]->getTimeInterval().getStart() ) {
1222    temperature = (*wrows)[nrow-1]->getTemperature().get() ;
1223    pressure = (*wrows)[nrow-1]->getPressure().get() ;
1224    humidity = (*wrows)[nrow-1]->getRelHumidity().get() ;
1225    windspeed = (*wrows)[nrow-1]->getWindSpeed().get() ;
1226    windaz = (*wrows)[nrow-1]->getWindDirection().get() ;
1227  }
1228  else {
1229    for ( unsigned int irow = 1 ; irow < wrows->size() ; irow++ ) {
1230      wrow = (*wrows)[irow-1] ;
1231      if ( startTime < (*wrows)[irow]->getTimeInterval().getStart() ) {
1232        temperature = wrow->getTemperature().get() ;
1233        pressure = wrow->getPressure().get() ;
1234        humidity = wrow->getRelHumidity().get() ;
1235        windspeed = wrow->getWindSpeed().get() ;
1236        windaz = wrow->getWindDirection().get() ;
1237        break ;
1238      }
1239    }
1240  }
1241
1242  return ;
1243}
1244
1245void ASDMReader::processStation()
1246{
1247  antennaPad_.resize(0) ;
1248  weatherStation_.resize(0) ;
1249  StationTable &stab = asdm_->getStation() ;
1250  vector<StationRow *> srows = stab.get() ;
1251  for ( unsigned int irow = 0 ; irow < srows.size() ; irow++ ) {
1252    StationType stype = srows[irow]->getType() ;
1253    Tag stag = srows[irow]->getStationId() ;
1254    if ( stype == ANTENNA_PAD )
1255      antennaPad_.push_back( stag ) ;
1256    else if ( stype == WEATHER_STATION )
1257      weatherStation_.push_back( stag ) ;
1258  }
1259
1260   weatherStationId_ = getClosestWeatherStation() ;
1261}
1262
1263int ASDMReader::getClosestWeatherStation()
1264{
1265  if ( weatherStation_.size() == 0 )
1266    return -1 ;
1267
1268  Tag atag( antennaId_, TagType::Antenna ) ;
1269  Tag stag = (asdm_->getAntenna().getRowByKey( atag ))->getStationId() ;
1270  vector<double> apos( 3 ) ;
1271  StationTable &stab = asdm_->getStation() ;
1272  StationRow *srow = stab.getRowByKey( stag ) ;
1273  vector<Length> pos = srow->getPosition() ;
1274  apos[0] = pos[0].get() ;
1275  apos[1] = pos[1].get() ;
1276  apos[2] = pos[2].get() ;
1277 
1278  double eps = 1.0e20 ;
1279  int retval = -1 ;
1280  for ( unsigned int ir = 0 ; ir < weatherStation_.size() ; ir++ ) {
1281    srow = stab.getRowByKey( weatherStation_[ir] ) ;
1282    vector<Length> wpos = srow->getPosition() ;
1283    double dist = (apos[0]-wpos[0].get())*(apos[0]-wpos[0].get())
1284      + (apos[1]-wpos[1].get())*(apos[1]-wpos[1].get())
1285      + (apos[2]-wpos[2].get())*(apos[2]-wpos[2].get()) ;
1286    if ( dist < eps ) {
1287      retval = (int)(weatherStation_[ir].getTagValue()) ;
1288    }
1289  }
1290
1291  return retval ;
1292}
1293
1294void ASDMReader::getPointingInfo( unsigned int idx,
1295                                  vector<double> &dir,
1296                                  double &az,
1297                                  double &el,
1298                                  vector<double> &srate )
1299{
1300  String funcName = "getPointingInfo" ;
1301
1302  dir.resize(0) ;
1303  az = -1.0 ;
1304  el = -1.0 ;
1305  srate.resize(0) ;
1306
1307  Tag atag( antennaId_, TagType::Antenna ) ;
1308  unsigned int index = dataIdList_[idx] ;
1309  vector<PointingRow *> *prows = asdm_->getPointing().getByContext( atag ) ;
1310  PointingRow *prow ;
1311  PointingRow *qrow ;
1312  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1313  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1314  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1315
1316  unsigned int nrow = prows->size() ;
1317  //logsink_->postLocally( LogMessage("There are " << nrow << " rows for given context: antennaId "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
1318
1319//   for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1320//     prow = (*prows)[irow] ;
1321//     ArrayTimeInterval ati = prow->getTimeInterval() ;
1322//     ArrayTime pst = ati.getStart() ;
1323//     ArrayTime pet( ati.getStartInMJD()+ati.getDurationInDays() ) ;
1324//     logsink_->postLocally( LogMessage("start: "+pst.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1325//     logsink_->postLocally( LogMessage("end: "+pet.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1326//   }
1327 
1328  if ( nrow == 0 )
1329    return ;
1330
1331  srate.resize(2) ;
1332  srate[0] = 0.0 ;
1333  srate[1] = 0.0 ;
1334  az = 0.0 ;
1335  el = 0.0 ;
1336  double tcen = 0.0 ;
1337
1338  //
1339  // shape of pointingDirection is (numSample,2) if usePolynomial = False, while it is
1340  // (numTerm,2) if usePolynomial = True.
1341  //
1342  // In the former case, typical sampled time interval is 48msec, which is very small
1343  // compared with typical integration time (~a few sec).
1344  // Scan rate for this case is always [0,0] (or get slope?).
1345  //
1346  // In the later case, scan rate is (pointingDirection[1][0],pointingDirection[1][1])
1347  //
1348  // PointingDirection seems to store AZEL
1349  //
1350  ArrayTimeInterval pTime0 = (*prows)[0]->getTimeInterval() ;
1351  ArrayTimeInterval pTime1 = (*prows)[nrow-1]->getTimeInterval() ;
1352  if ( tint.getStartInMJD()+tint.getDurationInDays() < pTime0.getStartInMJD() ) {
1353    logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint < ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1354    prow = (*prows)[0] ;
1355    //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1356    vector< vector<Angle> > dirA = prow->getTarget() ;
1357    vector< vector<Angle> > offA = prow->getOffset() ;
1358    //az = dirA[0][0].get() ;
1359    //el = dirA[0][1].get() ;
1360    az = dirA[0][0].get() + offA[0][0].get() ;
1361    el = dirA[0][1].get() + offA[0][1].get() ;
1362    if ( prow->getUsePolynomials() && prow->getNumTerm() > 1 ) {
1363      //srate[0] = dirA[1][0].get() ;
1364      //srate[1] = dirA[1][1].get() ;
1365      srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1366      srate[1] = dirA[1][1].get() + offA[1][1].get() ;
1367    }     
1368  }
1369  else if ( tint.getStartInMJD() > pTime1.getStartInMJD()+pTime1.getDurationInDays() ) {
1370    logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint > ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1371    prow = (*prows)[nrow-1] ;
1372    int numSample = prow->getNumSample() ;
1373    //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1374    vector< vector<Angle> > dirA = prow->getTarget() ;
1375    vector< vector<Angle> > offA = prow->getOffset() ;
1376    if ( prow->getUsePolynomials() ) {
1377      //az = dirA[0][0].get() ;
1378      //el = dirA[0][1].get() ;
1379      az = dirA[0][0].get() + offA[0][0].get() ;
1380      el = dirA[0][1].get() + offA[0][1].get() ;
1381      if ( prow->getNumTerm() > 1 ) {
1382        //srate[0] = dirA[1][0].get() ;
1383        //srate[1] = dirA[1][1].get() ;
1384        srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1385        srate[1] = dirA[1][1].get() + offA[1][0].get() ;
1386      }
1387    }
1388    else {
1389      //az = dirA[numSample-1][0].get() ;
1390      //el = dirA[numSample-1][1].get() ;
1391      az = dirA[numSample-1][0].get() + offA[numSample-1][0].get() ;
1392      el = dirA[numSample-1][1].get() + offA[numSample-1][1].get() ;
1393    }
1394  }
1395  else {
1396    ArrayTime startTime = tint.getStart() ;
1397    ArrayTime endTime( tint.getStartInMJD()+tint.getDurationInDays() ) ;
1398    int row0 = -1 ;
1399    int row1 = -1 ;
1400    int row2 = -1 ;
1401    double dt0 = getMidTime( tint ).getMJD() ;
1402    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1403      prow = (*prows)[irow] ;
1404      double dt = getMidTime( tint ).getMJD() - getMidTime( prow->getTimeInterval() ).getMJD() ;
1405      if ( dt > 0 && dt < dt0 ) {
1406        dt0 = dt ;
1407        row2 = irow ;
1408      }
1409      if ( prow->getTimeInterval().contains( startTime ) )
1410        row0 = irow ;
1411      else if ( prow->getTimeInterval().contains( endTime ) )
1412        row1 = irow ;
1413      if ( row0 != -1 && row1 != -1 )
1414        break ;
1415    }
1416    //logsink_->postLocally( LogMessage("row0 = "+String::toString(row0)+", row1 = "+String::toString(row1)+", row2 = "+String::toString(row2),LogOrigin(className_,funcName,WHERE)) ) ;
1417    if ( row0 == -1 && row1 == -1 ) {
1418      prow = (*prows)[row2] ;
1419      qrow = (*prows)[row2+1] ;
1420      double t0 = getEndTime( prow->getTimeInterval() ).getMJD() ;
1421      double t1 = qrow->getTimeInterval().getStartInMJD() ;
1422      double t = startTime.getMJD() ;
1423      //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1424      //vector< vector<Angle> > dirB = qrow->getPointingDirection() ;
1425      vector< vector<Angle> > dirA = prow->getTarget() ;
1426      vector< vector<Angle> > offA = prow->getOffset() ;
1427      vector< vector<Angle> > dirB = qrow->getTarget() ;
1428      vector< vector<Angle> > offB = qrow->getOffset() ;
1429      //double da0 = dirA[0][0].get() ;
1430      //double db0 = dirB[0][0].get() ;
1431      //double da1 = dirA[0][1].get() ;
1432      //double db1 = dirB[0][1].get() ;
1433      double da0 = dirA[0][0].get() + offA[0][0].get() ;
1434      double db0 = dirB[0][0].get() + offB[0][0].get() ;
1435      double da1 = dirA[0][1].get() + offA[0][1].get() ;
1436      double db1 = dirB[0][1].get() + offB[0][1].get() ;
1437      if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1438        double dt = ( t - t0 ) / ( t1 - t0 ) ;
1439        az = da0 + ( db0 - da0 ) * dt ;
1440        el = da1 + ( db1 - da1 ) * dt ;
1441        if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1442          //double ra0 = dirA[1][0].get() ;
1443          //double rb0 = dirB[1][0].get() ;
1444          //double ra1 = dirA[1][1].get() ;
1445          //double rb1 = dirB[1][1].get() ;
1446          double ra0 = dirA[1][0].get() + offA[1][0].get() ;
1447          double rb0 = dirB[1][0].get() + offB[1][0].get() ;
1448          double ra1 = dirA[1][1].get() + offA[1][1].get() ;
1449          double rb1 = dirB[1][1].get() + offB[1][1].get() ;         
1450          srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1451          srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1452        }
1453      }
1454      else if ( !(qrow->getUsePolynomials()) ) {
1455        double dt = ( t - t0 ) / ( t1 - t0 ) ;
1456        az = da0 + ( db0 - da0 ) * dt ;
1457        el = da1 + ( db1 - da1 ) * dt ;
1458      }
1459      //else {
1460      // nothing to do
1461      //}
1462    }
1463    else {
1464      int ndir = 0 ;
1465      if ( row0 == -1 ) {
1466        row0 = row1 ;
1467      }
1468      else if ( row1 == -1 ) {
1469        row1 = row0 ;
1470      }
1471      prow = (*prows)[row0] ;
1472      qrow = (*prows)[row1] ;
1473      if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1474        //logsink_->postLocally( LogMessage("usePolynomial = True",LogOrigin(className_,funcName,WHERE)) ) ;
1475        if ( row0 == row1 ) {
1476          prow = (*prows)[row0] ;
1477          //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1478          vector< vector<Angle> > dirA = prow->getTarget() ;
1479          vector< vector<Angle> > offA = prow->getOffset() ;
1480          //az = dirA[0][0].get() ;
1481          //el = dirA[0][1].get() ;
1482          az = dirA[0][0].get() + offA[0][0].get() ;
1483          el = dirA[0][1].get() + offA[0][1].get() ;
1484          if ( prow->getNumTerm() > 1 ) {
1485            //srate[0] = dirA[1][0].get() ;
1486            //srate[1] = dirA[1][1].get() ;
1487            srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1488            srate[1] = dirA[1][1].get() + offA[1][0].get() ;
1489          }
1490        }
1491        else {
1492          prow = (*prows)[row0] ;
1493          qrow = (*prows)[row1] ;
1494          //vector< vector<Angle> > dirA = qrow->getPointingDirection() ;
1495          //vector< vector<Angle> > dirB = qrow->getPointingDirection() ;
1496          vector< vector<Angle> > dirA = qrow->getTarget() ;
1497          vector< vector<Angle> > offA = prow->getOffset() ;
1498          vector< vector<Angle> > dirB = qrow->getTarget() ;
1499          vector< vector<Angle> > offB = qrow->getOffset() ;
1500          double t0 = getMidTime( prow->getTimeInterval() ).getMJD() ;
1501          double t1 = getMidTime( qrow->getTimeInterval() ).getMJD() ;
1502          double t = startTime.getMJD() ;
1503          double dt = ( t - t0 ) / ( t1 - t0 ) ;
1504          //double da0 = dirA[0][0].get() ;
1505          //double db0 = dirB[0][0].get() ;
1506          //double da1 = dirA[0][1].get() ;
1507          //double db1 = dirB[0][1].get() ;
1508          double da0 = dirA[0][0].get() + offA[0][0].get() ;
1509          double db0 = dirB[0][0].get() + offB[0][0].get() ;
1510          double da1 = dirA[0][1].get() + offA[0][1].get() ;
1511          double db1 = dirB[0][1].get() + offB[0][1].get() ;
1512          az = da0 + ( db0 - da0 ) * dt ;
1513          el = da1 + ( db1 - db0 ) * dt ;
1514          if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1515            //double ra0 = dirA[1][0].get() ;
1516            //double rb0 = dirB[1][0].get() ;
1517            //double ra1 = dirA[1][1].get() ;
1518            //double rb1 = dirB[1][1].get() ;
1519            double ra0 = dirA[1][0].get() + offA[1][0].get() ;
1520            double rb0 = dirB[1][0].get() + offB[1][0].get() ;
1521            double ra1 = dirA[1][1].get() + offA[1][1].get() ;
1522            double rb1 = dirB[1][1].get() + offB[1][1].get() ;         
1523            srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1524            srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1525          }
1526        }
1527      }
1528      else if ( prow->getUsePolynomials() == qrow->getUsePolynomials() ) {
1529        //logsink_->postLocally( LogMessage("numSample == numTerm",LogOrigin(className_,funcName,WHERE)) ) ;
1530        for ( int irow = row0 ; irow <= row1 ; irow++ ) {
1531          prow = (*prows)[irow] ;
1532          int numSample = prow->getNumSample() ;
1533          //logsink_->postLocally( LogMessage("numSample = "+String::toString(numSample),LogOrigin(className_,funcName,WHERE)) ) ;
1534          //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1535          vector< vector<Angle> > dirA = prow->getTarget() ;
1536          vector< vector<Angle> > offA = prow->getOffset() ;
1537          if ( prow->isSampledTimeIntervalExists() ) {
1538            //logsink_->postLocally( LogMessage("sampledTimeIntervalExists",LogOrigin(className_,funcName,WHERE)) ) ;
1539            vector<ArrayTimeInterval> stime = prow->getSampledTimeInterval() ;
1540            for ( int isam = 0 ; isam < numSample ; isam++ ) {
1541              //if ( tint.overlaps( stime[isam] ) ) {
1542              if ( tint.contains( stime[isam] ) ) {
1543                //az += dirA[isam][0].get() ;
1544                //el += dirA[isam][1].get() ;
1545                az += dirA[isam][0].get() + offA[isam][0].get() ;
1546                el += dirA[isam][1].get() + offA[isam][1].get() ;
1547                ndir++ ;
1548              }
1549            }
1550          }
1551          else {
1552            double sampleStart = prow->getTimeInterval().getStartInMJD() ;
1553            double sampleInterval = prow->getTimeInterval().getDurationInDays() / (double)numSample ;
1554            //logsink_->postLocally( LogMessage("sampleStart = "+String::toString(sampleStart),LogOrigin(className_,funcName,WHERE)) )
1555            //logsink_->postLocally( LogMessage("sampleInterval = "+String::toString(sampleInterval),LogOrigin(className_,funcName,WHERE)) ) ;
1556            //logsink_->postLocally( LogMessage("tint = "+tint.toString(),LogOrigin(className_,funcName,WHERE)) ) ;
1557            for ( int isam = 0 ; isam < numSample ; isam++ ) {
1558              ArrayTimeInterval stime( sampleStart+isam*sampleInterval, sampleInterval ) ;
1559              //if ( tint.overlaps( stime ) ) {
1560              if ( tint.contains( stime ) ) {
1561                //az += dirA[isam][0].get() ;
1562                //el += dirA[isam][1].get() ;
1563                az += dirA[isam][0].get() + offA[isam][0].get() ;
1564                el += dirA[isam][1].get() + offA[isam][1].get() ;
1565                tcen += getMidTime( stime ).getMJD() ;
1566                ndir++ ;
1567              }
1568            }
1569          }
1570        }
1571        if ( ndir > 1 ) {
1572          az /= (double)ndir ;
1573          el /= (double)ndir ;
1574          tcen /= (double)ndir ;
1575        }
1576      }
1577      //else {
1578      // nothing to do
1579      //}
1580    }
1581   
1582    AntennaRow *arow = asdm_->getAntenna().getRowByKey( Tag( antennaId_, TagType::Antenna ) ) ;
1583    StationRow *srow = arow->getStationUsingStationId() ;
1584    vector<Length> antposL = srow->getPosition() ;
1585    casa::Vector<casa::Double> antpos( 3 ) ;
1586    for ( int i = 0 ; i < 3 ; i++ )
1587      antpos[i] = antposL[i].get() ;
1588    //logsink_->postLocally( LogMessage("tcen = "+String::toString(tcen),LogOrigin(className_,funcName,WHERE)) ) ;
1589    //logsink_->postLocally( LogMessage("antpos = "+String::toString(antpos),LogOrigin(className_,funcName,WHERE)) ) ;
1590    toJ2000( dir, az, el, tcen, antpos ) ;
1591
1592  }
1593
1594  return ;
1595}
1596
1597ArrayTime ASDMReader::getMidTime( const ArrayTimeInterval &t )
1598{
1599  return ArrayTime( t.getStartInMJD() + 0.5 * t.getDurationInDays() ) ;
1600}
1601
1602ArrayTime ASDMReader::getEndTime( const ArrayTimeInterval &t )
1603{
1604  return ArrayTime( t.getStartInMJD() + t.getDurationInDays() ) ;
1605}
1606
1607ArrayTime ASDMReader::getStartTime( const ArrayTimeInterval &t )
1608{
1609  return ArrayTime( t.getStartInMJD() ) ;
1610}
1611
1612void ASDMReader::toJ2000( vector<double> &dir,
1613                          double az,
1614                          double el,
1615                          double mjd,
1616                          casa::Vector<casa::Double> antpos )
1617{
1618  String funcName = "toJ2000" ;
1619
1620  casa::Vector<casa::Double> azel( 2 ) ;
1621  azel[0] = az ;
1622  azel[1] = el ;
1623  casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
1624  casa::Vector<casa::Quantity> qantpos( 3 ) ;
1625  qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
1626  qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
1627  qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
1628  casa::MPosition mp( casa::MVPosition( qantpos ),
1629                      casa::MPosition::ITRF ) ;
1630  //ostringstream oss ;
1631  //mp.print( oss ) ;
1632  //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
1633
1634  casa::MeasFrame mf( me, mp ) ;
1635  casa::MDirection::Convert toj2000( casa::MDirection::AZELGEO,
1636  //MDirection::Convert toj2000( MDirection::AZEL,
1637  //MDirection::Convert toj2000( MDirection::AZELSW,
1638  //MDirection::Convert toj2000( MDirection::AZELSWGEO,
1639  //MDirection::Convert toj2000( MDirection::AZELNE,
1640  //MDirection::Convert toj2000( MDirection::AZELNEGEO,
1641                                     casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
1642  casa::Vector<casa::Double> cdir = toj2000( azel ).getAngle( "rad" ).getValue() ;
1643  //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
1644  dir.resize(2) ;
1645  dir[0] = (double)(cdir[0]) ;
1646  dir[1] = (double)(cdir[1]) ;
1647}
1648
1649void ASDMReader::setLogger( CountedPtr<LogSinkInterface> &logsink )
1650{
1651  logsink_ = logsink ;
1652}
1653
1654string ASDMReader::getFrame()
1655{
1656  String funcName = "getFrame" ;
1657 
1658  // default is TOPO
1659  string frame = "TOPO" ;
1660
1661  SpectralWindowTable &spwtab = asdm_->getSpectralWindow() ;
1662  vector<SpectralWindowRow *> rows = spwtab.get() ;
1663  vector<FrequencyReferenceCode> measFreqRef( rows.size() ) ;
1664  int nref = 0 ;
1665  for ( unsigned int irow = 0 ; irow < rows.size() ; irow++ ) {
1666    int nchan = rows[irow]->getNumChan() ;
1667    if ( nchan != 4 ) {
1668      if ( rows[irow]->isMeasFreqRefExists() ) {
1669        measFreqRef[nref] = rows[irow]->getMeasFreqRef() ;
1670        nref++ ;
1671      }
1672    }
1673  }
1674  if ( nref != 0 ) {
1675    frame = CFrequencyReferenceCode::toString( measFreqRef[0] ) ;
1676  }
1677 
1678  //logsink_->postLocally( LogMessage("frame = "+String::toString(frame),LogOrigin(className_,funcName,WHERE)) ) ;
1679
1680  return frame ;
1681}
1682
1683int ASDMReader::getNumIFs()
1684{
1685  String funcName = "getNumIFs" ;
1686
1687  int nif = 0 ;
1688  vector<SpectralWindowRow *> rows = asdm_->getSpectralWindow().get() ;
1689  unsigned int nrow = rows.size() ;
1690  // check if all rows have freqGroup attribute
1691  bool freqGroupExists = true ;
1692  bool countedWvr = false ;
1693  for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1694    freqGroupExists &= rows[irow]->isFreqGroupExists() ;
1695    if ( rows[irow]->getNumChan() == 4 ) {
1696      if ( !countedWvr ) {
1697        countedWvr = true ;
1698        nif++ ;
1699      }
1700    }
1701    else {
1702      nif++ ;
1703    }
1704  }
1705 
1706  if ( freqGroupExists ) {
1707    vector<int> freqGroup(0) ;
1708    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1709      int fg = rows[irow]->getFreqGroup() ;
1710      if ( (int)count( freqGroup.begin(), freqGroup.end(), fg ) == 0 ) {
1711        freqGroup.push_back( fg ) ;
1712      }
1713    }
1714    nif = freqGroup.size() ;
1715  }
1716
1717  //logsink_->postLocally( LogMessage("nif = "+String::toString(nif),LogOrigin(className_,funcName,WHERE)) ) ;
1718
1719  return nif ;
1720}
Note: See TracBrowser for help on using the repository browser.