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

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

New Development: No

JIRA Issue: Yes CAS-1913

Ready for Test: Yes/No

Interface Changes: Yes/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...

Accept multi spectral resolution type (e.g. fr,ca).
Separator changed from '+' to ',' for -srt and -corr-mode options.
File size: 61.9 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() ) ) ;
215 logsink_->postLocally( LogMessage( "antennaName_ = "+antennaName_, LogOrigin() ) ) ;
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
297 // nif
298 nif = spwrows.size() ;
299
300 // nchan
301 int refidx = -1 ;
302 vector<int> nchans ;
303 for ( int irow = 0 ; irow < nif ; irow++ ) {
304 nchans.push_back( spwrows[irow]->getNumChan() ) ;
305 if ( refidx == -1 && nchans[irow] != 1 && nchans[irow] != 4 )
306 refidx = irow ;
307 }
308 nchan = casa::Int( *max_element( nchans.begin(), nchans.end() ) ) ;
309
310 //logsink_->postLocally( LogMessage("refidx = "+String::toString(refidx),LogOrigin(className_,funcName,WHERE)) ) ;
311
312 // bandwidth
313 vector<double> bws ;
314 for ( int irow = 0 ; irow < nif ; irow++ ) {
315 if ( nchans[irow] != 4 ) { // exclude WVR data
316 bws.push_back( spwrows[irow]->getTotBandwidth().get() ) ;
317 }
318 }
319 bandwidth = casa::Double( *max_element( bws.begin(), bws.end() ) ) ;
320
321 // reffreq
322 reffreq = casa::Double( spwrows[refidx]->getRefFreq().get() ) ;
323
324 // freqref
325 if ( spwrows[refidx]->isMeasFreqRefExists() ) {
326 string mfr = CFrequencyReferenceCode::name( spwrows[refidx]->getMeasFreqRef() ) ;
327 if (mfr == "TOPO") {
328 freqref = "TOPOCENT";
329 } else if (mfr == "GEO") {
330 freqref = "GEOCENTR";
331 } else if (mfr == "BARY") {
332 freqref = "BARYCENT";
333 } else if (mfr == "GALACTO") {
334 freqref = "GALACTOC";
335 } else if (mfr == "LGROUP") {
336 freqref = "LOCALGRP";
337 } else if (mfr == "CMB") {
338 freqref = "CMBDIPOL";
339 } else if (mfr == "REST") {
340 freqref = "SOURCE";
341 }
342
343 }
344 else {
345 // frequency reference is TOPOCENT by default
346 freqref = "TOPOCENT" ;
347 }
348
349
350 PolarizationTable &ptab = asdm_->getPolarization() ;
351 vector<PolarizationRow *> prows = ptab.get() ;
352 vector<int> npols ;
353 refidx = -1 ;
354 for ( unsigned int irow = 0 ; irow < prows.size() ; irow++ ) {
355 npols.push_back( prows[irow]->getNumCorr() ) ;
356 if ( refidx == -1 && npols[irow] != 1 )
357 refidx = irow ;
358 }
359 if ( refidx == -1 )
360 refidx = 0 ;
361
362 // npol
363 npol = casa::Int( *max_element( npols.begin(), npols.end() ) ) ;
364
365 // poltype
366 vector<StokesParameter> corrType = prows[refidx]->getCorrType() ;
367 if ( corrType[0] == I ||
368 corrType[0] == Q ||
369 corrType[0] == U ||
370 corrType[0] == V )
371 poltype = "stokes" ;
372 else if ( corrType[0] == RR ||
373 corrType[0] == RL ||
374 corrType[0] == LR ||
375 corrType[0] == LL )
376 poltype = "circular" ;
377 else if ( corrType[0] == XX ||
378 corrType[0] == XY ||
379 corrType[0] == YX ||
380 corrType[0] == YY )
381 poltype = "linear" ;
382 else if ( corrType[0] == PLINEAR ||
383 corrType[0] == PANGLE ) {
384 poltype = "linpol" ;
385 }
386 else {
387 poltype = "notype" ;
388 }
389
390
391 FeedTable &ftab = asdm_->getFeed() ;
392 vector<FeedRow *> frows = ftab.get() ;
393 vector<int> nbeams ;
394 for ( unsigned int irow = 0 ; irow < frows.size() ; irow++ ) {
395 if ( frows[irow]->isFeedNumExists() )
396 nbeams.push_back( frows[irow]->getFeedNum() ) ;
397 else
398 nbeams.push_back( 1 ) ;
399 }
400
401 // nbeam
402 nbeam = casa::Int( *max_element( nbeams.begin(), nbeams.end() ) ) ;
403
404 // fluxunit
405 // tentatively set 'K'? or empty?
406 fluxunit = "K" ;
407
408 // equinox
409 // tentatively set 2000.0
410 equinox = 2000.0 ;
411
412 // epoch
413 // tentatively set "UTC"
414 epoch = "UTC" ;
415
416 // obstype
417 // at the moment take observingMode attribute in SBSummary table
418 SBSummaryRow *sbrow = ebrow->getSBSummaryUsingSBSummaryId() ;
419 vector<string> obsmode = sbrow->getObservingMode() ;
420 obstype = "" ;
421 for ( unsigned int imode = 0 ; imode < obsmode.size() ; imode++ ) {
422 obstype += casa::String(obsmode[imode]) ;
423 if ( imode != obsmode.size()-1 )
424 obstype += "#" ;
425 }
426}
427
428void ASDMReader::selectConfigDescription()
429{
430 String funcName = "selectConfigDescription" ;
431
432 vector<ConfigDescriptionRow *> cdrows = asdm_->getConfigDescription().get() ;
433 vector<Tag> cdidTags ;
434 for ( unsigned int irow = 0 ; irow < cdrows.size() ; irow++ ) {
435 //logsink_->postLocally( LogMessage("correlationMode["+String::toString(irow)+"] = "+String::toString(cdrows[irow]->getCorrelationMode()),LogOrigin(className_,funcName,WHERE)) ) ;
436 if ( cdrows[irow]->getCorrelationMode() != CROSS_ONLY ) {
437 cdidTags.push_back( cdrows[irow]->getConfigDescriptionId() ) ;
438 }
439 }
440
441 configDescIdList_.resize( cdidTags.size() ) ;
442 for ( unsigned int i = 0 ; i < cdidTags.size() ; i++ ) {
443 configDescIdList_[i] = casa::uInt( cdidTags[i].getTagValue() ) ;
444 }
445}
446
447void ASDMReader::selectFeed()
448{
449 feedIdList_.resize(0) ;
450 vector<FeedRow *> frows = asdm_->getFeed().get() ;
451 Tag atag( antennaId_, TagType::Antenna ) ;
452 for ( unsigned int irow = 0 ; irow < frows.size() ; irow++ ) {
453 casa::uInt feedId = (casa::uInt)(frows[irow]->getFeedId() ) ;
454 if ( casa::anyEQ( feedIdList_, feedId ) )
455 continue ;
456 if ( frows[irow]->getAntennaId() == atag ) {
457 unsigned int oldsize = feedIdList_.size() ;
458 feedIdList_.resize( oldsize+1, true ) ;
459 feedIdList_[oldsize] = feedId ;
460 }
461 }
462}
463
464casa::Vector<casa::uInt> ASDMReader::getFieldIdList()
465{
466 String funcName = "getFieldIdList" ;
467
468 vector<FieldRow *> frows = asdm_->getField().get() ;
469 fieldIdList_.resize( frows.size() ) ;
470 for ( unsigned int irow = 0 ; irow < frows.size() ; irow++ ) {
471 //logsink_->postLocally( LogMessage("fieldId["+String::toString(irow)+"]="+String(frows[irow]->getFieldId().toString()),LogOrigin(className_,funcName,WHERE)) ) ;
472 fieldIdList_[irow] = frows[irow]->getFieldId().getTagValue() ;
473 }
474
475 return fieldIdList_ ;
476}
477
478casa::uInt ASDMReader::getNumMainRow()
479{
480 casa::uInt nrow = casa::uInt( mainRow_.size() ) ;
481
482 return nrow ;
483}
484
485void ASDMReader::select()
486{
487 // selection by input CorrelationMode
488 EnumSet<CorrelationMode> esCorrs ;
489 sdmBin_->select( corrMode_ ) ;
490
491 // selection by TimeSampling
492 sdmBin_->select( timeSampling_ ) ;
493
494 // selection by SpectralResolutionType
495 sdmBin_->select( resolutionType_ ) ;
496
497 // selection by AtmPhaseCorrection and output CorrelationMode
498 EnumSet<AtmPhaseCorrection> esApcs ;
499 esApcs.set( apc_ ) ;
500 // always take only autocorrelation data
501 Enum<CorrelationMode> esCorr = AUTO_ONLY ;
502 sdmBin_->selectDataSubset( esCorr, esApcs ) ;
503
504 // select valid configDescriptionId
505 selectConfigDescription() ;
506
507 // select valid feedId
508 selectFeed() ;
509}
510
511casa::Bool ASDMReader::setMainRow( casa::uInt irow )
512{
513 casa::Bool status = true ;
514 row_ = irow ;
515
516 unsigned int cdid = mainRow_[row_]->getConfigDescriptionId().getTagValue() ;
517 if ( (int)count(configDescIdList_.begin(),configDescIdList_.end(),cdid) == 0 )
518 status = false ;
519 else {
520 status = (casa::Bool)(sdmBin_->acceptMainRow( mainRow_[row_] )) ;
521 }
522 return status ;
523}
524
525casa::Bool ASDMReader::setMainRow( casa::uInt configDescId, casa::uInt fieldId )
526{
527 clearMainRow() ;
528
529 Tag configDescTag( (unsigned int)configDescId, TagType::ConfigDescription ) ;
530 Tag fieldTag( (unsigned int)fieldId, TagType::Field ) ;
531 mainRow_ = casa::Vector<MainRow *>( *(asdm_->getMain().getByContext( configDescTag, fieldTag ) ) ) ;
532
533 return true ;
534}
535
536void ASDMReader::clearMainRow()
537{
538 mainRow_.resize(0) ;
539}
540
541void ASDMReader::setupIFNO()
542{
543 String funcName = "setupIFNO" ;
544
545 vector<SpectralWindowRow *> spwrows = asdm_->getSpectralWindow().get() ;
546 unsigned int nrow = spwrows.size() ;
547 ifno_.clear() ;
548 casa::uInt idx = 0 ;
549 casa::uInt wvridx = 0 ;
550 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
551 casa::uInt index ;
552 if ( isWVR( spwrows[irow] ) ) {
553 //logsink_->postLocally( LogMessage(spwrows[irow]->getSpectralWindowId().toString()+" is WVR",LogOrigin(className_,funcName,WHERE)) ) ;
554 index = wvridx ;
555 }
556 else {
557 index = ++idx ;
558 }
559 ifno_.insert( pair<Tag,casa::uInt>(spwrows[irow]->getSpectralWindowId(),index) ) ;
560 //logsink_->postLocally( LogMessage(spwrows[irow]->getSpectralWindowId().toString()+": IFNO="+String::toString(index),LogOrigin(className_,funcName,WHERE)) ) ;
561 }
562}
563
564bool ASDMReader::isWVR( SpectralWindowRow *row )
565{
566 BasebandName bbname = row->getBasebandName() ;
567 int nchan = row->getNumChan() ;
568 if ( bbname == NOBB && nchan == 4 )
569 return true ;
570 else
571 return false ;
572}
573
574// casa::Vector<casa::uInt> ASDMReader::getDataDescIdList( casa::uInt cdid )
575// {
576// Tag cdTag( (unsigned int)cdid, TagType::ConfigDescription ) ;
577// ConfigDescriptionRow *cdrow = asdm_->getConfigDescription().getRowByKey( cdTag ) ;
578// vector<Tag> ddTags = cdrow->getDataDescriptionId() ;
579// casa::Vector<casa::uInt> ddidList( ddTags.size() ) ;
580// for ( unsigned int idd = 0 ; idd < ddTags.size() ; idd++ ) {
581// ddidList[idd] = ddTags[idd].getTagValue() ;
582// }
583// return ddidList ;
584// }
585
586// casa::Vector<casa::uInt> ASDMReader::getSwitchCycleIdList( casa::uInt cdid )
587// {
588// Tag cdTag( (unsigned int)cdid, TagType::ConfigDescription ) ;
589// ConfigDescriptionRow *cdrow = asdm_->getConfigDescription().getRowByKey( cdTag ) ;
590// vector<Tag> scTags = cdrow->getSwitchCycleId() ;
591// casa::Vector<casa::uInt> scidList( scTags.size() ) ;
592// for ( unsigned int idd = 0 ; idd < scTags.size() ; idd++ ) {
593// scidList[idd] = scTags[idd].getTagValue() ;
594// }
595// return scidList ;
596// }
597
598// casa::Vector<casa::uInt> ASDMReader::getFeedIdList( casa::uInt cdid )
599// {
600// String funcName = "getFeedIdList" ;
601//
602// Tag cdTag( (unsigned int)cdid, TagType::ConfigDescription ) ;
603// ConfigDescriptionRow *cdrow = asdm_->getConfigDescription().getRowByKey( cdTag ) ;
604// casa::Vector<casa::uInt> feedIdList ;
605// vector<int> feedIds = cdrow->getFeedId() ;
606// for ( unsigned int ife = 0 ; ife < feedIds.size() ; ife++ ) {
607// logsink_->postLocally( LogMessage("feedIds["+String::toString(ife)+"]="+String::toString(feedIds[ife]),LogOrigin(className_,funcName,WHERE)) ) ;
608// if ( casa::anyEQ( feedIdList, casa::uInt( feedIds[ife] ) ) )
609// continue ;
610// if ( casa::anyEQ( feedIdList_, casa::uInt( feedIds[ife] ) ) ) {
611// casa::uInt oldsize = feedIdList.size() ;
612// feedIdList.resize( oldsize+1, true ) ;
613// feedIdList[oldsize] = casa::uInt( feedIds[ife] ) ;
614// }
615// }
616// logsink_->postLocally( LogMessage("feedIdList.size() = "+String::toString(feedIdList.size()),LogOrigin(className_,funcName,WHERE)) ) ;
617// return feedIdList ;
618// }
619
620casa::Bool ASDMReader::setData()
621{
622 String funcName = "setData" ;
623
624 //logsink_->postLocally( LogMessage("try to retrieve binary data",LogOrigin(className_,funcName,WHERE)) ) ;
625
626// EnumSet<AtmPhaseCorrection> esApcs ;
627// esApcs.set( apc_ ) ;
628// // always take only autocorrelation data
629// Enum<CorrelationMode> esCorr = AUTO_ONLY ;
630// vmsData_ = sdmBin_->getDataCols( esCorr, esApcs ) ;
631
632 // 2011/07/06 TN
633 // Workaround to avoid unwanted message from SDMBinData::getDataCols()
634 ostringstream oss ;
635 streambuf *buforg = cout.rdbuf(oss.rdbuf()) ;
636 vmsData_ = sdmBin_->getDataCols() ;
637 cout.rdbuf(buforg) ;
638
639// logsink_->postLocally( LogMessage("oss.str() = "+oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
640// cout << "This is test: oss.str()=" << oss.str() << endl ;
641
642
643 //logsink_->postLocally( LogMessage("processorId = "+String::toString(vmsData_->processorId),LogOrigin(className_,funcName,WHERE)) ) ;
644 //logsink_->postLocally( LogMessage("v_time.size() = "+String::toString(vmsData_->v_time.size()),LogOrigin(className_,funcName,WHERE)) ) ;
645 //logsink_->postLocally( LogMessage(" v_time[0] = "+String::toString(vmsData_->v_time[0]),LogOrigin(className_,funcName,WHERE)) ) ;
646 //logsink_->postLocally( LogMessage("v_interval.size() = "+String::toString(vmsData_->v_interval.size()),LogOrigin(className_,funcName,WHERE)) ) ;
647 //logsink_->postLocally( LogMessage(" v_interval[0] = "+String::toString(vmsData_->v_interval[0]),LogOrigin(className_,funcName,WHERE)) ) ;
648 //logsink_->postLocally( LogMessage("v_atmPhaseCorrection.size() = "+String::toString(vmsData_->v_atmPhaseCorrection.size()),LogOrigin(className_,funcName,WHERE)) ) ;
649 //logsink_->postLocally( LogMessage("binNum = "+String::toString(vmsData_->binNum),LogOrigin(className_,funcName,WHERE)) ) ;
650 //logsink_->postLocally( LogMessage("v_projectPath.size() = "+String::toString(vmsData_->v_projectPath.size()),LogOrigin(className_,funcName,WHERE)) ) ;
651 //logsink_->postLocally( LogMessage("v_antennaId1.size() = "+String::toString(vmsData_->v_antennaId1.size()),LogOrigin(className_,funcName,WHERE)) ) ;
652 //logsink_->postLocally( LogMessage("v_antennaId2.size() = "+String::toString(vmsData_->v_antennaId2.size()),LogOrigin(className_,funcName,WHERE)) ) ;
653 //logsink_->postLocally( LogMessage("v_feedId1.size() = "+String::toString(vmsData_->v_feedId1.size()),LogOrigin(className_,funcName,WHERE)) ) ;
654 //logsink_->postLocally( LogMessage("v_feedId2.size() = "+String::toString(vmsData_->v_feedId2.size()),LogOrigin(className_,funcName,WHERE)) ) ;
655 //logsink_->postLocally( LogMessage("v_dataDescId.size() = "+String::toString(vmsData_->v_dataDescId.size()),LogOrigin(className_,funcName,WHERE)) ) ;
656 //logsink_->postLocally( LogMessage("v_timeCentroid.size() = "+String::toString(vmsData_->v_timeCentroid.size()),LogOrigin(className_,funcName,WHERE)) ) ;
657 //logsink_->postLocally( LogMessage("v_exposure.size() = "+String::toString(vmsData_->v_exposure.size()),LogOrigin(className_,funcName,WHERE)) ) ;
658 //logsink_->postLocally( LogMessage("v_numData.size() = "+String::toString(vmsData_->v_numData.size()),LogOrigin(className_,funcName,WHERE)) ) ;
659 //logsink_->postLocally( LogMessage("vv_dataShape.size() = "+String::toString(vmsData_->vv_dataShape.size()),LogOrigin(className_,funcName,WHERE)) ) ;
660 //logsink_->postLocally( LogMessage("v_m_data.size() = "+String::toString(vmsData_->v_m_data.size()),LogOrigin(className_,funcName,WHERE)) ) ;
661 //logsink_->postLocally( LogMessage("v_phaseDir.size() = "+String::toString(vmsData_->v_phaseDir.size()),LogOrigin(className_,funcName,WHERE)) ) ;
662 //logsink_->postLocally( LogMessage("v_stateId.size() = "+String::toString(vmsData_->v_stateId.size()),LogOrigin(className_,funcName,WHERE)) ) ;
663 //logsink_->postLocally( LogMessage("v_msState.size() = "+String::toString(vmsData_->v_msState.size()),LogOrigin(className_,funcName,WHERE)) ) ;
664 //logsink_->postLocally( LogMessage("v_flag.size() = "+String::toString(vmsData_->v_flag.size()),LogOrigin(className_,funcName,WHERE)) ) ;
665
666 dataIdList_.clear() ;
667 unsigned int numTotalData = vmsData_->v_m_data.size() ;
668 for ( unsigned int idata = 0 ; idata < numTotalData ; idata++ ) {
669 if ( vmsData_->v_antennaId1[idata] == (int)antennaId_
670 && vmsData_->v_antennaId2[idata] == (int)antennaId_ )
671 dataIdList_.push_back( idata ) ;
672 }
673 numData_ = dataIdList_.size() ;
674 //logsink_->postLocally( LogMessage("numData_ = "+String::toString(numData_),LogOrigin(className_,funcName,WHERE)) ) ;
675 //logsink_->postLocally( LogMessage("dataSize = "+String::toString(mainRow_[row_]->getDataSize()),LogOrigin(className_,funcName,WHERE)) ) ;
676
677 return true ;
678}
679
680casa::uInt ASDMReader::getIFNo( unsigned int idx )
681{
682 Tag ddTag( vmsData_->v_dataDescId[dataIdList_[idx]], TagType::DataDescription ) ;
683 DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey( ddTag ) ;
684 Tag spwid = ddrow->getSpectralWindowId() ;
685 map<Tag,casa::uInt>::iterator iter = ifno_.find( spwid ) ;
686 if ( iter != ifno_.end() )
687 return iter->second ;
688 else {
689 return 0 ;
690 }
691}
692
693int ASDMReader::getNumPol( unsigned int idx )
694{
695 Tag ddTag( vmsData_->v_dataDescId[dataIdList_[idx]], TagType::DataDescription ) ;
696 DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey( ddTag ) ;
697 PolarizationRow *polrow = ddrow->getPolarizationUsingPolOrHoloId() ;
698 return polrow->getNumCorr() ;
699}
700
701void ASDMReader::getFrequency( unsigned int idx,
702 double &refpix,
703 double &refval,
704 double &incr )
705{
706 String funcName = "getFrequency" ;
707
708 Tag ddTag( vmsData_->v_dataDescId[dataIdList_[idx]], TagType::DataDescription ) ;
709 DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey( ddTag ) ;
710 //Tag spwid = ddrow->getSpectralWindowId() ;
711 SpectralWindowRow *spwrow = ddrow->getSpectralWindowUsingSpectralWindowId() ;
712 int nchan = spwrow->getNumChan() ;
713 if ( nchan == 1 ) {
714 //logsink_->postLocally( LogMessage("channel averaged data",LogOrigin(className_,funcName,WHERE)) ) ;
715 refpix = 0.0 ;
716 incr = spwrow->getTotBandwidth().get() ;
717 if ( spwrow->isChanFreqStartExists() ) {
718 refval = spwrow->getChanFreqStart().get() ;
719 }
720 else if ( spwrow->isChanFreqArrayExists() ) {
721 refval = spwrow->getChanFreqArray()[0].get() ;
722 }
723 else {
724 throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
725 }
726 }
727 else if ( nchan % 2 ) {
728 // odd
729 //logsink_->postLocally( LogMessage("odd case",LogOrigin(className_,funcName,WHERE)) ) ;
730 refpix = 0.5 * ( (double)nchan - 1.0 ) ;
731 int ic = ( nchan - 1 ) / 2 ;
732 if ( spwrow->isChanWidthExists() ) {
733 incr = spwrow->getChanWidth().get() ;
734 }
735 else if ( spwrow->isChanWidthArrayExists() ) {
736 incr = spwrow->getChanWidthArray()[0].get() ;
737 }
738 else {
739 throw (casa::AipsError( "Either chanWidthArray or chanWidth must exist." )) ;
740 }
741 if ( spwrow->isChanFreqStepExists() ) {
742 if ( spwrow->getChanFreqStep().get() < 0.0 )
743 incr *= -1.0 ;
744 }
745 else if ( spwrow->isChanFreqArrayExists() ) {
746 vector<Frequency> chanFreqArr = spwrow->getChanFreqArray() ;
747 if ( chanFreqArr[0].get() > chanFreqArr[1].get() )
748 incr *= -1.0 ;
749 }
750 else {
751 throw (casa::AipsError( "Either chanFreqArray or chanFreqStep must exist." )) ;
752 }
753 if ( spwrow->isChanFreqStartExists() ) {
754 refval = spwrow->getChanFreqStart().get() + refpix * incr ;
755 }
756 else if ( spwrow->isChanFreqArrayExists() ) {
757 refval = spwrow->getChanFreqArray()[ic].get() ;
758 }
759 else {
760 throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
761 }
762 }
763 else {
764 // even
765 //logsink_->postLocally( LogMessage("even case",LogOrigin(className_,funcName,WHERE)) ) ;
766 refpix = 0.5 * ( (double)nchan - 1.0 ) ;
767 int ic = nchan / 2 ;
768 if ( spwrow->isChanWidthExists() ) {
769 incr = spwrow->getChanWidth().get() ;
770 }
771 else if ( spwrow->isChanWidthArrayExists() ) {
772 incr = spwrow->getChanWidthArray()[0].get() ;
773 }
774 else {
775 throw (casa::AipsError( "Either chanWidthArray or chanWidth must exist." )) ;
776 }
777 if ( spwrow->isChanFreqStepExists() ) {
778 if ( spwrow->getChanFreqStep().get() < 0.0 )
779 incr *= -1.0 ;
780 }
781 else if ( spwrow->isChanFreqArrayExists() ) {
782 vector<Frequency> chanFreqArr = spwrow->getChanFreqArray() ;
783 if ( chanFreqArr[0].get() > chanFreqArr[1].get() )
784 incr *= -1.0 ;
785 }
786 else {
787 throw (casa::AipsError( "Either chanFreqArray or chanFreqStep must exist." )) ;
788 }
789 if ( spwrow->isChanFreqStartExists() ) {
790 refval = spwrow->getChanFreqStart().get() + refpix * incr ;
791 }
792 else if ( spwrow->isChanFreqArrayExists() ) {
793 vector<Frequency> freqs = spwrow->getChanFreqArray() ;
794 refval = 0.5 * ( freqs[ic-1].get() + freqs[ic].get() ) ;
795 }
796 else {
797 throw (casa::AipsError( "Either chanFreqArray or chanFreqStart must exist." )) ;
798 }
799 }
800}
801
802vector<double> ASDMReader::getRestFrequency( unsigned int idx )
803{
804 vector<double> rf( 0 ) ;
805 unsigned int index = dataIdList_[idx] ;
806 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
807 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
808 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
809 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
810 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
811 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
812 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
813 if ( frow->isSourceIdExists() ) {
814 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
815 int sid = frow->getSourceId() ;
816 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
817 if ( srow->isRestFrequencyExists() ) {
818 //logsink_->postLocally( LogMessage("restFrequency exists",LogOrigin(className_,funcName,WHERE)) ) ;
819 vector<Frequency> restfreq = srow->getRestFrequency() ;
820 rf.resize( restfreq.size() ) ;
821 for ( unsigned int i = 0 ; i < restfreq.size() ; i++ )
822 rf[i] = restfreq[i].get() ;
823 }
824 }
825 return rf ;
826}
827
828double ASDMReader::getTime( unsigned int idx )
829{
830 double tsec = vmsData_->v_time[dataIdList_[idx]] ;
831 return tsec * s2d ;
832}
833
834double ASDMReader::getInterval( unsigned int idx )
835{
836 return vmsData_->v_interval[dataIdList_[idx]] ;
837}
838
839string ASDMReader::getSourceName( unsigned int idx )
840{
841 unsigned int index = dataIdList_[idx] ;
842 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
843 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
844 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
845 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
846 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
847 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
848 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
849 string srcname ;
850 if ( frow->isSourceIdExists() ) {
851 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
852 int sid = frow->getSourceId() ;
853 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
854 srcname = srow->getSourceName() ;
855 }
856 else {
857 srcname = frow->getFieldName() ;
858 }
859 return srcname ;
860}
861
862string ASDMReader::getFieldName( unsigned int idx )
863{
864 int fid = vmsData_->v_fieldId[dataIdList_[idx]] ;
865 Tag ftag( fid, TagType::Field ) ;
866 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
867 ostringstream oss ;
868 oss << frow->getFieldName() << "__" << fid ;
869 return oss.str() ;
870}
871
872int ASDMReader::getSrcType( unsigned int scan,
873 unsigned int subscan )
874{
875 int srctype = SrcType::NOTYPE ;
876 Tag ebtag = mainRow_[row_]->getExecBlockId() ;
877 ScanRow *scanrow = asdm_->getScan().getRowByKey( ebtag, (int)scan ) ;
878 ScanIntent scanIntent = scanrow->getScanIntent()[0] ;
879 SubscanRow *subrow = asdm_->getSubscan().getRowByKey( ebtag, (int)scan, (int)subscan ) ;
880 SubscanIntent subIntent = subrow->getSubscanIntent() ;
881 SwitchingMode swmode = subrow->getSubscanMode() ;
882 if ( scanIntent == OBSERVE_TARGET ) {
883 // on sky scan
884 if ( swmode == NO_SWITCHING || swmode == POSITION_SWITCHING ) {
885 // position switching
886 // tentatively set NO_SWITCHING = POSITION_SWITCHING
887 if ( subIntent == ON_SOURCE ) {
888 srctype = SrcType::PSON ;
889 }
890 else if ( subIntent == OFF_SOURCE ) {
891 srctype = SrcType::PSOFF ;
892 }
893 }
894 else if ( swmode == FREQUENCY_SWITCHING ) {
895 // frequency switching
896 if ( subIntent == ON_SOURCE ) {
897 srctype = SrcType::FSON ;
898 }
899 else if ( subIntent == OFF_SOURCE ) {
900 srctype = SrcType::FSOFF ;
901 }
902 }
903 else if ( swmode == NUTATOR_SWITCHING ) {
904 // nutator switching
905 if ( subIntent == ON_SOURCE ) {
906 srctype = SrcType::PSON ;
907 }
908 else if ( subIntent == OFF_SOURCE ) {
909 srctype = SrcType::PSOFF ;
910 }
911 }
912 else {
913 // other switching mode
914 if ( subIntent == ON_SOURCE ) {
915 srctype = SrcType::SIG ;
916 }
917 else if ( subIntent == OFF_SOURCE ) {
918 srctype = SrcType::REF ;
919 }
920 }
921 }
922 else {
923 // off sky (e.g. calibrator device) scan
924 if ( subIntent == ON_SOURCE ) {
925 srctype = SrcType::SIG ;
926 }
927 else if ( subIntent == OFF_SOURCE ) {
928 srctype = SrcType::REF ;
929 }
930 else if ( subIntent == HOT ) {
931 srctype = SrcType::HOT ;
932 }
933 else if ( subIntent == AMBIENT ) {
934 srctype = SrcType::SKY ;
935 }
936 else {
937 srctype = SrcType::CAL ;
938 }
939 }
940
941 return srctype ;
942}
943
944unsigned int ASDMReader::getSubscanNo( unsigned int idx )
945{
946 //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)) ) ;
947 return vmsData_->v_msState[dataIdList_[idx]].subscanNum ;
948}
949
950vector<double> ASDMReader::getSourceDirection( unsigned int idx )
951{
952 vector<double> dir( 2, 0.0 ) ;
953 unsigned int index = dataIdList_[idx] ;
954 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
955 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
956 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
957 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
958 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
959 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
960 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
961 string srcname ;
962 if ( frow->isSourceIdExists() ) {
963 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
964 int sid = frow->getSourceId() ;
965 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
966 vector<Angle> srcdir = srow->getDirection() ;
967 if ( srow->isDirectionCodeExists() ) {
968 DirectionReferenceCode dircode = srow->getDirectionCode() ;
969 if ( dircode != J2000 ) {
970 // if not J2000, need direction conversion
971 }
972 }
973 dir[0] = srcdir[0].get() ;
974 dir[1] = srcdir[1].get() ;
975 }
976 return dir ;
977}
978
979vector<double> ASDMReader::getSourceProperMotion( unsigned int idx )
980{
981 vector<double> pm( 2, 0.0 ) ;
982 unsigned int index = dataIdList_[idx] ;
983 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
984 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
985 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
986 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
987 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
988 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
989 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
990 string srcname ;
991 if ( frow->isSourceIdExists() ) {
992 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
993 int sid = frow->getSourceId() ;
994 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
995 vector<AngularRate> srcpm = srow->getProperMotion() ;
996 pm[0] = srcpm[0].get() ;
997 pm[1] = srcpm[1].get() ;
998 }
999 return pm ;
1000}
1001
1002double ASDMReader::getSysVel( unsigned int idx )
1003{
1004 double sysvel = 0.0 ;
1005 unsigned int index = dataIdList_[idx] ;
1006 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1007 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1008 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1009 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1010 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1011 Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
1012 FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
1013 string srcname ;
1014 if ( frow->isSourceIdExists() ) {
1015 //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
1016 int sid = frow->getSourceId() ;
1017 SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
1018 if ( srow->isSysVelExists() ) {
1019 vector<Speed> sysvelV = srow->getSysVel() ;
1020 if ( sysvelV.size() > 0 )
1021 sysvel = sysvelV[0].get() ;
1022 }
1023 }
1024 return sysvel ;
1025}
1026
1027unsigned int ASDMReader::getFlagRow( unsigned int idx )
1028{
1029 return vmsData_->v_flag[dataIdList_[idx]] ;
1030}
1031
1032vector<unsigned int> ASDMReader::getDataShape( unsigned int idx )
1033{
1034 return vmsData_->vv_dataShape[dataIdList_[idx]] ;
1035}
1036
1037float * ASDMReader::getSpectrum( unsigned int idx )
1038{
1039 map<AtmPhaseCorrection, float*> data = vmsData_->v_m_data[dataIdList_[idx]] ;
1040 //map<AtmPhaseCorrection, float*>::iterator iter = data.find(AP_UNCORRECTED) ;
1041 map<AtmPhaseCorrection, float*>::iterator iter = data.find(apc_) ;
1042 float *autoCorr = iter->second ;
1043 return autoCorr ;
1044}
1045
1046// bool * ASDMReader::getFlagChannel( unsigned int idx )
1047// {
1048// return 0 ;
1049// }
1050
1051vector< vector<float> > ASDMReader::getTsys( unsigned int idx )
1052{
1053 unsigned int index = dataIdList_[idx] ;
1054 Tag anttag( antennaId_, TagType::Antenna ) ;
1055 int feedid = vmsData_->v_feedId1[index] ;
1056 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1057 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1058 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1059 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1060 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1061 //int nchan = asdm_->getSpectralWindow().getRowByKey(spwtag)->getNumChan() ;
1062 vector< vector<float> > defaultTsys( 1, vector<float>( 1, 1.0 ) ) ;
1063 SysCalTable &sctab = asdm_->getSysCal() ;
1064 if ( sctab.size() == 0 )
1065 return defaultTsys ;
1066 SysCalRow *scrow = sctab.getRowByKey( anttag, spwtag, tint, feedid ) ;
1067 if ( scrow->isTsysSpectrumExists() ) {
1068 vector< vector<Temperature> > tsysSpec = scrow->getTsysSpectrum() ;
1069 unsigned int numReceptor = tsysSpec.size() ;
1070 unsigned int numChan = tsysSpec[0].size() ;
1071 vector< vector<float> > tsys( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1072 for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1073 for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1074 tsys[ir][ic] = tsysSpec[ir][ic].get() ;
1075 }
1076 }
1077 return tsys ;
1078 }
1079// else if ( scrow->isTsysExists() ) {
1080// vector<Temperature> tsysScalar = scrow->getTsys() ;
1081// unsigned int numReceptor = tsysScalar.size() ;
1082// vector< vector<float> > tsys( numReceptor ) ;
1083// for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1084// tsys[ir] = vector<float>( 1, tsysScalar[ir].get() ) ;
1085// return tsys ;
1086// }
1087 else {
1088 return defaultTsys ;
1089 }
1090}
1091
1092vector< vector<float> > ASDMReader::getTcal( unsigned int idx )
1093{
1094 unsigned int index = dataIdList_[idx] ;
1095 Tag anttag( antennaId_, TagType::Antenna ) ;
1096 int feedid = vmsData_->v_feedId1[index] ;
1097 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1098 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1099 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1100 Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1101 Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
1102 //int nchan = asdm_->getSpectralWindow().getRowByKey(spwtag)->getNumChan() ;
1103 vector< vector<float> > defaultTcal( 1, vector<float>( 1, 1.0 ) ) ;
1104 SysCalTable &sctab = asdm_->getSysCal() ;
1105 if ( sctab.size() == 0 )
1106 return defaultTcal ;
1107 SysCalRow *scrow = sctab.getRowByKey( anttag, spwtag, tint, feedid ) ;
1108 if ( scrow->isTcalSpectrumExists() ) {
1109 vector< vector<Temperature> > tcalSpec = scrow->getTcalSpectrum() ;
1110 unsigned int numReceptor = tcalSpec.size() ;
1111 unsigned int numChan = tcalSpec[0].size() ;
1112 vector< vector<float> > tcal( numReceptor, vector<float>( numChan, 1.0 ) ) ;
1113 for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ ) {
1114 for ( unsigned int ic = 0 ; ic < numChan ; ic++ ) {
1115 tcal[ir][ic] = tcalSpec[ir][ic].get() ;
1116 }
1117 }
1118 return tcal ;
1119 }
1120// else if ( scrow->isTcalExists() ) {
1121// vector<Temperature> tcalScalar = scrow->getTcal() ;
1122// unsigned int numReceptor = tcalScalar.size() ;
1123// vector< vector<float> > tcal( numReceptor, vector<float>( 1, 1.0 ) ) ;
1124// for ( unsigned int ir = 0 ; ir < numReceptor ; ir++ )
1125// tcal[ir][0] = tcalScalar[ir][0].get() ;
1126// return tcal ;
1127// }
1128 else {
1129 return defaultTcal ;
1130 }
1131}
1132
1133vector<float> ASDMReader::getOpacity( unsigned int idx )
1134{
1135 vector<float> tau(0) ;
1136 CalAtmosphereTable &atmtab = asdm_->getCalAtmosphere() ;
1137 unsigned int nrow = atmtab.size() ;
1138 if ( nrow > 0 ) {
1139 unsigned int index = dataIdList_[idx] ;
1140 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1141 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1142 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1143 //int feedid = vmsData_->v_feedId1[index] ;
1144 //Tag atag( antennaId_, TagType::Antenna ) ;
1145 //Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
1146 //DataDescriptionRow *ddrow = asdm_->getDataDescription().getRowByKey(ddtag) ;
1147 //Tag spwtag = ddrow->getSpectralWindowId() ;
1148 //SpectralWindowRow *spwrow = ddrow->getSpectralWindowUsingSpectralWindowId() ;
1149 //BasebandName bbname = spwrow->getBasebandName() ;
1150 //FeedRow *frow = asdm_->getFeed().getRowByKey( atag, spwtag, tint, feedid ) ;
1151 //int nfeed = frow->getNumReceptor() ;
1152 //vector<ReceiverRow *> rrows = frow->getReceivers() ;
1153 vector<CalAtmosphereRow *> atmrows = atmtab.get() ;
1154 //ReceiverBand rb = rrows[0]->getFrequencyBand() ;
1155 int row0 = -1 ;
1156 double eps = tint.getStart().getMJD() ;
1157 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1158 CalAtmosphereRow *atmrow = atmrows[irow] ;
1159 if ( casa::String(atmrow->getAntennaName()) != antennaName_
1160 //|| atmrow->getReceiverBand() != rb
1161 //|| atmrow->getBasebandName() != bbname
1162 || atmrow->getCalDataUsingCalDataId()->getCalType() != CAL_ATMOSPHERE )
1163 continue ;
1164 else {
1165 double dt = tint.getStart().getMJD() - atmrow->getEndValidTime().getMJD() ;
1166 if ( dt >= 0 && dt < eps ) {
1167 eps = dt ;
1168 row0 = (int)irow ;
1169 }
1170 }
1171 }
1172 if ( row0 != -1 ) {
1173 CalAtmosphereRow *atmrow = atmrows[row0] ;
1174 tau = atmrow->getTau() ;
1175 }
1176 }
1177 return tau ;
1178}
1179
1180void ASDMReader::getWeatherInfo( unsigned int idx,
1181 float &temperature,
1182 float &pressure,
1183 float &humidity,
1184 float &windspeed,
1185 float &windaz )
1186{
1187 temperature = 0.0 ;
1188 pressure = 0.0 ;
1189 humidity = 0.0 ;
1190 windspeed = 0.0 ;
1191 windaz = 0.0 ;
1192
1193 //logsink_->postLocally( LogMessage("weatherStationId_ = "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1194
1195 WeatherTable &wtab = asdm_->getWeather() ;
1196 if ( wtab.size() == 0 || weatherStationId_ == -1 )
1197 return ;
1198
1199 unsigned int index = dataIdList_[idx] ;
1200 //Tag anttag( antennaId_, TagType::Antenna ) ;
1201 //Tag sttag = (asdm_->getAntenna().getRowByKey( anttag ))->getStationId() ;
1202 Tag sttag( (unsigned int)weatherStationId_, TagType::Station ) ;
1203 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1204 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1205 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1206 //WeatherRow *wrow = wtab.getRowByKey( sttag, tint ) ;
1207 vector<WeatherRow *> *wrows = wtab.getByContext( sttag ) ;
1208 WeatherRow *wrow = (*wrows)[0] ;
1209 unsigned int nrow = wrows->size() ;
1210 //logsink_->postLocally( LogMessage("There are "+String::toString(nrow)+" rows for given context: stationId "+String::toString(weatherStationId_),LogOrigin(className_,funcName,WHERE)) ) ;
1211 ArrayTime startTime = tint.getStart() ;
1212 if ( startTime < (*wrows)[0]->getTimeInterval().getStart() ) {
1213 temperature = (*wrows)[0]->getTemperature().get() ;
1214 pressure = (*wrows)[0]->getPressure().get() ;
1215 humidity = (*wrows)[0]->getRelHumidity().get() ;
1216 windspeed = (*wrows)[0]->getWindSpeed().get() ;
1217 windaz = (*wrows)[0]->getWindDirection().get() ;
1218 }
1219 else if ( startTime > (*wrows)[nrow-1]->getTimeInterval().getStart() ) {
1220 temperature = (*wrows)[nrow-1]->getTemperature().get() ;
1221 pressure = (*wrows)[nrow-1]->getPressure().get() ;
1222 humidity = (*wrows)[nrow-1]->getRelHumidity().get() ;
1223 windspeed = (*wrows)[nrow-1]->getWindSpeed().get() ;
1224 windaz = (*wrows)[nrow-1]->getWindDirection().get() ;
1225 }
1226 else {
1227 for ( unsigned int irow = 1 ; irow < wrows->size() ; irow++ ) {
1228 wrow = (*wrows)[irow-1] ;
1229 if ( startTime < (*wrows)[irow]->getTimeInterval().getStart() ) {
1230 temperature = wrow->getTemperature().get() ;
1231 pressure = wrow->getPressure().get() ;
1232 humidity = wrow->getRelHumidity().get() ;
1233 windspeed = wrow->getWindSpeed().get() ;
1234 windaz = wrow->getWindDirection().get() ;
1235 break ;
1236 }
1237 }
1238 }
1239
1240 return ;
1241}
1242
1243void ASDMReader::processStation()
1244{
1245 antennaPad_.resize(0) ;
1246 weatherStation_.resize(0) ;
1247 StationTable &stab = asdm_->getStation() ;
1248 vector<StationRow *> srows = stab.get() ;
1249 for ( unsigned int irow = 0 ; irow < srows.size() ; irow++ ) {
1250 StationType stype = srows[irow]->getType() ;
1251 Tag stag = srows[irow]->getStationId() ;
1252 if ( stype == ANTENNA_PAD )
1253 antennaPad_.push_back( stag ) ;
1254 else if ( stype == WEATHER_STATION )
1255 weatherStation_.push_back( stag ) ;
1256 }
1257
1258 weatherStationId_ = getClosestWeatherStation() ;
1259}
1260
1261int ASDMReader::getClosestWeatherStation()
1262{
1263 if ( weatherStation_.size() == 0 )
1264 return -1 ;
1265
1266 Tag atag( antennaId_, TagType::Antenna ) ;
1267 Tag stag = (asdm_->getAntenna().getRowByKey( atag ))->getStationId() ;
1268 vector<double> apos( 3 ) ;
1269 StationTable &stab = asdm_->getStation() ;
1270 StationRow *srow = stab.getRowByKey( stag ) ;
1271 vector<Length> pos = srow->getPosition() ;
1272 apos[0] = pos[0].get() ;
1273 apos[1] = pos[1].get() ;
1274 apos[2] = pos[2].get() ;
1275
1276 double eps = 1.0e20 ;
1277 int retval = -1 ;
1278 for ( unsigned int ir = 0 ; ir < weatherStation_.size() ; ir++ ) {
1279 srow = stab.getRowByKey( weatherStation_[ir] ) ;
1280 vector<Length> wpos = srow->getPosition() ;
1281 double dist = (apos[0]-wpos[0].get())*(apos[0]-wpos[0].get())
1282 + (apos[1]-wpos[1].get())*(apos[1]-wpos[1].get())
1283 + (apos[2]-wpos[2].get())*(apos[2]-wpos[2].get()) ;
1284 if ( dist < eps ) {
1285 retval = (int)(weatherStation_[ir].getTagValue()) ;
1286 }
1287 }
1288
1289 return retval ;
1290}
1291
1292void ASDMReader::getPointingInfo( unsigned int idx,
1293 vector<double> &dir,
1294 double &az,
1295 double &el,
1296 vector<double> &srate )
1297{
1298 String funcName = "getPointingInfo" ;
1299
1300 dir.resize(0) ;
1301 az = -1.0 ;
1302 el = -1.0 ;
1303 srate.resize(0) ;
1304
1305 Tag atag( antennaId_, TagType::Antenna ) ;
1306 unsigned int index = dataIdList_[idx] ;
1307 vector<PointingRow *> *prows = asdm_->getPointing().getByContext( atag ) ;
1308 PointingRow *prow ;
1309 PointingRow *qrow ;
1310 //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
1311 double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
1312 ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
1313
1314 unsigned int nrow = prows->size() ;
1315 //logsink_->postLocally( LogMessage("There are " << nrow << " rows for given context: antennaId "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
1316
1317// for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1318// prow = (*prows)[irow] ;
1319// ArrayTimeInterval ati = prow->getTimeInterval() ;
1320// ArrayTime pst = ati.getStart() ;
1321// ArrayTime pet( ati.getStartInMJD()+ati.getDurationInDays() ) ;
1322// logsink_->postLocally( LogMessage("start: "+pst.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1323// logsink_->postLocally( LogMessage("end: "+pet.toFITS(),LogOrigin(className_,funcName,WHERE)) ) ;
1324// }
1325
1326 if ( nrow == 0 )
1327 return ;
1328
1329 srate.resize(2) ;
1330 srate[0] = 0.0 ;
1331 srate[1] = 0.0 ;
1332 az = 0.0 ;
1333 el = 0.0 ;
1334 double tcen = 0.0 ;
1335
1336 //
1337 // shape of pointingDirection is (numSample,2) if usePolynomial = False, while it is
1338 // (numTerm,2) if usePolynomial = True.
1339 //
1340 // In the former case, typical sampled time interval is 48msec, which is very small
1341 // compared with typical integration time (~a few sec).
1342 // Scan rate for this case is always [0,0] (or get slope?).
1343 //
1344 // In the later case, scan rate is (pointingDirection[1][0],pointingDirection[1][1])
1345 //
1346 // PointingDirection seems to store AZEL
1347 //
1348 ArrayTimeInterval pTime0 = (*prows)[0]->getTimeInterval() ;
1349 ArrayTimeInterval pTime1 = (*prows)[nrow-1]->getTimeInterval() ;
1350 if ( tint.getStartInMJD()+tint.getDurationInDays() < pTime0.getStartInMJD() ) {
1351 logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint < ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1352 prow = (*prows)[0] ;
1353 //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1354 vector< vector<Angle> > dirA = prow->getTarget() ;
1355 vector< vector<Angle> > offA = prow->getOffset() ;
1356 //az = dirA[0][0].get() ;
1357 //el = dirA[0][1].get() ;
1358 az = dirA[0][0].get() + offA[0][0].get() ;
1359 el = dirA[0][1].get() + offA[0][1].get() ;
1360 if ( prow->getUsePolynomials() && prow->getNumTerm() > 1 ) {
1361 //srate[0] = dirA[1][0].get() ;
1362 //srate[1] = dirA[1][1].get() ;
1363 srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1364 srate[1] = dirA[1][1].get() + offA[1][1].get() ;
1365 }
1366 }
1367 else if ( tint.getStartInMJD() > pTime1.getStartInMJD()+pTime1.getDurationInDays() ) {
1368 logsink_->postLocally( LogMessage( "ArrayTimeInterval out of bounds: no data for given position (tint > ptime)", LogOrigin(className_,funcName,WHERE), LogMessage::WARN ) ) ;
1369 prow = (*prows)[nrow-1] ;
1370 int numSample = prow->getNumSample() ;
1371 //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1372 vector< vector<Angle> > dirA = prow->getTarget() ;
1373 vector< vector<Angle> > offA = prow->getOffset() ;
1374 if ( prow->getUsePolynomials() ) {
1375 //az = dirA[0][0].get() ;
1376 //el = dirA[0][1].get() ;
1377 az = dirA[0][0].get() + offA[0][0].get() ;
1378 el = dirA[0][1].get() + offA[0][1].get() ;
1379 if ( prow->getNumTerm() > 1 ) {
1380 //srate[0] = dirA[1][0].get() ;
1381 //srate[1] = dirA[1][1].get() ;
1382 srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1383 srate[1] = dirA[1][1].get() + offA[1][0].get() ;
1384 }
1385 }
1386 else {
1387 //az = dirA[numSample-1][0].get() ;
1388 //el = dirA[numSample-1][1].get() ;
1389 az = dirA[numSample-1][0].get() + offA[numSample-1][0].get() ;
1390 el = dirA[numSample-1][1].get() + offA[numSample-1][1].get() ;
1391 }
1392 }
1393 else {
1394 ArrayTime startTime = tint.getStart() ;
1395 ArrayTime endTime( tint.getStartInMJD()+tint.getDurationInDays() ) ;
1396 int row0 = -1 ;
1397 int row1 = -1 ;
1398 int row2 = -1 ;
1399 double dt0 = getMidTime( tint ).getMJD() ;
1400 for ( unsigned int irow = 0 ; irow < nrow ; irow++ ) {
1401 prow = (*prows)[irow] ;
1402 double dt = getMidTime( tint ).getMJD() - getMidTime( prow->getTimeInterval() ).getMJD() ;
1403 if ( dt > 0 && dt < dt0 ) {
1404 dt0 = dt ;
1405 row2 = irow ;
1406 }
1407 if ( prow->getTimeInterval().contains( startTime ) )
1408 row0 = irow ;
1409 else if ( prow->getTimeInterval().contains( endTime ) )
1410 row1 = irow ;
1411 if ( row0 != -1 && row1 != -1 )
1412 break ;
1413 }
1414 //logsink_->postLocally( LogMessage("row0 = "+String::toString(row0)+", row1 = "+String::toString(row1)+", row2 = "+String::toString(row2),LogOrigin(className_,funcName,WHERE)) ) ;
1415 if ( row0 == -1 && row1 == -1 ) {
1416 prow = (*prows)[row2] ;
1417 qrow = (*prows)[row2+1] ;
1418 double t0 = getEndTime( prow->getTimeInterval() ).getMJD() ;
1419 double t1 = qrow->getTimeInterval().getStartInMJD() ;
1420 double t = startTime.getMJD() ;
1421 //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1422 //vector< vector<Angle> > dirB = qrow->getPointingDirection() ;
1423 vector< vector<Angle> > dirA = prow->getTarget() ;
1424 vector< vector<Angle> > offA = prow->getOffset() ;
1425 vector< vector<Angle> > dirB = qrow->getTarget() ;
1426 vector< vector<Angle> > offB = qrow->getOffset() ;
1427 //double da0 = dirA[0][0].get() ;
1428 //double db0 = dirB[0][0].get() ;
1429 //double da1 = dirA[0][1].get() ;
1430 //double db1 = dirB[0][1].get() ;
1431 double da0 = dirA[0][0].get() + offA[0][0].get() ;
1432 double db0 = dirB[0][0].get() + offB[0][0].get() ;
1433 double da1 = dirA[0][1].get() + offA[0][1].get() ;
1434 double db1 = dirB[0][1].get() + offB[0][1].get() ;
1435 if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1436 double dt = ( t - t0 ) / ( t1 - t0 ) ;
1437 az = da0 + ( db0 - da0 ) * dt ;
1438 el = da1 + ( db1 - da1 ) * dt ;
1439 if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1440 //double ra0 = dirA[1][0].get() ;
1441 //double rb0 = dirB[1][0].get() ;
1442 //double ra1 = dirA[1][1].get() ;
1443 //double rb1 = dirB[1][1].get() ;
1444 double ra0 = dirA[1][0].get() + offA[1][0].get() ;
1445 double rb0 = dirB[1][0].get() + offB[1][0].get() ;
1446 double ra1 = dirA[1][1].get() + offA[1][1].get() ;
1447 double rb1 = dirB[1][1].get() + offB[1][1].get() ;
1448 srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1449 srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1450 }
1451 }
1452 else if ( !(qrow->getUsePolynomials()) ) {
1453 double dt = ( t - t0 ) / ( t1 - t0 ) ;
1454 az = da0 + ( db0 - da0 ) * dt ;
1455 el = da1 + ( db1 - da1 ) * dt ;
1456 }
1457 //else {
1458 // nothing to do
1459 //}
1460 }
1461 else {
1462 int ndir = 0 ;
1463 if ( row0 == -1 ) {
1464 row0 = row1 ;
1465 }
1466 else if ( row1 == -1 ) {
1467 row1 = row0 ;
1468 }
1469 prow = (*prows)[row0] ;
1470 qrow = (*prows)[row1] ;
1471 if ( prow->getUsePolynomials() && qrow->getUsePolynomials() ) {
1472 //logsink_->postLocally( LogMessage("usePolynomial = True",LogOrigin(className_,funcName,WHERE)) ) ;
1473 if ( row0 == row1 ) {
1474 prow = (*prows)[row0] ;
1475 //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1476 vector< vector<Angle> > dirA = prow->getTarget() ;
1477 vector< vector<Angle> > offA = prow->getOffset() ;
1478 //az = dirA[0][0].get() ;
1479 //el = dirA[0][1].get() ;
1480 az = dirA[0][0].get() + offA[0][0].get() ;
1481 el = dirA[0][1].get() + offA[0][1].get() ;
1482 if ( prow->getNumTerm() > 1 ) {
1483 //srate[0] = dirA[1][0].get() ;
1484 //srate[1] = dirA[1][1].get() ;
1485 srate[0] = dirA[1][0].get() + offA[1][0].get() ;
1486 srate[1] = dirA[1][1].get() + offA[1][0].get() ;
1487 }
1488 }
1489 else {
1490 prow = (*prows)[row0] ;
1491 qrow = (*prows)[row1] ;
1492 //vector< vector<Angle> > dirA = qrow->getPointingDirection() ;
1493 //vector< vector<Angle> > dirB = qrow->getPointingDirection() ;
1494 vector< vector<Angle> > dirA = qrow->getTarget() ;
1495 vector< vector<Angle> > offA = prow->getOffset() ;
1496 vector< vector<Angle> > dirB = qrow->getTarget() ;
1497 vector< vector<Angle> > offB = qrow->getOffset() ;
1498 double t0 = getMidTime( prow->getTimeInterval() ).getMJD() ;
1499 double t1 = getMidTime( qrow->getTimeInterval() ).getMJD() ;
1500 double t = startTime.getMJD() ;
1501 double dt = ( t - t0 ) / ( t1 - t0 ) ;
1502 //double da0 = dirA[0][0].get() ;
1503 //double db0 = dirB[0][0].get() ;
1504 //double da1 = dirA[0][1].get() ;
1505 //double db1 = dirB[0][1].get() ;
1506 double da0 = dirA[0][0].get() + offA[0][0].get() ;
1507 double db0 = dirB[0][0].get() + offB[0][0].get() ;
1508 double da1 = dirA[0][1].get() + offA[0][1].get() ;
1509 double db1 = dirB[0][1].get() + offB[0][1].get() ;
1510 az = da0 + ( db0 - da0 ) * dt ;
1511 el = da1 + ( db1 - db0 ) * dt ;
1512 if ( prow->getNumTerm() > 0 && qrow->getNumTerm() > 1 ) {
1513 //double ra0 = dirA[1][0].get() ;
1514 //double rb0 = dirB[1][0].get() ;
1515 //double ra1 = dirA[1][1].get() ;
1516 //double rb1 = dirB[1][1].get() ;
1517 double ra0 = dirA[1][0].get() + offA[1][0].get() ;
1518 double rb0 = dirB[1][0].get() + offB[1][0].get() ;
1519 double ra1 = dirA[1][1].get() + offA[1][1].get() ;
1520 double rb1 = dirB[1][1].get() + offB[1][1].get() ;
1521 srate[0] = ra0 + ( rb0 - ra0 ) * dt ;
1522 srate[1] = ra1 + ( rb1 - ra1 ) * dt ;
1523 }
1524 }
1525 }
1526 else if ( prow->getUsePolynomials() == qrow->getUsePolynomials() ) {
1527 //logsink_->postLocally( LogMessage("numSample == numTerm",LogOrigin(className_,funcName,WHERE)) ) ;
1528 for ( int irow = row0 ; irow <= row1 ; irow++ ) {
1529 prow = (*prows)[irow] ;
1530 int numSample = prow->getNumSample() ;
1531 //logsink_->postLocally( LogMessage("numSample = "+String::toString(numSample),LogOrigin(className_,funcName,WHERE)) ) ;
1532 //vector< vector<Angle> > dirA = prow->getPointingDirection() ;
1533 vector< vector<Angle> > dirA = prow->getTarget() ;
1534 vector< vector<Angle> > offA = prow->getOffset() ;
1535 if ( prow->isSampledTimeIntervalExists() ) {
1536 //logsink_->postLocally( LogMessage("sampledTimeIntervalExists",LogOrigin(className_,funcName,WHERE)) ) ;
1537 vector<ArrayTimeInterval> stime = prow->getSampledTimeInterval() ;
1538 for ( int isam = 0 ; isam < numSample ; isam++ ) {
1539 //if ( tint.overlaps( stime[isam] ) ) {
1540 if ( tint.contains( stime[isam] ) ) {
1541 //az += dirA[isam][0].get() ;
1542 //el += dirA[isam][1].get() ;
1543 az += dirA[isam][0].get() + offA[isam][0].get() ;
1544 el += dirA[isam][1].get() + offA[isam][1].get() ;
1545 ndir++ ;
1546 }
1547 }
1548 }
1549 else {
1550 double sampleStart = prow->getTimeInterval().getStartInMJD() ;
1551 double sampleInterval = prow->getTimeInterval().getDurationInDays() / (double)numSample ;
1552 //logsink_->postLocally( LogMessage("sampleStart = "+String::toString(sampleStart),LogOrigin(className_,funcName,WHERE)) )
1553 //logsink_->postLocally( LogMessage("sampleInterval = "+String::toString(sampleInterval),LogOrigin(className_,funcName,WHERE)) ) ;
1554 //logsink_->postLocally( LogMessage("tint = "+tint.toString(),LogOrigin(className_,funcName,WHERE)) ) ;
1555 for ( int isam = 0 ; isam < numSample ; isam++ ) {
1556 ArrayTimeInterval stime( sampleStart+isam*sampleInterval, sampleInterval ) ;
1557 //if ( tint.overlaps( stime ) ) {
1558 if ( tint.contains( stime ) ) {
1559 //az += dirA[isam][0].get() ;
1560 //el += dirA[isam][1].get() ;
1561 az += dirA[isam][0].get() + offA[isam][0].get() ;
1562 el += dirA[isam][1].get() + offA[isam][1].get() ;
1563 tcen += getMidTime( stime ).getMJD() ;
1564 ndir++ ;
1565 }
1566 }
1567 }
1568 }
1569 if ( ndir > 1 ) {
1570 az /= (double)ndir ;
1571 el /= (double)ndir ;
1572 tcen /= (double)ndir ;
1573 }
1574 }
1575 //else {
1576 // nothing to do
1577 //}
1578 }
1579
1580 AntennaRow *arow = asdm_->getAntenna().getRowByKey( Tag( antennaId_, TagType::Antenna ) ) ;
1581 StationRow *srow = arow->getStationUsingStationId() ;
1582 vector<Length> antposL = srow->getPosition() ;
1583 casa::Vector<casa::Double> antpos( 3 ) ;
1584 for ( int i = 0 ; i < 3 ; i++ )
1585 antpos[i] = antposL[i].get() ;
1586 //logsink_->postLocally( LogMessage("tcen = "+String::toString(tcen),LogOrigin(className_,funcName,WHERE)) ) ;
1587 //logsink_->postLocally( LogMessage("antpos = "+String::toString(antpos),LogOrigin(className_,funcName,WHERE)) ) ;
1588 toJ2000( dir, az, el, tcen, antpos ) ;
1589
1590 }
1591
1592 return ;
1593}
1594
1595ArrayTime ASDMReader::getMidTime( const ArrayTimeInterval &t )
1596{
1597 return ArrayTime( t.getStartInMJD() + 0.5 * t.getDurationInDays() ) ;
1598}
1599
1600ArrayTime ASDMReader::getEndTime( const ArrayTimeInterval &t )
1601{
1602 return ArrayTime( t.getStartInMJD() + t.getDurationInDays() ) ;
1603}
1604
1605ArrayTime ASDMReader::getStartTime( const ArrayTimeInterval &t )
1606{
1607 return ArrayTime( t.getStartInMJD() ) ;
1608}
1609
1610void ASDMReader::toJ2000( vector<double> &dir,
1611 double az,
1612 double el,
1613 double mjd,
1614 casa::Vector<casa::Double> antpos )
1615{
1616 String funcName = "toJ2000" ;
1617
1618 casa::Vector<casa::Double> azel( 2 ) ;
1619 azel[0] = az ;
1620 azel[1] = el ;
1621 casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
1622 casa::Vector<casa::Quantity> qantpos( 3 ) ;
1623 qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
1624 qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
1625 qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
1626 casa::MPosition mp( casa::MVPosition( qantpos ),
1627 casa::MPosition::ITRF ) ;
1628 //ostringstream oss ;
1629 //mp.print( oss ) ;
1630 //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
1631
1632 casa::MeasFrame mf( me, mp ) ;
1633 casa::MDirection::Convert toj2000( casa::MDirection::AZELGEO,
1634 //MDirection::Convert toj2000( MDirection::AZEL,
1635 //MDirection::Convert toj2000( MDirection::AZELSW,
1636 //MDirection::Convert toj2000( MDirection::AZELSWGEO,
1637 //MDirection::Convert toj2000( MDirection::AZELNE,
1638 //MDirection::Convert toj2000( MDirection::AZELNEGEO,
1639 casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
1640 casa::Vector<casa::Double> cdir = toj2000( azel ).getAngle( "rad" ).getValue() ;
1641 //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
1642 dir.resize(2) ;
1643 dir[0] = (double)(cdir[0]) ;
1644 dir[1] = (double)(cdir[1]) ;
1645}
1646
1647void ASDMReader::setLogger( CountedPtr<LogSinkInterface> &logsink )
1648{
1649 logsink_ = logsink ;
1650}
Note: See TracBrowser for help on using the repository browser.