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

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

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