Context Navigation

STMath.cpp@ 1711

Visit:

Last change on this file since 1711 was 1705, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-1799

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Fixed a bug that OSF data were not recognized as OTF scan.

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

File size: 162.6 KB

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

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

Context Navigation

source: branches/alma/src/STMath.cpp@ 1711

Download in other formats: