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

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

New Development: No

JIRA Issue: Yes CAS-1683

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: read/write GBT SDFITS data

Put in Release Notes: No

Module(s): atnf

Description: Describe your changes here...

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