1*34119Sbostic /* 2*34119Sbostic * Copyright (c) 1985 Regents of the University of California. 3*34119Sbostic * All rights reserved. The Berkeley software License Agreement 4*34119Sbostic * specifies the terms and conditions for redistribution. 5*34119Sbostic */ 6*34119Sbostic 724599Szliu #ifndef lint 8*34119Sbostic static char sccsid[] = "@(#)jn.c 5.2 (Berkeley) 04/29/88"; 9*34119Sbostic #endif /* not lint */ 1024599Szliu 1124599Szliu /* 1224599Szliu floating point Bessel's function of 1324599Szliu the first and second kinds and of 1424599Szliu integer order. 1524599Szliu 1624599Szliu int n; 1724599Szliu double x; 1824599Szliu jn(n,x); 1924599Szliu 2024599Szliu returns the value of Jn(x) for all 2124599Szliu integer values of n and all real values 2224599Szliu of x. 2324599Szliu 2424599Szliu There are no error returns. 2524599Szliu Calls j0, j1. 2624599Szliu 2724599Szliu For n=0, j0(x) is called, 2824599Szliu for n=1, j1(x) is called, 2924599Szliu for n<x, forward recursion us used starting 3024599Szliu from values of j0(x) and j1(x). 3124599Szliu for n>x, a continued fraction approximation to 3224599Szliu j(n,x)/j(n-1,x) is evaluated and then backward 3324599Szliu recursion is used starting from a supposed value 3424599Szliu for j(n,x). The resulting value of j(0,x) is 3524599Szliu compared with the actual value to correct the 3624599Szliu supposed value of j(n,x). 3724599Szliu 3824599Szliu yn(n,x) is similar in all respects, except 3924599Szliu that forward recursion is used for all 4024599Szliu values of n>1. 4124599Szliu */ 4224599Szliu 4324599Szliu #include <math.h> 4431853Szliu #if defined(vax)||defined(tahoe) 4524599Szliu #include <errno.h> 4631853Szliu #else /* defined(vax)||defined(tahoe) */ 4724599Szliu static double zero = 0.e0; 4831853Szliu #endif /* defined(vax)||defined(tahoe) */ 4924599Szliu 5024599Szliu double 5124599Szliu jn(n,x) int n; double x;{ 5224599Szliu int i; 5324599Szliu double a, b, temp; 5424599Szliu double xsq, t; 5524599Szliu double j0(), j1(); 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 double y0(), y1(); 9824599Szliu 9924599Szliu if (x <= 0) { 10031853Szliu #if defined(vax)||defined(tahoe) 10124599Szliu extern double infnan(); 10224599Szliu return(infnan(EDOM)); /* NaN */ 10331853Szliu #else /* defined(vax)||defined(tahoe) */ 10424599Szliu return(zero/zero); /* IEEE machines: invalid operation */ 10531853Szliu #endif /* defined(vax)||defined(tahoe) */ 10624599Szliu } 10724599Szliu sign = 1; 10824599Szliu if(n<0){ 10924599Szliu n = -n; 11024599Szliu if(n%2 == 1) sign = -1; 11124599Szliu } 11224599Szliu if(n==0) return(y0(x)); 11324599Szliu if(n==1) return(sign*y1(x)); 11424599Szliu 11524599Szliu a = y0(x); 11624599Szliu b = y1(x); 11724599Szliu for(i=1;i<n;i++){ 11824599Szliu temp = b; 11924599Szliu b = (2.*i/x)*b - a; 12024599Szliu a = temp; 12124599Szliu } 12224599Szliu return(sign*b); 12324599Szliu } 124