subroutine filename_inc ( filename ) c*********************************************************************72 c cc FILENAME_INC increments a partially numeric filename. c c Discussion: c c It is assumed that the digits in the name, whether scattered or c connected, represent a number that is to be increased by 1 on c each call. Non-numeric letters of the name are unaffected. c c If the name is empty, then the routine stops. c c If the name contains no digits, the empty string is returned. c c Example: c c Input Output c ----- ------ c a7to11.txt a7to12.txt c a7to99.txt a8to00.txt c a9to99.txt a0to00.txt c cat.txt ' ' c ' ' STOP! c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 08 September 2013 c c Author: c c John Burkardt c c Parameters: c c Input/output, character * ( * ) FILENAME. c On input, a character string to be incremented. c On output, the incremented string. c implicit none character c integer change integer digit character * ( * ) filename integer i integer lens lens = len_trim ( filename ) if ( lens .le. 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILENAME_INC - Fatal error!' write ( *, '(a)' ) ' The input string is empty.' stop 1 end if change = 0 do i = lens, 1, -1 c = filename(i:i) if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then change = change + 1 digit = ichar ( c ) - 48 digit = digit + 1 if ( digit .eq. 10 ) then digit = 0 end if c = char ( digit + 48 ) filename(i:i) = c if ( c .ne. '0' ) then return end if end if end do c c No digits were found. Return blank. c if ( change .eq. 0 ) then filename = ' ' return end if return end subroutine get_unit ( iunit ) c*********************************************************************72 c cc GET_UNIT returns a free FORTRAN unit number. c c Discussion: c c A "free" FORTRAN unit number is a value between 1 and 99 which c is not currently associated with an I/O device. A free FORTRAN unit c number is needed in order to open a file with the OPEN command. c c If IUNIT = 0, then no free FORTRAN unit could be found, although c all 99 units were checked (except for units 5, 6 and 9, which c are commonly reserved for console I/O). c c Otherwise, IUNIT is a value between 1 and 99, representing a c free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 c are special, and will never return those values. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 02 September 2013 c c Author: c c John Burkardt c c Parameters: c c Output, integer IUNIT, the free unit number. c implicit none integer i integer iunit logical value iunit = 0 do i = 1, 99 if ( i .ne. 5 .and. i .ne. 6 .and. i .ne. 9 ) then inquire ( unit = i, opened = value, err = 10 ) if ( .not. value ) then iunit = i return end if end if 10 continue end do return end subroutine padua_order ( l, n ) c*****************************************************************************80 c cc PADUA_ORDER returns the size of the Padua set of given level. c c Discussion: c c The Padua sets are indexed by a level that starts at 0. c This function returns the number of points in each level. c c Example: c c Level Size c ----- ---- c 0 1 c 1 3 c 2 6 c 3 10 c 4 15 c 5 21 c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 06 April 2014 c c Author: c c John Burkardt c c Reference: c c Marco Caliari, Stefano de Marchi, Marco Vianello, c Bivariate interpolation on the square at new nodal sets, c Applied Mathematics and Computation, c Volume 165, Number 2, 2005, pages 261-274. c c Parameters: c c Input, integer L, the level of the set. c 0 <= L c c Output, integer N, the order (number of points) in the set. c implicit none integer i integer l integer n n = 0 do i = 0, l n = n + ( l / 2 ) + 1 if ( mod ( l, 2 ) .eq. 1 .and. mod ( i, 2 ) .eq. 1 ) then n = n + 1 end if end do return end subroutine padua_plot ( l, filename ) c*****************************************************************************80 c cc PADUA_PLOT plots the Padua points of given level. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 06 April 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer L, the level of the set. c 0 <= L c c Input, character * ( * ) FILENAME, a prefix for the files to be created. c implicit none integer l character * 255 command_filename integer command_unit character * 255 data_filename integer data_unit character * ( * ) filename integer j integer n character * ( 255 ) plot_filename double precision xy(2,((l+1)*(l+2))/2) c c Count the points. c n = ( ( l + 1 ) * ( l + 2 ) ) / 2 c c Compute the points. c call padua_points ( l, xy ) c c Create graphics data file. c call get_unit ( data_unit ) data_filename = trim ( filename ) // '_data.txt' open ( unit = data_unit, file = data_filename, & status = 'replace' ) do j = 1, n write ( data_unit, '(2x,g24.16,2x,g24.16)' ) xy(1:2,j) end do close ( unit = data_unit ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Created data file "' // & trim ( data_filename ) // '".' c c Create graphics command file. c call get_unit ( command_unit ) command_filename = trim ( filename ) // '_commands.txt' open ( unit = command_unit, file = command_filename, & status = 'replace' ) write ( command_unit, '(a)' ) '# ' // trim ( command_filename ) write ( command_unit, '(a)' ) '#' write ( command_unit, '(a)' ) '# Usage:' write ( command_unit, '(a)' ) '# gnuplot < ' // & trim ( command_filename ) write ( command_unit, '(a)' ) '#' write ( command_unit, '(a)' ) 'set term png' plot_filename = trim ( filename ) // '.png' write ( command_unit, '(a)' ) 'set output "' & // trim ( plot_filename ) // '"' write ( command_unit, '(a)' ) 'set xlabel "<--- X --->"' write ( command_unit, '(a)' ) 'set ylabel "<--- Y --->"' write ( command_unit, '(a,i2,a)' ) & 'set title "Padua Points, Level ', l, '"' write ( command_unit, '(a)' ) 'set grid' write ( command_unit, '(a)' ) 'set key off' write ( command_unit, '(a)' ) 'set size ratio -1' write ( command_unit, '(a)' ) 'set style data lines' write ( command_unit, '(a)' ) 'set timestamp' write ( command_unit, '(a)' ) 'plot [-1:+1] [-1:+1] "' // & trim ( data_filename ) // & '" using 1:2 with points lt 3 pt 3' close ( unit = command_unit ) write ( *, '(a)' ) & ' Created command file "' // trim ( command_filename ) // '".' return end subroutine padua_points_set ( l, x, y ) c*********************************************************************72 c cc PADUA_POINTS_SET sets the Padua points. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 09 June 2014 c c Author: c c John Burkardt c c Reference: c c Marco Caliari, Stefano de Marchi, Marco Vianello, c Bivariate interpolation on the square at new nodal sets, c Applied Mathematics and Computation, c Volume 165, Number 2, 2005, pages 261-274. c c Parameters: c c Input, integer L, the level. c 0 <= L <= 10. c c Output, double precision X(N), Y(N), the Padua points. c implicit none integer l integer n double precision x(((l+1)*(l+2))/2) double precision y(((l+1)*(l+2))/2) if ( l .eq. 0 ) then x( 1) = 0.000000000000000D+00 y( 1) = 0.000000000000000D+00 else if ( l .eq. 1 ) then x( 1) = -1.000000000000000D+00 y( 1) = -1.000000000000000D+00 x( 2) = -1.000000000000000D+00 y( 2) = 1.000000000000000D+00 x( 3) = 1.000000000000000D+00 y( 3) = 0.000000000000000D+00 else if ( l .eq. 2 ) then x( 1) = -1.000000000000000D+00 y( 1) = -1.000000000000000D+00 x( 2) = -1.000000000000000D+00 y( 2) = 0.5000000000000001D+00 x( 3) = 0.000000000000000D+00 y( 3) = -0.4999999999999998D+00 x( 4) = 0.000000000000000D+00 y( 4) = 1.000000000000000D+00 x( 5) = 1.000000000000000D+00 y( 5) = -1.000000000000000D+00 x( 6) = 1.000000000000000D+00 y( 6) = 0.5000000000000001D+00 else if ( l .eq. 3 ) then x( 1) = -1.000000000000000D+00 y( 1) = -1.000000000000000D+00 x( 2) = -1.000000000000000D+00 y( 2) = 0.000000000000000D+00 x( 3) = -1.000000000000000D+00 y( 3) = 1.000000000000000D+00 x( 4) = -0.4999999999999998D+00 y( 4) = -0.7071067811865475D+00 x( 5) = -0.4999999999999998D+00 y( 5) = 0.7071067811865476D+00 x( 6) = 0.5000000000000001D+00 y( 6) = -1.000000000000000D+00 x( 7) = 0.5000000000000001D+00 y( 7) = 0.000000000000000D+00 x( 8) = 0.5000000000000001D+00 y( 8) = 1.000000000000000D+00 x( 9) = 1.000000000000000D+00 y( 9) = -0.7071067811865475D+00 x(10) = 1.000000000000000D+00 y(10) = 0.7071067811865476D+00 else if ( l .eq. 4 ) then x( 1) = -1.000000000000000D+00 y( 1) = -1.000000000000000D+00 x( 2) = -1.000000000000000D+00 y( 2) = -0.3090169943749473D+00 x( 3) = -1.000000000000000D+00 y( 3) = 0.8090169943749475D+00 x( 4) = -0.7071067811865475D+00 y( 4) = -0.8090169943749473D+00 x( 5) = -0.7071067811865475D+00 y( 5) = 0.3090169943749475D+00 x( 6) = -0.7071067811865475D+00 y( 6) = 1.000000000000000D+00 x( 7) = 0.000000000000000D+00 y( 7) = -1.000000000000000D+00 x( 8) = 0.000000000000000D+00 y( 8) = -0.3090169943749473D+00 x( 9) = 0.000000000000000D+00 y( 9) = 0.8090169943749475D+00 x(10) = 0.7071067811865476D+00 y(10) = -0.8090169943749473D+00 x(11) = 0.7071067811865476D+00 y(11) = 0.3090169943749475D+00 x(12) = 0.7071067811865476D+00 y(12) = 1.000000000000000D+00 x(13) = 1.000000000000000D+00 y(13) = -1.000000000000000D+00 x(14) = 1.000000000000000D+00 y(14) = -0.3090169943749473D+00 x(15) = 1.000000000000000D+00 y(15) = 0.8090169943749475D+00 else if ( l .eq. 5 ) then x( 1) = -1.000000000000000D+00 y( 1) = -1.000000000000000D+00 x( 2) = -1.000000000000000D+00 y( 2) = -0.4999999999999998D+00 x( 3) = -1.000000000000000D+00 y( 3) = 0.5000000000000001D+00 x( 4) = -1.000000000000000D+00 y( 4) = 1.000000000000000D+00 x( 5) = -0.8090169943749473D+00 y( 5) = -0.8660254037844387D+00 x( 6) = -0.8090169943749473D+00 y( 6) = 0.000000000000000D+00 x( 7) = -0.8090169943749473D+00 y( 7) = 0.8660254037844387D+00 x( 8) = -0.3090169943749473D+00 y( 8) = -1.000000000000000D+00 x( 9) = -0.3090169943749473D+00 y( 9) = -0.4999999999999998D+00 x(10) = -0.3090169943749473D+00 y(10) = 0.5000000000000001D+00 x(11) = -0.3090169943749473D+00 y(11) = 1.000000000000000D+00 x(12) = 0.3090169943749475D+00 y(12) = -0.8660254037844387D+00 x(13) = 0.3090169943749475D+00 y(13) = 0.000000000000000D+00 x(14) = 0.3090169943749475D+00 y(14) = 0.8660254037844387D+00 x(15) = 0.8090169943749475D+00 y(15) = -1.000000000000000D+00 x(16) = 0.8090169943749475D+00 y(16) = -0.4999999999999998D+00 x(17) = 0.8090169943749475D+00 y(17) = 0.5000000000000001D+00 x(18) = 0.8090169943749475D+00 y(18) = 1.000000000000000D+00 x(19) = 1.000000000000000D+00 y(19) = -0.8660254037844387D+00 x(20) = 1.000000000000000D+00 y(20) = 0.000000000000000D+00 x(21) = 1.000000000000000D+00 y(21) = 0.8660254037844387D+00 else if ( l .eq. 6 ) then x( 1) = -1.000000000000000D+00 y( 1) = -1.000000000000000D+00 x( 2) = -1.000000000000000D+00 y( 2) = -0.6234898018587335D+00 x( 3) = -1.000000000000000D+00 y( 3) = 0.2225209339563144D+00 x( 4) = -1.000000000000000D+00 y( 4) = 0.9009688679024191D+00 x( 5) = -0.8660254037844387D+00 y( 5) = -0.9009688679024190D+00 x( 6) = -0.8660254037844387D+00 y( 6) = -0.2225209339563143D+00 x( 7) = -0.8660254037844387D+00 y( 7) = 0.6234898018587336D+00 x( 8) = -0.8660254037844387D+00 y( 8) = 1.000000000000000D+00 x( 9) = -0.4999999999999998D+00 y( 9) = -1.000000000000000D+00 x(10) = -0.4999999999999998D+00 y(10) = -0.6234898018587335D+00 x(11) = -0.4999999999999998D+00 y(11) = 0.2225209339563144D+00 x(12) = -0.4999999999999998D+00 y(12) = 0.9009688679024191D+00 x(13) = 0.000000000000000D+00 y(13) = -0.9009688679024190D+00 x(14) = 0.000000000000000D+00 y(14) = -0.2225209339563143D+00 x(15) = 0.000000000000000D+00 y(15) = 0.6234898018587336D+00 x(16) = 0.000000000000000D+00 y(16) = 1.000000000000000D+00 x(17) = 0.5000000000000001D+00 y(17) = -1.000000000000000D+00 x(18) = 0.5000000000000001D+00 y(18) = -0.6234898018587335D+00 x(19) = 0.5000000000000001D+00 y(19) = 0.2225209339563144D+00 x(20) = 0.5000000000000001D+00 y(20) = 0.9009688679024191D+00 x(21) = 0.8660254037844387D+00 y(21) = -0.9009688679024190D+00 x(22) = 0.8660254037844387D+00 y(22) = -0.2225209339563143D+00 x(23) = 0.8660254037844387D+00 y(23) = 0.6234898018587336D+00 x(24) = 0.8660254037844387D+00 y(24) = 1.000000000000000D+00 x(25) = 1.000000000000000D+00 y(25) = -1.000000000000000D+00 x(26) = 1.000000000000000D+00 y(26) = -0.6234898018587335D+00 x(27) = 1.000000000000000D+00 y(27) = 0.2225209339563144D+00 x(28) = 1.000000000000000D+00 y(28) = 0.9009688679024191D+00 else if ( l .eq. 7 ) then x( 1) = -1.000000000000000D+00 y( 1) = -1.000000000000000D+00 x( 2) = -1.000000000000000D+00 y( 2) = -0.7071067811865475D+00 x( 3) = -1.000000000000000D+00 y( 3) = 0.000000000000000D+00 x( 4) = -1.000000000000000D+00 y( 4) = 0.7071067811865476D+00 x( 5) = -1.000000000000000D+00 y( 5) = 1.000000000000000D+00 x( 6) = -0.9009688679024190D+00 y( 6) = -0.9238795325112867D+00 x( 7) = -0.9009688679024190D+00 y( 7) = -0.3826834323650897D+00 x( 8) = -0.9009688679024190D+00 y( 8) = 0.3826834323650898D+00 x( 9) = -0.9009688679024190D+00 y( 9) = 0.9238795325112867D+00 x(10) = -0.6234898018587335D+00 y(10) = -1.000000000000000D+00 x(11) = -0.6234898018587335D+00 y(11) = -0.7071067811865475D+00 x(12) = -0.6234898018587335D+00 y(12) = 0.000000000000000D+00 x(13) = -0.6234898018587335D+00 y(13) = 0.7071067811865476D+00 x(14) = -0.6234898018587335D+00 y(14) = 1.000000000000000D+00 x(15) = -0.2225209339563143D+00 y(15) = -0.9238795325112867D+00 x(16) = -0.2225209339563143D+00 y(16) = -0.3826834323650897D+00 x(17) = -0.2225209339563143D+00 y(17) = 0.3826834323650898D+00 x(18) = -0.2225209339563143D+00 y(18) = 0.9238795325112867D+00 x(19) = 0.2225209339563144D+00 y(19) = -1.000000000000000D+00 x(20) = 0.2225209339563144D+00 y(20) = -0.7071067811865475D+00 x(21) = 0.2225209339563144D+00 y(21) = 0.000000000000000D+00 x(22) = 0.2225209339563144D+00 y(22) = 0.7071067811865476D+00 x(23) = 0.2225209339563144D+00 y(23) = 1.000000000000000D+00 x(24) = 0.6234898018587336D+00 y(24) = -0.9238795325112867D+00 x(25) = 0.6234898018587336D+00 y(25) = -0.3826834323650897D+00 x(26) = 0.6234898018587336D+00 y(26) = 0.3826834323650898D+00 x(27) = 0.6234898018587336D+00 y(27) = 0.9238795325112867D+00 x(28) = 0.9009688679024191D+00 y(28) = -1.000000000000000D+00 x(29) = 0.9009688679024191D+00 y(29) = -0.7071067811865475D+00 x(30) = 0.9009688679024191D+00 y(30) = 0.000000000000000D+00 x(31) = 0.9009688679024191D+00 y(31) = 0.7071067811865476D+00 x(32) = 0.9009688679024191D+00 y(32) = 1.000000000000000D+00 x(33) = 1.000000000000000D+00 y(33) = -0.9238795325112867D+00 x(34) = 1.000000000000000D+00 y(34) = -0.3826834323650897D+00 x(35) = 1.000000000000000D+00 y(35) = 0.3826834323650898D+00 x(36) = 1.000000000000000D+00 y(36) = 0.9238795325112867D+00 else if ( l .eq. 8 ) then x( 1) = -1.000000000000000D+00 y( 1) = -1.000000000000000D+00 x( 2) = -1.000000000000000D+00 y( 2) = -0.7660444431189779D+00 x( 3) = -1.000000000000000D+00 y( 3) = -0.1736481776669303D+00 x( 4) = -1.000000000000000D+00 y( 4) = 0.5000000000000001D+00 x( 5) = -1.000000000000000D+00 y( 5) = 0.9396926207859084D+00 x( 6) = -0.9238795325112867D+00 y( 6) = -0.9396926207859083D+00 x( 7) = -0.9238795325112867D+00 y( 7) = -0.4999999999999998D+00 x( 8) = -0.9238795325112867D+00 y( 8) = 0.1736481776669304D+00 x( 9) = -0.9238795325112867D+00 y( 9) = 0.7660444431189780D+00 x(10) = -0.9238795325112867D+00 y(10) = 1.000000000000000D+00 x(11) = -0.7071067811865475D+00 y(11) = -1.000000000000000D+00 x(12) = -0.7071067811865475D+00 y(12) = -0.7660444431189779D+00 x(13) = -0.7071067811865475D+00 y(13) = -0.1736481776669303D+00 x(14) = -0.7071067811865475D+00 y(14) = 0.5000000000000001D+00 x(15) = -0.7071067811865475D+00 y(15) = 0.9396926207859084D+00 x(16) = -0.3826834323650897D+00 y(16) = -0.9396926207859083D+00 x(17) = -0.3826834323650897D+00 y(17) = -0.4999999999999998D+00 x(18) = -0.3826834323650897D+00 y(18) = 0.1736481776669304D+00 x(19) = -0.3826834323650897D+00 y(19) = 0.7660444431189780D+00 x(20) = -0.3826834323650897D+00 y(20) = 1.000000000000000D+00 x(21) = 0.000000000000000D+00 y(21) = -1.000000000000000D+00 x(22) = 0.000000000000000D+00 y(22) = -0.7660444431189779D+00 x(23) = 0.000000000000000D+00 y(23) = -0.1736481776669303D+00 x(24) = 0.000000000000000D+00 y(24) = 0.5000000000000001D+00 x(25) = 0.000000000000000D+00 y(25) = 0.9396926207859084D+00 x(26) = 0.3826834323650898D+00 y(26) = -0.9396926207859083D+00 x(27) = 0.3826834323650898D+00 y(27) = -0.4999999999999998D+00 x(28) = 0.3826834323650898D+00 y(28) = 0.1736481776669304D+00 x(29) = 0.3826834323650898D+00 y(29) = 0.7660444431189780D+00 x(30) = 0.3826834323650898D+00 y(30) = 1.000000000000000D+00 x(31) = 0.7071067811865476D+00 y(31) = -1.000000000000000D+00 x(32) = 0.7071067811865476D+00 y(32) = -0.7660444431189779D+00 x(33) = 0.7071067811865476D+00 y(33) = -0.1736481776669303D+00 x(34) = 0.7071067811865476D+00 y(34) = 0.5000000000000001D+00 x(35) = 0.7071067811865476D+00 y(35) = 0.9396926207859084D+00 x(36) = 0.9238795325112867D+00 y(36) = -0.9396926207859083D+00 x(37) = 0.9238795325112867D+00 y(37) = -0.4999999999999998D+00 x(38) = 0.9238795325112867D+00 y(38) = 0.1736481776669304D+00 x(39) = 0.9238795325112867D+00 y(39) = 0.7660444431189780D+00 x(40) = 0.9238795325112867D+00 y(40) = 1.000000000000000D+00 x(41) = 1.000000000000000D+00 y(41) = -1.000000000000000D+00 x(42) = 1.000000000000000D+00 y(42) = -0.7660444431189779D+00 x(43) = 1.000000000000000D+00 y(43) = -0.1736481776669303D+00 x(44) = 1.000000000000000D+00 y(44) = 0.5000000000000001D+00 x(45) = 1.000000000000000D+00 y(45) = 0.9396926207859084D+00 else if ( l .eq. 9 ) then x( 1) = -1.000000000000000D+00 y( 1) = -1.000000000000000D+00 x( 2) = -1.000000000000000D+00 y( 2) = -0.8090169943749473D+00 x( 3) = -1.000000000000000D+00 y( 3) = -0.3090169943749473D+00 x( 4) = -1.000000000000000D+00 y( 4) = 0.3090169943749475D+00 x( 5) = -1.000000000000000D+00 y( 5) = 0.8090169943749475D+00 x( 6) = -1.000000000000000D+00 y( 6) = 1.000000000000000D+00 x( 7) = -0.9396926207859083D+00 y( 7) = -0.9510565162951535D+00 x( 8) = -0.9396926207859083D+00 y( 8) = -0.5877852522924730D+00 x( 9) = -0.9396926207859083D+00 y( 9) = 0.000000000000000D+00 x(10) = -0.9396926207859083D+00 y(10) = 0.5877852522924731D+00 x(11) = -0.9396926207859083D+00 y(11) = 0.9510565162951535D+00 x(12) = -0.7660444431189779D+00 y(12) = -1.000000000000000D+00 x(13) = -0.7660444431189779D+00 y(13) = -0.8090169943749473D+00 x(14) = -0.7660444431189779D+00 y(14) = -0.3090169943749473D+00 x(15) = -0.7660444431189779D+00 y(15) = 0.3090169943749475D+00 x(16) = -0.7660444431189779D+00 y(16) = 0.8090169943749475D+00 x(17) = -0.7660444431189779D+00 y(17) = 1.000000000000000D+00 x(18) = -0.4999999999999998D+00 y(18) = -0.9510565162951535D+00 x(19) = -0.4999999999999998D+00 y(19) = -0.5877852522924730D+00 x(20) = -0.4999999999999998D+00 y(20) = 0.000000000000000D+00 x(21) = -0.4999999999999998D+00 y(21) = 0.5877852522924731D+00 x(22) = -0.4999999999999998D+00 y(22) = 0.9510565162951535D+00 x(23) = -0.1736481776669303D+00 y(23) = -1.000000000000000D+00 x(24) = -0.1736481776669303D+00 y(24) = -0.8090169943749473D+00 x(25) = -0.1736481776669303D+00 y(25) = -0.3090169943749473D+00 x(26) = -0.1736481776669303D+00 y(26) = 0.3090169943749475D+00 x(27) = -0.1736481776669303D+00 y(27) = 0.8090169943749475D+00 x(28) = -0.1736481776669303D+00 y(28) = 1.000000000000000D+00 x(29) = 0.1736481776669304D+00 y(29) = -0.9510565162951535D+00 x(30) = 0.1736481776669304D+00 y(30) = -0.5877852522924730D+00 x(31) = 0.1736481776669304D+00 y(31) = 0.000000000000000D+00 x(32) = 0.1736481776669304D+00 y(32) = 0.5877852522924731D+00 x(33) = 0.1736481776669304D+00 y(33) = 0.9510565162951535D+00 x(34) = 0.5000000000000001D+00 y(34) = -1.000000000000000D+00 x(35) = 0.5000000000000001D+00 y(35) = -0.8090169943749473D+00 x(36) = 0.5000000000000001D+00 y(36) = -0.3090169943749473D+00 x(37) = 0.5000000000000001D+00 y(37) = 0.3090169943749475D+00 x(38) = 0.5000000000000001D+00 y(38) = 0.8090169943749475D+00 x(39) = 0.5000000000000001D+00 y(39) = 1.000000000000000D+00 x(40) = 0.7660444431189780D+00 y(40) = -0.9510565162951535D+00 x(41) = 0.7660444431189780D+00 y(41) = -0.5877852522924730D+00 x(42) = 0.7660444431189780D+00 y(42) = 0.000000000000000D+00 x(43) = 0.7660444431189780D+00 y(43) = 0.5877852522924731D+00 x(44) = 0.7660444431189780D+00 y(44) = 0.9510565162951535D+00 x(45) = 0.9396926207859084D+00 y(45) = -1.000000000000000D+00 x(46) = 0.9396926207859084D+00 y(46) = -0.8090169943749473D+00 x(47) = 0.9396926207859084D+00 y(47) = -0.3090169943749473D+00 x(48) = 0.9396926207859084D+00 y(48) = 0.3090169943749475D+00 x(49) = 0.9396926207859084D+00 y(49) = 0.8090169943749475D+00 x(50) = 0.9396926207859084D+00 y(50) = 1.000000000000000D+00 x(51) = 1.000000000000000D+00 y(51) = -0.9510565162951535D+00 x(52) = 1.000000000000000D+00 y(52) = -0.5877852522924730D+00 x(53) = 1.000000000000000D+00 y(53) = 0.000000000000000D+00 x(54) = 1.000000000000000D+00 y(54) = 0.5877852522924731D+00 x(55) = 1.000000000000000D+00 y(55) = 0.9510565162951535D+00 else if ( l .eq. 10 ) then x( 1) = -1.000000000000000D+00 y( 1) = -1.000000000000000D+00 x( 2) = -1.000000000000000D+00 y( 2) = -0.8412535328311811D+00 x( 3) = -1.000000000000000D+00 y( 3) = -0.4154150130018863D+00 x( 4) = -1.000000000000000D+00 y( 4) = 0.1423148382732851D+00 x( 5) = -1.000000000000000D+00 y( 5) = 0.6548607339452851D+00 x( 6) = -1.000000000000000D+00 y( 6) = 0.9594929736144974D+00 x( 7) = -0.9510565162951535D+00 y( 7) = -0.9594929736144974D+00 x( 8) = -0.9510565162951535D+00 y( 8) = -0.6548607339452850D+00 x( 9) = -0.9510565162951535D+00 y( 9) = -0.1423148382732850D+00 x(10) = -0.9510565162951535D+00 y(10) = 0.4154150130018864D+00 x(11) = -0.9510565162951535D+00 y(11) = 0.8412535328311812D+00 x(12) = -0.9510565162951535D+00 y(12) = 1.000000000000000D+00 x(13) = -0.8090169943749473D+00 y(13) = -1.000000000000000D+00 x(14) = -0.8090169943749473D+00 y(14) = -0.8412535328311811D+00 x(15) = -0.8090169943749473D+00 y(15) = -0.4154150130018863D+00 x(16) = -0.8090169943749473D+00 y(16) = 0.1423148382732851D+00 x(17) = -0.8090169943749473D+00 y(17) = 0.6548607339452851D+00 x(18) = -0.8090169943749473D+00 y(18) = 0.9594929736144974D+00 x(19) = -0.5877852522924730D+00 y(19) = -0.9594929736144974D+00 x(20) = -0.5877852522924730D+00 y(20) = -0.6548607339452850D+00 x(21) = -0.5877852522924730D+00 y(21) = -0.1423148382732850D+00 x(22) = -0.5877852522924730D+00 y(22) = 0.4154150130018864D+00 x(23) = -0.5877852522924730D+00 y(23) = 0.8412535328311812D+00 x(24) = -0.5877852522924730D+00 y(24) = 1.000000000000000D+00 x(25) = -0.3090169943749473D+00 y(25) = -1.000000000000000D+00 x(26) = -0.3090169943749473D+00 y(26) = -0.8412535328311811D+00 x(27) = -0.3090169943749473D+00 y(27) = -0.4154150130018863D+00 x(28) = -0.3090169943749473D+00 y(28) = 0.1423148382732851D+00 x(29) = -0.3090169943749473D+00 y(29) = 0.6548607339452851D+00 x(30) = -0.3090169943749473D+00 y(30) = 0.9594929736144974D+00 x(31) = 0.000000000000000D+00 y(31) = -0.9594929736144974D+00 x(32) = 0.000000000000000D+00 y(32) = -0.6548607339452850D+00 x(33) = 0.000000000000000D+00 y(33) = -0.1423148382732850D+00 x(34) = 0.000000000000000D+00 y(34) = 0.4154150130018864D+00 x(35) = 0.000000000000000D+00 y(35) = 0.8412535328311812D+00 x(36) = 0.000000000000000D+00 y(36) = 1.000000000000000D+00 x(37) = 0.3090169943749475D+00 y(37) = -1.000000000000000D+00 x(38) = 0.3090169943749475D+00 y(38) = -0.8412535328311811D+00 x(39) = 0.3090169943749475D+00 y(39) = -0.4154150130018863D+00 x(40) = 0.3090169943749475D+00 y(40) = 0.1423148382732851D+00 x(41) = 0.3090169943749475D+00 y(41) = 0.6548607339452851D+00 x(42) = 0.3090169943749475D+00 y(42) = 0.9594929736144974D+00 x(43) = 0.5877852522924731D+00 y(43) = -0.9594929736144974D+00 x(44) = 0.5877852522924731D+00 y(44) = -0.6548607339452850D+00 x(45) = 0.5877852522924731D+00 y(45) = -0.1423148382732850D+00 x(46) = 0.5877852522924731D+00 y(46) = 0.4154150130018864D+00 x(47) = 0.5877852522924731D+00 y(47) = 0.8412535328311812D+00 x(48) = 0.5877852522924731D+00 y(48) = 1.000000000000000D+00 x(49) = 0.8090169943749475D+00 y(49) = -1.000000000000000D+00 x(50) = 0.8090169943749475D+00 y(50) = -0.8412535328311811D+00 x(51) = 0.8090169943749475D+00 y(51) = -0.4154150130018863D+00 x(52) = 0.8090169943749475D+00 y(52) = 0.1423148382732851D+00 x(53) = 0.8090169943749475D+00 y(53) = 0.6548607339452851D+00 x(54) = 0.8090169943749475D+00 y(54) = 0.9594929736144974D+00 x(55) = 0.9510565162951535D+00 y(55) = -0.9594929736144974D+00 x(56) = 0.9510565162951535D+00 y(56) = -0.6548607339452850D+00 x(57) = 0.9510565162951535D+00 y(57) = -0.1423148382732850D+00 x(58) = 0.9510565162951535D+00 y(58) = 0.4154150130018864D+00 x(59) = 0.9510565162951535D+00 y(59) = 0.8412535328311812D+00 x(60) = 0.9510565162951535D+00 y(60) = 1.000000000000000D+00 x(61) = 1.000000000000000D+00 y(61) = -1.000000000000000D+00 x(62) = 1.000000000000000D+00 y(62) = -0.8412535328311811D+00 x(63) = 1.000000000000000D+00 y(63) = -0.4154150130018863D+00 x(64) = 1.000000000000000D+00 y(64) = 0.1423148382732851D+00 x(65) = 1.000000000000000D+00 y(65) = 0.6548607339452851D+00 x(66) = 1.000000000000000D+00 y(66) = 0.9594929736144974D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PADUA_POINTS_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal value of L = ', l write ( *, '(a)' ) ' Legal values are 0 through 10.' stop 1 end if ! ! Reverse the order to match the published data. ! n = ( ( l + 1 ) * ( l + 2 ) ) / 2 call r8vec_reverse ( n, x ) call r8vec_reverse ( n, y ) return end subroutine padua_points ( l, xy ) c*********************************************************************72 c cc PADUA_POINTS returns the Padua points of order N. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 09 June 2014 c c Author: c c John Burkardt c c Reference: c c Marco Caliari, Stefano de Marchi, Marco Vianello, c Bivariate interpolation on the square at new nodal sets, c Applied Mathematics and Computation, c Volume 165, Number 2, 2005, pages 261-274. c c Parameters: c c Input, integer L, the level of the set. c 0 <= L c c Output, double precision XY(2,((L+1)*(L+2))/2), the Padua points. c implicit none integer l double precision angle1 double precision angle2 integer i integer j integer j_hi integer k double precision r8_pi parameter ( r8_pi = 3.141592653589793D+00 ) double precision xy(2,((l+1)*(l+2))/2) if ( l .eq. 0 ) then xy(1,1) = 0.0D+00 xy(2,1) = 0.0D+00 return end if k = 0 do i = 0, l j_hi = ( l / 2 ) + 1 if ( mod ( l, 2 ) .eq. 1 .and. mod ( i, 2 ) .eq. 1 ) then j_hi = j_hi + 1 end if do j = 1, j_hi k = k + 1 if ( i * 2 .eq. l ) then xy(1,k) = 0.0D+00 else angle1 = dble ( i ) * r8_pi / dble ( l ) xy(1,k) = cos ( angle1 ) end if if ( mod ( i, 2 ) .eq. 0 ) then if ( 2 * ( 2 * j - 1 ) .eq. l + 1 ) then xy(2,k) = 0.0D+00 else angle2 = dble ( 2 * j - 1 ) * r8_pi / dble ( l + 1 ) xy(2,k) = cos ( angle2 ) end if else if ( 2 * ( 2 * j - 2 ) .eq. l + 1 ) then xy(2,k) = 0.0D+00 else angle2 = dble ( 2 * j - 2 ) * r8_pi / dble ( l + 1 ) xy(2,k) = cos ( angle2 ) end if end if end do end do return end subroutine padua_weights_set ( l, w ) c*********************************************************************72 c cc PADUA_WEIGHTS_SET sets quadrature weights for the Padua points. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 09 June 2014 c c Author: c c John Burkardt c c Reference: c c Marco Caliari, Stefano de Marchi, Marco Vianello, c Bivariate interpolation on the square at new nodal sets, c Applied Mathematics and Computation, c Volume 165, Number 2, 2005, pages 261-274. c c Parameters: c c Input, integer L, the level. c 0 <= L <= 10. c c Output, double precision W(N), the quadrature weights. c implicit none integer l integer n double precision w(((l+1)*(l+2))/2) if ( l .eq. 0 ) then w( 1) = 4.000000000000000D+00 else if ( l .eq. 1 ) then w( 1) = 1.000000000000000D+00 w( 2) = 1.000000000000000D+00 w( 3) = 2.000000000000000D+00 else if ( l .eq. 2 ) then w( 1) = 0.0D+00 w( 2) = 0.6666666666666663D+00 w( 3) = 2.222222222222222D+00 w( 4) = 0.4444444444444444D+00 w( 5) = 0.0D+00 w( 6) = 0.6666666666666664D+00 else if ( l .eq. 3 ) then w( 1) = -0.5555555555555480D-01 w( 2) = 0.3333333333333331D+00 w( 3) = -0.5555555555555580D-01 w( 4) = 0.8888888888888886D+00 w( 5) = 0.8888888888888893D+00 w( 6) = 0.2222222222222224D+00 w( 7) = 1.333333333333333D+00 w( 8) = 0.2222222222222220D+00 w( 9) = 0.1111111111111109D+00 w(10) = 0.1111111111111112D+00 else if ( l .eq. 4 ) then w( 1) = -0.8888888888888932D-02 w( 2) = 0.8104919101110961D-01 w( 3) = 0.6117303121111219D-01 w( 4) = 0.3874097078666789D+00 w( 5) = 0.6259236254666545D+00 w( 6) = 0.5333333333333362D-01 w( 7) = 0.7111111111111067D-01 w( 8) = 0.9830822022444241D+00 w( 9) = 0.5458066866444642D+00 w(10) = 0.3874097078666780D+00 w(11) = 0.6259236254666568D+00 w(12) = 0.5333333333333383D-01 w(13) = -0.8888888888888703D-02 w(14) = 0.8104919101110968D-01 w(15) = 0.6117303121111135D-01 else if ( l .eq. 5 ) then w( 1) = -0.1037037037037093D-01 w( 2) = 0.5037037037036911D-01 w( 3) = 0.5037037037037081D-01 w( 4) = -0.1037037037036947D-01 w( 5) = 0.1876963678740801D+00 w( 6) = 0.3460933466518654D+00 w( 7) = 0.1876963678740763D+00 w( 8) = 0.4514390511851724D-01 w( 9) = 0.5541130536814713D+00 w(10) = 0.5541130536814728D+00 w(11) = 0.4514390511851834D-01 w(12) = 0.2804517802740705D+00 w(13) = 0.6376103570518378D+00 w(14) = 0.2804517802740683D+00 w(15) = 0.3189313191851883D-01 w(16) = 0.3288499092814910D+00 w(17) = 0.3288499092814925D+00 w(18) = 0.3189313191851956D-01 w(19) = 0.2074074074074123D-01 w(20) = 0.3851851851851849D-01 w(21) = 0.2074074074074051D-01 else if ( l .eq. 6 ) then w( 1) = -0.3023431594858565D-02 w( 2) = 0.1957267632451884D-01 w( 3) = 0.2633929313290840D-01 w( 4) = 0.1425431928029237D-01 w( 5) = 0.1006383046329639D+00 w( 6) = 0.2208900184526934D+00 w( 7) = 0.1743144584714012D+00 w( 8) = 0.1209372637943976D-01 w( 9) = 0.1934996220710680D-01 w(10) = 0.3245064820875231D+00 w(11) = 0.4027058473592984D+00 w(12) = 0.1677234226317961D+00 w(13) = 0.1953319357827178D+00 w(14) = 0.4489633053035124D+00 w(15) = 0.3721824611057551D+00 w(16) = 0.2479213907785274D-01 w(17) = 0.1934996220710561D-01 w(18) = 0.3245064820875153D+00 w(19) = 0.4027058473592959D+00 w(20) = 0.1677234226317944D+00 w(21) = 0.1006383046329745D+00 w(22) = 0.2208900184526933D+00 w(23) = 0.1743144584714027D+00 w(24) = 0.1209372637944051D-01 w(25) = -0.3023431594861990D-02 w(26) = 0.1957267632451757D-01 w(27) = 0.2633929313290797D-01 w(28) = 0.1425431928029198D-01 else if ( l .eq. 7 ) then w( 1) = -0.3287981859413765D-02 w( 2) = 0.1337868480725671D-01 w( 3) = 0.2063492063491996D-01 w( 4) = 0.1337868480725546D-01 w( 5) = -0.3287981859408898D-02 w( 6) = 0.5949324721885513D-01 w( 7) = 0.1306477599993571D+00 w( 8) = 0.1306477599993581D+00 w( 9) = 0.5949324721885061D-01 w(10) = 0.1263869091685831D-01 w(11) = 0.1979944935601103D+00 w(12) = 0.2832184784823740D+00 w(13) = 0.1979944935601143D+00 w(14) = 0.1263869091685747D-01 w(15) = 0.1221817987389771D+00 w(16) = 0.3150266070593529D+00 w(17) = 0.3150266070593440D+00 w(18) = 0.1221817987389802D+00 w(19) = 0.1771365352315134D-01 w(20) = 0.2490926964598258D+00 w(21) = 0.3408041116306980D+00 w(22) = 0.2490926964598291D+00 w(23) = 0.1771365352314976D-01 w(24) = 0.9646986307476696D-01 w(25) = 0.2557725606433917D+00 w(26) = 0.2557725606433927D+00 w(27) = 0.9646986307476431D-01 w(28) = 0.8649923133686802D-02 w(29) = 0.1062007918394705D+00 w(30) = 0.1505805844901012D+00 w(31) = 0.1062007918394705D+00 w(32) = 0.8649923133690016D-02 w(33) = 0.6355881462931014D-02 w(34) = 0.1405228180237514D-01 w(35) = 0.1405228180237651D-01 w(36) = 0.6355881462928496D-02 else if ( l .eq. 8 ) then w( 1) = -0.1269841269835311D-02 w( 2) = 0.6706089639041270D-02 w( 3) = 0.1111455441352989D-01 w( 4) = 0.1026455026455282D-01 w( 5) = 0.4930678698742625D-02 w( 6) = 0.3633146869162523D-01 w( 7) = 0.8838322767333079D-01 w( 8) = 0.9965911758463214D-01 w( 9) = 0.6400185533755555D-01 w(10) = 0.4061629144893127D-02 w(11) = 0.6772486772485166D-02 w(12) = 0.1258344472781388D+00 w(13) = 0.1927501398511116D+00 w(14) = 0.1699470899470907D+00 w(15) = 0.6342599488133535D-01 w(16) = 0.8376332474107638D-01 w(17) = 0.2170841444607031D+00 w(18) = 0.2477307250801775D+00 w(19) = 0.1648098048612226D+00 w(20) = 0.1004771829779292D-01 w(21) = 0.1015873015872910D-01 w(22) = 0.1784328991205164D+00 w(23) = 0.2729409493576765D+00 w(24) = 0.2364021164021134D+00 w(25) = 0.8936689226256009D-01 w(26) = 0.8376332474107701D-01 w(27) = 0.2170841444607054D+00 w(28) = 0.2477307250801761D+00 w(29) = 0.1648098048612200D+00 w(30) = 0.1004771829779330D-01 w(31) = 0.6772486772485237D-02 w(32) = 0.1258344472781358D+00 w(33) = 0.1927501398511135D+00 w(34) = 0.1699470899470926D+00 w(35) = 0.6342599488133838D-01 w(36) = 0.3633146869162453D-01 w(37) = 0.8838322767332588D-01 w(38) = 0.9965911758463601D-01 w(39) = 0.6400185533755502D-01 w(40) = 0.4061629144888279D-02 w(41) = -0.1269841269836355D-02 w(42) = 0.6706089639046927D-02 w(43) = 0.1111455441352761D-01 w(44) = 0.1026455026454956D-01 w(45) = 0.4930678698747173D-02 else if ( l .eq. 9 ) then w( 1) = -0.1368606701945113D-02 w( 2) = 0.4837977417140975D-02 w( 3) = 0.8876308297144902D-02 w( 4) = 0.8876308297143068D-02 w( 5) = 0.4837977417150492D-02 w( 6) = -0.1368606701935084D-02 w( 7) = 0.2425285860992349D-01 w( 8) = 0.5727330842923516D-01 w( 9) = 0.7008257906578071D-01 w(10) = 0.5727330842922034D-01 w(11) = 0.2425285860989794D-01 w(12) = 0.4659404339099723D-02 w(13) = 0.8354521980498550D-01 w(14) = 0.1370796991940044D+00 w(15) = 0.1370796991940248D+00 w(16) = 0.8354521980500107D-01 w(17) = 0.4659404339109654D-02 w(18) = 0.5564545640233619D-01 w(19) = 0.1524391996823315D+00 w(20) = 0.1877107583774149D+00 w(21) = 0.1524391996823176D+00 w(22) = 0.5564545640232402D-01 w(23) = 0.8186176158691754D-02 w(24) = 0.1295355639606716D+00 w(25) = 0.2061407656847711D+00 w(26) = 0.2061407656847630D+00 w(27) = 0.1295355639606894D+00 w(28) = 0.8186176158692687D-02 w(29) = 0.6234969028097752D-01 w(30) = 0.1730419031522391D+00 w(31) = 0.2169418247419051D+00 w(32) = 0.1730419031522361D+00 w(33) = 0.6234969028097048D-01 w(34) = 0.7506172839505762D-02 w(35) = 0.1142161960569350D+00 w(36) = 0.1802176663769002D+00 w(37) = 0.1802176663769038D+00 w(38) = 0.1142161960569279D+00 w(39) = 0.7506172839512260D-02 w(40) = 0.4031900987631698D-01 w(41) = 0.1142976211857364D+00 w(42) = 0.1413353845521477D+00 w(43) = 0.1142976211857414D+00 w(44) = 0.4031900987631700D-01 w(45) = 0.3239075586856897D-02 w(46) = 0.4317587564913915D-01 w(47) = 0.7015250533601934D-01 w(48) = 0.7015250533601930D-01 w(49) = 0.4317587564913908D-01 w(50) = 0.3239075586852207D-02 w(51) = 0.2550690557469151D-02 w(52) = 0.6084230077461027D-02 w(53) = 0.7421516754852508D-02 w(54) = 0.6084230077458821D-02 w(55) = 0.2550690557473353D-02 else if ( l .eq. 10 ) then w( 1) = -0.6240762604463766D-03 w( 2) = 0.2843227149025789D-02 w( 3) = 0.5250031948150784D-02 w( 4) = 0.5891746241568810D-02 w( 5) = 0.4705736485964679D-02 w( 6) = 0.2135354637732944D-02 w( 7) = 0.1610939653924566D-01 w( 8) = 0.4099595211758227D-01 w( 9) = 0.5326500934654063D-01 w(10) = 0.4863338516658277D-01 w(11) = 0.2843474741781434D-01 w(12) = 0.1719619179693151D-02 w(13) = 0.2883769745121509D-02 w(14) = 0.5724711668876453D-01 w(15) = 0.9659872841640438D-01 w(16) = 0.1053210323353631D+00 w(17) = 0.8066212502628711D-01 w(18) = 0.2855765663647366D-01 w(19) = 0.3981286043310814D-01 w(20) = 0.1090390674981577D+00 w(21) = 0.1430169021081585D+00 w(22) = 0.1313686303763064D+00 w(23) = 0.7932850918298831D-01 w(24) = 0.4610696968783255D-02 w(25) = 0.5086495679684716D-02 w(26) = 0.9311356395361167D-01 w(27) = 0.1562320334111262D+00 w(28) = 0.1696057154254139D+00 w(29) = 0.1283581371975154D+00 w(30) = 0.4603059518094556D-01 w(31) = 0.4894888812994630D-01 w(32) = 0.1347281473526573D+00 w(33) = 0.1764193542601264D+00 w(34) = 0.1635037456303485D+00 w(35) = 0.9822749154565460D-01 w(36) = 0.5704840613923174D-02 w(37) = 0.5086495679679268D-02 w(38) = 0.9311356395362781D-01 w(39) = 0.1562320334111511D+00 w(40) = 0.1696057154253968D+00 w(41) = 0.1283581371975113D+00 w(42) = 0.4603059518094044D-01 w(43) = 0.3981286043311782D-01 w(44) = 0.1090390674981293D+00 w(45) = 0.1430169021081508D+00 w(46) = 0.1313686303763217D+00 w(47) = 0.7932850918299997D-01 w(48) = 0.4610696968790496D-02 w(49) = 0.2883769745110260D-02 w(50) = 0.5724711668875122D-01 w(51) = 0.9659872841642343D-01 w(52) = 0.1053210323353932D+00 w(53) = 0.8066212502626474D-01 w(54) = 0.2855765663644533D-01 w(55) = 0.1610939653928420D-01 w(56) = 0.4099595211758404D-01 w(57) = 0.5326500934649123D-01 w(58) = 0.4863338516656233D-01 w(59) = 0.2843474741784810D-01 w(60) = 0.1719619179720036D-02 w(61) = -0.6240762604606350D-03 w(62) = 0.2843227149011163D-02 w(63) = 0.5250031948172295D-02 w(64) = 0.5891746241587802D-02 w(65) = 0.4705736485965663D-02 w(66) = 0.2135354637703863D-02 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PADUA_WEIGHTS_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal value of L = ', l write ( *, '(a)' ) ' Legal values are 0 through 10.' stop 1 end if ! ! Reverse the order to match the published data. ! n = ( ( l + 1 ) * ( l + 2 ) ) / 2 call r8vec_reverse ( n, w ) return end subroutine padua_weights ( l, w ) c*********************************************************************72 c cc PADUA_WEIGHTS returns quadrature weights do Padua points. c c Discussion: c c The order of the weights corresponds to the ordering used c by the companion function padua_points(). c c Caliari, de Marchi and Vianello supplied a MATLAB code pdwtsMM c which carries out this same computation in a way that makes c more efficient use of MATLAB's vector and matrix capabilities. c This version of the computation was painfully rewritten to display c the individual scalar computations, so that it could be translated c into other languages. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 11 June 2014 c c Author: c c John Burkardt c c Reference: c c Marco Caliari, Stefano de Marchi, Marco Vianello, c Bivariate interpolation on the square at new nodal sets, c Applied Mathematics and Computation, c Volume 165, Number 2, 2005, pages 261-274. c c Parameters: c c Input, integer L, the level of the set. c 0 <= L c c Output, double precision W((L+1)*(L+2)/2), the quadrature weights. c implicit none integer l double precision angle integer i integer i2 integer j integer j2 integer lp1h integer lp2h integer lp3h double precision mi double precision mj double precision mom((l+2)/2,(l+2)/2) double precision r8_pi parameter ( r8_pi = 3.141592653589793D+00 ) double precision te1((l+2)/2,(l+2)/2) double precision te2((l+2)/2,(l+3)/2) double precision tmteo((l+3)/2,(l+1)/2) double precision tmtoe((l+2)/2,(l+2)/2) double precision to1((l+2)/2,(l+1)/2) double precision to2((l+2)/2,(l+2)/2) double precision w(((l+1)*(l+2))/2) double precision w1((l+2)/2,(l+2)/2) double precision w2((l+3)/2,(l+1)/2) if ( l .eq. 0 ) then w(1) = 4.0D+00 return end if c c Relatives of L/2: c lp1h = ( l + 1 ) / 2 lp2h = ( l + 2 ) / 2 lp3h = ( l + 3 ) / 2 c c TE1, TE2, TO1, TO2: c Even and odd Chebyshev polynomials on subgrids 1 and 2. c do j = 1, lp2h do i = 1, lp2h angle = r8_pi * dble ( 2 * ( i - 1 ) * ( 2 * j - 2 ) ) & / dble ( l ) te1(i,j) = cos ( angle ) end do end do do j = 1, lp2h do i = 2, lp2h te1(i,j) = te1(i,j) * sqrt ( 2.0D+00 ) end do end do do j = 1, lp1h do i = 1, lp2h angle = r8_pi * dble ( 2 * ( i - 1 ) * ( 2 * j - 1 ) ) & / dble ( l ) to1(i,j) = cos ( angle ) end do end do do j = 1, lp1h do i = 2, lp2h to1(i,j) = to1(i,j) * sqrt ( 2.0D+00 ) end do end do do j = 1, lp3h do i = 1, lp2h angle = r8_pi * dble ( 2 * ( i - 1 ) * ( 2 * j - 2 ) ) & / dble ( l + 1 ) te2(i,j) = cos ( angle ) end do end do do j = 1, lp3h do i = 2, lp2h te2(i,j) = te2(i,j) * sqrt ( 2.0D+00 ) end do end do do j = 1, lp2h do i = 1, lp2h angle = r8_pi * dble ( 2 * ( i - 1 ) * ( 2 * j - 1 ) ) & / dble ( l + 1 ) to2(i,j) = cos ( angle ) end do end do do j = 1, lp2h do i = 2, lp2h to2(i,j) = to2(i,j) * sqrt ( 2.0D+00 ) end do end do c c MOM: Moments matrix do even * even pairs. c do j = 1, lp2h mj = 2.0D+00 * sqrt ( 2.0D+00 ) & / dble ( 1 - ( 2 * j - 2 ) ** 2 ) do i = 1, lp2h + 1 - j mi = 2.0D+00 * sqrt ( 2.0D+00 ) & / dble ( 1 - ( 2 * i - 2 ) ** 2 ) mom(i,j) = mi * mj end do end do do j = 1, lp2h mom(1,j) = mom(1,j) / sqrt ( 2.0D+00 ) end do do i = 1, lp2h mom(i,1) = mom(i,1) / sqrt ( 2.0D+00 ) end do if ( mod ( l, 2 ) .eq. 0 ) then mom(lp2h,1) = mom(lp2h,1) / 2.0D+00 end if c c TMTOE and TMTEO: matrix products. c do j = 1, lp2h do i = 1, lp2h tmtoe(i,j) = 0.0D+00 end do end do do j2 = 1, lp2h do i2 = 1, lp2h + 1 - j2 do j = 1, lp2h do i = 1, lp2h tmtoe(i,j) = tmtoe(i,j) & + to2(i2,i) * mom(j2,i2) * te1(j2,j) end do end do end do end do do j = 1, lp1h do i = 1, lp3h tmteo(i,j) = 0.0D+00 end do end do do j2 = 1, lp2h do i2 = 1, lp2h + 1 - j2 do j = 1, lp1h do i = 1, lp3h tmteo(i,j) = tmteo(i,j) & + te2(i2,i) * mom(j2,i2) * to1(j2,j) end do end do end do end do c c W1 and W2: Interpolation weight matrices. c do j = 1, lp2h do i = 1, lp2h w1(i,j) = 2.0D+00 / dble ( l * ( l + 1 ) ) end do end do do i = 1, lp2h w1(i,1) = w1(i,1) / 2.0D+00 end do if ( mod ( l, 2 ) .eq. 0 ) then do i = 1, lp2h w1(i,lp2h) = w1(i,lp2h) / 2.0D+00 end do do j = 1, lp2h w1(lp2h,j) = w1(lp2h,j) / 2.0D+00 end do end if do j = 1, lp1h do i = 1, lp3h w2(i,j) = 2.0D+00 / dble ( l * ( l + 1 ) ) end do end do do j = 1, lp1h w2(1,j) = w2(1,j) / 2.0D+00 end do if ( mod ( l, 2 ) .eq. 1 ) then do j = 1, lp1h w2(lp3h,j) = w2(lp3h,j) / 2.0D+00 end do do i = 1, lp3h w2(i,lp1h) = w2(i,lp1h) / 2.0D+00 end do end if c c Cubature weights as matrices on the subgrids. c do j = 1, lp2h do i = 1, lp2h w1(i,j) = w1(i,j) * tmtoe(i,j) end do end do do j = 1, lp1h do i = 1, lp3h w2(i,j) = w2(i,j) * tmteo(i,j) end do end do c c Pack weight matrices W1 and W2 into the vector W. c if ( mod ( l, 2 ) .eq. 0 ) then do j = 1, lp2h do i = 1, lp2h w(i+(2*j-2)*lp2h) = w1(i,j) end do end do do j = 1, lp1h do i = 1, lp3h w(i+(2*j-1)*lp2h) = w2(i,j) end do end do else do j = 1, lp1h do i = 1, lp2h w(i+(j-1)*(l+2)) = w1(i,j) end do end do do j = 1, lp1h do i = 1, lp3h w(i+lp2h+(j-1)*(l+2)) = w2(i,j) end do end do end if return end subroutine r8mat_transpose_print ( m, n, a, title ) c*********************************************************************72 c cc R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. c c Discussion: c c An R8MAT is an array of R8's. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 28 April 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer M, N, the number of rows and columns. c c Input, double precision A(M,N), an M by N matrix to be printed. c c Input, character*(*) TITLE, a title. c implicit none integer m integer n double precision a(m,n) character*(*) 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 ) c*********************************************************************72 c cc R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT transposed. c c Discussion: c c An R8MAT is an array of R8's. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 28 April 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer M, N, the number of rows and columns. c c Input, double precision A(M,N), an M by N matrix to be printed. c c Input, integer ILO, JLO, the first row and column to print. c c Input, integer IHI, JHI, the last row and column to print. c c Input, character * ( * ) TITLE, a title. c implicit none integer incx parameter ( incx = 5 ) integer m integer n double precision a(m,n) character * ( 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 * ( * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) if ( m .le. 0 .or. n .le. 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' 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 ( *, '(2x,i8,a,5a14)' ) j, ':', ( ctemp(i), i = 1, inc ) end do end do return end subroutine r8vec_reverse ( n, a ) c*********************************************************************72 c cc R8VEC_REVERSE reverses the elements of an R8VEC. c c Discussion: c c An R8VEC is a vector of R8's. c c Example: c c Input: c c N = 5, c A = ( 11.0, 12.0, 13.0, 14.0, 15.0 ). c c Output: c c A = ( 15.0, 14.0, 13.0, 12.0, 11.0 ). c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 04 July 2009 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of entries in the array. c c Input/output, double precision A(N), the array to be reversed. c implicit none integer n double precision a(n) integer i integer i_hi double precision t i_hi = n / 2 do i = 1, i_hi t = a(i) a(i) = a(n+1-i) a(n+1-i) = t end do 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 GNU LGPL 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