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

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

Support latest ASDM definition in asdm2ASAP.
Now asdm2ASAP is linked to libalma_v3.so to be able to read latest ASDM.

The importasdm task invokes asdm2ASAP when singledish is True.

Older version of asdm2ASAP is renamec as oldasdm2ASAP.
The oldasdm2ASAP and related classes are built by linking to libalma.so.
It is installed but not used by any tasks.


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