home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / BESSELK.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  11.0 KB  |  367 lines
  1. function bk = besselk(alpha,xx,scale)
  2. %BESSELK Modified Bessel functions of the second kind.
  3. %    K = BESSELK(ALPHA,X) computes modified Bessel functions of the
  4. %    second kind, K_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 = BESSELK(ALPHA,X,1) computes K_sub_alpha(X)*EXP(X).
  14. %
  15. %    The relationship to unmodified Bessel functions with imaginary 
  16. %    argument:
  17. %
  18. %        K_sub_alpha(x) = pi/2 * i^(-alpha) * (J_sub_alpha(i*x) +
  19. %                         Y_sub_alpha(i*x))
  20. %
  21. %    Examples:
  22. %
  23. %        besselk(3:9,[0:.2:9.8 10:.5:20],1) generates the entire 
  24. %        71-by-7 table on page 424 of Abramowitz and Stegun,
  25. %        "Handbook of Mathematical Functions."
  26. %
  27. %    See also: BESSELJ, BESSELY, BESSELI.
  28.  
  29. %   Acknowledgement:
  30. %
  31. %    This program is based on a Fortran program by W. J. Cody and
  32. %    L. Stoltz, Applied Mathematics Division, Argonne National
  33. %    Laboratory, dated May 30, 1989.  Their references include:
  34. %
  35. %        "On Temme's algorithm for the modified Bessel functions of the
  36. %        third kind," Campbell, J. B., TOMS 6(4), Dec. 1980, pp. 581-586.
  37. %
  38. %        "A Fortran IV Subroutine for the modified Bessel functions of
  39. %        the third kind of real order and real argument," Campbell, J. B.,
  40. %        Report NRC/ERB-925, National Research Council, Canada.
  41. %
  42. %    MATLAB version: C. Moler, 10/9/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.    bk = NaN*ones(length(xx),nb);
  74.    alpha = alpha(1) - nfirst;
  75. %
  76. %  ln2mgamma = log(2) - Euler's constant
  77. %  p, q - approximation for log(gamma(1+alpha))/alpha
  78. %                                         + Euler's constant
  79. %  r, s - approximation for (1-alpha*pi/sin(alpha*pi))/(2.d0*alpha)
  80. %  t    - approximation for sinh(y)/y
  81.  
  82.    ln2mgamma = 0.11593151565841244881;
  83.    p = [0.805629875690432845e00, 0.204045500205365151e02, ...
  84.         0.157705605106676174e03, 0.536671116469207504e03, ...
  85.         0.900382759291288778e03, 0.730923886650660393e03, ...
  86.         0.229299301509425145e03, 0.822467033424113231e00];
  87.    q = [0.294601986247850434e02, 0.277577868510221208e03, ...
  88.         0.120670325591027438e04, 0.276291444159791519e04, ...
  89.         0.344374050506564618e04, 0.221063190113378647e04, ...
  90.         0.572267338359892221e03];
  91.    r = [-0.48672575865218401848e+0, 0.13079485869097804016e+2, ...
  92.         -0.10196490580880537526e+3, 0.34765409106507813131e+3, ...
  93.          0.34958981245219347820e-3];
  94.    s = [-0.25579105509976461286e+2, 0.21257260432226544008e+3, ...
  95.         -0.61069018684944109624e+3, 0.42269668805777760407e+3];
  96.    t = [ 0.16125990452916363814e-9, 0.25051878502858255354e-7, ...
  97.          0.27557319615147964774e-5, 0.19841269840928373686e-3, ...
  98.          0.83333333333334751799e-2, 0.16666666666666666446e+0];
  99. %
  100. %  This algorithm can not be conveniently vectorized in x.
  101. %  So, we do it the hard way.
  102. %
  103.    for v = 1:length(xx)
  104.       x = xx(v);
  105.       enu = alpha;
  106.       ncalc = 0;
  107.       k = 0;
  108.       if (enu < sqrt(realmin)), enu = 0; end
  109.       if (enu > 0.5);
  110.             k = 1;
  111.             enu = enu - 1;
  112.       end
  113.       twonu = enu+enu;
  114.       iend = nb+k-1;
  115.       c = enu*enu;
  116.       d3 = -c;
  117.       if (x == 0)
  118.          bk(v,:) = Inf*ones(1,nb);
  119.          ncalc = nb;
  120.       elseif (x <= 1);
  121. %
  122. %        Calculation of p0 = gamma(1+alpha) * (2/x)^alpha
  123. %                       q0 = gamma(1-alpha) * (x/2)^alpha
  124. %
  125.          d1 = 0;
  126.          d2 = p(1);
  127.          t1 = 1;
  128.          t2 = q(1);
  129.          for i = 2:2:7;
  130.             d1 = c*d1+p(i);
  131.             d2 = c*d2+p(i+1);
  132.             t1 = c*t1+q(i);
  133.             t2 = c*t2+q(i+1);
  134.          end
  135.          d1 = enu*d1;
  136.          t1 = enu*t1;
  137.          f1 = log(x);
  138.          f0 = ln2mgamma+enu*(p(8)-enu*(d1+d2)/(t1+t2))-f1;
  139.          q0 = exp(-enu*(ln2mgamma-enu*(p(8)+enu*(d1-d2)/(t1-t2))-f1));
  140.          f1 = enu*f0;
  141.          p0 = exp(f1);
  142.          d1 = r(5);
  143.          t1 = 1;
  144.          for i = 1:4;
  145.             d1 = c*d1+r(i);
  146.             t1 = c*t1+s(i);
  147.          end
  148.          if (abs(f1) <= 0.5)
  149.             f1 = f1*f1;
  150.             d2 = 0;
  151.             for i = 1:6;
  152.                d2 = f1*d2+t(i);
  153.             end
  154.             d2 = f0+f0*f1*d2;
  155.          else
  156.             d2 = sinh(f1)/enu;
  157.          end
  158.          f0 = d2-enu*d1/(t1*p0);
  159.          if (x <= 1.e-10)
  160. %
  161. %        Calculation of K(alpha,x) and x*K(alpha+1,x)/K(alpha,x)
  162. %
  163.             bk(v,1) = f0+x*f0;
  164.             if (~scale), bk(v,1) = bk(v,1)-x*bk(v,1); end
  165.             ratio = p0/f0;
  166.             if (k ~= 0)
  167. %
  168. %            Calculation of K(alpha,x) and x*K(alpha+1,x)/K(alpha,x),
  169. %            alpha >= 1/2
  170. %
  171.                bk(v,1) = ratio*bk(v,1)/x;
  172.                twonu = twonu+2;
  173.                ratio = twonu;
  174.             end
  175.             if (nb == 1)
  176.                ncalc = 1;
  177.             end
  178. %
  179. %           Calculate  K(alpha+l,x)/K(alpha+l-1,x), l = 1,:nb-1
  180. %
  181.             for i = 2:nb;
  182.                bk(v,i) = ratio/x;
  183.                twonu = twonu+2;
  184.                ratio = twonu;
  185.             end
  186.             ncalc = 1;
  187.          else
  188. %
  189. %        1.0e-10 < x <= 1.0
  190. %
  191.             c = 1;
  192.             x2by4 = x*x/4;
  193.             p0 = 0.5*p0;
  194.             q0 = 0.5*q0;
  195.             d1 = -1;
  196.             d2 = 0;
  197.             bk1 = 0;
  198.             bk2 = 0;
  199.             f1 = f0;
  200.             f2 = p0;
  201.             t1 = inf;
  202.             t2 = inf;
  203.             while ((abs(t1/(f1+bk1))>eps) | (abs(t2/(f2+bk2))>eps)) ;
  204.                d1 = d1+2;
  205.                d2 = d2+1;
  206.                d3 = d1+d3;
  207.                c = x2by4*c/d2;
  208.                f0 = (d2*f0+p0+q0)/d3;
  209.                p0 = p0/(d2-enu);
  210.                q0 = q0/(d2+enu);
  211.                t1 = c*f0;
  212.                t2 = c*(p0-d2*f0);
  213.                bk1 = bk1+t1;
  214.                bk2 = bk2+t2;
  215.             end
  216.             bk1 = f1+bk1;
  217.             bk2 = 2*(f2+bk2)/x;
  218.             if scale
  219.                   d1 = exp(x);
  220.                   bk1 = bk1*d1;
  221.                   bk2 = bk2*d1;
  222.             end
  223.             wminf = 41.8341*x+7.1075;
  224.          end
  225.       elseif (eps*x > 1);
  226. %
  227. %        1/eps < x
  228. %
  229.          bk1 = sqrt(2/(pi*x));
  230.          for i = 1:nb;
  231.             bk(v,i) = bk1;
  232.          end
  233.          ncalc = nb;
  234.       else
  235. %
  236. %        1.0 < x <= 1/eps
  237. %
  238.          twox = x+x;
  239.          blpha = 0;
  240.          ratio = 0;
  241.          if (x <= 4)
  242. %
  243. %        Calculation of K(alpha+1,x)/K(alpha,x),  1.0 <= x <= 4.0
  244. %
  245.             d2 = fix(52.0583/x+5.7607);
  246.             m = d2;
  247.             d1 = d2+d2;
  248.             d2 = d2-0.5;
  249.             d2 = d2*d2;
  250.             for i = 2:m;
  251.                d1 = d1-2;
  252.                d2 = d2-d1;
  253.                ratio = (d3+d2)/(twox+d1-ratio);
  254.             end
  255. %
  256. %           Calculation of I(alpha,x) and I(alpha+1,x) by backward
  257. %           recurrence and K(alpha,x) from the Wronskian
  258. %
  259.             d2 = fix(2.7782*x+14.4303);
  260.             m = d2;
  261.             c = abs(enu);
  262.             d3 = c+c;
  263.             d1 = d3-1;
  264.             f1 = realmin;
  265.             f0 = (2*(c+d2)/x+0.5*x/(c+d2+1))*realmin;
  266.             for i = 3:m;
  267.                d2 = d2-1;
  268.                f2 = (d3+d2+d2)*f0;
  269.                blpha = (1+d1/d2)*(f2+blpha);
  270.                f2 = f2/x+f1;
  271.                f1 = f0;
  272.                f0 = f2;
  273.             end
  274.             f1 = (d3+2)*f0/x+f1;
  275.             d1 = 0;
  276.             t1 = 1;
  277.             for i = 1:7;
  278.                d1 = c*d1+p(i);
  279.                t1 = c*t1+q(i);
  280.             end
  281.             p0 = exp(c*(ln2mgamma+c*(p(8)-c*d1/t1)-log(x)))/x;
  282.             f2 = (c+0.5-ratio)*f1/x;
  283.             bk1 = p0+(d3*f0-f2+f0+blpha)/(f2+f1+f0)*p0;
  284.             if (~scale), bk1 = bk1*exp(-x); end
  285.             wminf = 6.4306*x+42.511;
  286.          else
  287. %
  288. %        Calculation of K(alpha,x) and K(alpha+1,x)/K(alpha,x)
  289. %        by backward recurrence, for x > 4.0
  290. %
  291.             dm = fix(185.3004/x+9.3715);
  292.             m = dm;
  293.             d2 = dm-0.5;
  294.             d2 = d2*d2;
  295.             d1 = dm+dm;
  296.             for i = 2:m;
  297.                dm = dm-1;
  298.                d1 = d1-2;
  299.                d2 = d2-d1;
  300.                ratio = (d3+d2)/(twox+d1-ratio);
  301.                blpha = (ratio+ratio*blpha)/dm;
  302.             end
  303.             d = sqrt(2/pi);
  304.             bk1 = 1/((d+d*blpha)*sqrt(x));
  305.             if (~scale), bk1 = bk1*exp(-x); end
  306.             wminf = 1.35633*(x-abs(x-20))+84.5096;
  307.          end
  308. %
  309. %        Calculation of K(alpha+1,x) from K(alpha,x) and
  310. %        K(alpha+1,x)/K(alpha,x)
  311. %
  312.          bk2 = bk1+bk1*(enu+0.5-ratio)/x;
  313.       end
  314.       if ~ncalc
  315. %
  316. %        Calculation of ncalc, K(alpha+i,x), i = 0:ncalc-1,
  317. %        K(alpha+i,x)/K(alpha+i-1,x), i = ncalc:nb-1
  318. %
  319.          ncalc = nb;
  320.          bk(v,1) = bk1;
  321.          if (iend > 0)
  322.             j = 2-k;
  323.             if (j > 0), bk(v,j) = bk2; end
  324.             if (iend > 1)
  325.                m = iend;
  326.                if wminf-enu < m, m = fix(wminf-enu); end
  327.                for i = 2:m;
  328.                   t1 = bk1;
  329.                   bk1 = bk2;
  330.                   twonu = twonu+2;
  331.                   bk2 = twonu/x*bk1+t1;
  332.                   itemp = i;
  333.                   j = j+1;
  334.                   if (j > 0), bk(v,j) = bk2; end
  335.                end
  336.                m = itemp;
  337.                if (m ~= iend)
  338.                   ratio = bk2/bk1;
  339.                   mplus1 = m+1;
  340.                   for i = mplus1:iend;
  341.                      twonu = twonu+2;
  342.                      ratio = twonu/x+1/ratio;
  343.                      j = j+1;
  344.                      if (j > 1)
  345.                            bk(v,j) = ratio;
  346.                         else
  347.                            bk2 = ratio*bk2;
  348.                      end
  349.                   end
  350.                   ncalc = max(mplus1-k,1);
  351.                   if (ncalc == 1), bk(v,1) = bk2; end
  352.                end
  353.             end
  354.          end
  355.       end
  356.       for i = ncalc+1:nb
  357.          bk(v,i) = bk(v,i-1)*bk(v,i);
  358.       end   
  359.    end
  360. %
  361. %  Return the requested index range
  362. %
  363.    bk = bk(:,nfirst+1:nb);
  364.    if resize
  365.       bk = reshape(bk,resize(1),resize(2));
  366.    end
  367.