source: branches/parallelCasa3.3/src/STFiller.cpp@ 2493

Last change on this file since 2493 was 2263, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: No

Ready for Test: No

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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