ref: 4fedd7d4fe1c3b509c84d12e971aaa466b8b4b3d
dir: /sys/src/ape/lib/ap/math/jn.c/
#include <math.h>
#include <errno.h>
/*
floating point Bessel's function of
the first and second kinds and of
integer order.
int n;
double x;
jn(n,x);
returns the value of Jn(x) for all
integer values of n and all real values
of x.
There are no error returns.
Calls j0, j1.
For n=0, j0(x) is called,
for n=1, j1(x) is called,
for n<x, forward recursion us used starting
from values of j0(x) and j1(x).
for n>x, a continued fraction approximation to
j(n,x)/j(n-1,x) is evaluated and then backward
recursion is used starting from a supposed value
for j(n,x). The resulting value of j(0,x) is
compared with the actual value to correct the
supposed value of j(n,x).
yn(n,x) is similar in all respects, except
that forward recursion is used for all
values of n>1.
*/
double j0(double);
double j1(double);
double y0(double);
double y1(double);
double
jn(int n, double x)
{
int i;
double a, b, temp;
double xsq, t;
if(n < 0) {
n = -n;
x = -x;
}
if(n == 0)
return j0(x);
if(n == 1)
return j1(x);
if(x == 0)
return 0;
if(n > x)
goto recurs;
a = j0(x);
b = j1(x);
for(i=1; i<n; i++) {
temp = b;
b = (2*i/x)*b - a;
a = temp;
}
return b;
recurs:
xsq = x*x;
for(t=0,i=n+16; i>n; i--)
t = xsq/(2*i - t);
t = x/(2*n-t);
a = t;
b = 1;
for(i=n-1; i>0; i--) {
temp = b;
b = (2*i/x)*b - a;
a = temp;
}
return t*j0(x)/b;
}
double
yn(int n, double x)
{
int i;
int sign;
double a, b, temp;
if (x <= 0) {
errno = EDOM;
return -HUGE_VAL;
}
sign = 1;
if(n < 0) {
n = -n;
if(n%2 == 1)
sign = -1;
}
if(n == 0)
return y0(x);
if(n == 1)
return sign*y1(x);
a = y0(x);
b = y1(x);
for(i=1; i<n; i++) {
temp = b;
b = (2*i/x)*b - a;
a = temp;
}
return sign*b;
}