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

#ifndef _O_NAG_E02_H
#define _O_NAG_E02_H
#pragma dll(ONAG)
#include <NAG\nag_types.h>
/** e02adc
		computes weighted least-squares polynomial approximations to an arbitrary set of data points.

Example 1:	
	Assume we have "Data1" Worksheet and "Matrix1" Matrix.  "Data1" Worksheet has 4 columns, the first 
	3 columns contain 4 data in each column and we want to put result in the fourth column and at 
	"Matrix1" matrix.
	
	int tdim = 10, m = 4;
	//Attach four Datasets to these 4 columns
	Dataset xx("data1",0);
	Dataset yy("data1",1);
	Dataset ww("data1",2);
	Dataset ss("data1",3);
	xx.SetSize(m);
	yy.SetSize(m);
	ww.SetSize(m);
	ss.SetSize(m);
	//Because Dataset cannot be the parameter of function, but vector can be.
	vector x = xx;
	vector y = yy;
	vector w = ww;
	vector s = ss;
	//So you need 4 vectors for this function.
	//Because matrix can be the parameter of function, but Matrix cannot.
	matrix a;
	a.SetSize(20, 20);		//use a as a parameter of a function.
	
	//Call function here, x, y, w, s are four vectors, and a is a matrix.
	int success = nag_1d_cheb_fit(m, 4, tdim, x, y, w, a, s);
	
	//After call the function, let the fourth dataset eqauls to the fourth vector by using
	ss = s;	
	//After get result, let a Matrix equals to matrix by using
	Matrix aa("Matrix1");
	aa = a;	
			
Example 2:
	Determine weighted least-squares polynomial approximations of degrees 0, 1, 2 and 3 to 
	a set of 11 prescribed data points.

void test_nag_1d_cheb_fit()
{
	double d1;
	double xarg;
	matrix a;
	a.SetSize(50, 50);
	double s[200], xcapr;
	double x[200] = {1.00, 2.10, 3.10, 3.90, 4.90, 5.80, 6.50, 7.10, 7.80, 8.40, 9.00};
	double y[200] = {10.40, 7.90, 4.70, 2.50, 1.20, 2.20, 5.10, 9.20, 16.10, 24.50, 35.30};
	double w[200] = {1.00, 1.00, 1.00, 1.00, 1.00, 0.80, 0.80, 0.70, 0.50, 0.30, 0.20};
	double x1, ak[50], xm, fit;
	int i, j, k, m, r;
	int iwght = 2;
	int tdim = 50;
	int success;

	printf("e02adc Example Program Results \n");
	m = 11;
	k = 3;

	printf("the input data are as following:\n");
	printf("m = 11, k = 2, iwght = 3\n");
	printf("x:\n");
	for(i = 0; i < m; i++)
	{
		printf("%3.2f  ",x[i]);
		if((i + 1) % 5 == 0)
			printf("\n");			
	}
	printf("\ny:\n");
	for(i = 0; i < m; i++)
	{
		printf("%3.2f  ",y[i]);
		if((i + 1) % 5 == 0)
			printf("\n");			
	}
	printf("\nw:\n");
	for(i = 0; i < m; i++)
	{
		printf("%3.2f  ",w[i]);
		if((i + 1) % 5 == 0)
			printf("\n");			
	}
	
	success = nag_1d_cheb_fit(m, 4, tdim, x, y, w, a, s);
	
	if( success == 0)
	{
		for (i = 0; i <= k; ++i)
		{
			printf("\n");
			printf(" %s%4ld%s%12.2e\n","Degree",i," R.M.S. residual =",s[i]);
			printf("\n J Chebyshev coeff A(J) \n");
			for (j = 0; j < i+1; ++j)
				printf(" %3ld%15.4f\n",j+1,a[i][j]);
		}
		for (j = 0; j < k+1; ++j)
			ak[j] = a[k][j];
		x1 = x[0];
		xm = x[m-1];
		printf("\n %s%4ld\n","Polynomial approximation and residuals for degree",k);
		printf("\n R Abscissa Weight Ordinate Polynomial Residual\n");
		for (r = 1; r <= m; ++r)
		{
			xcapr = (x[r-1] - x1 - (xm - x[r-1])) / (xm - x1);
			nag_1d_cheb_eval(k+1, ak, xcapr, &fit);
			d1 = fit - y[r-1];
			printf(" %3ld%11.4f%11.4f%11.4f%11.4f%11.2e\n",r,x[r-1],w[r-1],y[r-1],fit,d1);
			if (r < m)
			{
				xarg = (x[r-1] + x[r]) * 0.5;
				xcapr = (xarg - x1 - (xm - xarg)) / (xm - x1);
				success = nag_1d_cheb_eval(k+1, ak, xcapr, &fit);
				if( success == 0)
					printf(" %11.4f %11.4f\n",xarg,fit);
				else
				{
					printf("there is some problem with the function");
					break;
				}					
			}
		}
	}
	else
		printf("there is some problem with this function");
}

	
	The output is following:
	
	Degree 0 R.M.S. residual = 4.07e+00
	J Chebyshev coeff A(J)
	1 		12.1740
	
	Degree 1 R.M.S. residual = 4.28e+00
	J 	Chebyshev coeff A(J)
	1 		12.2954
	2 	 	 0.2740
	
	Degree 2 R.M.S. residual = 1.69e+00
	J 	Chebyshev coeff A(J)
	1 		20.7345
	2 	 	 6.2016
	3 	 	 8.1876
	
	Degree 3 R.M.S. residual = 6.82e-02
	J 	Chebyshev coeff A(J)
	1 		24.1429
	2 	 	 9.4065
	3 		10.8400
	4 	 	 3.0589
	
	Polynomial approximation and residuals for degree 3
	
	R 	Abscissa 	Weight 	Ordinate 	Polynomial 		Residual
	
	1 	1.0000 		1.0000 	10.4000 	10.4461 		 4.61e-02 
		1.5500 							9.3106
	2 	2.1000 		1.0000 	7.9000 		7.7977 			-1.02e-01 
		2.6000 							6.2555
	3 	3.1000 		1.0000 	4.7000 		4.7025 			 2.52e-03 
		3.5000 							3.5488
	4 	3.9000 		1.0000 	2.5000 		2.5533 			 5.33e-02 
		4.4000 							1.6435
	5 	4.9000 		1.0000 	1.2000 		1.2390 			 3.90e-02 
		5.3500 							1.4257
	6 	5.8000 		0.8000 	2.2000 		2.2425 			 4.25e-02 
		6.1500 							3.3803
	7 	6.5000 		0.8000 	5.1000 		5.0116 			-8.84e-02 
		6.8000 							6.8400
	8 	7.1000 		0.7000 	9.2000 		9.0982 			-1.02e-01 
		7.4500 							12.3171
	9 	7.8000 		0.5000 	16.1000 	16.2123 		 1.12e-01 
		8.1000 							20.1266
	10 	8.4000 		0.3000 	24.5000 	24.6048 		 1.05e-01
	 	8.7000 							29.6779
	11 	9.0000 		0.2000 	35.3000 	35.3769 		 7.69e-02 


Return:
	This function returns NAG error code, 0 if no error.
	
	successfully call of the nag_1d_cheb_fit function.	
*/

	int  nag_1d_cheb_fit(
		int m, 		//the number m of data points.
		int kplus, 	//k + 1, where k is the maximum degree required
		int tda, 	//the second dimension of the array a.
		const double x[], //the values xr of the independent variable
		const double y[], //the values yr of the dependent variable.
		const double w[], //the set of weights, wr, for r = 1, 2, . . .,m.
		double a[], //the coefficients of in the approximating polynomial of degree i
		double s[] // the root mean square residual si, for i = 0, 1, . . . , k
); // least squares curve fit with polynomials and arbitrary data points


