source: branches/alma/src/STFiller.cpp@ 1603

Last change on this file since 1603 was 1603, checked in by TakTsutsumi, 15 years ago

New Development: No, merge with asap2.3.1

JIRA Issue: Yes CAS-1450

Ready to Release: Yes/No

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): single dish

Description: Upgrade of alma branch based on ASAP2.2.0

(rev.1562) to ASAP2.3.1 (rev.1561)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.8 KB
Line 
1//
2// C++ Implementation: STFiller
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006-2007
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <casa/iostream.h>
13#include <casa/iomanip.h>
14
15#include <casa/Exceptions.h>
16#include <casa/OS/Path.h>
17#include <casa/OS/File.h>
18#include <casa/Quanta/Unit.h>
19#include <casa/Arrays/ArrayMath.h>
20#include <casa/Arrays/ArrayLogical.h>
21#include <casa/Utilities/Regex.h>
22
23#include <casa/Containers/RecordField.h>
24
25#include <tables/Tables/TableRow.h>
26
27#include <measures/Measures/MDirection.h>
28#include <measures/Measures/MeasConvert.h>
29
30#include <atnf/PKSIO/PKSrecord.h>
31#include <atnf/PKSIO/PKSreader.h>
32#ifdef HAS_ALMA
33 #include <casa/System/ProgressMeter.h>
34#endif
35#include <casa/System/ProgressMeter.h>
36#include <atnf/PKSIO/NROReader.h>
37
38#include <time.h>
39
40
41#include "STDefs.h"
42#include "STAttr.h"
43
44#include "STFiller.h"
45#include "STHeader.h"
46
47using namespace casa;
48
49namespace asap {
50
51STFiller::STFiller() :
52 reader_(0),
53 header_(0),
54 table_(0),
55 refRx_(".*(e|w|_R)$"),
56 nreader_(0)
57{
58}
59
60STFiller::STFiller( CountedPtr< Scantable > stbl ) :
61 reader_(0),
62 header_(0),
63 table_(stbl),
64 refRx_(".*(e|w|_R)$"),
65 nreader_(0)
66{
67}
68
69STFiller::STFiller(const std::string& filename, int whichIF, int whichBeam ) :
70 reader_(0),
71 header_(0),
72 table_(0),
73 refRx_(".*(e|w|_R)$"),
74 nreader_(0)
75{
76 open(filename, whichIF, whichBeam);
77}
78
79STFiller::~STFiller()
80{
81 close();
82}
83
84void STFiller::open( const std::string& filename, int whichIF, int whichBeam, casa::Bool getPt )
85{
86 if (table_.null()) {
87 table_ = new Scantable();
88 }
89 if (reader_) { delete reader_; reader_ = 0; }
90 Bool haveBase, haveSpectra;
91
92 String inName(filename);
93 Path path(inName);
94 inName = path.expandedName();
95
96 File file(inName);
97 if ( !file.exists() ) {
98 throw(AipsError("File does not exist"));
99 }
100 filename_ = inName;
101
102 // Create reader and fill in values for arguments
103 String format;
104 Vector<Bool> beams, ifs;
105 Vector<uInt> nchans,npols;
106
107 //
108 // if isNRO_ is true, try NROReader
109 //
110 // 2008/11/11 Takeshi Nakazato
111 isNRO_ = fileCheck() ;
112 if ( isNRO_ ) {
113 if ( (nreader_ = getNROReader( inName, format )) == 0 ) {
114 throw(AipsError("Creation of NROReader failed")) ;
115 }
116 else {
117 openNRO( whichIF, whichBeam ) ;
118 return ;
119 }
120 }
121 //
122
123 if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, ifs,
124 nchans, npols, haveXPol_,haveBase, haveSpectra
125 )) == 0 ) {
126 throw(AipsError("Creation of PKSreader failed"));
127 }
128 if (!haveSpectra) {
129 delete reader_;
130 reader_ = 0;
131 throw(AipsError("No spectral data in file."));
132 return;
133 }
134 nBeam_ = beams.nelements();
135 nIF_ = ifs.nelements();
136 // Get basic parameters.
137 if ( anyEQ(haveXPol_, True) ) {
138 pushLog("Cross polarization present");
139 for (uInt i=0; i< npols.nelements();++i) {
140 if (npols[i] < 3) npols[i] += 2;// Convert Complex -> 2 Floats
141 }
142 }
143 if (header_) delete header_;
144 header_ = new STHeader();
145 header_->nchan = max(nchans);
146 header_->npol = max(npols);
147 header_->nbeam = nBeam_;
148
149 Int status = reader_->getHeader(header_->observer, header_->project,
150 header_->antennaname, header_->antennaposition,
151 header_->obstype,
152 header_->fluxunit,
153 header_->equinox,
154 header_->freqref,
155 header_->utc, header_->reffreq,
156 header_->bandwidth);
157
158 if (status) {
159 delete reader_;
160 reader_ = 0;
161 delete header_;
162 header_ = 0;
163 throw(AipsError("Failed to get header."));
164 }
165 if ((header_->obstype).matches("*SW*")) {
166 // need robust way here - probably read ahead of next timestamp
167 pushLog("Header indicates frequency switched observation.\n"
168 "setting # of IFs = 1 ");
169 nIF_ = 1;
170 header_->obstype = String("fswitch");
171 }
172 // Determine Telescope and set brightness unit
173
174
175 Bool throwIt = False;
176 Instrument inst = STAttr::convertInstrument(header_->antennaname, throwIt);
177
178 if (inst==ATMOPRA || inst==TIDBINBILLA) {
179 header_->fluxunit = "K";
180 } else {
181 // downcase for use with Quanta
182 if (header_->fluxunit == "JY") {
183 header_->fluxunit = "Jy";
184 }
185 }
186 STAttr stattr;
187 header_->poltype = stattr.feedPolType(inst);
188 header_->nif = nIF_;
189 header_->epoch = "UTC";
190 // *** header_->frequnit = "Hz"
191 // Apply selection criteria.
192 Vector<Int> ref;
193 ifOffset_ = 0;
194 if (whichIF>=0) {
195 if (whichIF>=0 && whichIF<nIF_) {
196 ifs = False;
197 ifs(whichIF) = True;
198 header_->nif = 1;
199 nIF_ = 1;
200 ifOffset_ = whichIF;
201 } else {
202 delete reader_;
203 reader_ = 0;
204 delete header_;
205 header_ = 0;
206 throw(AipsError("Illegal IF selection"));
207 }
208 }
209 beamOffset_ = 0;
210 if (whichBeam>=0) {
211 if (whichBeam>=0 && whichBeam<nBeam_) {
212 beams = False;
213 beams(whichBeam) = True;
214 header_->nbeam = 1;
215 nBeam_ = 1;
216 beamOffset_ = whichBeam;
217 } else {
218 delete reader_;
219 reader_ = 0;
220 delete header_;
221 header_ = 0;
222 throw(AipsError("Illegal Beam selection"));
223 }
224 }
225 Vector<Int> start(nIF_, 1);
226 Vector<Int> end(nIF_, 0);
227 reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False, getPt);
228 table_->setHeader(*header_);
229 //For MS, add the location of POINTING of the input MS so one get
230 //pointing data from there, if necessary.
231 //Also find nrow in MS
232 nInDataRow = 0;
233 if (format == "MS2") {
234 Path datapath(inName);
235 String ptTabPath = datapath.absoluteName();
236 Table inMS(ptTabPath);
237 nInDataRow = inMS.nrow();
238 ptTabPath.append("/POINTING");
239 table_->table().rwKeywordSet().define("POINTING", ptTabPath);
240 if ((header_->antennaname).matches("GBT")) {
241 String GOTabPath = datapath.absoluteName();
242 GOTabPath.append("/GBT_GO");
243 table_->table().rwKeywordSet().define("GBT_GO", GOTabPath);
244 }
245 }
246 String freqFrame = header_->freqref;
247 //translate frequency reference frame back to
248 //MS style (as PKSMS2reader converts the original frame
249 //in FITS standard style)
250 if (freqFrame == "TOPOCENT") {
251 freqFrame = "TOPO";
252 } else if (freqFrame == "GEOCENER") {
253 freqFrame = "GEO";
254 } else if (freqFrame == "BARYCENT") {
255 freqFrame = "BARY";
256 } else if (freqFrame == "GALACTOC") {
257 freqFrame = "GALACTO";
258 } else if (freqFrame == "LOCALGRP") {
259 freqFrame = "LGROUP";
260 } else if (freqFrame == "CMBDIPOL") {
261 freqFrame = "CMB";
262 } else if (freqFrame == "SOURCE") {
263 freqFrame = "REST";
264 }
265 table_->frequencies().setFrame(freqFrame);
266
267}
268
269void STFiller::close( )
270{
271 delete reader_;reader_=0;
272 delete nreader_;nreader_=0;
273 delete header_;header_=0;
274 table_ = 0;
275}
276
277int asap::STFiller::read( )
278{
279 int status = 0;
280
281 //
282 // for NRO data
283 //
284 // 2008/11/12 Takeshi Nakazato
285 if ( isNRO_ ) {
286 status = readNRO() ;
287 return status ;
288 }
289 //
290
291/**
292 Int beamNo, IFno, refBeam, scanNo, cycleNo;
293 Float azimuth, elevation, focusAxi, focusRot, focusTan,
294 humidity, parAngle, pressure, temperature, windAz, windSpeed;
295 Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
296 String fieldName, srcName, tcalTime, obsType;
297 Vector<Float> calFctr, sigma, tcal, tsys;
298 Matrix<Float> baseLin, baseSub;
299 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq(1);
300 Matrix<Float> spectra;
301 Matrix<uChar> flagtra;
302 Complex xCalFctr;
303 Vector<Complex> xPol;
304**/
305
306 Double min = 0.0;
307 Double max = nInDataRow;
308#ifdef HAS_ALMA
309 ProgressMeter fillpm(min, max, "Data importing progress");
310#endif
311 PKSrecord pksrec;
312 int n = 0;
313 while ( status == 0 ) {
314 status = reader_->read(pksrec);
315 if ( status != 0 ) break;
316 n += 1;
317
318 Regex filterrx(".*[SL|PA]$");
319 Regex obsrx("^AT.+");
320 if ( header_->antennaname.matches(obsrx) &&
321 pksrec.obsType.matches(filterrx)) {
322 //cerr << "ignoring paddle scan" << endl;
323 continue;
324 }
325 TableRow row(table_->table());
326 TableRecord& rec = row.record();
327 // fields that don't get used and are just passed through asap
328 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
329 // MRC changed type from double to float
330 Vector<Double> sratedbl(pksrec.scanRate.nelements());
331 convertArray(sratedbl, pksrec.scanRate);
332 *srateCol = sratedbl;
333 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
334 *spmCol = pksrec.srcPM;
335 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
336 *sdirCol = pksrec.srcDir;
337 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
338 *svelCol = pksrec.srcVel;
339 // the real stuff
340 RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
341 *fitCol = -1;
342 RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
343 *scanoCol = pksrec.scanNo-1;
344 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
345 *cyclenoCol = pksrec.cycleNo-1;
346 RecordFieldPtr<Double> mjdCol(rec, "TIME");
347 *mjdCol = pksrec.mjd;
348 RecordFieldPtr<Double> intCol(rec, "INTERVAL");
349 *intCol = pksrec.interval;
350 RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
351 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
352 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
353 *fieldnCol = pksrec.fieldName;
354 // try to auto-identify if it is on or off.
355 Regex rx(refRx_);
356 Regex rx2("_S$");
357 Int match = pksrec.srcName.matches(rx);
358 if (match) {
359 *srcnCol = pksrec.srcName;
360 } else {
361 *srcnCol = pksrec.srcName.before(rx2);
362 }
363 //*srcnCol = pksrec.srcName;//.before(rx2);
364 *srctCol = match;
365 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
366 *beamCol = pksrec.beamNo-beamOffset_-1;
367 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
368 Int rb = -1;
369 if (nBeam_ > 1 ) rb = pksrec.refBeam-1;
370 *rbCol = rb;
371 RecordFieldPtr<uInt> ifCol(rec, "IFNO");
372 *ifCol = pksrec.IFno-ifOffset_- 1;
373 uInt id;
374 /// @todo this has to change when nchan isn't global anymore
375 id = table_->frequencies().addEntry(Double(header_->nchan/2),
376 pksrec.refFreq, pksrec.freqInc);
377 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
378 *mfreqidCol = id;
379
380 id = table_->molecules().addEntry(pksrec.restFreq);
381 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
382 *molidCol = id;
383
384 id = table_->tcal().addEntry(pksrec.tcalTime, pksrec.tcal);
385 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
386 *mcalidCol = id;
387 id = table_->weather().addEntry(pksrec.temperature, pksrec.pressure,
388 pksrec.humidity, pksrec.windSpeed,
389 pksrec.windAz);
390 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
391 *mweatheridCol = id;
392 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
393 id = table_->focus().addEntry(pksrec.focusAxi, pksrec.focusTan,
394 pksrec.focusRot);
395 *mfocusidCol = id;
396 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
397 *dirCol = pksrec.direction;
398 RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
399 *azCol = pksrec.azimuth;
400 RecordFieldPtr<Float> elCol(rec, "ELEVATION");
401 *elCol = pksrec.elevation;
402
403 RecordFieldPtr<Float> parCol(rec, "PARANGLE");
404 *parCol = pksrec.parAngle;
405
406 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
407 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
408 RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
409
410 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
411 // Turn the (nchan,npol) matrix and possible complex xPol vector
412 // into 2-4 rows in the scantable
413 Vector<Float> tsysvec(1);
414 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
415 uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
416 for ( uInt i=0; i< npol; ++i ) {
417 tsysvec = pksrec.tsys(i);
418 *tsysCol = tsysvec;
419 *polnoCol = i;
420
421 *specCol = pksrec.spectra.column(i);
422 *flagCol = pksrec.flagged.column(i);
423 table_->table().addRow();
424 row.put(table_->table().nrow()-1, rec);
425 }
426 if ( haveXPol_[0] ) {
427 // no tsys given for xpol, so emulate it
428 tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
429 *tsysCol = tsysvec;
430 // add real part of cross pol
431 *polnoCol = 2;
432 Vector<Float> r(real(pksrec.xPol));
433 *specCol = r;
434 // make up flags from linears
435 /// @fixme this has to be a bitwise or of both pols
436 *flagCol = pksrec.flagged.column(0);// | pksrec.flagged.column(1);
437 table_->table().addRow();
438 row.put(table_->table().nrow()-1, rec);
439 // ad imaginary part of cross pol
440 *polnoCol = 3;
441 Vector<Float> im(imag(pksrec.xPol));
442 *specCol = im;
443 table_->table().addRow();
444 row.put(table_->table().nrow()-1, rec);
445 }
446#ifdef HAS_ALMA
447 fillpm._update(n);
448#endif
449 }
450 if (status > 0) {
451 close();
452 throw(AipsError("Reading error occured, data possibly corrupted."));
453 }
454#ifdef HAS_ALMA
455 fillpm.done();
456#endif
457 return status;
458}
459
460/**
461 * For NRO data
462 *
463 * 2008/11/11 Takeshi Nakazato
464 **/
465void STFiller::openNRO( int whichIF, int whichBeam )
466{
467 // open file
468 // DEBUG
469 time_t t0 ;
470 time( &t0 ) ;
471 tm *ttm = localtime( &t0 ) ;
472
473 cout << "STFiller::openNRO() Start time = " << t0
474 << " ("
475 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
476 << " "
477 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
478 << ")" << endl ;
479
480 // fill STHeader
481 header_ = new STHeader() ;
482
483 if ( nreader_->getHeaderInfo( header_->nchan,
484 header_->npol,
485 nIF_,
486 nBeam_,
487 header_->observer,
488 header_->project,
489 header_->obstype,
490 header_->antennaname,
491 header_->antennaposition,
492 header_->equinox,
493 header_->freqref,
494 header_->reffreq,
495 header_->bandwidth,
496 header_->utc,
497 header_->fluxunit,
498 header_->epoch,
499 header_->poltype ) ) {
500 cout << "STFiller::openNRO() Failed to get header information." << endl ;
501 return ;
502 }
503
504 // set frame keyword of FREQUENCIES table
505 if ( header_->freqref != "TOPO" ) {
506 table_->frequencies().setFrame( header_->freqref, false ) ;
507 }
508
509 ifOffset_ = 0;
510 vector<Bool> ifs = nreader_->getIFs() ;
511 if ( whichIF >= 0 ) {
512 if ( whichIF >= 0 && whichIF < nIF_ ) {
513 for ( int i = 0 ; i < nIF_ ; i++ )
514 ifs[i] = False ;
515 ifs[whichIF] = True ;
516 header_->nif = 1;
517 nIF_ = 1;
518 ifOffset_ = whichIF;
519 } else {
520 delete reader_;
521 reader_ = 0;
522 delete header_;
523 header_ = 0;
524 throw(AipsError("Illegal IF selection"));
525 }
526 }
527
528 // DEBUG
529 //cout << "STFiller::openNRO() nIF " << endl ;
530 //
531
532 beamOffset_ = 0;
533 vector<Bool> beams = nreader_->getBeams() ;
534 if (whichBeam>=0) {
535 if (whichBeam>=0 && whichBeam<nBeam_) {
536 for ( int i = 0 ; i < nBeam_ ; i++ )
537 beams[i] = False ;
538 beams[whichBeam] = True;
539 header_->nbeam = 1;
540 nBeam_ = 1;
541 beamOffset_ = whichBeam;
542 } else {
543 delete reader_;
544 reader_ = 0;
545 delete header_;
546 header_ = 0;
547 throw(AipsError("Illegal Beam selection"));
548 }
549 }
550
551 // DEBUG
552 //cout << "STFiller::openNRO() nBeam " << endl ;
553 //
554
555 header_->nbeam = nBeam_ ;
556 header_->nif = nIF_ ;
557
558 // set header
559 table_->setHeader( *header_ ) ;
560
561 // DEBUG
562 //cout << "STFiller::openNRO() Velocity Definition = " << nheader->getVDEF() << endl ;
563
564 // DEBUG
565 time_t t1 ;
566 time( &t1 ) ;
567 ttm = localtime( &t1 ) ;
568 cout << "STFiller::openNRO() End time = " << t1
569 << " ("
570 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
571 << " "
572 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
573 << ")" << endl ;
574 cout << "STFiller::openNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
575 //
576
577 return ;
578}
579
580int STFiller::readNRO()
581{
582 // DEBUG
583 time_t t0 ;
584 time( &t0 ) ;
585 tm *ttm = localtime( &t0 ) ;
586 cout << "STFiller::readNRO() Start time = " << t0
587 << " ("
588 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
589 << " "
590 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
591 << ")" << endl ;
592 //
593
594 // fill row
595 uInt id ;
596 uInt imax = nreader_->getRowNum() ;
597 vector< vector<double > > freqs ;
598 uInt i = 0 ;
599 int count = 0 ;
600 uInt scanno ;
601 uInt cycleno ;
602 uInt beamno ;
603 uInt polno ;
604 vector<double> fqs ;
605 Vector<Double> restfreq ;
606 uInt refbeamno ;
607 Double scantime ;
608 Double interval ;
609 String srcname ;
610 String fieldname ;
611 Array<Float> spectra ;
612 Array<uChar> flagtra ;
613 Array<Float> tsys ;
614 Array<Double> direction ;
615 Float azimuth ;
616 Float elevation ;
617 Float parangle ;
618 Float opacity ;
619 uInt tcalid ;
620 Int fitid ;
621 uInt focusid ;
622 Float temperature ;
623 Float pressure ;
624 Float humidity ;
625 Float windvel ;
626 Float winddir ;
627 Double srcvel ;
628 Array<Double> propermotion ;
629 Vector<Double> srcdir ;
630 Array<Double> scanrate ;
631 for ( i = 0 ; i < imax ; i++ ) {
632// if( nreader_->getDataset()->getRecord( i ) == NULL ) {
633// cerr << "STFiller::readNRO() error while reading row " << i << endl ;
634// return -1 ;
635// }
636
637 string scanType = nreader_->getScanType( i ) ;
638 Int srcType = -1 ;
639 if ( scanType.compare( 0, 2, "ON") == 0 ) {
640 // cout << "ON srcType: " << i << endl ;
641 srcType = 0 ;
642 }
643 else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
644 //cout << "OFF srcType: " << i << endl ;
645 srcType = 1 ;
646 }
647 else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
648 //cout << "ZERO srcType: " << i << endl ;
649 srcType = 2 ;
650 }
651 else {
652 //cout << "Undefined srcType: " << i << endl ;
653 srcType = 3 ;
654 }
655
656 // if srcType is 2 (ZERO scan), ignore scan
657 if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
658 TableRow row( table_->table() ) ;
659 TableRecord& rec = row.record();
660
661 if ( nreader_->getScanInfo( i,
662 scanno,
663 cycleno,
664 beamno,
665 polno,
666 fqs,
667 restfreq,
668 refbeamno,
669 scantime,
670 interval,
671 srcname,
672 fieldname,
673 spectra,
674 flagtra,
675 tsys,
676 direction,
677 azimuth,
678 elevation,
679 parangle,
680 opacity,
681 tcalid,
682 fitid,
683 focusid,
684 temperature,
685 pressure,
686 humidity,
687 windvel,
688 winddir,
689 srcvel,
690 propermotion,
691 srcdir,
692 scanrate ) ) {
693 cerr << "STFiller::readNRO() Failed to get scan information." << endl ;
694 return 1 ;
695 }
696
697 RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
698 *scannoCol = scanno ;
699 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ;
700 *cyclenoCol = cycleno ;
701 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
702 *beamCol = beamno ;
703 RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
704 RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
705 *polnoCol = polno ;
706 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
707 if ( freqs.size() == 0 ) {
708 id = table_->frequencies().addEntry( Double( fqs[0] ),
709 Double( fqs[1] ),
710 Double( fqs[2] ) ) ;
711 *mfreqidCol = id ;
712 *ifCol = id ;
713 freqs.push_back( fqs ) ;
714 }
715 else {
716 int iadd = -1 ;
717 for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
718 //cout << "STFiller::readNRO() freqs[" << iif << "][1] = " << freqs[iif][1] << endl ;
719 double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
720 //cout << "STFiller::readNRO() fdiff = " << fdiff << endl ;
721 if ( fdiff < 1.0e-8 ) {
722 iadd = iif ;
723 break ;
724 }
725 }
726 if ( iadd == -1 ) {
727 id = table_->frequencies().addEntry( Double( fqs[0] ),
728 Double( fqs[1] ),
729 Double( fqs[2] ) ) ;
730 *mfreqidCol = id ;
731 *ifCol = id ;
732 freqs.push_back( fqs ) ;
733 }
734 else {
735 *mfreqidCol = iadd ;
736 *ifCol = iadd ;
737 }
738 }
739 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
740 id = table_->molecules().addEntry( restfreq ) ;
741 *molidCol = id ;
742 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
743 *rbCol = refbeamno ;
744 RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
745 *mjdCol = scantime ;
746 RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
747 *intervalCol = interval ;
748 RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
749 *srcnCol = srcname ;
750 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
751 *srctCol = srcType ;
752 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
753 *fieldnCol = fieldname ;
754 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
755 *specCol = spectra ;
756 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
757 *flagCol = flagtra ;
758 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
759 *tsysCol = tsys ;
760 RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
761 *dirCol = direction ;
762 RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
763 *azCol = azimuth ;
764 RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
765 *elCol = elevation ;
766 RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
767 *parCol = parangle ;
768 RecordFieldPtr<Float> tauCol(rec, "OPACITY") ;
769 *tauCol = opacity ;
770 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
771 *mcalidCol = tcalid ;
772 RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
773 *fitCol = fitid ;
774 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
775 *mfocusidCol = focusid ;
776 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
777 id = table_->weather().addEntry( temperature,
778 pressure,
779 humidity,
780 windvel,
781 winddir ) ;
782 *mweatheridCol = id ;
783 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
784 *svelCol = srcvel ;
785 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
786 *spmCol = propermotion ;
787 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
788 *sdirCol = srcdir ;
789 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
790 *srateCol = scanrate ;
791
792 table_->table().addRow() ;
793 row.put(table_->table().nrow()-1, rec) ;
794 }
795 else {
796 count++ ;
797 }
798 // DEBUG
799 //int rownum = nreader_->getRowNum() ;
800 //cout << "STFiller::readNRO() Finished row " << i << "/" << rownum << endl ;
801 //
802 }
803
804 // DEBUG
805 time_t t1 ;
806 time( &t1 ) ;
807 ttm = localtime( &t1 ) ;
808 cout << "STFiller::readNRO() Processed " << i << " rows" << endl ;
809 cout << "STFiller::readNRO() Added " << i - count << " rows (ignored "
810 << count << " \"ZERO\" scans)" << endl ;
811 cout << "STFiller::readNRO() End time = " << t1
812 << " ("
813 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
814 << " "
815 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
816 << ")" << endl ;
817 cout << "STFiller::readNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
818 //
819
820 return 0 ;
821}
822
823Bool STFiller::fileCheck()
824{
825 bool bval = false ;
826
827 // if filename_ is directory, return false
828 File inFile( filename_ ) ;
829 if ( inFile.isDirectory() )
830 return bval ;
831
832 // if beginning of header data is "RW", return true
833 // otherwise, return false ;
834 FILE *fp = fopen( filename_.c_str(), "r" ) ;
835 char buf[9] ;
836 char buf2[80] ;
837 fread( buf, 4, 1, fp ) ;
838 buf[4] = '\0' ;
839 fseek( fp, 640, SEEK_SET ) ;
840 fread( buf2, 80, 1, fp ) ;
841 if ( ( strncmp( buf, "RW", 2 ) == 0 ) || ( strstr( buf2, "NRO45M" ) != NULL ) ) {
842 bval = true ;
843 }
844 fclose( fp ) ;
845 return bval ;
846}
847
848}//namespace asap
Note: See TracBrowser for help on using the repository browser.