Context Navigation

STMath.cpp @ 2950

Visit:

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

New Development: No

JIRA Issue: Yes CAS-6582

Ready for Test: Yes

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Smoothing functions are updated so that flagged rows are not processed.
Flag value for regridding function is changed from 1 to 128.

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

File size: 146.6 KB

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

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

Context Navigation

source: trunk/src/STMath.cpp @ 2950

Download in other formats: