134119Sbostic /* 234119Sbostic * Copyright (c) 1985 Regents of the University of California. 334119Sbostic * All rights reserved. The Berkeley software License Agreement 434119Sbostic * specifies the terms and conditions for redistribution. 534119Sbostic */ 634119Sbostic 724599Szliu #ifndef lint 8*35679Sbostic static char sccsid[] = "@(#)jn.c 5.3 (Berkeley) 09/22/88"; 934119Sbostic #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 43*35679Sbostic #include "mathimpl.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 5624599Szliu if(n<0){ 5724599Szliu n = -n; 5824599Szliu x = -x; 5924599Szliu } 6024599Szliu if(n==0) return(j0(x)); 6124599Szliu if(n==1) return(j1(x)); 6224599Szliu if(x == 0.) return(0.); 6324599Szliu if(n>x) goto recurs; 6424599Szliu 6524599Szliu a = j0(x); 6624599Szliu b = j1(x); 6724599Szliu for(i=1;i<n;i++){ 6824599Szliu temp = b; 6924599Szliu b = (2.*i/x)*b - a; 7024599Szliu a = temp; 7124599Szliu } 7224599Szliu return(b); 7324599Szliu 7424599Szliu recurs: 7524599Szliu xsq = x*x; 7624599Szliu for(t=0,i=n+16;i>n;i--){ 7724599Szliu t = xsq/(2.*i - t); 7824599Szliu } 7924599Szliu t = x/(2.*n-t); 8024599Szliu 8124599Szliu a = t; 8224599Szliu b = 1; 8324599Szliu for(i=n-1;i>0;i--){ 8424599Szliu temp = b; 8524599Szliu b = (2.*i/x)*b - a; 8624599Szliu a = temp; 8724599Szliu } 8824599Szliu return(t*j0(x)/b); 8924599Szliu } 9024599Szliu 9124599Szliu double 9224599Szliu yn(n,x) int n; double x;{ 9324599Szliu int i; 9424599Szliu int sign; 9524599Szliu double a, b, temp; 9624599Szliu 9724599Szliu if (x <= 0) { 9831853Szliu #if defined(vax)||defined(tahoe) 9924599Szliu return(infnan(EDOM)); /* NaN */ 10031853Szliu #else /* defined(vax)||defined(tahoe) */ 10124599Szliu return(zero/zero); /* IEEE machines: invalid operation */ 10231853Szliu #endif /* defined(vax)||defined(tahoe) */ 10324599Szliu } 10424599Szliu sign = 1; 10524599Szliu if(n<0){ 10624599Szliu n = -n; 10724599Szliu if(n%2 == 1) sign = -1; 10824599Szliu } 10924599Szliu if(n==0) return(y0(x)); 11024599Szliu if(n==1) return(sign*y1(x)); 11124599Szliu 11224599Szliu a = y0(x); 11324599Szliu b = y1(x); 11424599Szliu for(i=1;i<n;i++){ 11524599Szliu temp = b; 11624599Szliu b = (2.*i/x)*b - a; 11724599Szliu a = temp; 11824599Szliu } 11924599Szliu return(sign*b); 12024599Szliu } 121