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

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

New Development: No

JIRA Issue: Yes CAS-1810

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: Added 'antenna' parameter to scantable constructor

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): atnf

Description: Describe your changes here...

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