Context Navigation

source: trunk/src/STMath.cpp@ 1577

Visit:

Last change on this file since 1577 was 1571, checked in by Malte Marquarding, 15 years ago
Ticket #167: c++ part of running polynomial smoothing
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 65.9 KB

Line
1	//
2	// C++ Implementation: STMath
3	//
4	// Description:
5	//
6	//
7	// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8	//
9	// Copyright: See COPYING file that comes with this distribution
10	//
11	//
12
13	#include <casa/iomanip.h>
14	#include <casa/Exceptions/Error.h>
15	#include <casa/Containers/Block.h>
16	#include <casa/BasicSL/String.h>
17	#include <casa/Arrays/MaskArrLogi.h>
18	#include <casa/Arrays/MaskArrMath.h>
19	#include <casa/Arrays/ArrayLogical.h>
20	#include <casa/Arrays/ArrayMath.h>
21	#include <casa/Arrays/Slice.h>
22	#include <casa/Arrays/Slicer.h>
23	#include <casa/Containers/RecordField.h>
24	#include <tables/Tables/TableRow.h>
25	#include <tables/Tables/TableVector.h>
26	#include <tables/Tables/TabVecMath.h>
27	#include <tables/Tables/ExprNode.h>
28	#include <tables/Tables/TableRecord.h>
29	#include <tables/Tables/TableParse.h>
30	#include <tables/Tables/ReadAsciiTable.h>
31	#include <tables/Tables/TableIter.h>
32	#include <tables/Tables/TableCopy.h>
33	#include <scimath/Mathematics/FFTServer.h>
34
35	#include <lattices/Lattices/LatticeUtilities.h>
36
37	#include <coordinates/Coordinates/SpectralCoordinate.h>
38	#include <coordinates/Coordinates/CoordinateSystem.h>
39	#include <coordinates/Coordinates/CoordinateUtil.h>
40	#include <coordinates/Coordinates/FrequencyAligner.h>
41
42	#include <scimath/Mathematics/VectorKernel.h>
43	#include <scimath/Mathematics/Convolver.h>
44	#include <scimath/Functionals/Polynomial.h>
45
46	#include "MathUtils.h"
47	#include "RowAccumulator.h"
48	#include "STAttr.h"
49	#include "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);
1227	coeff = coeffs;
1228	msg = String("user");
1229	} else {
1230	STAttr sdAttr;
1231	coeff = sdAttr.gainElevationPoly(inst);
1232	ppoly = new Polynomial<Float>(3);
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	float tau )
1434	{
1435	CountedPtr< Scantable > out = getScantable(in, false);
1436
1437	Table tab = out->table();
1438	ROScalarColumn<Float> elev(tab, "ELEVATION");
1439	ArrayColumn<Float> specCol(tab, "SPECTRA");
1440	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1441	ArrayColumn<Float> tsysCol(tab, "TSYS");
1442	for ( uInt i=0; i<tab.nrow(); ++i) {
1443	Float zdist = Float(C::pi_2) - elev(i);
1444	Float factor = exp(tau/cos(zdist));
1445	MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1446	ma *= factor;
1447	specCol.put(i, ma.getArray());
1448	flagCol.put(i, flagsFromMA(ma));
1449	Vector<Float> tsys;
1450	tsysCol.get(i, tsys);
1451	tsys *= factor;
1452	tsysCol.put(i, tsys);
1453	}
1454	return out;
1455	}
1456
1457	CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
1458	const std::string& kernel,
1459	float width, int order)
1460	{
1461	CountedPtr< Scantable > out = getScantable(in, false);
1462	Table& table = out->table();
1463	ArrayColumn<Float> specCol(table, "SPECTRA");
1464	ArrayColumn<uChar> flagCol(table, "FLAGTRA");
1465	Vector<Float> spec;
1466	Vector<uChar> flag;
1467	for ( uInt i=0; i<table.nrow(); ++i) {
1468	specCol.get(i, spec);
1469	flagCol.get(i, flag);
1470	Vector<Bool> mask(flag.nelements());
1471	convertArray(mask, flag);
1472	Vector<Float> specout;
1473	Vector<Bool> maskout;
1474	if ( kernel == "hanning" ) {
1475	mathutil::hanning(specout, maskout, spec , !mask);
1476	convertArray(flag, !maskout);
1477	} else if ( kernel == "rmedian" ) {
1478	mathutil::runningMedian(specout, maskout, spec , mask, width);
1479	convertArray(flag, maskout);
1480	} else if ( kernel == "poly" ) {
1481	mathutil::polyfit(specout, maskout, spec, !mask, width, order);
1482	convertArray(flag, !maskout);
1483	}
1484	flagCol.put(i, flag);
1485	specCol.put(i, specout);
1486	}
1487	return out;
1488	}
1489
1490	CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
1491	const std::string& kernel, float width,
1492	int order)
1493	{
1494	if (kernel == "rmedian" \|\| kernel == "hanning" \|\| kernel == "poly") {
1495	return smoothOther(in, kernel, width, order);
1496	}
1497	CountedPtr< Scantable > out = getScantable(in, false);
1498	Table& table = out->table();
1499	VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
1500	// same IFNO should have same no of channels
1501	// this saves overhead
1502	TableIterator iter(table, "IFNO");
1503	while (!iter.pastEnd()) {
1504	Table tab = iter.table();
1505	ArrayColumn<Float> specCol(tab, "SPECTRA");
1506	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1507	Vector<Float> tmpspec; specCol.get(0, tmpspec);
1508	uInt nchan = tmpspec.nelements();
1509	Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
1510	Convolver<Float> conv(kvec, IPosition(1,nchan));
1511	Vector<Float> spec;
1512	Vector<uChar> flag;
1513	for ( uInt i=0; i<tab.nrow(); ++i) {
1514	specCol.get(i, spec);
1515	flagCol.get(i, flag);
1516	Vector<Bool> mask(flag.nelements());
1517	convertArray(mask, flag);
1518	Vector<Float> specout;
1519	mathutil::replaceMaskByZero(specout, mask);
1520	conv.linearConv(specout, spec);
1521	specCol.put(i, specout);
1522	}
1523	++iter;
1524	}
1525	return out;
1526	}
1527
1528	CountedPtr< Scantable >
1529	STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
1530	{
1531	if ( in.size() < 2 ) {
1532	throw(AipsError("Need at least two scantables to perform a merge."));
1533	}
1534	std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
1535	bool insitu = insitu_;
1536	setInsitu(false);
1537	CountedPtr< Scantable > out = getScantable(*it, false);
1538	setInsitu(insitu);
1539	Table& tout = out->table();
1540	ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
1541	ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
1542	// Renumber SCANNO to be 0-based
1543	Vector<uInt> scannos = scannocol.getColumn();
1544	uInt offset = min(scannos);
1545	scannos -= offset;
1546	scannocol.putColumn(scannos);
1547	uInt newscanno = max(scannos)+1;
1548	++it;
1549	while ( it != in.end() ){
1550	if ( ! (it)->conformant(out) ) {
1551	// non conformant.
1552	pushLog(String("Warning: Can't merge scantables as header info differs."));
1553	}
1554	out->appendToHistoryTable((*it)->history());
1555	const Table& tab = (*it)->table();
1556	TableIterator scanit(tab, "SCANNO");
1557	while (!scanit.pastEnd()) {
1558	TableIterator freqit(scanit.table(), "FREQ_ID");
1559	while ( !freqit.pastEnd() ) {
1560	Table thetab = freqit.table();
1561	uInt nrow = tout.nrow();
1562	tout.addRow(thetab.nrow());
1563	TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1564	ROTableRow row(thetab);
1565	for ( uInt i=0; i<thetab.nrow(); ++i) {
1566	uInt k = nrow+i;
1567	scannocol.put(k, newscanno);
1568	const TableRecord& rec = row.get(i);
1569	Double rv,rp,inc;
1570	(*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1571	uInt id;
1572	id = out->frequencies().addEntry(rp, rv, inc);
1573	freqidcol.put(k,id);
1574	String name,fname;Double rf;
1575	(*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1576	id = out->molecules().addEntry(rf, name, fname);
1577	molidcol.put(k, id);
1578	Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1579	(*it)->focus().getEntry(fax, ftan, frot, fhand,
1580	fmount,fuser, fxy, fxyp,
1581	rec.asuInt("FOCUS_ID"));
1582	id = out->focus().addEntry(fax, ftan, frot, fhand,
1583	fmount,fuser, fxy, fxyp);
1584	focusidcol.put(k, id);
1585	}
1586	++freqit;
1587	}
1588	++newscanno;
1589	++scanit;
1590	}
1591	++it;
1592	}
1593	return out;
1594	}
1595
1596	CountedPtr< Scantable >
1597	STMath::invertPhase( const CountedPtr < Scantable >& in )
1598	{
1599	return applyToPol(in, &STPol::invertPhase, Float(0.0));
1600	}
1601
1602	CountedPtr< Scantable >
1603	STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1604	{
1605	return applyToPol(in, &STPol::rotatePhase, Float(phase));
1606	}
1607
1608	CountedPtr< Scantable >
1609	STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1610	{
1611	return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1612	}
1613
1614	CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
1615	STPol::polOperation fptr,
1616	Float phase )
1617	{
1618	CountedPtr< Scantable > out = getScantable(in, false);
1619	Table& tout = out->table();
1620	Block<String> cols(4);
1621	cols[0] = String("SCANNO");
1622	cols[1] = String("BEAMNO");
1623	cols[2] = String("IFNO");
1624	cols[3] = String("CYCLENO");
1625	TableIterator iter(tout, cols);
1626	CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
1627	out->getPolType() );
1628	while (!iter.pastEnd()) {
1629	Table t = iter.table();
1630	ArrayColumn<Float> speccol(t, "SPECTRA");
1631	ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
1632	ScalarColumn<Float> parancol(t, "PARANGLE");
1633	Matrix<Float> pols(speccol.getColumn());
1634	try {
1635	stpol->setSpectra(pols);
1636	Float fang,fhand,parang;
1637	fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
1638	fhand = in->focusTable_.getFeedHand(focidcol(0));
1639	parang = parancol(0);
1640	/// @todo re-enable this
1641	// disable total feed angle to support paralactifying Caswell style
1642	stpol->setPhaseCorrections(parang, -parang, fhand);
1643	// use a member function pointer in STPol. This only works on
1644	// the STPol pointer itself, not the Counted Pointer so
1645	// derefernce it.
1646	(&((stpol))->fptr)(phase);
1647	speccol.putColumn(stpol->getSpectra());
1648	} catch (AipsError& e) {
1649	//delete stpol;stpol=0;
1650	throw(e);
1651	}
1652	++iter;
1653	}
1654	//delete stpol;stpol=0;
1655	return out;
1656	}
1657
1658	CountedPtr< Scantable >
1659	STMath::swapPolarisations( const CountedPtr< Scantable > & in )
1660	{
1661	CountedPtr< Scantable > out = getScantable(in, false);
1662	Table& tout = out->table();
1663	Table t0 = tout(tout.col("POLNO") == 0);
1664	Table t1 = tout(tout.col("POLNO") == 1);
1665	if ( t0.nrow() != t1.nrow() )
1666	throw(AipsError("Inconsistent number of polarisations"));
1667	ArrayColumn<Float> speccol0(t0, "SPECTRA");
1668	ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
1669	ArrayColumn<Float> speccol1(t1, "SPECTRA");
1670	ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
1671	Matrix<Float> s0 = speccol0.getColumn();
1672	Matrix<uChar> f0 = flagcol0.getColumn();
1673	speccol0.putColumn(speccol1.getColumn());
1674	flagcol0.putColumn(flagcol1.getColumn());
1675	speccol1.putColumn(s0);
1676	flagcol1.putColumn(f0);
1677	return out;
1678	}
1679
1680	CountedPtr< Scantable >
1681	STMath::averagePolarisations( const CountedPtr< Scantable > & in,
1682	const std::vector<bool>& mask,
1683	const std::string& weight )
1684	{
1685	if (in->npol() < 2 )
1686	throw(AipsError("averagePolarisations can only be applied to two or more"
1687	"polarisations"));
1688	bool insitu = insitu_;
1689	setInsitu(false);
1690	CountedPtr< Scantable > pols = getScantable(in, true);
1691	setInsitu(insitu);
1692	Table& tout = pols->table();
1693	std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
1694	Table tab = tableCommand(taql, in->table());
1695	if (tab.nrow() == 0 )
1696	throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
1697	TableCopy::copyRows(tout, tab);
1698	TableVector<uInt> vec(tout, "POLNO");
1699	vec = 0;
1700	pols->table_.rwKeywordSet().define("nPol", Int(1));
1701	//pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
1702	pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
1703	std::vector<CountedPtr<Scantable> > vpols;
1704	vpols.push_back(pols);
1705	CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
1706	return out;
1707	}
1708
1709	CountedPtr< Scantable >
1710	STMath::averageBeams( const CountedPtr< Scantable > & in,
1711	const std::vector<bool>& mask,
1712	const std::string& weight )
1713	{
1714	bool insitu = insitu_;
1715	setInsitu(false);
1716	CountedPtr< Scantable > beams = getScantable(in, false);
1717	setInsitu(insitu);
1718	Table& tout = beams->table();
1719	// give all rows the same BEAMNO
1720	TableVector<uInt> vec(tout, "BEAMNO");
1721	vec = 0;
1722	beams->table_.rwKeywordSet().define("nBeam", Int(1));
1723	std::vector<CountedPtr<Scantable> > vbeams;
1724	vbeams.push_back(beams);
1725	CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
1726	return out;
1727	}
1728
1729
1730	CountedPtr< Scantable >
1731	asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
1732	const std::string & refTime,
1733	const std::string & method)
1734	{
1735	// clone as this is not working insitu
1736	bool insitu = insitu_;
1737	setInsitu(false);
1738	CountedPtr< Scantable > out = getScantable(in, false);
1739	setInsitu(insitu);
1740	Table& tout = out->table();
1741	// Get reference Epoch to time of first row or given String
1742	Unit DAY(String("d"));
1743	MEpoch::Ref epochRef(in->getTimeReference());
1744	MEpoch refEpoch;
1745	if (refTime.length()>0) {
1746	Quantum<Double> qt;
1747	if (MVTime::read(qt,refTime)) {
1748	MVEpoch mv(qt);
1749	refEpoch = MEpoch(mv, epochRef);
1750	} else {
1751	throw(AipsError("Invalid format for Epoch string"));
1752	}
1753	} else {
1754	refEpoch = in->timeCol_(0);
1755	}
1756	MPosition refPos = in->getAntennaPosition();
1757
1758	InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1759	/*
1760	// Comment from MV.
1761	// the following code has been commented out because different FREQ_IDs have to be aligned together even
1762	// if the frame doesn't change. So far, lack of this check didn't cause any problems.
1763	// test if user frame is different to base frame
1764	if ( in->frequencies().getFrameString(true)
1765	== in->frequencies().getFrameString(false) ) {
1766	throw(AipsError("Can't convert as no output frame has been set"
1767	" (use set_freqframe) or it is aligned already."));
1768	}
1769	*/
1770	MFrequency::Types system = in->frequencies().getFrame();
1771	MVTime mvt(refEpoch.getValue());
1772	String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1773	ostringstream oss;
1774	oss << "Aligned at reference Epoch " << epochout
1775	<< " in frame " << MFrequency::showType(system);
1776	pushLog(String(oss));
1777	// set up the iterator
1778	Block<String> cols(4);
1779	// select by constant direction
1780	cols[0] = String("SRCNAME");
1781	cols[1] = String("BEAMNO");
1782	// select by IF ( no of channels varies over this )
1783	cols[2] = String("IFNO");
1784	// select by restfrequency
1785	cols[3] = String("MOLECULE_ID");
1786	TableIterator iter(tout, cols);
1787	while ( !iter.pastEnd() ) {
1788	Table t = iter.table();
1789	MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1790	TableIterator fiter(t, "FREQ_ID");
1791	// determine nchan from the first row. This should work as
1792	// we are iterating over BEAMNO and IFNO // we should have constant direction
1793
1794	ROArrayColumn<Float> sCol(t, "SPECTRA");
1795	const MDirection direction = dirCol(0);
1796	const uInt nchan = sCol(0).nelements();
1797
1798	// skip operations if there is nothing to align
1799	if (fiter.pastEnd()) {
1800	continue;
1801	}
1802
1803	Table ftab = fiter.table();
1804	// align all frequency ids with respect to the first encountered id
1805	ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1806	// get the SpectralCoordinate for the freqid, which we are iterating over
1807	SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1808	FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1809	direction, refPos, system );
1810	// realign the SpectralCoordinate and put into the output Scantable
1811	Vector<String> units(1);
1812	units = String("Hz");
1813	Bool linear=True;
1814	SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1815	sc2.setWorldAxisUnits(units);
1816	const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1817	sc2.referenceValue()[0],
1818	sc2.increment()[0]);
1819	while ( !fiter.pastEnd() ) {
1820	ftab = fiter.table();
1821	// spectral coordinate for the current FREQ_ID
1822	ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
1823	sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
1824	// create the "global" abcissa for alignment with same FREQ_ID
1825	Vector<Double> abc(nchan);
1826	for (uInt i=0; i<nchan; i++) {
1827	Double w;
1828	sC.toWorld(w,Double(i));
1829	abc[i] = w;
1830	}
1831	TableVector<uInt> tvec(ftab, "FREQ_ID");
1832	// assign new frequency id to all rows
1833	tvec = id;
1834	// cache abcissa for same time stamps, so iterate over those
1835	TableIterator timeiter(ftab, "TIME");
1836	while ( !timeiter.pastEnd() ) {
1837	Table tab = timeiter.table();
1838	ArrayColumn<Float> specCol(tab, "SPECTRA");
1839	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1840	MEpoch::ROScalarColumn timeCol(tab, "TIME");
1841	// use align abcissa cache after the first row
1842	// these rows should be just be POLNO
1843	bool first = true;
1844	for (int i=0; i<int(tab.nrow()); ++i) {
1845	// input values
1846	Vector<uChar> flag = flagCol(i);
1847	Vector<Bool> mask(flag.shape());
1848	Vector<Float> specOut, spec;
1849	spec = specCol(i);
1850	Vector<Bool> maskOut;Vector<uChar> flagOut;
1851	convertArray(mask, flag);
1852	// alignment
1853	Bool ok = fa.align(specOut, maskOut, abc, spec,
1854	mask, timeCol(i), !first,
1855	interp, False);
1856	// back into scantable
1857	flagOut.resize(maskOut.nelements());
1858	convertArray(flagOut, maskOut);
1859	flagCol.put(i, flagOut);
1860	specCol.put(i, specOut);
1861	// start abcissa caching
1862	first = false;
1863	}
1864	// next timestamp
1865	++timeiter;
1866	}
1867	// next FREQ_ID
1868	++fiter;
1869	}
1870	// next aligner
1871	++iter;
1872	}
1873	// set this afterwards to ensure we are doing insitu correctly.
1874	out->frequencies().setFrame(system, true);
1875	return out;
1876	}
1877
1878	CountedPtr<Scantable>
1879	asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
1880	const std::string & newtype )
1881	{
1882	if (in->npol() != 2 && in->npol() != 4)
1883	throw(AipsError("Can only convert two or four polarisations."));
1884	if ( in->getPolType() == newtype )
1885	throw(AipsError("No need to convert."));
1886	if ( ! in->selector_.empty() )
1887	throw(AipsError("Can only convert whole scantable. Unset the selection."));
1888	bool insitu = insitu_;
1889	setInsitu(false);
1890	CountedPtr< Scantable > out = getScantable(in, true);
1891	setInsitu(insitu);
1892	Table& tout = out->table();
1893	tout.rwKeywordSet().define("POLTYPE", String(newtype));
1894
1895	Block<String> cols(4);
1896	cols[0] = "SCANNO";
1897	cols[1] = "CYCLENO";
1898	cols[2] = "BEAMNO";
1899	cols[3] = "IFNO";
1900	TableIterator it(in->originalTable_, cols);
1901	String basetype = in->getPolType();
1902	STPol* stpol = STPol::getPolClass(in->factories_, basetype);
1903	try {
1904	while ( !it.pastEnd() ) {
1905	Table tab = it.table();
1906	uInt row = tab.rowNumbers()[0];
1907	stpol->setSpectra(in->getPolMatrix(row));
1908	Float fang,fhand,parang;
1909	fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
1910	fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
1911	parang = in->paraCol_(row);
1912	/// @todo re-enable this
1913	// disable total feed angle to support paralactifying Caswell style
1914	stpol->setPhaseCorrections(parang, -parang, fhand);
1915	Int npolout = 0;
1916	for (uInt i=0; i<tab.nrow(); ++i) {
1917	Vector<Float> outvec = stpol->getSpectrum(i, newtype);
1918	if ( outvec.nelements() > 0 ) {
1919	tout.addRow();
1920	TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
1921	ArrayColumn<Float> sCol(tout,"SPECTRA");
1922	ScalarColumn<uInt> pCol(tout,"POLNO");
1923	sCol.put(tout.nrow()-1 ,outvec);
1924	pCol.put(tout.nrow()-1 ,uInt(npolout));
1925	npolout++;
1926	}
1927	}
1928	tout.rwKeywordSet().define("nPol", npolout);
1929	++it;
1930	}
1931	} catch (AipsError& e) {
1932	delete stpol;
1933	throw(e);
1934	}
1935	delete stpol;
1936	return out;
1937	}
1938
1939	CountedPtr< Scantable >
1940	asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
1941	const std::string & scantype )
1942	{
1943	bool insitu = insitu_;
1944	setInsitu(false);
1945	CountedPtr< Scantable > out = getScantable(in, true);
1946	setInsitu(insitu);
1947	Table& tout = out->table();
1948	std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
1949	if (scantype == "on") {
1950	taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
1951	}
1952	Table tab = tableCommand(taql, in->table());
1953	TableCopy::copyRows(tout, tab);
1954	if (scantype == "on") {
1955	// re-index SCANNO to 0
1956	TableVector<uInt> vec(tout, "SCANNO");
1957	vec = 0;
1958	}
1959	return out;
1960	}
1961
1962	CountedPtr< Scantable >
1963	asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
1964	double frequency, double width )
1965	{
1966	CountedPtr< Scantable > out = getScantable(in, false);
1967	Table& tout = out->table();
1968	TableIterator iter(tout, "FREQ_ID");
1969	FFTServer<Float,Complex> ffts;
1970	while ( !iter.pastEnd() ) {
1971	Table tab = iter.table();
1972	Double rp,rv,inc;
1973	ROTableRow row(tab);
1974	const TableRecord& rec = row.get(0);
1975	uInt freqid = rec.asuInt("FREQ_ID");
1976	out->frequencies().getEntry(rp, rv, inc, freqid);
1977	ArrayColumn<Float> specCol(tab, "SPECTRA");
1978	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1979	for (int i=0; i<int(tab.nrow()); ++i) {
1980	Vector<Float> spec = specCol(i);
1981	Vector<uChar> flag = flagCol(i);
1982	Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
1983	Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
1984	for (int k=0; k < flag.nelements(); ++k ) {
1985	if (flag[k] > 0) {
1986	spec[k] = 0.0;
1987	}
1988	}
1989	Vector<Complex> lags;
1990	ffts.fft0(lags, spec);
1991	Int start = max(0, lag0);
1992	Int end = min(Int(lags.nelements()-1), lag1);
1993	if (start == end) {
1994	lags[start] = Complex(0.0);
1995	} else {
1996	for (int j=start; j <=end ;++j) {
1997	lags[j] = Complex(0.0);
1998	}
1999	}
2000	ffts.fft0(spec, lags);
2001	specCol.put(i, spec);
2002	}
2003	++iter;
2004	}
2005	return out;
2006	}

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

Download in other formats: