Context Navigation

STMath.cpp @ 2443

Visit:

Last change on this file since 2443 was 2443, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: No

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs:

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Test tuning ALMA calibration code.

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

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

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

Context Navigation

source: branches/hpc33/src/STMath.cpp @ 2443

Download in other formats: