source: branches/casa-release-4_3/src/STFiller.cpp@ 3090

Last change on this file since 3090 was 2658, checked in by Malte Marquarding, 12 years ago

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