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