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

#*****************************************************************************80
#
## cosine_integral_test() tests cosine_integral().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    08 August 2025
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'cosine_integral_test():' )
  print ( '  python version: ' + platform.python_version ( ) )
  print ( '  numpy version:  ' + np.version.version )
  print ( '  Test cosine_integral().' )

  ci_test ( );
  ci_plot ( );
#
#  Terminate.
#
  print ( '' )
  print ( 'cosine_integral_test():' )
  print ( '  Normal end of execution.' )

  return

def ci ( x ):

#*****************************************************************************80
#
## ci() computes the cosine integral Ci(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:
#
#    09 August 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.
#    0.0 <= X.
#
#  Output:
#
#    real value: the value of Ci(x).
#
  import numpy as np

  p2 = 1.570796326794897
  el = 0.5772156649015329
  eps = 1.0E-15
  x2 = x * x
  xabs = np.abs ( x )

  if ( xabs == 0.0 ):

    value = - float ( "inf" )

  elif ( xabs <= 16.0 ):

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

  elif ( xabs <= 32.0 ):

    bj = np.zeros ( 101 )
    m = int ( 47.2 + 0.82 * xabs )
    xa1 = 0.0
    xa0 = 1.0E-100
    for k in range ( m, 0, -1 ):
      xa = 4.0 * k * xa0 / xabs - 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 ) * xabs
      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 ) * xabs
      xg2 = xg2 + bj[k-1] * xr

    xcs = np.cos ( xabs / 2.0 )
    xss = np.sin ( xabs / 2.0 )
    value = el + np.log ( xabs ) - xabs * xss * xg1 + 2.0 * xcs * xg2 - 2.0 * xcs * xcs

  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 / xabs
    xg = xr
    for k in range ( 1, 9 ):
      xr = -2.0 * xr * ( 2 * k + 1 ) * k / x2
      xg = xg + xr
    value = xf * np.sin ( xabs ) / xabs - xg * np.cos ( xabs ) / xabs

  return value

def ci_plot ( ):

#*****************************************************************************80
#
## ci_plot() plots ci().
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    08 August 2025
#
#  Author:
#
#    John Burkardt
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'ci_plot():' )
  print ( '  ci() evaluates the cosine integral function.' )

  a = -10.0
  b =  10.0
  nplot = 101
  x = np.linspace ( a, b, nplot )
  y = np.zeros ( nplot )
  for i in range ( 0, nplot ):
    y[i] = ci ( x[i] )

  plt.plot ( x, y, 'b-', linewidth = 2 )
  plt.plot ( [a,b], [0,0], 'k-', linewidth = 2 )

  ymin = min ( y )
  ymax = max ( y )
  plt.plot ( [0, 0], [ymin,ymax], 'k-', linewidth = 2 )

  plt.grid ( True )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Ci(x) --->' )
  plt.title ( 'The cosine integral function' )
  filename = 'ci_plot.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( block = False )
  plt.close ( )

  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 = ci ( x )

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

  return

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 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 ( )
  cosine_integral_test ( )
  timestamp ( )