## Copyright (C) 2017 Michele Ginesi
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## Author: Michele Ginesi <address@hidden>

## -*- texinfo -*-
## @deftypefn {} {} sinint (@var{x})
## Compute the sine integral:
## @tex
## $$
## {\rm S_i} (x) = \int_0^z {\sin (t) \over t} dt
## $$
## @end tex
## @ifnottex
##
## @example
## @group
##            x
##           /
## S_i (x) = | sin (t) / t dt
##           /
##          0
## @end group
## @end example
##
## @end deftypefn

function [y] = sinint (x)

  if (nargin != 1)
    print_usage ();
  endif

  xx = x;
  ii_neg = (real (xx) < 0);
  xx(ii_neg) *= -1;
  y = -0.5i * (expint (1i * xx) - expint (-1i * xx)) + pi / 2;
  y(ii_neg) *= -1;

  # Trivial values
  y(x == 0) = 0;
  y(x == Inf) = pi / 2;
  y(x == - Inf) = - pi / 2;

endfunction

%!test
%! x = 1.1;
%! %y = sym(11)/10;
%! A = sinint (x);
%! %B = double (sinint (y));
%! B = 1.02868521867373;
%! assert (A, B, -5e-15);

%!test
%! %y = [2 3 sym(pi); exp(sym(1)) 5 6];
%! x = [2, 3, pi; exp(1), 5, 6];
%! A = sinint (x);
%! %B = double (sinint (y));
%! B = [1.60541297680269, 1.84865252799947, 1.85193705198247e+00; ...
%!      1.82104026914757, 1.54993124494467, 1.42468755128051e+00];
%! assert (A, B, -5e-15);

%!assert (sinint (0), 0)
%!assert (sinint (inf), pi/2)
%!assert (sinint (-inf), -pi/2)

%%tests against maple
%!assert (sinint (1), 0.9460830703671830149414, -2*eps)
%!assert (sinint (-1), -0.9460830703671830149414, -2*eps)
%!assert (sinint (pi), 1.851937051982466170361, -2*eps)
%!assert (sinint (-pi), -1.851937051982466170361, -2*eps)
%!assert (sinint (300), 1.5708810882137495193, -2*eps)
%!assert (sinint (1e4), 1.5708915453859619157, -2*eps)
%!assert (sinint (20i), 1.2807826332028294459e7*1i, -2*eps)

%!test
%! % maple:
%! % > A := [1+2*I, -2 + 5*I, 100, 10*I, -1e-4 + 1e-6*I, -20 + I];
%! % > for a in A do evalf(Si(a)) end do;
%! x = [1+2i; -2+5i; 100; 10i; -1e-4 + 1e-6*1i; -20-1i];
%! A = [ 1.6782404878293681180 + 2.0396845546022061045*1i
%!      -18.154174221650281533 + 1.6146414539230479060*1i
%!       1.5622254668890562934
%!       1246.1144901994233444*1i
%!      -0.000099999999944461111128 + 0.99999999833338888972e-6*1i
%!      -1.5386156269726011209 - 0.053969388020443786229*1i ];
%! B = sinint (x);
%! assert (A, B, -6e-16)