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

Last change on this file since 2209 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.