Context Navigation

source: trunk/src/STFiller.cpp @ 2657

Visit:

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

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

Download in other formats: