240 Chapter 6.Special Functions int ji float bi,bim,bip,tox,ans; if (n 2)nrerror("Index n less than 2 in bessi"); 1f(x==0.0) return 0.0; else tox=2.0/fabs(x); bip=ans=0.0; b1=1.0; for (j=2*(n+(int)sqrt(ACC*n));j>0;j--){Downward recurrence from even bim=bip+j*tox*bi; m. bip=bi; bi-bim; if (fabs(bi)>BIGNO){ Renormalize to prevent overflows. ans BIGNI: granted for 19881992 bi BIGNI; bip *BIGNI; 1.200 if (j =n)ans=bip; ans *bessio(x)/bi; Normalize with bessi0. to any from NUMERICAL RECIPESI return x <0.0&k(n&1)-ansans; (North America server computer, University Press. 令 one paper THE ART CITED REFERENCES AND FURTHER READING: Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- Programs matics Series,Volume 55(Washington:National Bureau of Standards;reprinted 1968 by Dover Publications,New York),89.8.[1] OF SCIENTIFIC Carrier,G.F.,Krook,M.and Pearson,C.E.1966,Functions of a Complex Variable (New York: McGraw-Hill),pp.220ff. 6.7 Bessel Functions of Fractional Order,Airy 1920 COMPUTING (ISBN Functions,Spherical Bessel Functions 10521 Many algorithms have been proposed for computing Bessel functions of fractional order numerically.Most of them are,in fact,not very good in practice.The routines given here are Numerical Recipes 43106 rather complicated,but they can be recommended wholeheartedly. Ordinary Bessel Functions (outside The basic idea is Steed's method,which was originally developed [1]for Coulomb wave North Software. functions.The method calculates,Y,and Y simultaneously,and so involves four relations among these functions.Three of the relations come from two continued fractions, one of which is complex.The fourth is provided by the Wronskian relation W三Y-上=名 (6.7.1) I The first continued fraction,CFl,is defined by (6.7.2) 王-20u+1)/x=2u+2/x=240 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). int j; float bi,bim,bip,tox,ans; if (n < 2) nrerror("Index n less than 2 in bessi"); if (x == 0.0) return 0.0; else { tox=2.0/fabs(x); bip=ans=0.0; bi=1.0; for (j=2*(n+(int) sqrt(ACC*n));j>0;j--) { Downward recurrence from even bim=bip+j*tox*bi; m. bip=bi; bi=bim; if (fabs(bi) > BIGNO) { Renormalize to prevent overflows. ans *= BIGNI; bi *= BIGNI; bip *= BIGNI; } if (j == n) ans=bip; } ans *= bessi0(x)/bi; Normalize with bessi0. return x < 0.0 && (n & 1) ? -ans : ans; } } CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathe￾matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), §9.8. [1] Carrier, G.F., Krook, M. and Pearson, C.E. 1966, Functions of a Complex Variable (New York: McGraw-Hill), pp. 220ff. 6.7 Bessel Functions of Fractional Order, Airy Functions, Spherical Bessel Functions Many algorithms have been proposed for computing Bessel functions of fractional order numerically. Most of them are, in fact, not very good in practice. The routines given here are rather complicated, but they can be recommended wholeheartedly. Ordinary Bessel Functions The basic idea is Steed’s method, which was originally developed [1] for Coulomb wave functions. The method calculates Jν , J
ν , Yν, and Y
ν simultaneously, and so involves four relations among these functions. Three of the relations come from two continued fractions, one of which is complex. The fourth is provided by the Wronskian relation W ≡ Jν Y
ν − YνJ
ν = 2 πx (6.7.1) The first continued fraction, CF1, is defined by fν ≡ J
ν Jν = ν x − Jν+1 Jν = ν x − 1 2(ν + 1)/x − 1 2(ν + 2)/x − ··· (6.7.2)