home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / math / w_jn.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-16  |  2.3 KB  |  89 lines
  1.  
  2. /* @(#)w_jn.c 1.3 95/01/18 */
  3. /*
  4.  * ====================================================
  5.  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  6.  *
  7.  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  8.  * Permission to use, copy, modify, and distribute this
  9.  * software is freely granted, provided that this notice 
  10.  * is preserved.
  11.  * ====================================================
  12.  */
  13.  
  14. /*
  15.  * wrapper jn(int n, double x), yn(int n, double x)
  16.  * floating point Bessel's function of the 1st and 2nd kind
  17.  * of order n
  18.  *          
  19.  * Special cases:
  20.  *    y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  21.  *    y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  22.  * Note 2. About jn(n,x), yn(n,x)
  23.  *    For n=0, j0(x) is called,
  24.  *    for n=1, j1(x) is called,
  25.  *    for n<x, forward recursion us used starting
  26.  *    from values of j0(x) and j1(x).
  27.  *    for n>x, a continued fraction approximation to
  28.  *    j(n,x)/j(n-1,x) is evaluated and then backward
  29.  *    recursion is used starting from a supposed value
  30.  *    for j(n,x). The resulting value of j(0,x) is
  31.  *    compared with the actual value to correct the
  32.  *    supposed value of j(n,x).
  33.  *
  34.  *    yn(n,x) is similar in all respects, except
  35.  *    that forward recursion is used for all
  36.  *    values of n>1.
  37.  *    
  38.  */
  39.  
  40. #include "fdlibm.h"
  41.  
  42. #ifdef __STDC__
  43.     double jn(int n, double x)    /* wrapper jn */
  44. #else
  45.     double jn(n,x)            /* wrapper jn */
  46.     double x; int n;
  47. #endif
  48. {
  49. #ifdef _IEEE_LIBM
  50.     return __ieee754_jn(n,x);
  51. #else
  52.     double z;
  53.     z = __ieee754_jn(n,x);
  54.     if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
  55.     if(fabs(x)>X_TLOSS) {
  56.         return __kernel_standard((double)n,x,38); /* jn(|x|>X_TLOSS,n) */
  57.     } else
  58.         return z;
  59. #endif
  60. }
  61.  
  62. #ifdef __STDC__
  63.     double yn(int n, double x)    /* wrapper yn */
  64. #else
  65.     double yn(n,x)            /* wrapper yn */
  66.     double x; int n;
  67. #endif
  68. {
  69. #ifdef _IEEE_LIBM
  70.     return __ieee754_yn(n,x);
  71. #else
  72.     double z;
  73.     z = __ieee754_yn(n,x);
  74.     if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z;
  75.         if(x <= 0.0){
  76.                 if(x==0.0)
  77.                     /* d= -one/(x-x); */
  78.                     return __kernel_standard((double)n,x,12);
  79.                 else
  80.                     /* d = zero/(x-x); */
  81.                     return __kernel_standard((double)n,x,13);
  82.         }
  83.     if(x>X_TLOSS) {
  84.         return __kernel_standard((double)n,x,39); /* yn(x>X_TLOSS,n) */
  85.     } else
  86.         return z;
  87. #endif
  88. }
  89.