1*24599Szliu #ifndef lint 2*24599Szliu static char sccsid[] = "@(#)jn.c 1.1 (ELEFUNT) 09/06/85"; 3*24599Szliu #endif not lint 4*24599Szliu 5*24599Szliu /* 6*24599Szliu floating point Bessel's function of 7*24599Szliu the first and second kinds and of 8*24599Szliu integer order. 9*24599Szliu 10*24599Szliu int n; 11*24599Szliu double x; 12*24599Szliu jn(n,x); 13*24599Szliu 14*24599Szliu returns the value of Jn(x) for all 15*24599Szliu integer values of n and all real values 16*24599Szliu of x. 17*24599Szliu 18*24599Szliu There are no error returns. 19*24599Szliu Calls j0, j1. 20*24599Szliu 21*24599Szliu For n=0, j0(x) is called, 22*24599Szliu for n=1, j1(x) is called, 23*24599Szliu for n<x, forward recursion us used starting 24*24599Szliu from values of j0(x) and j1(x). 25*24599Szliu for n>x, a continued fraction approximation to 26*24599Szliu j(n,x)/j(n-1,x) is evaluated and then backward 27*24599Szliu recursion is used starting from a supposed value 28*24599Szliu for j(n,x). The resulting value of j(0,x) is 29*24599Szliu compared with the actual value to correct the 30*24599Szliu supposed value of j(n,x). 31*24599Szliu 32*24599Szliu yn(n,x) is similar in all respects, except 33*24599Szliu that forward recursion is used for all 34*24599Szliu values of n>1. 35*24599Szliu */ 36*24599Szliu 37*24599Szliu #include <math.h> 38*24599Szliu #ifdef VAX 39*24599Szliu #include <errno.h> 40*24599Szliu #else /* IEEE double */ 41*24599Szliu static double zero = 0.e0; 42*24599Szliu #endif 43*24599Szliu 44*24599Szliu double 45*24599Szliu jn(n,x) int n; double x;{ 46*24599Szliu int i; 47*24599Szliu double a, b, temp; 48*24599Szliu double xsq, t; 49*24599Szliu double j0(), j1(); 50*24599Szliu 51*24599Szliu if(n<0){ 52*24599Szliu n = -n; 53*24599Szliu x = -x; 54*24599Szliu } 55*24599Szliu if(n==0) return(j0(x)); 56*24599Szliu if(n==1) return(j1(x)); 57*24599Szliu if(x == 0.) return(0.); 58*24599Szliu if(n>x) goto recurs; 59*24599Szliu 60*24599Szliu a = j0(x); 61*24599Szliu b = j1(x); 62*24599Szliu for(i=1;i<n;i++){ 63*24599Szliu temp = b; 64*24599Szliu b = (2.*i/x)*b - a; 65*24599Szliu a = temp; 66*24599Szliu } 67*24599Szliu return(b); 68*24599Szliu 69*24599Szliu recurs: 70*24599Szliu xsq = x*x; 71*24599Szliu for(t=0,i=n+16;i>n;i--){ 72*24599Szliu t = xsq/(2.*i - t); 73*24599Szliu } 74*24599Szliu t = x/(2.*n-t); 75*24599Szliu 76*24599Szliu a = t; 77*24599Szliu b = 1; 78*24599Szliu for(i=n-1;i>0;i--){ 79*24599Szliu temp = b; 80*24599Szliu b = (2.*i/x)*b - a; 81*24599Szliu a = temp; 82*24599Szliu } 83*24599Szliu return(t*j0(x)/b); 84*24599Szliu } 85*24599Szliu 86*24599Szliu double 87*24599Szliu yn(n,x) int n; double x;{ 88*24599Szliu int i; 89*24599Szliu int sign; 90*24599Szliu double a, b, temp; 91*24599Szliu double y0(), y1(); 92*24599Szliu 93*24599Szliu if (x <= 0) { 94*24599Szliu #ifdef VAX 95*24599Szliu extern double infnan(); 96*24599Szliu return(infnan(EDOM)); /* NaN */ 97*24599Szliu #else /* IEEE double */ 98*24599Szliu return(zero/zero); /* IEEE machines: invalid operation */ 99*24599Szliu #endif 100*24599Szliu } 101*24599Szliu sign = 1; 102*24599Szliu if(n<0){ 103*24599Szliu n = -n; 104*24599Szliu if(n%2 == 1) sign = -1; 105*24599Szliu } 106*24599Szliu if(n==0) return(y0(x)); 107*24599Szliu if(n==1) return(sign*y1(x)); 108*24599Szliu 109*24599Szliu a = y0(x); 110*24599Szliu b = y1(x); 111*24599Szliu for(i=1;i<n;i++){ 112*24599Szliu temp = b; 113*24599Szliu b = (2.*i/x)*b - a; 114*24599Szliu a = temp; 115*24599Szliu } 116*24599Szliu return(sign*b); 117*24599Szliu } 118