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