function value = r4_e1 ( x ) %*****************************************************************************80 % %% R4_E1 evaluates the exponential integral E1 for an R4 argument. % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 25 September 2011 % % Author: % % Original FORTRAN77 version by Wayne Fullerton. % MATLAB version by John Burkardt. % % Reference: % % Wayne Fullerton, % Portable Special Function Routines, % in Portability of Numerical Software, % edited by Wayne Cowell, % Lecture Notes in Computer Science, Volume 57, % Springer 1977, % ISBN: 978-3-540-08446-4, % LC: QA297.W65. % % Parameters: % % Input, real X, the argument. % % Output, real VALUE, the exponential integral E1 evaluated at X. % persistent ae11cs persistent ae12cs persistent ae13cs persistent ae14cs persistent e11cs persistent e12cs persistent ntae11 persistent ntae12 persistent ntae13 persistent ntae14 persistent nte11 persistent nte12 persistent xmax if ( isempty ( ntae11 ) ) ae11cs = [ ... 0.12150323971606579E+00, ... -0.065088778513550150E+00, ... 0.004897651357459670E+00, ... -0.000649237843027216E+00, ... 0.000093840434587471E+00, ... 0.000000420236380882E+00, ... -0.000008113374735904E+00, ... 0.000002804247688663E+00, ... 0.000000056487164441E+00, ... -0.000000344809174450E+00, ... 0.000000058209273578E+00, ... 0.000000038711426349E+00, ... -0.000000012453235014E+00, ... -0.000000005118504888E+00, ... 0.000000002148771527E+00, ... 0.000000000868459898E+00, ... -0.000000000343650105E+00, ... -0.000000000179796603E+00, ... 0.000000000047442060E+00, ... 0.000000000040423282E+00, ... -0.000000000003543928E+00, ... -0.000000000008853444E+00, ... -0.000000000000960151E+00, ... 0.000000000001692921E+00, ... 0.000000000000607990E+00, ... -0.000000000000224338E+00, ... -0.000000000000200327E+00, ... -0.000000000000006246E+00, ... 0.000000000000045571E+00, ... 0.000000000000016383E+00, ... -0.000000000000005561E+00, ... -0.000000000000006074E+00, ... -0.000000000000000862E+00, ... 0.000000000000001223E+00, ... 0.000000000000000716E+00, ... -0.000000000000000024E+00, ... -0.000000000000000201E+00, ... -0.000000000000000082E+00, ... 0.000000000000000017E+00 ]'; ae12cs = [ ... 0.58241749513472674E+00, ... -0.15834885090578275E+00, ... -0.006764275590323141E+00, ... 0.005125843950185725E+00, ... 0.000435232492169391E+00, ... -0.000143613366305483E+00, ... -0.000041801320556301E+00, ... -0.000002713395758640E+00, ... 0.000001151381913647E+00, ... 0.000000420650022012E+00, ... 0.000000066581901391E+00, ... 0.000000000662143777E+00, ... -0.000000002844104870E+00, ... -0.000000000940724197E+00, ... -0.000000000177476602E+00, ... -0.000000000015830222E+00, ... 0.000000000002905732E+00, ... 0.000000000001769356E+00, ... 0.000000000000492735E+00, ... 0.000000000000093709E+00, ... 0.000000000000010707E+00, ... -0.000000000000000537E+00, ... -0.000000000000000716E+00, ... -0.000000000000000244E+00, ... -0.000000000000000058E+00 ]'; ae13cs = [ ... -0.60577324664060346E+00, ... -0.11253524348366090E+00, ... 0.013432266247902779E+00, ... -0.001926845187381145E+00, ... 0.000309118337720603E+00, ... -0.000053564132129618E+00, ... 0.000009827812880247E+00, ... -0.000001885368984916E+00, ... 0.000000374943193568E+00, ... -0.000000076823455870E+00, ... 0.000000016143270567E+00, ... -0.000000003466802211E+00, ... 0.000000000758754209E+00, ... -0.000000000168864333E+00, ... 0.000000000038145706E+00, ... -0.000000000008733026E+00, ... 0.000000000002023672E+00, ... -0.000000000000474132E+00, ... 0.000000000000112211E+00, ... -0.000000000000026804E+00, ... 0.000000000000006457E+00, ... -0.000000000000001568E+00, ... 0.000000000000000383E+00, ... -0.000000000000000094E+00, ... 0.000000000000000023E+00 ]'; ae14cs = [ ... -0.1892918000753017E+00, ... -0.08648117855259871E+00, ... 0.00722410154374659E+00, ... -0.00080975594575573E+00, ... 0.00010999134432661E+00, ... -0.00001717332998937E+00, ... 0.00000298562751447E+00, ... -0.00000056596491457E+00, ... 0.00000011526808397E+00, ... -0.00000002495030440E+00, ... 0.00000000569232420E+00, ... -0.00000000135995766E+00, ... 0.00000000033846628E+00, ... -0.00000000008737853E+00, ... 0.00000000002331588E+00, ... -0.00000000000641148E+00, ... 0.00000000000181224E+00, ... -0.00000000000052538E+00, ... 0.00000000000015592E+00, ... -0.00000000000004729E+00, ... 0.00000000000001463E+00, ... -0.00000000000000461E+00, ... 0.00000000000000148E+00, ... -0.00000000000000048E+00, ... 0.00000000000000016E+00, ... -0.00000000000000005E+00 ]'; e11cs = [ ... -16.113461655571494026E+00, ... 7.7940727787426802769E+00, ... -1.9554058188631419507E+00, ... 0.37337293866277945612E+00, ... -0.05692503191092901938E+00, ... 0.00721107776966009185E+00, ... -0.00078104901449841593E+00, ... 0.00007388093356262168E+00, ... -0.00000620286187580820E+00, ... 0.00000046816002303176E+00, ... -0.00000003209288853329E+00, ... 0.00000000201519974874E+00, ... -0.00000000011673686816E+00, ... 0.00000000000627627066E+00, ... -0.00000000000031481541E+00, ... 0.00000000000001479904E+00, ... -0.00000000000000065457E+00, ... 0.00000000000000002733E+00, ... -0.00000000000000000108E+00 ]'; e12cs = [ ... -0.037390214792202795E+00, ... 0.042723986062209577E+00, ... -0.1303182079849700544E+00, ... 0.01441912402469889073E+00, ... -0.00134617078051068022E+00, ... 0.00010731029253063780E+00, ... -0.00000742999951611943E+00, ... 0.00000045377325690753E+00, ... -0.00000002476417211390E+00, ... 0.00000000122076581374E+00, ... -0.00000000005485141480E+00, ... 0.00000000000226362142E+00, ... -0.00000000000008635897E+00, ... 0.00000000000000306291E+00, ... -0.00000000000000010148E+00, ... 0.00000000000000000315E+00 ]'; eta = 0.1 * r4_mach ( 3 ); ntae11 = r4_inits ( ae11cs, 39, eta ); ntae12 = r4_inits ( ae12cs, 25, eta ); nte11 = r4_inits ( e11cs, 19, eta ); nte12 = r4_inits ( e12cs, 16, eta ); ntae13 = r4_inits ( ae13cs, 25, eta ); ntae14 = r4_inits ( ae14cs, 26, eta ); xmax = - log ( r4_mach ( 1 ) ); xmax = xmax - log ( xmax ); end if ( x <= - 10.0 ) value = exp ( - x ) / x * ( 1.0 ... + r4_csevl ( 20.0 / x + 1.0, ae11cs, ntae11 ) ); elseif ( x <= - 4.0 ) value = exp ( - x ) / x * ( 1.0 + r4_csevl ( ... ( 40.0 / x + 7.0 ) / 3.0, ae12cs, ntae12 ) ); elseif ( x <= - 1.0 ) value = - log ( abs ( x ) ) + r4_csevl ( ... ( 2.0 * x + 5.0 ) / 3.0, e11cs, nte11 ); elseif ( x == 0.0 ) fprintf ( 1, '\n' ); fprintf ( 1, 'R4_E1 - Fatal error!\n' ); fprintf ( 1, ' X is zero.\n' ); error ( 'R4_E1 - Fatal error!' ); elseif ( x <= 1.0 ) value = ( - log ( abs ( x ) ) - 0.6875 + x ) ... + r4_csevl ( x, e12cs, nte12 ); elseif ( x <= 4.0 ) value = exp ( - x ) / x * ( 1.0 + r4_csevl ( ... ( 8.0 / x - 5.0 ) / 3.0, ae13cs, ntae13 ) ); elseif ( x <= xmax ) value = exp ( - x ) / x * ( 1.0 + r4_csevl ( ... 8.0 / x - 1.0, ae14cs, ntae14 ) ); else value = 0.0; end return end