source: trunk/src/STFiller.cpp@ 2075

Last change on this file since 2075 was 1841, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: unit test for sdsave task (will be ready soon)

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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