#! /usr/bin/env python3
#
def sine_integral_test ( ):

#*****************************************************************************80
#
## sine_integral_test() tests sine_integral().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    27 March 2026
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'sine_integral_test():' )
  print ( '  python version: ' + platform.python_version ( ) )
  print ( '  numpy version:  ' + np.version.version )
  print ( '  Test sine_integral().' )

  ci_test ( );
  si_test ( );
#
#  Terminate.
#
  print ( '' )
  print ( 'sine_integral_test():' )
  print ( '  Normal end of execution.' )

  return

def ci_test ( ):

#*****************************************************************************80
#
## ci_test() tests ci().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    27 March 2025
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'ci_test():' )
  print ( '  ci() evaluates the Cosine Integral function CI(X).' )
  print ( '' )
  print ( '          X           CI(X)' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = ci_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2, _ = cisi ( x )

    print ( '  %14.6g  %24.16g  %24.16g' % ( x, fx1, fx2 ) )

  return

def si_test ( ):

#*****************************************************************************80
#
## si_test() tests si().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    27 March 2025
#
#  Author:
#
#    John Burkardt
#
  print ( '' )
  print ( 'si_test():' )
  print ( '  si() evalues the sine integral function.' )
  print ( '' )
  print ( '      X            SI(X)' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = si_values ( n_data )

    if ( n_data == 0 ):
      break

    _, fx2 = cisi ( x )

    print ( '  %14.6g  %24.16g  %24.16g' % ( x, fx1, fx2 ) )

  return

def cisi ( x ):

#*****************************************************************************80
#
## cisi() computes cosine Ci(x) and sine integrals 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:
#
#    26 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).
#
  import numpy as np

  p2 = 1.570796326794897
  el = 0.5772156649015329
  eps = 1.0E-15
  x2 = x * x

  if ( x == 0.0 ):

    ci = - float ( "inf" )
    si = 0.0

  elif ( x <= 16.0 ):

    xr = -0.25 * x2
    ci = el + np.log ( x ) + xr
    for k in range ( 2, 41 ):
      xr = -0.5 * xr * ( k - 1 ) / ( k * k * ( 2 * k - 1 ) ) * x2
      ci = ci + xr
      if ( np.abs ( xr ) < np.abs ( ci ) * eps ):
        break

    xr = x
    si = x
    for k in range ( 1, 41 ):
      xr = -0.5 * xr * ( 2 * k - 1 ) / k / ( 4 * k * k + 4 * k + 1 ) * x2
      si = si + xr
      if ( np.abs ( xr ) < np.abs ( si ) * eps ):
        return ci, si

  elif ( x <= 32.0 ):

    bj = np.zeros ( 101 )
    m = int ( 47.2 + 0.82 * x )
    xa1 = 0.0
    xa0 = 1.0E-100
    for k in range ( m, 0, -1 ):
      xa = 4.0 * k * xa0 / x - xa1
      bj[k-1] = xa
      xa1 = xa0
      xa0 = xa

    xs = bj[0]
    for k in range ( 2, m, 2 ):
      xs = xs + 2.0 * bj[k]

    bj[0] = bj[0] / xs
    for k in range ( 1, m ):
      bj[k] = bj[k] / xs

    xr = 1.0
    xg1 = bj[0]
    for k in range ( 2, m + 1 ):
      xr = 0.25 * xr * ( 2.0 * k - 3.0 )**2  \
        / ( ( k - 1.0 ) * ( 2.0 * k - 1.0 )**2 ) * x
      xg1 = xg1 + bj[k-1] * xr

    xr = 1.0
    xg2 = bj[0]
    for k in range ( 2, m + 1 ):
      xr = 0.25 * xr * ( 2.0 * k - 5.0 )**2 \
        / ( ( k - 1.0 ) * ( 2.0 * k - 3.0 )**2 ) * x
      xg2 = xg2 + bj[k-1] * xr

    xcs = np.cos ( x / 2.0 )
    xss = np.sin ( x / 2.0 )
    ci = el + np.log ( x ) - x * xss * xg1 + 2.0 * xcs * xg2 - 2.0 * xcs * xcs
    si = x * xcs * xg1 + 2.0 * xss * xg2 - np.sin ( x )

  else:

    xr = 1.0
    xf = 1.0
    for k in range ( 1, 10 ):
      xr = -2.0 * xr * k * ( 2 * k - 1 ) / x2
      xf = xf + xr
    xr = 1.0 / x
    xg = xr
    for k in range ( 1, 9 ):
      xr = -2.0 * xr * ( 2 * k + 1 ) * k / x2
      xg = xg + xr
    ci = xf * np.sin ( x ) / x - xg * np.cos ( x ) / x
    si = p2 - xf * np.cos ( x ) / x - xg * np.sin ( x ) / x

  return ci, si

def ci_values ( n_data ):

#*****************************************************************************80
#
## ci_values() returns some values of the cosine integral function.
#
#  Discussion:
#
#    The cosine integral is defined by
#
#      CI(X) = - integral ( X <= T < Infinity ) ( cos ( T ) ) / T  dT
#
#    In Mathematica, the function can be evaluated by:
#
#      CosIntegral[x]
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    29 January 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Milton Abramowitz and Irene Stegun,
#    Handbook of Mathematical Functions,
#    US Department of Commerce, 1964.
#
#    Stephen Wolfram,
#    The Mathematica Book,
#    Fourth Edition,
#    Wolfram Media / Cambridge University Press, 1999.
#
#  Input:
#
#    integer n_data.  The user sets n_data to 0 before the first call.  
#
#  Output:
#
#    integer n_data.  On each call, the routine increments n_data by 1, and
#    returns the corresponding data; when there is no more data, the
#    output value of n_data will be 0 again.
#
#    real X, the argument of the function.
#
#    real FX, the value of the function.
#
  import numpy as np

  n_max = 16

  fx_vec = np.array ( ( \
     -0.1777840788066129E+00, \
     -0.2227070695927976E-01, \
      0.1005147070088978E+00, \
      0.1982786159524672E+00, \
      0.2760678304677729E+00, \
      0.3374039229009681E+00, \
      0.4204591828942405E+00, \
      0.4620065850946773E+00, \
      0.4717325169318778E+00, \
      0.4568111294183369E+00, \
      0.4229808287748650E+00, \
      0.2858711963653835E+00, \
      0.1196297860080003E+00, \
     -0.3212854851248112E-01, \
     -0.1409816978869304E+00, \
     -0.1934911221017388E+00 ) )

  x_vec = np.array ( ( \
      0.5E+00, \
      0.6E+00, \
      0.7E+00, \
      0.8E+00, \
      0.9E+00, \
      1.0E+00, \
      1.2E+00, \
      1.4E+00, \
      1.6E+00, \
      1.8E+00, \
      2.0E+00, \
      2.5E+00, \
      3.0E+00, \
      3.5E+00, \
      4.0E+00, \
      4.5E+00 ) )

  if ( n_data < 0 ):
    n_data = 0

  if ( n_max <= n_data ):
    n_data = 0
    x = 0.0
    fx = 0.0
  else:
    x = x_vec[n_data]
    fx = fx_vec[n_data]
    n_data = n_data + 1

  return n_data, x, fx

def si_values ( n_data ):

#*****************************************************************************80
#
## si_values() returns some values of the sine integral function.
#
#  Discussion:
#
#    SI(X) = integral ( 0 <= T <= X ) sin ( T ) / T dt
#
#    In Mathematica, the function can be evaluated by:
#
#      SinIntegral[x]
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    21 February 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Milton Abramowitz and Irene Stegun,
#    Handbook of Mathematical Functions,
#    US Department of Commerce, 1964.
#
#    Stephen Wolfram,
#    The Mathematica Book,
#    Fourth Edition,
#    Wolfram Media / Cambridge University Press, 1999.
#
#  Input:
#
#    integer n_data.  The user sets n_data to 0 before the first call.  
#
#  Output:
#
#    integer n_data.  On each call, the routine increments n_data by 1, and
#    returns the corresponding data; when there is no more data, the
#    output value of n_data will be 0 again.
#
#    real X, the argument of the function.
#
#    real F, the value of the function.
#
  import numpy as np

  n_max = 16

  f_vec = np.array ( ( \
     0.4931074180430667E+00, \
     0.5881288096080801E+00, \
     0.6812222391166113E+00, \
     0.7720957854819966E+00, \
     0.8604707107452929E+00, \
     0.9460830703671830E+00, \
     0.1108047199013719E+01, \
     0.1256226732779218E+01, \
     0.1389180485870438E+01, \
     0.1505816780255579E+01, \
     0.1605412976802695E+01, \
     0.1778520173443827E+01, \
     0.1848652527999468E+01, \
     0.1833125398665997E+01, \
     0.1758203138949053E+01, \
     0.1654140414379244E+01 ))

  x_vec = np.array ( ( \
      0.5E+00, \
      0.6E+00, \
      0.7E+00, \
      0.8E+00, \
      0.9E+00, \
      1.0E+00, \
      1.2E+00, \
      1.4E+00, \
      1.6E+00, \
      1.8E+00, \
      2.0E+00, \
      2.5E+00, \
      3.0E+00, \
      3.5E+00, \
      4.0E+00, \
      4.5E+00 ))

  if ( n_data < 0 ):
    n_data = 0

  if ( n_max <= n_data ):
    n_data = 0
    x = 0.0
    f = 0.0
  else:
    x = x_vec[n_data]
    f = f_vec[n_data]
    n_data = n_data + 1

  return n_data, x, f

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  sine_integral_test ( )
  timestamp ( )