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

Last change on this file since 2965 was 2940, checked in by Malte Marquarding, 10 years ago

Finally fix aall warning: deprecated conversion from string constant to 'char*'

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