4.3BSD - /usr/src/usr.lib/libm/jn.c

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

Defined functions

Defined variables