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

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...


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

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