/*********************************************************************** * * * Program to calculate the first kind Bessel function of integer * * order N, for any REAL X, using the function BESSJ(N,X). * * * * -------------------------------------------------------------------- * * * * SAMPLE RUN: * * * * (Calculate Bessel function for N=2, X=0.75). * * * * Bessel function of order 2 for X = 0.7500: * * * * Y = 0.06707400 * * * * -------------------------------------------------------------------- * * Reference: From Numath Library By Tuan Dang Trong in Fortran 77. * * * * C++ Release 1.0 By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***********************************************************************/ #include #include double BESSJ0 (double X) { /*********************************************************************** This subroutine calculates the First Kind Bessel Function of order 0, for any real number X. The polynomial approximation by series of Chebyshev polynomials is used for 0 N. For X < N, the Miller's algorithm is used to avoid overflows. ----------------------------- REFERENCE: C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS, MATHEMATICAL TABLES, VOL.5, 1962. *************************************************************************/ const IACC = 40; const double BIGNO = 1e10, BIGNI = 1e-10; double TOX,BJM,BJ,BJP,SUM,TMP; int J, JSUM, M; if (N == 0) return BESSJ0(X); if (N == 1) return BESSJ1(X); if (X == 0.0) return 0.0; TOX = 2.0/X; if (X > 1.0*N) { BJM = BESSJ0(X); BJ = BESSJ1(X); for (J=1; J0; J--) { BJM = J*TOX*BJ-BJP; BJP = BJ; BJ = BJM; if (fabs(BJ) > BIGNO) { BJ = BJ*BIGNI; BJP = BJP*BIGNI; TMP = TMP*BIGNI; SUM = SUM*BIGNI; } if (JSUM != 0) SUM += BJ; JSUM = 1-JSUM; if (J == N) TMP = BJP; } SUM = 2.0*SUM-BJ; return (TMP/SUM); } } void main() { double X,Y; int N; N=2; X=0.75; Y = BESSJ(N,X); printf("\n Bessel Function of order %d, for X=%8.4f:\n\n", N, X); printf(" Y = %10.8f\n\n", Y); } // End of file Tbessj.cpp