source: trunk/external-alma/oldasdm2ASAP/OldASDMReader.cc @ 2299

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

New Development: No

JIRA Issue: Yes CAS-1913

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

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