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

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

Support latest ASDM definition in asdm2ASAP.
Now asdm2ASAP is linked to libalma_v3.so to be able to read latest ASDM.

The importasdm task invokes asdm2ASAP when singledish is True.

Older version of asdm2ASAP is renamec as oldasdm2ASAP.
The oldasdm2ASAP and related classes are built by linking to libalma.so.
It is installed but not used by any tasks.


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 = subrow->getSubscanMode() ;
894 if ( scanIntent == OBSERVE_TARGET ) {
895 // on sky scan
896 if ( swmode == NO_SWITCHING || swmode == POSITION_SWITCHING ) {
897 // position switching
898 // tentatively set NO_SWITCHING = POSITION_SWITCHING
899 if ( subIntent == ON_SOURCE ) {
900 srctype = SrcType::PSON ;
901 }
902 else if ( subIntent == OFF_SOURCE ) {
903 srctype = SrcType::PSOFF ;
904 }
905 }
906 else if ( swmode == FREQUENCY_SWITCHING ) {
907 // frequency switching
908 if ( subIntent == ON_SOURCE ) {
909 srctype = SrcType::FSON ;
910 }
911 else if ( subIntent == OFF_SOURCE ) {
912 srctype = SrcType::FSOFF ;
913 }
914 }
915 else if ( swmode == NUTATOR_SWITCHING ) {
916 // nutator switching
917 if ( subIntent == ON_SOURCE ) {
918 srctype = SrcType::PSON ;
919 }
920 else if ( subIntent == OFF_SOURCE ) {
921 srctype = SrcType::PSOFF ;
922 }
923 }
924 else {
925 // other switching mode
926 if ( subIntent == ON_SOURCE ) {
927 srctype = SrcType::SIG ;
928 }
929 else if ( subIntent == OFF_SOURCE ) {
930 srctype = SrcType::REF ;
931 }
932 }
933 }
934 else if ( scanIntent == CALIBRATE_ATMOSPHERE ) {
935 if ( swmode == NO_SWITCHING || swmode == POSITION_SWITCHING ) {
936 // position switching
937 // tentatively set NO_SWITCHING = POSITION_SWITCHING
938 if ( subIntent == ON_SOURCE ) {
939 srctype = SrcType::PONCAL ;
940 }
941 else if ( subIntent == OFF_SOURCE ) {
942 srctype = SrcType::POFFCAL ;
943 }
944 }
945 else if ( swmode == FREQUENCY_SWITCHING ) {
946 // frequency switching
947 if ( subIntent == ON_SOURCE ) {
948 srctype = SrcType::FONCAL ;
949 }
950 else if ( subIntent == OFF_SOURCE ) {
951 srctype = SrcType::FOFFCAL ;
952 }
953 }
954 else if ( swmode == NUTATOR_SWITCHING ) {
955 // nutator switching
956 if ( subIntent == ON_SOURCE ) {
957 srctype = SrcType::PONCAL ;
958 }
959 else if ( subIntent == OFF_SOURCE ) {
960 srctype = SrcType::POFFCAL ;
961 }
962 }
963 else {
964 // other switching mode
965 if ( subIntent == ON_SOURCE ) {
966 srctype = SrcType::CAL ;
967 }
968 else if ( subIntent == OFF_SOURCE ) {
969 srctype = SrcType::CAL ;
970 }
971 }
972 }
973 else {
974 // off sky (e.g. calibrator device) scan
975 if ( subIntent == ON_SOURCE ) {
976 srctype = SrcType::SIG ;
977 }
978 else if ( subIntent == OFF_SOURCE ) {
979 srctype = SrcType::REF ;
980 }
981 else if ( subIntent == HOT ) {
982 srctype = SrcType::HOT ;
983 }
984 else if ( subIntent == AMBIENT ) {
985 srctype = SrcType::SKY ;
986 }
987 else {
988 srctype = SrcType::CAL ;
989 }
990 }
991
992 return srctype ;
993}
994
995unsigned int ASDMReader::getSubscanNo( unsigned int idx )
996{
997 //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)) ) ;
998 return vmsData_->v_msState[dataIdList_[idx]].subscanNum ;
999}
1000
1001vector<double> ASDMReader::getSourceDirection( unsigned int idx )
1002{
1003 vector<double> dir( 2, 0.0 ) ;
1004 unsigned int index = dataIdList_[idx] ;
1005 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1006 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1007 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1008 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1009 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1010 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1011 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1012 string srcname ;
1013 if ( frow->isSourceIdExists() ) {
1014 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1015 int sid = frow->getSourceId() ;
1016 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1017 vector<Angle> srcdir = srow->getDirection() ;
1018 dir[0] = limitedAngle( srcdir[0].get() ) ;
1019 dir[1] = limitedAngle( srcdir[1].get() ) ;
1020 if ( srow->isDirectionCodeExists() ) {
1021 DirectionReferenceCode dircode = srow->getDirectionCode() ;
1022 //logsink_->postLocally( LogMessage("dircode="+CDirectionReferenceCode::toString(dircode),LogOrigin(className_,funcName,WHERE)) ) ;
1023 if ( dircode != J2000 ) {
1024 // if not J2000, need direction conversion
1025 string ref = CDirectionReferenceCode::toString( dircode ) ;
1026 double mjd = vmsData_->v_time[index] * s2d ;
1027 Tag atag( antennaId_, TagType::Antenna ) ;
1028 AntennaRow *arow = asdm_->getAntenna().getRowByKey( atag ) ;
1029 StationRow *srow = arow->getStationUsingStationId() ;
1030 vector<Length> antposL = srow->getPosition() ;
1031 casa::Vector<casa::Double> antpos( 3 ) ;
1032 for ( int i = 0 ; i < 3 ; i++ )
1033 antpos[i] = antposL[i].get() ;
1034 dir = toJ2000( dir, ref, mjd, antpos ) ;
1035 }
1036 }
1037 }
1038 return dir ;
1039}
1040
1041void ASDMReader::getSourceDirection( unsigned int idx,
1042 vector<double> &dir,
1043 string &ref )
1044{
1045 dir.resize( 2 ) ;
1046 dir[0] = 0.0 ;
1047 dir[1] = 0.0 ;
1048 ref = "J2000" ;
1049 unsigned int index = dataIdList_[idx] ;
1050 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1051 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1052 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1053 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1054 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1055 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1056 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1057 string srcname ;
1058 if ( frow->isSourceIdExists() ) {
1059 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1060 int sid = frow->getSourceId() ;
1061 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1062 vector<Angle> srcdir = srow->getDirection() ;
1063 if ( srow->isDirectionCodeExists() ) {
1064 ref = CDirectionReferenceCode::toString( srow->getDirectionCode() ) ;
1065 }
1066 dir[0] = limitedAngle( srcdir[0].get() ) ;
1067 dir[1] = limitedAngle( srcdir[1].get() ) ;
1068 }
1069}
1070
1071void ASDMReader::getSourceDirection( vector<double> &dir, string &ref )
1072{
1073 dir.resize( 2 ) ;
1074 ref = "J2000" ; // default is J2000
1075 SourceTable &tab = asdm_->getSource() ;
1076 SourceRow *row = tab.get()[0] ;
1077 vector<Angle> dirA = row->getDirection() ;
1078 dir[0] = dirA[0].get() ;
1079 dir[1] = dirA[1].get() ;
1080 if ( row->isDirectionCodeExists() ) {
1081 ref = CDirectionReferenceCode::toString( row->getDirectionCode() ) ;
1082 }
1083}
1084
1085vector<double> ASDMReader::getSourceProperMotion( unsigned int idx )
1086{
1087 vector<double> pm( 2, 0.0 ) ;
1088 unsigned int index = dataIdList_[idx] ;
1089 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1090 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1091 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1092 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1093 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1094 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1095 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1096 string srcname ;
1097 if ( frow->isSourceIdExists() ) {
1098 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1099 int sid = frow->getSourceId() ;
1100 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1101 vector<AngularRate> srcpm = srow->getProperMotion() ;
1102 pm[0] = srcpm[0].get() ;
1103 pm[1] = srcpm[1].get() ;
1104 }
1105 return pm ;
1106}
1107
1108double ASDMReader::getSysVel( unsigned int idx )
1109{
1110 double sysvel = 0.0 ;
1111 unsigned int index = dataIdList_[idx] ;
1112 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1113 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1114 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1115 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1116 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1117 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1118 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1119 string srcname ;
1120 if ( frow->isSourceIdExists() ) {
1121 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1122 int sid = frow->getSourceId() ;
1123 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1124 if ( srow->isSysVelExists() ) {
1125 vector<Speed> sysvelV = srow->getSysVel() ;
1126 if ( sysvelV.size() > 0 )
1127 sysvel = sysvelV[0].get() ;
1128 }
1129 }
1130 return sysvel ;
1131}
1132
1133unsigned int ASDMReader::getFlagRow( unsigned int idx )
1134{
1135 return vmsData_->v_flag[dataIdList_[idx]] ;
1136}
1137
1138vector<unsigned int> ASDMReader::getDataShape( unsigned int idx )
1139{
1140 return vmsData_->vv_dataShape[dataIdList_[idx]] ;
1141}
1142
1143float * ASDMReader::getSpectrum( unsigned int idx )
1144{
1145 map<AtmPhaseCorrection, float*> data = vmsData_->v_m_data[dataIdList_[idx]] ;
1146 //map<AtmPhaseCorrection, float*>::iterator iter = data.find(AP_UNCORRECTED) ;
1147 map<AtmPhaseCorrection, float*>::iterator iter = data.find(apc_) ;
1148 float *autoCorr = iter->second ;
1149 return autoCorr ;
1150}
1151
1152// bool * ASDMReader::getFlagChannel( unsigned int idx )
1153// {
1154// return 0 ;
1155// }
1156
1157vector< vector<float> > ASDMReader::getTsys( unsigned int idx )
1158{
1159 vector< vector<float> > defaultTsys( 1, vector<float>( 1, 1.0 ) ) ;
1160 SysCalRow *scrow = getSysCalRow( idx ) ;
1161 if ( scrow != 0 && scrow->isTsysSpectrumExists() ) {
1162 vector< vector<Temperature> > tsysSpec = scrow->getTsysSpectrum() ;
1163 unsigned int numReceptor = tsysSpec.size() ;
1164 unsigned int numChan = tsysSpec[0].size() ;
1165 vector< vector<float> > tsys( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1166 for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1167 for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1168 tsys[ir][ic] = tsysSpec[ir][ic].get() ;
1169 }
1170 }
1171 return tsys ;
1172 }
1173// else if ( scrow->isTsysExists() ) {
1174// vector<Temperature> tsysScalar = scrow->getTsys() ;
1175// unsigned int numReceptor = tsysScalar.size() ;
1176// vector< vector<float> > tsys( numReceptor ) ;
1177// for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1178// tsys[ir] = vector<float>( 1, tsysScalar[ir].get() ) ;
1179// return tsys ;
1180// }
1181 else {
1182 return defaultTsys ;
1183 }
1184}
1185
1186vector< vector<float> > ASDMReader::getTcal( unsigned int idx )
1187{
1188 vector< vector<float> > defaultTcal( 1, vector<float>( 1, 1.0 ) ) ;
1189 SysCalRow *scrow = getSysCalRow( idx ) ;
1190 if ( scrow != 0 && scrow->isTcalSpectrumExists() ) {
1191 vector< vector<Temperature> > tcalSpec = scrow->getTcalSpectrum() ;
1192 unsigned int numReceptor = tcalSpec.size() ;
1193 unsigned int numChan = tcalSpec[0].size() ;
1194 vector< vector<float> > tcal( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1195 for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1196 for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1197 tcal[ir][ic] = tcalSpec[ir][ic].get() ;
1198 }
1199 }
1200 return tcal ;
1201 }
1202// else if ( scrow->isTcalExists() ) {
1203// vector<Temperature> tcalScalar = scrow->getTcal() ;
1204// unsigned int numReceptor = tcalScalar.size() ;
1205// vector< vector<float> > tcal( numReceptor, vector<float>( 1, 1.0 ) ) ;
1206// for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1207// tcal[ir][0] = tcalScalar[ir][0].get() ;
1208// return tcal ;
1209// }
1210 else {
1211 return defaultTcal ;
1212 }
1213}
1214
1215void ASDMReader::getTcalAndTsys( unsigned int idx,
1216 vector< vector<float> > &tcal,
1217 vector< vector<float> > &tsys )
1218{
1219 String funcName = "getTcalAndTsys" ;
1220
1221 vector< vector<float> > defaultT( 1, vector<float>( 1, 1.0 ) ) ;
1222 SysCalRow *scrow = getSysCalRow( idx ) ;
1223 if ( scrow == 0 ) {
1224 tcal = defaultT ;
1225 tsys = defaultT ;
1226 }
1227 else {
1228 if ( scrow->isTsysSpectrumExists() ) {
1229 vector< vector<Temperature> > tsysSpec = scrow->getTsysSpectrum() ;
1230 unsigned int numReceptor = tsysSpec.size() ;
1231 unsigned int numChan = tsysSpec[0].size() ;
1232 tsys.resize( numReceptor ) ;
1233 for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1234 tsys[ir].resize( numChan ) ;
1235 for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1236 tsys[ir][ic] = tsysSpec[ir][ic].get() ;
1237 }
1238 }
1239 }
1240 else {
1241 tsys = defaultT ;
1242 }
1243 if ( scrow->isTcalSpectrumExists() ) {
1244 vector< vector<Temperature> > tcalSpec = scrow->getTcalSpectrum() ;
1245 unsigned int numReceptor = tcalSpec.size() ;
1246 unsigned int numChan = tcalSpec[0].size() ;
1247 tcal.resize( numReceptor ) ;
1248 for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1249 tcal[ir].resize( numReceptor ) ;
1250 for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1251 tcal[ir][ic] = tcalSpec[ir][ic].get() ;
1252 }
1253 }
1254 }
1255 else {
1256 tcal = defaultT ;
1257 }
1258 }
1259}
1260
1261vector<float> ASDMReader::getOpacity( unsigned int idx )
1262{
1263 vector<float> tau(0) ;
1264 CalAtmosphereTable &atmtab = asdm_->getCalAtmosphere() ;
1265 unsigned int nrow = atmtab.size() ;
1266 if ( nrow > 0 ) {
1267 unsigned int index = dataIdList_[idx] ;
1268 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1269 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1270 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1271 //int feedid = vmsData_->v_feedId1[index] ;
1272 //Tag atag( antennaId_, TagType::Antenna ) ;
1273 //Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1274 //DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey(ddtag) ;
1275 //Tag spwtag = ddrow->getSpectralWindowId() ;
1276 //SpectralWindowRow *spwrow = ddrow->getSpectralWindowUsingSpectralWindowId() ;
1277 //BasebandName bbname = spwrow->getBasebandName() ;
1278 //FeedRow *frow = asdm_->getFeed().getRowByKey( atag, spwtag, tint, feedid ) ;
1279 //int nfeed = frow->getNumReceptor() ;
1280 //vector<ReceiverRow *> rrows = frow->getReceivers() ;
1281 vector<CalAtmosphereRow *> atmrows = atmtab.get() ;
1282 //ReceiverBand rb = rrows[0]->getFrequencyBand() ;
1283 int row0 = -1 ;
1284 double eps = tint.getStart().getMJD() ;
1285 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1286 CalAtmosphereRow *atmrow = atmrows[irow] ;
1287 if ( casa::String(atmrow->getAntennaName()) != antennaName_
1288 //|| atmrow->getReceiverBand() != rb
1289 //|| atmrow->getBasebandName() != bbname
1290 || atmrow->getCalDataUsingCalDataId()->getCalType() != CAL_ATMOSPHERE )
1291 continue ;
1292 else {
1293 double dt = tint.getStart().getMJD() - atmrow->getEndValidTime().getMJD() ;
1294 if ( dt >= 0 && dt < eps ) {
1295 eps = dt ;
1296 row0 = (int)irow ;
1297 }
1298 }
1299 }
1300 if ( row0 != -1 ) {
1301 CalAtmosphereRow *atmrow = atmrows[row0] ;
1302 tau = atmrow->getTau() ;
1303 }
1304 }
1305 else {
1306 tau.resize( 1 ) ;
1307 tau[0] = 0.0 ;
1308 }
1309 return tau ;
1310}
1311
1312void ASDMReader::getWeatherInfo( unsigned int idx,
1313 float &temperature,
1314 float &pressure,
1315 float &humidity,
1316 float &windspeed,
1317 float &windaz )
1318{
1319 casa::String funcName = "getWeatherInfo" ;
1320
1321 temperature = 0.0 ;
1322 pressure = 0.0 ;
1323 humidity = 0.0 ;
1324 windspeed = 0.0 ;
1325 windaz = 0.0 ;
1326
1327 //logsink_->postLocally( LogMessage("weatherStationId_ = "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1328
1329 WeatherTable &wtab = asdm_->getWeather() ;
1330 if ( wtab.size() == 0 || weatherStationId_ == -1 )
1331 return ;
1332
1333 unsigned int index = dataIdList_[idx] ;
1334 //Tag anttag( antennaId_, TagType::Antenna ) ;
1335 //Tag sttag = (asdm_->getAntenna().getRowByKey( anttag ))->getStationId() ;
1336 Tag sttag( (unsigned int)weatherStationId_, TagType::Station ) ;
1337 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1338 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1339 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1340 //WeatherRow *wrow = wtab.getRowByKey( sttag, tint ) ;
1341 vector<WeatherRow *> *wrows = wtab.getByContext( sttag ) ;
1342 WeatherRow *wrow = (*wrows)[0] ;
1343 unsigned int nrow = wrows->size() ;
1344 //logsink_->postLocally( LogMessage("There are "+String::toString(nrow)+" rows for given context: stationId "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1345 ArrayTime startTime = getMidTime( tint ) ;
1346 if ( startTime < (*wrows)[0]->getTimeInterval().getStart() ) {
1347 temperature = (*wrows)[0]->getTemperature().get() ;
1348 pressure = (*wrows)[0]->getPressure().get() ;
1349 humidity = (*wrows)[0]->getRelHumidity().get() ;
1350 windspeed = (*wrows)[0]->getWindSpeed().get() ;
1351 windaz = (*wrows)[0]->getWindDirection().get() ;
1352 }
1353 else if ( startTime > getEndTime( (*wrows)[nrow-1]->getTimeInterval() ) ) {
1354 temperature = (*wrows)[nrow-1]->getTemperature().get() ;
1355 pressure = (*wrows)[nrow-1]->getPressure().get() ;
1356 humidity = (*wrows)[nrow-1]->getRelHumidity().get() ;
1357 windspeed = (*wrows)[nrow-1]->getWindSpeed().get() ;
1358 windaz = (*wrows)[nrow-1]->getWindDirection().get() ;
1359 }
1360 else {
1361 for ( unsigned int irow = 1 ; irow < wrows->size() ; irow++ ) {
1362 wrow = (*wrows)[irow-1] ;
1363 ArrayTime tStart = wrow->getTimeInterval().getStart() ;
1364 ArrayTime tEnd = (*wrows)[irow]->getTimeInterval().getStart() ;
1365 if ( startTime >= tStart && startTime <= tEnd ) {
1366 temperature = wrow->getTemperature().get() ;
1367 pressure = wrow->getPressure().get() ;
1368 humidity = wrow->getRelHumidity().get() ;
1369 windspeed = wrow->getWindSpeed().get() ;
1370 windaz = wrow->getWindDirection().get() ;
1371 break ;
1372 }
1373 }
1374 }
1375
1376 // Pa -> hPa
1377 pressure /= 100.0 ;
1378
1379 return ;
1380}
1381
1382void ASDMReader::processStation()
1383{
1384 antennaPad_.resize(0) ;
1385 weatherStation_.resize(0) ;
1386 StationTable &stab = asdm_->getStation() ;
1387 vector<StationRow *> srows = stab.get() ;
1388 for ( unsigned int irow = 0 ; irow < srows.size() ; irow++ ) {
1389 StationType stype = srows[irow]->getType() ;
1390 Tag stag = srows[irow]->getStationId() ;
1391 if ( stype == ANTENNA_PAD )
1392 antennaPad_.push_back( stag ) ;
1393 else if ( stype == WEATHER_STATION )
1394 weatherStation_.push_back( stag ) ;
1395 }
1396
1397 weatherStationId_ = getClosestWeatherStation() ;
1398}
1399
1400int ASDMReader::getClosestWeatherStation()
1401{
1402 if ( weatherStation_.size() == 0 )
1403 return -1 ;
1404
1405 Tag atag( antennaId_, TagType::Antenna ) ;
1406 Tag stag = (asdm_->getAntenna().getRowByKey( atag ))->getStationId() ;
1407 vector<double> apos( 3 ) ;
1408 StationTable &stab = asdm_->getStation() ;
1409 StationRow *srow = stab.getRowByKey( stag ) ;
1410 vector<Length> pos = srow->getPosition() ;
1411 apos[0] = pos[0].get() ;
1412 apos[1] = pos[1].get() ;
1413 apos[2] = pos[2].get() ;
1414
1415 double eps = 1.0e20 ;
1416 int retval = -1 ;
1417 for ( unsigned int ir = 0 ; ir < weatherStation_.size() ; ir++ ) {
1418 srow = stab.getRowByKey( weatherStation_[ir] ) ;
1419 vector<Length> wpos = srow->getPosition() ;
1420 double dist = (apos[0]-wpos[0].get())*(apos[0]-wpos[0].get())
1421 + (apos[1]-wpos[1].get())*(apos[1]-wpos[1].get())
1422 + (apos[2]-wpos[2].get())*(apos[2]-wpos[2].get()) ;
1423 if ( dist < eps ) {
1424 retval = (int)(weatherStation_[ir].getTagValue()) ;
1425 }
1426 }
1427
1428 return retval ;
1429}
1430
1431void ASDMReader::getPointingInfo( unsigned int idx,
1432 vector<double> &dir,
1433 double &az,
1434 double &el,
1435 vector<double> &srate )
1436{
1437 String funcName = "getPointingInfo" ;
1438
1439 dir.resize(0) ;
1440 az = -1.0 ;
1441 el = -1.0 ;
1442 srate.resize(0) ;
1443
1444 Tag atag( antennaId_, TagType::Antenna ) ;
1445 unsigned int index = dataIdList_[idx] ;
1446 vector<PointingRow *> *prows = asdm_->getPointing().getByContext( atag ) ;
1447
1448 if ( prows == 0 )
1449 return ;
1450
1451 PointingRow *prow ;
1452 PointingRow *qrow ;
1453 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1454 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1455 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1456
1457 unsigned int nrow = prows->size() ;
1458 //logsink_->postLocally( LogMessage("There are " << nrow << " rows for given context: antennaId "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
1459
1460// for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1461// prow = (*prows)[irow] ;
1462// ArrayTimeInterval ati = prow->getTimeInterval() ;
1463// ArrayTime pst = ati.getStart() ;
1464// ArrayTime pet( ati.getStartInMJD()+ati.getDurationInDays() ) ;
1465// logsink_->postLocally( LogMessage("start: "+pst.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1466// logsink_->postLocally( LogMessage("end: "+pet.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1467// }
1468
1469 srate.resize(2) ;
1470 srate[0] = 0.0 ;
1471 srate[1] = 0.0 ;
1472 az = 0.0 ;
1473 el = 0.0 ;
1474 //double tcen = 0.0 ;
1475 double tcen = getMidTime( tint ).getMJD() ;
1476
1477 //
1478 // shape of pointingDirection is (numSample,2) if usePolynomial = False, while it is
1479 // (numTerm,2) if usePolynomial = True.
1480 //
1481 // In the former case, typical sampled time interval is 48msec, which is very small
1482 // compared with typical integration time (~a few sec).
1483 // Scan rate for this case is always [0,0] (or get slope?).
1484 //
1485 // In the later case, scan rate is (pointingDirection[1][0],pointingDirection[1][1])
1486 //
1487 // PointingDirection seems to store AZEL
1488 //
1489 ArrayTimeInterval pTime0 = (*prows)[0]->getTimeInterval() ;
1490 ArrayTimeInterval pTime1 = (*prows)[nrow-1]->getTimeInterval() ;
1491 //if ( tint.getStartInMJD()+tint.getDurationInDays() < pTime0.getStartInMJD() ) {
1492 if ( getEndTime( tint ) < getStartTime( pTime0 ) ) {
1493 logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint < ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1494 prow = (*prows)[0] ;
1495 vector< vector<double> > dirA = pointingDir( prow ) ;
1496 az = dirA[0][0] ;
1497 el = dirA[0][1] ;
1498 if ( prow->getUsePolynomials() && prow->getNumTerm() > 1 ) {
1499 srate[0] = dirA[1][0] ;
1500 srate[1] = dirA[1][1] ;
1501 }
1502 }
1503 //else if ( tint.getStartInMJD() > pTime1.getStartInMJD()+pTime1.getDurationInDays() ) {
1504 else if ( getStartTime( tint ) > getEndTime( pTime1 ) ) {
1505 logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint > ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1506 prow = (*prows)[nrow-1] ;
1507 int numSample = prow->getNumSample() ;
1508 vector< vector<double> > dirA = pointingDir( prow ) ;
1509 if ( prow->getUsePolynomials() ) {
1510 az = dirA[0][0] ;
1511 el = dirA[0][1] ;
1512 if ( prow->getNumTerm() > 1 ) {
1513 srate[0] = dirA[1][0] ;
1514 srate[1] = dirA[1][1] ;
1515 }
1516 }
1517 else {
1518 az = dirA[numSample-1][0] ;
1519 el = dirA[numSample-1][1] ;
1520 }
1521 }
1522 else {
1523 ArrayTime startTime = tint.getStart() ;
1524 ArrayTime endTime = getEndTime( tint ) ;
1525 int row0 = -1 ;
1526 int row1 = -1 ;
1527 int row2 = -1 ;
1528 double dt0 = getMidTime( tint ).getMJD() ;
1529 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1530 prow = (*prows)[irow] ;
1531 double dt = getMidTime( tint ).getMJD() - getMidTime( prow->getTimeInterval() ).getMJD() ;
1532 if ( dt > 0 && dt < dt0 ) {
1533 dt0 = dt ;
1534 row2 = irow ;
1535 }
1536 if ( prow->getTimeInterval().contains( startTime ) )
1537 row0 = irow ;
1538 else if ( prow->getTimeInterval().contains( endTime ) )
1539 row1 = irow ;
1540 if ( row0 != -1 && row1 != -1 )
1541 break ;
1542 }
1543 //logsink_->postLocally( LogMessage("row0 = "+String::toString(row0)+", row1 = "+String::toString(row1)+", row2 = "+String::toString(row2),LogOrigin(className_,funcName,WHERE)) ) ;
1544 if ( row0 == -1 && row1 == -1 ) {
1545 prow = (*prows)[row2] ;
1546 qrow = (*prows)[row2+1] ;
1547 double t0 = getEndTime( prow->getTimeInterval() ).getMJD() ;
1548 double t1 = qrow->getTimeInterval().getStartInMJD() ;
1549 vector< vector<double> > dirA = pointingDir( prow ) ;
1550 vector< vector<double> > dirB = pointingDir( qrow ) ;
1551 double da0 = dirA[0][0] ;
1552 double db0 = dirB[0][0] ;
1553 double da1 = dirA[0][1] ;
1554 double db1 = dirB[0][1] ;
1555 if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1556 double dt = ( tcen - t0 ) / ( t1 - t0 ) ;
1557 az = da0 + ( db0 - da0 ) * dt ;
1558 el = da1 + ( db1 - da1 ) * dt ;
1559 if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1560 double ra0 = dirA[1][0] ;
1561 double rb0 = dirB[1][0] ;
1562 double ra1 = dirA[1][1] ;
1563 double rb1 = dirB[1][1] ;
1564 srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1565 srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1566 }
1567 }
1568 else if ( !(qrow->getUsePolynomials()) ) {
1569 double dt = ( tcen - t0 ) / ( t1 - t0 ) ;
1570 az = da0 + ( db0 - da0 ) * dt ;
1571 el = da1 + ( db1 - da1 ) * dt ;
1572 }
1573 //else {
1574 // nothing to do
1575 //}
1576 }
1577 else {
1578 int ndir = 0 ;
1579 if ( row0 == -1 ) {
1580 row0 = row1 ;
1581 }
1582 else if ( row1 == -1 ) {
1583 row1 = row0 ;
1584 }
1585 prow = (*prows)[row0] ;
1586 qrow = (*prows)[row1] ;
1587 if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1588 //logsink_->postLocally( LogMessage("usePolynomial = True",LogOrigin(className_,funcName,WHERE)) ) ;
1589 if ( row0 == row1 ) {
1590 prow = (*prows)[row0] ;
1591 vector< vector<double> > dirA = pointingDir( prow ) ;
1592 az = dirA[0][0] ;
1593 el = dirA[0][1] ;
1594 if ( prow->getNumTerm() > 1 ) {
1595 srate[0] = dirA[1][0] ;
1596 srate[1] = dirA[1][1] ;
1597 }
1598 }
1599 else {
1600 prow = (*prows)[row0] ;
1601 qrow = (*prows)[row1] ;
1602 vector< vector<double> > dirA = pointingDir( prow ) ;
1603 vector< vector<double> > dirB = pointingDir( qrow ) ;
1604 double t0 = getMidTime( prow->getTimeInterval() ).getMJD() ;
1605 double t1 = getMidTime( qrow->getTimeInterval() ).getMJD() ;
1606 double dt = ( tcen - t0 ) / ( t1 - t0 ) ;
1607 double da0 = dirA[0][0] ;
1608 double db0 = dirB[0][0] ;
1609 double da1 = dirA[0][1] ;
1610 double db1 = dirB[0][1] ;
1611 az = da0 + ( db0 - da0 ) * dt ;
1612 el = da1 + ( db1 - db0 ) * dt ;
1613 if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1614 double ra0 = dirA[1][0] ;
1615 double rb0 = dirB[1][0] ;
1616 double ra1 = dirA[1][1] ;
1617 double rb1 = dirB[1][1] ;
1618 srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1619 srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1620 }
1621 }
1622 }
1623 else if ( prow->getUsePolynomials() == qrow->getUsePolynomials() ) {
1624 //logsink_->postLocally( LogMessage("numSample == numTerm",LogOrigin(className_,funcName,WHERE)) ) ;
1625 tcen = 0.0 ;
1626 for ( int irow = row0 ; irow <= row1 ; irow++ ) {
1627 prow = (*prows)[irow] ;
1628 int numSample = prow->getNumSample() ;
1629 //logsink_->postLocally( LogMessage("numSample = "+String::toString(numSample),LogOrigin(className_,funcName,WHERE)) ) ;
1630 vector< vector<double> > dirA = pointingDir( prow ) ;
1631 if ( prow->isSampledTimeIntervalExists() ) {
1632 //logsink_->postLocally( LogMessage("sampledTimeIntervalExists",LogOrigin(className_,funcName,WHERE)) ) ;
1633 vector<ArrayTimeInterval> stime = prow->getSampledTimeInterval() ;
1634 for ( int isam = 0 ; isam < numSample ; isam++ ) {
1635 //if ( tint.overlaps( stime[isam] ) ) {
1636 if ( tint.contains( stime[isam] ) ) {
1637 az += dirA[isam][0] ;
1638 el += dirA[isam][1] ;
1639 tcen += getMidTime( stime[isam] ).getMJD() ;
1640 ndir++ ;
1641 }
1642 }
1643 }
1644 else {
1645 double sampleStart = prow->getTimeInterval().getStartInMJD() ;
1646 double sampleInterval = prow->getTimeInterval().getDurationInDays() / (double)numSample ;
1647 //logsink_->postLocally( LogMessage("sampleStart = "+String::toString(sampleStart),LogOrigin(className_,funcName,WHERE)) )
1648 //logsink_->postLocally( LogMessage("sampleInterval = "+String::toString(sampleInterval),LogOrigin(className_,funcName,WHERE)) ) ;
1649 //logsink_->postLocally( LogMessage("tint = "+tint.toString(),LogOrigin(className_,funcName,WHERE)) ) ;
1650 for ( int isam = 0 ; isam < numSample ; isam++ ) {
1651 ArrayTimeInterval stime( sampleStart+isam*sampleInterval, sampleInterval ) ;
1652 //if ( tint.overlaps( stime ) ) {
1653 if ( tint.contains( stime ) ) {
1654 az += dirA[isam][0] ;
1655 el += dirA[isam][1] ;
1656 tcen += getMidTime( stime ).getMJD() ;
1657 ndir++ ;
1658 }
1659 }
1660 }
1661 }
1662 if ( ndir > 1 ) {
1663 az /= (double)ndir ;
1664 el /= (double)ndir ;
1665 tcen /= (double)ndir ;
1666 }
1667 }
1668 //else {
1669 // nothing to do
1670 //}
1671 }
1672
1673 AntennaRow *arow = asdm_->getAntenna().getRowByKey( Tag( antennaId_, TagType::Antenna ) ) ;
1674 StationRow *srow = arow->getStationUsingStationId() ;
1675 vector<Length> antposL = srow->getPosition() ;
1676 casa::Vector<casa::Double> antpos( 3 ) ;
1677 for ( int i = 0 ; i < 3 ; i++ )
1678 antpos[i] = antposL[i].get() ;
1679 //logsink_->postLocally( LogMessage("tcen = "+String::toString(tcen),LogOrigin(className_,funcName,WHERE)) ) ;
1680 //logsink_->postLocally( LogMessage("antpos = "+String::toString(antpos),LogOrigin(className_,funcName,WHERE)) ) ;
1681 toJ2000( dir, az, el, tcen, antpos ) ;
1682
1683 }
1684
1685 return ;
1686}
1687
1688ArrayTime ASDMReader::getMidTime( const ArrayTimeInterval &t )
1689{
1690 return ArrayTime( t.getStartInMJD() + 0.5 * t.getDurationInDays() ) ;
1691}
1692
1693ArrayTime ASDMReader::getEndTime( const ArrayTimeInterval &t )
1694{
1695 return ArrayTime( t.getStartInMJD() + t.getDurationInDays() ) ;
1696}
1697
1698ArrayTime ASDMReader::getStartTime( const ArrayTimeInterval &t )
1699{
1700 return ArrayTime( t.getStartInMJD() ) ;
1701}
1702
1703void ASDMReader::toJ2000( vector<double> &dir,
1704 double az,
1705 double el,
1706 double mjd,
1707 casa::Vector<casa::Double> antpos )
1708{
1709 String funcName = "toJ2000" ;
1710
1711 vector<double> azel( 2 ) ;
1712 azel[0] = az ;
1713 azel[1] = el ;
1714 dir = toJ2000( azel, "AZELGEO", mjd, antpos ) ;
1715}
1716
1717vector<double> ASDMReader::toJ2000( vector<double> dir,
1718 casa::String dirref,
1719 double mjd,
1720 casa::Vector<casa::Double> antpos )
1721{
1722 casa::String funcName = "toJ2000" ;
1723
1724 vector<double> newd( dir ) ;
1725 if ( dirref != "J2000" ) {
1726 casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
1727 casa::Vector<casa::Quantity> qantpos( 3 ) ;
1728 qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
1729 qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
1730 qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
1731 casa::MPosition mp( casa::MVPosition( qantpos ),
1732 casa::MPosition::ITRF ) ;
1733 //ostringstream oss ;
1734 //mp.print( oss ) ;
1735 //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
1736
1737 casa::MeasFrame mf( me, mp ) ;
1738 casa::MDirection::Types dirtype ;
1739 casa::Bool b = casa::MDirection::getType( dirtype, dirref ) ;
1740 if ( b ) {
1741 casa::MDirection::Convert toj2000( dirtype,
1742 casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
1743 casa::Vector<casa::Double> cdir = toj2000( dir ).getAngle( "rad" ).getValue() ;
1744 //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
1745 newd[0] = (double)(cdir[0]) ;
1746 newd[1] = (double)(cdir[1]) ;
1747 }
1748 }
1749 return newd ;
1750}
1751
1752void ASDMReader::setLogger( CountedPtr<LogSinkInterface> &logsink )
1753{
1754 logsink_ = logsink ;
1755}
1756
1757string ASDMReader::getFrame()
1758{
1759 casa::String funcName = "getFrame" ;
1760
1761 // default is TOPO
1762 string frame = "TOPO" ;
1763
1764 SpectralWindowTable &spwtab = asdm_->getSpectralWindow() ;
1765 vector<SpectralWindowRow *> rows = spwtab.get() ;
1766 vector<FrequencyReferenceCode> measFreqRef( rows.size() ) ;
1767 int nref = 0 ;
1768 for ( unsigned int irow = 0 ; irow < rows.size() ; irow++ ) {
1769 int nchan = rows[irow]->getNumChan() ;
1770 if ( nchan != 4 ) {
1771 if ( rows[irow]->isMeasFreqRefExists() ) {
1772 measFreqRef[nref] = rows[irow]->getMeasFreqRef() ;
1773 nref++ ;
1774 }
1775 }
1776 }
1777 if ( nref != 0 ) {
1778 frame = CFrequencyReferenceCode::toString( measFreqRef[0] ) ;
1779 }
1780
1781 //logsink_->postLocally( LogMessage("frame = "+String::toString(frame),LogOrigin(className_,funcName,WHERE)) ) ;
1782
1783 return frame ;
1784}
1785
1786int ASDMReader::getNumIFs()
1787{
1788 casa::String funcName = "getNumIFs" ;
1789
1790 int nif = 0 ;
1791 vector<SpectralWindowRow *> rows = asdm_->getSpectralWindow().get() ;
1792 unsigned int nrow = rows.size() ;
1793 // check if all rows have freqGroup attribute
1794 bool freqGroupExists = true ;
1795 bool countedWvr = false ;
1796 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1797 freqGroupExists &= rows[irow]->isFreqGroupExists() ;
1798 if ( rows[irow]->getNumChan() == 4 ) {
1799 if ( !countedWvr ) {
1800 countedWvr = true ;
1801 nif++ ;
1802 }
1803 }
1804 else {
1805 nif++ ;
1806 }
1807 }
1808
1809 if ( freqGroupExists ) {
1810 vector<int> freqGroup(0) ;
1811 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1812 int fg = rows[irow]->getFreqGroup() ;
1813 if ( (int)count( freqGroup.begin(), freqGroup.end(), fg ) == 0 ) {
1814 freqGroup.push_back( fg ) ;
1815 }
1816 }
1817 nif = freqGroup.size() ;
1818 }
1819
1820 //logsink_->postLocally( LogMessage("nif = "+String::toString(nif),LogOrigin(className_,funcName,WHERE)) ) ;
1821
1822 return nif ;
1823}
1824
1825SysCalRow *ASDMReader::getSysCalRow( unsigned int idx )
1826{
1827 String funcName = "getSysCalRow" ;
1828
1829 SysCalRow *row = 0 ;
1830 unsigned int index = dataIdList_[idx] ;
1831 Tag anttag( antennaId_, TagType::Antenna ) ;
1832 int feedid = vmsData_->v_feedId1[index] ;
1833 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1834 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1835 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1836 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1837 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1838 //int nchan = asdm_->getSpectralWindow().getRowByKey(spwtag)->getNumChan() ;
1839 vector< vector<float> > defaultTsys( 1, vector<float>( 1, 1.0 ) ) ;
1840 SysCalTable &sctab = asdm_->getSysCal() ;
1841 //vector<SysCalRow *> rows = sctab.get() ;
1842 vector<SysCalRow *> *rows = sctab.getByContext( anttag, spwtag, feedid ) ;
1843 //if ( nrow == 0 ) {
1844 if ( rows == 0 ) {
1845 //logsink_->postLocally( LogMessage("no rows in SysCal table",LogOrigin(className_,funcName,WHERE)) ) ;
1846 row = 0 ;
1847 }
1848 else {
1849 unsigned int nrow = rows->size() ;
1850 //logsink_->postLocally( LogMessage("nrow = "+String::toString(nrow),LogOrigin(className_,funcName,WHERE)) ) ;
1851 int scindex = -1 ;
1852 if ( nrow == 1 ) {
1853 scindex = 0 ;
1854 }
1855 else if ( getEndTime( tint ) <= getStartTime( (*rows)[0]->getTimeInterval() ) )
1856 scindex = 0 ;
1857 else {
1858 for ( unsigned int irow = 0 ; irow < nrow-1 ; irow++ ) {
1859 ArrayTime t = getMidTime( tint ) ;
1860 if ( t > getStartTime( (*rows)[irow]->getTimeInterval() )
1861 && t <= getStartTime( (*rows)[irow+1]->getTimeInterval() ) ) {
1862 scindex = irow ;
1863 break ;
1864 }
1865 }
1866 if ( scindex == -1 )
1867 scindex = nrow-1 ;
1868 }
1869 //logsink_->postLocally( LogMessage("irow = "+String::toString(scindex),LogOrigin(className_,funcName,WHERE)) ) ;
1870 row = (*rows)[scindex] ;
1871 }
1872 return row ;
1873}
1874
1875double ASDMReader::limitedAngle( double angle )
1876{
1877 if ( angle > C::pi )
1878 while ( angle > C::pi ) angle -= C::_2pi ;
1879 else if ( angle < -C::pi )
1880 while ( angle < -C::pi ) angle += C::_2pi ;
1881 return angle ;
1882}
1883
1884vector< vector<double> > ASDMReader::pointingDir( PointingRow *row )
1885{
1886 //String funcName = "pointingDir" ;
1887 vector< vector<Angle> > aTar = row->getTarget() ;
1888 vector< vector<Angle> > aOff = row->getOffset() ;
1889 vector< vector<Angle> > aDir = row->getPointingDirection() ;
1890 vector< vector<Angle> > aEnc = row->getEncoder() ;
1891 unsigned int n = aTar.size() ;
1892 vector< vector<double> > dir( n ) ;
1893 double factor = 1.0 / cos( aTar[0][1].get() ) ;
1894 for ( unsigned int i = 0 ; i < n ; i++ ) {
1895 dir[i].resize( 2 ) ;
1896 /**
1897 * This is approximate way to add offset taking tracking error
1898 * into account
1899 *
1900 * az = dir[0][0] = target[0][0] + offset[0][0] / cos(el)
1901 * + encorder[0][0] - direction[0][0]
1902 * el = dir[0][1] = target[0][1] + offset[0][1]
1903 * + encorder[0][1] - direction[0][1]
1904 **/
1905 dir[i][0] = aTar[i][0].get() + factor * aOff[i][0].get()
1906 + aEnc[i][0].get() - aDir[i][0].get() ;
1907 dir[i][1] = aTar[i][1].get() + aOff[i][1].get()
1908 + aEnc[i][1].get() - aDir[i][1].get() ;
1909 //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)) ) ;
1910 }
1911 return dir ;
1912}
Note: See TracBrowser for help on using the repository browser.