#include <OC_nag.h> // this contains all the NAG headers, 
void nag_spline(string strYDataset, int nPts)
{
	// Check nPts
	if( nPts < 1 )
	{
		out_str("Number of points invalid!");
		return;
	}
	
	// Declare curve using Y dataset name
	Curve crvData(strYDataset);
	// Create a copy of the curve 
	Curve crvCopy(crvData);
	if( !crvCopy )
	{
		out_str("Failed to create curve from specified dataset name!");
		return;
	}
	
	// Copy Y, X values of curve to vectors
	vector vecX, vecY;
	Dataset dsX;
	crvCopy.AttachX(dsX);
	vecY = crvCopy;
	vecX = dsX;
	
	// Set up vectors for return values from NAG function
	vector vecXSpline(nPts), vecYSpline(nPts);
	// Fill X vector with nPts that cover the X range of the curve
	// Assume X is sorted ascending
	double dMin, dMax, dStep;
	int nDataSize = vecX.GetSize();
	dMin = vecX[0];
	dMax = vecX[nDataSize - 1];
	dStep = (dMax - dMin) / (nPts - 1);
	vecXSpline.Data(dMin, dMax, dStep);

	// First compute spline coefficients 
	Nag_Spline spline;
	int iRet;
	if( (iRet = nag_1d_spline_interpolant(nDataSize, vecX, vecY, &spline)) != 0 )
	{
		printf("NAG computation of spline coefficients returned error: %d\n", iRet);
		return;
	}
	
	
	// Now compute spline curve at each x point
	double fit;
	int iS = vecXSpline.GetSize();
	for(int ii = 0; ii < iS; ii++)
	{
		if( (iRet = nag_1d_spline_evaluate(vecXSpline[ii], &fit, &spline)) != 0 )
		{
			printf("NAG spline computation failed at x = %f with error code: %d\n", vecXSpline[ii], iRet);
			return;
		}
		vecYSpline[ii] = fit;
	}
	
	// Succeeded in computing spline - write results to new wks
	WorksheetPage wpg;
	wpg.Create("Origin");
	Worksheet wks = wpg.Layers(0);
	while( wks.DeleteCol(0) );
	wks.AddCol("SplineX");
	wks.Columns(0).SetType(OKDATAOBJ_DESIGNATION_X);
	wks.AddCol("SplineY");
	string strLabel;
	strLabel.Format("NAG Spline Interpolation for %s", strYDataset);
	wpg.Label = strLabel;
	Dataset dsXSpline(wks, 0);
	Dataset dsYSpline(wks, 1);
	dsXSpline = vecXSpline;
	dsYSpline = vecYSpline;

	// Clean up
	nag_free(spline.lamda);
	nag_free(spline.c);
}