Context Navigation

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

Visit:

Last change on this file since 1657 was 1657, checked in by WataruKawasaki, 15 years ago

New Development: Yes

JIRA Issue: Yes (CAS-1433)

Ready to Release: Yes

Interface Changes: No

Put in Release Notes: No

Module(s): sd.save()

Description: Modified to read/write row-based flagging info.

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

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

Download in other formats: