octave_packages/specfun-1.1.0/expint_Ei.m

   1 ## Copyright (C) 2006   Sylvain Pelissier   <sylvain.pelissier@gmail.com>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{y} =} expint_Ei (@var{x})
  18 ## Compute the exponential integral,
  19 ## @verbatim
  20 ##                      infinity
  21 ##                     /
  22 ##    expint_Ei(x) = - | exp(t)/t dt
  23 ##                     /
  24 ##                     -x
  25 ## @end verbatim
  26 ## @seealso{expint, expint_E1}
  27 ## @end deftypefn
  28
  29 function y = expint_Ei(x)
  30   if (nargin != 1)
  31     print_usage;
  32   endif
  33   y = zeros(size(x));
  34   F = @(x) exp(-x)./x;
  35   s = prod(size(x));
  36   for t = 1:s;
  37     if(x(t)<0 && imag(x(t)) == 0)
  38       y(t) = -quad(F,-x(t),Inf);
  39     else
  40       if(abs(x(t)) > 2 && imag(x(t)) == 0)
  41         y(t) = expint_Ei(2) - quad(F,-x(t),-2);
  42       else
  43         if(abs(x(t)) >= 10)
  44           if(imag(x(t)) <= 0)
  45             a1 = 4.03640;
  46             a2 = 1.15198;
  47             b1 = 5.03637;
  48             b2 = 4.19160;
  49             y(t) = -(x(t).^2 - a1.*x(t) + a2)./((x(t).^2-b1.*x(t)+b2).*(-x(t)).*exp(-x(t)))-i.*pi;
  50           else
  51             y(t) = conj(expint_Ei(conj(x(t))));
  52           endif;
  53         ## Serie Expansion
  54         else
  55           for k = 1:100;
  56             y(t) = y(t) + x(t).^k./(k.*factorial(k));
  57           endfor
  58           y(t) = 0.577215664901532860606512090082402431 + log(x(t)) + y(t);
  59         endif
  60       endif
  61     endif
  62   endfor
  63 endfunction