/** e02aec
		evaluates a polynomial from its Chebyshev-series representation

Example 1:
	Assume "Data1" Worksheet has three columns, the first column contains 5 data,  
	and we want to put result in the second and third column.

	int i, m = 11, n = 4, r, success;
	//Attach three Datasets to these 3 columns
	Dataset aa("data1",0);
	Dataset dxcap("data1",1);
	Dataset pp("data1",2);
	aa.SetSize(n+1);
	dxcap.SetSize(m);
	pp.SetSize(m);
	//Because Dataset cannot be the parameter of function, but vector can be.
	vector a = aa;
	vector xcap = dxcap;
	vector p = pp;
	for(r = 0; r < m; r++)
	{
		xcap[r] = 1.0 * (2*r - m + 1 ) / (1.0* (m - 1)) ;
		success = nag_1d_cheb_eval(5, a , xcap[r], &p[r]);
	}	
	//put the result to Worksheet columns.
	dxcap = xcap;
	pp = p;
	
Example 2:
	Evaluate at 11 equally-spaced points in the interval -1 = �x = 1 the polynomial of degree 4 with
	Chebyshev coe.cients, 2.0, 0.5, 0.25, 0.125, 0.0625.

void test_nag_1d_cheb_eval()
{
	double xcap;
	double a[200] = {2.0000, 0.5000, 0.2500, 0.1250, 0.0625};
	double p;
	int i, m = 11, n = 4;
	int r;
	int success;
	printf("e02aec Example Program Results \n");

	printf("the input data:\n");
	printf("m = 11, n = 4 ");
	printf("\na:\n");
	for(r = 0; r < 5; r++)
		printf("%10.4f",a[r]);	
		
	printf("\n R Argument Value of polynomial \n");
	for(r = 0; r < m; r++)
	{
		xcap = 1.0 * (2*r - m + 1 ) / (1.0* (m - 1)) ;
		success = nag_1d_cheb_eval(5, a , xcap, &p);
		if(success == 0)
			printf(" %3ld%14.4f %35.4f\n",r,xcap,p);
		else
		{
			printf("there is some problem with the function");
			break;
		}
	}
}
		
	The output is following:
	
	R 	Argument 	Value of polynomial
	
	1 	-1.0000 	0.6875
	2 	-0.8000 	0.6613
	3 	-0.6000 	0.6943
	4 	-0.4000 	0.7433
	5 	-0.2000 	0.7843
	6 	 0.0000 	0.8125
	7 	 0.2000 	0.8423
	8 	 0.4000 	0.9073
	9 	 0.6000 	1.0603
	10 	 0.8000 	1.3733
	11 	 1.0000 	1.9375

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, nplus1 must not be less than 1: nplus1 = _value_.
	271: On entry, abs(xcap) > 1.0 + 4 _, w here _ is the machine precision.  In this case the value of p is set arbitrarily to zero.
	
	successfully call of the nag_1d_cheb_eval function.	
*/

	int  nag_1d_cheb_eval(
		int nplus1,		// the number n + 1 of terms in the series
		const double a[], // the value of the ith coe.cient in the series, for i = 1, 2, . . ., n+1.
		double xcap, //the argument at which the polynomial is to be evaluated.
		double *p  // the value of the polynomial.
); // evaluation of polynomials in one variable


/** e02afc
		computes the coeffcients of a polynomial, in its Chebyshev series form, which 
		interpolates (passes exactly through) data at a special set of points. Least-squares
		polynomial approximations can also be obtained.

Example 1:
	Assume "Data1" Worksheet has two columns, the first column contains 11 data,  
	and we want to put result in the second column.
	
	int n = 10, success;
	//Attach two Datasets to these 2 columns
	Dataset ff("data1",0);
	Dataset dan("data1",1);
	
	ff.SetSize(n+1);
	dan.SetSize(n+1);
	
	//Because Dataset cannot be the parameter of function, but vector can be.
	vector f = ff;
	vector an = dan;
	
	success = nag_1d_cheb_interp_fit(n+1, f, an);
	
	//Put the result back to the Worksheet "data1" second column.
	dan = an;
	
Example2:	
	Determine the Chebyshev coeffcients of the polynomial which interpolates 
	the data �xr, fr, for r = 1, 2, . . . , 11, where �xr = cos((r - 1)p/10) 
	and fr = e�xr . Evaluate, for comparison with the values of fr, the 
	resulting Chebyshev series at �xr, for r = 1, 2, . . . , 11.  The example 
	program supplied is written in a general form that will enable polynomial
	interpolations of arbitrary data at the cosine points cos((r - 1)p/n), 
	for r = 1, 2, . . . , n + 1 to be obtained for any n (= nplus1-1). Note 
	that nag 1d cheb eval (e02aec) is used to evaluate the interpolating 
	polynomial. The program is self-starting in that any number of data sets 
	can be supplied.
	
void test_nag_1d_cheb_interp_fit()
{
	#define NMAX 199
	#define NP1MAX NMAX+1
	double d1;
	double xcap[NP1MAX];
	double f[NP1MAX] = {2.7182, 2.5884, 2.2456, 1.7999, 1.3620, 1.0000, 
	0.7341, 0.5555, 0.4452, 0.3863, 0.3678};
	double piby2n, pi1, an[NP1MAX], fit;
	int i, j, n = 10;
	int r,success;
		
	pi1 = 3.1415926;
	piby2n = pi1 * 0.5 / n;
		
	for (r = 0; r < n+1; ++r)
	{
		i = r;
		if (2*i <= n)
		{
			d1 = sin(piby2n * i);
			xcap[i] = 1.0 - d1 * d1 * 2.0;	
		}
		else if (2*i > n * 3)
		{
			d1 = sin(piby2n * (n - i));
			xcap[i] = d1 * d1 * 2.0 - 1.0;
		}
		else
		{
			xcap[i] = sin(piby2n * (n - 2*i));
		}
	}
		
	success = nag_1d_cheb_interp_fit(n+1, f, an);
		
	printf("\n Chebyshev \n");
	printf(" J coefficient A(J) \n");
	for (j = 0; j < n+1; ++j)
		printf(" %3ld%14.7f\n",j+1,an[j]);
	printf("\n R Abscissa Ordinate Fit \n");
	for (r = 0; r < n+1; ++r)
	{
		success =nag_1d_cheb_eval(n+1, an, xcap[r], &fit);
		printf(" %3ld%11.4f%11.4f%11.4f\n",r+1,xcap[r],f[r],fit);
	}		
}

	The output is following:
		
			Chebyshev
	J 		coefficient A(J)
	
	1 		 2.5320000
	2 		 1.1303095
	3 		 0.2714893
	4 		 0.0443462
	5 	 	 0.0055004
	6 		 0.0005400
	7 		 0.0000307
	8 		-0.0000006
	9 		-0.0000004
	10 		 0.0000049
	11 		-0.0000200
	
	R 	Abscissa 		Ordinate 	Fit
	
	1 	 1.0000 		2.7182 		2.7182
	2 	 0.9511 		2.5884 		2.5884
	3 	 0.8090 		2.2456 		2.2456
	4 	 0.5878 		1.7999 		1.7999
	5 	 0.3090 		1.3620 		1.3620
	6 	 0.0000 		1.0000 		1.0000
	7 	-0.3090 		0.7341 		0.7341
	8 	-0.5878 		0.5555 		0.5555
	9 	-0.8090 		0.4452 		0.4452
	10 	-0.9511 		0.3863 		0.3863
	11 	-1.0000 		0.3678 		0.3678
			
Return:
	This function returns NAG error code, 0 if no error
	
	11: On entry, nplus1 must not be less than 2: nplus1 = _value_.
	
	successfully call of the nag_1d_cheb_interp_fit function.	
*/

	int nag_1d_cheb_interp_fit(
		int nplus1,	//the number n + 1 of data points (one greater than the degree n of the interpolating polynomial.
		const double  f[],
		double a[] //the coefficient a[j] in the interpolating polynomial,for j = 1,2,.., n + 1.
); //Least-squares curve fit with polynomials,selected data points


