1 #include <math.h> 2 #include <errno.h> 3 4 /* 5 floating point Bessel's function of 6 the first and second kinds and of 7 integer order. 8 9 int n; 10 double x; 11 jn(n,x); 12 13 returns the value of Jn(x) for all 14 integer values of n and all real values 15 of x. 16 17 There are no error returns. 18 Calls j0, j1. 19 20 For n=0, j0(x) is called, 21 for n=1, j1(x) is called, 22 for n<x, forward recursion us used starting 23 from values of j0(x) and j1(x). 24 for n>x, a continued fraction approximation to 25 j(n,x)/j(n-1,x) is evaluated and then backward 26 recursion is used starting from a supposed value 27 for j(n,x). The resulting value of j(0,x) is 28 compared with the actual value to correct the 29 supposed value of j(n,x). 30 31 yn(n,x) is similar in all respects, except 32 that forward recursion is used for all 33 values of n>1. 34 */ 35 36 double j0(double); 37 double j1(double); 38 double y0(double); 39 double y1(double); 40 41 double 42 jn(int n, double x) 43 { 44 int i; 45 double a, b, temp; 46 double xsq, t; 47 48 if(n < 0) { 49 n = -n; 50 x = -x; 51 } 52 if(n == 0) 53 return j0(x); 54 if(n == 1) 55 return j1(x); 56 if(x == 0) 57 return 0; 58 if(n > x) 59 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 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(int n, double x) 88 { 89 int i; 90 int sign; 91 double a, b, temp; 92 93 if (x <= 0) { 94 errno = EDOM; 95 return -HUGE_VAL; 96 } 97 sign = 1; 98 if(n < 0) { 99 n = -n; 100 if(n%2 == 1) 101 sign = -1; 102 } 103 if(n == 0) 104 return y0(x); 105 if(n == 1) 106 return sign*y1(x); 107 108 a = y0(x); 109 b = y1(x); 110 for(i=1; i<n; i++) { 111 temp = b; 112 b = (2*i/x)*b - a; 113 a = temp; 114 } 115 return sign*b; 116 } 117