program main !*****************************************************************************80 ! !! MAIN is the main program for DELAUNAY_LMAP_2D. ! ! Discussion: ! ! DELAUNAY_LMAP_2D computes a Delaunay triangulation under a linear map. ! ! The dataset is simply a set of points in the plane. ! ! The linear map is applied using a matrix, which we denote by A. ! A embodies a linear transformation applied to the data. ! ! Thus, given a set of points V1, V2, ..., VN, we will actually ! be applying a standard Delaunay triangulation to A*V1, A*V2, ..., A*VN. ! ! In turn, the Delaunay triangulation is an organization of the ! data into triples, forming a triangulation of the data, with ! the property that the circumcircle of each triangle never contains ! another data point. ! ! For the linear mapped data, we can interpret this process in one ! of two ways: ! ! * Given V1, V2, ..., VN, we compute A*V1, A*V2, ..., A*VN, and ! carry out a standard Delaunay triangulation on that data, so that ! the linearly mapped Delaunay triangulation has the property that, ! if Va, Vb, and Vc form a linearly mapped Delaunay triangle, the ! circumcircle of (A*Va, A*Vb, A*Vc) does not contain the image ! A*Vd of any other data point Vd. ! ! * Given V1, V2, ..., VN, we compute a linear mapped Delaunay triangulation ! which has the property that, for any three points (Va,Vb,Vc) ! that form a linearly mapped Delaunay triangle, the "circumellipse" which ! passes through (Va, Vb, Vc) does not contain any other data point Vd. ! ! Thanks to Yoshimitsu Asahi for pointing out several errors in a previous ! version of this program, 08 June 2010. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 08 June 2010 ! ! Author: ! ! John Burkardt ! ! Usage: ! ! delaunay_lmap_2d node_file_name matrix_file_name ! implicit none integer ( kind = 4 ) arg_num character ( len = 255 ) :: eps_file_name = ' ' integer ( kind = 4 ) iarg integer ( kind = 4 ) iargc integer ( kind = 4 ) :: ierror = 0 character ( len = 255 ) :: node_file_name = ' ' integer ( kind = 4 ) m1 real ( kind = 8 ) matrix(2,2) character ( len = 255 ) :: matrix_file_name = ' ' integer ( kind = 4 ) node_num real ( kind = 8 ), allocatable, dimension ( :, : ) :: table integer ( kind = 4 ) triangle character ( len = 255 ) :: triangle_file_name = ' ' logical, allocatable, dimension ( : ) :: triangle_mask integer ( kind = 4 ), allocatable, dimension ( :, : ) :: triangle_neighbor integer ( kind = 4 ), allocatable, dimension ( :, : ) :: triangle_node integer ( kind = 4 ) triangle_num integer ( kind = 4 ), parameter :: triangle_order = 3 call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DELAUNAY_LMAP_2D' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Read a node file of N points in 2 dimensions,' write ( *, '(a)' ) ' Read a file defining a 2 by 2 matrix A, where' write ( *, '(a)' ) ' dot_product ( X, Y ) = (A*X)'' * (A*Y)' write ( *, '(a)' ) ' so that the norm is defined as' write ( *, '(a)' ) ' norm ( X ) = sqrt ( dot_product ( X, X ) )' write ( *, '(a)' ) ' so the distance from X to Y is defined by' write ( *, '(a)' ) ' dist ( X, Y ) = norm ( X - Y ).' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Under this modified distance, compute the ' write ( *, '(a)' ) ' corresponding Delaunay triangulation.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Write an file defining the triangulation of the nodes.' write ( *, '(a)' ) ' Write an EPS image of the triangulation.' ! ! Get the number of command line arguments. ! arg_num = iargc ( ) ! ! If at least one command line argument, it's the node file name. ! if ( 1 <= arg_num ) then iarg = 1 call getarg ( iarg, node_file_name ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DELAUNAY_LMAP_2D:' write ( *, '(a)' ) ' Please enter the name of the node file name.' read ( *, '(a)' ) node_file_name end if ! ! If at least two command line argument, it's the name of the matrix file. ! if ( 2 <= arg_num ) then iarg = 2 call getarg ( iarg, matrix_file_name ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DELAUNAY_LMAP_2D:' write ( *, '(a)' ) ' Please enter the name of the input matrix file.' read ( *, '(a)' ) matrix_file_name end if ! ! Create the output file name from the input file name. ! triangle_file_name = node_file_name call file_name_ext_swap ( triangle_file_name, 'delaunay.txt' ) ! ! Read the point coordinates. ! call r8mat_header_read ( node_file_name, m1, node_num ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Read the header of "' // trim ( node_file_name ) //'".' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Spatial dimension M = ', m1 write ( *, '(a,i8)' ) ' Number of points N = ', node_num allocate ( table(1:m1,1:node_num) ) call r8mat_data_read ( node_file_name, m1, node_num, table ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Read the data in "' // trim ( node_file_name ) //'".' call r8mat_transpose_print_some ( m1, node_num, table, 1, 1, m1, 5, & ' The first 5 nodes:' ) ! ! Read the matrix. ! call r8mat_data_read ( matrix_file_name, 2, 2, matrix ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Read the data in "' // trim ( matrix_file_name ) //'".' call r8mat_print ( 2, 2, matrix, ' Linear map matrix A:' ) ! ! Determine the Delaunay triangulation under the linear map. ! allocate ( triangle_node(triangle_order,3*node_num) ) allocate ( triangle_neighbor(triangle_order,3*node_num) ) call dtris2_lmap ( node_num, table, matrix, triangle_num, triangle_node, & triangle_neighbor ) ! ! Print a portion of the triangulation. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Computed the triangulation.' call i4mat_transpose_print_some ( triangle_order, triangle_num, & triangle_node, 1, 1, triangle_order, 5, ' The first 5 triangles:' ) ! ! Write the triangulation to a file. ! call i4mat_write ( triangle_file_name, triangle_order, triangle_num, & triangle_node ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Wrote the triangulation data to "' & // trim ( triangle_file_name ) //'".' ! ! Plot the triangulation. ! TRIANGLE_MASK selects the triangles we want to see. ! eps_file_name = node_file_name call file_name_ext_swap ( eps_file_name, 'delaunay.eps' ) allocate ( triangle_mask(1:triangle_num) ) triangle_mask(1:triangle_num) = .true. call triangulation_plot_eps ( eps_file_name, node_num, table(1,1:node_num), & table(2,1:node_num), triangle_num, triangle_mask, triangle_node, & ' Linear Map Delaunay triangulation' ) ! ! Free memory. ! deallocate ( table ) deallocate ( triangle_mask ) deallocate ( triangle_neighbor ) deallocate ( triangle_node ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DELAUNAY_LMAP_2D' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine ch_cap ( c ) !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none character c integer ( kind = 4 ) itemp itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end function ch_eqi ( c1, c2 ) !*****************************************************************************80 ! !! CH_EQI is a case insensitive comparison of two characters for equality. ! ! Example: ! ! CH_EQI ( 'A', 'a' ) is .TRUE. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 28 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C1, C2, the characters to compare. ! ! Output, logical CH_EQI, the result of the comparison. ! implicit none logical ch_eqi character c1 character c1_cap character c2 character c2_cap c1_cap = c1 c2_cap = c2 call ch_cap ( c1_cap ) call ch_cap ( c2_cap ) if ( c1_cap == c2_cap ) then ch_eqi = .true. else ch_eqi = .false. end if return end subroutine ch_to_digit ( c, digit ) !*****************************************************************************80 ! !! CH_TO_DIGIT returns the value of a base 10 digit. ! ! Example: ! ! C DIGIT ! --- ----- ! '0' 0 ! '1' 1 ! ... ... ! '9' 9 ! ' ' 0 ! 'X' -1 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, the decimal digit, '0' through '9' or blank ! are legal. ! ! Output, integer ( kind = 4 ) DIGIT, the corresponding value. ! If C was 'illegal', then DIGIT is -1. ! implicit none character c integer ( kind = 4 ) digit if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then digit = ichar ( c ) - 48 else if ( c == ' ' ) then digit = 0 else digit = - 1 end if return end function diaedg_lmap ( x0, y0, x1, y1, x2, y2, x3, y3, matrix ) !*****************************************************************************80 ! !! DIAEDG_LMAP chooses a diagonal edge under a linear map. ! ! Discussion: ! ! The routine determines whether 0--2 or 1--3 is the diagonal edge ! that should be chosen, based on the circumcircle criterion, where ! (X0,Y0), (X1,Y1), (X2,Y2), (X3,Y3) are the vertices of a simple ! quadrilateral in counterclockwise order. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 April 2004 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) X0, Y0, X1, Y1, X2, Y2, X3, Y3, the ! coordinates of the vertices of a quadrilateral, given in ! counter clockwise order. ! ! Input, real ( kind = 8 ) MATRIX(2,2), a linear map to be implicitly applied ! to the points. ! ! Output, integer ( kind = 4 ) DIAEDG_LMAP, chooses a diagonal: ! +1, if diagonal edge 02 is chosen; ! -1, if diagonal edge 13 is chosen; ! 0, if the four vertices are cocircular. ! implicit none real ( kind = 8 ) ca real ( kind = 8 ) cb integer ( kind = 4 ) diaedg_lmap real ( kind = 8 ) dx10 real ( kind = 8 ) dx12 real ( kind = 8 ) dx30 real ( kind = 8 ) dx32 real ( kind = 8 ) dy10 real ( kind = 8 ) dy12 real ( kind = 8 ) dy30 real ( kind = 8 ) dy32 real ( kind = 8 ) matrix(2,2) real ( kind = 8 ) s real ( kind = 8 ) tol real ( kind = 8 ) tola real ( kind = 8 ) tolb real ( kind = 8 ) x0 real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) y0 real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 tol = 100.0D+00 * epsilon ( tol ) dx10 = x1 - x0 dy10 = y1 - y0 dx12 = x1 - x2 dy12 = y1 - y2 dx30 = x3 - x0 dy30 = y3 - y0 dx32 = x3 - x2 dy32 = y3 - y2 call transform_lmap ( matrix, dx10, dy10 ) call transform_lmap ( matrix, dx12, dy12 ) call transform_lmap ( matrix, dx30, dy30 ) call transform_lmap ( matrix, dx32, dy32 ) tola = tol * max ( abs ( dx10 ), abs ( dy10 ), abs ( dx30 ), abs ( dy30 ) ) tolb = tol * max ( abs ( dx12 ), abs ( dy12 ), abs ( dx32 ), abs ( dy32 ) ) ca = dx10 * dx30 + dy10 * dy30 cb = dx12 * dx32 + dy12 * dy32 if ( tola < ca .and. tolb < cb ) then diaedg_lmap = -1 else if ( ca < -tola .and. cb < -tolb ) then diaedg_lmap = 1 else tola = max ( tola, tolb ) s = ( dx10 * dy30 - dx30 * dy10 ) * cb + ( dx32 * dy12 - dx12 * dy32 ) * ca if ( tola < s ) then diaedg_lmap = -1 else if ( s < -tola ) then diaedg_lmap = 1 else diaedg_lmap = 0 end if end if return end subroutine dtris2_lmap ( point_num, point_xy, matrix, tri_num, tri_vert, & tri_nabe ) !*****************************************************************************80 ! !! DTRIS2_LMAP constructs a Delaunay triangulation of 2D linear mapped vertices. ! ! Discussion: ! ! The routine constructs the Delaunay triangulation of a set of 2D vertices ! using an incremental approach and diagonal edge swaps. Vertices are ! first sorted in lexicographically increasing (X,Y) order, and ! then are inserted one at a time from outside the convex hull. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 April 2004 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, integer ( kind = 4 ) POINT_NUM, the number of vertices. ! ! Input, real ( kind = 8 ) POINT_XY(2,POINT_NUM), the coordinates ! of the vertices. ! ! Input, real ( kind = 8 ) MATRIX(2,2), the linear map which should ! be implicitly applied to the points before the triangulation is done. ! ! Output, integer ( kind = 4 ) TRI_NUM, the number of triangles in the ! triangulation; TRI_NUM is equal to 2*POINT_NUM - NB - 2, where NB is the ! number of boundary vertices. ! ! Output, integer ( kind = 4 ) TRI_VERT(3,TRI_NUM), the nodes that make up ! each triangle. The elements are indices of POINT_XY. The vertices of ! the triangles are in counter clockwise order. ! ! Output, integer ( kind = 4 ) TRI_NABE(3,TRI_NUM), the triangle neighbor ! list. Positive elements are indices of TIL; negative elements are used ! for links of a counter clockwise linked list of boundary edges; ! LINK = -(3*I + J-1) where I, J = triangle, edge index; TRI_NABE(J,I) ! refers to the neighbor along edge from vertex J to J+1 (mod 3). ! implicit none integer ( kind = 4 ) point_num real ( kind = 8 ) cmax integer ( kind = 4 ) e integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) indx(point_num) integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) ledg integer ( kind = 4 ) lr integer ( kind = 4 ) lrline integer ( kind = 4 ) ltri integer ( kind = 4 ) m integer ( kind = 4 ) m1 integer ( kind = 4 ) m2 real ( kind = 8 ) matrix(2,2) integer ( kind = 4 ) n real ( kind = 8 ) point_xy(2,point_num) integer ( kind = 4 ) redg integer ( kind = 4 ) rtri integer ( kind = 4 ) stack(point_num) integer ( kind = 4 ) t real ( kind = 8 ) tol integer ( kind = 4 ) top integer ( kind = 4 ) tri_nabe(3,point_num*2) integer ( kind = 4 ) tri_num integer ( kind = 4 ) tri_vert(3,point_num*2) tol = 100.0D+00 * epsilon ( tol ) ierr = 0 ! ! Sort the vertices by increasing (x,y). ! call r82vec_sort_heap_index_a ( point_num, point_xy, indx ) call r82vec_permute ( point_num, point_xy, indx ) ! ! Make sure that the data points are "reasonably" distinct. ! m1 = 1 do i = 2, point_num m = m1 m1 = i k = 0 do j = 1, 2 cmax = max ( abs ( point_xy(j,m) ), abs ( point_xy(j,m1) ) ) if ( tol * ( cmax + 1.0D+00 ) & < abs ( point_xy(j,m) - point_xy(j,m1) ) ) then k = j exit end if end do if ( k == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2_LMAP - Fatal error!' write ( *, '(a,i8)' ) ' Fails for point number I = ', i write ( *, '(a,i8)' ) ' M = ', m write ( *, '(a,i8)' ) ' M1 = ', m1 write ( *, '(a,2g14.6)' ) ' X,Y(M) = ', point_xy(1,m), point_xy(2,m) write ( *, '(a,2g14.6)' ) ' X,Y(M1) = ', point_xy(1,m1), point_xy(2,m1) ierr = 224 stop end if end do ! ! Starting from points M1 and M2, search for a third point M that ! makes a "healthy" triangle (M1,M2,M) ! m1 = 1 m2 = 2 j = 3 do if ( point_num < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2_LMAP - Fatal error!' ierr = 225 stop end if m = j lr = lrline ( point_xy(1,m), point_xy(2,m), point_xy(1,m1), & point_xy(2,m1), point_xy(1,m2), point_xy(2,m2), 0.0D+00 ) if ( lr /= 0 ) then exit end if j = j + 1 end do ! ! Set up the triangle information for (M1,M2,M), and for any other ! triangles you created because points were collinear with M1, M2. ! tri_num = j - 2 if ( lr == -1 ) then tri_vert(1,1) = m1 tri_vert(2,1) = m2 tri_vert(3,1) = m tri_nabe(3,1) = -3 do i = 2, tri_num m1 = m2 m2 = i+1 tri_vert(1,i) = m1 tri_vert(2,i) = m2 tri_vert(3,i) = m tri_nabe(1,i-1) = -3 * i tri_nabe(2,i-1) = i tri_nabe(3,i) = i - 1 end do tri_nabe(1,tri_num) = -3 * tri_num - 1 tri_nabe(2,tri_num) = -5 ledg = 2 ltri = tri_num else tri_vert(1,1) = m2 tri_vert(2,1) = m1 tri_vert(3,1) = m tri_nabe(1,1) = -4 do i = 2, tri_num m1 = m2 m2 = i+1 tri_vert(1,i) = m2 tri_vert(2,i) = m1 tri_vert(3,i) = m tri_nabe(3,i-1) = i tri_nabe(1,i) = -3 * i - 3 tri_nabe(2,i) = i - 1 end do tri_nabe(3,tri_num) = -3 * tri_num tri_nabe(2,1) = -3 * tri_num - 2 ledg = 2 ltri = 1 end if ! ! Insert the vertices one at a time from outside the convex hull, ! determine visible boundary edges, and apply diagonal edge swaps until ! Delaunay triangulation of vertices (so far) is obtained. ! top = 0 do i = j+1, point_num m = i m1 = tri_vert(ledg,ltri) if ( ledg <= 2 ) then m2 = tri_vert(ledg+1,ltri) else m2 = tri_vert(1,ltri) end if lr = lrline ( point_xy(1,m), point_xy(2,m), point_xy(1,m1), & point_xy(2,m1), point_xy(1,m2), point_xy(2,m2), 0.0D+00 ) if ( 0 < lr ) then rtri = ltri redg = ledg ltri = 0 else l = -tri_nabe(ledg,ltri) rtri = l / 3 redg = mod(l,3) + 1 end if call vbedg ( point_xy(1,m), point_xy(2,m), point_num, point_xy, tri_num, & tri_vert, tri_nabe, ltri, ledg, rtri, redg ) n = tri_num + 1 l = -tri_nabe(ledg,ltri) do t = l / 3 e = mod ( l, 3 ) + 1 l = -tri_nabe(e,t) m2 = tri_vert(e,t) if ( e <= 2 ) then m1 = tri_vert(e+1,t) else m1 = tri_vert(1,t) end if tri_num = tri_num + 1 tri_nabe(e,t) = tri_num tri_vert(1,tri_num) = m1 tri_vert(2,tri_num) = m2 tri_vert(3,tri_num) = m tri_nabe(1,tri_num) = t tri_nabe(2,tri_num) = tri_num - 1 tri_nabe(3,tri_num) = tri_num + 1 top = top + 1 if ( point_num < top ) then ierr = 8 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2_LMAP - Fatal error!' write ( *, '(a)' ) ' Stack overflow.' stop end if stack(top) = tri_num if ( t == rtri .and. e == redg ) then exit end if end do tri_nabe(ledg,ltri) = -3 * n - 1 tri_nabe(2,n) = -3 * tri_num - 2 tri_nabe(3,tri_num) = -l ltri = n ledg = 2 call swapec_lmap ( m, matrix, top, ltri, ledg, point_num, point_xy, & tri_num, tri_vert, tri_nabe, stack, ierr ) if ( ierr /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2_LMAP - Fatal error!' write ( *, '(a)' ) ' Error return from SWAPEC_LMAP.' stop end if end do ! ! Now account for the sorting that we did. ! do i = 1, 3 do j = 1, tri_num tri_vert(i,j) = indx ( tri_vert(i,j) ) end do end do call perm_inv ( point_num, indx ) call r82vec_permute ( point_num, point_xy, indx ) return end subroutine file_column_count ( input_filename, column_num ) !*****************************************************************************80 ! !! FILE_COLUMN_COUNT counts the number of columns in the first line of a file. ! ! Discussion: ! ! The file is assumed to be a simple text file. ! ! Most lines of the file is presumed to consist of COLUMN_NUM words, ! separated by spaces. There may also be some blank lines, and some ! comment lines, ! which have a "#" in column 1. ! ! The routine tries to find the first non-comment non-blank line and ! counts the number of words in that line. ! ! If all lines are blanks or comments, it goes back and tries to analyze ! a comment line. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 21 June 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the file. ! ! Output, integer ( kind = 4 ) COLUMN_NUM, the number of columns in the file. ! implicit none integer ( kind = 4 ) column_num logical got_one character ( len = * ) input_filename integer ( kind = 4 ) input_status integer ( kind = 4 ) input_unit character ( len = 255 ) line ! ! Open the file. ! call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & form = 'formatted', access = 'sequential', iostat = input_status ) if ( input_status /= 0 ) then column_num = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' & // trim ( input_filename ) // '" on unit ', input_unit stop end if ! ! Read one line, but skip blank lines and comment lines. ! got_one = .false. do read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if if ( line(1:1) == '#' ) then cycle end if got_one = .true. exit end do if ( .not. got_one ) then rewind ( input_unit ) do read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if got_one = .true. exit end do end if close ( unit = input_unit ) if ( .not. got_one ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Warning!' write ( *, '(a)' ) ' The file does not seem to contain any data.' column_num = -1 return end if call s_word_count ( line, column_num ) return end subroutine file_name_ext_get ( file_name, i, j ) !*****************************************************************************80 ! !! FILE_NAME_EXT_GET determines the "extension" of a file name. ! ! Discussion: ! ! The "extension" of a filename is the string of characters ! that appears after the LAST period in the name. A file ! with no period, or with a period as the last character ! in the name, has a "null" extension. ! ! Blanks are unusual in filenames. This routine ignores all ! trailing blanks, but will treat initial or internal blanks ! as regular characters acceptable in a file name. ! ! Example: ! ! FILE_NAME I J ! ! bob.for 4 7 ! N.B.C.D 6 7 ! Naomi. 6 6 ! Arthur 0 0 ! .com 1 1 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 17 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, a file name to be examined. ! ! Output, integer ( kind = 4 ) I, J, the indices of the first and last ! characters in the file extension. ! If no period occurs in FILE_NAME, then ! I = J = 0; ! Otherwise, ! I is the position of the LAST period in FILE_NAME, and J is the ! position of the last nonblank character following the period. ! implicit none character ( len = * ) file_name integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) s_index_last i = s_index_last ( file_name, '.' ) if ( i /= 0 ) then j = len_trim ( file_name ) else j = 0 end if return end subroutine file_name_ext_swap ( file_name, ext ) !*****************************************************************************80 ! !! FILE_NAME_EXT_SWAP replaces the current "extension" of a file name. ! ! Discussion: ! ! The "extension" of a filename is the string of characters ! that appears after the LAST period in the name. A file ! with no period, or with a period as the last character ! in the name, has a "null" extension. ! ! Example: ! ! Input Output ! ================ ========= ! FILE_NAME EXT FILE_NAME ! ! bob.for obj bob.obj ! bob.bob.bob txt bob.bob.txt ! bob yak bob.yak ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 09 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) FILE_NAME, a file name. ! On output, the extension of the file has been changed. ! ! Input, character ( len = * ) EXT, the extension to be used on the output ! copy of FILE_NAME, replacing the current extension if any. ! implicit none character ( len = * ) ext character ( len = * ) file_name integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) len_max integer ( kind = 4 ) len_name len_max = len ( file_name ) len_name = len_trim ( file_name ) call file_name_ext_get ( file_name, i, j ) if ( i == 0 ) then if ( len_max < len_name + 1 ) then return end if len_name = len_name + 1 file_name(len_name:len_name) = '.' i = len_name + 1 else i = i + 1 file_name(i:j) = ' ' end if file_name(i:) = ext return end subroutine file_row_count ( input_filename, row_num ) !*****************************************************************************80 ! !! FILE_ROW_COUNT counts the number of row records in a file. ! ! Discussion: ! ! It does not count lines that are blank, or that begin with a ! comment symbol '#'. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 06 March 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output, integer ( kind = 4 ) ROW_NUM, the number of rows found. ! implicit none integer ( kind = 4 ) bad_num integer ( kind = 4 ) comment_num integer ( kind = 4 ) ierror character ( len = * ) input_filename integer ( kind = 4 ) input_status integer ( kind = 4 ) input_unit character ( len = 255 ) line integer ( kind = 4 ) record_num integer ( kind = 4 ) row_num call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then row_num = -1; ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_ROW_COUNT - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_filename ) // '" on unit ', input_unit stop end if comment_num = 0 row_num = 0 record_num = 0 bad_num = 0 do read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then ierror = record_num exit end if record_num = record_num + 1 if ( line(1:1) == '#' ) then comment_num = comment_num + 1 cycle end if if ( len_trim ( line ) == 0 ) then comment_num = comment_num + 1 cycle end if row_num = row_num + 1 end do close ( unit = input_unit ) return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is a value between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is a value between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 October 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ( kind = 4 ) IUNIT, the free unit number. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ) ios integer ( kind = 4 ) iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end function i4_modp ( i, j ) !*****************************************************************************80 ! !! I4_MODP returns the nonnegative remainder of I4 division. ! ! Discussion: ! ! If ! NREM = I4_MODP ( I, J ) ! NMULT = ( I - NREM ) / J ! then ! I = J * NMULT + NREM ! where NREM is always nonnegative. ! ! The MOD function computes a result with the same sign as the ! quantity being divided. Thus, suppose you had an angle A, ! and you wanted to ensure that it was between 0 and 360. ! Then mod(A,360) would do, if A was positive, but if A ! was negative, your result would be between -360 and 0. ! ! On the other hand, I4_MODP(A,360) is between 0 and 360, always. ! ! Example: ! ! I J MOD I4_MODP Factorization ! ! 107 50 7 7 107 = 2 * 50 + 7 ! 107 -50 7 7 107 = -2 * -50 + 7 ! -107 50 -7 43 -107 = -3 * 50 + 43 ! -107 -50 -7 43 -107 = 3 * -50 + 43 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 02 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) I, the number to be divided. ! ! Input, integer ( kind = 4 ) J, the number that divides I. ! ! Output, integer ( kind = 4 ) I4_MODP, the nonnegative remainder when I is ! divided by J. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ) i4_modp integer ( kind = 4 ) j if ( j == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4_MODP - Fatal error!' write ( *, '(a,i8)' ) ' I4_MODP ( I, J ) called with J = ', j stop end if i4_modp = mod ( i, j ) if ( i4_modp < 0 ) then i4_modp = i4_modp + abs ( j ) end if return end function i4_wrap ( ival, ilo, ihi ) !*****************************************************************************80 ! !! I4_WRAP forces an I4 to lie between given limits by wrapping. ! ! Example: ! ! ILO = 4, IHI = 8 ! ! I I4_WRAP ! ! -2 8 ! -1 4 ! 0 5 ! 1 6 ! 2 7 ! 3 8 ! 4 4 ! 5 5 ! 6 6 ! 7 7 ! 8 8 ! 9 4 ! 10 5 ! 11 6 ! 12 7 ! 13 8 ! 14 4 ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 August 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) IVAL, a value. ! ! Input, integer ( kind = 4 ) ILO, IHI, the desired bounds for the value. ! ! Output, integer ( kind = 4 ) I4_WRAP, a "wrapped" version of IVAL. ! implicit none integer ( kind = 4 ) i4_modp integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) ival integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo integer ( kind = 4 ) wide jlo = min ( ilo, ihi ) jhi = max ( ilo, ihi ) wide = jhi - jlo + 1 if ( wide == 1 ) then i4_wrap = jlo else i4_wrap = jlo + i4_modp ( ival - jlo, wide ) end if return end subroutine i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT_SOME prints some of the transpose of an I4MAT. ! ! Discussion: ! ! An I4MAT is a rectangular array of I4 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 10 September 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, integer ( kind = 4 ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. ! ! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ), parameter :: incx = 10 integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) a(m,n) character ( len = 8 ) ctemp(incx) integer ( kind = 4 ) i integer ( kind = 4 ) i2 integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) 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)' ) i end do write ( *, '('' Row '',10a8)' ) 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), '(i8)' ) a(i,j) end do write ( *, '(i5,a,10a8)' ) j, ':', ( ctemp(i), i = 1, inc ) end do end do return end subroutine i4mat_write ( output_filename, m, n, table ) !*****************************************************************************80 ! !! I4MAT_WRITE writes an I4MAT file. ! ! Discussion: ! ! An I4MAT is an array of I4's. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 31 August 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! Input, integer ( kind = 4 ) M, the spatial dimension. ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Input, integer ( kind = 4 ) TABLE(M,N), the table data. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) j character ( len = * ) output_filename integer ( kind = 4 ) output_status integer ( kind = 4 ) output_unit character ( len = 30 ) string integer ( kind = 4 ) table(m,n) ! ! Open the file. ! call get_unit ( output_unit ) open ( unit = output_unit, file = output_filename, & status = 'replace', iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_WRITE - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the output file "' // & trim ( output_filename ) // '" on unit ', output_unit output_unit = -1 stop end if ! ! Create a format string. ! if ( 0 < m .and. 0 < n ) then write ( string, '(a1,i8,a4)' ) '(', m, 'i10)' ! ! Write the data. ! do j = 1, n write ( output_unit, string ) table(1:m,j) end do end if ! ! Close the file. ! close ( unit = output_unit ) return end subroutine i4vec_indicator ( n, a ) !*****************************************************************************80 ! !! I4VEC_INDICATOR sets an I4VEC to the indicator vector. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 09 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of elements of A. ! ! Output, integer ( kind = 4 ) A(N), the array to be initialized. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a(n) integer ( kind = 4 ) i do i = 1, n a(i) = i end do return end subroutine i4vec_sort_heap_index_a ( n, a, indx ) !*****************************************************************************80 ! !! I4VEC_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an I4VEC. ! ! Discussion: ! ! The sorting is not actually carried out. Rather an index array is ! created which defines the sorting. This array may be used to sort ! or index the array, or to sort or index related arrays keyed on the ! original array. ! ! Once the index array is computed, the sorting can be carried out ! "implicitly: ! ! A(INDX(I)), I = 1 to N is sorted, ! ! or explicitly, by the call ! ! call I4VEC_PERMUTE ( N, A, INDX ) ! ! after which A(I), I = 1 to N is sorted. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 September 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries in the array. ! ! Input, integer ( kind = 4 ) A(N), an array to be index-sorted. ! ! Output, integer ( kind = 4 ) INDX(N), the sort index. The ! I-th element of the sorted array is A(INDX(I)). ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a(n) integer ( kind = 4 ) aval integer ( kind = 4 ) i integer ( kind = 4 ) indx(n) integer ( kind = 4 ) indxt integer ( kind = 4 ) ir integer ( kind = 4 ) j integer ( kind = 4 ) l if ( n <= 1 ) then return end if do i = 1, n indx(i) = i end do l = n / 2 + 1 ir = n do if ( 1 < l ) then l = l - 1 indxt = indx(l) aval = a(indxt) else indxt = indx(ir) aval = a(indxt) indx(ir) = indx(1) ir = ir - 1 if ( ir == 1 ) then indx(1) = indxt exit end if end if i = l j = l + l do while ( j <= ir ) if ( j < ir ) then if ( a(indx(j)) < a(indx(j+1)) ) then j = j + 1 end if end if if ( aval < a(indx(j)) ) then indx(i) = indx(j) i = j j = j + j else j = ir + 1 end if end do indx(i) = indxt end do return end function lrline ( xu, yu, xv1, yv1, xv2, yv2, dv ) !*****************************************************************************80 ! !! LRLINE determines if a point is left of, right or, or on a directed line. ! ! Discussion: ! ! The directed line is paralled to, and at a signed distance DV from ! a directed base line from (XV1,YV1) to (XV2,YV2). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 July 2001 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) XU, YU, the coordinates of the point whose ! position relative to the directed line is to be determined. ! ! Input, real ( kind = 8 ) XV1, YV1, XV2, YV2, the coordinates of two points ! that determine the directed base line. ! ! Input, real ( kind = 8 ) DV, the signed distance of the directed line ! from the directed base line through the points (XV1,YV1) and (XV2,YV2). ! DV is positive for a line to the left of the base line. ! ! Output, integer ( kind = 4 ) LRLINE, the result: ! +1, the point is to the right of the directed line; ! 0, the point is on the directed line; ! -1, the point is to the left of the directed line. ! implicit none real ( kind = 8 ) dv real ( kind = 8 ) dx real ( kind = 8 ) dxu real ( kind = 8 ) dy real ( kind = 8 ) dyu integer ( kind = 4 ) lrline real ( kind = 8 ) t real ( kind = 8 ) tol real ( kind = 8 ) tolabs real ( kind = 8 ) xu real ( kind = 8 ) xv1 real ( kind = 8 ) xv2 real ( kind = 8 ) yu real ( kind = 8 ) yv1 real ( kind = 8 ) yv2 tol = 100.0D+00 * epsilon ( tol ) dx = xv2 - xv1 dy = yv2 - yv1 dxu = xu - xv1 dyu = yu - yv1 tolabs = tol * max ( abs ( dx ), abs ( dy ), abs ( dxu ), & abs ( dyu ), abs ( dv ) ) t = dy * dxu - dx * dyu + dv * sqrt ( dx * dx + dy * dy ) if ( tolabs < t ) then lrline = 1 else if ( -tolabs <= t ) then lrline = 0 else lrline = -1 end if return end subroutine perm_inv ( n, p ) !*****************************************************************************80 ! !! PERM_INV inverts a permutation "in place". ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of objects being permuted. ! ! Input/output, integer ( kind = 4 ) P(N), the permutation, in standard ! index form. On output, P describes the inverse permutation ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) i0 integer ( kind = 4 ) i1 integer ( kind = 4 ) i2 integer ( kind = 4 ) is integer ( kind = 4 ) p(n) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PERM_INV - Fatal error!' write ( *, '(a,i8)' ) ' Input value of N = ', n stop end if is = 1 do i = 1, n i1 = p(i) do while ( i < i1 ) i2 = p(i1) p(i1) = -i2 i1 = i2 end do is = -sign ( 1, p(i) ) p(i) = sign ( p(i), is ) end do do i = 1, n i1 = -p(i) if ( 0 <= i1 ) then i0 = i do i2 = p(i1) p(i1) = i0 if ( i2 < 0 ) then exit end if i0 = i1 i1 = i2 end do end if end do return end subroutine r82vec_permute ( n, a, p ) !*****************************************************************************80 ! !! R82VEC_PERMUTE permutes an R82VEC in place. ! ! Discussion: ! ! This routine permutes an array of real "objects", but the same ! logic can be used to permute an array of objects of any arithmetic ! type, or an array of objects of any complexity. The only temporary ! storage required is enough to store a single object. The number ! of data movements made is N + the number of cycles of order 2 or more, ! which is never more than N + N/2. ! ! Example: ! ! Input: ! ! N = 5 ! P = ( 2, 4, 5, 1, 3 ) ! A = ( 1.0, 2.0, 3.0, 4.0, 5.0 ) ! (11.0, 22.0, 33.0, 44.0, 55.0 ) ! ! Output: ! ! A = ( 2.0, 4.0, 5.0, 1.0, 3.0 ) ! ( 22.0, 44.0, 55.0, 11.0, 33.0 ). ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 11 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of objects. ! ! Input/output, real ( kind = 8 ) A(2,N), the array to be permuted. ! ! Input, integer ( kind = 4 ) P(N), the permutation. P(I) = J means ! that the I-th element of the output array should be the J-th ! element of the input array. P must be a legal permutation ! of the integers from 1 to N, otherwise the algorithm will ! fail catastrophically. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(2,n) real ( kind = 8 ) a_temp(2) integer ( kind = 4 ) iget integer ( kind = 4 ) iput integer ( kind = 4 ) istart integer ( kind = 4 ) p(n) ! ! Search for the next element of the permutation that has not been used. ! do istart = 1, n if ( p(istart) < 0 ) then cycle else if ( p(istart) == istart ) then p(istart) = - p(istart) cycle else a_temp(1:2) = a(1:2,istart) iget = istart ! ! Copy the new value into the vacated entry. ! do iput = iget iget = p(iget) p(iput) = - p(iput) if ( iget < 1 .or. n < iget ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R82VEC_PERMUTE - Fatal error!' stop end if if ( iget == istart ) then a(1:2,iput) = a_temp(1:2) exit end if a(1:2,iput) = a(1:2,iget) end do end if end do ! ! Restore the signs of the entries. ! p(1:n) = -p(1:n) return end subroutine r82vec_sort_heap_index_a ( n, a, indx ) !*****************************************************************************80 ! !! R82VEC_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R82VEC. ! ! Discussion: ! ! The sorting is not actually carried out. Rather an index array is ! created which defines the sorting. This array may be used to sort ! or index the array, or to sort or index related arrays keyed on the ! original array. ! ! Once the index array is computed, the sorting can be carried out ! "implicitly: ! ! A(1:2,INDX(I)), I = 1 to N is sorted, ! ! or explicitly, by the call ! ! call R82VEC_PERMUTE ( N, A, INDX ) ! ! after which A(1:2,I), I = 1 to N is sorted. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 11 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries in the array. ! ! Input, real ( kind = 8 ) A(2,N), an array to be index-sorted. ! ! Output, integer ( kind = 4 ) INDX(N), the sort index. The ! I-th element of the sorted array is A(1:2,INDX(I)). ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(2,n) real ( kind = 8 ) aval(2) integer ( kind = 4 ) i integer ( kind = 4 ) indx(n) integer ( kind = 4 ) indxt integer ( kind = 4 ) ir integer ( kind = 4 ) j integer ( kind = 4 ) l if ( n < 1 ) then return end if if ( n == 1 ) then indx(1) = 1 return end if call i4vec_indicator ( n, indx ) l = n / 2 + 1 ir = n do if ( 1 < l ) then l = l - 1 indxt = indx(l) aval(1:2) = a(1:2,indxt) else indxt = indx(ir) aval(1:2) = a(1:2,indxt) indx(ir) = indx(1) ir = ir - 1 if ( ir == 1 ) then indx(1) = indxt exit end if end if i = l j = l + l do while ( j <= ir ) if ( j < ir ) then if ( a(1,indx(j)) < a(1,indx(j+1)) .or. & ( a(1,indx(j)) == a(1,indx(j+1)) .and. & a(2,indx(j)) < a(2,indx(j+1)) ) ) then j = j + 1 end if end if if ( aval(1) < a(1,indx(j)) .or. & ( aval(1) == a(1,indx(j)) .and. & aval(2) < a(2,indx(j)) ) ) then indx(i) = indx(j) i = j j = j + j else j = ir + 1 end if end do indx(i) = indxt end do return end subroutine r8mat_data_read ( input_filename, m, n, table ) !*****************************************************************************80 ! !! R8MAT_DATA_READ reads data from an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Discussion: ! ! The file may contain more than N points, but this routine will ! return after reading N of them. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 18 October 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Input, integer ( kind = 4 ) M, the spatial dimension. ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Output, real ( kind = 8 ) TABLE(M,N), the table data. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) ierror character ( len = * ) input_filename integer ( kind = 4 ) input_status integer ( kind = 4 ) input_unit integer ( kind = 4 ) j character ( len = 255 ) line real ( kind = 8 ) table(m,n) real ( kind = 8 ) x(m) ierror = 0 call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_DATA_READ - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_filename ) // '" on unit ', input_unit stop end if j = 0 do while ( j < n ) read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_DATA_READ - Fatal error!' write ( *, '(a)' ) ' Error while reading lines of data.' write ( *, '(a,i8)' ) ' Number of values expected per line M = ', m write ( *, '(a,i8)' ) ' Number of data lines read, J = ', j write ( *, '(a,i8)' ) ' Number of data lines needed, N = ', n stop end if if ( line(1:1) == '#' .or. len_trim ( line ) == 0 ) then cycle end if call s_to_r8vec ( line, m, x, ierror ) if ( ierror /= 0 ) then cycle end if j = j + 1 table(1:m,j) = x(1:m) end do close ( unit = input_unit ) return end subroutine r8mat_header_read ( input_filename, m, n ) !*****************************************************************************80 ! !! R8MAT_HEADER_READ reads the header from an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output, integer ( kind = 4 ) M, spatial dimension. ! ! Output, integer ( kind = 4 ) N, the number of points. ! implicit none character ( len = * ) input_filename integer ( kind = 4 ) m integer ( kind = 4 ) n call file_column_count ( input_filename, m ) if ( m <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data columns in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if call file_row_count ( input_filename, n ) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data rows in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if return end subroutine r8mat_print ( m, n, a, title ) !*****************************************************************************80 ! !! R8MAT_PRINT prints an R8MAT. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of rows in A. ! ! Input, integer ( kind = 4 ) N, the number of columns in A. ! ! Input, real ( kind = 8 ) A(M,N), the matrix. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(m,n) character ( len = * ) title call r8mat_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8MAT_PRINT_SOME prints some of an R8MAT. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. ! ! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer ( kind = 4 ), parameter :: incx = 5 integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(m,n) character ( len = 14 ) ctemp(incx) integer ( kind = 4 ) i integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2 integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i8,6x)') j end do write ( *, '('' Col'',5a14)' ) ctemp(1:inc) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi do j2 = 1, inc j = j2lo - 1 + j2 write ( ctemp(j2), '(g14.6)' ) a(i,j) end do write ( *, '(2x,i8,5a14)' ) i, ( ctemp(j), j = 1, inc ) end do end do 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 array of R8 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 June 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. ! ! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. ! ! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer ( kind = 4 ), parameter :: incx = 5 integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(m,n) character ( len = 14 ) ctemp(incx) integer ( kind = 4 ) i integer ( kind = 4 ) i2 integer ( kind = 4 ) i2hi integer ( kind = 4 ) i2lo integer ( kind = 4 ) ihi integer ( kind = 4 ) ilo integer ( kind = 4 ) inc integer ( kind = 4 ) j integer ( kind = 4 ) j2hi integer ( kind = 4 ) j2lo integer ( kind = 4 ) jhi integer ( kind = 4 ) jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) 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,5a14)' ) j, ( ctemp(i), i = 1, inc ) end do end do return end subroutine s_blank_delete ( s ) !*****************************************************************************80 ! !! S_BLANK_DELETE removes blanks from a string, left justifying the remainder. ! ! Discussion: ! ! All TAB characters are also removed. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) S, the string to be transformed. ! implicit none character c integer ( kind = 4 ) get integer ( kind = 4 ) put integer ( kind = 4 ) nchar character ( len = * ) s character, parameter :: TAB = char ( 9 ) put = 0 nchar = len_trim ( s ) do get = 1, nchar c = s(get:get) if ( c /= ' ' .and. c /= TAB ) then put = put + 1 s(put:put) = c end if end do s(put+1:nchar) = ' ' return end function s_index_last ( s, sub ) !*****************************************************************************80 ! !! S_INDEX_LAST finds the LAST occurrence of a given substring. ! ! Discussion: ! ! It returns the location in the string at which the substring SUB is ! first found, or 0 if the substring does not occur at all. ! ! The routine is also trailing blank insensitive. This is very ! important for those cases where you have stored information in ! larger variables. If S is of length 80, and SUB is of ! length 80, then if S = 'FRED' and SUB = 'RED', a match would ! not be reported by the standard FORTRAN INDEX, because it treats ! both variables as being 80 characters long! This routine assumes that ! trailing blanks represent garbage! ! ! This means that this routine cannot be used to find, say, the last ! occurrence of a substring 'A ', since it assumes the blank space ! was not specified by the user, but is, rather, padding by the ! system. However, as a special case, this routine can properly handle ! the case where either S or SUB is all blanks. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be searched. ! ! Input, character ( len = * ) SUB, the substring to search for. ! ! Output, integer ( kind = 4 ) S_INDEX_LAST. 0 if SUB does not occur in ! the string. Otherwise S_INDEX_LAST = I, where S(I:I+LENS-1) = SUB, ! where LENS is the length of SUB, and is the last place ! this happens. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) llen1 integer ( kind = 4 ) llen2 character ( len = * ) s integer ( kind = 4 ) s_index_last character ( len = * ) sub s_index_last = 0 llen1 = len_trim ( s ) llen2 = len_trim ( sub ) ! ! In case S or SUB is blanks, use LEN ! if ( llen1 == 0 ) then llen1 = len ( s ) end if if ( llen2 == 0 ) then llen2 = len ( sub ) end if if ( llen1 < llen2 ) then return end if do j = 1, llen1+1-llen2 i = llen1 + 2 - llen2 - j if ( s(i:i+llen2-1) == sub ) then s_index_last = i return end if end do return end subroutine s_to_r8 ( s, dval, ierror, length ) !*****************************************************************************80 ! !! S_TO_R8 reads an R8 from a string. ! ! Discussion: ! ! The routine will read as many characters as possible until it reaches ! the end of the string, or encounters a character which cannot be ! part of the number. ! ! Legal input is: ! ! 1 blanks, ! 2 '+' or '-' sign, ! 2.5 blanks ! 3 integer part, ! 4 decimal point, ! 5 fraction part, ! 6 'E' or 'e' or 'D' or 'd', exponent marker, ! 7 exponent sign, ! 8 exponent integer part, ! 9 exponent decimal point, ! 10 exponent fraction part, ! 11 blanks, ! 12 final comma or semicolon, ! ! with most quantities optional. ! ! Example: ! ! S DVAL ! ! '1' 1.0 ! ' 1 ' 1.0 ! '1A' 1.0 ! '12,34,56' 12.0 ! ' 34 7' 34.0 ! '-1E2ABCD' -100.0 ! '-1X2ABCD' -1.0 ! ' 2E-1' 0.2 ! '23.45' 23.45 ! '-4.2E+2' -420.0 ! '17d2' 1700.0 ! '-14e-2' -0.14 ! 'e2' 100.0 ! '-12.73e-9.23' -12.73 * 10.0^(-9.23) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string containing the ! data to be read. Reading will begin at position 1 and ! terminate at the end of the string, or when no more ! characters can be read to form a legal real. Blanks, ! commas, or other nonnumeric data will, in particular, ! cause the conversion to halt. ! ! Output, real ( kind = 8 ) DVAL, the value read from the string. ! ! Output, integer ( kind = 4 ) IERROR, error flag. ! 0, no errors occurred. ! 1, 2, 6 or 7, the input number was garbled. The ! value of IERROR is the last type of input successfully ! read. For instance, 1 means initial blanks, 2 means ! a plus or minus sign, and so on. ! ! Output, integer ( kind = 4 ) LENGTH, the number of characters read ! to form the number, including any terminating ! characters such as a trailing comma or blanks. ! implicit none character c logical ch_eqi real ( kind = 8 ) dval integer ( kind = 4 ) ierror integer ( kind = 4 ) ihave integer ( kind = 4 ) isgn integer ( kind = 4 ) iterm integer ( kind = 4 ) jbot integer ( kind = 4 ) jsgn integer ( kind = 4 ) jtop integer ( kind = 4 ) length integer ( kind = 4 ) nchar integer ( kind = 4 ) ndig real ( kind = 8 ) rbot real ( kind = 8 ) rexp real ( kind = 8 ) rtop character ( len = * ) s nchar = len_trim ( s ) ierror = 0 dval = 0.0D+00 length = -1 isgn = 1 rtop = 0 rbot = 1 jsgn = 1 jtop = 0 jbot = 1 ihave = 1 iterm = 0 do length = length + 1 if ( nchar < length+1 ) then exit end if c = s(length+1:length+1) ! ! Blank character. ! if ( c == ' ' ) then if ( ihave == 2 ) then else if ( ihave == 6 .or. ihave == 7 ) then iterm = 1 else if ( 1 < ihave ) then ihave = 11 end if ! ! Comma. ! else if ( c == ',' .or. c == ';' ) then if ( ihave /= 1 ) then iterm = 1 ihave = 12 length = length + 1 end if ! ! Minus sign. ! else if ( c == '-' ) then if ( ihave == 1 ) then ihave = 2 isgn = -1 else if ( ihave == 6 ) then ihave = 7 jsgn = -1 else iterm = 1 end if ! ! Plus sign. ! else if ( c == '+' ) then if ( ihave == 1 ) then ihave = 2 else if ( ihave == 6 ) then ihave = 7 else iterm = 1 end if ! ! Decimal point. ! else if ( c == '.' ) then if ( ihave < 4 ) then ihave = 4 else if ( 6 <= ihave .and. ihave <= 8 ) then ihave = 9 else iterm = 1 end if ! ! Scientific notation exponent marker. ! else if ( ch_eqi ( c, 'E' ) .or. ch_eqi ( c, 'D' ) ) then if ( ihave < 6 ) then ihave = 6 else iterm = 1 end if ! ! Digit. ! else if ( ihave < 11 .and. lle ( '0', c ) .and. lle ( c, '9' ) ) then if ( ihave <= 2 ) then ihave = 3 else if ( ihave == 4 ) then ihave = 5 else if ( ihave == 6 .or. ihave == 7 ) then ihave = 8 else if ( ihave == 9 ) then ihave = 10 end if call ch_to_digit ( c, ndig ) if ( ihave == 3 ) then rtop = 10.0D+00 * rtop + real ( ndig, kind = 8 ) else if ( ihave == 5 ) then rtop = 10.0D+00 * rtop + real ( ndig, kind = 8 ) rbot = 10.0D+00 * rbot else if ( ihave == 8 ) then jtop = 10 * jtop + ndig else if ( ihave == 10 ) then jtop = 10 * jtop + ndig jbot = 10 * jbot end if ! ! Anything else is regarded as a terminator. ! else iterm = 1 end if ! ! If we haven't seen a terminator, and we haven't examined the ! entire string, go get the next character. ! if ( iterm == 1 ) then exit end if end do ! ! If we haven't seen a terminator, and we have examined the ! entire string, then we're done, and LENGTH is equal to NCHAR. ! if ( iterm /= 1 .and. length + 1 == nchar ) then length = nchar end if ! ! Number seems to have terminated. Have we got a legal number? ! Not if we terminated in states 1, 2, 6 or 7! ! if ( ihave == 1 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then ierror = ihave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'S_TO_R8 - Serious error!' write ( *, '(a)' ) ' Illegal or nonnumeric input:' write ( *, '(a)' ) ' ' // trim ( s ) return end if ! ! Number seems OK. Form it. ! if ( jtop == 0 ) then rexp = 1.0D+00 else if ( jbot == 1 ) then rexp = 10.0D+00 ** ( jsgn * jtop ) else rexp = 10.0D+00 ** ( real ( jsgn * jtop, kind = 8 ) & / real ( jbot, kind = 8 ) ) end if end if dval = real ( isgn, kind = 8 ) * rexp * rtop / rbot return end subroutine s_to_r8vec ( s, n, rvec, ierror ) !*****************************************************************************80 ! !! S_TO_R8VEC reads an R8VEC from a string. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 07 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be read. ! ! Input, integer ( kind = 4 ) N, the number of values expected. ! ! Output, real ( kind = 8 ) RVEC(N), the values read from the string. ! ! Output, integer ( kind = 4 ) IERROR, error flag. ! 0, no errors occurred. ! -K, could not read data for entries -K through N. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) ierror integer ( kind = 4 ) ilo integer ( kind = 4 ) lchar real ( kind = 8 ) rvec(n) character ( len = * ) s i = 0 ierror = 0 ilo = 1 do while ( i < n ) i = i + 1 call s_to_r8 ( s(ilo:), rvec(i), ierror, lchar ) if ( ierror /= 0 ) then ierror = -i exit end if ilo = ilo + lchar end do return end subroutine s_word_count ( s, nword ) !*****************************************************************************80 ! !! S_WORD_COUNT counts the number of "words" in a string. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be examined. ! ! Output, integer ( kind = 4 ) NWORD, the number of "words" in the string. ! Words are presumed to be separated by one or more blanks. ! implicit none logical blank integer ( kind = 4 ) i integer ( kind = 4 ) lens integer ( kind = 4 ) nword character ( len = * ) s nword = 0 lens = len ( s ) if ( lens <= 0 ) then return end if blank = .true. do i = 1, lens if ( s(i:i) == ' ' ) then blank = .true. else if ( blank ) then nword = nword + 1 blank = .false. end if end do return end subroutine swapec_lmap ( i, matrix, top, btri, bedg, point_num, point_xy, & tri_num, tri_vert, tri_nabe, stack, ierr ) !*****************************************************************************80 ! !! SWAPEC_LMAP swaps diagonal edges until all triangles are Delaunay. ! ! Discussion: ! ! The routine swaps diagonal edges in a 2D triangulation, based on ! the empty circumcircle criterion, until all triangles are Delaunay, ! given that I is the index of the new vertex added to the triangulation. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 July 2001 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, integer ( kind = 4 ) I, the index of the new vertex. ! ! Input, real ( kind = 8 ) MATRIX(2,2), the transformation matrix. ! ! Input/output, integer ( kind = 4 ) TOP, the index of the top of the stack. ! On output, TOP is zero. ! ! Input/output, integer ( kind = 4 ) BTRI, BEDG; on input, if positive, are ! the triangle and edge indices of a boundary edge whose updated indices ! must be recorded. On output, these may be updated because of swaps. ! ! Input, intger POINT_NUM, the number of points. ! ! Input, real ( kind = 8 ) POINT_XY(2,POINT_NUM), the coordinates ! of the points. ! ! Input, integer ( kind = 4 ) TRI_NUM, the number of triangles. ! ! Input/output, integer ( kind = 4 ) TRI_VERT(3,TRI_NUM), the triangle ! incidence list. May be updated on output because of swaps. ! ! Input/output, integer ( kind = 4 ) TRI_NABE(3,TRI_NUM), the triangle ! neighbor list; negative values are used for links of the counter-clockwise ! linked list of boundary edges; May be updated on output because of swaps. ! LINK = -(3*I + J-1) where I, J = triangle, edge index. ! ! Workspace, integer STACK(MAXST); on input, entries 1 through TOP ! contain the indices of initial triangles (involving vertex I) ! put in stack; the edges opposite I should be in interior; entries ! TOP+1 through MAXST are used as a stack. ! ! Output, integer ( kind = 4 ) IERR is set to 8 for abnormal return. ! implicit none integer ( kind = 4 ) point_num integer ( kind = 4 ) tri_num integer ( kind = 4 ) a integer ( kind = 4 ) b integer ( kind = 4 ) bedg integer ( kind = 4 ) btri integer ( kind = 4 ) c integer ( kind = 4 ) diaedg_lmap integer ( kind = 4 ) e integer ( kind = 4 ) ee integer ( kind = 4 ) em1 integer ( kind = 4 ) ep1 integer ( kind = 4 ) f integer ( kind = 4 ) fm1 integer ( kind = 4 ) fp1 integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) l real ( kind = 8 ) matrix(2,2) integer ( kind = 4 ) r integer ( kind = 4 ) s integer ( kind = 4 ) stack(point_num) integer ( kind = 4 ) swap integer ( kind = 4 ) t integer ( kind = 4 ) top integer ( kind = 4 ) tri_nabe(3,tri_num) integer ( kind = 4 ) tri_vert(3,tri_num) integer ( kind = 4 ) tt integer ( kind = 4 ) u real ( kind = 8 ) point_xy(2,point_num) real ( kind = 8 ) x real ( kind = 8 ) y ! ! Determine whether triangles in stack are Delaunay, and swap ! diagonal edge of convex quadrilateral if not. ! x = point_xy(1,i) y = point_xy(2,i) do if ( top <= 0 ) then exit end if t = stack(top) top = top - 1 if ( tri_vert(1,t) == i ) then e = 2 b = tri_vert(3,t) else if ( tri_vert(2,t) == i ) then e = 3 b = tri_vert(1,t) else e = 1 b = tri_vert(2,t) end if a = tri_vert(e,t) u = tri_nabe(e,t) if ( tri_nabe(1,u) == t ) then f = 1 c = tri_vert(3,u) else if ( tri_nabe(2,u) == t ) then f = 2 c = tri_vert(1,u) else f = 3 c = tri_vert(2,u) end if swap = diaedg_lmap ( x, y, point_xy(1,a), point_xy(2,a), point_xy(1,c), & point_xy(2,c), point_xy(1,b), point_xy(2,b), matrix ) if ( swap == 1 ) then em1 = i4_wrap ( e - 1, 1, 3 ) ep1 = i4_wrap ( e + 1, 1, 3 ) fm1 = i4_wrap ( f - 1, 1, 3 ) fp1 = i4_wrap ( f + 1, 1, 3 ) tri_vert(ep1,t) = c tri_vert(fp1,u) = i r = tri_nabe(ep1,t) s = tri_nabe(fp1,u) tri_nabe(ep1,t) = u tri_nabe(fp1,u) = t tri_nabe(e,t) = s tri_nabe(f,u) = r if ( 0 < tri_nabe(fm1,u) ) then top = top + 1 stack(top) = u end if if ( 0 < s ) then if ( tri_nabe(1,s) == u ) then tri_nabe(1,s) = t else if ( tri_nabe(2,s) == u ) then tri_nabe(2,s) = t else tri_nabe(3,s) = t end if top = top + 1 if ( point_num < top ) then ierr = 8 return end if stack(top) = t else if ( u == btri .and. fp1 == bedg ) then btri = t bedg = e end if l = - ( 3 * t + e - 1 ) tt = t ee = em1 do while ( 0 < tri_nabe(ee,tt) ) tt = tri_nabe(ee,tt) if ( tri_vert(1,tt) == a ) then ee = 3 else if ( tri_vert(2,tt) == a ) then ee = 1 else ee = 2 end if end do tri_nabe(ee,tt) = l end if if ( 0 < r ) then if ( tri_nabe(1,r) == t ) then tri_nabe(1,r) = u else if ( tri_nabe(2,r) == t ) then tri_nabe(2,r) = u else tri_nabe(3,r) = u end if else if ( t == btri .and. ep1 == bedg ) then btri = u bedg = f end if l = - ( 3 * u + f - 1 ) tt = u ee = fm1 do while ( 0 < tri_nabe(ee,tt) ) tt = tri_nabe(ee,tt) if ( tri_vert(1,tt) == b ) then ee = 3 else if ( tri_vert(2,tt) == b ) then ee = 1 else ee = 2 end if end do tri_nabe(ee,tt) = l end if end if 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 GNU LGPL license. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer ( kind = 4 ) d integer ( kind = 4 ) h integer ( kind = 4 ) m integer ( kind = 4 ) mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer ( kind = 4 ) n integer ( kind = 4 ) s integer ( kind = 4 ) values(8) integer ( kind = 4 ) 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 subroutine transform_lmap ( matrix, dx, dy ) !*****************************************************************************80 ! !! TRANSFORM_LMAP applies a transformation matrix to a vector. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 16 April 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) MATRIX(2,2), the linear map to be applied ! to the data. ! ! Input/output, real ( kind = 8 ) DX, DY, the components of the vector. ! implicit none real ( kind = 8 ) dx real ( kind = 8 ) dx2 real ( kind = 8 ) dy real ( kind = 8 ) dy2 real ( kind = 8 ) matrix(2,2) dx2 = matrix(1,1) * dx + matrix(1,2) * dy dy2 = matrix(2,1) * dx + matrix(2,2) * dy dx = dx2 dy = dy2 return end subroutine triangulation_plot_eps ( file_name, node_num, node_x, node_y, & element_num, element_mask, element_node, title ) !*****************************************************************************80 ! !! TRIANGULATION_PLOT_EPS creates an EPS file image of a triangulation. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 15 April 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the file to create. ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) NODE_X(NODE_NUM), NODE_Y(NODE_NUM), ! the coordinates of the nodes. ! ! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Input, logical ELEMENT_MASK(ELEMENT_NUM), a mask for the elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NODE(3,ELEMENT_NUM), the ! element->node data. ! ! Input, character ( len = * ) TITLE, a title for the plot. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) node_num real ( kind = 8 ) ave_x real ( kind = 8 ) ave_y integer ( kind = 4 ), parameter :: circle_size = 3 real ( kind = 8 ) dif integer ( kind = 4 ) element logical element_mask(element_num) integer ( kind = 4 ) element_node(3,element_num) integer ( kind = 4 ) eps_unit integer ( kind = 4 ) eps_x integer ( kind = 4 ) eps_y character ( len = * ) file_name integer ( kind = 4 ) i integer ( kind = 4 ) ios integer ( kind = 4 ) j integer ( kind = 4 ) local integer ( kind = 4 ) node logical node_mask(node_num) real ( kind = 8 ) node_x(node_num) real ( kind = 8 ) node_x_max real ( kind = 8 ) node_x_min real ( kind = 8 ) node_y(node_num) real ( kind = 8 ) node_y_max real ( kind = 8 ) node_y_min real ( kind = 8 ) scale character ( len = 40 ) string character ( len = * ) title ! ! Determine the range of the unmasked elements. ! node_x_min = huge ( node_x_min ) node_x_max = -huge ( node_x_max ) node_y_min = huge ( node_y_min ) node_y_max = -huge ( node_y_max ) node_mask(1:node_num) = .false. do element = 1, element_num if ( element_mask(element) ) then do j = 1, 3 node = element_node(j,element) node_mask(node) = .true. node_x_min = min ( node_x_min, node_x(node) ) node_x_max = max ( node_x_max, node_x(node) ) node_y_min = min ( node_y_min, node_y(node) ) node_y_max = max ( node_y_max, node_y(node) ) end do end if end do if ( node_y_max - node_y_min < node_x_max - node_x_min ) then scale = node_x_max - node_x_min dif = ( node_x_max - node_x_min ) - ( node_y_max - node_y_min ) node_y_max = node_y_max + 0.5D+00 * dif node_y_min = node_y_min - 0.5D+00 * dif else scale = node_y_max - node_y_min dif = ( node_y_max - node_y_min ) - ( node_x_max - node_x_min ) node_x_max = node_x_max + 0.5D+00 * dif node_x_min = node_x_min - 0.5D+00 * dif end if call get_unit ( eps_unit ) open ( unit = eps_unit, file = file_name, status = 'replace', iostat = ios ) if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TRIANGULATION_PLOT_EPS - Fatal error!' write ( *, '(a)' ) ' Could not open the output EPS file.' stop end if write ( eps_unit, '(a)' ) '%!PS-Adobe-3.0 EPSF-3.0' write ( eps_unit, '(a)' ) & '%%Creator: triangulation_plot_eps(delaunay_lmap_2d.f90)' write ( eps_unit, '(a)' ) '%%Title: ' // trim ( file_name ) write ( eps_unit, '(a)' ) '%%Pages: 1' write ( eps_unit, '(a)' ) '%%BoundingBox: 36 36 576 756' write ( eps_unit, '(a)' ) '%%Document-Fonts: Times-Roman' write ( eps_unit, '(a)' ) '%%LanguageLevel: 1' write ( eps_unit, '(a)' ) '%%EndComments' write ( eps_unit, '(a)' ) '%%BeginProlog' write ( eps_unit, '(a)' ) '/inch {72 mul} def' write ( eps_unit, '(a)' ) '%%EndProlog' write ( eps_unit, '(a)' ) '%%Page: 1 1' write ( eps_unit, '(a)' ) 'save' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Set RGB line color.' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.9000 0.9000 0.9000 setrgbcolor' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Draw a gray border around the page.' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) 'newpath' write ( eps_unit, '(a)' ) ' 36 126 moveto' write ( eps_unit, '(a)' ) ' 576 126 lineto' write ( eps_unit, '(a)' ) ' 576 666 lineto' write ( eps_unit, '(a)' ) ' 36 666 lineto' write ( eps_unit, '(a)' ) ' 36 126 lineto' write ( eps_unit, '(a)' ) 'stroke' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Set RGB line color.' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.0000 0.0000 0.0000 setrgbcolor' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Label the plot:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.0000 0.0000 0.0000 setrgbcolor' write ( eps_unit, '(a)' ) '/Times-Roman findfont 0.50 inch scalefont setfont' write ( eps_unit, '(a)' ) ' 36 666 moveto' write ( eps_unit, '(a)' ) '(' // trim ( title ) // ') show' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Define a clipping polygon' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 36 126 moveto' write ( eps_unit, '(a)' ) ' 576 126 lineto' write ( eps_unit, '(a)' ) ' 576 666 lineto' write ( eps_unit, '(a)' ) ' 36 666 lineto' write ( eps_unit, '(a)' ) ' 36 126 lineto' write ( eps_unit, '(a)' ) 'clip newpath' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Draw filled dots at each node:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.0000 0.0000 0.9000 setrgbcolor' do node = 1, node_num if ( node_mask(node) ) then eps_x = int & ( ( node_x_max - node_x(node) ) * 61.0D+00 & + ( + node_x(node) - node_x_min ) * 551.0D+00 ) & / scale eps_y = int & ( ( node_y_max - node_y(node) ) * 151.0D+00 & + ( node_y(node) - node_y_min ) * 641.0D+00 ) & / scale write ( eps_unit, '(a,i4,2x,i4,2x,i4,a)' ) & 'newpath ', eps_x, eps_y, circle_size, ' 0 360 arc closepath fill' end if end do ! ! Label the nodes. ! if ( node_num <= 200 ) then write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Label the nodes:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.0000 0.0000 1.0000 setrgbcolor' write ( eps_unit, '(a)' ) & '/Times-Roman findfont 0.20 inch scalefont setfont' do node = 1, node_num if ( node_mask(node) ) then eps_x = int & ( ( node_x_max - node_x(node) ) * 61.0D+00 & + ( + node_x(node) - node_x_min ) * 551.0D+00 ) & / scale eps_y = int & ( ( node_y_max - node_y(node) ) * 151.0D+00 & + ( node_y(node) - node_y_min ) * 641.0D+00 ) & / scale write ( string, '(i4)' ) node string = adjustl ( string ) write ( eps_unit, '(i4,2x,i4,a)' ) eps_x, eps_y+5, & ' moveto (' // trim ( string ) // ') show' end if end do end if write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Draw the element sides:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.9000 0.0000 0.0000 setrgbcolor' do element = 1, element_num if ( .not. element_mask(element) ) then cycle end if local = 1 node = element_node(local,element) eps_x = int & ( ( node_x_max - node_x(node) ) * 61.0D+00 & + ( + node_x(node) - node_x_min ) * 551.0D+00 ) & / scale eps_y = int & ( ( node_y_max - node_y(node) ) * 151.0D+00 & + ( node_y(node) - node_y_min ) * 641.0D+00 ) & / scale write ( eps_unit, '(a,i4,2x,i4,a)' ) 'newpath ', eps_x, eps_y, ' moveto' do local = mod ( local, 3 ) + 1 node = element_node(local,element) eps_x = int & ( ( node_x_max - node_x(node) ) * 61.0D+00 & + ( + node_x(node) - node_x_min ) * 551.0D+00 ) & / scale eps_y = int & ( ( node_y_max - node_y(node) ) * 151.0D+00 & + ( node_y(node) - node_y_min ) * 641.0D+00 ) & / scale write ( eps_unit, '(i4,2x,i4,a)' ) eps_x, eps_y, ' lineto' if ( local == 1 ) then exit end if end do write ( eps_unit, '(a)' ) 'stroke' end do ! ! Label the elements. ! (This is usually turned off here, because it clutters up the output.) ! if ( element_num <= 20 ) then write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Label the elements:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 1.0000 0.0000 0.0000 setrgbcolor' write ( eps_unit, '(a)' ) & '/Times-Roman findfont 0.30 inch scalefont setfont' do element = 1, element_num if ( .not. element_mask(element) ) then cycle end if ave_x = 0.0D+00 ave_y = 0.0D+00 do i = 1, 3 node = element_node(i,element) ave_x = ave_x + node_x(node) ave_y = ave_y + node_y(node) end do ave_x = ave_x / 3.0D+00 ave_y = ave_y / 3.0D+00 eps_x = int & ( ( node_x_max - ave_x ) * 61.0D+00 & + ( + ave_x - node_x_min ) * 551.0D+00 ) & / scale eps_y = int & ( ( node_y_max - ave_y ) * 151.0D+00 & + ( ave_y - node_y_min ) * 641.0D+00 ) & / scale write ( string, '(i4)' ) element string = adjustl ( string ) write ( eps_unit, '(i4,2x,i4,a)' ) eps_x, eps_y, ' moveto (' & // trim ( string ) // ') show' end do end if write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) 'restore showpage' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% End of page' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '%%Trailer' write ( eps_unit, '(a)' ) '%%EOF' close ( unit = eps_unit ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TRIANGULATION_PLOT_EPS:' write ( *, '(a)' ) ' An encapsulated PostScript file was created' write ( *, '(a)' ) ' containing an image of the triangulation.' write ( *, '(a)' ) ' The file is named "' // trim ( file_name ) // '".' return end subroutine vbedg ( x, y, point_num, point_xy, tri_num, tri_vert, tri_nabe, & ltri, ledg, rtri, redg ) !*****************************************************************************80 ! !! VBEDG determines which boundary edges are visible to a point. ! ! Discussion: ! ! The point (X,Y) is assumed to be outside the convex hull of the ! region covered by the 2D triangulation. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 19 January 2009 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Modified: ! ! 25 August 2001 ! ! Parameters: ! ! Input, real ( kind = 8 ) X, Y, the coordinates of a point outside ! the convex hull of the current triangulation. ! ! Input, integer ( kind = 4 ) POINT_NUM, the number of points. ! ! Input, real ( kind = 8 ) POINT_XY(2,POINT_NUM), the coordinates ! of the vertices. ! ! Input, integer ( kind = 4 ) TRI_NUM, the number of triangles. ! ! Input, integer ( kind = 4 ) TRI_VERT(3,TRI_NUM), the triangle ! incidence list. ! ! Input, integer ( kind = 4 ) TRI_NABE(3,TRI_NUM), the triangle neighbor ! list; negative values are used for links of a counter clockwise linked ! list of boundary edges; ! LINK = -(3*I + J-1) where I, J = triangle, edge index. ! ! Input/output, integer ( kind = 4 ) LTRI, LEDG. If LTRI /= 0 then these ! values are assumed to be already computed and are not changed, else they ! are updated. On output, LTRI is the index of boundary triangle to the left ! of the leftmost boundary triangle visible from (X,Y), and LEDG is the ! boundary edge of triangle LTRI to the left of the leftmost boundary edge ! visible from (X,Y). 1 <= LEDG <= 3. ! ! Input/output, integer ( kind = 4 ) RTRI. On input, the index of the ! boundary triangle to begin the search at. On output, the index of the ! rightmost boundary triangle visible from (X,Y). ! ! Input/output, integer ( kind = 4 ) REDG, the edge of triangle RTRI that ! is visible from (X,Y). 1 <= REDG <= 3. ! implicit none integer ( kind = 4 ) point_num integer ( kind = 4 ) tri_num integer ( kind = 4 ) a integer ( kind = 4 ) b integer ( kind = 4 ) e integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) l logical ldone integer ( kind = 4 ) ledg integer ( kind = 4 ) lr integer ( kind = 4 ) lrline integer ( kind = 4 ) ltri real ( kind = 8 ) point_xy(2,point_num) integer ( kind = 4 ) redg integer ( kind = 4 ) rtri integer ( kind = 4 ) t integer ( kind = 4 ) tri_nabe(3,tri_num) integer ( kind = 4 ) tri_vert(3,tri_num) real ( kind = 8 ) x real ( kind = 8 ) y ! ! Find the rightmost visible boundary edge using links, then possibly ! leftmost visible boundary edge using triangle neighbor information. ! if ( ltri == 0 ) then ldone = .false. ltri = rtri ledg = redg else ldone = .true. end if do l = -tri_nabe(redg,rtri) t = l / 3 e = mod ( l, 3 ) + 1 a = tri_vert(e,t) if ( e <= 2 ) then b = tri_vert(e+1,t) else b = tri_vert(1,t) end if lr = lrline ( x, y, point_xy(1,a), point_xy(2,a), point_xy(1,b), & point_xy(2,b), 0.0D+00 ) if ( lr <= 0 ) then exit end if rtri = t redg = e end do if ( ldone ) then return end if t = ltri e = ledg do b = tri_vert(e,t) e = i4_wrap ( e - 1, 1, 3 ) do while ( 0 < tri_nabe(e,t) ) t = tri_nabe(e,t) if ( tri_vert(1,t) == b ) then e = 3 else if ( tri_vert(2,t) == b ) then e = 1 else e = 2 end if end do a = tri_vert(e,t) lr = lrline ( x, y, point_xy(1,a), point_xy(2,a), point_xy(1,b), & point_xy(2,b), 0.0D+00 ) if ( lr <= 0 ) then exit end if end do ltri = t ledg = e return end