function alogam ( x, ifault ) !*****************************************************************************80 ! !! alogam() computes the logarithm of the Gamma function. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 March 1999 ! ! Author: ! ! Malcolm Pike, David Hill ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Malcolm Pike, David Hill, ! Algorithm 291: ! Logarithm of Gamma Function, ! Communications of the ACM, ! Volume 9, Number 9, September 1966, page 684. ! ! Parameters: ! ! Input, real ( kind = rk ) X, the argument of the Gamma function. ! X should be greater than 0. ! ! Output, integer IFAULT, error flag. ! 0, no error. ! 1, X <= 0. ! ! Output, real ( kind = rk ) ALOGAM, the logarithm of the Gamma ! function of X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) alogam real ( kind = rk ) f integer ifault real ( kind = rk ) x real ( kind = rk ) y real ( kind = rk ) z if ( x <= 0.0D+00 ) then ifault = 1 alogam = 0.0D+00 return end if ifault = 0 y = x if ( x < 7.0D+00 ) then f = 1.0D+00 z = y do while ( z < 7.0D+00 ) f = f * z z = z + 1.0D+00 end do y = z f = - log ( f ) else f = 0.0D+00 end if z = 1.0D+00 / y / y alogam = f + ( y - 0.5D+00 ) * log ( y ) - y & + 0.918938533204673D+00 + & ((( & - 0.000595238095238D+00 * z & + 0.000793650793651D+00 ) * z & - 0.002777777777778D+00 ) * z & + 0.083333333333333D+00 ) / y return end subroutine gamma_log_values ( n_data, x, fx ) !*****************************************************************************80 ! !! GAMMA_LOG_VALUES returns some values of the Log Gamma function. ! ! Discussion: ! ! In Mathematica, the function can be evaluated by: ! ! Log[Gamma[x]] ! ! Modified: ! ! 14 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. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 ! before the first call. 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. ! ! Output, real ( kind = rk ) X, the argument of the function. ! ! Output, real ( kind = rk ) FX, the value of the function. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n_max = 20 real ( kind = rk ) fx real ( kind = rk ), save, dimension ( n_max ) :: fx_vec = (/ & 0.1524063822430784D+01, & 0.7966778177017837D+00, & 0.3982338580692348D+00, & 0.1520596783998375D+00, & 0.0000000000000000D+00, & -0.4987244125983972D-01, & -0.8537409000331584D-01, & -0.1081748095078604D+00, & -0.1196129141723712D+00, & -0.1207822376352452D+00, & -0.1125917656967557D+00, & -0.9580769740706586D-01, & -0.7108387291437216D-01, & -0.3898427592308333D-01, & 0.00000000000000000D+00, & 0.69314718055994530D+00, & 0.17917594692280550D+01, & 0.12801827480081469D+02, & 0.39339884187199494D+02, & 0.71257038967168009D+02 /) integer n_data real ( kind = rk ) x real ( kind = rk ), save, dimension ( n_max ) :: x_vec = (/ & 0.20D+00, & 0.40D+00, & 0.60D+00, & 0.80D+00, & 1.00D+00, & 1.10D+00, & 1.20D+00, & 1.30D+00, & 1.40D+00, & 1.50D+00, & 1.60D+00, & 1.70D+00, & 1.80D+00, & 1.90D+00, & 2.00D+00, & 3.00D+00, & 4.00D+00, & 10.00D+00, & 20.00D+00, & 30.00D+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 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: ! ! 18 May 2013 ! ! 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,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