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