source: branches/parallelCasa3.3/external-alma/atnf/PKSIO/NRODataset.cc@ 2391

Last change on this file since 2391 was 2266, checked in by KohjiNakamura, 13 years ago

optimize NRODataset::getSpectrum()
File size: 22.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 // endian matches
74 same_ = -1 ;
75
76 // Record for frequency setting
77 frec_ = Record() ;
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(rowNum_);
346
347 for ( int i = 0 ; i < rowNum_ ; i++ ) {
348 spec[i] = 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 const int nchan = NUMCH ;
363 vector<double> bspec( nchan, 0.0 ) ; // spectrum "after" binding
364 // DEBUG
365 //cout << "NRODataset::getSpectrum() nchan = " << nchan << " chmax_ = " << chmax_ << endl ;
366 //
367
368 NRODataRecord *record = getRecord( i ) ;
369
370 const int bit = IBIT ; // fixed to 12 bit
371 double scale = record->SFCTR ;
372 // DEBUG
373 //cout << "NRODataset::getSpectrum() scale = " << scale << endl ;
374 //
375 double offset = record->ADOFF ;
376 // DEBUG
377 //cout << "NRODataset::getSpectrum() offset = " << offset << endl ;
378 //
379 if ( ( scale == 0.0 ) && ( offset == 0.0 ) ) {
380 //cerr << "NRODataset::getSpectrum() zero spectrum (" << i << ")" << endl ;
381 return bspec ;
382 }
383 unsigned char *cdata = (unsigned char *)record->LDATA ;
384 vector<double> mscale = MLTSCF ;
385 double dscale = mscale[getIndex( i )] ;
386 int cbind = CHBIND ;
387 int chmin = CHMIN ;
388
389 // char -> int
390 vector<int> ispec( chmax_, 0 ) ;
391
392 static const int shift_right[] = {
393 4, 0
394 };
395 static const int start_pos[] = {
396 0, 1
397 };
398 static const int incr[] = {
399 0, 3
400 };
401 int j = 0 ;
402 for ( int i = 0 ; i < chmax_ ; i++ ) {
403 int ivalue = 0 ;
404 if ( bit == 12 ) { // 12 bit qunatization
405 const int idx = j + start_pos[i & 1];
406 const unsigned tmp = unsigned(cdata[idx]) << 8 | cdata[idx + 1];
407 ivalue = int((tmp >> shift_right[i & 1]) & 0xFFF);
408 j += incr[i & 1];
409 }
410
411 ispec[i] = ivalue ;
412 if ( ( ispec[i] < 0 ) || ( ispec[i] > 4096 ) ) {
413 //cerr << "NRODataset::getSpectrum() ispec[" << i << "] is out of range" << endl ;
414 os << LogIO::SEVERE << "ispec[" << i << "] is out of range" << LogIO::EXCEPTION ;
415 return bspec ;
416 }
417 // DEBUG
418 //cout << "NRODataset::getSpectrum() ispec[" << i << "] = " << ispec[i] << endl ;
419 //
420 }
421
422 double *const spec = new double[ chmax_ ] ; // spectrum "before" binding
423 // int -> double
424 for ( int i = 0 ; i < chmax_ ; i++ ) {
425 spec[i] = (double)( ispec[i] * scale + offset ) * dscale ;
426 // DEBUG
427 //cout << "NRODataset::getSpectrum() spec[" << i << "] = " << spec[i] << endl ;
428 //
429 }
430
431 // channel binding
432 if ( cbind != 1 ) {
433 int k = chmin ;
434 double sum0 = 0 ;
435 double sum1 = 0 ;
436 for ( int i = 0 ; i < nchan ; i++ ) {
437 for ( int j = k ; j < k + cbind ; j++ ) {
438 sum0 += spec[k] ;
439 sum1++ ;
440 }
441 bspec[i] = sum0 / sum1 ;
442 k += cbind ;
443 // DEBUG
444 //cout << "NRODataset::getSpectrum() bspec[" << i << "] = " << bspec[i] << endl ;
445 //
446 }
447 }
448 else {
449 for ( int i = 0 ; i < nchan ; i++ )
450 bspec[i] = spec[i] ;
451 }
452 delete[] spec;
453
454 // DEBUG
455 //cout << "NRODataset::getSpectrum() end process" << endl ;
456 //
457
458 return bspec ;
459}
460
461int NRODataset::getIndex( int irow )
462{
463 // DEBUG
464 //cout << "NRODataset::getIndex() start" << endl ;
465 //
466 NRODataRecord *record = getRecord( irow ) ;
467 string str = record->ARRYT ;
468 // DEBUG
469 //cout << "NRODataset::getIndex() str = " << str << endl ;
470 //
471 string substr = str.substr( 1, 2 ) ;
472 unsigned int index = (unsigned int)(atoi( substr.c_str() ) - 1) ;
473 // DEBUG
474 //cout << "NRODataset::getIndex() irow = " << irow << " index = " << index << endl ;
475 //
476
477 // DEBUG
478 //cout << "NRODataset::getIndex() end" << endl ;
479 //
480 return index ;
481}
482
483int NRODataset::getPolarizationNum()
484{
485 // DEBUG
486 //cout << "NRODataset::getPolarizationNum() start process" << endl ;
487 //
488 int npol = 0 ;
489
490 vector<string> type( 2 ) ;
491 type[0] = "CIRC" ;
492 type[1] = "LINR" ;
493 vector<double> crot ;
494 vector<double> lagl ;
495 //vector<int> ntype( 2, 0 ) ;
496
497 unsigned int imax = rowNum_ ;
498 for ( unsigned int i = 0 ; i < imax ; i++ ) {
499 int index = getIndex( i ) ;
500 // DEBUG
501 //cout <<"NRODataset::getPolarizationNum() index = " << index << endl ;
502 //
503 if ( POLTP[index] == type[0] ) {
504 if( count( crot.begin(), crot.end(), POLDR[index] ) != 0 ) {
505 crot.push_back( POLDR[index] ) ;
506 npol++ ;
507 }
508 //ntype[0] = 1 ;
509 }
510 else if ( POLTP[index] == type[1] ) {
511 if ( count( lagl.begin(), lagl.end(), POLAN[index] ) != 0 ) {
512 lagl.push_back( POLAN[index] ) ;
513 npol++ ;
514 }
515 //ntype[1] = 1 ;
516 }
517 }
518
519 if ( npol == 0 )
520 npol = 1 ;
521
522 // DEBUG
523 //cout << "NRODataset::getPolarizationNum() end process" << endl ;
524 //
525
526 return npol ;
527}
528
529vector<double> NRODataset::getStartIntTime()
530{
531 vector<double> times ;
532 for ( int i = 0 ; i < rowNum_ ; i++ ) {
533 times.push_back( getStartIntTime( i ) ) ;
534 }
535 return times ;
536}
537
538double NRODataset::getStartIntTime( int i )
539{
540 NRODataRecord *record = getRecord( i ) ;
541
542 char *t = record->LAVST ;
543 return getMJD( t ) ;
544}
545
546double NRODataset::getMJD( char *time )
547{
548 // TODO: should be checked which time zone the time depends on
549 // 2008/11/14 Takeshi Nakazato
550 string strStartTime( time ) ;
551 string strYear = strStartTime.substr( 0, 4 ) ;
552 string strMonth = strStartTime.substr( 4, 2 ) ;
553 string strDay = strStartTime.substr( 6, 2 ) ;
554 string strHour = strStartTime.substr( 8, 2 ) ;
555 string strMinute = strStartTime.substr( 10, 2 ) ;
556 string strSecond = strStartTime.substr( 12, strStartTime.size() - 12 ) ;
557 unsigned int year = atoi( strYear.c_str() ) ;
558 unsigned int month = atoi( strMonth.c_str() ) ;
559 unsigned int day = atoi( strDay.c_str() ) ;
560 unsigned int hour = atoi( strHour.c_str() ) ;
561 unsigned int minute = atoi( strMinute.c_str() ) ;
562 double second = atof( strSecond.c_str() ) ;
563 Time t( year, month, day, hour, minute, second ) ;
564
565 return t.modifiedJulianDay() ;
566}
567
568double NRODataset::getScanTime( int i )
569{
570 double startTime = getStartIntTime( i ) ;
571 double interval = getIPTIM() / 86400.0 ; // [sec]->[day]
572 return startTime+0.5*interval ;
573}
574
575vector<bool> NRODataset::getIFs()
576{
577 vector<bool> v ;
578 vector< vector<double> > fref ;
579 vector< vector<double> > chcal = CHCAL ;
580 vector<double> f0cal = F0CAL ;
581 vector<double> beres = BERES ;
582 for ( int i = 0 ; i < rowNum_ ; i++ ) {
583 vector<double> f( 4, 0 ) ;
584 uInt index = getIndex( i ) ;
585 f[0] = chcal[index][0] ;
586 f[1] = f0cal[index] ;
587 f[2] = beres[index] ;
588 if ( f[0] != 0. ) {
589 f[1] = f[1] - f[0] * f[2] ;
590 }
591 NRODataRecord *record = getRecord( i ) ;
592 f[3] = record->FREQ0 ;
593 if ( v.size() == 0 ) {
594 v.push_back( True ) ;
595 fref.push_back( f ) ;
596 }
597 else {
598 bool b = true ;
599 int fsize = fref.size() ;
600 for ( int j = 0 ; j < fsize ; j++ ) {
601 if ( fref[j][1] == f[1] && fref[j][2] == f[2] && fref[j][3] == f[3] ) {
602 b = false ;
603 }
604 }
605 if ( b ) {
606 v.push_back( True ) ;
607 fref.push_back( f ) ;
608 }
609 }
610 }
611
612
613 // DEBUG
614 //cout << "NRODataset::getIFs() number of IF is " << v.size() << endl ;
615 //
616
617 return v ;
618}
619
620vector<double> NRODataset::getFrequencies( int i )
621{
622 // return value
623 // v[0] reference channel
624 // v[1] reference frequency
625 // v[2] frequency increment
626 vector<double> v( 3, 0.0 ) ;
627
628 NRODataRecord *record = getRecord( i ) ;
629 string arryt = string( record->ARRYT ) ;
630 uInt ib = getArrayId( arryt ) ;
631 string rxname = getRX()[0] ;
632 string key = arryt ;
633 if ( rxname.find("MULT2") != string::npos )
634 key = "BEARS" ;
635
636 if ( frec_.isDefined( key ) ) {
637 // frequency for the array is already calculated
638 Vector<Double> f = frec_.asArrayDouble( key ) ;
639 Double *f_p = f.data() ;
640 for ( int i = 0 ; i < 3 ; i++ )
641 v[i] = (double)(f_p[i]) ;
642 return v ;
643 }
644
645 //int ia = -1 ;
646 bool isAOS = false ;
647 //cout << "NRODataset::getFrequencies() record->ARRYT=" << record->ARRYT << endl ;
648 //cout << "NRODataset::getFrequencies() ib = " << ib << endl ;
649
650 if ( arryt[0] == 'W' || arryt[0] == 'U' || arryt[0] == 'H' )
651 isAOS = true ;
652
653 Bool isUSB ;
654 if ( record->FQIF1 > 0 )
655 isUSB = True ; // USB
656 else
657 isUSB = False ; // LSB
658
659 int ivdef = -1 ;
660 if ( (getVDEF()).compare( 0, 3, "RAD" ) == 0 )
661 ivdef = 0 ;
662 else if ( (getVDEF()).compare( 0, 3, "OPT" ) == 0 )
663 ivdef = 1 ;
664 // DEBUG
665 //cout << "NRODataset::getFrequencies() ivdef = " << ivdef << endl ;
666 //
667 double vel = getURVEL() + record->VRAD ;
668 double cvel = 2.99792458e8 ; // speed of light [m/s]
669 double fq0 = record->FREQ0 ;
670 //double fq0 = record->FQTRK ;
671
672 int ncal = getNFCAL()[ib] ;
673 double cw = 0.0 ;
674 vector<double> fqcal = getFQCAL()[ib] ;
675 vector<double> chcal = getCHCAL()[ib] ;
676 double f0cal = getF0CAL()[ib] ;
677 Vector<Double> freqs( ncal, fq0-f0cal ) ;
678
679 double factor = vel / cvel ;
680 if ( ivdef == 0 )
681 factor = 1.0 / ( 1.0 - factor ) ;
682 for ( int ii = 0 ; ii < ncal ; ii++ ) {
683 freqs[ii] += fqcal[ii] ;
684 if ( ivdef == 0 ) {
685 freqs[ii] = freqs[ii] * factor + record->FQTRK * ( 1.0 - factor ) ;
686 }
687 else if ( ivdef == 1 ) {
688 freqs[ii] = freqs[ii] * ( 1.0 + factor ) - record->FQTRK * factor ;
689 }
690
691 //ofstream ofs("freqs.txt",ios::out|ios::app) ;
692 //ofs << setprecision(16) ;
693 //ofs << i << " " << record->ARRYT << " " << chcal[ii] << " " << freqs[ii] << " " << record->ISCAN << " " << fqcal[ii] << " " << f0cal << " " << record->FQTRK << " " << vel << endl ;
694 //ofs.close() ;
695
696 }
697
698 if ( isAOS ) {
699 // regridding
700 while ( ncal < (int)chcal.size() ) {
701 chcal.pop_back() ;
702 }
703 Vector<Double> xin( chcal ) ;
704 Vector<Double> yin( freqs ) ;
705 int nchan = getNUMCH() ;
706 Vector<Double> xout( nchan ) ;
707 indgen( xout ) ;
708 Vector<Double> yout ;
709 InterpolateArray1D<Double, Double>::interpolate( yout, xout, xin, yin, InterpolateArray1D<Double,Double>::cubic ) ;
710 Double bw = abs( yout[nchan-1] - yout[0] ) ;
711 bw += 0.5 * abs( yout[nchan-1] - yout[nchan-2] + yout[1] - yout[0] ) ;
712 Double dz = bw / (Double) nchan ;
713 if ( yout[0] > yout[1] )
714 dz = - dz ;
715 v[0] = 0 ;
716 v[1] = yout[0] ;
717 v[2] = dz ;
718 }
719 else {
720
721 cw = ( freqs[1] - freqs[0] )
722 / ( chcal[1] - chcal[0] ) ;
723
724 if ( isUSB ) {
725 // channel frequency inversion
726 cw *= -1.0 ;
727 Double tmp = freqs[1] ;
728 freqs[1] = freqs[0] ;
729 freqs[0] = tmp ;
730 }
731
732 v[0] = chcal[0] - 1 ; // 0-base
733 v[1] = freqs[0] ;
734 v[2] = cw ;
735 }
736
737 if ( refFreq_[ib] == 0.0 )
738 refFreq_[ib] = v[1] ;
739
740 // register frequency setting to Record
741 Vector<Double> f( v ) ;
742 frec_.define( key, f ) ;
743
744 return v ;
745}
746
747uInt NRODataset::getArrayId( string type )
748{
749 string sbeamno = type.substr( 1, type.size()-1 ) ;
750 uInt ib = atoi( sbeamno.c_str() ) - 1 ;
751 return ib ;
752}
753
754void NRODataset::show()
755{
756 LogIO os( LogOrigin( "NRODataset", "show()", WHERE ) ) ;
757
758 os << LogIO::NORMAL << "------------------------------------------------------------" << endl ;
759 os << LogIO::NORMAL << "Number of scan = " << scanNum_ << endl ;
760 os << LogIO::NORMAL << "Number of data record = " << rowNum_ << endl ;
761 os << LogIO::NORMAL << "Length of data record = " << scanLen_ << " bytes" << endl ;
762 os << LogIO::NORMAL << "Allocated memory for spectral data = " << dataLen_ << " bytes" << endl ;
763 os << LogIO::NORMAL << "Max number of channel = " << chmax_ << endl ;
764 os << LogIO::NORMAL << "------------------------------------------------------------" << endl ;
765 os.post() ;
766
767}
768
769double NRODataset::toLSR( double v, double t, double x, double y )
770{
771 double vlsr ;
772
773 // get epoch
774 double tcent = t + 0.5*getIPTIM()/86400.0 ;
775 MEpoch me( Quantity( tcent, "d" ), MEpoch::UTC ) ;
776
777 // get antenna position
778 MPosition mp ;
779 if ( SITE.find( "45" ) != string::npos ) {
780 // 45m telescope
781 Double posx = -3.8710235e6 ;
782 Double posy = 3.4281068e6 ;
783 Double posz = 3.7240395e6 ;
784 mp = MPosition( MVPosition( posx, posy, posz ),
785 MPosition::ITRF ) ;
786 }
787 else {
788 // ASTE
789 Vector<Double> pos( 2 ) ;
790 pos[0] = -67.7031 ;
791 pos[1] = -22.9717 ;
792 Double sitealt = 4800.0 ;
793 mp = MPosition( MVPosition( Quantity( sitealt, "m" ),
794 Quantum< Vector<Double> >( pos, "deg" ) ),
795 MPosition::WGS84 ) ;
796 }
797
798 // get direction
799 MDirection md ;
800 if ( SCNCD == 0 ) {
801 // RADEC
802 if ( EPOCH == "B1950" ) {
803 md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
804 MDirection::B1950 ) ;
805 }
806 else {
807 md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
808 MDirection::J2000 ) ;
809 }
810 }
811 else if ( SCNCD == 1 ) {
812 // LB
813 md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
814 MDirection::GALACTIC ) ;
815 }
816 else {
817 // AZEL
818 md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
819 MDirection::AZEL ) ;
820 }
821
822 // to LSR
823 MeasFrame mf( me, mp, md ) ;
824 MFrequency::Convert tolsr( MFrequency::TOPO, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
825 vlsr = (double)(tolsr( Double(v) ).get( "Hz" ).getValue()) ;
826
827 return vlsr ;
828}
Note: See TracBrowser for help on using the repository browser.