/*------------------------------------------------------------------------------*
 * File Name: OCN_e01.h															*
 * Creation: TCZ 8/1/2001														*
 * Purpose: Origin C Header for NAG functions									*
 * Copyright (c) OriginLab Corp.	2001										*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *------------------------------------------------------------------------------*/


#ifndef NAGE01
#define NAGE01

  //#importdll "ONAG" // NAG DLL prepared by OriginLab
#pragma dll(ONAG)	
#include <NAG\nag_types.h>
 
/**	e01bac
		determines a cubic spline interpolant to a given set of data.

Example1:
	Assume "Data1" Worksheet has 2 columns, the first column contain 5 data 
	and we want to put result in the second column.
		
	int m = 5, i, j;
	Nag_Spline spline;
	//Attach two Datasets to these 2 columns
	Dataset xx("data1",0);
	xx.SetSize(m);
	Dataset yy("data1",1);
	yy.SetSize(m);
	//Because Dataset cannot be the parameter of function, but vector can be.
	vector x = xx, y = yy;	
	for (i=0; i<m; ++i)
		y[i] = exp(x[i]);
	nag_1d_spline_interpolant(m, x, y, &spline);
	//put the result back to the worksheet.
	yy = y;

Example2:
	The following example program sets up data from 7 values of the exponential 
	function in the interval 0 to 1. nag 1d spline interpolant is then called 
	to compute a spline interpolant to these data.  The spline is evaluated 
	by nag 1d spline evaluate (e02bbc), at the data points and at points halfway
	between each adjacent pair of data points, and the spline values and 
	the values of ex are printed out.
		
void test_nag_1d_spline_interpolant()
{
	int MMAX = 7;
	double x[] = {0.0, 0.2, 0.4, 0.6, 0.75, 0.9, 1.0};
	int i, j;
	double y[7], fit, xarg;
	Nag_Spline spline;
	int m = MMAX;
	
	for (i=0; i<m; ++i)
		y[i] = exp(x[i]);
	nag_1d_spline_interpolant(m, x, y, &spline);
	printf("\n Number of distinct knots = %ld\n\n", m-2);
	printf(" Distinct knots located at \n\n");
	for (j=3; j<m+1; j++)
	{
		printf("%7.4f",spline.lamda[j]);
		if((j-3)%5==4 || j==m )
			printf("\n");
	}

	printf("\n\n J B-spline coeff c\n\n");
	for (j=0; j<m; ++j)
		printf(" %ld %13.4f\n",j+1,spline.c[j]);
	printf("\n J Abscissa Ordinate Spline\n\n");

	for (j=0; j<m; ++j)
	{
		nag_1d_spline_evaluate(x[j], &fit, &spline);
		printf(" %ld %13.4f %13.4f %13.4f\n",j+1,x[j],y[j],fit);
		if (j < m-1)
		{
			xarg = (x[j] + x[j+1]) * 0.5;
			nag_1d_spline_evaluate(xarg, &fit, &spline);
			printf(" %15.4f %27.4f\n",xarg,fit);
		}
	}
	nag_free(spline.lamda);
	nag_free(spline.c);
	
}


	The output is as follows:
	
	 Number of distinct knots = 5

	 Distinct knots located at 

 	 0.0000 0.4000 0.6000 0.7500 1.0000


	 J B-spline coeff c

	 1        1.0000
	 2        1.1336
	 3        1.3726
	 4        1.7827
	 5        2.1744
	 6        2.4918
	 7        2.7183

	 J Abscissa Ordinate Spline

	 1        0.0000        1.0000        1.0000
    	      0.1000                      1.1052
	 2        0.2000        1.2214        1.2214
    	      0.3000                      1.3498
	 3        0.4000        1.4918        1.4918
    	      0.5000                      1.6487
	 4        0.6000        1.8221        1.8221
   		      0.6750                      1.9640
	 5        0.7500        2.1170        2.1170
    	      0.8250                      2.2819
	 6        0.9000        2.4596        2.4596
   	 	      0.9500                      2.5857
	 7        1.0000        2.7183        2.7183
	

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, m must not be less than 4: m = _value_.
	63: The sequence x is not strictly increasing: x[_value_] = _value_, x[_value_] = _value_.
	73: Memory allocation failed.
		
	successful call of the nag_1d_spline_interpolant function.

*/

	int  nag_1d_spline_interpolant(
	int m, 	// the number of data points.
	double x[], // the ith data value of the independent variable x, for i = 1, 2, . . .,m.
	double y[], // the ith data value of the dependent variable y, for i = 1, 2, . . .,m.
	Nag_Spline *spline
	);


/**	e01bec
		computes a monotonicity-preserving piecewise cubic Hermite interpolant to a set of data points.

Example1:
	Assume "Data1" Worksheet has 5 columns, the first 2 columns contain 9 data in each column 
	and we want to put result in the third, fourth and fifth column.
	
	int i, m=11, n = 9, r;
	//Attach four Datasets to these 5 columns
	Dataset dx("data1",0);
	dx.SetSize(n);
	Dataset df("data1",1);
	df.SetSize(n);
	Dataset dd("data1",2);
	dd.SetSize(m);
	Dataset dpf("data1",3);	
	dpf.SetSize(m);
	Dataset dpx("data1",4);
	dpx.SetSize(m);

	//vector can be a parameter, but Dataset cannot.
	vector x = dx, f = df, d = dd, pf = dpf, px = dpx;
	if(0 == nag_monotonic_interpolant(n, x, f, d))
	{
		double step = (x[n-1] - x[0]) / (1.0*(m-1));
		for (i = 0; i < m; i++)
		{
			if(x[0]+ i*step > x[n-1])
				px[i] =x[n-1];
			else
				px[i] = x[0]+ i*step;
		}
		nag_monotonic_evaluate(n, x, f, d, m, px, pf);
		//write the result back to the worksheet.
		dd = d;
		dpf = pf;
		dpx = px;
	}
	else
		printf("there is some problem with the function of nag_monotonic_interpolant");
		
Example2:

	This example program reads in a set of data points, calls nag 
	monotonic interpolant to compute a piecewise monotonic interpolant, 
	and then calls nag monotonic evaluate (e01bfc) to evaluate the
	interpolant at equally spaced points.

void test_nag_monotonic_interpolant()
{

	int i, m, n, r;
	double step, d[50], pf[50], px[50];

	n = 9;
	double x[50] = {7.99, 8.09, 8.19, 8.70, 9.20, 10.00, 12.00, 15.00, 20.00};
	double f[50] = {0.00000E+0, 0.27643E-4, 0.43750E-1, 0.16918E+0, 0.46943E+0,
					0.94374E+0, 0.99864E+0, 0.99992E+0, 0.99999E+0};
	
	if(0 == nag_monotonic_interpolant(n, x, f, d))
	{
		m = 11;
		step = (x[n-1] - x[0]) / (1.0*(m-1));
		for (i = 0; i < m; i++)
		{
			if(x[0]+ i*step > x[n-1])
				px[i] =x[n-1];
			else
				px[i] = x[0]+ i*step;
		}
		nag_monotonic_evaluate(n, x, f, d, m, px, pf);
	printf("                    Interpolated\n");
	printf("      Abscissa      Value\n");
		for (i = 0; i < m; i++)
			printf("%13.4f%13.4f\n", px[i], pf[i]);
	}
	else
		printf("there is some problem with the function of nag_monotonic_interpolant");
}

	The output is as follows:
	
					Interpolated
	Abscissa 		Value
	
	 7.9900 		0.0000
	 9.1910 		0.4640
	10.3920 		0.9645
	11.5930 		0.9965
	12.7940 		0.9992
	13.9950 		0.9998
	15.1960 		0.9999
	16.3970 		1.0000
	17.5980 		1.0000
	18.7990 		1.0000
	20.0000 		1.0000


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 2: n = _value_.
	236: On entry, x[r - 1] = x[r] for r = _value_: x[r - 1], x[r] = _values_.  The values of x[r], for r = 0, 1, . . . , n - 1, are not in strictly increasing order.
		
	successful call of the nag_monotonic_interpolant function.

*/
	int  nag_monotonic_interpolant(
	int n, 		// the number of data points.
	double x[], // the rth value of the independent variable (abscissa), for r = 0, 1 . . . , n - 1.
	double f[], // the rth value of the dependent variable (ordinate), for r = 0, 1, . . . , n - 1.
	double d[]  // d[r] contains the derivative at x[r].
	);

/**	e01bfc
		evaluates a piecewise cubic Hermite interpolant at a set of points.

Example1:
	Assume "Data1" Worksheet has 5 columns, the first 3 columns contain 9 data in each column 
	and we want to put result in the fourth and fifth column.
	
	int i, m = 11, n = 9, r;
	//Attach four Datasets to these 5 columns
	Dataset dx("data1",0);
	dx.SetSize(n);
	Dataset df("data1",1);
	df.SetSize(n);
	Dataset dd("data1",2);
	dd.SetSize(n);
	Dataset dpf("data1",3);	
	dpf.SetSize(m);
	Dataset dpx("data1",4);
	dpx.SetSize(m);
	//vector can be a parameter, but Dataset cannot.
	vector x = dx, f = df, d = dd, pf = dpf, px = dpx;	
	double step = (x[n-1] - x[0]) / (1.0*(m-1));
	for (i = 0; i < m; i++)
	{
		if(x[0]+ i*step > x[n-1])
			px[i] =x[n-1];
		else
			px[i] = x[0]+ i*step;
	}
	nag_monotonic_evaluate(n, x, f, d, m, px, pf);
	//put the result back to worksheet.
	dpf = pf;
	dpx = px;
			
Example2:
	This example program reads in values of n, x, f, d and m, and 
	then calls nag monotonic evaluate to evaluate the interpolant 
	at equally spaced points.

void test_nag_monotonic_evaluate()
{

	int i, m, n, r;
	double step, pf[50], px[50];
	
	n = 9;
	m = 11;
	
	double x[50] = {7.99, 8.09, 8.19, 8.70, 9.20, 10.00, 12.00, 15.00, 20.00};
	double f[50] = {0.00000E+0, 0.27643E-4, 0.43750E-1, 0.16918E+0, 0.46943E+0,
					0.94374E+0, 0.99864E+0, 0.99992E+0, 0.99999E+0};
	double d[50] = {0.00000E+0, 5.52510E-4, 0.33587E+0, 0.34944E+0, 0.59696E+0,
					6.03260E-2, 8.98335E-4, 2.93954E-5, 0.00000E+0};

	step = (x[n-1] - x[0]) / (1.0*(m-1));
	for (i = 0; i < m; i++)
	{
		if(x[0]+ i*step > x[n-1])
			px[i] =x[n-1];
		else
			px[i] = x[0]+ i*step;
	}
	nag_monotonic_evaluate(n, x, f, d, m, px, pf);
	printf("                    Interpolated\n");
	printf("      Abscissa      Value\n");
	for (i = 0; i < m; i++)
		printf("%13.4f%13.4f\n", px[i], pf[i]);
}


	The output is as follows:
	
                 Interpolated
   Abscissa      Value
    7.9900       0.0000
    9.1910       0.4640
   10.3920       0.9645
   11.5930       0.9965
   12.7940       0.9992
   13.9950       0.9998
   15.1960       0.9999
   16.3970       1.0000
   17.5980       1.0000
   18.7990       1.0000
   20.0000       1.0000


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 2: n = _value_.  On entry, m must not be less than 1: m = _value_.
	236: On entry, x[r - 1] = x[r] for r = _value_: x[r - 1], x[r] = _values_.  The values of x[r], for r = 0, 1, . . . , n - 1, are not in strictly increasing order.
	237: Warning - some points in array PX lie outside the range x[0]. . .x[n - 1]. Values at these points are unreliable as they have been computed by extrapolation.
		
	successful call of the nag_monotonic_evaluate function.

*/
	int  nag_monotonic_evaluate(
	int n, 
	double x[], 
	double f[],
	double d[], // Input: n, x, f and d must be unchanged from the previous call of nag monotonic interpolant (e01bec).
	int m, // the number of points at which the interpolant is to be evaluated.
	double px[], // Input: the m values of x at which the interpolant is to be evaluated.
	double pf[]  // Output: pf[i] contains the value of the interpolant evaluated at the point px[i], for i = 0, 1 . . .,m - 1.

	);
	
/**	e01bgc
		evaluates a piecewise cubic Hermite interpolant and its first derivative at a set of points.

Example1:
	Assume "Data1" Worksheet has 6 columns and 10 rows, the first 3 columns contain 5 data in each column 
	and we want to put result in the fourth to sixth column. 
	
	int i, m = 7, n = 5, r;
	double step;
	//Attach six Datasets to these 6 columns
	Dataset dx("data1",0);
	dx.SetSize(n);
	Dataset df("data1",1);
	df.SetSize(n);
	Dataset dd("data1",2);
	dd.SetSize(n);
	Dataset dpf("data1",3);	
	dpf.SetSize(m);
	Dataset dpx("data1",4);
	dpx.SetSize(m);
	Dataset dpd("data1",5);
	dpd.SetSize(m);
	//vector can be a parameter, but Dataset cannot.
	vector x = dx, f = df, d = dd, pf = dpf, px = dpx, pd = dpd;
	step = (x[n-1]-x[0]) / (double)(m-1);
	for (i = 0; i < m; i++)
	{
		if(x[0]+ i*step > x[n-1])
			px[i] =x[n-1];
		else
			px[i] = x[0]+ i*step;
	}
	nag_monotonic_deriv(n, x, f, d, m, px, pf, pd);
	//Put the results to the Worksheet.
	dpf = pf;
	dpx = px;
	dpd = pd;
	
					
Example2:
	This example program reads in values of n, x, f and d and calls 
	nag monotonic deriv to compute the values of the interpolant and 
	its derivative at equally spaced points.

void test_nag_monotonic_deriv()
{

	int i, m, n, r;
	double pd[21], pf[21], px[21], step;
	n = 9;
	m = 11;
	
	double x[50] = {7.99, 8.09, 8.19, 8.70, 9.20, 10.00, 12.00, 15.00, 20.00};
	double f[50] = {0.00000E+0, 0.27643E-4, 0.43750E-1, 0.16918E+0, 0.46943E+0,
					0.94374E+0, 0.99864E+0, 0.99992E+0, 0.99999E+0};
	double d[50] = {0.00000E+0, 5.52510E-4, 0.33587E+0, 0.34944E+0, 0.59696E+0,
					6.03260E-2, 8.98335E-4, 2.93954E-5, 0.00000E+0};

	
	step = (x[n-1]-x[0]) / (double)(m-1);
	for (i = 0; i < m; i++)
	{
		if(x[0]+ i*step > x[n-1])
			px[i] =x[n-1];
		else
			px[i] = x[0]+ i*step;
	}
	nag_monotonic_deriv(n, x, f, d, m, px, pf, pd);
	printf("                      Interpolated  Interpolated\n");
	printf("        Abscissa         Value       Derivative\n");
	for (i=0; i<m; i++)
		printf("%15.4f %15.4f %15.3e\n",px[i],pf[i],pd[i]);
}

	The output is as follows:
	
                  Interpolated  Interpolated
    Abscissa         Value       Derivative
     7.9900          0.0000      0.000e+000
     9.1910          0.4640      6.060e-001
    10.3920          0.9645      4.569e-002
    11.5930          0.9965      9.917e-003
    12.7940          0.9992      6.249e-004
    13.9950          0.9998      2.708e-004
    15.1960          0.9999      2.809e-005
    16.3970          1.0000      2.034e-005
    17.5980          1.0000      1.308e-005
    18.7990          1.0000      6.297e-006
    20.0000          1.0000      -3.388e-021


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 2: n = _value_.  On entry, m must not be less than 1: m = _value_.
	236: On entry, x[r - 1] = x[r] for r = _value_: x[r - 1 ] =_value_, x[r] = _value_.  The values of x[r], for r = 0, 1, . . . , n - 1a re not in strictly increasing order.
		
	successful call of the nag_monotonic_deriv function.

*/
	int  nag_monotonic_deriv(
	int n, 
	double x[], 
	double f[],
	double d[], // Input: n, x, f and d must be unchanged from the previous call of nag monotonic interpolant (e01bec).
	int m, // the number of points at which the interpolant is to be evaluated.
	double px[], // Input: the m values of x at which the interpolant is to be evaluated.
	double pf[], // Output: pf[i] contains the value of the interpolant evaluated at the point px[i], for i = 0, 1, . . .,m - 1.
	double pd[]  // Output: pd[i] contains the .rst derivative of the interpolant evaluated at the point px[i], for i = 0, 1, . . .,m - 1.
	);

