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

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

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