Context Navigation

STFiller.cpp@ 1674

Visit:

Last change on this file since 1674 was 1674, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-1799

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: srcType is added to PKSrecord class

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): casa/atnf

Description: Describe your changes here...

Calibration scheme is updated to use SRCTYPE as an identifier of scan intent
(ON, OFF, ...). Currently, updated calibration scheme is implemented only for
ALMA data. Calibration of other data is still use SRCNAME as the identifier.

casa/atnf must be updated.

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

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

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

Context Navigation

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

Download in other formats: