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

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

Defined functions

Defined variables