/** e02bac
		computes a weighted least-squares approximation to an arbitrary set of data points by a 
		cubic spline with knots prescribed by the user. Cubic spline interpolation can also be carried out.
		Return NAG error, 0 if no error.

Example1:
	Assume "Data1" Worksheet has three columns, every column has 14 data.  After call the function,
	we get the residual sum of squares and the result of a Nag_Spline structure variable.
	
	//Attach three Datasets to these 3 columns
	
	int m=14;
	Dataset xx("data1",0);
	Dataset yy("data1",1);
	Dataset dweights("data1",2);
	
	xx.SetSize(m);
	yy.SetSize(m);
	dweights.SetSize(m);
	
	
	double ss;
	Nag_Spline spline;
	vector x, y, weights;
	x = xx;
	y = yy;
	weights = dweights;
	BOOL success = nag_1d_spline_fit_knots(m, x, y, weights, &ss, &spline);	
	  		
Example2:
	Determine a weighted least-squares cubic spline approximation with five intervals (four interior knots) 
	to a set of 14 given data points. Tabulate the data and the corresponding values of the approximating 
	spline, together with the residual errors, and also the values of the approximating spline at points 
	half-way between each pair of adjacent data points.

void test_nag_1d_spline_fit_knots()
{
	int ncap, ncap7, j, m, r, wght;
	double weights[200] = {0.20, 0.20, 0.30, 0.70, 0.90, 1.00, 1.00, 1.00, 0.80, 0.50, 0.70, 1.00, 1.00, 1.00};;
	double xarg, ss, fit;
	double fg[12];
	double x[200] = {0.20, 0.47, 0.74, 1.09, 1.60, 1.90, 2.60, 3.10, 4.00, 5.15, 6.17, 8.00, 10.00, 12.00};
	double y[200] = {0.00, 2.00, 4.00, 6.00, 8.00, 8.62, 9.10, 8.90, 8.15, 7.00, 6.00, 4.54, 3.39, 2.56};
	Nag_Spline spline;
	int success;
	m = 14;
	ncap = 5;
	wght = 2;
	ncap7 = ncap+7;
	spline.n = ncap7;

	fg[4] = 1.50;
	fg[5] = 2.60;
	fg[6] = 4.00;
	fg[7] = 8.00;
	spline.lamda = fg;
	
	printf("spline = %d\n", spline.n);
	ASSERT(spline.n == 12);
	
	int nSize=sizeof(Nag_Spline);
	
 	success = nag_1d_spline_fit_knots(m, x, y, weights, &ss, &spline);
 	
	printf("\nNumber of distinct knots = %ld\n\n", ncap+1);
	printf("Distinct knots located at \n\n");
	for (j=3; j<ncap+4; j++)
	{
		printf("%8.4f",spline.lamda[j]);
		if((j-3)%6==5 || j==ncap+3)
			printf("\n");
	}
	printf("\n\n J B-spline coeff c\n\n");
	for (j=0; j<ncap+3; ++j)
	printf(" %ld %13.4f\n",j+1,spline.c[j]);
	printf("\nResidual sum of squares = ");
	printf("%9.2e\n\n",ss);
	printf("Cubic spline approximation and residuals\n");
	printf(" r Abscissa Weight Ordinate	Spline Residual\n\n");
	for (r=0; r<m; ++r)
	{
		nag_1d_spline_evaluate(x[r], &fit, &spline);
		printf("%3ld %11.4f %11.4f %11.4f %11.4f %10.1e\n",r+1,x[r],weights[r],y[r],fit,fit-y[r]);
		if (r<m-1)
		{
			xarg = (x[r] + x[r+1]) * 0.5;
			nag_1d_spline_evaluate(xarg, &fit, &spline);
			printf(" %14.4f %35.4f\n",xarg,fit);
		}
	}
	nag_free(spline.c);

}
	
	The output is following:
	spline = 12
	
	Number of distinct knots = 6
	
	Distinct knots located at
	
	0.2000 1.5000 2.6000 4.0000 8.0000 12.0000
	
	J 	B-spline coeff c
	
	1 	-0.0465
	2 	 3.6150
	3 	 8.5724
	4 	 9.4261
	5 	 7.2716
	6 	 4.1207
	7  	 3.0822
	8 	 2.5597
	
	Residual sum of squares = 1.78e-03
	
	Cubic spline approximation and residuals
	
	r 	Abscissa 	Weight 		Ordinate 	Spline 		Residual
	
	1 	0.2000 		0.2000 		0.0000 		-0.0465 	-4.7e-02
		0.3350 								 1.0622
	2 	0.4700 		0.2000 		2.0000 		 2.1057 	 1.1e-01
		0.6050 								 3.0817
	3 	0.7400 		0.3000 		4.0000 		 3.9880 	-1.2e-02
		0.9150 								 5.0558
	4 	1.0900 		0.7000 		6.0000 		 5.9983 	-1.7e-03
		1.3450 								 7.1376
	5 	1.6000 		0.9000 		8.0000 		 7.9872 	-1.3e-02
		1.7500 								 8.3544
	6 	1.9000 		1.0000 		8.6200 		 8.6348 	 1.5e-02
		2.2500 								 9.0076
	7 	2.6000 		1.0000 		9.1000 		 9.0896 	-1.0e-02
		2.8500 								 9.0353
	8 	3.1000 		1.0000 		8.9000 		 8.9125 	 1.2e-02
		3.5500 								 8.5660
	9 	4.0000 		0.8000 		8.1500 		 8.1321 	-1.8e-02
		4.5750 								 7.5592
	10 	5.1500 		0.5000 		7.0000 		 6.9925 	-7.5e-03
		5.6600 								 6.5010
	11 	6.1700 		0.7000 		6.0000 		 6.0255 	 2.6e-02
		7.0850 								 5.2292
	12  8.0000 		1.0000 		4.5400 		 4.5315 	-8.5e-03
		9.0000 								 3.9045
	13 10.0000 		1.0000 		3.3900 		 3.3928 	 2.8e-03
	   11.0000 								 2.9574
	14 12.0000 		1.0000 		2.5600 		 2.5597 	-3.5e-04

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, spline.n must not be less than 8: spline.n = _value_.
	73: Memory allocation failed.
	239: On entry, user-specified knots must be interior to the data interval, spline.lamda[4] must be greater than x[0] and spline.lamda[spline.n-5] must be less than x[m-1]: spline.lamda[4] = _value_, x[0] = _value_, spline.lamda[_value_] = _value_, x[_value_] = _value_.  
	62: The sequence spline.lamda is not increasing: spline.lamda[_value_] = _value_, spline.lamda[_value_] = _value_.  This condition on spline.lamda applies to user-speci.ed knots in the interval spline.lamda[4], spline.lamda[spline.n-5].  The sequence x is not increasing: x[_value_] = _value_, x[_value_] = _value_.  
	240: On entry, the weights are not strictly positive: weights[_value_] = _value_.
	244: Too many knots for the number of distinct abscissae, mdist: spline.n = _value_, mdist = _value_.  These must satisfy the constraint spline.n = mdist + 4.
	241: The conditions speci.ed by Schoenberg and Whitney fail.  The conditions speci.ed by Schoenberg andWhitney (1953) fail to hold for at least one subset of the distinct data abscissae. That is, there is no subset of spline.n-4 strictly increasing values, x[r0], x[r1],. . . ,x[rspline.n-5], among the abscissae such that x[r0] < spline.lamda[0] < x[r4], x[r1] < spline.lamda[1] < x[r5], . . . x[rspline.n-9] < spline.lamda[spline.n-9] < x[rspline.n-5].  This means that there is no unique solution: there are regions containing too many knots compared with the number of data points.  
		
	successfully call of the nag_1d_spline_fit_knots function.
*/

	int  nag_1d_spline_fit_knots(
		int m, 					//the number m of data points.
		const double x[], 		//the values xr of the independent variable (abscissa), for r = 1, 2,...,m.
		const double y[], 		//the values yr of the of the dependent variable (ordinate), for r = 1, 2,...,m.
		const double weights[], //the values wr of the weights, for r = 1, 2, . .,m. 
		double *ss, 			//the residual sum of squares
		Nag_Spline *spline  		//pointer to structure of type Nag_Spline.
); //Least-squares curve fit with cubic splines e02bac


/** e02bbc
		evaluates a cubic spline from its B-spline representation.
		
Example:
	Evaluate at 9 equally-spaced points in the interval [1.0, 9.0] the cubic spline with 
	(augmented) knots 1.0, 1.0, 1.0, 1.0, 3.0, 6.0, 8.0, 9.0, 9.0, 9.0, 9.0 and normalised 
	cubic B-spline coefficients 1.0, 2.0, 4.0, 7.0, 6.0, 4.0, 3.0.
	
void test_nag_1d_spline_evaluate()
{
	double a, b,ncap7;
	int m =9;
	int ncap = 4;
	int r;
	double lamda[11] = {1.00, 1.00, 1.00, 1.00, 3.00, 6.00, 8.00, 9.00, 9.00, 9.00, 9.00};
	double c[11] = {1.00, 2.00, 4.00, 7.00, 6.00, 4.00, 3.00};
	double s;
	double x;
	int i, j;
	Nag_Spline spline;
	ncap7 = ncap +7;
	spline.n = ncap7;	
	spline.lamda = lamda;
	spline.c = c;

	a = spline.lamda[3];
	printf("a = %f\n",a);
	b = spline.lamda[ncap+3];
	printf("b = %f\n",b);
	printf("Augmented set of knots stored in spline.lamda:\n");
	for(j = 0; j < ncap7-1; j++)
		printf("%10.4f",lamda[j]);	
	printf("\nB-spline coefficients stored in spline.c\n\n");
	for(j = 0; j< ncap+3; j++)
		printf("%10.4f",c[j]);		
	printf("\n x Value of cubic spline\n\n");
	for (r=1; r<=m; ++r)
	{
		x = 1.0*((m-r) * a + (r-1) * b) /(m-1);
		s= 0;
 		nag_1d_spline_evaluate(x, &s, &spline);
		printf("%10.4f%15.4f\n",x,s);
	}		
}

	The output is following:
	
	Augmented set of knots stored in spline.lamda:
	
	1.0000 	1.0000 	1.0000 	1.0000 	3.0000 	6.0000
	8.0000 	9.0000 	9.0000 	9.0000 	9.0000
	
	B-spline coefficients stored in spline.c:
	1.0000 	2.0000 	4.0000 	7.0000 	6.0000 	4.0000
	3.0000
	
	x 			Value of cubic spline
	1.0000 		1.0000
	2.0000 		2.3779
	3.0000 		3.6229
	4.0000 		4.8327
	5.0000 		5.8273
	6.0000 		6.3571
	7.0000 		6.1905
	8.0000 		5.1667
	9.0000 		3.0000

Return:
	This function returns NAG error, 0 if no error.
	
	11: On entry, spline.n must not be less than 8: spline.n = _value_.
	238: On entry, x must satisfy spline.lamda[3] = x = spline.lamda[spline.n-4]: spline.lamda[3] = _value_, x = _value_, spline.lamda[_value_] = _value_.  In this case s is set arbitrarily to zero.
	
	successfully call of the nag_1d_spline_evaluate function.
*/

	int  nag_1d_spline_evaluate(
		double x, 			//the argument x at which the cubic spline is to be evaluated.
		double *s, 			//the value of the spline, s(x).
		Nag_Spline *spline  //Pointer to structure of type Nag_Spline.
); //Evaluation of cubic splines.


/** e02bcc
		evaluates a cubic spline and its first three derivatives from its B-spline
		representation.
		
Example:
	Compute, at the 7 arguments x = 0, 1, 2 , 3, 4, 5, 6, the left- and right-hand 
	values and first 3 derivatives of the cubic spline defined over the interval 
	0 <= x <= 6 having the 6 interior knots x = 1, 3, 3, 3, 4, 4, the 8 additional 
	knots 0, 0, 0, 0, 6, 6, 6, 6, and the 10 B-spline coe.cients 10, 12, 13, 15, 
	22, 26, 24, 18, 14, 12.  The input data items (using the notation of Section 4) 
	comprise the following values in the order indicated: 
	
	�n 								m 
	spline.lamda[j], 				for j = 0, 1, . . . , �n + 6 
	spline.c[j], 					for j = 0, 1, . . . , �n + 2 
	x 								m values of x
	
	The example program is written in a general form that will enable the values 
	and derivatives of a cubic spline having an arbitrary number of knots to be 
	evaluated at a set of arbitrary points. Any number of data sets may be supplied.
	
void test_nag_1d_spline_deriv()
{
	int i, j, l, m, ncap, ncap7;
	double s[4];
	double lamd[14] = {0.0, 0.0, 0.0, 0.0, 1.0, 3.0, 3.0, 3.0, 4.0, 4.0, 6.0, 6.0, 6.0, 6.0};
	double c1[14] = {10.0, 12.0, 13.0, 15.0, 22.0, 26.0, 24.0, 18.0, 14.0, 12.0};
	double x[7] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
	Nag_Spline spline;
	Nag_DerivType derivs;

	printf("e02bcc Example Program Results\n");
	ncap = 7;
	m = 7;
	ncap7 = ncap+7;
   	spline.n = ncap7;
	spline.lamda = lamd;
	spline.c = c1;
	
	for( i = 0; i < 11; i++)
		printf("%f	",spline.lamda[i]);
	
	printf(" x Spline 1st deriv 2nd deriv 3rd deriv");
	for (i= 0; i< m; i++)
	{
		derivs = Nag_LeftDerivs;
		for (j=1; j<=2; j++)
		{
		 	nag_1d_spline_deriv(derivs, x[i], s, &spline);
		 	
		 	if(derivs ==Nag_LeftDerivs)
			{
				printf("\n\n%11.4f Left",x[i]);
				for (l=0; l<4; l++)
					printf("%11.4f",s[l]);
			}
			else
			{
				printf("\n%11.4f Right",x[i]);
				for (l=0; l<4; l++)
					printf("%11.4f",s[l]);
			}
			derivs = Nag_RightDerivs;
		}	
	}	
}		
		
	The output is following:	

	x 				Spline 		1st deriv 	2nd deriv 	3rd deriv
	
	0.0000 	Left 	10.0000 	6.0000 		-10.0000 	10.6667
	0.0000 	Right 	10.0000 	6.0000 		-10.0000 	10.6667

	1.0000 	Left 	12.7778 	1.3333 		  0.6667 	10.6667
	1.0000 	Right 	12.7778 	1.3333 		  0.6667 	 3.9167

	2.0000 	Left 	15.0972 	3.9583 	 	  4.5833 	 3.9167
	2.0000 	Right 	15.0972 	3.9583 		  4.5833 	 3.9167

	3.0000 	Left 	22.0000 	10.5000 	  8.5000 	 3.9167
	3.0000 	Right 	22.0000 	12.0000 	-36.0000 	36.0000

	4.0000 	Left 	22.0000 	-6.0000  	  0.0000 	36.0000
	4.0000 	Right 	22.0000 	-6.0000 	  0.0000 	 1.5000

	5.0000 	Left 	16.2500 	-5.2500 	  1.5000 	 1.5000
	5.0000 	Right 	16.2500 	-5.2500 	  1.5000	 1.5000

	6.0000 	Left 	12.0000 	-3.0000 	  3.0000	 1.5000
	6.0000 	Right 	12.0000 	-3.0000 	  3.0000	 1.5000
			
Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, spline.n must not be less than 8: spline.n = _value_.
	70: On entry, parameter derivs had an illegal value.
	238: On entry, x must satisfy spline.lamda[3] = x = spline.lamda[spline.n-4]: spline.lamda[3] = _value_, x = _value_, spline.lamda[_value_] = _value_.
	269: On entry, the cubic spline range is invalid: spline.lamda[3] = _value_ while spline.lamda[spline.n-4] = _value_.  These must satisfy spline.lamda[3] < spline.lamda[spline.n-4].
	
	successfully call of the nag_1d_spline_deriv function.		
*/

	int  nag_1d_spline_deriv(
	Nag_DerivType derivs, 
		double x, 	//the argument x at which the cubic spline and its derivatives are to be evaluated.
		double s[], // the value of the jth derivative of the spline at the argument x, for j = 0, 1, 2, 3.
	Nag_Spline *spline // Pointer to structure of type Nag_Spline
); //Evaluation of fitted cubic spline, function and derivatives


