subroutine r8vec_uniform_01 ( n, seed, r ) !*****************************************************************************80 ! !! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 13 August 2014 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! Bennett Fox, ! Algorithm 647: ! Implementation and Relative Efficiency of Quasirandom ! Sequence Generators, ! ACM Transactions on Mathematical Software, ! Volume 12, Number 4, pages 362-376, 1986. ! ! Peter Lewis, Allen Goodman, James Miller ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, pages 136-143, 1969. ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R(N), the vector of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer i integer, parameter :: i4_huge = 2147483647 integer k integer seed real ( kind = rk ) r(n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8VEC_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r(i) = real ( seed, kind = rk ) * 4.656612875D-10 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 function triangle_area ( t ) !*****************************************************************************80 ! !! TRIANGLE_AREA computes the area of a triangle in 2D. ! ! Discussion: ! ! If the triangle's vertices are given in counter clockwise order, ! the area will be positive. If the triangle's vertices are given ! in clockwise order, the area will be negative! ! ! An earlier version of this routine always returned the absolute ! value of the computed area. I am convinced now that that is ! a less useful result! For instance, by returning the signed ! area of a triangle, it is possible to easily compute the area ! of a nonconvex polygon as the sum of the (possibly negative) ! areas of triangles formed by node 1 and successive pairs of vertices. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) T(2,3), the triangle vertices. ! ! Output, real ( kind = rk ) TRIANGLE_AREA, the area of the triangle. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) triangle_area real ( kind = rk ) t(2,3) triangle_area = 0.5D+00 * ( & t(1,1) * ( t(2,2) - t(2,3) ) & + t(1,2) * ( t(2,3) - t(2,1) ) & + t(1,3) * ( t(2,1) - t(2,2) ) ) return end subroutine triangle_interpolate_linear ( m, n, p1, p2, p3, p, v1, v2, v3, v ) !*****************************************************************************80 ! !! TRIANGLE_INTERPOLATE_LINEAR interpolates data given on a triangle's vertices. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 January 2015 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the dimension of the value to be interpolated. ! ! Input, integer N, the number of points. ! ! Input, real P1(2), P2(2), P3(2), the vertices of the triangle, ! in counterclockwise order. ! ! Input, real P(2,N), the point at which the interpolant is desired. ! ! Input, real V1(M), V2(M), V3(M), the value at the vertices. ! ! Output, real V(M,N), the interpolated value of the quantity at P. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer j real ( kind = rk ) p(2,n) real ( kind = rk ) p1(2) real ( kind = rk ) p2(2) real ( kind = rk ) p3(2) real ( kind = rk ) triangle_area real ( kind = rk ) v(m,n) real ( kind = rk ) v1(m) real ( kind = rk ) v2(m) real ( kind = rk ) v3(m) v(1:m,1:n) = 0.0D+00 do j = 1, n v(1:m,j) = & ( triangle_area ( (/ p(1:2,j), p2, p3 /) ) * v1(1:m) & + triangle_area ( (/ p1, p(1:2,j), p3 /) ) * v2(1:m) & + triangle_area ( (/ p1, p2, p(1:2,j) /) ) * v3(1:m) ) & / triangle_area ( (/ p1, p2, p3 /) ) end do return end subroutine uniform_in_triangle_map1 ( v1, v2, v3, n, seed, x ) !*****************************************************************************80 ! !! UNIFORM_IN_TRIANGLE_MAP1 maps uniform points into a triangle. ! ! Discussion: ! ! The triangle is defined by three vertices. This routine ! uses Turk's rule #1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Greg Turk, ! Generating Random Points in a Triangle, ! in Graphics Gems, ! edited by Andrew Glassner, ! AP Professional, 1990, pages 24-28. ! ! Parameters: ! ! Input, real ( kind = rk ) V1(2), V2(2), V3(2), the vertices. ! ! Input, integer N, the number of points. ! ! Input/output, integer SEED, a seed for the random ! number generator. ! ! Output, real ( kind = rk ) X(2,N), the points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: dim_num = 2 integer n real ( kind = rk ) a real ( kind = rk ) b real ( kind = rk ) c integer j real ( kind = rk ) r(3) integer seed real ( kind = rk ) x(dim_num,n) real ( kind = rk ) v1(dim_num) real ( kind = rk ) v2(dim_num) real ( kind = rk ) v3(dim_num) ! ! Generate the points using Turk's rule 1. ! do j = 1, n call r8vec_uniform_01 ( 2, seed, r ) a = 1.0D+00 - sqrt ( r(2) ) b = ( 1.0D+00 - r(1) ) * sqrt ( r(2) ) c = r(1) * sqrt ( r(2) ) x(1:dim_num,j) = a * v1(1:dim_num) & + b * v2(1:dim_num) & + c * v3(1:dim_num) end do return end