Context Navigation

STFiller.cpp @ 2262

Visit:

Last change on this file since 2262 was 2262, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: No

Ready for Test: 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...

Preparation for check in tuning on NRO filler.
Merge changes in trunk.

r2154
r2156
r2158
r2198-2203
r2212
r2261

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

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

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

Context Navigation

source: branches/parallel/src/STFiller.cpp @ 2262

Download in other formats: