Context Navigation

source: trunk/src/STMath.cpp @ 3089

Visit:

Last change on this file since 3089 was 3089, checked in by Kana Sugimoto, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description: fixes to get rid of clang build warnings.

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

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

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

Download in other formats: