Context Navigation

STFiller.cpp @ 1671

Visit:

Last change on this file since 1671 was 1671, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: Yes CAS-1799

Ready to Release: No

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...

Wrong spectral window (IF) setting for ALMA TP (auto-correlation) data
were fixed.

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
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	bool isGBTFITS = false ;
315	if ((header_->antennaname.find( "GBT" ) != String::npos) && File(filename_).isRegular()) {
316	FILE *fp = fopen( filename_.c_str(), "r" ) ;
317	fseek( fp, 640, SEEK_SET ) ;
318	char buf[81] ;
319	fread( buf, 80, 1, fp ) ;
320	buf[80] = '\0' ;
321	if ( strstr( buf, "NRAO_GBT" ) != NULL ) {
322	isGBTFITS = true ;
323	}
324	fclose( fp ) ;
325	}
326	while ( status == 0 ) {
327	status = reader_->read(pksrec);
328	if ( status != 0 ) break;
329	n += 1;
330
331	Regex filterrx(".*[SL\|PA]$");
332	Regex obsrx("^AT.+");
333	if ( header_->antennaname.matches(obsrx) &&
334	pksrec.obsType.matches(filterrx)) {
335	//cerr << "ignoring paddle scan" << endl;
336	continue;
337	}
338	TableRow row(table_->table());
339	TableRecord& rec = row.record();
340	// fields that don't get used and are just passed through asap
341	RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
342	// MRC changed type from double to float
343	Vector<Double> sratedbl(pksrec.scanRate.nelements());
344	convertArray(sratedbl, pksrec.scanRate);
345	*srateCol = sratedbl;
346	RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
347	*spmCol = pksrec.srcPM;
348	RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
349	*sdirCol = pksrec.srcDir;
350	RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
351	*svelCol = pksrec.srcVel;
352	// the real stuff
353	RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
354	*fitCol = -1;
355	RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
356	*scanoCol = pksrec.scanNo-1;
357	RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
358	*cyclenoCol = pksrec.cycleNo-1;
359	RecordFieldPtr<Double> mjdCol(rec, "TIME");
360	*mjdCol = pksrec.mjd;
361	RecordFieldPtr<Double> intCol(rec, "INTERVAL");
362	*intCol = pksrec.interval;
363	RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
364	RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
365	RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
366	*fieldnCol = pksrec.fieldName;
367	// try to auto-identify if it is on or off.
368	Regex rx(refRx_);
369	Regex rx2("_S$");
370	Int match = pksrec.srcName.matches(rx);
371	if (match) {
372	*srcnCol = pksrec.srcName;
373	} else {
374	*srcnCol = pksrec.srcName.before(rx2);
375	}
376	//*srcnCol = pksrec.srcName;//.before(rx2);
377	*srctCol = match;
378	RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
379	*beamCol = pksrec.beamNo-beamOffset_-1;
380	RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
381	Int rb = -1;
382	if (nBeam_ > 1 ) rb = pksrec.refBeam-1;
383	*rbCol = rb;
384	RecordFieldPtr<uInt> ifCol(rec, "IFNO");
385	//*ifCol = pksrec.IFno-ifOffset_- 1;
386	uInt id;
387	/// @todo this has to change when nchan isn't global anymore
388	//id = table_->frequencies().addEntry(Double(header_->nchan/2),
389	// pksrec.refFreq, pksrec.freqInc);
390	if ( pksrec.nchan == 1 ) {
391	id = table_->frequencies().addEntry(Double(0),
392	pksrec.refFreq, pksrec.freqInc);
393	}
394	else {
395	id = table_->frequencies().addEntry(Double(pksrec.nchan/2),
396	pksrec.refFreq, pksrec.freqInc);
397	}
398	RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
399	*mfreqidCol = id;
400	*ifCol = id;
401
402	id = table_->molecules().addEntry(pksrec.restFreq);
403	RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
404	*molidCol = id;
405
406	id = table_->tcal().addEntry(pksrec.tcalTime, pksrec.tcal);
407	RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
408	*mcalidCol = id;
409	id = table_->weather().addEntry(pksrec.temperature, pksrec.pressure,
410	pksrec.humidity, pksrec.windSpeed,
411	pksrec.windAz);
412	RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
413	*mweatheridCol = id;
414	RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
415	id = table_->focus().addEntry(pksrec.focusAxi, pksrec.focusTan,
416	pksrec.focusRot);
417	*mfocusidCol = id;
418	RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
419	*dirCol = pksrec.direction;
420	RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
421	*azCol = pksrec.azimuth;
422	RecordFieldPtr<Float> elCol(rec, "ELEVATION");
423	*elCol = pksrec.elevation;
424
425	RecordFieldPtr<Float> parCol(rec, "PARANGLE");
426	*parCol = pksrec.parAngle;
427
428	RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
429	RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
430	RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
431
432	RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
433	// Turn the (nchan,npol) matrix and possible complex xPol vector
434	// into 2-4 rows in the scantable
435	Vector<Float> tsysvec(1);
436	// Why is spectra.ncolumn() == 3 for haveXPol_ == True
437	uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
438	for ( uInt i=0; i< npol; ++i ) {
439	tsysvec = pksrec.tsys(i);
440	*tsysCol = tsysvec;
441	if (isGBTFITS)
442	*polnoCol = pksrec.polNo ;
443	else
444	*polnoCol = i;
445
446	*specCol = pksrec.spectra.column(i);
447	*flagCol = pksrec.flagged.column(i);
448	table_->table().addRow();
449	row.put(table_->table().nrow()-1, rec);
450	}
451
452	RecordFieldPtr< uInt > flagrowCol(rec, "FLAGROW");
453	*flagrowCol = pksrec.flagrow;
454
455	if ( haveXPol_[0] ) {
456	// no tsys given for xpol, so emulate it
457	tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
458	*tsysCol = tsysvec;
459	// add real part of cross pol
460	*polnoCol = 2;
461	Vector<Float> r(real(pksrec.xPol));
462	*specCol = r;
463	// make up flags from linears
464	/// @fixme this has to be a bitwise or of both pols
465	*flagCol = pksrec.flagged.column(0);// \| pksrec.flagged.column(1);
466	table_->table().addRow();
467	row.put(table_->table().nrow()-1, rec);
468	// ad imaginary part of cross pol
469	*polnoCol = 3;
470	Vector<Float> im(imag(pksrec.xPol));
471	*specCol = im;
472	table_->table().addRow();
473	row.put(table_->table().nrow()-1, rec);
474	}
475	#ifdef HAS_ALMA
476	fillpm._update(n);
477	#endif
478	}
479	if (status > 0) {
480	close();
481	throw(AipsError("Reading error occured, data possibly corrupted."));
482	}
483	#ifdef HAS_ALMA
484	fillpm.done();
485	#endif
486	return status;
487	}
488
489	/**
490	* For NRO data
491	*
492	* 2008/11/11 Takeshi Nakazato
493	**/
494	void STFiller::openNRO( int whichIF, int whichBeam )
495	{
496	// open file
497	// DEBUG
498	time_t t0 ;
499	time( &t0 ) ;
500	tm *ttm = localtime( &t0 ) ;
501	LogIO os( LogOrigin( "STFiller", "openNRO()", WHERE ) ) ;
502	// cout << "STFiller::openNRO() Start time = " << t0
503	// << " ("
504	// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
505	// << " "
506	// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
507	// << ")" << endl ;
508	os << "Start time = " << t0
509	<< " ("
510	<< ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
511	<< " "
512	<< ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
513	<< ")" << LogIO::POST ;
514
515	// fill STHeader
516	header_ = new STHeader() ;
517
518	if ( nreader_->getHeaderInfo( header_->nchan,
519	header_->npol,
520	nIF_,
521	nBeam_,
522	header_->observer,
523	header_->project,
524	header_->obstype,
525	header_->antennaname,
526	header_->antennaposition,
527	header_->equinox,
528	header_->freqref,
529	header_->reffreq,
530	header_->bandwidth,
531	header_->utc,
532	header_->fluxunit,
533	header_->epoch,
534	header_->poltype ) ) {
535	// cout << "STFiller::openNRO() Failed to get header information." << endl ;
536	// return ;
537	throw( AipsError("Failed to get header information.") ) ;
538	}
539
540	// set frame keyword of FREQUENCIES table
541	if ( header_->freqref != "TOPO" ) {
542	table_->frequencies().setFrame( header_->freqref, false ) ;
543	}
544
545	ifOffset_ = 0;
546	vector<Bool> ifs = nreader_->getIFs() ;
547	if ( whichIF >= 0 ) {
548	if ( whichIF >= 0 && whichIF < nIF_ ) {
549	for ( int i = 0 ; i < nIF_ ; i++ )
550	ifs[i] = False ;
551	ifs[whichIF] = True ;
552	header_->nif = 1;
553	nIF_ = 1;
554	ifOffset_ = whichIF;
555	} else {
556	delete reader_;
557	reader_ = 0;
558	delete header_;
559	header_ = 0;
560	throw(AipsError("Illegal IF selection"));
561	}
562	}
563
564	beamOffset_ = 0;
565	vector<Bool> beams = nreader_->getBeams() ;
566	if (whichBeam>=0) {
567	if (whichBeam>=0 && whichBeam<nBeam_) {
568	for ( int i = 0 ; i < nBeam_ ; i++ )
569	beams[i] = False ;
570	beams[whichBeam] = True;
571	header_->nbeam = 1;
572	nBeam_ = 1;
573	beamOffset_ = whichBeam;
574	} else {
575	delete reader_;
576	reader_ = 0;
577	delete header_;
578	header_ = 0;
579	throw(AipsError("Illegal Beam selection"));
580	}
581	}
582
583	header_->nbeam = nBeam_ ;
584	header_->nif = nIF_ ;
585
586	// set header
587	table_->setHeader( *header_ ) ;
588
589	// DEBUG
590	time_t t1 ;
591	time( &t1 ) ;
592	ttm = localtime( &t1 ) ;
593	// cout << "STFiller::openNRO() End time = " << t1
594	// << " ("
595	// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
596	// << " "
597	// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
598	// << ")" << endl ;
599	// cout << "STFiller::openNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
600	os << "End time = " << t1
601	<< " ("
602	<< ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
603	<< " "
604	<< ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
605	<< ")" << endl ;
606	os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
607	os.post() ;
608	//
609
610	return ;
611	}
612
613	int STFiller::readNRO()
614	{
615	// DEBUG
616	time_t t0 ;
617	time( &t0 ) ;
618	tm *ttm = localtime( &t0 ) ;
619	LogIO os( LogOrigin( "STFiller", "readNRO()", WHERE ) ) ;
620	// cout << "STFiller::readNRO() Start time = " << t0
621	// << " ("
622	// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
623	// << " "
624	// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
625	// << ")" << endl ;
626	os << "Start time = " << t0
627	<< " ("
628	<< ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
629	<< " "
630	<< ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
631	<< ")" << LogIO::POST ;
632	//
633
634	// fill row
635	uInt id ;
636	uInt imax = nreader_->getRowNum() ;
637	vector< vector<double > > freqs ;
638	uInt i = 0 ;
639	int count = 0 ;
640	uInt scanno ;
641	uInt cycleno ;
642	uInt beamno ;
643	uInt polno ;
644	vector<double> fqs ;
645	Vector<Double> restfreq ;
646	uInt refbeamno ;
647	Double scantime ;
648	Double interval ;
649	String srcname ;
650	String fieldname ;
651	Array<Float> spectra ;
652	Array<uChar> flagtra ;
653	Array<Float> tsys ;
654	Array<Double> direction ;
655	Float azimuth ;
656	Float elevation ;
657	Float parangle ;
658	Float opacity ;
659	uInt tcalid ;
660	Int fitid ;
661	uInt focusid ;
662	Float temperature ;
663	Float pressure ;
664	Float humidity ;
665	Float windvel ;
666	Float winddir ;
667	Double srcvel ;
668	Array<Double> propermotion ;
669	Vector<Double> srcdir ;
670	Array<Double> scanrate ;
671	for ( i = 0 ; i < imax ; i++ ) {
672	string scanType = nreader_->getScanType( i ) ;
673	Int srcType = -1 ;
674	if ( scanType.compare( 0, 2, "ON") == 0 ) {
675	// os << "ON srcType: " << i << LogIO::POST ;
676	srcType = 0 ;
677	}
678	else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
679	//os << "OFF srcType: " << i << LogIO::POST ;
680	srcType = 1 ;
681	}
682	else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
683	//os << "ZERO srcType: " << i << LogIO::POST ;
684	srcType = 2 ;
685	}
686	else {
687	//os << "Undefined srcType: " << i << LogIO::POST ;
688	srcType = 3 ;
689	}
690
691	// if srcType is 2 (ZERO scan), ignore scan
692	if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
693	TableRow row( table_->table() ) ;
694	TableRecord& rec = row.record();
695
696	if ( nreader_->getScanInfo( i,
697	scanno,
698	cycleno,
699	beamno,
700	polno,
701	fqs,
702	restfreq,
703	refbeamno,
704	scantime,
705	interval,
706	srcname,
707	fieldname,
708	spectra,
709	flagtra,
710	tsys,
711	direction,
712	azimuth,
713	elevation,
714	parangle,
715	opacity,
716	tcalid,
717	fitid,
718	focusid,
719	temperature,
720	pressure,
721	humidity,
722	windvel,
723	winddir,
724	srcvel,
725	propermotion,
726	srcdir,
727	scanrate ) ) {
728	// cerr << "STFiller::readNRO() Failed to get scan information." << endl ;
729	// return 1 ;
730	throw( AipsError("Failed to get scan information.") ) ;
731	}
732
733	RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
734	*scannoCol = scanno ;
735	RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ;
736	*cyclenoCol = cycleno ;
737	RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
738	*beamCol = beamno ;
739	RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
740	RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
741	*polnoCol = polno ;
742	RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
743	if ( freqs.size() == 0 ) {
744	id = table_->frequencies().addEntry( Double( fqs[0] ),
745	Double( fqs[1] ),
746	Double( fqs[2] ) ) ;
747	*mfreqidCol = id ;
748	*ifCol = id ;
749	freqs.push_back( fqs ) ;
750	}
751	else {
752	int iadd = -1 ;
753	for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
754	//os << "freqs[" << iif << "][1] = " << freqs[iif][1] << LogIO::POST ;
755	double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
756	//os << "fdiff = " << fdiff << LogIO::POST ;
757	if ( fdiff < 1.0e-8 ) {
758	iadd = iif ;
759	break ;
760	}
761	}
762	if ( iadd == -1 ) {
763	id = table_->frequencies().addEntry( Double( fqs[0] ),
764	Double( fqs[1] ),
765	Double( fqs[2] ) ) ;
766	*mfreqidCol = id ;
767	*ifCol = id ;
768	freqs.push_back( fqs ) ;
769	}
770	else {
771	*mfreqidCol = iadd ;
772	*ifCol = iadd ;
773	}
774	}
775	RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
776	id = table_->molecules().addEntry( restfreq ) ;
777	*molidCol = id ;
778	RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
779	*rbCol = refbeamno ;
780	RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
781	*mjdCol = scantime ;
782	RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
783	*intervalCol = interval ;
784	RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
785	*srcnCol = srcname ;
786	RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
787	*srctCol = srcType ;
788	RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
789	*fieldnCol = fieldname ;
790	RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
791	*specCol = spectra ;
792	RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
793	*flagCol = flagtra ;
794	RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
795	*tsysCol = tsys ;
796	RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
797	*dirCol = direction ;
798	RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
799	*azCol = azimuth ;
800	RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
801	*elCol = elevation ;
802	RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
803	*parCol = parangle ;
804	RecordFieldPtr<Float> tauCol(rec, "OPACITY") ;
805	*tauCol = opacity ;
806	RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
807	*mcalidCol = tcalid ;
808	RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
809	*fitCol = fitid ;
810	RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
811	*mfocusidCol = focusid ;
812	RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
813	id = table_->weather().addEntry( temperature,
814	pressure,
815	humidity,
816	windvel,
817	winddir ) ;
818	*mweatheridCol = id ;
819	RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
820	*svelCol = srcvel ;
821	RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
822	*spmCol = propermotion ;
823	RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
824	*sdirCol = srcdir ;
825	RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
826	*srateCol = scanrate ;
827
828	table_->table().addRow() ;
829	row.put(table_->table().nrow()-1, rec) ;
830	}
831	else {
832	count++ ;
833	}
834	// DEBUG
835	//int rownum = nreader_->getRowNum() ;
836	//os << "Finished row " << i << "/" << rownum << LogIO::POST ;
837	//
838	}
839
840	// DEBUG
841	time_t t1 ;
842	time( &t1 ) ;
843	ttm = localtime( &t1 ) ;
844	// cout << "STFiller::readNRO() Processed " << i << " rows" << endl ;
845	// cout << "STFiller::readNRO() Added " << i - count << " rows (ignored "
846	// << count << " \"ZERO\" scans)" << endl ;
847	// cout << "STFiller::readNRO() End time = " << t1
848	// << " ("
849	// << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
850	// << " "
851	// << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
852	// << ")" << endl ;
853	// cout << "STFiller::readNRO() Elapsed time = " << t1 - t0 << " sec" << endl ;
854	os << "Processed " << i << " rows" << endl ;
855	os << "Added " << i - count << " rows (ignored "
856	<< count << " \"ZERO\" scans)" << endl ;
857	os.post() ;
858	os << "End time = " << t1
859	<< " ("
860	<< ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
861	<< " "
862	<< ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
863	<< ")" << endl ;
864	os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
865	os.post() ;
866	//
867
868	return 0 ;
869	}
870
871	Bool STFiller::fileCheck()
872	{
873	bool bval = false ;
874
875	// if filename_ is directory, return false
876	File inFile( filename_ ) ;
877	if ( inFile.isDirectory() )
878	return bval ;
879
880	// if beginning of header data is "RW", return true
881	// otherwise, return false ;
882	FILE *fp = fopen( filename_.c_str(), "r" ) ;
883	char buf[9] ;
884	char buf2[80] ;
885	fread( buf, 4, 1, fp ) ;
886	buf[4] = '\0' ;
887	fseek( fp, 640, SEEK_SET ) ;
888	fread( buf2, 80, 1, fp ) ;
889	if ( ( strncmp( buf, "RW", 2 ) == 0 ) \|\| ( strstr( buf2, "NRO45M" ) != NULL ) ) {
890	bval = true ;
891	}
892	fclose( fp ) ;
893	return bval ;
894	}
895
896	}//namespace asap

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

Context Navigation

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

Download in other formats: