Context Navigation

source: trunk/src/STFiller.cpp @ 2658

Visit:

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

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

Download in other formats: