# include <stdbool.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>

# include "sine_integral.h"

int main ( );
void ci_test ( );
void ci_values ( int *n_data, double *x, double *fx );
void si_test ( );
void si_values ( int *n_data, double *x, double *fx );
void timestamp ( );

/******************************************************************************/

int main ( )

/******************************************************************************/
/*
  Purpose:

    sine_integral_test() tests sine_integral().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    27 March 2025

  Author:

    John Burkardt
*/
{
  timestamp ( );
  printf ( "\n" );
  printf ( "sine_integral_test():\n" );
  printf ( "  C version\n" );
  printf ( "  Test sine_integral().\n" );

  ci_test ( );
  si_test ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "sine_integral_test():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

  return 0;
}
/******************************************************************************/

void ci_test ( )

/******************************************************************************/
/*
  Purpose:

    ci_test() tests ci().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    27 March 2025

  Author:

    John Burkardt
*/
{
  double fx1;
  double fx2;
  int n_data;
  double si;
  double x;

  printf ( "\n" );
  printf ( "ci_test():\n" );
  printf ( "  ci() evaluates the Cosine Integral function CI(X).\n" );
  printf ( "\n" );
  printf ( "          X           CI(X)\n" );
  printf ( "\n" );

  n_data = 0;

  while ( true )
  {
    ci_values ( &n_data, &x, &fx1 );

    if ( n_data == 0 )
    {
      break;
    }

    cisi ( x, &fx2, &si );

    printf ( "  %14.6g  %24.16g  %24.16g\n", x, fx1, fx2 );
  }

  return;
}
/******************************************************************************/

void ci_values ( int *n_data, double *x, double *fx )

/******************************************************************************/
/*
  Purpose:

    ci_values() returns some values of the cosine integral function.

  Discussion:

    The cosine integral is defined by

      CI(X) = - integral ( X <= T < +oo ) ( cos ( T ) ) / T  dT

    In Mathematica, the function can be evaluated by:

      CosIntegral[x]

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    12 August 2004

  Author:

    John Burkardt

  Reference:

    Milton Abramowitz, Irene Stegun,
    Handbook of Mathematical Functions,
    National Bureau of Standards, 1964,
    ISBN: 0-486-61272-4,
    LC: QA47.A34.

    Stephen Wolfram,
    The Mathematica Book,
    Fourth Edition,
    Cambridge University Press, 1999,
    ISBN: 0-521-64314-7,
    LC: QA76.95.W65.

  Input:

    int *N_DATA.  The user sets N_DATA to 0 before the first call.

  Output:

    int *N_DATA.  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.

    double *X, the argument of the function.

    double *FX, the value of the function.
*/
{
# define N_MAX 16

  static double fx_vec[N_MAX] = {
     -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 };

  static double x_vec[N_MAX] = {
      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;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

  return;
# undef N_MAX
}
/******************************************************************************/

void si_test ( )

/******************************************************************************/
/*
  Purpose:

    si_test() tests si().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    27 March 2025

  Author:

    John Burkardt
*/
{
  double ci;
  double fx1;
  double fx2;
  int n_data;
  double x;

  printf ( "\n" );
  printf ( "si_test():\n" );
  printf ( "  si() evaluates the sine integral function.\n" );
  printf ( "\n" );
  printf ( "      X            SI(X)\n" );
  printf ( "\n" );

  n_data = 0;

  while ( true )
  {
    si_values ( &n_data, &x, &fx1 );

    if ( n_data == 0 )
    {
      break;
    }

    cisi ( x, &ci, &fx2 );

    printf ( "  %14.6g  %24.16g  %24.16g\n", x, fx1, fx2 );

  }

  return;
}
/******************************************************************************/

void si_values ( int *n_data, double *x, double *fx )

/******************************************************************************/
/*
  Purpose:

    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:

    12 August 2004

  Author:

    John Burkardt

  Reference:

    Milton Abramowitz, Irene Stegun,
    Handbook of Mathematical Functions,
    National Bureau of Standards, 1964,
    ISBN: 0-486-61272-4,
    LC: QA47.A34.

    Stephen Wolfram,
    The Mathematica Book,
    Fourth Edition,
    Cambridge University Press, 1999,
    ISBN: 0-521-64314-7,
    LC: QA76.95.W65.

  Input:

    int *N_DATA.  The user sets N_DATA to 0 before the first call.

  Output:

    int *N_DATA.  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.

    double *X, the argument of the function.

    double *FX, the value of the function.
*/
{
# define N_MAX 16

  static double fx_vec[N_MAX] = {
     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 };

  static double x_vec[N_MAX] = {
      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;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

  return;
# undef N_MAX
}
/******************************************************************************/

void timestamp ( )

/******************************************************************************/
/*
  Purpose:
 
    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 May 2021

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}