Context Navigation

STGrid.cpp@ 3037

Visit:

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

New Development: No

JIRA Issue: No

Ready for Test: 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...

Flag out grids without data.

File size: 62.8 KB

Line
1	//
2	// C++ Implementation: STGrid
3	//
4	// Description:
5	//
6	//
7	// Author: Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2011
8	//
9	// Copyright: See COPYING file that comes with this distribution
10	//
11	//
12	#include <casa/BasicSL/String.h>
13	#include <casa/Arrays/Vector.h>
14	#include <casa/Arrays/ArrayMath.h>
15	#include <casa/Quanta/Quantum.h>
16	#include <casa/Quanta/QuantumHolder.h>
17	#include <casa/Utilities/CountedPtr.h>
18	#include <casa/Logging/LogIO.h>
19
20	#include <coordinates/Coordinates/DirectionCoordinate.h>
21
22	#include <tables/Tables/Table.h>
23	#include <tables/Tables/TableRecord.h>
24	#include <tables/Tables/TableRow.h>
25	#include <tables/Tables/ExprNode.h>
26	#include <tables/Tables/ScalarColumn.h>
27	#include <tables/Tables/ArrayColumn.h>
28	#include <tables/Tables/TableCopy.h>
29
30	#include <measures/Measures/MDirection.h>
31
32	#include "MathUtils.h"
33	#include <atnf/PKSIO/SrcType.h>
34
35	#include "STGrid.h"
36
37	using namespace std ;
38	using namespace concurrent ;
39	using namespace casa ;
40	using namespace asap ;
41
42	namespace asap {
43
44	// for performance check
45	double eToInt = 0.0 ;
46	double eGetWeight = 0.0 ;
47
48	// constructor
49	STGrid::STGrid()
50	: vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
51	{
52	init() ;
53	}
54
55	STGrid::STGrid( const string infile )
56	: vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
57	{
58	init() ;
59
60	setFileIn( infile ) ;
61	}
62
63	STGrid::STGrid( const vector<string> infile )
64	{
65	init() ;
66
67	setFileList( infile ) ;
68	}
69
70	void STGrid::init()
71	{
72	ifno_ = -1 ;
73	nx_ = -1 ;
74	ny_ = -1 ;
75	npol_ = 0 ;
76	nchan_ = 0 ;
77	nrow_ = 0 ;
78	cellx_ = 0.0 ;
79	celly_ = 0.0 ;
80	center_ = Vector<Double> ( 2, 0.0 ) ;
81	convType_ = "BOX" ;
82	wtype_ = "UNIFORM" ;
83	convSupport_ = -1 ;
84	userSupport_ = -1 ;
85	truncate_ = "";
86	gwidth_ = "";
87	jwidth_ = "";
88	convSampling_ = 100 ;
89	nprocessed_ = 0 ;
90	nchunk_ = 0 ;
91
92	// initialize user input
93	nxUI_ = -1 ;
94	nyUI_ = -1 ;
95	cellxUI_ = "" ;
96	cellyUI_ = "" ;
97	centerUI_ = "" ;
98	doclip_ = False ;
99	}
100
101	void STGrid::setFileIn( const string infile )
102	{
103	nfile_ = 1 ;
104	String name( infile ) ;
105	infileList_.resize( nfile_ ) ;
106	infileList_[0] = String(infile) ;
107	}
108
109	void STGrid::setFileList( const vector<string> infile )
110	{
111	nfile_ = infile.size() ;
112	infileList_.resize( nfile_ ) ;
113	for ( uInt i = 0 ; i < nfile_ ; i++ ) {
114	infileList_[i] = infile[i] ;
115	}
116	}
117
118	void STGrid::setPolList( vector<unsigned int> pols )
119	{
120	pollist_.assign( Vector<uInt>( pols ) ) ;
121	}
122
123	void STGrid::setScanList( vector<unsigned int> scans )
124	{
125	scanlist_.assign( Vector<uInt>( scans ) ) ;
126	}
127
128	void STGrid::setWeight( const string wType )
129	{
130	wtype_ = String( wType ) ;
131	wtype_.upcase() ;
132	}
133
134	void STGrid::defineImage( int nx,
135	int ny,
136	string scellx,
137	string scelly,
138	string scenter )
139	{
140	nxUI_ = (Int)nx ;
141	nyUI_ = (Int)ny ;
142	cellxUI_ = String( scellx ) ;
143	cellyUI_ = String( scelly ) ;
144	centerUI_ = String( scenter ) ;
145	}
146
147	void STGrid::setFunc( string convType,
148	int convSupport,
149	string truncate,
150	string gwidth,
151	string jwidth )
152	{
153	convType_ = String( convType ) ;
154	convType_.upcase() ;
155	userSupport_ = (Int)convSupport ;
156	truncate_ = String( truncate );
157	gwidth_ = String( gwidth );
158	jwidth_ = String( jwidth );
159	}
160
161	#define NEED_UNDERSCORES
162	#if defined(NEED_UNDERSCORES)
163	#define ggridsd ggridsd_
164	#endif
165	extern "C" {
166	void ggridsd(Double*,
167	const Complex*,
168	Int*,
169	Int*,
170	Int*,
171	const Int*,
172	const Int*,
173	const Float*,
174	Int*,
175	Int*,
176	Complex*,
177	Float*,
178	Int*,
179	Int*,
180	Int *,
181	Int *,
182	Int*,
183	Int*,
184	Float*,
185	Int*,
186	Int*,
187	Double*);
188	}
189	void STGrid::call_ggridsd( Array<Double> &xypos,
190	Array<Complex> &spectra,
191	Int &nvispol,
192	Int &nvischan,
193	Array<Int> &flagtra,
194	Array<Int> &flagrow,
195	Array<Float> &weight,
196	Int &nrow,
197	Int &irow,
198	Array<Complex> &gdata,
199	Array<Float> &gwgt,
200	Int &nx,
201	Int &ny,
202	Int &npol,
203	Int &nchan,
204	Int &support,
205	Int &sampling,
206	Vector<Float> &convFunc,
207	Int *chanMap,
208	Int *polMap )
209	{
210	// parameters for gridding
211	Int idopsf = 0 ;
212	Int len = npol*nchan ;
213	Double *sumw_p = new Double[len] ;
214	{
215	Double *work_p = sumw_p ;
216	for ( Int i = 0 ; i < len ; i++ ) {
217	*work_p = 0.0 ;
218	work_p++ ;
219	}
220	}
221
222	// prepare pointer
223	Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
224	Double *xy_p = xypos.getStorage( deletePos ) ;
225	const Complex *values_p = spectra.getStorage( deleteData ) ;
226	const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
227	const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
228	const Float *wgt_p = weight.getStorage( deleteWgt ) ;
229	Complex *grid_p = gdata.getStorage( deleteDataG ) ;
230	Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
231	Float *conv_p = convFunc.getStorage( deleteConv ) ;
232
233	// pass copy of irow to ggridsd since it will be modified in theroutine
234	Int irowCopy = irow ;
235
236	// call ggridsd
237	ggridsd( xy_p,
238	values_p,
239	&nvispol,
240	&nvischan,
241	&idopsf,
242	flag_p,
243	rflag_p,
244	wgt_p,
245	&nrow,
246	&irowCopy,
247	grid_p,
248	wgrid_p,
249	&nx,
250	&ny,
251	&npol,
252	&nchan,
253	&support,
254	&sampling,
255	conv_p,
256	chanMap,
257	polMap,
258	sumw_p ) ;
259
260	// finalization
261	xypos.putStorage( xy_p, deletePos ) ;
262	spectra.freeStorage( values_p, deleteData ) ;
263	flagtra.freeStorage( flag_p, deleteFlag ) ;
264	flagrow.freeStorage( rflag_p, deleteFlagR ) ;
265	weight.freeStorage( wgt_p, deleteWgt ) ;
266	gdata.putStorage( grid_p, deleteDataG ) ;
267	gwgt.putStorage( wgrid_p, deleteWgtG ) ;
268	convFunc.putStorage( conv_p, deleteConv ) ;
269	delete sumw_p ;
270	}
271
272	#define NEED_UNDERSCORES
273	#if defined(NEED_UNDERSCORES)
274	#define ggridsd2 ggridsd2_
275	#endif
276	extern "C" {
277	void ggridsd2(Double*,
278	const Complex*,
279	Int*,
280	Int*,
281	Int*,
282	const Int*,
283	const Int*,
284	const Float*,
285	Int*,
286	Int*,
287	Complex*,
288	Float*,
289	Int*,
290	Complex*,
291	Float*,
292	Float*,
293	Complex*,
294	Float*,
295	Float*,
296	Int*,
297	Int*,
298	Int *,
299	Int *,
300	Int*,
301	Int*,
302	Float*,
303	Int*,
304	Int*,
305	Double*);
306	}
307	void STGrid::call_ggridsd2( Array<Double> &xypos,
308	Array<Complex> &spectra,
309	Int &nvispol,
310	Int &nvischan,
311	Array<Int> &flagtra,
312	Array<Int> &flagrow,
313	Array<Float> &weight,
314	Int &nrow,
315	Int &irow,
316	Array<Complex> &gdata,
317	Array<Float> &gwgt,
318	Array<Int> &npoints,
319	Array<Complex> &clipmin,
320	Array<Float> &clipwmin,
321	Array<Float> &clipcmin,
322	Array<Complex> &clipmax,
323	Array<Float> &clipwmax,
324	Array<Float> &clipcmax,
325	Int &nx,
326	Int &ny,
327	Int &npol,
328	Int &nchan,
329	Int &support,
330	Int &sampling,
331	Vector<Float> &convFunc,
332	Int *chanMap,
333	Int *polMap )
334	{
335	// parameters for gridding
336	Int idopsf = 0 ;
337	Int len = npol*nchan ;
338	Double *sumw_p = new Double[len] ;
339	{
340	Double *work_p = sumw_p ;
341	for ( Int i = 0 ; i < len ; i++ ) {
342	*work_p = 0.0 ;
343	work_p++ ;
344	}
345	}
346
347	// prepare pointer
348	Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG, deleteNpts, deleteCMin, deleteCWMin, deleteCCMin, deleteCMax, deleteCWMax, deleteCCMax ;
349	Double *xy_p = xypos.getStorage( deletePos ) ;
350	const Complex *values_p = spectra.getStorage( deleteData ) ;
351	const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
352	const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
353	const Float *wgt_p = weight.getStorage( deleteWgt ) ;
354	Complex *grid_p = gdata.getStorage( deleteDataG ) ;
355	Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
356	Float *conv_p = convFunc.getStorage( deleteConv ) ;
357	Int *npts_p = npoints.getStorage( deleteNpts ) ;
358	Complex *cmin_p = clipmin.getStorage( deleteCMin ) ;
359	Float *cwmin_p = clipwmin.getStorage( deleteCWMin ) ;
360	Float *ccmin_p = clipcmin.getStorage( deleteCCMin ) ;
361	Complex *cmax_p = clipmax.getStorage( deleteCMax ) ;
362	Float *cwmax_p = clipwmax.getStorage( deleteCWMax ) ;
363	Float *ccmax_p = clipcmax.getStorage( deleteCCMax ) ;
364
365	// pass copy of irow to ggridsd since it will be modified in theroutine
366	Int irowCopy = irow ;
367
368	// call ggridsd
369	ggridsd2( xy_p,
370	values_p,
371	&nvispol,
372	&nvischan,
373	&idopsf,
374	flag_p,
375	rflag_p,
376	wgt_p,
377	&nrow,
378	&irowCopy,
379	grid_p,
380	wgrid_p,
381	npts_p,
382	cmin_p,
383	cwmin_p,
384	ccmin_p,
385	cmax_p,
386	cwmax_p,
387	ccmax_p,
388	&nx,
389	&ny,
390	&npol,
391	&nchan,
392	&support,
393	&sampling,
394	conv_p,
395	chanMap,
396	polMap,
397	sumw_p ) ;
398
399	// finalization
400	xypos.putStorage( xy_p, deletePos ) ;
401	spectra.freeStorage( values_p, deleteData ) ;
402	flagtra.freeStorage( flag_p, deleteFlag ) ;
403	flagrow.freeStorage( rflag_p, deleteFlagR ) ;
404	weight.freeStorage( wgt_p, deleteWgt ) ;
405	gdata.putStorage( grid_p, deleteDataG ) ;
406	gwgt.putStorage( wgrid_p, deleteWgtG ) ;
407	convFunc.putStorage( conv_p, deleteConv ) ;
408	clipmin.putStorage( cmin_p, deleteCMin ) ;
409	clipwmin.putStorage( cwmin_p, deleteCWMin ) ;
410	clipcmin.putStorage( ccmin_p, deleteCCMin ) ;
411	clipmax.putStorage( cmax_p, deleteCMax ) ;
412	clipwmax.putStorage( cwmax_p, deleteCWMax ) ;
413	clipcmax.putStorage( ccmax_p, deleteCCMax ) ;
414	delete sumw_p ;
415	}
416
417	void STGrid::grid()
418	{
419	LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
420	double t0,t1 ;
421
422	// data selection
423	t0 = mathutil::gettimeofday_sec() ;
424	selectData() ;
425	t1 = mathutil::gettimeofday_sec() ;
426	os << LogIO::DEBUGGING << "selectData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
427
428	setupGrid() ;
429	setupArray() ;
430
431	if ( wtype_.compare("UNIFORM") != 0 &&
432	wtype_.compare("TINT") != 0 &&
433	wtype_.compare("TSYS") != 0 &&
434	wtype_.compare("TINTSYS") != 0 ) {
435	LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
436	os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
437	wtype_ = "UNIFORM" ;
438	}
439
440	// Warn if gauss or gjinc gridding with non-square cells
441	if ((cellx_ != celly_) && (convType_=="GAUSS"\|\|convType_=="GJINC")) {
442	os << LogIO::WARN
443	<< "The " << convType_ << " gridding doesn't support non-square grid." << endl
444	<< "Result may be wrong." << LogIO::POST;
445	}
446
447	// grid parameter
448	os << LogIO::DEBUGGING ;
449	os << "----------" << endl ;
450	os << "Data selection summary" << endl ;
451	os << " ifno = " << ifno_ << endl ;
452	os << " pollist = " << pollist_ << endl ;
453	os << " scanlist = " << scanlist_ << endl ;
454	os << "----------" << endl ;
455	os << "Grid parameter summary" << endl ;
456	os << " (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
457	os << " (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
458	os << " center = " << center_ << endl ;
459	os << " weighting = " << wtype_ << endl ;
460	os << " convfunc = " << convType_ << endl;
461	if (convType_ == "GAUSS") {
462	os << " gwidth = " << gwidth_ << endl;
463	os << " truncate = " << truncate_ << endl;
464	}
465	else if (convType_ == "GJINC") {
466	os << " gwidth = " << gwidth_ << endl;
467	os << " jwidth = " << jwidth_ << endl;
468	os << " truncate = " << truncate_ << endl;
469	}
470	else {
471	os << " support = " << userSupport_ << endl;
472	}
473	os << " doclip = " << (doclip_?"True":"False") << endl ;
474	os << "----------" << LogIO::POST ;
475	os << LogIO::NORMAL ;
476
477	if ( doclip_ )
478	gridPerRowWithClipping() ;
479	else
480	gridPerRow() ;
481	}
482
483	void STGrid::updateChunkShape()
484	{
485	// TODO: nchunk_ must be determined from nchan_, npol_, and (nx_,ny_)
486	// by considering data size to be allocated for ggridsd input/output
487	nchunk_ = 400 ;
488	nchunk_ = min( nchunk_, nrow_ ) ;
489	vshape_ = IPosition( 1, nchunk_ ) ;
490	wshape_ = IPosition( 2, nchan_, nchunk_ ) ;
491	dshape_ = IPosition( 2, 2, nchunk_ ) ;
492	}
493
494	struct STGChunk {
495	Int nrow ;
496	Array<Complex> spectra;
497	Array<Int> flagtra;
498	Array<Int> rflag;
499	Array<Float> weight;
500	Array<Double> direction;
501	STGChunk(IPosition const &wshape, IPosition const &vshape,
502	IPosition const &dshape)
503	: spectra(wshape), flagtra(wshape), rflag(vshape), weight(wshape),
504	direction(dshape)
505	{ }
506	};
507
508	struct STCommonData {
509	Int gnx;
510	Int gny;
511	Int *chanMap;
512	Vector<Float> convFunc ;
513	Array<Complex> gdataArrC;
514	Array<Float> gwgtArr;
515	STCommonData(IPosition const &gshape, Array<Float> const &data)
516	: gdataArrC(gshape, 0.0), gwgtArr(data) {}
517	};
518
519	struct STCommonDataWithClipping {
520	Int gnx;
521	Int gny;
522	Int *chanMap;
523	Vector<Float> convFunc ;
524	Array<Complex> gdataArrC;
525	Array<Float> gwgtArr;
526	Array<Int> npoints ;
527	Array<Complex> clipMin ;
528	Array<Float> clipWMin ;
529	Array<Float> clipCMin ;
530	Array<Complex> clipMax ;
531	Array<Float> clipWMax ;
532	Array<Float> clipCMax ;
533	STCommonDataWithClipping(IPosition const &gshape,
534	IPosition const &pshape,
535	Array<Float> const &data)
536	: gdataArrC(gshape, 0.0),
537	gwgtArr(data),
538	npoints(pshape, 0),
539	clipMin(gshape, Complex(FLT_MAX,0.0)),
540	clipWMin(gshape, 0.0),
541	clipCMin(gshape, 0.0),
542	clipMax(gshape, Complex(-FLT_MAX,0.0)),
543	clipWMax(gshape, 0.0),
544	clipCMax(gshape, 0.0)
545	{}
546	};
547
548	#define DO_AHEAD 3
549
550	struct STContext {
551	STCommonData &common;
552	FIFO<STGChunk *, DO_AHEAD> queue;
553	STGrid *const self;
554	const Int pol;
555	STContext(STGrid *obj, STCommonData &common, Int pol)
556	: self(obj), common(common), pol(pol) {}
557	};
558
559	struct STContextWithClipping {
560	STCommonDataWithClipping &common;
561	FIFO<STGChunk *, DO_AHEAD> queue;
562	STGrid *const self;
563	const Int pol;
564	STContextWithClipping(STGrid *obj, STCommonDataWithClipping &common, Int pol)
565	: self(obj), common(common), pol(pol) {}
566	};
567
568
569	bool STGrid::produceChunk(void *ctx) throw(PCException)
570	{
571	STContext &context = (STContext )ctx;
572	if ( context.self->nprocessed_ >= context.self->nrow_ ) {
573	return false;
574	}
575	STGChunk *chunk = new STGChunk(context.self->wshape_,
576	context.self->vshape_,
577	context.self->dshape_);
578
579	double t0 = mathutil::gettimeofday_sec() ;
580	chunk->nrow = context.self->getDataChunk(
581	context.self->wshape_, context.self->vshape_, context.self->dshape_,
582	chunk->spectra, chunk->direction,
583	chunk->flagtra, chunk->rflag, chunk->weight);
584	double t1 = mathutil::gettimeofday_sec() ;
585	context.self->eGetData_ += t1-t0 ;
586
587	context.queue.lock();
588	context.queue.put(chunk);
589	context.queue.unlock();
590	return true;
591	}
592
593	void STGrid::consumeChunk(void *ctx) throw(PCException)
594	{
595	STContext &context = (STContext )ctx;
596	STGChunk *chunk = NULL;
597	try {
598	context.queue.lock();
599	chunk = context.queue.get();
600	context.queue.unlock();
601	} catch (FullException &e) {
602	context.queue.unlock();
603	// TODO: log error
604	throw PCException();
605	}
606
607	double t0, t1 ;
608	// world -> pixel
609	Array<Double> xypos( context.self->dshape_ ) ;
610	t0 = mathutil::gettimeofday_sec() ;
611	context.self->toPixel( chunk->direction, xypos ) ;
612	t1 = mathutil::gettimeofday_sec() ;
613	context.self->eToPixel_ += t1-t0 ;
614
615	// call ggridsd
616	Int nvispol = 1 ;
617	Int irow = -1 ;
618	t0 = mathutil::gettimeofday_sec() ;
619	context.self->call_ggridsd( xypos,
620	chunk->spectra,
621	nvispol,
622	context.self->nchan_,
623	chunk->flagtra,
624	chunk->rflag,
625	chunk->weight,
626	chunk->nrow,
627	irow,
628	context.common.gdataArrC,
629	context.common.gwgtArr,
630	context.common.gnx,
631	context.common.gny,
632	context.self->npol_,
633	context.self->nchan_,
634	context.self->convSupport_,
635	context.self->convSampling_,
636	context.common.convFunc,
637	context.common.chanMap,
638	(Int*)&context.pol ) ;
639	t1 = mathutil::gettimeofday_sec() ;
640	context.self->eGGridSD_ += t1-t0 ;
641
642	delete chunk;
643	}
644
645	void STGrid::gridPerRow()
646	{
647	LogIO os( LogOrigin("STGrid", "gridPerRow", WHERE) ) ;
648	double t0, t1 ;
649
650
651	// grid data
652	// Extend grid plane with convSupport_
653	// Int gnx = nx_+convSupport_*2 ;
654	// Int gny = ny_+convSupport_*2 ;
655	Int gnx = nx_;
656	Int gny = ny_;
657
658	IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
659	// 2011/12/20 TN
660	// data_ and gwgtArr share storage
661	data_.resize( gshape ) ;
662	data_ = 0.0 ;
663	STCommonData common = STCommonData(gshape, data_);
664	common.gnx = gnx ;
665	common.gny = gny ;
666
667	// parameters for gridding
668	Int *chanMap = new Int[nchan_] ;
669	for ( Int i = 0 ; i < nchan_ ; i++ ) {
670	chanMap[i] = i ;
671	}
672	common.chanMap = chanMap;
673
674	// convolution kernel
675	t0 = mathutil::gettimeofday_sec() ;
676	setConvFunc( common.convFunc ) ;
677	t1 = mathutil::gettimeofday_sec() ;
678	os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
679
680	// for performance check
681	eGetData_ = 0.0 ;
682	eToPixel_ = 0.0 ;
683	eGGridSD_ = 0.0 ;
684	double eInitPol = 0.0 ;
685
686	for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
687	initTable( ifile ) ;
688
689	os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;
690	Broker broker = Broker(produceChunk, consumeChunk);
691	for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
692	t0 = mathutil::gettimeofday_sec() ;
693	initPol( ipol ) ; // set ptab_ and attach()
694	t1 = mathutil::gettimeofday_sec() ;
695	eInitPol += t1-t0 ;
696
697	STContext context(this, common, ipol);
698
699	os << "start pol " << ipol << LogIO::POST ;
700
701	nprocessed_ = 0 ;
702	#if 1
703	broker.runProducerAsMasterThread(&context, DO_AHEAD);
704	#else
705	for (;;) {
706	bool produced = produceChunk(&context);
707	if (! produced) {
708	break;
709	}
710	consumeChunk(&context);
711	}
712	#endif
713
714	os << "end pol " << ipol << LogIO::POST ;
715
716	}
717	os << "end table " << ifile << LogIO::POST ;
718	}
719	os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
720	os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
721	os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
722	os << LogIO::DEBUGGING << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
723	os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
724	os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
725
726	delete chanMap ;
727
728	// set data
729	setData( common.gdataArrC, common.gwgtArr ) ;
730
731	}
732
733	void STGrid::consumeChunkWithClipping(void *ctx) throw(PCException)
734	{
735	STContextWithClipping &context = (STContextWithClipping )ctx;
736	STGChunk *chunk = NULL;
737	try {
738	context.queue.lock();
739	chunk = context.queue.get();
740	context.queue.unlock();
741	} catch (FullException &e) {
742	context.queue.unlock();
743	// TODO: log error
744	throw PCException();
745	}
746
747	double t0, t1 ;
748	// world -> pixel
749	Array<Double> xypos( context.self->dshape_ ) ;
750	t0 = mathutil::gettimeofday_sec() ;
751	context.self->toPixel( chunk->direction, xypos ) ;
752	t1 = mathutil::gettimeofday_sec() ;
753	context.self->eToPixel_ += t1-t0 ;
754
755	// call ggridsd
756	Int nvispol = 1 ;
757	Int irow = -1 ;
758	t0 = mathutil::gettimeofday_sec() ;
759	context.self->call_ggridsd2( xypos,
760	chunk->spectra,
761	nvispol,
762	context.self->nchan_,
763	chunk->flagtra,
764	chunk->rflag,
765	chunk->weight,
766	chunk->nrow,
767	irow,
768	context.common.gdataArrC,
769	context.common.gwgtArr,
770	context.common.npoints,
771	context.common.clipMin,
772	context.common.clipWMin,
773	context.common.clipCMin,
774	context.common.clipMax,
775	context.common.clipWMax,
776	context.common.clipCMax,
777	context.common.gnx,
778	context.common.gny,
779	context.self->npol_,
780	context.self->nchan_,
781	context.self->convSupport_,
782	context.self->convSampling_,
783	context.common.convFunc,
784	context.common.chanMap,
785	(Int*)&context.pol ) ;
786	t1 = mathutil::gettimeofday_sec() ;
787	context.self->eGGridSD_ += t1-t0 ;
788
789	delete chunk;
790	}
791
792	void STGrid::gridPerRowWithClipping()
793	{
794	LogIO os( LogOrigin("STGrid", "gridPerRowWithClipping", WHERE) ) ;
795	double t0, t1 ;
796
797
798	// grid data
799	// Extend grid plane with convSupport_
800	// Int gnx = nx_+convSupport_*2 ;
801	// Int gny = ny_+convSupport_*2 ;
802	Int gnx = nx_;
803	Int gny = ny_;
804
805	IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
806	IPosition pshape( 3, gnx, gny, npol_ ) ;
807	// 2011/12/20 TN
808	// data_ and gwgtArr share storage
809	data_.resize( gshape ) ;
810	data_ = 0.0 ;
811	STCommonDataWithClipping common = STCommonDataWithClipping( gshape,
812	pshape,
813	data_ ) ;
814	common.gnx = gnx ;
815	common.gny = gny ;
816
817	// parameters for gridding
818	Int *chanMap = new Int[nchan_] ;
819	for ( Int i = 0 ; i < nchan_ ; i++ ) {
820	chanMap[i] = i ;
821	}
822	common.chanMap = chanMap;
823
824	// convolution kernel
825	t0 = mathutil::gettimeofday_sec() ;
826	setConvFunc( common.convFunc ) ;
827	t1 = mathutil::gettimeofday_sec() ;
828	os << LogIO::DEBUGGING << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
829
830	// for performance check
831	eGetData_ = 0.0 ;
832	eToPixel_ = 0.0 ;
833	eGGridSD_ = 0.0 ;
834	double eInitPol = 0.0 ;
835
836	for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
837	initTable( ifile ) ;
838
839	os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;
840	Broker broker = Broker(produceChunk, consumeChunkWithClipping);
841	for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
842	t0 = mathutil::gettimeofday_sec() ;
843	initPol( ipol ) ; // set ptab_ and attach()
844	t1 = mathutil::gettimeofday_sec() ;
845	eInitPol += t1-t0 ;
846
847	STContextWithClipping context(this, common, ipol);
848
849	os << "start pol " << ipol << LogIO::POST ;
850
851	nprocessed_ = 0 ;
852	#if 1
853	broker.runProducerAsMasterThread(&context, DO_AHEAD);
854	#else
855	for (;;) {
856	bool produced = produceChunk(&context);
857	if (! produced) {
858	break;
859	}
860	consumeChunkWithClipping(&context);
861	}
862	#endif
863
864	os << "end pol " << ipol << LogIO::POST ;
865
866	}
867	os << "end table " << ifile << LogIO::POST ;
868	}
869	os << LogIO::DEBUGGING << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
870	os << LogIO::DEBUGGING << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
871	os << LogIO::DEBUGGING << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
872	os << LogIO::DEBUGGING << "ggridsd2: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
873	os << LogIO::DEBUGGING << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
874	os << LogIO::DEBUGGING << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
875
876	delete chanMap ;
877
878	// clip min and max in each grid
879	// os << "BEFORE CLIPPING" << LogIO::POST ;
880	// os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
881	// os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
882	t0 = mathutil::gettimeofday_sec() ;
883	clipMinMax( common.gdataArrC,
884	common.gwgtArr,
885	common.npoints,
886	common.clipMin,
887	common.clipWMin,
888	common.clipCMin,
889	common.clipMax,
890	common.clipWMax,
891	common.clipCMax ) ;
892	t1 = mathutil::gettimeofday_sec() ;
893	os << LogIO::DEBUGGING << "clipMinMax: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
894	// os << "AFTER CLIPPING" << LogIO::POST ;
895	// os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
896	// os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
897
898	// set data
899	setData( common.gdataArrC, common.gwgtArr ) ;
900
901	}
902
903	void STGrid::clipMinMax( Array<Complex> &grid,
904	Array<Float> &weight,
905	Array<Int> &npoints,
906	Array<Complex> &clipmin,
907	Array<Float> &clipwmin,
908	Array<Float> &clipcmin,
909	Array<Complex> &clipmax,
910	Array<Float> &clipwmax,
911	Array<Float> &clipcmax )
912	{
913	//LogIO os( LogOrigin("STGrid","clipMinMax",WHERE) ) ;
914
915	// prepare pointers
916	Bool delG, delW, delNP, delCMin, delCWMin, delCCMin, delCMax, delCWMax, delCCMax ;
917	Complex *grid_p = grid.getStorage( delG ) ;
918	Float *wgt_p = weight.getStorage( delW ) ;
919	const Int *npts_p = npoints.getStorage( delNP ) ;
920	const Complex *cmin_p = clipmin.getStorage( delCMin ) ;
921	const Float *cwmin_p = clipwmin.getStorage( delCWMin ) ;
922	const Float *ccmin_p = clipcmin.getStorage( delCCMin ) ;
923	const Complex *cmax_p = clipmax.getStorage( delCMax ) ;
924	const Float *cwmax_p = clipwmax.getStorage( delCWMax ) ;
925	const Float *ccmax_p = clipcmax.getStorage( delCCMax ) ;
926
927	const IPosition &gshape = grid.shape() ;
928	long offset = gshape[0] * gshape[1] * gshape[2] ; // nx * ny * npol
929	Int nchan = gshape[3] ;
930	long origin = nchan * offset ;
931	for ( long i = 0 ; i < offset ; i++ ) {
932	if ( *npts_p > 2 ) {
933	for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
934	// clip minimum and maximum
935	grid_p -= (cmin_p)(cwmin_p)(ccmin_p)
936	+ (cmax_p)(cwmax_p)(*ccmax_p) ;
937	wgt_p -= (cwmin_p)(ccmin_p)
938	+ (cwmax_p)(*ccmax_p) ;
939
940	grid_p += offset ;
941	wgt_p += offset ;
942	cmin_p += offset ;
943	cwmin_p += offset ;
944	ccmin_p += offset ;
945	cmax_p += offset ;
946	cwmax_p += offset ;
947	ccmax_p += offset ;
948	}
949	grid_p -= origin ;
950	wgt_p -= origin ;
951	cmin_p -= origin ;
952	cwmin_p -= origin ;
953	ccmin_p -= origin ;
954	cmax_p -= origin ;
955	cwmax_p -= origin ;
956	ccmax_p -= origin ;
957	}
958	grid_p++ ;
959	wgt_p++ ;
960	npts_p++ ;
961	cmin_p++ ;
962	cwmin_p++ ;
963	ccmin_p++ ;
964	cmax_p++ ;
965	cwmax_p++ ;
966	ccmax_p++ ;
967	}
968	grid_p -= offset ;
969	wgt_p -= offset ;
970	npts_p -= offset ;
971	cmin_p -= offset ;
972	cwmin_p -= offset ;
973	ccmin_p -= offset ;
974	cmax_p -= offset ;
975	cwmax_p -= offset ;
976	ccmax_p -= offset ;
977
978	// finalization
979	grid.putStorage( grid_p, delG ) ;
980	weight.putStorage( wgt_p, delW ) ;
981	npoints.freeStorage( npts_p, delNP ) ;
982	clipmin.freeStorage( cmin_p, delCMin ) ;
983	clipwmin.freeStorage( cwmin_p, delCWMin ) ;
984	clipcmin.freeStorage( ccmin_p, delCCMin ) ;
985	clipmax.freeStorage( cmax_p, delCMax ) ;
986	clipwmax.freeStorage( cwmax_p, delCWMax ) ;
987	clipcmax.freeStorage( ccmax_p, delCCMax ) ;
988	}
989
990	void STGrid::initPol( Int ipol )
991	{
992	LogIO os( LogOrigin("STGrid","initPol",WHERE) ) ;
993	if ( npolOrg_ == 1 ) {
994	os << "single polarization data." << LogIO::POST ;
995	ptab_ = tab_ ;
996	}
997	else
998	ptab_ = tab_( tab_.col("POLNO") == pollist_[ipol] ) ;
999
1000	attach( ptab_ ) ;
1001	}
1002
1003	void STGrid::initTable( uInt idx )
1004	{
1005	tab_ = tableList_[idx] ;
1006	nrow_ = rows_[idx] ;
1007	updateChunkShape() ;
1008	}
1009
1010	void STGrid::setData( Array<Complex> &gdata,
1011	Array<Float> &gwgt )
1012	{
1013	// 2011/12/20 TN
1014	// gwgt and data_ share storage
1015	LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
1016	double t0, t1 ;
1017	t0 = mathutil::gettimeofday_sec() ;
1018	uInt len = data_.nelements() ;
1019	const Complex *w1_p ;
1020	Float *w2_p ;
1021	Bool b1, b2 ;
1022	const Complex *gdata_p = gdata.getStorage( b1 ) ;
1023	Float *gwgt_p = data_.getStorage( b2 ) ;
1024	w1_p = gdata_p ;
1025	w2_p = gwgt_p ;
1026	for ( uInt i = 0 ; i < len ; i++ ) {
1027	if ( w2_p > 0.0 ) w2_p = (w1_p).real() / w2_p ;
1028	w1_p++ ;
1029	w2_p++ ;
1030	}
1031	gdata.freeStorage( gdata_p, b1 ) ;
1032	data_.putStorage( gwgt_p, b2 ) ;
1033	t1 = mathutil::gettimeofday_sec() ;
1034	os << LogIO::DEBUGGING << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1035	}
1036
1037	void STGrid::setupGrid()
1038	{
1039	Double xmin,xmax,ymin,ymax ;
1040	mapExtent( xmin, xmax, ymin, ymax ) ;
1041
1042	setupGrid( nxUI_, nyUI_, cellxUI_, cellyUI_,
1043	xmin, xmax, ymin, ymax, centerUI_ ) ;
1044	}
1045
1046	void STGrid::setupGrid( Int &nx,
1047	Int &ny,
1048	String &cellx,
1049	String &celly,
1050	Double &xmin,
1051	Double &xmax,
1052	Double &ymin,
1053	Double &ymax,
1054	String &center )
1055	{
1056	LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
1057	//cout << "nx=" << nx << ", ny=" << ny << endl ;
1058
1059	// center position
1060	if ( center.size() == 0 ) {
1061	center_(0) = 0.5 * ( xmin + xmax ) ;
1062	center_(1) = 0.5 * ( ymin + ymax ) ;
1063	}
1064	else {
1065	String::size_type pos0 = center.find( " " ) ;
1066	if ( pos0 == String::npos ) {
1067	throw AipsError( "bad string format in parameter center" ) ;
1068	}
1069	String::size_type pos1 = center.find( " ", pos0+1 ) ;
1070	String typestr, xstr, ystr ;
1071	if ( pos1 != String::npos ) {
1072	typestr = center.substr( 0, pos0 ) ;
1073	xstr = center.substr( pos0+1, pos1-pos0 ) ;
1074	ystr = center.substr( pos1+1 ) ;
1075	// todo: convert to J2000 (or direction ref for DIRECTION column)
1076	}
1077	else {
1078	typestr = "J2000" ;
1079	xstr = center.substr( 0, pos0 ) ;
1080	ystr = center.substr( pos0+1 ) ;
1081	}
1082	QuantumHolder qh ;
1083	String err ;
1084	qh.fromString( err, xstr ) ;
1085	Quantum<Double> xcen = qh.asQuantumDouble() ;
1086	qh.fromString( err, ystr ) ;
1087	Quantum<Double> ycen = qh.asQuantumDouble() ;
1088	center_(0) = xcen.getValue( "rad" ) ;
1089	center_(1) = ycen.getValue( "rad" ) ;
1090	double base = 0.5 * (xmin + xmax) ;
1091	int maxrotate = 1 ;
1092	int nelem = 2 * maxrotate + 1 ;
1093	double *sep = new double[nelem] ;
1094	for ( int i = 0 ; i < nelem ; i++ )
1095	sep[i] = abs(base - center_[0] - (i-maxrotate) * C::_2pi) ;
1096	// os << "sep[0]=" << sep[0] << endl
1097	// << "sep[1]=" << sep[1] << endl
1098	// << "sep[2]=" << sep[2] << LogIO::POST ;
1099	int idx = 0 ;
1100	base = sep[0] ;
1101	int nrotate = 0 ;
1102	while ( idx < nelem ) {
1103	if ( base > sep[idx] ) {
1104	base = sep[idx] ;
1105	nrotate = idx ;
1106	}
1107	idx++ ;
1108	}
1109	delete sep ;
1110	nrotate -= maxrotate ;
1111	// os << "nrotate = " << nrotate << LogIO::POST ;
1112	center_[0] += nrotate * C::_2pi ;
1113	}
1114	// os << "xmin=" << xmin << LogIO::POST ;
1115	// os << "center_=" << center_ << LogIO::POST ;
1116
1117	nx_ = nx ;
1118	ny_ = ny ;
1119	if ( nx < 0 && ny > 0 ) {
1120	nx_ = ny ;
1121	ny_ = ny ;
1122	}
1123	if ( ny < 0 && nx > 0 ) {
1124	nx_ = nx ;
1125	ny_ = nx ;
1126	}
1127
1128	//Double wx = xmax - xmin ;
1129	//Double wy = ymax - ymin ;
1130	Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
1131	Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
1132	// take 10% margin
1133	wx *= 1.10 ;
1134	wy *= 1.10 ;
1135
1136	Quantum<Double> qcellx ;
1137	Quantum<Double> qcelly ;
1138	//cout << "nx_ = " << nx_ << ", ny_ = " << ny_ << endl ;
1139	if ( cellx.size() != 0 && celly.size() != 0 ) {
1140	readQuantity( qcellx, cellx ) ;
1141	readQuantity( qcelly, celly ) ;
1142	}
1143	else if ( celly.size() != 0 ) {
1144	os << "Using celly to x-axis..." << LogIO::POST ;
1145	readQuantity( qcelly, celly ) ;
1146	qcellx = qcelly ;
1147	}
1148	else if ( cellx.size() != 0 ) {
1149	os << "Using cellx to y-axis..." << LogIO::POST ;
1150	readQuantity( qcellx, cellx ) ;
1151	qcelly = qcellx ;
1152	}
1153	else {
1154	if ( nx_ < 0 ) {
1155	os << "No user preference in grid setting. Using default..." << LogIO::POST ;
1156	readQuantity( qcellx, "1.0arcmin" ) ;
1157	qcelly = qcellx ;
1158	}
1159	else {
1160	if ( wx == 0.0 ) {
1161	os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1162	wx = 0.00290888 ;
1163	}
1164	if ( wy == 0.0 ) {
1165	os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1166	wy = 0.00290888 ;
1167	}
1168	qcellx = Quantum<Double>( wx/nx_*cos(center_[1]), "rad" ) ;
1169	qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
1170	}
1171	}
1172	cellx_ = qcellx.getValue( "rad" ) ;
1173	celly_ = qcelly.getValue( "rad" ) ;
1174	//os << "cellx_=" << cellx_ << ", celly_=" << celly_ << ", cos("<<center_(1)<<")=" << cos(center_(1)) << LogIO::POST ;
1175	if ( nx_ < 0 ) {
1176	if ( wx == 0.0 ) {
1177	os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1178	wx = 0.00290888 ;
1179	}
1180	if ( wy == 0.0 ) {
1181	os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1182	wy = 0.00290888 ;
1183	}
1184	nx_ = Int( ceil( wx/(cellx_/cos(center_[1])) ) ) ;
1185	ny_ = Int( ceil( wy/celly_ ) ) ;
1186	}
1187
1188	// create DirectionCoordinate
1189	Matrix<Double> xform(2,2) ;
1190	xform = 0.0 ;
1191	xform.diagonal() = 1.0 ;
1192	dircoord_ = new DirectionCoordinate(MDirection::J2000,
1193	Projection( Projection::SIN ),
1194	center_[0], center_[1],
1195	-cellx_, celly_,
1196	xform,
1197	0.5*Double(nx_-1),
1198	0.5*Double(ny_-1)) ;
1199	}
1200
1201	void STGrid::mapExtent( Double &xmin, Double &xmax,
1202	Double &ymin, Double &ymax )
1203	{
1204	//LogIO os( LogOrigin("STGrid","mapExtent",WHERE) ) ;
1205	directionCol_.attach( tableList_[0], "DIRECTION" ) ;
1206	Matrix<Double> direction = directionCol_.getColumn() ;
1207	//os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1208	Vector<Double> ra( direction.row(0) ) ;
1209	mathutil::rotateRA( ra ) ;
1210	minMax( xmin, xmax, direction.row( 0 ) ) ;
1211	minMax( ymin, ymax, direction.row( 1 ) ) ;
1212	Double amin, amax, bmin, bmax ;
1213	for ( uInt i = 1 ; i < nfile_ ; i++ ) {
1214	directionCol_.attach( tableList_[i], "DIRECTION" ) ;
1215	direction.assign( directionCol_.getColumn() ) ;
1216	//os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1217	// to make contiguous RA distribution (no 2pi jump)
1218	Vector<Double> ra( direction.row(0) ) ;
1219	mathutil::rotateRA( ra ) ;
1220	minMax( amin, amax, direction.row( 0 ) ) ;
1221	minMax( bmin, bmax, direction.row( 1 ) ) ;
1222	xmin = min( xmin, amin ) ;
1223	xmax = max( xmax, amax ) ;
1224	ymin = min( ymin, bmin ) ;
1225	ymax = max( ymax, bmax ) ;
1226	}
1227	//os << "(xmin,xmax)=(" << xmin << "," << xmax << ")" << LogIO::POST ;
1228	//os << "(ymin,ymax)=(" << ymin << "," << ymax << ")" << LogIO::POST ;
1229	}
1230
1231	void STGrid::table( Table &tab, uInt i )
1232	{
1233	if ( i >= 0 && i < nfile_ )
1234	tab = Table( infileList_[i] ) ;
1235	}
1236
1237	void STGrid::selectData()
1238	{
1239	LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
1240	Int ifno = ifno_ ;
1241	tableList_.resize( nfile_ ) ;
1242	if ( ifno_ == -1 ) {
1243	//Table taborg( infileList_[0] ) ;
1244	Table taborg ;
1245	table( taborg, 0 ) ;
1246	ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
1247	ifno_ = ifnoCol( 0 ) ;
1248	os << LogIO::WARN
1249	<< "IFNO is not given. Using default IFNO: " << ifno_ << LogIO::POST ;
1250	}
1251	for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1252	//Table taborg( infileList_[i] ) ;
1253	Table taborg ;
1254	table( taborg, i ) ;
1255	TableExprNode node ;
1256	if ( ifno != -1 \|\| isMultiIF( taborg ) ) {
1257	os << "apply selection on IFNO" << LogIO::POST ;
1258	node = taborg.col("IFNO") == ifno_ ;
1259	}
1260	if ( scanlist_.size() > 0 ) {
1261	os << "apply selection on SCANNO" << LogIO::POST ;
1262	node = node && taborg.col("SCANNO").in( scanlist_ ) ;
1263	}
1264	if ( node.isNull() ) {
1265	tableList_[i] = taborg ;
1266	}
1267	else {
1268	tableList_[i] = taborg( node ) ;
1269	}
1270	os << LogIO::DEBUGGING << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ;
1271	if ( tableList_[i].nrow() == 0 ) {
1272	os << LogIO::SEVERE
1273	<< "No corresponding rows for given selection: IFNO " << ifno_ ;
1274	if ( scanlist_.size() > 0 )
1275	os << " SCANNO " << scanlist_ ;
1276	os << LogIO::EXCEPTION ;
1277	}
1278	}
1279	}
1280
1281	Bool STGrid::isMultiIF( Table &tab )
1282	{
1283	ROScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
1284	Vector<uInt> ifnos = ifnoCol.getColumn() ;
1285	return anyNE( ifnos, ifnos[0] ) ;
1286	}
1287
1288	void STGrid::attach( Table &tab )
1289	{
1290	// attach to table
1291	spectraCol_.attach( tab, "SPECTRA" ) ;
1292	flagtraCol_.attach( tab, "FLAGTRA" ) ;
1293	directionCol_.attach( tab, "DIRECTION" ) ;
1294	flagRowCol_.attach( tab, "FLAGROW" ) ;
1295	tsysCol_.attach( tab, "TSYS" ) ;
1296	intervalCol_.attach( tab, "INTERVAL" ) ;
1297	}
1298
1299	Int STGrid::getDataChunk(
1300	IPosition const &wshape,
1301	IPosition const &vshape,
1302	IPosition const &dshape,
1303	Array<Complex> &spectra,
1304	Array<Double> &direction,
1305	Array<Int> &flagtra,
1306	Array<Int> &rflag,
1307	Array<Float> &weight )
1308	{
1309	LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1310
1311	Array<Float> spectraF_(wshape);
1312	Array<uChar> flagtraUC_(wshape);
1313	Array<uInt> rflagUI_(vshape);
1314	Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1315	if ( nrow < nchunk_ ) {
1316	spectra.resize( spectraF_.shape() ) ;
1317	flagtra.resize( flagtraUC_.shape() ) ;
1318	rflag.resize( rflagUI_.shape() ) ;
1319	}
1320	double t0, t1 ;
1321	t0 = mathutil::gettimeofday_sec() ;
1322	convertArray( spectra, spectraF_ ) ;
1323	toInt( flagtraUC_, flagtra ) ;
1324	toInt( rflagUI_, rflag ) ;
1325	t1 = mathutil::gettimeofday_sec() ;
1326	eToInt = t1 - t0 ;
1327
1328	return nrow ;
1329	}
1330
1331	#if 0
1332	Int STGrid::getDataChunk( Array<Complex> &spectra,
1333	Array<Double> &direction,
1334	Array<Int> &flagtra,
1335	Array<Int> &rflag,
1336	Array<Float> &weight )
1337	{
1338	LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1339	Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1340	if ( nrow < nchunk_ ) {
1341	spectra.resize( spectraF_.shape() ) ;
1342	flagtra.resize( flagtraUC_.shape() ) ;
1343	rflag.resize( rflagUI_.shape() ) ;
1344	}
1345	double t0, t1 ;
1346	t0 = mathutil::gettimeofday_sec() ;
1347	convertArray( spectra, spectraF_ ) ;
1348	toInt( flagtraUC_, flagtra ) ;
1349	toInt( rflagUI_, rflag ) ;
1350	t1 = mathutil::gettimeofday_sec() ;
1351	eToInt = t1 - t0 ;
1352
1353	return nrow ;
1354	}
1355	#endif
1356
1357	Int STGrid::getDataChunk( Array<Float> &spectra,
1358	Array<Double> &direction,
1359	Array<uChar> &flagtra,
1360	Array<uInt> &rflag,
1361	Array<Float> &weight )
1362	{
1363	LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1364	Int nrow = spectra.shape()[1] ;
1365	Int remainingRow = nrow_ - nprocessed_ ;
1366	if ( remainingRow < nrow ) {
1367	nrow = remainingRow ;
1368	IPosition mshape( 2, nchan_, nrow ) ;
1369	IPosition vshape( 1, nrow ) ;
1370	spectra.resize( mshape ) ;
1371	flagtra.resize( mshape ) ;
1372	direction.resize( IPosition(2,2,nrow) ) ;
1373	rflag.resize( vshape ) ;
1374	weight.resize( mshape ) ;
1375	}
1376	// 2011/12/22 TN
1377	// tsys shares its storage with weight
1378	Array<Float> tsys( weight ) ;
1379	Array<Double> tint( rflag.shape() ) ;
1380
1381	Vector<uInt> rflagVec( rflag ) ;
1382	Vector<Double> tintVec( tint ) ;
1383
1384	RefRows rows( nprocessed_, nprocessed_+nrow-1, 1 ) ;
1385	//os<<LogIO::DEBUGGING<<"nprocessed_="<<nprocessed_<<": rows.nrows()="<<rows.nrows()<<LogIO::POST ;
1386	spectraCol_.getColumnCells( rows, spectra ) ;
1387	flagtraCol_.getColumnCells( rows, flagtra ) ;
1388	directionCol_.getColumnCells( rows, direction ) ;
1389	// to make contiguous RA distribution (no 2pi jump)
1390	Vector<Double> v( Matrix<Double>(direction).row(0) ) ;
1391	mathutil::rotateRA( v ) ;
1392	flagRowCol_.getColumnCells( rows, rflagVec ) ;
1393	intervalCol_.getColumnCells( rows, tintVec ) ;
1394	Vector<Float> tsysTemp = tsysCol_( nprocessed_ ) ;
1395	if ( tsysTemp.nelements() == (uInt)nchan_ )
1396	tsysCol_.getColumnCells( rows, tsys ) ;
1397	else
1398	tsys = tsysTemp[0] ;
1399
1400	double t0,t1 ;
1401	t0 = mathutil::gettimeofday_sec() ;
1402	getWeight( weight, tsys, tint ) ;
1403	t1 = mathutil::gettimeofday_sec() ;
1404	eGetWeight += t1-t0 ;
1405
1406	nprocessed_ += nrow ;
1407
1408	return nrow ;
1409	}
1410
1411	void STGrid::setupArray()
1412	{
1413	LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
1414	ROScalarColumn<uInt> polnoCol( tableList_[0], "POLNO" ) ;
1415	Vector<uInt> pols = polnoCol.getColumn() ;
1416	//os << pols << LogIO::POST ;
1417	Vector<uInt> pollistOrg ;
1418	npolOrg_ = 0 ;
1419	uInt polno ;
1420	for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
1421	//polno = polnoCol( i ) ;
1422	polno = pols( i ) ;
1423	if ( allNE( pollistOrg, polno ) ) {
1424	pollistOrg.resize( npolOrg_+1, True ) ;
1425	pollistOrg[npolOrg_] = polno ;
1426	npolOrg_++ ;
1427	}
1428	}
1429	if ( pollist_.size() == 0 )
1430	pollist_ = pollistOrg ;
1431	else {
1432	Vector<uInt> newlist ;
1433	uInt newsize = 0 ;
1434	for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
1435	if ( anyEQ( pollistOrg, pollist_[i] ) ) {
1436	newlist.resize( newsize+1, True ) ;
1437	newlist[newsize] = pollist_[i] ;
1438	newsize++ ;
1439	}
1440	}
1441	pollist_.assign( newlist ) ;
1442	}
1443	npol_ = pollist_.size() ;
1444	if ( npol_ == 0 ) {
1445	os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
1446	}
1447	rows_.resize( nfile_ ) ;
1448	for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1449	rows_[i] = tableList_[i].nrow() / npolOrg_ ;
1450	//if ( nrow_ < rows_[i] )
1451	// nrow_ = rows_[i] ;
1452	}
1453	flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ;
1454	nchan_ = flagtraCol_( 0 ).nelements() ;
1455	// os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
1456	// << "nchan_ = " << nchan_ << endl
1457	// << "nrow_ = " << nrow_ << LogIO::POST ;
1458	}
1459
1460	void STGrid::getWeight( Array<Float> &w,
1461	Array<Float> &tsys,
1462	Array<Double> &tint )
1463	{
1464	LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
1465
1466	// 2011/12/22 TN
1467	// w (weight) and tsys share storage
1468	IPosition refShape = tsys.shape() ;
1469	Int nchan = refShape[0] ;
1470	Int nrow = refShape[1] ;
1471	// os << "nchan=" << nchan << ", nrow=" << nrow << LogIO::POST ;
1472	// os << "w.shape()=" << w.shape() << endl
1473	// << "tsys.shape()=" << tsys.shape() << endl
1474	// << "tint.shape()=" << tint.shape() << LogIO::POST ;
1475
1476	// set weight
1477	if ( wtype_.compare( "UNIFORM" ) == 0 ) {
1478	w = 1.0 ;
1479	}
1480	else if ( wtype_.compare( "TINT" ) == 0 ) {
1481	Bool b0, b1 ;
1482	Float *w_p = w.getStorage( b0 ) ;
1483	Float *w0_p = w_p ;
1484	const Double *ti_p = tint.getStorage( b1 ) ;
1485	const Double *w1_p = ti_p ;
1486	for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1487	for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1488	w0_p = w1_p ;
1489	w0_p++ ;
1490	}
1491	w1_p++ ;
1492	}
1493	w.putStorage( w_p, b0 ) ;
1494	tint.freeStorage( ti_p, b1 ) ;
1495	}
1496	else if ( wtype_.compare( "TSYS" ) == 0 ) {
1497	Bool b0 ;
1498	Float *w_p = w.getStorage( b0 ) ;
1499	Float *w0_p = w_p ;
1500	for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1501	for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1502	Float temp = *w0_p ;
1503	w0_p = 1.0 / ( temp temp ) ;
1504	w0_p++ ;
1505	}
1506	}
1507	w.putStorage( w_p, b0 ) ;
1508	}
1509	else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
1510	Bool b0, b1 ;
1511	Float *w_p = w.getStorage( b0 ) ;
1512	Float *w0_p = w_p ;
1513	const Double *ti_p = tint.getStorage( b1 ) ;
1514	const Double *w1_p = ti_p ;
1515	for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1516	Float interval = *w1_p ;
1517	for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1518	Float temp = *w0_p ;
1519	w0_p = interval / ( temp temp ) ;
1520	w0_p++ ;
1521	}
1522	w1_p++ ;
1523	}
1524	w.putStorage( w_p, b0 ) ;
1525	tint.freeStorage( ti_p, b1 ) ;
1526	}
1527	else {
1528	//LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
1529	//os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
1530	w = 1.0 ;
1531	}
1532	}
1533
1534	void STGrid::toInt( Array<uChar> &u, Array<Int> &v )
1535	{
1536	uInt len = u.nelements() ;
1537	Int *int_p = new Int[len] ;
1538	Bool deleteIt ;
1539	const uChar *data_p = u.getStorage( deleteIt ) ;
1540	Int *i_p = int_p ;
1541	const uChar *u_p = data_p ;
1542	for ( uInt i = 0 ; i < len ; i++ ) {
1543	i_p = ( u_p == 0 ) ? 0 : 1 ;
1544	i_p++ ;
1545	u_p++ ;
1546	}
1547	u.freeStorage( data_p, deleteIt ) ;
1548	v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1549	}
1550
1551	void STGrid::toInt( Array<uInt> &u, Array<Int> &v )
1552	{
1553	uInt len = u.nelements() ;
1554	Int *int_p = new Int[len] ;
1555	Bool deleteIt ;
1556	const uInt *data_p = u.getStorage( deleteIt ) ;
1557	Int *i_p = int_p ;
1558	const uInt *u_p = data_p ;
1559	for ( uInt i = 0 ; i < len ; i++ ) {
1560	i_p = ( u_p == 0 ) ? 0 : 1 ;
1561	i_p++ ;
1562	u_p++ ;
1563	}
1564	u.freeStorage( data_p, deleteIt ) ;
1565	v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1566	}
1567
1568	void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
1569	{
1570	uInt nrow = world.shape()[1] ;
1571	Bool bw, bp ;
1572	Double *w_p = world.getStorage( bw ) ;
1573	Double *p_p = pixel.getStorage( bp ) ;
1574	Double *ww_p = w_p ;
1575	Double *wp_p = p_p ;
1576	IPosition vshape( 1, 2 ) ;
1577	Vector<Double> _world, _pixel ;
1578	for ( uInt i = 0 ; i < nrow ; i++ ) {
1579	_world.takeStorage( vshape, ww_p, SHARE ) ;
1580	_pixel.takeStorage( vshape, wp_p, SHARE ) ;
1581	dircoord_->toPixel( _pixel, _world ) ;
1582	ww_p += 2 ;
1583	wp_p += 2 ;
1584	}
1585	world.putStorage( w_p, bw ) ;
1586	pixel.putStorage( p_p, bp ) ;
1587	}
1588
1589	void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
1590	{
1591	convFunc = 0.0 ;
1592	for ( Int i = 0 ; i < convSize/2 ; i++ )
1593	convFunc(i) = 1.0 ;
1594	}
1595
1596	#define NEED_UNDERSCORES
1597	#if defined(NEED_UNDERSCORES)
1598	#define grdsf grdsf_
1599	#define grdgauss grdgauss_
1600	#define grdjinc1 grdjinc1_
1601	#endif
1602	#if defined(USE_CASAPY)
1603	extern "C" {
1604	void grdsf(Double, Double);
1605	void grdgauss(Double, Double, Double*);
1606	void grdjinc1(Double, Double, Int, Double);
1607	}
1608	#else
1609	extern "C" {
1610	void grdsf(Double, Double);
1611	}
1612	void grdgauss(Double hwhm, Double val, Double *out)
1613	{
1614	out = exp(-log(2.0) (val / hwhm) * (val / hwhm));
1615	}
1616	void grdjinc1(Double c, Double val, Int normalize, Double out)
1617	{
1618	// Calculate J_1(x) using approximate formula
1619	Double x = C::pi * val / c;
1620	Double ax = fabs(x);
1621	Double ans;
1622	if ( ax < 8.0 ) {
1623	Double y = x * x;
1624	Double ans1 = x * (72362614232.0 + y * (-7895059235.0
1625	+ y * (242396853.1 + y * (-2972611.439
1626	+ y * (15704.48260 + y * (-30.16036606))))));
1627	Double ans2 = 144725228442.0 + y * (2300535178.0
1628	+ y * (18583304.74 + y * (99447.43394
1629	+ y * (376.9991397 + y * 1.0))));
1630	ans = ans1 / ans2;
1631	}
1632	else {
1633	Double z = 8.0 / ax;
1634	Double y = z * z;
1635	Double xx = ax - 2.356194491;
1636	Double ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
1637	+ y * (0.2457520174e-5 + y * (-0.240337019e-6))));
1638	Double ans2 = 0.04687499995 + y * (-0.2002690873e-3
1639	+ y * (0.8449199096e-5 + y * (-0.88228987e-6
1640	+ y * (0.105787412e-6))));
1641	ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1
1642	- z * sin(xx) * ans2);
1643	if (x < 0.0)
1644	ans = -ans;
1645	}
1646
1647	// Then, calculate Jinc
1648	if (x == 0.0) {
1649	*out = 0.5;
1650	}
1651	else {
1652	*out = ans / x;
1653	}
1654
1655	if (*normalize == 1)
1656	out = out / 0.5;
1657	}
1658	#endif
1659	void STGrid::spheroidalFunc( Vector<Float> &convFunc )
1660	{
1661	convFunc = 0.0 ;
1662	for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
1663	Double nu = Double(i) / Double(convSupport_*convSampling_) ;
1664	Double val ;
1665	grdsf( &nu, &val ) ;
1666	convFunc(i) = ( 1.0 - nu * nu ) * val ;
1667	}
1668	}
1669
1670	void STGrid::gaussFunc( Vector<Float> &convFunc, Double hwhm, Double truncate )
1671	{
1672	convFunc = 0.0 ;
1673	Int len = (Int)(truncate*Double(convSampling_)+0.5);
1674	Double out, val;
1675	for ( Int i = 0 ; i < len ; i++ ) {
1676	val = Double(i) / Double(convSampling_) ;
1677	grdgauss(&hwhm, &val, &out);
1678	convFunc(i) = out;
1679	}
1680	}
1681
1682	void STGrid::gjincFunc( Vector<Float> &convFunc, Double hwhm, Double c, Double truncate )
1683	{
1684	convFunc = 0.0;
1685	Double out1, out2, val;
1686	Int normalize = 1;
1687	if (truncate >= 0.0) {
1688	Int len = (Int)(truncate*Double(convSampling_)+0.5);
1689	for (Int i = 0 ; i < len ; i++) {
1690	val = Double(i) / Double(convSampling_);
1691	grdgauss(&hwhm, &val, &out1);
1692	grdjinc1(&c, &val, &normalize, &out2);
1693	convFunc(i) = out1 * out2;
1694	}
1695	}
1696	else {
1697	Int len = convFunc.nelements();
1698	for (Int i = 0 ; i < len ; i++) {
1699	val = Double(i) / Double(convSampling_);
1700	grdjinc1(&c, &val, &normalize, &out2);
1701	if (out2 <= 0.0) {
1702	LogIO os(LogOrigin("STGrid","gjincFunc",WHERE));
1703	os << LogIO::DEBUG1 << "convFunc is automatically truncated at radius " << val << LogIO::POST;
1704	break;
1705	}
1706	grdgauss(&hwhm, &val, &out1);
1707	convFunc(i) = out1 * out2;
1708	}
1709	}
1710	}
1711
1712	void STGrid::pbFunc( Vector<Float> &convFunc )
1713	{
1714	convFunc = 0.0 ;
1715	}
1716
1717	vector<float> STGrid::getConvFunc()
1718	{
1719	LogIO os(LogOrigin("STGrid","getConvFunc",WHERE));
1720	Vector<Float> convFunc;
1721	vector<float> out;
1722
1723	if (cellx_ <= 0.0 \|\| celly_ <= 0.0) {
1724	selectData();
1725	setupGrid();
1726	}
1727
1728	if (convType_ == "BOX" \|\| convType_ == "SF") {
1729	setConvFunc(convFunc);
1730	}
1731	else if (convType_ == "GAUSS") {
1732	Quantum<Double> q1,q2;
1733	readQuantity(q1,gwidth_);
1734	readQuantity(q2,truncate_);
1735	// if (celly_ <= 0.0
1736	// && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") \|\|
1737	// (!q2.getUnit().empty()&&q2.getUnit()!="pixel"))) {
1738	// throw AipsError("You have to call defineImage to get correct convFunc");
1739	// }
1740	setConvFunc(convFunc);
1741	}
1742	else if (convType_ == "GJINC") {
1743	Quantum<Double> q1,q2,q3;
1744	readQuantity(q1,gwidth_);
1745	readQuantity(q2,truncate_);
1746	readQuantity(q3,jwidth_);
1747	// if (celly_ <= 0.0
1748	// && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") \|\|
1749	// (!q2.getUnit().empty()&&q2.getUnit()!="pixel") \|\|
1750	// (!q3.getUnit().empty()&&q3.getUnit()!="pixel"))) {
1751	// throw AipsError("You have to call defineImage to get correct convFunc");
1752	// }
1753	setConvFunc(convFunc);
1754	}
1755	else if (convType_ == "PB") {
1756	throw AipsError("Grid function PB is not available");
1757	}
1758	else {
1759	throw AipsError("Unknown grid function: "+convType_);
1760	}
1761
1762	convFunc.tovector(out);
1763	return out;
1764	}
1765
1766	void STGrid::setConvFunc( Vector<Float> &convFunc )
1767	{
1768	LogIO os(LogOrigin("STGrid","setConvFunc",WHERE));
1769	convSupport_ = userSupport_ ;
1770	if ( convType_ == "BOX" ) {
1771	if ( convSupport_ < 0 )
1772	convSupport_ = 0 ;
1773	Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1774	convFunc.resize( convSize ) ;
1775	boxFunc( convFunc, convSize ) ;
1776	os << LogIO::DEBUGGING
1777	<< "convType_ = " << convType_ << endl
1778	<< "convSupport_ = " << convSupport_ << LogIO::POST;
1779	}
1780	else if ( convType_ == "SF" ) {
1781	if ( convSupport_ < 0 )
1782	convSupport_ = 3 ;
1783	Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1784	convFunc.resize( convSize ) ;
1785	spheroidalFunc( convFunc ) ;
1786	os << LogIO::DEBUGGING
1787	<< "convType_ = " << convType_ << endl
1788	<< "convSupport_ = " << convSupport_ << LogIO::POST;
1789	}
1790	else if ( convType_ == "GAUSS" ) {
1791	// determine pixel gwidth
1792	// default is HWHM corresponding to b = 1.0 (Mangum et al. 2007)
1793	Double pixelGW;
1794	Quantum<Double> q ;
1795	if (!gwidth_.empty()) {
1796	readQuantity( q, gwidth_ );
1797	if ( q.getUnit().empty() \|\| q.getUnit()=="pixel" ) {
1798	pixelGW = q.getValue();
1799	}
1800	else {
1801	pixelGW = q.getValue("rad")/celly_;
1802	}
1803	}
1804	pixelGW = (pixelGW >= 0.0) ? pixelGW : sqrt(log(2.0));
1805	if (pixelGW < 0.0) {
1806	os << LogIO::SEVERE
1807	<< "Negative width is specified for gaussian" << LogIO::EXCEPTION;
1808	}
1809	// determine truncation radius
1810	// default is 3 * HWHM
1811	Double truncate;
1812	if (!truncate_.empty()) {
1813	readQuantity( q, truncate_ );
1814	if ( q.getUnit().empty() \|\| q.getUnit()=="pixel" ) {
1815	truncate = q.getValue();
1816	}
1817	else {
1818	truncate = q.getValue("rad")/celly_;
1819	}
1820	}
1821	//convSupport_ = (Int)(truncate+0.5);
1822	truncate = (truncate >= 0.0) ? truncate : 3.0 * pixelGW;
1823	convSupport_ = Int(truncate);
1824	convSupport_ += (((truncate-(Double)convSupport_) > 0.0) ? 1 : 0);
1825	Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
1826	convFunc.resize( convSize ) ;
1827	gaussFunc( convFunc, pixelGW, truncate ) ;
1828	os << LogIO::DEBUGGING
1829	<< "convType_ = " << convType_ << endl
1830	<< "convSupport_ = " << convSupport_ << endl
1831	<< "truncate_ = " << truncate << "pixel" << endl
1832	<< "gwidth_ = " << pixelGW << "pixel" << LogIO::POST;
1833	}
1834	else if ( convType_ == "GJINC" ) {
1835	// determine pixel gwidth
1836	// default is HWHM corresponding to b = 2.52 (Mangum et al. 2007)
1837	Double pixelGW;
1838	Quantum<Double> q ;
1839	if (!gwidth_.empty()) {
1840	readQuantity( q, gwidth_ );
1841	if ( q.getUnit().empty() \|\| q.getUnit()=="pixel" ) {
1842	pixelGW = q.getValue();
1843	}
1844	else {
1845	pixelGW = q.getValue("rad")/celly_;
1846	}
1847	}
1848	pixelGW = (pixelGW >= 0.0) ? pixelGW : sqrt(log(2.0)) * 2.52;
1849	if (pixelGW < 0.0) {
1850	os << LogIO::SEVERE
1851	<< "Negative width is specified for gaussian" << LogIO::EXCEPTION;
1852	}
1853	// determine pixel c
1854	// default is c = 1.55 (Mangum et al. 2007)
1855	Double pixelJW;
1856	if (!jwidth_.empty()) {
1857	readQuantity( q, jwidth_ );
1858	if ( q.getUnit().empty() \|\| q.getUnit()=="pixel" ) {
1859	pixelJW = q.getValue();
1860	}
1861	else {
1862	pixelJW = q.getValue("rad")/celly_;
1863	}
1864	}
1865	pixelJW = (pixelJW >= 0.0) ? pixelJW : 1.55;
1866	if (pixelJW < 0.0) {
1867	os << LogIO::SEVERE
1868	<< "Negative width is specified for jinc" << LogIO::EXCEPTION;
1869	}
1870	// determine truncation radius
1871	// default is -1.0 (truncate at first null)
1872	Double truncate = -1.0;
1873	if (!truncate_.empty()) {
1874	readQuantity( q, truncate_ );
1875	if ( q.getUnit().empty() \|\| q.getUnit()=="pixel" ) {
1876	truncate = q.getValue();
1877	}
1878	else {
1879	truncate = q.getValue("rad")/celly_;
1880	}
1881	}
1882	//convSupport_ = (truncate >= 0.0) ? (Int)(truncate+0.5) : (Int)(2*pixelJW+0.5);
1883	Double convSupportF = (truncate >= 0.0) ? truncate : (2*pixelJW);
1884	convSupport_ = (Int)convSupportF;
1885	convSupport_ += (((convSupportF-(Double)convSupport_) > 0.0) ? 1 : 0);
1886	Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
1887	convFunc.resize( convSize ) ;
1888	gjincFunc( convFunc, pixelGW, pixelJW, truncate ) ;
1889	os << LogIO::DEBUGGING
1890	<< "convType_ = " << convType_ << endl
1891	<< "convSupport_ = " << convSupport_ << endl
1892	<< "truncate_ = " << truncate << "pixel" << endl
1893	<< "gwidth_ = " << pixelGW << "pixel" << endl
1894	<< "jwidth_ = " << pixelJW << "pixel" << LogIO::POST;
1895	}
1896	else if ( convType_ == "PB" ) {
1897	if ( convSupport_ < 0 )
1898	convSupport_ = 0 ;
1899	pbFunc( convFunc ) ;
1900	}
1901	else {
1902	throw AipsError( "Unsupported convolution function" ) ;
1903	}
1904	}
1905
1906	string STGrid::saveData( string outfile )
1907	{
1908	LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
1909	double t0, t1 ;
1910	t0 = mathutil::gettimeofday_sec() ;
1911
1912	//Int polno = 0 ;
1913	String outfile_ ;
1914	if ( outfile.size() == 0 ) {
1915	if ( infileList_[0].lastchar() == '/' ) {
1916	outfile_ = infileList_[0].substr( 0, infileList_[0].size()-1 ) ;
1917	}
1918	else {
1919	outfile_ = infileList_[0] ;
1920	}
1921	outfile_ += ".grid" ;
1922	}
1923	else {
1924	outfile_ = outfile ;
1925	}
1926	Table tab ;
1927	prepareTable( tab, outfile_ ) ;
1928	fillTable( tab ) ;
1929
1930	t1 = mathutil::gettimeofday_sec() ;
1931	os << LogIO::DEBUGGING << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1932
1933	return outfile_ ;
1934	}
1935
1936	void STGrid::prepareTable( Table &tab, String &name )
1937	{
1938	Table t( infileList_[0], Table::Old ) ;
1939	t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
1940	tab = Table( name, Table::Update ) ;
1941	// 2012/02/13 TN
1942	// explicitly copy subtables since no rows including subtables are
1943	// copied by Table::deepCopy with noRows=True
1944	//TableCopy::copySubTables( tab, t ) ;
1945	const TableRecord &inrec = t.keywordSet();
1946	TableRecord &outrec = tab.rwKeywordSet();
1947	for (uInt i = 0 ; i < inrec.nfields() ; i++) {
1948	if (inrec.type(i) == TpTable) {
1949	String name = inrec.name(i);
1950	Table intable = inrec.asTable(name);
1951	Table outtable = outrec.asTable(name);
1952	TableCopy::copyRows(outtable, intable);
1953	}
1954	}
1955	}
1956
1957	void STGrid::fillTable( Table &tab )
1958	{
1959	IPosition dshape = data_.shape() ;
1960	Int nrow = nx_ * ny_ * npol_ ;
1961	tab.rwKeywordSet().define( "nPol", npol_ ) ;
1962	tab.addRow( nrow ) ;
1963	Vector<Double> cpix( 2 ) ;
1964	cpix(0) = Double( nx_ - 1 ) * 0.5 ;
1965	cpix(1) = Double( ny_ - 1 ) * 0.5 ;
1966	Vector<Double> dir( 2 ) ;
1967	Vector<Double> pix( 2 );
1968	ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
1969	ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
1970	ArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
1971	ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
1972	ScalarColumn<uInt> scannoCol( tab, "SCANNO" ) ;
1973	Int irow = 0 ;
1974	Vector<Float> sp( nchan_ ) ;
1975	Vector<uChar> flag( nchan_, (uChar)1 ) ;
1976	Vector<uChar> unflag( nchan_, (uChar)0 ) ;
1977	Bool bsp, bdata ;
1978	const Float *data_p = data_.getStorage( bdata ) ;
1979	Float wsp_p, sp_p ;
1980	const Float *wdata_p = data_p ;
1981	long step = nx_ * ny_ * npol_ ;
1982	long offset ;
1983	uInt scanno = 0 ;
1984	for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
1985	pix(1) = (Double)(iy);
1986	for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
1987	pix(0) = (Double)(nx_-1-ix);
1988	dircoord_->toWorld(dir,pix);
1989	//os << "dir[" << ix << "," << iy << "]=" << dir << LogIO::POST;
1990	for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
1991	offset = ix + nx_ * (iy + ipol * ny_) ;
1992	//os << "offset = " << offset << LogIO::POST ;
1993	sp_p = sp.getStorage( bsp ) ;
1994	wsp_p = sp_p ;
1995	wdata_p = data_p + offset ;
1996	for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
1997	wsp_p = wdata_p ;
1998	wsp_p++ ;
1999	wdata_p += step ;
2000	}
2001	sp.putStorage( sp_p, bsp ) ;
2002	spectraCol.put( irow, sp ) ;
2003	if ( allEQ( sp, (Float)0.0 ) ) {
2004	flagtraCol.put( irow, flag ) ;
2005	}
2006	else {
2007	flagtraCol.put( irow, unflag ) ;
2008	}
2009	directionCol.put( irow, dir ) ;
2010	polnoCol.put( irow, pollist_[ipol] ) ;
2011	scannoCol.put( irow, scanno ) ;
2012	irow++ ;
2013	}
2014	scanno++ ;
2015	}
2016	}
2017	data_.freeStorage( data_p, bdata ) ;
2018
2019	fillMainColumns( tab ) ;
2020	}
2021
2022	void STGrid::fillMainColumns( Table &tab )
2023	{
2024	// values for fill
2025	//Table t( infileList_[0], Table::Old ) ;
2026	Table t ;
2027	table( t, 0 ) ;
2028	Table tsel = t( t.col( "IFNO" ) == (uInt)ifno_, 1 ) ;
2029	ROTableRow row( tsel ) ;
2030	row.get( 0 ) ;
2031	const TableRecord &rec = row.record() ;
2032	uInt freqId = rec.asuInt( "FREQ_ID" ) ;
2033	uInt molId = rec.asuInt( "MOLECULE_ID" ) ;
2034	uInt tcalId = rec.asuInt( "TCAL_ID" ) ;
2035	uInt focusId = rec.asuInt( "FOCUS_ID" ) ;
2036	uInt weatherId = rec.asuInt( "WEATHER_ID" ) ;
2037	String srcname = rec.asString( "SRCNAME" ) ;
2038	String fieldname = rec.asString( "FIELDNAME" ) ;
2039	Vector<Float> defaultTsys( 1, 1.0 ) ;
2040	// @todo how to set flagtra for gridded spectra?
2041	Vector<uChar> flagtra = rec.asArrayuChar( "FLAGTRA" ) ;
2042	flagtra = (uChar)0 ;
2043	Float opacity = rec.asFloat( "OPACITY" ) ;
2044	Double srcvel = rec.asDouble( "SRCVELOCITY" ) ;
2045	Vector<Double> srcpm = rec.asArrayDouble( "SRCPROPERMOTION" ) ;
2046	Vector<Double> srcdir = rec.asArrayDouble( "SRCDIRECTION" ) ;
2047	Vector<Double> scanrate = rec.asArrayDouble( "SCANRATE" ) ;
2048	Double time = rec.asDouble( "TIME" ) ;
2049	Double interval = rec.asDouble( "INTERVAL" ) ;
2050
2051	// fill columns
2052	Int nrow = tab.nrow() ;
2053	ScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
2054	ScalarColumn<uInt> beamnoCol(tab, "BEAMNO");
2055	ScalarColumn<uInt> freqIdCol( tab, "FREQ_ID" ) ;
2056	ScalarColumn<uInt> molIdCol( tab, "MOLECULE_ID" ) ;
2057	ScalarColumn<uInt> tcalidCol( tab, "TCAL_ID" ) ;
2058	ScalarColumn<Int> fitidCol( tab, "FIT_ID" ) ;
2059	ScalarColumn<uInt> focusidCol( tab, "FOCUS_ID" ) ;
2060	ScalarColumn<uInt> weatheridCol( tab, "WEATHER_ID" ) ;
2061	ArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
2062	ScalarColumn<uInt> rflagCol( tab, "FLAGROW" ) ;
2063	ArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
2064	ScalarColumn<String> srcnameCol( tab, "SRCNAME" ) ;
2065	ScalarColumn<String> fieldnameCol( tab, "FIELDNAME" ) ;
2066	ScalarColumn<Int> srctypeCol( tab, "SRCTYPE" ) ;
2067	ScalarColumn<Float> opacityCol( tab, "OPACITY" ) ;
2068	ScalarColumn<Double> srcvelCol( tab, "SRCVELOCITY" ) ;
2069	ArrayColumn<Double> srcpmCol( tab, "SRCPROPERMOTION" ) ;
2070	ArrayColumn<Double> srcdirCol( tab, "SRCDIRECTION" ) ;
2071	ArrayColumn<Double> scanrateCol( tab, "SCANRATE" ) ;
2072	ScalarColumn<Double> timeCol( tab, "TIME" ) ;
2073	ScalarColumn<Double> intervalCol( tab, "INTERVAL" ) ;
2074	for ( Int i = 0 ; i < nrow ; i++ ) {
2075	ifnoCol.put( i, (uInt)ifno_ ) ;
2076	beamnoCol.put(i, 0);
2077	freqIdCol.put( i, freqId ) ;
2078	molIdCol.put( i, molId ) ;
2079	tcalidCol.put( i, tcalId ) ;
2080	fitidCol.put( i, -1 ) ;
2081	focusidCol.put( i, focusId ) ;
2082	weatheridCol.put( i, weatherId ) ;
2083	//flagtraCol.put( i, flagtra ) ;
2084	rflagCol.put( i, 0 ) ;
2085	tsysCol.put( i, defaultTsys ) ;
2086	srcnameCol.put( i, srcname ) ;
2087	fieldnameCol.put( i, fieldname ) ;
2088	srctypeCol.put( i, (Int)SrcType::PSON ) ;
2089	opacityCol.put( i, opacity ) ;
2090	srcvelCol.put( i, srcvel ) ;
2091	srcpmCol.put( i, srcpm ) ;
2092	srcdirCol.put( i, srcdir ) ;
2093	scanrateCol.put( i, scanrate ) ;
2094	timeCol.put( i, time ) ;
2095	intervalCol.put( i, interval ) ;
2096	}
2097	}
2098
2099	vector<int> STGrid::getResultantMapSize()
2100	{
2101	vector<int> r(2);
2102	r[0] = nx_;
2103	r[1] = ny_;
2104	return r;
2105	}
2106
2107	vector<double> STGrid::getResultantCellSize()
2108	{
2109	vector<double> r(2);
2110	r[0] = cellx_;
2111	r[1] = celly_;
2112	return r;
2113	}
2114
2115	// STGrid2
2116	STGrid2::STGrid2()
2117	: STGrid()
2118	{
2119	}
2120
2121	STGrid2::STGrid2( const ScantableWrapper &s )
2122	: STGrid()
2123	{
2124	setScantable( s ) ;
2125	}
2126
2127	STGrid2::STGrid2( const vector<ScantableWrapper> &v )
2128	: STGrid()
2129	{
2130	setScantableList( v ) ;
2131	}
2132
2133	void STGrid2::setScantable( const ScantableWrapper &s )
2134	{
2135	nfile_ = 1 ;
2136	dataList_.resize( nfile_ ) ;
2137	dataList_[0] = s ;
2138	infileList_.resize( nfile_ ) ;
2139	infileList_[0] = s.getCP()->table().tableName() ;
2140	}
2141
2142	void STGrid2::setScantableList( const vector<ScantableWrapper> &v )
2143	{
2144	nfile_ = v.size() ;
2145	dataList_.resize( nfile_ ) ;
2146	infileList_.resize( nfile_ ) ;
2147	for ( uInt i = 0 ; i < nfile_ ; i++ ) {
2148	dataList_[i] = v[i] ;
2149	infileList_[i] = v[i].getCP()->table().tableName() ;
2150	}
2151	}
2152
2153	ScantableWrapper STGrid2::getResultAsScantable( int tp )
2154	{
2155	ScantableWrapper sw( tp ) ;
2156	CountedPtr<Scantable> s = sw.getCP() ;
2157	s->setHeader( dataList_[0].getCP()->getHeader() ) ;
2158	Table tout, tin ;
2159	String subt[] = { "FREQUENCIES", "FOCUS", "WEATHER",
2160	"TCAL", "MOLECULES", "HISTORY", "FIT" } ;
2161	for ( uInt i = 0 ; i < 7 ; i++ ) {
2162	tout = s->table().rwKeywordSet().asTable(subt[i]) ;
2163	tin = dataList_[0].getCP()->table().rwKeywordSet().asTable(subt[i]) ;
2164	TableCopy::copyRows( tout, tin ) ;
2165	tout.rwKeywordSet() = tin.rwKeywordSet();
2166	}
2167	fillTable( s->table() ) ;
2168	return sw ;
2169	}
2170
2171	void STGrid2::table( Table &tab, uInt i )
2172	{
2173	if ( i < nfile_ )
2174	tab = dataList_[i].getCP()->table() ;
2175	}
2176
2177	}

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

Context Navigation

source: branches/plotter2/src/STGrid.cpp@ 3037

Download in other formats: