Context Navigation

STMath.cpp @ 2928

Visit:

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: 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..

Merged IterationHelper? into STIdxIter2.

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

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

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

Context Navigation

source: trunk/src/STMath.cpp @ 2928

Download in other formats: