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

Last change on this file since 3150 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...


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