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

Last change on this file since 2756 was 2748, checked in by Takeshi Nakazato, 12 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...

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