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

Last change on this file since 2763 was 2761, 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...

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