source: trunk/external-alma/atnf/PKSIO/NRODataset.cc@ 2218

Last change on this file since 2218 was 2212, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: Yes CAS-2819

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...

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