1*48402Sbostic /*- 2*48402Sbostic * Copyright (c) 1985 The Regents of the University of California. 3*48402Sbostic * All rights reserved. 4*48402Sbostic * 5*48402Sbostic * %sccs.include.proprietary.c% 634119Sbostic */ 734119Sbostic 824599Szliu #ifndef lint 9*48402Sbostic static char sccsid[] = "@(#)jn.c 5.4 (Berkeley) 04/20/91"; 1034119Sbostic #endif /* not lint */ 1124599Szliu 1224599Szliu /* 1324599Szliu floating point Bessel's function of 1424599Szliu the first and second kinds and of 1524599Szliu integer order. 1624599Szliu 1724599Szliu int n; 1824599Szliu double x; 1924599Szliu jn(n,x); 2024599Szliu 2124599Szliu returns the value of Jn(x) for all 2224599Szliu integer values of n and all real values 2324599Szliu of x. 2424599Szliu 2524599Szliu There are no error returns. 2624599Szliu Calls j0, j1. 2724599Szliu 2824599Szliu For n=0, j0(x) is called, 2924599Szliu for n=1, j1(x) is called, 3024599Szliu for n<x, forward recursion us used starting 3124599Szliu from values of j0(x) and j1(x). 3224599Szliu for n>x, a continued fraction approximation to 3324599Szliu j(n,x)/j(n-1,x) is evaluated and then backward 3424599Szliu recursion is used starting from a supposed value 3524599Szliu for j(n,x). The resulting value of j(0,x) is 3624599Szliu compared with the actual value to correct the 3724599Szliu supposed value of j(n,x). 3824599Szliu 3924599Szliu yn(n,x) is similar in all respects, except 4024599Szliu that forward recursion is used for all 4124599Szliu values of n>1. 4224599Szliu */ 4324599Szliu 4435679Sbostic #include "mathimpl.h" 4531853Szliu #if defined(vax)||defined(tahoe) 4624599Szliu #include <errno.h> 4731853Szliu #else /* defined(vax)||defined(tahoe) */ 4824599Szliu static double zero = 0.e0; 4931853Szliu #endif /* defined(vax)||defined(tahoe) */ 5024599Szliu 5124599Szliu double 5224599Szliu jn(n,x) int n; double x;{ 5324599Szliu int i; 5424599Szliu double a, b, temp; 5524599Szliu double xsq, t; 5624599Szliu 5724599Szliu if(n<0){ 5824599Szliu n = -n; 5924599Szliu x = -x; 6024599Szliu } 6124599Szliu if(n==0) return(j0(x)); 6224599Szliu if(n==1) return(j1(x)); 6324599Szliu if(x == 0.) return(0.); 6424599Szliu if(n>x) goto recurs; 6524599Szliu 6624599Szliu a = j0(x); 6724599Szliu b = j1(x); 6824599Szliu for(i=1;i<n;i++){ 6924599Szliu temp = b; 7024599Szliu b = (2.*i/x)*b - a; 7124599Szliu a = temp; 7224599Szliu } 7324599Szliu return(b); 7424599Szliu 7524599Szliu recurs: 7624599Szliu xsq = x*x; 7724599Szliu for(t=0,i=n+16;i>n;i--){ 7824599Szliu t = xsq/(2.*i - t); 7924599Szliu } 8024599Szliu t = x/(2.*n-t); 8124599Szliu 8224599Szliu a = t; 8324599Szliu b = 1; 8424599Szliu for(i=n-1;i>0;i--){ 8524599Szliu temp = b; 8624599Szliu b = (2.*i/x)*b - a; 8724599Szliu a = temp; 8824599Szliu } 8924599Szliu return(t*j0(x)/b); 9024599Szliu } 9124599Szliu 9224599Szliu double 9324599Szliu yn(n,x) int n; double x;{ 9424599Szliu int i; 9524599Szliu int sign; 9624599Szliu double a, b, temp; 9724599Szliu 9824599Szliu if (x <= 0) { 9931853Szliu #if defined(vax)||defined(tahoe) 10024599Szliu return(infnan(EDOM)); /* NaN */ 10131853Szliu #else /* defined(vax)||defined(tahoe) */ 10224599Szliu return(zero/zero); /* IEEE machines: invalid operation */ 10331853Szliu #endif /* defined(vax)||defined(tahoe) */ 10424599Szliu } 10524599Szliu sign = 1; 10624599Szliu if(n<0){ 10724599Szliu n = -n; 10824599Szliu if(n%2 == 1) sign = -1; 10924599Szliu } 11024599Szliu if(n==0) return(y0(x)); 11124599Szliu if(n==1) return(sign*y1(x)); 11224599Szliu 11324599Szliu a = y0(x); 11424599Szliu b = y1(x); 11524599Szliu for(i=1;i<n;i++){ 11624599Szliu temp = b; 11724599Szliu b = (2.*i/x)*b - a; 11824599Szliu a = temp; 11924599Szliu } 12024599Szliu return(sign*b); 12124599Szliu } 122