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

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

Frequency is stored as LSRK value although raw frequency value is
in arbitrary frequency frame.
Also, direction is stored as J2000 value independently of original
direction reference code.
File size: 67.0 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    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  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  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  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//   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  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{
709  String funcName = "getFrequency" ;
710
711  Tag ddTag( vmsData_->v_dataDescId[dataIdList_[idx]], TagType::DataDescription ) ;
712  DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey( ddTag ) ;
713  //Tag spwid = ddrow->getSpectralWindowId() ;
714  SpectralWindowRow *spwrow = ddrow->getSpectralWindowUsingSpectralWindowId() ;
715  int nchan = spwrow->getNumChan() ;
716  if ( nchan == 1 ) {
717    //logsink_->postLocally( LogMessage("channel averaged data",LogOrigin(className_,funcName,WHERE)) ) ;
718    refpix = 0.0 ;
719    incr = spwrow->getTotBandwidth().get() ;
720    if ( spwrow->isChanFreqStartExists() ) {
721      refval = spwrow->getChanFreqStart().get() ;
722    }
723    else if ( spwrow->isChanFreqArrayExists() ) {
724      refval = spwrow->getChanFreqArray()[0].get() ;
725    }
726    else {
727      throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
728    }     
729  }
730  else if ( nchan % 2 ) {
731    // odd
732    //logsink_->postLocally( LogMessage("odd case",LogOrigin(className_,funcName,WHERE)) ) ;
733    refpix = 0.5 * ( (double)nchan - 1.0 ) ;
734    int ic = ( nchan - 1 ) / 2 ;
735    if ( spwrow->isChanWidthExists() ) {
736      incr = spwrow->getChanWidth().get() ;
737    }
738    else if ( spwrow->isChanWidthArrayExists() ) {
739      incr = spwrow->getChanWidthArray()[0].get() ;
740    }
741    else {
742      throw (casa::AipsError( "Either chanWidthArray or chanWidth must exist." )) ;
743    }     
744    if ( spwrow->isChanFreqStepExists() ) {
745      if ( spwrow->getChanFreqStep().get() < 0.0 )
746        incr *= -1.0 ;
747    }
748    else if ( spwrow->isChanFreqArrayExists() ) {
749      vector<Frequency> chanFreqArr = spwrow->getChanFreqArray() ;
750      if ( chanFreqArr[0].get() > chanFreqArr[1].get() )
751        incr *= -1.0 ;
752    }
753    else {
754      throw (casa::AipsError( "Either chanFreqArray or chanFreqStep must exist." )) ;
755    }         
756    if ( spwrow->isChanFreqStartExists() ) {
757      refval = spwrow->getChanFreqStart().get() + refpix * incr ;
758    }
759    else if ( spwrow->isChanFreqArrayExists() ) {
760      refval = spwrow->getChanFreqArray()[ic].get() ;
761    }
762    else {
763      throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
764    }     
765  }
766  else {
767    // even
768    //logsink_->postLocally( LogMessage("even case",LogOrigin(className_,funcName,WHERE)) ) ;
769    refpix = 0.5 * ( (double)nchan - 1.0 ) ;
770    int ic = nchan / 2 ;
771    if ( spwrow->isChanWidthExists() ) {
772      incr = spwrow->getChanWidth().get() ;
773    }
774    else if ( spwrow->isChanWidthArrayExists() ) {
775      incr = spwrow->getChanWidthArray()[0].get() ;
776    }
777    else {
778      throw (casa::AipsError( "Either chanWidthArray or chanWidth must exist." )) ;
779    }
780    if ( spwrow->isChanFreqStepExists() ) {
781      if ( spwrow->getChanFreqStep().get() < 0.0 )
782        incr *= -1.0 ;
783    }
784    else if ( spwrow->isChanFreqArrayExists() ) {
785      vector<Frequency> chanFreqArr = spwrow->getChanFreqArray() ;
786      if ( chanFreqArr[0].get() > chanFreqArr[1].get() )
787        incr *= -1.0 ;
788    }
789    else {
790      throw (casa::AipsError( "Either chanFreqArray or chanFreqStep must exist." )) ;
791    }         
792    if ( spwrow->isChanFreqStartExists() ) {
793      refval = spwrow->getChanFreqStart().get() + refpix * incr ;
794    }
795    else if ( spwrow->isChanFreqArrayExists() ) {
796      vector<Frequency> freqs = spwrow->getChanFreqArray() ;
797      refval = 0.5 * ( freqs[ic-1].get() + freqs[ic].get() ) ;     
798    }
799    else {
800      throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
801    }     
802  }
803}
804
805vector<double> ASDMReader::getRestFrequency( unsigned int idx )
806{
807  vector<double> rf( 0 ) ;
808  unsigned int index = dataIdList_[idx] ;
809  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
810  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
811  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
812  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
813  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
814  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
815  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
816  if ( frow->isSourceIdExists() ) {
817    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
818    int sid = frow->getSourceId() ;
819    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
820    if ( srow->isRestFrequencyExists() ) {
821      //logsink_->postLocally( LogMessage("restFrequency exists",LogOrigin(className_,funcName,WHERE)) ) ;
822      vector<Frequency> restfreq = srow->getRestFrequency() ;
823      rf.resize( restfreq.size() ) ;
824      for ( unsigned int i = 0 ; i < restfreq.size() ; i++ )
825        rf[i] = restfreq[i].get() ;
826    }
827  }
828  return rf ;
829}
830
831double ASDMReader::getTime( unsigned int idx )
832{
833  double tsec = vmsData_->v_time[dataIdList_[idx]] ;
834  return tsec * s2d ;
835}
836
837double ASDMReader::getInterval( unsigned int idx )
838{
839  return vmsData_->v_interval[dataIdList_[idx]] ;
840}
841
842string ASDMReader::getSourceName( unsigned int idx )
843{
844  unsigned int index = dataIdList_[idx] ;
845  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
846  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
847  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
848  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
849  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
850  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
851  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
852  string srcname ;
853  if ( frow->isSourceIdExists() ) {
854    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
855    int sid = frow->getSourceId() ;
856    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
857    srcname = srow->getSourceName() ;
858  }
859  else {
860    srcname = frow->getFieldName() ;
861  }
862  return srcname ;
863}
864
865string ASDMReader::getFieldName( unsigned int idx )
866{
867  int fid = vmsData_->v_fieldId[dataIdList_[idx]] ;
868  Tag ftag( fid, TagType::Field ) ;
869  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
870  ostringstream oss ;
871  oss << frow->getFieldName() << "__" << fid ;
872  return oss.str() ;
873}
874
875int ASDMReader::getSrcType( unsigned int scan,
876                            unsigned int subscan )
877{
878  int srctype = SrcType::NOTYPE ;
879  Tag ebtag = mainRow_[row_]->getExecBlockId() ;
880  ScanRow *scanrow = asdm_->getScan().getRowByKey( ebtag, (int)scan ) ;
881  ScanIntent scanIntent = scanrow->getScanIntent()[0] ;
882  SubscanRow *subrow = asdm_->getSubscan().getRowByKey( ebtag, (int)scan, (int)subscan ) ;
883  SubscanIntent subIntent = subrow->getSubscanIntent() ;
884  SwitchingMode swmode = subrow->getSubscanMode() ;
885  if ( scanIntent == OBSERVE_TARGET ) {
886    // on sky scan
887    if ( swmode == NO_SWITCHING || swmode == POSITION_SWITCHING ) {
888      // position switching
889      // tentatively set NO_SWITCHING = POSITION_SWITCHING
890      if ( subIntent == ON_SOURCE ) {
891        srctype = SrcType::PSON ;
892      }
893      else if ( subIntent == OFF_SOURCE ) {
894        srctype = SrcType::PSOFF ;
895      }
896    }
897    else if ( swmode == FREQUENCY_SWITCHING ) {
898      // frequency switching
899      if ( subIntent == ON_SOURCE ) {
900        srctype = SrcType::FSON ;
901      }
902      else if ( subIntent == OFF_SOURCE ) {
903        srctype = SrcType::FSOFF ;
904      }
905    }
906    else if ( swmode == NUTATOR_SWITCHING ) {
907      // nutator switching
908      if ( subIntent == ON_SOURCE ) {
909        srctype = SrcType::PSON ;
910      }
911      else if ( subIntent == OFF_SOURCE ) {
912        srctype = SrcType::PSOFF ;
913      }
914    }
915    else {
916      // other switching mode
917      if ( subIntent == ON_SOURCE ) {
918        srctype = SrcType::SIG ;
919      }
920      else if ( subIntent == OFF_SOURCE ) {
921        srctype = SrcType::REF ;
922      }
923    }
924  }
925  else {
926    // off sky (e.g. calibrator device) scan
927    if ( subIntent == ON_SOURCE ) {
928      srctype = SrcType::SIG ;
929    }
930    else if ( subIntent == OFF_SOURCE ) {
931      srctype = SrcType::REF ;
932    }
933    else if ( subIntent == HOT ) {
934      srctype = SrcType::HOT ;
935    }
936    else if ( subIntent == AMBIENT ) {
937      srctype = SrcType::SKY ;
938    }
939    else {
940      srctype = SrcType::CAL ;
941    }
942  }
943
944  return srctype ;
945}
946
947unsigned int ASDMReader::getSubscanNo( unsigned int idx )
948{
949  //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)) ) ;
950  return vmsData_->v_msState[dataIdList_[idx]].subscanNum ;
951}
952
953vector<double> ASDMReader::getSourceDirection( unsigned int idx )
954{
955  vector<double> dir( 2, 0.0 ) ;
956  unsigned int index = dataIdList_[idx] ;
957  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
958  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
959  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
960  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
961  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
962  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
963  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
964  string srcname ;
965  if ( frow->isSourceIdExists() ) {
966    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
967    int sid = frow->getSourceId() ;
968    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
969    vector<Angle> srcdir = srow->getDirection() ;
970    if ( srow->isDirectionCodeExists() ) {
971      DirectionReferenceCode dircode = srow->getDirectionCode() ;
972      if ( dircode != J2000 ) {
973        // if not J2000, need direction conversion
974      }
975    }
976    dir[0] = srcdir[0].get() ;
977    dir[1] = srcdir[1].get() ;
978  }
979  return dir ;
980}
981
982void ASDMReader::getSourceDirection( unsigned int idx,
983                                     vector<double> &dir,
984                                     string &ref )
985{
986  dir.resize( 2 ) ;
987  dir[0] = 0.0 ;
988  dir[1] = 0.0 ;
989  ref = "J2000" ;
990  unsigned int index = dataIdList_[idx] ;
991  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
992  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
993  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
994  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
995  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
996  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
997  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
998  string srcname ;
999  if ( frow->isSourceIdExists() ) {
1000    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1001    int sid = frow->getSourceId() ;
1002    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1003    vector<Angle> srcdir = srow->getDirection() ;
1004    if ( srow->isDirectionCodeExists() ) {
1005      ref = CDirectionReferenceCode::toString( srow->getDirectionCode() ) ;
1006    }
1007    dir[0] = srcdir[0].get() ;
1008    dir[1] = srcdir[1].get() ;
1009  }
1010}
1011
1012void ASDMReader::getSourceDirection( vector<double> &dir, string &ref )
1013{
1014  dir.resize( 2 ) ;
1015  ref = "J2000" ; // default is J2000
1016  SourceTable &tab = asdm_->getSource() ;
1017  SourceRow *row = tab.get()[0] ;
1018  vector<Angle> dirA = row->getDirection() ;
1019  dir[0] = dirA[0].get() ;
1020  dir[1] = dirA[1].get() ;
1021  if ( row->isDirectionCodeExists() ) {
1022    ref = CDirectionReferenceCode::toString( row->getDirectionCode() ) ;
1023  }
1024}
1025
1026vector<double> ASDMReader::getSourceProperMotion( unsigned int idx )
1027{
1028  vector<double> pm( 2, 0.0 ) ;
1029  unsigned int index = dataIdList_[idx] ;
1030  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1031  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1032  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1033  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1034  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1035  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1036  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1037  string srcname ;
1038  if ( frow->isSourceIdExists() ) {
1039    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1040    int sid = frow->getSourceId() ;
1041    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1042    vector<AngularRate> srcpm = srow->getProperMotion() ;
1043    pm[0] = srcpm[0].get() ;
1044    pm[1] = srcpm[1].get() ;
1045  }
1046  return pm ;
1047}
1048
1049double ASDMReader::getSysVel( unsigned int idx )
1050{
1051  double sysvel = 0.0 ;
1052  unsigned int index = dataIdList_[idx] ;
1053  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1054  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1055  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1056  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1057  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1058  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1059  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1060  string srcname ;
1061  if ( frow->isSourceIdExists() ) {
1062    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1063    int sid = frow->getSourceId() ;
1064    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1065    if ( srow->isSysVelExists() ) {
1066      vector<Speed> sysvelV = srow->getSysVel() ;
1067      if ( sysvelV.size() > 0 )
1068        sysvel = sysvelV[0].get() ;
1069    }
1070  }
1071  return sysvel ;
1072}
1073
1074unsigned int ASDMReader::getFlagRow( unsigned int idx )
1075{
1076  return vmsData_->v_flag[dataIdList_[idx]] ;
1077}
1078
1079vector<unsigned int> ASDMReader::getDataShape( unsigned int idx )
1080{
1081  return vmsData_->vv_dataShape[dataIdList_[idx]] ;
1082}
1083
1084float * ASDMReader::getSpectrum( unsigned int idx )
1085{
1086  map<AtmPhaseCorrection, float*> data = vmsData_->v_m_data[dataIdList_[idx]] ;
1087  //map<AtmPhaseCorrection, float*>::iterator iter = data.find(AP_UNCORRECTED) ;
1088  map<AtmPhaseCorrection, float*>::iterator iter = data.find(apc_) ;
1089  float *autoCorr = iter->second ;
1090  return autoCorr ;
1091}
1092
1093// bool * ASDMReader::getFlagChannel( unsigned int idx )
1094// {
1095//   return 0 ;
1096// }
1097
1098vector< vector<float> > ASDMReader::getTsys( unsigned int idx )
1099{
1100  unsigned int index = dataIdList_[idx] ;
1101  Tag anttag( antennaId_, TagType::Antenna ) ;
1102  int feedid = vmsData_->v_feedId1[index] ;
1103  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1104  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1105  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1106  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1107  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1108  //int nchan = asdm_->getSpectralWindow().getRowByKey(spwtag)->getNumChan() ;
1109  vector< vector<float> > defaultTsys( 1, vector<float>( 1, 1.0 ) ) ;
1110  SysCalTable &sctab = asdm_->getSysCal() ;
1111  if ( sctab.size() == 0 )
1112    return defaultTsys ;
1113  SysCalRow *scrow = sctab.getRowByKey( anttag, spwtag, tint, feedid ) ;
1114  if ( scrow->isTsysSpectrumExists() ) {
1115    vector< vector<Temperature> > tsysSpec = scrow->getTsysSpectrum() ;
1116    unsigned int numReceptor = tsysSpec.size() ;
1117    unsigned int numChan = tsysSpec[0].size() ;
1118    vector< vector<float> > tsys( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1119    for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1120      for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1121        tsys[ir][ic] = tsysSpec[ir][ic].get() ;
1122      }
1123    }
1124    return tsys ;
1125  }
1126//   else if ( scrow->isTsysExists() ) {
1127//     vector<Temperature> tsysScalar = scrow->getTsys() ;
1128//     unsigned int numReceptor = tsysScalar.size() ;
1129//     vector< vector<float> > tsys( numReceptor ) ;
1130//     for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1131//       tsys[ir] = vector<float>( 1, tsysScalar[ir].get() ) ;
1132//     return tsys ;
1133//   }
1134  else {
1135    return defaultTsys ;
1136  }
1137}
1138
1139vector< vector<float> > ASDMReader::getTcal( unsigned int idx )
1140{
1141  unsigned int index = dataIdList_[idx] ;
1142  Tag anttag( antennaId_, TagType::Antenna ) ;
1143  int feedid = vmsData_->v_feedId1[index] ;
1144  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1145  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1146  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1147  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1148  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1149  //int nchan = asdm_->getSpectralWindow().getRowByKey(spwtag)->getNumChan() ;
1150  vector< vector<float> > defaultTcal( 1, vector<float>( 1, 1.0 ) ) ;
1151  SysCalTable &sctab = asdm_->getSysCal() ;
1152  if ( sctab.size() == 0 )
1153    return defaultTcal ;
1154  SysCalRow *scrow = sctab.getRowByKey( anttag, spwtag, tint, feedid ) ;
1155  if ( scrow->isTcalSpectrumExists() ) {
1156    vector< vector<Temperature> > tcalSpec = scrow->getTcalSpectrum() ;
1157    unsigned int numReceptor = tcalSpec.size() ;
1158    unsigned int numChan = tcalSpec[0].size() ;
1159    vector< vector<float> > tcal( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1160    for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1161      for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1162        tcal[ir][ic] = tcalSpec[ir][ic].get() ;
1163      }
1164    }
1165    return tcal ;
1166  }
1167//   else if ( scrow->isTcalExists() ) {
1168//     vector<Temperature> tcalScalar = scrow->getTcal() ;
1169//     unsigned int numReceptor = tcalScalar.size() ;
1170//     vector< vector<float> > tcal( numReceptor, vector<float>( 1, 1.0 ) ) ;
1171//     for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1172//       tcal[ir][0] = tcalScalar[ir][0].get() ;
1173//     return tcal ;
1174//   }
1175  else {
1176    return defaultTcal ;
1177  }
1178}
1179
1180vector<float> ASDMReader::getOpacity( unsigned int idx )
1181{
1182  vector<float> tau(0) ;
1183  CalAtmosphereTable &atmtab = asdm_->getCalAtmosphere() ;
1184  unsigned int nrow = atmtab.size() ;
1185  if ( nrow > 0 ) {
1186    unsigned int index = dataIdList_[idx] ;
1187    //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1188    double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1189    ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1190    //int feedid = vmsData_->v_feedId1[index] ;
1191    //Tag atag( antennaId_, TagType::Antenna ) ;
1192    //Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1193    //DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey(ddtag) ;
1194    //Tag spwtag = ddrow->getSpectralWindowId() ;
1195    //SpectralWindowRow *spwrow = ddrow->getSpectralWindowUsingSpectralWindowId() ;
1196    //BasebandName bbname = spwrow->getBasebandName() ;
1197    //FeedRow *frow = asdm_->getFeed().getRowByKey( atag, spwtag, tint, feedid ) ;
1198    //int nfeed = frow->getNumReceptor() ;
1199    //vector<ReceiverRow *> rrows = frow->getReceivers() ;
1200    vector<CalAtmosphereRow *> atmrows = atmtab.get() ;
1201    //ReceiverBand rb = rrows[0]->getFrequencyBand() ;
1202    int row0 = -1 ;
1203    double eps = tint.getStart().getMJD() ;
1204    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1205      CalAtmosphereRow *atmrow = atmrows[irow] ;
1206      if ( casa::String(atmrow->getAntennaName()) != antennaName_
1207           //|| atmrow->getReceiverBand() != rb
1208           //|| atmrow->getBasebandName() != bbname
1209           || atmrow->getCalDataUsingCalDataId()->getCalType() != CAL_ATMOSPHERE )
1210        continue ;
1211      else {
1212        double dt = tint.getStart().getMJD() - atmrow->getEndValidTime().getMJD() ;
1213        if ( dt >= 0 && dt < eps ) {
1214          eps = dt ;
1215          row0 = (int)irow ;
1216        }
1217      }
1218    }
1219    if ( row0 != -1 ) {
1220      CalAtmosphereRow *atmrow = atmrows[row0] ;
1221      tau = atmrow->getTau() ;
1222    }
1223  }
1224  return tau ;
1225}
1226
1227void ASDMReader::getWeatherInfo( unsigned int idx,
1228                                 float &temperature,
1229                                 float &pressure,
1230                                 float &humidity,
1231                                 float &windspeed,
1232                                 float &windaz )
1233{
1234  temperature = 0.0 ;
1235  pressure = 0.0 ;
1236  humidity = 0.0 ;
1237  windspeed = 0.0 ;
1238  windaz = 0.0 ;
1239
1240  //logsink_->postLocally( LogMessage("weatherStationId_ = "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1241
1242  WeatherTable &wtab = asdm_->getWeather() ;
1243  if ( wtab.size() == 0 || weatherStationId_ == -1 )
1244    return ;
1245
1246  unsigned int index = dataIdList_[idx] ;
1247  //Tag anttag( antennaId_, TagType::Antenna ) ;
1248  //Tag sttag = (asdm_->getAntenna().getRowByKey( anttag ))->getStationId() ;
1249  Tag sttag( (unsigned int)weatherStationId_, TagType::Station ) ;
1250  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1251  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1252  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1253  //WeatherRow *wrow = wtab.getRowByKey( sttag, tint ) ;
1254  vector<WeatherRow *> *wrows = wtab.getByContext( sttag ) ;
1255  WeatherRow *wrow = (*wrows)[0] ;
1256  unsigned int nrow = wrows->size() ;
1257  //logsink_->postLocally( LogMessage("There are "+String::toString(nrow)+" rows for given context: stationId "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1258  ArrayTime startTime = tint.getStart() ;
1259  if ( startTime < (*wrows)[0]->getTimeInterval().getStart() ) {
1260    temperature = (*wrows)[0]->getTemperature().get() ;
1261    pressure = (*wrows)[0]->getPressure().get() ;
1262    humidity = (*wrows)[0]->getRelHumidity().get() ;
1263    windspeed = (*wrows)[0]->getWindSpeed().get() ;
1264    windaz = (*wrows)[0]->getWindDirection().get() ;
1265  }
1266  else if ( startTime > (*wrows)[nrow-1]->getTimeInterval().getStart() ) {
1267    temperature = (*wrows)[nrow-1]->getTemperature().get() ;
1268    pressure = (*wrows)[nrow-1]->getPressure().get() ;
1269    humidity = (*wrows)[nrow-1]->getRelHumidity().get() ;
1270    windspeed = (*wrows)[nrow-1]->getWindSpeed().get() ;
1271    windaz = (*wrows)[nrow-1]->getWindDirection().get() ;
1272  }
1273  else {
1274    for ( unsigned int irow = 1 ; irow < wrows->size() ; irow++ ) {
1275      wrow = (*wrows)[irow-1] ;
1276      if ( startTime < (*wrows)[irow]->getTimeInterval().getStart() ) {
1277        temperature = wrow->getTemperature().get() ;
1278        pressure = wrow->getPressure().get() ;
1279        humidity = wrow->getRelHumidity().get() ;
1280        windspeed = wrow->getWindSpeed().get() ;
1281        windaz = wrow->getWindDirection().get() ;
1282        break ;
1283      }
1284    }
1285  }
1286
1287  return ;
1288}
1289
1290void ASDMReader::processStation()
1291{
1292  antennaPad_.resize(0) ;
1293  weatherStation_.resize(0) ;
1294  StationTable &stab = asdm_->getStation() ;
1295  vector<StationRow *> srows = stab.get() ;
1296  for ( unsigned int irow = 0 ; irow < srows.size() ; irow++ ) {
1297    StationType stype = srows[irow]->getType() ;
1298    Tag stag = srows[irow]->getStationId() ;
1299    if ( stype == ANTENNA_PAD )
1300      antennaPad_.push_back( stag ) ;
1301    else if ( stype == WEATHER_STATION )
1302      weatherStation_.push_back( stag ) ;
1303  }
1304
1305   weatherStationId_ = getClosestWeatherStation() ;
1306}
1307
1308int ASDMReader::getClosestWeatherStation()
1309{
1310  if ( weatherStation_.size() == 0 )
1311    return -1 ;
1312
1313  Tag atag( antennaId_, TagType::Antenna ) ;
1314  Tag stag = (asdm_->getAntenna().getRowByKey( atag ))->getStationId() ;
1315  vector<double> apos( 3 ) ;
1316  StationTable &stab = asdm_->getStation() ;
1317  StationRow *srow = stab.getRowByKey( stag ) ;
1318  vector<Length> pos = srow->getPosition() ;
1319  apos[0] = pos[0].get() ;
1320  apos[1] = pos[1].get() ;
1321  apos[2] = pos[2].get() ;
1322 
1323  double eps = 1.0e20 ;
1324  int retval = -1 ;
1325  for ( unsigned int ir = 0 ; ir < weatherStation_.size() ; ir++ ) {
1326    srow = stab.getRowByKey( weatherStation_[ir] ) ;
1327    vector<Length> wpos = srow->getPosition() ;
1328    double dist = (apos[0]-wpos[0].get())*(apos[0]-wpos[0].get())
1329      + (apos[1]-wpos[1].get())*(apos[1]-wpos[1].get())
1330      + (apos[2]-wpos[2].get())*(apos[2]-wpos[2].get()) ;
1331    if ( dist < eps ) {
1332      retval = (int)(weatherStation_[ir].getTagValue()) ;
1333    }
1334  }
1335
1336  return retval ;
1337}
1338
1339void ASDMReader::getPointingInfo( unsigned int idx,
1340                                  vector<double> &dir,
1341                                  double &az,
1342                                  double &el,
1343                                  vector<double> &srate )
1344{
1345  String funcName = "getPointingInfo" ;
1346
1347  dir.resize(0) ;
1348  az = -1.0 ;
1349  el = -1.0 ;
1350  srate.resize(0) ;
1351
1352  Tag atag( antennaId_, TagType::Antenna ) ;
1353  unsigned int index = dataIdList_[idx] ;
1354  vector<PointingRow *> *prows = asdm_->getPointing().getByContext( atag ) ;
1355  PointingRow *prow ;
1356  PointingRow *qrow ;
1357  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1358  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1359  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1360
1361  unsigned int nrow = prows->size() ;
1362  //logsink_->postLocally( LogMessage("There are " << nrow << " rows for given context: antennaId "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
1363
1364//   for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1365//     prow = (*prows)[irow] ;
1366//     ArrayTimeInterval ati = prow->getTimeInterval() ;
1367//     ArrayTime pst = ati.getStart() ;
1368//     ArrayTime pet( ati.getStartInMJD()+ati.getDurationInDays() ) ;
1369//     logsink_->postLocally( LogMessage("start: "+pst.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1370//     logsink_->postLocally( LogMessage("end: "+pet.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1371//   }
1372 
1373  if ( nrow == 0 )
1374    return ;
1375
1376  srate.resize(2) ;
1377  srate[0] = 0.0 ;
1378  srate[1] = 0.0 ;
1379  az = 0.0 ;
1380  el = 0.0 ;
1381  double tcen = 0.0 ;
1382
1383  //
1384  // shape of pointingDirection is (numSample,2) if usePolynomial = False, while it is
1385  // (numTerm,2) if usePolynomial = True.
1386  //
1387  // In the former case, typical sampled time interval is 48msec, which is very small
1388  // compared with typical integration time (~a few sec).
1389  // Scan rate for this case is always [0,0] (or get slope?).
1390  //
1391  // In the later case, scan rate is (pointingDirection[1][0],pointingDirection[1][1])
1392  //
1393  // PointingDirection seems to store AZEL
1394  //
1395  ArrayTimeInterval pTime0 = (*prows)[0]->getTimeInterval() ;
1396  ArrayTimeInterval pTime1 = (*prows)[nrow-1]->getTimeInterval() ;
1397  if ( tint.getStartInMJD()+tint.getDurationInDays() < pTime0.getStartInMJD() ) {
1398    logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint < ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1399    prow = (*prows)[0] ;
1400    //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1401    vector< vector<Angle> > dirA = prow->getTarget() ;
1402    vector< vector<Angle> > offA = prow->getOffset() ;
1403    //az = dirA[0][0].get() ;
1404    //el = dirA[0][1].get() ;
1405    az = dirA[0][0].get() + offA[0][0].get() ;
1406    el = dirA[0][1].get() + offA[0][1].get() ;
1407    if ( prow->getUsePolynomials() && prow->getNumTerm() > 1 ) {
1408      //srate[0] = dirA[1][0].get() ;
1409      //srate[1] = dirA[1][1].get() ;
1410      srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1411      srate[1] = dirA[1][1].get() + offA[1][1].get() ;
1412    }     
1413  }
1414  else if ( tint.getStartInMJD() > pTime1.getStartInMJD()+pTime1.getDurationInDays() ) {
1415    logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint > ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1416    prow = (*prows)[nrow-1] ;
1417    int numSample = prow->getNumSample() ;
1418    //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1419    vector< vector<Angle> > dirA = prow->getTarget() ;
1420    vector< vector<Angle> > offA = prow->getOffset() ;
1421    if ( prow->getUsePolynomials() ) {
1422      //az = dirA[0][0].get() ;
1423      //el = dirA[0][1].get() ;
1424      az = dirA[0][0].get() + offA[0][0].get() ;
1425      el = dirA[0][1].get() + offA[0][1].get() ;
1426      if ( prow->getNumTerm() > 1 ) {
1427        //srate[0] = dirA[1][0].get() ;
1428        //srate[1] = dirA[1][1].get() ;
1429        srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1430        srate[1] = dirA[1][1].get() + offA[1][0].get() ;
1431      }
1432    }
1433    else {
1434      //az = dirA[numSample-1][0].get() ;
1435      //el = dirA[numSample-1][1].get() ;
1436      az = dirA[numSample-1][0].get() + offA[numSample-1][0].get() ;
1437      el = dirA[numSample-1][1].get() + offA[numSample-1][1].get() ;
1438    }
1439  }
1440  else {
1441    ArrayTime startTime = tint.getStart() ;
1442    ArrayTime endTime( tint.getStartInMJD()+tint.getDurationInDays() ) ;
1443    int row0 = -1 ;
1444    int row1 = -1 ;
1445    int row2 = -1 ;
1446    double dt0 = getMidTime( tint ).getMJD() ;
1447    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1448      prow = (*prows)[irow] ;
1449      double dt = getMidTime( tint ).getMJD() - getMidTime( prow->getTimeInterval() ).getMJD() ;
1450      if ( dt > 0 && dt < dt0 ) {
1451        dt0 = dt ;
1452        row2 = irow ;
1453      }
1454      if ( prow->getTimeInterval().contains( startTime ) )
1455        row0 = irow ;
1456      else if ( prow->getTimeInterval().contains( endTime ) )
1457        row1 = irow ;
1458      if ( row0 != -1 && row1 != -1 )
1459        break ;
1460    }
1461    //logsink_->postLocally( LogMessage("row0 = "+String::toString(row0)+", row1 = "+String::toString(row1)+", row2 = "+String::toString(row2),LogOrigin(className_,funcName,WHERE)) ) ;
1462    if ( row0 == -1 && row1 == -1 ) {
1463      prow = (*prows)[row2] ;
1464      qrow = (*prows)[row2+1] ;
1465      double t0 = getEndTime( prow->getTimeInterval() ).getMJD() ;
1466      double t1 = qrow->getTimeInterval().getStartInMJD() ;
1467      double t = startTime.getMJD() ;
1468      //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1469      //vector< vector<Angle> > dirB = qrow->getPointingDirection() ;
1470      vector< vector<Angle> > dirA = prow->getTarget() ;
1471      vector< vector<Angle> > offA = prow->getOffset() ;
1472      vector< vector<Angle> > dirB = qrow->getTarget() ;
1473      vector< vector<Angle> > offB = qrow->getOffset() ;
1474      //double da0 = dirA[0][0].get() ;
1475      //double db0 = dirB[0][0].get() ;
1476      //double da1 = dirA[0][1].get() ;
1477      //double db1 = dirB[0][1].get() ;
1478      double da0 = dirA[0][0].get() + offA[0][0].get() ;
1479      double db0 = dirB[0][0].get() + offB[0][0].get() ;
1480      double da1 = dirA[0][1].get() + offA[0][1].get() ;
1481      double db1 = dirB[0][1].get() + offB[0][1].get() ;
1482      if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1483        double dt = ( t - t0 ) / ( t1 - t0 ) ;
1484        az = da0 + ( db0 - da0 ) * dt ;
1485        el = da1 + ( db1 - da1 ) * dt ;
1486        if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1487          //double ra0 = dirA[1][0].get() ;
1488          //double rb0 = dirB[1][0].get() ;
1489          //double ra1 = dirA[1][1].get() ;
1490          //double rb1 = dirB[1][1].get() ;
1491          double ra0 = dirA[1][0].get() + offA[1][0].get() ;
1492          double rb0 = dirB[1][0].get() + offB[1][0].get() ;
1493          double ra1 = dirA[1][1].get() + offA[1][1].get() ;
1494          double rb1 = dirB[1][1].get() + offB[1][1].get() ;         
1495          srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1496          srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1497        }
1498      }
1499      else if ( !(qrow->getUsePolynomials()) ) {
1500        double dt = ( t - t0 ) / ( t1 - t0 ) ;
1501        az = da0 + ( db0 - da0 ) * dt ;
1502        el = da1 + ( db1 - da1 ) * dt ;
1503      }
1504      //else {
1505      // nothing to do
1506      //}
1507    }
1508    else {
1509      int ndir = 0 ;
1510      if ( row0 == -1 ) {
1511        row0 = row1 ;
1512      }
1513      else if ( row1 == -1 ) {
1514        row1 = row0 ;
1515      }
1516      prow = (*prows)[row0] ;
1517      qrow = (*prows)[row1] ;
1518      if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1519        //logsink_->postLocally( LogMessage("usePolynomial = True",LogOrigin(className_,funcName,WHERE)) ) ;
1520        if ( row0 == row1 ) {
1521          prow = (*prows)[row0] ;
1522          //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1523          vector< vector<Angle> > dirA = prow->getTarget() ;
1524          vector< vector<Angle> > offA = prow->getOffset() ;
1525          //az = dirA[0][0].get() ;
1526          //el = dirA[0][1].get() ;
1527          az = dirA[0][0].get() + offA[0][0].get() ;
1528          el = dirA[0][1].get() + offA[0][1].get() ;
1529          if ( prow->getNumTerm() > 1 ) {
1530            //srate[0] = dirA[1][0].get() ;
1531            //srate[1] = dirA[1][1].get() ;
1532            srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1533            srate[1] = dirA[1][1].get() + offA[1][0].get() ;
1534          }
1535        }
1536        else {
1537          prow = (*prows)[row0] ;
1538          qrow = (*prows)[row1] ;
1539          //vector< vector<Angle> > dirA = qrow->getPointingDirection() ;
1540          //vector< vector<Angle> > dirB = qrow->getPointingDirection() ;
1541          vector< vector<Angle> > dirA = qrow->getTarget() ;
1542          vector< vector<Angle> > offA = prow->getOffset() ;
1543          vector< vector<Angle> > dirB = qrow->getTarget() ;
1544          vector< vector<Angle> > offB = qrow->getOffset() ;
1545          double t0 = getMidTime( prow->getTimeInterval() ).getMJD() ;
1546          double t1 = getMidTime( qrow->getTimeInterval() ).getMJD() ;
1547          double t = startTime.getMJD() ;
1548          double dt = ( t - t0 ) / ( t1 - t0 ) ;
1549          //double da0 = dirA[0][0].get() ;
1550          //double db0 = dirB[0][0].get() ;
1551          //double da1 = dirA[0][1].get() ;
1552          //double db1 = dirB[0][1].get() ;
1553          double da0 = dirA[0][0].get() + offA[0][0].get() ;
1554          double db0 = dirB[0][0].get() + offB[0][0].get() ;
1555          double da1 = dirA[0][1].get() + offA[0][1].get() ;
1556          double db1 = dirB[0][1].get() + offB[0][1].get() ;
1557          az = da0 + ( db0 - da0 ) * dt ;
1558          el = da1 + ( db1 - db0 ) * dt ;
1559          if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1560            //double ra0 = dirA[1][0].get() ;
1561            //double rb0 = dirB[1][0].get() ;
1562            //double ra1 = dirA[1][1].get() ;
1563            //double rb1 = dirB[1][1].get() ;
1564            double ra0 = dirA[1][0].get() + offA[1][0].get() ;
1565            double rb0 = dirB[1][0].get() + offB[1][0].get() ;
1566            double ra1 = dirA[1][1].get() + offA[1][1].get() ;
1567            double rb1 = dirB[1][1].get() + offB[1][1].get() ;         
1568            srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1569            srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1570          }
1571        }
1572      }
1573      else if ( prow->getUsePolynomials() == qrow->getUsePolynomials() ) {
1574        //logsink_->postLocally( LogMessage("numSample == numTerm",LogOrigin(className_,funcName,WHERE)) ) ;
1575        for ( int irow = row0 ; irow <= row1 ; irow++ ) {
1576          prow = (*prows)[irow] ;
1577          int numSample = prow->getNumSample() ;
1578          //logsink_->postLocally( LogMessage("numSample = "+String::toString(numSample),LogOrigin(className_,funcName,WHERE)) ) ;
1579          //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1580          vector< vector<Angle> > dirA = prow->getTarget() ;
1581          vector< vector<Angle> > offA = prow->getOffset() ;
1582          if ( prow->isSampledTimeIntervalExists() ) {
1583            //logsink_->postLocally( LogMessage("sampledTimeIntervalExists",LogOrigin(className_,funcName,WHERE)) ) ;
1584            vector<ArrayTimeInterval> stime = prow->getSampledTimeInterval() ;
1585            for ( int isam = 0 ; isam < numSample ; isam++ ) {
1586              //if ( tint.overlaps( stime[isam] ) ) {
1587              if ( tint.contains( stime[isam] ) ) {
1588                //az += dirA[isam][0].get() ;
1589                //el += dirA[isam][1].get() ;
1590                az += dirA[isam][0].get() + offA[isam][0].get() ;
1591                el += dirA[isam][1].get() + offA[isam][1].get() ;
1592                ndir++ ;
1593              }
1594            }
1595          }
1596          else {
1597            double sampleStart = prow->getTimeInterval().getStartInMJD() ;
1598            double sampleInterval = prow->getTimeInterval().getDurationInDays() / (double)numSample ;
1599            //logsink_->postLocally( LogMessage("sampleStart = "+String::toString(sampleStart),LogOrigin(className_,funcName,WHERE)) )
1600            //logsink_->postLocally( LogMessage("sampleInterval = "+String::toString(sampleInterval),LogOrigin(className_,funcName,WHERE)) ) ;
1601            //logsink_->postLocally( LogMessage("tint = "+tint.toString(),LogOrigin(className_,funcName,WHERE)) ) ;
1602            for ( int isam = 0 ; isam < numSample ; isam++ ) {
1603              ArrayTimeInterval stime( sampleStart+isam*sampleInterval, sampleInterval ) ;
1604              //if ( tint.overlaps( stime ) ) {
1605              if ( tint.contains( stime ) ) {
1606                //az += dirA[isam][0].get() ;
1607                //el += dirA[isam][1].get() ;
1608                az += dirA[isam][0].get() + offA[isam][0].get() ;
1609                el += dirA[isam][1].get() + offA[isam][1].get() ;
1610                tcen += getMidTime( stime ).getMJD() ;
1611                ndir++ ;
1612              }
1613            }
1614          }
1615        }
1616        if ( ndir > 1 ) {
1617          az /= (double)ndir ;
1618          el /= (double)ndir ;
1619          tcen /= (double)ndir ;
1620        }
1621      }
1622      //else {
1623      // nothing to do
1624      //}
1625    }
1626   
1627    AntennaRow *arow = asdm_->getAntenna().getRowByKey( Tag( antennaId_, TagType::Antenna ) ) ;
1628    StationRow *srow = arow->getStationUsingStationId() ;
1629    vector<Length> antposL = srow->getPosition() ;
1630    casa::Vector<casa::Double> antpos( 3 ) ;
1631    for ( int i = 0 ; i < 3 ; i++ )
1632      antpos[i] = antposL[i].get() ;
1633    //logsink_->postLocally( LogMessage("tcen = "+String::toString(tcen),LogOrigin(className_,funcName,WHERE)) ) ;
1634    //logsink_->postLocally( LogMessage("antpos = "+String::toString(antpos),LogOrigin(className_,funcName,WHERE)) ) ;
1635    toJ2000( dir, az, el, tcen, antpos ) ;
1636
1637  }
1638
1639  return ;
1640}
1641
1642ArrayTime ASDMReader::getMidTime( const ArrayTimeInterval &t )
1643{
1644  return ArrayTime( t.getStartInMJD() + 0.5 * t.getDurationInDays() ) ;
1645}
1646
1647ArrayTime ASDMReader::getEndTime( const ArrayTimeInterval &t )
1648{
1649  return ArrayTime( t.getStartInMJD() + t.getDurationInDays() ) ;
1650}
1651
1652ArrayTime ASDMReader::getStartTime( const ArrayTimeInterval &t )
1653{
1654  return ArrayTime( t.getStartInMJD() ) ;
1655}
1656
1657void ASDMReader::toJ2000( vector<double> &dir,
1658                          double az,
1659                          double el,
1660                          double mjd,
1661                          casa::Vector<casa::Double> antpos )
1662{
1663  String funcName = "toJ2000" ;
1664
1665//   casa::Vector<casa::Double> azel( 2 ) ;
1666//   azel[0] = az ;
1667//   azel[1] = el ;
1668//   casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
1669//   casa::Vector<casa::Quantity> qantpos( 3 ) ;
1670//   qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
1671//   qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
1672//   qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
1673//   casa::MPosition mp( casa::MVPosition( qantpos ),
1674//                       casa::MPosition::ITRF ) ;
1675//   //ostringstream oss ;
1676//   //mp.print( oss ) ;
1677//   //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
1678
1679//   casa::MeasFrame mf( me, mp ) ;
1680//   casa::MDirection::Convert toj2000( casa::MDirection::AZELGEO,
1681//                                      casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
1682//   casa::Vector<casa::Double> cdir = toj2000( azel ).getAngle( "rad" ).getValue() ;
1683//   //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
1684//   dir.resize(2) ;
1685//   dir[0] = (double)(cdir[0]) ;
1686//   dir[1] = (double)(cdir[1]) ;
1687  vector<double> azel( 2 ) ;
1688  azel[0] = az ;
1689  azel[1] = el ;
1690  dir = toJ2000( azel, "AZELGEO", mjd, antpos ) ;
1691}
1692
1693vector<double> ASDMReader::toJ2000( vector<double> dir,
1694                                    casa::String dirref,
1695                                    double mjd,
1696                                    casa::Vector<casa::Double> antpos )
1697{
1698  casa::String funcname = "toJ2000" ;
1699
1700  vector<double> newd( dir ) ;
1701  if ( dirref != "J2000" ) {
1702    casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
1703    casa::Vector<casa::Quantity> qantpos( 3 ) ;
1704    qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
1705    qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
1706    qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
1707    casa::MPosition mp( casa::MVPosition( qantpos ),
1708                        casa::MPosition::ITRF ) ;
1709    //ostringstream oss ;
1710    //mp.print( oss ) ;
1711    //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
1712   
1713    casa::MeasFrame mf( me, mp ) ;
1714    casa::MDirection::Types dirtype ;
1715    casa::Bool b = casa::MDirection::getType( dirtype, dirref ) ;
1716    if ( b ) {
1717      casa::MDirection::Convert toj2000( dirtype,
1718                                         casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
1719      casa::Vector<casa::Double> cdir = toj2000( dir ).getAngle( "rad" ).getValue() ;
1720      //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
1721      newd[0] = (double)(cdir[0]) ;
1722      newd[1] = (double)(cdir[1]) ;
1723    }
1724  }
1725  return newd ;
1726}
1727
1728void ASDMReader::setLogger( CountedPtr<LogSinkInterface> &logsink )
1729{
1730  logsink_ = logsink ;
1731}
1732
1733string ASDMReader::getFrame()
1734{
1735  String funcName = "getFrame" ;
1736 
1737  // default is TOPO
1738  string frame = "TOPO" ;
1739
1740  SpectralWindowTable &spwtab = asdm_->getSpectralWindow() ;
1741  vector<SpectralWindowRow *> rows = spwtab.get() ;
1742  vector<FrequencyReferenceCode> measFreqRef( rows.size() ) ;
1743  int nref = 0 ;
1744  for ( unsigned int irow = 0 ; irow < rows.size() ; irow++ ) {
1745    int nchan = rows[irow]->getNumChan() ;
1746    if ( nchan != 4 ) {
1747      if ( rows[irow]->isMeasFreqRefExists() ) {
1748        measFreqRef[nref] = rows[irow]->getMeasFreqRef() ;
1749        nref++ ;
1750      }
1751    }
1752  }
1753  if ( nref != 0 ) {
1754    frame = CFrequencyReferenceCode::toString( measFreqRef[0] ) ;
1755  }
1756 
1757  //logsink_->postLocally( LogMessage("frame = "+String::toString(frame),LogOrigin(className_,funcName,WHERE)) ) ;
1758
1759  return frame ;
1760}
1761
1762int ASDMReader::getNumIFs()
1763{
1764  String funcName = "getNumIFs" ;
1765
1766  int nif = 0 ;
1767  vector<SpectralWindowRow *> rows = asdm_->getSpectralWindow().get() ;
1768  unsigned int nrow = rows.size() ;
1769  // check if all rows have freqGroup attribute
1770  bool freqGroupExists = true ;
1771  bool countedWvr = false ;
1772  for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1773    freqGroupExists &= rows[irow]->isFreqGroupExists() ;
1774    if ( rows[irow]->getNumChan() == 4 ) {
1775      if ( !countedWvr ) {
1776        countedWvr = true ;
1777        nif++ ;
1778      }
1779    }
1780    else {
1781      nif++ ;
1782    }
1783  }
1784 
1785  if ( freqGroupExists ) {
1786    vector<int> freqGroup(0) ;
1787    for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1788      int fg = rows[irow]->getFreqGroup() ;
1789      if ( (int)count( freqGroup.begin(), freqGroup.end(), fg ) == 0 ) {
1790        freqGroup.push_back( fg ) ;
1791      }
1792    }
1793    nif = freqGroup.size() ;
1794  }
1795
1796  //logsink_->postLocally( LogMessage("nif = "+String::toString(nif),LogOrigin(className_,funcName,WHERE)) ) ;
1797
1798  return nif ;
1799}
Note: See TracBrowser for help on using the repository browser.