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

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

Frequency is stored as LSRK value although raw frequency value is
in arbitrary frequency frame.
Also, direction is stored as J2000 value independently of original
direction reference code.
File size: 67.0 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 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 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 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 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 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// 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 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{
709 String funcName = "getFrequency" ;
710
711 Tag ddTag( vmsData_->v_dataDescId[dataIdList_[idx]], TagType::DataDescription ) ;
712 DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey( ddTag ) ;
713 //Tag spwid = ddrow->getSpectralWindowId() ;
714 SpectralWindowRow *spwrow = ddrow->getSpectralWindowUsingSpectralWindowId() ;
715 int nchan = spwrow->getNumChan() ;
716 if ( nchan == 1 ) {
717 //logsink_->postLocally( LogMessage("channel averaged data",LogOrigin(className_,funcName,WHERE)) ) ;
718 refpix = 0.0 ;
719 incr = spwrow->getTotBandwidth().get() ;
720 if ( spwrow->isChanFreqStartExists() ) {
721 refval = spwrow->getChanFreqStart().get() ;
722 }
723 else if ( spwrow->isChanFreqArrayExists() ) {
724 refval = spwrow->getChanFreqArray()[0].get() ;
725 }
726 else {
727 throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
728 }
729 }
730 else if ( nchan % 2 ) {
731 // odd
732 //logsink_->postLocally( LogMessage("odd case",LogOrigin(className_,funcName,WHERE)) ) ;
733 refpix = 0.5 * ( (double)nchan - 1.0 ) ;
734 int ic = ( nchan - 1 ) / 2 ;
735 if ( spwrow->isChanWidthExists() ) {
736 incr = spwrow->getChanWidth().get() ;
737 }
738 else if ( spwrow->isChanWidthArrayExists() ) {
739 incr = spwrow->getChanWidthArray()[0].get() ;
740 }
741 else {
742 throw (casa::AipsError( "Either chanWidthArray or chanWidth must exist." )) ;
743 }
744 if ( spwrow->isChanFreqStepExists() ) {
745 if ( spwrow->getChanFreqStep().get() < 0.0 )
746 incr *= -1.0 ;
747 }
748 else if ( spwrow->isChanFreqArrayExists() ) {
749 vector<Frequency> chanFreqArr = spwrow->getChanFreqArray() ;
750 if ( chanFreqArr[0].get() > chanFreqArr[1].get() )
751 incr *= -1.0 ;
752 }
753 else {
754 throw (casa::AipsError( "Either chanFreqArray or chanFreqStep must exist." )) ;
755 }
756 if ( spwrow->isChanFreqStartExists() ) {
757 refval = spwrow->getChanFreqStart().get() + refpix * incr ;
758 }
759 else if ( spwrow->isChanFreqArrayExists() ) {
760 refval = spwrow->getChanFreqArray()[ic].get() ;
761 }
762 else {
763 throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
764 }
765 }
766 else {
767 // even
768 //logsink_->postLocally( LogMessage("even case",LogOrigin(className_,funcName,WHERE)) ) ;
769 refpix = 0.5 * ( (double)nchan - 1.0 ) ;
770 int ic = nchan / 2 ;
771 if ( spwrow->isChanWidthExists() ) {
772 incr = spwrow->getChanWidth().get() ;
773 }
774 else if ( spwrow->isChanWidthArrayExists() ) {
775 incr = spwrow->getChanWidthArray()[0].get() ;
776 }
777 else {
778 throw (casa::AipsError( "Either chanWidthArray or chanWidth must exist." )) ;
779 }
780 if ( spwrow->isChanFreqStepExists() ) {
781 if ( spwrow->getChanFreqStep().get() < 0.0 )
782 incr *= -1.0 ;
783 }
784 else if ( spwrow->isChanFreqArrayExists() ) {
785 vector<Frequency> chanFreqArr = spwrow->getChanFreqArray() ;
786 if ( chanFreqArr[0].get() > chanFreqArr[1].get() )
787 incr *= -1.0 ;
788 }
789 else {
790 throw (casa::AipsError( "Either chanFreqArray or chanFreqStep must exist." )) ;
791 }
792 if ( spwrow->isChanFreqStartExists() ) {
793 refval = spwrow->getChanFreqStart().get() + refpix * incr ;
794 }
795 else if ( spwrow->isChanFreqArrayExists() ) {
796 vector<Frequency> freqs = spwrow->getChanFreqArray() ;
797 refval = 0.5 * ( freqs[ic-1].get() + freqs[ic].get() ) ;
798 }
799 else {
800 throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
801 }
802 }
803}
804
805vector<double> ASDMReader::getRestFrequency( unsigned int idx )
806{
807 vector<double> rf( 0 ) ;
808 unsigned int index = dataIdList_[idx] ;
809 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
810 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
811 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
812 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
813 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
814 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
815 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
816 if ( frow->isSourceIdExists() ) {
817 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
818 int sid = frow->getSourceId() ;
819 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
820 if ( srow->isRestFrequencyExists() ) {
821 //logsink_->postLocally( LogMessage("restFrequency exists",LogOrigin(className_,funcName,WHERE)) ) ;
822 vector<Frequency> restfreq = srow->getRestFrequency() ;
823 rf.resize( restfreq.size() ) ;
824 for ( unsigned int i = 0 ; i < restfreq.size() ; i++ )
825 rf[i] = restfreq[i].get() ;
826 }
827 }
828 return rf ;
829}
830
831double ASDMReader::getTime( unsigned int idx )
832{
833 double tsec = vmsData_->v_time[dataIdList_[idx]] ;
834 return tsec * s2d ;
835}
836
837double ASDMReader::getInterval( unsigned int idx )
838{
839 return vmsData_->v_interval[dataIdList_[idx]] ;
840}
841
842string ASDMReader::getSourceName( unsigned int idx )
843{
844 unsigned int index = dataIdList_[idx] ;
845 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
846 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
847 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
848 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
849 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
850 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
851 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
852 string srcname ;
853 if ( frow->isSourceIdExists() ) {
854 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
855 int sid = frow->getSourceId() ;
856 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
857 srcname = srow->getSourceName() ;
858 }
859 else {
860 srcname = frow->getFieldName() ;
861 }
862 return srcname ;
863}
864
865string ASDMReader::getFieldName( unsigned int idx )
866{
867 int fid = vmsData_->v_fieldId[dataIdList_[idx]] ;
868 Tag ftag( fid, TagType::Field ) ;
869 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
870 ostringstream oss ;
871 oss << frow->getFieldName() << "__" << fid ;
872 return oss.str() ;
873}
874
875int ASDMReader::getSrcType( unsigned int scan,
876 unsigned int subscan )
877{
878 int srctype = SrcType::NOTYPE ;
879 Tag ebtag = mainRow_[row_]->getExecBlockId() ;
880 ScanRow *scanrow = asdm_->getScan().getRowByKey( ebtag, (int)scan ) ;
881 ScanIntent scanIntent = scanrow->getScanIntent()[0] ;
882 SubscanRow *subrow = asdm_->getSubscan().getRowByKey( ebtag, (int)scan, (int)subscan ) ;
883 SubscanIntent subIntent = subrow->getSubscanIntent() ;
884 SwitchingMode swmode = subrow->getSubscanMode() ;
885 if ( scanIntent == OBSERVE_TARGET ) {
886 // on sky scan
887 if ( swmode == NO_SWITCHING || swmode == POSITION_SWITCHING ) {
888 // position switching
889 // tentatively set NO_SWITCHING = POSITION_SWITCHING
890 if ( subIntent == ON_SOURCE ) {
891 srctype = SrcType::PSON ;
892 }
893 else if ( subIntent == OFF_SOURCE ) {
894 srctype = SrcType::PSOFF ;
895 }
896 }
897 else if ( swmode == FREQUENCY_SWITCHING ) {
898 // frequency switching
899 if ( subIntent == ON_SOURCE ) {
900 srctype = SrcType::FSON ;
901 }
902 else if ( subIntent == OFF_SOURCE ) {
903 srctype = SrcType::FSOFF ;
904 }
905 }
906 else if ( swmode == NUTATOR_SWITCHING ) {
907 // nutator switching
908 if ( subIntent == ON_SOURCE ) {
909 srctype = SrcType::PSON ;
910 }
911 else if ( subIntent == OFF_SOURCE ) {
912 srctype = SrcType::PSOFF ;
913 }
914 }
915 else {
916 // other switching mode
917 if ( subIntent == ON_SOURCE ) {
918 srctype = SrcType::SIG ;
919 }
920 else if ( subIntent == OFF_SOURCE ) {
921 srctype = SrcType::REF ;
922 }
923 }
924 }
925 else {
926 // off sky (e.g. calibrator device) scan
927 if ( subIntent == ON_SOURCE ) {
928 srctype = SrcType::SIG ;
929 }
930 else if ( subIntent == OFF_SOURCE ) {
931 srctype = SrcType::REF ;
932 }
933 else if ( subIntent == HOT ) {
934 srctype = SrcType::HOT ;
935 }
936 else if ( subIntent == AMBIENT ) {
937 srctype = SrcType::SKY ;
938 }
939 else {
940 srctype = SrcType::CAL ;
941 }
942 }
943
944 return srctype ;
945}
946
947unsigned int ASDMReader::getSubscanNo( unsigned int idx )
948{
949 //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)) ) ;
950 return vmsData_->v_msState[dataIdList_[idx]].subscanNum ;
951}
952
953vector<double> ASDMReader::getSourceDirection( unsigned int idx )
954{
955 vector<double> dir( 2, 0.0 ) ;
956 unsigned int index = dataIdList_[idx] ;
957 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
958 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
959 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
960 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
961 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
962 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
963 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
964 string srcname ;
965 if ( frow->isSourceIdExists() ) {
966 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
967 int sid = frow->getSourceId() ;
968 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
969 vector<Angle> srcdir = srow->getDirection() ;
970 if ( srow->isDirectionCodeExists() ) {
971 DirectionReferenceCode dircode = srow->getDirectionCode() ;
972 if ( dircode != J2000 ) {
973 // if not J2000, need direction conversion
974 }
975 }
976 dir[0] = srcdir[0].get() ;
977 dir[1] = srcdir[1].get() ;
978 }
979 return dir ;
980}
981
982void ASDMReader::getSourceDirection( unsigned int idx,
983 vector<double> &dir,
984 string &ref )
985{
986 dir.resize( 2 ) ;
987 dir[0] = 0.0 ;
988 dir[1] = 0.0 ;
989 ref = "J2000" ;
990 unsigned int index = dataIdList_[idx] ;
991 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
992 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
993 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
994 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
995 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
996 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
997 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
998 string srcname ;
999 if ( frow->isSourceIdExists() ) {
1000 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1001 int sid = frow->getSourceId() ;
1002 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1003 vector<Angle> srcdir = srow->getDirection() ;
1004 if ( srow->isDirectionCodeExists() ) {
1005 ref = CDirectionReferenceCode::toString( srow->getDirectionCode() ) ;
1006 }
1007 dir[0] = srcdir[0].get() ;
1008 dir[1] = srcdir[1].get() ;
1009 }
1010}
1011
1012void ASDMReader::getSourceDirection( vector<double> &dir, string &ref )
1013{
1014 dir.resize( 2 ) ;
1015 ref = "J2000" ; // default is J2000
1016 SourceTable &tab = asdm_->getSource() ;
1017 SourceRow *row = tab.get()[0] ;
1018 vector<Angle> dirA = row->getDirection() ;
1019 dir[0] = dirA[0].get() ;
1020 dir[1] = dirA[1].get() ;
1021 if ( row->isDirectionCodeExists() ) {
1022 ref = CDirectionReferenceCode::toString( row->getDirectionCode() ) ;
1023 }
1024}
1025
1026vector<double> ASDMReader::getSourceProperMotion( unsigned int idx )
1027{
1028 vector<double> pm( 2, 0.0 ) ;
1029 unsigned int index = dataIdList_[idx] ;
1030 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1031 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1032 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1033 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1034 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1035 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1036 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1037 string srcname ;
1038 if ( frow->isSourceIdExists() ) {
1039 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1040 int sid = frow->getSourceId() ;
1041 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1042 vector<AngularRate> srcpm = srow->getProperMotion() ;
1043 pm[0] = srcpm[0].get() ;
1044 pm[1] = srcpm[1].get() ;
1045 }
1046 return pm ;
1047}
1048
1049double ASDMReader::getSysVel( unsigned int idx )
1050{
1051 double sysvel = 0.0 ;
1052 unsigned int index = dataIdList_[idx] ;
1053 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1054 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1055 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1056 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1057 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1058 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1059 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1060 string srcname ;
1061 if ( frow->isSourceIdExists() ) {
1062 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1063 int sid = frow->getSourceId() ;
1064 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1065 if ( srow->isSysVelExists() ) {
1066 vector<Speed> sysvelV = srow->getSysVel() ;
1067 if ( sysvelV.size() > 0 )
1068 sysvel = sysvelV[0].get() ;
1069 }
1070 }
1071 return sysvel ;
1072}
1073
1074unsigned int ASDMReader::getFlagRow( unsigned int idx )
1075{
1076 return vmsData_->v_flag[dataIdList_[idx]] ;
1077}
1078
1079vector<unsigned int> ASDMReader::getDataShape( unsigned int idx )
1080{
1081 return vmsData_->vv_dataShape[dataIdList_[idx]] ;
1082}
1083
1084float * ASDMReader::getSpectrum( unsigned int idx )
1085{
1086 map<AtmPhaseCorrection, float*> data = vmsData_->v_m_data[dataIdList_[idx]] ;
1087 //map<AtmPhaseCorrection, float*>::iterator iter = data.find(AP_UNCORRECTED) ;
1088 map<AtmPhaseCorrection, float*>::iterator iter = data.find(apc_) ;
1089 float *autoCorr = iter->second ;
1090 return autoCorr ;
1091}
1092
1093// bool * ASDMReader::getFlagChannel( unsigned int idx )
1094// {
1095// return 0 ;
1096// }
1097
1098vector< vector<float> > ASDMReader::getTsys( unsigned int idx )
1099{
1100 unsigned int index = dataIdList_[idx] ;
1101 Tag anttag( antennaId_, TagType::Antenna ) ;
1102 int feedid = vmsData_->v_feedId1[index] ;
1103 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1104 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1105 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1106 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1107 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1108 //int nchan = asdm_->getSpectralWindow().getRowByKey(spwtag)->getNumChan() ;
1109 vector< vector<float> > defaultTsys( 1, vector<float>( 1, 1.0 ) ) ;
1110 SysCalTable &sctab = asdm_->getSysCal() ;
1111 if ( sctab.size() == 0 )
1112 return defaultTsys ;
1113 SysCalRow *scrow = sctab.getRowByKey( anttag, spwtag, tint, feedid ) ;
1114 if ( scrow->isTsysSpectrumExists() ) {
1115 vector< vector<Temperature> > tsysSpec = scrow->getTsysSpectrum() ;
1116 unsigned int numReceptor = tsysSpec.size() ;
1117 unsigned int numChan = tsysSpec[0].size() ;
1118 vector< vector<float> > tsys( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1119 for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1120 for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1121 tsys[ir][ic] = tsysSpec[ir][ic].get() ;
1122 }
1123 }
1124 return tsys ;
1125 }
1126// else if ( scrow->isTsysExists() ) {
1127// vector<Temperature> tsysScalar = scrow->getTsys() ;
1128// unsigned int numReceptor = tsysScalar.size() ;
1129// vector< vector<float> > tsys( numReceptor ) ;
1130// for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1131// tsys[ir] = vector<float>( 1, tsysScalar[ir].get() ) ;
1132// return tsys ;
1133// }
1134 else {
1135 return defaultTsys ;
1136 }
1137}
1138
1139vector< vector<float> > ASDMReader::getTcal( unsigned int idx )
1140{
1141 unsigned int index = dataIdList_[idx] ;
1142 Tag anttag( antennaId_, TagType::Antenna ) ;
1143 int feedid = vmsData_->v_feedId1[index] ;
1144 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1145 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1146 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1147 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1148 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1149 //int nchan = asdm_->getSpectralWindow().getRowByKey(spwtag)->getNumChan() ;
1150 vector< vector<float> > defaultTcal( 1, vector<float>( 1, 1.0 ) ) ;
1151 SysCalTable &sctab = asdm_->getSysCal() ;
1152 if ( sctab.size() == 0 )
1153 return defaultTcal ;
1154 SysCalRow *scrow = sctab.getRowByKey( anttag, spwtag, tint, feedid ) ;
1155 if ( scrow->isTcalSpectrumExists() ) {
1156 vector< vector<Temperature> > tcalSpec = scrow->getTcalSpectrum() ;
1157 unsigned int numReceptor = tcalSpec.size() ;
1158 unsigned int numChan = tcalSpec[0].size() ;
1159 vector< vector<float> > tcal( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1160 for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1161 for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1162 tcal[ir][ic] = tcalSpec[ir][ic].get() ;
1163 }
1164 }
1165 return tcal ;
1166 }
1167// else if ( scrow->isTcalExists() ) {
1168// vector<Temperature> tcalScalar = scrow->getTcal() ;
1169// unsigned int numReceptor = tcalScalar.size() ;
1170// vector< vector<float> > tcal( numReceptor, vector<float>( 1, 1.0 ) ) ;
1171// for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1172// tcal[ir][0] = tcalScalar[ir][0].get() ;
1173// return tcal ;
1174// }
1175 else {
1176 return defaultTcal ;
1177 }
1178}
1179
1180vector<float> ASDMReader::getOpacity( unsigned int idx )
1181{
1182 vector<float> tau(0) ;
1183 CalAtmosphereTable &atmtab = asdm_->getCalAtmosphere() ;
1184 unsigned int nrow = atmtab.size() ;
1185 if ( nrow > 0 ) {
1186 unsigned int index = dataIdList_[idx] ;
1187 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1188 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1189 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1190 //int feedid = vmsData_->v_feedId1[index] ;
1191 //Tag atag( antennaId_, TagType::Antenna ) ;
1192 //Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1193 //DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey(ddtag) ;
1194 //Tag spwtag = ddrow->getSpectralWindowId() ;
1195 //SpectralWindowRow *spwrow = ddrow->getSpectralWindowUsingSpectralWindowId() ;
1196 //BasebandName bbname = spwrow->getBasebandName() ;
1197 //FeedRow *frow = asdm_->getFeed().getRowByKey( atag, spwtag, tint, feedid ) ;
1198 //int nfeed = frow->getNumReceptor() ;
1199 //vector<ReceiverRow *> rrows = frow->getReceivers() ;
1200 vector<CalAtmosphereRow *> atmrows = atmtab.get() ;
1201 //ReceiverBand rb = rrows[0]->getFrequencyBand() ;
1202 int row0 = -1 ;
1203 double eps = tint.getStart().getMJD() ;
1204 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1205 CalAtmosphereRow *atmrow = atmrows[irow] ;
1206 if ( casa::String(atmrow->getAntennaName()) != antennaName_
1207 //|| atmrow->getReceiverBand() != rb
1208 //|| atmrow->getBasebandName() != bbname
1209 || atmrow->getCalDataUsingCalDataId()->getCalType() != CAL_ATMOSPHERE )
1210 continue ;
1211 else {
1212 double dt = tint.getStart().getMJD() - atmrow->getEndValidTime().getMJD() ;
1213 if ( dt >= 0 && dt < eps ) {
1214 eps = dt ;
1215 row0 = (int)irow ;
1216 }
1217 }
1218 }
1219 if ( row0 != -1 ) {
1220 CalAtmosphereRow *atmrow = atmrows[row0] ;
1221 tau = atmrow->getTau() ;
1222 }
1223 }
1224 return tau ;
1225}
1226
1227void ASDMReader::getWeatherInfo( unsigned int idx,
1228 float &temperature,
1229 float &pressure,
1230 float &humidity,
1231 float &windspeed,
1232 float &windaz )
1233{
1234 temperature = 0.0 ;
1235 pressure = 0.0 ;
1236 humidity = 0.0 ;
1237 windspeed = 0.0 ;
1238 windaz = 0.0 ;
1239
1240 //logsink_->postLocally( LogMessage("weatherStationId_ = "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1241
1242 WeatherTable &wtab = asdm_->getWeather() ;
1243 if ( wtab.size() == 0 || weatherStationId_ == -1 )
1244 return ;
1245
1246 unsigned int index = dataIdList_[idx] ;
1247 //Tag anttag( antennaId_, TagType::Antenna ) ;
1248 //Tag sttag = (asdm_->getAntenna().getRowByKey( anttag ))->getStationId() ;
1249 Tag sttag( (unsigned int)weatherStationId_, TagType::Station ) ;
1250 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1251 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1252 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1253 //WeatherRow *wrow = wtab.getRowByKey( sttag, tint ) ;
1254 vector<WeatherRow *> *wrows = wtab.getByContext( sttag ) ;
1255 WeatherRow *wrow = (*wrows)[0] ;
1256 unsigned int nrow = wrows->size() ;
1257 //logsink_->postLocally( LogMessage("There are "+String::toString(nrow)+" rows for given context: stationId "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1258 ArrayTime startTime = tint.getStart() ;
1259 if ( startTime < (*wrows)[0]->getTimeInterval().getStart() ) {
1260 temperature = (*wrows)[0]->getTemperature().get() ;
1261 pressure = (*wrows)[0]->getPressure().get() ;
1262 humidity = (*wrows)[0]->getRelHumidity().get() ;
1263 windspeed = (*wrows)[0]->getWindSpeed().get() ;
1264 windaz = (*wrows)[0]->getWindDirection().get() ;
1265 }
1266 else if ( startTime > (*wrows)[nrow-1]->getTimeInterval().getStart() ) {
1267 temperature = (*wrows)[nrow-1]->getTemperature().get() ;
1268 pressure = (*wrows)[nrow-1]->getPressure().get() ;
1269 humidity = (*wrows)[nrow-1]->getRelHumidity().get() ;
1270 windspeed = (*wrows)[nrow-1]->getWindSpeed().get() ;
1271 windaz = (*wrows)[nrow-1]->getWindDirection().get() ;
1272 }
1273 else {
1274 for ( unsigned int irow = 1 ; irow < wrows->size() ; irow++ ) {
1275 wrow = (*wrows)[irow-1] ;
1276 if ( startTime < (*wrows)[irow]->getTimeInterval().getStart() ) {
1277 temperature = wrow->getTemperature().get() ;
1278 pressure = wrow->getPressure().get() ;
1279 humidity = wrow->getRelHumidity().get() ;
1280 windspeed = wrow->getWindSpeed().get() ;
1281 windaz = wrow->getWindDirection().get() ;
1282 break ;
1283 }
1284 }
1285 }
1286
1287 return ;
1288}
1289
1290void ASDMReader::processStation()
1291{
1292 antennaPad_.resize(0) ;
1293 weatherStation_.resize(0) ;
1294 StationTable &stab = asdm_->getStation() ;
1295 vector<StationRow *> srows = stab.get() ;
1296 for ( unsigned int irow = 0 ; irow < srows.size() ; irow++ ) {
1297 StationType stype = srows[irow]->getType() ;
1298 Tag stag = srows[irow]->getStationId() ;
1299 if ( stype == ANTENNA_PAD )
1300 antennaPad_.push_back( stag ) ;
1301 else if ( stype == WEATHER_STATION )
1302 weatherStation_.push_back( stag ) ;
1303 }
1304
1305 weatherStationId_ = getClosestWeatherStation() ;
1306}
1307
1308int ASDMReader::getClosestWeatherStation()
1309{
1310 if ( weatherStation_.size() == 0 )
1311 return -1 ;
1312
1313 Tag atag( antennaId_, TagType::Antenna ) ;
1314 Tag stag = (asdm_->getAntenna().getRowByKey( atag ))->getStationId() ;
1315 vector<double> apos( 3 ) ;
1316 StationTable &stab = asdm_->getStation() ;
1317 StationRow *srow = stab.getRowByKey( stag ) ;
1318 vector<Length> pos = srow->getPosition() ;
1319 apos[0] = pos[0].get() ;
1320 apos[1] = pos[1].get() ;
1321 apos[2] = pos[2].get() ;
1322
1323 double eps = 1.0e20 ;
1324 int retval = -1 ;
1325 for ( unsigned int ir = 0 ; ir < weatherStation_.size() ; ir++ ) {
1326 srow = stab.getRowByKey( weatherStation_[ir] ) ;
1327 vector<Length> wpos = srow->getPosition() ;
1328 double dist = (apos[0]-wpos[0].get())*(apos[0]-wpos[0].get())
1329 + (apos[1]-wpos[1].get())*(apos[1]-wpos[1].get())
1330 + (apos[2]-wpos[2].get())*(apos[2]-wpos[2].get()) ;
1331 if ( dist < eps ) {
1332 retval = (int)(weatherStation_[ir].getTagValue()) ;
1333 }
1334 }
1335
1336 return retval ;
1337}
1338
1339void ASDMReader::getPointingInfo( unsigned int idx,
1340 vector<double> &dir,
1341 double &az,
1342 double &el,
1343 vector<double> &srate )
1344{
1345 String funcName = "getPointingInfo" ;
1346
1347 dir.resize(0) ;
1348 az = -1.0 ;
1349 el = -1.0 ;
1350 srate.resize(0) ;
1351
1352 Tag atag( antennaId_, TagType::Antenna ) ;
1353 unsigned int index = dataIdList_[idx] ;
1354 vector<PointingRow *> *prows = asdm_->getPointing().getByContext( atag ) ;
1355 PointingRow *prow ;
1356 PointingRow *qrow ;
1357 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1358 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1359 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1360
1361 unsigned int nrow = prows->size() ;
1362 //logsink_->postLocally( LogMessage("There are " << nrow << " rows for given context: antennaId "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
1363
1364// for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1365// prow = (*prows)[irow] ;
1366// ArrayTimeInterval ati = prow->getTimeInterval() ;
1367// ArrayTime pst = ati.getStart() ;
1368// ArrayTime pet( ati.getStartInMJD()+ati.getDurationInDays() ) ;
1369// logsink_->postLocally( LogMessage("start: "+pst.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1370// logsink_->postLocally( LogMessage("end: "+pet.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1371// }
1372
1373 if ( nrow == 0 )
1374 return ;
1375
1376 srate.resize(2) ;
1377 srate[0] = 0.0 ;
1378 srate[1] = 0.0 ;
1379 az = 0.0 ;
1380 el = 0.0 ;
1381 double tcen = 0.0 ;
1382
1383 //
1384 // shape of pointingDirection is (numSample,2) if usePolynomial = False, while it is
1385 // (numTerm,2) if usePolynomial = True.
1386 //
1387 // In the former case, typical sampled time interval is 48msec, which is very small
1388 // compared with typical integration time (~a few sec).
1389 // Scan rate for this case is always [0,0] (or get slope?).
1390 //
1391 // In the later case, scan rate is (pointingDirection[1][0],pointingDirection[1][1])
1392 //
1393 // PointingDirection seems to store AZEL
1394 //
1395 ArrayTimeInterval pTime0 = (*prows)[0]->getTimeInterval() ;
1396 ArrayTimeInterval pTime1 = (*prows)[nrow-1]->getTimeInterval() ;
1397 if ( tint.getStartInMJD()+tint.getDurationInDays() < pTime0.getStartInMJD() ) {
1398 logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint < ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1399 prow = (*prows)[0] ;
1400 //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1401 vector< vector<Angle> > dirA = prow->getTarget() ;
1402 vector< vector<Angle> > offA = prow->getOffset() ;
1403 //az = dirA[0][0].get() ;
1404 //el = dirA[0][1].get() ;
1405 az = dirA[0][0].get() + offA[0][0].get() ;
1406 el = dirA[0][1].get() + offA[0][1].get() ;
1407 if ( prow->getUsePolynomials() && prow->getNumTerm() > 1 ) {
1408 //srate[0] = dirA[1][0].get() ;
1409 //srate[1] = dirA[1][1].get() ;
1410 srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1411 srate[1] = dirA[1][1].get() + offA[1][1].get() ;
1412 }
1413 }
1414 else if ( tint.getStartInMJD() > pTime1.getStartInMJD()+pTime1.getDurationInDays() ) {
1415 logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint > ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1416 prow = (*prows)[nrow-1] ;
1417 int numSample = prow->getNumSample() ;
1418 //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1419 vector< vector<Angle> > dirA = prow->getTarget() ;
1420 vector< vector<Angle> > offA = prow->getOffset() ;
1421 if ( prow->getUsePolynomials() ) {
1422 //az = dirA[0][0].get() ;
1423 //el = dirA[0][1].get() ;
1424 az = dirA[0][0].get() + offA[0][0].get() ;
1425 el = dirA[0][1].get() + offA[0][1].get() ;
1426 if ( prow->getNumTerm() > 1 ) {
1427 //srate[0] = dirA[1][0].get() ;
1428 //srate[1] = dirA[1][1].get() ;
1429 srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1430 srate[1] = dirA[1][1].get() + offA[1][0].get() ;
1431 }
1432 }
1433 else {
1434 //az = dirA[numSample-1][0].get() ;
1435 //el = dirA[numSample-1][1].get() ;
1436 az = dirA[numSample-1][0].get() + offA[numSample-1][0].get() ;
1437 el = dirA[numSample-1][1].get() + offA[numSample-1][1].get() ;
1438 }
1439 }
1440 else {
1441 ArrayTime startTime = tint.getStart() ;
1442 ArrayTime endTime( tint.getStartInMJD()+tint.getDurationInDays() ) ;
1443 int row0 = -1 ;
1444 int row1 = -1 ;
1445 int row2 = -1 ;
1446 double dt0 = getMidTime( tint ).getMJD() ;
1447 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1448 prow = (*prows)[irow] ;
1449 double dt = getMidTime( tint ).getMJD() - getMidTime( prow->getTimeInterval() ).getMJD() ;
1450 if ( dt > 0 && dt < dt0 ) {
1451 dt0 = dt ;
1452 row2 = irow ;
1453 }
1454 if ( prow->getTimeInterval().contains( startTime ) )
1455 row0 = irow ;
1456 else if ( prow->getTimeInterval().contains( endTime ) )
1457 row1 = irow ;
1458 if ( row0 != -1 && row1 != -1 )
1459 break ;
1460 }
1461 //logsink_->postLocally( LogMessage("row0 = "+String::toString(row0)+", row1 = "+String::toString(row1)+", row2 = "+String::toString(row2),LogOrigin(className_,funcName,WHERE)) ) ;
1462 if ( row0 == -1 && row1 == -1 ) {
1463 prow = (*prows)[row2] ;
1464 qrow = (*prows)[row2+1] ;
1465 double t0 = getEndTime( prow->getTimeInterval() ).getMJD() ;
1466 double t1 = qrow->getTimeInterval().getStartInMJD() ;
1467 double t = startTime.getMJD() ;
1468 //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1469 //vector< vector<Angle> > dirB = qrow->getPointingDirection() ;
1470 vector< vector<Angle> > dirA = prow->getTarget() ;
1471 vector< vector<Angle> > offA = prow->getOffset() ;
1472 vector< vector<Angle> > dirB = qrow->getTarget() ;
1473 vector< vector<Angle> > offB = qrow->getOffset() ;
1474 //double da0 = dirA[0][0].get() ;
1475 //double db0 = dirB[0][0].get() ;
1476 //double da1 = dirA[0][1].get() ;
1477 //double db1 = dirB[0][1].get() ;
1478 double da0 = dirA[0][0].get() + offA[0][0].get() ;
1479 double db0 = dirB[0][0].get() + offB[0][0].get() ;
1480 double da1 = dirA[0][1].get() + offA[0][1].get() ;
1481 double db1 = dirB[0][1].get() + offB[0][1].get() ;
1482 if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1483 double dt = ( t - t0 ) / ( t1 - t0 ) ;
1484 az = da0 + ( db0 - da0 ) * dt ;
1485 el = da1 + ( db1 - da1 ) * dt ;
1486 if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1487 //double ra0 = dirA[1][0].get() ;
1488 //double rb0 = dirB[1][0].get() ;
1489 //double ra1 = dirA[1][1].get() ;
1490 //double rb1 = dirB[1][1].get() ;
1491 double ra0 = dirA[1][0].get() + offA[1][0].get() ;
1492 double rb0 = dirB[1][0].get() + offB[1][0].get() ;
1493 double ra1 = dirA[1][1].get() + offA[1][1].get() ;
1494 double rb1 = dirB[1][1].get() + offB[1][1].get() ;
1495 srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1496 srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1497 }
1498 }
1499 else if ( !(qrow->getUsePolynomials()) ) {
1500 double dt = ( t - t0 ) / ( t1 - t0 ) ;
1501 az = da0 + ( db0 - da0 ) * dt ;
1502 el = da1 + ( db1 - da1 ) * dt ;
1503 }
1504 //else {
1505 // nothing to do
1506 //}
1507 }
1508 else {
1509 int ndir = 0 ;
1510 if ( row0 == -1 ) {
1511 row0 = row1 ;
1512 }
1513 else if ( row1 == -1 ) {
1514 row1 = row0 ;
1515 }
1516 prow = (*prows)[row0] ;
1517 qrow = (*prows)[row1] ;
1518 if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1519 //logsink_->postLocally( LogMessage("usePolynomial = True",LogOrigin(className_,funcName,WHERE)) ) ;
1520 if ( row0 == row1 ) {
1521 prow = (*prows)[row0] ;
1522 //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1523 vector< vector<Angle> > dirA = prow->getTarget() ;
1524 vector< vector<Angle> > offA = prow->getOffset() ;
1525 //az = dirA[0][0].get() ;
1526 //el = dirA[0][1].get() ;
1527 az = dirA[0][0].get() + offA[0][0].get() ;
1528 el = dirA[0][1].get() + offA[0][1].get() ;
1529 if ( prow->getNumTerm() > 1 ) {
1530 //srate[0] = dirA[1][0].get() ;
1531 //srate[1] = dirA[1][1].get() ;
1532 srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1533 srate[1] = dirA[1][1].get() + offA[1][0].get() ;
1534 }
1535 }
1536 else {
1537 prow = (*prows)[row0] ;
1538 qrow = (*prows)[row1] ;
1539 //vector< vector<Angle> > dirA = qrow->getPointingDirection() ;
1540 //vector< vector<Angle> > dirB = qrow->getPointingDirection() ;
1541 vector< vector<Angle> > dirA = qrow->getTarget() ;
1542 vector< vector<Angle> > offA = prow->getOffset() ;
1543 vector< vector<Angle> > dirB = qrow->getTarget() ;
1544 vector< vector<Angle> > offB = qrow->getOffset() ;
1545 double t0 = getMidTime( prow->getTimeInterval() ).getMJD() ;
1546 double t1 = getMidTime( qrow->getTimeInterval() ).getMJD() ;
1547 double t = startTime.getMJD() ;
1548 double dt = ( t - t0 ) / ( t1 - t0 ) ;
1549 //double da0 = dirA[0][0].get() ;
1550 //double db0 = dirB[0][0].get() ;
1551 //double da1 = dirA[0][1].get() ;
1552 //double db1 = dirB[0][1].get() ;
1553 double da0 = dirA[0][0].get() + offA[0][0].get() ;
1554 double db0 = dirB[0][0].get() + offB[0][0].get() ;
1555 double da1 = dirA[0][1].get() + offA[0][1].get() ;
1556 double db1 = dirB[0][1].get() + offB[0][1].get() ;
1557 az = da0 + ( db0 - da0 ) * dt ;
1558 el = da1 + ( db1 - db0 ) * dt ;
1559 if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1560 //double ra0 = dirA[1][0].get() ;
1561 //double rb0 = dirB[1][0].get() ;
1562 //double ra1 = dirA[1][1].get() ;
1563 //double rb1 = dirB[1][1].get() ;
1564 double ra0 = dirA[1][0].get() + offA[1][0].get() ;
1565 double rb0 = dirB[1][0].get() + offB[1][0].get() ;
1566 double ra1 = dirA[1][1].get() + offA[1][1].get() ;
1567 double rb1 = dirB[1][1].get() + offB[1][1].get() ;
1568 srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1569 srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1570 }
1571 }
1572 }
1573 else if ( prow->getUsePolynomials() == qrow->getUsePolynomials() ) {
1574 //logsink_->postLocally( LogMessage("numSample == numTerm",LogOrigin(className_,funcName,WHERE)) ) ;
1575 for ( int irow = row0 ; irow <= row1 ; irow++ ) {
1576 prow = (*prows)[irow] ;
1577 int numSample = prow->getNumSample() ;
1578 //logsink_->postLocally( LogMessage("numSample = "+String::toString(numSample),LogOrigin(className_,funcName,WHERE)) ) ;
1579 //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1580 vector< vector<Angle> > dirA = prow->getTarget() ;
1581 vector< vector<Angle> > offA = prow->getOffset() ;
1582 if ( prow->isSampledTimeIntervalExists() ) {
1583 //logsink_->postLocally( LogMessage("sampledTimeIntervalExists",LogOrigin(className_,funcName,WHERE)) ) ;
1584 vector<ArrayTimeInterval> stime = prow->getSampledTimeInterval() ;
1585 for ( int isam = 0 ; isam < numSample ; isam++ ) {
1586 //if ( tint.overlaps( stime[isam] ) ) {
1587 if ( tint.contains( stime[isam] ) ) {
1588 //az += dirA[isam][0].get() ;
1589 //el += dirA[isam][1].get() ;
1590 az += dirA[isam][0].get() + offA[isam][0].get() ;
1591 el += dirA[isam][1].get() + offA[isam][1].get() ;
1592 ndir++ ;
1593 }
1594 }
1595 }
1596 else {
1597 double sampleStart = prow->getTimeInterval().getStartInMJD() ;
1598 double sampleInterval = prow->getTimeInterval().getDurationInDays() / (double)numSample ;
1599 //logsink_->postLocally( LogMessage("sampleStart = "+String::toString(sampleStart),LogOrigin(className_,funcName,WHERE)) )
1600 //logsink_->postLocally( LogMessage("sampleInterval = "+String::toString(sampleInterval),LogOrigin(className_,funcName,WHERE)) ) ;
1601 //logsink_->postLocally( LogMessage("tint = "+tint.toString(),LogOrigin(className_,funcName,WHERE)) ) ;
1602 for ( int isam = 0 ; isam < numSample ; isam++ ) {
1603 ArrayTimeInterval stime( sampleStart+isam*sampleInterval, sampleInterval ) ;
1604 //if ( tint.overlaps( stime ) ) {
1605 if ( tint.contains( stime ) ) {
1606 //az += dirA[isam][0].get() ;
1607 //el += dirA[isam][1].get() ;
1608 az += dirA[isam][0].get() + offA[isam][0].get() ;
1609 el += dirA[isam][1].get() + offA[isam][1].get() ;
1610 tcen += getMidTime( stime ).getMJD() ;
1611 ndir++ ;
1612 }
1613 }
1614 }
1615 }
1616 if ( ndir > 1 ) {
1617 az /= (double)ndir ;
1618 el /= (double)ndir ;
1619 tcen /= (double)ndir ;
1620 }
1621 }
1622 //else {
1623 // nothing to do
1624 //}
1625 }
1626
1627 AntennaRow *arow = asdm_->getAntenna().getRowByKey( Tag( antennaId_, TagType::Antenna ) ) ;
1628 StationRow *srow = arow->getStationUsingStationId() ;
1629 vector<Length> antposL = srow->getPosition() ;
1630 casa::Vector<casa::Double> antpos( 3 ) ;
1631 for ( int i = 0 ; i < 3 ; i++ )
1632 antpos[i] = antposL[i].get() ;
1633 //logsink_->postLocally( LogMessage("tcen = "+String::toString(tcen),LogOrigin(className_,funcName,WHERE)) ) ;
1634 //logsink_->postLocally( LogMessage("antpos = "+String::toString(antpos),LogOrigin(className_,funcName,WHERE)) ) ;
1635 toJ2000( dir, az, el, tcen, antpos ) ;
1636
1637 }
1638
1639 return ;
1640}
1641
1642ArrayTime ASDMReader::getMidTime( const ArrayTimeInterval &t )
1643{
1644 return ArrayTime( t.getStartInMJD() + 0.5 * t.getDurationInDays() ) ;
1645}
1646
1647ArrayTime ASDMReader::getEndTime( const ArrayTimeInterval &t )
1648{
1649 return ArrayTime( t.getStartInMJD() + t.getDurationInDays() ) ;
1650}
1651
1652ArrayTime ASDMReader::getStartTime( const ArrayTimeInterval &t )
1653{
1654 return ArrayTime( t.getStartInMJD() ) ;
1655}
1656
1657void ASDMReader::toJ2000( vector<double> &dir,
1658 double az,
1659 double el,
1660 double mjd,
1661 casa::Vector<casa::Double> antpos )
1662{
1663 String funcName = "toJ2000" ;
1664
1665// casa::Vector<casa::Double> azel( 2 ) ;
1666// azel[0] = az ;
1667// azel[1] = el ;
1668// casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
1669// casa::Vector<casa::Quantity> qantpos( 3 ) ;
1670// qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
1671// qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
1672// qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
1673// casa::MPosition mp( casa::MVPosition( qantpos ),
1674// casa::MPosition::ITRF ) ;
1675// //ostringstream oss ;
1676// //mp.print( oss ) ;
1677// //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
1678
1679// casa::MeasFrame mf( me, mp ) ;
1680// casa::MDirection::Convert toj2000( casa::MDirection::AZELGEO,
1681// casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
1682// casa::Vector<casa::Double> cdir = toj2000( azel ).getAngle( "rad" ).getValue() ;
1683// //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
1684// dir.resize(2) ;
1685// dir[0] = (double)(cdir[0]) ;
1686// dir[1] = (double)(cdir[1]) ;
1687 vector<double> azel( 2 ) ;
1688 azel[0] = az ;
1689 azel[1] = el ;
1690 dir = toJ2000( azel, "AZELGEO", mjd, antpos ) ;
1691}
1692
1693vector<double> ASDMReader::toJ2000( vector<double> dir,
1694 casa::String dirref,
1695 double mjd,
1696 casa::Vector<casa::Double> antpos )
1697{
1698 casa::String funcname = "toJ2000" ;
1699
1700 vector<double> newd( dir ) ;
1701 if ( dirref != "J2000" ) {
1702 casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
1703 casa::Vector<casa::Quantity> qantpos( 3 ) ;
1704 qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
1705 qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
1706 qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
1707 casa::MPosition mp( casa::MVPosition( qantpos ),
1708 casa::MPosition::ITRF ) ;
1709 //ostringstream oss ;
1710 //mp.print( oss ) ;
1711 //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
1712
1713 casa::MeasFrame mf( me, mp ) ;
1714 casa::MDirection::Types dirtype ;
1715 casa::Bool b = casa::MDirection::getType( dirtype, dirref ) ;
1716 if ( b ) {
1717 casa::MDirection::Convert toj2000( dirtype,
1718 casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
1719 casa::Vector<casa::Double> cdir = toj2000( dir ).getAngle( "rad" ).getValue() ;
1720 //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
1721 newd[0] = (double)(cdir[0]) ;
1722 newd[1] = (double)(cdir[1]) ;
1723 }
1724 }
1725 return newd ;
1726}
1727
1728void ASDMReader::setLogger( CountedPtr<LogSinkInterface> &logsink )
1729{
1730 logsink_ = logsink ;
1731}
1732
1733string ASDMReader::getFrame()
1734{
1735 String funcName = "getFrame" ;
1736
1737 // default is TOPO
1738 string frame = "TOPO" ;
1739
1740 SpectralWindowTable &spwtab = asdm_->getSpectralWindow() ;
1741 vector<SpectralWindowRow *> rows = spwtab.get() ;
1742 vector<FrequencyReferenceCode> measFreqRef( rows.size() ) ;
1743 int nref = 0 ;
1744 for ( unsigned int irow = 0 ; irow < rows.size() ; irow++ ) {
1745 int nchan = rows[irow]->getNumChan() ;
1746 if ( nchan != 4 ) {
1747 if ( rows[irow]->isMeasFreqRefExists() ) {
1748 measFreqRef[nref] = rows[irow]->getMeasFreqRef() ;
1749 nref++ ;
1750 }
1751 }
1752 }
1753 if ( nref != 0 ) {
1754 frame = CFrequencyReferenceCode::toString( measFreqRef[0] ) ;
1755 }
1756
1757 //logsink_->postLocally( LogMessage("frame = "+String::toString(frame),LogOrigin(className_,funcName,WHERE)) ) ;
1758
1759 return frame ;
1760}
1761
1762int ASDMReader::getNumIFs()
1763{
1764 String funcName = "getNumIFs" ;
1765
1766 int nif = 0 ;
1767 vector<SpectralWindowRow *> rows = asdm_->getSpectralWindow().get() ;
1768 unsigned int nrow = rows.size() ;
1769 // check if all rows have freqGroup attribute
1770 bool freqGroupExists = true ;
1771 bool countedWvr = false ;
1772 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1773 freqGroupExists &= rows[irow]->isFreqGroupExists() ;
1774 if ( rows[irow]->getNumChan() == 4 ) {
1775 if ( !countedWvr ) {
1776 countedWvr = true ;
1777 nif++ ;
1778 }
1779 }
1780 else {
1781 nif++ ;
1782 }
1783 }
1784
1785 if ( freqGroupExists ) {
1786 vector<int> freqGroup(0) ;
1787 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1788 int fg = rows[irow]->getFreqGroup() ;
1789 if ( (int)count( freqGroup.begin(), freqGroup.end(), fg ) == 0 ) {
1790 freqGroup.push_back( fg ) ;
1791 }
1792 }
1793 nif = freqGroup.size() ;
1794 }
1795
1796 //logsink_->postLocally( LogMessage("nif = "+String::toString(nif),LogOrigin(className_,funcName,WHERE)) ) ;
1797
1798 return nif ;
1799}
Note: See TracBrowser for help on using the repository browser.