/* @(#)jn.c 4.1 12/25/82 */ /* floating point Bessel's function of the first and second kinds and of integer order. int n; double x; jn(n,x); returns the value of Jn(x) for all integer values of n and all real values of x. There are no error returns. Calls j0, j1. For n=0, j0(x) is called, for n=1, j1(x) is called, for nx, a continued fraction approximation to j(n,x)/j(n-1,x) is evaluated and then backward recursion is used starting from a supposed value for j(n,x). The resulting value of j(0,x) is compared with the actual value to correct the supposed value of j(n,x). yn(n,x) is similar in all respects, except that forward recursion is used for all values of n>1. */ #include #include int errno; double jn(n,x) int n; double x;{ int i; double a, b, temp; double xsq, t; double j0(), j1(); if(n<0){ n = -n; x = -x; } if(n==0) return(j0(x)); if(n==1) return(j1(x)); if(x == 0.) return(0.); if(n>x) goto recurs; a = j0(x); b = j1(x); for(i=1;in;i--){ t = xsq/(2.*i - t); } t = x/(2.*n-t); a = t; b = 1; for(i=n-1;i>0;i--){ temp = b; b = (2.*i/x)*b - a; a = temp; } return(t*j0(x)/b); } double yn(n,x) int n; double x;{ int i; int sign; double a, b, temp; double y0(), y1(); if (x <= 0) { errno = EDOM; return(-HUGE); } sign = 1; if(n<0){ n = -n; if(n%2 == 1) sign = -1; } if(n==0) return(y0(x)); if(n==1) return(sign*y1(x)); a = y0(x); b = y1(x); for(i=1;i