Context Navigation

source: trunk/src/STMath.cpp@ 1400

Visit:

Last change on this file since 1400 was 1391, checked in by Malte Marquarding, 17 years ago
merge from alma branch to get alma/GBT support. Commented out fluxUnit changes as they are using a chnaged interface to PKSreader/writer. Also commented out progress meter related code.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 64.9 KB

Line
1	//
2	// C++ Implementation: STMath
3	//
4	// Description:
5	//
6	//
7	// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8	//
9	// Copyright: See COPYING file that comes with this distribution
10	//
11	//
12
13	#include <casa/iomanip.h>
14	#include <casa/Exceptions/Error.h>
15	#include <casa/Containers/Block.h>
16	#include <casa/BasicSL/String.h>
17	#include <casa/Arrays/MaskArrLogi.h>
18	#include <casa/Arrays/MaskArrMath.h>
19	#include <casa/Arrays/ArrayLogical.h>
20	#include <casa/Arrays/ArrayMath.h>
21	#include <casa/Arrays/Slice.h>
22	#include <casa/Arrays/Slicer.h>
23	#include <casa/Containers/RecordField.h>
24	#include <tables/Tables/TableRow.h>
25	#include <tables/Tables/TableVector.h>
26	#include <tables/Tables/TabVecMath.h>
27	#include <tables/Tables/ExprNode.h>
28	#include <tables/Tables/TableRecord.h>
29	#include <tables/Tables/TableParse.h>
30	#include <tables/Tables/ReadAsciiTable.h>
31	#include <tables/Tables/TableIter.h>
32	#include <tables/Tables/TableCopy.h>
33	#include <scimath/Mathematics/FFTServer.h>
34
35	#include <lattices/Lattices/LatticeUtilities.h>
36
37	#include <coordinates/Coordinates/SpectralCoordinate.h>
38	#include <coordinates/Coordinates/CoordinateSystem.h>
39	#include <coordinates/Coordinates/CoordinateUtil.h>
40	#include <coordinates/Coordinates/FrequencyAligner.h>
41
42	#include <scimath/Mathematics/VectorKernel.h>
43	#include <scimath/Mathematics/Convolver.h>
44	#include <scimath/Functionals/Polynomial.h>
45
46	#include "MathUtils.h"
47	#include "RowAccumulator.h"
48	#include "STAttr.h"
49	#include "STMath.h"
50	#include "STSelector.h"
51
52	using namespace casa;
53
54	using namespace asap;
55
56	STMath::STMath(bool insitu) :
57	insitu_(insitu)
58	{
59	}
60
61
62	STMath::~STMath()
63	{
64	}
65
66	CountedPtr<Scantable>
67	STMath::average( const std::vector<CountedPtr<Scantable> >& in,
68	const std::vector<bool>& mask,
69	const std::string& weight,
70	const std::string& avmode)
71	{
72	if ( avmode == "SCAN" && in.size() != 1 )
73	throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
74	"Use merge first."));
75	WeightType wtype = stringToWeight(weight);
76
77	// output
78	// clone as this is non insitu
79	bool insitu = insitu_;
80	setInsitu(false);
81	CountedPtr< Scantable > out = getScantable(in[0], true);
82	setInsitu(insitu);
83	std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
84	++stit;
85	while ( stit != in.end() ) {
86	out->appendToHistoryTable((*stit)->history());
87	++stit;
88	}
89
90	Table& tout = out->table();
91
92	/// @todo check if all scantables are conformant
93
94	ArrayColumn<Float> specColOut(tout,"SPECTRA");
95	ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
96	ArrayColumn<Float> tsysColOut(tout,"TSYS");
97	ScalarColumn<Double> mjdColOut(tout,"TIME");
98	ScalarColumn<Double> intColOut(tout,"INTERVAL");
99	ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
100	ScalarColumn<uInt> scanColOut(tout,"SCANNO");
101
102	// set up the output table rows. These are based on the structure of the
103	// FIRST scantable in the vector
104	const Table& baset = in[0]->table();
105
106	Block<String> cols(3);
107	cols[0] = String("BEAMNO");
108	cols[1] = String("IFNO");
109	cols[2] = String("POLNO");
110	if ( avmode == "SOURCE" ) {
111	cols.resize(4);
112	cols[3] = String("SRCNAME");
113	}
114	if ( avmode == "SCAN" && in.size() == 1) {
115	cols.resize(4);
116	cols[3] = String("SCANNO");
117	}
118	uInt outrowCount = 0;
119	TableIterator iter(baset, cols);
120	while (!iter.pastEnd()) {
121	Table subt = iter.table();
122	// copy the first row of this selection into the new table
123	tout.addRow();
124	TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
125	// re-index to 0
126	if ( avmode != "SCAN" && avmode != "SOURCE" ) {
127	scanColOut.put(outrowCount, uInt(0));
128	}
129	++outrowCount;
130	++iter;
131	}
132	RowAccumulator acc(wtype);
133	Vector<Bool> cmask(mask);
134	acc.setUserMask(cmask);
135	ROTableRow row(tout);
136	ROArrayColumn<Float> specCol, tsysCol;
137	ROArrayColumn<uChar> flagCol;
138	ROScalarColumn<Double> mjdCol, intCol;
139	ROScalarColumn<Int> scanIDCol;
140
141	Vector<uInt> rowstodelete;
142
143	for (uInt i=0; i < tout.nrow(); ++i) {
144	for ( int j=0; j < int(in.size()); ++j ) {
145	const Table& tin = in[j]->table();
146	const TableRecord& rec = row.get(i);
147	ROScalarColumn<Double> tmp(tin, "TIME");
148	Double td;tmp.get(0,td);
149	Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
150	&& tin.col("IFNO") == Int(rec.asuInt("IFNO"))
151	&& tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
152	Table subt;
153	if ( avmode == "SOURCE") {
154	subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
155	} else if (avmode == "SCAN") {
156	subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
157	} else {
158	subt = basesubt;
159	}
160	specCol.attach(subt,"SPECTRA");
161	flagCol.attach(subt,"FLAGTRA");
162	tsysCol.attach(subt,"TSYS");
163	intCol.attach(subt,"INTERVAL");
164	mjdCol.attach(subt,"TIME");
165	Vector<Float> spec,tsys;
166	Vector<uChar> flag;
167	Double inter,time;
168	for (uInt k = 0; k < subt.nrow(); ++k ) {
169	flagCol.get(k, flag);
170	Vector<Bool> bflag(flag.shape());
171	convertArray(bflag, flag);
172	/*
173	if ( allEQ(bflag, True) ) {
174	continue;//don't accumulate
175	}
176	*/
177	specCol.get(k, spec);
178	tsysCol.get(k, tsys);
179	intCol.get(k, inter);
180	mjdCol.get(k, time);
181	// spectrum has to be added last to enable weighting by the other values
182	acc.add(spec, !bflag, tsys, inter, time);
183	}
184	}
185	const Vector<Bool>& msk = acc.getMask();
186	if ( allEQ(msk, False) ) {
187	uint n = rowstodelete.nelements();
188	rowstodelete.resize(n+1, True);
189	rowstodelete[n] = i;
190	continue;
191	}
192	//write out
193	Vector<uChar> flg(msk.shape());
194	convertArray(flg, !msk);
195	flagColOut.put(i, flg);
196	specColOut.put(i, acc.getSpectrum());
197	tsysColOut.put(i, acc.getTsys());
198	intColOut.put(i, acc.getInterval());
199	mjdColOut.put(i, acc.getTime());
200	// we should only have one cycle now -> reset it to be 0
201	// frequency switched data has different CYCLENO for different IFNO
202	// which requires resetting this value
203	cycColOut.put(i, uInt(0));
204	acc.reset();
205	}
206	if (rowstodelete.nelements() > 0) {
207	cout << rowstodelete << endl;
208	tout.removeRow(rowstodelete);
209	if (tout.nrow() == 0) {
210	throw(AipsError("Can't average fully flagged data."));
211	}
212	}
213	return out;
214	}
215
216	CountedPtr< Scantable >
217	STMath::averageChannel( const CountedPtr < Scantable > & in,
218	const std::string & mode,
219	const std::string& avmode )
220	{
221	// clone as this is non insitu
222	bool insitu = insitu_;
223	setInsitu(false);
224	CountedPtr< Scantable > out = getScantable(in, true);
225	setInsitu(insitu);
226	Table& tout = out->table();
227	ArrayColumn<Float> specColOut(tout,"SPECTRA");
228	ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
229	ArrayColumn<Float> tsysColOut(tout,"TSYS");
230	ScalarColumn<uInt> scanColOut(tout,"SCANNO");
231	ScalarColumn<Double> intColOut(tout, "INTERVAL");
232	Table tmp = in->table().sort("BEAMNO");
233	Block<String> cols(3);
234	cols[0] = String("BEAMNO");
235	cols[1] = String("IFNO");
236	cols[2] = String("POLNO");
237	if ( avmode == "SCAN") {
238	cols.resize(4);
239	cols[3] = String("SCANNO");
240	}
241	uInt outrowCount = 0;
242	uChar userflag = 1 << 7;
243	TableIterator iter(tmp, cols);
244	while (!iter.pastEnd()) {
245	Table subt = iter.table();
246	ROArrayColumn<Float> specCol, tsysCol;
247	ROArrayColumn<uChar> flagCol;
248	ROScalarColumn<Double> intCol(subt, "INTERVAL");
249	specCol.attach(subt,"SPECTRA");
250	flagCol.attach(subt,"FLAGTRA");
251	tsysCol.attach(subt,"TSYS");
252	tout.addRow();
253	TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
254	if ( avmode != "SCAN") {
255	scanColOut.put(outrowCount, uInt(0));
256	}
257	Vector<Float> tmp;
258	specCol.get(0, tmp);
259	uInt nchan = tmp.nelements();
260	// have to do channel by channel here as MaskedArrMath
261	// doesn't have partialMedians
262	Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
263	Vector<Float> outspec(nchan);
264	Vector<uChar> outflag(nchan,0);
265	Vector<Float> outtsys(1);/// @fixme when tsys is channel based
266	for (uInt i=0; i<nchan; ++i) {
267	Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
268	MaskedArray<Float> ma = maskedArray(specs,flags);
269	outspec[i] = median(ma);
270	if ( allEQ(ma.getMask(), False) )
271	outflag[i] = userflag;// flag data
272	}
273	outtsys[0] = median(tsysCol.getColumn());
274	specColOut.put(outrowCount, outspec);
275	flagColOut.put(outrowCount, outflag);
276	tsysColOut.put(outrowCount, outtsys);
277	Double intsum = sum(intCol.getColumn());
278	intColOut.put(outrowCount, intsum);
279	++outrowCount;
280	++iter;
281	}
282	return out;
283	}
284
285	CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
286	bool droprows)
287	{
288	if (insitu_) return in;
289	else {
290	// clone
291	Scantable* tabp = new Scantable(*in, Bool(droprows));
292	return CountedPtr<Scantable>(tabp);
293	}
294	}
295
296	CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
297	float val,
298	const std::string& mode,
299	bool tsys )
300	{
301	CountedPtr< Scantable > out = getScantable(in, false);
302	Table& tab = out->table();
303	ArrayColumn<Float> specCol(tab,"SPECTRA");
304	ArrayColumn<Float> tsysCol(tab,"TSYS");
305	for (uInt i=0; i<tab.nrow(); ++i) {
306	Vector<Float> spec;
307	Vector<Float> ts;
308	specCol.get(i, spec);
309	tsysCol.get(i, ts);
310	if (mode == "MUL" \|\| mode == "DIV") {
311	if (mode == "DIV") val = 1.0/val;
312	spec *= val;
313	specCol.put(i, spec);
314	if ( tsys ) {
315	ts *= val;
316	tsysCol.put(i, ts);
317	}
318	} else if ( mode == "ADD" \|\| mode == "SUB") {
319	if (mode == "SUB") val *= -1.0;
320	spec += val;
321	specCol.put(i, spec);
322	if ( tsys ) {
323	ts += val;
324	tsysCol.put(i, ts);
325	}
326	}
327	}
328	return out;
329	}
330
331	CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
332	const CountedPtr<Scantable>& right,
333	const std::string& mode)
334	{
335	bool insitu = insitu_;
336	if ( ! left->conformant(*right) ) {
337	throw(AipsError("'left' and 'right' scantables are not conformant."));
338	}
339	setInsitu(false);
340	CountedPtr< Scantable > out = getScantable(left, false);
341	setInsitu(insitu);
342	Table& tout = out->table();
343	Block<String> coln(5);
344	coln[0] = "SCANNO"; coln[1] = "CYCLENO"; coln[2] = "BEAMNO";
345	coln[3] = "IFNO"; coln[4] = "POLNO";
346	Table tmpl = tout.sort(coln);
347	Table tmpr = right->table().sort(coln);
348	ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
349	ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
350	ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
351	ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
352
353	for (uInt i=0; i<tout.nrow(); ++i) {
354	Vector<Float> lspecvec, rspecvec;
355	Vector<uChar> lflagvec, rflagvec;
356	lspecvec = lspecCol(i); rspecvec = rspecCol(i);
357	lflagvec = lflagCol(i); rflagvec = rflagCol(i);
358	MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
359	MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
360	if (mode == "ADD") {
361	mleft += mright;
362	} else if ( mode == "SUB") {
363	mleft -= mright;
364	} else if ( mode == "MUL") {
365	mleft *= mright;
366	} else if ( mode == "DIV") {
367	mleft /= mright;
368	} else {
369	throw(AipsError("Illegal binary operator"));
370	}
371	lspecCol.put(i, mleft.getArray());
372	}
373	return out;
374	}
375
376
377
378	MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
379	const Vector<uChar>& f)
380	{
381	Vector<Bool> mask;
382	mask.resize(f.shape());
383	convertArray(mask, f);
384	return MaskedArray<Float>(s,!mask);
385	}
386
387	Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
388	{
389	const Vector<Bool>& m = ma.getMask();
390	Vector<uChar> flags(m.shape());
391	convertArray(flags, !m);
392	return flags;
393	}
394
395	CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
396	const std::string & mode,
397	bool preserve )
398	{
399	/// @todo make other modes available
400	/// modes should be "nearest", "pair"
401	// make this operation non insitu
402	const Table& tin = in->table();
403	Table ons = tin(tin.col("SRCTYPE") == Int(0));
404	Table offs = tin(tin.col("SRCTYPE") == Int(1));
405	if ( offs.nrow() == 0 )
406	throw(AipsError("No 'off' scans present."));
407	// put all "on" scans into output table
408
409	bool insitu = insitu_;
410	setInsitu(false);
411	CountedPtr< Scantable > out = getScantable(in, true);
412	setInsitu(insitu);
413	Table& tout = out->table();
414
415	TableCopy::copyRows(tout, ons);
416	TableRow row(tout);
417	ROScalarColumn<Double> offtimeCol(offs, "TIME");
418	ArrayColumn<Float> outspecCol(tout, "SPECTRA");
419	ROArrayColumn<Float> outtsysCol(tout, "TSYS");
420	ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
421	for (uInt i=0; i < tout.nrow(); ++i) {
422	const TableRecord& rec = row.get(i);
423	Double ontime = rec.asDouble("TIME");
424	Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
425	&& offs.col("IFNO") == Int(rec.asuInt("IFNO"))
426	&& offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
427	ROScalarColumn<Double> offtimeCol(presel, "TIME");
428
429	Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
430	// Timestamp may vary within a cycle ???!!!
431	// increase this by 0.01 sec in case of rounding errors...
432	// There might be a better way to do this.
433	// fix to this fix. TIME is MJD, so 1.0d not 1.0s
434	mindeltat += 0.01/24./60./60.;
435	Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
436
437	if ( sel.nrow() < 1 ) {
438	throw(AipsError("No closest in time found... This could be a rounding "
439	"issue. Try quotient instead."));
440	}
441	TableRow offrow(sel);
442	const TableRecord& offrec = offrow.get(0);//should only be one row
443	RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
444	RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
445	RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
446	/// @fixme this assumes tsys is a scalar not vector
447	Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
448	Vector<Float> specon, tsyson;
449	outtsysCol.get(i, tsyson);
450	outspecCol.get(i, specon);
451	Vector<uChar> flagon;
452	outflagCol.get(i, flagon);
453	MaskedArray<Float> mon = maskedArray(specon, flagon);
454	MaskedArray<Float> moff = maskedArray(specoff, flagoff);
455	MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
456	if (preserve) {
457	quot -= tsysoffscalar;
458	} else {
459	quot -= tsyson[0];
460	}
461	outspecCol.put(i, quot.getArray());
462	outflagCol.put(i, flagsFromMA(quot));
463	}
464	// renumber scanno
465	TableIterator it(tout, "SCANNO");
466	uInt i = 0;
467	while ( !it.pastEnd() ) {
468	Table t = it.table();
469	TableVector<uInt> vec(t, "SCANNO");
470	vec = i;
471	++i;
472	++it;
473	}
474	return out;
475	}
476
477
478	CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
479	const CountedPtr< Scantable > & off,
480	bool preserve )
481	{
482	bool insitu = insitu_;
483	if ( ! on->conformant(*off) ) {
484	throw(AipsError("'on' and 'off' scantables are not conformant."));
485	}
486	setInsitu(false);
487	CountedPtr< Scantable > out = getScantable(on, false);
488	setInsitu(insitu);
489	Table& tout = out->table();
490	const Table& toff = off->table();
491	TableIterator sit(tout, "SCANNO");
492	TableIterator s2it(toff, "SCANNO");
493	while ( !sit.pastEnd() ) {
494	Table ton = sit.table();
495	TableRow row(ton);
496	Table t = s2it.table();
497	ArrayColumn<Float> outspecCol(ton, "SPECTRA");
498	ROArrayColumn<Float> outtsysCol(ton, "TSYS");
499	ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
500	for (uInt i=0; i < ton.nrow(); ++i) {
501	const TableRecord& rec = row.get(i);
502	Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
503	&& t.col("IFNO") == Int(rec.asuInt("IFNO"))
504	&& t.col("POLNO") == Int(rec.asuInt("POLNO")) );
505	if ( offsel.nrow() == 0 )
506	throw AipsError("STMath::quotient: no matching off");
507	TableRow offrow(offsel);
508	const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
509	RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
510	RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
511	RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
512	Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
513	Vector<Float> specon, tsyson;
514	outtsysCol.get(i, tsyson);
515	outspecCol.get(i, specon);
516	Vector<uChar> flagon;
517	outflagCol.get(i, flagon);
518	MaskedArray<Float> mon = maskedArray(specon, flagon);
519	MaskedArray<Float> moff = maskedArray(specoff, flagoff);
520	MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
521	if (preserve) {
522	quot -= tsysoffscalar;
523	} else {
524	quot -= tsyson[0];
525	}
526	outspecCol.put(i, quot.getArray());
527	outflagCol.put(i, flagsFromMA(quot));
528	}
529	++sit;
530	++s2it;
531	// take the first off for each on scan which doesn't have a
532	// matching off scan
533	// non <= noff: matching pairs, non > noff matching pairs then first off
534	if ( s2it.pastEnd() ) s2it.reset();
535	}
536	return out;
537	}
538
539	// dototalpower (migration of GBTIDL procedure dototalpower.pro)
540	// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
541	// do it for each cycles in a specific scan.
542	CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
543	const CountedPtr< Scantable >& caloff, Float tcal )
544	{
545	if ( ! calon->conformant(*caloff) ) {
546	throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
547	}
548	setInsitu(false);
549	CountedPtr< Scantable > out = getScantable(caloff, false);
550	Table& tout = out->table();
551	const Table& tcon = calon->table();
552	Vector<Float> tcalout;
553	Vector<Float> tcalout2; //debug
554
555	if ( tout.nrow() != tcon.nrow() ) {
556	throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
557	}
558	// iteration by scanno or cycle no.
559	TableIterator sit(tout, "SCANNO");
560	TableIterator s2it(tcon, "SCANNO");
561	while ( !sit.pastEnd() ) {
562	Table toff = sit.table();
563	TableRow row(toff);
564	Table t = s2it.table();
565	ScalarColumn<Double> outintCol(toff, "INTERVAL");
566	ArrayColumn<Float> outspecCol(toff, "SPECTRA");
567	ArrayColumn<Float> outtsysCol(toff, "TSYS");
568	ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
569	ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
570	ROScalarColumn<uInt> outpolCol(toff, "POLNO");
571	ROScalarColumn<Double> onintCol(t, "INTERVAL");
572	ROArrayColumn<Float> onspecCol(t, "SPECTRA");
573	ROArrayColumn<Float> ontsysCol(t, "TSYS");
574	ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
575	//ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
576
577	for (uInt i=0; i < toff.nrow(); ++i) {
578	//skip these checks -> assumes the data order are the same between the cal on off pairs
579	//
580	Vector<Float> specCalon, specCaloff;
581	// to store scalar (mean) tsys
582	Vector<Float> tsysout(1);
583	uInt tcalId, polno;
584	Double offint, onint;
585	outpolCol.get(i, polno);
586	outspecCol.get(i, specCaloff);
587	onspecCol.get(i, specCalon);
588	Vector<uChar> flagCaloff, flagCalon;
589	outflagCol.get(i, flagCaloff);
590	onflagCol.get(i, flagCalon);
591	outtcalIdCol.get(i, tcalId);
592	outintCol.get(i, offint);
593	onintCol.get(i, onint);
594	// caluculate mean Tsys
595	uInt nchan = specCaloff.nelements();
596	// percentage of edge cut off
597	uInt pc = 10;
598	uInt bchan = nchan/pc;
599	uInt echan = nchan-bchan;
600
601	Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
602	Vector<Float> testsubsp = specCaloff(chansl);
603	MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
604	MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
605	MaskedArray<Float> spdiff = spon-spoff;
606	uInt noff = spoff.nelementsValid();
607	//uInt non = spon.nelementsValid();
608	uInt ndiff = spdiff.nelementsValid();
609	Float meantsys;
610
611	/**
612	Double subspec, subdiff;
613	uInt usednchan;
614	subspec = 0;
615	subdiff = 0;
616	usednchan = 0;
617	for(uInt k=(bchan-1); k<echan; k++) {
618	subspec += specCaloff[k];
619	subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
620	++usednchan;
621	}
622	**/
623	// get tcal if input tcal <= 0
624	String tcalt;
625	Float tcalUsed;
626	tcalUsed = tcal;
627	if ( tcal <= 0.0 ) {
628	caloff->tcal().getEntry(tcalt, tcalout, tcalId);
629	if (polno<=3) {
630	tcalUsed = tcalout[polno];
631	}
632	else {
633	tcalUsed = tcalout[0];
634	}
635	}
636
637	Float meanoff;
638	Float meandiff;
639	if (noff && ndiff) {
640	//Debug
641	//if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
642	meanoff = sum(spoff)/noff;
643	meandiff = sum(spdiff)/ndiff;
644	meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
645	}
646	else {
647	meantsys=1;
648	}
649
650	tsysout[0] = Float(meantsys);
651	MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
652	MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
653	MaskedArray<Float> sig = Float(0.5) * (mcaloff + mcalon);
654	//uInt ncaloff = mcaloff.nelementsValid();
655	//uInt ncalon = mcalon.nelementsValid();
656
657	outintCol.put(i, offint+onint);
658	outspecCol.put(i, sig.getArray());
659	outflagCol.put(i, flagsFromMA(sig));
660	outtsysCol.put(i, tsysout);
661	}
662	++sit;
663	++s2it;
664	}
665	return out;
666	}
667
668	//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
669	// observatiions.
670	// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
671	// no smoothing).
672	// output: resultant scantable [= (sig-ref/ref)*tsys]
673	CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
674	const CountedPtr < Scantable >& ref,
675	int smoothref,
676	casa::Float tsysv,
677	casa::Float tau )
678	{
679	if ( ! ref->conformant(*sig) ) {
680	throw(AipsError("'sig' and 'ref' scantables are not conformant."));
681	}
682	setInsitu(false);
683	CountedPtr< Scantable > out = getScantable(sig, false);
684	CountedPtr< Scantable > smref;
685	if ( smoothref > 1 ) {
686	float fsmoothref = static_cast<float>(smoothref);
687	std::string inkernel = "boxcar";
688	smref = smooth(ref, inkernel, fsmoothref );
689	ostringstream oss;
690	oss<<"Applied smoothing of "<<fsmoothref<<" on the reference."<<endl;
691	pushLog(String(oss));
692	}
693	else {
694	smref = ref;
695	}
696	Table& tout = out->table();
697	const Table& tref = smref->table();
698	if ( tout.nrow() != tref.nrow() ) {
699	throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
700	}
701	// iteration by scanno? or cycle no.
702	TableIterator sit(tout, "SCANNO");
703	TableIterator s2it(tref, "SCANNO");
704	while ( !sit.pastEnd() ) {
705	Table ton = sit.table();
706	Table t = s2it.table();
707	ScalarColumn<Double> outintCol(ton, "INTERVAL");
708	ArrayColumn<Float> outspecCol(ton, "SPECTRA");
709	ArrayColumn<Float> outtsysCol(ton, "TSYS");
710	ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
711	ArrayColumn<Float> refspecCol(t, "SPECTRA");
712	ROScalarColumn<Double> refintCol(t, "INTERVAL");
713	ROArrayColumn<Float> reftsysCol(t, "TSYS");
714	ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
715	ROScalarColumn<Float> refelevCol(t, "ELEVATION");
716	for (uInt i=0; i < ton.nrow(); ++i) {
717
718	Double onint, refint;
719	Vector<Float> specon, specref;
720	// to store scalar (mean) tsys
721	Vector<Float> tsysref;
722	outintCol.get(i, onint);
723	refintCol.get(i, refint);
724	outspecCol.get(i, specon);
725	refspecCol.get(i, specref);
726	Vector<uChar> flagref, flagon;
727	outflagCol.get(i, flagon);
728	refflagCol.get(i, flagref);
729	reftsysCol.get(i, tsysref);
730
731	Float tsysrefscalar;
732	if ( tsysv > 0.0 ) {
733	ostringstream oss;
734	Float elev;
735	refelevCol.get(i, elev);
736	oss << "user specified Tsys = " << tsysv;
737	// do recalc elevation if EL = 0
738	if ( elev == 0 ) {
739	throw(AipsError("EL=0, elevation data is missing."));
740	} else {
741	if ( tau <= 0.0 ) {
742	throw(AipsError("Valid tau is not supplied."));
743	} else {
744	tsysrefscalar = tsysv * exp(tau/elev);
745	}
746	}
747	oss << ", corrected (for El) tsys= "<<tsysrefscalar;
748	pushLog(String(oss));
749	}
750	else {
751	tsysrefscalar = tsysref[0];
752	}
753	//get quotient spectrum
754	MaskedArray<Float> mref = maskedArray(specref, flagref);
755	MaskedArray<Float> mon = maskedArray(specon, flagon);
756	MaskedArray<Float> specres = tsysrefscalar*((mon - mref)/mref);
757	Double resint = onintrefintsmoothref/(onint+refint*smoothref);
758
759	//Debug
760	//cerr<<"Tsys used="<<tsysrefscalar<<endl;
761	// fill the result, replay signal tsys by reference tsys
762	outintCol.put(i, resint);
763	outspecCol.put(i, specres.getArray());
764	outflagCol.put(i, flagsFromMA(specres));
765	outtsysCol.put(i, tsysref);
766	}
767	++sit;
768	++s2it;
769	}
770	return out;
771	}
772
773	CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
774	const std::vector<int>& scans,
775	int smoothref,
776	casa::Float tsysv,
777	casa::Float tau,
778	casa::Float tcal )
779
780	{
781	setInsitu(false);
782	STSelector sel;
783	std::vector<int> scan1, scan2, beams;
784	std::vector< vector<int> > scanpair;
785	std::vector<string> calstate;
786	String msg;
787
788	CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
789	CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
790
791	std::vector< CountedPtr< Scantable > > sctables;
792	sctables.push_back(s1b1on);
793	sctables.push_back(s1b1off);
794	sctables.push_back(s1b2on);
795	sctables.push_back(s1b2off);
796	sctables.push_back(s2b1on);
797	sctables.push_back(s2b1off);
798	sctables.push_back(s2b2on);
799	sctables.push_back(s2b2off);
800
801	//check scanlist
802	int n=s->checkScanInfo(scans);
803	if (n==1) {
804	throw(AipsError("Incorrect scan pairs. "));
805	}
806
807	// Assume scans contain only a pair of consecutive scan numbers.
808	// It is assumed that first beam, b1, is on target.
809	// There is no check if the first beam is on or not.
810	if ( scans.size()==1 ) {
811	scan1.push_back(scans[0]);
812	scan2.push_back(scans[0]+1);
813	} else if ( scans.size()==2 ) {
814	scan1.push_back(scans[0]);
815	scan2.push_back(scans[1]);
816	} else {
817	if ( scans.size()%2 == 0 ) {
818	for (uInt i=0; i<scans.size(); i++) {
819	if (i%2 == 0) {
820	scan1.push_back(scans[i]);
821	}
822	else {
823	scan2.push_back(scans[i]);
824	}
825	}
826	} else {
827	throw(AipsError("Odd numbers of scans, cannot form pairs."));
828	}
829	}
830	scanpair.push_back(scan1);
831	scanpair.push_back(scan2);
832	calstate.push_back("*calon");
833	calstate.push_back("*[^calon]");
834	CountedPtr< Scantable > ws = getScantable(s, false);
835	uInt l=0;
836	while ( l < sctables.size() ) {
837	for (uInt i=0; i < 2; i++) {
838	for (uInt j=0; j < 2; j++) {
839	for (uInt k=0; k < 2; k++) {
840	sel.reset();
841	sel.setScans(scanpair[i]);
842	sel.setName(calstate[k]);
843	beams.clear();
844	beams.push_back(j);
845	sel.setBeams(beams);
846	ws->setSelection(sel);
847	sctables[l]= getScantable(ws, false);
848	l++;
849	}
850	}
851	}
852	}
853
854	// replace here by splitData or getData functionality
855	CountedPtr< Scantable > sig1;
856	CountedPtr< Scantable > ref1;
857	CountedPtr< Scantable > sig2;
858	CountedPtr< Scantable > ref2;
859	CountedPtr< Scantable > calb1;
860	CountedPtr< Scantable > calb2;
861
862	msg=String("Processing dototalpower for subset of the data");
863	ostringstream oss1;
864	oss1 << msg << endl;
865	pushLog(String(oss1));
866	// Debug for IRC CS data
867	//float tcal1=7.0;
868	//float tcal2=4.0;
869	sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
870	ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
871	ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
872	sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
873
874	// correction of user-specified tsys for elevation here
875
876	// dosigref calibration
877	msg=String("Processing dosigref for subset of the data");
878	ostringstream oss2;
879	oss2 << msg << endl;
880	pushLog(String(oss2));
881	calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
882	calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
883
884	// iteration by scanno or cycle no.
885	Table& tcalb1 = calb1->table();
886	Table& tcalb2 = calb2->table();
887	TableIterator sit(tcalb1, "SCANNO");
888	TableIterator s2it(tcalb2, "SCANNO");
889	while ( !sit.pastEnd() ) {
890	Table t1 = sit.table();
891	Table t2= s2it.table();
892	ArrayColumn<Float> outspecCol(t1, "SPECTRA");
893	ArrayColumn<Float> outtsysCol(t1, "TSYS");
894	ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
895	ScalarColumn<Double> outintCol(t1, "INTERVAL");
896	ArrayColumn<Float> t2specCol(t2, "SPECTRA");
897	ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
898	ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
899	ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
900	for (uInt i=0; i < t1.nrow(); ++i) {
901	Vector<Float> spec1, spec2;
902	// to store scalar (mean) tsys
903	Vector<Float> tsys1, tsys2;
904	Vector<uChar> flag1, flag2;
905	Double tint1, tint2;
906	outspecCol.get(i, spec1);
907	t2specCol.get(i, spec2);
908	outflagCol.get(i, flag1);
909	t2flagCol.get(i, flag2);
910	outtsysCol.get(i, tsys1);
911	t2tsysCol.get(i, tsys2);
912	outintCol.get(i, tint1);
913	t2intCol.get(i, tint2);
914	// average
915	// assume scalar tsys for weights
916	Float wt1, wt2, tsyssq1, tsyssq2;
917	tsyssq1 = tsys1[0]*tsys1[0];
918	tsyssq2 = tsys2[0]*tsys2[0];
919	wt1 = Float(tint1)/tsyssq1;
920	wt2 = Float(tint2)/tsyssq2;
921	Float invsumwt=1/(wt1+wt2);
922	MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
923	MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
924	MaskedArray<Float> avspec = invsumwt * (wt1mspec1 + wt2mspec2);
925	//Array<Float> avtsys = Float(0.5) * (tsys1 + tsys2);
926	// cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
927	tsys1[0] = sqrt(tsyssq1 + tsyssq2);
928	Array<Float> avtsys = tsys1;
929
930	outspecCol.put(i, avspec.getArray());
931	outflagCol.put(i, flagsFromMA(avspec));
932	outtsysCol.put(i, avtsys);
933	}
934	++sit;
935	++s2it;
936	}
937	return calb1;
938	}
939
940	//GBTIDL version of frequency switched data calibration
941	CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
942	const std::vector<int>& scans,
943	int smoothref,
944	casa::Float tsysv,
945	casa::Float tau,
946	casa::Float tcal )
947	{
948
949
950	STSelector sel;
951	CountedPtr< Scantable > ws = getScantable(s, false);
952	CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
953	CountedPtr< Scantable > calsig, calref, out;
954
955	//split the data
956	sel.setName("*_fs");
957	ws->setSelection(sel);
958	sig = getScantable(ws,false);
959	sel.reset();
960	sel.setName("*_fs_calon");
961	ws->setSelection(sel);
962	sigwcal = getScantable(ws,false);
963	sel.reset();
964	sel.setName("*_fsr");
965	ws->setSelection(sel);
966	ref = getScantable(ws,false);
967	sel.reset();
968	sel.setName("*_fsr_calon");
969	ws->setSelection(sel);
970	refwcal = getScantable(ws,false);
971
972	calsig = dototalpower(sigwcal, sig, tcal=tcal);
973	calref = dototalpower(refwcal, ref, tcal=tcal);
974
975	out=dosigref(calsig,calref,smoothref,tsysv,tau);
976
977	return out;
978	}
979
980
981	CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
982	{
983	// make copy or reference
984	CountedPtr< Scantable > out = getScantable(in, false);
985	Table& tout = out->table();
986	Block<String> cols(4);
987	cols[0] = String("SCANNO");
988	cols[1] = String("CYCLENO");
989	cols[2] = String("BEAMNO");
990	cols[3] = String("POLNO");
991	TableIterator iter(tout, cols);
992	while (!iter.pastEnd()) {
993	Table subt = iter.table();
994	// this should leave us with two rows for the two IFs....if not ignore
995	if (subt.nrow() != 2 ) {
996	continue;
997	}
998	ArrayColumn<Float> specCol(subt, "SPECTRA");
999	ArrayColumn<Float> tsysCol(subt, "TSYS");
1000	ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1001	Vector<Float> onspec,offspec, ontsys, offtsys;
1002	Vector<uChar> onflag, offflag;
1003	tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
1004	specCol.get(0, onspec); specCol.get(1, offspec);
1005	flagCol.get(0, onflag); flagCol.get(1, offflag);
1006	MaskedArray<Float> on = maskedArray(onspec, onflag);
1007	MaskedArray<Float> off = maskedArray(offspec, offflag);
1008	MaskedArray<Float> oncopy = on.copy();
1009
1010	on /= off; on -= 1.0f;
1011	on *= ontsys[0];
1012	off /= oncopy; off -= 1.0f;
1013	off *= offtsys[0];
1014	specCol.put(0, on.getArray());
1015	const Vector<Bool>& m0 = on.getMask();
1016	Vector<uChar> flags0(m0.shape());
1017	convertArray(flags0, !m0);
1018	flagCol.put(0, flags0);
1019
1020	specCol.put(1, off.getArray());
1021	const Vector<Bool>& m1 = off.getMask();
1022	Vector<uChar> flags1(m1.shape());
1023	convertArray(flags1, !m1);
1024	flagCol.put(1, flags1);
1025	++iter;
1026	}
1027
1028	return out;
1029	}
1030
1031	std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1032	const std::vector< bool > & mask,
1033	const std::string& which )
1034	{
1035
1036	Vector<Bool> m(mask);
1037	const Table& tab = in->table();
1038	ROArrayColumn<Float> specCol(tab, "SPECTRA");
1039	ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1040	std::vector<float> out;
1041	for (uInt i=0; i < tab.nrow(); ++i ) {
1042	Vector<Float> spec; specCol.get(i, spec);
1043	Vector<uChar> flag; flagCol.get(i, flag);
1044	MaskedArray<Float> ma = maskedArray(spec, flag);
1045	float outstat = 0.0;
1046	if ( spec.nelements() == m.nelements() ) {
1047	outstat = mathutil::statistics(which, ma(m));
1048	} else {
1049	outstat = mathutil::statistics(which, ma);
1050	}
1051	out.push_back(outstat);
1052	}
1053	return out;
1054	}
1055
1056	CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
1057	int width )
1058	{
1059	if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
1060	CountedPtr< Scantable > out = getScantable(in, false);
1061	Table& tout = out->table();
1062	out->frequencies().rescale(width, "BIN");
1063	ArrayColumn<Float> specCol(tout, "SPECTRA");
1064	ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
1065	for (uInt i=0; i < tout.nrow(); ++i ) {
1066	MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
1067	MaskedArray<Float> maout;
1068	LatticeUtilities::bin(maout, main, 0, Int(width));
1069	/// @todo implement channel based tsys binning
1070	specCol.put(i, maout.getArray());
1071	flagCol.put(i, flagsFromMA(maout));
1072	// take only the first binned spectrum's length for the deprecated
1073	// global header item nChan
1074	if (i==0) tout.rwKeywordSet().define(String("nChan"),
1075	Int(maout.getArray().nelements()));
1076	}
1077	return out;
1078	}
1079
1080	CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
1081	const std::string& method,
1082	float width )
1083	//
1084	// Should add the possibility of width being specified in km/s. This means
1085	// that for each freqID (SpectralCoordinate) we will need to convert to an
1086	// average channel width (say at the reference pixel). Then we would need
1087	// to be careful to make sure each spectrum (of different freqID)
1088	// is the same length.
1089	//
1090	{
1091	//InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1092	Int interpMethod(stringToIMethod(method));
1093
1094	CountedPtr< Scantable > out = getScantable(in, false);
1095	Table& tout = out->table();
1096
1097	// Resample SpectralCoordinates (one per freqID)
1098	out->frequencies().rescale(width, "RESAMPLE");
1099	TableIterator iter(tout, "IFNO");
1100	TableRow row(tout);
1101	while ( !iter.pastEnd() ) {
1102	Table tab = iter.table();
1103	ArrayColumn<Float> specCol(tab, "SPECTRA");
1104	//ArrayColumn<Float> tsysCol(tout, "TSYS");
1105	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1106	Vector<Float> spec;
1107	Vector<uChar> flag;
1108	specCol.get(0,spec); // the number of channels should be constant per IF
1109	uInt nChanIn = spec.nelements();
1110	Vector<Float> xIn(nChanIn); indgen(xIn);
1111	Int fac = Int(nChanIn/width);
1112	Vector<Float> xOut(fac+10); // 10 to be safe - resize later
1113	uInt k = 0;
1114	Float x = 0.0;
1115	while (x < Float(nChanIn) ) {
1116	xOut(k) = x;
1117	k++;
1118	x += width;
1119	}
1120	uInt nChanOut = k;
1121	xOut.resize(nChanOut, True);
1122	// process all rows for this IFNO
1123	Vector<Float> specOut;
1124	Vector<Bool> maskOut;
1125	Vector<uChar> flagOut;
1126	for (uInt i=0; i < tab.nrow(); ++i) {
1127	specCol.get(i, spec);
1128	flagCol.get(i, flag);
1129	Vector<Bool> mask(flag.nelements());
1130	convertArray(mask, flag);
1131
1132	IPosition shapeIn(spec.shape());
1133	//sh.nchan = nChanOut;
1134	InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
1135	xIn, spec, mask,
1136	interpMethod, True, True);
1137	/// @todo do the same for channel based Tsys
1138	flagOut.resize(maskOut.nelements());
1139	convertArray(flagOut, maskOut);
1140	specCol.put(i, specOut);
1141	flagCol.put(i, flagOut);
1142	}
1143	++iter;
1144	}
1145
1146	return out;
1147	}
1148
1149	STMath::imethod STMath::stringToIMethod(const std::string& in)
1150	{
1151	static STMath::imap lookup;
1152
1153	// initialize the lookup table if necessary
1154	if ( lookup.empty() ) {
1155	lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
1156	lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
1157	lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
1158	lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
1159	}
1160
1161	STMath::imap::const_iterator iter = lookup.find(in);
1162
1163	if ( lookup.end() == iter ) {
1164	std::string message = in;
1165	message += " is not a valid interpolation mode";
1166	throw(AipsError(message));
1167	}
1168	return iter->second;
1169	}
1170
1171	WeightType STMath::stringToWeight(const std::string& in)
1172	{
1173	static std::map<std::string, WeightType> lookup;
1174
1175	// initialize the lookup table if necessary
1176	if ( lookup.empty() ) {
1177	lookup["NONE"] = asap::NONE;
1178	lookup["TINT"] = asap::TINT;
1179	lookup["TINTSYS"] = asap::TINTSYS;
1180	lookup["TSYS"] = asap::TSYS;
1181	lookup["VAR"] = asap::VAR;
1182	}
1183
1184	std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
1185
1186	if ( lookup.end() == iter ) {
1187	std::string message = in;
1188	message += " is not a valid weighting mode";
1189	throw(AipsError(message));
1190	}
1191	return iter->second;
1192	}
1193
1194	CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
1195	const vector< float > & coeff,
1196	const std::string & filename,
1197	const std::string& method)
1198	{
1199	// Get elevation data from Scantable and convert to degrees
1200	CountedPtr< Scantable > out = getScantable(in, false);
1201	Table& tab = out->table();
1202	ROScalarColumn<Float> elev(tab, "ELEVATION");
1203	Vector<Float> x = elev.getColumn();
1204	x *= Float(180 / C::pi); // Degrees
1205
1206	Vector<Float> coeffs(coeff);
1207	const uInt nc = coeffs.nelements();
1208	if ( filename.length() > 0 && nc > 0 ) {
1209	throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1210	}
1211
1212	// Correct
1213	if ( nc > 0 \|\| filename.length() == 0 ) {
1214	// Find instrument
1215	Bool throwit = True;
1216	Instrument inst =
1217	STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1218	throwit);
1219
1220	// Set polynomial
1221	Polynomial<Float>* ppoly = 0;
1222	Vector<Float> coeff;
1223	String msg;
1224	if ( nc > 0 ) {
1225	ppoly = new Polynomial<Float>(nc);
1226	coeff = coeffs;
1227	msg = String("user");
1228	} else {
1229	STAttr sdAttr;
1230	coeff = sdAttr.gainElevationPoly(inst);
1231	ppoly = new Polynomial<Float>(3);
1232	msg = String("built in");
1233	}
1234
1235	if ( coeff.nelements() > 0 ) {
1236	ppoly->setCoefficients(coeff);
1237	} else {
1238	delete ppoly;
1239	throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
1240	}
1241	ostringstream oss;
1242	oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
1243	oss << " " << coeff;
1244	pushLog(String(oss));
1245	const uInt nrow = tab.nrow();
1246	Vector<Float> factor(nrow);
1247	for ( uInt i=0; i < nrow; ++i ) {
1248	factor[i] = 1.0 / (*ppoly)(x[i]);
1249	}
1250	delete ppoly;
1251	scaleByVector(tab, factor, true);
1252
1253	} else {
1254	// Read and correct
1255	pushLog("Making correction from ascii Table");
1256	scaleFromAsciiTable(tab, filename, method, x, true);
1257	}
1258	return out;
1259	}
1260
1261	void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
1262	const std::string& method,
1263	const Vector<Float>& xout, bool dotsys)
1264	{
1265
1266	// Read gain-elevation ascii file data into a Table.
1267
1268	String formatString;
1269	Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
1270	scaleFromTable(in, tbl, method, xout, dotsys);
1271	}
1272
1273	void STMath::scaleFromTable(Table& in,
1274	const Table& table,
1275	const std::string& method,
1276	const Vector<Float>& xout, bool dotsys)
1277	{
1278
1279	ROScalarColumn<Float> geElCol(table, "ELEVATION");
1280	ROScalarColumn<Float> geFacCol(table, "FACTOR");
1281	Vector<Float> xin = geElCol.getColumn();
1282	Vector<Float> yin = geFacCol.getColumn();
1283	Vector<Bool> maskin(xin.nelements(),True);
1284
1285	// Interpolate (and extrapolate) with desired method
1286
1287	InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1288
1289	Vector<Float> yout;
1290	Vector<Bool> maskout;
1291	InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
1292	xin, yin, maskin, interp,
1293	True, True);
1294
1295	scaleByVector(in, Float(1.0)/yout, dotsys);
1296	}
1297
1298	void STMath::scaleByVector( Table& in,
1299	const Vector< Float >& factor,
1300	bool dotsys )
1301	{
1302	uInt nrow = in.nrow();
1303	if ( factor.nelements() != nrow ) {
1304	throw(AipsError("factors.nelements() != table.nelements()"));
1305	}
1306	ArrayColumn<Float> specCol(in, "SPECTRA");
1307	ArrayColumn<uChar> flagCol(in, "FLAGTRA");
1308	ArrayColumn<Float> tsysCol(in, "TSYS");
1309	for (uInt i=0; i < nrow; ++i) {
1310	MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1311	ma *= factor[i];
1312	specCol.put(i, ma.getArray());
1313	flagCol.put(i, flagsFromMA(ma));
1314	if ( dotsys ) {
1315	Vector<Float> tsys = tsysCol(i);
1316	tsys *= factor[i];
1317	tsysCol.put(i,tsys);
1318	}
1319	}
1320	}
1321
1322	CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
1323	float d, float etaap,
1324	float jyperk )
1325	{
1326	CountedPtr< Scantable > out = getScantable(in, false);
1327	Table& tab = in->table();
1328	Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
1329	Unit K(String("K"));
1330	Unit JY(String("Jy"));
1331
1332	bool tokelvin = true;
1333	Double cfac = 1.0;
1334
1335	if ( fluxUnit == JY ) {
1336	pushLog("Converting to K");
1337	Quantum<Double> t(1.0,fluxUnit);
1338	Quantum<Double> t2 = t.get(JY);
1339	cfac = (t2 / t).getValue(); // value to Jy
1340
1341	tokelvin = true;
1342	out->setFluxUnit("K");
1343	} else if ( fluxUnit == K ) {
1344	pushLog("Converting to Jy");
1345	Quantum<Double> t(1.0,fluxUnit);
1346	Quantum<Double> t2 = t.get(K);
1347	cfac = (t2 / t).getValue(); // value to K
1348
1349	tokelvin = false;
1350	out->setFluxUnit("Jy");
1351	} else {
1352	throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
1353	}
1354	// Make sure input values are converted to either Jy or K first...
1355	Float factor = cfac;
1356
1357	// Select method
1358	if (jyperk > 0.0) {
1359	factor *= jyperk;
1360	if ( tokelvin ) factor = 1.0 / jyperk;
1361	ostringstream oss;
1362	oss << "Jy/K = " << jyperk;
1363	pushLog(String(oss));
1364	Vector<Float> factors(tab.nrow(), factor);
1365	scaleByVector(tab,factors, false);
1366	} else if ( etaap > 0.0) {
1367	if (d < 0) {
1368	Instrument inst =
1369	STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1370	True);
1371	STAttr sda;
1372	d = sda.diameter(inst);
1373	}
1374	jyperk = STAttr::findJyPerK(etaap, d);
1375	ostringstream oss;
1376	oss << "Jy/K = " << jyperk;
1377	pushLog(String(oss));
1378	factor *= jyperk;
1379	if ( tokelvin ) {
1380	factor = 1.0 / factor;
1381	}
1382	Vector<Float> factors(tab.nrow(), factor);
1383	scaleByVector(tab, factors, False);
1384	} else {
1385
1386	// OK now we must deal with automatic look up of values.
1387	// We must also deal with the fact that the factors need
1388	// to be computed per IF and may be different and may
1389	// change per integration.
1390
1391	pushLog("Looking up conversion factors");
1392	convertBrightnessUnits(out, tokelvin, cfac);
1393	}
1394
1395	return out;
1396	}
1397
1398	void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
1399	bool tokelvin, float cfac )
1400	{
1401	Table& table = in->table();
1402	Instrument inst =
1403	STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
1404	TableIterator iter(table, "FREQ_ID");
1405	STFrequencies stfreqs = in->frequencies();
1406	STAttr sdAtt;
1407	while (!iter.pastEnd()) {
1408	Table tab = iter.table();
1409	ArrayColumn<Float> specCol(tab, "SPECTRA");
1410	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1411	ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
1412	MEpoch::ROScalarColumn timeCol(tab, "TIME");
1413
1414	uInt freqid; freqidCol.get(0, freqid);
1415	Vector<Float> tmpspec; specCol.get(0, tmpspec);
1416	// STAttr.JyPerK has a Vector interface... change sometime.
1417	Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
1418	for ( uInt i=0; i<tab.nrow(); ++i) {
1419	Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
1420	Float factor = cfac * jyperk;
1421	if ( tokelvin ) factor = Float(1.0) / factor;
1422	MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1423	ma *= factor;
1424	specCol.put(i, ma.getArray());
1425	flagCol.put(i, flagsFromMA(ma));
1426	}
1427	++iter;
1428	}
1429	}
1430
1431	CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
1432	float tau )
1433	{
1434	CountedPtr< Scantable > out = getScantable(in, false);
1435
1436	Table tab = out->table();
1437	ROScalarColumn<Float> elev(tab, "ELEVATION");
1438	ArrayColumn<Float> specCol(tab, "SPECTRA");
1439	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1440	for ( uInt i=0; i<tab.nrow(); ++i) {
1441	Float zdist = Float(C::pi_2) - elev(i);
1442	Float factor = exp(tau/cos(zdist));
1443	MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1444	ma *= factor;
1445	specCol.put(i, ma.getArray());
1446	flagCol.put(i, flagsFromMA(ma));
1447	}
1448	return out;
1449	}
1450
1451	CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
1452	const std::string& kernel,
1453	float width )
1454	{
1455	CountedPtr< Scantable > out = getScantable(in, false);
1456	Table& table = out->table();
1457	ArrayColumn<Float> specCol(table, "SPECTRA");
1458	ArrayColumn<uChar> flagCol(table, "FLAGTRA");
1459	Vector<Float> spec;
1460	Vector<uChar> flag;
1461	for ( uInt i=0; i<table.nrow(); ++i) {
1462	specCol.get(i, spec);
1463	flagCol.get(i, flag);
1464	Vector<Bool> mask(flag.nelements());
1465	convertArray(mask, flag);
1466	Vector<Float> specout;
1467	Vector<Bool> maskout;
1468	if ( kernel == "hanning" ) {
1469	mathutil::hanning(specout, maskout, spec , !mask);
1470	convertArray(flag, !maskout);
1471	} else if ( kernel == "rmedian" ) {
1472	mathutil::runningMedian(specout, maskout, spec , mask, width);
1473	convertArray(flag, maskout);
1474	}
1475	flagCol.put(i, flag);
1476	specCol.put(i, specout);
1477	}
1478	return out;
1479	}
1480
1481	CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
1482	const std::string& kernel, float width )
1483	{
1484	if (kernel == "rmedian" \|\| kernel == "hanning") {
1485	return smoothOther(in, kernel, width);
1486	}
1487	CountedPtr< Scantable > out = getScantable(in, false);
1488	Table& table = out->table();
1489	VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
1490	// same IFNO should have same no of channels
1491	// this saves overhead
1492	TableIterator iter(table, "IFNO");
1493	while (!iter.pastEnd()) {
1494	Table tab = iter.table();
1495	ArrayColumn<Float> specCol(tab, "SPECTRA");
1496	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1497	Vector<Float> tmpspec; specCol.get(0, tmpspec);
1498	uInt nchan = tmpspec.nelements();
1499	Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
1500	Convolver<Float> conv(kvec, IPosition(1,nchan));
1501	Vector<Float> spec;
1502	Vector<uChar> flag;
1503	for ( uInt i=0; i<tab.nrow(); ++i) {
1504	specCol.get(i, spec);
1505	flagCol.get(i, flag);
1506	Vector<Bool> mask(flag.nelements());
1507	convertArray(mask, flag);
1508	Vector<Float> specout;
1509	mathutil::replaceMaskByZero(specout, mask);
1510	conv.linearConv(specout, spec);
1511	specCol.put(i, specout);
1512	}
1513	++iter;
1514	}
1515	return out;
1516	}
1517
1518	CountedPtr< Scantable >
1519	STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
1520	{
1521	if ( in.size() < 2 ) {
1522	throw(AipsError("Need at least two scantables to perform a merge."));
1523	}
1524	std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
1525	bool insitu = insitu_;
1526	setInsitu(false);
1527	CountedPtr< Scantable > out = getScantable(*it, false);
1528	setInsitu(insitu);
1529	Table& tout = out->table();
1530	ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
1531	ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
1532	// Renumber SCANNO to be 0-based
1533	Vector<uInt> scannos = scannocol.getColumn();
1534	uInt offset = min(scannos);
1535	scannos -= offset;
1536	scannocol.putColumn(scannos);
1537	uInt newscanno = max(scannos)+1;
1538	++it;
1539	while ( it != in.end() ){
1540	if ( ! (it)->conformant(out) ) {
1541	// log message: "ignoring scantable i, as it isn't
1542	// conformant with the other(s)"
1543	cerr << "oh oh" << endl;
1544	++it;
1545	continue;
1546	}
1547	out->appendToHistoryTable((*it)->history());
1548	const Table& tab = (*it)->table();
1549	TableIterator scanit(tab, "SCANNO");
1550	while (!scanit.pastEnd()) {
1551	TableIterator freqit(scanit.table(), "FREQ_ID");
1552	while ( !freqit.pastEnd() ) {
1553	Table thetab = freqit.table();
1554	uInt nrow = tout.nrow();
1555	//tout.addRow(thetab.nrow());
1556	TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1557	ROTableRow row(thetab);
1558	for ( uInt i=0; i<thetab.nrow(); ++i) {
1559	uInt k = nrow+i;
1560	scannocol.put(k, newscanno);
1561	const TableRecord& rec = row.get(i);
1562	Double rv,rp,inc;
1563	(*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1564	uInt id;
1565	id = out->frequencies().addEntry(rp, rv, inc);
1566	freqidcol.put(k,id);
1567	String name,fname;Double rf;
1568	(*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1569	id = out->molecules().addEntry(rf, name, fname);
1570	molidcol.put(k, id);
1571	Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1572	(*it)->focus().getEntry(fax, ftan, frot, fhand,
1573	fmount,fuser, fxy, fxyp,
1574	rec.asuInt("FOCUS_ID"));
1575	id = out->focus().addEntry(fax, ftan, frot, fhand,
1576	fmount,fuser, fxy, fxyp);
1577	focusidcol.put(k, id);
1578	}
1579	++freqit;
1580	}
1581	++newscanno;
1582	++scanit;
1583	}
1584	++it;
1585	}
1586	return out;
1587	}
1588
1589	CountedPtr< Scantable >
1590	STMath::invertPhase( const CountedPtr < Scantable >& in )
1591	{
1592	return applyToPol(in, &STPol::invertPhase, Float(0.0));
1593	}
1594
1595	CountedPtr< Scantable >
1596	STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1597	{
1598	return applyToPol(in, &STPol::rotatePhase, Float(phase));
1599	}
1600
1601	CountedPtr< Scantable >
1602	STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1603	{
1604	return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1605	}
1606
1607	CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
1608	STPol::polOperation fptr,
1609	Float phase )
1610	{
1611	CountedPtr< Scantable > out = getScantable(in, false);
1612	Table& tout = out->table();
1613	Block<String> cols(4);
1614	cols[0] = String("SCANNO");
1615	cols[1] = String("BEAMNO");
1616	cols[2] = String("IFNO");
1617	cols[3] = String("CYCLENO");
1618	TableIterator iter(tout, cols);
1619	CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
1620	out->getPolType() );
1621	while (!iter.pastEnd()) {
1622	Table t = iter.table();
1623	ArrayColumn<Float> speccol(t, "SPECTRA");
1624	ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
1625	ScalarColumn<Float> parancol(t, "PARANGLE");
1626	Matrix<Float> pols(speccol.getColumn());
1627	try {
1628	stpol->setSpectra(pols);
1629	Float fang,fhand,parang;
1630	fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
1631	fhand = in->focusTable_.getFeedHand(focidcol(0));
1632	parang = parancol(0);
1633	/// @todo re-enable this
1634	// disable total feed angle to support paralactifying Caswell style
1635	stpol->setPhaseCorrections(parang, -parang, fhand);
1636	// use a member function pointer in STPol. This only works on
1637	// the STPol pointer itself, not the Counted Pointer so
1638	// derefernce it.
1639	(&((stpol))->fptr)(phase);
1640	speccol.putColumn(stpol->getSpectra());
1641	} catch (AipsError& e) {
1642	//delete stpol;stpol=0;
1643	throw(e);
1644	}
1645	++iter;
1646	}
1647	//delete stpol;stpol=0;
1648	return out;
1649	}
1650
1651	CountedPtr< Scantable >
1652	STMath::swapPolarisations( const CountedPtr< Scantable > & in )
1653	{
1654	CountedPtr< Scantable > out = getScantable(in, false);
1655	Table& tout = out->table();
1656	Table t0 = tout(tout.col("POLNO") == 0);
1657	Table t1 = tout(tout.col("POLNO") == 1);
1658	if ( t0.nrow() != t1.nrow() )
1659	throw(AipsError("Inconsistent number of polarisations"));
1660	ArrayColumn<Float> speccol0(t0, "SPECTRA");
1661	ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
1662	ArrayColumn<Float> speccol1(t1, "SPECTRA");
1663	ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
1664	Matrix<Float> s0 = speccol0.getColumn();
1665	Matrix<uChar> f0 = flagcol0.getColumn();
1666	speccol0.putColumn(speccol1.getColumn());
1667	flagcol0.putColumn(flagcol1.getColumn());
1668	speccol1.putColumn(s0);
1669	flagcol1.putColumn(f0);
1670	return out;
1671	}
1672
1673	CountedPtr< Scantable >
1674	STMath::averagePolarisations( const CountedPtr< Scantable > & in,
1675	const std::vector<bool>& mask,
1676	const std::string& weight )
1677	{
1678	if (in->npol() < 2 )
1679	throw(AipsError("averagePolarisations can only be applied to two or more"
1680	"polarisations"));
1681	bool insitu = insitu_;
1682	setInsitu(false);
1683	CountedPtr< Scantable > pols = getScantable(in, true);
1684	setInsitu(insitu);
1685	Table& tout = pols->table();
1686	std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
1687	Table tab = tableCommand(taql, in->table());
1688	if (tab.nrow() == 0 )
1689	throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
1690	TableCopy::copyRows(tout, tab);
1691	TableVector<uInt> vec(tout, "POLNO");
1692	vec = 0;
1693	pols->table_.rwKeywordSet().define("nPol", Int(1));
1694	//pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
1695	pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
1696	std::vector<CountedPtr<Scantable> > vpols;
1697	vpols.push_back(pols);
1698	CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
1699	return out;
1700	}
1701
1702	CountedPtr< Scantable >
1703	STMath::averageBeams( const CountedPtr< Scantable > & in,
1704	const std::vector<bool>& mask,
1705	const std::string& weight )
1706	{
1707	bool insitu = insitu_;
1708	setInsitu(false);
1709	CountedPtr< Scantable > beams = getScantable(in, false);
1710	setInsitu(insitu);
1711	Table& tout = beams->table();
1712	// give all rows the same BEAMNO
1713	TableVector<uInt> vec(tout, "BEAMNO");
1714	vec = 0;
1715	beams->table_.rwKeywordSet().define("nBeam", Int(1));
1716	std::vector<CountedPtr<Scantable> > vbeams;
1717	vbeams.push_back(beams);
1718	CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
1719	return out;
1720	}
1721
1722
1723	CountedPtr< Scantable >
1724	asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
1725	const std::string & refTime,
1726	const std::string & method)
1727	{
1728	// clone as this is not working insitu
1729	bool insitu = insitu_;
1730	setInsitu(false);
1731	CountedPtr< Scantable > out = getScantable(in, false);
1732	setInsitu(insitu);
1733	Table& tout = out->table();
1734	// Get reference Epoch to time of first row or given String
1735	Unit DAY(String("d"));
1736	MEpoch::Ref epochRef(in->getTimeReference());
1737	MEpoch refEpoch;
1738	if (refTime.length()>0) {
1739	Quantum<Double> qt;
1740	if (MVTime::read(qt,refTime)) {
1741	MVEpoch mv(qt);
1742	refEpoch = MEpoch(mv, epochRef);
1743	} else {
1744	throw(AipsError("Invalid format for Epoch string"));
1745	}
1746	} else {
1747	refEpoch = in->timeCol_(0);
1748	}
1749	MPosition refPos = in->getAntennaPosition();
1750
1751	InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1752	// test if user frame is different to base frame
1753	if ( in->frequencies().getFrameString(true)
1754	== in->frequencies().getFrameString(false) ) {
1755	throw(AipsError("Can't convert as no output frame has been set"
1756	" (use set_freqframe) or it is aligned already."));
1757	}
1758	MFrequency::Types system = in->frequencies().getFrame();
1759	MVTime mvt(refEpoch.getValue());
1760	String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1761	ostringstream oss;
1762	oss << "Aligned at reference Epoch " << epochout
1763	<< " in frame " << MFrequency::showType(system);
1764	pushLog(String(oss));
1765	// set up the iterator
1766	Block<String> cols(4);
1767	// select by constant direction
1768	cols[0] = String("SRCNAME");
1769	cols[1] = String("BEAMNO");
1770	// select by IF ( no of channels varies over this )
1771	cols[2] = String("IFNO");
1772	// select by restfrequency
1773	cols[3] = String("MOLECULE_ID");
1774	TableIterator iter(tout, cols);
1775	while ( !iter.pastEnd() ) {
1776	Table t = iter.table();
1777	MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1778	TableIterator fiter(t, "FREQ_ID");
1779	// determine nchan from the first row. This should work as
1780	// we are iterating over BEAMNO and IFNO // we should have constant direction
1781
1782	ROArrayColumn<Float> sCol(t, "SPECTRA");
1783	MDirection direction = dirCol(0);
1784	uInt nchan = sCol(0).nelements();
1785	while ( !fiter.pastEnd() ) {
1786	Table ftab = fiter.table();
1787	ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1788	// get the SpectralCoordinate for the freqid, which we are iterating over
1789	SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1790	FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1791	direction, refPos, system );
1792	// realign the SpectralCoordinate and put into the output Scantable
1793	Vector<String> units(1);
1794	units = String("Hz");
1795	Bool linear=True;
1796	SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1797	sc2.setWorldAxisUnits(units);
1798	uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1799	sc2.referenceValue()[0],
1800	sc2.increment()[0]);
1801	TableVector<uInt> tvec(ftab, "FREQ_ID");
1802	tvec = id;
1803	// create the "global" abcissa for alignment with same FREQ_ID
1804	Vector<Double> abc(nchan);
1805	Double w;
1806	for (uInt i=0; i<nchan; i++) {
1807	sC.toWorld(w,Double(i));
1808	abc[i] = w;
1809	}
1810	// cache abcissa for same time stamps, so iterate over those
1811	TableIterator timeiter(ftab, "TIME");
1812	while ( !timeiter.pastEnd() ) {
1813	Table tab = timeiter.table();
1814	ArrayColumn<Float> specCol(tab, "SPECTRA");
1815	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1816	MEpoch::ROScalarColumn timeCol(tab, "TIME");
1817	// use align abcissa cache after the first row
1818	bool first = true;
1819	// these rows should be just be POLNO
1820	for (int i=0; i<int(tab.nrow()); ++i) {
1821	// input values
1822	Vector<uChar> flag = flagCol(i);
1823	Vector<Bool> mask(flag.shape());
1824	Vector<Float> specOut, spec;
1825	spec = specCol(i);
1826	Vector<Bool> maskOut;Vector<uChar> flagOut;
1827	convertArray(mask, flag);
1828	// alignment
1829	Bool ok = fa.align(specOut, maskOut, abc, spec,
1830	mask, timeCol(i), !first,
1831	interp, False);
1832	// back into scantable
1833	flagOut.resize(maskOut.nelements());
1834	convertArray(flagOut, maskOut);
1835	flagCol.put(i, flagOut);
1836	specCol.put(i, specOut);
1837	// start abcissa caching
1838	first = false;
1839	}
1840	// next timestamp
1841	++timeiter;
1842	}
1843	// next FREQ_ID
1844	++fiter;
1845	}
1846	// next aligner
1847	++iter;
1848	}
1849	// set this afterwards to ensure we are doing insitu correctly.
1850	out->frequencies().setFrame(system, true);
1851	return out;
1852	}
1853
1854	CountedPtr<Scantable>
1855	asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
1856	const std::string & newtype )
1857	{
1858	if (in->npol() != 2 && in->npol() != 4)
1859	throw(AipsError("Can only convert two or four polarisations."));
1860	if ( in->getPolType() == newtype )
1861	throw(AipsError("No need to convert."));
1862	if ( ! in->selector_.empty() )
1863	throw(AipsError("Can only convert whole scantable. Unset the selection."));
1864	bool insitu = insitu_;
1865	setInsitu(false);
1866	CountedPtr< Scantable > out = getScantable(in, true);
1867	setInsitu(insitu);
1868	Table& tout = out->table();
1869	tout.rwKeywordSet().define("POLTYPE", String(newtype));
1870
1871	Block<String> cols(4);
1872	cols[0] = "SCANNO";
1873	cols[1] = "CYCLENO";
1874	cols[2] = "BEAMNO";
1875	cols[3] = "IFNO";
1876	TableIterator it(in->originalTable_, cols);
1877	String basetype = in->getPolType();
1878	STPol* stpol = STPol::getPolClass(in->factories_, basetype);
1879	try {
1880	while ( !it.pastEnd() ) {
1881	Table tab = it.table();
1882	uInt row = tab.rowNumbers()[0];
1883	stpol->setSpectra(in->getPolMatrix(row));
1884	Float fang,fhand,parang;
1885	fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
1886	fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
1887	parang = in->paraCol_(row);
1888	/// @todo re-enable this
1889	// disable total feed angle to support paralactifying Caswell style
1890	stpol->setPhaseCorrections(parang, -parang, fhand);
1891	Int npolout = 0;
1892	for (uInt i=0; i<tab.nrow(); ++i) {
1893	Vector<Float> outvec = stpol->getSpectrum(i, newtype);
1894	if ( outvec.nelements() > 0 ) {
1895	tout.addRow();
1896	TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
1897	ArrayColumn<Float> sCol(tout,"SPECTRA");
1898	ScalarColumn<uInt> pCol(tout,"POLNO");
1899	sCol.put(tout.nrow()-1 ,outvec);
1900	pCol.put(tout.nrow()-1 ,uInt(npolout));
1901	npolout++;
1902	}
1903	}
1904	tout.rwKeywordSet().define("nPol", npolout);
1905	++it;
1906	}
1907	} catch (AipsError& e) {
1908	delete stpol;
1909	throw(e);
1910	}
1911	delete stpol;
1912	return out;
1913	}
1914
1915	CountedPtr< Scantable >
1916	asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
1917	const std::string & scantype )
1918	{
1919	bool insitu = insitu_;
1920	setInsitu(false);
1921	CountedPtr< Scantable > out = getScantable(in, true);
1922	setInsitu(insitu);
1923	Table& tout = out->table();
1924	std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
1925	if (scantype == "on") {
1926	taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
1927	}
1928	Table tab = tableCommand(taql, in->table());
1929	TableCopy::copyRows(tout, tab);
1930	if (scantype == "on") {
1931	// re-index SCANNO to 0
1932	TableVector<uInt> vec(tout, "SCANNO");
1933	vec = 0;
1934	}
1935	return out;
1936	}
1937
1938	CountedPtr< Scantable >
1939	asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
1940	double frequency, double width )
1941	{
1942	CountedPtr< Scantable > out = getScantable(in, false);
1943	Table& tout = out->table();
1944	TableIterator iter(tout, "FREQ_ID");
1945	FFTServer<Float,Complex> ffts;
1946	while ( !iter.pastEnd() ) {
1947	Table tab = iter.table();
1948	Double rp,rv,inc;
1949	ROTableRow row(tab);
1950	const TableRecord& rec = row.get(0);
1951	uInt freqid = rec.asuInt("FREQ_ID");
1952	out->frequencies().getEntry(rp, rv, inc, freqid);
1953	ArrayColumn<Float> specCol(tab, "SPECTRA");
1954	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1955	for (int i=0; i<int(tab.nrow()); ++i) {
1956	Vector<Float> spec = specCol(i);
1957	Vector<uChar> flag = flagCol(i);
1958	Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
1959	Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
1960	for (int k=0; k < flag.nelements(); ++k ) {
1961	if (flag[k] > 0) {
1962	spec[k] = 0.0;
1963	}
1964	}
1965	Vector<Complex> lags;
1966	ffts.fft0(lags, spec);
1967	Int start = max(0, lag0);
1968	Int end = min(Int(lags.nelements()-1), lag1);
1969	if (start == end) {
1970	lags[start] = Complex(0.0);
1971	} else {
1972	for (int j=start; j <=end ;++j) {
1973	lags[j] = Complex(0.0);
1974	}
1975	}
1976	ffts.fft0(spec, lags);
1977	specCol.put(i, spec);
1978	}
1979	++iter;
1980	}
1981	return out;
1982	}

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

Download in other formats: