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

Last change on this file since 1628 was 1616, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-729, CAS-1147

Ready to Release: Yes

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...

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