Context Navigation

source: trunk/src/SDMemTable.cc@ 192

Visit:

Last change on this file since 192 was 164, checked in by kil064, 20 years ago
make a few more functions 'const'
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.3 KB

Rev	Line
[2]	1	//#---------------------------------------------------------------------------
	2	//# SDMemTable.cc: A MemoryTable container for single dish integrations
	3	//#---------------------------------------------------------------------------
	4	//# Copyright (C) 2004
[125]	5	//# ATNF
[2]	6	//#
	7	//# This program is free software; you can redistribute it and/or modify it
	8	//# under the terms of the GNU General Public License as published by the Free
	9	//# Software Foundation; either version 2 of the License, or (at your option)
	10	//# any later version.
	11	//#
	12	//# This program is distributed in the hope that it will be useful, but
	13	//# WITHOUT ANY WARRANTY; without even the implied warranty of
	14	//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	15	//# Public License for more details.
	16	//#
	17	//# You should have received a copy of the GNU General Public License along
	18	//# with this program; if not, write to the Free Software Foundation, Inc.,
	19	//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
	20	//#
	21	//# Correspondence concerning this software should be addressed as follows:
	22	//# Internet email: Malte.Marquarding@csiro.au
	23	//# Postal address: Malte Marquarding,
	24	//# Australia Telescope National Facility,
	25	//# P.O. Box 76,
	26	//# Epping, NSW, 2121,
	27	//# AUSTRALIA
	28	//#
	29	//# $Id:
	30	//#---------------------------------------------------------------------------
	31
[125]	32	#include <casa/aips.h>
[80]	33	#include <casa/iostream.h>
	34	#include <casa/iomanip.h>
	35	#include <casa/Arrays/Array.h>
	36	#include <casa/Arrays/ArrayMath.h>
	37	#include <casa/Arrays/MaskArrMath.h>
	38	#include <casa/Arrays/ArrayLogical.h>
	39	#include <casa/Arrays/ArrayAccessor.h>
[2]	40
[80]	41	#include <tables/Tables/TableParse.h>
	42	#include <tables/Tables/TableDesc.h>
	43	#include <tables/Tables/SetupNewTab.h>
	44	#include <tables/Tables/ScaColDesc.h>
	45	#include <tables/Tables/ArrColDesc.h>
[2]	46
[80]	47	#include <tables/Tables/ExprNode.h>
	48	#include <tables/Tables/ScalarColumn.h>
	49	#include <tables/Tables/ArrayColumn.h>
	50	#include <tables/Tables/TableRecord.h>
	51	#include <measures/Measures/MFrequency.h>
	52	#include <measures/Measures/MeasTable.h>
[105]	53	#include <coordinates/Coordinates/CoordinateUtil.h>
[80]	54	#include <casa/Quanta/MVTime.h>
[2]	55
	56	#include "SDMemTable.h"
	57	#include "SDContainer.h"
	58
[125]	59	using namespace casa;
[83]	60	using namespace asap;
[2]	61
[18]	62	SDMemTable::SDMemTable() :
	63	IFSel_(0),
	64	beamSel_(0),
	65	polSel_(0) {
	66	setup();
	67	}
[2]	68	SDMemTable::SDMemTable(const std::string& name) :
	69	IFSel_(0),
	70	beamSel_(0),
	71	polSel_(0) {
[50]	72	Table tab(name);
[22]	73	table_ = tab.copyToMemoryTable("dummy");
[2]	74	}
	75
[16]	76	SDMemTable::SDMemTable(const SDMemTable& other, Bool clear) {
[148]	77	IFSel_= other.IFSel_;
	78	beamSel_= other.beamSel_;
	79	polSel_= other.polSel_;
	80	chanMask_ = other.chanMask_;
	81	table_ = other.table_.copyToMemoryTable(String("dummy"));
[2]	82	// clear all rows()
[16]	83	if (clear) {
[148]	84	table_.removeRow(this->table_.rowNumbers());
[16]	85	} else {
[148]	86	IFSel_ = other.IFSel_;
	87	beamSel_ = other.beamSel_;
	88	polSel_ = other.polSel_;
[16]	89	}
[2]	90	}
	91
[80]	92	SDMemTable::SDMemTable(const Table& tab, const std::string& exprs) :
[2]	93	IFSel_(0),
	94	beamSel_(0),
	95	polSel_(0) {
	96	Table t = tableCommand(exprs,tab);
[89]	97	if (t.nrow() == 0)
	98	throw(AipsError("Query unsuccessful."));
[22]	99	table_ = t.copyToMemoryTable("dummy");
[2]	100	}
	101
	102	SDMemTable::~SDMemTable(){
[105]	103	//cerr << "goodbye from SDMemTable @ " << this << endl;
[2]	104	}
	105
[164]	106	SDMemTable SDMemTable::getScan(Int scanID) const {
[80]	107	String cond("SELECT * from $1 WHERE SCANID == ");
	108	cond += String::toString(scanID);
	109	return SDMemTable(table_, cond);
[2]	110	}
	111
[138]	112	SDMemTable &SDMemTable::operator=(const SDMemTable& other) {
[148]	113	if (this != &other) {
	114	IFSel_= other.IFSel_;
	115	beamSel_= other.beamSel_;
	116	polSel_= other.polSel_;
	117	chanMask_.resize(0);
	118	chanMask_ = other.chanMask_;
	119	table_ = other.table_.copyToMemoryTable(String("dummy"));
	120	}
[138]	121	return *this;
	122	}
	123
