Context Navigation

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

Visit:

Last change on this file since 1491 was 1491, 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: Fixed bug in frequency conversion again.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: