Context Navigation

source: trunk/src/STMath.cpp @ 2974

Visit:

Last change on this file since 2974 was 2974, checked in by WataruKawasaki, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6594

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs: test_sdmath

Put in Release Notes: No

Module(s): sd

Description: modified STMath::binaryOperate() so that:

(1) If one of or both spectra are row-flagged,

operation will be skipped so the 'left'
spectrum values are not modified, while
the FLAGROW of the output spectrum will
be set flagged.

(2) Otherwise, operation is executed for all

channels even if there are masked ones.
If a channel is flagged in one of or both
spectra, operation is done for that channel
but the channel is flagged.

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

File size: 148.1 KB

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

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

Download in other formats: