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

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

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