/**	e01bhc
		evaluates the definite integral of a piecewise cubic Hermite interpolant over the interval [a, b].

Example1:
	Assume "Data1" Worksheet has 3 columns and 9 rows, and contain 9 data in each column. 
	"Data2" Worksheet has 3 columns and 4 rows, the first two columns contain 4 data in 
	each column.  This piece of code reads in values of "Data1", and then reads in pair
	of "Data2" first two columns, and evaluates the definite integral of the intepolant 
	over the interval [a, b].
	
	int n = 9, nn = 4;
	//Attach six Datasets to these 6 columns
	Dataset dx("data1",0);
	dx.SetSize(n);
	Dataset df("data1",1);
	df.SetSize(n);
	Dataset dd("data1",2);
	dd.SetSize(n);
	Dataset da("data2",0);
	da.SetSize(nn);
	Dataset db("data2",1);
	db.SetSize(nn);
	Dataset dintegral("data2",2);
	dintegral.SetSize(nn);
	//vector can be a parameter, but Dataset cannot.
	vector x = dx, f = df, d = dd, a = da, b = db, integral = dintegral;
	for(int i = 0; i < nn; i++)
	{
		nag_monotonic_intg(n, x, f, d, a[i], b[i], &integral[i]);
	}
	//put the result to the "data2" third column.
	dintegral = integral;
	
	
Example2:
	This example program reads in values of n, x, f and d. It then reads 
	in pairs of values for a and b, and evaluates the definite integral 
	of the interpolant over the interval [a, b] until end-of-file is reached.

void test_nag_monotonic_intg()
{
	double integral;
	int n, r;

	n = 9;
	double x[50] = {7.99, 8.09, 8.19, 8.70, 9.20, 10.00, 12.00, 15.00, 20.00};
	double f[50] = {0.00000E+0, 0.27643E-4, 0.43750E-1, 0.16918E+0, 0.46943E+0,
					0.94374E+0, 0.99864E+0, 0.99992E+0, 0.99999E+0};
	double d[50] = {0.00000E+0, 5.52510E-4, 0.33587E+0, 0.34944E+0, 0.59696E+0,
					6.03260E-2, 8.98335E-4, 2.93954E-5, 0.00000E+0};
 
	double a[4] = {7.99, 10.0, 12.0, 15.0};
	double b[4] = {20.0, 12.0, 10.0, 15.0};
	printf("					    Integral\n");
	printf("	a 		  b 		    over (a,b)\n");
	for(int i = 0; i < 4; i++)
	{
		nag_monotonic_intg(n, x, f, d, a[i], b[i], &integral);
		printf("%13.4f %13.4f %13.4f\n",a[i],b[i],integral);
	}
}


	The output is as follows:
	
						    Integral
	a 		  b	 		    over (a,b)
     7.9900   20.0000       10.7648
    10.0000   12.0000        1.9622
    12.0000   10.0000       -1.9622
    15.0000   15.0000        0.0000
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, n must not be less than 2: n = _value_.
	236: On entry, x[r - 1] = x[r] for r = _value_: x[r - 1] = _value_, x[r] = _value_.  The values of x[r], for r = 0, 1, . . . , n - 1 are not in strictly increasing order.
	237: On entry, limits a, b must not be outside interval [x[0], x[n - 1]], a = _value_, b = _value_, x[0] = _value_, x[_value_] = _value_. Extrapolation was performed to compute the integral.  The value returned is therefore unreliable.
		
	successful call of the nag_monotonic_intg function.

*/
	int  nag_monotonic_intg(
	int n, 
	double x[], 
	double f[],
	double d[], // Input: n, x, f and d must be unchanged from the previous call of nag monotonic interpolant (e01bec).
	double a,   
	double b,	// Input: the interval [a, b] over which integration is to be performed.
	double *integral // Output: the value of the definite integral of the interpolant over the interval [a, b].
	);

/**	e01dac
		computes a bicubic spline interpolating surface through a set 
		of data values, given on a rectangular grid in the x-y plane.

Example1:
	Assume we have a Worksheet "data1" which has three columns and 20
	rows.  The first column shows the x axis value, and there are 5 data 
	in this column.  The second column shows the y axis value, and there
	are 4 data in this column.  The third column contains 20 data.  This 
	piece of code computes a bicubic spline interpolating surface through
	a set of data values.
	
	int mx = 5, my = 4;
	Nag_2dSpline spline;
	//Attach three Datasets to these 3 columns
	Dataset dx("data1",0);
	dx.SetSize(mx);
	Dataset dy("data1",1);
	dy.SetSize(my);
	Dataset df("data1",2);
	df.SetSize(20);
	//vector can be a parameter, but Dataset cannot.
	vector x = dx, y = dy, f = df;
	nag_2d_spline_interpolant(mx, my, x, y, f, &spline);
			
Example2:
   	This program reads in values of mx, xq for q = 1, 2, . . .,mx, my 
	and yr for r = 1, 2, . . .,my, followed by values of the 
	ordinates fq,r defined at the grid points (xq, yr). It then calls 
	nag 2d spline interpolant to compute a bicubic spline interpolant 
	of the data values, and prints the values of the knots and 
	B-spline coefficients. Finally it evaluates the spline at a 
	small sample of points on a rectangular grid.

void test_nag_2d_spline_interpolant()
{
	int i, j, mx, my, npx, npy;
	double fg[400], tx[20], ty[20];
	double xhi, yhi, xlo, ylo, step;
	Nag_2dSpline spline;

	mx = 7;
	my = 6;
	double x[20] = {1.0, 1.10, 1.30, 1.50, 1.60, 1.80, 2.00};

	double y[20] = {0.00, 0.10, 0.40, 0.70, 0.90, 1.00};

	double f[400] = {1.00, 1.10, 1.40, 1.70, 1.90, 2.00, 1.21, 1.31, 1.61, 1.91,
					 2.11, 2.21, 1.69, 1.79, 2.09, 2.39, 2.59, 2.69, 2.25, 2.35,
					 2.65, 2.95, 3.15, 3.25, 2.56, 2.66, 2.96, 3.26, 3.46, 3.56,
					 3.24, 3.34, 3.64, 3.94, 4.14, 4.24, 4.00, 4.10, 4.40, 4.70,
					 4.90, 5.00};
					 

	nag_2d_spline_interpolant(mx, my, x, y, f, &spline);
	
	printf("Distinct knots in x direction located at\n");
	for (j=3; j<spline.nx-3; j++)
	{
		printf("%12.4f",spline.lamda[j]);
		if((j-3)%5==4 || j==spline.nx - 4)
			printf("\n");
	}

	printf("\nDistinct knots in y direction located at\n");	
	for (j=3; j<spline.ny-3; j++)
	{
		printf("%12.4f",spline.mu[j]);
		if((j-3)%5==4 || j==spline.ny-4)
			printf("\n");
	}

	printf("\nThe B-Spline coefficients:\n");
	for (i=0; i<mx; i++)
	{
		for (j=0; j<my; j++)
			printf("%9.4f",spline.c[my*i+j]);
		printf("\n");
	}

	npx = 6;
	npy = 6;
	xlo = 1.0;
	ylo = 0.0;
	xhi = 2.0;
	yhi = 1.0;

	step = (xhi-xlo)/(1.0*(npx-1));
	printf("\nSpline evaluated on a regular mesh (x across, y down): \n ");
	
	for (i=0; i<npx; i++)
	{
		if(xlo + i * step > xhi)
			tx[i] = xhi;
		else
			tx[i] = xlo + i * step;
		printf(" %5.2f ",tx[i]);
	}
	
	step = (yhi-ylo)/(npy-1);
	for (i=0; i<npy; i++)
	{
		if(ylo + i * step > yhi)
			ty[i] = yhi;
		else
			ty[i] = ylo + i * step;
	}
	nag_2d_spline_eval_rect(npx, npy, tx, ty, fg, &spline);

	printf("\n");
	for (j=0; j<npy; j++)
	{
		printf("%5.2f",ty[j]);
		for (i=0; i<npx; i++)
			printf("%8.3f ",fg[npy*i+j]);
		printf("\n");
	}
	nag_free(spline.lamda);
	nag_free(spline.mu);
	nag_free(spline.c);

}


	The output is as follows:
	
	
	Distinct knots in x direction located at
      1.0000      1.3000      1.5000      1.6000      2.0000

	Distinct knots in y direction located at
      0.0000      0.4000      0.7000      1.0000

	The B-Spline coefficients:
	   1.0000   1.1333   1.3667   1.7000   1.9000   2.0000
	   1.2000   1.3333   1.5667   1.9000   2.1000   2.2000
	   1.5833   1.7167   1.9500   2.2833   2.4833   2.5833
	   2.1433   2.2767   2.5100   2.8433   3.0433   3.1433
	   2.8667   3.0000   3.2333   3.5667   3.7667   3.8667
	   3.4667   3.6000   3.8333   4.1667   4.3667   4.4667
	   4.0000   4.1333   4.3667   4.7000   4.9000   5.0000

	Spline evaluated on a regular mesh (x across, y down): 
	   1.00   1.20   1.40   1.60   1.80   2.00 
	 0.00   1.000    1.440    1.960    2.560    3.240    4.000 
	 0.20   1.200    1.640    2.160    2.760    3.440    4.200 
	 0.40   1.400    1.840    2.360    2.960    3.640    4.400 
	 0.60   1.600    2.040    2.560    3.160    3.840    4.600 
	 0.80   1.800    2.240    2.760    3.360    4.040    4.800 
	 1.00   2.000    2.440    2.960    3.560    4.240    5.000 

	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, mx must not be less than 4: mx = _value_.  On entry, my must not be less than 4: my = _value_.
	63: The sequence x is not strictly increasing: x[_value_] = _value_, x[_value_] = _value_.  The sequence y is not strictly increasing: y[_value_] = _value_, y[_value_] = _value_.
	73: Memory allocation failed.
	243: An intermediate set of linear equations is singular, the data is too ill-conditioned to compute B-spline coefficients.
		
	successful call of the nag_2d_spline_interpolant function.
	
*/
	int  nag_2d_spline_interpolant(
	int mx, 
	int my, // Input: mx and my must specify mx and my respectively, the number of points along the x and y axis that define the rectangular grid.
	double x[], 
	double y[], // Input: x[q-1] and y[r-1] must contain xq, for q = 1, 2, . . .,mx, and yr, for r = 1, 2, . . .,my, respectively.
	double f[], // Input: f[my × (q - 1) + r - 1] must contain fq,r, for q = 1, 2, . . .,mx; r = 1, 2, . . .,my.
	Nag_2dSpline *spline
	);
	
