特殊函数的C#实现

源代码来源于以下网站，并添加了第一类修正贝塞尔函数。

Special Function(s) for C# - CodeProjecthttps://www.codeproject.com/articles/11647/special-function-s-for-c

using System;
/*
	**************************************************************************
	**
	**    Class  SpecialFunction (C#)
	**
	**************************************************************************
	**    Copyright (C) 1984 Stephen L. Moshier (original C version - Cephes Math Library)
	**    Copyright (C) 1996 Leigh Brookshaw	(Java version)
	**    Copyright (C) 2005 Miroslav Stampar	(C# version [->this<-])
	**
	**    This program is free software; you can redistribute it and/or modify
	**    it under the terms of the GNU General Public License as published by
	**    the Free Software Foundation; either version 2 of the License, or
	**    (at your option) any later version.
	**
	**    This program is distributed in the hope that it will be useful,
	**    but WITHOUT ANY WARRANTY; without even the implied warranty of
	**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	**    GNU General Public License for more details.
	**
	**    You should have received a copy of the GNU General Public License
	**    along with this program; if not, write to the Free Software
	**    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
	**************************************************************************
	**
	**    This class is an extension of System.Math. It includes a number
	**    of special functions not found in the Math class.
	**
	*************************************************************************/


	/**
	 * This class contains physical constants and special functions not found
	 * in the System.Math class.
	 * Like the System.Math class this class is final and cannot be
	 * subclassed.
	 * All physical constants are in cgs units.
	 * NOTE: These special functions do not necessarily use the fastest
	 * or most accurate algorithms.
	 *
	 * @version $Revision: 1.8 $, $Date: 2005/09/12 09:52:34 $
	 */


public class SpecialFunction
{
	// Machine constants

	private const double MACHEP = 1.11022302462515654042E-16;
	private const double MAXLOG = 7.09782712893383996732E2;
	private const double MINLOG = -7.451332191019412076235E2;
	private const double MAXGAM = 171.624376956302725;
	private const double SQTPI = 2.50662827463100050242E0;
	private const double SQRTH = 7.07106781186547524401E-1;
	private const double LOGPI = 1.14472988584940017414;


	// Physical Constants in cgs Units

	/// <summary>
	/// Boltzman Constant. Units erg/deg(K) 
	/// </summary>
	public const double BOLTZMAN = 1.3807e-16;

	/// <summary>
	/// Elementary Charge. Units statcoulomb 
	/// </summary>
	public const double ECHARGE = 4.8032e-10;

	/// <summary>
	/// Electron Mass. Units g 
	/// </summary>
	public const double EMASS = 9.1095e-28;

	/// <summary>
	/// Proton Mass. Units g 
	/// </summary>
	public const double PMASS = 1.6726e-24;

	/// <summary>
	/// Gravitational Constant. Units dyne-cm^2/g^2
	/// </summary>
	public const double GRAV = 6.6720e-08;

	/// <summary>
	/// Planck constant. Units erg-sec 
	/// </summary>
	public const double PLANCK = 6.6262e-27;

	/// <summary>
	/// Speed of Light in a Vacuum. Units cm/sec 
	/// </summary>
	public const double LIGHTSPEED = 2.9979e10;

	/// <summary>
	/// Stefan-Boltzman Constant. Units erg/cm^2-sec-deg^4 
	/// </summary>
	public const double STEFANBOLTZ = 5.6703e-5;

	/// <summary>
	/// Avogadro Number. Units  1/mol 
	/// </summary>
	public const double AVOGADRO = 6.0220e23;

	/// <summary>
	/// Gas Constant. Units erg/deg-mol 
	/// </summary>
	public const double GASCONSTANT = 8.3144e07;

	/// <summary>
	/// Gravitational Acceleration at the Earths surface. Units cm/sec^2 
	/// </summary>
	public const double GRAVACC = 980.67;

	/// <summary>
	/// Solar Mass. Units g 
	/// </summary>
	public const double SOLARMASS = 1.99e33;

	/// <summary>
	/// Solar Radius. Units cm
	/// </summary>
	public const double SOLARRADIUS = 6.96e10;

	/// <summary>
	/// Solar Luminosity. Units erg/sec
	/// </summary>
	public const double SOLARLUM = 3.90e33;

	/// <summary>
	/// Solar Flux. Units erg/cm^2-sec
	/// </summary>
	public const double SOLARFLUX = 6.41e10;

	/// <summary>
	/// Astronomical Unit (radius of the Earth's orbit). Units cm
	/// </summary>
	public const double AU = 1.50e13;

	/// <summary>
	/// Don't let anyone instantiate this class. 
	/// </summary>
	private SpecialFunction()
	{
	}


	// Function Methods

	/// <summary>
	/// Returns the base 10 logarithm of the specified number.
	/// </summary>
	/// <param name="x"></param>
	/// <returns></returns>
	public static double log10(double x)
	{
		if (x <= 0.0) throw new ArithmeticException("range exception");
		return Math.Log(x)/2.30258509299404568401;
	}

		
	/// <summary>
	/// Returns the hyperbolic cosine of the specified number.
	/// </summary>
	/// <param name="x"></param>
	/// <returns></returns>
	public static double cosh(double x)
	{
		double a;
		a = x;
		if (a < 0.0) a = Math.Abs(x);
		a = Math.Exp(a);
		return 0.5*(a + 1/a);
	}

		
	/// <summary>
	/// Returns the hyperbolic sine of the specified number.
	/// </summary>
	/// <param name="x"></param>
	/// <returns></returns>
	public static double sinh(double x)
	{
		double a;
		if (x == 0.0) return x;
		a = x;
		if (a < 0.0) a = Math.Abs(x);
		a = Math.Exp(a);
		if (x < 0.0) return -0.5*(a - 1/a);
		else return 0.5*(a - 1/a);
	}

		
	/// <summary>
	/// Returns the hyperbolic tangent of the specified number.
	/// </summary>
	/// <param name="x"></param>
	/// <returns></returns>
	public static double tanh(double x)
	{
		double a;
		if (x == 0.0) return x;
		a = x;
		if (a < 0.0) a = Math.Abs(x);
		a = Math.Exp(2.0*a);
		if (x < 0.0) return -(1.0 - 2.0/(a + 1.0));
		else return (1.0 - 2.0/(a + 1.0));
	}

		
	/// <summary>
	/// Returns the hyperbolic arc cosine of the specified number.
	/// </summary>
	/// <param name="x"></param>
	/// <returns></returns>
	public static double acosh(double x)
	{
		if (x < 1.0) throw new ArithmeticException("range exception");
		return Math.Log(x + Math.Sqrt(x*x - 1));
	}

		
	/// <summary>
	/// Returns the hyperbolic arc sine of the specified number.
	/// </summary>
	/// <param name="xx"></param>
	/// <returns></returns>
	public static double asinh(double xx)
	{
		double x;
		int sign;
		if (xx == 0.0) return xx;
		if (xx < 0.0)
		{
			sign = -1;
			x = -xx;
		}
		else
		{
			sign = 1;
			x = xx;
		}
		return sign*Math.Log(x + Math.Sqrt(x*x + 1));
	}

		
	/// <summary>
	/// Returns the hyperbolic arc tangent of the specified number.
	/// </summary>
	/// <param name="x"></param>
	/// <re