subroutine monomial_value ( m, n, e, x, v ) !*****************************************************************************80 ! !! monomial_value() evaluates a monomial. ! ! Discussion: ! ! F(X) = product ( 1 <= I <= M ) X(I)^E(I) ! ! with the convention that 0^0 = 1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 13 April 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the spatial dimension. ! ! Input, integer N, the number of points. ! ! Input, integer E(M), the exponents. ! ! Input, real ( kind = rk ) X(M,N), the evaluation points. ! ! Output, real ( kind = rk ) V(N), the monomial values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer i integer e(m) real ( kind = rk ) v(n) real ( kind = rk ) x(m,n) v(1:n) = 1.0D+00 do i = 1, m if ( e(i) /= 0.0D+00 ) then v(1:n) = v(1:n) * x(i,1:n) ** e(i) end if end do return end subroutine pyramid01_integral ( expon, value ) !*****************************************************************************80 ! !! PYRAMID01_INTEGRAL: monomial integral in a unit pyramid. ! ! Discussion: ! ! This routine returns the integral of ! ! product ( 1 <= I <= 3 ) X(I)^EXPON(I) ! ! over the unit pyramid. ! ! The integration region is: ! ! - ( 1 - Z ) <= X <= 1 - Z ! - ( 1 - Z ) <= Y <= 1 - Z ! 0 <= Z <= 1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 March 2008 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Arthur Stroud, ! Approximate Calculation of Multiple Integrals, ! Prentice Hall, 1971, ! ISBN: 0130438936, ! LC: QA311.S85. ! ! Parameters: ! ! Input, integer EXPON(3), the exponents. ! ! Output, real ( kind = rk ) VALUE, the integral of the monomial. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer expon(3) integer i integer i_hi real ( kind = rk ) r8_choose real ( kind = rk ) r8_mop real ( kind = rk ) value value = 0.0D+00 if ( mod ( expon(1), 2 ) == 0 .and. mod ( expon(2), 2 ) == 0 ) then i_hi = 2 + expon(1) + expon(2) do i = 0, i_hi value = value + r8_mop ( i ) * r8_choose ( i_hi, i ) & / real ( i + expon(3) + 1, kind = rk ) end do value = value & * 2.0D+00 / real ( expon(1) + 1, kind = rk ) & * 2.0D+00 / real ( expon(2) + 1, kind = rk ) end if return end subroutine pyramid01_sample ( n, x ) !*****************************************************************************80 ! !! PYRAMID01_SAMPLE: sample the unit pyramid. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 13 April 2014 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of samples desired. ! ! Output, real ( kind = rk ) X(3,N), the sample values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ), parameter :: one_third = 1.0D+00 / 3.0D+00 real ( kind = rk ) x(3,n) call random_number ( harvest = x(1:3,1:n) ) x(3,1:n) = 1.0D+00 - x(3,1:n) ** one_third x(2,1:n) = ( 1.0D+00 - x(3,1:n) ) * ( 2.0D+00 * x(2,1:n) - 1.0D+00 ) x(1,1:n) = ( 1.0D+00 - x(3,1:n) ) * ( 2.0D+00 * x(1,1:n) - 1.0D+00 ) return end function pyramid01_volume ( ) !*****************************************************************************80 ! !! PYRAMID01_VOLUME: volume of a unit pyramid with square base. ! ! Discussion: ! ! The volume of this unit pyramid is 4/3. ! ! The integration region is: ! ! - ( 1 - Z ) <= X <= 1 - Z ! - ( 1 - Z ) <= Y <= 1 - Z ! 0 <= Z <= 1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 March 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) PYRAMID01_VOLUME, the volume. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) pyramid01_volume pyramid01_volume = 4.0D+00 / 3.0D+00 return end function r8_choose ( n, k ) !*****************************************************************************80 ! !! R8_CHOOSE computes the binomial coefficient C(N,K) as an R8. ! ! Discussion: ! ! The value is calculated in such a way as to avoid overflow and ! roundoff. The calculation is done in R8 arithmetic. ! ! The formula used is: ! ! C(N,K) = N! / ( K! * (N-K)! ) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 March 2008 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! ML Wolfson, HV Wright, ! Algorithm 160: ! Combinatorial of M Things Taken N at a Time, ! Communications of the ACM, ! Volume 6, Number 4, April 1963, page 161. ! ! Parameters: ! ! Input, integer N, K, are the values of N and K. ! ! Output, real ( kind = rk ) R8_CHOOSE, the number of combinations of N ! things taken K at a time. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer i integer k integer mn integer mx integer n real ( kind = rk ) r8_choose real ( kind = rk ) value mn = min ( k, n - k ) if ( mn < 0 ) then value = 0.0D+00 else if ( mn == 0 ) then value = 1.0D+00 else mx = max ( k, n - k ) value = real ( mx + 1, kind = rk ) do i = 2, mn value = ( value * real ( mx + i, kind = rk ) ) / real ( i, kind = rk ) end do end if r8_choose = value return end function r8_mop ( i ) !*****************************************************************************80 ! !! R8_MOP returns the I-th power of -1 as an R8. ! ! Discussion: ! ! An R8 is a real ( kind = rk ) value. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 07 November 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the power of -1. ! ! Output, real ( kind = rk ) R8_MOP, the I-th power of -1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer i real ( kind = rk ) r8_mop if ( mod ( i, 2 ) == 0 ) then r8_mop = + 1.0D+00 else r8_mop = - 1.0D+00 end if return end subroutine r8mat_transpose_print ( m, n, a, title ) !*****************************************************************************80 ! !! R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. ! ! Discussion: ! ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 June 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), an M by N matrix to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) character ( len = * ) title call r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. ! ! Discussion: ! ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 10 September 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: incx = 5 integer m integer n real ( kind = rk ) a(m,n) character ( len = 14 ) ctemp(incx) integer i integer i2 integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) if ( m <= 0 .or. n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (None)' return end if do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i8,6x)' ) i end do write ( *, '('' Row '',5a14)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' write ( *, '(a)' ) ' ' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(g14.6)' ) a(i,j) end do write ( *, '(i5,a,5a14)' ) j, ':', ( ctemp(i), i = 1, inc ) end do end do return end subroutine r8vec_print ( n, a, title ) !*****************************************************************************80 ! !! R8VEC_PRINT prints an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 August 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, real ( kind = rk ) A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a(n) integer i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,a,1x,g16.8)' ) i, ':', a(i) end do 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 ! ! Parameters: ! ! None ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) 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