source: branches/hpc34/external-alma/atnf/PKSIO/NRODataset.cc@ 2596

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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