Context Navigation

STMath.cpp@ 2526

Visit:

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

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

Context Navigation

source: branches/asap-3.x/src/STMath.cpp@ 2526

Download in other formats: