program main !*****************************************************************************80 ! !! MAIN is the main program for STRI_QUAD_PRB. ! ! Discussion: ! ! STRI_QUAD_PRB tests STRI_QUAD. ! ! Modified: ! ! 22 October 2006 ! ! Author: ! ! John Burkardt ! implicit none call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'STRI_QUAD_PRB' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' Test the routines in STRI_QUAD.' call test00 call test01 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'STRI_QUAD_PRB' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine test00 !*****************************************************************************80 ! !! TEST00 takes care of initializing the random number seed for us. ! implicit none integer :: seed = 0 call random_initialize ( seed ) return end subroutine test01 !*****************************************************************************80 ! !! TEST01 tests SPHERE_QUAD_**. ! implicit none integer e(3) real ( kind = 8 ), external :: polyterm_value_3d real ( kind = 8 ) h integer h_test integer i real ( kind = 8 ), parameter :: r = 1.0E+00 real ( kind = 8 ) result_00 real ( kind = 8 ) result_01 real ( kind = 8 ) result_02 real ( kind = 8 ) result_03 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Approximate the integral of a function on a sphere.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' SPHERE_QUAD_00 uses a Monte Carlo method.' write ( *, '(a)' ) ' SPHERE_QUAD_01 uses centroids of spherical triangles.' write ( *, '(a)' ) ' SPHERE_QUAD_02 uses vertices of spherical triangles.' write ( *, '(a)' ) ' SPHERE_QUAD_03 uses midsides of spherical triangles.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' H QUAD_00 QUAD_01 QUAD_02 QUAD_03' write ( *, '(a)' ) ' ' do i = 1, 17 if ( i == 1 ) then e(1:3) = (/ 0, 0, 0 /) else if ( i == 2 ) then e(1:3) = (/ 1, 0, 0 /) else if ( i == 3 ) then e(1:3) = (/ 0, 1, 0 /) else if ( i == 4 ) then e(1:3) = (/ 0, 0, 1 /) else if ( i == 5 ) then e(1:3) = (/ 2, 0, 0 /) else if ( i == 6 ) then e(1:3) = (/ 0, 2, 2 /) else if ( i == 7 ) then e(1:3) = (/ 2, 2, 2 /) else if ( i == 8 ) then e(1:3) = (/ 0, 2, 4 /) else if ( i == 9 ) then e(1:3) = (/ 0, 0, 6 /) else if ( i == 10 ) then e(1:3) = (/ 1, 2, 4 /) else if ( i == 11 ) then e(1:3) = (/ 2, 4, 2 /) else if ( i == 12 ) then e(1:3) = (/ 6, 2, 0 /) else if ( i == 13 ) then e(1:3) = (/ 0, 0, 8 /) else if ( i == 14 ) then e(1:3) = (/ 6, 0, 4 /) else if ( i == 15 ) then e(1:3) = (/ 4, 6, 2 /) else if ( i == 16 ) then e(1:3) = (/ 2, 4, 8 /) else if ( i == 17 ) then e(1:3) = (/ 16, 0, 0 /) end if call polyterm_exponent ( 'SET', e ) call polyterm_exponent ( 'PRINT', e ) do h_test = 1, 2 if ( h_test == 1 ) then h = real ( 1.0, kind = 8 ) else if ( h_test == 2 ) then h = real ( 0.1, kind = 8 ) end if call sphere_quad_00 ( r, polyterm_value_3d, h, result_00 ) call sphere_quad_01 ( r, polyterm_value_3d, h, result_01 ) call sphere_quad_02 ( r, polyterm_value_3d, h, result_02 ) call sphere_quad_03 ( r, polyterm_value_3d, h, result_03 ) write ( *, '(5g14.6)' ) h, result_00, result_01, result_02, result_03 end do end do return end subroutine polyterm_exponent ( action, e ) !*****************************************************************************80 ! !! POLYTERM_EXPONENT gets or sets the exponents for the polynomial term. ! ! Modified: ! ! 08 June 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = 3 ) ACTION. ! 'GET' asks the routine to return the current values in E. ! 'SET' asks the routine to set the current values to E. ! ! Input/output, integer E(3), storage used to set or get values. ! implicit none character ( len = * ) action integer e(3) integer, save, dimension ( 3 ) :: e_save = (/ 0, 0, 0 /) character ( len = 80 ) text character ( len = 80 ) text2 if ( action(1:1) == 'G' ) then e(1:3) = e_save(1:3) else if ( action(1:1) == 'P' ) then write ( *, * ) ' ' if ( all ( e_save(1:3) == 0 ) ) then text = 'P(X,Y,Z) = 1' else text = 'P(X,Y,Z) = ' if ( e_save(1) == 0 ) then else if ( e_save(1) == 1 ) then call s_cat ( text, ' X', text ) else call s_cat ( text, ' X^', text ) write ( text2, '(i2)' ) e_save(1) text2 = adjustl ( text2 ) call s_cat ( text, text2, text ) end if if ( e_save(2) == 0 ) then else if ( e_save(2) == 1 ) then call s_cat ( text, ' Y', text ) else call s_cat ( text, ' Y^', text ) write ( text2, '(i2)' ) e_save(2) text2 = adjustl ( text2 ) call s_cat ( text, text2, text ) end if if ( e_save(3) == 0 ) then else if ( e_save(3) == 1 ) then call s_cat ( text, ' Z', text ) else call s_cat ( text, ' Z^', text ) write ( text2, '(i2)' ) e_save(3) text2 = adjustl ( text2 ) call s_cat ( text, text2, text ) end if end if write ( *, '(a)' ) trim ( text ) else if ( action(1:1) == 'S' ) then e_save(1:3) = e(1:3) end if return end function polyterm_value_3d ( x ) !*****************************************************************************80 ! !! POLYTERM_VALUE_3D evaluates a single polynomial term in 3D. ! ! Discussion: ! ! The polynomial term has the form: ! ! X(1)**E(1) * X(2)**E(2) * X(3)**E(3) ! ! The exponents E(1:3) are set by calling POLYTERM_EXPONENT_SET. ! ! Modified: ! ! 08 June 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X(3), the point where the polynomial term ! is to be evaluated. ! ! Output, real ( kind = 8 ) POLYTERM_VALUE_3D, the value of the ! polynomial term. ! implicit none integer e(3) real factor integer i real ( kind = 8 ) polyterm_value_3d real value real ( kind = 8 ) x(3) call polyterm_exponent ( 'GET', e ) value = 1.0 do i = 1, 3 if ( e(i) == 0 ) then factor = 1.0 else if ( x(i) == 0.0 ) then factor = 0.0 else factor = x(i)**e(i) end if value = value * factor end do polyterm_value_3d = value return end