/** e02bdc
		computes the definite integral of a cubic spline from its B-spline
		representation.
	
Example:
	Determine the definite integral over the interval [0,6] of a cubic spline having 6 interior 
	knots at the positions lamda= 1, 3, 3, 3, 4, 4, the 8 additional knots 0, 0, 0, 0, 6, 6, 6, 6, 
	and the 10 B-spline coefficients 10, 12, 13, 15, 22, 26, 24, 18, 14, 12.

void test_nag_1d_spline_intg()
{
	int j;
	double integral;
	Nag_Spline spline;
	spline.n = 14;
	double lamda[] = {0.0, 0.0, 0.0, 0.0, 1.0, 3.0, 3.0, 3.0,4.0, 4.0, 6.0, 6.0, 6.0, 6.0};
	double c[] = {10.0, 12.0, 13.0, 15.0, 22.0, 26.0, 24.0, 18.0,14.0, 12.0, 0, 0, 0}
	spline.c = c;
	spline.lamda = lamda;
	
	nag_1d_spline_intg(&spline, &integral);
	
	printf("Definite integral = %11.3e\n",integral);
}

	The output is following:

Definite integral = 1.000e+02

Return:
	The function returns NAG error, 0 if no error.
	
	11: On entry, spline.n must not be less than 8: spline.n = _value_.
	253: On entry, the knots must satisfy the following constraints: spline.lamda[spline.n-4] > spline.lamda[3], spline.lamda[j] = spline.lamda[j - 1], for j = 1, 2, . . . ,spline.n-1, with equality in the cases j = 1, 2, 3,spline.n-3,spline.n-2 and spline.n-1.
		
	successfully call of the nag_1d_spline_intg function.
*/
	
	int  nag_1d_spline_intg(
	 	Nag_Spline *spline, //pointer to structure of type Nag_Spline
		double *integral	//the value of the definite integral of s(x) between the limits x = a and x = b,
); // Evaluation of fitted cubic spline, definite integral


/** e02bec
		computes a cubic spline approximation to an arbitrary set of data points.
		
Example:
	This example program reads in a set of data values, followed by a set of values of S. For each value
	of S it calls nag_1d_spline_fit to compute a spline approximation, and prints the values of the knots
	and the B-spline coefficients ci.

void test_nag_1d_spline_fit()
{
	int NEST = 54;
	int m = 15, r, j;
	double weights[50] = {1.00, 2.00, 1.50, 1.00, 3.00, 1.00, 0.50, 1.00, 2.00, 2.50, 1.00, 3.00, 1.00, 2.00, 1.00};  
	double x[50] = {0.0000E+00, 5.0000E-01, 1.0000E+00, 1.5000E+00, 2.0000E+00, 2.5000E+00, 3.0000E+00, 4.0000E+00, 
	4.5000E+00, 5.0000E+00, 5.5000E+00, 6.0000E+00, 7.0000E+00, 7.5000E+00, 8.0000E+00};
	double y[50] = {-1.1000E+00, -3.7200E-01, 4.3100E-01, 1.6900E+00, 2.1100E+00, 3.1000E+00, 4.2300E+00, 4.3500E+00,
	4.8100E+00, 4.6100E+00, 4.7900E+00, 5.2300E+00, 6.3500E+00, 7.1900E+00, 7.9700E+00};
	double s[3] = {1.0, 0.5, 0.1};
	double fp, sp[99], txr;
	Nag_Start start;
	Nag_Comm warmstartinf;
	Nag_Spline spline;
	
	start = Nag_Cold;

	for(int i = 0; i < 3; i ++)
	{
		nag_1d_spline_fit(start, m, x, y, weights, s[i], NEST, &fp, &warmstartinf, &spline);
		
		for (r=0; r<m; r++)
			nag_1d_spline_evaluate(x[r], &sp[r*2], &spline);
		for (r=0; r<m-1; r++)
		{
			txr = (x[r] + x[r+1]) / 2;
			nag_1d_spline_evaluate(txr, &sp[r*2], &spline);
		}
		
		printf("\nCalling with smoothing factor s = %12.3e\n",s);
		printf("\nNumber of distinct knots = %ld\n\n", spline.n-6);
		printf("Distinct knots located at \n\n");
		for (j=3; j<spline.n-3; j++)
		{
			printf("%8.4f",spline.lamda[j]);
			if((j-3)%6==5 || j==spline.n-4)
				printf("\n");
		}
		printf("\n\n J B-spline coeff c\n\n");
		for (j=0; j<spline.n-4; ++j)
			printf(" %3ld %13.4f\n",j+1,spline.c[j]);
		printf("\nWeighted sum of squared residuals fp = %12.3e\n",fp);
		if (fp == 0.0)
			printf("The spline is an interpolating spline\n");
		if (spline.n == 8)
			printf("The spline is the weighted least-squares cubic	polynomial\n");
		start = Nag_Warm;	
	}
	nag_free(spline.c);
	nag_free(spline.lamda);
	nag_free(warmstartinf.nag_w);
	nag_free(warmstartinf.nag_iw);	

}

	The output is following:
	
	Calling with smoothing factor s = 1.000e+00
	
	Number of distinct knots = 3
	
	Distinct knots located at
	
	0.0000 4.0000 8.0000
	
	J 	B-spline coeff c
	
	1 	-1.3201
	2 	 1.3542
	3 	 5.5510
	4 	 4.7031
	5 	 8.2277
	
	Weighted sum of squared residuals fp = 1.000e+00
	
	Calling with smoothing factor s = 5.000e-01
	
	Number of distinct knots = 7
	
	Distinct knots located at
	
	0.0000 1.0000 2.0000 4.0000 5.0000 6.0000	8.0000
	
	J 	B-spline coeff c
	
	1 	-1.1072
	2 	-0.6571
	3 	 0.4350
	4 	 2.8061
	5 	 4.6824
	6 	 4.6416
	7 	 5.1976
	8 	 6.9008
	9 	 7.9979
	
	Weighted sum of squared residuals fp = 5.001e-01
	
	Calling with smoothing factor s = 1.000e-01
	
	Number of distinct knots = 10
	
	Distinct knots located at
	
	0.0000 1.0000 1.5000 2.0000 3.0000 4.0000	4.5000 5.0000 6.0000 8.0000
	
	J 	B-spline coeff c
	1 	-1.0900
	2 	-0.6422
	3 	 0.0369
	4 	 1.6353
	5 	 2.1274
	6 	 4.5526
	7 	 4.2225
	8 	 4.9108
	9 	 4.4159
	10 	 5.4794
	11 	 6.8308
	12   7.9935
	
	Weighted sum of squared residuals fp = 1.000e-01
	
Return:
	The function returns NAG error, 0 if no error.
	
	70: On entry, parameter start had an illegal value.
	11: On entry, m must not be less than 4: m = _value_.  On entry, nest must not be less than 8: nest = _value_.
	5: On entry, s must not be less than 0.0: s = _value_.
	240: On entry, the weights are not strictly positive: weights[_value_] = _value_.
	63: The sequence x is not strictly increasing: x[_value_] = _value_ x[_value_] = _value_.
	268: On entry, nest = _value_, s = _value_, m = _value_.  Constraint: nest = m + 4 when s = 0.0.
	73: Memory allocation failed.
	266: Start has been set to Nag Warm at the .rst call of this function. It must be set to Nag Cold at the FIrst call. 
	264: The number of knots needed is greater than nest, nest = _value_. If nest is already large, say nest > m/2, this may indicate that possibly s is too small: s = _value_.  
	254: The iterative process has failed to converge. Possibly s is too small: s = _value_.  If the function fails with an error exit of NE SPLINE COEFF CONV or NE NUM KNOTS ID GT, a spline approximation is returned, but it fails to satisfy the .tting criterion (see (2) and (3) in Section 3) - perhaps by only a small amount, however. 
		
	successfully call of the nag_1d_spline_fit function.	
*/

	int  nag_1d_spline_fit(
		Nag_Start start, // start must be set to Nag Cold or Nag Warm.
		int m, // the number of data points
		const double x[], // the value x[r] of the independent variable (abscissa) x, for r = 1, 2, . . .,m.
		const double y[], // the value y[r] of the dependent variable (ordinate) y, for r = 1, 2, . . .,m.
		const double weights[], // the value w[r] of the weights,for r = 1, 2, . . .,m.
		double s, //the smoothing factor
		int nest, // an over-estimate for the number, n, of knots required.
		double *fp, //the sum of the squared weighted residuals of the computed spline approximation
		Nag_Comm *warmstartinf, // Pointer to structure of type Nag_Comm
		Nag_Spline *spline //Pointer to structure of type Nag_Spline
);//Least-squares cubic spline curve fit, automatic knot placement, one variable


/** e02dcc
		computes a bicubic spline approximation to a set of data values,
		given on a rectangular grid in the x-y plane.

Example:

	This example program reads in values of mx, my, xq, for q = 1, 2, . . . ,mx, 
	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 fit
	grid to compute a bicubic spline approximation for one specified value of s,
	and prints the values of the computed 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_fit_grid()
{
	int NXEST = 15;
	int NYEST = 13;
	int mx = 11, my = 9, nx, ny, npx, npy, i, j;
	double fg[49], px[7], py[7];
	double xhi, yhi, xlo, ylo, delta, s, fp;
	Nag_Start start;
	Nag_2dSpline spline;
	Nag_Comm warmstartinf;
	double x[] ={0.0000E+00, 5.0000E-01, 1.0000E+00, 1.5000E+00, 2.0000E+00,2.5000E+00,
		 	3.0000E+00, 3.5000E+00, 4.0000E+00, 4.5000E+00,5.0000E+00};
	double y[] ={0.0000E+00, 5.0000E-01, 1.0000E+00, 1.5000E+00, 2.0000E+00,2.5000E+00,
				 3.0000E+00, 3.5000E+00, 4.0000E+00}	
	double f[] ={1.0000E+00, 8.8758E-01, 5.4030E-01, 7.0737E-02, -4.1515E-01,
			-8.0114E-01, -9.7999E-01, -9.3446E-01, -6.5664E-01, 1.5000E+00,
			1.3564E+00, 8.2045E-01, 1.0611E-01, -6.2422E-01, -1.2317E+00,
			-1.4850E+00, -1.3047E+00, -9.8547E-01, 2.0600E+00, 1.7552E+00,
			1.0806E+00, 1.5147E-01, -8.3229E-01, -1.6023E+00, -1.9700E+00,
			-1.8729E+00, -1.4073E+00, 2.5700E+00, 2.1240E+00, 1.3508E+00,
			1.7684E-01, -1.0404E+00, -2.0029E+00, -2.4750E+00, -2.3511E+00,
			-1.6741E+00, 3.0000E+00, 2.6427E+00, 1.6309E+00, 2.1221E-01,
			-1.2484E+00, -2.2034E+00, -2.9700E+00, -2.8094E+00, -1.9809E+00,
			3.5000E+00, 3.1715E+00, 1.8611E+00, 2.4458E-01, -1.4565E+00,
			-2.8640E+00, -3.2650E+00, -3.2776E+00, -2.2878E+00, 4.0400E+00,
			3.5103E+00, 2.0612E+00, 2.8595E-01, -1.6946E+00, -3.2046E+00,
			-3.9600E+00, -3.7958E+00, -2.6146E+00, 4.5000E+00, 3.9391E+00,
			2.4314E+00, 3.1632E-01, -1.8627E+00, -3.6351E+00, -4.4550E+00,
			-4.2141E+00, -2.9314E+00, 5.0400E+00, 4.3879E+00, 2.7515E+00,
			3.5369E-01, -2.0707E+00, -4.0057E+00, -4.9700E+00, -4.6823E+00,
			-3.2382E+00, 5.5050E+00, 4.8367E+00, 2.9717E+00, 3.8505E-01,
			-2.2888E+00, -4.4033E+00, -5.4450E+00, -5.1405E+00, -3.5950E+00,
			6.0000E+00, 5.2755E+00, 3.2418E+00, 4.2442E-01, -2.4769E+00,
			-4.8169E+00, -5.9300E+00, -5.6387E+00, -3.9319E+00}	
	
	start = Nag_Cold;
	s = 0.1;

	nag_2d_spline_fit_grid(start, mx, x, my, y, f, s, NXEST, NYEST,
						&fp, &warmstartinf, &spline);
	nx = spline.nx;
	ny = spline.ny;
	printf("\nCalling with smoothing factor s = %11.4e:	spline.nx = %2ld, spline.ny = %2ld.\n\n",s,nx,ny);
	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\n");
	for (i=0; i<ny-4; i++)
	{
		for (j=0; j<nx-4; j++)
			printf("%9.4f",spline.c[i+j*(ny-4)]);
		printf("\n");
	}
	printf("\nSum of squared residuals fp = %11.4e\n",fp);
	if (fp==0.0)
		printf("\nThe spline is an interpolating spline\n");
	if (nx==8 && ny==8)
		printf("\nThe spline is the least-squares bi-cubic polynomial\n");

	npx = 6;
	xlo = 0.0;
	xhi = 5.0;
	npy = 5;
	ylo = 0.0;
	yhi = 4.0;
	delta = (xhi-xlo) / (npx-1);
	for (i=0; i<npx; i++)
	{
		if(xlo+i*delta > xhi)
			px[i] = xhi;
		else
			px[i] =	xlo+i*delta; 		
	}
	for (i=0; i<npy; i++)
	{
		if(ylo+i*delta > yhi)
			py[i] = yhi;
		else
			py[i] =	ylo+i*delta; 		
	}
	nag_2d_spline_eval_rect(npx, npy, px, py, fg, &spline);
	printf("\nValues of computed spline:\n");
	printf("\n x");
	for (i=0; i<npx; i++)
		printf("%7.2f ",px[i]);
	printf("\n y\n");
	for (i=npy-1; i>=0; i--)
	{
		printf("%7.2f ",py[i]);
		for (j=0; j<npx; ++j)
			printf("%8.2f ",fg[npy*j+i]);
		printf("\n");
	}
	nag_free(spline.mu);
	nag_free(spline.c);
	nag_free(spline.lamda);
	nag_free(warmstartinf.nag_w);
	nag_free(warmstartinf.nag_iw);
}

	The output is following:
	
	Calling with smoothing factor s = 1.0000e-01: spline.nx = 10, spline.ny = 13.
	
	Distinct knots in x direction located at
	
	0.0000 1.5000 2.5000 5.0000
	
	Distinct knots in y direction located at
	
	0.0000 1.0000 2.0000 2.5000 3.0000	3.5000 4.0000

	The B-spline coefficients:

	 0.9918  1.5381  2.3913  3.9845  5.2138  5.9965
	 1.0546  1.5270  2.2441  4.2217  5.0860  6.0821
	 0.6098  0.9557  1.5587  2.3458  3.3860  3.7716
	-0.2915 -0.4199 -0.7399 -1.1763 -1.5527 -1.7775
	-0.8476 -1.3296 -1.8521 -3.3468 -4.3628 -5.0085
	-1.0168 -1.5952 -2.4022 -3.9390 -5.4680 -6.1656
	-0.9529 -1.3381 -2.2844 -3.9559 -5.0032 -5.8709
	-0.7711 -1.0914 -1.8488 -3.2549 -3.9444 -4.7297
	-0.6476 -1.0373 -1.5936 -2.5887 -3.3485 -3.9330
	
	Sum of squared residuals fp = 1.0004e-01
	
	Values of computed spline:
	
	x 	0.00 1.00 2.00 3.00 4.00 5.00
	
	y
	4.00 -0.65 -1.36 -1.99 -2.61 -3.25 -3.93
	3.00 -0.98 -1.97 -2.91 -3.91 -4.97 -5.92
	2.00 -0.42 -0.83 -1.24 -1.66 -2.08 -2.48
	1.00  0.54  1.09  1.61  2.14  2.71  3.24
	0.00  0.99  2.04  3.03  4.01  5.02  6.00

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

	70: On entry, parameter start had an illegal value.
	11: On entry, mx must not be less than 4: mx = _value_.  On entry, my must not be less than 4: my = _value_.  On entry, nxest must not be less than 8: nxest = _value_.  On entry, nyest must not be less than 8: nyest = _value_.
	5: On entry, s must not be less than 0.0: s = _value_.
	268: On entry, s = _value_, nxest = _value_, mx = _value_. Constraint: nxest = mx+4 when s = 0.0.  On entry, s = _value_, nyest = _value_, my = _value_. Constraint: nyest = mx+4 when s = 0.0.
	73: Memory allocation failed.
	266: start has been set to Nag Warm at the .rst call of this function. It must be set to Nag Cold at the .rst call.
	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_.
	261: The number of knots required is greater than allowed by nxest or nyest, nxest = _value_, nyest = _value_. Possibly s is too small, especially if nxest, nyest > mx/2, my/2. s = _value_, mx = _value_, my = _value_.
	254: The iterative process has failed to converge. Possibly s is too small: s = _value_.  If the function fails with an error exit of NE NUM KNOTS 2D GT RECT or NE SPLINE COEFF CONV, then a spline approximation is returned, but it fails to satisfy the .tting criterion (see (2) and (3) in Section 3) - perhaps by only a small amount, however.
		
	successfully call of the nag_2d_spline_fit_grid function.
*/
	
	int nag_2d_spline_fit_grid(
		Nag_Start start, // start must be set to Nag Cold or Nag Warm.
		int mx, // the number of grid points along the x axis
		const double x[],//the x co-ordinate of the qth grid point along the x axis, for q = 1, 2, . . .,mx.
		int my, //the number of grid points along the y axis.
		const double y[], // the y co-ordinate of the rth grid point along the y axis, for r = 1, 2, . . .,my.
		const double f[],// contain the data value fq,r, for q = 1, 2, . . .,mx and r = 1, 2, . . .,my.
		double s, // the smoothing factor
		int nxest,
		int nyest, // an upper bound for the number of knots nx and ny required in the x and y directions respectively.
		double *fp, // the sum of squared residuals of the computed spline approximation
		Nag_Comm *warmstartinf, // Pointer to structure of type Nag_Comm
		Nag_2dSpline *spline // Pointer to structure of type Nag_2dSpline
);// Least-squares bicubic spline fit with automatic knot placement, two variables (rectangular grid).


/** e02ddc
		computes a bicubic spline approximation to a set of scattered data.
		
Example:
	This example program reads in a value of m, followed by a set of m data 
	points (xr, yr, fr) and their weights wr. It then calls nag 2d spline fit scat 
	to compute a bicubic spline approximation for one specified value of S, 
	and prints the values of the computed 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_fit_scat()
{
	double  weights[30], fg[49],px[7], py[7];
	int i, j, m, rank, nx, ny, npx, npy, nxest, nyest;
	double s, delta, fp, xhi, yhi, xlo, ylo;
	Nag_Start start;
	Nag_2dSpline spline;
	double warmstartinf;
	double x[] ={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[] ={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[] ={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};
			
	nxest = 14;
	nyest = 14;			
	for(i = 0; i < 30; i++)
		weights[i] = 1.00;
	m = 30;
	start = Nag_Cold;
	s = 10.0;
	
	nag_2d_spline_fit_scat(start, m, x, y, f, weights, s, nxest, nyest, 
							&fp,&rank, &warmstartinf, &spline);
	
	nx= spline.nx;
	ny = spline.ny;
	printf("\nCalling with smoothing factor s = %11.4e, nx= %1ld,ny = %1ld\n",s,nx,ny);
	printf("rank deficiency = %1ld\n\n",(nx-4)*(ny-4)-rank);
	printf("Distinct knots in xdirection 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\n");
	for (i=0; i <ny-4; i++)
	{
		for (j=0; j < nx-4; j++)
			printf("%9.2f",spline.c[i+j*(ny-4)]);
		printf("\n");
	}
	printf("\n Sum of squared residuals fp = %11.4e\n",fp);
	if (nx==8 && ny==8)
		printf("The spline is the least-squares bi-cubic polynomial\n");
	
	npx = 7;
	xlo = 3.0;
	xhi = 21.0;

	npy = 6;
	ylo = 2.0;
	yhi = 17.0;
	delta = (xhi-xlo)/(npx-1);
	for (i=0; i<npx; i++)
	{
		if(xlo+i*delta > xhi)
			px[i] = xhi;
		else
			px[i] =	xlo+i*delta; 		
	}
	for (i=0; i<npy; i++)
	{
		if(ylo+i*delta > yhi)
			py[i] = yhi;
		else
			py[i] =	ylo+i*delta; 		
	}
	nag_2d_spline_eval_rect(npx, npy, px, py, fg, &spline);
	
	printf("\nValues of computed spline:\n\n");
	printf(" x");
	for (i=0; i< npx;  i++)
		printf("%8.2f ",px[i]);
	printf("\n y\n");
	for (i=npy-1; i>=0; i--)
	{
		printf("%8.2f ",py[i]);
		for (j=0; j < npx; ++j)
			printf("%8.2f ",fg[npy*j+i]);
		printf("\n");
	}
	nag_free(spline.lamda);
	nag_free(spline.c);
	nag_free(spline.mu);
}


	The output is following:
	
	
	Calling with smoothing factor s = 1.0000e+01, nx= 10, ny = 9 rank deficiency = 0
	
	Distinct knots in xdirection located at
	
	0.0000 9.7575 18.2582 25.0000
	
	Distinct knots in y direction located at
	
	0.0000 9.0008 20.0000
	
	The B-spline coefficients:
	
	58.16  46.31  6.01  32.00  5.86 -23.78
	63.78  46.74 33.37  18.30 14.36  15.95
	40.84 -33.79  5.17  13.10 -4.13  19.37
	75.44 111.92  6.94  17.33  7.09 -13.24
	34.61 -42.61 25.20  -1.96 10.37  -9.09
	
	Sum of squared residuals fp = 1.0002e+01
	
	Values of computed spline:
	
	x	3.00 6.00 9.00 12.00 15.00 18.00 21.00
	
	y
	
	17.00 40.74 28.62 19.84 14.29 11.21 9.46 7.09
	14.00 48.34 33.97 21.56 14.71 12.32 10.82 7.15
	11.00 37.26 24.46 17.21 14.14 13.02 11.23 7.29
	8.00 30.25 19.66 16.90 16.28 15.21 12.71 8.99
	5.00 36.64 26.75 23.07 21.13 18.97 15.90 11.98
	2.00 45.04 33.70 26.25 22.88 21.62 19.39 13.40


Return:
	This function returns NAG error code, 0 if no error.
	
	70: On entry, parameter start had an illegal value.
	11: On entry, nxest must not be less than 8: nxest = _value_.  On entry, nyest must not be less than 8: nyest = _value_.
	5: On entry, s must not be less than or equal to 0.0: s = _value_.
	270: On entry, the number of data points with non-zero weights = _value_. Constraint: the number of non-zero weights = 16.
	73: Memory allocation failed.
	266: start has been set to Nag Warm at the .rst call of this function. It must be set to Nag Cold at the .rst call.
	260: On entry, all the values in the array x must not be equal.  On entry, all the values in the array y must not be equal.
	262: The number of knots required is greater than allowed by nxest or nyest, nxest = _value_, nyest = _value_. Possibly s is too small, especially if nxest, nyest > 5 + ��m. s = _value_, m = _value_.
	263: No more knots can be added because the number of B-spline coe.cients already exceeds m.  Either m or s is probably too small: m = _value_, s = _value_.
	265: No more knots added; the additional knot would coincide with an old one. Possibly an inaccurate data point has too large a weight, or s is too small. s = _value_.
	254: The iterative process has failed to converge. Possibly s is too small: s = _value_.  If the function fails with an error exit of NE NUM KNOTS 2D GT SCAT, NE NUM COEFF GT, NE NO ADDITIONAL KNOTS or NE SPLINE COEFF CONV, then a spline approximation is returned, but it fails to satisfy the setting criterion (see (2) and (3) in Section 3) - perhaps by only a small amount, however. 
		
	successfully call of the  nag_1d_cheb_fit function.	
*/
	
