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