# include # include # include # include # include # include "cosine_integral.h" int main ( ); void ci_plot ( ); void ci_test ( ); double *r8vec_linspace_new ( int n, double a, double b ); void ci_values ( int *n_data, double *x, double *fx ); void timestamp ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: cosine_integral_test() tests cosine_integral(). Licensing: This code is distributed under the MIT license. Modified: 11 August 2025 Author: John Burkardt */ { timestamp ( ); printf ( "\n" ); printf ( "cosine_integral_test():\n" ); printf ( " C version\n" ); printf ( " Test cosine_integral().\n" ); ci_plot ( ); ci_test ( ); /* Terminate. */ printf ( "\n" ); printf ( "cosine_integral_test():\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void ci_plot ( ) /******************************************************************************/ /* Purpose: ci_plot() plots ci(x). Licensing: This code is distributed under the MIT license. Modified: 11 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 = "ci"; double *xplot; double *yplot; /* Idiotic gnuplot can't handle plot with bad data (0,-oo) so replace by (0,-10) and restrict yrange of plot. */ nplot = 101; xplot = r8vec_linspace_new ( nplot, -10.0, +10.0 ); yplot = ( double * ) malloc ( nplot * sizeof ( double ) ); for ( i = 0; i < nplot; i++ ) { if ( fabs ( xplot[i] ) < 0.0001 ) { yplot[i] = -10.0; } else { yplot[i] = ci ( 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 '<---Ci(X)--->'\n" ); fprintf ( command_unit, "set yrange [ -1.0 : +1.0 ]\n" ); fprintf ( command_unit, "set title 'Cosine integral Ci(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 ci_test ( ) /******************************************************************************/ /* Purpose: ci_test() tests ci(). 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 ( "ci_test():\n" ); printf ( " ci() evaluates the cosine integral function.\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; } fx2 = ci ( x ); 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 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 }