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

#ifndef _O_NAG_S_H
#define _O_NAG_S_H

#importdll "ONAG" // NAG DLL prepared by OriginLab

#include <NAG\nag_types.h>

/* begin proto */
/** s10aac
		returns a value for the hyperbolic tangent, tanh x.

Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_tanh()
{

	double  y;
	double x[] = { -20.0, -5.0,	0.5, 5.0};
	printf(" 	   x       y\n");
	for(int i = 0; i < 4; i++)
	{
		y = nag_tanh(x[i]);
		printf("%12.1f%12.5f\n", x[i], y);
	}
	
}

	The output if following:
	
		
		
	  x 			y
	
	-20.0 		-1.00000
	-5.0 		-0.99991
	 0.5  		 0.46212
	 5.0 		 0.99991
	
Parameters:

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

*/
double	nag_tanh(
		double x // the argument x of the function.
);


/**	s10abc
		returns the value of the hyperbolic sine, sinh x.
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	

void test_nag_sinh()
{
	double  y;
	double x[] = {-10.0, -0.5, 0.0, 0.5, 25.0};
	printf("   x		  y\n");
	for(int i = 0; i < 5; i++)
	{
		y = nag_sinh(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}

	The output is following:
	
		x 			y
	
	-1.000e+01 	-1.101e+04
	-5.000e-01 	-5.211e-01
	 0.000e+00 	 0.000e+00
	 5.000e-01 	 5.211e-01
	 2.500e+01 	 3.600e+10
	 
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	5: On entry, |x| must not be greater than _value_: x = _value_.  The function has been called with an argument too large in absolute magnitude. There is a danger of setting over.ow. The result is the value of sinh at the closest argument for which a valid call could be made. 
		
	successfully call of the nag_sinh function.

*/
double	nag_sinh(
		double x, // the argument x of the function.
		NagError* fail = NULL 
);

/**	s10acc
		returns the value of the hyperbolic cosine, cosh x.
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_cosh()
{
	double  y;
	double x[] = {-10.0, -0.5, 0.0, 0.5, 25.0};
	printf("   x		  y\n");
	for(int i = 0; i < 5; i++)
	{
		y = nag_cosh(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}


	The output is following:
	
		x 			y
	
	-1.000e+01 	1.101e+04
	-5.000e-01 	1.128e+00
	 0.000e+00 	1.000e+00
	 5.000e-01 	1.128e+00
	 2.500e+01 	3.600e+10
	 
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	7: On entry, |x| must not be greater than _value_: x = _value_.  The function has been called with an argument too large in absolute magnitude. There is a danger of over.ow. The result returned is the value of cosh x at the nearest valid argument.
	
	successfully call of the nag_cosh function.
	
*/

double	nag_cosh(
		double x, // the argument x of the function.
		NagError* fail = NULL 
);
/**	s11aac
		returns the value of the inverse hyperbolic tangent, arctanh x.
	
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_arctanh()
{
	double  y;
	double x[] = {-0.5, 0.0, 0.5, -0.9999};
	printf("   x		  y\n");
	for(int i = 0; i < 4; i++)
	{
		y = nag_arctanh(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}	

	

	The output is following:
	
		x 			y
	
	-5.000e-01 	-5.493e-01
	 0.000e+00 	 0.000e+00
	 5.000e-01 	 5.493e-01
	-9.999e-01 	-4.952e+00
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	8: On entry, |x| must not be greater than or equal to 1.0: x = _value_.  The function has been called with an argument greater than or equal to 1.0in magnitude, for which arctanh is not de.ned. The result is returned as zero.
	
	successfully call of the nag_arctanh function.

*/
double	nag_arctanh(
		double x, // the argument x of the function.
		NagError* fail = NULL 
);
/**	s11abc
		returns the value of the inverse hyperbolic sine, arcsinh x.
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_arcsinh()
{
	double  y;
	double x[] = {-2.0, -0.5, 1.0, 6.0};
	printf("   x		  y\n");
	for(int i = 0; i < 4; i++)
	{
		y = nag_arcsinh(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}




	The output is following:
	
		x 			y
	
	-2.000e+00 	-1.444e+00
	-5.000e-01 	-4.812e-01
	 1.000e+00 	 8.814e-01
	 6.000e+00 	 2.492e+00

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
		
	sucessfully call of the nag_arcsinh function.

				
*/

double	nag_arcsinh(
		double x // the argument x of the function.
);
/**	s11acc
		returns the value of the inverse hyperbolic cosine, arccoshx. 
		The result is in the principal positive branch.

Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_arccosh()
{
	double  y;
	double x[] = {1.00, 2.0, 5.0, 10.0};
	printf("   x		  y\n");
	for(int i = 0; i < 4; i++)
	{
		y = nag_arccosh(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}


	The output is following:
	
		x 			y
	
	1.000e+00 	0.000e+00
	2.000e+00 	1.317e+00
	5.000e+00 	2.292e+00
	1.000e+01 	2.993e+00

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	5: On entry, x must not be less than 1.0: x = _value_.  arccosh x is not de.ned and the result returned is zero.
		
	successfully call of the nag_arccosh function.

*/

double	 nag_arccosh(
		 double x, // the argument x of the function.
		 NagError* fail = NULL 
);
/**	s13aac
		returns the value of the exponential integral E1(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_exp_integral()
{
	double  y;
	double x[] = {2.0, 9.0};
	printf("   x		  y\n");
	for(int i = 0; i < 2; i++)
	{
		y = nag_exp_integral(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}


	The output is following:
	
		x 			y
	
	2.000e+00 	4.890e-02
	9.000e+00 	1.245e-05
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	6: On entry, x must not be less than or equal to 0.0: x = _value_.  The function is not de.ned for this value and the result returned is zero.
	
	successfully call of the nag_exp_integral function.

*/

double		nag_exp_integral(
		double x, // the argument x of the function.
		NagError *fail = NULL
);
/**	s13acc
		returns the value of the cosine integral Ci(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_cos_integral()
{
	double  y;
	double x[] = {0.2, 0.4, 0.6, 0.8, 1.0};
	printf("   x		  y\n");
	for(int i = 0; i < 5; i++)
	{
		y = nag_cos_integral(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}

		
		
	The output is following:
	
		x 			y
		
	2.000e-01 	-1.042e+00
	4.000e-01 	-3.788e-01
	6.000e-01 	-2.227e-02
	8.000e-01 	 1.983e-01
	1.000e+00 	 3.374e-01
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	6: On entry, x must not be less than or equal to 0.0: x = _value_.  The function is not de.ned for this value and the result returned is zero.

	successfully call of the nag_cos_integral function.

*/

double		nag_cos_integral(
		double x, // the argument x of the function.
		NagError *fail = NULL
);

/**	s13adc
		returns the value of the sine integral Si(x).

Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_sin_integral()
{
	double  y;
	double x[] = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0};
	printf("   x		  y\n");
	for(int i = 0; i < 6; i++)
	{
		y = nag_sin_integral(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}

	
	
	
	The output is following:
	
		x 			y
	
	0.000e+00 	0.000e+00
	2.000e-01 	1.996e-01
	4.000e-01 	3.965e-01
	6.000e-01 	5.881e-01
	8.000e-01 	7.721e-01
	1.000e+00 	9.461e-01

Parameters:

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

double		nag_sin_integral(
		double x // the argument x of the function.
);

/**	s14aac
		returns the value of the Gamma function G(x).

Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_gamma()
{
	double  y;
	double x[] = {1.0, 1.25, 1.5, 1.75, 2.0, 5.0, 10.0, -1.5};
	printf("   x		  y\n");
	for(int i = 0; i < 8; i++)
	{
		y = nag_gamma(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}


	The output is following:
	
		x 				y
		
	1.000e+00 		1.000e+00
	1.250e+00 		9.064e-01
	1.500e+00 		8.862e-01
	1.750e+00 		9.191e-01
	2.000e+00 		1.000e+00
	5.000e+00 		2.400e+01
	1.000e+01 		3.629e+05
	-1.500e+00 		2.363e+00

Parameters:

Return:
	This funtion returns NAG error code, 0 if no error.
	
	7: On entry, x must not be greater than _value_: x = _value_.  The argument is too large, the function returns the approximate value of G(x) at the nearest valid argument.
	5: On entry, x must not be less than _value_: x = _value_.  The argument is too large and negative, the function returns zero.
	692: On entry, x must be greater than _value_: x = _value_.  The argument is too close to zero, the function returns the approximate value of G(x) at the nearest valid argument.
	693: On entry, x must not be e.ectively a negative integer: x = _value_.  The argument is a negative integer, at which values G(x) is in.nite. The function returns a large positive value.
	
	successfully call of the nag_gamma function.
*/

double		nag_gamma(
		double x, // the argument x of the function.
		NagError *fail = NULL
);

/**	s14abc
		returns a value for the logarithm of the Gamma function, ln G(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_log_gamma()
{
	double  y;
	double x[] = {1.0, 1.25, 1.5, 1.75, 2.0, 5.0, 10.0, 20.0, 1000.0};
	printf("   x		  y\n");
	for(int i = 0; i < 9; i++)
	{
		y = nag_log_gamma(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}




	The output is following:
	
		x 			y
		
	1.000e+00 	 0.000e+00
	1.250e+00 	-9.827e-02
	1.500e+00 	-1.208e-01
	1.750e+00 	-8.440e-02
	2.000e+00 	 0.000e+00
	5.000e+00 	 3.178e+00
	1.000e+01    1.280e+01
	2.000e+01 	 3.934e+01
	1.000e+03 	 5.905e+03
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	6: On entry, x must not be less than or equal to 0.0: x = _value_.  For zero and negative values the function is undefined, the function returns zero.
	7: On entry, x must not be greater than _value_: x = _value_.  x is too large, the function would over.ow, the function returns the value of the function at the largest permissible argument.
	
	successfully call of the nag_log_gamma function.
	
*/

double		nag_log_gamma(
		double x, // the argument x of the function.
		NagError *fail = NULL
);

/**	s14aec
		returns the value of the kth derivative of the psi function psi(x) for real x and k = 0, 1, ..., 6.
		
Example:
	The example program evaluates \psi(2)(x) at  x = 2.5, and prints the results

void test_nag_real_polygamma()
{
	double y;
	double x[] ={1.0, 0.5, -3.6, 8.0, 2.9, -4.7, -5.4};
	int k[] = {0, 1, 2, 3, 4, 5, 6};
	printf("  x	    k	    (D^K/DX^K)psi(X)\n");
	for( int i = 0; i < 6; i++)
	{
		y = nag_real_polygamma(x[i], k[i]);
		printf("%5.1f %5ld     %12.4e\n", x[i], k[i], y);
	}	
}


The output is following:

 X		    K		(D^K/DX^K)psi(X)
 
 1.0		0		-5.772 e-01
 0.5		1		 4.9348e+00
-3.6		2		-2.2335e+01
 8.0		3		 4.6992e-03
 2.9		4		-1.5897e-01
-4.7		5		 1.6566e+05
-5.4		6		 4.1378e+05

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	90: On entry, k = <value>.  Constraint: 0 <= k <= 6.
	96: On entry, x = <value>.  Constraint: x must not be 'too close' to a non-positive integer.  That is, |x - nint(x)| >= machine precision x nint(x).
	694: The evaluation has been abandoned due to the likelihood of underflow.  The result is returned as zero.
	695: The evaluation has been abandoned due to the likelihood of overflow.  The result is returned as zero.
	74:  An internal error has occurred in this function.  Check the function call and any array sizes.  If the call is correct the please consult NAG for assistance.
		
	successfully call of the nag_real_polygamma.		

*/

double		nag_real_polygamma(
		double x, // the argument x of the function.
		int k, 	  // the function psi-k(z) to be evaluated.
		NagError *fail = NULL
);

/**	s14afc
		returns the value of the kth derivative of the psi function psi(z), for complex z and k = 0, 1, ..., 4.
		
Example:
	The example program evaluates the psi(trigamma) function at z = -1.5 +2.5i, 
	and prints the results.

void test_nag_complex_polygamma()
{
	complex z, z_l;
	double re[5] = {1.2, 0.5, -1.5, 8.0, 2.9};
	double im[5] = {5.0, -0.2, 2.5, 3.3, 7.5};
	int k[5] = {0, 1, 1, 3, 4};
	printf("          z        k          (d^K/DZ^k)psi(z)\n\n");

	for(int i = 0; i < 5; i++)
	{
		z.m_re = re[i];
		z.m_im = im[i];
		z_l = nag_complex_polygamma(z, k[i]);
		printf("(%5.1f, %5.1f) %6ld (%12.4e, %12.4e) \n",z.m_re, z.m_im, k[i], z_l.m_re, z_l.m_im);
	}
}


	The output is following:
	
	z 					k				(D^K/DZ^K)psi(z)
	
	( 1.2,  5.0)		0			( 1.6176e+00,  1.4312e+00)
	( 0.5, -0.2)		1			( 3.4044e+00,  2.5394e+00)
	(-1.5,  2.5)		2			(-1.9737e-01, -2.4271e-01)
	( 8.0,  3.3)		3			( 1.1814e-03, -3.4188e-03)
	( 2.9,  7.5)		4			(-5.0227e-04, -1.4955e-03)
	    
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	90: On entry, k = <value>.  Constraint: 0 <= k <= 4.	    
	102: On entry, z = (<value>, <value>).  Constraint: z.re must not be ' too close' to a non-positive integer when z.im = 0.0.  That is, |z.re - nint(z.re)| >= machine precision x |nint(z.re)|
	695: The evaluation has been abandoned due to the likelihood of overflow.  The result is returned as zero.
	74: An internal error has occurred in the function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.

	successfully call of the nag_complex_polygamma function.	
		
*/	

complex  nag_complex_polygamma(
		complex z, // the argument z of the function. 
		int k, // the function psi-k(z) to be evaluated.
		NagError *fail = NULL
);

/**	s14bac
		computes values for the incomplete gamma functions P(a, x) and Q(a, x).
		
Example:
	The following program reads values of the argument a and x from a file, 
	evaluates the function and prints the results.


void test_nag_incomplete_gamma()
{
	double q, tol, p;
	double a[] = {2.0,7.0, 0.5, 20.0, 21.0};
	double x[] = {3.0, 1.0, 99.0, 21.0, 20.0}
	tol = 1e-5;
	printf("       a            x           p           q\n");
	
	for(int i = 0; i < 5; i++)
	{
		nag_incomplete_gamma(a[i], x[i], tol, &p, &q);
		printf("%12.4f%12.4f%12.4f%12.4f\n", a[i], x[i], p, q);
	}
}



	The output is following:
			
		a 		x 		p 	 q
	
	 2.0000	 3.0000 0.8009 0.1991
	 7.0000	 1.0000 0.0001 0.9999
	 0.5000	99.0000 1.0000 0.0000
	20.0000 21.0000 0.6157 0.3843
	21.0000 20.0000 0.4409 0.5591

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	On error nag incomplete gamma returns with a value of 0.0 for p and q.
	6: On entry, a must not be less than or equal to 0.0: a = _value_.
	11: On entry, x must not be less than 0.0: x = _value_.
	405: The algorithm has failed to converge in _value_ iterations.  Convergence of the Taylor series or Legendre continued fraction has failed within the speci.ed number of iterations. This error is extremely unlikely to occur; if it does, contact NAG.
	
	successfully call of the nag_incomplete_gamma function.

*/

int	nag_incomplete_gamma(
	double a, // the argument a of the functions.
	double x, // the argument x of the functions.
	double tol, // the relative accuracy required by the user in the results. If nag incomplete gamma is entered with tol greater than 1.0 or less than machine precision, then the value of machine precision is used instead.
	double *p, // the values of the functions P(a, x) and Q(a, x) respectively.
	double *q  // the values of the functions P(a, x) and Q(a, x) respectively.
);

/**	s15abc
		returns the value of the cumulative Normal distribution function P(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_cumul_normal()
{
	double y;
	double x[] = { -20.0, -1.0, 0.0, 1.0, 2.0, 20.0};

	printf("   x		  y\n");
	for(int i = 0; i < 6; i++)
	{
		y= nag_cumul_normal(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}



	The output is following:
	
		x 			y
		
	-2.000e+01 	2.754e-89
	-1.000e+00 	1.587e-01
	 0.000e+00 	5.000e-01
	 1.000e+00 	8.413e-01
	 2.000e+00 	9.772e-01
	 2.000e+01 	1.000e+00
	 
Parameters:

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

*/

double 	nag_cumul_normal(
		double x // the argument x of the function.
);

/**	s15acc
		returns the value of the complement of the cumulative normal distribution function Q(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_cumul_normal_complem()
{
	double y;
	double x[] = { -20.0, -1.0, 0.0, 1.0, 2.0, 20.0};
	
	printf("   x		  y\n");
	for(int i = 0; i < 6; i++)
	{
		y= nag_cumul_normal_complem(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}

}

	
	

	The output is following:
	
	x 				y
	
	-2.000e+01 	1.000e+00
	-1.000e+00 	8.413e-01
	 0.000e+00 	5.000e-01
	 1.000e+00 	1.587e-01
	 2.000e+00 	2.275e-02
	 2.000e+01 	2.754e-89
	 
Parameters:

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

*/

double		nag_cumul_normal_complem(
		double x // the argument x of the function.
);

/**	s15adc
		returns the value of the complementary error function, erfc x.
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_erfc()
{
	double  y;
	double x[] = { -10.0, -1.0, 0.0, 1.0, 10.0}
	printf("   x		  y\n");
	for(int i = 0; i < 5; i++)
	{
		y= nag_erfc(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}

}

	The output is following:
	
	x 					y
	
	-1.000e+01 		2.000e+00
	-1.000e+00 		1.843e+00
	 0.000e+00 		1.000e+00
	 1.000e+00 		1.573e-01
	 1.000e+01 		2.088e-45
	 
Parameters:

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

*/
double	nag_erfc(
		double x // the argument x of the function.
);

/**	s15aec
		returns the value of the error function, erf x.
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_erf()
{
	double y;
	double x[] = {-6.0, -4.5, -1.0, 1.0, 4.5, 6.0};
	printf("   x		  y\n");
	for(int i = 0; i < 6; i++)
	{
		y= nag_erf(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
	
}



	The output is following:
	
	    x			 y
	
	-6.000e+00 	-1.000e+00
	-4.500e+00 	-1.000e+00
	-1.000e+00 	-8.427e-01
	 1.000e+00 	 8.427e-01
	 4.500e+00 	 1.000e+00
	 6.000e+00 	 1.000e+00

Parameters:

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

*/

	double		nag_erf(
	double x // the argument x of the function.
	);

/**	s17acc
		returns the value of the Bessel function Y0(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_bessel_y0()
{
	double y;
	double x[] = {0.5, 1.0, 3.0, 6.0, 8.0, 10.0, 1000.0};
	printf("   x		  y\n");
	for(int i = 0; i < 7; i++)
	{
		y= nag_bessel_y0(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}


	The output is following:
	
		x 				y
	5.000e-01 		-4.445e-01
	1.000e+00 		 8.826e-02
	3.000e+00 		 3.769e-01
	6.000e+00 		-2.882e-01
	8.000e+00 		 2.235e-01
	1.000e+01 		 5.567e-02
	1.000e+03 		 4.716e-03

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	7: On entry, x must not be greater than _value_: x = _value_.  x is too large, the function returns the amplitude of the Y0 oscillation, _2/px.
	6: On entry, x must not be less than or equal to 0.0: x = _value_.  Y0 is undefined, the function returns zero.

	successfully call of the nag_bessel_y0 function.
*/
double	nag_bessel_y0(
		double x, 	//the argument x of the function.
		NagError *fail = NULL
);

/**	s17adc
		returns the value of the Bessel function Y1(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_bessel_y1()
{
	double y;
	double x[] = {0.5, 1.0, 3.0, 6.0, 8.0, 10.0, 1000.0};
	printf("   x		  y\n");
	for(int i = 0; i < 7; i++)
	{
		y= nag_bessel_y1(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}



	The output is following:
	
	x 				y

	5.000e-01 	-1.471e+00
	1.000e+00 	-7.812e-01
	3.000e+00 	 3.247e-01
	6.000e+00 	-1.750e-01
	8.000e+00 	-1.581e-01
	1.000e+01 	 2.490e-01
	1.000e+03 	-2.478e-02
	
Parameters:

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

	7: On entry, x must not be greater than _value_: x = _value_.  x is too large, the function returns the amplitude of the Y1 oscillation, _2/px.
	6: On entry, x must not be less than or equal to 0.0: x = _value_.  Y1 is undefined, the function returns zero.
	692: On entry, x must be greater than _value_: x = _value_.  x is too close to zero, there is a danger of over.ow, the function returns the value of Y1(x) at the smallest valid argument.
		
	successfully call of the nag_bessel_y1 function.

*/
double nag_bessel_y1(
	   double x, // the argument x of the function.
	   NagError *fail = NULL
);

/**	s17aec
		returns the value of the Bessel function J0(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results

void test_nag_bessel_j0()
{
	double y;
	double x[] ={0.0, 0.5, 1.0, 3.0, 6.0, 8.0, 10.0, -1.0, 1000.0}
	
	printf("   x		  y\n");
	for(int i = 0; i < 9; i++)
	{
		y= nag_bessel_j0(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}




	The output is following:
	
	x 					y
	
	0.000e+00 		1.000e+00
	5.000e-01 		9.385e-01
	1.000e+00 		7.652e-01
	3.000e+00 	   -2.601e-01
	6.000e+00 		1.506e-01
	8.000e+00 		1.717e-01
	1.000e+01 	   -2.459e-01
	-1.000e+00 		7.652e-01
	1.000e+03 		2.479e-02
	
Parameters:

Return:
	This function return NAG error code, 0 if no error.
	
	7: On entry, x must not be greater than _value_: x = _value_.  x is too large, the function returns the amplitude of the J0 oscillation, _2/p|x|.
		
	successfully call of the nag_bessel_j0 function.

*/
double 	nag_bessel_j0(
	   double x, // the argument x of the function.
	   NagError *fail = NULL
);

/**	s17afc
		returns the value of the Bessel function J1(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_bessel_j1()
{
	double y;
	double x[] ={0.0, 0.5, 1.0, 3.0, 6.0, 8.0, 10.0, -1.0, 1000.0}
	printf("   x		  y\n");
	for(int i = 0; i < 9; i++)
	{
		y= nag_bessel_j1(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}




	The output is following:
	
		x 			y
		
	0.000e+00 	0.000e+00
	5.000e-01 	2.423e-01
	1.000e+00 	4.401e-01
	3.000e+00 	3.391e-01
	6.000e+00  -2.767e-01
	8.000e+00	2.346e-01
	1.000e+01 	4.347e-02
	-1.000e+00 -4.401e-01
	1.000e+03 	4.728e-03
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	7: On entry, x must not be greater than _value_: x = _value_.  x is too large. The function returns the amplitude of the J1 oscillation, _2/p|x|.

	successfully call of the 	nag_bessel_j1.
	
*/
double 	nag_bessel_j1(
	   double x, // the argument x of the function.
	   NagError *fail = NULL
);

/**	s17agc
		returns a value for the Airy function Ai(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_airy_ai()
{
	double y;
	double x[] ={-10.0, -1.0, 0.0, 1.0, 5.0, 10.0, 20.0};
	
	printf("   x		  y\n");
	for(int i = 0; i < 7; i++)
	{
		y= nag_airy_ai(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}



	The output is following:
	
		x 				y
		
	-1.000e+01 		4.024e-02
	-1.000e+00 		5.356e-01
	 0.000e+00 		3.550e-01
	 1.000e+00 		1.353e-01
	 5.000e+00 		1.083e-04
	 1.000e+01 		1.105e-10
	 2.000e+01 	  	1.692e-27
	 
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	7: On entry, x must not be greater than _value_: x = _value_.  x is too large and positive. The function returns zero.
	5: On entry, x must not be less than _value_: x = _value_.  x is too large and negative. The function returns zero.
		
	successfully call of the nag_airy_bi function.

*/
double 	nag_airy_ai(
	   double x, // the argument x of the function.
	   NagError *fail = NULL
);

/**	s17ahc
		returns a value of the Airy function Bi(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_airy_bi()
{
	double y;
	double x[] ={-10.0, -1.0, 0.0, 1.0, 5.0, 10.0, 20.0};
	
	printf("   x		  y\n");
	for(int i = 0; i < 7; i++)
	{
		y= nag_airy_bi(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}




	The output is following:
	
		x 			y
	
	-1.000e+01 	-3.147e-01
	-1.000e+00 	 1.040e-01
	 0.000e+00   6.149e-01
	 1.000e+00 	 1.207e+00
	 5.000e+00 	 6.578e+02
	 1.000e+01 	 4.556e+08
	 2.000e+01 	 2.104e+25
	 
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	7: On entry, x must not be greater than _value_: x = _value_.  x is too large and positive. The function returns zero.
	5: On entry, x must not be less than _value_: x = _value_.  x is too large and negative. The function returns zero.

	successfully call of the nag_airy_bi function.

*/
double 	nag_airy_bi(
	   double x, // the argument x of the function.
	   NagError *fail = NULL
);

/**	s17ajc
		returns a value of the derivative of the Airy function Ai(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_airy_ai_deriv()
{	
	double y;
	double x[] ={-10.0, -1.0, 0.0, 1.0, 5.0, 10.0, 20.0};
		
	printf("   x		  y\n");
	for(int i = 0; i < 7; i++)
	{
		y= nag_airy_ai_deriv(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}


	The output is following:
	
		x 			y
		
	-1.000e+01 	 9.963e-01
	-1.000e+00 	-1.016e-02
	 0.000e+00 	-2.588e-01
	 1.000e+00 	-1.591e-01
	 5.000e+00 	-2.474e-04
	 1.000e+01 	-3.521e-10
	 2.000e+01 	-7.586e-27
	 
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	7: On entry, x must not be greater than _value_: x = _value_.  x is too large and positive. The function returns zero.
	5: On entry, x must not be less than _value_: x = _value_.  x is too large and negative. The function returns zero.
	
	successfully call of the nag_airy_ai_deriv function.

*/
double 	nag_airy_ai_deriv(
	   double x, // the argument x of the function.
	   NagError *fail = NULL
);

/**	s17akc
		returns a value of the derivative of the Airy function Bi(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_bi_deriv()
{
	double y;
	double x[] ={-10.0, -1.0, 0.0, 1.0, 5.0, 10.0, 20.0};
		
	printf("   x		  y\n");
	for(int i = 0; i < 7; i++)
	{
		y= nag_airy_bi_deriv(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}		
}



	The output is following:
	
		x 			y
	
	-1.000e+01 	1.194e-01
	-1.000e+00 	5.924e-01
	 0.000e+00  4.483e-01
	 1.000e+00  9.324e-01
	 5.000e+00  1.436e+03
	 1.000e+01  1.429e+09
	 2.000e+01  9.382e+25
	 
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	7: On entry, x must not be greater than _value_: x = _value_.  x is too large and positive. The function returns zero.
	5: On entry, x must not be less than _value_: x = _value_.  x is too large and negative. The function returns zero.

	successfully call of the nag_airy_bi_deriv function.

*/
double 	nag_airy_bi_deriv(
	   double x, // the argument x of the function.
	   NagError *fail = NULL
);

