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

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

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