source: branches/plotter2/external-alma/atnf/PKSIO/NRODataset.cc@ 2887

Last change on this file since 2887 was 2785, 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...

Candidate bug fix for seg. fault on Mac.

File size: 44.7 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#define STRING2CHAR(s) const_cast<char *>((s).c_str())
49
50//#include <casa/namespace.h>
51
52using namespace std ;
53
54//
55// NRODataset
56//
57// Base class for NRO dataset.
58//
59
60// constructor
61NRODataset::NRODataset( string name )
62 : scanNum_(0),
63 rowNum_(0),
64 scanLen_(0),
65 dataLen_(0),
66 dataid_(-1),
67 filename_(name),
68 fp_(NULL),
69 same_(-1),
70 frec_()
71{
72 // size for common part of data
73 datasize_ = sizeof( char ) * 8 // LOFIL
74 + sizeof( char ) * 8 // VER
75 + sizeof( char ) * 16 // GROUP
76 + sizeof( char ) * 16 // PROJ
77 + sizeof( char ) * 24 // SCHED
78 + sizeof( char ) * 40 // OBSVR
79 + sizeof( char ) * 16 // LOSTM
80 + sizeof( char ) * 16 // LOETM
81 + sizeof( int ) * 2 // ARYNM, NSCAN
82 + sizeof( char ) * 120 // TITLE
83 + sizeof( char ) * 16 // OBJ
84 + sizeof( char ) * 8 // EPOCH
85 + sizeof( double ) * 4 // RA0, DEC0, GLNG0, GLAT0
86 + sizeof( int ) * 2 // NCALB, SCNCD
87 + sizeof( char ) * 120 // SCMOD
88 + sizeof( double ) // URVEL
89 + sizeof( char ) * 4 // VREF
90 + sizeof( char ) * 4 // VDEF
91 + sizeof( char ) * 8 // SWMOD
92 + sizeof( double ) * 8 // FRQSW, DBEAM, MLTOF, CMTQ, CMTE, CMTSOM, CMTNODE, CMTI
93 + sizeof( char ) * 24 // CMTTM
94 + sizeof( double ) * 6 // SBDX, SBDY, SBDZ1, SBDZ2, DAZP, DELP
95 + sizeof( int ) * 4 // CHBIND, NUMCH, CHMIN, CHMAX
96 + sizeof( double ) * 3 // ALCTM, IPTIM, PA
97 + sizeof( int ) * 3 // SCNLEN, SBIND, IBIT
98 + sizeof( char ) * 8 ; // SITE
99}
100
101// destructor
102NRODataset::~NRODataset()
103{
104 // release memory
105 releaseRecord() ;
106
107 // close file
108 close() ;
109}
110
111// data initialization
112void NRODataset::initializeCommon()
113{
114 LogIO os( LogOrigin( "NRODataset", "initialize()", WHERE ) ) ;
115
116 int arymax = arrayMax() ;
117
118 // check endian
119 open() ;
120 fseek( fp_, 144, SEEK_SET ) ;
121 int tmp ;
122 if( fread( &tmp, 1, sizeof(int), fp_ ) != sizeof(int) ) {
123 os << LogIO::SEVERE << "Error while checking endian of the file. " << LogIO::EXCEPTION ;
124 return ;
125 }
126 if ( ( 0 < tmp ) && ( tmp <= arymax ) ) {
127 same_ = 1 ;
128 os << LogIO::NORMAL << "same endian " << LogIO::POST ;
129 }
130 else {
131 same_ = 0 ;
132 os << LogIO::NORMAL << "different endian " << LogIO::POST ;
133 }
134 fseek( fp_, 0, SEEK_SET ) ;
135
136 // common part of calculation of data size and memory allocation
137 RX.resize( arymax ) ;
138 HPBW.resize( arymax ) ;
139 EFFA.resize( arymax ) ;
140 EFFB.resize( arymax ) ;
141 EFFL.resize( arymax ) ;
142 EFSS.resize( arymax ) ;
143 GAIN.resize( arymax ) ;
144 HORN.resize( arymax ) ;
145 POLTP.resize( arymax ) ;
146 POLDR.resize( arymax ) ;
147 POLAN.resize( arymax ) ;
148 DFRQ.resize( arymax ) ;
149 SIDBD.resize( arymax ) ;
150 REFN.resize( arymax ) ;
151 IPINT.resize( arymax ) ;
152 MULTN.resize( arymax ) ;
153 MLTSCF.resize( arymax ) ;
154 LAGWIND.resize( arymax ) ;
155 BEBW.resize( arymax ) ;
156 BERES.resize( arymax ) ;
157 CHWID.resize( arymax ) ;
158 ARRY.resize( arymax ) ;
159 NFCAL.resize( arymax ) ;
160 F0CAL.resize( arymax ) ;
161 FQCAL.resize( arymax ) ;
162 CHCAL.resize( arymax ) ;
163 CWCAL.resize( arymax ) ;
164 DSBFC.resize( arymax ) ;
165
166 for ( int i = 0 ; i < arymax ; i++ ) {
167 FQCAL[i].resize( 10 ) ;
168 CHCAL[i].resize( 10 ) ;
169 CWCAL[i].resize( 10 ) ;
170 }
171
172 // NRODataRecord
173 record_ = new NRODataRecord() ;
174 record_->LDATA = NULL ;
175
176 // reference frequency
177 refFreq_.resize( arymax, 0.0 ) ;
178}
179
180// fill data header
181int NRODataset::fillHeader()
182{
183 LogIO os( LogOrigin( "NRODataset", "fillHeader()", WHERE ) ) ;
184
185 // open file
186 if ( open() ) {
187 os << LogIO::SEVERE << "Error opening file " << filename_ << "." << LogIO::EXCEPTION ;
188 return -1 ;
189 }
190
191 // fill
192 int status = fillHeader( same_ ) ;
193
194 if ( status != 0 ) {
195 os << LogIO::SEVERE << "Error while reading header " << filename_ << "." << LogIO::EXCEPTION ;
196 return status ;
197 }
198
199 initArray();
200
201 show() ;
202
203 return status ;
204}
205
206int NRODataset::fillHeaderCommon( int sameEndian )
207{
208 LogIO os( LogOrigin( "NRODataset", "fillHeader()", WHERE ) ) ;
209
210 int arymax = arrayMax() ;
211
212 // make sure file pointer points a beginning of the file
213 fseek( fp_, 0, SEEK_SET ) ;
214
215 // read data header
216 LOFIL.resize(8);
217 if ( readHeader( STRING2CHAR(LOFIL), 8 ) == -1 ) {
218 os << LogIO::WARN << "Error while reading data LOFIL." << LogIO::POST ;
219 return -1 ;
220 }
221 // DEBUG
222 //cout << "LOFIL = " << LOFIL << endl ;
223 //
224 VER.resize(8);
225 if ( readHeader( STRING2CHAR(VER), 8 ) == -1 ) {
226 os << LogIO::WARN << "Error while reading data VER." << LogIO::POST ;
227 return -1 ;
228 }
229 // DEBUG
230 //cout << "VER = " << VER << endl ;
231 //
232 GROUP.resize(16);
233 if ( readHeader( STRING2CHAR(GROUP), 16 ) == -1 ) {
234 os << LogIO::WARN << "Error while reading data GROUP." << LogIO::POST ;
235 return -1 ;
236 }
237 // DEBUG
238 //cout << "GROUP = " << GROUP << endl ;
239 //
240 PROJ.resize(16);
241 if ( readHeader( STRING2CHAR(PROJ), 16 ) == -1 ) {
242 os << LogIO::WARN << "Error while reading data PROJ." << LogIO::POST ;
243 return -1 ;
244 }
245 // DEBUG
246 //cout << "PROJ = " << PROJ << endl ;
247 //
248 SCHED.resize(24);
249 if ( readHeader( STRING2CHAR(SCHED), 24 ) == -1 ) {
250 os << LogIO::WARN << "Error while reading data SCHED." << LogIO::POST ;
251 return -1 ;
252 }
253 // DEBUG
254 //cout << "SCHED = " << SCHED << endl ;
255 //
256 OBSVR.resize(40);
257 if ( readHeader( STRING2CHAR(OBSVR), 40 ) == -1 ) {
258 os << LogIO::WARN << "Error while reading data OBSVR." << LogIO::POST ;
259 return -1 ;
260 }
261 // DEBUG
262 //cout << "OBSVR = " << OBSVR << endl ;
263 //
264 LOSTM.resize(16);
265 if ( readHeader( STRING2CHAR(LOSTM), 16 ) == -1 ) {
266 os << LogIO::WARN << "Error while reading data LOSTM." << LogIO::POST ;
267 return -1 ;
268 }
269 // DEBUG
270 //cout << "LOSTM = " << LOSTM << endl ;
271 //
272 LOETM.resize(16);
273 if ( readHeader( STRING2CHAR(LOETM), 16 ) == -1 ) {
274 os << LogIO::WARN << "Error while reading data LOETM." << LogIO::POST ;
275 return -1 ;
276 }
277 // DEBUG
278 //cout << "LOETM = " << LOETM << endl ;
279 //
280 if ( readHeader( ARYNM, sameEndian ) == -1 ) {
281 os << LogIO::WARN << "Error while reading data ARYNM." << LogIO::POST ;
282 return -1 ;
283 }
284 // DEBUG
285 //cout << "ARYNM = " << ARYNM << endl ;
286 //
287 if ( readHeader( NSCAN, sameEndian ) == -1 ) {
288 os << LogIO::WARN << "Error while reading data NSCAN." << LogIO::POST ;
289 return -1 ;
290 }
291 // DEBUG
292 //cout << "NSCAN = " << NSCAN << endl ;
293 //
294 TITLE.resize(120);
295 if ( readHeader( STRING2CHAR(TITLE), 120 ) == -1 ) {
296 os << LogIO::WARN << "Error while reading data TITLE." << LogIO::POST ;
297 return -1 ;
298 }
299 // DEBUG
300 //cout << "TITLE = " << TITLE << endl ;
301 //
302 OBJ.resize(16);
303 if ( readHeader( STRING2CHAR(OBJ), 16 ) == -1 ) {
304 os << LogIO::WARN << "Error while reading data OBJ." << LogIO::POST ;
305 return -1 ;
306 }
307 // DEBUG
308 //cout << "OBJ = " << OBJ << endl ;
309 //
310 EPOCH.resize(8);
311 if ( readHeader( STRING2CHAR(EPOCH), 8 ) == -1 ) {
312 os << LogIO::WARN << "Error while reading data EPOCH." << LogIO::POST ;
313 return -1 ;
314 }
315 // DEBUG
316 //cout << "EPOCH = " << EPOCH << endl ;
317 //
318 if ( readHeader( RA0, sameEndian ) == -1 ) {
319 os << LogIO::WARN << "Error while reading data RA0." << LogIO::POST ;
320 return -1 ;
321 }
322 // DEBUG
323 //cout << "RA0 = " << RA0 << endl ;
324 //
325 if ( readHeader( DEC0, sameEndian ) == -1 ) {
326 os << LogIO::WARN << "Error while reading data DEC0." << LogIO::POST ;
327 return -1 ;
328 }
329 // DEBUG
330 //cout << "DEC0 = " << DEC0 << endl ;
331 //
332 if ( readHeader( GLNG0, sameEndian ) == -1 ) {
333 os << LogIO::WARN << "Error while reading data GLNG0." << LogIO::POST ;
334 return -1 ;
335 }
336 // DEBUG
337 //cout << "GLNG0 = " << GLNG0 << endl ;
338 //
339 if ( readHeader( GLAT0, sameEndian ) == -1 ) {
340 os << LogIO::WARN << "Error while reading data GLAT0." << LogIO::POST ;
341 return -1 ;
342 }
343 // DEBUG
344 //cout << "GLAT0 = " << GLAT0 << endl ;
345 //
346 if ( readHeader( NCALB, sameEndian ) == -1 ) {
347 os << LogIO::WARN << "Error while reading data NCALB." << LogIO::POST ;
348 return -1 ;
349 }
350 // DEBUG
351 //cout << "NCALB = " << NCALB << endl ;
352 //
353 if ( readHeader( SCNCD, sameEndian ) == -1 ) {
354 os << LogIO::WARN << "Error while reading data SCNCD." << LogIO::POST ;
355 return -1 ;
356 }
357 // DEBUG
358 //cout << "SCNCD = " << SCNCD << endl ;
359 //
360 SCMOD.resize(120);
361 if ( readHeader( STRING2CHAR(SCMOD), 120 ) == -1 ) {
362 os << LogIO::WARN << "Error while reading data SCMOD." << LogIO::POST ;
363 return -1 ;
364 }
365 // DEBUG
366 //cout << "SCMOD = " << SCMOD << endl ;
367 //
368 if ( readHeader( URVEL, sameEndian ) == -1 ) {
369 os << LogIO::WARN << "Error while reading data URVEL." << LogIO::POST ;
370 return -1 ;
371 }
372 // DEBUG
373 //cout << "URVEL = " << URVEL << endl ;
374 //
375 VREF.resize(4);
376 if ( readHeader( STRING2CHAR(VREF), 4 ) == -1 ) {
377 os << LogIO::WARN << "Error while reading data VREF." << LogIO::POST ;
378 return -1 ;
379 }
380 // DEBUG
381 //cout << "VREF = " << VREF << endl ;
382 //
383 VDEF.resize(4);
384 if ( readHeader( STRING2CHAR(VDEF), 4 ) == -1 ) {
385 os << LogIO::WARN << "Error while reading data VDEF." << LogIO::POST ;
386 return -1 ;
387 }
388 // DEBUG
389 //cout << "VDEF = " << VDEF << endl ;
390 //
391 SWMOD.resize(8);
392 if ( readHeader( STRING2CHAR(SWMOD), 8 ) == -1 ) {
393 os << LogIO::WARN << "Error while reading data SWMOD." << LogIO::POST ;
394 return -1 ;
395 }
396 SWMOD += "::OTF" ;
397 // DEBUG
398 //cout << "SWMOD = " << SWMOD << endl ;
399 //
400 if ( readHeader( FRQSW, sameEndian ) == -1 ) {
401 os << LogIO::WARN << "Error while reading data FRQSW." << LogIO::POST ;
402 return -1 ;
403 }
404 // DEBUG
405 //cout << "FRQSW = " << FRQSW << endl ;
406 //
407 if ( readHeader( DBEAM, sameEndian ) == -1 ) {
408 os << LogIO::WARN << "Error while reading data DBEAM." << LogIO::POST ;
409 return -1 ;
410 }
411 // DEBUG
412 //cout << "DBEAM = " << DBEAM << endl ;
413 //
414 if ( readHeader( MLTOF, sameEndian ) == -1 ) {
415 os << LogIO::WARN << "Error while reading data MLTOF." << LogIO::POST ;
416 return -1 ;
417 }
418 // DEBUG
419 //cout << "MLTOF = " << MLTOF << endl ;
420 //
421 if ( readHeader( CMTQ, sameEndian ) == -1 ) {
422 os << LogIO::WARN << "Error while reading data CMTQ." << LogIO::POST ;
423 return -1 ;
424 }
425 // DEBUG
426 //cout << "CMTQ = " << CMTQ << endl ;
427 //
428 if ( readHeader( CMTE, sameEndian ) == -1 ) {
429 os << LogIO::WARN << "Error while reading data CMTE." << LogIO::POST ;
430 return -1 ;
431 }
432 // DEBUG
433 //cout << "CMTE = " << CMTE << endl ;
434 //
435 if ( readHeader( CMTSOM, sameEndian ) == -1 ) {
436 os << LogIO::WARN << "Error while reading data CMTSOM." << LogIO::POST ;
437 return -1 ;
438 }
439 // DEBUG
440 //cout << "CMTSOM = " << CMTSOM << endl ;
441 //
442 if ( readHeader( CMTNODE, sameEndian ) == -1 ) {
443 os << LogIO::WARN << "Error while reading data CMTNODE." << LogIO::POST ;
444 return -1 ;
445 }
446 // DEBUG
447 //cout << "CMTNODE = " << CMTNODE << endl ;
448 //
449 if ( readHeader( CMTI, sameEndian ) == -1 ) {
450 os << LogIO::WARN << "Error while reading data CMTI." << LogIO::POST ;
451 return -1 ;
452 }
453 // DEBUG
454 //cout << "CMTI = " << CMTI << endl ;
455 //
456 CMTTM.resize(24);
457 if ( readHeader( STRING2CHAR(CMTTM), 24 ) == -1 ) {
458 os << LogIO::WARN << "Error while reading data CMTTM." << LogIO::POST ;
459 return -1 ;
460 }
461 // DEBUG
462 //cout << "CMTTM = " << CMTTM << endl ;
463 //
464 if ( readHeader( SBDX, sameEndian ) == -1 ) {
465 os << LogIO::WARN << "Error while reading data SBDX." << LogIO::POST ;
466 return -1 ;
467 }
468 // DEBUG
469 //cout << "SBDX = " << SBDX << endl ;
470 //
471 if ( readHeader( SBDY, sameEndian ) == -1 ) {
472 os << LogIO::WARN << "Error while reading data SBDY." << LogIO::POST ;
473 return -1 ;
474 }
475 // DEBUG
476 //cout << "SBDY = " << SBDY << endl ;
477 //
478 if ( readHeader( SBDZ1, sameEndian ) == -1 ) {
479 os << LogIO::WARN << "Error while reading data SBDZ1." << LogIO::POST ;
480 return -1 ;
481 }
482 // DEBUG
483 //cout << "SBDZ1 = " << SBDZ1 << endl ;
484 //
485 if ( readHeader( SBDZ2, sameEndian ) == -1 ) {
486 os << LogIO::WARN << "Error while reading data SBDZ2." << LogIO::POST ;
487 return -1 ;
488 }
489 // DEBUG
490 //cout << "SBDZ2 = " << SBDZ2 << endl ;
491 //
492 if ( readHeader( DAZP, sameEndian ) == -1 ) {
493 os << LogIO::WARN << "Error while reading data DAZP." << LogIO::POST ;
494 return -1 ;
495 }
496 // DEBUG
497 //cout << "DAZP = " << DAZP << endl ;
498 //
499 if ( readHeader( DELP, sameEndian ) == -1 ) {
500 os << LogIO::WARN << "Error while reading data DELP." << LogIO::POST ;
501 return -1 ;
502 }
503 // DEBUG
504 //cout << "DELP = " << DELP << endl ;
505 //
506 if ( readHeader( CHBIND, sameEndian ) == -1 ) {
507 os << LogIO::WARN << "Error while reading data CHBIND." << LogIO::POST ;
508 return -1 ;
509 }
510 // DEBUG
511 //cout << "CHBIND = " << CHBIND << endl ;
512 //
513 if ( readHeader( NUMCH, sameEndian ) == -1 ) {
514 os << LogIO::WARN << "Error while reading data NUMCH." << LogIO::POST ;
515 return -1 ;
516 }
517 // DEBUG
518 //cout << "NUMCH = " << NUMCH << endl ;
519 //
520 if ( readHeader( CHMIN, sameEndian ) == -1 ) {
521 os << LogIO::WARN << "Error while reading data CHMIN." << LogIO::POST ;
522 return -1 ;
523 }
524 // DEBUG
525 //cout << "CHMIN = " << CHMIN << endl ;
526 //
527 if ( readHeader( CHMAX, sameEndian ) == -1 ) {
528 os << LogIO::WARN << "Error while reading data CHMAX." << LogIO::POST ;
529 return -1 ;
530 }
531 // DEBUG
532 //cout << "CHMAX = " << CHMAX << endl ;
533 //
534 if ( readHeader( ALCTM, sameEndian ) == -1 ) {
535 os << LogIO::WARN << "Error while reading data ALCTM." << LogIO::POST ;
536 return -1 ;
537 }
538 // DEBUG
539 //cout << "ALCTM = " << ALCTM << endl ;
540 //
541 if ( readHeader( IPTIM, sameEndian ) == -1 ) {
542 os << LogIO::WARN << "Error while reading data IPTIM." << LogIO::POST ;
543 return -1 ;
544 }
545 // DEBUG
546 //cout << "IPTIM = " << IPTIM << endl ;
547 //
548 if ( readHeader( PA, sameEndian ) == -1 ) {
549 os << LogIO::WARN << "Error while reading data PA." << LogIO::POST ;
550 return -1 ;
551 }
552 // DEBUG
553 //cout << "PA = " << PA << endl ;
554 //
555 for ( int i = 0 ; i < arymax ; i++ ) {
556 RX[i].resize(16);
557 if ( readHeader( STRING2CHAR(RX[i]), 16 ) == -1 ) {
558 os << LogIO::WARN << "Error while reading data RX[" << i << "]." << LogIO::POST ;
559 return -1 ;
560 }
561 }
562 // DEBUG
563// nro_debug_output( "RX", arymax, RX ) ;
564 //
565 for ( int i = 0 ; i < arymax ; i++ ) {
566 if ( readHeader( HPBW[i], sameEndian ) == -1 ) {
567 os << LogIO::WARN << "Error while reading data HPBW[" << i << "]." << LogIO::POST ;
568 return -1 ;
569 }
570 }
571 // DEBUG
572// nro_debug_output( "HPBW", arymax, HPBW ) ;
573 //
574 for ( int i = 0 ; i < arymax ; i++ ) {
575 if ( readHeader( EFFA[i], sameEndian ) == -1 ) {
576 os << LogIO::WARN << "Error while reading data EFFA[" << i << "]." << LogIO::POST ;
577 return -1 ;
578 }
579 }
580 // DEBUG
581// nro_debug_output( "EFFA", arymax, EFFA ) ;
582 //
583 for ( int i = 0 ; i < arymax ; i++ ) {
584 if ( readHeader( EFFB[i], sameEndian ) == -1 ) {
585 os << LogIO::WARN << "Error while reading data EFFB[" << i << "]." << LogIO::POST ;
586 return -1 ;
587 }
588 }
589 // DEBUG
590// nro_debug_output( "EFFB", arymax, EFFB ) ;
591 //
592 for ( int i = 0 ; i < arymax ; i++ ) {
593 if ( readHeader( EFFL[i], sameEndian ) == -1 ) {
594 os << LogIO::WARN << "Error while reading data EFFL[" << i << "]." << LogIO::POST ;
595 return -1 ;
596 }
597 }
598 // DEBUG
599// nro_debug_output( "EFFL", arymax, EFFL ) ;
600 //
601 for ( int i = 0 ; i < arymax ; i++ ) {
602 if ( readHeader( EFSS[i], sameEndian ) == -1 ) {
603 os << LogIO::WARN << "Error while reading data EFSS[" << i << "]." << LogIO::POST ;
604 return -1 ;
605 }
606 }
607 // DEBUG
608// nro_debug_output( "EFSS", arymax, EFSS ) ;
609 //
610 for ( int i = 0 ; i < arymax ; i++) {
611 if ( readHeader( GAIN[i], sameEndian ) == -1 ) {
612 os << LogIO::WARN << "Error while reading data GAIN[" << i << "]." << LogIO::POST ;
613 return -1 ;
614 }
615 }
616 // DEBUG
617// nro_debug_output( "GAIN", arymax, GAIN ) ;
618 //
619 for ( int i = 0 ; i < arymax ; i++) {
620 HORN[i].resize(4);
621 if ( readHeader( STRING2CHAR(HORN[i]), 4 ) == -1 ) {
622 os << LogIO::WARN << "Error while reading data HORN[" << i << "]." << LogIO::POST ;
623 return -1 ;
624 }
625 }
626 // DEBUG
627// nro_debug_output( "HORN", arymax, HORN ) ;
628 //
629 for ( int i = 0 ; i < arymax ; i++) {
630 POLTP[i].resize(4);
631 if ( readHeader( STRING2CHAR(POLTP[i]), 4 ) == -1 ) {
632 os << LogIO::WARN << "Error while reading data POLTP[" << i << "]." << LogIO::POST ;
633 return -1 ;
634 }
635 }
636 // DEBUG
637// nro_debug_output( "POLTP", arymax, POLTP ) ;
638 //
639 for ( int i = 0 ; i < arymax ; i++) {
640 if ( readHeader( POLDR[i], sameEndian ) == -1 ) {
641 os << LogIO::WARN << "Error while reading data POLDR[" << i << "]." << LogIO::POST ;
642 return -1 ;
643 }
644 }
645 // DEBUG
646// nro_debug_output( "POLDR", arymax, POLDR ) ;
647 //
648 for ( int i = 0 ; i < arymax ; i++) {
649 if ( readHeader( POLAN[i], sameEndian ) == -1 ) {
650 os << LogIO::WARN << "Error while reading data POLAN[" << i << "]." << LogIO::POST ;
651 return -1 ;
652 }
653 }
654 // DEBUG
655// nro_debug_output( "POLAN", arymax, POLAN ) ;
656 //
657 for ( int i = 0 ; i < arymax ; i++) {
658 if ( readHeader( DFRQ[i], sameEndian ) == -1 ) {
659 os << LogIO::WARN << "Error while reading data DFRQ[" << i << "]." << LogIO::POST ;
660 return -1 ;
661 }
662 }
663 // DEBUG
664// nro_debug_output( "DFRQ", arymax, DFRQ ) ;
665 //
666 for ( int i = 0 ; i < arymax ; i++) {
667 SIDBD[i].resize(4);
668 if ( readHeader( STRING2CHAR(SIDBD[i]), 4 ) == -1 ) {
669 os << LogIO::WARN << "Error while reading data SIDBD[" << i << "]." << LogIO::POST ;
670 return -1 ;
671 }
672 }
673 // DEBUG
674// nro_debug_output( "SIDBD", arymax, SIDBD ) ;
675 //
676 for ( int i = 0 ; i < arymax ; i++) {
677 if ( readHeader( REFN[i], sameEndian ) == -1 ) {
678 os << LogIO::WARN << "Error while reading data REFN[" << i << "]." << LogIO::POST ;
679 return -1 ;
680 }
681 }
682 // DEBUG
683// nro_debug_output( "REFN", arymax, REFN ) ;
684 //
685 for ( int i = 0 ; i < arymax ; i++) {
686 if ( readHeader( IPINT[i], sameEndian ) == -1 ) {
687 os << LogIO::WARN << "Error while reading data IPINT[" << i << "]." << LogIO::POST ;
688 return -1 ;
689 }
690 }
691 // DEBUG
692// nro_debug_output( "IPINT", arymax, IPINT ) ;
693 //
694 for ( int i = 0 ; i < arymax ; i++) {
695 if ( readHeader( MULTN[i], sameEndian ) == -1 ) {
696 os << LogIO::WARN << "Error while reading data MULTN[" << i << "]." << LogIO::POST ;
697 return -1 ;
698 }
699 }
700 // DEBUG
701// nro_debug_output( "MULTN", arymax, MULTN ) ;
702 //
703 for ( int i = 0 ; i < arymax ; i++) {
704 if ( readHeader( MLTSCF[i], sameEndian ) == -1 ) {
705 os << LogIO::WARN << "Error while reading data MLTSCF[" << i << "]." << LogIO::POST ;
706 return -1 ;
707 }
708 }
709 // DEBUG
710// nro_debug_output( "MLTSCF", arymax, MLTSCF ) ;
711 //
712 for ( int i = 0 ; i < arymax ; i++) {
713 LAGWIND[i].resize(8);
714 if ( readHeader( STRING2CHAR(LAGWIND[i]), 8 ) == -1 ) {
715 os << LogIO::WARN << "Error while reading data LAGWIND[" << i << "]." << LogIO::POST ;
716 return -1 ;
717 }
718 }
719 // DEBUG
720// nro_debug_output( "LAGWIND", arymax, LAGWIND ) ;
721 //
722 for ( int i = 0 ; i < arymax ; i++) {
723 if ( readHeader( BEBW[i], sameEndian ) == -1 ) {
724 os << LogIO::WARN << "Error while reading data BEBW[" << i << "]." << LogIO::POST ;
725 return -1 ;
726 }
727 }
728 // DEBUG
729// nro_debug_output( "BEBW", arymax, BEBW ) ;
730 //
731 for ( int i = 0 ; i < arymax ; i++) {
732 if ( readHeader( BERES[i], sameEndian ) == -1 ) {
733 os << LogIO::WARN << "Error while reading data BERES[" << i << "]." << LogIO::POST ;
734 return -1 ;
735 }
736 }
737 // DEBUG
738// nro_debug_output( "BERES", arymax, BERES ) ;
739 //
740 for ( int i = 0 ; i < arymax ; i++) {
741 if ( readHeader( CHWID[i], sameEndian ) == -1 ) {
742 os << LogIO::WARN << "Error while reading data CHWID[" << i << "]." << LogIO::POST ;
743 return -1 ;
744 }
745 }
746 // DEBUG
747// nro_debug_output( "CHWID", arymax, CHWID ) ;
748 //
749 for ( int i = 0 ; i < arymax ; i++) {
750 if ( readHeader( ARRY[i], sameEndian ) == -1 ) {
751 os << LogIO::WARN << "Error while reading data ARRY[" << i << "]." << LogIO::POST ;
752 return -1 ;
753 }
754 }
755 // DEBUG
756// nro_debug_output( "ARRY", arymax, ARRY ) ;
757 //
758 for ( int i = 0 ; i < arymax ; i++) {
759 if ( readHeader( NFCAL[i], sameEndian ) == -1 ) {
760 os << LogIO::WARN << "Error while reading data NFCAL[" << i << "]." << LogIO::POST ;
761 return -1 ;
762 }
763 }
764 // DEBUG
765// nro_debug_output( "NFCAL", arymax, NFCAL ) ;
766 //
767 for ( int i = 0 ; i < arymax ; i++) {
768 if ( readHeader( F0CAL[i], sameEndian ) == -1 ) {
769 os << LogIO::WARN << "Error while reading data F0CAL[" << i << "]." << LogIO::POST ;
770 return -1 ;
771 }
772 }
773 // DEBUG
774// nro_debug_output( "F0CAL", arymax, F0CAL ) ;
775 //
776 for ( int i = 0 ; i < arymax ; i++) {
777 for ( int j = 0 ; j < 10 ; j++ ) {
778 if ( readHeader( FQCAL[i][j], sameEndian ) == -1 ) {
779 os << LogIO::WARN << "Error while reading data FQCAL[" << i << "][" << j << "]." << LogIO::POST ;
780 return -1 ;
781 }
782 }
783 }
784 // DEBUG
785// nro_debug_output( "FQCAL", arymax, 10, FQCAL ) ;
786 //
787 for ( int i = 0 ; i < arymax ; i++) {
788 for ( int j = 0 ; j < 10 ; j++ ) {
789 if ( readHeader( CHCAL[i][j], sameEndian ) == -1 ) {
790 os << LogIO::WARN << "Error while reading data CHCAL[" << i << "][" << j << "]." << LogIO::POST ;
791 return -1 ;
792 }
793 }
794 }
795 // DEBUG
796// nro_debug_output( "CHCAL", arymax, 10, CHCAL ) ;
797 //
798 for ( int i = 0 ; i < arymax ; i++) {
799 for ( int j = 0 ; j < 10 ; j++ ) {
800 if ( readHeader( CWCAL[i][j], sameEndian ) == -1 ) {
801 os << LogIO::WARN << "Error while reading data CWCAL[" << i << "][" << j << "]." << LogIO::POST ;
802 return -1 ;
803 }
804 }
805 }
806 // DEBUG
807// nro_debug_output( "CWCAL", arymax, 10, CWCAL ) ;
808 //
809 if ( readHeader( SCNLEN, sameEndian ) == -1 ) {
810 os << LogIO::WARN << "Error while reading data SCNLEN." << LogIO::POST ;
811 return -1 ;
812 }
813 // DEBUG
814 //cout << "SCNLEN = " << SCNLEN << endl ;
815 //
816 if ( readHeader( SBIND, sameEndian ) == -1 ) {
817 os << LogIO::WARN << "Error while reading data SBIND." << LogIO::POST ;
818 return -1 ;
819 }
820 // DEBUG
821 //cout << "SBIND = " << SBIND << endl ;
822 //
823 if ( readHeader( IBIT, sameEndian ) == -1 ) {
824 os << LogIO::WARN << "Error while reading data IBIT." << LogIO::POST ;
825 return -1 ;
826 }
827 // DEBUG
828 //cout << "IBIT = " << IBIT << endl ;
829 //
830 SITE.resize(8);
831 if ( readHeader( STRING2CHAR(SITE), 8 ) == -1 ) {
832 os << LogIO::WARN << "Error while reading data SITE." << LogIO::POST ;
833 return -1 ;
834 }
835 // DEBUG
836 //cout << "SITE = " << SITE << endl ;
837 //
838
839 //scanNum_ = NSCAN + 1 ; // includes ZERO scan
840 scanLen_ = SCNLEN ;
841 dataLen_ = scanLen_ - SCAN_HEADER_SIZE ;
842 scanNum_ = getScanNum();
843 rowNum_ = scanNum_ * ARYNM ;
844 chmax_ = (int) ( dataLen_ * 8 / IBIT ) ;
845 record_->LDATA = new char[dataLen_] ;
846
847 return 0 ;
848}
849
850void NRODataset::convertEndian( int &value )
851{
852 char volatile *first = reinterpret_cast<char volatile *>( &value ) ;
853 char volatile *last = first + sizeof( int ) ;
854 std::reverse( first, last ) ;
855}
856
857void NRODataset::convertEndian( float &value )
858{
859 char volatile *first = reinterpret_cast<char volatile *>( &value ) ;
860 char volatile *last = first + sizeof( float ) ;
861 std::reverse( first, last ) ;
862}
863
864void NRODataset::convertEndian( double &value )
865{
866 char volatile *first = reinterpret_cast<char volatile *>( &value ) ;
867 char volatile *last = first + sizeof( double ) ;
868 std::reverse( first, last ) ;
869}
870
871int NRODataset::readHeader( char *v, int size )
872{
873 if ( (int)( fread( v, 1, size, fp_ ) ) != size ) {
874 return -1 ;
875 }
876 return 0 ;
877}
878
879int NRODataset::readHeader( int &v, int b )
880{
881 if ( fread( &v, 1, sizeof(int), fp_ ) != sizeof(int) ) {
882 return -1 ;
883 }
884
885 if ( b == 0 )
886 convertEndian( v ) ;
887
888 return 0 ;
889}
890
891int NRODataset::readHeader( float &v, int b )
892{
893 if ( fread( &v, 1, sizeof(float), fp_ ) != sizeof(float) ) {
894 return -1 ;
895 }
896
897 if ( b == 0 )
898 convertEndian( v ) ;
899
900 return 0 ;
901}
902
903int NRODataset::readHeader( double &v, int b )
904{
905 if ( fread( &v, 1, sizeof(double), fp_ ) != sizeof(double) ) {
906 return -1 ;
907 }
908
909 if ( b == 0 )
910 convertEndian( v ) ;
911
912 return 0 ;
913}
914
915void NRODataset::convertEndian( NRODataRecord &r )
916{
917 convertEndian( r.ISCAN ) ;
918 convertEndian( r.DSCX ) ;
919 convertEndian( r.DSCY ) ;
920 convertEndian( r.SCX ) ;
921 convertEndian( r.SCY ) ;
922 convertEndian( r.PAZ ) ;
923 convertEndian( r.PEL ) ;
924 convertEndian( r.RAZ ) ;
925 convertEndian( r.REL ) ;
926 convertEndian( r.XX ) ;
927 convertEndian( r.YY ) ;
928 convertEndian( r.TEMP ) ;
929 convertEndian( r.PATM ) ;
930 convertEndian( r.PH2O ) ;
931 convertEndian( r.VWIND ) ;
932 convertEndian( r.DWIND ) ;
933 convertEndian( r.TAU ) ;
934 convertEndian( r.TSYS ) ;
935 convertEndian( r.BATM ) ;
936 convertEndian( r.LINE ) ;
937 for ( int i = 0 ; i < 4 ; i++ )
938 convertEndian( r.IDMY1[i] ) ;
939 convertEndian( r.VRAD ) ;
940 convertEndian( r.FREQ0 ) ;
941 convertEndian( r.FQTRK ) ;
942 convertEndian( r.FQIF1 ) ;
943 convertEndian( r.ALCV ) ;
944 for ( int i = 0 ; i < 2 ; i++ )
945 for ( int j = 0 ; j < 2 ; j++ )
946 convertEndian( r.OFFCD[i][j] ) ;
947 convertEndian( r.IDMY0 ) ;
948 convertEndian( r.IDMY2 ) ;
949 convertEndian( r.DPFRQ ) ;
950 convertEndian( r.SFCTR ) ;
951 convertEndian( r.ADOFF ) ;
952}
953
954void NRODataset::releaseRecord()
955{
956 if ( !record_.null() ) {
957 record_ = NULL ;
958 }
959 dataid_ = -1 ;
960}
961
962// Get specified scan
963NRODataRecord *NRODataset::getRecord( int i )
964{
965 // DEBUG
966 //cout << "NRODataset::getRecord() Start " << i << endl ;
967 //
968 if ( i < 0 || i >= rowNum_ ) {
969 LogIO os( LogOrigin( "NRODataset", "getRecord()", WHERE ) ) ;
970 //cerr << "NRODataset::getRecord() data index out of range." << endl ;
971 os << LogIO::SEVERE << "data index " << i << " out of range. return NULL." << LogIO::POST ;
972 return NULL ;
973 }
974
975 if ( i == dataid_ ) {
976 return &(*record_) ;
977 }
978
979 // DEBUG
980 //cout << "NRODataset::getData() Get new dataset" << endl ;
981 //
982 // read data
983 int status = fillRecord( i ) ;
984 if ( status == 0 ) {
985 dataid_ = i ;
986 }
987 else {
988 LogIO os( LogOrigin( "NRODataset", "getRecord()", WHERE ) ) ;
989 //cerr << "NRODataset::getRecord() error while reading data " << i << endl ;
990 os << LogIO::SEVERE << "error while reading data " << i << ". return NULL." << LogIO::EXCEPTION ;
991 dataid_ = -1 ;
992 return NULL ;
993 }
994
995 return &(*record_) ;
996}
997
998int NRODataset::fillRecord( int i )
999{
1000 int status = 0 ;
1001
1002 status = open() ;
1003 if ( status != 0 )
1004 return status ;
1005
1006
1007 // fill NRODataset
1008 long offset = getDataSize() + scanLen_ * i ;
1009 // DEBUG
1010 //cout << "NRODataset::fillRecord() offset (header) = " << offset << endl ;
1011 //cout << "NRODataset::fillRecord() sizeof(NRODataRecord) = " << sizeof( NRODataRecord ) << " byte" << endl ;
1012 fseek( fp_, offset, SEEK_SET ) ;
1013 if ( (int)fread( &(*record_), 1, SCAN_HEADER_SIZE, fp_ ) != SCAN_HEADER_SIZE ) {
1014 //cerr << "Failed to read scan header: " << i << endl ;
1015 LogIO os( LogOrigin( "NRODataset", "fillRecord()", WHERE ) ) ;
1016 os << LogIO::SEVERE << "Failed to read scan header for " << i << "th row." << LogIO::POST ;
1017 return -1 ;
1018 }
1019 if ( (int)fread( &(*record_->LDATA), 1, dataLen_, fp_ ) != dataLen_ ) {
1020 //cerr << "Failed to read spectral data: " << i << endl ;
1021 LogIO os( LogOrigin( "NRODataset", "fillRecord()", WHERE ) ) ;
1022 os << LogIO::SEVERE << "Failed to read spectral data for " << i << "th row." << LogIO::POST ;
1023 return -1 ;
1024 }
1025
1026 if ( same_ == 0 ) {
1027 convertEndian( *record_ ) ;
1028 }
1029
1030 // DWIND unit conversion (deg -> rad)
1031 record_->DWIND = record_->DWIND * M_PI / 180.0 ;
1032
1033 return status ;
1034}
1035
1036// open
1037int NRODataset::open()
1038{
1039 int status = 0 ;
1040
1041 if ( fp_ == NULL ) {
1042 if ( (fp_ = fopen( filename_.c_str(), "rb" )) == NULL )
1043 status = -1 ;
1044 else
1045 status = 0 ;
1046 }
1047
1048 return status ;
1049}
1050
1051// close
1052void NRODataset::close()
1053{
1054 // DEBUG
1055 //cout << "NRODataset::close() close file" << endl ;
1056 //
1057 if ( fp_ != NULL )
1058 fclose( fp_ ) ;
1059 fp_ = NULL ;
1060}
1061
1062// get spectrum
1063vector< vector<double> > NRODataset::getSpectrum()
1064{
1065 vector< vector<double> > spec(rowNum_);
1066
1067 for ( int i = 0 ; i < rowNum_ ; i++ ) {
1068 spec[i] = getSpectrum( i ) ;
1069 }
1070
1071 return spec ;
1072}
1073
1074vector<double> NRODataset::getSpectrum( int i )
1075{
1076 // DEBUG
1077 //cout << "NRODataset::getSpectrum() start process (" << i << ")" << endl ;
1078 //
1079 // size of spectrum is not chmax_ but dataset_->getNCH() after binding
1080 const int nchan = NUMCH ;
1081 vector<double> spec( chmax_ ) ; // spectrum "before" binding
1082 // DEBUG
1083 //cout << "NRODataset::getSpectrum() nchan = " << nchan << " chmax_ = " << chmax_ << endl ;
1084 //
1085
1086 const NRODataRecord *record = getRecord( i ) ;
1087
1088 const int bit = IBIT ; // fixed to 12 bit
1089 double scale = record->SFCTR ;
1090 // DEBUG
1091 //cout << "NRODataset::getSpectrum() scale = " << scale << endl ;
1092 //
1093 double offset = record->ADOFF ;
1094 // DEBUG
1095 //cout << "NRODataset::getSpectrum() offset = " << offset << endl ;
1096 //
1097 if ( ( scale == 0.0 ) && ( offset == 0.0 ) ) {
1098 //cerr << "NRODataset::getSpectrum() zero spectrum (" << i << ")" << endl ;
1099 LogIO os( LogOrigin("NRODataset","getSpectrum",WHERE) ) ;
1100 os << LogIO::WARN << "zero spectrum for row " << i << LogIO::POST ;
1101 if ( spec.size() != (unsigned int)nchan )
1102 spec.resize( nchan ) ;
1103 for ( vector<double>::iterator i = spec.begin() ;
1104 i != spec.end() ; i++ )
1105 *i = 0.0 ;
1106 return spec ;
1107 }
1108 unsigned char *cdata = (unsigned char *)&(*record->LDATA) ;
1109 vector<double> mscale = MLTSCF ;
1110 double dscale = mscale[getIndex( i )] ;
1111 int cbind = CHBIND ;
1112 int chmin = CHMIN ;
1113
1114 // char -> int -> double
1115 vector<double>::iterator iter = spec.begin() ;
1116
1117 static const int shift_right[] = {
1118 4, 0
1119 };
1120 static const int start_pos[] = {
1121 0, 1
1122 };
1123 static const int incr[] = {
1124 0, 3
1125 };
1126 int j = 0 ;
1127 for ( int i = 0 ; i < chmax_ ; i++ ) {
1128 // char -> int
1129 int ivalue = 0 ;
1130 if ( bit == 12 ) { // 12 bit qunatization
1131 const int ialt = i & 1 ;
1132 const int idx = j + start_pos[ialt];
1133 const unsigned tmp = unsigned(cdata[idx]) << 8 | cdata[idx + 1];
1134 ivalue = int((tmp >> shift_right[ialt]) & 0xFFF);
1135 j += incr[ialt];
1136 }
1137
1138 if ( ( ivalue < 0 ) || ( ivalue > 4096 ) ) {
1139 //cerr << "NRODataset::getSpectrum() ispec[" << i << "] is out of range" << endl ;
1140 LogIO os( LogOrigin( "NRODataset", "getSpectrum", WHERE ) ) ;
1141 os << LogIO::SEVERE << "ivalue for row " << i << " is out of range" << LogIO::EXCEPTION ;
1142 if ( spec.size() != (unsigned int)nchan )
1143 spec.resize( nchan ) ;
1144 for ( vector<double>::iterator i = spec.begin() ;
1145 i != spec.end() ; i++ )
1146 *i = 0.0 ;
1147 return spec ;
1148 }
1149 // DEBUG
1150 //cout << "NRODataset::getSpectrum() ispec[" << i << "] = " << ispec[i] << endl ;
1151 //
1152
1153 // int -> double
1154 *iter = (double)( ivalue * scale + offset ) * dscale ;
1155 // DEBUG
1156 //cout << "NRODataset::getSpectrum() spec[" << i << "] = " << *iter << endl ;
1157 //
1158 iter++ ;
1159 }
1160
1161 // channel binding if necessary
1162 if ( cbind != 1 ) {
1163 iter = spec.begin() ;
1164 advance( iter, chmin ) ;
1165 vector<double>::iterator iter2 = spec.begin() ;
1166 for ( int i = 0 ; i < nchan ; i++ ) {
1167 double sum0 = 0 ;
1168 double sum1 = 0 ;
1169 for ( int j = 0 ; j < cbind ; j++ ) {
1170 sum0 += *iter ;
1171 sum1 += 1.0 ;
1172 iter++ ;
1173 }
1174 *iter2 = sum0 / sum1 ;
1175 iter2++ ;
1176 // DEBUG
1177 //cout << "NRODataset::getSpectrum() bspec[" << i << "] = " << bspec[i] << endl ;
1178 //
1179 }
1180 spec.resize( nchan ) ;
1181 }
1182
1183 // DEBUG
1184 //cout << "NRODataset::getSpectrum() end process" << endl ;
1185 //
1186
1187 return spec ;
1188}
1189
1190int NRODataset::getIndex( int irow )
1191{
1192 // DEBUG
1193 //cout << "NRODataset::getIndex() start" << endl ;
1194 //
1195 const NRODataRecord *record = getRecord( irow ) ;
1196
1197 const string str = record->ARRYT ;
1198 // DEBUG
1199 //cout << "NRODataset::getIndex() str = " << str << endl ;
1200 //
1201 int index = (int)getArrayId(str);
1202 // DEBUG
1203 //cout << "NRODataset::getIndex() irow = " << irow << "str = " << str << " index = " << index << endl ;
1204 //
1205
1206 // DEBUG
1207 //cout << "NRODataset::getIndex() end" << endl ;
1208 //
1209 return index ;
1210}
1211
1212int NRODataset::getPolarizationNum()
1213{
1214 // DEBUG
1215 //cout << "NRODataset::getPolarizationNum() start process" << endl ;
1216 //
1217 int npol = 1;
1218 Regex reRx2("(.*V|H20ch2)$");
1219 Regex reRx1("(.*H|H20ch1)$");
1220 Bool match1 = false;
1221 Bool match2 = false;
1222 for (int i = 0; i < arrayMax(); i++) {
1223 //cout << "RX[" << i << "]=" << RX[i] << endl;
1224 if (!match1) {
1225 match1 = (reRx1.match(RX[i].c_str(), RX[i].size()) != String::npos);
1226 }
1227 if (!match2) {
1228 match2 = (reRx2.match(RX[i].c_str(), RX[i].size()) != String::npos);
1229 }
1230 }
1231
1232 if (match1 && match2)
1233 npol = 2;
1234
1235 //cout << "npol = " << npol << endl;
1236
1237 // DEBUG
1238 //cout << "NRODataset::getPolarizationNum() end process" << endl ;
1239 //
1240
1241 return npol ;
1242}
1243
1244vector<double> NRODataset::getStartIntTime()
1245{
1246 vector<double> times ;
1247 for ( int i = 0 ; i < rowNum_ ; i++ ) {
1248 times.push_back( getStartIntTime( i ) ) ;
1249 }
1250 return times ;
1251}
1252
1253double NRODataset::getStartIntTime( int i )
1254{
1255 const NRODataRecord *record = getRecord( i ) ;
1256
1257 const char *t = record->LAVST ;
1258 return getMJD( t ) ;
1259}
1260
1261double NRODataset::getMJD( const char *time )
1262{
1263 // TODO: should be checked which time zone the time depends on
1264 // 2008/11/14 Takeshi Nakazato
1265 string strStartTime( time ) ;
1266 string strYear = strStartTime.substr( 0, 4 ) ;
1267 string strMonth = strStartTime.substr( 4, 2 ) ;
1268 string strDay = strStartTime.substr( 6, 2 ) ;
1269 string strHour = strStartTime.substr( 8, 2 ) ;
1270 string strMinute = strStartTime.substr( 10, 2 ) ;
1271 string strSecond = strStartTime.substr( 12, strStartTime.size() - 12 ) ;
1272 unsigned int year = atoi( strYear.c_str() ) ;
1273 unsigned int month = atoi( strMonth.c_str() ) ;
1274 unsigned int day = atoi( strDay.c_str() ) ;
1275 unsigned int hour = atoi( strHour.c_str() ) ;
1276 unsigned int minute = atoi( strMinute.c_str() ) ;
1277 double second = atof( strSecond.c_str() ) ;
1278 Time t( year, month, day, hour, minute, second ) ;
1279
1280 return t.modifiedJulianDay() ;
1281}
1282
1283double NRODataset::getScanTime( int i )
1284{
1285 double startTime = getStartIntTime( i ) ;
1286 double interval = getIPTIM() / 86400.0 ; // [sec]->[day]
1287 return startTime+0.5*interval ;
1288}
1289
1290vector<bool> NRODataset::getIFs()
1291{
1292 vector<bool> v ;
1293 vector< vector<double> > fref ;
1294 vector< vector<double> > chcal = CHCAL ;
1295 vector<double> f0cal = F0CAL ;
1296 vector<double> beres = BERES ;
1297 for ( int i = 0 ; i < rowNum_ ; i++ ) {
1298 vector<double> f( 4, 0 ) ;
1299 uInt index = getIndex( i ) ;
1300 f[0] = chcal[index][0] ;
1301 f[1] = f0cal[index] ;
1302 f[2] = beres[index] ;
1303 if ( f[0] != 0. ) {
1304 f[1] = f[1] - f[0] * f[2] ;
1305 }
1306 const NRODataRecord *record = getRecord( i ) ;
1307 f[3] = record->FREQ0 ;
1308 if ( v.size() == 0 ) {
1309 v.push_back( True ) ;
1310 fref.push_back( f ) ;
1311 }
1312 else {
1313 bool b = true ;
1314 int fsize = fref.size() ;
1315 for ( int j = 0 ; j < fsize ; j++ ) {
1316 if ( fref[j][1] == f[1] && fref[j][2] == f[2] && fref[j][3] == f[3] ) {
1317 b = false ;
1318 }
1319 }
1320 if ( b ) {
1321 v.push_back( True ) ;
1322 fref.push_back( f ) ;
1323 }
1324 }
1325 }
1326
1327
1328 // DEBUG
1329 //cout << "NRODataset::getIFs() number of IF is " << v.size() << endl ;
1330 //
1331
1332 return v ;
1333}
1334
1335vector<double> NRODataset::getFrequencies( int i )
1336{
1337 // return value
1338 // v[0] reference channel
1339 // v[1] reference frequency
1340 // v[2] frequency increment
1341 vector<double> v( 3, 0.0 ) ;
1342
1343 const NRODataRecord *record = getRecord( i ) ;
1344 string arryt = string( record->ARRYT ) ;
1345 uInt ib = getArrayId( arryt ) ;
1346 string rxname = getRX()[0] ;
1347 string key = arryt ;
1348 if ( rxname.find("MULT2") != string::npos )
1349 key = "BEARS" ;
1350
1351 if ( frec_.isDefined( key ) ) {
1352 // frequency for the array is already calculated
1353 Vector<Double> f = frec_.asArrayDouble( key ) ;
1354 Double *f_p = f.data() ;
1355 for ( int i = 0 ; i < 3 ; i++ )
1356 v[i] = (double)(f_p[i]) ;
1357 return v ;
1358 }
1359
1360 //int ia = -1 ;
1361 bool isAOS = false ;
1362 //cout << "NRODataset::getFrequencies() record->ARRYT=" << record->ARRYT << endl ;
1363 //cout << "NRODataset::getFrequencies() ib = " << ib << endl ;
1364
1365 if ( arryt[0] == 'W' || arryt[0] == 'U' || arryt[0] == 'H' )
1366 isAOS = true ;
1367
1368 Bool isUSB ;
1369 if ( record->FQIF1 > 0 )
1370 isUSB = True ; // USB
1371 else
1372 isUSB = False ; // LSB
1373
1374 int ivdef = -1 ;
1375 if ( (getVDEF()).compare( 0, 3, "RAD" ) == 0 )
1376 ivdef = 0 ;
1377 else if ( (getVDEF()).compare( 0, 3, "OPT" ) == 0 )
1378 ivdef = 1 ;
1379 // DEBUG
1380 //cout << "NRODataset::getFrequencies() ivdef = " << ivdef << endl ;
1381 //
1382 double vel = getURVEL() + record->VRAD ;
1383 double cvel = 2.99792458e8 ; // speed of light [m/s]
1384 double fq0 = record->FREQ0 ;
1385 //double fq0 = record->FQTRK ;
1386
1387 int ncal = getNFCAL()[ib] ;
1388 double cw = 0.0 ;
1389 vector<double> fqcal = getFQCAL()[ib] ;
1390 vector<double> chcal = getCHCAL()[ib] ;
1391 double f0cal = getF0CAL()[ib] ;
1392 Vector<Double> freqs( ncal, fq0-f0cal ) ;
1393
1394 double factor = vel / cvel ;
1395 if ( ivdef == 0 )
1396 factor = 1.0 / ( 1.0 - factor ) ;
1397 for ( int ii = 0 ; ii < ncal ; ii++ ) {
1398 freqs[ii] += fqcal[ii] ;
1399 if ( ivdef == 0 ) {
1400 freqs[ii] = freqs[ii] * factor + record->FQTRK * ( 1.0 - factor ) ;
1401 }
1402 else if ( ivdef == 1 ) {
1403 freqs[ii] = freqs[ii] * ( 1.0 + factor ) - record->FQTRK * factor ;
1404 }
1405
1406 //ofstream ofs("freqs.txt",ios::out|ios::app) ;
1407 //ofs << setprecision(16) ;
1408 //ofs << i << " " << record->ARRYT << " " << chcal[ii] << " " << freqs[ii] << " " << record->ISCAN << " " << fqcal[ii] << " " << f0cal << " " << record->FQTRK << " " << vel << endl ;
1409 //ofs.close() ;
1410
1411 }
1412
1413 if ( isAOS ) {
1414 // regridding
1415 while ( ncal < (int)chcal.size() ) {
1416 chcal.pop_back() ;
1417 }
1418 Vector<Double> xin( chcal ) ;
1419 Vector<Double> yin( freqs ) ;
1420 int nchan = getNUMCH() ;
1421 Vector<Double> xout( nchan ) ;
1422 indgen( xout ) ;
1423 Vector<Double> yout ;
1424 InterpolateArray1D<Double, Double>::interpolate( yout, xout, xin, yin, InterpolateArray1D<Double,Double>::cubic ) ;
1425 Double bw = abs( yout[nchan-1] - yout[0] ) ;
1426 bw += 0.5 * abs( yout[nchan-1] - yout[nchan-2] + yout[1] - yout[0] ) ;
1427 Double dz = bw / (Double) nchan ;
1428 if ( yout[0] > yout[1] )
1429 dz = - dz ;
1430 v[0] = 0 ;
1431 v[1] = yout[0] ;
1432 v[2] = dz ;
1433 }
1434 else {
1435
1436 cw = ( freqs[1] - freqs[0] )
1437 / ( chcal[1] - chcal[0] ) ;
1438
1439 if ( isUSB ) {
1440 // channel frequency inversion
1441 cw *= -1.0 ;
1442 Double tmp = freqs[1] ;
1443 freqs[1] = freqs[0] ;
1444 freqs[0] = tmp ;
1445 }
1446
1447 v[0] = chcal[0] - 1 ; // 0-base
1448 v[1] = freqs[0] ;
1449 v[2] = cw ;
1450 }
1451
1452 if ( refFreq_[ib] == 0.0 )
1453 refFreq_[ib] = v[1] ;
1454
1455 // register frequency setting to Record
1456 Vector<Double> f( v ) ;
1457 frec_.define( key, f ) ;
1458
1459 return v ;
1460}
1461
1462uInt NRODataset::getArrayId( string type )
1463{
1464 string sbeamno = type.substr( 1, type.size()-1 ) ;
1465 uInt ib = atoi( sbeamno.c_str() ) - 1 ;
1466 return ib ;
1467}
1468
1469uInt NRODataset::getSortedArrayId( string type )
1470{
1471 uInt index = 0;
1472 while (arrayNames_[index] != type && index < (uInt)ARYNM)
1473 ++index;
1474 return index;
1475}
1476
1477void NRODataset::show()
1478{
1479 LogIO os( LogOrigin( "NRODataset", "show()", WHERE ) ) ;
1480
1481 os << LogIO::NORMAL << "------------------------------------------------------------" << endl ;
1482 os << LogIO::NORMAL << "Number of scan = " << scanNum_ << endl ;
1483 os << LogIO::NORMAL << "Number of data record = " << rowNum_ << endl ;
1484 os << LogIO::NORMAL << "Length of data record = " << scanLen_ << " bytes" << endl ;
1485 os << LogIO::NORMAL << "Allocated memory for spectral data = " << dataLen_ << " bytes" << endl ;
1486 os << LogIO::NORMAL << "Max number of channel = " << chmax_ << endl ;
1487 os << LogIO::NORMAL << "------------------------------------------------------------" << endl ;
1488 os.post() ;
1489
1490}
1491
1492double NRODataset::toLSR( double v, double t, double x, double y )
1493{
1494 double vlsr ;
1495
1496 // get epoch
1497 double tcent = t + 0.5*getIPTIM()/86400.0 ;
1498 MEpoch me( Quantity( tcent, "d" ), MEpoch::UTC ) ;
1499
1500 // get antenna position
1501 MPosition mp ;
1502 if ( SITE.find( "45" ) != string::npos ) {
1503 // 45m telescope
1504 Double posx = -3.8710235e6 ;
1505 Double posy = 3.4281068e6 ;
1506 Double posz = 3.7240395e6 ;
1507 mp = MPosition( MVPosition( posx, posy, posz ),
1508 MPosition::ITRF ) ;
1509 }
1510 else {
1511 // ASTE
1512 Vector<Double> pos( 2 ) ;
1513 pos[0] = -67.7031 ;
1514 pos[1] = -22.9717 ;
1515 Double sitealt = 4800.0 ;
1516 mp = MPosition( MVPosition( Quantity( sitealt, "m" ),
1517 Quantum< Vector<Double> >( pos, "deg" ) ),
1518 MPosition::WGS84 ) ;
1519 }
1520
1521 // get direction
1522 MDirection md ;
1523 if ( SCNCD == 0 ) {
1524 // RADEC
1525 if ( EPOCH == "B1950" ) {
1526 md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
1527 MDirection::B1950 ) ;
1528 }
1529 else {
1530 md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
1531 MDirection::J2000 ) ;
1532 }
1533 }
1534 else if ( SCNCD == 1 ) {
1535 // LB
1536 md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
1537 MDirection::GALACTIC ) ;
1538 }
1539 else {
1540 // AZEL
1541 md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
1542 MDirection::AZEL ) ;
1543 }
1544
1545 // to LSR
1546 MeasFrame mf( me, mp, md ) ;
1547 MFrequency::Convert tolsr( MFrequency::TOPO, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
1548 vlsr = (double)(tolsr( Double(v) ).get( "Hz" ).getValue()) ;
1549
1550 return vlsr ;
1551}
1552
1553uInt NRODataset::getPolNo( int i )
1554{
1555 int idx = getIndex( i ) ;
1556// cout << "HORN[" << idx << "]=" << HORN[idx]
1557// << ", RX[" << idx << "]=" << RX[idx] << endl ;
1558 return polNoFromRX( RX[idx] ) ;
1559}
1560
1561uInt NRODataset::polNoFromRX( const string &rx )
1562{
1563 uInt polno = 0 ;
1564 // 2013/01/23 TN
1565 // In NRO 45m telescope, naming convension for dual-polarization
1566 // receiver is as follows:
1567 //
1568 // xxxH for horizontal component,
1569 // xxxV for vertical component.
1570 //
1571 // Exception is H20ch1/ch2.
1572 // Here, POLNO is assigned as follows:
1573 //
1574 // POLNO=0: xxxH or H20ch1
1575 // 1: xxxV or H20ch2
1576 //
1577 // For others, POLNO is always 0.
1578 String rxString(rx);
1579 rxString.trim();
1580 //cout << "rx='" << rxString << "' (size " << rxString.size() << ")" << endl;
1581 Regex reRx("(.*V|H20ch2)$");
1582 if (reRx.match(rxString.c_str(), rxString.size()) != String::npos) {
1583 //cout << "match!" << endl;
1584 polno = 1;
1585 }
1586 return polno ;
1587}
1588
1589void NRODataset::initArray()
1590{
1591 if (ARYNM <= 0)
1592 throw AipsError("ARYNM must be greater than zero.");
1593
1594 int numArray = 0;
1595 arrayNames_.resize(ARYNM);
1596 for (int irow = 0; numArray < ARYNM && irow < rowNum_; irow++) {
1597 //cout << "irow " << irow << endl;
1598 const NRODataRecord *record = getRecord( irow ) ;
1599 const string str = record->ARRYT ;
1600 if (find(arrayNames_.begin(), arrayNames_.end(), str) == arrayNames_.end()) {
1601 arrayNames_[numArray] = str;
1602 //cout << "arrayNames_[" << numArray << "]=" << str << endl;
1603 ++numArray;
1604 }
1605 }
1606 //cout << "numArray=" << numArray << endl;
1607}
1608
1609int NRODataset::getScanNum()
1610{
1611 long offset = getDataSize() + scanLen_ * NSCAN * ARYNM ;
1612 fseek( fp_, offset, SEEK_SET ) ;
1613 // try to read data
1614 fgetc( fp_ ) ;
1615 int eof = feof( fp_ ) ;
1616 //cout << "eof=" << eof << endl;
1617 // reset file pointer
1618 fseek( fp_, 0, SEEK_SET ) ;
1619 return NSCAN + (eof > 0 ? 0 : 1) ;
1620}
Note: See TracBrowser for help on using the repository browser.