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

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

New Development: No

JIRA Issue: Yes CAS-1913

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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