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

Last change on this file since 1488 was 1482, checked in by Takeshi Nakazato, 16 years ago

New Development: No

JIRA Issue: Yes CAS-1043

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: asap_init()

s=sd.scantable('NRO_DATA_FILENAME',False)

Put in Release Notes: No

Module(s): No

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