Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

source: branches/mergetest/src/STMath.cpp@ 1777

Visit:

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

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

Download in other formats: