# include # include # include # include # include "sine_integral.h" int main ( ); void si_plot ( ); void si_test ( ); double *r8vec_linspace_new ( int n, double a, double b ); 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: 10 August 2025 Author: John Burkardt */ { timestamp ( ); printf ( "\n" ); printf ( "sine_integral_test():\n" ); printf ( " C version\n" ); printf ( " Test sine_integral().\n" ); si_plot ( ); si_test ( ); /* Terminate. */ printf ( "\n" ); printf ( "sine_integral_test():\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void si_plot ( ) /******************************************************************************/ /* Purpose: si_plot() plots si(x). Licensing: This code is distributed under the MIT license. Modified: 10 August 2025 Author: John Burkardt */ { char command_filename[255]; FILE *command_unit; char data_filename[255]; FILE *data_unit; int i; int nplot; char output_filename[255]; char *prefix = "si"; double *xplot; double *yplot; nplot = 101; xplot = r8vec_linspace_new ( nplot, -10.0, +10.0 ); yplot = ( double * ) malloc ( nplot * sizeof ( double ) ); for ( i = 0; i < nplot; i++ ) { yplot[i] = si ( xplot[i] ); } /* Create the data file. */ sprintf ( data_filename, "%s_data.txt", prefix ); data_unit = fopen ( data_filename, "wt" ); for ( i = 0; i < nplot; i++ ) { fprintf ( data_unit, "%g %g\n", xplot[i], yplot[i] ); } fclose ( data_unit ); printf ( "\n" ); printf ( " Created graphics data file '%s'.\n", data_filename ); /* Plot the selected data. */ sprintf ( command_filename, "%s_commands.txt", prefix ); command_unit = fopen ( command_filename, "wt" ); fprintf ( command_unit, "# %s\n", command_filename ); fprintf ( command_unit, "#\n" ); fprintf ( command_unit, "# Usage:\n" ); fprintf ( command_unit, "# gnuplot < %s\n", command_filename ); fprintf ( command_unit, "#\n" ); fprintf ( command_unit, "set term png\n" ); fprintf ( command_unit, "set nokey\n" ); sprintf ( output_filename, "%s_plot.png", prefix ); fprintf ( command_unit, "set output '%s'\n", output_filename ); fprintf ( command_unit, "set xlabel '<---X--->'\n" ); fprintf ( command_unit, "set ylabel '<---Si(X)--->'\n" ); fprintf ( command_unit, "set title 'Sine integral Si(x)'\n" ); fprintf ( command_unit, "set grid\n" ); fprintf ( command_unit, "set style data lines\n" ); fprintf ( command_unit, "plot '%s' using 1:2 lw 3 linecolor rgb 'blue'\n", data_filename ); fclose ( command_unit ); printf ( " Created graphics command file '%s'.\n", command_filename ); free ( xplot ); free ( yplot ); return; } /******************************************************************************/ double *r8vec_linspace_new ( int n, double a, double b ) /******************************************************************************/ /* Purpose: r8vec_linspace_new() creates a vector of linearly spaced values. Discussion: An R8VEC is a vector of R8's. 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12. In other words, the interval is divided into N-1 even subintervals, and the endpoints of intervals are used as the points. Licensing: This code is distributed under the MIT license. Modified: 29 March 2011 Author: John Burkardt Input: int N, the number of entries in the vector. double A, B, the first and last entries. Output: double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data. */ { int i; double *x; x = ( double * ) malloc ( n * sizeof ( double ) ); if ( n == 1 ) { x[0] = ( a + b ) / 2.0; } else { for ( i = 0; i < n; i++ ) { x[i] = ( ( double ) ( n - 1 - i ) * a + ( double ) ( i ) * b ) / ( double ) ( n - 1 ); } } return x; } /******************************************************************************/ void si_test ( ) /******************************************************************************/ /* Purpose: si_test() tests si(). Licensing: This code is distributed under the MIT license. Modified: 10 August 2025 Author: John Burkardt */ { 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; } fx2 = si ( x ); 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 }