source: branches/casa-prerelease/pre-asap/external-alma/asdm2ASAP/ASDMReader.cc@ 2300

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

Do not access Subscan/SubscanMode if it doesn't exist.


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