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

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

New Development: No

JIRA Issue: Yes CAS-1913

Ready for Test: No

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...

Update of asdm2ASAP and related classes.

  • replaced LogIO with StreamLogSink?
  • all log messages are written to logger via StreamLogSink?
  • no output to stdout/stderr
  • commented out junk log
  • added time_sampling selection
  • supported wvr_corrected_data='both'
  • added logfile specification
  • added corr_mode selection


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