source: trunk/external-alma/oldasdm2ASAP/OldASDMReader.cc@ 2299

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

Do not access Subscan/SubscanMode if it doesn't exist.


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