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

Last change on this file since 1544 was 1533, 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_FILENAME',False)

Put in Release Notes: No

Module(s): atnf

Description:

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