Context Navigation

STGrid.cpp@ 2383

Visit:

Last change on this file since 2383 was 2383, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: Yes CAS-2816

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...

Re-defined row-by-row gridding such that each polarization component
is processed separately. After this change, sorting table rows is
no more necessary.

File size: 35.3 KB

Line
1	#include <iostream>
2	#include <fstream>
3
4	#include <casa/BasicSL/String.h>
5	#include <casa/Arrays/Vector.h>
6	#include <casa/Arrays/Matrix.h>
7	#include <casa/Arrays/Cube.h>
8	#include <casa/Arrays/ArrayMath.h>
9	#include <casa/Arrays/ArrayIter.h>
10	#include <casa/Quanta/Quantum.h>
11	#include <casa/Quanta/QuantumHolder.h>
12	#include <casa/Utilities/CountedPtr.h>
13	#include <casa/Logging/LogIO.h>
14
15	#include <tables/Tables/Table.h>
16	#include <tables/Tables/TableRecord.h>
17	#include <tables/Tables/ExprNode.h>
18	#include <tables/Tables/ScalarColumn.h>
19	#include <tables/Tables/ArrayColumn.h>
20	#include <tables/Tables/TableIter.h>
21
22	#include <measures/Measures/MDirection.h>
23
24	#include <MathUtils.h>
25
26	#include "STGrid.h"
27
28	using namespace std ;
29	using namespace casa ;
30	using namespace asap ;
31
32	namespace asap {
33
34	// constructor
35	STGrid::STGrid()
36	{
37	init() ;
38	}
39
40	STGrid::STGrid( const string infile )
41	{
42	init() ;
43
44	setFileIn( infile ) ;
45	}
46
47	void STGrid::init()
48	{
49	ifno_ = -1 ;
50	nx_ = -1 ;
51	ny_ = -1 ;
52	npol_ = 0 ;
53	nchan_ = 0 ;
54	nrow_ = 0 ;
55	cellx_ = 0.0 ;
56	celly_ = 0.0 ;
57	center_ = Vector<Double> ( 2, 0.0 ) ;
58	convType_ = "BOX" ;
59	wtype_ = "UNIFORM" ;
60	convSupport_ = -1 ;
61	userSupport_ = -1 ;
62	convSampling_ = 100 ;
63	nprocessed_ = 0 ;
64	nchunk_ = 0 ;
65	}
66
67	void STGrid::setFileIn( const string infile )
68	{
69	String name( infile ) ;
70	if ( infile_.compare( name ) != 0 ) {
71	infile_ = String( infile ) ;
72	tab_ = Table( infile_ ) ;
73	}
74	}
75
76	void STGrid::setPolList( vector<unsigned int> pols )
77	{
78	pollist_.assign( Vector<uInt>( pols ) ) ;
79	cout << "pollist_ = " << pollist_ << endl ;
80	}
81
82	void STGrid::setScanList( vector<unsigned int> scans )
83	{
84	scanlist_.assign( Vector<uInt>( scans ) ) ;
85	cout << "scanlist_ = " << scanlist_ << endl ;
86	}
87
88	void STGrid::setWeight( const string wType )
89	{
90	wtype_ = String( wType ) ;
91	wtype_.upcase() ;
92	cout << "wtype_ = " << wtype_ << endl ;
93	}
94
95	void STGrid::defineImage( int nx,
96	int ny,
97	string scellx,
98	string scelly,
99	string scenter )
100	{
101	ROArrayColumn<Double> dirCol( tab_, "DIRECTION" ) ;
102	Matrix<Double> direction = dirCol.getColumn() ;
103	Double rmax, rmin, dmax, dmin ;
104	minMax( rmin, rmax, direction.row( 0 ) ) ;
105	minMax( dmin, dmax, direction.row( 1 ) ) ;
106
107	Int npx = (Int)nx ;
108	Int npy = (Int)ny ;
109	String cellx( scellx ) ;
110	String celly( scelly ) ;
111	String center( scenter ) ;
112	setupGrid( npx, npy,
113	cellx, celly,
114	rmin, rmax,
115	dmin, dmax,
116	center ) ;
117	}
118
119	void STGrid::setFunc( string convType,
120	int convSupport )
121	{
122	convType_ = String( convType ) ;
123	convType_.upcase() ;
124	userSupport_ = (Int)convSupport ;
125	}
126
127	#define NEED_UNDERSCORES
128	#if defined(NEED_UNDERSCORES)
129	#define ggridsd ggridsd_
130	#endif
131	extern "C" {
132	void ggridsd(Double*,
133	const Complex*,
134	Int*,
135	Int*,
136	Int*,
137	const Int*,
138	const Int*,
139	const Float*,
140	Int*,
141	Int*,
142	Complex*,
143	Float*,
144	Int*,
145	Int*,
146	Int *,
147	Int *,
148	Int*,
149	Int*,
150	Float*,
151	Int*,
152	Int*,
153	Double*);
154	}
155	void STGrid::call_ggridsd( Double* xy,
156	const Complex* values,
157	Int* nvispol,
158	Int* nvischan,
159	const Int* flag,
160	const Int* rflag,
161	const Float* weight,
162	Int* nrow,
163	Int* irow,
164	Complex* grid,
165	Float* wgrid,
166	Int* nx,
167	Int* ny,
168	Int * npol,
169	Int * nchan,
170	Int* support,
171	Int* sampling,
172	Float* convFunc,
173	Int *chanMap,
174	Int *polMap )
175	{
176	// parameters for gridding
177	Int idopsf = 0 ;
178	Int len = (npol)(*nchan) ;
179	Double *sumw_p = new Double[len] ;
180	{
181	Double *work_p = sumw_p ;
182	for ( Int i = 0 ; i < len ; i++ ) {
183	*work_p = 0.0 ;
184	work_p++ ;
185	}
186	}
187
188	// call ggridsd
189	ggridsd( xy,
190	values,
191	nvispol,
192	nvischan,
193	&idopsf,
194	flag,
195	rflag,
196	weight,
197	nrow,
198	irow,
199	grid,
200	wgrid,
201	nx,
202	ny,
203	npol,
204	nchan,
205	support,
206	sampling,
207	convFunc,
208	chanMap,
209	polMap,
210	sumw_p ) ;
211
212	// finalization
213	delete sumw_p ;
214	}
215
216	Bool STGrid::pastEnd()
217	{
218	LogIO os( LogOrigin("STGrid","pastEnd",WHERE) ) ;
219	Bool b = nprocessed_ >= nrow_ ;
220	return b ;
221	}
222
223	void STGrid::grid()
224	{
225	// data selection
226	selectData() ;
227	setupArray() ;
228
229	if ( npol_ > 1 ) {
230	LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
231	os << LogIO::WARN << "STGrid doesn't support assigning polarization-dependent weight. Use averaged weight over polarization." << LogIO::POST ;
232	}
233
234	if ( wtype_.compare("UNIFORM") != 0 &&
235	wtype_.compare("TINT") != 0 &&
236	wtype_.compare("TSYS") != 0 &&
237	wtype_.compare("TINTSYS") != 0 ) {
238	LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
239	os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
240	wtype_ = "UNIFORM" ;
241	}
242
243	Bool doAll = examine() ;
244
245	if ( doAll )
246	gridPerPol() ;
247	else
248	gridPerRow() ;
249	}
250
251	Bool STGrid::examine()
252	{
253	// TODO: nchunk_ must be determined from nchan_, npol_, and (nx_,ny_)
254	// by considering data size to be allocated for ggridsd input/output
255	nchunk_ = 100 ;
256	Bool b = nchunk_ >= nrow_ ;
257	nchunk_ = min( nchunk_, nrow_ ) ;
258	return b ;
259	}
260
261	void STGrid::gridPerRow()
262	{
263	LogIO os( LogOrigin("STGrid", "gridPerRow", WHERE) ) ;
264	double t0, t1 ;
265
266	// grid parameter
267	os << LogIO::DEBUGGING ;
268	os << "----------" << endl ;
269	os << "Grid parameter summary" << endl ;
270	os << " (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
271	os << " (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
272	os << " center = " << center_ << endl ;
273	os << "----------" << LogIO::POST ;
274	os << LogIO::NORMAL ;
275
276	// boolean for getStorage
277	Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
278
279	// convolution kernel
280	Vector<Float> convFunc ;
281	t0 = mathutil::gettimeofday_sec() ;
282	setConvFunc( convFunc ) ;
283	t1 = mathutil::gettimeofday_sec() ;
284	os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
285	//cout << "convSupport=" << convSupport_ << endl ;
286	//cout << "convFunc=" << convFunc << endl ;
287
288	// grid data
289	Int gnx = nx_ ;
290	Int gny = ny_ ;
291	// Extend grid plane with convSupport_
292	// Int gnx = nx_+convSupport_*2 ;
293	// Int gny = ny_+convSupport_*2 ;
294	IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
295	Array<Complex> gdataArrC( gshape, 0.0 ) ;
296	// 2011/12/20 TN
297	// data_ and gwgtArr share storage
298	data_.resize( gshape ) ;
299	data_ = 0.0 ;
300	Array<Float> gwgtArr( data_ ) ;
301
302	// data storage
303	Int irow = -1 ;
304	IPosition cshape( 3, npol_, nchan_, nchunk_ ) ;
305	IPosition mshape( 2, npol_, nchunk_ ) ;
306	IPosition vshape( 1, nchunk_ ) ;
307	IPosition dshape( 2, 2, nchunk_ ) ;
308	IPosition wshape( 2, nchan_, nchunk_ ) ;
309	Array<Complex> spectra( wshape ) ;
310	Array<Int> flagtra( wshape ) ;
311	Array<Int> rflag( vshape ) ;
312	Array<Float> weight( wshape ) ;
313	Array<Double> direction( dshape ) ;
314	Array<Double> xypos( dshape ) ;
315
316	spectraF_.resize( wshape ) ;
317	flagtraUC_.resize( wshape ) ;
318	rflagUI_.resize( vshape ) ;
319
320	// parameters for gridding
321	Int *chanMap = new Int[nchan_] ;
322	{
323	Int *work_p = chanMap ;
324	for ( Int i = 0 ; i < nchan_ ; i++ ) {
325	*work_p = i ;
326	work_p++ ;
327	}
328	}
329	Int *polMap = new Int[1] ;
330	Int nvispol = 1 ;
331
332	// for performance check
333	double eGetDataChunk = 0.0 ;
334	double eToPixel = 0.0 ;
335	double eGGridSD = 0.0 ;
336	double eToInt = 0.0 ;
337
338	for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
339	initPol( ipol ) ;
340
341	polMap[0] = ipol ;
342
343	while( !pastEnd() ) {
344	// retrieve data
345	t0 = mathutil::gettimeofday_sec() ;
346	Int nrow = getDataChunk( spectra, direction, flagtra, rflag, weight ) ;
347	//os << "nrow = " << nrow << LogIO::POST ;
348	t1 = mathutil::gettimeofday_sec() ;
349	eGetDataChunk += t1-t0 ;
350	eToInt += subtime_ ;
351
352	// world -> pixel
353	t0 = mathutil::gettimeofday_sec() ;
354	toPixel( direction, xypos ) ;
355	t1 = mathutil::gettimeofday_sec() ;
356	eToPixel += t1-t0 ;
357
358	// prepare pointer
359	Float *conv_p = convFunc.getStorage( deleteConv ) ;
360	Complex *gdata_p = gdataArrC.getStorage( deleteDataG ) ;
361	Float *wdata_p = gwgtArr.getStorage( deleteWgtG ) ;
362	const Complex *data_p = spectra.getStorage( deleteData ) ;
363	const Float *wgt_p = weight.getStorage( deleteWgt ) ;
364	const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
365	const Int *rflag_p = rflag.getStorage( deleteFlagR ) ;
366	Double *xypos_p = xypos.getStorage( deletePos ) ;
367
368	// call ggridsd
369	irow = -1 ;
370	t0 = mathutil::gettimeofday_sec() ;
371	call_ggridsd( xypos_p,
372	data_p,
373	//&npol_,
374	&nvispol,
375	&nchan_,
376	flag_p,
377	rflag_p,
378	wgt_p,
379	&nrow,
380	&irow,
381	gdata_p,
382	wdata_p,
383	&gnx,
384	&gny,
385	&npol_,
386	&nchan_,
387	&convSupport_,
388	&convSampling_,
389	conv_p,
390	chanMap,
391	polMap ) ;
392	t1 = mathutil::gettimeofday_sec() ;
393	eGGridSD += t1-t0 ;
394
395	// finalization
396	convFunc.putStorage( conv_p, deleteConv ) ;
397	gdataArrC.putStorage( gdata_p, deleteDataG ) ;
398	gwgtArr.putStorage( wdata_p, deleteWgtG ) ;
399	xypos.putStorage( xypos_p, deletePos ) ;
400	spectra.freeStorage( data_p, deleteData ) ;
401	weight.freeStorage( wgt_p, deleteWgt ) ;
402	flagtra.freeStorage( flag_p, deleteFlag ) ;
403	rflag.freeStorage( rflag_p, deleteFlagR ) ;
404
405	}
406
407	nprocessed_ = 0 ;
408	}
409	os << "getDataChunk: elapsed time is " << eGetDataChunk << " sec." << LogIO::POST ;
410	os << "toPixel: elapsed time is " << eToPixel << " sec." << LogIO::POST ;
411	os << "ggridsd: elapsed time is " << eGGridSD << " sec." << LogIO::POST ;
412	os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
413
414	delete polMap ;
415	delete chanMap ;
416
417	// set data
418	setData( gdataArrC, gwgtArr ) ;
419
420	}
421
422	void STGrid::gridPerPol()
423	{
424	LogIO os( LogOrigin("STGrid", "gridPerPol", WHERE) ) ;
425	double t0, t1 ;
426
427	// grid parameter
428	os << LogIO::DEBUGGING ;
429	os << "----------" << endl ;
430	os << "Grid parameter summary" << endl ;
431	os << " (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
432	os << " (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
433	os << " center = " << center_ << endl ;
434	os << "----------" << LogIO::POST ;
435	os << LogIO::NORMAL ;
436
437	// convolution kernel
438	Vector<Float> convFunc ;
439	t0 = mathutil::gettimeofday_sec() ;
440	setConvFunc( convFunc ) ;
441	t1 = mathutil::gettimeofday_sec() ;
442	os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
443	//cout << "convSupport=" << convSupport_ << endl ;
444	//cout << "convFunc=" << convFunc << endl ;
445
446	// prepare grid data storage
447	Int gnx = nx_ ;
448	Int gny = ny_ ;
449	// // Extend grid plane with convSupport_
450	// Int gnx = nx_+convSupport_*2 ;
451	// Int gny = ny_+convSupport_*2 ;
452	IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
453	Array<Complex> gdataArrC( gshape, 0.0 ) ;
454	//Array<Float> gwgtArr( gshape, 0.0 ) ;
455	// 2011/12/20 TN
456	// data_ and weight array shares storage
457	data_.resize( gshape ) ;
458	data_ = 0.0 ;
459	Array<Float> gwgtArr( data_ ) ;
460
461	// to read data from the table
462	IPosition mshape( 2, nchan_, nrow_ ) ;
463	IPosition dshape( 2, 2, nrow_ ) ;
464	Array<Complex> spectra( mshape ) ;
465	Array<Double> direction( dshape ) ;
466	Array<Int> flagtra( mshape ) ;
467	Array<Int> rflag( IPosition(1,nrow_) ) ;
468	Array<Float> weight( mshape ) ;
469	Array<Double> xypos( dshape ) ;
470
471	// maps
472	Int *chanMap = new Int[nchan_] ;
473	{
474	Int *work_p = chanMap ;
475	for ( Int i = 0 ; i < nchan_ ; i++ ) {
476	*work_p = i ;
477	work_p++ ;
478	}
479	}
480	Int *polMap = new Int[1] ;
481	Int nvispol = 1 ;
482
483	// for performance check
484	double eInitPol = 0.0 ;
485	double eGetData = 0.0 ;
486	double eToPixel = 0.0 ;
487	double eGGridSD = 0.0 ;
488	double eToInt = 0.0 ;
489
490	for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
491	t0 = mathutil::gettimeofday_sec() ;
492	initPol( ipol ) ;
493	t1 = mathutil::gettimeofday_sec() ;
494	eInitPol += t1-t0 ;
495
496	// retrieve data
497	t0 = mathutil::gettimeofday_sec() ;
498	getDataPerPol( spectra, direction, flagtra, rflag, weight ) ;
499	t1 = mathutil::gettimeofday_sec() ;
500	eGetData += t1-t0 ;
501
502	IPosition sshape = spectra.shape() ;
503
504	// world -> pixel
505	t0 = mathutil::gettimeofday_sec() ;
506	toPixel( direction, xypos ) ;
507	t1 = mathutil::gettimeofday_sec() ;
508	//os << "xypos=" << xypos << LogIO::POST ;
509	os << "toPixel: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
510
511	// call ggridsd
512	Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
513	Double *xypos_p = xypos.getStorage( deletePos ) ;
514	const Complex *data_p = spectra.getStorage( deleteData ) ;
515	const Float *wgt_p = weight.getStorage( deleteWgt ) ;
516	const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
517	const Int *rflag_p = rflag.getStorage( deleteFlagR ) ;
518	Float *conv_p = convFunc.getStorage( deleteConv ) ;
519	Complex *gdata_p = gdataArrC.getStorage( deleteDataG ) ;
520	Float *wdata_p = gwgtArr.getStorage( deleteWgtG ) ;
521	polMap[0] = ipol ;
522	t0 = mathutil::gettimeofday_sec() ;
523	Int irow = -1 ;
524	call_ggridsd( xypos_p,
525	data_p,
526	&nvispol,
527	&nchan_,
528	flag_p,
529	rflag_p,
530	wgt_p,
531	&nrow_,
532	&irow,
533	gdata_p,
534	wdata_p,
535	&gnx,
536	&gny,
537	&npol_,
538	&nchan_,
539	&convSupport_,
540	&convSampling_,
541	conv_p,
542	chanMap,
543	polMap ) ;
544	t1 = mathutil::gettimeofday_sec() ;
545	os << "ggridsd: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
546	xypos.putStorage( xypos_p, deletePos ) ;
547	spectra.freeStorage( data_p, deleteData ) ;
548	weight.freeStorage( wgt_p, deleteWgt ) ;
549	flagtra.freeStorage( flag_p, deleteFlag ) ;
550	rflag.freeStorage( rflag_p, deleteFlagR ) ;
551	convFunc.putStorage( conv_p, deleteConv ) ;
552	gdataArrC.putStorage( gdata_p, deleteDataG ) ;
553	gwgtArr.putStorage( wdata_p, deleteWgtG ) ;
554	}
555
556	os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
557	os << "getData: elapsed time is " << eGetData << " sec." << LogIO::POST ;
558	// delete maps
559	delete polMap ;
560	delete chanMap ;
561
562	setData( gdataArrC, gwgtArr ) ;
563	// os << "gdataArr = " << gdataArr << LogIO::POST ;
564	// os << "gwgtArr = " << gwgtArr << LogIO::POST ;
565	// os << "data_ " << data_ << LogIO::POST ;
566	}
567
568	void STGrid::initPol( Int ipol )
569	{
570	ptab_ = tab_( tab_.col("POLNO") == (uInt)ipol ) ;
571
572	attach( ptab_ ) ;
573	}
574
575	void STGrid::setData( Array<Complex> &gdata,
576	Array<Float> &gwgt )
577	{
578	// 2011/12/20 TN
579	// gwgt and data_ share storage
580	LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
581	double t0, t1 ;
582	t0 = mathutil::gettimeofday_sec() ;
583	uInt len = data_.nelements() ;
584	const Complex *w1_p ;
585	Float *w2_p ;
586	Bool b1, b2 ;
587	const Complex *gdata_p = gdata.getStorage( b1 ) ;
588	Float *gwgt_p = data_.getStorage( b2 ) ;
589	w1_p = gdata_p ;
590	w2_p = gwgt_p ;
591	for ( uInt i = 0 ; i < len ; i++ ) {
592	if ( w2_p > 0.0 ) w2_p = (w1_p).real() / w2_p ;
593	w1_p++ ;
594	w2_p++ ;
595	}
596	gdata.freeStorage( gdata_p, b1 ) ;
597	data_.putStorage( gwgt_p, b2 ) ;
598	t1 = mathutil::gettimeofday_sec() ;
599	os << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
600	}
601
602	void STGrid::setupGrid( Int &nx,
603	Int &ny,
604	String &cellx,
605	String &celly,
606	Double &xmin,
607	Double &xmax,
608	Double &ymin,
609	Double &ymax,
610	String &center )
611	{
612	LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
613	//cout << "nx=" << nx << ", ny=" << ny << endl ;
614
615	// center position
616	if ( center.size() == 0 ) {
617	center_(0) = 0.5 * ( xmin + xmax ) ;
618	center_(1) = 0.5 * ( ymin + ymax ) ;
619	}
620	else {
621	String::size_type pos0 = center.find( " " ) ;
622	if ( pos0 == String::npos ) {
623	throw AipsError( "bad string format in parameter center" ) ;
624	}
625	String::size_type pos1 = center.find( " ", pos0+1 ) ;
626	String typestr, xstr, ystr ;
627	if ( pos1 != String::npos ) {
628	typestr = center.substr( 0, pos0 ) ;
629	xstr = center.substr( pos0+1, pos1-pos0 ) ;
630	ystr = center.substr( pos1+1 ) ;
631	// todo: convert to J2000 (or direction ref for DIRECTION column)
632	}
633	else {
634	typestr = "J2000" ;
635	xstr = center.substr( 0, pos0 ) ;
636	ystr = center.substr( pos0+1 ) ;
637	}
638	QuantumHolder qh ;
639	String err ;
640	qh.fromString( err, xstr ) ;
641	Quantum<Double> xcen = qh.asQuantumDouble() ;
642	qh.fromString( err, ystr ) ;
643	Quantum<Double> ycen = qh.asQuantumDouble() ;
644	center_(0) = xcen.getValue( "rad" ) ;
645	center_(1) = ycen.getValue( "rad" ) ;
646	}
647
648
649	nx_ = nx ;
650	ny_ = ny ;
651	if ( nx < 0 && ny > 0 ) {
652	nx_ = ny ;
653	ny_ = ny ;
654	}
655	if ( ny < 0 && nx > 0 ) {
656	nx_ = nx ;
657	ny_ = nx ;
658	}
659
660	//Double wx = xmax - xmin ;
661	//Double wy = ymax - ymin ;
662	Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
663	Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
664	// take 10% margin
665	wx *= 1.10 ;
666	wy *= 1.10 ;
667
668	Quantum<Double> qcellx ;
669	Quantum<Double> qcelly ;
670	//cout << "nx_ = " << nx_ << ", ny_ = " << ny_ << endl ;
671	if ( cellx.size() != 0 && celly.size() != 0 ) {
672	readQuantity( qcellx, cellx ) ;
673	readQuantity( qcelly, celly ) ;
674	}
675	else if ( celly.size() != 0 ) {
676	os << "Using celly to x-axis..." << LogIO::POST ;
677	readQuantity( qcelly, celly ) ;
678	qcellx = qcelly ;
679	}
680	else if ( cellx.size() != 0 ) {
681	os << "Using cellx to y-axis..." << LogIO::POST ;
682	readQuantity( qcellx, cellx ) ;
683	qcelly = qcellx ;
684	}
685	else {
686	if ( nx_ < 0 ) {
687	os << "No user preference in grid setting. Using default..." << LogIO::POST ;
688	readQuantity( qcellx, "1.0arcmin" ) ;
689	qcelly = qcellx ;
690	}
691	else {
692	if ( wx == 0.0 ) {
693	os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
694	wx = 0.00290888 ;
695	}
696	if ( wy == 0.0 ) {
697	os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
698	wy = 0.00290888 ;
699	}
700	qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
701	qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
702	}
703	}
704	cellx_ = qcellx.getValue( "rad" ) ;
705	celly_ = qcelly.getValue( "rad" ) ;
706	if ( nx_ < 0 ) {
707	if ( wx == 0.0 ) {
708	os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
709	wx = 0.00290888 ;
710	}
711	if ( wy == 0.0 ) {
712	os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
713	wy = 0.00290888 ;
714	}
715	nx_ = Int( ceil( wx/cellx_ ) ) ;
716	ny_ = Int( ceil( wy/celly_ ) ) ;
717	}
718	}
719
720	void STGrid::selectData()
721	{
722	Int ifno = ifno_ ;
723	Table taborg( infile_ ) ;
724	if ( ifno == -1 ) {
725	LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
726	// os << LogIO::SEVERE
727	// << "Please set IFNO before actual gridding"
728	// << LogIO::EXCEPTION ;
729	ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
730	ifno = ifnoCol( 0 ) ;
731	os << LogIO::WARN
732	<< "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
733	}
734	// tab = taborg( taborg.col("IFNO") == ifno ) ;
735	TableExprNode node ;
736	node = taborg.col("IFNO") == ifno ;
737	if ( scanlist_.size() > 0 ) {
738	node = node && taborg.col("SCANNO").in( scanlist_ ) ;
739	}
740	tab_ = taborg( node ) ;
741	if ( tab_.nrow() == 0 ) {
742	LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
743	os << LogIO::SEVERE
744	<< "No corresponding rows for given selection: IFNO " << ifno ;
745	if ( scanlist_.size() > 0 )
746	os << " SCANNO " << scanlist_ ;
747	os << LogIO::EXCEPTION ;
748	}
749	//attach( tab_ ) ;
750	}
751
752	void STGrid::attach( Table &tab )
753	{
754	// attach to table
755	spectraCol_.attach( tab, "SPECTRA" ) ;
756	flagtraCol_.attach( tab, "FLAGTRA" ) ;
757	directionCol_.attach( tab, "DIRECTION" ) ;
758	flagRowCol_.attach( tab, "FLAGROW" ) ;
759	tsysCol_.attach( tab, "TSYS" ) ;
760	intervalCol_.attach( tab, "INTERVAL" ) ;
761	// Vector<String> colnames( 6 ) ;
762	// colnames[0] = "SPECTRA" ;
763	// colnames[1] = "FLAGTRA" ;
764	// colnames[2] = "DIRECTION" ;
765	// colnames[3] = "FLAGROW" ;
766	// colnames[4] = "TSYS" ;
767	// colnames[5] = "INTERVAL" ;
768	// row_ = ROTableRow( tab_, colnames ) ;
769	// const TableRecord &rec = row_.record() ;
770	// spectraRF_.attachToRecord( rec, colnames[0] ) ;
771	// flagtraRF_.attachToRecord( rec, colnames[1] ) ;
772	// directionRF_.attachToRecord( rec, colnames[2] ) ;
773	// flagRowRF_.attachToRecord( rec, colnames[3] ) ;
774	// tsysRF_.attachToRecord( rec, colnames[4] ) ;
775	// intervalRF_.attachToRecord( rec, colnames[5] ) ;
776	}
777
778	void STGrid::getDataPerPol( Array<Float> &spectra,
779	Array<Double> &direction,
780	Array<uChar> &flagtra,
781	Array<uInt> &rflag,
782	Array<Float> &weight )
783	{
784	LogIO os( LogOrigin("STGrid","getDataPerPol",WHERE) ) ;
785
786	// 2011/12/22 TN
787	// weight and tsys shares its storage
788	Array<Float> tsys( weight ) ;
789	Array<Double> tint( rflag.shape() ) ;
790
791	Vector<uInt> rflagVec( rflag ) ;
792	Vector<Double> tintVec( tint ) ;
793
794	spectraCol_.getColumn( spectra ) ;
795	flagtraCol_.getColumn( flagtra ) ;
796	flagRowCol_.getColumn( rflagVec ) ;
797	directionCol_.getColumn( direction ) ;
798	intervalCol_.getColumn( tintVec ) ;
799	Vector<Float> tsysTemp = tsysCol_( 0 ) ;
800	if ( tsysTemp.nelements() == (uInt)nchan_ ) {
801	tsysCol_.getColumn( tsys ) ;
802	}
803	else {
804	tsys = tsysTemp[0] ;
805	}
806
807	getWeight( weight, tsys, tint ) ;
808	}
809
810	void STGrid::getDataPerPol( Array<Complex> &spectra,
811	Array<Double> &direction,
812	Array<Int> &flagtra,
813	Array<Int> &rflag,
814	Array<Float> &weight )
815	{
816	LogIO os( LogOrigin("STGrid","getDataPerPol",WHERE) ) ;
817	double t0, t1 ;
818
819	Array<uChar> flagUC( flagtra.shape() ) ;
820	Array<uInt> rflagUI( rflag.shape() ) ;
821	Array<Float> spectraF( spectra.shape() ) ;
822	getDataPerPol( spectraF, direction, flagUC, rflagUI, weight ) ;
823
824	convertArray( spectra, spectraF ) ;
825	t0 = mathutil::gettimeofday_sec() ;
826	toInt( flagUC, flagtra ) ;
827	toInt( rflagUI, rflag ) ;
828	t1 = mathutil::gettimeofday_sec() ;
829	//os << "toInt: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
830	}
831
832	Int STGrid::getDataChunk( Array<Complex> &spectra,
833	Array<Double> &direction,
834	Array<Int> &flagtra,
835	Array<Int> &rflag,
836	Array<Float> &weight )
837	{
838	Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
839	if ( nrow < nchunk_ ) {
840	spectra.resize( spectraF_.shape() ) ;
841	flagtra.resize( flagtraUC_.shape() ) ;
842	rflag.resize( rflagUI_.shape() ) ;
843	}
844	double t0, t1 ;
845	t0 = mathutil::gettimeofday_sec() ;
846	convertArray( spectra, spectraF_ ) ;
847	toInt( flagtraUC_, flagtra ) ;
848	toInt( rflagUI_, rflag ) ;
849	t1 = mathutil::gettimeofday_sec() ;
850	subtime_ = t1 - t0 ;
851
852	return nrow ;
853	}
854
855	Int STGrid::getDataChunk( Array<Float> &spectra,
856	Array<Double> &direction,
857	Array<uChar> &flagtra,
858	Array<uInt> &rflag,
859	Array<Float> &weight )
860	{
861	LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
862	Int nrow = spectra.shape()[1] ;
863	Int remainingRow = nrow_ - nprocessed_ ;
864	if ( remainingRow < nrow ) {
865	nrow = remainingRow ;
866	IPosition mshape( 2, nchan_, nrow ) ;
867	IPosition vshape( 1, nrow ) ;
868	spectra.resize( mshape ) ;
869	flagtra.resize( mshape ) ;
870	direction.resize( IPosition(2,2,nrow) ) ;
871	rflag.resize( vshape ) ;
872	weight.resize( mshape ) ;
873	}
874	// 2011/12/22 TN
875	// tsys shares its storage with weight
876	Array<Float> tsys( weight ) ;
877	Array<Double> tint( rflag.shape() ) ;
878
879	Vector<uInt> rflagVec( rflag ) ;
880	Vector<Double> tintVec( tint ) ;
881
882	RefRows rows( nprocessed_, nprocessed_+nrow-1, 1 ) ;
883	os << "rows.nrows()=" << rows.nrows() << LogIO::POST ;
884	spectraCol_.getColumnCells( rows, spectra ) ;
885	flagtraCol_.getColumnCells( rows, flagtra ) ;
886	directionCol_.getColumnCells( rows, direction ) ;
887	flagRowCol_.getColumnCells( rows, rflagVec ) ;
888	intervalCol_.getColumnCells( rows, tintVec ) ;
889	Vector<Float> tsysTemp = tsysCol_( nprocessed_ ) ;
890	if ( tsysTemp.nelements() == (uInt)nchan_ )
891	tsysCol_.getColumnCells( rows, tsys ) ;
892	else
893	tsys = tsysTemp[0] ;
894
895	getWeight( weight, tsys, tint ) ;
896
897	nprocessed_ += nrow ;
898
899	return nrow ;
900	}
901
902	void STGrid::setupArray()
903	{
904	LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
905	ROScalarColumn<uInt> polnoCol( tab_, "POLNO" ) ;
906	Vector<uInt> pols = polnoCol.getColumn() ;
907	//os << pols << LogIO::POST ;
908	Vector<uInt> pollistOrg ;
909	uInt npolOrg = 0 ;
910	uInt polno ;
911	for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
912	//polno = polnoCol( i ) ;
913	polno = pols( i ) ;
914	if ( allNE( pollistOrg, polno ) ) {
915	pollistOrg.resize( npolOrg+1, True ) ;
916	pollistOrg[npolOrg] = polno ;
917	npolOrg++ ;
918	}
919	}
920	if ( pollist_.size() == 0 )
921	pollist_ = pollistOrg ;
922	else {
923	Vector<uInt> newlist ;
924	uInt newsize = 0 ;
925	for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
926	if ( anyEQ( pollistOrg, pollist_[i] ) ) {
927	newlist.resize( newsize+1, True ) ;
928	newlist[newsize] = pollist_[i] ;
929	newsize++ ;
930	}
931	}
932	pollist_.assign( newlist ) ;
933	}
934	npol_ = pollist_.size() ;
935	if ( npol_ == 0 ) {
936	os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
937	}
938	nrow_ = tab_.nrow() / npolOrg ;
939	ROArrayColumn<uChar> tmpCol( tab_, "FLAGTRA" ) ;
940	nchan_ = tmpCol( 0 ).nelements() ;
941	// os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
942	// << "nchan_ = " << nchan_ << endl
943	// << "nrow_ = " << nrow_ << LogIO::POST ;
944	}
945
946	void STGrid::getWeight( Array<Float> &w,
947	Array<Float> &tsys,
948	Array<Double> &tint )
949	{
950	LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
951	//os << "start getWeight" << LogIO::POST ;
952	double t0, t1 ;
953	t0 = mathutil::gettimeofday_sec() ;
954	// 2011/12/22 TN
955	// w (weight) and tsys share storage
956	IPosition refShape = tsys.shape() ;
957	Int nchan = refShape[0] ;
958	Int nrow = refShape[1] ;
959	// os << "nchan=" << nchan << ", nrow=" << nrow << LogIO::POST ;
960	// os << "w.shape()=" << w.shape() << endl
961	// << "tsys.shape()=" << tsys.shape() << endl
962	// << "tint.shape()=" << tint.shape() << LogIO::POST ;
963
964	// set weight
965	if ( wtype_.compare( "UNIFORM" ) == 0 ) {
966	w = 1.0 ;
967	}
968	else if ( wtype_.compare( "TINT" ) == 0 ) {
969	Bool b0, b1 ;
970	Float *w_p = w.getStorage( b0 ) ;
971	Float *w0_p = w_p ;
972	const Double *ti_p = tint.getStorage( b1 ) ;
973	const Double *w1_p = ti_p ;
974	for ( Int irow = 0 ; irow < nrow ; irow++ ) {
975	for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
976	w0_p = w1_p ;
977	w0_p++ ;
978	}
979	w1_p++ ;
980	}
981	w.putStorage( w_p, b0 ) ;
982	tint.freeStorage( ti_p, b1 ) ;
983	}
984	else if ( wtype_.compare( "TSYS" ) == 0 ) {
985	Bool b0 ;
986	Float *w_p = w.getStorage( b0 ) ;
987	Float *w0_p = w_p ;
988	for ( Int irow = 0 ; irow < nrow ; irow++ ) {
989	for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
990	Float temp = *w0_p ;
991	w0_p = 1.0 / ( temp temp ) ;
992	w0_p++ ;
993	}
994	}
995	w.putStorage( w_p, b0 ) ;
996	}
997	else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
998	Bool b0, b1 ;
999	Float *w_p = w.getStorage( b0 ) ;
1000	Float *w0_p = w_p ;
1001	const Double *ti_p = tint.getStorage( b1 ) ;
1002	const Double *w1_p = ti_p ;
1003	for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1004	Float interval = *w1_p ;
1005	for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1006	Float temp = *w0_p ;
1007	w0_p = interval / ( temp temp ) ;
1008	w0_p++ ;
1009	}
1010	w1_p++ ;
1011	}
1012	w.putStorage( w_p, b0 ) ;
1013	tint.freeStorage( ti_p, b1 ) ;
1014	}
1015	else {
1016	//LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
1017	//os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
1018	w = 1.0 ;
1019	}
1020
1021	t1 = mathutil::gettimeofday_sec() ;
1022	//os << "getWeight: elapsed time is " << t1-t0 << " sec" << LogIO::POST ;
1023	subtime_ += t1-t0 ;
1024	}
1025
1026	void STGrid::toInt( Array<uChar> &u, Array<Int> &v )
1027	{
1028	uInt len = u.nelements() ;
1029	Int *int_p = new Int[len] ;
1030	Bool deleteIt ;
1031	const uChar *data_p = u.getStorage( deleteIt ) ;
1032	Int *i_p = int_p ;
1033	const uChar *u_p = data_p ;
1034	for ( uInt i = 0 ; i < len ; i++ ) {
1035	i_p = ( u_p == 0 ) ? 0 : 1 ;
1036	i_p++ ;
1037	u_p++ ;
1038	}
1039	u.freeStorage( data_p, deleteIt ) ;
1040	v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1041	}
1042
1043	void STGrid::toInt( Array<uInt> &u, Array<Int> &v )
1044	{
1045	uInt len = u.nelements() ;
1046	Int *int_p = new Int[len] ;
1047	Bool deleteIt ;
1048	const uInt *data_p = u.getStorage( deleteIt ) ;
1049	Int *i_p = int_p ;
1050	const uInt *u_p = data_p ;
1051	for ( uInt i = 0 ; i < len ; i++ ) {
1052	i_p = ( u_p == 0 ) ? 0 : 1 ;
1053	i_p++ ;
1054	u_p++ ;
1055	}
1056	u.freeStorage( data_p, deleteIt ) ;
1057	v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1058	}
1059
1060	void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
1061	{
1062	// gridding will be done on (nx_+2convSupport_) x (ny_+2convSupport_)
1063	// grid plane to avoid unexpected behavior on grid edge
1064	Block<Double> pixc( 2 ) ;
1065	pixc[0] = Double( nx_-1 ) * 0.5 ;
1066	pixc[1] = Double( ny_-1 ) * 0.5 ;
1067	// pixc[0] = Double( nx_+2convSupport_-1 ) 0.5 ;
1068	// pixc[1] = Double( ny_+2convSupport_-1 ) 0.5 ;
1069	uInt nrow = world.shape()[1] ;
1070	Bool bw, bp ;
1071	const Double *w_p = world.getStorage( bw ) ;
1072	Double *p_p = pixel.getStorage( bp ) ;
1073	const Double *ww_p = w_p ;
1074	Double *wp_p = p_p ;
1075	for ( uInt i = 0 ; i < nrow ; i++ ) {
1076	wp_p = pixc[0] + ( ww_p - center_[0] ) / cellx_ ;
1077	wp_p++ ;
1078	ww_p++ ;
1079	wp_p = pixc[1] + ( ww_p - center_[1] ) / celly_ ;
1080	wp_p++ ;
1081	ww_p++ ;
1082	}
1083	world.freeStorage( w_p, bw ) ;
1084	pixel.putStorage( p_p, bp ) ;
1085	// String gridfile = "grid."+convType_+"."+String::toString(convSupport_)+".dat" ;
1086	// ofstream ofs( gridfile.c_str(), ios::out ) ;
1087	// ofs << "center " << center_(0) << " " << pixc(0)
1088	// << " " << center_(1) << " " << pixc(1) << endl ;
1089	// for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1090	// ofs << irow ;
1091	// for ( uInt i = 0 ; i < 2 ; i++ ) {
1092	// ofs << " " << world(i, irow) << " " << pixel(i, irow) ;
1093	// }
1094	// ofs << endl ;
1095	// }
1096	// ofs.close() ;
1097	}
1098
1099	void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
1100	{
1101	convFunc = 0.0 ;
1102	for ( Int i = 0 ; i < convSize/2 ; i++ )
1103	convFunc(i) = 1.0 ;
1104	}
1105
1106	#define NEED_UNDERSCORES
1107	#if defined(NEED_UNDERSCORES)
1108	#define grdsf grdsf_
1109	#endif
1110	extern "C" {
1111	void grdsf(Double, Double);
1112	}
1113	void STGrid::spheroidalFunc( Vector<Float> &convFunc )
1114	{
1115	convFunc = 0.0 ;
1116	for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
1117	Double nu = Double(i) / Double(convSupport_*convSampling_) ;
1118	Double val ;
1119	grdsf( &nu, &val ) ;
1120	convFunc(i) = ( 1.0 - nu * nu ) * val ;
1121	}
1122	}
1123
1124	void STGrid::gaussFunc( Vector<Float> &convFunc )
1125	{
1126	convFunc = 0.0 ;
1127	// HWHM of the Gaussian is convSupport_ / 4
1128	// To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
1129	Int len = convSampling_ * convSupport_ ;
1130	Double hwhm = len * 0.25 ;
1131	for ( Int i = 0 ; i < len ; i++ ) {
1132	Double val = Double(i) / hwhm ;
1133	convFunc(i) = exp( -log(2)valval ) ;
1134	}
1135	}
1136
1137	void STGrid::pbFunc( Vector<Float> &convFunc )
1138	{
1139	convFunc = 0.0 ;
1140	}
1141
1142	void STGrid::setConvFunc( Vector<Float> &convFunc )
1143	{
1144	convSupport_ = userSupport_ ;
1145	if ( convType_ == "BOX" ) {
1146	if ( convSupport_ < 0 )
1147	convSupport_ = 0 ;
1148	Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1149	convFunc.resize( convSize ) ;
1150	boxFunc( convFunc, convSize ) ;
1151	}
1152	else if ( convType_ == "SF" ) {
1153	if ( convSupport_ < 0 )
1154	convSupport_ = 3 ;
1155	Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1156	convFunc.resize( convSize ) ;
1157	spheroidalFunc( convFunc ) ;
1158	}
1159	else if ( convType_ == "GAUSS" ) {
1160	// to take into account Gaussian tail
1161	if ( convSupport_ < 0 )
1162	convSupport_ = 12 ; // 3 * 4
1163	else {
1164	convSupport_ = userSupport_ * 4 ;
1165	}
1166	Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1167	convFunc.resize( convSize ) ;
1168	gaussFunc( convFunc ) ;
1169	}
1170	else if ( convType_ == "PB" ) {
1171	if ( convSupport_ < 0 )
1172	convSupport_ = 0 ;
1173	pbFunc( convFunc ) ;
1174	}
1175	else {
1176	throw AipsError( "Unsupported convolution function" ) ;
1177	}
1178	}
1179
1180	string STGrid::saveData( string outfile )
1181	{
1182	LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
1183	double t0, t1 ;
1184	t0 = mathutil::gettimeofday_sec() ;
1185
1186	//Int polno = 0 ;
1187	String outfile_ ;
1188	if ( outfile.size() == 0 ) {
1189	if ( infile_.lastchar() == '/' ) {
1190	outfile_ = infile_.substr( 0, infile_.size()-1 ) ;
1191	}
1192	else {
1193	outfile_ = infile_ ;
1194	}
1195	outfile_ += ".grid" ;
1196	}
1197	else {
1198	outfile_ = outfile ;
1199	}
1200	Table tab ;
1201	prepareTable( tab, outfile_ ) ;
1202	IPosition dshape = data_.shape() ;
1203	Int nrow = nx_ * ny_ * npol_ ;
1204	tab.rwKeywordSet().define( "nPol", npol_ ) ;
1205	tab.addRow( nrow ) ;
1206	Vector<Double> cpix( 2 ) ;
1207	cpix(0) = Double( nx_ - 1 ) * 0.5 ;
1208	cpix(1) = Double( ny_ - 1 ) * 0.5 ;
1209	Vector<Double> dir( 2 ) ;
1210	ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
1211	ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
1212	ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
1213	Int irow = 0 ;
1214	Vector<Float> sp( nchan_ ) ;
1215	Bool bsp, bdata ;
1216	const Float *data_p = data_.getStorage( bdata ) ;
1217	Float wsp_p, sp_p ;
1218	const Float *wdata_p = data_p ;
1219	long step = nx_ * ny_ * npol_ ;
1220	long offset ;
1221	for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
1222	dir(1) = center_(1) - ( cpix(1) - (Double)iy ) * celly_ ;
1223	for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
1224	dir(0) = center_(0) - ( cpix(0) - (Double)ix ) * cellx_ ;
1225	for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
1226	offset = ix + iy * nx_ + ipol * nx_ * ny_ ;
1227	//os << "offset = " << offset << LogIO::POST ;
1228	sp_p = sp.getStorage( bsp ) ;
1229	wsp_p = sp_p ;
1230	wdata_p = data_p + offset ;
1231	for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
1232	wsp_p = wdata_p ;
1233	wsp_p++ ;
1234	wdata_p += step ;
1235	}
1236	sp.putStorage( sp_p, bsp ) ;
1237	spectraCol.put( irow, sp ) ;
1238	directionCol.put( irow, dir ) ;
1239	polnoCol.put( irow, pollist_[ipol] ) ;
1240	irow++ ;
1241	}
1242	}
1243	}
1244	data_.freeStorage( data_p, bdata ) ;
1245
1246	t1 = mathutil::gettimeofday_sec() ;
1247	os << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1248
1249	return outfile_ ;
1250	}
1251
1252	void STGrid::prepareTable( Table &tab, String &name )
1253	{
1254	Table t( infile_, Table::Old ) ;
1255	t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
1256	tab = Table( name, Table::Update ) ;
1257	}
1258
1259	}

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

Context Navigation

source: trunk/src/STGrid.cpp@ 2383

Download in other formats: