source: branches/hpc34/external-alma/asdm2ASAP/ASDMReader.cc@ 3019

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

New Development: No

JIRA Issue: No

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

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