/*------------------------------------------------------------------------------*
* File Name: NAGSpecialFunctions.c												*
* Creation: ER, 04/25/03														*
* Purpose: OriginC file to expose NAG special functions to LabTalk				*
* Copyright (c) OriginLab Corp.	2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007	*
* All Rights Reserved															*
* 																				*
* Modification Log:															    *
*------------------------------------------------------------------------------*/

////////////////////////////////////////////////////////////////////////////////////
#include <origin.h>
#include <NAG/OCN_s.h> 
////////////////////////////////////////////////////////////////////////////////////


////////////////////////////////////////////////////////////////////////////////////
// This file contains functions that expose the NAG Special functions to LabTalk
// so that they can be accessed from the script window, and in the GUI in places
// such as the Set Column Values dialog, the Function Plotting dialog etc.
//
// Add this file to the system folder of the Code Builder workspace in order to
// have these functions available in all Origin sessions.
//
// You can add your own custom math functions to this file and have them available
// in all sessions as well.
//
// Note: Not all functions from the NAG Special functions chapter are included here.
// Those that use complex variables or arrays for input/ouput are not included.
// For a complete listing of all NAG Special functions and for syntax/documentation
// on these functions, refer to the Special Functions section under NAG, under the
// Origin C Programming Help file.
//
//
// exponential integral
double exp_integral(double x)
{
	NagError fail;
	double y = nag_exp_integral(x, &fail);
	if(fail.code != 0) report_nag_error(fail);
	return y;
}

// cosine integral
double cos_integral(double x)
{
	NagError fail;
	double y = nag_cos_integral(x, &fail);
	if(fail.code != 0) report_nag_error(fail);
	return y;
}


// sine integral
double sin_integral(double x)
{
	double y = nag_sin_integral(x);
	return y;
}

// Gamma function
double gamma(double x)
{
	NagError fail;
	double y = nag_gamma(x, &fail);
	if(fail.code != 0) report_nag_error(fail);
	return y;
}

// Gamma function
double log_gamma(double x)
{
	NagError fail;
	double y = nag_log_gamma(x, &fail);
	if(fail.code != 0) report_nag_error(fail);
	return y;
}

// Polygamma function
double real_polygamma(double x, int k)
{
	NagError fail;
	double y = nag_real_polygamma(x, k, &fail);
	if(fail.code != 0) report_nag_error(fail);
	return y;
}


// Incomplete gamma function
double incomplete_gamma(double a, double x)
{
	double p, q;
	double tol = 1e-5;
	int ierr = nag_incomplete_gamma(a, x, tol, &p, &q);
	return p;
}

// Cumulative normal distribution function
double cumul_normal(double x)
{
	double y = nag_cumul_normal(x);
	return y;
}

// Complement of the cumulative normal distribution function
double cumul_normal_complem(double x)
{
	double y = nag_cumul_normal_complem(x);
	return y;
}

// Complement of the error function
double erfc(double x)
{
	double y = nag_erfc(x);
	return y;
}

// Error function
// Use Erf instead as erf for name - to avoid conflict with built-in error function
double Erf(double x)
{
	double y = nag_erf(x);
	return y;
}


// Bessel Y0
double bessel_y0(double x)
{
	NagError fail;
	double y = nag_bessel_y0(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}


// Bessel Y1
double bessel_y1(double x)
{
	NagError fail;
	double y = nag_bessel_y1(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Bessel J0
double bessel_j0(double x)
{
	NagError fail;
	double y = nag_bessel_j0(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Bessel J1
double bessel_j1(double x)
{
	NagError fail;
	double y = nag_bessel_j1(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Airy function Ai
double airy_ai(double x)
{
	NagError fail;
	double y = nag_airy_ai(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Airy function Bi
double airy_bi(double x)
{
	NagError fail;
	double y = nag_airy_bi(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Derivative of Airy function Ai
double airy_ai_deriv(double x)
{
	NagError fail;
	double y = nag_airy_ai_deriv(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Derivative of Airy function Bi
double airy_bi_deriv(double x)
{
	NagError fail;
	double y = nag_airy_bi_deriv(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Modified Bessel Function K0
double bessel_k0(double x)
{
	NagError fail;
	double y = nag_bessel_k0(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Modified Bessel Function K1
double bessel_k1(double x)
{
	NagError fail;
	double y = nag_bessel_k1(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Modified Bessel Function I0
double bessel_i0(double x)
{
	NagError fail;
	double y = nag_bessel_i0(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Modified Bessel Function I1
double bessel_i1(double x)
{
	NagError fail;
	double y = nag_bessel_i1(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Scaled Modified Bessel Function K0
double bessel_k0_scaled(double x)
{
	NagError fail;
	double y = nag_bessel_k0_scaled(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Scaled Modified Bessel Function K1
double bessel_k1_scaled(double x)
{
	NagError fail;
	double y = nag_bessel_k1_scaled(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Scaled Modified Bessel Function I0
double bessel_i0_scaled(double x)
{
	double y = nag_bessel_i0_scaled(x);
	return y;
}

// Scaled Modified Bessel Function I1
double bessel_i1_scaled(double x)
{
	double y = nag_bessel_i1_scaled(x);
	return y;
}


// Scaled Modified Bessel Function I_nu
double bessel_i_nu_scaled(double x, int nu)
{
	NagError fail;
	double y = nag_bessel_i_nu_scaled(x, nu, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Scaled Modified Bessel Function K_nu
double bessel_k_nu_scaled(double x, int nu)
{
	NagError fail;
	double y = nag_bessel_k_nu_scaled(x, nu, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Modified Bessel Function I_nu
double bessel_i_nu(double x, int nu)
{
	NagError fail;
	double y = nag_bessel_i_nu(x, nu, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Modified Bessel Function K_nu
double bessel_k_nu(double x, int nu)
{
	NagError fail;
	double y = nag_bessel_k_nu(x, nu, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Kelvin Function berx
double kelvin_ber(double x)
{
	NagError fail;
	double y = nag_kelvin_ber(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Kelvin Function beix
double kelvin_bei(double x)
{
	NagError fail;
	double y = nag_kelvin_bei(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Kelvin Function kerx
double kelvin_ker(double x)
{
	NagError fail;
	double y = nag_kelvin_ker(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Kelvin Function keix
double kelvin_kei(double x)
{
	NagError fail;
	double y = nag_kelvin_kei(x, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// Fresnel Integral S
double fresnel_s(double x)
{
	double y = nag_fresnel_s(x);
	return y;
}

// Fresnel Integral C
double fresnel_c(double x)
{
	double y = nag_fresnel_c(x);
	return y;
}

// Elliptical Integral rc
double elliptic_integral_rc(double x, double y)
{
	NagError fail;
	double yret = nag_elliptic_integral_rc(x, y, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return yret;
}

// Elliptical Integral rf
double elliptic_integral_rf(double x, double y, double z)
{
	NagError fail;
	double yret = nag_elliptic_integral_rf(x, y, z, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return yret;
}

// Elliptical Integral rd
double elliptic_integral_rd(double x, double y, double z)
{
	NagError fail;
	double yret = nag_elliptic_integral_rd(x, y, z, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return yret;
}

// Elliptical Integral rj
double elliptic_integral_rj(double x, double y, double z, double r)
{
	NagError fail;
	double yret = nag_elliptic_integral_rj(x, y, z, r, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return yret;
}

// Jacobian Theta Function
double jacobian_theta(int k, double x, double q)
{
	NagError fail;
	double y = nag_jacobian_theta(k, x, q, &fail);
	if(fail.code !=0) report_nag_error(fail);
	return y;
}

// This function prints error messages from these math functions, to the script window.
// Error messages are suppressed if the LabTalk variable report_nag_sf_errors is set to zero.
static void report_nag_error(NagError& fail)
{
	double dErrFlag;
	LT_get_var("REPORT_NAG_SF_ERRORS", &dErrFlag);
	if(dErrFlag)
	{
		using type = LabTalk.type;
		double dOldRedir = type.redirection;
		type.redirection = 5;
		printf("NAG Special Function Error - Error Code: %d\n", fail.code);
		string strErrMsg = fail.message;
		printf("Error message:\n%s\n", strErrMsg);
		type.redirection = dOldRedir;
	}
}


// This function types a list of all the NAG functions in this file, to the script window.
void List_NAG_SF()
{
	printf("----------------------------------------------------------------------------------------\n");
	printf("List of NAG Special Functions:\n");
	printf("----------------------------------------------------------------------------------------\n");
	printf("  exp_integral(x): \tExponential Integral.\n");
	printf("  cos_integral(x): \tCosine Integral.\n");
	printf("  sin_integral(x): \tSine Integral.\n");
	printf("\n");
	printf("  gamma(x): \t\tGamma Function.\n");
	printf("  log_gamma(x): \tLogarithm of Gamma Function.\n");
	printf("  real_polygamma(x,k): \tkth derivative of the psi function psi(x); k = 0 to 6.\n");
	printf("  incomplete_gamma(a, x): \tLower Incomplete Gamma Function (upper = 1 - lower).\n");
	printf("\n");
	printf("  cumul_normal(x): \tCumulative Normal Distribution Function.\n");
	printf("  cumul_normal_complem(x): \tComplement of the Cumulative Normal Distribution Function.\n");
	printf("  erfc(x): \t\tComplementary Error Function.\n");
	printf("  erf(x): \t\tError Function.\n");
	printf("  airy_ai(x): \t\tAiry Function Ai.\n");
	printf("  airy_bi(x): \t\tAiry Function Bi.\n");
	printf("  airy_ai_deriv(x): \tDerivative of Airy Function Ai.\n");
	printf("  airy_bi_deriv(x): \tDerivative of Airy Function Bi.\n");
	printf("\n");
	printf("  bessel_y0(x): \tBessel Function Y0.\n");
	printf("  bessel_y1(x): \tBessel Function Y1.\n");
	printf("  bessel_j0(x): \tBessel Function J0.\n");
	printf("  bessel_j1(x): \tBessel Function J1.\n");
	printf("  bessel_k0(x): \tModified Bessel Function K0.\n");
	printf("  bessel_k1(x): \tModified Bessel Function K1.\n");
	printf("  bessel_i0(x): \tModified Bessel Function I0.\n");
	printf("  bessel_i1(x): \tModified Bessel Function I1.\n");
	printf("  bessel_i_nu(x, nu): \tModified Bessel Function I_nu.\n");
	printf("  bessel_k_nu(x, nu): \tModified Bessel Function K_nu.\n");
	printf("  bessel_k0_scaled(x): \tScaled Modified Bessel Function K0.\n");
	printf("  bessel_k1_scaled(x): \tScaled Modified Bessel Function K1.\n");
	printf("  bessel_i0_scaled(x): \tScaled Modified Bessel Function I0.\n");
	printf("  bessel_i1_scaled(x): \tScaled Modified Bessel Function I1.\n");
	printf("  bessel_i_nu_scaled(x, nu): \tScaled Modified Bessel Function I_nu.\n");
	printf("  bessel_k_nu_scaled(x, nu): \tScaled Modified Bessel Function K_nu.\n");
	printf("  kelvin_ber(x): \tKelvin Function ber.\n");
	printf("  kelvin_bei(x): \tKelvin Function bei.\n");
	printf("  kelvin_ker(x): \tKelvin Function ker.\n");
	printf("  kelvin_kei(x): \tKelvin Function kei.\n");
	printf("  fresnel_s(x): \tFresnel Integral S.\n");
	printf("  fresnel_c(x): \tFresnel Integral C.\n");
	printf("  jacobian_theta(k, x, q): \tJacobian Theta Function.\n");
	printf("\n");
	printf("  elliptic_integral_rc(x, y): \tEllipitcal Integral rc.\n");
	printf("  elliptic_integral_rf(x, y, z): \tEllipitcal Integral rf.\n");
	printf("  elliptic_integral_rd(x, y, z): \tEllipitcal Integral rd.\n");
	printf("  elliptic_integral_rj(x, y, z, r): \tEllipitcal Integral rj.\n");
	printf("\n  For more information, please refer to NAG - Approximation of Special Functions\n");
	printf("  chapter of the Origin C Language Reference Help File.\n");
	printf("----------------------------------------------------------------------------------------\n");
	printf("\n");
}

// end of file.