program main c*********************************************************************72 c cc polygon_monte_carlo_test() tests polygon_monte_carlo(). c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 06 May 2014 c c Author: c c John Burkardt c implicit none integer nv1 parameter ( nv1 = 4 ) double precision v1(2,nv1) save v1 data v1 / & -1.0, -1.0, & 1.0, -1.0, & 1.0, 1.0, & -1.0, 1.0 / call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'polygon_monte_carlo_test():' write ( *, '(a)' ) ' FORTRAN77 version' write ( *, '(a)' ) ' Test polygon_monte_carlo().' call test01 ( nv1, v1 ) c c Terminate. c write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'polygon_monte_carlo_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine test01 ( nv, v ) c*****************************************************************************80 c cc TEST01 estimates integrals over a polygon in 2D. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 06 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 65536 ) integer nv integer e(2) integer e_test(2,7) integer j integer n double precision polygon_area double precision r8vec_sum double precision result(7) integer seed double precision v(2,nv) double precision value(n_max) double precision x(2,n_max) save e_test data e_test / & 0, 0, & 2, 0, & 0, 2, & 4, 0, & 2, 2, & 0, 4, & 6, 0 / write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) ' Use POLYGON_SAMPLE to estimate integrals ' write ( *, '(a)' ) ' over the interior of a polygon in 2D.' seed = 123456789 write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' N' // & ' 1' // & ' X^2 ' // & ' Y^2' // & ' X^4' // & ' X^2Y^2' // & ' Y^4' // & ' X^6' write ( *, '(a)' ) ' ' n = 1 10 continue if ( n .le. 65536 ) then call polygon_sample ( nv, v, n, seed, x ) do j = 1, 7 e(1:2) = e_test(1:2,j) call monomial_value ( 2, n, e, x, value ) result(j) = polygon_area ( nv, v ) * r8vec_sum ( n, value ) & / dble ( n ) end do write ( *, '(2x,i8,7(2x,g14.6))' ) n, result(1:7) n = 2 * n go to 10 end if write ( *, '(a)' ) ' ' do j = 1, 7 e(1:2) = e_test(1:2,j) call polygon_monomial_integral ( nv, v, e, result(j) ) end do write ( *, '(2x,a8,7(2x,g14.6))' ) ' Exact', result(1:7) return end subroutine timestamp ( ) c*********************************************************************72 c cc TIMESTAMP prints out the current YMDHMS date as a timestamp. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 January 2007 c c Author: c c John Burkardt c c Parameters: c c None c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 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, month(m), y, h, ':', n, ':', s, '.', mm, ampm return end