/**	s17alc
		determines the leading N zeros of one of the Bessel functions Ja(x), 
		Ya(x), J'a(x), Y'a(x) for real x and non-nagative a (alpha).

Example:

	to determine the leading five positive zeros of the Bessel function j0(x)

void test_nag_bessel_zeros()
{

	double a = 0.0;
	int n = 5;
	int mode = 1;
	double rel = 1.1e-8;
	double x[100];
	nag_bessel_zeros(a, n, mode, rel, x);
	printf(" a,       n,   mode,      rel\n");
	printf("                     (machine-dependent)\n");
	printf("%4.1f   %3ld %6ld        %9.1e\n\n", a, n, mode, rel);

	if(mode == 1)
		printf("Leading N positive zeros of J\n");
	if(mode == 2)
		printf("Leading N positive zeros of Y \n");
	if(mode == 3)
	{
		if(a == 0)
			printf("Leading N non_negative zeros of J'\n");
		else 
			printf("Leading N positive zeros of J'\n");
	}
	if(mode == 4)
		printf("Leading N positive zeros of Y'\n\n");
	
	for(int i = 0; i <= n - 1; i++)
		printf("x = %12.4e\n", x[i]);
	printf("\n");		

}	



	The output is following:
	
	a 		n		mode			rel
								(machine-dependent)
	0.0		5		1				1.1e-08
	
	Leading N positive zeros of J
	
	x = 2.4048 e+00
	x = 5.5201e+00
	x = 8.6537e+00
	x = 1.1792e+01
	x = 1.4931e+01			

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	90: On entry, n = <value>.  Constraint: n >= 1.  On entry, mode = <value>.  Constraint: 1<= mode <= 4.
	96: On entry, a = <value>.  Constraint: 0.0 <= a <= 100000.0.  On entry, rel = <value>.  Constraint: rel > 0.0
	74: An internal error has occurred in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	
	successfully call of the nag_bessel_zeros function.
							
*/
int 	nag_bessel_zeros(
	double a, // the order alpha of fuction.
	int n, 	  // the number n of zeros required.
	int mode, // specifies the form of the function whose zeros are required. 
	double rel, // the relative accuracy to which the zeros are required.
	double x[], // the n required zeros of the function specified by mode.
	NagError *fail = NULL
);

/**	s18acc
		returns the value of the modified Bessel function K0(x).

Example:
	The following program reads values of the argument x from a .le, evaluates the function at each
	value of x and prints the results.

void test_nag_bessel_k0()
{
	double y;
	double x[] ={0.4, 0.6, 1.4, 1.6, 2.5, 3.5, 6.0, 8.0, 10.0, 1000.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 10; i++)
	{
		y=nag_bessel_k0(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}


	The output is following:
	
		x 			y
		
	4.000e-01 	1.115e+00
	6.000e-01 	7.775e-01
	1.400e+00 	2.437e-01
	1.600e+00 	1.880e-01
	2.500e+00 	6.235e-02
	3.500e+00 	1.960e-02
	6.000e+00 	1.244e-03
	8.000e+00 	1.465e-04
	1.000e+01 	1.778e-05
	1.000e+03 	0.000e+00
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	6: On entry, x must not be less than or equal to 0.0: x = _value_.  K0 is undefined and the function returns zero.

	successfully call of the nag_bessel_k0 function.


*/
double 	nag_bessel_k0(
	double x, // the argument x of the function.
	NagError *fail = NULL
);

/**	s18adc
		returns the value of the modified Bessel function, K1(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_bessel_k1()
{
	double y;
	double x[] ={0.4, 0.6, 1.4, 1.6, 2.5, 3.5, 6.0, 8.0, 10.0, 1000.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 10; i++)
	{
		y=nag_bessel_k1(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}

	The output is following:
	
		x 			y
		
	4.000e-01 	2.184e+00
	6.000e-01 	1.303e+00
	1.400e+00 	3.208e-01
	1.600e+00 	2.406e-01
	2.500e+00 	7.389e-02
	3.500e+00 	2.224e-02
	6.000e+00 	1.344e-03
	8.000e+00 	1.554e-04
	1.000e+01 	1.865e-05
	1.000e+03 	0.000e+00

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	6: On entry, x must not be less than or equal to 0.0: x = _value_.  K0 is undefined and the function returns zero.
	692: On entry, x must be greater than _value_: x = _value_.  x is too small, there is a danger of over.ow and the function returns approximately the largest representable value.
	
	successfully call of the the nag_bessel_k1 function.
*/
double		nag_bessel_k1(
		double x, // the argument x of the function.
		NagError *fail = NULL
);

/**	s18aec
		returns the value of the modified Bessel function I0(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.

void test_nag_bessel_i0()
{
	double y;
	double x[] ={0.0, 0.5, 1.0, 3.0, 6.0, 8.0, 10.0, 15.0, 20.0, -1.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 10; i++)
	{
		y=nag_bessel_i0(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}


	The output is following:
	
		x 			y
		
	0.000e+00 	1.000e+00
	5.000e-01 	1.063e+00
	1.000e+00 	1.266e+00
	3.000e+00 	4.881e+00
	6.000e+00 	6.723e+01
	8.000e+00 	4.276e+02
	1.000e+01 	2.816e+03
	1.500e+01 	3.396e+05
	2.000e+01 	4.356e+07
	-1.000e+00 	1.266e+00
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	7: On entry, |x| must not be greater than _value_: x = _value_.  |x| is too large and the function returns the approximate value of I0(x) at the nearest valid argument.
		
	successfully call of the nag_bessel_i1 function.

*/

double 	nag_bessel_i0(
	   double x, // the argument x of the function.
	   NagError *fail = NULL
);

/** s18afc
		returns a value for the modified Bessel function I1(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_bessel_i1()
{
	double y;
	double x[] ={0.0, 0.5, 1.0, 3.0, 6.0, 8.0, 10.0, 15.0, 20.0, -1.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 10; i++)
	{
		y=nag_bessel_i1(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}		
}


	The output is following:
	
		x 			y
		
	0.000e+00 	0.000e+00
	5.000e-01 	2.579e-01
	1.000e+00 	5.652e-01
	3.000e+00 	3.953e+00
	6.000e+00 	6.134e+01
	8.000e+00 	3.999e+02
	1.000e+01 	2.671e+03
	1.500e+01 	3.281e+05
	2.000e+01 	4.245e+07
	-1.000e+00 -5.652e-01

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	7: On entry, |x| must not be greater than _value_: x = _value_.  |x| is too large and the function returns the approximate value of I0(x) at the nearest valid argument.
	
	successfully call of the nag_bessel_i1 function.

*/
double		nag_bessel_i1(
		double x, // the argument x of the function.
		NagError *fail = NULL
);

/**	s18ccc
		returns a value of the scaled modified Bessel function exK0(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_bessel_k0_scaled()
{
	double y;
	double x[] ={0.4, 0.6, 1.4, 2.5, 10.0, 1000.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 6; i++)
	{
		y=nag_bessel_k0_scaled(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}		
}


	The output is following:
	
		x 		y
	
	4.000e-01 1.663e+00
	6.000e-01 1.417e+00
	1.400e+00 9.881e-01
	2.500e+00 7.595e-01
	1.000e+01 3.916e-01
	1.000e+03 3.963e-02
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	6: On entry, x must not be less than or equal to 0.0: x = _value_.  K0 is undefined and the function returns zero.

	successfully call of the nag_bessel_k0_scaled function.

*/
double 	nag_bessel_k0_scaled(
	   double x, // the argument x of the function.
	   NagError *fail = NULL
);

/**	s18cdc
		returns a value of the scaled modified Bessel function exK1(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_bessel_k1_scaled()
{
	double y;
	double x[] ={0.4, 0.6, 1.4, 2.5, 10.0, 1000.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 6; i++)
	{
		y=nag_bessel_k1_scaled(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}		
}


	The output is following:
	
		x 		y
		
	4.000e-01 3.259e+00
	6.000e-01 2.374e+00
	1.400e+00 1.301e+00
	2.500e+00 9.002e-01
	1.000e+01 4.108e-01
	1.000e+03 3.965e-02

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	6: On entry, x must not be less than or equal to 0.0: x = _value_.  K1 is undefined and the function returns zero.
	692: On entry, x must be greater than _value_ : x = _value_.  The function returns the value of the function at the smallest permitted value of the argument.

	successfully call of the nag_bessel_k1_scaled.
	
*/
double	 nag_bessel_k1_scaled(
		double x, // the argument x of the function.
		NagError *fail = NULL
);

/**	s18cec
		returns a value of the scaled modified Bessel function e-|x|I0(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_bessel_i0_scaled()
{
	double y;
	double x[] ={0.0, 0.5, 1.0, 3.0, 6.0, 10.0, 1000.0, -1.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 8; i++)
	{
		y=nag_bessel_i0_scaled(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}


	The output is following:
	
		x 			y
		
	0.000e+00 	1.000e+00
	5.000e-01 	6.450e-01
	1.000e+00 	4.658e-01
	3.000e+00 	2.430e-01
	6.000e+00 	1.667e-01
	1.000e+01 	1.278e-01
	1.000e+03 	1.262e-02
	-1.000e+00 	4.658e-01

Parameters:

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

*/
double		nag_bessel_i0_scaled(
		double x // the argument x of the function.
);

/**	s18cfc
		returns a value of the scaled modified Bessel function e-|x|I1(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_bessel_i1_scaled()
{
	double y;
	double x[] ={0.0, 0.5, 1.0, 3.0, 6.0, 10.0, 1000.0, -1.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 8; i++)
	{
		y=nag_bessel_i1_scaled(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}	
}


	The output is following:
	
		x			y
		
	0.000e+00 	1.000e+00
	5.000e-01 	1.564e-01
	1.000e+00 	2.079e-01
	3.000e+00 	1.968e-01
	6.000e+00 	1.521e-01
	1.000e+01 	1.213e-01
	1.000e+03 	1.261e-02
	-1.000e+00 	-2.079e-01

Parameters:

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

*/

