Context Navigation

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

Visit:

Last change on this file since 1603 was 1603, checked in by TakTsutsumi, 15 years ago

New Development: No, merge with asap2.3.1

JIRA Issue: Yes CAS-1450

Ready to Release: Yes/No?

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): single dish

Description: Upgrade of alma branch based on ASAP2.2.0

(rev.1562) to ASAP2.3.1 (rev.1561)

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

File size: 24.8 KB

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

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

Download in other formats: