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

Last change on this file since 3112 was 3111, checked in by Kana Sugimoto, 8 years ago

New Development: No

JIRA Issue: Yes (CAS-9502)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sd.scantable, sdsave

Description: Cleared up commented codes.


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