1 /*- 2 * Copyright (c) 1985 The Regents of the University of California. 3 * All rights reserved. 4 * 5 * %sccs.include.proprietary.c% 6 */ 7 8 #ifndef lint 9 static char sccsid[] = "@(#)jn.c 5.4 (Berkeley) 04/20/91"; 10 #endif /* not lint */ 11 12 /* 13 floating point Bessel's function of 14 the first and second kinds and of 15 integer order. 16 17 int n; 18 double x; 19 jn(n,x); 20 21 returns the value of Jn(x) for all 22 integer values of n and all real values 23 of x. 24 25 There are no error returns. 26 Calls j0, j1. 27 28 For n=0, j0(x) is called, 29 for n=1, j1(x) is called, 30 for n<x, forward recursion us used starting 31 from values of j0(x) and j1(x). 32 for n>x, a continued fraction approximation to 33 j(n,x)/j(n-1,x) is evaluated and then backward 34 recursion is used starting from a supposed value 35 for j(n,x). The resulting value of j(0,x) is 36 compared with the actual value to correct the 37 supposed value of j(n,x). 38 39 yn(n,x) is similar in all respects, except 40 that forward recursion is used for all 41 values of n>1. 42 */ 43 44 #include "mathimpl.h" 45 #if defined(vax)||defined(tahoe) 46 #include <errno.h> 47 #else /* defined(vax)||defined(tahoe) */ 48 static double zero = 0.e0; 49 #endif /* defined(vax)||defined(tahoe) */ 50 51 double 52 jn(n,x) int n; double x;{ 53 int i; 54 double a, b, temp; 55 double xsq, t; 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 98 if (x <= 0) { 99 #if defined(vax)||defined(tahoe) 100 return(infnan(EDOM)); /* NaN */ 101 #else /* defined(vax)||defined(tahoe) */ 102 return(zero/zero); /* IEEE machines: invalid operation */ 103 #endif /* defined(vax)||defined(tahoe) */ 104 } 105 sign = 1; 106 if(n<0){ 107 n = -n; 108 if(n%2 == 1) sign = -1; 109 } 110 if(n==0) return(y0(x)); 111 if(n==1) return(sign*y1(x)); 112 113 a = y0(x); 114 b = y1(x); 115 for(i=1;i<n;i++){ 116 temp = b; 117 b = (2.*i/x)*b - a; 118 a = temp; 119 } 120 return(sign*b); 121 } 122