source: trunk/src/STFiller.cpp@ 2328

Last change on this file since 2328 was 2289, checked in by ShinnosukeKawakami, 13 years ago

merged parallel branch to trunk
  • 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 iftmp ;
638 uInt beamno ;
639 uInt polno ;
640 vector<double> fqs ;
641 Vector<Double> restfreq ;
642 uInt refbeamno ;
643 Double scantime ;
644 Double interval ;
645 String srcname ;
646 String fieldname ;
647 Vector<Float> spectra ;
648 Vector<uChar> flagtra ;
649 Vector<Float> tsys ;
650 Vector<Double> direction ;
651 Float azimuth ;
652 Float elevation ;
653 Float parangle ;
654 Float opacity ;
655 uInt tcalid ;
656 Int fitid ;
657 uInt focusid ;
658 Float temperature ;
659 Float pressure ;
660 Float humidity ;
661 Float windvel ;
662 Float winddir ;
663 Double srcvel ;
664 Vector<Double> propermotion( 2, 0.0 ) ;
665 Vector<Double> srcdir ;
666 Vector<Double> scanrate( 2, 0.0 ) ;
667 for ( i = 0 ; i < imax ; i++ ) {
668 string scanType = nreader_->getScanType( i ) ;
669 Int srcType = -1 ;
670 if ( scanType.compare( 0, 2, "ON") == 0 ) {
671 // os << "ON srcType: " << i << LogIO::POST ;
672 srcType = 0 ;
673 }
674 else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
675 //os << "OFF srcType: " << i << LogIO::POST ;
676 srcType = 1 ;
677 }
678 else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
679 //os << "ZERO srcType: " << i << LogIO::POST ;
680 srcType = 2 ;
681 }
682 else {
683 //os << "Undefined srcType: " << i << LogIO::POST ;
684 srcType = 3 ;
685 }
686
687 // if srcType is 2 (ZERO scan), ignore scan
688 if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
689 TableRow row( table_->table() ) ;
690 TableRecord& rec = row.record();
691
692 if ( nreader_->getScanInfo( i,
693 scanno,
694 cycleno,
695 iftmp,
696 beamno,
697 polno,
698 fqs,
699 restfreq,
700 refbeamno,
701 scantime,
702 interval,
703 srcname,
704 fieldname,
705 spectra,
706 flagtra,
707 tsys,
708 direction,
709 azimuth,
710 elevation,
711 parangle,
712 opacity,
713 tcalid,
714 fitid,
715 focusid,
716 temperature,
717 pressure,
718 humidity,
719 windvel,
720 winddir,
721 srcvel,
722 propermotion,
723 srcdir,
724 scanrate ) ) {
725// cerr << "STFiller::readNRO() Failed to get scan information." << endl ;
726// return 1 ;
727 throw( AipsError("Failed to get scan information.") ) ;
728 }
729
730 RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
731 *scannoCol = scanno ;
732 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ;
733 *cyclenoCol = cycleno ;
734 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
735 *beamCol = beamno ;
736 RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
737 RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
738 *polnoCol = polno ;
739 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
740 if ( freqs.size() == 0 ) {
741 id = table_->frequencies().addEntry( Double( fqs[0] ),
742 Double( fqs[1] ),
743 Double( fqs[2] ) ) ;
744 *mfreqidCol = id ;
745 *ifCol = id ;
746 freqs.push_back( fqs ) ;
747 }
748 else {
749 int iadd = -1 ;
750 for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
751 //os << "freqs[" << iif << "][1] = " << freqs[iif][1] << LogIO::POST ;
752 double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
753 //os << "fdiff = " << fdiff << LogIO::POST ;
754 if ( fdiff < 1.0e-8 ) {
755 iadd = iif ;
756 break ;
757 }
758 }
759 if ( iadd == -1 ) {
760 id = table_->frequencies().addEntry( Double( fqs[0] ),
761 Double( fqs[1] ),
762 Double( fqs[2] ) ) ;
763 *mfreqidCol = id ;
764 *ifCol = id ;
765 freqs.push_back( fqs ) ;
766 }
767 else {
768 *mfreqidCol = iadd ;
769 *ifCol = iadd ;
770 }
771 }
772 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
773 id = table_->molecules().addEntry( restfreq ) ;
774 *molidCol = id ;
775 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
776 *rbCol = refbeamno ;
777 RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
778 *mjdCol = scantime ;
779 RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
780 *intervalCol = interval ;
781 RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
782 *srcnCol = srcname ;
783 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
784 *srctCol = srcType ;
785 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
786 *fieldnCol = fieldname ;
787 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
788 *specCol = spectra ;
789 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
790 *flagCol = flagtra ;
791 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
792 *tsysCol = tsys ;
793 RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
794 *dirCol = direction ;
795 RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
796 *azCol = azimuth ;
797 RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
798 *elCol = elevation ;
799 //RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
800 //*parCol = parangle ;
801 RecordFieldPtr<Float> tauCol(rec, "OPACITY") ;
802 *tauCol = opacity ;
803 id = table_->tcal().addEntry("",Vector<Float>(1,0.0));
804 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
805 *mcalidCol = id;
806 //RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
807 //*mcalidCol = tcalid ;
808 RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
809 *fitCol = fitid ;
810 //RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
811 //*mfocusidCol = focusid ;
812 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
813 id = table_->focus().addEntry(parangle, 0.0, 0.0, 0.0 ) ;
814 *mfocusidCol = id;
815 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
816 id = table_->weather().addEntry( temperature,
817 pressure,
818 humidity,
819 windvel,
820 winddir ) ;
821 *mweatheridCol = id ;
822 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
823 *svelCol = srcvel ;
824 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
825 *spmCol = propermotion ;
826 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
827 *sdirCol = srcdir ;
828 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
829 *srateCol = scanrate ;
830
831 table_->table().addRow() ;
832 row.put(table_->table().nrow()-1, rec) ;
833 }
834 else {
835 count++ ;
836 }
837 // DEBUG
838 //int rownum = nreader_->getRowNum() ;
839 //os << "Finished row " << i << "/" << rownum << LogIO::POST ;
840 //
841 }
842
843 // DEBUG
844 time_t t1 ;
845 time( &t1 ) ;
846 ttm = localtime( &t1 ) ;
847// cout << "STFiller::readNRO() Processed " << i << " rows" << endl ;
848// cout << "STFiller::readNRO() Added " << i - count << " rows (ignored "
849// << count << " \"ZERO\" scans)" << endl ;
850// cout << "STFiller::readNRO() End time = " << t1
851// << " ("
852// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
853// << " "
854// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
855// << ")" << endl ;
856// cout << "STFiller::readNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
857 os << "Processed " << i << " rows" << endl ;
858 os << "Added " << i - count << " rows (ignored "
859 << count << " \"ZERO\" scans)" << endl ;
860 os.post() ;
861 os << "End time = " << t1
862 << " ("
863 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
864 << " "
865 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
866 << ")" << endl ;
867 os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
868 os.post() ;
869 //
870
871 return 0 ;
872}
873
874Bool STFiller::fileCheck()
875{
876 bool bval = false ;
877
878 // if filename_ is directory, return false
879 File inFile( filename_ ) ;
880 if ( inFile.isDirectory() )
881 return bval ;
882
883 // if beginning of header data is "RW", return true
884 // otherwise, return false ;
885 FILE *fp = fopen( filename_.c_str(), "r" ) ;
886 char buf[9] ;
887 char buf2[80] ;
888 size_t tmpret;
889 tmpret = fread( buf, 4, 1, fp ) ;
890 buf[4] = '\0' ;
891 fseek( fp, 640, SEEK_SET ) ;
892 tmpret = fread( buf2, 80, 1, fp ) ;
893 (void) tmpret; //suppress unused warning
894 if ( ( strncmp( buf, "RW", 2 ) == 0 ) || ( strstr( buf2, "NRO45M" ) != NULL ) ) {
895 bval = true ;
896 }
897 fclose( fp ) ;
898 return bval ;
899}
900
901}//namespace asap
Note: See TracBrowser for help on using the repository browser.