source: trunk/src/STFiller.cpp@ 2176

Last change on this file since 2176 was 2163, checked in by Malte Marquarding, 14 years ago

Remove various compiler warnings
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.0 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, const std::string& antenna, 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, antenna, 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 // set both "FRAME" and "BASEFRAME"
267 table_->frequencies().setFrame(freqFrame, false);
268 table_->frequencies().setFrame(freqFrame,true);
269 //table_->focus().setParallactify(true);
270}
271
272void STFiller::close( )
273{
274 delete reader_;reader_=0;
275 delete nreader_;nreader_=0;
276 delete header_;header_=0;
277 table_ = 0;
278}
279
280int asap::STFiller::read( )
281{
282 int status = 0;
283
284 //
285 // for NRO data
286 //
287 // 2008/11/12 Takeshi Nakazato
288 if ( isNRO_ ) {
289 status = readNRO() ;
290 return status ;
291 }
292 //
293
294/**
295 Int beamNo, IFno, refBeam, scanNo, cycleNo;
296 Float azimuth, elevation, focusAxi, focusRot, focusTan,
297 humidity, parAngle, pressure, temperature, windAz, windSpeed;
298 Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
299 String fieldName, srcName, tcalTime, obsType;
300 Vector<Float> calFctr, sigma, tcal, tsys;
301 Matrix<Float> baseLin, baseSub;
302 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq(1);
303 Matrix<Float> spectra;
304 Matrix<uChar> flagtra;
305 Complex xCalFctr;
306 Vector<Complex> xPol;
307**/
308
309#ifdef HAS_ALMA
310 Double min = 0.0;
311 Double max = nInDataRow;
312 ProgressMeter fillpm(min, max, "Data importing progress");
313#endif
314 PKSrecord pksrec;
315 pksrec.srcType=-1;
316 int n = 0;
317 bool isGBTFITS = false ;
318 if ((header_->antennaname.find( "GBT" ) != String::npos) && File(filename_).isRegular()) {
319 FILE *fp = fopen( filename_.c_str(), "r" ) ;
320 fseek( fp, 640, SEEK_SET ) ;
321 char buf[81] ;
322 size_t tmpret = fread( buf, 80, 1, fp ) ;
323 (void) tmpret;
324 buf[80] = '\0' ;
325 if ( strstr( buf, "NRAO_GBT" ) != NULL ) {
326 isGBTFITS = true ;
327 }
328 fclose( fp ) ;
329 }
330 while ( status == 0 ) {
331 status = reader_->read(pksrec);
332 if ( status != 0 ) break;
333 n += 1;
334 Regex filterrx(".*[SL|PA]$");
335 Regex obsrx("^AT.+");
336 if ( header_->antennaname.matches(obsrx) &&
337 pksrec.obsType.matches(filterrx)) {
338 //cerr << "ignoring paddle scan" << endl;
339 continue;
340 }
341 TableRow row(table_->table());
342 TableRecord& rec = row.record();
343 // fields that don't get used and are just passed through asap
344 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
345 // MRC changed type from double to float
346 Vector<Double> sratedbl(pksrec.scanRate.nelements());
347 convertArray(sratedbl, pksrec.scanRate);
348 *srateCol = sratedbl;
349 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
350 *spmCol = pksrec.srcPM;
351 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
352 *sdirCol = pksrec.srcDir;
353 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
354 *svelCol = pksrec.srcVel;
355 // the real stuff
356 RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
357 *fitCol = -1;
358 RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
359 *scanoCol = pksrec.scanNo-1;
360 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
361 *cyclenoCol = pksrec.cycleNo-1;
362 RecordFieldPtr<Double> mjdCol(rec, "TIME");
363 *mjdCol = pksrec.mjd;
364 RecordFieldPtr<Double> intCol(rec, "INTERVAL");
365 *intCol = pksrec.interval;
366 RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
367 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
368 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
369 *fieldnCol = pksrec.fieldName;
370 // try to auto-identify if it is on or off.
371 Regex rx(refRx_);
372 Regex rx2("_S$");
373 Int match = pksrec.srcName.matches(rx);
374 if (match) {
375 *srcnCol = pksrec.srcName;
376 } else {
377 *srcnCol = pksrec.srcName.before(rx2);
378 }
379 //*srcnCol = pksrec.srcName;//.before(rx2);
380 *srctCol = match;
381 if ( pksrec.srcType != -1 ) {
382 *srctCol = pksrec.srcType ;
383 }
384 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
385 *beamCol = pksrec.beamNo-beamOffset_-1;
386 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
387 Int rb = -1;
388 if (nBeam_ > 1 ) rb = pksrec.refBeam-1;
389 *rbCol = rb;
390 RecordFieldPtr<uInt> ifCol(rec, "IFNO");
391 *ifCol = pksrec.IFno-ifOffset_- 1;
392 uInt id = table_->frequencies().addEntry(Double(pksrec.spectra.nrow()/2),
393 pksrec.refFreq, pksrec.freqInc);
394 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
395 *mfreqidCol = id;
396 //*ifCol = id;
397
398 id = table_->molecules().addEntry(pksrec.restFreq);
399 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
400 *molidCol = id;
401
402 id = table_->tcal().addEntry(pksrec.tcalTime, pksrec.tcal);
403 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
404 *mcalidCol = id;
405 id = table_->weather().addEntry(pksrec.temperature, pksrec.pressure,
406 pksrec.humidity, pksrec.windSpeed,
407 pksrec.windAz);
408 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
409 *mweatheridCol = id;
410
411 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
412 id = table_->focus().addEntry(pksrec.parAngle, pksrec.focusAxi,
413 pksrec.focusTan, pksrec.focusRot);
414 *mfocusidCol = id;
415 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
416 *dirCol = pksrec.direction;
417 RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
418 *azCol = pksrec.azimuth;
419 RecordFieldPtr<Float> elCol(rec, "ELEVATION");
420 *elCol = pksrec.elevation;
421
422 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
423 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
424 RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
425
426 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
427 // Turn the (nchan,npol) matrix and possible complex xPol vector
428 // into 2-4 rows in the scantable
429 Vector<Float> tsysvec(1);
430 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
431 uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
432 for ( uInt i=0; i< npol; ++i ) {
433 tsysvec = pksrec.tsys(i);
434 *tsysCol = tsysvec;
435 if (isGBTFITS)
436 *polnoCol = pksrec.polNo ;
437 else
438 *polnoCol = i;
439
440 *specCol = pksrec.spectra.column(i);
441 *flagCol = pksrec.flagged.column(i);
442 table_->table().addRow();
443 row.put(table_->table().nrow()-1, rec);
444 }
445
446 RecordFieldPtr< uInt > flagrowCol(rec, "FLAGROW");
447 *flagrowCol = pksrec.flagrow;
448
449 if ( haveXPol_[0] ) {
450 // no tsys given for xpol, so emulate it
451 tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
452 *tsysCol = tsysvec;
453 // add real part of cross pol
454 *polnoCol = 2;
455 Vector<Float> r(real(pksrec.xPol));
456 *specCol = r;
457 // make up flags from linears
458 /// @fixme this has to be a bitwise or of both pols
459 *flagCol = pksrec.flagged.column(0);// | pksrec.flagged.column(1);
460 table_->table().addRow();
461 row.put(table_->table().nrow()-1, rec);
462 // ad imaginary part of cross pol
463 *polnoCol = 3;
464 Vector<Float> im(imag(pksrec.xPol));
465 *specCol = im;
466 table_->table().addRow();
467 row.put(table_->table().nrow()-1, rec);
468 }
469#ifdef HAS_ALMA
470 fillpm._update(n);
471#endif
472 }
473 if (status > 0) {
474 close();
475 throw(AipsError("Reading error occured, data possibly corrupted."));
476 }
477#ifdef HAS_ALMA
478 fillpm.done();
479#endif
480 return status;
481}
482
483/**
484 * For NRO data
485 *
486 * 2008/11/11 Takeshi Nakazato
487 **/
488void STFiller::openNRO( int whichIF, int whichBeam )
489{
490 // open file
491 // DEBUG
492 time_t t0 ;
493 time( &t0 ) ;
494 tm *ttm = localtime( &t0 ) ;
495 LogIO os( LogOrigin( "STFiller", "openNRO()", WHERE ) ) ;
496// cout << "STFiller::openNRO() Start time = " << t0
497// << " ("
498// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
499// << " "
500// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
501// << ")" << endl ;
502 os << "Start time = " << t0
503 << " ("
504 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
505 << " "
506 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
507 << ")" << LogIO::POST ;
508
509 // fill STHeader
510 header_ = new STHeader() ;
511
512 if ( nreader_->getHeaderInfo( header_->nchan,
513 header_->npol,
514 nIF_,
515 nBeam_,
516 header_->observer,
517 header_->project,
518 header_->obstype,
519 header_->antennaname,
520 header_->antennaposition,
521 header_->equinox,
522 header_->freqref,
523 header_->reffreq,
524 header_->bandwidth,
525 header_->utc,
526 header_->fluxunit,
527 header_->epoch,
528 header_->poltype ) ) {
529// cout << "STFiller::openNRO() Failed to get header information." << endl ;
530// return ;
531 throw( AipsError("Failed to get header information.") ) ;
532 }
533
534 // set FRAME and BASEFRAME keyword of FREQUENCIES table
535 if ( header_->freqref != "TOPO" ) {
536 table_->frequencies().setFrame( header_->freqref, false ) ;
537 table_->frequencies().setFrame( header_->freqref, true ) ;
538 }
539
540 ifOffset_ = 0;
541 vector<Bool> ifs = nreader_->getIFs() ;
542 if ( whichIF >= 0 ) {
543 if ( whichIF >= 0 && whichIF < nIF_ ) {
544 for ( int i = 0 ; i < nIF_ ; i++ )
545 ifs[i] = False ;
546 ifs[whichIF] = True ;
547 header_->nif = 1;
548 nIF_ = 1;
549 ifOffset_ = whichIF;
550 } else {
551 delete reader_;
552 reader_ = 0;
553 delete header_;
554 header_ = 0;
555 throw(AipsError("Illegal IF selection"));
556 }
557 }
558
559 beamOffset_ = 0;
560 vector<Bool> beams = nreader_->getBeams() ;
561 if (whichBeam>=0) {
562 if (whichBeam>=0 && whichBeam<nBeam_) {
563 for ( int i = 0 ; i < nBeam_ ; i++ )
564 beams[i] = False ;
565 beams[whichBeam] = True;
566 header_->nbeam = 1;
567 nBeam_ = 1;
568 beamOffset_ = whichBeam;
569 } else {
570 delete reader_;
571 reader_ = 0;
572 delete header_;
573 header_ = 0;
574 throw(AipsError("Illegal Beam selection"));
575 }
576 }
577
578 header_->nbeam = nBeam_ ;
579 header_->nif = nIF_ ;
580
581 // set header
582 table_->setHeader( *header_ ) ;
583
584 // DEBUG
585 time_t t1 ;
586 time( &t1 ) ;
587 ttm = localtime( &t1 ) ;
588// cout << "STFiller::openNRO() End time = " << t1
589// << " ("
590// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
591// << " "
592// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
593// << ")" << endl ;
594// cout << "STFiller::openNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
595 os << "End time = " << t1
596 << " ("
597 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
598 << " "
599 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
600 << ")" << endl ;
601 os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
602 os.post() ;
603 //
604
605 return ;
606}
607
608int STFiller::readNRO()
609{
610 // DEBUG
611 time_t t0 ;
612 time( &t0 ) ;
613 tm *ttm = localtime( &t0 ) ;
614 LogIO os( LogOrigin( "STFiller", "readNRO()", WHERE ) ) ;
615// cout << "STFiller::readNRO() Start time = " << t0
616// << " ("
617// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
618// << " "
619// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
620// << ")" << endl ;
621 os << "Start time = " << t0
622 << " ("
623 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
624 << " "
625 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
626 << ")" << LogIO::POST ;
627 //
628
629 // fill row
630 uInt id ;
631 uInt imax = nreader_->getRowNum() ;
632 vector< vector<double > > freqs ;
633 uInt i = 0 ;
634 int count = 0 ;
635 uInt scanno ;
636 uInt cycleno ;
637 uInt beamno ;
638 uInt polno ;
639 vector<double> fqs ;
640 Vector<Double> restfreq ;
641 uInt refbeamno ;
642 Double scantime ;
643 Double interval ;
644 String srcname ;
645 String fieldname ;
646 Array<Float> spectra ;
647 Array<uChar> flagtra ;
648 Array<Float> tsys ;
649 Array<Double> direction ;
650 Float azimuth ;
651 Float elevation ;
652 Float parangle ;
653 Float opacity ;
654 uInt tcalid ;
655 Int fitid ;
656 uInt focusid ;
657 Float temperature ;
658 Float pressure ;
659 Float humidity ;
660 Float windvel ;
661 Float winddir ;
662 Double srcvel ;
663 Array<Double> propermotion ;
664 Vector<Double> srcdir ;
665 Array<Double> scanrate ;
666 for ( i = 0 ; i < imax ; i++ ) {
667 string scanType = nreader_->getScanType( i ) ;
668 Int srcType = -1 ;
669 if ( scanType.compare( 0, 2, "ON") == 0 ) {
670 // os << "ON srcType: " << i << LogIO::POST ;
671 srcType = 0 ;
672 }
673 else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
674 //os << "OFF srcType: " << i << LogIO::POST ;
675 srcType = 1 ;
676 }
677 else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
678 //os << "ZERO srcType: " << i << LogIO::POST ;
679 srcType = 2 ;
680 }
681 else {
682 //os << "Undefined srcType: " << i << LogIO::POST ;
683 srcType = 3 ;
684 }
685
686 // if srcType is 2 (ZERO scan), ignore scan
687 if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
688 TableRow row( table_->table() ) ;
689 TableRecord& rec = row.record();
690
691 if ( nreader_->getScanInfo( i,
692 scanno,
693 cycleno,
694 beamno,
695 polno,
696 fqs,
697 restfreq,
698 refbeamno,
699 scantime,
700 interval,
701 srcname,
702 fieldname,
703 spectra,
704 flagtra,
705 tsys,
706 direction,
707 azimuth,
708 elevation,
709 parangle,
710 opacity,
711 tcalid,
712 fitid,
713 focusid,
714 temperature,
715 pressure,
716 humidity,
717 windvel,
718 winddir,
719 srcvel,
720 propermotion,
721 srcdir,
722 scanrate ) ) {
723// cerr << "STFiller::readNRO() Failed to get scan information." << endl ;
724// return 1 ;
725 throw( AipsError("Failed to get scan information.") ) ;
726 }
727
728 RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
729 *scannoCol = scanno ;
730 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ;
731 *cyclenoCol = cycleno ;
732 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
733 *beamCol = beamno ;
734 RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
735 RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
736 *polnoCol = polno ;
737 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
738 if ( freqs.size() == 0 ) {
739 id = table_->frequencies().addEntry( Double( fqs[0] ),
740 Double( fqs[1] ),
741 Double( fqs[2] ) ) ;
742 *mfreqidCol = id ;
743 *ifCol = id ;
744 freqs.push_back( fqs ) ;
745 }
746 else {
747 int iadd = -1 ;
748 for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
749 //os << "freqs[" << iif << "][1] = " << freqs[iif][1] << LogIO::POST ;
750 double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
751 //os << "fdiff = " << fdiff << LogIO::POST ;
752 if ( fdiff < 1.0e-8 ) {
753 iadd = iif ;
754 break ;
755 }
756 }
757 if ( iadd == -1 ) {
758 id = table_->frequencies().addEntry( Double( fqs[0] ),
759 Double( fqs[1] ),
760 Double( fqs[2] ) ) ;
761 *mfreqidCol = id ;
762 *ifCol = id ;
763 freqs.push_back( fqs ) ;
764 }
765 else {
766 *mfreqidCol = iadd ;
767 *ifCol = iadd ;
768 }
769 }
770 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
771 id = table_->molecules().addEntry( restfreq ) ;
772 *molidCol = id ;
773 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
774 *rbCol = refbeamno ;
775 RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
776 *mjdCol = scantime ;
777 RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
778 *intervalCol = interval ;
779 RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
780 *srcnCol = srcname ;
781 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
782 *srctCol = srcType ;
783 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
784 *fieldnCol = fieldname ;
785 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
786 *specCol = spectra ;
787 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
788 *flagCol = flagtra ;
789 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
790 *tsysCol = tsys ;
791 RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
792 *dirCol = direction ;
793 RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
794 *azCol = azimuth ;
795 RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
796 *elCol = elevation ;
797 //RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
798 //*parCol = parangle ;
799 RecordFieldPtr<Float> tauCol(rec, "OPACITY") ;
800 *tauCol = opacity ;
801 id = table_->tcal().addEntry("",Vector<Float>(1,0.0));
802 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
803 *mcalidCol = id;
804 //RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
805 //*mcalidCol = tcalid ;
806 RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
807 *fitCol = fitid ;
808 //RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
809 //*mfocusidCol = focusid ;
810 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
811 id = table_->focus().addEntry(parangle, 0.0, 0.0, 0.0 ) ;
812 *mfocusidCol = id;
813 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
814 id = table_->weather().addEntry( temperature,
815 pressure,
816 humidity,
817 windvel,
818 winddir ) ;
819 *mweatheridCol = id ;
820 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
821 *svelCol = srcvel ;
822 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
823 *spmCol = propermotion ;
824 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
825 *sdirCol = srcdir ;
826 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
827 *srateCol = scanrate ;
828
829 table_->table().addRow() ;
830 row.put(table_->table().nrow()-1, rec) ;
831 }
832 else {
833 count++ ;
834 }
835 // DEBUG
836 //int rownum = nreader_->getRowNum() ;
837 //os << "Finished row " << i << "/" << rownum << LogIO::POST ;
838 //
839 }
840
841 // DEBUG
842 time_t t1 ;
843 time( &t1 ) ;
844 ttm = localtime( &t1 ) ;
845// cout << "STFiller::readNRO() Processed " << i << " rows" << endl ;
846// cout << "STFiller::readNRO() Added " << i - count << " rows (ignored "
847// << count << " \"ZERO\" scans)" << endl ;
848// cout << "STFiller::readNRO() End time = " << t1
849// << " ("
850// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
851// << " "
852// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
853// << ")" << endl ;
854// cout << "STFiller::readNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
855 os << "Processed " << i << " rows" << endl ;
856 os << "Added " << i - count << " rows (ignored "
857 << count << " \"ZERO\" scans)" << endl ;
858 os.post() ;
859 os << "End time = " << t1
860 << " ("
861 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
862 << " "
863 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
864 << ")" << endl ;
865 os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
866 os.post() ;
867 //
868
869 return 0 ;
870}
871
872Bool STFiller::fileCheck()
873{
874 bool bval = false ;
875
876 // if filename_ is directory, return false
877 File inFile( filename_ ) ;
878 if ( inFile.isDirectory() )
879 return bval ;
880
881 // if beginning of header data is "RW", return true
882 // otherwise, return false ;
883 FILE *fp = fopen( filename_.c_str(), "r" ) ;
884 char buf[9] ;
885 char buf2[80] ;
886 size_t tmpret;
887 tmpret = fread( buf, 4, 1, fp ) ;
888 buf[4] = '\0' ;
889 fseek( fp, 640, SEEK_SET ) ;
890 tmpret = fread( buf2, 80, 1, fp ) ;
891 (void) tmpret; //suppress unused warning
892 if ( ( strncmp( buf, "RW", 2 ) == 0 ) || ( strstr( buf2, "NRO45M" ) != NULL ) ) {
893 bval = true ;
894 }
895 fclose( fp ) ;
896 return bval ;
897}
898
899}//namespace asap
Note: See TracBrowser for help on using the repository browser.