Context Navigation

STMath.cpp @ 3106

Visit:

Last change on this file since 3106 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No?

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Check-in asap modifications from Jim regarding casacore namespace conversion.

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

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

Context Navigation

source: trunk/src/STMath.cpp @ 3106

Download in other formats: