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

Last change on this file since 1490 was 1490, 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: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Changed criteria if new row is added to the

FREQUENCIES subtable or not.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.7 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 //cout << "STFiller::openNRO() Velocity Definition = " << nheader->getVDEF() << endl ;
555
556 // DEBUG
557 time_t t1 ;
558 time( &t1 ) ;
559 ttm = localtime( &t1 ) ;
560 cout << "STFiller::openNRO() End time = " << t1
561 << " ("
562 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
563 << " "
564 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
565 << ")" << endl ;
566 cout << "STFiller::openNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
567 //
568
569 return ;
570}
571
572int STFiller::readNRO()
573{
574 // DEBUG
575 time_t t0 ;
576 time( &t0 ) ;
577 tm *ttm = localtime( &t0 ) ;
578 cout << "STFiller::readNRO() Start time = " << t0
579 << " ("
580 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
581 << " "
582 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
583 << ")" << endl ;
584 //
585
586 // get header
587 //vector<NRODataset *> data = nreader_->getData() ;
588 NROHeader *h = nreader_->getHeader() ;
589
590 // fill row
591 uInt id ;
592 uInt imax = nreader_->getRowNum() ;
593 vector< vector<double > > freqs ;
594 //uInt imax = 30 ;
595 uInt i = 0 ;
596 int count = 0 ;
597 for ( i = 0 ; i < imax ; i++ ) {
598 if( nreader_->getData( i ) != 0 ) {
599 cerr << "STFiller::readNRO() error while reading row " << i << endl ;
600 return -1 ;
601 }
602 NRODataset *d = nreader_->getData() ;
603
604 //char *scanType = d->getSCANTP() ;
605 char *scanType = d->SCANTP ;
606 Int srcType = -1 ;
607 if ( strncmp( scanType, "ON", 2 ) == 0 ) {
608 srcType = 0 ;
609 }
610 else if ( strncmp( scanType, "OFF", 3 ) == 0 ) {
611 //cout << "OFF srcType: " << i << endl ;
612 srcType = 1 ;
613 }
614 else if ( strncmp( scanType, "ZERO", 4 ) == 0 ) {
615 //cout << "ZERO srcType: " << i << endl ;
616 srcType = 2 ;
617 }
618 else {
619 //cout << "Undefined srcType: " << i << endl ;
620 srcType = 3 ;
621 }
622
623 // if srcType is 2 (ZERO scan), ignore scan
624 if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
625 TableRow row( table_->table() ) ;
626 TableRecord& rec = row.record();
627
628 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
629 *srctCol = srcType ;
630 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
631 Array<Double> srcarr( IPosition( 1, 2 ) ) ;
632 srcarr = 0.0 ;
633 *srateCol = srcarr ;
634 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
635 *spmCol = srcarr ;
636 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
637 *sdirCol = nreader_->getSourceDirection() ;
638 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
639 //*svelCol = h->getURVEL() ; // [m/s]
640 *svelCol = h->getURVEL() + d->VRAD ; // [m/s]
641 RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
642 *fitCol = -1 ;
643 RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
644 //*scannoCol = (uInt)(d->getISCAN()) ;
645 *scannoCol = (uInt)(d->ISCAN) ;
646 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
647 RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
648 *mjdCol = Double( nreader_->getStartIntTime( i ) ) ;
649 RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
650 *intervalCol = Double( h->getIPTIM() ) ;
651 RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
652 *srcnCol = String( h->getOBJ() ) ;
653 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
654 *fieldnCol = String( h->getOBJ() ) ;
655 // BEAMNO is 0-base
656 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
657 string arryt = string( d->ARRYT ) ;
658 string sbeamno = arryt.substr( 1, arryt.size()-1 ) ;
659 uInt ibeamno = atoi( sbeamno.c_str() ) ;
660 *beamCol = ibeamno - 1 ;
661 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
662 RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
663 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
664 vector<double> fqs = nreader_->getFrequencies( i ) ;
665 //cout << "STFiller::readNRO() fqs[1] = " << fqs[1] << endl ;
666 if ( freqs.size() == 0 ) {
667 id = table_->frequencies().addEntry( Double( fqs[0] ),
668 Double( fqs[1] ),
669 Double( fqs[2] ) ) ;
670 *mfreqidCol = id ;
671 *ifCol = id ;
672 freqs.push_back( fqs ) ;
673 }
674 else {
675 int iadd = -1 ;
676 for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
677 //cout << "STFiller::readNRO() freqs[" << iif << "][1] = " << freqs[iif][1] << endl ;
678 double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
679 //cout << "STFiller::readNRO() fdiff = " << fdiff << endl ;
680 if ( fdiff < 1.0e-8 ) {
681 iadd = iif ;
682 break ;
683 }
684 }
685 if ( iadd == -1 ) {
686 id = table_->frequencies().addEntry( Double( fqs[0] ),
687 Double( fqs[1] ),
688 Double( fqs[2] ) ) ;
689 *mfreqidCol = id ;
690 *ifCol = id ;
691 freqs.push_back( fqs ) ;
692 }
693 else {
694 *mfreqidCol = iadd ;
695 *ifCol = iadd ;
696 }
697 }
698 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
699 Vector<Double> restfreq( IPosition( 1, 1 ) ) ;
700 restfreq( 0 ) = d->FREQ0 ;
701 id = table_->molecules().addEntry( restfreq ) ;
702 *molidCol = id ;
703 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
704 //
705 // No Tcal information in the data
706 //
707 // 2008/11/20 Takeshi Nakazato
708 //
709 *mcalidCol = 0 ;
710 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
711
712 id = table_->weather().addEntry( Float( d->TEMP ),
713 Float( d->PATM ),
714 Float( d->PH2O ),
715 Float( d->VWIND ),
716 Float( d->DWIND ) ) ;
717
718 *mweatheridCol = id ;
719 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
720 RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
721 *dirCol = nreader_->getDirection( i ) ;
722 RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
723 *azCol = d->RAZ ;
724 RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
725 *elCol = d->REL ;
726 RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
727 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
728 vector<double> spec = nreader_->getSpectrum( i ) ;
729 Array<Float> sp( IPosition( 1, spec.size() ) ) ;
730 int index = 0 ;
731 for ( Array<Float>::iterator itr = sp.begin() ; itr != sp.end() ; itr++ ) {
732 *itr = spec[index++] ;
733 }
734 *specCol = sp ;
735 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
736 Array<uChar> flag( sp.shape() ) ;
737 flag.set( 0 ) ;
738 *flagCol = flag ;
739 RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
740 *polnoCol = 0 ;
741 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
742 Array<Float> tsys( IPosition( 1, 1 ), d->TSYS ) ;
743 *tsysCol = tsys ;
744
745 table_->table().addRow() ;
746 row.put(table_->table().nrow()-1, rec) ;
747 }
748 else {
749 count++ ;
750 }
751 // DEBUG
752 //int rownum = nreader_->getRowNum() ;
753 //cout << "STFiller::readNRO() Finished row " << i << "/" << rownum << endl ;
754 //
755 }
756
757 // DEBUG
758 time_t t1 ;
759 time( &t1 ) ;
760 ttm = localtime( &t1 ) ;
761 cout << "STFiller::readNRO() Processed " << i << " rows" << endl ;
762 cout << "STFiller::readNRO() Added " << i - count << " rows (ignored "
763 << count << " \"ZERO\" scans)" << endl ;
764 cout << "STFiller::readNRO() End time = " << t1
765 << " ("
766 << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
767 << " "
768 << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
769 << ")" << endl ;
770 cout << "STFiller::readNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
771 //
772
773 return 0 ;
774}
775
776Bool STFiller::fileCheck()
777{
778 // if filename_ is directory, return false
779 File inFile( filename_ ) ;
780 if ( inFile.isDirectory() )
781 return false ;
782
783 // if beginning of header data is "RW", return true
784 // otherwise, return false ;
785 FILE *fp = fopen( filename_.c_str(), "r" ) ;
786 char buf[9] ;
787 fread( buf, 4, 1, fp ) ;
788 buf[4] = '\0' ;
789 if ( ( strncmp( buf, "RW", 2 ) == 0 ) )
790 return true ;
791
792 return false ;
793}
794
795}//namespace asap
Note: See TracBrowser for help on using the repository browser.