source: trunk/src/STFiller.cpp@ 2657

Last change on this file since 2657 was 2657, checked in by Malte Marquarding, 13 years ago

Ticket #280: Use FITSSpectralUtil::specsysFromFrame to write correct FITS output frame
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.7 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#include <fits/FITS/FITSSpectralUtil.h>
39
40#include <time.h>
41
42
43#include "STDefs.h"
44#include "STAttr.h"
45
46#include "STFiller.h"
47#include "STHeader.h"
48
49using namespace casa;
50
51namespace asap {
52
53STFiller::STFiller() :
54 reader_(0),
55 header_(0),
56 table_(0),
57 refRx_(".*(e|w|_R)$"),
58 nreader_(0)
59{
60}
61
62STFiller::STFiller( CountedPtr< Scantable > stbl ) :
63 reader_(0),
64 header_(0),
65 table_(stbl),
66 refRx_(".*(e|w|_R)$"),
67 nreader_(0)
68{
69}
70
71STFiller::STFiller(const std::string& filename, int whichIF, int whichBeam ) :
72 reader_(0),
73 header_(0),
74 table_(0),
75 refRx_(".*(e|w|_R)$"),
76 nreader_(0)
77{
78 open(filename, "", whichIF, whichBeam);
79}
80
81STFiller::~STFiller()
82{
83 close();
84}
85
86void STFiller::open( const std::string& filename, const std::string& antenna, int whichIF, int whichBeam, casa::Bool getPt )
87{
88 if (table_.null()) {
89 table_ = new Scantable();
90 }
91 if (reader_) { delete reader_; reader_ = 0; }
92 Bool haveBase, haveSpectra;
93
94 String inName(filename);
95 Path path(inName);
96 inName = path.expandedName();
97
98 File file(inName);
99 if ( !file.exists() ) {
100 throw(AipsError("File does not exist"));
101 }
102 filename_ = inName;
103
104 // Create reader and fill in values for arguments
105 String format;
106 Vector<Bool> beams, ifs;
107 Vector<uInt> nchans,npols;
108
109 //
110 // if isNRO_ is true, try NROReader
111 //
112 // 2008/11/11 Takeshi Nakazato
113 isNRO_ = fileCheck() ;
114 if ( isNRO_ ) {
115 if ( (nreader_ = getNROReader( inName, format )) == 0 ) {
116 throw(AipsError("Creation of NROReader failed")) ;
117 }
118 else {
119 openNRO( whichIF, whichBeam ) ;
120 return ;
121 }
122 }
123 //
124
125 if ( (reader_ = getPKSreader(inName, antenna, 0, 0, format, beams, ifs,
126 nchans, npols, haveXPol_,haveBase, haveSpectra
127 )) == 0 ) {
128 throw(AipsError("Creation of PKSreader failed"));
129 }
130 if (!haveSpectra) {
131 delete reader_;
132 reader_ = 0;
133 throw(AipsError("No spectral data in file."));
134 return;
135 }
136 nBeam_ = beams.nelements();
137 nIF_ = ifs.nelements();
138 // Get basic parameters.
139 if ( anyEQ(haveXPol_, True) ) {
140 pushLog("Cross polarization present");
141 for (uInt i=0; i< npols.nelements();++i) {
142 if (npols[i] < 3) npols[i] += 2;// Convert Complex -> 2 Floats
143 }
144 }
145 if (header_) delete header_;
146 header_ = new STHeader();
147 header_->nchan = max(nchans);
148 header_->npol = max(npols);
149 header_->nbeam = nBeam_;
150
151 Int status = reader_->getHeader(header_->observer, header_->project,
152 header_->antennaname, header_->antennaposition,
153 header_->obstype,
154 header_->fluxunit,
155 header_->equinox,
156 header_->freqref,
157 header_->utc, header_->reffreq,
158 header_->bandwidth);
159
160 if (status) {
161 delete reader_;
162 reader_ = 0;
163 delete header_;
164 header_ = 0;
165 throw(AipsError("Failed to get header."));
166 }
167 if ((header_->obstype).matches("*SW*")) {
168 // need robust way here - probably read ahead of next timestamp
169 pushLog("Header indicates frequency switched observation.\n"
170 "setting # of IFs = 1 ");
171 nIF_ = 1;
172 header_->obstype = String("fswitch");
173 }
174 // Determine Telescope and set brightness unit
175
176
177 Bool throwIt = False;
178 Instrument inst = STAttr::convertInstrument(header_->antennaname, throwIt);
179
180 if (inst==ATMOPRA || inst==TIDBINBILLA) {
181 header_->fluxunit = "K";
182 } else {
183 // downcase for use with Quanta
184 if (header_->fluxunit == "JY") {
185 header_->fluxunit = "Jy";
186 }
187 }
188 STAttr stattr;
189 header_->poltype = stattr.feedPolType(inst);
190 header_->nif = nIF_;
191 header_->epoch = "UTC";
192 // *** header_->frequnit = "Hz"
193 // Apply selection criteria.
194 Vector<Int> ref;
195 ifOffset_ = 0;
196 if (whichIF>=0) {
197 if (whichIF>=0 && whichIF<nIF_) {
198 ifs = False;
199 ifs(whichIF) = True;
200 header_->nif = 1;
201 nIF_ = 1;
202 ifOffset_ = whichIF;
203 } else {
204 delete reader_;
205 reader_ = 0;
206 delete header_;
207 header_ = 0;
208 throw(AipsError("Illegal IF selection"));
209 }
210 }
211 beamOffset_ = 0;
212 if (whichBeam>=0) {
213 if (whichBeam>=0 && whichBeam<nBeam_) {
214 beams = False;
215 beams(whichBeam) = True;
216 header_->nbeam = 1;
217 nBeam_ = 1;
218 beamOffset_ = whichBeam;
219 } else {
220 delete reader_;
221 reader_ = 0;
222 delete header_;
223 header_ = 0;
224 throw(AipsError("Illegal Beam selection"));
225 }
226 }
227 Vector<Int> start(nIF_, 1);
228 Vector<Int> end(nIF_, 0);
229 reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False,
230 getPt);
231 table_->setHeader(*header_);
232 //For MS, add the location of POINTING of the input MS so one get
233 //pointing data from there, if necessary.
234 //Also find nrow in MS
235 nInDataRow = 0;
236 if (format == "MS2") {
237 Path datapath(inName);
238 String ptTabPath = datapath.absoluteName();
239 Table inMS(ptTabPath);
240 nInDataRow = inMS.nrow();
241 ptTabPath.append("/POINTING");
242 table_->table().rwKeywordSet().define("POINTING", ptTabPath);
243 if ((header_->antennaname).matches("GBT")) {
244 String GOTabPath = datapath.absoluteName();
245 GOTabPath.append("/GBT_GO");
246 table_->table().rwKeywordSet().define("GBT_GO", GOTabPath);
247 }
248 }
249 String freqFrame = header_->freqref;
250 //translate frequency reference frame back to
251 //MS style (as PKSMS2reader converts the original frame
252 //in FITS standard style)
253 MFrequency::Types fframe;
254 FITSSpectralUtil::frameFromSpecsys(fframe, freqFrame);
255 // set both "FRAME" and "BASEFRAME"
256 table_->frequencies().setFrame(fframe, false);
257 table_->frequencies().setFrame(fframe,true);
258
259 //table_->focus().setParallactify(true);
260}
261
262void STFiller::close( )
263{
264 delete reader_;reader_=0;
265 delete nreader_;nreader_=0;
266 delete header_;header_=0;
267 table_ = 0;
268}
269
270int asap::STFiller::read( )
271{
272 int status = 0;
273
274 //
275 // for NRO data
276 //
277 // 2008/11/12 Takeshi Nakazato
278 if ( isNRO_ ) {
279 status = readNRO() ;
280 return status ;
281 }
282 //
283
284/**
285 Int beamNo, IFno, refBeam, scanNo, cycleNo;
286 Float azimuth, elevation, focusAxi, focusRot, focusTan,
287 humidity, parAngle, pressure, temperature, windAz, windSpeed;
288 Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
289 String fieldName, srcName, tcalTime, obsType;
290 Vector<Float> calFctr, sigma, tcal, tsys;
291 Matrix<Float> baseLin, baseSub;
292 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq(1);
293 Matrix<Float> spectra;
294 Matrix<uChar> flagtra;
295 Complex xCalFctr;
296 Vector<Complex> xPol;
297**/
298
299#ifdef HAS_ALMA
300 Double min = 0.0;
301 Double max = nInDataRow;
302 ProgressMeter fillpm(min, max, "Data importing progress");
303#endif
304 PKSrecord pksrec;
305 pksrec.srcType=-1;
306 int n = 0;
307 bool isGBTFITS = false ;
308 if ((header_->antennaname.find( "GBT" ) != String::npos) && File(filename_).isRegular()) {
309 FILE *fp = fopen( filename_.c_str(), "r" ) ;
310 fseek( fp, 640, SEEK_SET ) ;
311 char buf[81] ;
312 size_t tmpret = fread( buf, 80, 1, fp ) ;
313 (void) tmpret;
314 buf[80] = '\0' ;
315 if ( strstr( buf, "NRAO_GBT" ) != NULL ) {
316 isGBTFITS = true ;
317 }
318 fclose( fp ) ;
319 }
320 while ( status == 0 ) {
321 status = reader_->read(pksrec);
322 if ( status != 0 ) break;
323 n += 1;
324 Regex filterrx(".*[SL|PA]$");
325 Regex obsrx("^AT.+");
326 if ( header_->antennaname.matches(obsrx) &&
327 pksrec.obsType.matches(filterrx)) {
328 //cerr << "ignoring paddle scan" << endl;
329 continue;
330 }
331 TableRow row(table_->table());
332 TableRecord& rec = row.record();
333 // fields that don't get used and are just passed through asap
334 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
335 // MRC changed type from double to float
336 Vector<Double> sratedbl(pksrec.scanRate.nelements());
337 convertArray(sratedbl, pksrec.scanRate);
338 *srateCol = sratedbl;
339 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
340 *spmCol = pksrec.srcPM;
341 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
342 *sdirCol = pksrec.srcDir;
343 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
344 *svelCol = pksrec.srcVel;
345 // the real stuff
346 RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
347 *fitCol = -1;
348 RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
349 *scanoCol = pksrec.scanNo-1;
350 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
351 *cyclenoCol = pksrec.cycleNo-1;
352 RecordFieldPtr<Double> mjdCol(rec, "TIME");
353 *mjdCol = pksrec.mjd;
354 RecordFieldPtr<Double> intCol(rec, "INTERVAL");
355 *intCol = pksrec.interval;
356 RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
357 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
358 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
359 *fieldnCol = pksrec.fieldName;
360 // try to auto-identify if it is on or off.
361 Regex rx(refRx_);
362 Regex rx2("_S$");
363 Int match = pksrec.srcName.matches(rx);
364 if (match) {
365 *srcnCol = pksrec.srcName;
366 } else {
367 *srcnCol = pksrec.srcName.before(rx2);
368 }
369 //*srcnCol = pksrec.srcName;//.before(rx2);
370 *srctCol = match;
371 if ( pksrec.srcType != -1 ) {
372 *srctCol = pksrec.srcType ;
373 }
374 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
375 *beamCol = pksrec.beamNo-beamOffset_-1;
376 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
377 Int rb = -1;
378 if (nBeam_ > 1 ) rb = pksrec.refBeam-1;
379 *rbCol = rb;
380 RecordFieldPtr<uInt> ifCol(rec, "IFNO");
381 *ifCol = pksrec.IFno-ifOffset_- 1;
382 uInt id = table_->frequencies().addEntry(Double(pksrec.spectra.nrow()/2),
383 pksrec.refFreq, pksrec.freqInc);
384 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
385 *mfreqidCol = id;
386 //*ifCol = id;
387
388 id = table_->molecules().addEntry(pksrec.restFreq);
389 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
390 *molidCol = id;
391
392 id = table_->tcal().addEntry(pksrec.tcalTime, pksrec.tcal);
393 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
394 *mcalidCol = id;
395 id = table_->weather().addEntry(pksrec.temperature, pksrec.pressure,
396 pksrec.humidity, pksrec.windSpeed,
397 pksrec.windAz);
398 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
399 *mweatheridCol = id;
400
401 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
402 id = table_->focus().addEntry(pksrec.parAngle, pksrec.focusAxi,
403 pksrec.focusTan, pksrec.focusRot);
404 *mfocusidCol = id;
405 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
406 *dirCol = pksrec.direction;
407 RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
408 *azCol = pksrec.azimuth;
409 RecordFieldPtr<Float> elCol(rec, "ELEVATION");
410 *elCol = pksrec.elevation;
411
412 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
413 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
414 RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
415
416 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
417 // Turn the (nchan,npol) matrix and possible complex xPol vector
418 // into 2-4 rows in the scantable
419 Vector<Float> tsysvec(1);
420 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
421 uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
422 for ( uInt i=0; i< npol; ++i ) {
423 tsysvec = pksrec.tsys(i);
424 *tsysCol = tsysvec;
425 if (isGBTFITS)
426 *polnoCol = pksrec.polNo ;
427 else
428 *polnoCol = i;
429
430 *specCol = pksrec.spectra.column(i);
431 *flagCol = pksrec.flagged.column(i);
432 table_->table().addRow();
433 row.put(table_->table().nrow()-1, rec);
434 }
435
436 RecordFieldPtr< uInt > flagrowCol(rec, "FLAGROW");
437 *flagrowCol = pksrec.flagrow;
438
439 if ( haveXPol_[0] ) {
440 // no tsys given for xpol, so emulate it
441 tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
442 *tsysCol = tsysvec;
443 // add real part of cross pol
444 *polnoCol = 2;
445 Vector<Float> r(real(pksrec.xPol));
446 *specCol = r;
447 // make up flags from linears
448 /// @fixme this has to be a bitwise or of both pols
449 *flagCol = pksrec.flagged.column(0);// | pksrec.flagged.column(1);
450 table_->table().addRow();
451 row.put(table_->table().nrow()-1, rec);
452 // ad imaginary part of cross pol
453 *polnoCol = 3;
454 Vector<Float> im(imag(pksrec.xPol));
455 *specCol = im;
456 table_->table().addRow();
457 row.put(table_->table().nrow()-1, rec);
458 }
459#ifdef HAS_ALMA
460 fillpm._update(n);
461#endif
462 }
463 if (status > 0) {
464 close();
465 throw(AipsError("Reading error occured, data possibly corrupted."));
466 }
467#ifdef HAS_ALMA
468 fillpm.done();
469#endif
470 return status;
471}
472
473/**
474 * For NRO data
475 *
476 * 2008/11/11 Takeshi Nakazato
477 **/
478void STFiller::openNRO( int whichIF, int whichBeam )
479{
480 // open file
481 // DEBUG
482 time_t t0 ;
483 time( &t0 ) ;
484 tm *ttm = localtime( &t0 ) ;
485 LogIO os( LogOrigin( "STFiller", "openNRO()", WHERE ) ) ;
486// cout << "STFiller::openNRO() Start time = " << t0
487// << " ("
488// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
489// << " "
490// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
491// << ")" << endl ;
492 os << "Start time = " << t0
493 << " ("
494 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
495 << " "
496 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
497 << ")" << LogIO::POST ;
498
499 // fill STHeader
500 header_ = new STHeader() ;
501
502 if ( nreader_->getHeaderInfo( header_->nchan,
503 header_->npol,
504 nIF_,
505 nBeam_,
506 header_->observer,
507 header_->project,
508 header_->obstype,
509 header_->antennaname,
510 header_->antennaposition,
511 header_->equinox,
512 header_->freqref,
513 header_->reffreq,
514 header_->bandwidth,
515 header_->utc,
516 header_->fluxunit,
517 header_->epoch,
518 header_->poltype ) ) {
519// cout << "STFiller::openNRO() Failed to get header information." << endl ;
520// return ;
521 throw( AipsError("Failed to get header information.") ) ;
522 }
523
524 // set FRAME and BASEFRAME keyword of FREQUENCIES table
525 if ( header_->freqref != "TOPO" ) {
526 table_->frequencies().setFrame( header_->freqref, false ) ;
527 table_->frequencies().setFrame( header_->freqref, true ) ;
528 }
529
530 ifOffset_ = 0;
531 vector<Bool> ifs = nreader_->getIFs() ;
532 if ( whichIF >= 0 ) {
533 if ( whichIF >= 0 && whichIF < nIF_ ) {
534 for ( int i = 0 ; i < nIF_ ; i++ )
535 ifs[i] = False ;
536 ifs[whichIF] = True ;
537 header_->nif = 1;
538 nIF_ = 1;
539 ifOffset_ = whichIF;
540 } else {
541 delete reader_;
542 reader_ = 0;
543 delete header_;
544 header_ = 0;
545 throw(AipsError("Illegal IF selection"));
546 }
547 }
548
549 beamOffset_ = 0;
550 vector<Bool> beams = nreader_->getBeams() ;
551 if (whichBeam>=0) {
552 if (whichBeam>=0 && whichBeam<nBeam_) {
553 for ( int i = 0 ; i < nBeam_ ; i++ )
554 beams[i] = False ;
555 beams[whichBeam] = True;
556 header_->nbeam = 1;
557 nBeam_ = 1;
558 beamOffset_ = whichBeam;
559 } else {
560 delete reader_;
561 reader_ = 0;
562 delete header_;
563 header_ = 0;
564 throw(AipsError("Illegal Beam selection"));
565 }
566 }
567
568 header_->nbeam = nBeam_ ;
569 header_->nif = nIF_ ;
570
571 // set header
572 table_->setHeader( *header_ ) ;
573
574 // DEBUG
575 time_t t1 ;
576 time( &t1 ) ;
577 ttm = localtime( &t1 ) ;
578// cout << "STFiller::openNRO() End time = " << t1
579// << " ("
580// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
581// << " "
582// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
583// << ")" << endl ;
584// cout << "STFiller::openNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
585 os << "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 os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
592 os.post() ;
593 //
594
595 return ;
596}
597
598int STFiller::readNRO()
599{
600 // DEBUG
601 time_t t0 ;
602 time( &t0 ) ;
603 tm *ttm = localtime( &t0 ) ;
604 LogIO os( LogOrigin( "STFiller", "readNRO()", WHERE ) ) ;
605// cout << "STFiller::readNRO() Start time = " << t0
606// << " ("
607// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
608// << " "
609// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
610// << ")" << endl ;
611 os << "Start time = " << t0
612 << " ("
613 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
614 << " "
615 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
616 << ")" << LogIO::POST ;
617 //
618
619 // fill row
620 uInt id ;
621 uInt imax = nreader_->getRowNum() ;
622 vector< vector<double > > freqs ;
623 uInt i = 0 ;
624 int count = 0 ;
625 uInt scanno ;
626 uInt cycleno ;
627 uInt iftmp ;
628 uInt beamno ;
629 uInt polno ;
630 vector<double> fqs ;
631 Vector<Double> restfreq ;
632 uInt refbeamno ;
633 Double scantime ;
634 Double interval ;
635 String srcname ;
636 String fieldname ;
637 Vector<Float> spectra ;
638 Vector<uChar> flagtra ;
639 Vector<Float> tsys ;
640 Vector<Double> direction ;
641 Float azimuth ;
642 Float elevation ;
643 Float parangle ;
644 Float opacity ;
645 uInt tcalid ;
646 Int fitid ;
647 uInt focusid ;
648 Float temperature ;
649 Float pressure ;
650 Float humidity ;
651 Float windvel ;
652 Float winddir ;
653 Double srcvel ;
654 Vector<Double> propermotion( 2, 0.0 ) ;
655 Vector<Double> srcdir ;
656 Vector<Double> scanrate( 2, 0.0 ) ;
657 for ( i = 0 ; i < imax ; i++ ) {
658 string scanType = nreader_->getScanType( i ) ;
659 Int srcType = -1 ;
660 if ( scanType.compare( 0, 2, "ON") == 0 ) {
661 // os << "ON srcType: " << i << LogIO::POST ;
662 srcType = 0 ;
663 }
664 else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
665 //os << "OFF srcType: " << i << LogIO::POST ;
666 srcType = 1 ;
667 }
668 else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
669 //os << "ZERO srcType: " << i << LogIO::POST ;
670 srcType = 2 ;
671 }
672 else {
673 //os << "Undefined srcType: " << i << LogIO::POST ;
674 srcType = 3 ;
675 }
676
677 // if srcType is 2 (ZERO scan), ignore scan
678 if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
679 TableRow row( table_->table() ) ;
680 TableRecord& rec = row.record();
681
682 if ( nreader_->getScanInfo( i,
683 scanno,
684 cycleno,
685 iftmp,
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 id = table_->tcal().addEntry("",Vector<Float>(1,0.0));
794 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
795 *mcalidCol = id;
796 //RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
797 //*mcalidCol = tcalid ;
798 RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
799 *fitCol = fitid ;
800 //RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
801 //*mfocusidCol = focusid ;
802 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
803 id = table_->focus().addEntry(parangle, 0.0, 0.0, 0.0 ) ;
804 *mfocusidCol = id;
805 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
806 id = table_->weather().addEntry( temperature,
807 pressure,
808 humidity,
809 windvel,
810 winddir ) ;
811 *mweatheridCol = id ;
812 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
813 *svelCol = srcvel ;
814 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
815 *spmCol = propermotion ;
816 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
817 *sdirCol = srcdir ;
818 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
819 *srateCol = scanrate ;
820
821 table_->table().addRow() ;
822 row.put(table_->table().nrow()-1, rec) ;
823 }
824 else {
825 count++ ;
826 }
827 // DEBUG
828 //int rownum = nreader_->getRowNum() ;
829 //os << "Finished row " << i << "/" << rownum << LogIO::POST ;
830 //
831 }
832
833 // DEBUG
834 time_t t1 ;
835 time( &t1 ) ;
836 ttm = localtime( &t1 ) ;
837// cout << "STFiller::readNRO() Processed " << i << " rows" << endl ;
838// cout << "STFiller::readNRO() Added " << i - count << " rows (ignored "
839// << count << " \"ZERO\" scans)" << endl ;
840// cout << "STFiller::readNRO() End time = " << t1
841// << " ("
842// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
843// << " "
844// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
845// << ")" << endl ;
846// cout << "STFiller::readNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
847 os << "Processed " << i << " rows" << endl ;
848 os << "Added " << i - count << " rows (ignored "
849 << count << " \"ZERO\" scans)" << endl ;
850 os.post() ;
851 os << "End time = " << t1
852 << " ("
853 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
854 << " "
855 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
856 << ")" << endl ;
857 os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
858 os.post() ;
859 //
860
861 return 0 ;
862}
863
864Bool STFiller::fileCheck()
865{
866 bool bval = false ;
867
868 // if filename_ is directory, return false
869 File inFile( filename_ ) ;
870 if ( inFile.isDirectory() )
871 return bval ;
872
873 // if beginning of header data is "RW", return true
874 // otherwise, return false ;
875 FILE *fp = fopen( filename_.c_str(), "r" ) ;
876 char buf[9] ;
877 char buf2[80] ;
878 size_t tmpret;
879 tmpret = fread( buf, 4, 1, fp ) ;
880 buf[4] = '\0' ;
881 fseek( fp, 640, SEEK_SET ) ;
882 tmpret = fread( buf2, 80, 1, fp ) ;
883 (void) tmpret; //suppress unused warning
884 if ( ( strncmp( buf, "RW", 2 ) == 0 ) || ( strstr( buf2, "NRO45M" ) != NULL ) ) {
885 bval = true ;
886 }
887 fclose( fp ) ;
888 return bval ;
889}
890
891}//namespace asap
Note: See TracBrowser for help on using the repository browser.