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

Last change on this file since 1502 was 1495, checked in by TakTsutsumi, 16 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: a new parameter in open method

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Added a boolean, getPt in open()

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