source: branches/newfiller/src/STFiller.cpp@ 1784

Last change on this file since 1784 was 1757, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.0 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;
392 /// @todo this has to change when nchan isn't global anymore
393 //id = table_->frequencies().addEntry(Double(header_->nchan/2),
394 // pksrec.refFreq, pksrec.freqInc);
395 if ( pksrec.nchan == 1 ) {
396 id = table_->frequencies().addEntry(Double(0),
397 pksrec.refFreq, pksrec.freqInc);
398 }
399 else {
400 id = table_->frequencies().addEntry(Double(pksrec.nchan/2),
401 pksrec.refFreq, pksrec.freqInc);
402 }
403 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
404 *mfreqidCol = id;
405 //*ifCol = id;
406
407 id = table_->molecules().addEntry(pksrec.restFreq);
408 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
409 *molidCol = id;
410
411 id = table_->tcal().addEntry(pksrec.tcalTime, pksrec.tcal);
412 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
413 *mcalidCol = id;
414 id = table_->weather().addEntry(pksrec.temperature, pksrec.pressure,
415 pksrec.humidity, pksrec.windSpeed,
416 pksrec.windAz);
417 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
418 *mweatheridCol = id;
419
420 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
421 id = table_->focus().addEntry(pksrec.parAngle, pksrec.focusAxi,
422 pksrec.focusTan, pksrec.focusRot);
423 *mfocusidCol = id;
424 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
425 *dirCol = pksrec.direction;
426 RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
427 *azCol = pksrec.azimuth;
428 RecordFieldPtr<Float> elCol(rec, "ELEVATION");
429 *elCol = pksrec.elevation;
430
431 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
432 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
433 RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
434
435 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
436 // Turn the (nchan,npol) matrix and possible complex xPol vector
437 // into 2-4 rows in the scantable
438 Vector<Float> tsysvec(1);
439 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
440 uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
441 for ( uInt i=0; i< npol; ++i ) {
442 tsysvec = pksrec.tsys(i);
443 *tsysCol = tsysvec;
444 if (isGBTFITS)
445 *polnoCol = pksrec.polNo ;
446 else
447 *polnoCol = i;
448
449 *specCol = pksrec.spectra.column(i);
450 *flagCol = pksrec.flagged.column(i);
451 table_->table().addRow();
452 row.put(table_->table().nrow()-1, rec);
453 }
454
455 RecordFieldPtr< uInt > flagrowCol(rec, "FLAGROW");
456 *flagrowCol = pksrec.flagrow;
457
458 if ( haveXPol_[0] ) {
459 // no tsys given for xpol, so emulate it
460 tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
461 *tsysCol = tsysvec;
462 // add real part of cross pol
463 *polnoCol = 2;
464 Vector<Float> r(real(pksrec.xPol));
465 *specCol = r;
466 // make up flags from linears
467 /// @fixme this has to be a bitwise or of both pols
468 *flagCol = pksrec.flagged.column(0);// | pksrec.flagged.column(1);
469 table_->table().addRow();
470 row.put(table_->table().nrow()-1, rec);
471 // ad imaginary part of cross pol
472 *polnoCol = 3;
473 Vector<Float> im(imag(pksrec.xPol));
474 *specCol = im;
475 table_->table().addRow();
476 row.put(table_->table().nrow()-1, rec);
477 }
478#ifdef HAS_ALMA
479 fillpm._update(n);
480#endif
481 }
482 if (status > 0) {
483 close();
484 throw(AipsError("Reading error occured, data possibly corrupted."));
485 }
486#ifdef HAS_ALMA
487 fillpm.done();
488#endif
489 return status;
490}
491
492/**
493 * For NRO data
494 *
495 * 2008/11/11 Takeshi Nakazato
496 **/
497void STFiller::openNRO( int whichIF, int whichBeam )
498{
499 // open file
500 // DEBUG
501 time_t t0 ;
502 time( &t0 ) ;
503 tm *ttm = localtime( &t0 ) ;
504 LogIO os( LogOrigin( "STFiller", "openNRO()", WHERE ) ) ;
505// cout << "STFiller::openNRO() Start time = " << t0
506// << " ("
507// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
508// << " "
509// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
510// << ")" << endl ;
511 os << "Start time = " << t0
512 << " ("
513 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
514 << " "
515 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
516 << ")" << LogIO::POST ;
517
518 // fill STHeader
519 header_ = new STHeader() ;
520
521 if ( nreader_->getHeaderInfo( header_->nchan,
522 header_->npol,
523 nIF_,
524 nBeam_,
525 header_->observer,
526 header_->project,
527 header_->obstype,
528 header_->antennaname,
529 header_->antennaposition,
530 header_->equinox,
531 header_->freqref,
532 header_->reffreq,
533 header_->bandwidth,
534 header_->utc,
535 header_->fluxunit,
536 header_->epoch,
537 header_->poltype ) ) {
538// cout << "STFiller::openNRO() Failed to get header information." << endl ;
539// return ;
540 throw( AipsError("Failed to get header information.") ) ;
541 }
542
543 // set FRAME and BASEFRAME keyword of FREQUENCIES table
544 if ( header_->freqref != "TOPO" ) {
545 table_->frequencies().setFrame( header_->freqref, false ) ;
546 table_->frequencies().setFrame( header_->freqref, true ) ;
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.