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

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

New Development: Yes

JIRA Issue: Yes CAS-1913

Ready for Test: No

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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