124599Szliu #ifndef lint 224706Selefunt static char sccsid[] = 3*31853Szliu "@(#)jn.c 4.2 (Berkeley) 8/21/85; 1.4 (ucb.elefunt) 07/13/87"; 4*31853Szliu #endif /* not lint */ 524599Szliu 624599Szliu /* 724599Szliu floating point Bessel's function of 824599Szliu the first and second kinds and of 924599Szliu integer order. 1024599Szliu 1124599Szliu int n; 1224599Szliu double x; 1324599Szliu jn(n,x); 1424599Szliu 1524599Szliu returns the value of Jn(x) for all 1624599Szliu integer values of n and all real values 1724599Szliu of x. 1824599Szliu 1924599Szliu There are no error returns. 2024599Szliu Calls j0, j1. 2124599Szliu 2224599Szliu For n=0, j0(x) is called, 2324599Szliu for n=1, j1(x) is called, 2424599Szliu for n<x, forward recursion us used starting 2524599Szliu from values of j0(x) and j1(x). 2624599Szliu for n>x, a continued fraction approximation to 2724599Szliu j(n,x)/j(n-1,x) is evaluated and then backward 2824599Szliu recursion is used starting from a supposed value 2924599Szliu for j(n,x). The resulting value of j(0,x) is 3024599Szliu compared with the actual value to correct the 3124599Szliu supposed value of j(n,x). 3224599Szliu 3324599Szliu yn(n,x) is similar in all respects, except 3424599Szliu that forward recursion is used for all 3524599Szliu values of n>1. 3624599Szliu */ 3724599Szliu 3824599Szliu #include <math.h> 39*31853Szliu #if defined(vax)||defined(tahoe) 4024599Szliu #include <errno.h> 41*31853Szliu #else /* defined(vax)||defined(tahoe) */ 4224599Szliu static double zero = 0.e0; 43*31853Szliu #endif /* defined(vax)||defined(tahoe) */ 4424599Szliu 4524599Szliu double 4624599Szliu jn(n,x) int n; double x;{ 4724599Szliu int i; 4824599Szliu double a, b, temp; 4924599Szliu double xsq, t; 5024599Szliu double j0(), j1(); 5124599Szliu 5224599Szliu if(n<0){ 5324599Szliu n = -n; 5424599Szliu x = -x; 5524599Szliu } 5624599Szliu if(n==0) return(j0(x)); 5724599Szliu if(n==1) return(j1(x)); 5824599Szliu if(x == 0.) return(0.); 5924599Szliu if(n>x) goto recurs; 6024599Szliu 6124599Szliu a = j0(x); 6224599Szliu b = j1(x); 6324599Szliu for(i=1;i<n;i++){ 6424599Szliu temp = b; 6524599Szliu b = (2.*i/x)*b - a; 6624599Szliu a = temp; 6724599Szliu } 6824599Szliu return(b); 6924599Szliu 7024599Szliu recurs: 7124599Szliu xsq = x*x; 7224599Szliu for(t=0,i=n+16;i>n;i--){ 7324599Szliu t = xsq/(2.*i - t); 7424599Szliu } 7524599Szliu t = x/(2.*n-t); 7624599Szliu 7724599Szliu a = t; 7824599Szliu b = 1; 7924599Szliu for(i=n-1;i>0;i--){ 8024599Szliu temp = b; 8124599Szliu b = (2.*i/x)*b - a; 8224599Szliu a = temp; 8324599Szliu } 8424599Szliu return(t*j0(x)/b); 8524599Szliu } 8624599Szliu 8724599Szliu double 8824599Szliu yn(n,x) int n; double x;{ 8924599Szliu int i; 9024599Szliu int sign; 9124599Szliu double a, b, temp; 9224599Szliu double y0(), y1(); 9324599Szliu 9424599Szliu if (x <= 0) { 95*31853Szliu #if defined(vax)||defined(tahoe) 9624599Szliu extern double infnan(); 9724599Szliu return(infnan(EDOM)); /* NaN */ 98*31853Szliu #else /* defined(vax)||defined(tahoe) */ 9924599Szliu return(zero/zero); /* IEEE machines: invalid operation */ 100*31853Szliu #endif /* defined(vax)||defined(tahoe) */ 10124599Szliu } 10224599Szliu sign = 1; 10324599Szliu if(n<0){ 10424599Szliu n = -n; 10524599Szliu if(n%2 == 1) sign = -1; 10624599Szliu } 10724599Szliu if(n==0) return(y0(x)); 10824599Szliu if(n==1) return(sign*y1(x)); 10924599Szliu 11024599Szliu a = y0(x); 11124599Szliu b = y1(x); 11224599Szliu for(i=1;i<n;i++){ 11324599Szliu temp = b; 11424599Szliu b = (2.*i/x)*b - a; 11524599Szliu a = temp; 11624599Szliu } 11724599Szliu return(sign*b); 11824599Szliu } 119