Context Navigation

source: trunk/src/STMath.cpp @ 3010

Visit:

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

New Development: No

JIRA Issue: Yes CAS-6599

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sd

Description: modified so that sd.scantable.stats() returns None for flagged rows and rows with all channels flagged.

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

File size: 148.5 KB

Line
1	//
2	// C++ Implementation: STMath
3	//
4	// Description:
5	//
6	//
7	// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8	//
9	// Copyright: See COPYING file that comes with this distribution
10	//
11	//
12
13	//#include <math.h>
14	//#include <bits/nan.h>
15	#include <sstream>
16	#include <iostream>
17
18	#include <casa/iomanip.h>
19	#include <casa/Arrays/MaskArrLogi.h>
20	#include <casa/Arrays/MaskArrMath.h>
21	#include <casa/Arrays/ArrayLogical.h>
22	#include <casa/Arrays/ArrayMath.h>
23	#include <casa/Arrays/Slice.h>
24	#include <casa/Arrays/Slicer.h>
25	#include <casa/BasicSL/String.h>
26	#include <casa/Containers/Block.h>
27	#include <casa/Containers/RecordField.h>
28	#include <casa/Exceptions/Error.h>
29	#include <casa/Logging/LogIO.h>
30	#include <casa/Quanta/Quantum.h>
31
32	#include <coordinates/Coordinates/CoordinateSystem.h>
33	#include <coordinates/Coordinates/CoordinateUtil.h>
34	#include <coordinates/Coordinates/FrequencyAligner.h>
35	#include <coordinates/Coordinates/SpectralCoordinate.h>
36
37	#include <lattices/Lattices/LatticeUtilities.h>
38
39	#include <scimath/Functionals/Polynomial.h>
40	#include <scimath/Mathematics/Convolver.h>
41	#include <scimath/Mathematics/VectorKernel.h>
42
43	#include <tables/Tables/ExprNode.h>
44	#include <tables/Tables/ReadAsciiTable.h>
45	#include <tables/Tables/TableCopy.h>
46	#include <tables/Tables/TableIter.h>
47	#include <tables/Tables/TableParse.h>
48	#include <tables/Tables/TableRecord.h>
49	#include <tables/Tables/TableRow.h>
50	#include <tables/Tables/TableVector.h>
51	#include <tables/Tables/TabVecMath.h>
52
53	#include <atnf/PKSIO/SrcType.h>
54
55	#include "RowAccumulator.h"
56	#include "STAttr.h"
57	#include "STMath.h"
58	#include "STSelector.h"
59	#include "Accelerator.h"
60	#include "STIdxIter.h"
61
62	#include "CalibrationHelper.h"
63
64	using namespace casa;
65	using namespace asap;
66
67	// 2012/02/17 TN
68	// Since STGrid is implemented, average doesn't consider direction
69	// when accumulating
70	// tolerance for direction comparison (rad)
71	// #define TOL_OTF 1.0e-15
72	// #define TOL_POINT 2.9088821e-4 // 1 arcmin
73
74	STMath::STMath(bool insitu) :
75	insitu_(insitu)
76	{
77	}
78
79
80	STMath::~STMath()
81	{
82	}
83
84	CountedPtr<Scantable>
85	STMath::average( const std::vector<CountedPtr<Scantable> >& in,
86	const std::vector<bool>& mask,
87	const std::string& weight,
88	const std::string& avmode)
89	{
90	// double t0, t1 ;
91	// t0 = mathutil::gettimeofday_sec() ;
92
93	LogIO os( LogOrigin( "STMath", "average()", WHERE ) ) ;
94	if ( avmode == "SCAN" && in.size() != 1 )
95	throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
96	"Use merge first."));
97	WeightType wtype = stringToWeight(weight);
98
99	// 2012/02/17 TN
100	// Since STGrid is implemented, average doesn't consider direction
101	// when accumulating
102	// check if OTF observation
103	// String obstype = in[0]->getHeader().obstype ;
104	// Double tol = 0.0 ;
105	// if ( (obstype.find( "OTF" ) != String::npos) \|\| (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
106	// tol = TOL_OTF ;
107	// }
108	// else {
109	// tol = TOL_POINT ;
110	// }
111
112	// output
113	// clone as this is non insitu
114	bool insitu = insitu_;
115	setInsitu(false);
116	CountedPtr< Scantable > out = getScantable(in[0], true);
117	setInsitu(insitu);
118	std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
119	++stit;
120	while ( stit != in.end() ) {
121	out->appendToHistoryTable((*stit)->history());
122	++stit;
123	}
124
125	Table& tout = out->table();
126
127	/// @todo check if all scantables are conformant
128
129	ArrayColumn<Float> specColOut(tout,"SPECTRA");
130	ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
131	ArrayColumn<Float> tsysColOut(tout,"TSYS");
132	ScalarColumn<Double> mjdColOut(tout,"TIME");
133	ScalarColumn<Double> intColOut(tout,"INTERVAL");
134	ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
135	ScalarColumn<uInt> scanColOut(tout,"SCANNO");
136	ScalarColumn<uInt> flagRowColOut(tout,"FLAGROW");
137
138	// set up the output table rows. These are based on the structure of the
139	// FIRST scantable in the vector
140	const Table& baset = in[0]->table();
141
142	RowAccumulator acc(wtype);
143	Vector<Bool> cmask(mask);
144	acc.setUserMask(cmask);
145	// ROTableRow row(tout);
146	ROArrayColumn<Float> specCol, tsysCol;
147	ROArrayColumn<uChar> flagCol;
148	ROScalarColumn<Double> mjdCol, intCol;
149	ROScalarColumn<Int> scanIDCol;
150	ROScalarColumn<uInt> flagRowCol;
151
152	//Vector<uInt> rowstodelete;
153	Block<uInt> rowstodelB( in[0]->nrow() ) ;
154	uInt nrowdel = 0 ;
155
156	// Block<String> cols(3);
157	vector<string> cols(3) ;
158	cols[0] = String("BEAMNO");
159	cols[1] = String("IFNO");
160	cols[2] = String("POLNO");
161	if ( avmode == "SOURCE" ) {
162	cols.resize(4);
163	cols[3] = String("SRCNAME");
164	}
165	if ( avmode == "SCAN" && in.size() == 1) {
166	//cols.resize(4);
167	//cols[3] = String("SCANNO");
168	cols.resize(5);
169	cols[3] = String("SRCNAME");
170	cols[4] = String("SCANNO");
171	}
172	uInt outrowCount = 0;
173	// use STIdxIter2 instead of TableIterator
174	STIdxIter2 iter( in[0], cols ) ;
175	// double t2 = 0 ;
176	// double t3 = 0 ;
177	// double t4 = 0 ;
178	// double t5 = 0 ;
179	// TableIterator iter(baset, cols);
180	// int count = 0 ;
181	while (!iter.pastEnd()) {
182	Vector<uInt> rows = iter.getRows( SHARE ) ;
183	if ( rows.nelements() == 0 ) {
184	iter.next() ;
185	continue ;
186	}
187	Record current = iter.currentValue() ;
188	//Table subt = iter.table();
189	// copy the first row of this selection into the new table
190	tout.addRow();
191	// t4 = mathutil::gettimeofday_sec() ;
192	// skip to copy SPECTRA, FLAGTRA, and TSYS since those heavy columns are
193	// overwritten in the following process
194	copyRows( tout, baset, outrowCount, rows[0], 1, False, False, False ) ;
195	// t5 += mathutil::gettimeofday_sec() - t4 ;
196	// re-index to 0
197	if ( avmode != "SCAN" && avmode != "SOURCE" ) {
198	scanColOut.put(outrowCount, uInt(0));
199	}
200
201	// 2012/02/17 TN
202	// Since STGrid is implemented, average doesn't consider direction
203	// when accumulating
204	// MDirection::ScalarColumn dircol ;
205	// dircol.attach( subt, "DIRECTION" ) ;
206	// Int length = subt.nrow() ;
207	// vector< Vector<Double> > dirs ;
208	// vector<int> indexes ;
209	// for ( Int i = 0 ; i < length ; i++ ) {
210	// Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
211	// //os << << count++ << ": " ;
212	// //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
213	// bool adddir = true ;
214	// for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
215	// //if ( allTrue( t == dirs[j] ) ) {
216	// Double dx = t[0] - dirs[j][0] ;
217	// Double dy = t[1] - dirs[j][1] ;
218	// Double dd = sqrt( dx * dx + dy * dy ) ;
219	// //if ( allNearAbs( t, dirs[j], tol ) ) {
220	// if ( dd <= tol ) {
221	// adddir = false ;
222	// break ;
223	// }
224	// }
225	// if ( adddir ) {
226	// dirs.push_back( t ) ;
227	// indexes.push_back( i ) ;
228	// }
229	// }
230	// uInt rowNum = dirs.size() ;
231	// tout.addRow( rowNum ) ;
232	// for ( uInt i = 0 ; i < rowNum ; i++ ) {
233	// TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
234	// // re-index to 0
235	// if ( avmode != "SCAN" && avmode != "SOURCE" ) {
236	// scanColOut.put(outrowCount+i, uInt(0));
237	// }
238	// }
239	// outrowCount += rowNum ;
240
241	// merge loop
242	uInt i = outrowCount ;
243	// in[0] is already selected by iterator
244	specCol.attach(baset,"SPECTRA");
245	flagCol.attach(baset,"FLAGTRA");
246	tsysCol.attach(baset,"TSYS");
247	intCol.attach(baset,"INTERVAL");
248	mjdCol.attach(baset,"TIME");
249	flagRowCol.attach(baset,"FLAGROW");
250	Vector<Float> spec,tsys;
251	Vector<uChar> flag;
252	Double inter,time;
253	uInt flagRow;
254
255	for (uInt l = 0; l < rows.nelements(); ++l ) {
256	uInt k = rows[l] ;
257	flagCol.get(k, flag);
258	Vector<Bool> bflag(flag.shape());
259	flagRowCol.get(k, flagRow);
260	if (flagRow > 0)
261	bflag = true;
262	else
263	convertArray(bflag, flag);
264	/*
265	if ( allEQ(bflag, True) ) {
266	continue;//don't accumulate
267	}
268	*/
269	specCol.get(k, spec);
270	tsysCol.get(k, tsys);
271	intCol.get(k, inter);
272	mjdCol.get(k, time);
273	// spectrum has to be added last to enable weighting by the other values
274	// t2 = mathutil::gettimeofday_sec() ;
275	acc.add(spec, !bflag, tsys, inter, time);
276	// t3 += mathutil::gettimeofday_sec() - t2 ;
277
278	}
279
280
281	// in[0] is already selected by TableIterator so that index is
282	// started from 1
283	for ( int j=1; j < int(in.size()); ++j ) {
284	const Table& tin = in[j]->table();
285	//const TableRecord& rec = row.get(i);
286	ROScalarColumn<Double> tmp(tin, "TIME");
287	Double td;tmp.get(0,td);
288
289	#if 1
290	static char const*const colNames1[] = { "IFNO", "BEAMNO", "POLNO" };
291	//uInt const values1[] = { rec.asuInt("IFNO"), rec.asuInt("BEAMNO"), rec.asuInt("POLNO") };
292	uInt const values1[] = { current.asuInt("IFNO"), current.asuInt("BEAMNO"), current.asuInt("POLNO") };
293	SingleTypeEqPredicate<uInt, 3> myPred(tin, colNames1, values1);
294	CustomTableExprNodeRep myNodeRep(tin, myPred);
295	myNodeRep.link(); // to avoid automatic delete when myExpr is destructed.
296	CustomTableExprNode myExpr(myNodeRep);
297	Table basesubt = tin(myExpr);
298	#else
299	// Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
300	// && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
301	// && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
302	Table basesubt = tin( tin.col("BEAMNO") == current.asuInt("BEAMNO")
303	&& tin.col("IFNO") == current.asuInt("IFNO")
304	&& tin.col("POLNO") == current.asuInt("POLNO") );
305	#endif
306	Table subt;
307	if ( avmode == "SOURCE") {
308	// subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME"));
309	subt = basesubt( basesubt.col("SRCNAME") == current.asString("SRCNAME") );
310
311	} else if (avmode == "SCAN") {
312	// subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME")
313	// && basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
314	subt = basesubt( basesubt.col("SRCNAME") == current.asString("SRCNAME")
315	&& basesubt.col("SCANNO") == current.asuInt("SCANNO") );
316	} else {
317	subt = basesubt;
318	}
319
320	// 2012/02/17 TN
321	// Since STGrid is implemented, average doesn't consider direction
322	// when accumulating
323	// vector<uInt> removeRows ;
324	// uInt nrsubt = subt.nrow() ;
325	// for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
326	// //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
327	// Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
328	// Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
329	// double dx = x0[0] - x1[0];
330	// double dy = x0[1] - x1[1];
331	// Double dd = sqrt( dx * dx + dy * dy ) ;
332	// //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
333	// if ( dd > tol ) {
334	// removeRows.push_back( irow ) ;
335	// }
336	// }
337	// if ( removeRows.size() != 0 ) {
338	// subt.removeRow( removeRows ) ;
339	// }
340
341	// if ( nrsubt == removeRows.size() )
342	// throw(AipsError("Averaging data is empty.")) ;
343
344	specCol.attach(subt,"SPECTRA");
345	flagCol.attach(subt,"FLAGTRA");
346	tsysCol.attach(subt,"TSYS");
347	intCol.attach(subt,"INTERVAL");
348	mjdCol.attach(subt,"TIME");
349	flagRowCol.attach(subt,"FLAGROW");
350	for (uInt k = 0; k < subt.nrow(); ++k ) {
351	flagCol.get(k, flag);
352	Vector<Bool> bflag(flag.shape());
353	flagRowCol.get(k, flagRow);
354	if (flagRow > 0)
355	bflag = true;
356	else
357	convertArray(bflag, flag);
358	/*
359	if ( allEQ(bflag, True) ) {
360	continue;//don't accumulate
361	}
362	*/
363	specCol.get(k, spec);
364	//tsysCol.get(k, tsys);
365	tsys.assign( tsysCol(k) );
366	intCol.get(k, inter);
367	mjdCol.get(k, time);
368	// spectrum has to be added last to enable weighting by the other values
369	// t2 = mathutil::gettimeofday_sec() ;
370	acc.add(spec, !bflag, tsys, inter, time);
371	// t3 += mathutil::gettimeofday_sec() - t2 ;
372	}
373
374	}
375	const Vector<Bool>& msk = acc.getMask();
376	uInt outFlagRow = 0;
377	if ( allEQ(msk, False) ) {
378	rowstodelB[nrowdel] = i ;
379	nrowdel++ ;
380	outFlagRow = 1;
381	}
382	//write out
383	if (acc.state()) {
384	// If there exists a channel at which all the input spectra are masked,
385	// spec has 'nan' values for that channel and it may affect the following
386	// processes. To avoid this, replacing 'nan' values in spec with
387	// weighted-mean of all spectra in the following line.
388	// (done for CAS-2776, 2011/04/07 by Wataru Kawasaki)
389	acc.replaceNaN();
390
391	Vector<uChar> flg(msk.shape());
392	convertArray(flg, !msk);
393	for (uInt k = 0; k < flg.nelements(); ++k) {
394	uChar userFlag = 1 << 7;
395	if (msk[k]==True) userFlag = 0 << 7;
396	flg(k) = userFlag;
397	}
398
399	flagColOut.put(i, flg);
400	specColOut.put(i, acc.getSpectrum());
401	tsysColOut.put(i, acc.getTsys());
402	intColOut.put(i, acc.getInterval());
403	mjdColOut.put(i, acc.getTime());
404	// we should only have one cycle now -> reset it to be 0
405	// frequency switched data has different CYCLENO for different IFNO
406	// which requires resetting this value
407	cycColOut.put(i, uInt(0));
408	flagRowColOut.put(i, outFlagRow);
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	if (nrowdel == tout.nrow()) {
422	os << LogIO::WARN << "Output data are fully flagged." << LogIO::POST;
423	}
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	Vector<Bool> m(mask);
1998	const Table& tab = in->table();
1999	ROArrayColumn<Float> specCol(tab, "SPECTRA");
2000	ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2001	std::vector<float> out;
2002	for (uInt i = 0; i < tab.nrow(); ++i) {
2003	if (in->isAllChannelsFlagged(i)) {
2004	out.push_back(NAN);
2005	} else {
2006	float outstat = 0.0;
2007	Vector<Float> spec; specCol.get(i, spec);
2008	Vector<uChar> flag; flagCol.get(i, flag);
2009	MaskedArray<Float> ma = maskedArray(spec, flag);
2010	if (spec.nelements() == m.nelements()) {
2011	outstat = mathutil::statistics(which, ma(m));
2012	} else {
2013	outstat = mathutil::statistics(which, ma);
2014	}
2015	out.push_back(outstat);
2016	}
2017	}
2018	return out;
2019	}
2020
2021	std::vector< float > STMath::statisticRow( const CountedPtr< Scantable > & in,
2022	const std::vector< bool > & mask,
2023	const std::string& which,
2024	int row )
2025	{
2026
2027	Vector<Bool> m(mask);
2028	const Table& tab = in->table();
2029	ROArrayColumn<Float> specCol(tab, "SPECTRA");
2030	ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2031	std::vector<float> out;
2032
2033	Vector<Float> spec; specCol.get(row, spec);
2034	Vector<uChar> flag; flagCol.get(row, flag);
2035	MaskedArray<Float> ma = maskedArray(spec, flag);
2036	float outstat = 0.0;
2037	if ( spec.nelements() == m.nelements() ) {
2038	outstat = mathutil::statistics(which, ma(m));
2039	} else {
2040	outstat = mathutil::statistics(which, ma);
2041	}
2042	out.push_back(outstat);
2043
2044	return out;
2045	}
2046
2047	std::vector< int > STMath::minMaxChan( const CountedPtr< Scantable > & in,
2048	const std::vector< bool > & mask,
2049	const std::string& which )
2050	{
2051
2052	Vector<Bool> m(mask);
2053	const Table& tab = in->table();
2054	ROArrayColumn<Float> specCol(tab, "SPECTRA");
2055	ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2056	std::vector<int> out;
2057	for (uInt i=0; i < tab.nrow(); ++i ) {
2058	if (in->isAllChannelsFlagged(i)) {
2059	out.push_back(NAN);
2060	} else {
2061	Vector<Float> spec; specCol.get(i, spec);
2062	Vector<uChar> flag; flagCol.get(i, flag);
2063	MaskedArray<Float> ma = maskedArray(spec, flag);
2064	if (ma.ndim() != 1) {
2065	throw (ArrayError(
2066	"std::vector<int> STMath::minMaxChan("
2067	"ContedPtr<Scantable> &in, std::vector<bool> &mask, "
2068	" std::string &which)"
2069	" - MaskedArray is not 1D"));
2070	}
2071	IPosition outpos(1,0);
2072	if ( spec.nelements() == m.nelements() ) {
2073	outpos = mathutil::minMaxPos(which, ma(m));
2074	} else {
2075	outpos = mathutil::minMaxPos(which, ma);
2076	}
2077	out.push_back(outpos[0]);
2078	}
2079	}
2080	return out;
2081	}
2082
2083	CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
2084	int width )
2085	{
2086	if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
2087	CountedPtr< Scantable > out = getScantable(in, false);
2088	Table& tout = out->table();
2089	out->frequencies().rescale(width, "BIN");
2090	ArrayColumn<Float> specCol(tout, "SPECTRA");
2091	ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
2092	ArrayColumn<Float> tsysCol(tout, "TSYS");
2093
2094	for (uInt i=0; i < tout.nrow(); ++i ) {
2095	MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
2096	MaskedArray<Float> maout;
2097	LatticeUtilities::bin(maout, main, 0, Int(width));
2098	specCol.put(i, maout.getArray());
2099	flagCol.put(i, flagsFromMA(maout));
2100	if (tsysCol(i).nelements() == specCol(i).nelements()) {
2101	MaskedArray<Float> matsysin = maskedArray(tsysCol(i), flagCol(i));
2102	MaskedArray<Float> matsysout;
2103	LatticeUtilities::bin(matsysout, matsysin, 0, Int(width));
2104	tsysCol.put(i, matsysout.getArray());
2105	}
2106	// take only the first binned spectrum's length for the deprecated
2107	// global header item nChan
2108	if (i==0) tout.rwKeywordSet().define(String("nChan"),
2109	Int(maout.getArray().nelements()));
2110	}
2111	return out;
2112	}
2113
2114	CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
2115	const std::string& method,
2116	float width )
2117	//
2118	// Should add the possibility of width being specified in km/s. This means
2119	// that for each freqID (SpectralCoordinate) we will need to convert to an
2120	// average channel width (say at the reference pixel). Then we would need
2121	// to be careful to make sure each spectrum (of different freqID)
2122	// is the same length.
2123	//
2124	{
2125	//InterpolateArray1D<Double,Float>::InterpolationMethod interp;
2126	Int interpMethod(stringToIMethod(method));
2127
2128	CountedPtr< Scantable > out = getScantable(in, false);
2129	Table& tout = out->table();
2130
2131	// Resample SpectralCoordinates (one per freqID)
2132	out->frequencies().rescale(width, "RESAMPLE");
2133	TableIterator iter(tout, "IFNO");
2134	TableRow row(tout);
2135	while ( !iter.pastEnd() ) {
2136	Table tab = iter.table();
2137	ArrayColumn<Float> specCol(tab, "SPECTRA");
2138	//ArrayColumn<Float> tsysCol(tout, "TSYS");
2139	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2140	Vector<Float> spec;
2141	Vector<uChar> flag;
2142	specCol.get(0,spec); // the number of channels should be constant per IF
2143	uInt nChanIn = spec.nelements();
2144	Vector<Float> xIn(nChanIn); indgen(xIn);
2145	Int fac = Int(nChanIn/width);
2146	Vector<Float> xOut(fac+10); // 10 to be safe - resize later
2147	uInt k = 0;
2148	Float x = 0.0;
2149	while (x < Float(nChanIn) ) {
2150	xOut(k) = x;
2151	k++;
2152	x += width;
2153	}
2154	uInt nChanOut = k;
2155	xOut.resize(nChanOut, True);
2156	// process all rows for this IFNO
2157	Vector<Float> specOut;
2158	Vector<Bool> maskOut;
2159	Vector<uChar> flagOut;
2160	for (uInt i=0; i < tab.nrow(); ++i) {
2161	specCol.get(i, spec);
2162	flagCol.get(i, flag);
2163	Vector<Bool> mask(flag.nelements());
2164	convertArray(mask, flag);
2165
2166	IPosition shapeIn(spec.shape());
2167	//sh.nchan = nChanOut;
2168	InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
2169	xIn, spec, mask,
2170	interpMethod, True, True);
2171	/// @todo do the same for channel based Tsys
2172	flagOut.resize(maskOut.nelements());
2173	convertArray(flagOut, maskOut);
2174	specCol.put(i, specOut);
2175	flagCol.put(i, flagOut);
2176	}
2177	++iter;
2178	}
2179
2180	return out;
2181	}
2182
2183	STMath::imethod STMath::stringToIMethod(const std::string& in)
2184	{
2185	static STMath::imap lookup;
2186
2187	// initialize the lookup table if necessary
2188	if ( lookup.empty() ) {
2189	lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
2190	lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
2191	lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
2192	lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
2193	}
2194
2195	STMath::imap::const_iterator iter = lookup.find(in);
2196
2197	if ( lookup.end() == iter ) {
2198	std::string message = in;
2199	message += " is not a valid interpolation mode";
2200	throw(AipsError(message));
2201	}
2202	return iter->second;
2203	}
2204
2205	WeightType STMath::stringToWeight(const std::string& in)
2206	{
2207	static std::map<std::string, WeightType> lookup;
2208
2209	// initialize the lookup table if necessary
2210	if ( lookup.empty() ) {
2211	lookup["NONE"] = asap::W_NONE;
2212	lookup["TINT"] = asap::W_TINT;
2213	lookup["TINTSYS"] = asap::W_TINTSYS;
2214	lookup["TSYS"] = asap::W_TSYS;
2215	lookup["VAR"] = asap::W_VAR;
2216	}
2217
2218	std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
2219
2220	if ( lookup.end() == iter ) {
2221	std::string message = in;
2222	message += " is not a valid weighting mode";
2223	throw(AipsError(message));
2224	}
2225	return iter->second;
2226	}
2227
2228	CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
2229	const vector< float > & coeff,
2230	const std::string & filename,
2231	const std::string& method)
2232	{
2233	LogIO os( LogOrigin( "STMath", "gainElevation", WHERE ) ) ;
2234	// Get elevation data from Scantable and convert to degrees
2235	CountedPtr< Scantable > out = getScantable(in, false);
2236	Table& tab = out->table();
2237	ROScalarColumn<Float> elev(tab, "ELEVATION");
2238	Vector<Float> x = elev.getColumn();
2239	x *= Float(180 / C::pi); // Degrees
2240
2241	Vector<Float> coeffs(coeff);
2242	const uInt nc = coeffs.nelements();
2243	if ( filename.length() > 0 && nc > 0 ) {
2244	throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
2245	}
2246
2247	// Correct
2248	if ( nc > 0 \|\| filename.length() == 0 ) {
2249	// Find instrument
2250	Bool throwit = True;
2251	Instrument inst =
2252	STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2253	throwit);
2254
2255	// Set polynomial
2256	Polynomial<Float>* ppoly = 0;
2257	Vector<Float> coeff;
2258	String msg;
2259	if ( nc > 0 ) {
2260	ppoly = new Polynomial<Float>(nc-1);
2261	coeff = coeffs;
2262	msg = String("user");
2263	} else {
2264	STAttr sdAttr;
2265	coeff = sdAttr.gainElevationPoly(inst);
2266	ppoly = new Polynomial<Float>(coeff.nelements()-1);
2267	msg = String("built in");
2268	}
2269
2270	if ( coeff.nelements() > 0 ) {
2271	ppoly->setCoefficients(coeff);
2272	} else {
2273	delete ppoly;
2274	throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
2275	}
2276	os << "Making polynomial correction with " << msg << " coefficients:" << endl;
2277	os << " " << coeff << LogIO::POST;
2278	const uInt nrow = tab.nrow();
2279	Vector<Float> factor(nrow);
2280	for ( uInt i=0; i < nrow; ++i ) {
2281	factor[i] = 1.0 / (*ppoly)(x[i]);
2282	}
2283	delete ppoly;
2284	scaleByVector(tab, factor, true);
2285
2286	} else {
2287	// Read and correct
2288	os << "Making correction from ascii Table" << LogIO::POST;
2289	scaleFromAsciiTable(tab, filename, method, x, true);
2290	}
2291	return out;
2292	}
2293
2294	void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
2295	const std::string& method,
2296	const Vector<Float>& xout, bool dotsys)
2297	{
2298
2299	// Read gain-elevation ascii file data into a Table.
2300
2301	String formatString;
2302	Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
2303	scaleFromTable(in, tbl, method, xout, dotsys);
2304	}
2305
2306	void STMath::scaleFromTable(Table& in,
2307	const Table& table,
2308	const std::string& method,
2309	const Vector<Float>& xout, bool dotsys)
2310	{
2311
2312	ROScalarColumn<Float> geElCol(table, "ELEVATION");
2313	ROScalarColumn<Float> geFacCol(table, "FACTOR");
2314	Vector<Float> xin = geElCol.getColumn();
2315	Vector<Float> yin = geFacCol.getColumn();
2316	Vector<Bool> maskin(xin.nelements(),True);
2317
2318	// Interpolate (and extrapolate) with desired method
2319
2320	InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2321
2322	Vector<Float> yout;
2323	Vector<Bool> maskout;
2324	InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
2325	xin, yin, maskin, interp,
2326	True, True);
2327
2328	scaleByVector(in, Float(1.0)/yout, dotsys);
2329	}
2330
2331	void STMath::scaleByVector( Table& in,
2332	const Vector< Float >& factor,
2333	bool dotsys )
2334	{
2335	uInt nrow = in.nrow();
2336	if ( factor.nelements() != nrow ) {
2337	throw(AipsError("factors.nelements() != table.nelements()"));
2338	}
2339	ArrayColumn<Float> specCol(in, "SPECTRA");
2340	ArrayColumn<uChar> flagCol(in, "FLAGTRA");
2341	ArrayColumn<Float> tsysCol(in, "TSYS");
2342	for (uInt i=0; i < nrow; ++i) {
2343	MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2344	ma *= factor[i];
2345	specCol.put(i, ma.getArray());
2346	flagCol.put(i, flagsFromMA(ma));
2347	if ( dotsys ) {
2348	Vector<Float> tsys = tsysCol(i);
2349	tsys *= factor[i];
2350	tsysCol.put(i,tsys);
2351	}
2352	}
2353	}
2354
2355	CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
2356	float d, float etaap,
2357	float jyperk )
2358	{
2359	LogIO os( LogOrigin( "STMath", "convertFlux", WHERE ) ) ;
2360
2361	CountedPtr< Scantable > out = getScantable(in, false);
2362	Table& tab = in->table();
2363	Table& outtab = out->table();
2364	Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
2365	Unit K(String("K"));
2366	Unit JY(String("Jy"));
2367
2368	bool tokelvin = true;
2369	Double cfac = 1.0;
2370
2371	if ( fluxUnit == JY ) {
2372	os << "Converting to K" << LogIO::POST;
2373	Quantum<Double> t(1.0,fluxUnit);
2374	Quantum<Double> t2 = t.get(JY);
2375	cfac = (t2 / t).getValue(); // value to Jy
2376
2377	tokelvin = true;
2378	out->setFluxUnit("K");
2379	} else if ( fluxUnit == K ) {
2380	os << "Converting to Jy" << LogIO::POST;
2381	Quantum<Double> t(1.0,fluxUnit);
2382	Quantum<Double> t2 = t.get(K);
2383	cfac = (t2 / t).getValue(); // value to K
2384
2385	tokelvin = false;
2386	out->setFluxUnit("Jy");
2387	} else {
2388	throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
2389	}
2390	// Make sure input values are converted to either Jy or K first...
2391	Float factor = cfac;
2392
2393	// Select method
2394	if (jyperk > 0.0) {
2395	factor *= jyperk;
2396	if ( tokelvin ) factor = 1.0 / jyperk;
2397	os << "Jy/K = " << jyperk << LogIO::POST;
2398	Vector<Float> factors(outtab.nrow(), factor);
2399	scaleByVector(outtab,factors, false);
2400	} else if ( etaap > 0.0) {
2401	if (d < 0) {
2402	Instrument inst =
2403	STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2404	True);
2405	STAttr sda;
2406	d = sda.diameter(inst);
2407	}
2408	jyperk = STAttr::findJyPerK(etaap, d);
2409	os << "Jy/K = " << jyperk << LogIO::POST;
2410	factor *= jyperk;
2411	if ( tokelvin ) {
2412	factor = 1.0 / factor;
2413	}
2414	Vector<Float> factors(outtab.nrow(), factor);
2415	scaleByVector(outtab, factors, False);
2416	} else {
2417
2418	// OK now we must deal with automatic look up of values.
2419	// We must also deal with the fact that the factors need
2420	// to be computed per IF and may be different and may
2421	// change per integration.
2422
2423	os <<"Looking up conversion factors" << LogIO::POST;
2424	convertBrightnessUnits(out, tokelvin, cfac);
2425	}
2426
2427	return out;
2428	}
2429
2430	void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
2431	bool tokelvin, float cfac )
2432	{
2433	Table& table = in->table();
2434	Instrument inst =
2435	STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
2436	TableIterator iter(table, "FREQ_ID");
2437	STFrequencies stfreqs = in->frequencies();
2438	STAttr sdAtt;
2439	while (!iter.pastEnd()) {
2440	Table tab = iter.table();
2441	ArrayColumn<Float> specCol(tab, "SPECTRA");
2442	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2443	ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
2444	MEpoch::ROScalarColumn timeCol(tab, "TIME");
2445
2446	uInt freqid; freqidCol.get(0, freqid);
2447	Vector<Float> tmpspec; specCol.get(0, tmpspec);
2448	// STAttr.JyPerK has a Vector interface... change sometime.
2449	Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
2450	for ( uInt i=0; i<tab.nrow(); ++i) {
2451	Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
2452	Float factor = cfac * jyperk;
2453	if ( tokelvin ) factor = Float(1.0) / factor;
2454	MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2455	ma *= factor;
2456	specCol.put(i, ma.getArray());
2457	flagCol.put(i, flagsFromMA(ma));
2458	}
2459	++iter;
2460	}
2461	}
2462
2463	CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
2464	const std::vector<float>& tau )
2465	{
2466	CountedPtr< Scantable > out = getScantable(in, false);
2467
2468	Table outtab = out->table();
2469
2470	const Int ntau = uInt(tau.size());
2471	std::vector<float>::const_iterator tauit = tau.begin();
2472	AlwaysAssert((ntau == 1 \|\| ntau == in->nif() \|\| ntau == in->nif() * in->npol()),
2473	AipsError);
2474	TableIterator iiter(outtab, "IFNO");
2475	while ( !iiter.pastEnd() ) {
2476	Table itab = iiter.table();
2477	TableIterator piter(itab, "POLNO");
2478	while ( !piter.pastEnd() ) {
2479	Table tab = piter.table();
2480	ROScalarColumn<Float> elev(tab, "ELEVATION");
2481	ArrayColumn<Float> specCol(tab, "SPECTRA");
2482	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2483	ArrayColumn<Float> tsysCol(tab, "TSYS");
2484	for ( uInt i=0; i<tab.nrow(); ++i) {
2485	Float zdist = Float(C::pi_2) - elev(i);
2486	Float factor = exp(*tauit/cos(zdist));
2487	MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2488	ma *= factor;
2489	specCol.put(i, ma.getArray());
2490	flagCol.put(i, flagsFromMA(ma));
2491	Vector<Float> tsys;
2492	tsysCol.get(i, tsys);
2493	tsys *= factor;
2494	tsysCol.put(i, tsys);
2495	}
2496	if (ntau == in->nif()*in->npol() ) {
2497	tauit++;
2498	}
2499	piter++;
2500	}
2501	if (ntau >= in->nif() ) {
2502	tauit++;
2503	}
2504	iiter++;
2505	}
2506	return out;
2507	}
2508
2509	CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
2510	const std::string& kernel,
2511	float width, int order)
2512	{
2513	CountedPtr< Scantable > out = getScantable(in, false);
2514	Table table = out->table();
2515
2516	TableIterator iter(table, "IFNO");
2517	while (!iter.pastEnd()) {
2518	Table tab = iter.table();
2519	ArrayColumn<Float> specCol(tab, "SPECTRA");
2520	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2521	ScalarColumn<uInt> flagrowCol(tab, "FLAGROW");
2522	Vector<Float> spec;
2523	Vector<uChar> flag;
2524	Vector<uInt> flagrow = flagrowCol.getColumn();
2525	for (uInt i = 0; i < tab.nrow(); ++i) {
2526	if (flagrow[i] != 0) {
2527	// do not process flagged row
2528	continue;
2529	}
2530	specCol.get(i, spec);
2531	flagCol.get(i, flag);
2532	Vector<Bool> mask(flag.nelements());
2533	convertArray(mask, flag);
2534	Vector<Float> specout;
2535	Vector<Bool> maskout;
2536	if (kernel == "hanning") {
2537	mathutil::hanning(specout, maskout, spec, !mask);
2538	} else if (kernel == "rmedian") {
2539	mathutil::runningMedian(specout, maskout, spec , mask, width);
2540	} else if (kernel == "poly") {
2541	mathutil::polyfit(specout, maskout, spec, !mask, width, order);
2542	}
2543
2544	for (uInt j = 0; j < flag.nelements(); ++j) {
2545	uChar userFlag = 1 << 7;
2546	if (maskout[j]==True) userFlag = 0 << 7;
2547	flag(j) = userFlag;
2548	}
2549
2550	flagCol.put(i, flag);
2551	specCol.put(i, specout);
2552	}
2553	++iter;
2554	}
2555	return out;
2556	}
2557
2558	CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
2559	const std::string& kernel, float width,
2560	int order)
2561	{
2562	if (kernel == "rmedian" \|\| kernel == "hanning" \|\| kernel == "poly") {
2563	return smoothOther(in, kernel, width, order);
2564	}
2565	CountedPtr< Scantable > out = getScantable(in, false);
2566	Table& table = out->table();
2567	VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
2568	// same IFNO should have same no of channels
2569	// this saves overhead
2570	TableIterator iter(table, "IFNO");
2571	while (!iter.pastEnd()) {
2572	Table tab = iter.table();
2573	ArrayColumn<Float> specCol(tab, "SPECTRA");
2574	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2575	ScalarColumn<uInt> flagrowCol(tab, "FLAGROW");
2576	Vector<Float> spec = specCol( 0 );
2577	uInt nchan = spec.nelements();
2578	Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
2579	Convolver<Float> conv(kvec, IPosition(1,nchan));
2580	Vector<uChar> flag;
2581	Vector<Bool> mask(nchan);
2582	Vector<uInt> flagrow = flagrowCol.getColumn();
2583	for ( uInt i=0; i<tab.nrow(); ++i) {
2584	if (flagrow[i] != 0) {
2585	// do not process flagged row
2586	continue;
2587	}
2588
2589	specCol.get(i, spec);
2590	flagCol.get(i, flag);
2591	convertArray(mask, flag);
2592	Vector<Float> specout;
2593	//mathutil::replaceMaskByZero(specout, mask);
2594	mathutil::replaceMaskByZero(spec, !mask);
2595	//std::vector<bool> vmask;
2596	//(!mask).tovector(vmask);
2597	//mathutil::doZeroOrderInterpolation(spec, vmask);
2598	conv.linearConv(specout, spec);
2599	specCol.put(i, specout);
2600	}
2601	++iter;
2602	}
2603	return out;
2604	}
2605
2606	CountedPtr< Scantable >
2607	STMath::merge( const std::vector< CountedPtr < Scantable > >& in,
2608	const std::string &freqTol )
2609	{
2610	Double freqTolInHz = 0.0; // default is 0.0Hz (merge only when exact match)
2611	if (freqTol.size() > 0) {
2612	Quantum<Double> freqTolInQuantity;
2613	if (!Quantum<Double>::read(freqTolInQuantity, freqTol)) {
2614	throw(AipsError("Failed to convert freqTol string to quantity"));
2615	}
2616	if (!freqTolInQuantity.isConform("Hz")) {
2617	throw(AipsError("Invalid freqTol string"));
2618	}
2619	freqTolInHz = freqTolInQuantity.getValue("Hz");
2620	LogIO os(LogOrigin("STMath", "merge", WHERE));
2621	os << "frequency tolerance = " << freqTolInHz << "Hz" << LogIO::POST;
2622	}
2623
2624	if ( in.size() < 2 ) {
2625	throw(AipsError("Need at least two scantables to perform a merge."));
2626	}
2627	std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
2628	bool insitu = insitu_;
2629	setInsitu(false);
2630	CountedPtr< Scantable > out = getScantable(*it, false);
2631	setInsitu(insitu);
2632	Table& tout = out->table();
2633	ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
2634	ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
2635	ScalarColumn<uInt> ifnocol(tout, "IFNO");
2636	// Renumber SCANNO to be 0-based
2637	Vector<uInt> scannos = scannocol.getColumn();
2638	uInt offset = min(scannos);
2639	scannos -= offset;
2640	scannocol.putColumn(scannos);
2641	uInt newscanno = max(scannos)+1;
2642	++it;
2643
2644	// new IFNO
2645	uInt ifnoCounter = max(ifnocol.getColumn()) + 1;
2646
2647	// Here we assume that each IFNO has unique MOLECULE_ID
2648	// molIdMap:
2649	// KEY: IFNO
2650	// VALUE: MOLECULE_ID
2651	map<uInt, uInt> molIdMap;
2652	{
2653	TableIterator ifit(tout, "IFNO");
2654	while (!ifit.pastEnd()) {
2655	ROTableRow row(ifit.table());
2656	const TableRecord& rec = row.get(0);
2657	molIdMap[rec.asuInt("IFNO")] = rec.asuInt("MOLECULE_ID");
2658	ifit.next();
2659	}
2660	}
2661
2662	while ( it != in.end() ){
2663	// Check FREQUENCIES/BASEFRAME
2664	if ( out->frequencies().getFrame(true) != (*it)->frequencies().getFrame(true) ) {
2665	throw(AipsError("BASEFRAME is not identical"));
2666	}
2667
2668	if ( ! (it)->conformant(out) ) {
2669	// non conformant.
2670	LogIO os( LogOrigin( "STMath", "merge()", WHERE ) ) ;
2671	os << LogIO::SEVERE << "Can't merge scantables as header informations (any one of AntennaName, Equinox, and FluxUnit) differ." << LogIO::EXCEPTION ;
2672	}
2673	out->appendToHistoryTable((*it)->history());
2674	const Table& tab = (*it)->table();
2675
2676	Block<String> cols(3);
2677	cols[0] = String("FREQ_ID");
2678	cols[1] = String("MOLECULE_ID");
2679	cols[2] = String("FOCUS_ID");
2680
2681	TableIterator scanit(tab, "SCANNO");
2682	while (!scanit.pastEnd()) {
2683	TableIterator subit(scanit.table(), cols);
2684	while ( !subit.pastEnd() ) {
2685	uInt nrow = tout.nrow();
2686	Table thetab = subit.table();
2687	ROTableRow row(thetab);
2688	Vector<uInt> thecolvals(thetab.nrow());
2689	// The selected subset of table should have
2690	// the equal FREQ_ID, MOLECULE_ID, and FOCUS_ID values.
2691	const TableRecord& rec = row.get(0);
2692
2693	// Set the proper FREQ_ID
2694	Double rv,rp,inc;
2695	(*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
2696	uInt id;
2697
2698	// default value is new unique IFNO
2699	uInt newifno = ifnoCounter;
2700	Bool isIfMerged = False;
2701	uInt nchan = rec.asArrayFloat("SPECTRA").shape()[0];
2702	//id = out->frequencies().addEntry(rp, rv, inc);
2703	if ( !out->frequencies().match(rp, rv, inc, freqTolInHz, id) ) {
2704	// add new entry to FREQUENCIES table
2705	id = out->frequencies().addEntry(rp, rv, inc);
2706
2707	// increment counter for IFNO
2708	ifnoCounter++;
2709	}
2710	else {
2711	// should renumber IFNO to be same as existing rows that have same FREQ_ID
2712	LogIO os(LogOrigin("STMath", "merge", WHERE));
2713	Table outFreqIdSelected = tout(tout.col("FREQ_ID") == id);
2714	TableIterator _iter(outFreqIdSelected, "IFNO");
2715	map<uInt, uInt> nchanMap;
2716	while (!_iter.pastEnd()) {
2717	const Table _table = _iter.table();
2718	ROTableRow _row(_table);
2719	const TableRecord &_rec = _row.get(0);
2720	uInt nchan = _rec.asArrayFloat("SPECTRA").shape()[0];
2721	if (nchanMap.find(nchan) != nchanMap.end()) {
2722	throw(AipsError("There are non-unique IFNOs assigned to spectra that have same FREQ_ID and same nchan. Something wrong."));
2723	}
2724	nchanMap[nchan] = _rec.asuInt("IFNO");
2725	_iter.next();
2726	}
2727
2728	os << LogIO::DEBUGGING << "nchanMap for " << id << ":" << LogIO::POST;
2729	for (map<uInt, uInt>::iterator i = nchanMap.begin(); i != nchanMap.end(); ++i) {
2730	os << LogIO::DEBUGGING << "nchanMap[" << i->first << "] = " << i->second << LogIO::POST;
2731	}
2732
2733	if (nchanMap.find(nchan) == nchanMap.end()) {
2734	// increment counter for IFNO
2735	ifnoCounter++;
2736	}
2737	else {
2738	// renumber IFNO to be same as existing value that corresponds to nchan
2739	newifno = nchanMap[nchan];
2740	isIfMerged = True;
2741	}
2742	os << LogIO::DEBUGGING << "newifno = " << newifno << LogIO::POST;
2743	}
2744	// copy rows to output table
2745	tout.addRow(thetab.nrow());
2746	TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
2747
2748	Slicer slice(IPosition(1, nrow), IPosition(1, thetab.nrow()),
2749	Slicer::endIsLength);
2750
2751	thecolvals = id;
2752	freqidcol.putColumnRange(slice, thecolvals);
2753
2754	thecolvals = newifno;
2755	ifnocol.putColumnRange(slice, thecolvals);
2756
2757	// Set the proper MOLECULE_ID
2758	Vector<String> name,fname;Vector<Double> rf;
2759	(*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
2760	id = out->molecules().addEntry(rf, name, fname);
2761	if (molIdMap.find(newifno) == molIdMap.end()) {
2762	// add new entry to molIdMap
2763	molIdMap[newifno] = id;
2764	}
2765	if (isIfMerged) {
2766	id = molIdMap[newifno];
2767	}
2768	thecolvals = id;
2769	molidcol.putColumnRange(slice, thecolvals);
2770
2771	// Set the proper FOCUS_ID
2772	Float fpa,frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
2773	(*it)->focus().getEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2774	fxy, fxyp, rec.asuInt("FOCUS_ID"));
2775	id = out->focus().addEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2776	fxy, fxyp);
2777	thecolvals = id;
2778	focusidcol.putColumnRange(slice, thecolvals);
2779
2780	// Set the proper SCANNO
2781	thecolvals = newscanno;
2782	scannocol.putColumnRange(slice, thecolvals);
2783
2784	++subit;
2785	}
2786	++newscanno;
2787	++scanit;
2788	}
2789	++it;
2790	}
2791	return out;
2792	}
2793
2794	CountedPtr< Scantable >
2795	STMath::invertPhase( const CountedPtr < Scantable >& in )
2796	{
2797	return applyToPol(in, &STPol::invertPhase, Float(0.0));
2798	}
2799
2800	CountedPtr< Scantable >
2801	STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
2802	{
2803	return applyToPol(in, &STPol::rotatePhase, Float(phase));
2804	}
2805
2806	CountedPtr< Scantable >
2807	STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
2808	{
2809	return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
2810	}
2811
2812	CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
2813	STPol::polOperation fptr,
2814	Float phase )
2815	{
2816	CountedPtr< Scantable > out = getScantable(in, false);
2817	Table& tout = out->table();
2818	Block<String> cols(4);
2819	cols[0] = String("SCANNO");
2820	cols[1] = String("BEAMNO");
2821	cols[2] = String("IFNO");
2822	cols[3] = String("CYCLENO");
2823	TableIterator iter(tout, cols);
2824	CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
2825	out->getPolType() );
2826	while (!iter.pastEnd()) {
2827	Table t = iter.table();
2828	ArrayColumn<Float> speccol(t, "SPECTRA");
2829	ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
2830	Matrix<Float> pols(speccol.getColumn());
2831	try {
2832	stpol->setSpectra(pols);
2833	Float fang,fhand;
2834	fang = in->focusTable_.getTotalAngle(focidcol(0));
2835	fhand = in->focusTable_.getFeedHand(focidcol(0));
2836	stpol->setPhaseCorrections(fang, fhand);
2837	// use a member function pointer in STPol. This only works on
2838	// the STPol pointer itself, not the Counted Pointer so
2839	// derefernce it.
2840	(&((stpol))->fptr)(phase);
2841	speccol.putColumn(stpol->getSpectra());
2842	} catch (AipsError& e) {
2843	//delete stpol;stpol=0;
2844	throw(e);
2845	}
2846	++iter;
2847	}
2848	//delete stpol;stpol=0;
2849	return out;
2850	}
2851
2852	CountedPtr< Scantable >
2853	STMath::swapPolarisations( const CountedPtr< Scantable > & in )
2854	{
2855	CountedPtr< Scantable > out = getScantable(in, false);
2856	Table& tout = out->table();
2857	Table t0 = tout(tout.col("POLNO") == 0);
2858	Table t1 = tout(tout.col("POLNO") == 1);
2859	if ( t0.nrow() != t1.nrow() )
2860	throw(AipsError("Inconsistent number of polarisations"));
2861	ArrayColumn<Float> speccol0(t0, "SPECTRA");
2862	ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
2863	ArrayColumn<Float> speccol1(t1, "SPECTRA");
2864	ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
2865	Matrix<Float> s0 = speccol0.getColumn();
2866	Matrix<uChar> f0 = flagcol0.getColumn();
2867	speccol0.putColumn(speccol1.getColumn());
2868	flagcol0.putColumn(flagcol1.getColumn());
2869	speccol1.putColumn(s0);
2870	flagcol1.putColumn(f0);
2871	return out;
2872	}
2873
2874	CountedPtr< Scantable >
2875	STMath::averagePolarisations( const CountedPtr< Scantable > & in,
2876	const std::vector<bool>& mask,
2877	const std::string& weight )
2878	{
2879	if (in->npol() < 2 )
2880	throw(AipsError("averagePolarisations can only be applied to two or more"
2881	"polarisations"));
2882	bool insitu = insitu_;
2883	setInsitu(false);
2884	CountedPtr< Scantable > pols = getScantable(in, true);
2885	setInsitu(insitu);
2886	Table& tout = pols->table();
2887	std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
2888	Table tab = tableCommand(taql, in->table());
2889	if (tab.nrow() == 0 )
2890	throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
2891	TableCopy::copyRows(tout, tab);
2892	TableVector<uInt> vec(tout, "POLNO");
2893	vec = 0;
2894	pols->table_.rwKeywordSet().define("nPol", Int(1));
2895	pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
2896	//pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
2897	std::vector<CountedPtr<Scantable> > vpols;
2898	vpols.push_back(pols);
2899	CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
2900	return out;
2901	}
2902
2903	CountedPtr< Scantable >
2904	STMath::averageBeams( const CountedPtr< Scantable > & in,
2905	const std::vector<bool>& mask,
2906	const std::string& weight )
2907	{
2908	bool insitu = insitu_;
2909	setInsitu(false);
2910	CountedPtr< Scantable > beams = getScantable(in, false);
2911	setInsitu(insitu);
2912	Table& tout = beams->table();
2913	// give all rows the same BEAMNO
2914	TableVector<uInt> vec(tout, "BEAMNO");
2915	vec = 0;
2916	beams->table_.rwKeywordSet().define("nBeam", Int(1));
2917	std::vector<CountedPtr<Scantable> > vbeams;
2918	vbeams.push_back(beams);
2919	CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
2920	return out;
2921	}
2922
2923
2924	CountedPtr< Scantable >
2925	asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
2926	const std::string & refTime,
2927	const std::string & method)
2928	{
2929	LogIO os( casa::LogOrigin("STMath", "frequencyAlign()", WHERE));
2930	// clone as this is not working insitu
2931	bool insitu = insitu_;
2932	setInsitu(false);
2933	CountedPtr< Scantable > out = getScantable(in, false);
2934	setInsitu(insitu);
2935	Table& tout = out->table();
2936	// Get reference Epoch to time of first row or given String
2937	Unit DAY(String("d"));
2938	MEpoch::Ref epochRef(in->getTimeReference());
2939	MEpoch refEpoch;
2940	if (refTime.length()>0) {
2941	Quantum<Double> qt;
2942	if (MVTime::read(qt,refTime)) {
2943	MVEpoch mv(qt);
2944	refEpoch = MEpoch(mv, epochRef);
2945	} else {
2946	throw(AipsError("Invalid format for Epoch string"));
2947	}
2948	} else {
2949	refEpoch = in->timeCol_(0);
2950	}
2951	MPosition refPos = in->getAntennaPosition();
2952
2953	InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2954	/*
2955	// Comment from MV.
2956	// the following code has been commented out because different FREQ_IDs have to be aligned together even
2957	// if the frame doesn't change. So far, lack of this check didn't cause any problems.
2958	// test if user frame is different to base frame
2959	if ( in->frequencies().getFrameString(true)
2960	== in->frequencies().getFrameString(false) ) {
2961	throw(AipsError("Can't convert as no output frame has been set"
2962	" (use set_freqframe) or it is aligned already."));
2963	}
2964	*/
2965	MFrequency::Types system = in->frequencies().getFrame();
2966	MVTime mvt(refEpoch.getValue());
2967	String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
2968	os << "Aligning at reference Epoch " << epochout
2969	<< " in frame " << MFrequency::showType(system) << LogIO::POST;
2970	// set up the iterator
2971	Block<String> cols(4);
2972	// select by constant direction
2973	cols[0] = String("SRCNAME");
2974	cols[1] = String("BEAMNO");
2975	// select by IF ( no of channels varies over this )
2976	cols[2] = String("IFNO");
2977	// select by restfrequency
2978	cols[3] = String("MOLECULE_ID");
2979	TableIterator iter(tout, cols);
2980	while ( !iter.pastEnd() ) {
2981	Table t = iter.table();
2982	ROScalarColumn<String> snCol(t, "SRCNAME");
2983	os << "Aligning to position of source '" << snCol(0) << "'" << LogIO::POST;
2984	MDirection::ROScalarColumn dirCol(t, "DIRECTION");
2985	TableIterator fiter(t, "FREQ_ID");
2986	// determine nchan from the first row. This should work as
2987	// we are iterating over BEAMNO and IFNO
2988	// we should have constant direction
2989
2990	ROArrayColumn<Float> sCol(t, "SPECTRA");
2991	const MDirection direction = dirCol(0);
2992	const uInt nchan = sCol(0).nelements();
2993
2994	// skip operations if there is nothing to align
2995	if (fiter.pastEnd()) {
2996	continue;
2997	}
2998
2999	Table ftab = fiter.table();
3000	// align all frequency ids with respect to the first encountered id
3001	ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
3002	// get the SpectralCoordinate for the freqid, which we are iterating over
3003	SpectralCoordinate sC = \
3004	in->frequencies().getSpectralCoordinate(freqidCol(0));
3005	FrequencyAligner<Float> fa( sC, nchan, refEpoch,
3006	direction, refPos, system );
3007	// realign the SpectralCoordinate and put into the output Scantable
3008	Vector<String> units(1);
3009	units = String("Hz");
3010	Bool linear=True;
3011	SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
3012	sc2.setWorldAxisUnits(units);
3013	const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
3014	sc2.referenceValue()[0],
3015	sc2.increment()[0]);
3016	while ( !fiter.pastEnd() ) {
3017
3018	ftab = fiter.table();
3019	// spectral coordinate for the current FREQ_ID
3020	ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
3021	sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
3022	// create the "global" abcissa for alignment with same FREQ_ID
3023	Vector<Double> abc(nchan);
3024	for (uInt i=0; i<nchan; i++) {
3025	Double w;
3026	sC.toWorld(w,Double(i));
3027	abc[i] = w;
3028	}
3029	TableVector<uInt> tvec(ftab, "FREQ_ID");
3030	// assign new frequency id to all rows
3031	tvec = id;
3032	// cache abcissa for same time stamps, so iterate over those
3033	TableIterator timeiter(ftab, "TIME");
3034	while ( !timeiter.pastEnd() ) {
3035	Table tab = timeiter.table();
3036	ArrayColumn<Float> specCol(tab, "SPECTRA");
3037	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
3038	MEpoch::ROScalarColumn timeCol(tab, "TIME");
3039	// use align abcissa cache after the first row
3040	// these rows should be just be POLNO
3041	bool first = true;
3042	for (int i=0; i<int(tab.nrow()); ++i) {
3043	// input values
3044	Vector<uChar> flag = flagCol(i);
3045	Vector<Bool> mask(flag.shape());
3046	Vector<Float> specOut, spec;
3047	spec = specCol(i);
3048	Vector<Bool> maskOut;Vector<uChar> flagOut;
3049	convertArray(mask, flag);
3050	// alignment
3051	Bool ok = fa.align(specOut, maskOut, abc, spec,
3052	mask, timeCol(i), !first,
3053	interp, False);
3054	(void) ok; // unused stop compiler nagging
3055	// back into scantable
3056	flagOut.resize(maskOut.nelements());
3057	convertArray(flagOut, maskOut);
3058	flagCol.put(i, flagOut);
3059	specCol.put(i, specOut);
3060	// start abcissa caching
3061	first = false;
3062	}
3063	// next timestamp
3064	++timeiter;
3065	}
3066	// next FREQ_ID
3067	++fiter;
3068	}
3069	// next aligner
3070	++iter;
3071	}
3072	// set this afterwards to ensure we are doing insitu correctly.
3073	out->frequencies().setFrame(system, true);
3074	return out;
3075	}
3076
3077	CountedPtr<Scantable>
3078	asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
3079	const std::string & newtype )
3080	{
3081	if (in->npol() != 2 && in->npol() != 4)
3082	throw(AipsError("Can only convert two or four polarisations."));
3083	if ( in->getPolType() == newtype )
3084	throw(AipsError("No need to convert."));
3085	if ( ! in->selector_.empty() )
3086	throw(AipsError("Can only convert whole scantable. Unset the selection."));
3087	bool insitu = insitu_;
3088	setInsitu(false);
3089	CountedPtr< Scantable > out = getScantable(in, true);
3090	setInsitu(insitu);
3091	Table& tout = out->table();
3092	tout.rwKeywordSet().define("POLTYPE", String(newtype));
3093
3094	Block<String> cols(4);
3095	cols[0] = "SCANNO";
3096	cols[1] = "CYCLENO";
3097	cols[2] = "BEAMNO";
3098	cols[3] = "IFNO";
3099	TableIterator it(in->originalTable_, cols);
3100	String basetype = in->getPolType();
3101	STPol* stpol = STPol::getPolClass(in->factories_, basetype);
3102	try {
3103	while ( !it.pastEnd() ) {
3104	Table tab = it.table();
3105	uInt row = tab.rowNumbers()[0];
3106	stpol->setSpectra(in->getPolMatrix(row));
3107	Float fang,fhand;
3108	fang = in->focusTable_.getTotalAngle(in->mfocusidCol_(row));
3109	fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
3110	stpol->setPhaseCorrections(fang, fhand);
3111	Int npolout = 0;
3112	for (uInt i=0; i<tab.nrow(); ++i) {
3113	Vector<Float> outvec = stpol->getSpectrum(i, newtype);
3114	if ( outvec.nelements() > 0 ) {
3115	tout.addRow();
3116	TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
3117	ArrayColumn<Float> sCol(tout,"SPECTRA");
3118	ScalarColumn<uInt> pCol(tout,"POLNO");
3119	sCol.put(tout.nrow()-1 ,outvec);
3120	pCol.put(tout.nrow()-1 ,uInt(npolout));
3121	npolout++;
3122	}
3123	}
3124	tout.rwKeywordSet().define("nPol", npolout);
3125	++it;
3126	}
3127	} catch (AipsError& e) {
3128	delete stpol;
3129	throw(e);
3130	}
3131	delete stpol;
3132	return out;
3133	}
3134
3135	CountedPtr< Scantable >
3136	asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
3137	const std::string & scantype )
3138	{
3139	bool insitu = insitu_;
3140	setInsitu(false);
3141	CountedPtr< Scantable > out = getScantable(in, true);
3142	setInsitu(insitu);
3143	Table& tout = out->table();
3144	std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
3145	if (scantype == "on") {
3146	taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
3147	}
3148	Table tab = tableCommand(taql, in->table());
3149	TableCopy::copyRows(tout, tab);
3150	if (scantype == "on") {
3151	// re-index SCANNO to 0
3152	TableVector<uInt> vec(tout, "SCANNO");
3153	vec = 0;
3154	}
3155	return out;
3156	}
3157
3158	std::vector<float>
3159	asap::STMath::fft( const casa::CountedPtr< Scantable > & in,
3160	const std::vector<int>& whichrow,
3161	bool getRealImag )
3162	{
3163	std::vector<float> res;
3164	Table tab = in->table();
3165	std::vector<bool> mask;
3166
3167	if (whichrow.size() < 1) { // for all rows (by default)
3168	int nrow = int(tab.nrow());
3169	for (int i = 0; i < nrow; ++i) {
3170	res = in->execFFT(i, mask, getRealImag);
3171	}
3172	} else { // for specified rows
3173	for (uInt i = 0; i < whichrow.size(); ++i) {
3174	res = in->execFFT(i, mask, getRealImag);
3175	}
3176	}
3177
3178	return res;
3179	}
3180
3181
3182	CountedPtr<Scantable>
3183	asap::STMath::lagFlag( const CountedPtr<Scantable>& in,
3184	double start, double end,
3185	const std::string& mode )
3186	{
3187	CountedPtr<Scantable> out = getScantable(in, false);
3188	Table& tout = out->table();
3189	TableIterator iter(tout, "FREQ_ID");
3190	FFTServer<Float,Complex> ffts;
3191
3192	while ( !iter.pastEnd() ) {
3193	Table tab = iter.table();
3194	Double rp,rv,inc;
3195	ROTableRow row(tab);
3196	const TableRecord& rec = row.get(0);
3197	uInt freqid = rec.asuInt("FREQ_ID");
3198	out->frequencies().getEntry(rp, rv, inc, freqid);
3199	ArrayColumn<Float> specCol(tab, "SPECTRA");
3200	ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
3201
3202	for (int i=0; i<int(tab.nrow()); ++i) {
3203	Vector<Float> spec = specCol(i);
3204	Vector<uChar> flag = flagCol(i);
3205	std::vector<bool> mask;
3206	for (uInt j = 0; j < flag.nelements(); ++j) {
3207	mask.push_back(!(flag[j]>0));
3208	}
3209	mathutil::doZeroOrderInterpolation(spec, mask);
3210
3211	Vector<Complex> lags;
3212	ffts.fft0(lags, spec);
3213
3214	Int lag0(start+0.5);
3215	Int lag1(end+0.5);
3216	if (mode == "frequency") {
3217	lag0 = Int(spec.nelements()*abs(inc)/(start)+0.5);
3218	lag1 = Int(spec.nelements()*abs(inc)/(end)+0.5);
3219	}
3220	Int lstart = max(0, lag0);
3221	Int lend = min(Int(lags.nelements()-1), lag1);
3222	if (lstart == lend) {
3223	lags[lstart] = Complex(0.0);
3224	} else {
3225	if (lstart > lend) {
3226	Int tmp = lend;
3227	lend = lstart;
3228	lstart = tmp;
3229	}
3230	for (int j=lstart; j <=lend ;++j) {
3231	lags[j] = Complex(0.0);
3232	}
3233	}
3234
3235	ffts.fft0(spec, lags);
3236
3237	specCol.put(i, spec);
3238	}
3239	++iter;
3240	}
3241	return out;
3242	}
3243
3244	// Averaging spectra with different channel/resolution
3245	CountedPtr<Scantable>
3246	STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
3247	const bool& compel,
3248	const std::vector<bool>& mask,
3249	const std::string& weight,
3250	const std::string& avmode )
3251	throw ( casa::AipsError )
3252	{
3253	LogIO os( LogOrigin( "STMath", "new_average()", WHERE ) ) ;
3254	if ( avmode == "SCAN" && in.size() != 1 )
3255	throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
3256	"Use merge first."));
3257
3258	CountedPtr<Scantable> out ; // processed result
3259	if ( compel ) {
3260	std::vector< CountedPtr<Scantable> > newin ; // input for average process
3261	uInt insize = in.size() ; // number of input scantables
3262
3263	// setup newin
3264	bool oldInsitu = insitu_ ;
3265	setInsitu( false ) ;
3266	newin.resize( insize ) ;
3267	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3268	newin[itable] = getScantable( in[itable], false ) ;
3269	}
3270	setInsitu( oldInsitu ) ;
3271
3272	// warning
3273	os << "Average spectra with different spectral resolution" << LogIO::POST ;
3274
3275	// temporarily set coordinfo
3276	vector<string> oldinfo( insize ) ;
3277	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3278	vector<string> coordinfo = in[itable]->getCoordInfo() ;
3279	oldinfo[itable] = coordinfo[0] ;
3280	coordinfo[0] = "Hz" ;
3281	newin[itable]->setCoordInfo( coordinfo ) ;
3282	}
3283
3284	ostringstream oss ;
3285
3286	// check IF frequency coverage
3287	// freqid: list of FREQ_ID, which is used, in each table
3288	// iffreq: list of minimum and maximum frequency for each FREQ_ID in
3289	// each table
3290	// freqid[insize][numIF]
3291	// freqid: [[id00, id01, ...],
3292	// [id10, id11, ...],
3293	// ...
3294	// [idn0, idn1, ...]]
3295	// iffreq[insize][numIF*2]
3296	// iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
3297	// [min_id10, max_id10, min_id11, max_id11, ...],
3298	// ...
3299	// [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
3300	//os << "Check IF settings in each table" << LogIO::POST ;
3301	vector< vector<uInt> > freqid( insize );
3302	vector< vector<double> > iffreq( insize ) ;
3303	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3304	Vector<uInt> freqIds = newin[itable]->mfreqidCol_.getColumn() ;
3305	vector<uInt> uniqueFreqId = newin[itable]->getNumbers(newin[itable]->mfreqidCol_) ;
3306	for ( vector<uInt>::iterator i = uniqueFreqId.begin() ;
3307	i != uniqueFreqId.end() ; i++ ) {
3308	//os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
3309	uInt target = 0 ;
3310	while ( freqIds[target] != *i )
3311	target++ ;
3312	vector<double> abcissa = newin[itable]->getAbcissa( target ) ;
3313	freqid[itable].push_back( *i ) ;
3314	double incr = abs( abcissa[1] - abcissa[0] ) ;
3315	iffreq[itable].push_back( (min_element(abcissa.begin(),abcissa.end()))-0.5incr ) ;
3316	iffreq[itable].push_back( (max_element(abcissa.begin(),abcissa.end()))+0.5incr ) ;
3317	}
3318	}
3319
3320	// debug
3321	// os << "IF settings summary:" << endl ;
3322	// for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
3323	// os << " Table" << i << endl ;
3324	// for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
3325	// os << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2j] << "," << iffreq[i][2j+1] << ")" << endl ;
3326	// }
3327	// }
3328	// os << endl ;
3329	// os.post() ;
3330
3331	// IF grouping based on their frequency coverage
3332	// ifgrp: number of member in each IF group
3333	// ifgrp[numgrp]
3334	// ifgrp: [n0, n1, ...]
3335	//os << "IF grouping based on their frequency coverage" << LogIO::POST ;
3336
3337	// parameter for IF grouping
3338	// groupmode = OR retrieve all region
3339	// AND only retrieve overlaped region
3340	//string groupmode = "AND" ;
3341	string groupmode = "OR" ;
3342	uInt sizecr = 0 ;
3343	if ( groupmode == "AND" )
3344	sizecr = 1 ;
3345	else if ( groupmode == "OR" )
3346	sizecr = 0 ;
3347
3348	vector<double> sortedfreq ;
3349	for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
3350	for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
3351	if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
3352	sortedfreq.push_back( iffreq[i][j] ) ;
3353	}
3354	}
3355	sort( sortedfreq.begin(), sortedfreq.end() ) ;
3356	vector<uInt> ifgrp( sortedfreq.size()-1, 0 ) ;
3357	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3358	for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
3359	double range0 = iffreq[itable][2*iif] ;
3360	double range1 = iffreq[itable][2*iif+1] ;
3361	for ( uInt j = 0 ; j < sortedfreq.size()-1 ; j++ ) {
3362	double fmin = max( range0, sortedfreq[j] ) ;
3363	double fmax = min( range1, sortedfreq[j+1] ) ;
3364	if ( fmin < fmax ) {
3365	ifgrp[j]++ ;
3366	}
3367	}
3368	}
3369	}
3370
3371	// Grouping continuous IF groups (without frequency gap)
3372	// freqgrp: list of IF group indexes in each frequency group
3373	// freqgrp[numgrp][nummember]
3374	// freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
3375	// [ifgrp10, ifgrp11, ifgrp12, ...],
3376	// ...
3377	// [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
3378	// grprange[2*numgrp]
3379	// grprange: [fmin0,fmax0,fmin1,fmax1,...]
3380	vector< vector<uInt> > freqgrp ;
3381	vector<double> grprange ;
3382	vector<uInt> grpedge ;
3383	for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
3384	if ( ifgrp[igrp] <= sizecr ) {
3385	grpedge.push_back( igrp ) ;
3386	}
3387	}
3388	grpedge.push_back( ifgrp.size() ) ;
3389	uInt itmp = 0 ;
3390	for ( uInt i = 0 ; i < grpedge.size() ; i++ ) {
3391	int n = grpedge[i] - itmp ;
3392	if ( n > 0 ) {
3393	vector<uInt> members( n ) ;
3394	for ( int j = 0 ; j < n ; j++ ) {
3395	members[j] = itmp+j ;
3396	}
3397	freqgrp.push_back( members ) ;
3398	grprange.push_back( sortedfreq[itmp] ) ;
3399	grprange.push_back( sortedfreq[grpedge[i]] ) ;
3400	}
3401	itmp += n + 1 ;
3402	}
3403
3404	// print frequency group
3405	oss.str("") ;
3406	oss << "Frequency Group summary: " << endl ;
3407	oss << " GROUP_ID: [FREQ_MIN, FREQ_MAX]" << endl ;
3408	for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
3409	oss << " GROUP " << setw( 2 ) << i << ": [" << grprange[2i] << "," << grprange[2i+1] << "]" ;
3410	oss << endl ;
3411	}
3412	oss << endl ;
3413	os << oss.str() << LogIO::POST ;
3414
3415	// groups: list of frequency group index whose frequency range overlaps
3416	// with that of each table and IF
3417	// groups[numtable][numIF]
3418	// groups: [[grpx, grpy,...],
3419	// [grpa, grpb,...],
3420	// ...
3421	// [grpk, grpm,...]]
3422	vector< vector<uInt> > groups( insize ) ;
3423	for ( uInt i = 0 ; i < insize ; i++ ) {
3424	groups[i].resize( freqid[i].size() ) ;
3425	}
3426	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3427	for ( uInt ifreq = 0 ; ifreq < freqid[itable].size() ; ifreq++ ) {
3428	double minf = iffreq[itable][2*ifreq] ;
3429	uInt groupid ;
3430	for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3431	vector<uInt> memberlist = freqgrp[igrp] ;
3432	if ( (minf >= grprange[2igrp]) && (minf <= grprange[2igrp+1]) ) {
3433	groupid = igrp ;
3434	break ;
3435	}
3436	}
3437	groups[itable][ifreq] = groupid ;
3438	}
3439	}
3440
3441
3442	// print membership
3443	oss.str("") ;
3444	for ( uInt i = 0 ; i < insize ; i++ ) {
3445	oss << "Table " << i << endl ;
3446	for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
3447	oss << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ;
3448	oss << setw( 2 ) << groups[i][j] ;
3449	oss << endl ;
3450	}
3451	}
3452	os << oss.str() << LogIO::POST ;
3453
3454	// reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
3455	//os << "All IF number is set to IF group index" << LogIO::POST ;
3456	// reset SCANNO only when avmode != "SCAN"
3457	if ( avmode != "SCAN" ) {
3458	os << "All scan number is set to 0" << LogIO::POST ;
3459	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3460	uInt nrow = newin[itable]->nrow() ;
3461	Vector<uInt> resetScan( nrow, 0 ) ;
3462	newin[itable]->scanCol_.putColumn( resetScan ) ;
3463	}
3464	}
3465
3466	// reset spectra and flagtra: align spectral resolution
3467	//os << "Align spectral resolution" << LogIO::POST ;
3468	// gmaxdnu: the coarsest frequency resolution in the frequency group
3469	// gminfreq: lower frequency edge of the frequency group
3470	// gnchan: number of channels for the frequency group
3471	vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
3472	vector<double> gminfreq( freqgrp.size() ) ;
3473	vector<double> gnchan( freqgrp.size() ) ;
3474	for ( uInt i = 0 ; i < insize ; i++ ) {
3475	vector<uInt> members = groups[i] ;
3476	for ( uInt j = 0 ; j < members.size() ; j++ ) {
3477	uInt groupid = members[j] ;
3478	Double rp,rv,ic ;
3479	newin[i]->frequencies().getEntry( rp, rv, ic, j ) ;
3480	if ( abs(ic) > abs(gmaxdnu[groupid]) )
3481	gmaxdnu[groupid] = ic ;
3482	}
3483	}
3484	for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3485	gminfreq[igrp] = grprange[2*igrp] ;
3486	double maxfreq = grprange[2*igrp+1] ;
3487	gnchan[igrp] = (int)(abs((maxfreq-gminfreq[igrp])/gmaxdnu[igrp])+0.9) ;
3488	}
3489
3490	// regrid spectral data and update frequency info
3491	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3492	Vector<uInt> oldFreqId = newin[itable]->mfreqidCol_.getColumn() ;
3493	Vector<uInt> newFreqId( oldFreqId.nelements() ) ;
3494
3495	// update MAIN
3496	for ( uInt irow = 0 ; irow < newin[itable]->nrow() ; irow++ ) {
3497	uInt groupid = groups[itable][oldFreqId[irow]] ;
3498	newFreqId[irow] = groupid ;
3499	newin[itable]->regridChannel( gnchan[groupid],
3500	gmaxdnu[groupid],
3501	gminfreq[groupid],
3502	irow ) ;
3503	}
3504	newin[itable]->mfreqidCol_.putColumn( newFreqId ) ;
3505	newin[itable]->ifCol_.putColumn( newFreqId ) ;
3506
3507	// update FREQUENCIES
3508	Table tab = newin[itable]->frequencies().table() ;
3509	ScalarColumn<uInt> fIdCol( tab, "ID" ) ;
3510	ScalarColumn<Double> fRefPixCol( tab, "REFPIX" ) ;
3511	ScalarColumn<Double> fRefValCol( tab, "REFVAL" ) ;
3512	ScalarColumn<Double> fIncrCol( tab, "INCREMENT" ) ;
3513	if ( freqgrp.size() > tab.nrow() ) {
3514	tab.addRow( freqgrp.size()-tab.nrow() ) ;
3515	}
3516	for ( uInt irow = 0 ; irow < freqgrp.size() ; irow++ ) {
3517	Double refval = gminfreq[irow] + 0.5 * abs(gmaxdnu[irow]) ;
3518	Double refpix = (gmaxdnu[irow] > 0.0) ? 0 : gnchan[irow]-1 ;
3519	Double increment = gmaxdnu[irow] ;
3520	fIdCol.put( irow, irow ) ;
3521	fRefPixCol.put( irow, refpix ) ;
3522	fRefValCol.put( irow, refval ) ;
3523	fIncrCol.put( irow, increment ) ;
3524	}
3525	}
3526
3527	// set back coordinfo
3528	for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3529	vector<string> coordinfo = newin[itable]->getCoordInfo() ;
3530	coordinfo[0] = oldinfo[itable] ;
3531	newin[itable]->setCoordInfo( coordinfo ) ;
3532	}
3533
3534	// average
3535	out = average( newin, mask, weight, avmode ) ;
3536	}
3537	else {
3538	// simple average
3539	out = average( in, mask, weight, avmode ) ;
3540	}
3541
3542	return out;
3543	}
3544
3545	CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s,
3546	const String calmode,
3547	const String antname )
3548	{
3549	// frequency switch
3550	if ( calmode == "fs" ) {
3551	return cwcalfs( s, antname ) ;
3552	}
3553	else {
3554	vector<bool> masks = s->getMask( 0 ) ;
3555	vector<int> types ;
3556
3557	// save original table selection
3558	Table torg = s->table_ ;
3559
3560	// sky scan
3561	bool insitu = insitu_ ;
3562	insitu_ = false ;
3563	// share calibration scans before average with out
3564	CountedPtr<Scantable> out = getScantable( s, true ) ;
3565	insitu_ = insitu ;
3566	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
3567	out->attach() ;
3568	CountedPtr<Scantable> asky = averageWithinSession( out,
3569	masks,
3570	"TINT" ) ;
3571	// hot scan
3572	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
3573	out->attach() ;
3574	CountedPtr<Scantable> ahot = averageWithinSession( out,
3575	masks,
3576	"TINT" ) ;
3577	// cold scan
3578	CountedPtr<Scantable> acold ;
3579	// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
3580	// out->attach() ;
3581	// CountedPtr<Scantable> acold = averageWithinSession( out,
3582	// masks,
3583	// "TINT" ) ;
3584
3585	// off scan
3586	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3587	out->attach() ;
3588	CountedPtr<Scantable> aoff = averageWithinSession( out,
3589	masks,
3590	"TINT" ) ;
3591
3592	// on scan
3593	s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3594	s->attach() ;
3595	out->table_ = out->originalTable_ ;
3596	out->attach() ;
3597	out->table().addRow( s->nrow() ) ;
3598	copyRows( out->table(), s->table(), 0, 0, s->nrow(), False, True, False ) ;
3599
3600	// iterate throgh STIdxIter2
3601	ChopperWheelCalibrator cal(out, s, aoff, asky, ahot, acold);
3602	STIdxIter2::Iterate<ChopperWheelCalibrator>(cal, "BEAMNO,POLNO,IFNO");
3603
3604	s->table_ = torg ;
3605	s->attach() ;
3606
3607	// flux unit
3608	out->setFluxUnit( "K" ) ;
3609
3610	return out ;
3611	}
3612	}
3613
3614	CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s,
3615	const String calmode )
3616	{
3617	// frequency switch
3618	if ( calmode == "fs" ) {
3619	return almacalfs( s ) ;
3620	}
3621	else {
3622	// double t0, t1 ;
3623	// t0 = mathutil::gettimeofday_sec() ;
3624	vector<bool> masks = s->getMask( 0 ) ;
3625
3626	// save original table selection
3627	Table torg = s->table_ ;
3628
3629	// off scan
3630	// TODO 2010/01/08 TN
3631	// Grouping by time should be needed before averaging.
3632	// Each group must have own unique SCANNO (should be renumbered).
3633	// See PIPELINE/SDCalibration.py
3634	bool insitu = insitu_ ;
3635	insitu_ = false ;
3636	// share off scan before average with out
3637	CountedPtr<Scantable> out = getScantable( s, true ) ;
3638	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3639	out->attach() ;
3640	insitu_ = insitu ;
3641	CountedPtr<Scantable> aoff = averageWithinSession( out,
3642	masks,
3643	"TINT" ) ;
3644
3645	// on scan
3646	// t0 = mathutil::gettimeofday_sec() ;
3647	s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3648	s->attach() ;
3649	out->table_ = out->originalTable_ ;
3650	out->attach() ;
3651	out->table().addRow( s->nrow() ) ;
3652	copyRows( out->table(), s->table(), 0, 0, s->nrow(), False ) ;
3653	// t1 = mathutil::gettimeofday_sec() ;
3654	// cout << "elapsed time for preparing output table: " << t1-t0 << " sec" << endl ;
3655
3656	// process each on scan
3657	// t0 = mathutil::gettimeofday_sec() ;
3658
3659	// iterate throgh STIdxIter2
3660	AlmaCalibrator cal(out, s, aoff);
3661	STIdxIter2::Iterate<AlmaCalibrator>(cal, "BEAMNO,POLNO,IFNO");
3662
3663	// finalize
3664	s->table_ = torg ;
3665	s->attach() ;
3666
3667	// t1 = mathutil::gettimeofday_sec() ;
3668	// cout << "elapsed time for calibration: " << t1-t0 << " sec" << endl ;
3669
3670	// flux unit
3671	out->setFluxUnit( "K" ) ;
3672
3673	return out ;
3674	}
3675	}
3676
3677	CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
3678	const String antname )
3679	{
3680	vector<int> types ;
3681
3682	// APEX calibration mode
3683	int apexcalmode = 1 ;
3684
3685	if ( antname.find( "APEX" ) != string::npos ) {
3686	// check if off scan exists or not
3687	STSelector sel = STSelector() ;
3688	//sel.setName( offstr1 ) ;
3689	types.push_back( SrcType::FLOOFF ) ;
3690	sel.setTypes( types ) ;
3691	try {
3692	s->setSelection( sel ) ;
3693	}
3694	catch ( AipsError &e ) {
3695	apexcalmode = 0 ;
3696	}
3697	sel.reset() ;
3698	}
3699	s->unsetSelection() ;
3700	types.clear() ;
3701
3702	vector<bool> masks = s->getMask( 0 ) ;
3703	CountedPtr<Scantable> ssig, sref ;
3704	//CountedPtr<Scantable> out ;
3705	bool insitu = insitu_ ;
3706	insitu_ = False ;
3707	CountedPtr<Scantable> out = getScantable( s, true ) ;
3708	insitu_ = insitu ;
3709
3710	if ( antname.find( "APEX" ) != string::npos ) {
3711	// APEX calibration
3712	// sky scan
3713	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOSKY ) ;
3714	out->attach() ;
3715	CountedPtr<Scantable> askylo = averageWithinSession( out,
3716	masks,
3717	"TINT" ) ;
3718	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHISKY ) ;
3719	out->attach() ;
3720	CountedPtr<Scantable> askyhi = averageWithinSession( out,
3721	masks,
3722	"TINT" ) ;
3723
3724	// hot scan
3725	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOHOT ) ;
3726	out->attach() ;
3727	CountedPtr<Scantable> ahotlo = averageWithinSession( out,
3728	masks,
3729	"TINT" ) ;
3730	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIHOT ) ;
3731	out->attach() ;
3732	CountedPtr<Scantable> ahothi = averageWithinSession( out,
3733	masks,
3734	"TINT" ) ;
3735
3736	// cold scan
3737	CountedPtr<Scantable> acoldlo, acoldhi ;
3738	// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOCOLD ) ;
3739	// out->attach() ;
3740	// CountedPtr<Scantable> acoldlo = averageWithinSession( out,
3741	// masks,
3742	// "TINT" ) ;
3743	// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHICOLD ) ;
3744	// out->attach() ;
3745	// CountedPtr<Scantable> acoldhi = averageWithinSession( out,
3746	// masks,
3747	// "TINT" ) ;
3748
3749	// ref scan
3750	insitu_ = false ;
3751	sref = getScantable( s, true ) ;
3752	CountedPtr<Scantable> rref = getScantable( s, true ) ;
3753	insitu_ = insitu ;
3754	rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSLO ) ;
3755	rref->attach() ;
3756	copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
3757
3758	// sig scan
3759	insitu_ = false ;
3760	ssig = getScantable( s, true ) ;
3761	CountedPtr<Scantable> rsig = getScantable( s, true ) ;
3762	insitu_ = insitu ;
3763	rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSHI ) ;
3764	rsig->attach() ;
3765	copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
3766
3767	if ( apexcalmode == 0 ) {
3768	// using STIdxIter2
3769	vector<string> cols( 3 ) ;
3770	cols[0] = "BEAMNO" ;
3771	cols[1] = "POLNO" ;
3772	cols[2] = "IFNO" ;
3773	STIdxIter2 iter(ssig, cols) ;
3774	STSelector sel ;
3775	vector< CountedPtr<Scantable> > on( 2 ) ;
3776	on[0] = rsig ;
3777	on[1] = rref ;
3778	vector< CountedPtr<Scantable> > sky( 2 ) ;
3779	sky[0] = askylo ;
3780	sky[1] = askyhi ;
3781	vector< CountedPtr<Scantable> > hot( 2 ) ;
3782	hot[0] = ahotlo ;
3783	hot[1] = ahothi ;
3784	vector< CountedPtr<Scantable> > cold( 2 ) ;
3785	while ( !iter.pastEnd() ) {
3786	Record ids = iter.currentValue() ;
3787	stringstream ss ;
3788	ss << "SELECT FROM $1 WHERE "
3789	<< "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3790	<< "POLNO==" << ids.asuInt("POLNO") << "&&"
3791	<< "IFNO==" << ids.asuInt("IFNO") ;
3792	//cout << "TaQL string: " << ss.str() << endl ;
3793	sel.setTaQL( ss.str() ) ;
3794	sky[0]->setSelection( sel ) ;
3795	sky[1]->setSelection( sel ) ;
3796	hot[0]->setSelection( sel ) ;
3797	hot[1]->setSelection( sel ) ;
3798	Vector<uInt> rows = iter.getRows( SHARE ) ;
3799	calibrateAPEXFS( ssig, sref, on, sky, hot, cold, rows ) ;
3800	sky[0]->unsetSelection() ;
3801	sky[1]->unsetSelection() ;
3802	hot[0]->unsetSelection() ;
3803	hot[1]->unsetSelection() ;
3804	sel.reset() ;
3805	iter.next() ;
3806	}
3807
3808	}
3809	else if ( apexcalmode == 1 ) {
3810	// APEX fs data with off scan
3811	// off scan
3812	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOOFF ) ;
3813	out->attach() ;
3814	CountedPtr<Scantable> aofflo = averageWithinSession( out,
3815	masks,
3816	"TINT" ) ;
3817	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIOFF ) ;
3818	out->attach() ;
3819	CountedPtr<Scantable> aoffhi = averageWithinSession( out,
3820	masks,
3821	"TINT" ) ;
3822
3823	// process each sig and ref scan
3824	// iterate throgh STIdxIter2
3825	ChopperWheelCalibrator cal_sig(ssig, rsig, aofflo, askylo, ahotlo, acoldlo);
3826	STIdxIter2::Iterate<ChopperWheelCalibrator>(cal_sig, "BEAMNO,POLNO,IFNO");
3827	ChopperWheelCalibrator cal_ref(sref, rref, aoffhi, askyhi, ahothi, acoldhi);
3828	STIdxIter2::Iterate<ChopperWheelCalibrator>(cal_ref, "BEAMNO,POLNO,IFNO");
3829	}
3830	}
3831	else {
3832	// non-APEX fs data
3833	// sky scan
3834	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
3835	out->attach() ;
3836	CountedPtr<Scantable> asky = averageWithinSession( out,
3837	masks,
3838	"TINT" ) ;
3839	STSelector sel = STSelector() ;
3840
3841	// hot scan
3842	out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
3843	out->attach() ;
3844	CountedPtr<Scantable> ahot = averageWithinSession( out,
3845	masks,
3846	"TINT" ) ;
3847
3848	// cold scan
3849	CountedPtr<Scantable> acold ;
3850	// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
3851	// out->attach() ;
3852	// CountedPtr<Scantable> acold = averageWithinSession( out,
3853	// masks,
3854	// "TINT" ) ;
3855
3856	// ref scan
3857	bool insitu = insitu_ ;
3858	insitu_ = false ;
3859	sref = getScantable( s, true ) ;
3860	CountedPtr<Scantable> rref = getScantable( s, true ) ;
3861	insitu_ = insitu ;
3862	rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3863	rref->attach() ;
3864	copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
3865
3866	// sig scan
3867	insitu_ = false ;
3868	ssig = getScantable( s, true ) ;
3869	CountedPtr<Scantable> rsig = getScantable( s, true ) ;
3870	insitu_ = insitu ;
3871	rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3872	rsig->attach() ;
3873	copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
3874
3875	// process each sig and ref scan
3876	vector<string> cols( 3 ) ;
3877	cols[0] = "BEAMNO" ;
3878	cols[1] = "POLNO" ;
3879	cols[2] = "IFNO" ;
3880	STIdxIter2 iter(ssig, cols);
3881	while ( !iter.pastEnd() ) {
3882	Record ids = iter.currentValue() ;
3883	stringstream ss ;
3884	ss << "SELECT FROM $1 WHERE "
3885	<< "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3886	<< "POLNO==" << ids.asuInt("POLNO") << "&&"
3887	<< "IFNO==" << ids.asuInt("IFNO") ;
3888	//cout << "TaQL string: " << ss.str() << endl ;
3889	sel.setTaQL( ss.str() ) ;
3890	ahot->setSelection( sel ) ;
3891	asky->setSelection( sel ) ;
3892	Vector<uInt> rows = iter.getRows( SHARE ) ;
3893	// out should be an exact copy of s except that SPECTRA column is empty
3894	calibrateFS( ssig, sref, rsig, rref, asky, ahot, acold, rows ) ;
3895	ahot->unsetSelection() ;
3896	asky->unsetSelection() ;
3897	sel.reset() ;
3898	iter.next() ;
3899	}
3900	}
3901
3902	// do folding if necessary
3903	Table sigtab = ssig->table() ;
3904	Table reftab = sref->table() ;
3905	ScalarColumn<uInt> reffidCol ;
3906	Int nchan = (Int)ssig->nchan() ;
3907	reffidCol.attach( reftab, "FREQ_ID" ) ;
3908	Vector<uInt> sfids = ssig->mfreqidCol_.getColumn() ;
3909	Vector<uInt> rfids = sref->mfreqidCol_.getColumn() ;
3910	vector<uInt> sfids_unique ;
3911	vector<uInt> rfids_unique ;
3912	vector<uInt> sifno_unique ;
3913	vector<uInt> rifno_unique ;
3914	for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) {
3915	if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) {
3916	sfids_unique.push_back( sfids[i] ) ;
3917	sifno_unique.push_back( ssig->getIF( i ) ) ;
3918	}
3919	if ( count( rfids_unique.begin(), rfids_unique.end(), rfids[i] ) == 0 ) {
3920	rfids_unique.push_back( rfids[i] ) ;
3921	rifno_unique.push_back( sref->getIF( i ) ) ;
3922	}
3923	}
3924	double refpix_sig, refval_sig, increment_sig ;
3925	double refpix_ref, refval_ref, increment_ref ;
3926	vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ;
3927	for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) {
3928	ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ;
3929	sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ;
3930	if ( refpix_sig == refpix_ref ) {
3931	double foffset = refval_ref - refval_sig ;
3932	int choffset = static_cast<int>(foffset/increment_sig) ;
3933	double doffset = foffset / increment_sig ;
3934	if ( abs(choffset) >= nchan ) {
3935	LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ;
3936	os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ;
3937	os << "Just return signal data" << LogIO::POST ;
3938	//std::vector< CountedPtr<Scantable> > tabs ;
3939	//tabs.push_back( ssig ) ;
3940	//tabs.push_back( sref ) ;
3941	//out = merge( tabs ) ;
3942	tmp[i] = ssig ;
3943	}
3944	else {
3945	STSelector sel = STSelector() ;
3946	vector<int> v( 1, sifno_unique[i] ) ;
3947	sel.setIFs( v ) ;
3948	ssig->setSelection( sel ) ;
3949	sel.reset() ;
3950	v[0] = rifno_unique[i] ;
3951	sel.setIFs( v ) ;
3952	sref->setSelection( sel ) ;
3953	sel.reset() ;
3954	if ( antname.find( "APEX" ) != string::npos ) {
3955	tmp[i] = dofold( ssig, sref, 0.5doffset, -0.5doffset ) ;
3956	//tmp[i] = dofold( ssig, sref, doffset ) ;
3957	}
3958	else {
3959	tmp[i] = dofold( ssig, sref, doffset ) ;
3960	}
3961	ssig->unsetSelection() ;
3962	sref->unsetSelection() ;
3963	}
3964	}
3965	}
3966
3967	if ( tmp.size() > 1 ) {
3968	out = merge( tmp ) ;
3969	}
3970	else {
3971	out = tmp[0] ;
3972	}
3973
3974	// flux unit
3975	out->setFluxUnit( "K" ) ;
3976
3977	return out ;
3978	}
3979
3980	CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s )
3981	{
3982	(void) s; //currently unused
3983	CountedPtr<Scantable> out ;
3984
3985	return out ;
3986	}
3987
3988	void STMath::calibrateAPEXFS( CountedPtr<Scantable> &sig,
3989	CountedPtr<Scantable> &ref,
3990	const vector< CountedPtr<Scantable> >& on,
3991	const vector< CountedPtr<Scantable> >& sky,
3992	const vector< CountedPtr<Scantable> >& hot,
3993	const vector< CountedPtr<Scantable> >& cold,
3994	const Vector<uInt> &rows )
3995	{
3996	// if rows is empty, just return
3997	if ( rows.nelements() == 0 )
3998	return ;
3999	ROScalarColumn<Double> timeCol( sky[0]->table(), "TIME" ) ;
4000	Vector<Double> timeSkyS = timeCol.getColumn() ;
4001	timeCol.attach( sky[1]->table(), "TIME" ) ;
4002	Vector<Double> timeSkyR = timeCol.getColumn() ;
4003	timeCol.attach( hot[0]->table(), "TIME" ) ;
4004	Vector<Double> timeHotS = timeCol.getColumn() ;
4005	timeCol.attach( hot[1]->table(), "TIME" ) ;
4006	Vector<Double> timeHotR = timeCol.getColumn() ;
4007	timeCol.attach( sig->table(), "TIME" ) ;
4008	ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
4009	ROArrayColumn<Float> arrayFloatCol( sky[0]->table(), "SPECTRA" ) ;
4010	SpectralData skyspectraS(arrayFloatCol.getColumn());
4011	arrayFloatCol.attach( sky[1]->table(), "SPECTRA" ) ;
4012	SpectralData skyspectraR(arrayFloatCol.getColumn());
4013	arrayFloatCol.attach( hot[0]->table(), "SPECTRA" ) ;
4014	SpectralData hotspectraS(arrayFloatCol.getColumn());
4015	arrayFloatCol.attach( hot[1]->table(), "SPECTRA" ) ;
4016	SpectralData hotspectraR(arrayFloatCol.getColumn());
4017	TcalData tcaldataS(sky[0]);
4018	TsysData tsysdataS(sky[0]);
4019	TcalData tcaldataR(sky[1]);
4020	TsysData tsysdataR(sky[1]);
4021	unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
4022	Vector<Float> spec( spsize ) ;
4023	// I know that the data is contiguous
4024	const uInt *p = rows.data() ;
4025	vector<int> ids( 2 ) ;
4026	Block<uInt> flagchan( spsize ) ;
4027	uInt nflag = 0 ;
4028	for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4029	double reftime = timeCol.asdouble(*p) ;
4030	ids = getRowIdFromTime( reftime, timeSkyS ) ;
4031	Vector<Float> spskyS = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSkyS, ids, skyspectraS, "linear");
4032	Vector<Float> tcalS = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSkyS, ids, tcaldataS, "linear");
4033	Vector<Float> tsysS = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSkyS, ids, tsysdataS, "linear");
4034	ids = getRowIdFromTime( reftime, timeHotS ) ;
4035	Vector<Float> sphotS = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHotS, ids, hotspectraS, "linear");
4036	reftime = timeCol2.asdouble(*p) ;
4037	ids = getRowIdFromTime( reftime, timeSkyR ) ;
4038	Vector<Float> spskyR = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSkyR, ids, skyspectraR, "linear");
4039	Vector<Float> tcalR = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSkyR, ids, tcaldataR, "linear");
4040	Vector<Float> tsysR = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSkyR, ids, tsysdataR, "linear");
4041	ids = getRowIdFromTime( reftime, timeHotR ) ;
4042	Vector<Float> sphotR = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHotR, ids, hotspectraR, "linear");
4043	Vector<Float> spsig = on[0]->specCol_( *p ) ;
4044	Vector<Float> spref = on[1]->specCol_( *p ) ;
4045	for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4046	if ( (sphotS[j]-spskyS[j]) == 0.0 \|\| (sphotR[j]-spskyR[j]) == 0.0 ) {
4047	spec[j] = 0.0 ;
4048	flagchan[nflag++] = j ;
4049	}
4050	else {
4051	spec[j] = tcalS[j] * spsig[j] / ( sphotS[j] - spskyS[j] )
4052	- tcalR[j] * spref[j] / ( sphotR[j] - spskyR[j] ) ;
4053	}
4054	}
4055	sig->specCol_.put( *p, spec ) ;
4056	sig->tsysCol_.put( *p, tsysS ) ;
4057	spec *= (Float)-1.0 ;
4058	ref->specCol_.put( *p, spec ) ;
4059	ref->tsysCol_.put( *p, tsysR ) ;
4060	if ( nflag > 0 ) {
4061	Vector<uChar> flsig = sig->flagsCol_( *p ) ;
4062	Vector<uChar> flref = ref->flagsCol_( *p ) ;
4063	for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4064	flsig[flagchan[j]] = (uChar)True ;
4065	flref[flagchan[j]] = (uChar)True ;
4066	}
4067	sig->flagsCol_.put( *p, flsig ) ;
4068	ref->flagsCol_.put( *p, flref ) ;
4069	}
4070	nflag = 0 ;
4071	p++ ;
4072	}
4073	}
4074
4075	void STMath::calibrateFS( CountedPtr<Scantable> &sig,
4076	CountedPtr<Scantable> &ref,
4077	const CountedPtr<Scantable>& rsig,
4078	const CountedPtr<Scantable>& rref,
4079	const CountedPtr<Scantable>& sky,
4080	const CountedPtr<Scantable>& hot,
4081	const CountedPtr<Scantable>& cold,
4082	const Vector<uInt> &rows )
4083	{
4084	// if rows is empty, just return
4085	if ( rows.nelements() == 0 )
4086	return ;
4087	ROScalarColumn<Double> timeCol( sky->table(), "TIME" ) ;
4088	Vector<Double> timeSky = timeCol.getColumn() ;
4089	timeCol.attach( hot->table(), "TIME" ) ;
4090	Vector<Double> timeHot = timeCol.getColumn() ;
4091	timeCol.attach( sig->table(), "TIME" ) ;
4092	ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
4093	ROArrayColumn<Float> arrayFloatCol( sky->table(), "SPECTRA" ) ;
4094	SpectralData skyspectra(arrayFloatCol.getColumn());
4095	arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
4096	SpectralData hotspectra(arrayFloatCol.getColumn());
4097	TcalData tcaldata(sky);
4098	TsysData tsysdata(sky);
4099	unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
4100	Vector<Float> spec( spsize ) ;
4101	// I know that the data is contiguous
4102	const uInt *p = rows.data() ;
4103	vector<int> ids( 2 ) ;
4104	Block<uInt> flagchan( spsize ) ;
4105	uInt nflag = 0 ;
4106	for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4107	double reftime = timeCol.asdouble(*p) ;
4108	ids = getRowIdFromTime( reftime, timeSky ) ;
4109	Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
4110	Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
4111	Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
4112	ids = getRowIdFromTime( reftime, timeHot ) ;
4113	Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
4114	Vector<Float> spsig = rsig->specCol_( *p ) ;
4115	Vector<Float> spref = rref->specCol_( *p ) ;
4116	// using gain array
4117	for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4118	if ( spref[j] == 0.0 \|\| (sphot[j]-spsky[j]) == 0.0 ) {
4119	spec[j] = 0.0 ;
4120	flagchan[nflag++] = j ;
4121	}
4122	else {
4123	spec[j] = ( ( spsig[j] - spref[j] ) / spref[j] )
4124	* ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4125	}
4126	}
4127	sig->specCol_.put( *p, spec ) ;
4128	sig->tsysCol_.put( *p, tsys ) ;
4129	if ( nflag > 0 ) {
4130	Vector<uChar> fl = sig->flagsCol_( *p ) ;
4131	for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4132	fl[flagchan[j]] = (uChar)True ;
4133	}
4134	sig->flagsCol_.put( *p, fl ) ;
4135	}
4136	nflag = 0 ;
4137
4138	reftime = timeCol2.asdouble(*p) ;
4139	ids = getRowIdFromTime( reftime, timeSky ) ;
4140	spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
4141	tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
4142	tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
4143	ids = getRowIdFromTime( reftime, timeHot ) ;
4144	sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
4145	// using gain array
4146	for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4147	if ( spsig[j] == 0.0 \|\| (sphot[j]-spsky[j]) == 0.0 ) {
4148	spec[j] = 0.0 ;
4149	flagchan[nflag++] = j ;
4150	}
4151	else {
4152	spec[j] = ( ( spref[j] - spsig[j] ) / spsig[j] )
4153	* ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4154	}
4155	}
4156	ref->specCol_.put( *p, spec ) ;
4157	ref->tsysCol_.put( *p, tsys ) ;
4158	if ( nflag > 0 ) {
4159	Vector<uChar> fl = ref->flagsCol_( *p ) ;
4160	for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4161	fl[flagchan[j]] = (uChar)True ;
4162	}
4163	ref->flagsCol_.put( *p, fl ) ;
4164	}
4165	nflag = 0 ;
4166	p++ ;
4167	}
4168	}
4169
4170	void STMath::copyRows( Table &out,
4171	const Table &in,
4172	uInt startout,
4173	uInt startin,
4174	uInt nrow,
4175	Bool copySpectra,
4176	Bool copyFlagtra,
4177	Bool copyTsys )
4178	{
4179	uInt nexclude = 0 ;
4180	Block<String> excludeColsBlock( 3 ) ;
4181	if ( !copySpectra ) {
4182	excludeColsBlock[nexclude] = "SPECTRA" ;
4183	nexclude++ ;
4184	}
4185	if ( !copyFlagtra ) {
4186	excludeColsBlock[nexclude] = "FLAGTRA" ;
4187	nexclude++ ;
4188	}
4189	if ( !copyTsys ) {
4190	excludeColsBlock[nexclude] = "TSYS" ;
4191	nexclude++ ;
4192	}
4193	// if ( nexclude < 3 ) {
4194	// excludeCols.resize( nexclude, True ) ;
4195	// }
4196	Vector<String> excludeCols( IPosition(1,nexclude),
4197	excludeColsBlock.storage(),
4198	SHARE ) ;
4199	// cout << "excludeCols=" << excludeCols << endl ;
4200	TableRow rowout( out, excludeCols, True ) ;
4201	ROTableRow rowin( in, excludeCols, True ) ;
4202	uInt rin = startin ;
4203	uInt rout = startout ;
4204	for ( uInt i = 0 ; i < nrow ; i++ ) {
4205	rowin.get( rin ) ;
4206	rowout.putMatchingFields( rout, rowin.record() ) ;
4207	rin++ ;
4208	rout++ ;
4209	}
4210	}
4211
4212	CountedPtr<Scantable> STMath::averageWithinSession( CountedPtr<Scantable> &s,
4213	vector<bool> &mask,
4214	string weight )
4215	{
4216	// prepare output table
4217	bool insitu = insitu_ ;
4218	insitu_ = false ;
4219	CountedPtr<Scantable> a = getScantable( s, true ) ;
4220	insitu_ = insitu ;
4221	Table &atab = a->table() ;
4222	ScalarColumn<Double> timeColOut( atab, "TIME" ) ;
4223
4224	if ( s->nrow() == 0 )
4225	return a ;
4226
4227	// setup RowAccumulator
4228	WeightType wtype = stringToWeight( weight ) ;
4229	RowAccumulator acc( wtype ) ;
4230	Vector<Bool> cmask( mask ) ;
4231	acc.setUserMask( cmask ) ;
4232
4233	vector<string> cols( 3 ) ;
4234	cols[0] = "IFNO" ;
4235	cols[1] = "POLNO" ;
4236	cols[2] = "BEAMNO" ;
4237	STIdxIter2 iter( s, cols ) ;
4238
4239	Table ttab = s->table() ;
4240	ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ;
4241	Vector<Double> timeVec = timeCol->getColumn() ;
4242	delete timeCol ;
4243	Vector<Double> interval = s->integrCol_.getColumn() ;
4244	uInt nrow = timeVec.nelements() ;
4245	uInt outrow = 0 ;
4246
4247	while( !iter.pastEnd() ) {
4248
4249	Vector<uInt> rows = iter.getRows( SHARE ) ;
4250	uInt len = rows.nelements() ;
4251
4252	if ( len == 0 ) {
4253	iter.next() ;
4254	continue ;
4255	}
4256
4257	uInt nchan = s->nchan(s->getIF(rows[0])) ;
4258	Vector<uChar> flag( nchan ) ;
4259	Vector<Bool> bflag( nchan ) ;
4260	Vector<Float> spec( nchan ) ;
4261	Vector<Float> tsys( nchan ) ;
4262
4263	Vector<Double> timeSep( len-1 ) ;
4264	for ( uInt i = 0 ; i < len-1 ; i++ ) {
4265	timeSep[i] = timeVec[rows[i+1]] - timeVec[rows[i]] ;
4266	}
4267
4268	uInt irow ;
4269	uInt jrow ;
4270	for ( uInt i = 0 ; i < len-1 ; i++ ) {
4271	irow = rows[i] ;
4272	jrow = rows[i+1] ;
4273	// accumulate data
4274	s->flagsCol_.get( irow, flag ) ;
4275	//if row-flagged, all channels set flagged
4276	if (s->getFlagRow(irow)) {
4277	for (uInt k = 0; k < nchan; ++k) {
4278	flag(k) = 1 << 7;
4279	}
4280	}
4281	convertArray( bflag, flag ) ;
4282	s->specCol_.get( irow, spec ) ;
4283	tsys.assign( s->tsysCol_( irow ) ) ;
4284	//if ( !allEQ(bflag,True) )
4285	acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
4286	double gap = 2.0 * 86400.0 * timeSep[i] / ( interval[jrow] + interval[irow] ) ;
4287	//cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
4288	if ( gap > 1.1 ) {
4289	//cout << "detected gap between " << i << " and " << i+1 << endl ;
4290	// put data to output table
4291	// reset RowAccumulator
4292	if ( acc.state() ) {
4293	atab.addRow() ;
4294	copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
4295	acc.replaceNaN() ;
4296	const Vector<Bool> &msk = acc.getMask() ;
4297	convertArray( flag, !msk ) ;
4298	for (uInt k = 0; k < nchan; ++k) {
4299	uChar userFlag = 1 << 7;
4300	if (msk[k]==True) userFlag = 0 << 7;
4301	flag(k) = userFlag;
4302	}
4303	a->flagsCol_.put( outrow, flag ) ;
4304	a->specCol_.put( outrow, acc.getSpectrum() ) ;
4305	a->tsysCol_.put( outrow, acc.getTsys() ) ;
4306	a->integrCol_.put( outrow, acc.getInterval() ) ;
4307	timeColOut.put( outrow, acc.getTime() ) ;
4308	a->cycleCol_.put( outrow, 0 ) ;
4309	}
4310	acc.reset() ;
4311	outrow++ ;
4312	}
4313	}
4314
4315	// accumulate and add last data
4316	irow = rows[len-1] ;
4317	s->flagsCol_.get( irow, flag ) ;
4318	//if row-flagged, all channels set flagged
4319	if (s->getFlagRow(irow)) {
4320	for (uInt k = 0; k < nchan; ++k) {
4321	flag(k) = 1 << 7;
4322	}
4323	}
4324	convertArray( bflag, flag ) ;
4325	s->specCol_.get( irow, spec ) ;
4326	tsys.assign( s->tsysCol_( irow ) ) ;
4327	//if (!allEQ(bflag,True) )
4328	acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
4329	if ( acc.state() ) {
4330	atab.addRow() ;
4331	copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
4332	acc.replaceNaN() ;
4333	const Vector<Bool> &msk = acc.getMask() ;
4334	convertArray( flag, !msk ) ;
4335	for (uInt k = 0; k < nchan; ++k) {
4336	uChar userFlag = 1 << 7;
4337	if (msk[k]==True) userFlag = 0 << 7;
4338	flag(k) = userFlag;
4339	}
4340	a->flagsCol_.put( outrow, flag ) ;
4341	a->specCol_.put( outrow, acc.getSpectrum() ) ;
4342	a->tsysCol_.put( outrow, acc.getTsys() ) ;
4343	a->integrCol_.put( outrow, acc.getInterval() ) ;
4344	timeColOut.put( outrow, acc.getTime() ) ;
4345	a->cycleCol_.put( outrow, 0 ) ;
4346	}
4347	acc.reset() ;
4348	outrow++ ;
4349
4350	iter.next() ;
4351	}
4352
4353	return a ;
4354	}

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

Download in other formats: