source: branches/newfiller/external-alma/atnf/PKSIO/NRODataset.cc@ 1800

Last change on this file since 1800 was 1757, checked in by Kana Sugimoto, 15 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.


File size: 21.5 KB
Line 
1//#---------------------------------------------------------------------------
2//# NRODataset.cc: Base class for NRO dataset.
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: 2009/02/27, Takeshi Nakazato, NAOJ
31//#---------------------------------------------------------------------------
32
33#include <atnf/PKSIO/NRODataset.h>
34#include <casa/OS/Time.h>
35#include <scimath/Mathematics/InterpolateArray1D.h>
36
37#include <math.h>
38#include <fstream>
39
40//#include <casa/namespace.h>
41
42using namespace std ;
43
44//
45// NRODataset
46//
47// Base class for NRO dataset.
48//
49
50// constructor
51NRODataset::NRODataset( string name )
52{
53 LogIO os( LogOrigin( "NRODataset", "NRODataset()", WHERE ) ) ;
54
55 // memory allocation
56 initialize() ;
57
58 filename_ = name ;
59 fp_ = NULL ;
60 scanNum_ = 0 ;
61 rowNum_ = 0 ;
62 scanLen_ = 0 ;
63 dataLen_ = 0 ;
64 dataid_ = -1 ;
65
66 // OS endian
67 int i = 1 ;
68 endian_ = -1 ;
69 if ( *reinterpret_cast<char *>(&i) == 1 ) {
70 endian_ = LITTLE_ENDIAN ;
71 os << LogIO::NORMAL << "LITTLE_ENDIAN " << LogIO::POST ;
72 }
73 else {
74 endian_ = BIG_ENDIAN ;
75 os << LogIO::NORMAL << "BIG_ENDIAN " << LogIO::POST ;
76 }
77 same_ = -1 ;
78}
79
80// destructor
81NRODataset::~NRODataset()
82{
83 // release memory
84 releaseRecord() ;
85
86 // close file
87 close() ;
88}
89
90// data initialization
91void NRODataset::initialize()
92{
93 datasize_ = sizeof( char ) * 8 // LOFIL
94 + sizeof( char ) * 8 // VER
95 + sizeof( char ) * 16 // GROUP
96 + sizeof( char ) * 16 // PROJ
97 + sizeof( char ) * 24 // SCHED
98 + sizeof( char ) * 40 // OBSVR
99 + sizeof( char ) * 16 // LOSTM
100 + sizeof( char ) * 16 // LOETM
101 + sizeof( int ) * 2 // ARYNM, NSCAN
102 + sizeof( char ) * 120 // TITLE
103 + sizeof( char ) * 16 // OBJ
104 + sizeof( char ) * 8 // EPOCH
105 + sizeof( double ) * 4 // RA0, DEC0, GLNG0, GLAT0
106 + sizeof( int ) * 2 // NCALB, SCNCD
107 + sizeof( char ) * 120 // SCMOD
108 + sizeof( double ) // URVEL
109 + sizeof( char ) * 4 // VREF
110 + sizeof( char ) * 4 // VDEF
111 + sizeof( char ) * 8 // SWMOD
112 + sizeof( double ) * 8 // FRQSW, DBEAM, MLTOF, CMTQ, CMTE, CMTSOM, CMTNODE, CMTI
113 + sizeof( char ) * 24 // CMTTM
114 + sizeof( double ) * 6 // SBDX, SBDY, SBDZ1, SBDZ2, DAZP, DELP
115 + sizeof( int ) * 4 // CHBIND, NUMCH, CHMIN, CHMAX
116 + sizeof( double ) * 3 // ALCTM, IPTIM, PA
117 + sizeof( int ) * 3 // SCNLEN, SBIND, IBIT
118 + sizeof( char ) * 8 ; // SITE
119
120 // NRODataRecord
121 record_ = new NRODataRecord() ;
122 record_->LDATA = NULL ;
123}
124
125void NRODataset::convertEndian( int &value )
126{
127 char volatile *first = reinterpret_cast<char volatile *>( &value ) ;
128 char volatile *last = first + sizeof( int ) ;
129 std::reverse( first, last ) ;
130}
131
132void NRODataset::convertEndian( float &value )
133{
134 char volatile *first = reinterpret_cast<char volatile *>( &value ) ;
135 char volatile *last = first + sizeof( float ) ;
136 std::reverse( first, last ) ;
137}
138
139void NRODataset::convertEndian( double &value )
140{
141 char volatile *first = reinterpret_cast<char volatile *>( &value ) ;
142 char volatile *last = first + sizeof( double ) ;
143 std::reverse( first, last ) ;
144}
145
146int NRODataset::readHeader( char *v, int size )
147{
148 if ( (int)( fread( v, 1, size, fp_ ) ) != size ) {
149 return -1 ;
150 }
151 return 0 ;
152}
153
154int NRODataset::readHeader( int &v, int b )
155{
156 if ( fread( &v, 1, sizeof(int), fp_ ) != sizeof(int) ) {
157 return -1 ;
158 }
159
160 if ( b == 0 )
161 convertEndian( v ) ;
162
163 return 0 ;
164}
165
166int NRODataset::readHeader( float &v, int b )
167{
168 if ( fread( &v, 1, sizeof(float), fp_ ) != sizeof(float) ) {
169 return -1 ;
170 }
171
172 if ( b == 0 )
173 convertEndian( v ) ;
174
175 return 0 ;
176}
177
178int NRODataset::readHeader( double &v, int b )
179{
180 if ( fread( &v, 1, sizeof(double), fp_ ) != sizeof(double) ) {
181 return -1 ;
182 }
183
184 if ( b == 0 )
185 convertEndian( v ) ;
186
187 return 0 ;
188}
189
190void NRODataset::convertEndian( NRODataRecord *r )
191{
192 convertEndian( r->ISCAN ) ;
193 convertEndian( r->DSCX ) ;
194 convertEndian( r->DSCY ) ;
195 convertEndian( r->SCX ) ;
196 convertEndian( r->SCY ) ;
197 convertEndian( r->PAZ ) ;
198 convertEndian( r->PEL ) ;
199 convertEndian( r->RAZ ) ;
200 convertEndian( r->REL ) ;
201 convertEndian( r->XX ) ;
202 convertEndian( r->YY ) ;
203 convertEndian( r->TEMP ) ;
204 convertEndian( r->PATM ) ;
205 convertEndian( r->PH2O ) ;
206 convertEndian( r->VWIND ) ;
207 convertEndian( r->DWIND ) ;
208 convertEndian( r->TAU ) ;
209 convertEndian( r->TSYS ) ;
210 convertEndian( r->BATM ) ;
211 convertEndian( r->LINE ) ;
212 for ( int i = 0 ; i < 4 ; i++ )
213 convertEndian( r->IDMY1[i] ) ;
214 convertEndian( r->VRAD ) ;
215 convertEndian( r->FREQ0 ) ;
216 convertEndian( r->FQTRK ) ;
217 convertEndian( r->FQIF1 ) ;
218 convertEndian( r->ALCV ) ;
219 for ( int i = 0 ; i < 2 ; i++ )
220 for ( int j = 0 ; j < 2 ; j++ )
221 convertEndian( r->OFFCD[i][j] ) ;
222 convertEndian( r->IDMY0 ) ;
223 convertEndian( r->IDMY2 ) ;
224 convertEndian( r->DPFRQ ) ;
225 convertEndian( r->SFCTR ) ;
226 convertEndian( r->ADOFF ) ;
227}
228
229void NRODataset::releaseRecord()
230{
231 if ( record_ ) {
232 if ( record_->LDATA != NULL ) {
233 delete record_->LDATA ;
234 record_->LDATA = NULL ;
235 }
236 delete record_ ;
237 record_ = NULL ;
238 }
239 dataid_ = -1 ;
240}
241
242// Get specified scan
243NRODataRecord *NRODataset::getRecord( int i )
244{
245 LogIO os( LogOrigin( "NRODataset", "getRecord()", WHERE ) ) ;
246
247 // DEBUG
248 //cout << "NRODataset::getRecord() Start " << i << endl ;
249 //
250 if ( i < 0 || i >= rowNum_ ) {
251 //cerr << "NRODataset::getRecord() data index out of range." << endl ;
252 os << LogIO::SEVERE << "data index " << i << " out of range. return NULL." << LogIO::POST ;
253 return NULL ;
254 }
255
256 if ( i == dataid_ ) {
257 return record_ ;
258 }
259
260 // DEBUG
261 //cout << "NRODataset::getData() Get new dataset" << endl ;
262 //
263 // read data
264 int status = fillRecord( i ) ;
265 if ( status == 0 ) {
266 dataid_ = i ;
267 }
268 else {
269 //cerr << "NRODataset::getRecord() error while reading data " << i << endl ;
270 os << LogIO::SEVERE << "error while reading data " << i << ". return NULL." << LogIO::POST ;
271 dataid_ = -1 ;
272 return NULL ;
273 }
274
275 return record_ ;
276}
277
278int NRODataset::fillRecord( int i )
279{
280 LogIO os( LogOrigin( "NRODataset", "fillRecord()", WHERE ) ) ;
281
282 int status = 0 ;
283
284 status = open() ;
285 if ( status != 0 )
286 return status ;
287
288
289 // fill NRODataset
290 int offset = getDataSize() + scanLen_ * i ;
291 // DEBUG
292 //cout << "NRODataset::fillRecord() offset (header) = " << offset << endl ;
293 //cout << "NRODataset::fillRecord() sizeof(NRODataRecord) = " << sizeof( NRODataRecord ) << " byte" << endl ;
294 fseek( fp_, offset, SEEK_SET ) ;
295 if ( (int)fread( record_, 1, SCAN_HEADER_SIZE, fp_ ) != SCAN_HEADER_SIZE ) {
296 //cerr << "Failed to read scan header: " << i << endl ;
297 os << LogIO::SEVERE << "Failed to read scan header for " << i << "th row." << LogIO::POST ;
298 return -1 ;
299 }
300 if ( (int)fread( record_->LDATA, 1, dataLen_, fp_ ) != dataLen_ ) {
301 //cerr << "Failed to read spectral data: " << i << endl ;
302 os << LogIO::SEVERE << "Failed to read spectral data for " << i << "th row." << LogIO::POST ;
303 return -1 ;
304 }
305
306 if ( same_ == 0 ) {
307 convertEndian( record_ ) ;
308 }
309
310 // DWIND unit conversion (deg -> rad)
311 record_->DWIND = record_->DWIND * M_PI / 180.0 ;
312
313 return status ;
314}
315
316// open
317int NRODataset::open()
318{
319 int status = 0 ;
320
321 if ( fp_ == NULL ) {
322 if ( (fp_ = fopen( filename_.c_str(), "rb" )) == NULL )
323 status = -1 ;
324 else
325 status = 0 ;
326 }
327
328 return status ;
329}
330
331// close
332void NRODataset::close()
333{
334 // DEBUG
335 //cout << "NRODataset::close() close file" << endl ;
336 //
337 if ( fp_ != NULL )
338 fclose( fp_ ) ;
339 fp_ = NULL ;
340}
341
342// get spectrum
343vector< vector<double> > NRODataset::getSpectrum()
344{
345 vector< vector<double> > spec;
346
347 for ( int i = 0 ; i < rowNum_ ; i++ ) {
348 spec.push_back( getSpectrum( i ) ) ;
349 }
350
351 return spec ;
352}
353
354vector<double> NRODataset::getSpectrum( int i )
355{
356 LogIO os( LogOrigin( "NRODataset", "getSpectrum", WHERE ) ) ;
357
358 // DEBUG
359 //cout << "NRODataset::getSpectrum() start process (" << i << ")" << endl ;
360 //
361 // size of spectrum is not chmax_ but dataset_->getNCH() after binding
362 int nchan = NUMCH ;
363 vector<double> spec( chmax_, 0.0 ) ; // spectrum "before" binding
364 vector<double> bspec( nchan, 0.0 ) ; // spectrum "after" binding
365 // DEBUG
366 //cout << "NRODataset::getSpectrum() nchan = " << nchan << " chmax_ = " << chmax_ << endl ;
367 //
368
369 NRODataRecord *record = getRecord( i ) ;
370
371 int bit = IBIT ; // fixed to 12 bit
372 double scale = record->SFCTR ;
373 // DEBUG
374 //cout << "NRODataset::getSpectrum() scale = " << scale << endl ;
375 //
376 double offset = record->ADOFF ;
377 // DEBUG
378 //cout << "NRODataset::getSpectrum() offset = " << offset << endl ;
379 //
380 if ( ( scale == 0.0 ) && ( offset == 0.0 ) ) {
381 //cerr << "NRODataset::getSpectrum() zero spectrum (" << i << ")" << endl ;
382 return bspec ;
383 }
384 char *cdata = record->LDATA ;
385 vector<double> mscale = MLTSCF ;
386 double dscale = mscale[getIndex( i )] ;
387 int cbind = CHBIND ;
388 int chmin = CHMIN ;
389
390 // char -> int
391 vector<int> ispec( chmax_, 0 ) ;
392 union SharedMemory {
393 int ivalue ;
394 unsigned char cbuf[4] ;
395 } ;
396 SharedMemory u ;
397 int j = 0 ;
398 char ctmp = 0x00 ;
399 int sw = 0 ;
400 for ( int i = 0 ; i < chmax_ ; i++ ) {
401 if ( bit == 12 ) { // 12 bit qunatization
402 u.ivalue = 0 ;
403
404 if ( endian_ == BIG_ENDIAN ) {
405 // big endian
406 if ( sw == 0 ) {
407 char c0 = (cdata[j] >> 4) & 0x0f ;
408 char c1 = ((cdata[j] << 4) & 0xf0) | ((cdata[j+1] >> 4) & 0x0f) ;
409 ctmp = cdata[j+1] & 0x0f ;
410 u.cbuf[2] = c0 ;
411 u.cbuf[3] = c1 ;
412 j += 2 ;
413 sw = 1 ;
414 }
415 else if ( sw == 1 ) {
416 u.cbuf[2] = ctmp ;
417 u.cbuf[3] = cdata[j] ;
418 j++ ;
419 sw = 0 ;
420 }
421 }
422 else if ( endian_ == LITTLE_ENDIAN ) {
423 // little endian
424 if ( sw == 0 ) {
425 char c0 = (cdata[j] >> 4) & 0x0f ;
426 char c1 = ((cdata[j] << 4) & 0xf0) | ((cdata[j+1] >> 4) & 0x0f) ;
427 ctmp = cdata[j+1] & 0x0f ;
428 u.cbuf[1] = c0 ;
429 u.cbuf[0] = c1 ;
430 j += 2 ;
431 sw = 1 ;
432 }
433 else if ( sw == 1 ) {
434 u.cbuf[1] = ctmp ;
435 u.cbuf[0] = cdata[j] ;
436 j++ ;
437 sw = 0 ;
438 }
439 }
440 }
441
442 ispec[i] = u.ivalue ;
443 if ( ( ispec[i] < 0 ) || ( ispec[i] > 4096 ) ) {
444 //cerr << "NRODataset::getSpectrum() ispec[" << i << "] is out of range" << endl ;
445 os << LogIO::SEVERE << "ispec[" << i << "] is out of range" << LogIO::EXCEPTION ;
446 return bspec ;
447 }
448 // DEBUG
449 //cout << "NRODataset::getSpectrum() ispec[" << i << "] = " << ispec[i] << endl ;
450 //
451 }
452
453 // int -> double
454 for ( int i = 0 ; i < chmax_ ; i++ ) {
455 spec[i] = (double)( ispec[i] * scale + offset ) * dscale ;
456 // DEBUG
457 //cout << "NRODataset::getSpectrum() spec[" << i << "] = " << spec[i] << endl ;
458 //
459 }
460
461 // channel binding
462 if ( cbind != 1 ) {
463 int k = chmin ;
464 double sum0 = 0 ;
465 double sum1 = 0 ;
466 for ( int i = 0 ; i < nchan ; i++ ) {
467 for ( int j = k ; j < k + cbind ; j++ ) {
468 sum0 += spec[k] ;
469 sum1++ ;
470 }
471 bspec[i] = sum0 / sum1 ;
472 k += cbind ;
473 // DEBUG
474 //cout << "NRODataset::getSpectrum() bspec[" << i << "] = " << bspec[i] << endl ;
475 //
476 }
477 }
478 else {
479 for ( int i = 0 ; i < nchan ; i++ )
480 bspec[i] = spec[i] ;
481 }
482
483 // DEBUG
484 //cout << "NRODataset::getSpectrum() end process" << endl ;
485 //
486
487 return bspec ;
488}
489
490int NRODataset::getIndex( int irow )
491{
492 // DEBUG
493 //cout << "NRODataset::getIndex() start" << endl ;
494 //
495 NRODataRecord *record = getRecord( irow ) ;
496 string str = record->ARRYT ;
497 // DEBUG
498 //cout << "NRODataset::getIndex() str = " << str << endl ;
499 //
500 string substr = str.substr( 1, 2 ) ;
501 unsigned int index = (unsigned int)(atoi( substr.c_str() ) - 1) ;
502 // DEBUG
503 //cout << "NRODataset::getIndex() irow = " << irow << " index = " << index << endl ;
504 //
505
506 // DEBUG
507 //cout << "NRODataset::getIndex() end" << endl ;
508 //
509 return index ;
510}
511
512int NRODataset::getPolarizationNum()
513{
514 // DEBUG
515 //cout << "NRODataset::getPolarizationNum() start process" << endl ;
516 //
517 int npol = 0 ;
518
519 vector<string> type( 2 ) ;
520 type[0] = "CIRC" ;
521 type[1] = "LINR" ;
522 vector<double> crot ;
523 vector<double> lagl ;
524 //vector<int> ntype( 2, 0 ) ;
525
526 unsigned int imax = rowNum_ ;
527 for ( unsigned int i = 0 ; i < imax ; i++ ) {
528 int index = getIndex( i ) ;
529 // DEBUG
530 //cout <<"NRODataset::getPolarizationNum() index = " << index << endl ;
531 //
532 if ( POLTP[index] == type[0] ) {
533 if( count( crot.begin(), crot.end(), POLDR[index] ) != 0 ) {
534 crot.push_back( POLDR[index] ) ;
535 npol++ ;
536 }
537 //ntype[0] = 1 ;
538 }
539 else if ( POLTP[index] == type[1] ) {
540 if ( count( lagl.begin(), lagl.end(), POLAN[index] ) != 0 ) {
541 lagl.push_back( POLAN[index] ) ;
542 npol++ ;
543 }
544 //ntype[1] = 1 ;
545 }
546 }
547
548 if ( npol == 0 )
549 npol = 1 ;
550
551 // DEBUG
552 //cout << "NRODataset::getPolarizationNum() end process" << endl ;
553 //
554
555 return npol ;
556}
557
558vector<double> NRODataset::getStartIntTime()
559{
560 vector<double> times ;
561 for ( int i = 0 ; i < rowNum_ ; i++ ) {
562 times.push_back( getStartIntTime( i ) ) ;
563 }
564 return times ;
565}
566
567double NRODataset::getStartIntTime( int i )
568{
569 NRODataRecord *record = getRecord( i ) ;
570
571 char *t = record->LAVST ;
572 return getMJD( t ) ;
573}
574
575double NRODataset::getMJD( char *time )
576{
577 // TODO: should be checked which time zone the time depends on
578 // 2008/11/14 Takeshi Nakazato
579 string strStartTime( time ) ;
580 string strYear = strStartTime.substr( 0, 4 ) ;
581 string strMonth = strStartTime.substr( 4, 2 ) ;
582 string strDay = strStartTime.substr( 6, 2 ) ;
583 string strHour = strStartTime.substr( 8, 2 ) ;
584 string strMinute = strStartTime.substr( 10, 2 ) ;
585 string strSecond = strStartTime.substr( 12, strStartTime.size() - 12 ) ;
586 unsigned int year = atoi( strYear.c_str() ) ;
587 unsigned int month = atoi( strMonth.c_str() ) ;
588 unsigned int day = atoi( strDay.c_str() ) ;
589 unsigned int hour = atoi( strHour.c_str() ) ;
590 unsigned int minute = atoi( strMinute.c_str() ) ;
591 double second = atof( strSecond.c_str() ) ;
592 Time t( year, month, day, hour, minute, second ) ;
593
594 return t.modifiedJulianDay() ;
595}
596
597vector<bool> NRODataset::getIFs()
598{
599 vector<bool> v ;
600 vector< vector<double> > fref ;
601 vector< vector<double> > chcal = CHCAL ;
602 vector<double> f0cal = F0CAL ;
603 vector<double> beres = BERES ;
604 for ( int i = 0 ; i < rowNum_ ; i++ ) {
605 vector<double> f( 4, 0 ) ;
606 uInt index = getIndex( i ) ;
607 f[0] = chcal[index][0] ;
608 f[1] = f0cal[index] ;
609 f[2] = beres[index] ;
610 if ( f[0] != 0. ) {
611 f[1] = f[1] - f[0] * f[2] ;
612 }
613 NRODataRecord *record = getRecord( i ) ;
614 f[3] = record->FREQ0 ;
615 if ( v.size() == 0 ) {
616 v.push_back( True ) ;
617 fref.push_back( f ) ;
618 }
619 else {
620 bool b = true ;
621 int fsize = fref.size() ;
622 for ( int j = 0 ; j < fsize ; j++ ) {
623 if ( fref[j][1] == f[1] && fref[j][2] == f[2] && fref[j][3] == f[3] ) {
624 b = false ;
625 }
626 }
627 if ( b ) {
628 v.push_back( True ) ;
629 fref.push_back( f ) ;
630 }
631 }
632 }
633
634
635 // DEBUG
636 //cout << "NRODataset::getIFs() number of IF is " << v.size() << endl ;
637 //
638
639 return v ;
640}
641
642vector<double> NRODataset::getFrequencies( int i )
643{
644 // return value
645 // v[0] reference channel
646 // v[1] reference frequency
647 // v[2] frequency increment
648 vector<double> v( 3, 0.0 ) ;
649
650 NRODataRecord *record = getRecord( i ) ;
651 string arryt = string( record->ARRYT ) ;
652 //string sbeamno = arryt.substr( 1, arryt.size()-1 ) ;
653 //uInt ib = atoi( sbeamno.c_str() ) - 1 ;
654 uInt ib = getArrayId( arryt ) ;
655
656 int ia = -1 ;
657 bool isAOS = false ;
658 //cout << "NRODataset::getFrequencies() record->ARRYT=" << record->ARRYT << endl ;
659 //cout << "NRODataset::getFrequencies() ib = " << ib << endl ;
660
661 if ( strncmp( record->ARRYT, "X", 1 ) == 0 ) {
662 // FX
663 if ( strncmp( (record->ARRYT)+1, "1", 1 ) == 0
664 || strncmp( (record->ARRYT)+1, "3", 1 ) ) {
665 // FX1, 3
666 ia = 2 ;
667 }
668 else {
669 // FX2, 4
670 ia = 1 ;
671 }
672 }
673 else if ( strncmp( record->ARRYT, "A", 1 ) == 0 )
674 ia = 2 ; // AC
675 else if ( strncmp( record->ARRYT, "W", 1 ) == 0 ) {
676 // AOS-W
677 ia = 2 ;
678 isAOS = true ;
679 }
680 else if ( strncmp( record->ARRYT, "U", 1 ) == 0 ) {
681 // AOS-U
682 ia = 2 ;
683 isAOS = true ;
684 }
685 else if ( strncmp( record->ARRYT, "H", 1 ) == 0 ) {
686 // AOS-H
687 isAOS = true ;
688 //cout << record->ARRYT << " " << strlen(record->ARRYT) << endl ;
689 //cout << (record->ARRYT)+1 << endl ;
690 if ( strncmp( (record->ARRYT)+2, " ", 1 ) == 0 ) {
691 // H1-9
692 if ( strncmp( (record->ARRYT)+1, "9", 1 ) == 0 ) {
693 // H9
694 ia = 2 ;
695 }
696 else {
697 // H1-8
698 ia = 1 ;
699 }
700 }
701 else {
702 // H10-16
703 ia = 2 ;
704 }
705 }
706
707 int iu ;
708 if ( record->FQIF1 > 0 )
709 iu = 1 ; // USB
710 else
711 iu = 2 ; // LSB
712
713 int ivdef = -1 ;
714 //if ( strncmp( (dataset_->getVDEF()).c_str(), "RAD", 3 ) == 0 )
715 if ( (getVDEF()).compare( 0, 3, "RAD" ) == 0 )
716 ivdef = 0 ;
717 //else if ( strncmp( dataset_->getVDEF(), "OPT", 3 ) == 0 )
718 else if ( (getVDEF()).compare( 0, 3, "OPT" ) == 0 )
719 ivdef = 1 ;
720 // DEBUG
721 //cout << "NRODataset::getFrequencies() ivdef = " << ivdef << endl ;
722 //
723 double vel = getURVEL() + record->VRAD ;
724 double cvel = 2.99792458e8 ; // speed of light [m/s]
725 double fq0 = record->FREQ0 ;
726 //double fq0 = record->FQTRK ;
727
728 int ncal = getNFCAL()[ib] ;
729 vector<double> freqs( ncal ) ;
730 double cw = 0.0 ;
731 vector<double> fqcal = getFQCAL()[ib] ;
732 vector<double> chcal = getCHCAL()[ib] ;
733
734 for ( int ii = 0 ; ii < ncal ; ii++ ) {
735 freqs[ii] = fqcal[ii] ;
736 freqs[ii] -= getF0CAL()[ib] ;
737 if ( ia == 1 ) {
738 if ( iu == 1 ) {
739 freqs[ii] = fq0 + freqs[ii] ;
740 }
741 else if ( iu == 2 ) {
742 freqs[ii] = fq0 - freqs[ii] ;
743 }
744 }
745 else if ( ia == 2 ) {
746 if ( iu == 1 ) {
747 freqs[ii] = fq0 - freqs[ii] ;
748 }
749 else if ( iu == 2 ) {
750 freqs[ii] = fq0 + freqs[ii] ;
751 }
752 }
753// if ( ivdef == 0 ) {
754// double factor = 1.0 / ( 1.0 - vel / cvel ) ;
755// freqs[ii] = freqs[ii] * factor - record->FQTRK * ( factor - 1.0 ) ;
756// }
757// else if ( ivdef == 1 ) {
758// double factor = vel / cvel ;
759// freqs[ii] = freqs[ii] * ( 1.0 + factor ) - record->FQTRK * factor ;
760// }
761 }
762
763 if ( isAOS ) {
764 // regridding
765 while ( ncal < (int)chcal.size() ) {
766 chcal.pop_back() ;
767 }
768 Vector<Double> xin( chcal ) ;
769 Vector<Double> yin( freqs ) ;
770 int nchan = getNUMCH() ;
771 Vector<Double> xout( nchan ) ;
772 indgen( xout ) ;
773 Vector<Double> yout ;
774 InterpolateArray1D<Double, Double>::interpolate( yout, xout, xin, yin, InterpolateArray1D<Double,Double>::cubic ) ;
775 Double bw = abs( yout[nchan-1] - yout[0] ) ;
776 bw += 0.5 * abs( yout[nchan-1] - yout[nchan-2] + yout[1] - yout[0] ) ;
777 Double dz = bw / (Double) nchan ;
778 if ( yout[0] > yout[1] )
779 dz = - dz ;
780 v[0] = 0 ;
781 v[1] = yout[0] ;
782 v[2] = dz ;
783 }
784 else {
785 cw = getBERES()[ib] ;
786
787 if ( cw == 0.0 ) {
788 cw = ( freqs[1] - freqs[0] )
789 / ( chcal[1] - chcal[0] ) ;
790// if ( cw < 0.0 )
791// cw = - cw ;
792 }
793 v[0] = chcal[0] - 1 ; // 0-base
794 v[1] = freqs[0] ;
795 v[2] = cw ;
796 }
797
798 return v ;
799}
800
801uInt NRODataset::getArrayId( string type )
802{
803 string sbeamno = type.substr( 1, type.size()-1 ) ;
804 uInt ib = atoi( sbeamno.c_str() ) - 1 ;
805 return ib ;
806}
807
808void NRODataset::show()
809{
810 LogIO os( LogOrigin( "NRODataset", "show()", WHERE ) ) ;
811
812 os << LogIO::NORMAL << "------------------------------------------------------------" << endl ;
813 os << LogIO::NORMAL << "Number of scan = " << scanNum_ << endl ;
814 os << LogIO::NORMAL << "Number of data record = " << rowNum_ << endl ;
815 os << LogIO::NORMAL << "Length of data record = " << scanLen_ << " bytes" << endl ;
816 os << LogIO::NORMAL << "Allocated memory for spectral data = " << dataLen_ << " bytes" << endl ;
817 os << LogIO::NORMAL << "Max number of channel = " << chmax_ << endl ;
818 os << LogIO::NORMAL << "------------------------------------------------------------" << endl ;
819 os.post() ;
820
821}
Note: See TracBrowser for help on using the repository browser.