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

Last change on this file since 2872 was 2804, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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