Context Navigation

STFiller.cpp @ 1616

Visit:

Last change on this file since 1616 was 1616, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-729, CAS-1147

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: Describe your changes here...

Use LogIO instead of std::cout and std::cerr.

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

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

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

Context Navigation

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

Download in other formats: