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

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

New Development: No

JIRA Issue: Yes CAS-1799

Ready to Release: No

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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