/**	e01sac
		generates a two-dimensional surface interpolating a set of 
		scattered data points, using either the method of Renka and 
		Cline or a modification of Shepard's method.

Example1:
	Assume we have Worksheet "data1" which has 9 columns and 30 rows.
	The first three columns contain 30 data in each column.  This example
	shows user to use nag_2d_scat_interpolant(e01sac), nag_2d_scat_eval(e01sbc),
	and nag_ad_scat_free(e01szc).  The result will be put the in the rest
	six columns of "data1".
		
	int i, isel, j, m = 30, n = 0, nx = 7, ny = 6;
	double xhi = 21.0, xlo = 3.0, yhi = 17.0, ylo = 2.0;
	Nag_Scat_Struct comm;
	Nag_2d_Scat_Method method;
	Nag_E01_Opt optional;
	//Attach nine Datasets to these 9 columns
	Dataset dx("data1",0);
	dx.SetSize(m);
	Dataset dy("data1",1);
	dy.SetSize(m);
	Dataset df("data1",2);
	df.SetSize(m);
	Dataset dpx("data1",3);
	dpx.SetSize(100);
	Dataset dpy("data1",4);
	dpy.SetSize(100);
	Dataset dpf("data1",5);
	dpf.SetSize(100);
	Dataset dmpx("data1",6);
	dmpx.SetSize(100);
	Dataset dmpy("data1",7);
	dmpy.SetSize(100);
	Dataset dmpf("data1",8);
	dmpf.SetSize(100);
	//vector can be a parameter, but Dataset cannot.
	vector x = dx, y = dy, f = df, px = dpx, py = dpy, pf = dpf;
	
	for (isel=1; isel<=2; ++isel)
	{
		if (isel==1)
			method = Nag_RC;
		else if (isel==2)
			method = Nag_Shep;
	
		if (method==Nag_RC)
			nag_2d_scat_interpolant(method, m, x, y, f, &comm, NULL);		
		if (method==Nag_Shep)
		{		
			optional.nq = 24;
			optional.nw = 12;
			optional.rnq = -1.0;
			nag_2d_scat_interpolant(method, m, x, y, f, &comm, &optional);
			printf("   optional.rnw =%8.2f optional.rnq =%8.2f\n\n",	optional.rnw, optional.rnq);
			printf("   minimum number of data points that lie within radius 	optional.rnq =%3ld\n", optional.minnq);
		}
		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;
			}
		}
					
		if (method == Nag_RC)
		{
			nag_2d_scat_eval(&comm, n, px, py, pf);
			//put the result to "data1" the fourth column to the sixth column.
			dpx = px;
			dpy = py;
			dpf = pf;
		}
		if (method == Nag_Shep)
		{
			nag_2d_scat_eval(&comm, n, px, py, pf);
			//put the result to "data1" the seventh column to ninth column.
			dmpx = px;
			dmpy = py;
			dmpf = pf;
		}
		
		nag_2d_scat_free(&comm);
	}
		
				
Example2:
	This program reads in a set of 30 data points and calls 
	nag 2d scat interpolant to construct an interpolating surface. 
	It then calls nag 2d scat eval (e01sbc) to evaluate the 
	interpolant at a sample of points on a rectangular grid. 
	The two methods described in Section 3 are used in the 
	construction of the interpolating surface and the subsequent 
	evaluation of the interpolant.  Note that this example is not 
	typical of a realistic problem: the number of data points would 
	normally be larger, and the interpolant would need to be evaluated 
	on a finer grid to obtain an accurate plot, say.
	
void test_nag_2d_scat_interpolant()
{

	int i, isel, j, m, n, nx, ny;
	double xhi, xlo, yhi, ylo;
	double pf[100], px[100], py[100];
	Nag_Scat_Struct comm;
	Nag_2d_Scat_Method method;
	Nag_E01_Opt optional;

	nx = 7;
	ny = 6;
	xlo = 3.0;
	ylo = 2.0;
	xhi = 21.0;
	yhi = 17.0;
	m = 30;	
	double x[100] = {11.16, 12.85, 19.85, 19.72, 15.91, 0.00, 20.87, 3.45, 14.26, 17.43,
					 22.80, 7.58, 25.00, 0.00, 9.66, 5.22, 17.25, 25.00, 12.13, 22.23,
					 11.52, 15.20, 7.54, 17.32, 2.14, 0.51, 22.69, 5.47, 21.67, 3.31};
	
  	double y[100] = {1.24, 3.06, 10.72, 1.39, 7.74, 20.00, 20.00, 12.78, 17.87, 3.46,
  					 12.39, 1.98, 11.87, 0.00, 20.00, 14.66, 19.57, 3.87, 10.79, 6.21,
  					 8.53, 0.00, 10.69, 13.78, 15.03, 8.37, 19.63, 17.13, 14.36, 0.33};

  
	double f[100] = {22.15, 22.11, 7.97, 16.83, 15.30, 34.60, 5.74, 41.24, 10.74, 18.60,
					 5.47, 29.87, 4.40, 58.20, 4.73, 40.36, 6.43, 8.74, 13.71, 10.25,
					 15.74, 21.60, 19.31, 12.11, 53.10, 49.43, 3.25, 28.63, 5.52, 44.08};
	
	for (isel=1; isel<=2; ++isel)
	{
		if (isel==1)
			method = Nag_RC;
		else if (isel==2)
			method = Nag_Shep;
	
		if (method==Nag_RC)
			nag_2d_scat_interpolant(method, m, x, y, f, &comm, NULL);
		
		if (method==Nag_Shep)
		{
			
			optional.nq = 24;
			optional.nw = 12;
			optional.rnq = -1.0;
			nag_2d_scat_interpolant(method, m, x, y, f, &comm, &optional);
			printf("   optional.rnw =%8.2f optional.rnq =%8.2f\n\n",	optional.rnw, optional.rnq);
			printf("   minimum number of data points that lie within radius 	optional.rnq =%3ld\n", optional.minnq);
		}
	
		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;
			}
		}
		
		if (method == Nag_RC)
			nag_2d_scat_eval(&comm, n, px, py, pf);
		if (method == Nag_Shep)
			nag_2d_scat_eval(&comm, n, px, py, pf);
			
		printf("\n   px\n");
		for (i = 0; i < nx; i++)
			printf("%8.2f", px[i]);
		printf("\n   py       pf\n");
		for (i = ny-1; i >= 0; --i)
		{
			printf("%8.2f ", py[nx * i]);
			for (j = 0; j < nx; j++)
				printf("%8.2f", pf[nx * i + j]);
			printf("\n");
		}
		nag_2d_scat_free(&comm);
	}
}



	The output is as follows:
	
	
	px
   	 3.00    6.00    9.00   12.00   15.00   18.00   21.00
   	py       pf
  	17.00    41.25   27.62   18.03   12.29   11.68    9.09    5.37
  	14.00    47.61   36.66   22.87   14.02   13.44   11.20    6.46
  	11.00    38.55   25.25   16.72   13.83   13.08   10.71    6.88
   	8.00    37.90   23.97   16.79   16.43   15.46   13.02    9.30
    5.00    40.49   29.26   22.51   20.72   19.30   16.72   12.87
    2.00    43.52   33.91   26.59   22.23   21.15   18.67   14.88
    optional.rnw =    9.49 optional.rnq =   13.42

    minimum number of data points that lie within radius 	optional.rnq =  7

    px
     3.00    6.00    9.00   12.00   15.00   18.00   21.00
    py       pf
    17.00    40.23   27.72   21.23   14.59   12.00    9.43    5.46
    14.00    46.96   37.37   23.74   14.67   13.25   11.29    6.26
    11.00    39.42   25.42   16.32   13.78   12.60   10.39    7.03
     8.00    37.50   22.36   18.57   15.63   15.55   13.05    9.69
     5.00    41.25   31.76   24.74   21.17   18.93   16.83   12.65
     2.00    44.58   34.35   26.47   22.27   20.98   18.69   15.06

	 
Parameters:

Return:
	This function returns NAG error code, 0 if no error.

	70: On entry, parameter method had an illegal value.
	11: On entry, m must not be less than 3: m = _value_.
	73: Memory allocation failed.
	245: On entry, all the (x,y) pairs are collinear. Consider specifying method = Nag Shep or usinga one-dimensional interpolating function nag 1d spline interpolant (e01bac).
	246: On entry, each data pair is not unique since data points (x[_value_],y[_value_]) and (x[_value_],y[_value_]) are identical and equal to (_value_,_value_).
	247: On entry, either or both of optional.nq and optional.nw are invalid, optional.nq = _value_ and optional.nw = _value_. optional.nq and optional.nw must satisfy the following constraints: 0 < optional.nw = optional.nq. 
	248: On entry, either or both of optional.rnq and optional.rnw are invalid, optional.rnq = _value_ and optional.rnw = _value_. optional.rnq and optional.rnw must satisfy the following constraints: 0.0 < optional.rnw = optional.rnq.
	249: The minimum number of data points _value_ that lie within the radius optional.rnq of any node is small enough to indicate that the interpolant may be unsatisfactory in regions where the data points are sparse. Current values of other relevant parameters (available as members of the structure optional, if this has been defined) are rnq = _value_, rnw = _value_, nq = _value_, nw = _value_.  
		
	successful call of the nag_2d_scat_interpolant function.
	
*/
	int  nag_2d_scat_interpolant(
	Nag_2d_Scat_Method method,
	int m, 		// the number of data points, m.
	double x[], // Input: the co-ordinates of the rth data point, for r = 1, 2, . . .,m. The data points are accepted in any order, but see Section 6.
	double y[], // Input: the co-ordinates of the rth data point, for r = 1, 2, . . .,m. The data points are accepted in any order, but see Section 6.
	double f[], // Input: the co-ordinates of the rth data point, for r = 1, 2, . . .,m. The data points are accepted in any order, but see Section 6.
	Nag_Scat_Struct *comm, // Pointer to a communication structure of type Nag Scat Struct. For method = Nag RC, this structure contains the computed triangulation and the estimated partial derivatives at the nodes. For method = Nag Shep, this structure contains the coefficients of the constructed nodal functions.  This structure must be passed unchanged to the interpolant evaluating function nag 2d scat eval (e01sbc).
	Nag_E01_Opt *optional  // Pointer to structure of type Nag E01 Opt, which may be used only if method = Nag Shep.
	);

/**	e01sbc
		evaluates at given points the two-dimensional interpolant function
		computed by nag 2d scat interpolant (e01sac).
		
Example:

	See the example program for nag 2d scat interpolant (e01sac).

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	178: The setup function nag 2d scat interpolant (e01sac) has not been called.
	179: The call to setup function nag 2d scat interpolant (e01sac) produced an error.
	81: The structure _value_ has been corrupted since the previous call to _value_.
	11: On entry, n must not be less than 1: n = _value_.
	251: The evaluation point (_value_, _value_) of (px, py) lies outside the triangulation boundary.  The returned value, _value_,of pf was computed by extrapolation.
	252: On entry,the interpolant cannot be evaluted because the evaluation point (px, py) of (_value_, _value_) is outside the support region of the input data points defined by optional.rnw= _value_ as set in nag 2d scat interpolant (e01sac).  
	
	successful call of the nag_2d_scat_eval function.

*/
	int  nag_2d_scat_eval(
	Nag_Scat_Struct *comm, // Pointer to a communication structure of type Nag Scat Struct which must be unchanged from the previous call of nag 2d scat interpolant (e01sac).
	int n, // the number of points at which the evaluation of the interpolant is required.
	double px[],  
	double py[], // Input: the x- and y-coordinates of the kth point (pxk, pyk),fo r k = 1, 2, . . . ,n,at which the interpolant is to be evaluated.
	double pf[]  // Output: the values of the interpolant evaluated at the points (pxk, pyk),fo r k = 1, 2, . . . ,n.
	);
/**	e01szc
		frees the memory allocated by NAG function nag 2d scat interpolant
		(e01sac) to the communication structure comm of type Nag Scat Struct.
		
Example:

	See nag 2d scat interpolant (e01sac) for an example of hown ag 2d scat free is used.

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	successful call of the nag_2d_scat_free function.
	
*/
	void nag_2d_scat_free(
	Nag_Scat_Struct *comm // Input: the communication structure that was used in the call to nag 2d scat interpolant (e01sac).  Output: pointers which pointed to NAG allocated memory will have been freed and set to NULL.
	);
   
#endif /* not NAGE01 */