source: tags/asap-4.1.0/external-alma/atnf/PKSIO/NRODataset.cc@ 2752

Last change on this file since 2752 was 2643, checked in by ShinnosukeKawakami, 12 years ago

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