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