2.9BSD - /usr/src/lib/m/jn.c

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

Defined functions

Defined variables