source: branches/casa-release-4_3/external-alma/asdm2ASAP/ASDMReader.cc@ 3074

Last change on this file since 3074 was 2755, checked in by Takeshi Nakazato, 12 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...

Speed up asdm2ASAP by reusing MDirection::Convert object.


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