/*------------------------------------------------------------------------------*
 * File Name: NAG 2D Grid RC				 									*
 * Creation: ER, 08/17/2001														*
 * Purpose: OriginC program to demonstrate NAG's 2D Renka-Cline Gridding		*
 * Copyright (c) OriginLab Corp.2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007	*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:  LAS, 9/10/01 pass wksData and wksCopy to function			*
 *                    ER, 7/30/02 introduced new parameters, qil, wfl for       *
 *                                controlling gridding parameters for Shepard   *
/*------------------------------------------------------------------------------*/
 
////////////////////////////////////////////////////////////////////////////////////
// you can include just this typical header file for most Origin built-in functions and classes
// and it takes a reasonable amount of time to compile, 
#include <origin.h>
//#include <data.h>
//#include <math.h>
//#include <nag\OCN_e01.h>
//#include <nag\nag_types.h>

////////////////////////////////////////////////////////////////////////////////////


////////////////////////////////////////////////////////////////////////////////////
// This function accepts random xyz data and then performs 2D gridding on the data
// using the NAG library functions. Two NAG algorithms for gridding are available:
// Renka-Cline, and Shepard. For more documentation, view Chapter e01 of the 
// NAG PDF files
// 
// Parameters:
//				wksData: 	name of worksheet with data
//				wksCopy: 	name of worksheet with regular/output xyz data
//				mthd:  		conversion method:
//							0: Renka-Cline
//							non-zero: Shepard 
//				qil:		scale factor for qudratic interpolant locality, for Shepard method
//				wfl:		scale factor for weight function locality, for Shepard method
//	
// Return:		0 = success
//				
bool Grid2D(string wksData, string wksCopy, int mthd, double qil=0.0, double wfl=0.0)
{
	int i, isel, j, m, n, nx, ny;
	double xhi, xlo, yhi, ylo;
	
	// datasets for the scatter points
	Dataset xData(wksData,0);
	Dataset yData(wksData,1);
	Dataset zData(wksData,2);
		
	m = xData.GetSize();			// get number of scatter points
	
	nx = 2 * sqrt(m);				// set grid x,y size accordingly
	ny = nx;
	
	xlo = min(xData);				// get lo and hi values of x and y scatter
	ylo = min(yData);
	xhi = max(xData);
	yhi = max(yData);
	

	// datasets for grid points
	Dataset pxData(wksCopy,0);
	Dataset pyData(wksCopy,1);
	Dataset pzData(wksCopy,2);
	
	// set their size to equal grid dimensions
	pxData.SetSize(nx*ny);
	pyData.SetSize(nx*ny);
	pzData.SetSize(nx*ny);
	

	// Create vectors to hold scatter and grid points
	// We cannot pass datasets to NAG functions - need to use vectors
	vector<double> x, y, z, px, py, pz;

	// set vector dimensions accordingly
	x.SetSize(m);
	y.SetSize(m);
	z.SetSize(m);
	px.SetSize(nx*ny);
	py.SetSize(nx*ny);
	pz.SetSize(nx*ny);

	// copy scatter points into vectors
	x = xData;
	y = yData;
	z = zData;		

	// set up NAG structures for calling gridding routine
	Nag_Scat_Struct comm;
	Nag_2d_Scat_Method method;	
	Nag_E01_Opt optional;
	
	// Set method based on flag passed to this routine
	if (mthd == 0) method = Nag_RC;
	else method = Nag_Shep;
	
	// call the appropriate interpolation routine
	if (method == Nag_RC)
	{
		nag_2d_scat_interpolant(method, m, x, y, z, &comm, NULL);
	}
	else
	{
		// default values of parameters for Shepard's method
		double nq = 24.0;
		double nw = 12.0;
		double rnq = -1.0;
		// scale defaults with values passed by user; leave rnq as is
		nq *= qil;			// quadratic interpolant locality
		nw *= wfl;			// weight function locality
		// now set these values for the conversion method
		optional.nq = nq;
		optional.nw = nw;
		optional.rnq = rnq;
		nag_2d_scat_interpolant(method, m, x, y, z, &comm, &optional);
	}
			
	// now compute positions of the grid points using the hi and lo scatter values
	n = 0;		
	for (j=0; j<ny; ++j)
	{
		for (i=0; i<nx; ++i)
		{
			px[i+nx*j] = (1.0 * (nx-i-1) / (nx-1)) * xlo + (1.0*i / (nx-1)) * xhi;
			py[i+nx*j] = (1.0 * (ny-j-1) / (ny-1)) * ylo +	(1.0* j / (ny-1)) * yhi;
			++n;
		}
	}
		

	// 	evaluate two-dimensional interpolant function computed by the interpolant function call
	nag_2d_scat_eval(&comm, n, px, py, pz);		
			
	// write results to datasets
	pxData = px;
	pyData = py;
	pzData = pz;
	
	// clean up 
	nag_2d_scat_free(&comm);
	
	// success; return 0
	return 0;	
}
// end
