sys/src/ape/lib/ap/math/jn.c

   1 #include <math.h>
   2 #include <errno.h>
   3
   4 /*
   5         floating point Bessel's function of
   6         the first and second kinds and of
   7         integer order.
   8
   9         int n;
  10         double x;
  11         jn(n,x);
  12
  13         returns the value of Jn(x) for all
  14         integer values of n and all real values
  15         of x.
  16
  17         There are no error returns.
  18         Calls j0, j1.
  19
  20         For n=0, j0(x) is called,
  21         for n=1, j1(x) is called,
  22         for n<x, forward recursion us used starting
  23         from values of j0(x) and j1(x).
  24         for n>x, a continued fraction approximation to
  25         j(n,x)/j(n-1,x) is evaluated and then backward
  26         recursion is used starting from a supposed value
  27         for j(n,x). The resulting value of j(0,x) is
  28         compared with the actual value to correct the
  29         supposed value of j(n,x).
  30
  31         yn(n,x) is similar in all respects, except
  32         that forward recursion is used for all
  33         values of n>1.
  34 */
  35
  36 double  j0(double);
  37 double  j1(double);
  38 double  y0(double);
  39 double  y1(double);
  40
  41 double
  42 jn(int n, double x)
  43 {
  44         int i;
  45         double a, b, temp;
  46         double xsq, t;
  47
  48         if(n < 0) {
  49                 n = -n;
  50                 x = -x;
  51         }
  52         if(n == 0)
  53                 return j0(x);
  54         if(n == 1)
  55                 return j1(x);
  56         if(x == 0)
  57                 return 0;
  58         if(n > x)
  59                 goto recurs;
  60
  61         a = j0(x);
  62         b = j1(x);
  63         for(i=1; i<n; i++) {
  64                 temp = b;
  65                 b = (2*i/x)*b - a;
  66                 a = temp;
  67         }
  68         return b;
  69
  70 recurs:
  71         xsq = x*x;
  72         for(t=0,i=n+16; i>n; i--)
  73                 t = xsq/(2*i - t);
  74         t = x/(2*n-t);
  75
  76         a = t;
  77         b = 1;
  78         for(i=n-1; i>0; i--) {
  79                 temp = b;
  80                 b = (2*i/x)*b - a;
  81                 a = temp;
  82         }
  83         return t*j0(x)/b;
  84 }
  85
  86 double
  87 yn(int n, double x)
  88 {
  89         int i;
  90         int sign;
  91         double a, b, temp;
  92
  93         if (x <= 0) {
  94                 errno = EDOM;
  95                 return -HUGE_VAL;
  96         }
  97         sign = 1;
  98         if(n < 0) {
  99                 n = -n;
 100                 if(n%2 == 1)
 101                         sign = -1;
 102         }
 103         if(n == 0)
 104                 return y0(x);
 105         if(n == 1)
 106                 return sign*y1(x);
 107
 108         a = y0(x);
 109         b = y1(x);
 110         for(i=1; i<n; i++) {
 111                 temp = b;
 112                 b = (2*i/x)*b - a;
 113                 a = temp;
 114         }
 115         return sign*b;
 116 }