int  nag_2d_spline_fit_scat(
		Nag_Start start, 
		int m,
		const double x[],
		const double y[],
		const double f[],
		const double weights[], 
		double s,
		int nxest,
		int nyest,
		double *fp,
		int *rank,
		double *warmstartinf, 
		Nag_2dSpline *spline);


/** e02dec
		calculates values of a bicubic spline from its B-spline representation.
		
Example:
	This program reads in knot sets spline.lamda[0],. . . ,spline.lamda[spline.nx-1] 
	and spline.mu[0],. . . ,spline.mu[spline.ny and a set of bicubic spline coefficients cij. 
	Following these are a value for m and the co-ordinates (xr, yr), for r = 1, 2, . . .,m, 
	at which the spline is to be evaluated.

void test_nag_2d_spline_eval()
{
	int i, m = 7;
	double x[20] = {1.0, 1.1, 1.5, 1.6, 1.9, 1.9, 2.0}; 
	double y[20] = {0.0, 0.1, 0.7, 0.4, 0.3, 0.8, 1.0};
	double ff[20];
	Nag_2dSpline spline;
	double lamda[] = {1.0,1.0,1.0,1.0,1.3,1.5,1.6,2.0,2.0,2.0,2.0};
	double mu[] = {0.0,0.0, 0.0, 0.0, 0.4, 0.7, 1.0, 1.0, 1.0, 1.0}; 
	double c[] ={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.nx = 11;
	spline.ny = 10;
	spline.c = c;
	spline.lamda =lamda;
	spline.mu = mu;

	nag_2d_spline_eval(m, x, y, ff, &spline);

	printf(" i x[i] y[i] ff[i]\n");
	for (i=0; i<m; i++)
		printf("%7ld %11.3f%11.3f%11.3f\n",i,x[i],y[i],ff[i]);
}

	The output is following:
	
	
	i 	x[i] 	y[i] 	ff[i]
	
	0 	1.000 	0.000 	1.000
	1 	1.100 	0.100 	1.310
	2 	1.500 	0.700 	2.950
	3 	1.600 	0.400 	2.960
	4 	1.900 	0.300 	3.910
	5 	1.900 	0.800 	4.410
	6 	2.000 	1.000 	5.000

Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, m must not be less than 1: m = _value_.  On entry, spline.nx must not be less than 8: spline.nx = _value_.  On entry, spline.ny must not be less than 8: spline.ny = _value_.  
	73: Memory allocation failed.
	255: On entry, the end knots must satisfy _value_, _value_ = _value_, _value_ = _value_.
	62: The sequence spline.lamda is not increasing: spline.lamda [_value_] = _value_, spline.lamda [_value_] = _value_. The sequence spline.mu is not increasing: spline.mu[_value_] = _value_, spline.mu[_value_] = _value_.
	257: On entry, point (x[_value_] = _value_, y[_value_] = _value_) lies outside the rectangle bounded by spline.lamda[3] = _value_, spline.lamda[_value_] = _value_, spline.mu[3] = _value_, spline.mu[_value_] = _value_.  
		
	successfully call of the nag_2d_spline_eval.	
*/

	int nag_2d_spline_eval(
		int m , // the number of points at which values of the spline are required
		const double x[],
		const double y[], // These are the coordinates of the points at which values of the spline are required
		double ff[], // contains the value of the spline at the point (xr, yr),
		Nag_2dSpline *spline // Pointer to structure of type Nag_2dSpline
  ); // Evaluation of bicubic spline, at a set of points


/** e02dfc
		calculates values of a bicubic spline from its B-spline representation.
		The spline is evaluated at all points on a rectangular grid.
		
Example:
	This program reads in knot sets spline.lamda[0],. . . ,spline.lamda[spline.nx-1]
	and spline.mu[0],. . ., spline.mu[spline.ny-1], and a set of bicubic spline 
	coefficients cij . Following these are values for mx and the x co-ordinates x[q], 
	for q = 1, 2, . . .,mx, and values for my and the y co-ordinates y[r],
	for r = 1, 2, . . .,my, defining the grid of points on which the spline 
	is to be evaluated.
		
void test_nag_2d_spline_eval_rect()
{
	int i, j, mx = 7, my = 6;
	double x[20] = {1.0, 1.1, 1.3, 1.4, 1.5, 1.7, 2.0};
	double y[20] = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0};
	double ff[400];
	Nag_2dSpline spline;
	double lamda[] = {1.0,1.0,1.0,1.0,1.3,1.5,1.6,2.0,2.0,2.0,2.0};
	double mu[] = {0.0,0.0, 0.0, 0.0, 0.4, 0.7, 1.0, 1.0, 1.0, 1.0}; 
	double c[] ={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.nx = 11;
	spline.ny = 10;
	spline.c = c;
	spline.lamda =lamda;
	spline.mu = mu;

	nag_2d_spline_eval_rect(mx, my, x, y, ff, &spline);
	
	printf("Spline evaluated on x-y grid (x across, y down):\n   ");
	for(i=0; i<mx; i++)
		printf(" %8.3f ",x[i]);
	printf("\n");
	for (j=0; j<my; j++)
	{
		printf("%2.1f",y[j]);
		for (i=0; i<mx; i++)
		{
			printf(" %8.3f ", ff[my*i+j]);
			if(i%7==6 || i==mx-1) 
				printf("\n");
		}
	}
}

	The output is following:
		
	Spline evaluated on x-y grid (x across, y down):
	
			1.0 	1.1 	1.3 	1.4 	1.5 	1.7 	2.0
	0.0 	1.000 	1.210 	1.690 	1.960 	2.250 	2.890 	4.000
	0.2 	1.200 	1.410 	1.890 	2.160 	2.450 	3.090 	4.200
	0.4 	1.400 	1.610 	2.090 	2.360 	2.650 	3.290 	4.400
	0.6 	1.600 	1.810 	2.290 	2.560 	2.850 	3.490 	4.600
	0.8 	1.800 	2.010 	2.490 	2.760 	3.050 	3.690 	4.800
	1.0 	2.000 	2.210 	2.690 	2.960 	3.250 	3.890 	5.000
	
Return:
	This function returns NAG error code, 0 if no error.
	
	11: On entry, mx must not be less than 1: mx = _value.  On entry, my must not be less than 1: my = _value.  On entry, spline.nx must not be less than 8: spline.nx = _value.  On entry, spline.ny must not be less than 8: spline.ny = _value.
	73: Memory allocation failed.
	255: On entry, the end knots must satisfy _value	, _value	 = _value	, _value	 = _value	.
	62: The sequence spline.lamda is not increasing: spline.lamda[_value	] = _value	, spline.lamda [_value	] = _value	.  The sequence spline.mu is not increasing: spline.mu[_value	] = _value	, spline.mu[_value	] = _value	.  
	256: On entry, the end knots and coordinates must satisfy spline.lamda[3] = x[0] and x[mx]-1] = spline.lamda[spline.nx-4]. spline.lamda[3] = _value	, x[0] = _value	, x[_value	] = _value	, spline.lamda[_value	] = _value	.  On entry, the end knots and coordinates must satisfy spline.mu[3] = y[0] and y[my-1] = spline.mu[spline.ny-4]. spline.mu[3] = _value	, y[0] = _value	, y[_value	] = _value	, spline.mu[_value	] = _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	.
		
	successfully call of the nag_2d_spline_eval_rect function.	
*/

	int  nag_2d_spline_eval_rect(
		int mx,
		int my, //the number of points along the x and y axes that define the rectangular grid.
		const double x[],
		const double y[], // These are the x and y co-ordinates that define the rectangular grid of points at which values of the spline are required
		double ff[], // contains the value of the spline at the point (xq, yr),
		Nag_2dSpline *spline // Pointer to structure of type Nag_2dSpline
  ); // Evaluation of bicubic spline, at a mesh of points