Context Navigation

NROReader.cc @ 3106

Visit:

Last change on this file since 3106 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No?

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Check-in asap modifications from Jim regarding casacore namespace conversion.

File size: 24.2 KB

Line
1	//#---------------------------------------------------------------------------
2	//# NROReader.cc: Base class for NRO headerdata.
3	//#---------------------------------------------------------------------------
4	//# Copyright (C) 2000-2006
5	//# Associated Universities, Inc. Washington DC, USA.
6	//#
7	//# This library is free software; you can redistribute it and/or modify it
8	//# under the terms of the GNU Library General Public License as published by
9	//# the Free Software Foundation; either version 2 of the License, or (at your
10	//# option) any later version.
11	//#
12	//# This library is distributed in the hope that it will be useful, but WITHOUT
13	//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14	//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15	//# License for more details.
16	//#
17	//# You should have received a copy of the GNU Library General Public License
18	//# along with this library; if not, write to the Free Software Foundation,
19	//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20	//#
21	//# Correspondence concerning AIPS++ should be addressed as follows:
22	//# Internet email: aips2-request@nrao.edu.
23	//# Postal address: AIPS++ Project Office
24	//# National Radio Astronomy Observatory
25	//# 520 Edgemont Road
26	//# Charlottesville, VA 22903-2475 USA
27	//#
28	//# $Id$
29	//#---------------------------------------------------------------------------
30	//# Original: 2008/10/30, Takeshi Nakazato, NAOJ
31	//#---------------------------------------------------------------------------
32
33	#include <atnf/PKSIO/NROReader.h>
34	#include <atnf/PKSIO/NRO45Reader.h>
35	#include <atnf/PKSIO/ASTEReader.h>
36	#include <atnf/PKSIO/ASTEFXReader.h>
37	#include <atnf/PKSIO/NRO45FITSReader.h>
38	#include <atnf/PKSIO/NROOTFDataset.h>
39	#include <atnf/PKSIO/ASTEDataset.h>
40
41	#include <measures/Measures/MEpoch.h>
42	#include <measures/Measures/MPosition.h>
43	#include <measures/Measures/MDirection.h>
44	#include <measures/Measures/MCDirection.h>
45	#include <measures/Measures/MeasFrame.h>
46	#include <measures/Measures/MeasConvert.h>
47
48	#include <casa/Logging/LogIO.h>
49	#include <casa/IO/RegularFileIO.h>
50	#include <casa/OS/File.h>
51	#include <casa/OS/Time.h>
52
53	#include <stdio.h>
54	#include <string>
55	#include <iomanip>
56	#include <cassert>
57
58	using namespace std ;
59
60	namespace {
61	// to limit template argument to specified type
62	template<class T>
63	inline bool type_guard() {
64	return false;
65	}
66	template<>
67	inline bool type_guard<string>() {
68	return true;
69	}
70	template<>
71	inline bool type_guard<String>() {
72	return true;
73	}
74
75	// T must be std::string or casacore::String
76	template<class T>
77	inline String trim_nro_string(T const &s) {
78	assert(type_guard<T>());
79	return s.substr(0, s.find_first_of('\0'));
80	}
81	}
82
83	//
84	// getNROReader
85	//
86	// Return an appropriate NROReader for a NRO 45m and ASTE dataset.
87	//
88	NROReader *getNROReader( const String filename,
89	String &datatype )
90	{
91	LogIO os( LogOrigin( "", "getNROReader()", WHERE ) ) ;
92
93	// Check accessibility of the input.
94	File inFile( filename ) ;
95	if ( !inFile.exists() ) {
96	datatype = filename + " not found." ;
97	return 0 ;
98	}
99
100	if ( !inFile.isReadable() ) {
101	datatype = filename + " is not readable." ;
102	return 0 ;
103	}
104
105	// Determine the type of input.
106	NROReader *reader = 0;
107	if ( inFile.isRegular() ) {
108	FILE *file ;
109	file = fopen( filename.c_str(), "r" ) ;
110	// read LOFIL0
111	char buf[9];
112	size_t retval = fread( buf, 4, 1, file ) ;
113	if (retval < 1) {
114	os << LogIO::SEVERE << "Failed to read data" << LogIO::POST;
115	fclose(file);
116	return 0;
117	}
118	buf[4] = '\0' ;
119	// DEBUG
120	//os << LogIO::NORMAL << "getNROReader:: buf = " << String(buf) << LogIO::POST ;
121	//
122	if ( string( buf ) == "XTEN" ) {
123	// FITS data
124	datatype = "NRO 45m FITS" ;
125	reader = new NRO45FITSReader( filename ) ;
126	}
127	else if ( string( buf ) == "RW-F") {
128	// ASTE-FX data
129	datatype = "ASTE-FX";
130	reader = new ASTEFXReader( filename );
131	} else {
132	// otherwise, read SITE0
133	NRODataset *d = new NROOTFDataset( filename ) ;
134	d->initialize() ;
135	int size = d->getDataSize() - 188 ;
136	delete d ;
137	fseek( file, size, SEEK_SET ) ;
138	size_t retval = fread( buf, 8, 1, file ) ;
139	if (retval < 1) {
140	os << LogIO::SEVERE << "Failed to read data" << LogIO::POST;
141	fclose(file);
142	return 0;
143	}
144	buf[8] = '\0' ;
145	// DEBUG
146	//cout << "getNROReader:: buf = " << buf << endl ;
147	//
148	if ( string( buf ) == "NRO" ) {
149	// NRO 45m data
150	datatype = "NRO 45m OTF" ;
151	reader = new NRO45Reader( filename );
152	}
153	else {
154	d = new ASTEDataset( filename ) ;
155	d->initialize() ;
156	size = d->getDataSize() - 188 ;
157	delete d ;
158	fseek( file, size, SEEK_SET ) ;
159	size_t retval = fread( buf, 8, 1, file ) ;
160	if (retval < 1) {
161	os << LogIO::SEVERE << "Failed to read data" << LogIO::POST;
162	fclose(file);
163	return 0;
164	}
165	buf[8] = '\0' ;
166	// DEBUG
167	//cout << "getNROReader:: buf = " << buf << endl ;
168	//
169	if ( string( buf ) == "ASTE" ) {
170	// ASTE data
171	datatype = "ASTE" ;
172	reader = new ASTEReader( filename ) ;
173	}
174	else {
175	datatype = "UNRECOGNIZED INPUT FORMAT";
176	}
177	}
178	}
179	fclose( file ) ;
180	} else {
181	datatype = "UNRECOGNIZED INPUT FORMAT";
182	}
183
184	// DEBUG
185	os << LogIO::NORMAL << "Data format of " << filename << ": " << datatype << LogIO::POST ;
186	//
187
188	// return reader if exists
189	if ( reader ) {
190	reader->read() ;
191	return reader ;
192	}
193
194	return 0 ;
195	}
196
197
198	//
199	// getNROReader
200	//
201	// Search a list of directories for a NRO 45m and ASTE dataset and return an
202	// appropriate NROReader for it.
203	//
204	NROReader* getNROReader( const String name,
205	const Vector<String> directories,
206	int &iDir,
207	String &datatype )
208	{
209	int nDir = directories.size();
210	for ( iDir = 0; iDir < nDir; iDir++ ) {
211	string inName = directories[iDir] + "/" + name;
212	NROReader *reader = getNROReader( inName, datatype ) ;
213
214	if (reader != 0) {
215	return reader;
216	}
217	}
218
219	iDir = -1;
220	return 0;
221	}
222
223
224	//----------------------------------------------------------------------
225	// constructor
226	NROReader::NROReader( string name )
227	: dataset_( NULL ),
228	srcdir_( 0 ),
229	msrcdir_( 0 ),
230	freqRefFromVREF_( false ) // import frequency as REST
231	{
232	// initialization
233	filename_ = name ;
234	coord_ = -1 ;
235	}
236
237	// destructor
238	NROReader::~NROReader()
239	{
240	}
241
242	// Read data header
243	Int NROReader::read()
244	{
245	LogIO os( LogOrigin( "NROReader", "read()", WHERE ) ) ;
246
247	int status = 0 ;
248
249	// initialize Dataset object
250	initDataset();
251
252	// fill Dataset
253	status = dataset_->fillHeader() ;
254
255	if ( status != 0 ) {
256	os << LogIO::SEVERE << "Failed to fill data header." << LogIO::EXCEPTION ;
257	}
258
259	return status ;
260	}
261
262	// set frequency reference frame
263	void NROReader::setFreqRefFromVREF( bool fromVREF )
264	{
265	os_.origin( LogOrigin( "NROReader", "setFreqRefFromVREF", WHERE ) ) ;
266	os_ << ((fromVREF) ? "Take frequency reference frame from VREF" :
267	"Use frequency reference frame REST") << LogIO::POST ;
268
269	freqRefFromVREF_ = fromVREF ;
270	}
271
272	// get spectrum
273	vector< vector<double> > NROReader::getSpectrum()
274	{
275	return dataset_->getSpectrum() ;
276	}
277
278	Int NROReader::getPolarizationNum()
279	{
280	return dataset_->getPolarizationNum() ;
281	}
282
283	double NROReader::getStartTime()
284	{
285	//char *startTime = dataset_->getLOSTM() ;
286	string startTime = dataset_->getLOSTM() ;
287	//cout << "getStartTime() startTime = " << startTime << endl ;
288	return getMJD( startTime ) ;
289	}
290
291	double NROReader::getEndTime()
292	{
293	//char *endTime = dataset_->getLOETM() ;
294	string endTime = dataset_->getLOETM() ;
295	return getMJD( endTime ) ;
296	}
297
298	vector<double> NROReader::getStartIntTime()
299	{
300	return dataset_->getStartIntTime() ;
301	}
302
303	double NROReader::getMJD( char *time )
304	{
305	// TODO: should be checked which time zone the time depends on
306	// 2008/11/14 Takeshi Nakazato
307	string strStartTime = string( time ) ;
308	string strYear = strStartTime.substr( 0, 4 ) ;
309	string strMonth = strStartTime.substr( 4, 2 ) ;
310	string strDay = strStartTime.substr( 6, 2 ) ;
311	string strHour = strStartTime.substr( 8, 2 ) ;
312	string strMinute = strStartTime.substr( 10, 2 ) ;
313	string strSecond = strStartTime.substr( 12, strStartTime.size() - 12 ) ;
314	uInt year = atoi( strYear.c_str() ) ;
315	uInt month = atoi( strMonth.c_str() ) ;
316	uInt day = atoi( strDay.c_str() ) ;
317	uInt hour = atoi( strHour.c_str() ) ;
318	uInt minute = atoi( strMinute.c_str() ) ;
319	double second = atof( strSecond.c_str() ) ;
320	Time t( year, month, day, hour, minute, second ) ;
321
322	return t.modifiedJulianDay() ;
323	}
324
325	double NROReader::getMJD( string strStartTime )
326	{
327	// TODO: should be checked which time zone the time depends on
328	// 2008/11/14 Takeshi Nakazato
329	string strYear = strStartTime.substr( 0, 4 ) ;
330	string strMonth = strStartTime.substr( 4, 2 ) ;
331	string strDay = strStartTime.substr( 6, 2 ) ;
332	string strHour = strStartTime.substr( 8, 2 ) ;
333	string strMinute = strStartTime.substr( 10, 2 ) ;
334	string strSecond = strStartTime.substr( 12, strStartTime.size() - 12 ) ;
335	uInt year = atoi( strYear.c_str() ) ;
336	uInt month = atoi( strMonth.c_str() ) ;
337	uInt day = atoi( strDay.c_str() ) ;
338	uInt hour = atoi( strHour.c_str() ) ;
339	uInt minute = atoi( strMinute.c_str() ) ;
340	double second = atof( strSecond.c_str() ) ;
341	Time t( year, month, day, hour, minute, second ) ;
342
343	return t.modifiedJulianDay() ;
344	}
345
346	vector<Bool> NROReader::getIFs()
347	{
348	return dataset_->getIFs() ;
349	}
350
351	vector<Bool> NROReader::getBeams()
352	{
353	vector<Bool> v ;
354	vector<int> arry = dataset_->getARRY() ;
355	for ( uInt i = 0 ; i < arry.size() ; i++ ) {
356	if ( arry[i] != 0 ) {
357	v.push_back( True ) ;
358	}
359	}
360
361	// DEBUG
362	//cout << "NROReader::getBeams() number of beam is " << v.size() << endl ;
363	//
364
365	return v ;
366	}
367
368	// Get SRCDIRECTION in RADEC(J2000)
369	Vector<Double> NROReader::getSourceDirection()
370	{
371	if ( msrcdir_.nelements() == 2 )
372	return msrcdir_ ;
373
374	srcdir_.resize( 2 ) ;
375	srcdir_[0] = Double( dataset_->getRA0() ) ;
376	srcdir_[1] = Double( dataset_->getDEC0() ) ;
377	char epoch[5] ;
378	strncpy( epoch, (dataset_->getEPOCH()).c_str(), 5 ) ;
379	if ( strncmp( epoch, "B1950", 5 ) == 0 ) {
380	// convert to J2000 value
381	MDirection result =
382	MDirection::Convert( MDirection( Quantity( srcdir_[0], "rad" ),
383	Quantity( srcdir_[1], "rad" ),
384	MDirection::Ref( MDirection::B1950 ) ),
385	MDirection::Ref( MDirection::J2000 ) ) () ;
386	msrcdir_ = result.getAngle().getValue() ;
387	if ( msrcdir_[0] < 0.0 && srcdir_[0] >= 0.0 )
388	msrcdir_[0] = 2.0 * M_PI + msrcdir_[0] ;
389	//cout << "NROReader::getSourceDirection() SRCDIRECTION convert from ("
390	//<< srcra << "," << srcdec << ") B1950 to ("
391	//<< v( 0 ) << ","<< v( 1 ) << ") J2000" << endl ;
392	//LogIO os( LogOrigin( "NROReader", "getSourceDirection()", WHERE ) ) ;
393	//os << LogIO::NORMAL << "SRCDIRECTION convert from ("
394	// << srcra << "," << srcdec << ") B1950 to ("
395	// << v( 0 ) << ","<< v( 1 ) << ") J2000" << LogIO::POST ;
396	}
397	else if ( strncmp( epoch, "J2000", 5 ) == 0 ) {
398	msrcdir_.reference( srcdir_ ) ;
399	}
400
401	return msrcdir_ ;
402	}
403
404	// Get DIRECTION in RADEC(J2000)
405	Vector<Double> NROReader::getDirection( int i )
406	{
407	Vector<Double> v( 2 ) ;
408	const NRODataRecord *record = dataset_->getRecord( i ) ;
409	char epoch[5] ;
410	strncpy( epoch, (dataset_->getEPOCH()).c_str(), 5 ) ;
411	int icoord = dataset_->getSCNCD() ;
412	double scantime = dataset_->getScanTime( i ) ;
413	initConvert( icoord, scantime, epoch ) ;
414	v[0]= Double( record->SCX ) ;
415	v[1] = Double( record->SCY ) ;
416	if ( converter_.null() )
417	// no conversion
418	return v ;
419	else {
420	Vector<Double> v2 = (*converter_)( v ).getAngle().getValue() ;
421	if ( v2[0] < 0.0 && v[0] >= 0.0 )
422	v2[0] = 2.0 * M_PI + v2[0] ;
423	return v2 ;
424	}
425	}
426
427	void NROReader::initConvert( int icoord, double t, char *epoch )
428	{
429	if ( icoord == 0 && strncmp( epoch, "J2000", 5 ) == 0 )
430	// no conversion
431	return ;
432
433	if ( converter_.null() \|\| icoord != coord_ ) {
434	LogIO os( LogOrigin( "NROReader", "initConvert()", WHERE ) ) ;
435	coord_ = icoord ;
436	if ( coord_ == 0 ) {
437	// RADEC (B1950) -> RADEC (J2000)
438	os << "Creating converter from RADEC (B1950) to RADEC (J2000)" << LogIO::POST ;
439	converter_ = new MDirection::Convert( MDirection::B1950,
440	MDirection::J2000 ) ;
441	}
442	else if ( coord_ == 1 ) {
443	// LB -> RADEC (J2000)
444	os << "Creating converter from GALACTIC to RADEC (J2000)" << LogIO::POST ;
445	converter_ = new MDirection::Convert( MDirection::GALACTIC,
446	MDirection::J2000 ) ;
447	}
448	else {
449	// coord_ == 2
450	// AZEL -> RADEC (J2000)
451	os << "Creating converter from AZEL to RADEC (J2000)" << LogIO::POST ;
452	if ( mf_.null() ) {
453	mf_ = new MeasFrame() ;
454	vector<double> antpos = getAntennaPosition() ;
455	Vector<Quantity> qantpos( 3 ) ;
456	for ( int ip = 0 ; ip < 3 ; ip++ )
457	qantpos[ip] = Quantity( antpos[ip], "m" ) ;
458	mp_ = MPosition( MVPosition( qantpos ), MPosition::ITRF ) ;
459	mf_->set( mp_ ) ;
460	}
461	converter_ = new MDirection::Convert( MDirection::AZEL,
462	MDirection::Ref(MDirection::J2000,
463	*mf_ ) ) ;
464	}
465	}
466
467	if ( coord_ == 2 ) {
468	me_ = MEpoch( Quantity( t, "d" ), MEpoch::UTC ) ;
469	mf_->set( me_ ) ;
470	}
471	}
472
473	int NROReader::getHeaderInfo( Int &nchan,
474	Int &npol,
475	Int &nif,
476	Int &nbeam,
477	String &observer,
478	String &project,
479	String &obstype,
480	String &antname,
481	Vector<Double> &antpos,
482	Float &equinox,
483	String &freqref,
484	Double &reffreq,
485	Double &bw,
486	Double &utc,
487	String &fluxunit,
488	String &epoch,
489	String &poltype )
490	{
491	nchan = dataset_->getNUMCH() ;
492	//cout << "nchan = " << nchan << endl ;
493	npol = getPolarizationNum() ;
494	//cout << "npol = " << npol << endl ;
495	observer = trim_nro_string(dataset_->getOBSVR()) ;
496	//cout << "observer = \"" << observer << "\" (size " << observer.size() << ")" << endl ;
497	project = trim_nro_string(dataset_->getPROJ()) ;
498	//cout << "project = \"" << project << "\" (size " << project.size() << ")" << endl ;
499	obstype = trim_nro_string(dataset_->getSWMOD()) ;
500	//cout << "obstype = \"" << obstype << "\" (size " << obstype.size() << ")" << endl ;
501	antname = trim_nro_string(dataset_->getSITE()) ;
502	//cout << "antname = \"" << antname << "\" (size " << antname.size() << ")" << endl ;
503	// TODO: should be investigated antenna position since there are
504	// no corresponding information in the header
505	// 2008/11/13 Takeshi Nakazato
506	//
507	// INFO: tentative antenna posiiton is obtained for NRO 45m from ITRF website
508	// 2008/11/26 Takeshi Nakazato
509	vector<double> pos = getAntennaPosition() ;
510	antpos = pos ;
511	//cout << "antpos = " << antpos << endl ;
512	string eq = dataset_->getEPOCH() ;
513	// if ( eq.compare( 0, 5, "B1950" ) == 0 )
514	// equinox = 1950.0 ;
515	// else if ( eq.compare( 0, 5, "J2000" ) == 0 )
516	// equinox = 2000.0 ;
517	// equinox is always 2000.0
518	equinox = 2000.0 ;
519	//cout << "equinox = " << equinox << endl ;
520	string vref = dataset_->getVREF() ;
521	if ( vref.compare( 0, 3, "LSR" ) == 0 ) {
522	if ( vref.size() == 3 ) {
523	vref.append( "K" ) ;
524	}
525	else {
526	vref[3] = 'K' ;
527	}
528	}
529	else if ( vref.compare( 0, 3, "GAL" ) == 0 ) {
530	vref = "GALACTO" ;
531	}
532	else if (vref.compare( 0, 3, "HEL" ) == 0 ) {
533	os_.origin( LogOrigin( "NROReader", "getHeaderInfo", WHERE ) ) ;
534	os_ << LogIO::WARN << "Heliocentric frame is not supported. Use Barycentric frame instead." << LogIO::POST ;
535	vref = "BARY" ;
536	}
537	//freqref = vref ;
538	//freqref = "LSRK" ;
539	//freqref = "REST" ;
540	freqref = freqRefFromVREF_ ? trim_nro_string(vref) : "REST" ;
541	//cout << "freqref = \"" << freqref << "\" (size " << freqref.size() << ")" << endl ;
542	const NRODataRecord *record = dataset_->getRecord( 0 ) ;
543	reffreq = record->FREQ0 ; // rest frequency
544
545	//cout << "reffreq = " << reffreq << endl ;
546	bw = dataset_->getBEBW()[0] ;
547	//cout << "bw = " << bw << endl ;
548	utc = getStartTime() ;
549	//cout << "utc = " << utc << endl ;
550	fluxunit = "K" ;
551	//cout << "fluxunit = " << fluxunit << endl ;
552	epoch = "UTC" ;
553	//cout << "epoch = " << epoch << endl ;
554	string poltp = trim_nro_string(dataset_->getPOLTP()[0]) ;
555	//cout << "poltp = \"" << poltp << "\" (size " << poltp.size() << ")" << endl ;
556	if ( poltp.empty() \|\| poltp[0] == ' ' \|\| poltp[0] == '\0' )
557	//poltp = "None" ;
558	poltp = "linear" ; // if no polarization type specified, set to "linear"
559	//else if ( strcmp( poltp, "LINR" ) == 0 )
560	else if ( poltp.compare( 0, 1, "LINR", 0, 1 ) == 0 )
561	poltp = "linear" ;
562	//else if ( strcmp( poltp, "CIRL" ) == 0 )
563	else if ( poltp.compare( 0, 1, "CIRL", 0, 1 ) == 0 )
564	poltp = "circular" ;
565	poltype = poltp ;
566	//cout << "poltype = '" << poltype << "'" << endl ;
567
568	//vector<Bool> ifs = getIFs() ;
569	//nif = ifs.size() ;
570	nif = getNumIF() ;
571	//cout << "nif = " << nif << endl ;
572
573	//vector<Bool> beams = getBeams() ;
574	//nbeam = beams.size() ;
575	nbeam = getNumBeam() ;
576	//cout << "nbeam = " << nbeam << endl ;
577
578	return 0 ;
579	}
580
581	string NROReader::getScanType( int i )
582	{
583	const NRODataRecord *record = dataset_->getRecord( i ) ;
584	string s = record->SCANTP ;
585
586	return s ;
587	}
588
589	int NROReader::getScanInfo( int irow,
590	uInt &scanno,
591	uInt &cycleno,
592	uInt &ifno,
593	uInt &beamno,
594	uInt &polno,
595	vector<double> &freqs,
596	Vector<Double> &restfreq,
597	uInt &refbeamno,
598	Double &scantime,
599	Double &interval,
600	String &srcname,
601	String &fieldname,
602	Vector<Float> &spectra,
603	Vector<uChar> &flagtra,
604	Vector<Float> &tsys,
605	Vector<Double> &direction,
606	Float &azimuth,
607	Float &elevation,
608	Float &parangle,
609	Float &opacity,
610	uInt &tcalid,
611	Int &fitid,
612	uInt &focusid,
613	Float &temperature,
614	Float &pressure,
615	Float &humidity,
616	Float &windvel,
617	Float &winddir,
618	Double &srcvel,
619	Vector<Double> &/propermotion/,
620	Vector<Double> &srcdir,
621	Vector<Double> &/scanrate/ )
622	{
623	static const IPosition oneByOne( 1, 1 );
624
625	// DEBUG
626	//cout << "NROReader::getScanInfo() irow = " << irow << endl ;
627	//
628	NRODataRecord *record = dataset_->getRecord( irow ) ;
629
630	// scanno
631	scanno = (uInt)(record->ISCAN) ;
632	//cout << "scanno = " << scanno << endl ;
633
634	// cycleno
635	cycleno = 0 ;
636	//cout << "cycleno = " << cycleno << endl ;
637
638	// beamno and ifno
639	string rxname = dataset_->getRX()[0] ;
640	if ( rxname.find("MULT2") != string::npos ) {
641	string arryt = string( record->ARRYT ) ;
642	beamno = dataset_->getSortedArrayId( arryt ) ;
643	ifno = 0 ;
644	}
645	else {
646	beamno = 0 ;
647	string arryt = string( record->ARRYT ) ;
648	ifno = dataset_->getSortedArrayId( arryt ) ;
649	}
650	//cout << "beamno = " << beamno << endl ;
651
652	// polno
653	polno = dataset_->getPolNo( irow ) ;
654	//cout << "polno = " << polno << endl ;
655
656	// freqs (for IFNO and FREQ_ID)
657	//freqs = getFrequencies( irow ) ;
658	freqs = dataset_->getFrequencies( irow ) ;
659	//cout << "freqs = [" << freqs[0] << ", " << freqs[1] << ", " << freqs[2] << "]" << endl ;
660
661	if ( freqRefFromVREF_ ) {
662	freqs = shiftFrequency( freqs,
663	dataset_->getURVEL(),
664	dataset_->getVDEF() ) ;
665	}
666
667	// restfreq (for MOLECULE_ID)
668	restfreq.resize( oneByOne ) ;
669	restfreq[0] = record->FREQ0 ;
670	//cout << "restfreq = " << restfreq[0] << endl ;
671
672	// refbeamno
673	refbeamno = 0 ;
674	//cout << "refbeamno = " << refbeamno << endl ;
675
676	// scantime
677	//scantime = Double( dataset_->getStartIntTime( irow ) ) ;
678	scantime = Double( dataset_->getScanTime( irow ) ) ;
679	//cout << "scantime = " << scantime << endl ;
680
681	// interval
682	interval = Double( dataset_->getIPTIM() ) ;
683	//cout << "interval = " << interval << endl ;
684
685	// srcname
686	srcname = trim_nro_string( dataset_->getOBJ() ) ;
687	//cout << "srcname = \"" << srcname << "\" (size " << srcname.size() << ")" << endl ;
688
689	// fieldname
690	fieldname = trim_nro_string( dataset_->getOBJ() ) ;
691	//cout << "fieldname = \"" << fieldname << "\" (size " << fieldname.size() << ")" << endl ;
692
693	// spectra
694	vector<double> spec = dataset_->getSpectrum( irow ) ;
695	spectra.resize( spec.size() ) ;
696	//int index = 0 ;
697	Bool b ;
698	Float *fp = spectra.getStorage( b ) ;
699	Float *wp = fp ;
700	for ( vector<double>::iterator i = spec.begin() ;
701	i != spec.end() ; i++ ) {
702	wp = i ;
703	wp++ ;
704	}
705	spectra.putStorage( fp, b ) ;
706	//cout << "spec.size() = " << spec.size() << endl ;
707
708	// flagtra
709	bool setValue = !( flagtra.nelements() == spectra.nelements() ) ;
710	if ( setValue ) {
711	//cout << "flagtra resized. reset values..." << endl ;
712	flagtra.resize( spectra.nelements() ) ;
713	flagtra.set( 0 ) ;
714	}
715	//cout << "flag.size() = " << flagtra.size() << endl ;
716
717	// tsys
718	tsys.resize( oneByOne ) ;
719	tsys[0] = record->TSYS ;
720	//cout << "tsys[0] = " << tsys[0] << endl ;
721
722	// direction
723	direction = getDirection( irow ) ;
724	//cout << "direction = [" << direction[0] << ", " << direction[1] << "]" << endl ;
725
726	// azimuth
727	azimuth = record->RAZ ;
728	//cout << "azimuth = " << azimuth << endl ;
729
730	// elevation
731	elevation = record->REL ;
732	//cout << "elevation = " << elevation << endl ;
733
734	// parangle
735	parangle = 0.0 ;
736	//cout << "parangle = " << parangle << endl ;
737
738	// opacity
739	opacity = 0.0 ;
740	//cout << "opacity = " << opacity << endl ;
741
742	// tcalid
743	tcalid = 0 ;
744	//cout << "tcalid = " << tcalid << endl ;
745
746	// fitid
747	fitid = -1 ;
748	//cout << "fitid = " << fitid << endl ;
749
750	// focusid
751	focusid = 0 ;
752	//cout << "focusid = " << focusid << endl ;
753
754	// temperature (for WEATHER_ID)
755	temperature = Float( record->TEMP ) ;
756	//cout << "temperature = " << temperature << endl ;
757
758	// pressure (for WEATHER_ID)
759	pressure = Float( record->PATM ) ;
760	//cout << "pressure = " << pressure << endl ;
761
762	// humidity (for WEATHER_ID)
763	humidity = Float( record->PH2O ) ;
764	//cout << "humidity = " << humidity << endl ;
765
766	// windvel (for WEATHER_ID)
767	windvel = Float( record->VWIND ) ;
768	//cout << "windvel = " << windvel << endl ;
769
770	// winddir (for WEATHER_ID)
771	winddir = Float( record->DWIND ) ;
772	//cout << "winddir = " << winddir << endl ;
773
774	// srcvel
775	srcvel = dataset_->getURVEL() ;
776	//cout << "srcvel = " << srcvel << endl ;
777
778	// propermotion
779	// do nothing
780
781	// srcdir
782	srcdir = getSourceDirection() ;
783	//cout << "srcdir = [" << srcdir[0] << ", " << srcdir[1] << endl ;
784
785	// scanrate
786	// do nothing
787
788	return 0 ;
789	}
790
791	Int NROReader::getRowNum()
792	{
793	return dataset_->getRowNum() ;
794	}
795
796	vector<double> NROReader::shiftFrequency( const vector<double> &f,
797	const double &v,
798	const string &vdef )
799	{
800	vector<double> r( f ) ;
801	double factor = v / 2.99792458e8 ;
802	if ( vdef.compare( 0, 3, "RAD" ) == 0 ) {
803	factor = 1.0 / ( 1.0 + factor ) ;
804	r[1] *= factor ;
805	r[2] *= factor ;
806	}
807	else if ( vdef.compare( 0, 3, "OPT" ) == 0 ) {
808	factor += 1.0 ;
809	r[1] *= factor ;
810	r[2] *= factor ;
811	}
812	else {
813	cout << "vdef=" << vdef << " is not supported." << endl;
814	}
815	return r ;
816	}

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

Context Navigation

source: trunk/external-alma/atnf/PKSIO/NROReader.cc @ 3106

Download in other formats: