home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / BESSELI.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  9.3 KB  |  330 lines
  1. function b = besseli(alpha,xx,scale)
  2. %BESSELI Modified Bessel functions of the first kind.
  3. %    I = BESSELI(ALPHA,X) computes modified Bessel functions of the
  4. %    first kind, I_sub_alpha(X) for real, non-negative order ALPHA
  5. %    and argument X.  In general, both ALPHA and X may be vectors.
  6. %    The output I is m-by-n with m = lenth(X), n = length(ALPHA) and
  7. %        I(k,j) = I_sub_alpha(j)(X(k)).
  8. %    The elements of X can be any nonnegative real values in any order.
  9. %    For ALPHA there are two important restrictions: the increment
  10. %    in ALPHA must be one, i.e. ALPHA = alpha:1:alpha+n-1, and the
  11. %    elements must also be in the range 0 <= ALPHA(j) <= 1000.
  12. %
  13. %    E = BESSELI(ALPHA,X,1) computes I_sub_alpha(X)*EXP(-X).
  14. %
  15. %    The relationship between the Bessel functions of the first kind is
  16. %
  17. %        I_sub_alpha(x) = i^(-alpha) * J_sub_alpha(i*x)
  18. %
  19. %    Examples:
  20. %
  21. %        besseli(3:9,[0:.2:9.8 10:.5:20],1) generates the entire 
  22. %        71-by-7 table on page 423 of Abramowitz and Stegun,
  23. %        "Handbook of Mathematical Functions."
  24. %
  25. %    See also: BESSELJ, BESSELY, BESSELK.
  26.  
  27. %   Acknowledgement:
  28. %
  29. %    This program is based on a Fortran program by W. J. Cody and
  30. %    L. Stoltz, Applied Mathematics Division, Argonne National
  31. %    Laboratory, dated May 30, 1989.  Their references include:
  32. %
  33. %       "A Note on backward recurrence algorithms," Olver, F. W. J.,
  34. %       and Sookne, D. J., Math. Comp. 26, 1972, pp 941-947.
  35. %
  36. %       "Bessel functions of real argument and integer order,"
  37. %       Sookne, D. J., NBS Jour. of Res. B. 77b, 1973, pp. 125-132.
  38. %
  39. %       "Algorithm 597, Sequence of modified Bessel functions of the
  40. %       "first kind," Cody, W. J., Trans. Math. Soft., 1983, pp. 242-245.
  41. %
  42. %    MATLAB version: C. Moler, 10/6/92.
  43. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  44.  
  45. %
  46. %  Check for real, non-negative arguments.
  47. %
  48.    if nargin < 3, scale = 0; end
  49.    if any(imag(xx)) | any(xx < 0) | any(imag(alpha)) | any(alpha < 0)
  50.       error('Input arguments must be real and nonnegative.')
  51.    end
  52.    if isempty(alpha) | isempty(xx)
  53.       bk = []; 
  54.       return
  55.    end
  56. %
  57. %  Break alpha into integer and fractional parts,
  58. %  and initialize result array.
  59. %
  60.    nfirst = fix(alpha(1));
  61.    nb = fix(alpha(length(alpha))) + 1;
  62.    if ~(nb <= 1001)
  63.       error('Alpha must be <= 1000.')
  64.    end
  65.    if length(alpha) > 1
  66.       if any(abs(diff(alpha)-1) > 4*nb*eps)
  67.          error('Increment in alpha must be 1.')
  68.       end
  69.    end
  70.    resize = (length(alpha) == 1);
  71.    if resize, resize = size(xx); end
  72.    xx = xx(:);
  73.    b = NaN*ones(length(xx),nb);
  74.    alpha = alpha(1) - nfirst;
  75. %
  76. %  Asymptotic expansion for large x.
  77. %
  78.    v = find(1.e3 < xx);
  79.    if any(v)
  80.       x = xx(v);
  81.       if scale
  82.          z = 1./(8*x);
  83.          mu = 4*(nfirst+alpha:1:nb-1+alpha).^2;
  84.          z = (1 - z*(mu-1).*(1 - (z/2)*(mu-9).*(1 - (z/3)*(mu - 25))));
  85.          b(v,nfirst+1:nb) = ((1./sqrt(2*pi*x))*ones(1,nb-nfirst)).*z;
  86.       else
  87.          b(v,:) = Inf*ones(length(v),nb);
  88.       end
  89.    end
  90. %
  91. %  Three-term backward recurrence for most x.
  92. %
  93.    v = find((1.e-4 <= xx) & (xx <= 1.e3));
  94.    if any(v)
  95.       x = xx(v);
  96.       magx = max(fix(x));
  97. %
  98. %     Initialize the forward sweep, the p-sequence of Olver.
  99. %
  100.       nbmx = nb-magx;
  101.       n = magx+1;
  102.       en = (n+n) + (alpha+alpha);
  103.       plast = 1;
  104.       p = en ./ x;
  105. %
  106. %     Calculate general significance test.
  107. %
  108.       test = 2/eps;
  109.       if (2*magx > 80);
  110.          test = sqrt(test*max(p));
  111.       else
  112.          test = test / 1.585^magx;
  113.       end
  114.       if (nbmx >= 3)
  115. %
  116. %     Calculate p-sequence until n = nb-1.
  117. %
  118.          tover = eps*realmax;
  119.          nstart = magx+2;
  120.          nend = nb - 1;
  121.          for k = nstart:nend
  122.             n = k;
  123.             en = en + 2;
  124.             pold = plast;
  125.             plast = p;
  126.             p = en * plast./x + pold;
  127.             if any(p > tover)
  128. %
  129. %           To avoid overflow, divide p-sequence by tover.
  130. %           Calculate p-sequence until abs(p) > 1.
  131. %
  132.                tover = realmax;
  133.                p = p / tover;
  134.                plast = plast / tover;
  135.                psave = p;
  136.                psavel = plast;
  137.                nstart = n + 1;
  138.                while any(p <= 1);
  139.                   n = n + 1;
  140.                   en = en + 2;
  141.                   pold = plast;
  142.                   plast = p;
  143.                   p = en * plast./x + pold;
  144.                end
  145.                tempb = en ./ x;
  146. %
  147. %              Calculate backward test, and find ncalc,
  148. %              the highest n such that the test is passed.
  149. %
  150.                test = pold*plast*eps;
  151.                test = test*(0.5-0.5/min(tempb*tempb));
  152.                p = plast * tover;
  153.                n = n - 1;
  154.                en = en - 2;
  155.                nend = min0(nb,n);
  156.                ncalc = nend + 1;
  157.                for l = nstart:nend
  158.                   pold = psavel;
  159.                   psavel = psave;
  160.                   psave = en * psavel./x + pold;
  161.                   if any(psave.*psavel > test);
  162.                      ncalc = l-1;
  163.                      test = 0;
  164.                      break
  165.                   end
  166.                end
  167.                break
  168.             end
  169.          end
  170.          if ~test
  171.             n = nend;
  172.             en = (n+n) + (alpha+alpha);
  173. %
  174. %           Calculate special significance test for nbmx > 2.
  175. %
  176.             test = max(test,max(sqrt(plast/eps).*sqrt(p+p)));
  177.          end
  178.       end
  179. %
  180. %     Calculate p-sequence until significance test passed.
  181. %
  182.       while any(p < test)
  183.          n = n + 1;
  184.          en = en + 2;
  185.          pold = plast;
  186.          plast = p;
  187.          p = en * plast./x + pold;
  188.       end
  189. %
  190. %     Initialize the backward recursion and the normalization sum.
  191. %
  192.       n = n + 1;
  193.       en = en + 2;
  194.       tempb = 0;
  195.       tempa = 1 ./ p;
  196.       em = n - 1;
  197.       empal = em + alpha;
  198.       emp2al = (em - 1) + (alpha + alpha);
  199.       sum = tempa * empal * emp2al / em;
  200.       nend = n - nb;
  201.       skip = 0;
  202.       if (nend < 0)
  203.          b(v,n) = tempa;
  204.          nend = -nend;
  205.          b(v,n+1:n+nend) = zeros(length(v),nend);
  206.       else
  207.          if (nend > 0)
  208. %
  209. %           Recur backward via difference equation, calculating 
  210. %           (but not storing) b(n), until n = nb.
  211. %
  212.             for l = 1:nend
  213.                n = n - 1;
  214.                en = en - 2;
  215.                tempc = tempb;
  216.                tempb = tempa;
  217.                tempa = (en*tempb) ./ x + tempc;
  218.                em = em - 1;
  219.                emp2al = emp2al - 1;
  220.                if (n == 1), break, end
  221.                if (n == 2), emp2al = 1; end
  222.                empal = empal - 1;
  223.                sum = (sum + tempa*empal) * emp2al / em;
  224.             end
  225.          end
  226. %
  227. %        Store b(nb).
  228. %
  229.          b(v,n) = tempa;
  230.          if (nb <= 1)
  231.             sum = (sum + sum) + tempa;
  232.             skip = -1;
  233.          else
  234. %
  235. %           Calculate and store b(nb-1).
  236. %
  237.             n = n - 1;
  238.             en = en - 2;
  239.             b(v,n)  = (en*tempa) ./ x + tempb;
  240.             if (n == 1)
  241.                skip = 1;
  242.             else
  243.                em = em - 1;
  244.                emp2al = emp2al - 1;
  245.                if (n == 2), emp2al = 1; end
  246.                empal = empal - 1;
  247.                sum = (sum + b(v,n)*empal) * emp2al / em;
  248.             end
  249.          end
  250.       end
  251.       if skip == 0
  252.          nend = n - 2;
  253.          if (nend > 0)
  254. %
  255. %        Calculate via difference equation and store b(n), until n = 2.
  256. %
  257.             for l = 1:nend
  258.                n = n - 1;
  259.                en = en - 2;
  260.                b(v,n) = (en*b(v,n+1)) ./ x + b(v,n+2);
  261.                em = em - 1;
  262.                emp2al = emp2al - 1;
  263.                if (n == 2), emp2al = 1; end
  264.                empal = empal - 1;
  265.                sum = (sum + b(v,n)*empal) * emp2al / em;
  266.             end
  267.          end
  268.       end
  269. %
  270. %     Calculate b(1)
  271. %
  272.       if skip == 0, b(v,1) = 2*empal*b(v,2) ./ x + b(v,3); end
  273.       if skip >= 0, sum = (sum + sum) + b(v,1); end
  274. %
  275. %     Normalize.  Divide all b(n) by sum.
  276. %
  277.       if (alpha ~= 0)
  278.          sum = sum * gamma(1+alpha) .* (x*0.5).^(-alpha);
  279.       end
  280.       if ~scale, sum = sum .* exp(-x); end
  281.       for n = 1:nb
  282.          b(v,n) = b(v,n) ./ sum;
  283.       end
  284.    end
  285. %
  286. %  Two-term ascending series for small x.
  287. %
  288.    v = find((0 < xx) & (xx < 1.e-4));
  289.    if any(v)
  290.       x = xx(v);
  291.       tempa = ones(length(x),1);
  292.       empal = 1 + alpha;
  293.       halfx = 0.5 * x;
  294.       if (alpha ~= 0),  tempa = halfx.^alpha /gamma(empal); end
  295.       if scale, tempa = tempa .* exp(-x); end
  296.       tempb = halfx .* halfx;
  297.       b(v,1) = tempa + tempa.*tempb / empal;
  298.       if (nb > 1)
  299. %
  300. %        Calculate higher-order functions.
  301. %
  302.          tempc = halfx;
  303.          tover = 8*realmin / x;
  304.          if (tempb ~= 0), tover = 4*realmin / tempb; end
  305.          for n = 2:nb;
  306.             tempa = tempa / empal;
  307.             empal = empal + 1;
  308.             tempa = tempa * tempc;
  309.             if (tempa <= tover*empal),  tempa = 0; end
  310.             b(v,n) = tempa + tempa*tempb / empal;
  311.          end
  312.       end
  313.    end
  314. %
  315. %  x == 0
  316. %
  317.    v = find(xx == 0);
  318.    if any(v)
  319.       % if alpha>0, I(alpha,0) = 0; I(0,0) = 1.
  320.       b(v,:) = zeros(length(v),nb);
  321.       if alpha == 0, b(v,1) = ones(length(v),1); end
  322.    end
  323. %
  324. %  Return the requested index range
  325. %
  326.    b = b(:,nfirst+1:nb);
  327.    if resize
  328.       b = reshape(b,resize(1),resize(2));
  329.    end
  330.