[164]	124	SDMemTable SDMemTable::getSource(const std::string& source) const {
[80]	125	String cond("SELECT * from $1 WHERE SRCNAME == ");
	126	cond += source;
	127	return SDMemTable(table_, cond);
	128	}
	129
[2]	130	void SDMemTable::setup() {
	131	TableDesc td("", "1", TableDesc::Scratch);
	132	td.comment() = "A SDMemTable";
	133	td.addColumn(ScalarColumnDesc<Double>("TIME"));
	134	td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
	135	td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
	136	td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
[89]	137	td.addColumn(ArrayColumnDesc<Float>("TSYS"));
	138	td.addColumn(ScalarColumnDesc<Int>("SCANID"));
	139	td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
[39]	140	td.addColumn(ArrayColumnDesc<uInt>("FREQID"));
[78]	141	td.addColumn(ArrayColumnDesc<Double>("DIRECTION"));
[105]	142	td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
	143	td.addColumn(ScalarColumnDesc<String>("TCALTIME"));
	144	td.addColumn(ArrayColumnDesc<Float>("TCAL"));
	145	td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
	146	td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
	147	td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
	148	td.addColumn(ScalarColumnDesc<Int>("REFBEAM"));
	149
[2]	150	// Now create a new table from the description.
[18]	151
[22]	152	SetupNewTable aNewTab("dummy", td, Table::New);
[89]	153	table_ = Table(aNewTab, Table::Memory, 0);
[2]	154	}
	155
	156	std::string SDMemTable::getSourceName(Int whichRow) const {
	157	ROScalarColumn<String> src(table_, "SRCNAME");
	158	String name;
	159	src.get(whichRow, name);
	160	return name;
	161	}
	162
[50]	163	std::string SDMemTable::getTime(Int whichRow) const {
[2]	164	ROScalarColumn<Double> src(table_, "TIME");
	165	Double tm;
	166	src.get(whichRow, tm);
[50]	167	MVTime mvt(tm);
[105]	168	mvt.setFormat(MVTime::TIME);
[50]	169	ostringstream oss;
	170	oss << mvt;
	171	String str(oss);
	172	return str;
[2]	173	}
[50]	174	double SDMemTable::getInterval(Int whichRow) const {
	175	ROScalarColumn<Double> src(table_, "INTERVAL");
	176	Double intval;
	177	src.get(whichRow, intval);
	178	return intval;
	179	}
[2]	180
	181	bool SDMemTable::setIF(Int whichIF) {
[50]	182	if ( whichIF >= 0 && whichIF < nIF()) {
[2]	183	IFSel_ = whichIF;
	184	return true;
[50]	185	}
	186	return false;
[2]	187	}
[50]	188
[2]	189	bool SDMemTable::setBeam(Int whichBeam) {
[50]	190	if ( whichBeam >= 0 && whichBeam < nBeam()) {
[2]	191	beamSel_ = whichBeam;
	192	return true;
[50]	193	}
	194	return false;
	195	}
[2]	196
	197	bool SDMemTable::setPol(Int whichPol) {
[50]	198	if ( whichPol >= 0 && whichPol < nPol()) {
[2]	199	polSel_ = whichPol;
	200	return true;
[50]	201	}
	202	return false;
[2]	203	}
	204
[68]	205	bool SDMemTable::setMask(std::vector<int> whichChans) {
[16]	206	ROArrayColumn<uChar> spec(table_, "FLAGTRA");
	207	std::vector<int>::iterator it;
	208	uInt n = spec.shape(0)(3);
[105]	209	if (whichChans.empty()) {
	210	chanMask_ = std::vector<bool>(n,true);
	211	return true;
	212	}
[16]	213	chanMask_.resize(n,true);
[39]	214	for (it = whichChans.begin(); it != whichChans.end(); ++it) {
[105]	215	if (*it < n) {
[16]	216	chanMask_[*it] = false;
[105]	217	}
[16]	218	}
[2]	219	return true;
	220	}
	221
[16]	222	std::vector<bool> SDMemTable::getMask(Int whichRow) const {
	223	std::vector<bool> mask;
	224	ROArrayColumn<uChar> spec(table_, "FLAGTRA");
	225	Array<uChar> arr;
	226	spec.get(whichRow, arr);
	227	ArrayAccessor<uChar, Axis<0> > aa0(arr);
	228	aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
	229	ArrayAccessor<uChar, Axis<1> > aa1(aa0);
	230	aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
	231	ArrayAccessor<uChar, Axis<2> > aa2(aa1);
	232	aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
	233
	234	Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
	235
	236	std::vector<bool> tmp;
	237	tmp = chanMask_; // WHY the fxxx do I have to make a copy here
	238	std::vector<bool>::iterator miter;
	239	miter = tmp.begin();
	240
	241	for (ArrayAccessor<uChar, Axis<3> > i(aa2); i != i.end(); ++i) {
	242	bool out =!static_cast<bool>(*i);
	243	if (useUserMask) {
	244	out = out && (*miter);
	245	miter++;
	246	}
	247	mask.push_back(out);
[89]	248	}
[16]	249	return mask;
[2]	250	}
[50]	251	std::vector<float> SDMemTable::getSpectrum(Int whichRow) const {
[2]	252
	253	std::vector<float> spectrum;
	254	ROArrayColumn<Float> spec(table_, "SPECTRA");
	255	Array<Float> arr;
	256	spec.get(whichRow, arr);
	257	ArrayAccessor<Float, Axis<0> > aa0(arr);
	258	aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
	259	ArrayAccessor<Float, Axis<1> > aa1(aa0);
	260	aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
[16]	261	ArrayAccessor<Float, Axis<2> > aa2(aa1);
	262	aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
	263	for (ArrayAccessor<Float, Axis<3> > i(aa2); i != i.end(); ++i) {
	264	spectrum.push_back(*i);
[2]	265	}
	266	return spectrum;
	267	}
[105]	268	std::vector<string> SDMemTable::getCoordInfo() const {
	269	String un;
	270	Table t = table_.keywordSet().asTable("FREQUENCIES");
	271	String sunit;
	272	t.keywordSet().get("UNIT",sunit);
	273	String dpl;
	274	t.keywordSet().get("DOPPLER",dpl);
	275	if (dpl == "") dpl = "RADIO";
	276	String rfrm;
	277	t.keywordSet().get("REFFRAME",rfrm);
	278	std::vector<string> inf;
	279	inf.push_back(sunit);
	280	inf.push_back(rfrm);
	281	inf.push_back(dpl);
	282	return inf;
	283	}
[39]	284
[105]	285	void SDMemTable::setCoordInfo(std::vector<string> theinfo) {
	286
	287	std::vector<string>::iterator it;
	288	String un,rfrm,dpl;
	289	un = theinfo[0];
	290	rfrm = theinfo[1];
	291	dpl = theinfo[2];
	292
	293	//String un(theunit);
	294	Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
	295	Vector<Double> rstf;
	296	t.keywordSet().get("RESTFREQS",rstf);
	297	Bool canDo = True;
	298	Unit u1("km/s");Unit u2("Hz");
	299	if (Unit(un) == u1) {
	300	Vector<Double> rstf;
	301	t.keywordSet().get("RESTFREQS",rstf);
	302	if (rstf.nelements() == 0) {
	303	throw(AipsError("Can't set unit to km/s if no restfrequencies are specified"));
	304	}
	305	} else if (Unit(un) != u2 && un != "") {
	306	throw(AipsError("Unit not conformant with Spectral Coordinates"));
	307	}
	308	t.rwKeywordSet().define("UNIT", un);
	309
	310	MFrequency::Types mdr;
	311	if (!MFrequency::getType(mdr, rfrm)) {
	312
	313	Int a,b;const uInt* c;
	314	const String* valid = MFrequency::allMyTypes(a, b, c);
	315	String pfix = "Please specify a legal frame type. Types are\n";
	316	throw(AipsError(pfix+(*valid)));
	317	} else {
	318	t.rwKeywordSet().define("REFFRAME",rfrm);
	319	}
	320
	321	}
	322
[164]	323	std::vector<double> SDMemTable::getAbcissa(Int whichRow) const {
[39]	324	std::vector<double> absc(nChan());
[89]	325	Vector<Double> absc1(nChan());
[39]	326	indgen(absc1);
	327	ROArrayColumn<uInt> fid(table_, "FREQID");
	328	Vector<uInt> v;
	329	fid.get(whichRow, v);
	330	uInt specidx = v(IFSel_);
	331	SpectralCoordinate spc = getCoordinate(specidx);
[78]	332	Table t = table_.keywordSet().asTable("FREQUENCIES");
	333	String rf;
[105]	334	//t.keywordSet().get("EQUINOX",rf);
[78]	335	MDirection::Types mdr;
[105]	336	//if (!MDirection::getType(mdr, rf)) {
	337	mdr = MDirection::J2000;
	338	//cout << "Unknown equinox using J2000" << endl;
	339	//}
[78]	340	ROArrayColumn<Double> dir(table_, "DIRECTION");
	341	Array<Double> posit;
	342	dir.get(whichRow,posit);
	343	Vector<Double> wpos(2);
	344	wpos[0] = posit(IPosition(2,beamSel_,0));
	345	wpos[1] = posit(IPosition(2,beamSel_,1));
	346	Quantum<Double> lon(wpos[0],Unit(String("rad")));
	347	Quantum<Double> lat(wpos[1],Unit(String("rad")));
	348	MDirection direct(lon, lat, mdr);
	349	ROScalarColumn<Double> tme(table_, "TIME");
	350	Double obstime;
	351	tme.get(whichRow,obstime);
[89]	352	MVEpoch tm2(Quantum<Double>(obstime, Unit(String("d"))));
[78]	353	MEpoch epoch(tm2);
	354
	355	Vector<Double> antpos;
	356	table_.keywordSet().get("AntennaPosition", antpos);
	357	MVPosition mvpos(antpos(0),antpos(1),antpos(2));
	358	MPosition pos(mvpos);
[105]	359	String sunit;
	360	t.keywordSet().get("UNIT",sunit);
	361	if (sunit == "") sunit = "pixel";
	362	Unit u(sunit);
	363	String frm;
	364	t.keywordSet().get("REFFRAME",frm);
	365	if (frm == "") frm = "TOPO";
	366	String dpl;
	367	t.keywordSet().get("DOPPLER",dpl);
	368	if (dpl == "") dpl = "RADIO";
[78]	369	MFrequency::Types mtype;
[105]	370	if (!MFrequency::getType(mtype, frm)) {
[89]	371	cout << "Frequency type unknown assuming TOPO" << endl;
[78]	372	mtype = MFrequency::TOPO;
	373	}
[105]	374
	375	if (!spc.setReferenceConversion(mtype,epoch,pos,direct)) {
	376	throw(AipsError("Couldn't convert frequency frame."));
	377	}
[78]	378
[39]	379	if ( u == Unit("km/s") ) {
[89]	380	Vector<Double> rstf;
	381	t.keywordSet().get("RESTFREQS",rstf);
	382	if (rstf.nelements() > 0) {
[121]	383	if (rstf.nelements() >= nIF())
[125]	384	spc.selectRestFrequency(uInt(IFSel_));
[89]	385	spc.setVelocity(u.getName());
[39]	386	Vector<Double> wrld;
	387	spc.pixelToVelocity(wrld,absc1);
	388	std::vector<double>::iterator it;
	389	uInt i = 0;
	390	for (it = absc.begin(); it != absc.end(); ++it) {
[89]	391	(*it) = wrld[i];
	392	i++;
[39]	393	}
	394	}
	395	} else if (u == Unit("Hz")) {
	396	Vector<String> wau(1); wau = u.getName();
	397	spc.setWorldAxisUnits(wau);
	398	std::vector<double>::iterator it;
	399	Double tmp;
	400	uInt i = 0;
[89]	401	for (it = absc.begin(); it != absc.end(); ++it) {
[39]	402	spc.toWorld(tmp,absc1[i]);
	403	(*it) = tmp;
	404	i++;
	405	}
[105]	406
	407	} else {
	408	// assume channels/pixels
	409	std::vector<double>::iterator it;
	410	uInt i=0;
	411	for (it = absc.begin(); it != absc.end(); ++it) {
	412	(*it) = Double(i++);
	413	}
[39]	414	}
	415	return absc;
	416	}
	417
[164]	418	std::string SDMemTable::getAbcissaString(Int whichRow) const
[105]	419	{
	420	ROArrayColumn<uInt> fid(table_, "FREQID");
	421	Table t = table_.keywordSet().asTable("FREQUENCIES");
	422	String sunit;
	423	t.keywordSet().get("UNIT",sunit);
	424	if (sunit == "") sunit = "pixel";
	425	Unit u(sunit);
	426	Vector<uInt> v;
	427	fid.get(whichRow, v);
	428	uInt specidx = v(IFSel_);
	429	SpectralCoordinate spc = getCoordinate(specidx);
	430	String frm;
	431	t.keywordSet().get("REFFRAME",frm);
	432	MFrequency::Types mtype;
	433	if (!MFrequency::getType(mtype, frm)) {
	434	cout << "Frequency type unknown assuming TOPO" << endl;
	435	mtype = MFrequency::TOPO;
	436	}
	437	spc.setFrequencySystem(mtype);
	438	String s = "Channel";
	439	if (u == Unit("km/s")) {
	440	spc.setVelocity(u.getName());
	441	s = CoordinateUtil::axisLabel(spc,0,True,True,True);
	442	} else if (u == Unit("Hz")) {
	443	Vector<String> wau(1);wau = u.getName();
	444	spc.setWorldAxisUnits(wau);
	445	s = CoordinateUtil::axisLabel(spc);
	446	}
	447	return s;
	448	}
	449
[89]	450	void SDMemTable::setSpectrum(std::vector<float> spectrum, int whichRow) {
	451	ArrayColumn<Float> spec(table_, "SPECTRA");
	452	Array<Float> arr;
	453	spec.get(whichRow, arr);
	454	if (spectrum.size() != arr.shape()(3)) {
	455	throw(AipsError("Attempting to set spectrum with incorrect length."));
	456	}
	457
	458	ArrayAccessor<Float, Axis<0> > aa0(arr);
	459	aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
	460	ArrayAccessor<Float, Axis<1> > aa1(aa0);
	461	aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
	462	ArrayAccessor<Float, Axis<2> > aa2(aa1);
	463	aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
	464
	465	std::vector<float>::iterator it = spectrum.begin();
	466	for (ArrayAccessor<Float, Axis<3> > i(aa2); i != i.end(); ++i) {
	467	(i) = Float(it);
	468	it++;
	469	}
	470	spec.put(whichRow, arr);
	471	}
	472
[164]	473	void SDMemTable::getSpectrum(Vector<Float>& spectrum, Int whichRow) const {
[21]	474	ROArrayColumn<Float> spec(table_, "SPECTRA");
	475	Array<Float> arr;
	476	spec.get(whichRow, arr);
	477	spectrum.resize(arr.shape()(3));
	478	ArrayAccessor<Float, Axis<0> > aa0(arr);
	479	aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
	480	ArrayAccessor<Float, Axis<1> > aa1(aa0);
	481	aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
	482	ArrayAccessor<Float, Axis<2> > aa2(aa1);
	483	aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
[2]	484
[21]	485	ArrayAccessor<Float, Axis<0> > va(spectrum);
	486	for (ArrayAccessor<Float, Axis<3> > i(aa2); i != i.end(); ++i) {
	487	(va) = (i);
	488	va++;
	489	}
	490	}
[89]	491	/*
[68]	492	void SDMemTable::getMask(Vector<Bool>& mask, Int whichRow) const {
[21]	493	ROArrayColumn<uChar> spec(table_, "FLAGTRA");
	494	Array<uChar> arr;
	495	spec.get(whichRow, arr);
	496	mask.resize(arr.shape()(3));
	497
	498	ArrayAccessor<uChar, Axis<0> > aa0(arr);
	499	aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
	500	ArrayAccessor<uChar, Axis<1> > aa1(aa0);
	501	aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
	502	ArrayAccessor<uChar, Axis<2> > aa2(aa1);
	503	aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
	504
	505	Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
[89]	506
[21]	507	ArrayAccessor<Bool, Axis<0> > va(mask);
	508	std::vector<bool> tmp;
	509	tmp = chanMask_; // WHY the fxxx do I have to make a copy here. The
[89]	510	// iterator should work on chanMask_??
[21]	511	std::vector<bool>::iterator miter;
	512	miter = tmp.begin();
	513
	514	for (ArrayAccessor<uChar, Axis<3> > i(aa2); i != i.end(); ++i) {
	515	bool out =!static_cast<bool>(*i);
	516	if (useUserMask) {
	517	out = out && (*miter);
	518	miter++;
	519	}
	520	(*va) = out;
	521	va++;
[89]	522	}
[21]	523	}
[89]	524	*/
[16]	525	MaskedArray<Float> SDMemTable::rowAsMaskedArray(uInt whichRow,
[164]	526	Bool useSelection) const
	527	{
[2]	528	ROArrayColumn<Float> spec(table_, "SPECTRA");
	529	Array<Float> arr;
	530	ROArrayColumn<uChar> flag(table_, "FLAGTRA");
	531	Array<uChar> farr;
	532	spec.get(whichRow, arr);
	533	flag.get(whichRow, farr);
	534	Array<Bool> barr(farr.shape());convertArray(barr, farr);
	535	MaskedArray<Float> marr;
[16]	536	if (useSelection) {
	537	ArrayAccessor<Float, Axis<0> > aa0(arr);
	538	aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
	539	ArrayAccessor<Float, Axis<1> > aa1(aa0);
	540	aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
	541	ArrayAccessor<Float, Axis<2> > aa2(aa1);
	542	aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
[89]	543
[16]	544	ArrayAccessor<Bool, Axis<0> > baa0(barr);
	545	baa0.reset(baa0.begin(uInt(beamSel_)));//go to beam
	546	ArrayAccessor<Bool, Axis<1> > baa1(baa0);
	547	baa1.reset(baa1.begin(uInt(IFSel_)));// go to IF
	548	ArrayAccessor<Bool, Axis<2> > baa2(baa1);
	549	baa2.reset(baa2.begin(uInt(polSel_)));// go to pol
	550
	551	Vector<Float> a(arr.shape()(3));
	552	Vector<Bool> b(barr.shape()(3));
	553	ArrayAccessor<Float, Axis<0> > a0(a);
	554	ArrayAccessor<Bool, Axis<0> > b0(b);
	555
	556	ArrayAccessor<Bool, Axis<3> > j(baa2);
	557	for (ArrayAccessor<Float, Axis<3> > i(aa2); i != i.end(); ++i) {
	558	(a0) = (i);
	559	(b0) = !(j);
	560	j++;
	561	a0++;
	562	b0++;
	563	}
	564	marr.setData(a,b);
	565	} else {
	566	marr.setData(arr,!barr);
	567	}
[2]	568	return marr;
	569	}
	570
	571	Float SDMemTable::getTsys(Int whichRow) const {
	572	ROArrayColumn<Float> ts(table_, "TSYS");
	573	Array<Float> arr;
	574	ts.get(whichRow, arr);
	575	Float out;
	576	IPosition ip(arr.shape());
[16]	577	ip(0) = beamSel_;ip(1) = IFSel_;ip(2) = polSel_;ip(3)=0;
[2]	578	out = arr(ip);
	579	return out;
	580	}
	581
[39]	582	SpectralCoordinate SDMemTable::getCoordinate(uInt whichIdx) const {
	583
	584	Table t = table_.keywordSet().asTable("FREQUENCIES");
	585	if (whichIdx > t.nrow() ) {
	586	cerr << "SDMemTable::getCoordinate - whichIdx out of range" << endl;
[68]	587	return SpectralCoordinate();
[39]	588	}
[89]	589
[39]	590	Double rp,rv,inc;
	591	String rf;
[89]	592	Vector<Double> vec;
[39]	593	ROScalarColumn<Double> rpc(t, "REFPIX");
	594	ROScalarColumn<Double> rvc(t, "REFVAL");
	595	ROScalarColumn<Double> incc(t, "INCREMENT");
[89]	596	t.keywordSet().get("RESTFREQS",vec);
[105]	597	t.keywordSet().get("BASEREFFRAME",rf);
[89]	598
[39]	599	MFrequency::Types mft;
	600	if (!MFrequency::getType(mft, rf)) {
	601	cerr << "Frequency type unknown assuming TOPO" << endl;
[89]	602	mft = MFrequency::TOPO;
	603	}
[39]	604	rpc.get(whichIdx, rp);
	605	rvc.get(whichIdx, rv);
	606	incc.get(whichIdx, inc);
	607	SpectralCoordinate spec(mft,rv,inc,rp);
[89]	608	if (vec.nelements() > 0)
	609	spec.setRestFrequencies(vec);
[39]	610	return spec;
	611	}
	612
[89]	613	Bool SDMemTable::setCoordinate(const SpectralCoordinate& speccord,
	614	uInt whichIdx) {
[50]	615	Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
	616	if (whichIdx > t.nrow() ) {
[89]	617	throw(AipsError("SDMemTable::setCoordinate - coord no out of range"));
[50]	618	}
	619	ScalarColumn<Double> rpc(t, "REFPIX");
	620	ScalarColumn<Double> rvc(t, "REFVAL");
	621	ScalarColumn<Double> incc(t, "INCREMENT");
[89]	622
[50]	623	rpc.put(whichIdx, speccord.referencePixel()[0]);
	624	rvc.put(whichIdx, speccord.referenceValue()[0]);
	625	incc.put(whichIdx, speccord.increment()[0]);
	626
	627	return True;
	628	}
	629
[89]	630	Int SDMemTable::nCoordinates() const
	631	{
	632	return table_.keywordSet().asTable("FREQUENCIES").nrow();
	633	}
[50]	634
[89]	635	void SDMemTable::setRestFreqs(std::vector<double> freqs, const std::string& theunit)
	636	{
	637	Vector<Double> tvec(freqs);
	638	Quantum<Vector<Double> > q(tvec, String(theunit));
	639	tvec.resize();
	640	tvec = q.getValue("Hz");
	641	Table t = table_.keywordSet().asTable("FREQUENCIES");
	642	t.rwKeywordSet().define("RESTFREQS",tvec);
	643	}
	644
[39]	645	bool SDMemTable::putSDFreqTable(const SDFrequencyTable& sdft) {
	646	TableDesc td("", "1", TableDesc::Scratch);
	647	td.addColumn(ScalarColumnDesc<Double>("REFPIX"));
	648	td.addColumn(ScalarColumnDesc<Double>("REFVAL"));
	649	td.addColumn(ScalarColumnDesc<Double>("INCREMENT"));
	650	SetupNewTable aNewTab("freqs", td, Table::New);
	651	Table aTable (aNewTab, Table::Memory, sdft.length());
	652	ScalarColumn<Double> sc0(aTable, "REFPIX");
	653	ScalarColumn<Double> sc1(aTable, "REFVAL");
	654	ScalarColumn<Double> sc2(aTable, "INCREMENT");
	655	for (uInt i=0; i < sdft.length(); ++i) {
	656	sc0.put(i,sdft.referencePixel(i));
	657	sc1.put(i,sdft.referenceValue(i));
	658	sc2.put(i,sdft.increment(i));
	659	}
[105]	660	String rf = sdft.refFrame();
	661	if (rf.contains("TOPO")) rf = "TOPO";
	662
	663	aTable.rwKeywordSet().define("BASEREFFRAME", rf);
	664	aTable.rwKeywordSet().define("REFFRAME", rf);
[39]	665	aTable.rwKeywordSet().define("EQUINOX", sdft.equinox());
[105]	666	aTable.rwKeywordSet().define("UNIT", String(""));
	667	aTable.rwKeywordSet().define("DOPPLER", String("RADIO"));
[89]	668	Vector<Double> rfvec;
	669	aTable.rwKeywordSet().define("RESTFREQS", rfvec);
[39]	670	table_.rwKeywordSet().defineTable ("FREQUENCIES", aTable);
	671	return True;
	672	}
	673
	674	SDFrequencyTable SDMemTable::getSDFreqTable() const {
	675	SDFrequencyTable sdft;
[89]	676
[39]	677	return sdft;
	678	}
	679
[2]	680	bool SDMemTable::putSDContainer(const SDContainer& sdc) {
	681	ScalarColumn<Double> mjd(table_, "TIME");
	682	ScalarColumn<String> srcn(table_, "SRCNAME");
[105]	683	ScalarColumn<String> fldn(table_, "FIELDNAME");
[2]	684	ArrayColumn<Float> spec(table_, "SPECTRA");
	685	ArrayColumn<uChar> flags(table_, "FLAGTRA");
	686	ArrayColumn<Float> ts(table_, "TSYS");
	687	ScalarColumn<Int> scan(table_, "SCANID");
[16]	688	ScalarColumn<Double> integr(table_, "INTERVAL");
[39]	689	ArrayColumn<uInt> freqid(table_, "FREQID");
[78]	690	ArrayColumn<Double> dir(table_, "DIRECTION");
[105]	691	ScalarColumn<Int> rbeam(table_, "REFBEAM");
	692	ArrayColumn<Float> tcal(table_, "TCAL");
	693	ScalarColumn<String> tcalt(table_, "TCALTIME");
	694	ScalarColumn<Float> az(table_, "AZIMUTH");
	695	ScalarColumn<Float> el(table_, "ELEVATION");
	696	ScalarColumn<Float> para(table_, "PARANGLE");
[2]	697
	698	uInt rno = table_.nrow();
	699	table_.addRow();
[89]	700
[2]	701	mjd.put(rno, sdc.timestamp);
	702	srcn.put(rno, sdc.sourcename);
[105]	703	fldn.put(rno, sdc.fieldname);
[2]	704	spec.put(rno, sdc.getSpectrum());
	705	flags.put(rno, sdc.getFlags());
	706	ts.put(rno, sdc.getTsys());
	707	scan.put(rno, sdc.scanid);
[16]	708	integr.put(rno, sdc.interval);
[39]	709	freqid.put(rno, sdc.getFreqMap());
[78]	710	dir.put(rno, sdc.getDirection());
[105]	711	rbeam.put(rno, sdc.refbeam);
	712	tcal.put(rno, sdc.tcal);
	713	tcalt.put(rno, sdc.tcaltime);
	714	az.put(rno, sdc.azimuth);
	715	el.put(rno, sdc.elevation);
	716	para.put(rno, sdc.parangle);
[89]	717
[2]	718	return true;
	719	}
[18]	720
[21]	721	SDContainer SDMemTable::getSDContainer(uInt whichRow) const {
	722	ROScalarColumn<Double> mjd(table_, "TIME");
	723	ROScalarColumn<String> srcn(table_, "SRCNAME");
[105]	724	ROScalarColumn<String> fldn(table_, "FIELDNAME");
[21]	725	ROArrayColumn<Float> spec(table_, "SPECTRA");
	726	ROArrayColumn<uChar> flags(table_, "FLAGTRA");
	727	ROArrayColumn<Float> ts(table_, "TSYS");
	728	ROScalarColumn<Int> scan(table_, "SCANID");
	729	ROScalarColumn<Double> integr(table_, "INTERVAL");
[39]	730	ROArrayColumn<uInt> freqid(table_, "FREQID");
[78]	731	ROArrayColumn<Double> dir(table_, "DIRECTION");
[105]	732	ROScalarColumn<Int> rbeam(table_, "REFBEAM");
	733	ROArrayColumn<Float> tcal(table_, "TCAL");
	734	ROScalarColumn<String> tcalt(table_, "TCALTIME");
	735	ROScalarColumn<Float> az(table_, "AZIMUTH");
	736	ROScalarColumn<Float> el(table_, "ELEVATION");
	737	ROScalarColumn<Float> para(table_, "PARANGLE");
[21]	738
	739	SDContainer sdc(nBeam(),nIF(),nPol(),nChan());
	740	mjd.get(whichRow, sdc.timestamp);
	741	srcn.get(whichRow, sdc.sourcename);
	742	integr.get(whichRow, sdc.interval);
	743	scan.get(whichRow, sdc.scanid);
[105]	744	fldn.get(whichRow, sdc.fieldname);
	745	rbeam.get(whichRow, sdc.refbeam);
	746	az.get(whichRow, sdc.azimuth);
	747	el.get(whichRow, sdc.elevation);
	748	para.get(whichRow, sdc.parangle);
	749	Vector<Float> tc;
	750	tcal.get(whichRow, tc);
	751	sdc.tcal[0] = tc[0];sdc.tcal[1] = tc[1];
	752	tcalt.get(whichRow, sdc.tcaltime);
[21]	753	Array<Float> spectrum;
	754	Array<Float> tsys;
	755	Array<uChar> flagtrum;
[39]	756	Vector<uInt> fmap;
[78]	757	Array<Double> direction;
[21]	758	spec.get(whichRow, spectrum);
	759	sdc.putSpectrum(spectrum);
	760	flags.get(whichRow, flagtrum);
	761	sdc.putFlags(flagtrum);
	762	ts.get(whichRow, tsys);
	763	sdc.putTsys(tsys);
[39]	764	freqid.get(whichRow, fmap);
	765	sdc.putFreqMap(fmap);
[78]	766	dir.get(whichRow, direction);
	767	sdc.putDirection(direction);
[21]	768	return sdc;
	769	}
[78]	770
[18]	771	bool SDMemTable::putSDHeader(const SDHeader& sdh) {
	772	table_.lock();
	773	table_.rwKeywordSet().define("nIF", sdh.nif);
	774	table_.rwKeywordSet().define("nBeam", sdh.nbeam);
	775	table_.rwKeywordSet().define("nPol", sdh.npol);
	776	table_.rwKeywordSet().define("nChan", sdh.nchan);
	777	table_.rwKeywordSet().define("Observer", sdh.observer);
	778	table_.rwKeywordSet().define("Project", sdh.project);
	779	table_.rwKeywordSet().define("Obstype", sdh.obstype);
	780	table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
	781	table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
	782	table_.rwKeywordSet().define("Equinox", sdh.equinox);
	783	table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
	784	table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
	785	table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
	786	table_.rwKeywordSet().define("UTC", sdh.utc);
	787	table_.unlock();
	788	return true;
[50]	789	}
[21]	790
	791	SDHeader SDMemTable::getSDHeader() const {
	792	SDHeader sdh;
	793	table_.keywordSet().get("nBeam",sdh.nbeam);
	794	table_.keywordSet().get("nIF",sdh.nif);
	795	table_.keywordSet().get("nPol",sdh.npol);
	796	table_.keywordSet().get("nChan",sdh.nchan);
	797	table_.keywordSet().get("Observer", sdh.observer);
	798	table_.keywordSet().get("Project", sdh.project);
	799	table_.keywordSet().get("Obstype", sdh.obstype);
	800	table_.keywordSet().get("AntennaName", sdh.antennaname);
	801	table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
	802	table_.keywordSet().get("Equinox", sdh.equinox);
	803	table_.keywordSet().get("FreqRefFrame", sdh.freqref);
	804	table_.keywordSet().get("FreqRefVal", sdh.reffreq);
	805	table_.keywordSet().get("Bandwidth", sdh.bandwidth);
	806	table_.keywordSet().get("UTC", sdh.utc);
	807	return sdh;
[18]	808	}
[2]	809	void SDMemTable::makePersistent(const std::string& filename) {
	810	table_.deepCopy(filename,Table::New);
	811	}
	812
[89]	813	Int SDMemTable::nScan() const {
[50]	814	Int n = 0;
	815	ROScalarColumn<Int> scans(table_, "SCANID");
	816	Int previous = -1;Int current=0;
	817	for (uInt i=0; i< scans.nrow();i++) {
	818	scans.getScalar(i,current);
	819	if (previous != current) {
[89]	820	previous = current;
[50]	821	n++;
	822	}
	823	}
	824	return n;
	825	}
	826
[105]	827	String SDMemTable::formatSec(Double x) {
	828	Double xcop = x;
	829	MVTime mvt(xcop/24./3600.); // make days
	830	if (x < 59.95)
	831	return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
	832	return mvt.string(MVTime::TIME_CLEAN_NO_H, 7)+" ";
	833	};
	834
[164]	835	std::string SDMemTable::summary() {
[2]	836	ROScalarColumn<Int> scans(table_, "SCANID");
	837	ROScalarColumn<String> srcs(table_, "SRCNAME");
[89]	838	ostringstream oss;
	839	oss << endl;
	840	oss << "--------------------------------------------------" << endl;
	841	oss << " Scan Table Summary" << endl;
	842	oss << "--------------------------------------------------" << endl;
	843	oss.flags(std::ios_base::left);
	844	oss << setw(15) << "Beams:" << setw(4) << nBeam() << endl
	845	<< setw(15) << "IFs:" << setw(4) << nIF() << endl
	846	<< setw(15) << "Polarisations:" << setw(4) << nPol() << endl
	847	<< setw(15) << "Channels:" << setw(4) << nChan() << endl;
	848	oss << endl;
	849	String tmp;
	850	table_.keywordSet().get("Observer", tmp);
	851	oss << setw(15) << "Observer:" << tmp << endl;
	852	table_.keywordSet().get("Project", tmp);
	853	oss << setw(15) << "Project:" << tmp << endl;
	854	table_.keywordSet().get("Obstype", tmp);
	855	oss << setw(15) << "Obs. Type:" << tmp << endl;
	856	table_.keywordSet().get("AntennaName", tmp);
	857	oss << setw(15) << "Antenna Name:" << tmp << endl;
[105]	858	Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
	859	Vector<Double> vec;
	860	t.keywordSet().get("RESTFREQS",vec);
	861	oss << setw(15) << "Rest Freqs:";
	862	if (vec.nelements() > 0) {
	863	oss << setprecision(0) << vec << " [Hz]" << endl;
	864	} else {
	865	oss << "None set" << endl;
	866	}
[158]	867	oss << setw(15) << "Abcissa:" << getAbcissaString() << endl;
[105]	868	oss << setw(15) << "Cursor:" << "Beam[" << getBeam() << "] "
	869	<< "IF[" << getIF() << "] " << "Pol[" << getPol() << "]" << endl;
[89]	870	oss << endl;
[2]	871	uInt count = 0;
	872	String name;
	873	Int previous = -1;Int current=0;
[89]	874	Int integ = 0;
	875	oss << setw(6) << "Scan"
	876	<< setw(12) << "Source"
	877	<< setw(21) << "Time"
	878	<< setw(11) << "Integration" << endl;
[105]	879	oss << "--------------------------------------------------" << endl;
[2]	880	for (uInt i=0; i< scans.nrow();i++) {
	881	scans.getScalar(i,current);
	882	if (previous != current) {
	883	srcs.getScalar(i,name);
[50]	884	previous = current;
[105]	885	String t = formatSec(Double(getInterval(i)));
[89]	886	oss << setw(6) << count << setw(12) << name << setw(21) << getTime(i)
[105]	887	<< setw(2) << setprecision(1)
	888	<< t << endl;
[50]	889	count++;
[89]	890	} else {
	891	integ++;
[2]	892	}
	893	}
[89]	894	oss << endl;
	895	oss << "Table contains " << table_.nrow() << " integration(s)." << endl;
	896	oss << "in " << count << " scan(s)." << endl;
	897	oss << "--------------------------------------------------";
	898	return String(oss);
[2]	899	}
[18]	900
	901	Int SDMemTable::nBeam() const {
	902	Int n;
	903	table_.keywordSet().get("nBeam",n);
	904	return n;
	905	}
	906	Int SDMemTable::nIF() const {
	907	Int n;
	908	table_.keywordSet().get("nIF",n);
	909	return n;
	910	}
	911	Int SDMemTable::nPol() const {
	912	Int n;
	913	table_.keywordSet().get("nPol",n);
	914	return n;
	915	}
	916	Int SDMemTable::nChan() const {
	917	Int n;
	918	table_.keywordSet().get("nChan",n);
	919	return n;
	920	}
[16]	921	/*
[18]	922	void SDMemTable::maskChannels(const std::vector<Int>& whichChans ) {
[89]	923
[16]	924	std::vector<int>::iterator it;
	925	ArrayAccessor<uChar, Axis<2> > j(flags_);
	926	for (it = whichChans.begin(); it != whichChans.end(); it++) {
	927	j.reset(j.begin(uInt(*it)));
	928	for (ArrayAccessor<uChar, Axis<0> > i(j); i != i.end(); ++i) {
	929	for (ArrayAccessor<uChar, Axis<1> > ii(i); ii != ii.end(); ++ii) {
[89]	930	for (ArrayAccessor<uChar, Axis<3> > iii(ii);
	931	iii != iii.end(); ++iii) {
	932	(*iii) =
	933	}
[16]	934	}
	935	}
	936	}
[89]	937
[16]	938	}
	939	*/
[89]	940	void SDMemTable::flag(int whichRow) {
	941	ArrayColumn<uChar> spec(table_, "FLAGTRA");
	942	Array<uChar> arr;
	943	spec.get(whichRow, arr);
	944
	945	ArrayAccessor<uChar, Axis<0> > aa0(arr);
	946	aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
	947	ArrayAccessor<uChar, Axis<1> > aa1(aa0);
	948	aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
	949	ArrayAccessor<uChar, Axis<2> > aa2(aa1);
	950	aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
	951
	952	for (ArrayAccessor<uChar, Axis<3> > i(aa2); i != i.end(); ++i) {
	953	(*i) = uChar(True);
	954	}
	955
	956	spec.put(whichRow, arr);
	957	}

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

Download in other formats: