Context Navigation

source: trunk/src/STMath.cpp@ 1480

Visit:

Last change on this file since 1480 was 1475, checked in by Max Voronkov, 16 years ago
fixed frequency alignment problem (#142). Different FREQ_IDs are now gridded together. \ Different IFs are left intact therefore causing no problem for the Tid case. Tested on Tid and Mopra data successfully.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 65.5 KB

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

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

Download in other formats: