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

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

New Development: No

JIRA Issue: Yes CAS-1683

Ready to Release: No

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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