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

Last change on this file since 1524 was 1521, checked in by Takeshi Nakazato, 16 years ago

New Development: No

JIRA Issue: Yes CAS-1043

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: No

Module(s): atnf

Description:

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