source: trunk/external-alma/atnf/PKSIO/NROReader.cc@ 3103

Last change on this file since 3103 was 3097, checked in by Takeshi Nakazato, 9 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

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...


Bug fix on string handling in NRO data filler. Remove trailing null character ('\0').

File size: 24.2 KB
Line 
1//#---------------------------------------------------------------------------
2//# NROReader.cc: Base class for NRO headerdata.
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2000-2006
5//# Associated Universities, Inc. Washington DC, USA.
6//#
7//# This library is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU Library General Public License as published by
9//# the Free Software Foundation; either version 2 of the License, or (at your
10//# option) any later version.
11//#
12//# This library is distributed in the hope that it will be useful, but WITHOUT
13//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15//# License for more details.
16//#
17//# You should have received a copy of the GNU Library General Public License
18//# along with this library; if not, write to the Free Software Foundation,
19//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning AIPS++ should be addressed as follows:
22//# Internet email: aips2-request@nrao.edu.
23//# Postal address: AIPS++ Project Office
24//# National Radio Astronomy Observatory
25//# 520 Edgemont Road
26//# Charlottesville, VA 22903-2475 USA
27//#
28//# $Id$
29//#---------------------------------------------------------------------------
30//# Original: 2008/10/30, Takeshi Nakazato, NAOJ
31//#---------------------------------------------------------------------------
32
33#include <atnf/PKSIO/NROReader.h>
34#include <atnf/PKSIO/NRO45Reader.h>
35#include <atnf/PKSIO/ASTEReader.h>
36#include <atnf/PKSIO/ASTEFXReader.h>
37#include <atnf/PKSIO/NRO45FITSReader.h>
38#include <atnf/PKSIO/NROOTFDataset.h>
39#include <atnf/PKSIO/ASTEDataset.h>
40
41#include <measures/Measures/MEpoch.h>
42#include <measures/Measures/MPosition.h>
43#include <measures/Measures/MDirection.h>
44#include <measures/Measures/MCDirection.h>
45#include <measures/Measures/MeasFrame.h>
46#include <measures/Measures/MeasConvert.h>
47
48#include <casa/Logging/LogIO.h>
49#include <casa/IO/RegularFileIO.h>
50#include <casa/OS/File.h>
51#include <casa/OS/Time.h>
52
53#include <stdio.h>
54#include <string>
55#include <iomanip>
56
57using namespace std ;
58
59namespace {
60// to limit template argument to specified type
61template<class T>
62constexpr inline bool type_guard() {
63 return false;
64}
65template<>
66constexpr inline bool type_guard<string>() {
67 return true;
68}
69template<>
70constexpr inline bool type_guard<String>() {
71 return true;
72}
73
74// T must be std::string or casa::String
75template<class T>
76inline String trim_nro_string(T const &s) {
77 static_assert(type_guard<T>(), "");
78 return s.substr(0, s.find_first_of('\0'));
79}
80}
81
82//
83// getNROReader
84//
85// Return an appropriate NROReader for a NRO 45m and ASTE dataset.
86//
87NROReader *getNROReader( const String filename,
88 String &datatype )
89{
90 LogIO os( LogOrigin( "", "getNROReader()", WHERE ) ) ;
91
92 // Check accessibility of the input.
93 File inFile( filename ) ;
94 if ( !inFile.exists() ) {
95 datatype = filename + " not found." ;
96 return 0 ;
97 }
98
99 if ( !inFile.isReadable() ) {
100 datatype = filename + " is not readable." ;
101 return 0 ;
102 }
103
104 // Determine the type of input.
105 NROReader *reader = 0;
106 if ( inFile.isRegular() ) {
107 FILE *file ;
108 file = fopen( filename.c_str(), "r" ) ;
109 // read LOFIL0
110 char buf[9];
111 size_t retval = fread( buf, 4, 1, file ) ;
112 if (retval < 1) {
113 os << LogIO::SEVERE << "Failed to read data" << LogIO::POST;
114 fclose(file);
115 return 0;
116 }
117 buf[4] = '\0' ;
118 // DEBUG
119 //os << LogIO::NORMAL << "getNROReader:: buf = " << String(buf) << LogIO::POST ;
120 //
121 if ( string( buf ) == "XTEN" ) {
122 // FITS data
123 datatype = "NRO 45m FITS" ;
124 reader = new NRO45FITSReader( filename ) ;
125 }
126 else if ( string( buf ) == "RW-F") {
127 // ASTE-FX data
128 datatype = "ASTE-FX";
129 reader = new ASTEFXReader( filename );
130 } else {
131 // otherwise, read SITE0
132 NRODataset *d = new NROOTFDataset( filename ) ;
133 d->initialize() ;
134 int size = d->getDataSize() - 188 ;
135 delete d ;
136 fseek( file, size, SEEK_SET ) ;
137 size_t retval = fread( buf, 8, 1, file ) ;
138 if (retval < 1) {
139 os << LogIO::SEVERE << "Failed to read data" << LogIO::POST;
140 fclose(file);
141 return 0;
142 }
143 buf[8] = '\0' ;
144 // DEBUG
145 //cout << "getNROReader:: buf = " << buf << endl ;
146 //
147 if ( string( buf ) == "NRO" ) {
148 // NRO 45m data
149 datatype = "NRO 45m OTF" ;
150 reader = new NRO45Reader( filename );
151 }
152 else {
153 d = new ASTEDataset( filename ) ;
154 d->initialize() ;
155 size = d->getDataSize() - 188 ;
156 delete d ;
157 fseek( file, size, SEEK_SET ) ;
158 size_t retval = fread( buf, 8, 1, file ) ;
159 if (retval < 1) {
160 os << LogIO::SEVERE << "Failed to read data" << LogIO::POST;
161 fclose(file);
162 return 0;
163 }
164 buf[8] = '\0' ;
165 // DEBUG
166 //cout << "getNROReader:: buf = " << buf << endl ;
167 //
168 if ( string( buf ) == "ASTE" ) {
169 // ASTE data
170 datatype = "ASTE" ;
171 reader = new ASTEReader( filename ) ;
172 }
173 else {
174 datatype = "UNRECOGNIZED INPUT FORMAT";
175 }
176 }
177 }
178 fclose( file ) ;
179 } else {
180 datatype = "UNRECOGNIZED INPUT FORMAT";
181 }
182
183 // DEBUG
184 os << LogIO::NORMAL << "Data format of " << filename << ": " << datatype << LogIO::POST ;
185 //
186
187 // return reader if exists
188 if ( reader ) {
189 reader->read() ;
190 return reader ;
191 }
192
193 return 0 ;
194}
195
196
197//
198// getNROReader
199//
200// Search a list of directories for a NRO 45m and ASTE dataset and return an
201// appropriate NROReader for it.
202//
203NROReader* getNROReader( const String name,
204 const Vector<String> directories,
205 int &iDir,
206 String &datatype )
207{
208 int nDir = directories.size();
209 for ( iDir = 0; iDir < nDir; iDir++ ) {
210 string inName = directories[iDir] + "/" + name;
211 NROReader *reader = getNROReader( inName, datatype ) ;
212
213 if (reader != 0) {
214 return reader;
215 }
216 }
217
218 iDir = -1;
219 return 0;
220}
221
222
223//----------------------------------------------------------------------
224// constructor
225NROReader::NROReader( string name )
226 : dataset_( NULL ),
227 srcdir_( 0 ),
228 msrcdir_( 0 ),
229 freqRefFromVREF_( false ) // import frequency as REST
230{
231 // initialization
232 filename_ = name ;
233 coord_ = -1 ;
234}
235
236// destructor
237NROReader::~NROReader()
238{
239}
240
241// Read data header
242Int NROReader::read()
243{
244 LogIO os( LogOrigin( "NROReader", "read()", WHERE ) ) ;
245
246 int status = 0 ;
247
248 // initialize Dataset object
249 initDataset();
250
251 // fill Dataset
252 status = dataset_->fillHeader() ;
253
254 if ( status != 0 ) {
255 os << LogIO::SEVERE << "Failed to fill data header." << LogIO::EXCEPTION ;
256 }
257
258 return status ;
259}
260
261// set frequency reference frame
262void NROReader::setFreqRefFromVREF( bool fromVREF )
263{
264 os_.origin( LogOrigin( "NROReader", "setFreqRefFromVREF", WHERE ) ) ;
265 os_ << ((fromVREF) ? "Take frequency reference frame from VREF" :
266 "Use frequency reference frame REST") << LogIO::POST ;
267
268 freqRefFromVREF_ = fromVREF ;
269}
270
271// get spectrum
272vector< vector<double> > NROReader::getSpectrum()
273{
274 return dataset_->getSpectrum() ;
275}
276
277Int NROReader::getPolarizationNum()
278{
279 return dataset_->getPolarizationNum() ;
280}
281
282double NROReader::getStartTime()
283{
284 //char *startTime = dataset_->getLOSTM() ;
285 string startTime = dataset_->getLOSTM() ;
286 //cout << "getStartTime() startTime = " << startTime << endl ;
287 return getMJD( startTime ) ;
288}
289
290double NROReader::getEndTime()
291{
292 //char *endTime = dataset_->getLOETM() ;
293 string endTime = dataset_->getLOETM() ;
294 return getMJD( endTime ) ;
295}
296
297vector<double> NROReader::getStartIntTime()
298{
299 return dataset_->getStartIntTime() ;
300}
301
302double NROReader::getMJD( char *time )
303{
304 // TODO: should be checked which time zone the time depends on
305 // 2008/11/14 Takeshi Nakazato
306 string strStartTime = string( time ) ;
307 string strYear = strStartTime.substr( 0, 4 ) ;
308 string strMonth = strStartTime.substr( 4, 2 ) ;
309 string strDay = strStartTime.substr( 6, 2 ) ;
310 string strHour = strStartTime.substr( 8, 2 ) ;
311 string strMinute = strStartTime.substr( 10, 2 ) ;
312 string strSecond = strStartTime.substr( 12, strStartTime.size() - 12 ) ;
313 uInt year = atoi( strYear.c_str() ) ;
314 uInt month = atoi( strMonth.c_str() ) ;
315 uInt day = atoi( strDay.c_str() ) ;
316 uInt hour = atoi( strHour.c_str() ) ;
317 uInt minute = atoi( strMinute.c_str() ) ;
318 double second = atof( strSecond.c_str() ) ;
319 Time t( year, month, day, hour, minute, second ) ;
320
321 return t.modifiedJulianDay() ;
322}
323
324double NROReader::getMJD( string strStartTime )
325{
326 // TODO: should be checked which time zone the time depends on
327 // 2008/11/14 Takeshi Nakazato
328 string strYear = strStartTime.substr( 0, 4 ) ;
329 string strMonth = strStartTime.substr( 4, 2 ) ;
330 string strDay = strStartTime.substr( 6, 2 ) ;
331 string strHour = strStartTime.substr( 8, 2 ) ;
332 string strMinute = strStartTime.substr( 10, 2 ) ;
333 string strSecond = strStartTime.substr( 12, strStartTime.size() - 12 ) ;
334 uInt year = atoi( strYear.c_str() ) ;
335 uInt month = atoi( strMonth.c_str() ) ;
336 uInt day = atoi( strDay.c_str() ) ;
337 uInt hour = atoi( strHour.c_str() ) ;
338 uInt minute = atoi( strMinute.c_str() ) ;
339 double second = atof( strSecond.c_str() ) ;
340 Time t( year, month, day, hour, minute, second ) ;
341
342 return t.modifiedJulianDay() ;
343}
344
345vector<Bool> NROReader::getIFs()
346{
347 return dataset_->getIFs() ;
348}
349
350vector<Bool> NROReader::getBeams()
351{
352 vector<Bool> v ;
353 vector<int> arry = dataset_->getARRY() ;
354 for ( uInt i = 0 ; i < arry.size() ; i++ ) {
355 if ( arry[i] != 0 ) {
356 v.push_back( True ) ;
357 }
358 }
359
360 // DEBUG
361 //cout << "NROReader::getBeams() number of beam is " << v.size() << endl ;
362 //
363
364 return v ;
365}
366
367// Get SRCDIRECTION in RADEC(J2000)
368Vector<Double> NROReader::getSourceDirection()
369{
370 if ( msrcdir_.nelements() == 2 )
371 return msrcdir_ ;
372
373 srcdir_.resize( 2 ) ;
374 srcdir_[0] = Double( dataset_->getRA0() ) ;
375 srcdir_[1] = Double( dataset_->getDEC0() ) ;
376 char epoch[5] ;
377 strncpy( epoch, (dataset_->getEPOCH()).c_str(), 5 ) ;
378 if ( strncmp( epoch, "B1950", 5 ) == 0 ) {
379 // convert to J2000 value
380 MDirection result =
381 MDirection::Convert( MDirection( Quantity( srcdir_[0], "rad" ),
382 Quantity( srcdir_[1], "rad" ),
383 MDirection::Ref( MDirection::B1950 ) ),
384 MDirection::Ref( MDirection::J2000 ) ) () ;
385 msrcdir_ = result.getAngle().getValue() ;
386 if ( msrcdir_[0] < 0.0 && srcdir_[0] >= 0.0 )
387 msrcdir_[0] = 2.0 * M_PI + msrcdir_[0] ;
388 //cout << "NROReader::getSourceDirection() SRCDIRECTION convert from ("
389 //<< srcra << "," << srcdec << ") B1950 to ("
390 //<< v( 0 ) << ","<< v( 1 ) << ") J2000" << endl ;
391 //LogIO os( LogOrigin( "NROReader", "getSourceDirection()", WHERE ) ) ;
392 //os << LogIO::NORMAL << "SRCDIRECTION convert from ("
393 // << srcra << "," << srcdec << ") B1950 to ("
394 // << v( 0 ) << ","<< v( 1 ) << ") J2000" << LogIO::POST ;
395 }
396 else if ( strncmp( epoch, "J2000", 5 ) == 0 ) {
397 msrcdir_.reference( srcdir_ ) ;
398 }
399
400 return msrcdir_ ;
401}
402
403// Get DIRECTION in RADEC(J2000)
404Vector<Double> NROReader::getDirection( int i )
405{
406 Vector<Double> v( 2 ) ;
407 const NRODataRecord *record = dataset_->getRecord( i ) ;
408 char epoch[5] ;
409 strncpy( epoch, (dataset_->getEPOCH()).c_str(), 5 ) ;
410 int icoord = dataset_->getSCNCD() ;
411 double scantime = dataset_->getScanTime( i ) ;
412 initConvert( icoord, scantime, epoch ) ;
413 v[0]= Double( record->SCX ) ;
414 v[1] = Double( record->SCY ) ;
415 if ( converter_.null() )
416 // no conversion
417 return v ;
418 else {
419 Vector<Double> v2 = (*converter_)( v ).getAngle().getValue() ;
420 if ( v2[0] < 0.0 && v[0] >= 0.0 )
421 v2[0] = 2.0 * M_PI + v2[0] ;
422 return v2 ;
423 }
424}
425
426void NROReader::initConvert( int icoord, double t, char *epoch )
427{
428 if ( icoord == 0 && strncmp( epoch, "J2000", 5 ) == 0 )
429 // no conversion
430 return ;
431
432 if ( converter_.null() || icoord != coord_ ) {
433 LogIO os( LogOrigin( "NROReader", "initConvert()", WHERE ) ) ;
434 coord_ = icoord ;
435 if ( coord_ == 0 ) {
436 // RADEC (B1950) -> RADEC (J2000)
437 os << "Creating converter from RADEC (B1950) to RADEC (J2000)" << LogIO::POST ;
438 converter_ = new MDirection::Convert( MDirection::B1950,
439 MDirection::J2000 ) ;
440 }
441 else if ( coord_ == 1 ) {
442 // LB -> RADEC (J2000)
443 os << "Creating converter from GALACTIC to RADEC (J2000)" << LogIO::POST ;
444 converter_ = new MDirection::Convert( MDirection::GALACTIC,
445 MDirection::J2000 ) ;
446 }
447 else {
448 // coord_ == 2
449 // AZEL -> RADEC (J2000)
450 os << "Creating converter from AZEL to RADEC (J2000)" << LogIO::POST ;
451 if ( mf_.null() ) {
452 mf_ = new MeasFrame() ;
453 vector<double> antpos = getAntennaPosition() ;
454 Vector<Quantity> qantpos( 3 ) ;
455 for ( int ip = 0 ; ip < 3 ; ip++ )
456 qantpos[ip] = Quantity( antpos[ip], "m" ) ;
457 mp_ = MPosition( MVPosition( qantpos ), MPosition::ITRF ) ;
458 mf_->set( mp_ ) ;
459 }
460 converter_ = new MDirection::Convert( MDirection::AZEL,
461 MDirection::Ref(MDirection::J2000,
462 *mf_ ) ) ;
463 }
464 }
465
466 if ( coord_ == 2 ) {
467 me_ = MEpoch( Quantity( t, "d" ), MEpoch::UTC ) ;
468 mf_->set( me_ ) ;
469 }
470}
471
472int NROReader::getHeaderInfo( Int &nchan,
473 Int &npol,
474 Int &nif,
475 Int &nbeam,
476 String &observer,
477 String &project,
478 String &obstype,
479 String &antname,
480 Vector<Double> &antpos,
481 Float &equinox,
482 String &freqref,
483 Double &reffreq,
484 Double &bw,
485 Double &utc,
486 String &fluxunit,
487 String &epoch,
488 String &poltype )
489{
490 nchan = dataset_->getNUMCH() ;
491 //cout << "nchan = " << nchan << endl ;
492 npol = getPolarizationNum() ;
493 //cout << "npol = " << npol << endl ;
494 observer = trim_nro_string(dataset_->getOBSVR()) ;
495 //cout << "observer = \"" << observer << "\" (size " << observer.size() << ")" << endl ;
496 project = trim_nro_string(dataset_->getPROJ()) ;
497 //cout << "project = \"" << project << "\" (size " << project.size() << ")" << endl ;
498 obstype = trim_nro_string(dataset_->getSWMOD()) ;
499 //cout << "obstype = \"" << obstype << "\" (size " << obstype.size() << ")" << endl ;
500 antname = trim_nro_string(dataset_->getSITE()) ;
501 //cout << "antname = \"" << antname << "\" (size " << antname.size() << ")" << endl ;
502 // TODO: should be investigated antenna position since there are
503 // no corresponding information in the header
504 // 2008/11/13 Takeshi Nakazato
505 //
506 // INFO: tentative antenna posiiton is obtained for NRO 45m from ITRF website
507 // 2008/11/26 Takeshi Nakazato
508 vector<double> pos = getAntennaPosition() ;
509 antpos = pos ;
510 //cout << "antpos = " << antpos << endl ;
511 string eq = dataset_->getEPOCH() ;
512// if ( eq.compare( 0, 5, "B1950" ) == 0 )
513// equinox = 1950.0 ;
514// else if ( eq.compare( 0, 5, "J2000" ) == 0 )
515// equinox = 2000.0 ;
516 // equinox is always 2000.0
517 equinox = 2000.0 ;
518 //cout << "equinox = " << equinox << endl ;
519 string vref = dataset_->getVREF() ;
520 if ( vref.compare( 0, 3, "LSR" ) == 0 ) {
521 if ( vref.size() == 3 ) {
522 vref.append( "K" ) ;
523 }
524 else {
525 vref[3] = 'K' ;
526 }
527 }
528 else if ( vref.compare( 0, 3, "GAL" ) == 0 ) {
529 vref = "GALACTO" ;
530 }
531 else if (vref.compare( 0, 3, "HEL" ) == 0 ) {
532 os_.origin( LogOrigin( "NROReader", "getHeaderInfo", WHERE ) ) ;
533 os_ << LogIO::WARN << "Heliocentric frame is not supported. Use Barycentric frame instead." << LogIO::POST ;
534 vref = "BARY" ;
535 }
536 //freqref = vref ;
537 //freqref = "LSRK" ;
538 //freqref = "REST" ;
539 freqref = freqRefFromVREF_ ? trim_nro_string(vref) : "REST" ;
540 //cout << "freqref = \"" << freqref << "\" (size " << freqref.size() << ")" << endl ;
541 const NRODataRecord *record = dataset_->getRecord( 0 ) ;
542 reffreq = record->FREQ0 ; // rest frequency
543
544 //cout << "reffreq = " << reffreq << endl ;
545 bw = dataset_->getBEBW()[0] ;
546 //cout << "bw = " << bw << endl ;
547 utc = getStartTime() ;
548 //cout << "utc = " << utc << endl ;
549 fluxunit = "K" ;
550 //cout << "fluxunit = " << fluxunit << endl ;
551 epoch = "UTC" ;
552 //cout << "epoch = " << epoch << endl ;
553 string poltp = trim_nro_string(dataset_->getPOLTP()[0]) ;
554 //cout << "poltp = \"" << poltp << "\" (size " << poltp.size() << ")" << endl ;
555 if ( poltp.empty() || poltp[0] == ' ' || poltp[0] == '\0' )
556 //poltp = "None" ;
557 poltp = "linear" ; // if no polarization type specified, set to "linear"
558 //else if ( strcmp( poltp, "LINR" ) == 0 )
559 else if ( poltp.compare( 0, 1, "LINR", 0, 1 ) == 0 )
560 poltp = "linear" ;
561 //else if ( strcmp( poltp, "CIRL" ) == 0 )
562 else if ( poltp.compare( 0, 1, "CIRL", 0, 1 ) == 0 )
563 poltp = "circular" ;
564 poltype = poltp ;
565 //cout << "poltype = '" << poltype << "'" << endl ;
566
567 //vector<Bool> ifs = getIFs() ;
568 //nif = ifs.size() ;
569 nif = getNumIF() ;
570 //cout << "nif = " << nif << endl ;
571
572 //vector<Bool> beams = getBeams() ;
573 //nbeam = beams.size() ;
574 nbeam = getNumBeam() ;
575 //cout << "nbeam = " << nbeam << endl ;
576
577 return 0 ;
578}
579
580string NROReader::getScanType( int i )
581{
582 const NRODataRecord *record = dataset_->getRecord( i ) ;
583 string s = record->SCANTP ;
584
585 return s ;
586}
587
588int NROReader::getScanInfo( int irow,
589 uInt &scanno,
590 uInt &cycleno,
591 uInt &ifno,
592 uInt &beamno,
593 uInt &polno,
594 vector<double> &freqs,
595 Vector<Double> &restfreq,
596 uInt &refbeamno,
597 Double &scantime,
598 Double &interval,
599 String &srcname,
600 String &fieldname,
601 Vector<Float> &spectra,
602 Vector<uChar> &flagtra,
603 Vector<Float> &tsys,
604 Vector<Double> &direction,
605 Float &azimuth,
606 Float &elevation,
607 Float &parangle,
608 Float &opacity,
609 uInt &tcalid,
610 Int &fitid,
611 uInt &focusid,
612 Float &temperature,
613 Float &pressure,
614 Float &humidity,
615 Float &windvel,
616 Float &winddir,
617 Double &srcvel,
618 Vector<Double> &/*propermotion*/,
619 Vector<Double> &srcdir,
620 Vector<Double> &/*scanrate*/ )
621{
622 static const IPosition oneByOne( 1, 1 );
623
624 // DEBUG
625 //cout << "NROReader::getScanInfo() irow = " << irow << endl ;
626 //
627 NRODataRecord *record = dataset_->getRecord( irow ) ;
628
629 // scanno
630 scanno = (uInt)(record->ISCAN) ;
631 //cout << "scanno = " << scanno << endl ;
632
633 // cycleno
634 cycleno = 0 ;
635 //cout << "cycleno = " << cycleno << endl ;
636
637 // beamno and ifno
638 string rxname = dataset_->getRX()[0] ;
639 if ( rxname.find("MULT2") != string::npos ) {
640 string arryt = string( record->ARRYT ) ;
641 beamno = dataset_->getSortedArrayId( arryt ) ;
642 ifno = 0 ;
643 }
644 else {
645 beamno = 0 ;
646 string arryt = string( record->ARRYT ) ;
647 ifno = dataset_->getSortedArrayId( arryt ) ;
648 }
649 //cout << "beamno = " << beamno << endl ;
650
651 // polno
652 polno = dataset_->getPolNo( irow ) ;
653 //cout << "polno = " << polno << endl ;
654
655 // freqs (for IFNO and FREQ_ID)
656 //freqs = getFrequencies( irow ) ;
657 freqs = dataset_->getFrequencies( irow ) ;
658 //cout << "freqs = [" << freqs[0] << ", " << freqs[1] << ", " << freqs[2] << "]" << endl ;
659
660 if ( freqRefFromVREF_ ) {
661 freqs = shiftFrequency( freqs,
662 dataset_->getURVEL(),
663 dataset_->getVDEF() ) ;
664 }
665
666 // restfreq (for MOLECULE_ID)
667 restfreq.resize( oneByOne ) ;
668 restfreq[0] = record->FREQ0 ;
669 //cout << "restfreq = " << restfreq[0] << endl ;
670
671 // refbeamno
672 refbeamno = 0 ;
673 //cout << "refbeamno = " << refbeamno << endl ;
674
675 // scantime
676 //scantime = Double( dataset_->getStartIntTime( irow ) ) ;
677 scantime = Double( dataset_->getScanTime( irow ) ) ;
678 //cout << "scantime = " << scantime << endl ;
679
680 // interval
681 interval = Double( dataset_->getIPTIM() ) ;
682 //cout << "interval = " << interval << endl ;
683
684 // srcname
685 srcname = trim_nro_string( dataset_->getOBJ() ) ;
686 //cout << "srcname = \"" << srcname << "\" (size " << srcname.size() << ")" << endl ;
687
688 // fieldname
689 fieldname = trim_nro_string( dataset_->getOBJ() ) ;
690 //cout << "fieldname = \"" << fieldname << "\" (size " << fieldname.size() << ")" << endl ;
691
692 // spectra
693 vector<double> spec = dataset_->getSpectrum( irow ) ;
694 spectra.resize( spec.size() ) ;
695 //int index = 0 ;
696 Bool b ;
697 Float *fp = spectra.getStorage( b ) ;
698 Float *wp = fp ;
699 for ( vector<double>::iterator i = spec.begin() ;
700 i != spec.end() ; i++ ) {
701 *wp = *i ;
702 wp++ ;
703 }
704 spectra.putStorage( fp, b ) ;
705 //cout << "spec.size() = " << spec.size() << endl ;
706
707 // flagtra
708 bool setValue = !( flagtra.nelements() == spectra.nelements() ) ;
709 if ( setValue ) {
710 //cout << "flagtra resized. reset values..." << endl ;
711 flagtra.resize( spectra.nelements() ) ;
712 flagtra.set( 0 ) ;
713 }
714 //cout << "flag.size() = " << flagtra.size() << endl ;
715
716 // tsys
717 tsys.resize( oneByOne ) ;
718 tsys[0] = record->TSYS ;
719 //cout << "tsys[0] = " << tsys[0] << endl ;
720
721 // direction
722 direction = getDirection( irow ) ;
723 //cout << "direction = [" << direction[0] << ", " << direction[1] << "]" << endl ;
724
725 // azimuth
726 azimuth = record->RAZ ;
727 //cout << "azimuth = " << azimuth << endl ;
728
729 // elevation
730 elevation = record->REL ;
731 //cout << "elevation = " << elevation << endl ;
732
733 // parangle
734 parangle = 0.0 ;
735 //cout << "parangle = " << parangle << endl ;
736
737 // opacity
738 opacity = 0.0 ;
739 //cout << "opacity = " << opacity << endl ;
740
741 // tcalid
742 tcalid = 0 ;
743 //cout << "tcalid = " << tcalid << endl ;
744
745 // fitid
746 fitid = -1 ;
747 //cout << "fitid = " << fitid << endl ;
748
749 // focusid
750 focusid = 0 ;
751 //cout << "focusid = " << focusid << endl ;
752
753 // temperature (for WEATHER_ID)
754 temperature = Float( record->TEMP ) ;
755 //cout << "temperature = " << temperature << endl ;
756
757 // pressure (for WEATHER_ID)
758 pressure = Float( record->PATM ) ;
759 //cout << "pressure = " << pressure << endl ;
760
761 // humidity (for WEATHER_ID)
762 humidity = Float( record->PH2O ) ;
763 //cout << "humidity = " << humidity << endl ;
764
765 // windvel (for WEATHER_ID)
766 windvel = Float( record->VWIND ) ;
767 //cout << "windvel = " << windvel << endl ;
768
769 // winddir (for WEATHER_ID)
770 winddir = Float( record->DWIND ) ;
771 //cout << "winddir = " << winddir << endl ;
772
773 // srcvel
774 srcvel = dataset_->getURVEL() ;
775 //cout << "srcvel = " << srcvel << endl ;
776
777 // propermotion
778 // do nothing
779
780 // srcdir
781 srcdir = getSourceDirection() ;
782 //cout << "srcdir = [" << srcdir[0] << ", " << srcdir[1] << endl ;
783
784 // scanrate
785 // do nothing
786
787 return 0 ;
788}
789
790Int NROReader::getRowNum()
791{
792 return dataset_->getRowNum() ;
793}
794
795vector<double> NROReader::shiftFrequency( const vector<double> &f,
796 const double &v,
797 const string &vdef )
798{
799 vector<double> r( f ) ;
800 double factor = v / 2.99792458e8 ;
801 if ( vdef.compare( 0, 3, "RAD" ) == 0 ) {
802 factor = 1.0 / ( 1.0 + factor ) ;
803 r[1] *= factor ;
804 r[2] *= factor ;
805 }
806 else if ( vdef.compare( 0, 3, "OPT" ) == 0 ) {
807 factor += 1.0 ;
808 r[1] *= factor ;
809 r[2] *= factor ;
810 }
811 else {
812 cout << "vdef=" << vdef << " is not supported." << endl;
813 }
814 return r ;
815}
Note: See TracBrowser for help on using the repository browser.