function value = bessel_i0 ( arg ) %*****************************************************************************80 % %% BESSEL_I0 evaluates the modified Bessel function of the first kind and order zero. % % Discussion: % % The main computation evaluates slightly modified forms of % minimax approximations generated by Blair and Edwards, Chalk % River (Atomic Energy of Canada Limited) Report AECL-4928, % October, 1974. This transportable program is patterned after % the machine dependent FUNPACK packet NATSI0, but cannot match % that version for efficiency or accuracy. This version uses % rational functions that theoretically approximate I-SUB-0(X) % to at least 18 significant decimal digits. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 28 October 2004 % % Author: % % Original FORTRAN77 version by William Cody, Laura Stoltz. % MATLAB version by John Burkardt. % % Parameters: % % Input, real ARG, the argument. % % Output, real VALUE, the value of the modified Bessel function % of the first kind. % exp40 = 2.353852668370199854E+17; p = [ ... -5.2487866627945699800E-18, ... -1.5982226675653184646E-14, ... -2.6843448573468483278E-11, ... -3.0517226450451067446E-08, ... -2.5172644670688975051E-05, ... -1.5453977791786851041E-02, ... -7.0935347449210549190E+00, ... -2.4125195876041896775E+03, ... -5.9545626019847898221E+05, ... -1.0313066708737980747E+08, ... -1.1912746104985237192E+10, ... -8.4925101247114157499E+11, ... -3.2940087627407749166E+13, ... -5.5050369673018427753E+14, ... -2.2335582639474375249E+15 ]; pp = [ ... -3.9843750000000000000E-01, ... 2.9205384596336793945E+00, ... -2.4708469169133954315E+00, ... 4.7914889422856814203E-01, ... -3.7384991926068969150E-03, ... -2.6801520353328635310E-03, ... 9.9168777670983678974E-05, ... -2.1877128189032726730E-06 ]; q = [ ... -3.7277560179962773046E+03, ... 6.5158506418655165707E+06, ... -6.5626560740833869295E+09, ... 3.7604188704092954661E+12, ... -9.7087946179594019126E+14 ]; qq = [ ... -3.1446690275135491500E+01, ... 8.5539563258012929600E+01, ... -6.0228002066743340583E+01, ... 1.3982595353892851542E+01, ... -1.1151759188741312645E+00, ... 3.2547697594819615062E-02, ... -5.5194330231005480228E-04 ]; rec15 = 6.6666666666666666666E-02; xmax = 91.9E+00; xsmall = 2.98E-08; x = abs ( arg ); if ( x < xsmall ) value = 1.0; elseif ( x < 15.0 ) % % XSMALL <= ABS(ARG) < 15.0 % xx = x * x; sump = p(1); for i = 2 : 15 sump = sump * xx + p(i); end xx = xx - 225.0; sumq = (((( ... xx + q(1) ) ... * xx + q(2) ) ... * xx + q(3) ) ... * xx + q(4) ) ... * xx + q(5); value = sump / sumq; elseif ( 15.0 <= x ) if ( xmax < x ) value = r8_huge ( ); else % % 15.0 <= ABS(ARG) % xx = 1.0 / x - rec15; sump = (((((( ... pp(1) ... * xx + pp(2) ) ... * xx + pp(3) ) ... * xx + pp(4) ) ... * xx + pp(5) ) ... * xx + pp(6) ) ... * xx + pp(7) ) ... * xx + pp(8); sumq = (((((( ... xx + qq(1) ) ... * xx + qq(2) ) ... * xx + qq(3) ) ... * xx + qq(4) ) ... * xx + qq(5) ) ... * xx + qq(6) ) ... * xx + qq(7); value = sump / sumq; % % Calculation reformulated to avoid premature overflow. % if ( x <= xmax - 15.0 ) a = exp ( x ); b = 1.0; else a = exp ( x - 40.0 ); b = exp40; end value = ( ( value * a - pp(1) * a ) / sqrt ( x ) ) * b; end end return end