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

Last change on this file since 1520 was 1519, 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: asap_init()

s=sd.scantable('NRO_DATA_FILE',False)

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Get all necessary informations for ASAP table

using NROReader::getScanInfo() method.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.1 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 // DEBUG
483 //cout << "STFiller::openNRO() getHeaderInfo " << endl ;
484 //
485
486 ifOffset_ = 0;
487 vector<Bool> ifs = nreader_->getIFs() ;
488 if ( whichIF >= 0 ) {
489 if ( whichIF >= 0 && whichIF < nIF_ ) {
490 for ( int i = 0 ; i < nIF_ ; i++ )
491 ifs[i] = False ;
492 ifs[whichIF] = True ;
493 header_->nif = 1;
494 nIF_ = 1;
495 ifOffset_ = whichIF;
496 } else {
497 delete reader_;
498 reader_ = 0;
499 delete header_;
500 header_ = 0;
501 throw(AipsError("Illegal IF selection"));
502 }
503 }
504
505 // DEBUG
506 //cout << "STFiller::openNRO() nIF " << endl ;
507 //
508
509 beamOffset_ = 0;
510 vector<Bool> beams = nreader_->getBeams() ;
511 if (whichBeam>=0) {
512 if (whichBeam>=0 && whichBeam<nBeam_) {
513 for ( int i = 0 ; i < nBeam_ ; i++ )
514 beams[i] = False ;
515 beams[whichBeam] = True;
516 header_->nbeam = 1;
517 nBeam_ = 1;
518 beamOffset_ = whichBeam;
519 } else {
520 delete reader_;
521 reader_ = 0;
522 delete header_;
523 header_ = 0;
524 throw(AipsError("Illegal Beam selection"));
525 }
526 }
527
528 // DEBUG
529 //cout << "STFiller::openNRO() nBeam " << endl ;
530 //
531
532 header_->nbeam = nBeam_ ;
533 header_->nif = nIF_ ;
534
535 // set header
536 table_->setHeader( *header_ ) ;
537
538 // DEBUG
539 //cout << "STFiller::openNRO() Velocity Definition = " << nheader->getVDEF() << endl ;
540
541 // DEBUG
542 time_t t1 ;
543 time( &t1 ) ;
544 ttm = localtime( &t1 ) ;
545 cout << "STFiller::openNRO() End time = " << t1
546 << " ("
547 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
548 << " "
549 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
550 << ")" << endl ;
551 cout << "STFiller::openNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
552 //
553
554 return ;
555}
556
557int STFiller::readNRO()
558{
559 // DEBUG
560 time_t t0 ;
561 time( &t0 ) ;
562 tm *ttm = localtime( &t0 ) ;
563 cout << "STFiller::readNRO() Start time = " << t0
564 << " ("
565 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
566 << " "
567 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
568 << ")" << endl ;
569 //
570
571 // fill row
572 uInt id ;
573 uInt imax = nreader_->getRowNum() ;
574 vector< vector<double > > freqs ;
575 uInt i = 0 ;
576 int count = 0 ;
577 uInt scanno ;
578 uInt cycleno ;
579 uInt beamno ;
580 uInt polno ;
581 vector<double> fqs ;
582 Vector<Double> restfreq ;
583 uInt refbeamno ;
584 Double scantime ;
585 Double interval ;
586 String srcname ;
587 String fieldname ;
588 Array<Float> spectra ;
589 Array<uChar> flagtra ;
590 Array<Float> tsys ;
591 Array<Double> direction ;
592 Float azimuth ;
593 Float elevation ;
594 Float parangle ;
595 Float opacity ;
596 uInt tcalid ;
597 Int fitid ;
598 uInt focusid ;
599 Float temperature ;
600 Float pressure ;
601 Float humidity ;
602 Float windvel ;
603 Float winddir ;
604 Double srcvel ;
605 Array<Double> propermotion ;
606 Vector<Double> srcdir ;
607 Array<Double> scanrate ;
608 for ( i = 0 ; i < imax ; i++ ) {
609 if( nreader_->getData( i ) != 0 ) {
610 cerr << "STFiller::readNRO() error while reading row " << i << endl ;
611 return -1 ;
612 }
613
614 string scanType = nreader_->getScanType( i ) ;
615 Int srcType = -1 ;
616 if ( scanType.compare( 0, 2, "ON") == 0 ) {
617 srcType = 0 ;
618 }
619 else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
620 //cout << "OFF srcType: " << i << endl ;
621 srcType = 1 ;
622 }
623 else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
624 //cout << "ZERO srcType: " << i << endl ;
625 srcType = 2 ;
626 }
627 else {
628 //cout << "Undefined srcType: " << i << endl ;
629 srcType = 3 ;
630 }
631
632 // if srcType is 2 (ZERO scan), ignore scan
633 if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
634 TableRow row( table_->table() ) ;
635 TableRecord& rec = row.record();
636
637 if ( nreader_->getScanInfo( i,
638 scanno,
639 cycleno,
640 beamno,
641 polno,
642 fqs,
643 restfreq,
644 refbeamno,
645 scantime,
646 interval,
647 srcname,
648 fieldname,
649 spectra,
650 flagtra,
651 tsys,
652 direction,
653 azimuth,
654 elevation,
655 parangle,
656 opacity,
657 tcalid,
658 fitid,
659 focusid,
660 temperature,
661 pressure,
662 humidity,
663 windvel,
664 winddir,
665 srcvel,
666 propermotion,
667 srcdir,
668 scanrate ) ) {
669 cerr << "STFiller::readNRO() Failed to get scan information." << endl ;
670 return 1 ;
671 }
672
673 RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
674 *scannoCol = scanno ;
675 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ;
676 *cyclenoCol = cycleno ;
677 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
678 *beamCol = beamno ;
679 RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
680 RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
681 *polnoCol = polno ;
682 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
683 if ( freqs.size() == 0 ) {
684 id = table_->frequencies().addEntry( Double( fqs[0] ),
685 Double( fqs[1] ),
686 Double( fqs[2] ) ) ;
687 *mfreqidCol = id ;
688 *ifCol = id ;
689 freqs.push_back( fqs ) ;
690 }
691 else {
692 int iadd = -1 ;
693 for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
694 //cout << "STFiller::readNRO() freqs[" << iif << "][1] = " << freqs[iif][1] << endl ;
695 double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
696 //cout << "STFiller::readNRO() fdiff = " << fdiff << endl ;
697 if ( fdiff < 1.0e-8 ) {
698 iadd = iif ;
699 break ;
700 }
701 }
702 if ( iadd == -1 ) {
703 id = table_->frequencies().addEntry( Double( fqs[0] ),
704 Double( fqs[1] ),
705 Double( fqs[2] ) ) ;
706 *mfreqidCol = id ;
707 *ifCol = id ;
708 freqs.push_back( fqs ) ;
709 }
710 else {
711 *mfreqidCol = iadd ;
712 *ifCol = iadd ;
713 }
714 }
715 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
716 id = table_->molecules().addEntry( restfreq ) ;
717 *molidCol = id ;
718 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
719 *rbCol = refbeamno ;
720 RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
721 *mjdCol = scantime ;
722 RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
723 *intervalCol = interval ;
724 RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
725 *srcnCol = srcname ;
726 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
727 *srctCol = srcType ;
728 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
729 *fieldnCol = fieldname ;
730 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
731 *specCol = spectra ;
732 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
733 *flagCol = flagtra ;
734 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
735 *tsysCol = tsys ;
736 RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
737 *dirCol = direction ;
738 RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
739 *azCol = azimuth ;
740 RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
741 *elCol = elevation ;
742 RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
743 *parCol = parangle ;
744 RecordFieldPtr<Float> tauCol(rec, "OPACITY") ;
745 *tauCol = opacity ;
746 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
747 *mcalidCol = tcalid ;
748 RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
749 *fitCol = fitid ;
750 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
751 *mfocusidCol = focusid ;
752 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
753 id = table_->weather().addEntry( temperature,
754 pressure,
755 humidity,
756 windvel,
757 winddir ) ;
758 *mweatheridCol = id ;
759 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
760 *svelCol = srcvel ;
761 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
762 *spmCol = propermotion ;
763 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
764 *sdirCol = srcdir ;
765 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
766 *srateCol = scanrate ;
767
768 table_->table().addRow() ;
769 row.put(table_->table().nrow()-1, rec) ;
770 }
771 else {
772 count++ ;
773 }
774 // DEBUG
775 //int rownum = nreader_->getRowNum() ;
776 //cout << "STFiller::readNRO() Finished row " << i << "/" << rownum << endl ;
777 //
778 }
779
780 // DEBUG
781 time_t t1 ;
782 time( &t1 ) ;
783 ttm = localtime( &t1 ) ;
784 cout << "STFiller::readNRO() Processed " << i << " rows" << endl ;
785 cout << "STFiller::readNRO() Added " << i - count << " rows (ignored "
786 << count << " \"ZERO\" scans)" << endl ;
787 cout << "STFiller::readNRO() End time = " << t1
788 << " ("
789 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
790 << " "
791 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
792 << ")" << endl ;
793 cout << "STFiller::readNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
794 //
795
796 return 0 ;
797}
798
799Bool STFiller::fileCheck()
800{
801 // if filename_ is directory, return false
802 File inFile( filename_ ) ;
803 if ( inFile.isDirectory() )
804 return false ;
805
806 // if beginning of header data is "RW", return true
807 // otherwise, return false ;
808 FILE *fp = fopen( filename_.c_str(), "r" ) ;
809 char buf[9] ;
810 fread( buf, 4, 1, fp ) ;
811 buf[4] = '\0' ;
812 if ( ( strncmp( buf, "RW", 2 ) == 0 ) )
813 return true ;
814
815 return false ;
816}
817
818}//namespace asap
Note: See TracBrowser for help on using the repository browser.