Context Navigation

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

Visit:

Last change on this file since 1490 was 1490, checked in by Takeshi Nakazato, 15 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
42	using namespace casa;
43
44	namespace asap {
45
46	STFiller::STFiller() :
47	reader_(0),
48	header_(0),
49	table_(0),
50	nreader_(0)
51	{
52	}
53
54	STFiller::STFiller( CountedPtr< Scantable > stbl ) :
55	reader_(0),
56	header_(0),
57	table_(stbl),
58	nreader_(0)
59	{
60	}
61
62	STFiller::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
71	STFiller::~STFiller()
72	{
73	close();
74	}
75
76	void 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
255	void STFiller::close( )
256	{
257	delete reader_;reader_=0;
258	delete nreader_;nreader_=0;
259	delete header_;header_=0;
260	table_ = 0;
261	}
262
263	int 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	**/
443	void 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
572	int 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
776	Bool 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.

Download in other formats: