projects / unix-history / blame
| Commit | Line | Data |
|---|---|---|
| 0bdde995 SL |
1 | /* @(#)jn.c 4.1 %G% */ |
| 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 | } |