Context Navigation

STFiller.cpp @ 1668

Visit:

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

New Development: No

JIRA Issue: Yes CAS-1683

Ready to Release: No

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Modification to support GBT SDFITS data.
IFNO is same as FREQ_ID.

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

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

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

Context Navigation

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

Download in other formats: