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.3 (Berkeley) 09/22/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 "mathimpl.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 56 if(n<0){ 57 n = -n; 58 x = -x; 59 } 60 if(n==0) return(j0(x)); 61 if(n==1) return(j1(x)); 62 if(x == 0.) return(0.); 63 if(n>x) goto recurs; 64 65 a = j0(x); 66 b = j1(x); 67 for(i=1;i<n;i++){ 68 temp = b; 69 b = (2.*i/x)*b - a; 70 a = temp; 71 } 72 return(b); 73 74 recurs: 75 xsq = x*x; 76 for(t=0,i=n+16;i>n;i--){ 77 t = xsq/(2.*i - t); 78 } 79 t = x/(2.*n-t); 80 81 a = t; 82 b = 1; 83 for(i=n-1;i>0;i--){ 84 temp = b; 85 b = (2.*i/x)*b - a; 86 a = temp; 87 } 88 return(t*j0(x)/b); 89 } 90 91 double 92 yn(n,x) int n; double x;{ 93 int i; 94 int sign; 95 double a, b, temp; 96 97 if (x <= 0) { 98 #if defined(vax)||defined(tahoe) 99 return(infnan(EDOM)); /* NaN */ 100 #else /* defined(vax)||defined(tahoe) */ 101 return(zero/zero); /* IEEE machines: invalid operation */ 102 #endif /* defined(vax)||defined(tahoe) */ 103 } 104 sign = 1; 105 if(n<0){ 106 n = -n; 107 if(n%2 == 1) sign = -1; 108 } 109 if(n==0) return(y0(x)); 110 if(n==1) return(sign*y1(x)); 111 112 a = y0(x); 113 b = y1(x); 114 for(i=1;i<n;i++){ 115 temp = b; 116 b = (2.*i/x)*b - a; 117 a = temp; 118 } 119 return(sign*b); 120 } 121