Context Navigation

STMath.cpp @ 2959

Visit:

Last change on this file since 2959 was 2959, checked in by Takeshi Nakazato, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6582

Ready for Test: Yes

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix for handling of data that one scan is fully flagged but other
scans are valid.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 147.4 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 <sstream>
14
15	#include <casa/iomanip.h>
16	#include <casa/Arrays/MaskArrLogi.h>
17	#include <casa/Arrays/MaskArrMath.h>
18	#include <casa/Arrays/ArrayLogical.h>
19	#include <casa/Arrays/ArrayMath.h>
20	#include <casa/Arrays/Slice.h>
21	#include <casa/Arrays/Slicer.h>
22	#include <casa/BasicSL/String.h>
23	#include <casa/Containers/Block.h>
24	#include <casa/Containers/RecordField.h>
25	#include <casa/Exceptions/Error.h>
26	#include <casa/Logging/LogIO.h>
27	#include <casa/Quanta/Quantum.h>
28
29	#include <coordinates/Coordinates/CoordinateSystem.h>
30	#include <coordinates/Coordinates/CoordinateUtil.h>
31	#include <coordinates/Coordinates/FrequencyAligner.h>
32	#include <coordinates/Coordinates/SpectralCoordinate.h>
33
34	#include <lattices/Lattices/LatticeUtilities.h>
35
36	#include <scimath/Functionals/Polynomial.h>
37	#include <scimath/Mathematics/Convolver.h>
38	#include <scimath/Mathematics/VectorKernel.h>
39
40	#include <tables/Tables/ExprNode.h>
41	#include <tables/Tables/ReadAsciiTable.h>
42	#include <tables/Tables/TableCopy.h>
43	#include <tables/Tables/TableIter.h>
44	#include <tables/Tables/TableParse.h>
45	#include <tables/Tables/TableRecord.h>
46	#include <tables/Tables/TableRow.h>
47	#include <tables/Tables/TableVector.h>
48	#include <tables/Tables/TabVecMath.h>
49
50	#include <atnf/PKSIO/SrcType.h>
51
52	#include "RowAccumulator.h"
53	#include "STAttr.h"
54	#include "STMath.h"
55	#include "STSelector.h"
56	#include "Accelerator.h"
57	#include "STIdxIter.h"
58
59	#include "CalibrationHelper.h"
60
61	using namespace casa;
62	using namespace asap;
63
64	// 2012/02/17 TN
65	// Since STGrid is implemented, average doesn't consider direction
66	// when accumulating
67	// tolerance for direction comparison (rad)
68	// #define TOL_OTF 1.0e-15
69	// #define TOL_POINT 2.9088821e-4 // 1 arcmin
70
71	STMath::STMath(bool insitu) :
72	insitu_(insitu)
73	{
74	}
75
76
77	STMath::~STMath()
78	{
79	}
80
81	CountedPtr<Scantable>
82	STMath::average( const std::vector<CountedPtr<Scantable> >& in,
83	const std::vector<bool>& mask,
84	const std::string& weight,
85	const std::string& avmode)
86	{
87	// double t0, t1 ;
88	// t0 = mathutil::gettimeofday_sec() ;
89
90	LogIO os( LogOrigin( "STMath", "average()", WHERE ) ) ;
91	if ( avmode == "SCAN" && in.size() != 1 )
92	throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
93	"Use merge first."));
94	WeightType wtype = stringToWeight(weight);
95
96	// 2012/02/17 TN
97	// Since STGrid is implemented, average doesn't consider direction
98	// when accumulating
99	// check if OTF observation
100	// String obstype = in[0]->getHeader().obstype ;
101	// Double tol = 0.0 ;
102	// if ( (obstype.find( "OTF" ) != String::npos) \|\| (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
103	// tol = TOL_OTF ;
104	// }
105	// else {
106	// tol = TOL_POINT ;
107	// }
108
109	// output
110	// clone as this is non insitu
111	bool insitu = insitu_;
112	setInsitu(false);
113	CountedPtr< Scantable > out = getScantable(in[0], true);
114	setInsitu(insitu);
115	std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
116	++stit;
117	while ( stit != in.end() ) {
118	out->appendToHistoryTable((*stit)->history());
119	++stit;
120	}
121
122	Table& tout = out->table();
123
124	/// @todo check if all scantables are conformant
125
126	ArrayColumn<Float> specColOut(tout,"SPECTRA");
127	ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
128	ArrayColumn<Float> tsysColOut(tout,"TSYS");
129	ScalarColumn<Double> mjdColOut(tout,"TIME");
130	ScalarColumn<Double> intColOut(tout,"INTERVAL");
131	ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
132	ScalarColumn<uInt> scanColOut(tout,"SCANNO");
133	ScalarColumn<uInt> flagRowColOut(tout,"FLAGROW");
134
135	// set up the output table rows. These are based on the structure of the
136	// FIRST scantable in the vector
137	const Table& baset = in[0]->table();
138
139	RowAccumulator acc(wtype);
140	Vector<Bool> cmask(mask);
141	acc.setUserMask(cmask);
142	// ROTableRow row(tout);
143	ROArrayColumn<Float> specCol, tsysCol;
144	ROArrayColumn<uChar> flagCol;
145	ROScalarColumn<Double> mjdCol, intCol;
146	ROScalarColumn<Int> scanIDCol;
147	ROScalarColumn<uInt> flagRowCol;
148
149	//Vector<uInt> rowstodelete;
150	Block<uInt> rowstodelB( in[0]->nrow() ) ;
151	uInt nrowdel = 0 ;
152
153	// Block<String> cols(3);
154	vector<string> cols(3) ;
155	cols[0] = String("BEAMNO");
156	cols[1] = String("IFNO");
157	cols[2] = String("POLNO");
158	if ( avmode == "SOURCE" ) {
159	cols.resize(4);
160	cols[3] = String("SRCNAME");
161	}
162	if ( avmode == "SCAN" && in.size() == 1) {
163	//cols.resize(4);
164	//cols[3] = String("SCANNO");
165	cols.resize(5);
166	cols[3] = String("SRCNAME");
167	cols[4] = String("SCANNO");
168	}
169	uInt outrowCount = 0;
170	// use STIdxIter2 instead of TableIterator
171	STIdxIter2 iter( in[0], cols ) ;
172	// double t2 = 0 ;
173	// double t3 = 0 ;
174	// double t4 = 0 ;
175	// double t5 = 0 ;
176	// TableIterator iter(baset, cols);
177	// int count = 0 ;
178	while (!iter.pastEnd()) {
179	Vector<uInt> rows = iter.getRows( SHARE ) ;
180	if ( rows.nelements() == 0 ) {
181	iter.next() ;
182	continue ;
183	}
184	Record current = iter.currentValue() ;
185	//Table subt = iter.table();
186	// copy the first row of this selection into the new table
187	tout.addRow();
188	// t4 = mathutil::gettimeofday_sec() ;
189	// skip to copy SPECTRA, FLAGTRA, and TSYS since those heavy columns are
190	// overwritten in the following process
191	copyRows( tout, baset, outrowCount, rows[0], 1, False, False, False ) ;
192	// t5 += mathutil::gettimeofday_sec() - t4 ;
193	// re-index to 0
194	if ( avmode != "SCAN" && avmode != "SOURCE" ) {
195	scanColOut.put(outrowCount, uInt(0));
196	}
197
198	// 2012/02/17 TN
199	// Since STGrid is implemented, average doesn't consider direction
200	// when accumulating
201	// MDirection::ScalarColumn dircol ;
202	// dircol.attach( subt, "DIRECTION" ) ;
203	// Int length = subt.nrow() ;
204	// vector< Vector<Double> > dirs ;
205	// vector<int> indexes ;
206	// for ( Int i = 0 ; i < length ; i++ ) {
207	// Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
208	// //os << << count++ << ": " ;
209	// //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
210	// bool adddir = true ;
211	// for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
212	// //if ( allTrue( t == dirs[j] ) ) {
213	// Double dx = t[0] - dirs[j][0] ;
214	// Double dy = t[1] - dirs[j][1] ;
215	// Double dd = sqrt( dx * dx + dy * dy ) ;
216	// //if ( allNearAbs( t, dirs[j], tol ) ) {
217	// if ( dd <= tol ) {
218	// adddir = false ;
219	// break ;
220	// }
221	// }
222	// if ( adddir ) {
223	// dirs.push_back( t ) ;
224	// indexes.push_back( i ) ;
225	// }
226	// }
227	// uInt rowNum = dirs.size() ;
228	// tout.addRow( rowNum ) ;
229	// for ( uInt i = 0 ; i < rowNum ; i++ ) {
230	// TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
231	// // re-index to 0
232	// if ( avmode != "SCAN" && avmode != "SOURCE" ) {
233	// scanColOut.put(outrowCount+i, uInt(0));
234	// }
235	// }
236	// outrowCount += rowNum ;
237
238	// merge loop
239	uInt i = outrowCount ;
240	// in[0] is already selected by iterator
241	specCol.attach(baset,"SPECTRA");
242	flagCol.attach(baset,"FLAGTRA");
243	tsysCol.attach(baset,"TSYS");
244	intCol.attach(baset,"INTERVAL");
245	mjdCol.attach(baset,"TIME");
246	flagRowCol.attach(baset,"FLAGROW");
247	Vector<Float> spec,tsys;
248	Vector<uChar> flag;
249	Double inter,time;
250	uInt flagRow;
251
252	for (uInt l = 0; l < rows.nelements(); ++l ) {
253	uInt k = rows[l] ;
254	flagCol.get(k, flag);
255	Vector<Bool> bflag(flag.shape());
256	flagRowCol.get(k, flagRow);
257	if (flagRow > 0)
258	bflag = true;
259	else
260	convertArray(bflag, flag);
261	/*
262	if ( allEQ(bflag, True) ) {
263	continue;//don't accumulate
264	}
265	*/
266	specCol.get(k, spec);
267	tsysCol.get(k, tsys);
268	intCol.get(k, inter);
269	mjdCol.get(k, time);
270	// spectrum has to be added last to enable weighting by the other values
271	// t2 = mathutil::gettimeofday_sec() ;
272	acc.add(spec, !bflag, tsys, inter, time);
273	// t3 += mathutil::gettimeofday_sec() - t2 ;
274
275	}
276
277
278	// in[0] is already selected by TableIterator so that index is
279	// started from 1
280	for ( int j=1; j < int(in.size()); ++j ) {
281	const Table& tin = in[j]->table();
282	//const TableRecord& rec = row.get(i);
283	ROScalarColumn<Double> tmp(tin, "TIME");
284	Double td;tmp.get(0,td);
285
286	#if 1
287	static char const*const colNames1[] = { "IFNO", "BEAMNO", "POLNO" };
288	//uInt const values1[] = { rec.asuInt("IFNO"), rec.asuInt("BEAMNO"), rec.asuInt("POLNO") };
289	uInt const values1[] = { current.asuInt("IFNO"), current.asuInt("BEAMNO"), current.asuInt("POLNO") };
290	SingleTypeEqPredicate<uInt, 3> myPred(tin, colNames1, values1);
291	CustomTableExprNodeRep myNodeRep(tin, myPred);
292	myNodeRep.link(); // to avoid automatic delete when myExpr is destructed.
293	CustomTableExprNode myExpr(myNodeRep);
294	Table basesubt = tin(myExpr);
295	#else
296	// Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
297	// && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
298	// && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
299	Table basesubt = tin( tin.col("BEAMNO") == current.asuInt("BEAMNO")
300	&& tin.col("IFNO") == current.asuInt("IFNO")
301	&& tin.col("POLNO") == current.asuInt("POLNO") );
302	#endif
303	Table subt;
304	if ( avmode == "SOURCE") {
305	// subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME"));
306	subt = basesubt( basesubt.col("SRCNAME") == current.asString("SRCNAME") );
307
308	} else if (avmode == "SCAN") {
309	// subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME")
310	// && basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
311	subt = basesubt( basesubt.col("SRCNAME") == current.asString("SRCNAME")
312	&& basesubt.col("SCANNO") == current.asuInt("SCANNO") );
313	} else {
314	subt = basesubt;
315	}
316
317	// 2012/02/17 TN
318	// Since STGrid is implemented, average doesn't consider direction
319	// when accumulating
320	// vector<uInt> removeRows ;
321	// uInt nrsubt = subt.nrow() ;
322	// for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
323	// //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
324	// Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
325	// Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
326	// double dx = x0[0] - x1[0];
327	// double dy = x0[1] - x1[1];
328	// Double dd = sqrt( dx * dx + dy * dy ) ;
329	// //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
330	// if ( dd > tol ) {
331	// removeRows.push_back( irow ) ;
332	// }
333	// }
334	// if ( removeRows.size() != 0 ) {
335	// subt.removeRow( removeRows ) ;
336	// }
337
338	// if ( nrsubt == removeRows.size() )
339	// throw(AipsError("Averaging data is empty.")) ;
340
341	specCol.attach(subt,"SPECTRA");
342	flagCol.attach(subt,"FLAGTRA");
343	tsysCol.attach(subt,"TSYS");
344	intCol.attach(subt,"INTERVAL");
345	mjdCol.attach(subt,"TIME");
346	flagRowCol.attach(subt,"FLAGROW");
347	for (uInt k = 0; k < subt.nrow(); ++k ) {
348	flagCol.get(k, flag);
349	Vector<Bool> bflag(flag.shape());
350	flagRowCol.get(k, flagRow);
351	if (flagRow > 0)
352	bflag = true;
353	else
354	convertArray(bflag, flag);
355	/*
356	if ( allEQ(bflag, True) ) {
357	continue;//don't accumulate
358	}
359	*/
360	specCol.get(k, spec);
361	//tsysCol.get(k, tsys);
362	tsys.assign( tsysCol(k) );
363	intCol.get(k, inter);
364	mjdCol.get(k, time);
365	// spectrum has to be added last to enable weighting by the other values
366	// t2 = mathutil::gettimeofday_sec() ;
367	acc.add(spec, !bflag, tsys, inter, time);
368	// t3 += mathutil::gettimeofday_sec() - t2 ;
369	}
370
371	}
372	const Vector<Bool>& msk = acc.getMask();
373	if ( allEQ(msk, False) ) {
374	rowstodelB[nrowdel] = i ;
375	nrowdel++ ;
376	outrowCount++;
377	acc.reset();
378	iter.next();
379	continue;
380	}
381	//write out
382	if (acc.state()) {
383	// If there exists a channel at which all the input spectra are masked,
384	// spec has 'nan' values for that channel and it may affect the following
385	// processes. To avoid this, replacing 'nan' values in spec with
386	// weighted-mean of all spectra in the following line.
387	// (done for CAS-2776, 2011/04/07 by Wataru Kawasaki)
388	acc.replaceNaN();
389
390	Vector<uChar> flg(msk.shape());
391	convertArray(flg, !msk);
392	for (uInt k = 0; k < flg.nelements(); ++k) {
393	uChar userFlag = 1 << 7;
394	if (msk[k]==True) userFlag = 0 << 7;
395	flg(k) = userFlag;
396	}
397
398	flagColOut.put(i, flg);
399	specColOut.put(i, acc.getSpectrum());
400	tsysColOut.put(i, acc.getTsys());
401	intColOut.put(i, acc.getInterval());
402	mjdColOut.put(i, acc.getTime());
403	// we should only have one cycle now -> reset it to be 0
404	// frequency switched data has different CYCLENO for different IFNO
405	// which requires resetting this value
406	cycColOut.put(i, uInt(0));
407	// completely flagged rows are removed anyway
408	flagRowColOut.put(i, uInt(0));
409	} else {
410	os << "For output row="<<i<<", all input rows of data are flagged. no averaging" << LogIO::POST;
411	}
412	acc.reset();
413
414	// merge with while loop for preparing out table
415	++outrowCount;
416	// ++iter ;
417	iter.next() ;
418	}
419
420	if ( nrowdel > 0 ) {
421	Vector<uInt> rowstodelete( IPosition(1,nrowdel), rowstodelB.storage(), SHARE ) ;
422	//os << rowstodelete << LogIO::POST ;
423	tout.removeRow(rowstodelete);
424	if (tout.nrow() == 0) {
425	throw(AipsError("Can't average fully flagged data."));
426	}
427	}
428
429	// t1 = mathutil::gettimeofday_sec() ;
430	// cout << "elapsed time for average(): " << t1-t0 << " sec" << endl ;
431	// cout << " elapsed time for acc.add(): " << t3 << " sec" << endl ;
432	// cout << " elapsed time for copyRows(): " << t5 << " sec" << endl ;
433
434	return out;
435	}
436
437	CountedPtr< Scantable >
438	STMath::averageChannel( const CountedPtr < Scantable > & in,
439	const std::string & mode,
440	const std::string& avmode )
441	{
442	(void) mode; // currently unused
443	// 2012/02/17 TN
444	// Since STGrid is implemented, average doesn't consider direction
445	// when accumulating
446	// check if OTF observation
447	// String obstype = in->getHeader().obstype ;
448	// Double tol = 0.0 ;
449	// if ( obstype.find( "OTF" ) != String::npos ) {
450	// tol = TOL_OTF ;
451	// }
452	// else {
453	// tol = TOL_POINT ;
454	// }
455
456	// clone as this is non insitu
457	bool insitu = insitu_;
458	setInsitu(false);
459	CountedPtr< Scantable > out = getScantable(in, true);
460	setInsitu(insitu);
461	Table& tout = out->table();
462	ArrayColumn<Float> specColOut(tout,"SPECTRA");
463	ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
464	ArrayColumn<Float> tsysColOut(tout,"TSYS");
465	ScalarColumn<uInt> scanColOut(tout,"SCANNO");
466	ScalarColumn<Double> intColOut(tout, "INTERVAL");
467	Table tmp = in->table().sort("BEAMNO");
468	Block<String> cols(3);
469	cols[0] = String("BEAMNO");
470	cols[1] = String("IFNO");
471	cols[2] = String("POLNO");
472	if ( avmode == "SCAN") {
473	cols.resize(4);
474	cols[3] = String("SCANNO");
475	}
476	uInt outrowCount = 0;
477	uChar userflag = 1 << 7;
478	TableIterator iter(tmp, cols);
479	while (!iter.pastEnd()) {
480	Table subt = iter.table();
481	ROArrayColumn<Float> specCol, tsysCol;
482	ROArrayColumn<uChar> flagCol;
483	ROScalarColumn<Double> intCol(subt, "INTERVAL");
484	specCol.attach(subt,"SPECTRA");
485	flagCol.attach(subt,"FLAGTRA");
486	tsysCol.attach(subt,"TSYS");
487
488	tout.addRow();
489	TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
490	if ( avmode != "SCAN") {
491	scanColOut.put(outrowCount, uInt(0));
492	}
493	Vector<Float> tmp;
494	specCol.get(0, tmp);
495	uInt nchan = tmp.nelements();
496	// have to do channel by channel here as MaskedArrMath
497	// doesn't have partialMedians
498	Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
499	Vector<Float> outspec(nchan);
500	Vector<uChar> outflag(nchan,0);
501	Vector<Float> outtsys(1);/// @fixme when tsys is channel based
502	for (uInt i=0; i<nchan; ++i) {
503	Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
504	MaskedArray<Float> ma = maskedArray(specs,flags);
505	outspec[i] = median(ma);
506	if ( allEQ(ma.getMask(), False) )
507	outflag[i] = userflag;// flag data
508	}
509	outtsys[0] = median(tsysCol.getColumn());
510	specColOut.put(outrowCount, outspec);
511	flagColOut.put(outrowCount, outflag);
512	tsysColOut.put(outrowCount, outtsys);
513	Double intsum = sum(intCol.getColumn());
514	intColOut.put(outrowCount, intsum);
515	++outrowCount;
516	++iter;
517
518	// 2012/02/17 TN
519	// Since STGrid is implemented, average doesn't consider direction
520	// when accumulating
521	// MDirection::ScalarColumn dircol ;
522	// dircol.attach( subt, "DIRECTION" ) ;
523	// Int length = subt.nrow() ;
524	// vector< Vector<Double> > dirs ;
525	// vector<int> indexes ;
526	// // Handle MX mode averaging
527	// if (in->nbeam() > 1 ) {
528	// length = 1;
529	// }
530	// for ( Int i = 0 ; i < length ; i++ ) {
531	// Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
532	// bool adddir = true ;
533	// for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
534	// //if ( allTrue( t == dirs[j] ) ) {
535	// Double dx = t[0] - dirs[j][0] ;
536	// Double dy = t[1] - dirs[j][1] ;
537	// Double dd = sqrt( dx * dx + dy * dy ) ;
538	// //if ( allNearAbs( t, dirs[j], tol ) ) {
539	// if ( dd <= tol ) {
540	// adddir = false ;
541	// break ;
542	// }
543	// }
544	// if ( adddir ) {
545	// dirs.push_back( t ) ;
546	// indexes.push_back( i ) ;
547	// }
548	// }
549	// uInt rowNum = dirs.size() ;
550	// tout.addRow( rowNum );
551	// for ( uInt i = 0 ; i < rowNum ; i++ ) {
552	// TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
553	// // Handle MX mode averaging
554	// if ( avmode != "SCAN") {
555	// scanColOut.put(outrowCount+i, uInt(0));
556	// }
557	// }
558	// MDirection::ScalarColumn dircolOut ;
559	// dircolOut.attach( tout, "DIRECTION" ) ;
560	// for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
561	// Vector<Double> t = \
562	// dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
563	// Vector<Float> tmp;
564	// specCol.get(0, tmp);
565	// uInt nchan = tmp.nelements();
566	// // have to do channel by channel here as MaskedArrMath
567	// // doesn't have partialMedians
568	// Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
569	// // mask spectra for different DIRECTION
570	// for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
571	// Vector<Double> direction = \
572	// dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
573	// //if ( t[0] != direction[0] \|\| t[1] != direction[1] ) {
574	// Double dx = t[0] - direction[0];
575	// Double dy = t[1] - direction[1];
576	// Double dd = sqrt(dxdx + dydy);
577	// //if ( !allNearAbs( t, direction, tol ) ) {
578	// if ( dd > tol && in->nbeam() < 2 ) {
579	// flags[jrow] = userflag ;
580	// }
581	// }
582	// Vector<Float> outspec(nchan);
583	// Vector<uChar> outflag(nchan,0);
584	// Vector<Float> outtsys(1);/// @fixme when tsys is channel based
585	// for (uInt i=0; i<nchan; ++i) {
586	// Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
587	// MaskedArray<Float> ma = maskedArray(specs,flags);
588	// outspec[i] = median(ma);
589	// if ( allEQ(ma.getMask(), False) )
590	// outflag[i] = userflag;// flag data
591	// }
592	// outtsys[0] = median(tsysCol.getColumn());
593	// specColOut.put(outrowCount+irow, outspec);
594	// flagColOut.put(outrowCount+irow, outflag);
595	// tsysColOut.put(outrowCount+irow, outtsys);
596	// Vector<Double> integ = intCol.getColumn() ;
597	// MaskedArray<Double> mi = maskedArray( integ, flags ) ;
598	// Double intsum = sum(mi);
599	// intColOut.put(outrowCount+irow, intsum);
600	// }
601	// outrowCount += rowNum ;
602	// ++iter;
603	}
604	return out;
605	}
606
607	CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
608	bool droprows)
609	{
610	if (insitu_) {
611	return in;
612	}
613	else {
614	// clone
615	return CountedPtr<Scantable>(new Scantable(*in, Bool(droprows)));
616	}
617	}
618
619	CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
620	float val,
621	const std::string& mode,
622	bool tsys,
623	bool skip_flaggedrow )
624	{
625	CountedPtr< Scantable > out = getScantable(in, false);
626	Table& tab = out->table();
627	ArrayColumn<Float> specCol(tab,"SPECTRA");
628	ArrayColumn<Float> tsysCol(tab,"TSYS");
629	ScalarColumn<uInt> flagrowCol(tab, "FLAGROW");
630	if (mode=="DIV") val = 1.0/val ;
631	else if (mode=="SUB") val *= -1.0 ;
632	for (uInt i=0; i<tab.nrow(); ++i) {
633	Vector<Float> spec;
634	Vector<Float> ts;
635	uInt flagrow;
636	specCol.get(i, spec);
637	tsysCol.get(i, ts);
638	flagrowCol.get(i, flagrow);
639	if (skip_flaggedrow && (flagrow > 0)) continue;
640	if (mode == "MUL" \|\| mode == "DIV") {
641	//if (mode == "DIV") val = 1.0/val;
642	spec *= val;
643	specCol.put(i, spec);
644	if ( tsys ) {
645	ts *= val;
646	tsysCol.put(i, ts);
647	}
648	} else if ( mode == "ADD" \|\| mode == "SUB") {
649	//if (mode == "SUB") val *= -1.0;
650	spec += val;
651	specCol.put(i, spec);
652	if ( tsys ) {
653	ts += val;
654	tsysCol.put(i, ts);
655	}
656	}
657	}
658	return out;
659	}
660
661	CountedPtr< Scantable > STMath::arrayOperate( const CountedPtr< Scantable >& in,
662	const std::vector<float> val,
663	const std::string& mode,
664	const std::string& opmode,
665	bool tsys,
666	bool skip_flaggedrow )
667	{
668	CountedPtr< Scantable > out ;
669	if ( opmode == "channel" ) {
670	out = arrayOperateChannel( in, val, mode, tsys, skip_flaggedrow ) ;
671	}
672	else if ( opmode == "row" ) {
673	out = arrayOperateRow( in, val, mode, tsys, skip_flaggedrow ) ;
674	}
675	else {
676	throw( AipsError( "Unknown array operation mode." ) ) ;
677	}
678	return out ;
679	}
680
681	CountedPtr< Scantable > STMath::arrayOperateChannel( const CountedPtr< Scantable >& in,
682	const std::vector<float> val,
683	const std::string& mode,
684	bool tsys,
685	bool skip_flaggedrow )
686	{
687	if ( val.size() == 1 ){
688	return unaryOperate( in, val[0], mode, tsys ) ;
689	}
690
691	// conformity of SPECTRA and TSYS
692	if ( tsys ) {
693	TableIterator titer(in->table(), "IFNO");
694	while ( !titer.pastEnd() ) {
695	ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
696	ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
697	Array<Float> spec = specCol.getColumn() ;
698	Array<Float> ts = tsysCol.getColumn() ;
699	if ( !spec.conform( ts ) ) {
700	throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
701	}
702	titer.next() ;
703	}
704	}
705
706	// check if all spectra in the scantable have the same number of channel
707	vector<uInt> nchans;
708	vector<uInt> ifnos = in->getIFNos() ;
709	for ( uInt i = 0 ; i < ifnos.size() ; i++ ) {
710	nchans.push_back( in->nchan( ifnos[i] ) ) ;
711	}
712	Vector<uInt> mchans( nchans ) ;
713	if ( anyNE( mchans, mchans[0] ) ) {
714	throw( AipsError("All spectra in the input scantable must have the same number of channel for vector operation." ) ) ;
715	}
716
717	// check if vector size is equal to nchan
718	Vector<Float> fact( val ) ;
719	if ( fact.nelements() != mchans[0] ) {
720	throw( AipsError("Vector size must be 1 or be same as number of channel.") ) ;
721	}
722
723	// check divided by zero
724	if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
725	throw( AipsError("Divided by zero is not recommended." ) ) ;
726	}
727
728	CountedPtr< Scantable > out = getScantable(in, false);
729	Table& tab = out->table();
730	ArrayColumn<Float> specCol(tab,"SPECTRA");
731	ArrayColumn<Float> tsysCol(tab,"TSYS");
732	ScalarColumn<uInt> flagrowCol(tab, "FLAGROW");
733	if (mode == "DIV") fact = (float)1.0 / fact;
734	else if (mode == "SUB") fact *= (float)-1.0 ;
735	for (uInt i=0; i<tab.nrow(); ++i) {
736	Vector<Float> spec;
737	Vector<Float> ts;
738	uInt flagrow;
739	specCol.get(i, spec);
740	tsysCol.get(i, ts);
741	flagrowCol.get(i, flagrow);
742	if (skip_flaggedrow && (flagrow > 0)) continue;
743	if (mode == "MUL" \|\| mode == "DIV") {
744	//if (mode == "DIV") fact = (float)1.0 / fact;
745	spec *= fact;
746	specCol.put(i, spec);
747	if ( tsys ) {
748	ts *= fact;
749	tsysCol.put(i, ts);
750	}
751	} else if ( mode == "ADD" \|\| mode == "SUB") {
752	//if (mode == "SUB") fact *= (float)-1.0 ;
753	spec += fact;
754	specCol.put(i, spec);
755	if ( tsys ) {
756	ts += fact;
757	tsysCol.put(i, ts);
758	}
759	}
760	}
761	return out;
762	}
763
764	CountedPtr< Scantable > STMath::arrayOperateRow( const CountedPtr< Scantable >& in,
765	const std::vector<float> val,
766	const std::string& mode,
767	bool tsys,
768	bool skip_flaggedrow )
769	{
770	if ( val.size() == 1 ) {
771	return unaryOperate( in, val[0], mode, tsys ) ;
772	}
773
774	// conformity of SPECTRA and TSYS
775	if ( tsys ) {
776	TableIterator titer(in->table(), "IFNO");
777	while ( !titer.pastEnd() ) {
778	ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
779	ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
780	Array<Float> spec = specCol.getColumn() ;
781	Array<Float> ts = tsysCol.getColumn() ;
782	if ( !spec.conform( ts ) ) {
783	throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
784	}
785	titer.next() ;
786	}
787	}
788
789	// check if vector size is equal to nrow
790	Vector<Float> fact( val ) ;
791	if (fact.nelements() != uInt(in->nrow())) {
792	throw( AipsError("Vector size must be 1 or be same as number of row.") ) ;
793	}
794
795	// check divided by zero
796	if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
797	throw( AipsError("Divided by zero is not recommended." ) ) ;
798	}
799
800	CountedPtr< Scantable > out = getScantable(in, false);
801	Table& tab = out->table();
802	ArrayColumn<Float> specCol(tab,"SPECTRA");
803	ArrayColumn<Float> tsysCol(tab,"TSYS");
804	ScalarColumn<uInt> flagrowCol(tab, "FLAGROW");
805	if (mode == "DIV") fact = (float)1.0 / fact;
806	if (mode == "SUB") fact *= (float)-1.0 ;
807	for (uInt i=0; i<tab.nrow(); ++i) {
808	Vector<Float> spec;
809	Vector<Float> ts;
810	uInt flagrow;
811	specCol.get(i, spec);
812	tsysCol.get(i, ts);
813	flagrowCol.get(i, flagrow);
814	if (skip_flaggedrow && (flagrow > 0)) continue;
815	if (mode == "MUL" \|\| mode == "DIV") {
816	spec *= fact[i];
817	specCol.put(i, spec);
818	if ( tsys ) {
819	ts *= fact[i];
820	tsysCol.put(i, ts);
821	}
822	} else if ( mode == "ADD" \|\| mode == "SUB") {
823	spec += fact[i];
824	specCol.put(i, spec);
825	if ( tsys ) {
826	ts += fact[i];
827	tsysCol.put(i, ts);
828	}
829	}
830	}
831	return out;
832	}
833
834	CountedPtr< Scantable > STMath::array2dOperate( const CountedPtr< Scantable >& in,
835	const std::vector< std::vector<float> > val,
836	const std::string& mode,
837	bool tsys )
838	{
839	// conformity of SPECTRA and TSYS
840	if ( tsys ) {
841	TableIterator titer(in->table(), "IFNO");
842	while ( !titer.pastEnd() ) {
843	ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
844	ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
845	Array<Float> spec = specCol.getColumn() ;
846	Array<Float> ts = tsysCol.getColumn() ;
847	if ( !spec.conform( ts ) ) {
848	throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
849	}
850	titer.next() ;
851	}
852	}
853
854	// some checks
855	vector<uInt> nchans;
856	for (Int i = 0 ; i < in->nrow() ; i++) {
857	nchans.push_back((in->getSpectrum(i)).size());
858	}
859	//Vector<uInt> mchans( nchans ) ;
860	vector< Vector<Float> > facts ;
861	for ( uInt i = 0 ; i < nchans.size() ; i++ ) {
862	Vector<Float> tmp( val[i] ) ;
863	// check divided by zero
864	if ( ( mode == "DIV" ) && anyEQ( tmp, (float)0.0 ) ) {
865	throw( AipsError("Divided by zero is not recommended." ) ) ;
866	}
867	// conformity check
868	if ( tmp.nelements() != nchans[i] ) {
869	stringstream ss ;
870	ss << "Row " << i << ": Vector size must be same as number of channel." ;
871	throw( AipsError( ss.str() ) ) ;
872	}
873	facts.push_back( tmp ) ;
874	}
875
876
877	CountedPtr< Scantable > out = getScantable(in, false);
878	Table& tab = out->table();
879	ArrayColumn<Float> specCol(tab,"SPECTRA");
880	ArrayColumn<Float> tsysCol(tab,"TSYS");
881	for (uInt i=0; i<tab.nrow(); ++i) {
882	Vector<Float> fact = facts[i] ;
883	Vector<Float> spec;
884	Vector<Float> ts;
885	specCol.get(i, spec);
886	tsysCol.get(i, ts);
887	if (mode == "MUL" \|\| mode == "DIV") {
888	if (mode == "DIV") fact = (float)1.0 / fact;
889	spec *= fact;
890	specCol.put(i, spec);
891	if ( tsys ) {
892	ts *= fact;
893	tsysCol.put(i, ts);
894	}
895	} else if ( mode == "ADD" \|\| mode == "SUB") {
896	if (mode == "SUB") fact *= (float)-1.0 ;
897	spec += fact;
898	specCol.put(i, spec);
899	if ( tsys ) {
900	ts += fact;
901	tsysCol.put(i, ts);
902	}
903	}
904	}
905	return out;
906	}
907
908	CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
909	const CountedPtr<Scantable>& right,
910	const std::string& mode)
911	{
912	bool insitu = insitu_;
913	if ( ! left->conformant(*right) ) {
914	throw(AipsError("'left' and 'right' scantables are not conformant."));
915	}
916	setInsitu(false);
917	CountedPtr< Scantable > out = getScantable(left, false);
918	setInsitu(insitu);
919	Table& tout = out->table();
920	Block<String> coln(5);
921	coln[0] = "SCANNO"; coln[1] = "CYCLENO"; coln[2] = "BEAMNO";
922	coln[3] = "IFNO"; coln[4] = "POLNO";
923	Table tmpl = tout.sort(coln);
924	Table tmpr = right->table().sort(coln);
925	ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
926	ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
927	ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
928	ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
929
930	for (uInt i=0; i<tout.nrow(); ++i) {
931	Vector<Float> lspecvec, rspecvec;
932	Vector<uChar> lflagvec, rflagvec;
933	lspecvec = lspecCol(i); rspecvec = rspecCol(i);
934	lflagvec = lflagCol(i); rflagvec = rflagCol(i);
935	MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
936	MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
937	if (mode == "ADD") {
938	mleft += mright;
939	} else if ( mode == "SUB") {
940	mleft -= mright;
941	} else if ( mode == "MUL") {
942	mleft *= mright;
943	} else if ( mode == "DIV") {
944	mleft /= mright;
945	} else {
946	throw(AipsError("Illegal binary operator"));
947	}
948	lspecCol.put(i, mleft.getArray());
949	}
950	return out;
951	}
952
953
954
955	MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
956	const Vector<uChar>& f)
957	{
958	Vector<Bool> mask;
959	mask.resize(f.shape());
960	convertArray(mask, f);
961	return MaskedArray<Float>(s,!mask);
962	}
963
964	MaskedArray<Double> STMath::maskedArray( const Vector<Double>& s,
965	const Vector<uChar>& f)
966	{
967	Vector<Bool> mask;
968	mask.resize(f.shape());
969	convertArray(mask, f);
970	return MaskedArray<Double>(s,!mask);
971	}
972
973	Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
974	{
975	const Vector<Bool>& m = ma.getMask();
976	Vector<uChar> flags(m.shape());
977	convertArray(flags, !m);
978	return flags;
979	}
980
981	CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
982	const std::string & mode,
983	bool preserve )
984	{
985	/// @todo make other modes available
986	/// modes should be "nearest", "pair"
987	// make this operation non insitu
988	(void) mode; //currently unused
989	const Table& tin = in->table();
990	Table ons = tin(tin.col("SRCTYPE") == Int(SrcType::PSON));
991	Table offs = tin(tin.col("SRCTYPE") == Int(SrcType::PSOFF));
992	if ( offs.nrow() == 0 )
993	throw(AipsError("No 'off' scans present."));
994	// put all "on" scans into output table
995
996	bool insitu = insitu_;
997	setInsitu(false);
998	CountedPtr< Scantable > out = getScantable(in, true);
999	setInsitu(insitu);
1000	Table& tout = out->table();
1001
1002	TableCopy::copyRows(tout, ons);
1003	TableRow row(tout);
1004	ROScalarColumn<Double> offtimeCol(offs, "TIME");
1005	ArrayColumn<Float> outspecCol(tout, "SPECTRA");
1006	ROArrayColumn<Float> outtsysCol(tout, "TSYS");
1007	ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
1008	for (uInt i=0; i < tout.nrow(); ++i) {
1009	const TableRecord& rec = row.get(i);
1010	Double ontime = rec.asDouble("TIME");
1011	Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1012	&& offs.col("IFNO") == Int(rec.asuInt("IFNO"))
1013	&& offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
1014	ROScalarColumn<Double> offtimeCol(presel, "TIME");
1015
1016	Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
1017	// Timestamp may vary within a cycle ???!!!
1018	// increase this by 0.01 sec in case of rounding errors...
1019	// There might be a better way to do this.
1020	// fix to this fix. TIME is MJD, so 1.0d not 1.0s
1021	mindeltat += 0.01/24./60./60.;
1022	Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
1023
1024	if ( sel.nrow() < 1 ) {
1025	throw(AipsError("No closest in time found... This could be a rounding "
1026	"issue. Try quotient instead."));
1027	}
1028	TableRow offrow(sel);
1029	const TableRecord& offrec = offrow.get(0);//should only be one row
1030	RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
1031	RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
1032	RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
1033	/// @fixme this assumes tsys is a scalar not vector
1034	Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
1035	Vector<Float> specon, tsyson;
1036	outtsysCol.get(i, tsyson);
1037	outspecCol.get(i, specon);
1038	Vector<uChar> flagon;
1039	outflagCol.get(i, flagon);
1040	MaskedArray<Float> mon = maskedArray(specon, flagon);
1041	MaskedArray<Float> moff = maskedArray(specoff, flagoff);
1042	MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
1043	if (preserve) {
1044	quot -= tsysoffscalar;
1045	} else {
1046	quot -= tsyson[0];
1047	}
1048	outspecCol.put(i, quot.getArray());
1049	outflagCol.put(i, flagsFromMA(quot));
1050	}
1051	// renumber scanno
1052	TableIterator it(tout, "SCANNO");
1053	uInt i = 0;
1054	while ( !it.pastEnd() ) {
1055	Table t = it.table();
1056	TableVector<uInt> vec(t, "SCANNO");
1057	vec = i;
1058	++i;
1059	++it;
1060	}
1061	return out;
1062	}
1063
1064
1065	CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
1066	const CountedPtr< Scantable > & off,
1067	bool preserve )
1068	{
1069	bool insitu = insitu_;
1070	if ( ! on->conformant(*off) ) {
1071	throw(AipsError("'on' and 'off' scantables are not conformant."));
1072	}
1073	setInsitu(false);
1074	CountedPtr< Scantable > out = getScantable(on, false);
1075	setInsitu(insitu);
1076	Table& tout = out->table();
1077	const Table& toff = off->table();
1078	TableIterator sit(tout, "SCANNO");
1079	TableIterator s2it(toff, "SCANNO");
1080	while ( !sit.pastEnd() ) {
1081	Table ton = sit.table();
1082	TableRow row(ton);
1083	Table t = s2it.table();
1084	ArrayColumn<Float> outspecCol(ton, "SPECTRA");
1085	ROArrayColumn<Float> outtsysCol(ton, "TSYS");
1086	ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
1087	for (uInt i=0; i < ton.nrow(); ++i) {
1088	const TableRecord& rec = row.get(i);
1089	Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1090	&& t.col("IFNO") == Int(rec.asuInt("IFNO"))
1091	&& t.col("POLNO") == Int(rec.asuInt("POLNO")) );
1092	if ( offsel.nrow() == 0 )
1093	throw AipsError("STMath::quotient: no matching off");
1094	TableRow offrow(offsel);
1095	const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
1096	RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
1097	RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
1098	RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
1099	Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
1100	Vector<Float> specon, tsyson;
1101	outtsysCol.get(i, tsyson);
1102	outspecCol.get(i, specon);
1103	Vector<uChar> flagon;
1104	outflagCol.get(i, flagon);
1105	MaskedArray<Float> mon = maskedArray(specon, flagon);
1106	MaskedArray<Float> moff = maskedArray(specoff, flagoff);
1107	MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
1108	if (preserve) {
1109	quot -= tsysoffscalar;
1110	} else {
1111	quot -= tsyson[0];
1112	}
1113	outspecCol.put(i, quot.getArray());
1114	outflagCol.put(i, flagsFromMA(quot));
1115	}
1116	++sit;
1117	++s2it;
1118	// take the first off for each on scan which doesn't have a
1119	// matching off scan
1120	// non <= noff: matching pairs, non > noff matching pairs then first off
1121	if ( s2it.pastEnd() ) s2it.reset();
1122	}
1123	return out;
1124	}
1125
1126	// dototalpower (migration of GBTIDL procedure dototalpower.pro)
1127	// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
1128	// do it for each cycles in a specific scan.
1129	CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
1130	const CountedPtr< Scantable >& caloff, Float tcal )
1131	{
1132	if ( ! calon->conformant(*caloff) ) {
1133	throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
1134	}
1135	setInsitu(false);
1136	CountedPtr< Scantable > out = getScantable(caloff, false);
1137	Table& tout = out->table();
1138	const Table& tcon = calon->table();
1139	Vector<Float> tcalout;
1140
1141	std::map<uInt,uInt> tcalIdToRecNoMap;
1142	const Table& calOffTcalTable = caloff->tcal().table();
1143	{
1144	ROScalarColumn<uInt> calOffTcalTable_IDcol(calOffTcalTable, "ID");
1145	const Vector<uInt> tcalIds(calOffTcalTable_IDcol.getColumn());
1146	size_t tcalIdsEnd = tcalIds.nelements();
1147	for (uInt i = 0; i < tcalIdsEnd; i++) {
1148	tcalIdToRecNoMap[tcalIds[i]] = i;
1149	}
1150	}
1151	ROArrayColumn<Float> calOffTcalTable_TCALcol(calOffTcalTable, "TCAL");
1152
1153	if ( tout.nrow() != tcon.nrow() ) {
1154	throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
1155	}
1156	// iteration by scanno or cycle no.
1157	TableIterator sit(tout, "SCANNO");
1158	TableIterator s2it(tcon, "SCANNO");
1159	while ( !sit.pastEnd() ) {
1160	Table toff = sit.table();
1161	TableRow row(toff);
1162	Table t = s2it.table();
1163	ScalarColumn<Double> outintCol(toff, "INTERVAL");
1164	ArrayColumn<Float> outspecCol(toff, "SPECTRA");
1165	ArrayColumn<Float> outtsysCol(toff, "TSYS");
1166	ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
1167	ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
1168	ROScalarColumn<uInt> outpolCol(toff, "POLNO");
1169	ROScalarColumn<Double> onintCol(t, "INTERVAL");
1170	ROArrayColumn<Float> onspecCol(t, "SPECTRA");
1171	ROArrayColumn<Float> ontsysCol(t, "TSYS");
1172	ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
1173	//ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
1174
1175	for (uInt i=0; i < toff.nrow(); ++i) {
1176	//skip these checks -> assumes the data order are the same between the cal on off pairs
1177	//
1178	Vector<Float> specCalon, specCaloff;
1179	// to store scalar (mean) tsys
1180	Vector<Float> tsysout(1);
1181	uInt tcalId, polno;
1182	Double offint, onint;
1183	outpolCol.get(i, polno);
1184	outspecCol.get(i, specCaloff);
1185	onspecCol.get(i, specCalon);
1186	Vector<uChar> flagCaloff, flagCalon;
1187	outflagCol.get(i, flagCaloff);
1188	onflagCol.get(i, flagCalon);
1189	outtcalIdCol.get(i, tcalId);
1190	outintCol.get(i, offint);
1191	onintCol.get(i, onint);
1192	// caluculate mean Tsys
1193	uInt nchan = specCaloff.nelements();
1194	// percentage of edge cut off
1195	uInt pc = 10;
1196	uInt bchan = nchan/pc;
1197	uInt echan = nchan-bchan;
1198
1199	Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
1200	Vector<Float> testsubsp = specCaloff(chansl);
1201	MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
1202	MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
1203	MaskedArray<Float> spdiff = spon-spoff;
1204	uInt noff = spoff.nelementsValid();
1205	//uInt non = spon.nelementsValid();
1206	uInt ndiff = spdiff.nelementsValid();
1207	Float meantsys;
1208
1209	/**
1210	Double subspec, subdiff;
1211	uInt usednchan;
1212	subspec = 0;
1213	subdiff = 0;
1214	usednchan = 0;
1215	for(uInt k=(bchan-1); k<echan; k++) {
1216	subspec += specCaloff[k];
1217	subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
1218	++usednchan;
1219	}
1220	**/
1221	// get tcal if input tcal <= 0
1222	Float tcalUsed;
1223	tcalUsed = tcal;
1224	if ( tcal <= 0.0 ) {
1225	uInt tcalRecNo = tcalIdToRecNoMap[tcalId];
1226	calOffTcalTable_TCALcol.get(tcalRecNo, tcalout);
1227	// if (polno<=3) {
1228	// tcalUsed = tcalout[polno];
1229	// }
1230	// else {
1231	// tcalUsed = tcalout[0];
1232	// }
1233	if ( tcalout.size() == 1 )
1234	tcalUsed = tcalout[0] ;
1235	else if ( tcalout.size() == nchan )
1236	tcalUsed = mean(tcalout) ;
1237	else {
1238	uInt ipol = polno ;
1239	if ( ipol > 3 ) ipol = 0 ;
1240	tcalUsed = tcalout[ipol] ;
1241	}
1242	}
1243
1244	Float meanoff;
1245	Float meandiff;
1246	if (noff && ndiff) {
1247	//Debug
1248	//if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
1249	//LogIO os( LogOrigin( "STMath", "dototalpower()", WHERE ) ) ;
1250	//if(noff!=ndiff) os<<"noff and ndiff is not equal"<<LogIO::POST;
1251	meanoff = sum(spoff)/noff;
1252	meandiff = sum(spdiff)/ndiff;
1253	meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
1254	}
1255	else {
1256	meantsys=1;
1257	}
1258
1259	tsysout[0] = Float(meantsys);
1260	MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
1261	MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
1262	MaskedArray<Float> sig = Float(0.5) * (mcaloff + mcalon);
1263	//uInt ncaloff = mcaloff.nelementsValid();
1264	//uInt ncalon = mcalon.nelementsValid();
1265
1266	outintCol.put(i, offint+onint);
1267	outspecCol.put(i, sig.getArray());
1268	outflagCol.put(i, flagsFromMA(sig));
1269	outtsysCol.put(i, tsysout);
1270	}
1271	++sit;
1272	++s2it;
1273	}
1274	return out;
1275	}
1276
1277	//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
1278	// observatiions.
1279	// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
1280	// no smoothing).
1281	// output: resultant scantable [= (sig-ref/ref)*tsys]
1282	CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
1283	const CountedPtr < Scantable >& ref,
1284	int smoothref,
1285	casa::Float tsysv,
1286	casa::Float tau )
1287	{
1288	LogIO os( casa::LogOrigin( "STMath", "dosigref()"));
1289	if ( ! ref->conformant(*sig) ) {
1290	throw(AipsError("'sig' and 'ref' scantables are not conformant."));
1291	}
1292	setInsitu(false);
1293	CountedPtr< Scantable > out = getScantable(sig, false);
1294	CountedPtr< Scantable > smref;
1295	if ( smoothref > 1 ) {
1296	float fsmoothref = static_cast<float>(smoothref);
1297	std::string inkernel = "boxcar";
1298	smref = smooth(ref, inkernel, fsmoothref );
1299	ostringstream oss;
1300	os <<"Applied smoothing of "<<fsmoothref<<" on the reference."
1301	<< LogIO::POST;
1302	}
1303	else {
1304	smref = ref;
1305	}
1306	Table& tout = out->table();
1307	const Table& tref = smref->table();
1308	if ( tout.nrow() != tref.nrow() ) {
1309	throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
1310	}
1311	// iteration by scanno? or cycle no.
1312	TableIterator sit(tout, "SCANNO");
1313	TableIterator s2it(tref, "SCANNO");
1314	while ( !sit.pastEnd() ) {
1315	Table ton = sit.table();
1316	Table t = s2it.table();
1317	ScalarColumn<Double> outintCol(ton, "INTERVAL");
1318	ArrayColumn<Float> outspecCol(ton, "SPECTRA");
1319	ArrayColumn<Float> outtsysCol(ton, "TSYS");
1320	ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
1321	ArrayColumn<Float> refspecCol(t, "SPECTRA");
1322	ROScalarColumn<Double> refintCol(t, "INTERVAL");
1323	ROArrayColumn<Float> reftsysCol(t, "TSYS");
1324	ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
1325	ROScalarColumn<Float> refelevCol(t, "ELEVATION");
1326	for (uInt i=0; i < ton.nrow(); ++i) {
1327
1328	Double onint, refint;
1329	Vector<Float> specon, specref;
1330	// to store scalar (mean) tsys
1331	Vector<Float> tsysref;
1332	outintCol.get(i, onint);
1333	refintCol.get(i, refint);
1334	outspecCol.get(i, specon);
1335	refspecCol.get(i, specref);
1336	Vector<uChar> flagref, flagon;
1337	outflagCol.get(i, flagon);
1338	refflagCol.get(i, flagref);
1339	reftsysCol.get(i, tsysref);
1340
1341	Float tsysrefscalar;
1342	if ( tsysv > 0.0 ) {
1343	ostringstream oss;
1344	Float elev;
1345	refelevCol.get(i, elev);
1346	os << "user specified Tsys = " << tsysv;
1347	// do recalc elevation if EL = 0
1348	if ( elev == 0 ) {
1349	throw(AipsError("EL=0, elevation data is missing."));
1350	} else {
1351	if ( tau <= 0.0 ) {
1352	throw(AipsError("Valid tau is not supplied."));
1353	} else {
1354	tsysrefscalar = tsysv * exp(tau/elev);
1355	}
1356	}
1357	os << ", corrected (for El) tsys= "<<tsysrefscalar;
1358	}
1359	else {
1360	tsysrefscalar = tsysref[0];
1361	}
1362	//get quotient spectrum
1363	MaskedArray<Float> mref = maskedArray(specref, flagref);
1364	MaskedArray<Float> mon = maskedArray(specon, flagon);
1365	MaskedArray<Float> specres = tsysrefscalar*((mon - mref)/mref);
1366	Double resint = onintrefintsmoothref/(onint+refint*smoothref);
1367
1368	//Debug
1369	//cerr<<"Tsys used="<<tsysrefscalar<<endl;
1370	//LogIO os( LogOrigin( "STMath", "dosigref", WHERE ) ) ;
1371	//os<<"Tsys used="<<tsysrefscalar<<LogIO::POST;
1372	// fill the result, replay signal tsys by reference tsys
1373	outintCol.put(i, resint);
1374	outspecCol.put(i, specres.getArray());
1375	outflagCol.put(i, flagsFromMA(specres));
1376	outtsysCol.put(i, tsysref);
1377	}
1378	++sit;
1379	++s2it;
1380	}
1381
1382	out->setFluxUnit("K");
1383
1384	return out;
1385	}
1386
1387	CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
1388	const std::vector<int>& scans,
1389	int smoothref,
1390	casa::Float tsysv,
1391	casa::Float tau,
1392	casa::Float tcal )
1393
1394	{
1395	setInsitu(false);
1396	LogIO os( casa::LogOrigin( "STMath", "donod()"));
1397	STSelector sel;
1398	std::vector<int> scan1, scan2, beams, types;
1399	std::vector< vector<int> > scanpair;
1400	//std::vector<string> calstate;
1401	std::vector<int> calstate;
1402	String msg;
1403
1404	CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
1405	CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
1406
1407	std::vector< CountedPtr< Scantable > > sctables;
1408	sctables.push_back(s1b1on);
1409	sctables.push_back(s1b1off);
1410	sctables.push_back(s1b2on);
1411	sctables.push_back(s1b2off);
1412	sctables.push_back(s2b1on);
1413	sctables.push_back(s2b1off);
1414	sctables.push_back(s2b2on);
1415	sctables.push_back(s2b2off);
1416
1417	//check scanlist
1418	int n=s->checkScanInfo(scans);
1419	if (n==1) {
1420	throw(AipsError("Incorrect scan pairs. "));
1421	}
1422
1423	// Assume scans contain only a pair of consecutive scan numbers.
1424	// It is assumed that first beam, b1, is on target.
1425	// There is no check if the first beam is on or not.
1426	if ( scans.size()==1 ) {
1427	scan1.push_back(scans[0]);
1428	scan2.push_back(scans[0]+1);
1429	} else if ( scans.size()==2 ) {
1430	scan1.push_back(scans[0]);
1431	scan2.push_back(scans[1]);
1432	} else {
1433	if ( scans.size()%2 == 0 ) {
1434	for (uInt i=0; i<scans.size(); i++) {
1435	if (i%2 == 0) {
1436	scan1.push_back(scans[i]);
1437	}
1438	else {
1439	scan2.push_back(scans[i]);
1440	}
1441	}
1442	} else {
1443	throw(AipsError("Odd numbers of scans, cannot form pairs."));
1444	}
1445	}
1446	scanpair.push_back(scan1);
1447	scanpair.push_back(scan2);
1448	//calstate.push_back("*calon");
1449	//calstate.push_back("*[^calon]");
1450	calstate.push_back(SrcType::NODCAL);
1451	calstate.push_back(SrcType::NOD);
1452	CountedPtr< Scantable > ws = getScantable(s, false);
1453	uInt l=0;
1454	while ( l < sctables.size() ) {
1455	for (uInt i=0; i < 2; i++) {
1456	for (uInt j=0; j < 2; j++) {
1457	for (uInt k=0; k < 2; k++) {
1458	sel.reset();
1459	sel.setScans(scanpair[i]);
1460	//sel.setName(calstate[k]);
1461	types.clear();
1462	types.push_back(calstate[k]);
1463	sel.setTypes(types);
1464	beams.clear();
1465	beams.push_back(j);
1466	sel.setBeams(beams);
1467	ws->setSelection(sel);
1468	sctables[l]= getScantable(ws, false);
1469	l++;
1470	}
1471	}
1472	}
1473	}
1474
1475	// replace here by splitData or getData functionality
1476	CountedPtr< Scantable > sig1;
1477	CountedPtr< Scantable > ref1;
1478	CountedPtr< Scantable > sig2;
1479	CountedPtr< Scantable > ref2;
1480	CountedPtr< Scantable > calb1;
1481	CountedPtr< Scantable > calb2;
1482
1483	msg=String("Processing dototalpower for subset of the data");
1484	os << msg << LogIO::POST;
1485	// Debug for IRC CS data
1486	//float tcal1=7.0;
1487	//float tcal2=4.0;
1488	sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
1489	ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
1490	ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
1491	sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
1492
1493	// correction of user-specified tsys for elevation here
1494
1495	// dosigref calibration
1496	msg=String("Processing dosigref for subset of the data");
1497	os << msg << endl;
1498	calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
1499	calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
1500
1501	// iteration by scanno or cycle no.
1502	Table& tcalb1 = calb1->table();
1503	Table& tcalb2 = calb2->table();
1504	TableIterator sit(tcalb1, "SCANNO");
1505	TableIterator s2it(tcalb2, "SCANNO");
1506	while ( !sit.pastEnd() ) {
1507	Table t1 = sit.table();
1508	Table t2= s2it.table();
1509	ArrayColumn<Float> outspecCol(t1, "SPECTRA");
1510	ArrayColumn<Float> outtsysCol(t1, "TSYS");
1511	ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
1512	ScalarColumn<Double> outintCol(t1, "INTERVAL");
1513	ArrayColumn<Float> t2specCol(t2, "SPECTRA");
1514	ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
1515	ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
1516	ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
1517	for (uInt i=0; i < t1.nrow(); ++i) {
1518	Vector<Float> spec1, spec2;
1519	// to store scalar (mean) tsys
1520	Vector<Float> tsys1, tsys2;
1521	Vector<uChar> flag1, flag2;
1522	Double tint1, tint2;
1523	outspecCol.get(i, spec1);
1524	t2specCol.get(i, spec2);
1525	outflagCol.get(i, flag1);
1526	t2flagCol.get(i, flag2);
1527	outtsysCol.get(i, tsys1);
1528	t2tsysCol.get(i, tsys2);
1529	outintCol.get(i, tint1);
1530	t2intCol.get(i, tint2);
1531	// average
1532	// assume scalar tsys for weights
1533	Float wt1, wt2, tsyssq1, tsyssq2;
1534	tsyssq1 = tsys1[0]*tsys1[0];
1535	tsyssq2 = tsys2[0]*tsys2[0];
1536	wt1 = Float(tint1)/tsyssq1;
1537	wt2 = Float(tint2)/tsyssq2;
1538	Float invsumwt=1/(wt1+wt2);
1539	MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
1540	MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
1541	MaskedArray<Float> avspec = invsumwt * (wt1mspec1 + wt2mspec2);
1542	//Array<Float> avtsys = Float(0.5) * (tsys1 + tsys2);
1543	// cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
1544	// LogIO os( LogOrigin( "STMath", "donod", WHERE ) ) ;
1545	// os<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<LogIO::POST;
1546	tsys1[0] = sqrt(tsyssq1 + tsyssq2);
1547	Array<Float> avtsys = tsys1;
1548
1549	outspecCol.put(i, avspec.getArray());
1550	outflagCol.put(i, flagsFromMA(avspec));
1551	outtsysCol.put(i, avtsys);
1552	}
1553	++sit;
1554	++s2it;
1555	}
1556
1557	calb1->setFluxUnit("K");
1558
1559	return calb1;
1560	}
1561
1562	//GBTIDL version of frequency switched data calibration
1563	CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
1564	const std::vector<int>& scans,
1565	int smoothref,
1566	casa::Float tsysv,
1567	casa::Float tau,
1568	casa::Float tcal )
1569	{
1570
1571
1572	(void) scans; //currently unused
1573	STSelector sel;
1574	CountedPtr< Scantable > ws = getScantable(s, false);
1575	CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
1576	CountedPtr< Scantable > calsig, calref, out, out1, out2;
1577	Bool nofold=False;
1578	vector<int> types ;
1579
1580	//split the data
1581	//sel.setName("*_fs");
1582	types.push_back( SrcType::FSON ) ;
1583	sel.setTypes( types ) ;
1584	ws->setSelection(sel);
1585	sig = getScantable(ws,false);
1586	sel.reset();
1587	types.clear() ;
1588	//sel.setName("*_fs_calon");
1589	types.push_back( SrcType::FONCAL ) ;
1590	sel.setTypes( types ) ;
1591	ws->setSelection(sel);
1592	sigwcal = getScantable(ws,false);
1593	sel.reset();
1594	types.clear() ;
1595	//sel.setName("*_fsr");
1596	types.push_back( SrcType::FSOFF ) ;
1597	sel.setTypes( types ) ;
1598	ws->setSelection(sel);
1599	ref = getScantable(ws,false);
1600	sel.reset();
1601	types.clear() ;
1602	//sel.setName("*_fsr_calon");
1603	types.push_back( SrcType::FOFFCAL ) ;
1604	sel.setTypes( types ) ;
1605	ws->setSelection(sel);
1606	refwcal = getScantable(ws,false);
1607	sel.reset() ;
1608	types.clear() ;
1609
1610	calsig = dototalpower(sigwcal, sig, tcal=tcal);
1611	calref = dototalpower(refwcal, ref, tcal=tcal);
1612
1613	out1=dosigref(calsig,calref,smoothref,tsysv,tau);
1614	out2=dosigref(calref,calsig,smoothref,tsysv,tau);
1615
1616	Table& tabout1=out1->table();
1617	Table& tabout2=out2->table();
1618	ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
1619	ScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
1620	ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
1621	Vector<Float> spec; specCol.get(0, spec);
1622	uInt nchan = spec.nelements();
1623	uInt freqid1; freqidCol1.get(0,freqid1);
1624	uInt freqid2; freqidCol2.get(0,freqid2);
1625	Double rp1, rp2, rv1, rv2, inc1, inc2;
1626	out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
1627	out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
1628	//cerr << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << endl ;
1629	//LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
1630	//os << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << LogIO::POST ;
1631	if (rp1==rp2) {
1632	Double foffset = rv1 - rv2;
1633	uInt choffset = static_cast<uInt>(foffset/abs(inc2));
1634	if (choffset >= nchan) {
1635	//cerr<<"out-band frequency switching, no folding"<<endl;
1636	LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
1637	os<<"out-band frequency switching, no folding"<<LogIO::POST;
1638	nofold = True;
1639	}
1640	}
1641
1642	if (nofold) {
1643	std::vector< CountedPtr< Scantable > > tabs;
1644	tabs.push_back(out1);
1645	tabs.push_back(out2);
1646	out = merge(tabs);
1647	}
1648	else {
1649	//out = out1;
1650	Double choffset = ( rv1 - rv2 ) / inc2 ;
1651	out = dofold( out1, out2, choffset ) ;
1652	}
1653
1654	out->setFluxUnit("K");
1655
1656	return out;
1657	}
1658
1659	CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig,
1660	const CountedPtr<Scantable> &ref,
1661	Double choffset,
1662	Double choffset2 )
1663	{
1664	LogIO os( LogOrigin( "STMath", "dofold", WHERE ) ) ;
1665	os << "choffset=" << choffset << " choffset2=" << choffset2 << LogIO::POST ;
1666
1667	// output scantable
1668	CountedPtr<Scantable> out = getScantable( sig, false ) ;
1669
1670	// separate choffset to integer part and decimal part
1671	Int ioffset = (Int)choffset ;
1672	Double doffset = choffset - ioffset ;
1673	Int ioffset2 = (Int)choffset2 ;
1674	Double doffset2 = choffset2 - ioffset2 ;
1675	os << "ioffset=" << ioffset << " doffset=" << doffset << LogIO::POST ;
1676	os << "ioffset2=" << ioffset2 << " doffset2=" << doffset2 << LogIO::POST ;
1677
1678	// get column
1679	ROArrayColumn<Float> specCol1( sig->table(), "SPECTRA" ) ;
1680	ROArrayColumn<Float> specCol2( ref->table(), "SPECTRA" ) ;
1681	ROArrayColumn<Float> tsysCol1( sig->table(), "TSYS" ) ;
1682	ROArrayColumn<Float> tsysCol2( ref->table(), "TSYS" ) ;
1683	ROArrayColumn<uChar> flagCol1( sig->table(), "FLAGTRA" ) ;
1684	ROArrayColumn<uChar> flagCol2( ref->table(), "FLAGTRA" ) ;
1685	ROScalarColumn<Double> mjdCol1( sig->table(), "TIME" ) ;
1686	ROScalarColumn<Double> mjdCol2( ref->table(), "TIME" ) ;
1687	ROScalarColumn<Double> intervalCol1( sig->table(), "INTERVAL" ) ;
1688	ROScalarColumn<Double> intervalCol2( ref->table(), "INTERVAL" ) ;
1689
1690	// check
1691	if ( ioffset == 0 ) {
1692	LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
1693	os << "channel offset is zero, no folding" << LogIO::POST ;
1694	return out ;
1695	}
1696	int nchan = ref->nchan() ;
1697	if ( abs(ioffset) >= nchan ) {
1698	LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
1699	os << "out-band frequency switching, no folding" << LogIO::POST ;
1700	return out ;
1701	}
1702
1703	// attach column for output scantable
1704	ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ;
1705	ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ;
1706	ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ;
1707	ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ;
1708	ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ;
1709	ScalarColumn<uInt> fidColOut( out->table(), "FREQ_ID" ) ;
1710
1711	// for each row
1712	// assume that the data order are same between sig and ref
1713	RowAccumulator acc( asap::W_TINTSYS ) ;
1714	for ( int i = 0 ; i < sig->nrow() ; i++ ) {
1715	// get values
1716	Vector<Float> spsig ;
1717	specCol1.get( i, spsig ) ;
1718	Vector<Float> spref ;
1719	specCol2.get( i, spref ) ;
1720	Vector<Float> tsyssig ;
1721	tsysCol1.get( i, tsyssig ) ;
1722	Vector<Float> tsysref ;
1723	tsysCol2.get( i, tsysref ) ;
1724	Vector<uChar> flagsig ;
1725	flagCol1.get( i, flagsig ) ;
1726	Vector<uChar> flagref ;
1727	flagCol2.get( i, flagref ) ;
1728	Double timesig ;
1729	mjdCol1.get( i, timesig ) ;
1730	Double timeref ;
1731	mjdCol2.get( i, timeref ) ;
1732	Double intsig ;
1733	intervalCol1.get( i, intsig ) ;
1734	Double intref ;
1735	intervalCol2.get( i, intref ) ;
1736
1737	// shift reference spectra
1738	int refchan = spref.nelements() ;
1739	Vector<Float> sspref( spref.nelements() ) ;
1740	Vector<Float> stsysref( tsysref.nelements() ) ;
1741	Vector<uChar> sflagref( flagref.nelements() ) ;
1742	if ( ioffset > 0 ) {
1743	// SPECTRA and FLAGTRA
1744	for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
1745	sspref[j] = spref[j+ioffset] ;
1746	sflagref[j] = flagref[j+ioffset] ;
1747	}
1748	for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
1749	sspref[j] = spref[j-refchan+ioffset] ;
1750	sflagref[j] = flagref[j-refchan+ioffset] ;
1751	}
1752	spref = sspref.copy() ;
1753	flagref = sflagref.copy() ;
1754	for ( int j = 0 ; j < refchan - 1 ; j++ ) {
1755	sspref[j] = doffset * spref[j+1] + ( 1.0 - doffset ) * spref[j] ;
1756	sflagref[j] = flagref[j+1] + flagref[j] ;
1757	}
1758	sspref[refchan-1] = doffset * spref[0] + ( 1.0 - doffset ) * spref[refchan-1] ;
1759	sflagref[refchan-1] = flagref[0] + flagref[refchan-1] ;
1760
1761	// TSYS
1762	if ( spref.nelements() == tsysref.nelements() ) {
1763	for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
1764	stsysref[j] = tsysref[j+ioffset] ;
1765	}
1766	for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
1767	stsysref[j] = tsysref[j-refchan+ioffset] ;
1768	}
1769	tsysref = stsysref.copy() ;
1770	for ( int j = 0 ; j < refchan - 1 ; j++ ) {
1771	stsysref[j] = doffset * tsysref[j+1] + ( 1.0 - doffset ) * tsysref[j] ;
1772	}
1773	stsysref[refchan-1] = doffset * tsysref[0] + ( 1.0 - doffset ) * tsysref[refchan-1] ;
1774	}
1775	}
1776	else {
1777	// SPECTRA and FLAGTRA
1778	for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
1779	sspref[j] = spref[refchan+ioffset+j] ;
1780	sflagref[j] = flagref[refchan+ioffset+j] ;
1781	}
1782	for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
1783	sspref[j] = spref[j+ioffset] ;
1784	sflagref[j] = flagref[j+ioffset] ;
1785	}
1786	spref = sspref.copy() ;
1787	flagref = sflagref.copy() ;
1788	sspref[0] = doffset * spref[refchan-1] + ( 1.0 - doffset ) * spref[0] ;
1789	sflagref[0] = flagref[0] + flagref[refchan-1] ;
1790	for ( int j = 1 ; j < refchan ; j++ ) {
1791	sspref[j] = doffset * spref[j-1] + ( 1.0 - doffset ) * spref[j] ;
1792	sflagref[j] = flagref[j-1] + flagref[j] ;
1793	}
1794	// TSYS
1795	if ( spref.nelements() == tsysref.nelements() ) {
1796	for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
1797	stsysref[j] = tsysref[refchan+ioffset+j] ;
1798	}
1799	for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
1800	stsysref[j] = tsysref[j+ioffset] ;
1801	}
1802	tsysref = stsysref.copy() ;
1803	stsysref[0] = doffset * tsysref[refchan-1] + ( 1.0 - doffset ) * tsysref[0] ;
1804	for ( int j = 1 ; j < refchan ; j++ ) {
1805	stsysref[j] = doffset * tsysref[j-1] + ( 1.0 - doffset ) * tsysref[j] ;
1806	}
1807	}
1808	}
1809
1810	// shift signal spectra if necessary (only for APEX?)
1811	if ( choffset2 != 0.0 ) {
1812	int sigchan = spsig.nelements() ;
1813	Vector<Float> sspsig( spsig.nelements() ) ;
1814	Vector<Float> stsyssig( tsyssig.nelements() ) ;
1815	Vector<uChar> sflagsig( flagsig.nelements() ) ;
1816	if ( ioffset2 > 0 ) {
1817	// SPECTRA and FLAGTRA
1818	for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
1819	sspsig[j] = spsig[j+ioffset2] ;
1820	sflagsig[j] = flagsig[j+ioffset2] ;
1821	}
1822	for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
1823	sspsig[j] = spsig[j-sigchan+ioffset2] ;
1824	sflagsig[j] = flagsig[j-sigchan+ioffset2] ;
1825	}
1826	spsig = sspsig.copy() ;
1827	flagsig = sflagsig.copy() ;
1828	for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
1829	sspsig[j] = doffset2 * spsig[j+1] + ( 1.0 - doffset2 ) * spsig[j] ;
1830	sflagsig[j] = flagsig[j+1] \|\| flagsig[j] ;
1831	}
1832	sspsig[sigchan-1] = doffset2 * spsig[0] + ( 1.0 - doffset2 ) * spsig[sigchan-1] ;
1833	sflagsig[sigchan-1] = flagsig[0] \|\| flagsig[sigchan-1] ;
1834	// TSTS
1835	if ( spsig.nelements() == tsyssig.nelements() ) {
1836	for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
1837	stsyssig[j] = tsyssig[j+ioffset2] ;
1838	}
1839	for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
1840	stsyssig[j] = tsyssig[j-sigchan+ioffset2] ;
1841	}
1842	tsyssig = stsyssig.copy() ;
1843	for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
1844	stsyssig[j] = doffset2 * tsyssig[j+1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
1845	}
1846	stsyssig[sigchan-1] = doffset2 * tsyssig[0] + ( 1.0 - doffset2 ) * tsyssig[sigchan-1] ;
1847	}
1848	}
1849	else {
1850	// SPECTRA and FLAGTRA
1851	for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
1852	sspsig[j] = spsig[sigchan+ioffset2+j] ;
1853	sflagsig[j] = flagsig[sigchan+ioffset2+j] ;
1854	}
1855	for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
1856	sspsig[j] = spsig[j+ioffset2] ;
1857	sflagsig[j] = flagsig[j+ioffset2] ;
1858	}
1859	spsig = sspsig.copy() ;
1860	flagsig = sflagsig.copy() ;
1861	sspsig[0] = doffset2 * spsig[sigchan-1] + ( 1.0 - doffset2 ) * spsig[0] ;
1862	sflagsig[0] = flagsig[0] + flagsig[sigchan-1] ;
1863	for ( int j = 1 ; j < sigchan ; j++ ) {
1864	sspsig[j] = doffset2 * spsig[j-1] + ( 1.0 - doffset2 ) * spsig[j] ;
1865	sflagsig[j] = flagsig[j-1] + flagsig[j] ;
1866	}
1867	// TSYS
1868	if ( spsig.nelements() == tsyssig.nelements() ) {
1869	for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
1870	stsyssig[j] = tsyssig[sigchan+ioffset2+j] ;
1871	}
1872	for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
1873	stsyssig[j] = tsyssig[j+ioffset2] ;
1874	}
1875	tsyssig = stsyssig.copy() ;
1876	stsyssig[0] = doffset2 * tsyssig[sigchan-1] + ( 1.0 - doffset2 ) * tsyssig[0] ;
1877	for ( int j = 1 ; j < sigchan ; j++ ) {
1878	stsyssig[j] = doffset2 * tsyssig[j-1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
1879	}
1880	}
1881	}
1882	}
1883
1884	// folding
1885	acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ;
1886	acc.add( sspref, !sflagref, stsysref, intref, timeref ) ;
1887
1888	// put result
1889	specColOut.put( i, acc.getSpectrum() ) ;
1890	const Vector<Bool> &msk = acc.getMask() ;
1891	Vector<uChar> flg( msk.shape() ) ;
1892	convertArray( flg, !msk ) ;
1893	flagColOut.put( i, flg ) ;
1894	tsysColOut.put( i, acc.getTsys() ) ;
1895	intervalColOut.put( i, acc.getInterval() ) ;
1896	mjdColOut.put( i, acc.getTime() ) ;
1897	// change FREQ_ID to unshifted IF setting (only for APEX?)
1898	if ( choffset2 != 0.0 ) {
1899	uInt freqid = fidColOut( 0 ) ; // assume single-IF data
1900	double refpix, refval, increment ;
1901	out->frequencies().getEntry( refpix, refval, increment, freqid ) ;
1902	refval -= choffset * increment ;
1903	uInt newfreqid = out->frequencies().addEntry( refpix, refval, increment ) ;
1904	Vector<uInt> freqids = fidColOut.getColumn() ;
1905	for ( uInt j = 0 ; j < freqids.nelements() ; j++ ) {
1906	if ( freqids[j] == freqid )
1907	freqids[j] = newfreqid ;
1908	}
1909	fidColOut.putColumn( freqids ) ;
1910	}
1911
1912	acc.reset() ;
1913	}
1914
1915	return out ;
1916	}
1917
1918
1919	CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
1920	{
1921	// make copy or reference
1922	CountedPtr< Scantable > out = getScantable(in, false);
1923	Table& tout = out->table();
1924	Block<String> cols(4);
1925	cols[0] = String("SCANNO");
1926	cols[1] = String("CYCLENO");
1927	cols[2] = String("BEAMNO");
1928	cols[3] = String("POLNO");
1929	TableIterator iter(tout, cols);
1930	while (!iter.pastEnd()) {
1931	Table subt = iter.table();
1932	// this should leave us with two rows for the two IFs....if not ignore
1933	if (subt.nrow() != 2 ) {
1934	continue;
1935	}
1936	ArrayColumn<Float> specCol(subt, "SPECTRA");
1937	ArrayColumn<Float> tsysCol(subt, "TSYS");
1938	ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1939	Vector<Float> onspec,offspec, ontsys, offtsys;
1940	Vector<uChar> onflag, offflag;
1941	tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
1942	specCol.get(0, onspec); specCol.get(1, offspec);
1943	flagCol.get(0, onflag); flagCol.get(1, offflag);
1944	MaskedArray<Float> on = maskedArray(onspec, onflag);
1945	MaskedArray<Float> off = maskedArray(offspec, offflag);
1946	MaskedArray<Float> oncopy = on.copy();
1947
1948	on /= off; on -= 1.0f;
1949	on *= ontsys[0];
1950	off /= oncopy; off -= 1.0f;
1951	off *= offtsys[0];
1952	specCol.put(0, on.getArray());
1953	const Vector<Bool>& m0 = on.getMask();
1954	Vector<uChar> flags0(m0.shape());
1955	convertArray(flags0, !m0);
1956	flagCol.put(0, flags0);
1957
1958	specCol.put(1, off.getArray());
1959	const Vector<Bool>& m1 = off.getMask();
1960	Vector<uChar> flags1(m1.shape());
1961	convertArray(flags1, !m1);
1962	flagCol.put(1, flags1);
1963	++iter;
1964	}
1965
1966	return out;
1967	}
1968
1969	std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1970	const std::vector< bool > & mask,
1971	const std::string& which )
1972	{
1973
1974	Vector<Bool> m(mask);
1975	const Table& tab = in->table();
1976	ROArrayColumn<Float> specCol(tab, "SPECTRA");
1977	ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1978	std::vector<float> out;
1979	for (uInt i=0; i < tab.nrow(); ++i ) {
1980	Vector<Float> spec; specCol.get(i, spec);
1981	Vector<uChar> flag; flagCol.get(i, flag);
1982	MaskedArray<Float> ma = maskedArray(spec, flag);
1983	float outstat = 0.0;
1984	if ( spec.nelements() == m.nelements() ) {
1985	outstat = mathutil::statistics(which, ma(m));
1986	} else {
1987	outstat = mathutil::statistics(which, ma);
1988	}
1989	out.push_back(outstat);
1990	}
1991	return out;
1992	}
1993
1994	std::vector< float > STMath::statisticRow( const CountedPtr< Scantable > & in,
1995	const std::vector< bool > & mask,
1996	const std::string& which,
1997	int row )
1998	{
1999
2000	Vector<Bool> m(mask);
2001	const Table& tab = in->table();
2002	ROArrayColumn<Float> specCol(tab, "SPECTRA");
2003	ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2004	std::vector<float> out;
2005
2006	Vector<Float> spec; specCol.get(row, spec);
2007	Vector<uChar> flag; flagCol.get(row, flag);
2008	MaskedArray<Float> ma = maskedArray(spec, flag);
2009	float outstat = 0.0;
2010	if ( spec.nelements() == m.nelements() ) {
2011	outstat = mathutil::statistics(which, ma(m));
2012	} else {
2013	outstat = mathutil::statistics(which, ma);
2014	}
2015	out.push_back(outstat);
2016
2017	return out;
2018	}
2019
2020	std::vector< int > STMath::minMaxChan( const CountedPtr< Scantable > & in,
2021	const std::vector< bool > & mask,
2022	const std::string& which )
2023	{
2024
2025	Vector<Bool> m(mask);
2026	const Table& tab = in->table();
2027	ROArrayColumn<Float> specCol(tab, "SPECTRA");
2028	ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2029	std::vector<int> out;
2030	for (uInt i=0; i < tab.nrow(); ++i ) {
2031	Vector<Float> spec; specCol.get(i, spec);
2032	Vector<uChar> flag; flagCol.get(i, flag);
2033	MaskedArray<Float> ma = maskedArray(spec, flag);
2034	if (ma.ndim() != 1) {
2035	throw (ArrayError(
2036	"std::vector<int> STMath::minMaxChan("
2037	"ContedPtr<Scantable> &in, std::vector<bool> &mask, "
2038	" std::string &which)"
2039	" - MaskedArray is not 1D"));
2040	}
2041	IPosition outpos(1,0);
2042	if ( spec.nelements() == m.nelements() ) {
2043	outpos = mathutil::minMaxPos(which, ma(m));
2044	} else {
2045	outpos = mathutil::minMaxPos(which, ma);
2046	}
2047	out.push_back(outpos[0]);
2048	}
2049	return out;
2050	}
2051
2052	CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
2053	int width )
2054	{
2055	if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
2056	CountedPtr< Scantable > out = getScantable(in, false);
2057	Table& tout = out->table();
2058	out->frequencies().rescale(width, "BIN");
2059	ArrayColumn<Float> specCol(tout, "SPECTRA");
2060	ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
2061	ArrayColumn<Float> tsysCol(tout, "TSYS");
2062
2063	for (uInt i=0; i < tout.nrow(); ++i ) {
2064	MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
2065	MaskedArray<Float> maout;
2066	LatticeUtilities::bin(maout, main, 0, Int(width));
2067	specCol.put(i, maout.getArray());
2068	flagCol.put(i, flagsFromMA(maout));
2069	if (tsysCol(i).nelements() == specCol(i).nelements()) {
2070	MaskedArray<Float> matsysin = maskedArray(tsysCol(i), flagCol(i));
2071	MaskedArray<Float> matsysout;
2072	LatticeUtilities::bin(matsysout, matsysin, 0, Int(width));
2073	tsysCol.put(i, matsysout.getArray());
2074	}
2075	// take only the first binned spectrum's length for the deprecated
2076	// global header item nChan
2077	if (i==0) tout.rwKeywordSet().define(String("nChan"),
2078	Int(maout.getArray().nelements()));
2079	}
2080	return out;
2081	}
2082
2083	CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
2084	const std::string& method,
2085	float width )
2086	//
2087	// Should add the possibility of width being specified in km/s. This means
2088	// that for each freqID (SpectralCoordinate) we will need to convert to an
2089	// average channel width (say at the reference pixel). Then we would need
2090	// to be careful to make sure each spectrum (of different freqID)
2091	// is the same length.
2092	//
2093	{
2094	//InterpolateArray1D<Double,Float>::InterpolationMethod interp;
2095	Int interpMethod(stringToIMethod(method));
2096
2097	CountedPtr< Scantable > out = getScantable(in, false);
2098	Table& tout = out->table();
2099
2100	// Resample SpectralCoordinates (one per freqID)
2101	out->frequencies().rescale(width, "RESAMPLE");
2102	TableIterator iter(tout, "IFNO");
2103	TableRow row(tout);
2104	while ( !iter.pastEnd() ) {
2105	Table tab = iter.table();
2106	ArrayColumn<Float> specCol(tab, "SPECTRA");
2107	//ArrayColumn<Float> tsysCol(tout, "TSYS");
2108	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2109	Vector<Float> spec;
2110	Vector<uChar> flag;
2111	specCol.get(0,spec); // the number of channels should be constant per IF
2112	uInt nChanIn = spec.nelements();
2113	Vector<Float> xIn(nChanIn); indgen(xIn);
2114	Int fac = Int(nChanIn/width);
2115	Vector<Float> xOut(fac+10); // 10 to be safe - resize later
2116	uInt k = 0;
2117	Float x = 0.0;
2118	while (x < Float(nChanIn) ) {
2119	xOut(k) = x;
2120	k++;
2121	x += width;
2122	}
2123	uInt nChanOut = k;
2124	xOut.resize(nChanOut, True);
2125	// process all rows for this IFNO
2126	Vector<Float> specOut;
2127	Vector<Bool> maskOut;
2128	Vector<uChar> flagOut;
2129	for (uInt i=0; i < tab.nrow(); ++i) {
2130	specCol.get(i, spec);
2131	flagCol.get(i, flag);
2132	Vector<Bool> mask(flag.nelements());
2133	convertArray(mask, flag);
2134
2135	IPosition shapeIn(spec.shape());
2136	//sh.nchan = nChanOut;
2137	InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
2138	xIn, spec, mask,
2139	interpMethod, True, True);
2140	/// @todo do the same for channel based Tsys
2141	flagOut.resize(maskOut.nelements());
2142	convertArray(flagOut, maskOut);
2143	specCol.put(i, specOut);
2144	flagCol.put(i, flagOut);
2145	}
2146	++iter;
2147	}
2148
2149	return out;
2150	}
2151
2152	STMath::imethod STMath::stringToIMethod(const std::string& in)
2153	{
2154	static STMath::imap lookup;
2155
2156	// initialize the lookup table if necessary
2157	if ( lookup.empty() ) {
2158	lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
2159	lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
2160	lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
2161	lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
2162	}
2163
2164	STMath::imap::const_iterator iter = lookup.find(in);
2165
2166	if ( lookup.end() == iter ) {
2167	std::string message = in;
2168	message += " is not a valid interpolation mode";
2169	throw(AipsError(message));
2170	}
2171	return iter->second;
2172	}
2173
2174	WeightType STMath::stringToWeight(const std::string& in)
2175	{
2176	static std::map<std::string, WeightType> lookup;
2177
2178	// initialize the lookup table if necessary
2179	if ( lookup.empty() ) {
2180	lookup["NONE"] = asap::W_NONE;
2181	lookup["TINT"] = asap::W_TINT;
2182	lookup["TINTSYS"] = asap::W_TINTSYS;
2183	lookup["TSYS"] = asap::W_TSYS;
2184	lookup["VAR"] = asap::W_VAR;
2185	}
2186
2187	std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
2188
2189	if ( lookup.end() == iter ) {
2190	std::string message = in;
2191	message += " is not a valid weighting mode";
2192	throw(AipsError(message));
2193	}
2194	return iter->second;
2195	}
2196
2197	CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
2198	const vector< float > & coeff,
2199	const std::string & filename,
2200	const std::string& method)
2201	{
2202	LogIO os( LogOrigin( "STMath", "gainElevation", WHERE ) ) ;
2203	// Get elevation data from Scantable and convert to degrees
2204	CountedPtr< Scantable > out = getScantable(in, false);
2205	Table& tab = out->table();
2206	ROScalarColumn<Float> elev(tab, "ELEVATION");
2207	Vector<Float> x = elev.getColumn();
2208	x *= Float(180 / C::pi); // Degrees
2209
2210	Vector<Float> coeffs(coeff);
2211	const uInt nc = coeffs.nelements();
2212	if ( filename.length() > 0 && nc > 0 ) {
2213	throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
2214	}
2215
2216	// Correct
2217	if ( nc > 0 \|\| filename.length() == 0 ) {
2218	// Find instrument
2219	Bool throwit = True;
2220	Instrument inst =
2221	STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2222	throwit);
2223
2224	// Set polynomial
2225	Polynomial<Float>* ppoly = 0;
2226	Vector<Float> coeff;
2227	String msg;
2228	if ( nc > 0 ) {
2229	ppoly = new Polynomial<Float>(nc-1);
2230	coeff = coeffs;
2231	msg = String("user");
2232	} else {
2233	STAttr sdAttr;
2234	coeff = sdAttr.gainElevationPoly(inst);
2235	ppoly = new Polynomial<Float>(coeff.nelements()-1);
2236	msg = String("built in");
2237	}
2238
2239	if ( coeff.nelements() > 0 ) {
2240	ppoly->setCoefficients(coeff);
2241	} else {
2242	delete ppoly;
2243	throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
2244	}
2245	os << "Making polynomial correction with " << msg << " coefficients:" << endl;
2246	os << " " << coeff << LogIO::POST;
2247	const uInt nrow = tab.nrow();
2248	Vector<Float> factor(nrow);
2249	for ( uInt i=0; i < nrow; ++i ) {
2250	factor[i] = 1.0 / (*ppoly)(x[i]);
2251	}
2252	delete ppoly;
2253	scaleByVector(tab, factor, true);
2254
2255	} else {
2256	// Read and correct
2257	os << "Making correction from ascii Table" << LogIO::POST;
2258	scaleFromAsciiTable(tab, filename, method, x, true);
2259	}
2260	return out;
2261	}
2262
2263	void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
2264	const std::string& method,
2265	const Vector<Float>& xout, bool dotsys)
2266	{
2267
2268	// Read gain-elevation ascii file data into a Table.
2269
2270	String formatString;
2271	Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
2272	scaleFromTable(in, tbl, method, xout, dotsys);
2273	}
2274
2275	void STMath::scaleFromTable(Table& in,
2276	const Table& table,
2277	const std::string& method,
2278	const Vector<Float>& xout, bool dotsys)
2279	{
2280
2281	ROScalarColumn<Float> geElCol(table, "ELEVATION");
2282	ROScalarColumn<Float> geFacCol(table, "FACTOR");
2283	Vector<Float> xin = geElCol.getColumn();
2284	Vector<Float> yin = geFacCol.getColumn();
2285	Vector<Bool> maskin(xin.nelements(),True);
2286
2287	// Interpolate (and extrapolate) with desired method
2288
2289	InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2290
2291	Vector<Float> yout;
2292	Vector<Bool> maskout;
2293	InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
2294	xin, yin, maskin, interp,
2295	True, True);
2296
2297	scaleByVector(in, Float(1.0)/yout, dotsys);
2298	}
2299
2300	void STMath::scaleByVector( Table& in,
2301	const Vector< Float >& factor,
2302	bool dotsys )
2303	{
2304	uInt nrow = in.nrow();
2305	if ( factor.nelements() != nrow ) {
2306	throw(AipsError("factors.nelements() != table.nelements()"));
2307	}
2308	ArrayColumn<Float> specCol(in, "SPECTRA");
2309	ArrayColumn<uChar> flagCol(in, "FLAGTRA");
2310	ArrayColumn<Float> tsysCol(in, "TSYS");
2311	for (uInt i=0; i < nrow; ++i) {
2312	MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2313	ma *= factor[i];
2314	specCol.put(i, ma.getArray());
2315	flagCol.put(i, flagsFromMA(ma));
2316	if ( dotsys ) {
2317	Vector<Float> tsys = tsysCol(i);
2318	tsys *= factor[i];
2319	tsysCol.put(i,tsys);
2320	}
2321	}
2322	}
2323
2324	CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
2325	float d, float etaap,
2326	float jyperk )
2327	{
2328	LogIO os( LogOrigin( "STMath", "convertFlux", WHERE ) ) ;
2329
2330	CountedPtr< Scantable > out = getScantable(in, false);
2331	Table& tab = in->table();
2332	Table& outtab = out->table();
2333	Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
2334	Unit K(String("K"));
2335	Unit JY(String("Jy"));
2336
2337	bool tokelvin = true;
2338	Double cfac = 1.0;
2339
2340	if ( fluxUnit == JY ) {
2341	os << "Converting to K" << LogIO::POST;
2342	Quantum<Double> t(1.0,fluxUnit);
2343	Quantum<Double> t2 = t.get(JY);
2344	cfac = (t2 / t).getValue(); // value to Jy
2345
2346	tokelvin = true;
2347	out->setFluxUnit("K");
2348	} else if ( fluxUnit == K ) {
2349	os << "Converting to Jy" << LogIO::POST;
2350	Quantum<Double> t(1.0,fluxUnit);
2351	Quantum<Double> t2 = t.get(K);
2352	cfac = (t2 / t).getValue(); // value to K
2353
2354	tokelvin = false;
2355	out->setFluxUnit("Jy");
2356	} else {
2357	throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
2358	}
2359	// Make sure input values are converted to either Jy or K first...
2360	Float factor = cfac;
2361
2362	// Select method
2363	if (jyperk > 0.0) {
2364	factor *= jyperk;
2365	if ( tokelvin ) factor = 1.0 / jyperk;
2366	os << "Jy/K = " << jyperk << LogIO::POST;
2367	Vector<Float> factors(outtab.nrow(), factor);
2368	scaleByVector(outtab,factors, false);
2369	} else if ( etaap > 0.0) {
2370	if (d < 0) {
2371	Instrument inst =
2372	STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2373	True);
2374	STAttr sda;
2375	d = sda.diameter(inst);
2376	}
2377	jyperk = STAttr::findJyPerK(etaap, d);
2378	os << "Jy/K = " << jyperk << LogIO::POST;
2379	factor *= jyperk;
2380	if ( tokelvin ) {
2381	factor = 1.0 / factor;
2382	}
2383	Vector<Float> factors(outtab.nrow(), factor);
2384	scaleByVector(outtab, factors, False);
2385	} else {
2386
2387	// OK now we must deal with automatic look up of values.
2388	// We must also deal with the fact that the factors need
2389	// to be computed per IF and may be different and may
2390	// change per integration.
2391
2392	os <<"Looking up conversion factors" << LogIO::POST;
2393	convertBrightnessUnits(out, tokelvin, cfac);
2394	}
2395
2396	return out;
2397	}
2398
2399	void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
2400	bool tokelvin, float cfac )
2401	{
2402	Table& table = in->table();
2403	Instrument inst =
2404	STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
2405	TableIterator iter(table, "FREQ_ID");
2406	STFrequencies stfreqs = in->frequencies();
2407	STAttr sdAtt;
2408	while (!iter.pastEnd()) {
2409	Table tab = iter.table();
2410	ArrayColumn<Float> specCol(tab, "SPECTRA");
2411	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2412	ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
2413	MEpoch::ROScalarColumn timeCol(tab, "TIME");
2414
2415	uInt freqid; freqidCol.get(0, freqid);
2416	Vector<Float> tmpspec; specCol.get(0, tmpspec);
2417	// STAttr.JyPerK has a Vector interface... change sometime.
2418	Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
2419	for ( uInt i=0; i<tab.nrow(); ++i) {
2420	Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
2421	Float factor = cfac * jyperk;
2422	if ( tokelvin ) factor = Float(1.0) / factor;
2423	MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2424	ma *= factor;
2425	specCol.put(i, ma.getArray());
2426	flagCol.put(i, flagsFromMA(ma));
2427	}
2428	++iter;
2429	}
2430	}
2431
2432	CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
2433	const std::vector<float>& tau )
2434	{
2435	CountedPtr< Scantable > out = getScantable(in, false);
2436
2437	Table outtab = out->table();
2438
2439	const Int ntau = uInt(tau.size());
2440	std::vector<float>::const_iterator tauit = tau.begin();
2441	AlwaysAssert((ntau == 1 \|\| ntau == in->nif() \|\| ntau == in->nif() * in->npol()),
2442	AipsError);
2443	TableIterator iiter(outtab, "IFNO");
2444	while ( !iiter.pastEnd() ) {
2445	Table itab = iiter.table();
2446	TableIterator piter(itab, "POLNO");
2447	while ( !piter.pastEnd() ) {
2448	Table tab = piter.table();
2449	ROScalarColumn<Float> elev(tab, "ELEVATION");
2450	ArrayColumn<Float> specCol(tab, "SPECTRA");
2451	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2452	ArrayColumn<Float> tsysCol(tab, "TSYS");
2453	for ( uInt i=0; i<tab.nrow(); ++i) {
2454	Float zdist = Float(C::pi_2) - elev(i);
2455	Float factor = exp(*tauit/cos(zdist));
2456	MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2457	ma *= factor;
2458	specCol.put(i, ma.getArray());
2459	flagCol.put(i, flagsFromMA(ma));
2460	Vector<Float> tsys;
2461	tsysCol.get(i, tsys);
2462	tsys *= factor;
2463	tsysCol.put(i, tsys);
2464	}
2465	if (ntau == in->nif()*in->npol() ) {
2466	tauit++;
2467	}
2468	piter++;
2469	}
2470	if (ntau >= in->nif() ) {
2471	tauit++;
2472	}
2473	iiter++;
2474	}
2475	return out;
2476	}
2477
2478	CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
2479	const std::string& kernel,
2480	float width, int order)
2481	{
2482	CountedPtr< Scantable > out = getScantable(in, false);
2483	Table table = out->table();
2484
2485	TableIterator iter(table, "IFNO");
2486	while (!iter.pastEnd()) {
2487	Table tab = iter.table();
2488	ArrayColumn<Float> specCol(tab, "SPECTRA");
2489	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2490	ScalarColumn<uInt> flagrowCol(tab, "FLAGROW");
2491	Vector<Float> spec;
2492	Vector<uChar> flag;
2493	Vector<uInt> flagrow = flagrowCol.getColumn();
2494	for (uInt i = 0; i < tab.nrow(); ++i) {
2495	if (flagrow[i] != 0) {
2496	// do not process flagged row
2497	continue;
2498	}
2499	specCol.get(i, spec);
2500	flagCol.get(i, flag);
2501	Vector<Bool> mask(flag.nelements());
2502	convertArray(mask, flag);
2503	Vector<Float> specout;
2504	Vector<Bool> maskout;
2505	if (kernel == "hanning") {
2506	mathutil::hanning(specout, maskout, spec, !mask);
2507	} else if (kernel == "rmedian") {
2508	mathutil::runningMedian(specout, maskout, spec , mask, width);
2509	} else if (kernel == "poly") {
2510	mathutil::polyfit(specout, maskout, spec, !mask, width, order);
2511	}
2512
2513	for (uInt j = 0; j < flag.nelements(); ++j) {
2514	uChar userFlag = 1 << 7;
2515	if (maskout[j]==True) userFlag = 0 << 7;
2516	flag(j) = userFlag;
2517	}
2518
2519	flagCol.put(i, flag);
2520	specCol.put(i, specout);
2521	}
2522	++iter;
2523	}
2524	return out;
2525	}
2526
2527	CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
2528	const std::string& kernel, float width,
2529	int order)
2530	{
2531	if (kernel == "rmedian" \|\| kernel == "hanning" \|\| kernel == "poly") {
2532	return smoothOther(in, kernel, width, order);
2533	}
2534	CountedPtr< Scantable > out = getScantable(in, false);
2535	Table& table = out->table();
2536	VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
2537	// same IFNO should have same no of channels
2538	// this saves overhead
2539	TableIterator iter(table, "IFNO");
2540	while (!iter.pastEnd()) {
2541	Table tab = iter.table();
2542	ArrayColumn<Float> specCol(tab, "SPECTRA");
2543	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2544	ScalarColumn<uInt> flagrowCol(tab, "FLAGROW");
2545	Vector<Float> spec = specCol( 0 );
2546	uInt nchan = spec.nelements();
2547	Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
2548	Convolver<Float> conv(kvec, IPosition(1,nchan));
2549	Vector<uChar> flag;
2550	Vector<Bool> mask(nchan);
2551	Vector<uInt> flagrow = flagrowCol.getColumn();
2552	for ( uInt i=0; i<tab.nrow(); ++i) {
2553	if (flagrow[i] != 0) {
2554	// do not process flagged row
2555	continue;
2556	}
2557
2558	specCol.get(i, spec);
2559	flagCol.get(i, flag);
2560	convertArray(mask, flag);
2561	Vector<Float> specout;
2562	//mathutil::replaceMaskByZero(specout, mask);
2563	mathutil::replaceMaskByZero(spec, !mask);
2564	//std::vector<bool> vmask;
2565	//(!mask).tovector(vmask);
2566	//mathutil::doZeroOrderInterpolation(spec, vmask);
2567	conv.linearConv(specout, spec);
2568	specCol.put(i, specout);
2569	}
2570	++iter;
2571	}
2572	return out;
2573	}
2574
2575	CountedPtr< Scantable >
2576	STMath::merge( const std::vector< CountedPtr < Scantable > >& in,
2577	const std::string &freqTol )
2578	{
2579	Double freqTolInHz = 0.0; // default is 0.0Hz (merge only when exact match)
2580	if (freqTol.size() > 0) {
2581	Quantum<Double> freqTolInQuantity;
2582	if (!Quantum<Double>::read(freqTolInQuantity, freqTol)) {
2583	throw(AipsError("Failed to convert freqTol string to quantity"));
2584	}
2585	if (!freqTolInQuantity.isConform("Hz")) {
2586	throw(AipsError("Invalid freqTol string"));
2587	}
2588	freqTolInHz = freqTolInQuantity.getValue("Hz");
2589	LogIO os(LogOrigin("STMath", "merge", WHERE));
2590	os << "frequency tolerance = " << freqTolInHz << "Hz" << LogIO::POST;
2591	}
2592
2593	if ( in.size() < 2 ) {
2594	throw(AipsError("Need at least two scantables to perform a merge."));
2595	}
2596	std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
2597	bool insitu = insitu_;
2598	setInsitu(false);
2599	CountedPtr< Scantable > out = getScantable(*it, false);
2600	setInsitu(insitu);
2601	Table& tout = out->table();
2602	ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
2603	ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
2604	ScalarColumn<uInt> ifnocol(tout, "IFNO");
2605	// Renumber SCANNO to be 0-based
2606	Vector<uInt> scannos = scannocol.getColumn();
2607	uInt offset = min(scannos);
2608	scannos -= offset;
2609	scannocol.putColumn(scannos);
2610	uInt newscanno = max(scannos)+1;
2611	++it;
2612
2613	// new IFNO
2614	uInt ifnoCounter = max(ifnocol.getColumn()) + 1;
2615
2616	// Here we assume that each IFNO has unique MOLECULE_ID
2617	// molIdMap:
2618	// KEY: IFNO
2619	// VALUE: MOLECULE_ID
2620	map<uInt, uInt> molIdMap;
2621	{
2622	TableIterator ifit(tout, "IFNO");
2623	while (!ifit.pastEnd()) {
2624	ROTableRow row(ifit.table());
2625	const TableRecord& rec = row.get(0);
2626	molIdMap[rec.asuInt("IFNO")] = rec.asuInt("MOLECULE_ID");
2627	ifit.next();
2628	}
2629	}
2630
2631	while ( it != in.end() ){
2632	// Check FREQUENCIES/BASEFRAME
2633	if ( out->frequencies().getFrame(true) != (*it)->frequencies().getFrame(true) ) {
2634	throw(AipsError("BASEFRAME is not identical"));
2635	}
2636
2637	if ( ! (it)->conformant(out) ) {
2638	// non conformant.
2639	LogIO os( LogOrigin( "STMath", "merge()", WHERE ) ) ;
2640	os << LogIO::SEVERE << "Can't merge scantables as header informations (any one of AntennaName, Equinox, and FluxUnit) differ." << LogIO::EXCEPTION ;
2641	}
2642	out->appendToHistoryTable((*it)->history());
2643	const Table& tab = (*it)->table();
2644
2645	Block<String> cols(3);
2646	cols[0] = String("FREQ_ID");
2647	cols[1] = String("MOLECULE_ID");
2648	cols[2] = String("FOCUS_ID");
2649
2650	TableIterator scanit(tab, "SCANNO");
2651	while (!scanit.pastEnd()) {
2652	TableIterator subit(scanit.table(), cols);
2653	while ( !subit.pastEnd() ) {
2654	uInt nrow = tout.nrow();
2655	Table thetab = subit.table();
2656	ROTableRow row(thetab);
2657	Vector<uInt> thecolvals(thetab.nrow());
2658	// The selected subset of table should have
2659	// the equal FREQ_ID, MOLECULE_ID, and FOCUS_ID values.
2660	const TableRecord& rec = row.get(0);
2661
2662	// Set the proper FREQ_ID
2663	Double rv,rp,inc;
2664	(*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
2665	uInt id;
2666
2667	// default value is new unique IFNO
2668	uInt newifno = ifnoCounter;
2669	Bool isIfMerged = False;
2670	uInt nchan = rec.asArrayFloat("SPECTRA").shape()[0];
2671	//id = out->frequencies().addEntry(rp, rv, inc);
2672	if ( !out->frequencies().match(rp, rv, inc, freqTolInHz, id) ) {
2673	// add new entry to FREQUENCIES table
2674	id = out->frequencies().addEntry(rp, rv, inc);
2675
2676	// increment counter for IFNO
2677	ifnoCounter++;
2678	}
2679	else {
2680	// should renumber IFNO to be same as existing rows that have same FREQ_ID
2681	LogIO os(LogOrigin("STMath", "merge", WHERE));
2682	Table outFreqIdSelected = tout(tout.col("FREQ_ID") == id);
2683	TableIterator _iter(outFreqIdSelected, "IFNO");
2684	map<uInt, uInt> nchanMap;
2685	while (!_iter.pastEnd()) {
2686	const Table _table = _iter.table();
2687	ROTableRow _row(_table);
2688	const TableRecord &_rec = _row.get(0);
2689	uInt nchan = _rec.asArrayFloat("SPECTRA").shape()[0];
2690	if (nchanMap.find(nchan) != nchanMap.end()) {
2691	throw(AipsError("There are non-unique IFNOs assigned to spectra that have same FREQ_ID and same nchan. Something wrong."));
2692	}
2693	nchanMap[nchan] = _rec.asuInt("IFNO");
2694	_iter.next();
2695	}
2696
2697	os << LogIO::DEBUGGING << "nchanMap for " << id << ":" << LogIO::POST;
2698	for (map<uInt, uInt>::iterator i = nchanMap.begin(); i != nchanMap.end(); ++i) {
2699	os << LogIO::DEBUGGING << "nchanMap[" << i->first << "] = " << i->second << LogIO::POST;
2700	}
2701
2702	if (nchanMap.find(nchan) == nchanMap.end()) {
2703	// increment counter for IFNO
2704	ifnoCounter++;
2705	}
2706	else {
2707	// renumber IFNO to be same as existing value that corresponds to nchan
2708	newifno = nchanMap[nchan];
2709	isIfMerged = True;
2710	}
2711	os << LogIO::DEBUGGING << "newifno = " << newifno << LogIO::POST;
2712	}
2713	// copy rows to output table
2714	tout.addRow(thetab.nrow());
2715	TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
2716
2717	Slicer slice(IPosition(1, nrow), IPosition(1, thetab.nrow()),
2718	Slicer::endIsLength);
2719
2720	thecolvals = id;
2721	freqidcol.putColumnRange(slice, thecolvals);
2722
2723	thecolvals = newifno;
2724	ifnocol.putColumnRange(slice, thecolvals);
2725
2726	// Set the proper MOLECULE_ID
2727	Vector<String> name,fname;Vector<Double> rf;
2728	(*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
2729	id = out->molecules().addEntry(rf, name, fname);
2730	if (molIdMap.find(newifno) == molIdMap.end()) {
2731	// add new entry to molIdMap
2732	molIdMap[newifno] = id;
2733	}
2734	if (isIfMerged) {
2735	id = molIdMap[newifno];
2736	}
2737	thecolvals = id;
2738	molidcol.putColumnRange(slice, thecolvals);
2739
2740	// Set the proper FOCUS_ID
2741	Float fpa,frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
2742	(*it)->focus().getEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2743	fxy, fxyp, rec.asuInt("FOCUS_ID"));
2744	id = out->focus().addEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2745	fxy, fxyp);
2746	thecolvals = id;
2747	focusidcol.putColumnRange(slice, thecolvals);
2748
2749	// Set the proper SCANNO
2750	thecolvals = newscanno;
2751	scannocol.putColumnRange(slice, thecolvals);
2752
2753	++subit;
2754	}
2755	++newscanno;
2756	++scanit;
2757	}
2758	++it;
2759	}
2760	return out;
2761	}
2762
2763	CountedPtr< Scantable >
2764	STMath::invertPhase( const CountedPtr < Scantable >& in )
2765	{
2766	return applyToPol(in, &STPol::invertPhase, Float(0.0));
2767	}
2768
2769	CountedPtr< Scantable >
2770	STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
2771	{
2772	return applyToPol(in, &STPol::rotatePhase, Float(phase));
2773	}
2774
2775	CountedPtr< Scantable >
2776	STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
2777	{
2778	return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
2779	}
2780
2781	CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
2782	STPol::polOperation fptr,
2783	Float phase )
2784	{
2785	CountedPtr< Scantable > out = getScantable(in, false);
2786	Table& tout = out->table();
2787	Block<String> cols(4);
2788	cols[0] = String("SCANNO");
2789	cols[1] = String("BEAMNO");
2790	cols[2] = String("IFNO");
2791	cols[3] = String("CYCLENO");
2792	TableIterator iter(tout, cols);
2793	CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
2794	out->getPolType() );
2795	while (!iter.pastEnd()) {
2796	Table t = iter.table();
2797	ArrayColumn<Float> speccol(t, "SPECTRA");
2798	ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
2799	Matrix<Float> pols(speccol.getColumn());
2800	try {
2801	stpol->setSpectra(pols);
2802	Float fang,fhand;
2803	fang = in->focusTable_.getTotalAngle(focidcol(0));
2804	fhand = in->focusTable_.getFeedHand(focidcol(0));
2805	stpol->setPhaseCorrections(fang, fhand);
2806	// use a member function pointer in STPol. This only works on
2807	// the STPol pointer itself, not the Counted Pointer so
2808	// derefernce it.
2809	(&((stpol))->fptr)(phase);
2810	speccol.putColumn(stpol->getSpectra());
2811	} catch (AipsError& e) {
2812	//delete stpol;stpol=0;
2813	throw(e);
2814	}
2815	++iter;
2816	}
2817	//delete stpol;stpol=0;
2818	return out;
2819	}
2820
2821	CountedPtr< Scantable >
2822	STMath::swapPolarisations( const CountedPtr< Scantable > & in )
2823	{
2824	CountedPtr< Scantable > out = getScantable(in, false);
2825	Table& tout = out->table();
2826	Table t0 = tout(tout.col("POLNO") == 0);
2827	Table t1 = tout(tout.col("POLNO") == 1);
2828	if ( t0.nrow() != t1.nrow() )
2829	throw(AipsError("Inconsistent number of polarisations"));
2830	ArrayColumn<Float> speccol0(t0, "SPECTRA");
2831	ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
2832	ArrayColumn<Float> speccol1(t1, "SPECTRA");
2833	ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
2834	Matrix<Float> s0 = speccol0.getColumn();
2835	Matrix<uChar> f0 = flagcol0.getColumn();
2836	speccol0.putColumn(speccol1.getColumn());
2837	flagcol0.putColumn(flagcol1.getColumn());
2838	speccol1.putColumn(s0);
2839	flagcol1.putColumn(f0);
2840	return out;
2841	}
2842
2843	CountedPtr< Scantable >
2844	STMath::averagePolarisations( const CountedPtr< Scantable > & in,
2845	const std::vector<bool>& mask,
2846	const std::string& weight )
2847	{
2848	if (in->npol() < 2 )
2849	throw(AipsError("averagePolarisations can only be applied to two or more"
2850	"polarisations"));
2851	bool insitu = insitu_;
2852	setInsitu(false);
2853	CountedPtr< Scantable > pols = getScantable(in, true);
2854	setInsitu(insitu);
2855	Table& tout = pols->table();
2856	std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
2857	Table tab = tableCommand(taql, in->table());
2858	if (tab.nrow() == 0 )
2859	throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
2860	TableCopy::copyRows(tout, tab);
2861	TableVector<uInt> vec(tout, "POLNO");
2862	vec = 0;
2863	pols->table_.rwKeywordSet().define("nPol", Int(1));
2864	pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
2865	//pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
2866	std::vector<CountedPtr<Scantable> > vpols;
2867	vpols.push_back(pols);
2868	CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
2869	return out;
2870	}
2871
2872	CountedPtr< Scantable >
2873	STMath::averageBeams( const CountedPtr< Scantable > & in,
2874	const std::vector<bool>& mask,
2875	const std::string& weight )
2876	{
2877	bool insitu = insitu_;
2878	setInsitu(false);
2879	CountedPtr< Scantable > beams = getScantable(in, false);
2880	setInsitu(insitu);
2881	Table& tout = beams->table();
2882	// give all rows the same BEAMNO
2883	TableVector<uInt> vec(tout, "BEAMNO");
2884	vec = 0;
2885	beams->table_.rwKeywordSet().define("nBeam", Int(1));
2886	std::vector<CountedPtr<Scantable> > vbeams;
2887	vbeams.push_back(beams);
2888	CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
2889	return out;
2890	}
2891
2892
2893	CountedPtr< Scantable >
2894	asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
2895	const std::string & refTime,
2896	const std::string & method)
2897	{
2898	LogIO os( casa::LogOrigin("STMath", "frequencyAlign()", WHERE));
2899	// clone as this is not working insitu
2900	bool insitu = insitu_;
2901	setInsitu(false);
2902	CountedPtr< Scantable > out = getScantable(in, false);
2903	setInsitu(insitu);
2904	Table& tout = out->table();
2905	// Get reference Epoch to time of first row or given String
2906	Unit DAY(String("d"));
2907	MEpoch::Ref epochRef(in->getTimeReference());
2908	MEpoch refEpoch;
2909	if (refTime.length()>0) {
2910	Quantum<Double> qt;
2911	if (MVTime::read(qt,refTime)) {
2912	MVEpoch mv(qt);
2913	refEpoch = MEpoch(mv, epochRef);
2914	} else {
2915	throw(AipsError("Invalid format for Epoch string"));
2916	}
2917	} else {
2918	refEpoch = in->timeCol_(0);
2919	}
2920	MPosition refPos = in->getAntennaPosition();
2921
2922	InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2923	/*
2924	// Comment from MV.
2925	// the following code has been commented out because different FREQ_IDs have to be aligned together even
2926	// if the frame doesn't change. So far, lack of this check didn't cause any problems.
2927	// test if user frame is different to base frame
2928	if ( in->frequencies().getFrameString(true)
2929	== in->frequencies().getFrameString(false) ) {
2930	throw(AipsError("Can't convert as no output frame has been set"
2931	" (use set_freqframe) or it is aligned already."));
2932	}
2933	*/
2934	MFrequency::Types system = in->frequencies().getFrame();
2935	MVTime mvt(refEpoch.getValue());
2936	String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
2937	os << "Aligning at reference Epoch " << epochout
2938	<< " in frame " << MFrequency::showType(system) << LogIO::POST;
2939	// set up the iterator
2940	Block<String> cols(4);
2941	// select by constant direction
2942	cols[0] = String("SRCNAME");
2943	cols[1] = String("BEAMNO");
2944	// select by IF ( no of channels varies over this )
2945	cols[2] = String("IFNO");
2946	// select by restfrequency
2947	cols[3] = String("MOLECULE_ID");
2948	TableIterator iter(tout, cols);
2949	while ( !iter.pastEnd() ) {
2950	Table t = iter.table();
2951	ROScalarColumn<String> snCol(t, "SRCNAME");
2952	os << "Aligning to position of source '" << snCol(0) << "'" << LogIO::POST;
2953	MDirection::ROScalarColumn dirCol(t, "DIRECTION");
2954	TableIterator fiter(t, "FREQ_ID");
2955	// determine nchan from the first row. This should work as
2956	// we are iterating over BEAMNO and IFNO
2957	// we should have constant direction
2958
2959	ROArrayColumn<Float> sCol(t, "SPECTRA");
2960	const MDirection direction = dirCol(0);
2961	const uInt nchan = sCol(0).nelements();
2962
2963	// skip operations if there is nothing to align
2964	if (fiter.pastEnd()) {
2965	continue;
2966	}
2967
2968	Table ftab = fiter.table();
2969	// align all frequency ids with respect to the first encountered id
2970	ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
2971	// get the SpectralCoordinate for the freqid, which we are iterating over
2972	SpectralCoordinate sC = \
2973	in->frequencies().getSpectralCoordinate(freqidCol(0));
2974	FrequencyAligner<Float> fa( sC, nchan, refEpoch,
2975	direction, refPos, system );
2976	// realign the SpectralCoordinate and put into the output Scantable
2977	Vector<String> units(1);
2978	units = String("Hz");
2979	Bool linear=True;
2980	SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
2981	sc2.setWorldAxisUnits(units);
2982	const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
2983	sc2.referenceValue()[0],
2984	sc2.increment()[0]);
2985	while ( !fiter.pastEnd() ) {
2986
2987	ftab = fiter.table();
2988	// spectral coordinate for the current FREQ_ID
2989	ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
2990	sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
2991	// create the "global" abcissa for alignment with same FREQ_ID
2992	Vector<Double> abc(nchan);
2993	for (uInt i=0; i<nchan; i++) {
2994	Double w;
2995	sC.toWorld(w,Double(i));
2996	abc[i] = w;
2997	}
2998	TableVector<uInt> tvec(ftab, "FREQ_ID");
2999	// assign new frequency id to all rows
3000	tvec = id;
3001	// cache abcissa for same time stamps, so iterate over those
3002	TableIterator timeiter(ftab, "TIME");
3003	while ( !timeiter.pastEnd() ) {
3004	Table tab = timeiter.table();
3005	ArrayColumn<Float> specCol(tab, "SPECTRA");
3006	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
3007	MEpoch::ROScalarColumn timeCol(tab, "TIME");
3008	// use align abcissa cache after the first row
3009	// these rows should be just be POLNO
3010	bool first = true;
3011	for (int i=0; i<int(tab.nrow()); ++i) {
3012	// input values
3013	Vector<uChar> flag = flagCol(i);
3014	Vector<Bool> mask(flag.shape());
3015	Vector<Float> specOut, spec;
3016	spec = specCol(i);
3017	Vector<Bool> maskOut;Vector<uChar> flagOut;
3018	convertArray(mask, flag);
3019	// alignment
3020	Bool ok = fa.align(specOut, maskOut, abc, spec,
3021	mask, timeCol(i), !first,
3022	interp, False);
3023	(void) ok; // unused stop compiler nagging
3024	// back into scantable
3025	flagOut.resize(maskOut.nelements());
3026	convertArray(flagOut, maskOut);
3027	flagCol.put(i, flagOut);
3028	specCol.put(i, specOut);
3029	// start abcissa caching
3030	first = false;
3031	}
3032	// next timestamp
3033	++timeiter;
3034	}
3035	// next FREQ_ID
3036	++fiter;
3037	}
3038	// next aligner
3039	++iter;
3040	}
3041	// set this afterwards to ensure we are doing insitu correctly.
3042	out->frequencies().setFrame(system, true);
3043	return out;
3044	}
3045
3046	CountedPtr<Scantable>
3047	asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
3048	const std::string & newtype )
3049	{
3050	if (in->npol() != 2 && in->npol() != 4)
3051	throw(AipsError("Can only convert two or four polarisations."));
3052	if ( in->getPolType() == newtype )
3053	throw(AipsError("No need to convert."));
3054	if ( ! in->selector_.empty() )
3055	throw(AipsError("Can only convert whole scantable. Unset the selection."));
3056	bool insitu = insitu_;
3057	setInsitu(false);
3058	CountedPtr< Scantable > out = getScantable(in, true);
3059	setInsitu(insitu);
3060	Table& tout = out->table();
3061	tout.rwKeywordSet().define("POLTYPE", String(newtype));
3062
3063	Block<String> cols(4);
3064	cols[0] = "SCANNO";
3065	cols[1] = "CYCLENO";
3066	cols[2] = "BEAMNO";
3067	cols[3] = "IFNO";
3068	TableIterator it(in->originalTable_, cols);
3069	String basetype = in->getPolType();
3070	STPol* stpol = STPol::getPolClass(in->factories_, basetype);
3071	try {
3072	while ( !it.pastEnd() ) {
3073	Table tab = it.table();
3074	uInt row = tab.rowNumbers()[0];
3075	stpol->setSpectra(in->getPolMatrix(row));
3076	Float fang,fhand;
3077	fang = in->focusTable_.getTotalAngle(in->mfocusidCol_(row));
3078	fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
3079	stpol->setPhaseCorrections(fang, fhand);
3080	Int npolout = 0;
3081	for (uInt i=0; i<tab.nrow(); ++i) {
3082	Vector<Float> outvec = stpol->getSpectrum(i, newtype);
3083	if ( outvec.nelements() > 0 ) {
3084	tout.addRow();
3085	TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
3086	ArrayColumn<Float> sCol(tout,"SPECTRA");
3087	ScalarColumn<uInt> pCol(tout,"POLNO");
3088	sCol.put(tout.nrow()-1 ,outvec);
3089	pCol.put(tout.nrow()-1 ,uInt(npolout));
3090	npolout++;
3091	}
3092	}
3093	tout.rwKeywordSet().define("nPol", npolout);
3094	++it;
3095	}
3096	} catch (AipsError& e) {
3097	delete stpol;
3098	throw(e);
3099	}
3100	delete stpol;
3101	return out;
3102	}
3103
3104	CountedPtr< Scantable >
3105	asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
3106	const std::string & scantype )
3107	{
3108	bool insitu = insitu_;
3109	setInsitu(false);
3110	CountedPtr< Scantable > out = getScantable(in, true);
3111	setInsitu(insitu);
3112	Table& tout = out->table();
3113	std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
3114	if (scantype == "on") {
3115	taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
3116	}
3117	Table tab = tableCommand(taql, in->table());
3118	TableCopy::copyRows(tout, tab);
3119	if (scantype == "on") {
3120	// re-index SCANNO to 0
3121	TableVector<uInt> vec(tout, "SCANNO");
3122	vec = 0;
3123	}
3124	return out;
3125	}
3126
3127	std::vector<float>
3128	asap::STMath::fft( const casa::CountedPtr< Scantable > & in,
3129	const std::vector<int>& whichrow,
3130	bool getRealImag )
3131	{
3132	std::vector<float> res;
3133	Table tab = in->table();
3134	std::vector<bool> mask;
3135
3136	if (whichrow.size() < 1) { // for all rows (by default)
3137	int nrow = int(tab.nrow());
3138	for (int i = 0; i < nrow; ++i) {
3139	res = in->execFFT(i, mask, getRealImag);
3140	}
3141	} else { // for specified rows
3142	for (uInt i = 0; i < whichrow.size(); ++i) {
3143	res = in->execFFT(i, mask, getRealImag);
3144	}
3145	}
3146
3147	return res;
3148	}
3149
3150
3151	CountedPtr<Scantable>
3152	asap::STMath::lagFlag( const CountedPtr<Scantable>& in,
3153	double start, double end,
3154	const std::string& mode )
3155	{
3156	CountedPtr<Scantable> out = getScantable(in, false);
3157	Table& tout = out->table();
3158	TableIterator iter(tout, "FREQ_ID");
3159	FFTServer<Float,Complex> ffts;
3160
3161	while ( !iter.pastEnd() ) {
3162	Table tab = iter.table();
3163	Double rp,rv,inc;
3164	ROTableRow row(tab);
3165	const TableRecord& rec = row.get(0);
3166	uInt freqid = rec.asuInt("FREQ_ID");
3167	out->frequencies().getEntry(rp, rv, inc, freqid);
3168	ArrayColumn<Float> specCol(tab, "SPECTRA");
3169	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
3170
3171	for (int i=0; i<int(tab.nrow()); ++i) {
3172	Vector<Float> spec = specCol(i);
3173	Vector<uChar> flag = flagCol(i);
3174	std::vector<bool> mask;
3175	for (uInt j = 0; j < flag.nelements(); ++j) {
3176	mask.push_back(!(flag[j]>0));
3177	}
3178	mathutil::doZeroOrderInterpolation(spec, mask);
3179
3180	Vector<Complex> lags;
3181	ffts.fft0(lags, spec);
3182
3183	Int lag0(start+0.5);
3184	Int lag1(end+0.5);
3185	if (mode == "frequency") {
3186	lag0 = Int(spec.nelements()*abs(inc)/(start)+0.5);
3187	lag1 = Int(spec.nelements()*abs(inc)/(end)+0.5);
3188	}
3189	Int lstart = max(0, lag0);
3190	Int lend = min(Int(lags.nelements()-1), lag1);
3191	if (lstart == lend) {
3192	lags[lstart] = Complex(0.0);
3193	} else {
3194	if (lstart > lend) {
3195	Int tmp = lend;
3196	lend = lstart;
3197	lstart = tmp;
3198	}
3199	for (int j=lstart; j <=lend ;++j) {
3200	lags[j] = Complex(0.0);
3201	}
3202	}
3203
3204	ffts.fft0(spec, lags);
3205
3206	specCol.put(i, spec);
3207	}
3208	++iter;
3209	}
3210	return out;
3211	}
3212
3213	// Averaging spectra with different channel/resolution
3214	CountedPtr<Scantable>
3215	STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
3216	const bool& compel,
3217	const std::vector<bool>& mask,
3218	const std::string& weight,
3219	const std::string& avmode )
3220	throw ( casa::AipsError )
3221	{
3222	LogIO os( LogOrigin( "STMath", "new_average()", WHERE ) ) ;
3223	if ( avmode == "SCAN" && in.size() != 1 )
3224	throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
3225	"Use merge first."));
3226
3227	CountedPtr<Scantable> out ; // processed result
3228	if ( compel ) {
3229	std::vector< CountedPtr<Scantable> > newin ; // input for average process
3230	uInt insize = in.size() ; // number of input scantables
3231
3232	// setup newin
3233	bool oldInsitu = insitu_ ;
3234	setInsitu( false ) ;
3235	newin.resize( insize ) ;
3236	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3237	newin[itable] = getScantable( in[itable], false ) ;
3238	}
3239	setInsitu( oldInsitu ) ;
3240
3241	// warning
3242	os << "Average spectra with different spectral resolution" << LogIO::POST ;
3243
3244	// temporarily set coordinfo
3245	vector<string> oldinfo( insize ) ;
3246	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3247	vector<string> coordinfo = in[itable]->getCoordInfo() ;
3248	oldinfo[itable] = coordinfo[0] ;
3249	coordinfo[0] = "Hz" ;
3250	newin[itable]->setCoordInfo( coordinfo ) ;
3251	}
3252
3253	ostringstream oss ;
3254
3255	// check IF frequency coverage
3256	// freqid: list of FREQ_ID, which is used, in each table
3257	// iffreq: list of minimum and maximum frequency for each FREQ_ID in
3258	// each table
3259	// freqid[insize][numIF]
3260	// freqid: [[id00, id01, ...],
3261	// [id10, id11, ...],
3262	// ...
3263	// [idn0, idn1, ...]]
3264	// iffreq[insize][numIF*2]
3265	// iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
3266	// [min_id10, max_id10, min_id11, max_id11, ...],
3267	// ...
3268	// [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
3269	//os << "Check IF settings in each table" << LogIO::POST ;
3270	vector< vector<uInt> > freqid( insize );
3271	vector< vector<double> > iffreq( insize ) ;
3272	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3273	Vector<uInt> freqIds = newin[itable]->mfreqidCol_.getColumn() ;
3274	vector<uInt> uniqueFreqId = newin[itable]->getNumbers(newin[itable]->mfreqidCol_) ;
3275	for ( vector<uInt>::iterator i = uniqueFreqId.begin() ;
3276	i != uniqueFreqId.end() ; i++ ) {
3277	//os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
3278	uInt target = 0 ;
3279	while ( freqIds[target] != *i )
3280	target++ ;
3281	vector<double> abcissa = newin[itable]->getAbcissa( target ) ;
3282	freqid[itable].push_back( *i ) ;
3283	double incr = abs( abcissa[1] - abcissa[0] ) ;
3284	iffreq[itable].push_back( (min_element(abcissa.begin(),abcissa.end()))-0.5incr ) ;
3285	iffreq[itable].push_back( (max_element(abcissa.begin(),abcissa.end()))+0.5incr ) ;
3286	}
3287	}
3288
3289	// debug
3290	// os << "IF settings summary:" << endl ;
3291	// for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
3292	// os << " Table" << i << endl ;
3293	// for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
3294	// os << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2j] << "," << iffreq[i][2j+1] << ")" << endl ;
3295	// }
3296	// }
3297	// os << endl ;
3298	// os.post() ;
3299
3300	// IF grouping based on their frequency coverage
3301	// ifgrp: number of member in each IF group
3302	// ifgrp[numgrp]
3303	// ifgrp: [n0, n1, ...]
3304	//os << "IF grouping based on their frequency coverage" << LogIO::POST ;
3305
3306	// parameter for IF grouping
3307	// groupmode = OR retrieve all region
3308	// AND only retrieve overlaped region
3309	//string groupmode = "AND" ;
3310	string groupmode = "OR" ;
3311	uInt sizecr = 0 ;
3312	if ( groupmode == "AND" )
3313	sizecr = 1 ;
3314	else if ( groupmode == "OR" )
3315	sizecr = 0 ;
3316
3317	vector<double> sortedfreq ;
3318	for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
3319	for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
3320	if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
3321	sortedfreq.push_back( iffreq[i][j] ) ;
3322	}
3323	}
3324	sort( sortedfreq.begin(), sortedfreq.end() ) ;
3325	vector<uInt> ifgrp( sortedfreq.size()-1, 0 ) ;
3326	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3327	for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
3328	double range0 = iffreq[itable][2*iif] ;
3329	double range1 = iffreq[itable][2*iif+1] ;
3330	for ( uInt j = 0 ; j < sortedfreq.size()-1 ; j++ ) {
3331	double fmin = max( range0, sortedfreq[j] ) ;
3332	double fmax = min( range1, sortedfreq[j+1] ) ;
3333	if ( fmin < fmax ) {
3334	ifgrp[j]++ ;
3335	}
3336	}
3337	}
3338	}
3339
3340	// Grouping continuous IF groups (without frequency gap)
3341	// freqgrp: list of IF group indexes in each frequency group
3342	// freqgrp[numgrp][nummember]
3343	// freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
3344	// [ifgrp10, ifgrp11, ifgrp12, ...],
3345	// ...
3346	// [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
3347	// grprange[2*numgrp]
3348	// grprange: [fmin0,fmax0,fmin1,fmax1,...]
3349	vector< vector<uInt> > freqgrp ;
3350	vector<double> grprange ;
3351	vector<uInt> grpedge ;
3352	for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
3353	if ( ifgrp[igrp] <= sizecr ) {
3354	grpedge.push_back( igrp ) ;
3355	}
3356	}
3357	grpedge.push_back( ifgrp.size() ) ;
3358	uInt itmp = 0 ;
3359	for ( uInt i = 0 ; i < grpedge.size() ; i++ ) {
3360	int n = grpedge[i] - itmp ;
3361	if ( n > 0 ) {
3362	vector<uInt> members( n ) ;
3363	for ( int j = 0 ; j < n ; j++ ) {
3364	members[j] = itmp+j ;
3365	}
3366	freqgrp.push_back( members ) ;
3367	grprange.push_back( sortedfreq[itmp] ) ;
3368	grprange.push_back( sortedfreq[grpedge[i]] ) ;
3369	}
3370	itmp += n + 1 ;
3371	}
3372
3373	// print frequency group
3374	oss.str("") ;
3375	oss << "Frequency Group summary: " << endl ;
3376	oss << " GROUP_ID: [FREQ_MIN, FREQ_MAX]" << endl ;
3377	for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
3378	oss << " GROUP " << setw( 2 ) << i << ": [" << grprange[2i] << "," << grprange[2i+1] << "]" ;
3379	oss << endl ;
3380	}
3381	oss << endl ;
3382	os << oss.str() << LogIO::POST ;
3383
3384	// groups: list of frequency group index whose frequency range overlaps
3385	// with that of each table and IF
3386	// groups[numtable][numIF]
3387	// groups: [[grpx, grpy,...],
3388	// [grpa, grpb,...],
3389	// ...
3390	// [grpk, grpm,...]]
3391	vector< vector<uInt> > groups( insize ) ;
3392	for ( uInt i = 0 ; i < insize ; i++ ) {
3393	groups[i].resize( freqid[i].size() ) ;
3394	}
3395	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3396	for ( uInt ifreq = 0 ; ifreq < freqid[itable].size() ; ifreq++ ) {
3397	double minf = iffreq[itable][2*ifreq] ;
3398	uInt groupid ;
3399	for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3400	vector<uInt> memberlist = freqgrp[igrp] ;
3401	if ( (minf >= grprange[2igrp]) && (minf <= grprange[2igrp+1]) ) {
3402	groupid = igrp ;
3403	break ;
3404	}
3405	}
3406	groups[itable][ifreq] = groupid ;
3407	}
3408	}
3409
3410
3411	// print membership
3412	oss.str("") ;
3413	for ( uInt i = 0 ; i < insize ; i++ ) {
3414	oss << "Table " << i << endl ;
3415	for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
3416	oss << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ;
3417	oss << setw( 2 ) << groups[i][j] ;
3418	oss << endl ;
3419	}
3420	}
3421	os << oss.str() << LogIO::POST ;
3422
3423	// reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
3424	//os << "All IF number is set to IF group index" << LogIO::POST ;
3425	// reset SCANNO only when avmode != "SCAN"
3426	if ( avmode != "SCAN" ) {
3427	os << "All scan number is set to 0" << LogIO::POST ;
3428	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3429	uInt nrow = newin[itable]->nrow() ;
3430	Vector<uInt> resetScan( nrow, 0 ) ;
3431	newin[itable]->scanCol_.putColumn( resetScan ) ;
3432	}
3433	}
3434
3435	// reset spectra and flagtra: align spectral resolution
3436	//os << "Align spectral resolution" << LogIO::POST ;
3437	// gmaxdnu: the coarsest frequency resolution in the frequency group
3438	// gminfreq: lower frequency edge of the frequency group
3439	// gnchan: number of channels for the frequency group
3440	vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
3441	vector<double> gminfreq( freqgrp.size() ) ;
3442	vector<double> gnchan( freqgrp.size() ) ;
3443	for ( uInt i = 0 ; i < insize ; i++ ) {
3444	vector<uInt> members = groups[i] ;
3445	for ( uInt j = 0 ; j < members.size() ; j++ ) {
3446	uInt groupid = members[j] ;
3447	Double rp,rv,ic ;
3448	newin[i]->frequencies().getEntry( rp, rv, ic, j ) ;
3449	if ( abs(ic) > abs(gmaxdnu[groupid]) )
3450	gmaxdnu[groupid] = ic ;
3451	}
3452	}
3453	for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3454	gminfreq[igrp] = grprange[2*igrp] ;
3455	double maxfreq = grprange[2*igrp+1] ;
3456	gnchan[igrp] = (int)(abs((maxfreq-gminfreq[igrp])/gmaxdnu[igrp])+0.9) ;
3457	}
3458
3459	// regrid spectral data and update frequency info
3460	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3461	Vector<uInt> oldFreqId = newin[itable]->mfreqidCol_.getColumn() ;
3462	Vector<uInt> newFreqId( oldFreqId.nelements() ) ;
3463
3464	// update MAIN
3465	for ( uInt irow = 0 ; irow < newin[itable]->nrow() ; irow++ ) {
3466	uInt groupid = groups[itable][oldFreqId[irow]] ;
3467	newFreqId[irow] = groupid ;
3468	newin[itable]->regridChannel( gnchan[groupid],
3469	gmaxdnu[groupid],
3470	gminfreq[groupid],
3471	irow ) ;
3472	}
3473	newin[itable]->mfreqidCol_.putColumn( newFreqId ) ;
3474	newin[itable]->ifCol_.putColumn( newFreqId ) ;
3475
3476	// update FREQUENCIES
3477	Table tab = newin[itable]->frequencies().table() ;
3478	ScalarColumn<uInt> fIdCol( tab, "ID" ) ;
3479	ScalarColumn<Double> fRefPixCol( tab, "REFPIX" ) ;
3480	ScalarColumn<Double> fRefValCol( tab, "REFVAL" ) ;
3481	ScalarColumn<Double> fIncrCol( tab, "INCREMENT" ) ;
3482	if ( freqgrp.size() > tab.nrow() ) {
3483	tab.addRow( freqgrp.size()-tab.nrow() ) ;
3484	}
3485	for ( uInt irow = 0 ; irow < freqgrp.size() ; irow++ ) {
3486	Double refval = gminfreq[irow] + 0.5 * abs(gmaxdnu[irow]) ;
3487	Double refpix = (gmaxdnu[irow] > 0.0) ? 0 : gnchan[irow]-1 ;
3488	Double increment = gmaxdnu[irow] ;
3489	fIdCol.put( irow, irow ) ;
3490	fRefPixCol.put( irow, refpix ) ;
3491	fRefValCol.put( irow, refval ) ;
3492	fIncrCol.put( irow, increment ) ;
3493	}
3494	}
3495
3496	// set back coordinfo
3497	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3498	vector<string> coordinfo = newin[itable]->getCoordInfo() ;
3499	coordinfo[0] = oldinfo[itable] ;
3500	newin[itable]->setCoordInfo( coordinfo ) ;
3501	}
3502
3503	// average
3504	out = average( newin, mask, weight, avmode ) ;
3505	}
3506	else {
3507	// simple average
3508	out = average( in, mask, weight, avmode ) ;
3509	}
3510
3511	return out;
3512	}
3513
3514	CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s,
3515	const String calmode,
3516	const String antname )
3517	{
3518	// frequency switch
3519	if ( calmode == "fs" ) {
3520	return cwcalfs( s, antname ) ;
3521	}
3522	else {
3523	vector<bool> masks = s->getMask( 0 ) ;
3524	vector<int> types ;
3525
3526	// save original table selection
3527	Table torg = s->table_ ;
3528
3529	// sky scan
3530	bool insitu = insitu_ ;
3531	insitu_ = false ;
3532	// share calibration scans before average with out
3533	CountedPtr<Scantable> out = getScantable( s, true ) ;
3534	insitu_ = insitu ;
3535	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
3536	out->attach() ;
3537	CountedPtr<Scantable> asky = averageWithinSession( out,
3538	masks,
3539	"TINT" ) ;
3540	// hot scan
3541	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
3542	out->attach() ;
3543	CountedPtr<Scantable> ahot = averageWithinSession( out,
3544	masks,
3545	"TINT" ) ;
3546	// cold scan
3547	CountedPtr<Scantable> acold ;
3548	// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
3549	// out->attach() ;
3550	// CountedPtr<Scantable> acold = averageWithinSession( out,
3551	// masks,
3552	// "TINT" ) ;
3553
3554	// off scan
3555	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3556	out->attach() ;
3557	CountedPtr<Scantable> aoff = averageWithinSession( out,
3558	masks,
3559	"TINT" ) ;
3560
3561	// on scan
3562	s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3563	s->attach() ;
3564	out->table_ = out->originalTable_ ;
3565	out->attach() ;
3566	out->table().addRow( s->nrow() ) ;
3567	copyRows( out->table(), s->table(), 0, 0, s->nrow(), False, True, False ) ;
3568
3569	// iterate throgh STIdxIter2
3570	ChopperWheelCalibrator cal(out, s, aoff, asky, ahot, acold);
3571	STIdxIter2::Iterate<ChopperWheelCalibrator>(cal, "BEAMNO,POLNO,IFNO");
3572
3573	s->table_ = torg ;
3574	s->attach() ;
3575
3576	// flux unit
3577	out->setFluxUnit( "K" ) ;
3578
3579	return out ;
3580	}
3581	}
3582
3583	CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s,
3584	const String calmode )
3585	{
3586	// frequency switch
3587	if ( calmode == "fs" ) {
3588	return almacalfs( s ) ;
3589	}
3590	else {
3591	// double t0, t1 ;
3592	// t0 = mathutil::gettimeofday_sec() ;
3593	vector<bool> masks = s->getMask( 0 ) ;
3594
3595	// save original table selection
3596	Table torg = s->table_ ;
3597
3598	// off scan
3599	// TODO 2010/01/08 TN
3600	// Grouping by time should be needed before averaging.
3601	// Each group must have own unique SCANNO (should be renumbered).
3602	// See PIPELINE/SDCalibration.py
3603	bool insitu = insitu_ ;
3604	insitu_ = false ;
3605	// share off scan before average with out
3606	CountedPtr<Scantable> out = getScantable( s, true ) ;
3607	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3608	out->attach() ;
3609	insitu_ = insitu ;
3610	CountedPtr<Scantable> aoff = averageWithinSession( out,
3611	masks,
3612	"TINT" ) ;
3613
3614	// on scan
3615	// t0 = mathutil::gettimeofday_sec() ;
3616	s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3617	s->attach() ;
3618	out->table_ = out->originalTable_ ;
3619	out->attach() ;
3620	out->table().addRow( s->nrow() ) ;
3621	copyRows( out->table(), s->table(), 0, 0, s->nrow(), False ) ;
3622	// t1 = mathutil::gettimeofday_sec() ;
3623	// cout << "elapsed time for preparing output table: " << t1-t0 << " sec" << endl ;
3624
3625	// process each on scan
3626	// t0 = mathutil::gettimeofday_sec() ;
3627
3628	// iterate throgh STIdxIter2
3629	AlmaCalibrator cal(out, s, aoff);
3630	STIdxIter2::Iterate<AlmaCalibrator>(cal, "BEAMNO,POLNO,IFNO");
3631
3632	// finalize
3633	s->table_ = torg ;
3634	s->attach() ;
3635
3636	// t1 = mathutil::gettimeofday_sec() ;
3637	// cout << "elapsed time for calibration: " << t1-t0 << " sec" << endl ;
3638
3639	// flux unit
3640	out->setFluxUnit( "K" ) ;
3641
3642	return out ;
3643	}
3644	}
3645
3646	CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
3647	const String antname )
3648	{
3649	vector<int> types ;
3650
3651	// APEX calibration mode
3652	int apexcalmode = 1 ;
3653
3654	if ( antname.find( "APEX" ) != string::npos ) {
3655	// check if off scan exists or not
3656	STSelector sel = STSelector() ;
3657	//sel.setName( offstr1 ) ;
3658	types.push_back( SrcType::FLOOFF ) ;
3659	sel.setTypes( types ) ;
3660	try {
3661	s->setSelection( sel ) ;
3662	}
3663	catch ( AipsError &e ) {
3664	apexcalmode = 0 ;
3665	}
3666	sel.reset() ;
3667	}
3668	s->unsetSelection() ;
3669	types.clear() ;
3670
3671	vector<bool> masks = s->getMask( 0 ) ;
3672	CountedPtr<Scantable> ssig, sref ;
3673	//CountedPtr<Scantable> out ;
3674	bool insitu = insitu_ ;
3675	insitu_ = False ;
3676	CountedPtr<Scantable> out = getScantable( s, true ) ;
3677	insitu_ = insitu ;
3678
3679	if ( antname.find( "APEX" ) != string::npos ) {
3680	// APEX calibration
3681	// sky scan
3682	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOSKY ) ;
3683	out->attach() ;
3684	CountedPtr<Scantable> askylo = averageWithinSession( out,
3685	masks,
3686	"TINT" ) ;
3687	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHISKY ) ;
3688	out->attach() ;
3689	CountedPtr<Scantable> askyhi = averageWithinSession( out,
3690	masks,
3691	"TINT" ) ;
3692
3693	// hot scan
3694	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOHOT ) ;
3695	out->attach() ;
3696	CountedPtr<Scantable> ahotlo = averageWithinSession( out,
3697	masks,
3698	"TINT" ) ;
3699	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIHOT ) ;
3700	out->attach() ;
3701	CountedPtr<Scantable> ahothi = averageWithinSession( out,
3702	masks,
3703	"TINT" ) ;
3704
3705	// cold scan
3706	CountedPtr<Scantable> acoldlo, acoldhi ;
3707	// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOCOLD ) ;
3708	// out->attach() ;
3709	// CountedPtr<Scantable> acoldlo = averageWithinSession( out,
3710	// masks,
3711	// "TINT" ) ;
3712	// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHICOLD ) ;
3713	// out->attach() ;
3714	// CountedPtr<Scantable> acoldhi = averageWithinSession( out,
3715	// masks,
3716	// "TINT" ) ;
3717
3718	// ref scan
3719	insitu_ = false ;
3720	sref = getScantable( s, true ) ;
3721	CountedPtr<Scantable> rref = getScantable( s, true ) ;
3722	insitu_ = insitu ;
3723	rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSLO ) ;
3724	rref->attach() ;
3725	copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
3726
3727	// sig scan
3728	insitu_ = false ;
3729	ssig = getScantable( s, true ) ;
3730	CountedPtr<Scantable> rsig = getScantable( s, true ) ;
3731	insitu_ = insitu ;
3732	rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSHI ) ;
3733	rsig->attach() ;
3734	copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
3735
3736	if ( apexcalmode == 0 ) {
3737	// using STIdxIter2
3738	vector<string> cols( 3 ) ;
3739	cols[0] = "BEAMNO" ;
3740	cols[1] = "POLNO" ;
3741	cols[2] = "IFNO" ;
3742	STIdxIter2 iter(ssig, cols) ;
3743	STSelector sel ;
3744	vector< CountedPtr<Scantable> > on( 2 ) ;
3745	on[0] = rsig ;
3746	on[1] = rref ;
3747	vector< CountedPtr<Scantable> > sky( 2 ) ;
3748	sky[0] = askylo ;
3749	sky[1] = askyhi ;
3750	vector< CountedPtr<Scantable> > hot( 2 ) ;
3751	hot[0] = ahotlo ;
3752	hot[1] = ahothi ;
3753	vector< CountedPtr<Scantable> > cold( 2 ) ;
3754	while ( !iter.pastEnd() ) {
3755	Record ids = iter.currentValue() ;
3756	stringstream ss ;
3757	ss << "SELECT FROM $1 WHERE "
3758	<< "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3759	<< "POLNO==" << ids.asuInt("POLNO") << "&&"
3760	<< "IFNO==" << ids.asuInt("IFNO") ;
3761	//cout << "TaQL string: " << ss.str() << endl ;
3762	sel.setTaQL( ss.str() ) ;
3763	sky[0]->setSelection( sel ) ;
3764	sky[1]->setSelection( sel ) ;
3765	hot[0]->setSelection( sel ) ;
3766	hot[1]->setSelection( sel ) ;
3767	Vector<uInt> rows = iter.getRows( SHARE ) ;
3768	calibrateAPEXFS( ssig, sref, on, sky, hot, cold, rows ) ;
3769	sky[0]->unsetSelection() ;
3770	sky[1]->unsetSelection() ;
3771	hot[0]->unsetSelection() ;
3772	hot[1]->unsetSelection() ;
3773	sel.reset() ;
3774	iter.next() ;
3775	}
3776
3777	}
3778	else if ( apexcalmode == 1 ) {
3779	// APEX fs data with off scan
3780	// off scan
3781	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOOFF ) ;
3782	out->attach() ;
3783	CountedPtr<Scantable> aofflo = averageWithinSession( out,
3784	masks,
3785	"TINT" ) ;
3786	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIOFF ) ;
3787	out->attach() ;
3788	CountedPtr<Scantable> aoffhi = averageWithinSession( out,
3789	masks,
3790	"TINT" ) ;
3791
3792	// process each sig and ref scan
3793	// iterate throgh STIdxIter2
3794	ChopperWheelCalibrator cal_sig(ssig, rsig, aofflo, askylo, ahotlo, acoldlo);
3795	STIdxIter2::Iterate<ChopperWheelCalibrator>(cal_sig, "BEAMNO,POLNO,IFNO");
3796	ChopperWheelCalibrator cal_ref(sref, rref, aoffhi, askyhi, ahothi, acoldhi);
3797	STIdxIter2::Iterate<ChopperWheelCalibrator>(cal_ref, "BEAMNO,POLNO,IFNO");
3798	}
3799	}
3800	else {
3801	// non-APEX fs data
3802	// sky scan
3803	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
3804	out->attach() ;
3805	CountedPtr<Scantable> asky = averageWithinSession( out,
3806	masks,
3807	"TINT" ) ;
3808	STSelector sel = STSelector() ;
3809
3810	// hot scan
3811	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
3812	out->attach() ;
3813	CountedPtr<Scantable> ahot = averageWithinSession( out,
3814	masks,
3815	"TINT" ) ;
3816
3817	// cold scan
3818	CountedPtr<Scantable> acold ;
3819	// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
3820	// out->attach() ;
3821	// CountedPtr<Scantable> acold = averageWithinSession( out,
3822	// masks,
3823	// "TINT" ) ;
3824
3825	// ref scan
3826	bool insitu = insitu_ ;
3827	insitu_ = false ;
3828	sref = getScantable( s, true ) ;
3829	CountedPtr<Scantable> rref = getScantable( s, true ) ;
3830	insitu_ = insitu ;
3831	rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3832	rref->attach() ;
3833	copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
3834
3835	// sig scan
3836	insitu_ = false ;
3837	ssig = getScantable( s, true ) ;
3838	CountedPtr<Scantable> rsig = getScantable( s, true ) ;
3839	insitu_ = insitu ;
3840	rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3841	rsig->attach() ;
3842	copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
3843
3844	// process each sig and ref scan
3845	vector<string> cols( 3 ) ;
3846	cols[0] = "BEAMNO" ;
3847	cols[1] = "POLNO" ;
3848	cols[2] = "IFNO" ;
3849	STIdxIter2 iter(ssig, cols);
3850	while ( !iter.pastEnd() ) {
3851	Record ids = iter.currentValue() ;
3852	stringstream ss ;
3853	ss << "SELECT FROM $1 WHERE "
3854	<< "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3855	<< "POLNO==" << ids.asuInt("POLNO") << "&&"
3856	<< "IFNO==" << ids.asuInt("IFNO") ;
3857	//cout << "TaQL string: " << ss.str() << endl ;
3858	sel.setTaQL( ss.str() ) ;
3859	ahot->setSelection( sel ) ;
3860	asky->setSelection( sel ) ;
3861	Vector<uInt> rows = iter.getRows( SHARE ) ;
3862	// out should be an exact copy of s except that SPECTRA column is empty
3863	calibrateFS( ssig, sref, rsig, rref, asky, ahot, acold, rows ) ;
3864	ahot->unsetSelection() ;
3865	asky->unsetSelection() ;
3866	sel.reset() ;
3867	iter.next() ;
3868	}
3869	}
3870
3871	// do folding if necessary
3872	Table sigtab = ssig->table() ;
3873	Table reftab = sref->table() ;
3874	ScalarColumn<uInt> reffidCol ;
3875	Int nchan = (Int)ssig->nchan() ;
3876	reffidCol.attach( reftab, "FREQ_ID" ) ;
3877	Vector<uInt> sfids = ssig->mfreqidCol_.getColumn() ;
3878	Vector<uInt> rfids = sref->mfreqidCol_.getColumn() ;
3879	vector<uInt> sfids_unique ;
3880	vector<uInt> rfids_unique ;
3881	vector<uInt> sifno_unique ;
3882	vector<uInt> rifno_unique ;
3883	for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) {
3884	if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) {
3885	sfids_unique.push_back( sfids[i] ) ;
3886	sifno_unique.push_back( ssig->getIF( i ) ) ;
3887	}
3888	if ( count( rfids_unique.begin(), rfids_unique.end(), rfids[i] ) == 0 ) {
3889	rfids_unique.push_back( rfids[i] ) ;
3890	rifno_unique.push_back( sref->getIF( i ) ) ;
3891	}
3892	}
3893	double refpix_sig, refval_sig, increment_sig ;
3894	double refpix_ref, refval_ref, increment_ref ;
3895	vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ;
3896	for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) {
3897	ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ;
3898	sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ;
3899	if ( refpix_sig == refpix_ref ) {
3900	double foffset = refval_ref - refval_sig ;
3901	int choffset = static_cast<int>(foffset/increment_sig) ;
3902	double doffset = foffset / increment_sig ;
3903	if ( abs(choffset) >= nchan ) {
3904	LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ;
3905	os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ;
3906	os << "Just return signal data" << LogIO::POST ;
3907	//std::vector< CountedPtr<Scantable> > tabs ;
3908	//tabs.push_back( ssig ) ;
3909	//tabs.push_back( sref ) ;
3910	//out = merge( tabs ) ;
3911	tmp[i] = ssig ;
3912	}
3913	else {
3914	STSelector sel = STSelector() ;
3915	vector<int> v( 1, sifno_unique[i] ) ;
3916	sel.setIFs( v ) ;
3917	ssig->setSelection( sel ) ;
3918	sel.reset() ;
3919	v[0] = rifno_unique[i] ;
3920	sel.setIFs( v ) ;
3921	sref->setSelection( sel ) ;
3922	sel.reset() ;
3923	if ( antname.find( "APEX" ) != string::npos ) {
3924	tmp[i] = dofold( ssig, sref, 0.5doffset, -0.5doffset ) ;
3925	//tmp[i] = dofold( ssig, sref, doffset ) ;
3926	}
3927	else {
3928	tmp[i] = dofold( ssig, sref, doffset ) ;
3929	}
3930	ssig->unsetSelection() ;
3931	sref->unsetSelection() ;
3932	}
3933	}
3934	}
3935
3936	if ( tmp.size() > 1 ) {
3937	out = merge( tmp ) ;
3938	}
3939	else {
3940	out = tmp[0] ;
3941	}
3942
3943	// flux unit
3944	out->setFluxUnit( "K" ) ;
3945
3946	return out ;
3947	}
3948
3949	CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s )
3950	{
3951	(void) s; //currently unused
3952	CountedPtr<Scantable> out ;
3953
3954	return out ;
3955	}
3956
3957	void STMath::calibrateAPEXFS( CountedPtr<Scantable> &sig,
3958	CountedPtr<Scantable> &ref,
3959	const vector< CountedPtr<Scantable> >& on,
3960	const vector< CountedPtr<Scantable> >& sky,
3961	const vector< CountedPtr<Scantable> >& hot,
3962	const vector< CountedPtr<Scantable> >& cold,
3963	const Vector<uInt> &rows )
3964	{
3965	// if rows is empty, just return
3966	if ( rows.nelements() == 0 )
3967	return ;
3968	ROScalarColumn<Double> timeCol( sky[0]->table(), "TIME" ) ;
3969	Vector<Double> timeSkyS = timeCol.getColumn() ;
3970	timeCol.attach( sky[1]->table(), "TIME" ) ;
3971	Vector<Double> timeSkyR = timeCol.getColumn() ;
3972	timeCol.attach( hot[0]->table(), "TIME" ) ;
3973	Vector<Double> timeHotS = timeCol.getColumn() ;
3974	timeCol.attach( hot[1]->table(), "TIME" ) ;
3975	Vector<Double> timeHotR = timeCol.getColumn() ;
3976	timeCol.attach( sig->table(), "TIME" ) ;
3977	ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
3978	ROArrayColumn<Float> arrayFloatCol( sky[0]->table(), "SPECTRA" ) ;
3979	SpectralData skyspectraS(arrayFloatCol.getColumn());
3980	arrayFloatCol.attach( sky[1]->table(), "SPECTRA" ) ;
3981	SpectralData skyspectraR(arrayFloatCol.getColumn());
3982	arrayFloatCol.attach( hot[0]->table(), "SPECTRA" ) ;
3983	SpectralData hotspectraS(arrayFloatCol.getColumn());
3984	arrayFloatCol.attach( hot[1]->table(), "SPECTRA" ) ;
3985	SpectralData hotspectraR(arrayFloatCol.getColumn());
3986	TcalData tcaldataS(sky[0]);
3987	TsysData tsysdataS(sky[0]);
3988	TcalData tcaldataR(sky[1]);
3989	TsysData tsysdataR(sky[1]);
3990	unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
3991	Vector<Float> spec( spsize ) ;
3992	// I know that the data is contiguous
3993	const uInt *p = rows.data() ;
3994	vector<int> ids( 2 ) ;
3995	Block<uInt> flagchan( spsize ) ;
3996	uInt nflag = 0 ;
3997	for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
3998	double reftime = timeCol.asdouble(*p) ;
3999	ids = getRowIdFromTime( reftime, timeSkyS ) ;
4000	Vector<Float> spskyS = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSkyS, ids, skyspectraS, "linear");
4001	Vector<Float> tcalS = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSkyS, ids, tcaldataS, "linear");
4002	Vector<Float> tsysS = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSkyS, ids, tsysdataS, "linear");
4003	ids = getRowIdFromTime( reftime, timeHotS ) ;
4004	Vector<Float> sphotS = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHotS, ids, hotspectraS, "linear");
4005	reftime = timeCol2.asdouble(*p) ;
4006	ids = getRowIdFromTime( reftime, timeSkyR ) ;
4007	Vector<Float> spskyR = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSkyR, ids, skyspectraR, "linear");
4008	Vector<Float> tcalR = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSkyR, ids, tcaldataR, "linear");
4009	Vector<Float> tsysR = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSkyR, ids, tsysdataR, "linear");
4010	ids = getRowIdFromTime( reftime, timeHotR ) ;
4011	Vector<Float> sphotR = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHotR, ids, hotspectraR, "linear");
4012	Vector<Float> spsig = on[0]->specCol_( *p ) ;
4013	Vector<Float> spref = on[1]->specCol_( *p ) ;
4014	for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4015	if ( (sphotS[j]-spskyS[j]) == 0.0 \|\| (sphotR[j]-spskyR[j]) == 0.0 ) {
4016	spec[j] = 0.0 ;
4017	flagchan[nflag++] = j ;
4018	}
4019	else {
4020	spec[j] = tcalS[j] * spsig[j] / ( sphotS[j] - spskyS[j] )
4021	- tcalR[j] * spref[j] / ( sphotR[j] - spskyR[j] ) ;
4022	}
4023	}
4024	sig->specCol_.put( *p, spec ) ;
4025	sig->tsysCol_.put( *p, tsysS ) ;
4026	spec *= (Float)-1.0 ;
4027	ref->specCol_.put( *p, spec ) ;
4028	ref->tsysCol_.put( *p, tsysR ) ;
4029	if ( nflag > 0 ) {
4030	Vector<uChar> flsig = sig->flagsCol_( *p ) ;
4031	Vector<uChar> flref = ref->flagsCol_( *p ) ;
4032	for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4033	flsig[flagchan[j]] = (uChar)True ;
4034	flref[flagchan[j]] = (uChar)True ;
4035	}
4036	sig->flagsCol_.put( *p, flsig ) ;
4037	ref->flagsCol_.put( *p, flref ) ;
4038	}
4039	nflag = 0 ;
4040	p++ ;
4041	}
4042	}
4043
4044	void STMath::calibrateFS( CountedPtr<Scantable> &sig,
4045	CountedPtr<Scantable> &ref,
4046	const CountedPtr<Scantable>& rsig,
4047	const CountedPtr<Scantable>& rref,
4048	const CountedPtr<Scantable>& sky,
4049	const CountedPtr<Scantable>& hot,
4050	const CountedPtr<Scantable>& cold,
4051	const Vector<uInt> &rows )
4052	{
4053	// if rows is empty, just return
4054	if ( rows.nelements() == 0 )
4055	return ;
4056	ROScalarColumn<Double> timeCol( sky->table(), "TIME" ) ;
4057	Vector<Double> timeSky = timeCol.getColumn() ;
4058	timeCol.attach( hot->table(), "TIME" ) ;
4059	Vector<Double> timeHot = timeCol.getColumn() ;
4060	timeCol.attach( sig->table(), "TIME" ) ;
4061	ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
4062	ROArrayColumn<Float> arrayFloatCol( sky->table(), "SPECTRA" ) ;
4063	SpectralData skyspectra(arrayFloatCol.getColumn());
4064	arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
4065	SpectralData hotspectra(arrayFloatCol.getColumn());
4066	TcalData tcaldata(sky);
4067	TsysData tsysdata(sky);
4068	unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
4069	Vector<Float> spec( spsize ) ;
4070	// I know that the data is contiguous
4071	const uInt *p = rows.data() ;
4072	vector<int> ids( 2 ) ;
4073	Block<uInt> flagchan( spsize ) ;
4074	uInt nflag = 0 ;
4075	for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4076	double reftime = timeCol.asdouble(*p) ;
4077	ids = getRowIdFromTime( reftime, timeSky ) ;
4078	Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
4079	Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
4080	Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
4081	ids = getRowIdFromTime( reftime, timeHot ) ;
4082	Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
4083	Vector<Float> spsig = rsig->specCol_( *p ) ;
4084	Vector<Float> spref = rref->specCol_( *p ) ;
4085	// using gain array
4086	for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4087	if ( spref[j] == 0.0 \|\| (sphot[j]-spsky[j]) == 0.0 ) {
4088	spec[j] = 0.0 ;
4089	flagchan[nflag++] = j ;
4090	}
4091	else {
4092	spec[j] = ( ( spsig[j] - spref[j] ) / spref[j] )
4093	* ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4094	}
4095	}
4096	sig->specCol_.put( *p, spec ) ;
4097	sig->tsysCol_.put( *p, tsys ) ;
4098	if ( nflag > 0 ) {
4099	Vector<uChar> fl = sig->flagsCol_( *p ) ;
4100	for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4101	fl[flagchan[j]] = (uChar)True ;
4102	}
4103	sig->flagsCol_.put( *p, fl ) ;
4104	}
4105	nflag = 0 ;
4106
4107	reftime = timeCol2.asdouble(*p) ;
4108	ids = getRowIdFromTime( reftime, timeSky ) ;
4109	spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
4110	tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
4111	tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
4112	ids = getRowIdFromTime( reftime, timeHot ) ;
4113	sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
4114	// using gain array
4115	for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4116	if ( spsig[j] == 0.0 \|\| (sphot[j]-spsky[j]) == 0.0 ) {
4117	spec[j] = 0.0 ;
4118	flagchan[nflag++] = j ;
4119	}
4120	else {
4121	spec[j] = ( ( spref[j] - spsig[j] ) / spsig[j] )
4122	* ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4123	}
4124	}
4125	ref->specCol_.put( *p, spec ) ;
4126	ref->tsysCol_.put( *p, tsys ) ;
4127	if ( nflag > 0 ) {
4128	Vector<uChar> fl = ref->flagsCol_( *p ) ;
4129	for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4130	fl[flagchan[j]] = (uChar)True ;
4131	}
4132	ref->flagsCol_.put( *p, fl ) ;
4133	}
4134	nflag = 0 ;
4135	p++ ;
4136	}
4137	}
4138
4139	void STMath::copyRows( Table &out,
4140	const Table &in,
4141	uInt startout,
4142	uInt startin,
4143	uInt nrow,
4144	Bool copySpectra,
4145	Bool copyFlagtra,
4146	Bool copyTsys )
4147	{
4148	uInt nexclude = 0 ;
4149	Block<String> excludeColsBlock( 3 ) ;
4150	if ( !copySpectra ) {
4151	excludeColsBlock[nexclude] = "SPECTRA" ;
4152	nexclude++ ;
4153	}
4154	if ( !copyFlagtra ) {
4155	excludeColsBlock[nexclude] = "FLAGTRA" ;
4156	nexclude++ ;
4157	}
4158	if ( !copyTsys ) {
4159	excludeColsBlock[nexclude] = "TSYS" ;
4160	nexclude++ ;
4161	}
4162	// if ( nexclude < 3 ) {
4163	// excludeCols.resize( nexclude, True ) ;
4164	// }
4165	Vector<String> excludeCols( IPosition(1,nexclude),
4166	excludeColsBlock.storage(),
4167	SHARE ) ;
4168	// cout << "excludeCols=" << excludeCols << endl ;
4169	TableRow rowout( out, excludeCols, True ) ;
4170	ROTableRow rowin( in, excludeCols, True ) ;
4171	uInt rin = startin ;
4172	uInt rout = startout ;
4173	for ( uInt i = 0 ; i < nrow ; i++ ) {
4174	rowin.get( rin ) ;
4175	rowout.putMatchingFields( rout, rowin.record() ) ;
4176	rin++ ;
4177	rout++ ;
4178	}
4179	}
4180
4181	CountedPtr<Scantable> STMath::averageWithinSession( CountedPtr<Scantable> &s,
4182	vector<bool> &mask,
4183	string weight )
4184	{
4185	// prepare output table
4186	bool insitu = insitu_ ;
4187	insitu_ = false ;
4188	CountedPtr<Scantable> a = getScantable( s, true ) ;
4189	insitu_ = insitu ;
4190	Table &atab = a->table() ;
4191	ScalarColumn<Double> timeColOut( atab, "TIME" ) ;
4192
4193	if ( s->nrow() == 0 )
4194	return a ;
4195
4196	// setup RowAccumulator
4197	WeightType wtype = stringToWeight( weight ) ;
4198	RowAccumulator acc( wtype ) ;
4199	Vector<Bool> cmask( mask ) ;
4200	acc.setUserMask( cmask ) ;
4201
4202	vector<string> cols( 3 ) ;
4203	cols[0] = "IFNO" ;
4204	cols[1] = "POLNO" ;
4205	cols[2] = "BEAMNO" ;
4206	STIdxIter2 iter( s, cols ) ;
4207
4208	Table ttab = s->table() ;
4209	ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ;
4210	Vector<Double> timeVec = timeCol->getColumn() ;
4211	delete timeCol ;
4212	Vector<Double> interval = s->integrCol_.getColumn() ;
4213	uInt nrow = timeVec.nelements() ;
4214	uInt outrow = 0 ;
4215
4216	while( !iter.pastEnd() ) {
4217
4218	Vector<uInt> rows = iter.getRows( SHARE ) ;
4219	uInt len = rows.nelements() ;
4220
4221	if ( len == 0 ) {
4222	iter.next() ;
4223	continue ;
4224	}
4225
4226	uInt nchan = s->nchan(s->getIF(rows[0])) ;
4227	Vector<uChar> flag( nchan ) ;
4228	Vector<Bool> bflag( nchan ) ;
4229	Vector<Float> spec( nchan ) ;
4230	Vector<Float> tsys( nchan ) ;
4231
4232	Vector<Double> timeSep( len-1 ) ;
4233	for ( uInt i = 0 ; i < len-1 ; i++ ) {
4234	timeSep[i] = timeVec[rows[i+1]] - timeVec[rows[i]] ;
4235	}
4236
4237	uInt irow ;
4238	uInt jrow ;
4239	for ( uInt i = 0 ; i < len-1 ; i++ ) {
4240	irow = rows[i] ;
4241	jrow = rows[i+1] ;
4242	// accumulate data
4243	s->flagsCol_.get( irow, flag ) ;
4244	convertArray( bflag, flag ) ;
4245	s->specCol_.get( irow, spec ) ;
4246	tsys.assign( s->tsysCol_( irow ) ) ;
4247	if ( !allEQ(bflag,True) )
4248	acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
4249	double gap = 2.0 * 86400.0 * timeSep[i] / ( interval[jrow] + interval[irow] ) ;
4250	//cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
4251	if ( gap > 1.1 ) {
4252	//cout << "detected gap between " << i << " and " << i+1 << endl ;
4253	// put data to output table
4254	// reset RowAccumulator
4255	if ( acc.state() ) {
4256	atab.addRow() ;
4257	copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
4258	acc.replaceNaN() ;
4259	const Vector<Bool> &msk = acc.getMask() ;
4260	convertArray( flag, !msk ) ;
4261	for (uInt k = 0; k < nchan; ++k) {
4262	uChar userFlag = 1 << 7;
4263	if (msk[k]==True) userFlag = 0 << 7;
4264	flag(k) = userFlag;
4265	}
4266	a->flagsCol_.put( outrow, flag ) ;
4267	a->specCol_.put( outrow, acc.getSpectrum() ) ;
4268	a->tsysCol_.put( outrow, acc.getTsys() ) ;
4269	a->integrCol_.put( outrow, acc.getInterval() ) ;
4270	timeColOut.put( outrow, acc.getTime() ) ;
4271	a->cycleCol_.put( outrow, 0 ) ;
4272	}
4273	acc.reset() ;
4274	outrow++ ;
4275	}
4276	}
4277
4278	// accumulate and add last data
4279	irow = rows[len-1] ;
4280	s->flagsCol_.get( irow, flag ) ;
4281	convertArray( bflag, flag ) ;
4282	s->specCol_.get( irow, spec ) ;
4283	tsys.assign( s->tsysCol_( irow ) ) ;
4284	if (!allEQ(bflag,True) )
4285	acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
4286	if ( acc.state() ) {
4287	atab.addRow() ;
4288	copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
4289	acc.replaceNaN() ;
4290	const Vector<Bool> &msk = acc.getMask() ;
4291	convertArray( flag, !msk ) ;
4292	for (uInt k = 0; k < nchan; ++k) {
4293	uChar userFlag = 1 << 7;
4294	if (msk[k]==True) userFlag = 0 << 7;
4295	flag(k) = userFlag;
4296	}
4297	a->flagsCol_.put( outrow, flag ) ;
4298	a->specCol_.put( outrow, acc.getSpectrum() ) ;
4299	a->tsysCol_.put( outrow, acc.getTsys() ) ;
4300	a->integrCol_.put( outrow, acc.getInterval() ) ;
4301	timeColOut.put( outrow, acc.getTime() ) ;
4302	a->cycleCol_.put( outrow, 0 ) ;
4303	}
4304	acc.reset() ;
4305	outrow++ ;
4306
4307	iter.next() ;
4308	}
4309
4310	return a ;
4311	}

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

Context Navigation

source: trunk/src/STMath.cpp @ 2959

Download in other formats: