function [ ci, si ] = cisi ( x )

%*****************************************************************************80
%
%% cisi() computes cosine integral Ci(x) and sine integral Si(x).
%
%  Licensing:
%
%    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However,
%    they give permission to incorporate this routine into a user program
%    provided that the copyright is acknowledged.
%
%  Modified:
%
%    27 March 2025
%
%  Author:
%
%    Original FORTRAN77 version by Shanjie Zhang, Jianming Jin.
%    This version by John Burkardt.
%
%  Reference:
%
%    Shanjie Zhang, Jianming Jin,
%    Computation of Special Functions,
%    Wiley, 1996,
%    ISBN: 0-471-11963-6,
%    LC: QA351.C45.
%
%  Input:
%
%    real x: the argument of Ci(x) and Si(x).
%
%  Output:
%
%    real ci, si: the values of Ci(x) and Si(x).
%
  p2 = 1.570796326794897E+00;
  el = 0.5772156649015329E+00;
  epsilon = 1.0E-15;
  x2 = x * x;

  if ( x == 0.0E+00 )

    ci = - Inf;
    si = 0.0E+00;

  elseif ( x <= 16.0E+00 )

    xr = -0.25E+00 * x2;
    ci = el + log ( x ) + xr;
    for k = 2 : 40
      xr = -0.5E+00 * xr * ( k - 1 ) / ( k * k * ( 2 * k - 1 ) ) * x2;
      ci = ci + xr;
      if ( abs ( xr ) < abs ( ci ) * epsilon )
        break;
      end
    end

    xr = x;
    si = x;
    for k = 1 : 40
      xr = -0.5E+00 * xr * ( 2 * k - 1 ) / k / ( 4 * k * k + 4 * k + 1 ) * x2;
      si = si + xr;
      if ( abs ( xr ) < abs ( si ) * epsilon )
        return
      end
    end

  elseif ( x <= 32.0E+00 )

    m = floor ( 47.2E+00 + 0.82E+00 * x );
    bj = zeros ( m, 1 );
    xa1 = 0.0E+00;
    xa0 = 1.0E-100;
    for k = m : -1 : 1
      xa = 4.0E+00 * k * xa0 / x - xa1;
      bj(k) = xa;
      xa1 = xa0;
      xa0 = xa;
    end
    xs = bj(1);
    for k = 3 : 2 : m
      xs = xs + 2.0E+00 * bj(k);
    end
    bj(1) = bj(1) / xs;
    for k = 2 : m
      bj(k) = bj(k) / xs;
    end
    xr = 1.0E+00;
    xg1 = bj(1);
    for k = 2 : m
      xr = 0.25E+00 * xr * ( 2.0E+00 * k - 3.0E+00 )^2 ...
        / ( ( k - 1.0E+00 ) * ( 2.0E+00 * k - 1.0E+00 )^2 ) * x;
      xg1 = xg1 + bj(k) * xr;
    end

    xr = 1.0E+00;
    xg2 = bj(1);
    for k = 2 : m
      xr = 0.25E+00 * xr * ( 2.0E+00 * k - 5.0E+00 )^2 ...
        / ( ( k-1.0E+00 ) * ( 2.0E+00 * k - 3.0E+00 )^2 ) * x;
      xg2 = xg2 + bj(k) * xr;
    end

    xcs = cos ( x / 2.0E+00 );
    xss = sin ( x / 2.0E+00 );
    ci = el + log ( x ) - x * xss * xg1 + 2.0 * xcs * xg2 - 2.0 * xcs * xcs;
    si = x * xcs * xg1 + 2.0 * xss * xg2 - sin ( x );

  else

    xr = 1.0E+00;
    xf = 1.0E+00;
    for k = 1 : 9
      xr = -2.0E+00 * xr * k * ( 2 * k - 1 ) / x2;
      xf = xf + xr;
    end
    xr = 1.0E+00 / x;
    xg = xr;
    for k = 1 : 8
      xr = -2.0E+00 * xr * ( 2 * k + 1 ) * k / x2;
      xg = xg + xr;
    end
    ci = xf * sin ( x ) / x - xg * cos ( x ) / x;
    si = p2 - xf * cos ( x ) / x - xg * sin ( x ) / x;

  end

  return
end
