program main !*****************************************************************************80 ! !! clausen_test() tests clausen(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 06 September 2021 ! ! Author: ! ! John Burkardt ! implicit none call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'clausen_test():' write ( *, '(a)' ) ' Fortran90 version' write ( *, '(a)' ) ' Test clausen().' call clausen_values_test ( ) call clausen_plot_test ( ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'clausen_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop 0 end subroutine clausen_values_test ( ) !*****************************************************************************80 ! !! clausen_values_test() tests clausen() against a table of values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 06 September 2021 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) real ( kind = rk8 ) clausen real ( kind = rk8 ) diff real ( kind = rk8 ) fx1 real ( kind = rk8 ) fx2 integer n_data real ( kind = rk8 ) x write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'clausen_values_test():' write ( *, '(a)' ) ' clausen() evaluates the Clausen function.' write ( *, '(a)' ) ' Compare its results to tabulated data.' write ( *, '(a)' ) '' write ( *, '(a)' ) ' Tabulated Computed' write ( *, '(a)' ) ' X FX FX Diff' write ( *, '(a)' ) '' n_data = 0 do call clausen_values ( n_data, x, fx1 ) if ( n_data == 0 ) then exit end if fx2 = clausen ( x ) diff = abs ( fx1 - fx2 ) write ( *, '(2x,g14.6,2x,g24.16,2x,g24.16,2x,e7.1)' ) x, fx1, fx2, diff end do return end subroutine clausen_plot_test ( ) !*****************************************************************************80 ! !! clausen_plot_test() plots the clausen function. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 March 2025 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) real ( kind = rk8 ) a real ( kind = rk8 ) b real ( kind = rk8 ) clausen real ( kind = rk8 ), allocatable :: fx(:) character ( len = 50 ) header integer i integer nx real ( kind = rk8 ), parameter :: r8_pi = 3.141592653589793D+00 real ( kind = rk8 ), allocatable :: x(:) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'clausen_plot_test():' write ( *, '(a)' ) ' Plot the clausen() function.' nx = 201 a = 0.0D+00 b = 4.0D+00 * r8_pi allocate ( x(1:nx) ) call r8vec_linspace ( nx, a, b, x ) allocate ( fx(1:nx) ) do i = 1, nx fx(i) = clausen ( x(i) ) end do header = 'clausen_test' call gnuplot_fx ( nx, x, fx, trim ( header ) ) deallocate ( fx ) deallocate ( x ) return end subroutine clausen_values ( n_data, x, fx ) !*****************************************************************************80 ! !! clausen_values() returns some values of the Clausen's integral. ! ! Discussion: ! ! The function is defined by: ! ! CLAUSEN(x) = integral ( 0 <= t <= x ) -ln ( 2 * sin ( t / 2 ) ) dt ! ! The data was reported by McLeod. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 06 September 2021 ! ! 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. ! ! Allan McLeod, ! Algorithm 757: ! MISCFUN: A software package to compute uncommon special functions, ! ACM Transactions on Mathematical Software, ! Volume 22, Number 3, September 1996, pages 288-301. ! ! 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 ( kind = rk8 ) X, the argument of the function. ! ! real ( kind = rk8 ) FX, the value of the function. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: n_max = 20 real ( kind = rk8 ) fx real ( kind = rk8 ), save, dimension ( n_max ) :: fx_vec = (/ & 0.14137352886760576684D-01, & 0.13955467081981281934D+00, & -0.38495732156574238507D+00, & 0.84831187770367927099D+00, & 0.10139591323607685043D+01, & -0.93921859275409211003D+00, & 0.72714605086327924743D+00, & 0.43359820323553277936D+00, & -0.98026209391301421161D-01, & -0.56814394442986978080D+00, & -0.70969701784448921625D+00, & 0.99282013254695671871D+00, & -0.98127747477447367875D+00, & -0.64078266570172320959D+00, & 0.86027963733231192456D+00, & 0.39071647608680211043D+00, & 0.47574793926539191502D+00, & 0.10105014481412878253D+01, & 0.96332089044363075154D+00, & -0.61782699481929311757D+00 /) integer n_data real ( kind = rk8 ) x real ( kind = rk8 ), save, dimension ( n_max ) :: x_vec = (/ & 0.0019531250D+00, & 0.0312500000D+00, & -0.1250000000D+00, & 0.5000000000D+00, & 1.0000000000D+00, & -1.5000000000D+00, & 2.0000000000D+00, & 2.5000000000D+00, & -3.0000000000D+00, & 4.0000000000D+00, & 4.2500000000D+00, & -5.0000000000D+00, & 5.5000000000D+00, & 6.0000000000D+00, & 8.0000000000D+00, & -10.0000000000D+00, & 15.0000000000D+00, & 20.0000000000D+00, & -30.0000000000D+00, & 50.0000000000D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! get_unit() returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is a value between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is a value between 1 and 99, representing a ! free FORTRAN unit. The code assumes that units 5 and 6 ! are special, and will never return those values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 October 2008 ! ! Author: ! ! John Burkardt ! ! Output: ! ! integer IUNIT, the free unit number. ! implicit none integer i integer ios integer iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end subroutine gnuplot_fx ( n, x, y, header ) !*****************************************************************************80 ! !! gnuplot_fx() plots x, f(x) data. ! ! Discussion: ! ! Actually, we simply create two files for processing by gnuplot(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 March 2025 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer n character ( len = 255 ) command_filename integer command_unit character ( len = 255 ) data_filename integer data_unit character ( len = * ) header integer i character ( len = 255 ) png_filename real ( kind = rk8 ) x(n) real ( kind = rk8 ) y(n) ! ! Create a graphics data file. ! data_filename = header // '_data.txt' call get_unit ( data_unit ) open ( unit = data_unit, file = data_filename, status = 'replace' ) do i = 1, n write ( data_unit, '(2g14.6)' ) x(i), y(i) end do close ( unit = data_unit ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Created graphics data file "' & // trim ( data_filename ) // '".' png_filename = header // '.png' ! ! Create graphics command file. ! call get_unit ( command_unit ) command_filename = header // '_commands.txt' open ( unit = command_unit, file = command_filename, status = 'replace' ) write ( command_unit, '(a)' ) '# ' // trim ( command_filename ) write ( command_unit, '(a)' ) '#' write ( command_unit, '(a)' ) '# Usage:' write ( command_unit, '(a)' ) '# gnuplot < ' // trim ( command_filename ) write ( command_unit, '(a)' ) '#' write ( command_unit, '(a)' ) 'set term png' write ( command_unit, '(a)' ) 'set output "' // trim ( png_filename ) // '"' write ( command_unit, '(a)' ) 'set xlabel "<--- X --->"' write ( command_unit, '(a)' ) 'set ylabel "<-- Y(X) -->"' write ( command_unit, '(a)' ) 'set title "' // header // '"' write ( command_unit, '(a)' ) 'set grid' write ( command_unit, '(a)' ) 'set style data lines' write ( command_unit, '(a)' ) 'plot "' // trim ( data_filename ) // & '" using 1:2 lw 3 linecolor rgb "blue"' close ( unit = command_unit ) write ( *, '(a)' ) & ' Created command file "' // trim ( command_filename ) // '".' return end subroutine r8vec_linspace ( n, a, b, x ) !*****************************************************************************80 ! !! r8vec_linspace() 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: ! ! 14 March 2011 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer N, the number of entries in the vector. ! ! real ( kind = rk8 ) A, B, the first and last entries. ! ! Output: ! ! real ( kind = rk8 ) X(N), a vector of linearly spaced data. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer n real ( kind = rk8 ) a real ( kind = rk8 ) b integer i real ( kind = rk8 ) x(n) if ( n == 1 ) then x(1) = ( a + b ) / 2.0D+00 else do i = 1, n x(i) = ( real ( n - i, kind = rk8 ) * a & + real ( i - 1, kind = rk8 ) * b ) & / real ( n - 1, kind = rk8 ) end do end if return end subroutine timestamp ( ) !*****************************************************************************80 ! !! timestamp() prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 06 September 2021 ! ! Author: ! ! John Burkardt ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end