home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / BESSELY.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  8.9 KB  |  318 lines
  1. function y = bessely(alpha,xx)
  2. %BESSELY Bessel functions of the second kind.
  3. %    Y = BESSELY(ALPHA,X) calculates Bessel functions of the second kind,
  4. %    Y_sub_alpha(X) for real, non-negative order ALPHA and argument X.
  5. %    In general, both ALPHA and X may be vectors.  The output Y is
  6. %    an m-y-n matrix with m = lenth(X), n = length(ALPHA) and
  7. %        Y(i,k) = Y_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. %        bessely(3:9,(10:.2:20)') generates the 51-by-7 table on page 401
  16. %        of Abramowitz and Stegun, "Handbook of Mathematical Functions."
  17. %
  18. %    See also: BESSELJ, BESSELI, BESSELK.
  19.  
  20. %   Acknowledgement:
  21. %
  22. %    This program is based on a Fortran program by W. J. Cody,
  23. %    Applied Mathematics Division, Argonne National Laboratory, 
  24. %    dated March 19, 1990.  Cody references earlier work by
  25. %    N. M. Temme and J. B. Campbell.
  26. %
  27. %    References: "Bessel functions J_nu(x) and Y_nu(x) of real
  28. %       order and real argument," Campbell, J. B.,
  29. %       Comp. Phy. Comm. 18, 1979, pp. 133-142.
  30. %
  31. %       "On the numerical evaluation of the ordinary
  32. %       Bessel function of the second kind," Temme,
  33. %       N. M., J. Comput. Phys. 21, 1976, pp. 343-350.
  34. %
  35. %    MATLAB version: C. Moler, 10/2/92, 12/11/92.
  36. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  37.  
  38. %
  39. %  Check for real, non-negative arguments.
  40. %
  41.    if any(imag(xx)) | any(xx < 0) | any(imag(alpha)) | any(alpha < 0)
  42.       error('Input arguments must be real and nonnegative.')
  43.    end
  44.    if isempty(alpha) | isempty(xx)
  45.       bk = []; 
  46.       return
  47.    end
  48. %
  49. %  Break alpha into integer and fractional parts,
  50. %  and initialize result array.
  51. %
  52.    nfirst = fix(alpha(1));
  53.    nb = fix(alpha(length(alpha))) + 1;
  54.    if ~(nb <= 1001)
  55.       error('Alpha must be <= 1000.')
  56.    end
  57.    if length(alpha) > 1
  58.       if any(abs(diff(alpha)-1) > 4*nb*eps)
  59.          error('Increment in alpha must be 1.')
  60.       end
  61.    end
  62.    resize = (length(alpha) == 1);
  63.    if resize, resize = size(xx); end
  64.    xx = xx(:);
  65.    zz = (xx==0);  % Remember location of x == 0.
  66.    if any(zz), xx = xx + zz; end
  67.    y = NaN*ones(length(xx),nb);
  68.    y0 = NaN*ones(length(xx),1);
  69.    y1 = NaN*ones(length(xx),1);
  70.    alpha = alpha(1) - nfirst;
  71. %
  72.    enu = alpha;
  73.    na = round(enu);
  74.    if (na == 1), enu = enu - na; end
  75.    if (enu == -0.5)
  76.       p = sqrt(2./(pi*xx));
  77.       y0 = p .* sin(xx);
  78.       y1 = -p .* cos(xx);
  79.    else
  80. %
  81. %     Use Temme's scheme for small x.
  82. %
  83.       v = find(xx < 3);
  84.       if any(v)
  85.          x = xx(v);
  86.          b = x * 0.5;
  87.          d = -log(b);
  88.          f = enu * d;
  89.          e = b.^(-enu);
  90.          if (abs(enu) < 1.e-8)
  91.             c = 1/pi;
  92.          else
  93.             c = enu / sin(enu*pi);
  94.          end
  95. %
  96. %        Computation of sinh(f)/f
  97. %
  98.          if any(abs(f) < 1)
  99.             x2 = f.*f;
  100.             en = 19;
  101.             s = 1;
  102.             for i = 1:9
  103.                s = s.*x2/en/(en-1)+1;
  104.                en = en - 2;
  105.             end
  106.          else
  107.             s = (e - 1./e) * 0.5 ./ f;
  108.          end
  109. %
  110. %        Computation of 1/gamma(1-a) using Chebyshev polynomials.
  111. %
  112.          ch = [-0.67735241822398840964d-23,-0.61455180116049879894d-22, ...
  113.                 0.29017595056104745456d-20, 0.13639417919073099464d-18, ...
  114.                 0.23826220476859635824d-17,-0.90642907957550702534d-17, ...
  115.                -0.14943667065169001769d-14,-0.33919078305362211264d-13, ...
  116.                -0.17023776642512729175d-12, 0.91609750938768647911d-11, ...
  117.                 0.24230957900482704055d-09, 0.17451364971382984243d-08, ...
  118.                -0.33126119768180852711d-07,-0.86592079961391259661d-06, ...
  119.                -0.49717367041957398581d-05, 0.76309597585908126618d-04, ...
  120.                 0.12719271366545622927d-02, 0.17063050710955562222d-02, ...
  121.                -0.76852840844786673690d-01,-0.28387654227602353814d+00, ...
  122.                 0.92187029365045265648d+00];
  123.  
  124.          x2 = enu*enu*8;
  125.          aye = ch(1);
  126.          even = 0;
  127.          alfa = ch(2);
  128.          odd = 0;
  129.          for i = 3:2:19
  130.             even = -(aye+aye+even);
  131.             aye = -even*x2 - aye + ch(i);
  132.             odd = -(alfa+alfa+odd);
  133.             alfa = -odd*x2 - alfa + ch(i+1);
  134.          end
  135.          even = (even*0.5+aye)*x2 - aye + ch(21);
  136.          odd = (odd+alfa)*2;
  137.          gamma = odd*enu + even;
  138. %
  139. %        End of computation of 1/gamma(1-a)
  140. %
  141.          g = e * gamma;
  142.          e = (e + 1./e) * 0.5;
  143.          f = 2*c*(odd.*e+even.*s.*d);
  144.          e = enu*enu;
  145.          p = g.*c;
  146.          q = 1./(pi*g);
  147.          c = enu*pi/2;
  148.          if (abs(c) < 1.e-8)
  149.             r = 1;
  150.          else
  151.             r = sin(c)/c;
  152.          end
  153.          r = pi*c*r*r;
  154.          c = 1;
  155.          d = -b.*b;
  156.          h = zeros(size(x));
  157.          y0(v) = f + r*q;
  158.          y1(v) = p;
  159.          en = 1;
  160.          while (abs(g./(1+abs(y0(v)))) + abs(h./(1+abs(y1(v)))) > eps)
  161.             f = (f*en+p+q)/(en*en-e);
  162.             c = c.*d/en;
  163.             p = p/(en-enu);
  164.             q = q/(en+enu);
  165.             g = c.*(f+r.*q);
  166.             h = c.*p - en*g;
  167.             y0(v) = y0(v) + g;
  168.             y1(v) = y1(v) + h;
  169.             en = en + 1;
  170.          end
  171.          y0(v) = -y0(v);
  172.          y1(v) = -y1(v)./b;
  173.       end
  174. %
  175. %     Use Temme's scheme for moderate x
  176. %
  177.       v = find((xx >= 3) & (xx < 16));
  178.       if any(v)
  179.          x = xx(v);
  180.          c = (0.5-enu)*(0.5+enu);
  181.          b = x + x;
  182.          e = (x/pi*cos(enu*pi)/eps);
  183.          e = e.*e;
  184.          p = ones(size(x));
  185.          q = -x;
  186.          r = 1 + x.*x;
  187.          s = r;
  188.          en = 2;
  189.          while any(r*en*en < e)
  190.             en1 = en+1;
  191.             d = (en-1+c/en)./s;
  192.             p = (en+en-p.*d)/en1;
  193.             q = (-b+q.*d)/en1;
  194.             s = p.*p + q.*q;
  195.             r = r.*s;
  196.             en = en1;
  197.          end
  198.          f = p./s;
  199.          p = f;
  200.          g = -q./s;
  201.          q = g;
  202.          en = en - 1;
  203.          while any(en > 0);
  204.             r = en1*(2-p)-2;
  205.             s = b + en1*q;
  206.             d = (en-1+c/en)./(r.*r+s.*s);
  207.             p = d.*r;
  208.             q = d.*s;
  209.             e = f + 1;
  210.             f = p.*e - g.*q;
  211.             g = q.*e + p.*g;
  212.             en1 = en;
  213.             en = en - 1;
  214.          end
  215.          f = 1 + f;
  216.          d = f.*f + g.*g;
  217.          pa = f./d;
  218.          qa = -g./d;
  219.          d = enu + 0.5 -p;
  220.          q = q + x;
  221.          pa1 = (pa.*q-qa.*d)./x;
  222.          qa1 = (qa.*q+pa.*d)./x;
  223.          b = x - (pi/2)*(enu+0.5);
  224.          c = cos(b);
  225.          s = sin(b);
  226.          d = sqrt(2./(pi*x));
  227.          y0(v) = d.*(pa.*s+qa.*c);
  228.          y1(v) = d.*(qa1.*s-pa1.*c);
  229.       end
  230. %
  231. %     Use Campbell's asymptotic scheme for large x.
  232. %
  233.       v = find(xx >= 16);
  234.       if any(v)
  235.          x = xx(v);
  236.          d1 = fix(x/(5*pi));
  237.          dmu = ((x-15*d1)-d1*(5*pi-15))-(alpha+0.5)*pi/2;
  238.          sig = 1-2*rem(d1,2);
  239.          cosmu = sig.*cos(dmu);
  240.          sinmu = sig.*sin(dmu);
  241.          ddiv = 8 * x;
  242.          dmu = alpha;
  243.          den = sqrt(x);
  244.          for k = 1:2
  245.             p = cosmu;
  246.             cosmu = sinmu;
  247.             sinmu = -p;
  248.             d1 = (2*dmu-1)*(2*dmu+1);
  249.             d2 = 0;
  250.             div = ddiv;
  251.             p = 0;
  252.             q = 0;
  253.             q0 = d1./div;
  254.             term = q0;
  255.             for i = 2:20
  256.                d2 = d2 + 8;
  257.                d1 = d1 - d2;
  258.                div = div + ddiv;
  259.                term = -term*d1./div;
  260.                p = p + term;
  261.                d2 = d2 + 8;
  262.                d1 = d1 - d2;
  263.                div = div + ddiv;
  264.                term = term*d1./div;
  265.                q = q + term;
  266.                if (abs(term) <= eps), break, end
  267.             end
  268.             p = p + 1;
  269.             q = q + q0;
  270.             if (k == 1)
  271.                y0(v) = sqrt(2/pi) * (p.*cosmu-q.*sinmu) ./ den;
  272.             else
  273.                y1(v) = sqrt(2/pi) * (p.*cosmu-q.*sinmu) ./ den;
  274.             end
  275.             dmu = dmu + 1;
  276.          end
  277.       end
  278.    end
  279.    if (na == 1)
  280.       v = find(xx < 16);
  281.       if any(v)
  282.          h = 2*(enu+1)./xx(v);
  283.          h = h.*y1(v) - y0(v);
  284.          y0(v) = y1(v);
  285.          y1(v) = h;
  286.       end
  287.    end
  288. %
  289. %  Now have first two Y's.
  290. %  Use three-term forward recurrence for the rest.
  291. %
  292.    y(:,1) = y0;
  293.    if nb > 1
  294.       y(:,2) = y1;
  295.       if nb > 2
  296.          aye = 1 + alpha;
  297.          twobyx = 2./xx;
  298.          for i = 3:nb
  299.             y(:,i) = twobyx.*aye.*y(:,i-1) - y(:,i-2);
  300.             aye = aye + 1;
  301.          end
  302.       end
  303.    end
  304. %
  305. %  Y(alpha,0) = -Inf.
  306. %
  307.    v = find(zz);
  308.    if any(v)
  309.       y(v,:) = -inf*ones(length(v),nb);
  310.    end
  311. %
  312. %  Return the requested index range
  313. %
  314.    y = y(:,nfirst+1:nb);
  315.    if resize
  316.       y = reshape(y,resize(1),resize(2));
  317.    end
  318.