double 	nag_bessel_i1_scaled(
	   double x // the argument x of the function.
);

/**	s18ecc
		returns the value of the scaled modified Bessel function (e^-x)I_nu/4(x).

Example:
	The example program reads values of the arguments x and v from a file, 
	evaluates the function and prints the results
	
void test_nag_bessel_i_nu_scaled()
{
	double x[] = {3.9, 1.4, 8.2, 6.7, 0.5, 2.3};
	int nu[] = {-3, -2, -1, 1, 2, 3};
	printf("   x	   nu	 y\n");
	for(int i = 0; i < 6; i++)
	{
		double y = nag_bessel_i_nu_scaled(x[i], nu[i]);
		printf("%4.1f %6ld %12.4e\n", x[i], nu[i],  y);
	}	
}


	The output is following:
	
	x		nu			y
	
	3.9		-3		1.9272e-01
	1.4		-2		3.5767e-01
	8.2		-1		1.4103e-01
	6.7		 1		1.5649e-01
	0.5		 2		3.5664e-01
	2.3		 3		2.3748e-01


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	118: On entry, x = <value>, nu = <value>.  Constraint: x > 0.0 when nu < 0.  On entry, x = <value>, nu = <value>.  Constraint: x >= 0 when nu > 0 .
	90: On entry, nu = <value>.  Constraint: 1<= abs(nu) <=3.
	695: The evaluation has been abandoned due to the likelihood of overflow.  The result is returned as zero.
	696: The evaluation has been completed but some precision has been lost.
	697: The evaluation has been abandoned due to total loss of precision.  The result is returned as zero.
	698: The evaluation has been abandoned due to failure to satisfy the termination condition.  The result is returned as zero.
	74: An internal error has occured in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	
	successfully call of the nag_bessel_i_nu_scaled function.			


*/

double		nag_bessel_i_nu_scaled(
		double x, // the argument x of the function.
		int nu,  // the argument nu of the function.
		NagError *fail = NULL
);

/**	s18edc
		returns the value of the scaled modified Bessel function (e^x)K_nu/4(x) for real x > 0.

Example:
	The example program reads values of the arguments x and v from a file, 
	evaluates the function and prints the results.
	 
void test_nag_bessel_k_nu_scaled()
{
	double x[] = {3.9, 1.4, 8.2, 6.7, 0.5, 2.3};
	int nu[] = {-3, -2, -1, 1, 2, 3};
	printf("   x	   nu	 y\n");
	for(int i = 0; i < 6; i++)
	{
		double y = nag_bessel_k_nu_scaled(x[i], nu[i]);
		printf("%4.1f %6ld %12.4e\n", x[i], nu[i],  y);
	}	
}


	The output is following:
	
	x		nu			y
		
	3.9		-3		6.5781e-01
	1.4		-2		1.0592e+00
	8.2		-1		4.3297e-01
	6.7		 1		4.7791e-01
	0.5		 2		1.7725e+00
	2.3		 3		8.7497e-01


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	96: On entry, nu = <value>.  Constraint: 1 <= abs(nu) <= 3.
	90: On entry, nu = <value>.  Constraint: 1<= abs(nu) <=3.
	695: The evaluation has been abandoned due to the likelihood of overflow.  The result is returned as zero.
	696: The evaluation has been completed but some precision has been lost.
	697: The evaluation has been abandoned due to total loss of precision.  The result is returned as zero.
	698: The evaluation has been abandoned due to failure to satisfy the termination condition.  The result is returned as zero.
	74: An internal error has occured in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
		
	successfully call of the nag_bessel_k_nu_scaled function.
*/

double		nag_bessel_k_nu_scaled(
		double x, // the argument x of the function.
		int nu,  // the argument nu of the function.
		NagError *fail = NULL
);

/**	s18eec
		returns the value of the modified Bessel function I_nu/4(x) for real x > 0.

Example:
	The example program reads values of the arguments x and v from a file, 
	evaluats the function and prints the results.
	
void test_nag_bessel_i_nu()
{
	double x[] = {3.9, 1.4, 8.2, 6.7, 0.5, 2.3};
	int nu[] = {-3, -2, -1, 1, 2, 3};
	printf("   x	   nu	 y\n");
	for(int i = 0; i < 6; i++)
	{
		double y = nag_bessel_i_nu(x[i], nu[i]);
		printf("%4.1f %6ld %12.4e\n", x[i], nu[i],  y);
	}	
}


	The output is following:

	x		nu			y
		
	3.9		-3		9.5207e+00
	1.4		-2		1.4504e+00
	8.2		-1		5.1349e+02
	6.7		 1		1.2714e+02
	0.5		 2		5.8799e-01
	2.3		 3		2.3687e+00	


Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	118: On entry, x = <value>, nu = <value>.  Constraint: x > 0 when nu < 0.  On entry, x = <value>, nu = <value>.  Constraint: 1 <= abs(nu) <= 3.
	90: On entry, nu = <value>.  Constraint: 1<= abs(nu) <=3.
	695: The evaluation has been abandoned due to the likelihood of overflow.  The result is returned as zero.
	696: The evaluation has been completed but some precision has been lost.
	697: The evaluation has been abandoned due to total loss of precision.  The result is returned as zero.
	698: The evaluation has been abandoned due to failure to satisfy the termination condition.  The result is returned as zero.
	74: An internal error has occured in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	
	successfully call of the nag_bessel_i_nu function.
	

*/

double 	nag_bessel_i_nu(
		double x, // the argument x of the function.
		int nu,  // the argument nu of the function.
		NagError *fail = NULL
);

/**	s18efc
		returns the value of the modified Bessel function K_nu/4(x) for real x > 0.
	
Example:
	The example program reads values of the arguments x and v from a file, 
	evaluates the function and prints the results.
	
void test_nag_bessel_k_nu()
{
	double x[] = {3.9, 1.4, 8.2, 6.7, 0.5, 2.3};
	int nu[] = {-3, -2, -1, 1, 2, 3};
	printf("   x	   nu	 y\n");
	for(int i = 0; i < 6; i++)
	{
		double y = nag_bessel_k_nu(x[i], nu[i]);
		printf("%4.1f %6ld %12.4e\n", x[i], nu[i],  y);
	}	
}


	The output is following:

	x		nu			y
		
	3.9		-3		1.3315e-02
	1.4		-2		2.6121e-01
	8.2		-1		1.1892e-04
	6.7		 1		5.8826e-04
	0.5		 2		1.0750e+00
	2.3		 3		8.7724e-02	

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	96: On entry, nu = <value>.  Constraint: 1 <= abs(nu) <= 3.
	90: On entry, nu = <value>.  Constraint: 1<= abs(nu) <=3.
	695: The evaluation has been abandoned due to the likelihood of overflow.  The result is returned as zero.
	696: The evaluation has been completed but some precision has been lost.
	697: The evaluation has been abandoned due to total loss of precision.  The result is returned as zero.
	698: The evaluation has been abandoned due to failure to satisfy the termination condition.  The result is returned as zero.
	74: An internal error has occured in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	
	successfully call of the nag_bessel_k_nu function.
	
*/

double 	nag_bessel_k_nu(
	   double x, // the argument x of the function.
	   int nu,  // the argument nu of the function.
	   NagError *fail = NULL
);

/**	s18egc
		returns a sequence of values for the modified Bessel functions K_alpha+n(x) 
		for real x > 0, selected values of apha >= 0 and n = 0, 1, 2, ..., N.

Example:
	The example program evaluates k0(x), k1(x), k2(x) and k3(x) 
	at x = 0.5 and prints the results.

void test_nag_bessel_k_alpha()
{
	double b[101];
	double x = 0.5;
	int ia = 0;
	int ja = 1;
	int nl = 3;
	printf("x = %4.1f ia = %6ld ja = %6ld nl = %6ld\n\n",x, ia, ja, nl);
	int success;
	nag_bessel_k_alpha(x, ia, ja, nl, b);
	if(success == 0)
	{
		printf("Requested values of K_alpha(X)\n\n");
		double alpha = 1.0 * ia / ja;
		printf("	alpha		K_alpha(X)\n");
		for(int i = 0; i <= nl; i++, alpha+=1.0)
			printf("%12.4e %12.4e \n",alpha, b[i]);
	}
	else 
		printf("nag_bessel_k_alpha has some problem");
}

	The output is following:
	
	x		ia		ja		nl
	
	0.5		0		1		3
	
	Requested values of K_alpha(x)
	
	alpha			K_alpha(x)
	
	0.0000e+00		9.2552e-01
	1.0000e+00		1.6554e+00
	2.0000e+00		7.5502e+00
	3.0000e+00		6.2058e+01
	

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	96: On entry, nu = <value>.  Constraint: 1 <= abs(nu) <= 3.
	90: On entry, nl = <value>.  Constraint: 0<= abs(nl) <= 100.
	91: On entry, ia = <value>, ja = <value>.  Constraint: ia and ja must constitue a valid pair (ia, ja).
	695: The evaluation has been abandoned due to the likelihood of overflow.  
	696: The evaluation has been completed but some precision has been lost.
	697: The evaluation has been abandoned due to total loss of precision.  .
	698: The evaluation has been abandoned due to failure to satisfy the termination condition.  
	74: An internal error has occured in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	
	successfully call of the nag_bessel_k_alpha function.


*/

int 	nag_bessel_k_alpha(
	double x, // the argument x of the function.
	int ia, 
	int ja, // the numerator i and denominator j, respectively, of the order alpha = i/j of the first member in the rquered sequence of function values.  Only the following combinations of pairs of values of i and j are allowed: (ia,ja) = (0,1), (1,2), (1,3), (1,4), (2,3) or (3,4).
	int nl, // the value of N.  Note that the order of the last member in the required sequence of function values is given by alpha + N.
	double b[] // with fail.code = NE_NOERROR or fail.code = NW_SOME_PRECISION_LOSS, the required sequence of function values: b(n) contains K_alpha+n(x) n = 0, 1, ..., N.

);

/**	s18ehc
		eturns a sequence of values for the modified Bessel functions (e^x)K_alpha+n(x) 
		for real x > 0, selected values of apha >= 0 and n = 0, 1, 2, ..., N.


Example:
	The example program evaluates e^xk0(x), e^xk1(x),e^xk2(x) and e^xk3(x) 
	at x = 0.5 and prints the results.

void test_nag_bessel_kalpha_scaled()
{
	double b[101];
	double x = 0.5;
	int ia = 0;
	int ja = 1;
	int nl = 3;
	printf("x = %4.1f ia = %6ld ja = %6ld nl = %6ld\n\n",x, ia, ja, nl);
	int success;
	success = nag_bessel_kalpha_scaled(x, ia, ja, nl, b);
	if(success == 0)
	{
		printf("Requested values of K_alpha(X)\n\n");
		double alpha = 1.0 * ia / ja;
		printf("	alpha		K_alpha(X)\n");
		for(int i = 0; i <= nl; i++, alpha+=1.0)
			printf("%12.4e %12.4e \n",alpha, b[i]);
	}
	else 
		printf("nag_bessel_k_alpha has some problem");
}

The output is following:
	
	x		ia		ja		nl
	
	0.5		0		1		3
	
	Requested values of K_alpha(x)
	
	alpha			exp(x)K_alpha(x)
	
	0.0000e+00		1.5241e-01
	1.0000e+00		2.7310e+00
	2.0000e+00		1.2448e+00
	3.0000e+00		1.0232e+01
	

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	96: On entry, nu = <value>.  Constraint: 1 <= abs(nu) <= 3.
	90: On entry, nl = <value>.  Constraint: 0<= abs(nl) <= 100.
	91: On entry, ia = <value>, ja = <value>.  Constraint: ia and ja must constitue a valid pair (ia, ja).
	695: The evaluation has been abandoned due to the likelihood of overflow.  
	696: The evaluation has been completed but some precision has been lost.
	697: The evaluation has been abandoned due to total loss of precision.  .
	698: The evaluation has been abandoned due to failure to satisfy the termination condition.  
	74: An internal error has occured in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	
	successfully call of the nag_bessel_kalpha_scaled function.

*/

int 	nag_bessel_kalpha_scaled(
	double x, // the argument x of the function.
	int ia, 
	int ja, // the numerator i and denominator j, respectively, of the order alpha = i/j of the first member in the rquered sequence of function values.  Only the following combinations of pairs of values of i and j are allowed: (ia,ja) = (0,1), (1,2), (1,3), (1,4), (2,3) or (3,4).
	int nl, // the value of N.  Note that the order of the last member in the required sequence of function values is given by alpha + N.
	double b[] // with fail.code = NE_NOERROR or fail.code = NW_SOME_PRECISION_LOSS, the required sequence of function values: b(n) contains K_alpha+n(x) n = 0, 1, ..., N.

);

/**	s18ejc
		returns a sequence of values for the modified Bessel functions I_alpha+n-1(x) or I_alpha-n+1(x) for real x, non-negative alpha < 1 and n = 1, 2, ..., |N| + 1.

Example:
	The example program evalutes j0(x), j1(x), j2(x) and j3(x) at x = 0.5, and prints the results.

**note this function requires the support of the complex array 
*/

int 	nag_bessel_i_alpha(
	double x, // the argument x of the functions.
	double a, // the order alpha of the first member in the required sequence of function values.
	int nl,  //  the value of N.
	complex *b // with fail.code = NE_NOERROR or fail.code = NW_SOME_PRECISION_LOSS, the required sequence of function values: b(n) contains I_alpha+n-1(x) if nl >= 1 and I_alpha-n+1(x) otherwise, for n = 0, 1, ..., abs(nl) + 1.

	
);

/**	s18ekc
		returns a sequence of values for the Bessel functions J_alpha+n-1(x) or 
		J_alpha-n+1(x) for real x, non-negative alpha < 1 and n = 1, 2, .. , abs(N) +1.
		
Example:
	The example program evalutes j0(x), j1(x), j2(x) and j3(x) at x = 0.5, and prints the results.

void test_nag_bessel_j_alpha()
{
	double c_1 = 1.0;
	matrix <complex> b;
	b.SetSize(1,101);
	double a;
	double alpha;
	double d_1;
	double x;
	int i, nl;
	a = 0.0;
	x = 0.5;
	nl = 3;
	printf("%4.1f	%4.1f	%6ld\n\n",x, a, nl);
	if(nag_bessel_j_alpha(x, a, nl, b) == 0)
	{
		printf("\n Requested values of I_alpha(x)\n\n");
		alpha = a;
		printf("	alpha				I_alpha(x)\n");
		for(i = 1; i <= fabs(nl) + 1; i++)
		{
			complex cmTemp = b[0][i-1];
			printf("%12.4e	(%12.4e, %12.4e)\n", alpha, cmTemp.m_re, cmTemp.m_im);
			d_1 = 1.0 * nl;
			if(d_1 > 0)
				alpha = fabs(c_1) + alpha;
			else
				alpha  =  -fabs(c_1) + alpha;
		}
	}
	else
		printf("there is some problem with the fuction");	
}

	The output is following:
	
	x		a		nl
	0.5		0.0		3
	
	Requested values of J_alpha(x)
	
		alpha					J_alpha(x)

		0.0000e+00		(9.3847e-01, 0.0000e+00)
		1.0000e+00		(2.4227e-01, 0.0000e+00)
		2.0000e+00		(3.0604e-02, 0.0000e+00)
		3.0000e+00		(2.5637e-03, 0.0000e+00)
		
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	118: On entry, x = <value>, nl = <value>.  Constraint x not equal 0 when nl < 0.
	96: On entry, a = <value>.  Constraint 0 <= a < 1.
	90: On entry, nl = <value>. Constraint: abs(nl) <= 101.
	695: The evaluation has been abandoned due to the likelihood of overflow.
	696: The evaluation has been completed but some precision has been lost.
	697: The evaluation has been abandoned due to total loss of precision.
	698: The evaluation has been abandoned due to failure to satisfy the termination condition.
	74: An internal error has occured in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	
	successfully call of the nag_bessel_j_alpha function.

*/
	int nag_bessel_j_alpha(
	 	double x, // the argument x of the function.
	 	double a, // the order alpha of the first member in the required sequence of function values.
	 	int nl, // the value of N
	 	complex b[], // the dimension, dim1, of the array b must be at least abs(nl) + 1.
	 NagError *fail = NULL
);	



/**	s19aac
		returns a value for the Kelvin function berx.

Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_kelvin_ber()
{
	double y;
	double x[] ={0.1, 1.0, 2.5, 5.0, 10.0, 15.0, -1.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 7; i++)
	{
		y = nag_kelvin_ber(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}		
}


	The output is following:
	
		x 				y
	1.000e-01 		1.000e+00
	1.000e+00 		9.844e-01
	2.500e+00 		4.000e-01
	5.000e+00 	   -6.230e+00
	1.000e+01 		1.388e+02
	1.500e+01 	   -2.967e+03
	-1.000e+00 		9.844e-01

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	7: On entry, |x| must not be greater than _value_: x = _value_.  |x| is too large for an accurate result to be returned and the function returns zero.

	successfully call of the nag_kelvin_ber function.
*/

double 	nag_kelvin_ber(
		double x, // the argument x of the function.
		NagError *fail = NULL
);

/**	s19abc
		returns a value for the Kelvin function bei x.
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
void test_nag_kelvin_bei()
{
	double y;
	double x[] ={0.1, 1.0, 2.5, 5.0, 10.0, 15.0, -1.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 7; i++)
	{
		y = nag_kelvin_bei(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}		
}


	The output is following:
	
	x 				y
	
	1.000e-01 	2.500e-03
	1.000e+00 	2.496e-01
	2.500e+00 	1.457e+00
	5.000e+00	1.160e-01
	1.000e+01 	5.637e+01
	1.500e+01 	-2.953e+03
	-1.000e+00 	2.496e-01
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
		
	7: On entry, |x| must not be greater than _value_: x = _value_.  |x| is too large for an accurate result to be returned and the function returns zero.

	successfully call of the nag_kelvin_bei function.

*/


double 	nag_kelvin_bei(
		double x, // the argument x of the function.
		NagError *fail = NULL
);

/**	s19acc
		returns a value for the Kelvin function ker x.

Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_kelvin_ker()
{
	double y;
	double x[] ={0.1, 1.0, 2.5, 5.0, 10.0, 15.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 6; i++)
	{
		y = nag_kelvin_ker(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}		
}


	The output is following:
	
	x			 y
	
	1.000e-01	 2.420e+00
	1.000e+00	 2.867e-01
	2.500e+00	-6.969e-02
	5.000e+00	-1.151e-02
	1.000e+01	 1.295e-04
	1.500e+01	-1.514e-08

Parameters:

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

	7: On entry, x must not be greater than _value_: x = _value_.  x is too large, the result under.ows and the function returns zero.
	12: On entry, x must not be less than or equal to 0.0: x = _value_.  The function is undefined and returns zero.

	successfully call of the nag_kelvin_ker function.

*/

double 	nag_kelvin_ker(
		double x, // the argument x of the function.
		NagError *fail = NULL
);

/**	s19adc
		returns a value for the Kelvin function kei x.
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_kelvin_kei()
{
	double y;
	double x[] ={0.1, 1.0, 2.5, 5.0, 10.0, 15.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 6; i++)
	{
		y = nag_kelvin_kei(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}		
}

	The output is following:
	
	x			 	y
	
	0.000e+00	 -7.854e-01
	1.000e-01	 -7.769e-01
	1.000e+00	 -4.950e-01
	2.500e+00	 -1.107e-01
	5.000e+00	  1.119e-02
	1.000e+01	 -3.075e-04
	1.500e+01	  7.963e-06

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	7: On entry, x must not be greater than _value_: x = _value_.  x is too large, and the result under.ows and the function returns zero.
	11: On entry, x must not be less than 0.0: x = _value_.  The function is undefined and returns zero.
	
	successfully call of the nag_kelvin_kei function.
*/

double 	nag_kelvin_kei(
		double x, // the argument x of the function.
		NagError *fail = NULL
);

/**	s20acc
		returns a value for the Fresnel Integral S(x).

Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_fresnel_s()
{
	double y;
	double x[] ={0.1, 1.0, 2.5, 5.0, 10.0, 15.0};		
	printf("   x		  y\n");
	for(int i = 0; i < 6; i++)
	{
		y = nag_fresnel_s(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}		
}

	The output is following:
	
		x			 y
	
     1.000e-001  5.236e-004
     1.000e+000  4.383e-001
     2.500e+000  6.192e-001
     5.000e+000  4.992e-001
     1.000e+001  4.682e-001
     1.500e+001  5.000e-001

	 
Parameters:

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

*/

double 	nag_fresnel_s(
		double x // the argument x of the function.
);

/**	s20adc
		returns a value for the Fresnel Integral C(x).
		
Example:
	The following program reads values of the argument x from a file, 
	evaluates the function at each value of x and prints the results.
	
void test_nag_fresnel_c()
{
	double  y;
	double x[] ={0.0, 0.5, 1.0, 2.0, 4.0, 5.0, 6.0, 8.0, 10.0, -1.0, 1000.0};
 
	printf("   x		  y\n");
	for(int i = 0; i < 11; i++)
	{
		y = nag_fresnel_c(x[i]);
		printf("%12.3e%12.3e\n", x[i], y);
	}
}

	The output is following:
	
	 	x 			y
	 	
	 0.000e+00 	  0.000e+00
	 5.000e-01	  4.923e-01
	 1.000e+00	  7.799e-01
	 2.000e+00	  4.883e-01
	 4.000e+00	  4.984e-01
	 5.000e+00	  5.636e-01
	 6.000e+00	  4.995e-01
	 8.000e+00	  4.998e-01
	 1.000e+01 	  4.999e-01
	-1.000e+00 	 -7.799e-01
	 1.000e+03 	  5.000e-01

Parameters:

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

*/

double 	nag_fresnel_c(
		double x // the argument x of the function.
);

/**	s21bac
		returns a value of an elementary integral, which occurs as a
		degenerate case of an elliptic integral of the first kind.

Example:
	This example program simply generates a small set of non-extreme arguments 
	which are used with the function to produce the table of low accuracy results.
	
void test_nag_elliptic_integral_rc()
{
	double rc, x, y;
	int ix;
	printf(" x y s21bac \n");
	for (ix=1; ix<=3; ix++)
	{
		x = ix*0.5;
		y = 1.0;
		rc = nag_elliptic_integral_rc(x, y);
		printf("%7.2f%7.2f%12.4f\n", x, y, rc);
	}
}

	The output is following:
	
	x 	   	  y 	s21bac
	
	0.50 	1.00 	1.1107
	1.00 	1.00	1.0000
	1.50 	1.00 	0.9312

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	5: On entry, x must not be less than 0.0: x = _value_.  The function is undefined.
	10: On entry, y must not be equal to 0.0: y = _value_.  The function is undefined and returns zero.

	successfully call of the nag_elliptic_integral_rc function.


*/

double 	nag_elliptic_integral_rc(
	   double x, // the argument x of the function.
	   double y, // the argument y of the function.
	   NagError *fail = NULL
);

/**	s21bbc
		returns a value of the symmetrised elliptic integral of the first kind.

Example:
	This example program simply generates a small set of non-extreme arguments 
	which are used with the function to produce the table of low accuracy results

void test_nag_elliptic_integral_rf()
{
	double rf, x, y, z;
	int ix;
	printf(" x y z s21bbc\n");
	for (ix=1; ix<=3; ix++)
	{
		x = ix*0.5;
		y = (ix+1)*0.5;
		z = (ix+2)*0.5;
		rf = nag_elliptic_integral_rf(x, y, z);
		printf("%7.2f%7.2f%7.2f%12.4f\n", x, y, z, rf);
	}
}

	The output is following:
	
	x 		y 		z 		s21bbc
	
	0.50 	1.00 	1.50 	1.0281
	1.00 	1.50 	2.00 	0.8260
	1.50 	2.00 	2.50 	0.7116
	

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	5: On entry, x must not be less than 0.0: x = _value_.  On entry, y must not be less than 0.0: y = _value_.  On entry, z must not be less than 0.0: z = _value_. The function is undefined.
	10: On entry, _parameters_ must not be equal to 0.0: _parameters_ = _value_.  On entry,t wo or more of x, y and z are zero; the function is undefined and the function returns zero.
	
	successfully call of the nag_elliptic_integral_rf function.


*/

double 	nag_elliptic_integral_rf(
	   double x, // the argument x of the function.
	   double y, // the argument y of the function.
	   double z, // the argument z of the function.
	   NagError *fail = NULL
);

/**	s21bcc
		returns a value of the symmetrised elliptic integral of the second kind.

Example:
	This example program simply generates a small set of non-extreme arguments 
	which are used with the function to produce the table of low accuracy results.

void test_nag_elliptic_integral_rd()
{
	double rd, x, y, z;
	int ix, iy;
	printf (" x y z s21bcc \n");
	for ( ix=1; ix<=3; ix++)
	{
		x = ix*0.5;
		for (iy=ix; iy<=3; iy++)
		{
			y = iy*0.5;
			z = 1.0;
			rd = nag_elliptic_integral_rd(x, y, z);
			printf (" %7.2f%7.2f%7.2f%12.4f\n", x, y, z, rd);
		}
	}
}

	The output is following:
	
	x 		y 		z 		s21bcc
	
	0.50 	0.50 	1.00 	1.4787
	0.50 	1.00 	1.00 	1.2108
	0.50 	1.50 	1.00 	1.0611
	1.00 	1.00 	1.00 	1.0000
	1.00 	1.50 	1.00 	0.8805
	1.50 	1.50 	1.00 	0.7775

Parameters:

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

	5: On entry, x must not be less than 0.0: x = _value_.  On entry, y must not be less than 0.0: y = _value_.  The function is undefined. On entry, _parameters_ must not be less than _value_: _parameters_ = _value_.  On entry, either z is too close to zero or both x and y are too close to zero: there is a danger of setting overflow.
	10: On entry, x+y must not be equal to 0.0: x+y = _value_.  Both x and y are zero and the function is undefined.
	6: On entry, z must not be less than or equal to 0.0: z = _value_.  The function is undefined.
	8: On entry, x must not be greater than or equal to _value_: x = _value_.  On entry, y must not be greater than or equal to _value_: y = _value_.  On entry, z must not be greater than or equal to _value_: z = _value_.  There is a danger of setting under.ow and the function returns zero.
	
	successfully call of the nag_elliptic_integral_rd function.


*/

double 	nag_elliptic_integral_rd(
	   double x, // the argument x of the function.
	   double y, // the argument y of the function.
	   double z, // the argument z of the function.
	   NagError *fail = NULL
);

/**	s21bdc
		returns a value of the symmetrised elliptic integral of the third kind.

Example:
	This example program simply generates a small set of non-extreme arguments 
	which are used with the function to produce the table of low accuracy results.

void test_nag_elliptic_integral_rj()
{
	double r, rj, x, y, z;
	int ix, iy, iz;
	printf(" x y z r s21bdc \n");
	for (ix = 1; ix <= 3; ix++)
	{
		x = ix*0.5;
		for ( iy=ix; iy<=3; iy++)
		{
			y = iy*0.5;
			for ( iz=iy; iz<=3; iz++)
			{
				z = iz*0.5;
				r = 2.0;
				rj = nag_elliptic_integral_rj(x, y, z, r);
				printf (" %7.2f%7.2f%7.2f%7.2f%12.4f\n", x, y, z, r, rj);
			}
		}
	}
}


	The output is following:
	
	x 		y 		z 		r 		s21bdc
	
	0.50 	0.50 	0.50 	2.00 	1.1184
	0.50 	0.50 	1.00 	2.00 	0.9221
	0.50 	0.50 	1.50 	2.00 	0.8115
	0.50 	1.00 	1.00 	2.00 	0.7671
	0.50 	1.00 	1.50 	2.00 	0.6784
	0.50 	1.50 	1.50 	2.00 	0.6017
	1.00 	1.00 	1.00 	2.00 	0.6438
	1.00 	1.00 	1.50 	2.00 	0.5722
	1.00 	1.50 	1.50 	2.00 	0.5101
	1.50 	1.50 	1.50 	2.00 	0.4561

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	5: On entry, x must not be less than 0.0: x = _value_.  On entry, y must not be less than 0.0: y = _value_.  On entry, z must not be less than 0.0: z = _value_.  The function is undefined.  On entry, _parameters_ must not be less than _value_: _parameters_ = _value_.  On entry, either r is too close to zero, or any two of x, y and z are too close to zero; there is a danger of setting overflow.
	10: On entry, _parameters_ must not be equal to 0.0: _parameters_ = _value_.  At least two of x, y and z are zero and the function is undefined.  On entry, r must not be equal to 0.0: r = _value_.  The function is undefined.
	7: On entry, x must not be greater than _value_: x = _value_.  On entry, y must not be greater than _value_: y = _value_.  On entry, z must not be greater than _value_: z = _value_.  On entry, |r| must not be greater than _value_: |r| = _value_.  
	
	successfully call of the nag_elliptic_integral_rj function.


*/

double 	nag_elliptic_integral_rj(
	   double x, // the argument x of the function.
	   double y, // the argument y of the function.
	   double z, // the argument z of the function.
	   double r, // the argument rho of the function.
	   NagError *fail = NULL
);

/**	s21cbc
		evaluates the Jacobian elliptic functions snz, cnz, and dnz for a complex argument z.

Example:
	The program evaluates snz, cnz, and dnz at z = -2.0+3.0i 
	when k = 0.5 and prints the results

void test_nag_jacobian_elliptic()
{
	complex cn, dn, sn, z;
	double ak2;

	z = -2.0 + 3.0i;
	ak2 = 0.25;
	nag_jacobian_elliptic(z, ak2, &sn, &cn, &dn);
	printf("		z			ak2\n");
	printf(" (%8.4f,%8.4f)	%10.2f\n\n",z.m_re,z.m_im, ak2);
	printf("	sn				cn				dn\n ");
	printf("(%8.4f, %8.4f)	",sn.m_re, sn.m_im);
	printf("(%8.4f, %8.4f)	",cn.m_re, cn.m_im);
	printf("(%8.4f, %8.4f)	",dn.m_re, dn.m_im);
	printf("\n");	
}

	The output is following:
	
		z							ak2
		
	(-2.0000, 3.0000)				0.25
	
		sn						cn				dn
	
	(-1.5865, 0.2456)	(0.3125, 1.2468)	(-0.6395, -0.1523)
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	96: On entry, ak2 = <value>.  Constraint: 0 <= ak2 <= 1
	102: On entry, z = (<value>, <value>).  Constraint: abs(z.ro) <= lamda and abs(z.im) <= lamda, where lamda = (1/X02AMC)^(1/2).
	74: An internal error has occurred in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	
	successfully call of the nag_jacobian_elliptic function.

	
*/

int		nag_jacobian_elliptic(
	complex z, // the argument z of the function.
	double ak2, // the value of k^2.
	complex *sn, // the value of the functions snz.
	complex *cn, // the value of the functions cnz.
	complex *dn // the value of the functions dnz.
);

/**	s21ccc
		returns the value of one of the Jacobian theta functions theta_0(x,q), 
		theta_0(x,q), theta_1(x,q), theta_2(x,q), theta_3(x,q), or 
		theta_0(x,q) for a real argument x and non-negative q <= 1.
	
Example:
	The example program evaluates theta(x,q) at x =0.7 
	when q = 0.4, and prints the results.

void test_nag_jacobian_theta()
{
	double q, x, y;
	int k;

	k = 2;
	x = 0.7;
	q = 0.4;
	y  = nag_jacobian_theta(k, x, q);
	printf("%d	%f	%f	%10.4e", k, x, q, y);	
}


	The output is following:

	k		x		q			y
	
	2		0.7		0.4		-6.9289e-01

Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	90: On entry, k = <value>.  Constraint: 0 <= k <= 4.
	96: On entry, q = <value>.  Constraint: 0 <= q <= 1.  On entry, x = <value>.  Constraint: (x - 0.5) must not be an integer when q = 1 and k = 1. On entry, x = <value>.  Constraint: x must not be an integer when q = 1 and k = 2.
	699: The evaluation has been abandoned because the function value is infinite.
	74: An internal error has occurred in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
	 
	successfully call of the nag_jacobian_theta function.	


*/

double		nag_jacobian_theta(
		int k, // the function theta_k(x,q) to be evaluated.  Note that k = 4 is equivalent to k = 0.
		double x, // the argument x of the function.
		double q, // the argument q of the function.
		NagError *fail = NULL
);

/**	s21dac
		returns the value of the general elliptic integral of the second kind F(z, k',a,b) for a complex argument z.

Example:
	The example program evaluates the elliptic integral of the first kind 
	F1(z,k') given by where z = 1.2 + 3.7i; and k' = 0.5, and prints the results.
	
void test_nag_elliptic_integral_f()
{
 	complex y;
 	matrix<complex> z;
 	z.SetSize(1,2);
 	double a[2], akp[2], b[2];
 	
 	z[0][0] = 1.2 + 3.7i;
 	z[0][1] = 9.2 - 3.4i;
 	akp[0] = 0.5;
 	akp[1] = 0.8;
 	a[0] = 1.0;
 	a[1] = 0.4;
 	
 	b[0] = 1.0;
 	b[1] = 2.7;
 	for(int i = 0; i < 2; i++)
 	{
 		complex temp = z[0][i];
 	
  		y =  nag_elliptic_integral_f(temp, akp[i], a[i], b[i]);
  	
		printf("(%4.1f,	%4.1f)	%7.1f	%7.1f	%7.1f	(%12.4e,	%12.4e)\n",temp.m_re, temp.m_im, akp[i], a[i], b[i], y.m_re, y.m_im);
 	}
}
 
	
	The output is following:
	
	z				akp		a		b			y
	
	(1.2, 3.7)		0.5		1.0		1.0		(1.9713e+00, 5.0538e-01)
	(9.2, -3.4)		0.8		0.4		2.7		(2.5042e+00,-1.1709e-01)
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	102: On entry, z = (<value>,<value>). Constraints: 0 <= zire <= lamda, |z.im| <= lamda, where lamda^6 = 1/X02AMC.
	700: The iterative procedure used to evaluate the integral has failed to converge.  The result is returned as zero.
	74: An internal error has occurred in this function.  Check the function call and any array sizes.  If the call is correct then please consult NAG for assistance.
			 
	successfully call of the nag_elliptic_integral_f function.	


*/

complex	nag_elliptic_integral_f(
		
		complex z, // the argument z of the function.
		double akp, // the argument k' of the function.
		double a, // the argument a of the function.
		double b, // the argument b of the function.
		NagError *fail = NULL
);

/**	s22aac	
		returns a sequence of values for either the unnormalized or 
		normalized Legendre functions of the first kind P_m_n(x) or _P_m_n(x)
		for real x of a given order m and degree n = 0, 1, ..., N.
		

Example:
	The following program reads the values of the arguments x, m, N, calculates the sequence 
	of unnormalized assocaited Legendre function values P and prints out the results.

void test_nag_legendre_p()
{
	int m, n, nl, mode;
	double x;
	double p[101];
	mode = 1;
	x = 0.5;
	m = 2;
	nl = 3;
	nag_legendre_p(mode, x, m, nl, p);
	printf("mode	x	m	nl\n\n");
	printf("%3ld	%5.2f%6ld%6ld\n\n",mode, x, m, nl);
	printf("Unnormalized associated Legendre function values\n");
	printf("\n");
	printf("\n	P(n)\n");
	for(n = 0; n <= nl; n++)
		printf("%2ld %12.4e\n", n, p[n]);	
}


	The output is following:
	
	mode	x		m		nl
	
	1		0.5		2		3

	Unnormalized associated Legendre function values
	
	n			P(n)
	
	0		0.0000e+00
	1		0.0000e+00
	2		2.2500e+00
	3		5.6250e+00
	
Parameters:

Return:
	This function returns NAG error code, 0 if no error.
	
	96: On entry, x = <value>.  Constraint: abs(x) <= 1.
	90: On entry, mode = <value>.  Constraint: mode = 1 or 2.  On entry, nl = <value>.  Constraint: nl >= 0.  On entry, m = <value>.  Constraint: abs(m) <= 27.
	91: On entry, nl = <value>, m = <value>.  Constraint: nl <= 100 when m = 0. On entry, nl = <value>, m = <value>.  Constraint: nl <= 55 - abs(m) when m not equal 0.
		 
	successfully call of the nag_legendre_p function.	

*/

int		nag_legendre_p(
	int mode, // indicates whether the sequence of function values is to be returned unnormalized or normalized as follows: if mode = 1, then the sequence of function values is returned unnormalized; if mode = 2, then the sequence of function values is returned normalized.
	double x, // the argument x of the function.
	int m,  // the order m of the function.
	int nl, // the degree N of the last function required in the sequence.
	double p[] // the required sequence of function values as follows: if mode = 1, p(n) contains P_m_n(x), for n = 0, 1, ..., N.  if mode = 2, p(n) contains _P_m_n(x),for n = 0, 1, ..., N.
);

/* end proto */
#endif //!_O_NAG_S_H