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

Last change on this file since 1661 was 1657, checked in by WataruKawasaki, 15 years ago

New Development: Yes

JIRA Issue: Yes (CAS-1433)

Ready to Release: Yes

Interface Changes: No

Put in Release Notes: No

Module(s): sd.save()

Description: Modified to read/write row-based flagging info.


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