home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / BESSJ.C < prev    next >
Encoding:
Text File  |  1984-08-29  |  2.5 KB  |  93 lines
  1.    bessj(x1,b1,y1)
  2.  
  3.       /* function to find the Bessel J function of integer or  */
  4.       /* fractional order for any x greater than zero.         */
  5.  
  6.       float x1,b1,*y1;
  7.  
  8.     {
  9.  
  10.      int ib1,iflag;
  11.      float c[6],d[6],a,aa,bb,bbo,a1,b,c1,p1,t,t2,term1,x,xx1,y,yy1,z1;
  12.      static float
  13.            s[26]={-3.375,54.4921875,-272.953125,354.375,-2.0625,
  14.                   13.7578125,-19.6875,-1.0625,1.875,-.375,.0030381944,
  15.                   -.4861111111,10.28645833333,-58.,78.75,
  16.                   .0005580357,-.4241071428,3.6026785714,-5.625,.0125,
  17.                   -.35,.75,.0416666666,.25,3.1415926535,1.5707963267};
  18.      extern double cos(),fabs(),pow(),sqrt();
  19.  
  20.  
  21.       *y1 = 0.0;
  22.       if(b1 <  0.0) return;
  23.       if(x1 < 8.0)
  24.       {
  25.        c1 = b1 + 1.;
  26.        gamma(c1,&z1);
  27.        term1 = pow((x1*0.5),b1)/z1;
  28.        p1 = 1.;
  29.        *y1 = term1;
  30.        while((fabs(*y1) * .00000000001) <= fabs(term1))
  31.        {
  32.         term1 = -term1 * x1 * x1 / (4.0 * (b1 + p1) * p1);
  33.         *y1 = *y1 + term1;
  34.         p1 = p1 + 1.;
  35.        }
  36.        return;
  37.       }
  38.       if(b1 <= 1.)
  39.       {
  40.       iflag = 1;
  41.       x = x1;
  42.       b = b1;
  43.       y = *y1;
  44.       }
  45.        else
  46.        {
  47.         ib1 = b1;
  48.         a1 = b1 - ib1;
  49.         iflag = 2;
  50.         x = x1;
  51.         b = a1;
  52.         xx1 = 0.;
  53.         y = xx1;
  54.        }
  55. reit: ;
  56.       a = b * b - 0.25;
  57.       t = 1./x;
  58.       t2 = t * t;
  59.       c[0] = (((s[0] * a + s[1]) * a + s[2]) * a +s[3]) * a;
  60.       c[1] = ((s[4] * a + s[5]) * a + s[6]) * a;
  61.       c[2] = (s[7] * a + s[8]) * a;
  62.       c[3] = s[9] * a;
  63.       bbo = (((c[0] * t2 + c[1]) * t2 + c[2]) * t2 + c[3]) * t2 * t2 + 1.;
  64.       bb = bbo * sqrt(2. * t / (s[24] * sqrt(1. - a * t2)));
  65.       d[0] = ((((s[10] * a + s[11]) * a + s[12]) * a + s[13]) * a + s[14])*a;
  66.       d[1] = (((s[15] * a + s[16]) * a + s[17]) * a + s[18]) * a;
  67.       d[2] = ((s[19] * a + s[20]) * a + s[21]) * a;
  68.       d[3] = (s[22] * a + s[23]) * a;
  69.       d[4] = .5 * a;
  70.       aa = ((((d[0]* t2+ d[1])* t2+ d[2])* t2+ d[3])* t2 +d[4])* t2+1.;
  71.       aa = aa * x - (b + .5) * s[25];
  72.       y = bb * cos(aa);
  73.       if(iflag == 1) return;
  74.       if(iflag == 2)
  75.       {
  76.        p1 = a1 - 1.;
  77.        iflag = 3;
  78.        x = x1;
  79.        b = p1;
  80.        yy1 = 0.;
  81.        y = yy1;
  82.        goto reit;
  83.       }
  84.       while(a1 != b1)
  85.       {
  86.        *y1 = 2. * pow(a1,x1) / x1 - yy1;
  87.        if(a1 == 0) *y1 = 0;
  88.        a1 = a1 + 1.;
  89.        yy1 = xx1;
  90.        xx1 = *y1;
  91.       }
  92.     }
  93.