home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / BESSELJ.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  10.4 KB  |  370 lines
  1. function b = besselj(alpha,xx)
  2. %BESSELJ Bessel functions of the first kind.
  3. %    J = BESSELJ(ALPHA,X) computes Bessel functions of the first kind,
  4. %    J_sub_alpha(X) for real, non-negative order ALPHA and argument X.
  5. %    In general, both ALPHA and X may be vectors.  The output J is
  6. %    an m-by-n matrix with m = lenth(X), n = length(ALPHA) and
  7. %        J(i,k) = J_sub_alpha(k)(X(i)).
  8. %    The elements of X can be any nonnegative real values in any order.
  9. %    For ALPHA, however, there are two important restrictions: the
  10. %    increment in ALPHA must be one, i.e. ALPHA = alpha:1:alpha+n-1,
  11. %    and the values must satisfy 0 <= alpha(k) <= 1000.
  12. %
  13. %    Examples:
  14. %
  15. %        besselj(3:9,(10:.2:20)') generates the 51-by-7 table on page 400
  16. %        of Abramowitz and Stegun, "Handbook of Mathematical Functions."
  17. %
  18. %        besselj(2/3:1:98/3,r) generates the fractional order Bessel
  19. %        functions used by the MathWorks Logo, the L-shaped membrane.
  20. %        J_sub_2/3(r) matches the singularity at the interior corner
  21. %        where the angle is pi/(2/3).
  22. %
  23. %    See also: BESSELY, BESSELI, BESSELK.
  24.  
  25. %   Acknowledgement:
  26. %
  27. %    This program is based on a Fortran program by W. J. Cody,
  28. %    Applied Mathematics Division, Argonne National Laboratory, 
  29. %    dated March 19, 1990.  Cody references earlier work by
  30. %    David J. Sookne and F. W. J. Olver:
  31. %
  32. %    "A Note on Backward Recurrence Algorithms," Olver, F. W. J.,
  33. %    and Sookne, D. J., Math. Comp. 26, 1972, pp 941-947.
  34. %
  35. %    "Bessel Functions of Real Argument and Integer Order,"
  36. %    Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp 125-132.
  37. %
  38. %    Note:
  39. %
  40. %    This version supercedes earlier versions described in the reference
  41. %    manuals for MATLAB 3.5 and 4.0.  Earlier usage is still acceptable,
  42. %    but has been generalized to allow vector ALPHA.  The algorithm has
  43. %    been substantially improved.  The auxilliary routines BESSELN and
  44. %    BESSELA are not longer used.  BESSELH has been superceded by BESSELY.
  45. %
  46. %    MATLAB version: C. Moler, 10/2/92.
  47. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  48.  
  49. %
  50. %  Check for real, non-negative arguments.
  51. %
  52.    if any(imag(xx)) | any(xx < 0) | any(imag(alpha)) | any(alpha < 0)
  53.       error('Input arguments must be real and nonnegative.')
  54.    end
  55.    if isempty(alpha) | isempty(xx)
  56.       bk = []; 
  57.       return
  58.    end
  59. %
  60. %  Break alpha into integer and fractional parts,
  61. %  and initialize result array.
  62. %
  63.    nfirst = fix(alpha(1));
  64.    nb = fix(alpha(length(alpha))) + 1;
  65.    if ~(nb <= 1001)
  66.       error('Alpha must be <= 1000.')
  67.    end
  68.    if length(alpha) > 1
  69.       if any(abs(diff(alpha)-1) > 4*nb*eps)
  70.          error('Increment in alpha must be 1.')
  71.       end
  72.    end
  73.    resize = (length(alpha) == 1);
  74.    if resize, resize = size(xx); end
  75.    xx = xx(:);
  76.    b = NaN*ones(length(xx),nb);
  77.    alpha = alpha(1) - nfirst;
  78. %
  79. %  Two-term ascending series for small x.
  80. %
  81.    v = find(xx < 1.e-4);
  82.    if any(v)
  83.       x = xx(v);
  84.       tempa = ones(size(x));
  85.       alpem = 1 + alpha;
  86.       halfx = 0.5*x;
  87.       if (alpha ~= 0)
  88.          tempa = halfx.^alpha/(alpha*gamma(alpha));
  89.       end
  90.       tempb = -halfx.*halfx;
  91.       b(v,1) = tempa + tempa.*tempb/alpem;
  92.       if (nb ~= 1)
  93. %
  94. %        Calculate higher order functions.
  95. %
  96.          tempc = halfx;
  97.          for n = 2:nb
  98.             tempa = tempa/alpem;
  99.             alpem = alpem + 1;
  100.             tempa = tempa.*tempc;
  101.             b(v,n) = tempa + tempa.*tempb/alpem;
  102.          end
  103.       end
  104.    end
  105. %
  106. %  Aymptotic series for large x when nb is not too large.
  107. %
  108.    v = find((xx > 25) & (nb <= xx+1));
  109.    if any(v)
  110.       x = xx(v);
  111.       xc = sqrt(2./(pi*x));
  112.       xin = (1./(8*x)).^2;
  113.       m = 11;
  114.       if all(x >= 35), m = 8; end
  115.       if all(x >= 130), m = 4; end
  116.       xm = 4*m;
  117. %
  118. %     Argument reduction for sin and cos routines.
  119. %     twopi1 + twopi2 = 2 * pi to extra precision.
  120. %
  121.       twopi1 = 201/32;
  122.       twopi2 = 0.001935307179586476925286767;
  123.       t = round(x/(2*pi));
  124.       z = (x-t*twopi1) - t*twopi2 - (alpha+0.5)*pi/2;
  125.       vsin = sin(z);
  126.       vcos = cos(z);
  127.       gnu = alpha + alpha;
  128.       for i = 1:2
  129.          s = ((xm-1)-gnu)*((xm-1)+gnu)*xin*0.5;
  130.          t = (gnu-(xm-3)).*(gnu+(xm-3));
  131.          capp = s.*t/prod(1:2*m);
  132.          t1 = (gnu-(xm+1)).*(gnu+(xm+1));
  133.          capq = s.*t1/prod(1:2*m+1);
  134.          xk = xm;
  135.          k = m + m;
  136.          t1 = t;
  137.          for j = 2:m
  138.             xk = xk - 4;
  139.             s = ((xk-1)-gnu).*((xk-1)+gnu);
  140.             t = (gnu-(xk-3)).*(gnu+(xk-3));
  141.             capp = (capp+1/prod(1:k-2)).*s.*t.*xin;
  142.             capq = (capq+1/prod(1:k-1)).*s.*t1.*xin;
  143.             k = k - 2;
  144.             t1 = t;
  145.          end
  146.          capp = capp + 1;
  147.          capq = (capq+1).*(gnu*gnu-1).*(1./(8*x));
  148.          b(v,i) = xc.*(capp.*vcos-capq.*vsin);
  149.          if (nb == 1), break, end
  150.          t = vsin;
  151.          vsin = -vcos;
  152.          vcos = t;
  153.          gnu = gnu + 2;
  154.       end
  155. %
  156.       gnu = alpha + alpha + 2;
  157.       for j = 3:nb
  158.          b(v,j) = gnu*b(v,j-1)./x - b(v,j-2);
  159.          gnu = gnu + 2;
  160.       end
  161.    end
  162. %
  163. %  For most x, use three-term recurrence.
  164. %
  165.    v = find((xx >= 1.e-4) & ((xx <= 25) | (nb > xx+1)));
  166.    if any(v)
  167.       x = xx(v);
  168.       ncalc = 0;
  169.       magx = max(x);
  170.       nbmx = nb - floor(magx);
  171.       n = floor(magx) + 1;
  172.       en = n+n + (alpha+alpha);
  173.       plast = 1;
  174.       p = en/magx;
  175. %
  176. %     Calculate general significance test.
  177. %
  178.       test = 2/eps;
  179.       if (nbmx >= 3)
  180. %
  181. %     Calculate p*s until n = nb-1.  Check for possible overflow.
  182. %
  183.          tover = eps*realmax;
  184.          nstart = floor(magx) + 2;
  185.          nend = nb - 1;
  186.          en = nstart+nstart - 2 + (alpha+alpha);
  187.          for k = nstart:nend
  188.             n = k;
  189.             en = en + 2;
  190.             pold = plast;
  191.             plast = p;
  192.             p = en*plast/magx - pold;
  193.             if p > tover
  194. %
  195. %           To avoid overflow, divide p*s by tover.
  196. %           Calculate p*s until abs(p) > 1.
  197. %
  198.                tover = realmax;
  199.                p = p/tover;
  200.                plast = plast/tover;
  201.                psave = p;
  202.                psavel = plast;
  203.                nstart = n + 1;
  204.                while p <= 1
  205.                   n = n + 1;
  206.                   en = en + 2;
  207.                   pold = plast;
  208.                   plast = p;
  209.                   p = en*plast/magx - pold;
  210.                end
  211.                tempb = en/magx;
  212. %
  213. %              Calculate backward test and find ncalc, the highest n
  214. %              such that the test is passed.
  215. %
  216.                test = pold*plast*(0.5-0.5/(tempb*tempb));
  217.                test = eps*test;
  218.                p = plast*tover;
  219.                n = n - 1;
  220.                en = en - 2;
  221.                nend = min(nb,n);
  222.                ncalc = nend;
  223.                for l = nstart:nend
  224.                   pold = psavel;
  225.                   psavel = psave;
  226.                   psave = en*psavel/magx - pold;
  227.                   if psave*psavel > test
  228.                      ncalc = l - 1;
  229.                      break
  230.                   end
  231.                end
  232.                break
  233.             end
  234.          end
  235.          if ncalc == 0
  236.             n = nend;
  237.             en = n+n + (alpha+alpha);
  238. %
  239. %           Calculate special significance test for nbmx > 2.
  240. %
  241.             test = max(test,max(sqrt(plast/eps)*sqrt(p+p)));
  242.          end
  243.       end
  244. %
  245. %     Calculate p*s until significance test passes.
  246. %
  247.       if ncalc == 0
  248.          ncalc = nb;
  249.          while p < test
  250.             n = n + 1;
  251.             en = en + 2;
  252.             pold = plast;
  253.             plast = p;
  254.             p = en*plast/magx - pold;
  255.          end
  256.       end
  257. %
  258. %     Initialize the backward recursion and the normalization sum.
  259. %
  260.       n = n + 1;
  261.       en = en + 2;
  262.       tempb = zeros(size(x));
  263.       tempa = ones(size(x))/p;
  264.       m = 2*n - 4*fix(n/2);
  265.       sum = zeros(size(x));
  266.       em = fix(n/2);
  267.       alpem = (em-1) + alpha;
  268.       alp2em = (em+em) + alpha;
  269.       if (m ~= 0), sum = tempa*alpem*alp2em/em; end
  270.       nend = n - nb;
  271.       if (nend > 0)
  272. %
  273. %     Recur backward via difference equation, calculating (but not
  274. %     storing) b(:,n), until n = nb.
  275. %
  276.          for l = 1:nend
  277.             n = n - 1;
  278.             en = en - 2;
  279.             tempc = tempb;
  280.             tempb = tempa;
  281.             tempa = (en*tempb)./x - tempc;
  282.             m = 2 - m;
  283.             if (m ~= 0) 
  284.                em = em - 1;
  285.                alp2em = (em+em) + alpha;
  286.                if (n == 1), break, end
  287.                alpem = (em-1) + alpha;
  288.                if (alpem == 0), alpem = 1; end
  289.                sum = (sum+tempa*alp2em)*alpem/em;
  290.             end
  291.          end
  292.       end
  293. %
  294. %     Store b(:,nb).
  295. %
  296.       b(v,n) = tempa;
  297.       if (nend >= 0)
  298.          if (nb <= 1)
  299.             alp2em = alpha;
  300.             if ((alpha+1) == 1), alp2em = 1; end
  301.             sum = sum + b(v,1)*alp2em;
  302.          else
  303. %
  304. %        Calculate and store b(:,nb-1).
  305. %
  306.             n = n - 1;
  307.             en = en - 2;
  308.             b(v,n) = (en*tempa)./x - tempb;
  309.             if (n == 1)
  310.                m = 2;
  311.             end
  312.             m = 2 - m;
  313.             if (m ~= 0)
  314.                em = em - 1;
  315.                alp2em = (em+em) + alpha;
  316.                alpem = (em-1) + alpha;
  317.                if (alpem == 0), alpem = 1; end
  318.                sum = (sum+b(v,n)*alp2em)*alpem/em;
  319.             end
  320.          end
  321.       end
  322.       nend = n - 2;
  323.       if (nend > 0)
  324. %
  325. %     Calculate via difference equation and store b(:,n), until n = 2.
  326. %
  327.          for l = 1:nend
  328.             n = n - 1;
  329.             en = en - 2;
  330.             b(v,n) = (en*b(v,n+1))./x - b(v,n+2);
  331.             m = 2 - m;
  332.             if (m ~= 0)
  333.                em = em - 1;
  334.                alp2em = (em+em) + alpha;
  335.                alpem = (em-1) + alpha;
  336.                if (alpem == 0), alpem = 1; end
  337.                sum = (sum+b(v,n)*alp2em)*alpem/em;
  338.             end
  339.          end
  340.       end
  341. %
  342. %     Calculate b(:,1), i necessary.
  343. %
  344.       if nb > 1, 
  345.          if n > 1
  346.             b(v,1) = 2*(alpha+1)*b(v,2)./x - b(v,3);
  347.          end
  348.          em = em - 1;
  349.          alp2em = (em+em) + alpha;
  350.          if (alp2em == 0), alp2em = 1; end
  351.          sum = sum + b(v,1)*alp2em;
  352.       end
  353. %
  354. %     Normalize.  Divide all b(:,n) by sum.
  355. %
  356.       if alpha > eps
  357.          sum = gamma(alpha)*sum.*(x*0.5).^(-alpha);
  358.       end
  359.       for n = 1:nb
  360.          b(v,n) = b(v,n)./sum;
  361.       end
  362.    end
  363. %
  364. %  Return the requested index range
  365. %
  366.    b = b(:,nfirst+1:nb);
  367.    if resize
  368.       b = reshape(b,resize(1),resize(2));
  369.    end
  370.