program plot_points ! !******************************************************************************* ! !! PLOT_POINTS plots the points in a file. ! ! ! Discussion: ! ! The program can be invoked by: ! ! plot_points point_file_name ps_file_name ! ! or: ! ! plot_points plot_file_name ! ! in which case the output file name will be constructed by replacing ! the extension of PLOT_FILE_NAME by ".ps", or: ! ! plot_points ! ! in which case the user will be asked to supply both file names. ! ! Modified: ! ! 06 March 2002 ! ! Author: ! ! John Burkardt ! implicit none ! character ( len = 100 ) command integer, allocatable, dimension ( : ) :: connect real, allocatable, dimension ( :, : ) :: coord real, allocatable, dimension ( : ) :: coord_max real, allocatable, dimension ( : ) :: coord_min logical, parameter :: delaunay = .true. character ( len = 256 ) :: file_in_name = ' ' integer file_in_unit character ( len = 256 ) :: file_out_name = 'plot000.eps' real, parameter :: grace = 0.05E+00 integer i integer iarg integer iargc integer :: ierror = 0 integer ilen integer ios integer ipxfargc integer lens integer arg_num integer dim_num integer point_num integer point_num2 character ( len = 80 ) point_file integer, allocatable, dimension ( : ) :: tag real tol logical, parameter :: voronoi = .true. real width integer x_index integer y_index ! call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_POINTS' write ( *, '(a)' ) ' Given a file of 2D points, make a dot plot, and' write ( *, '(a)' ) ' a Delaunay triangulation plot.' ! ! Get the number of command line arguments. ! ! Old style: ! arg_num = iargc ( ) ! ! New style: ! ! arg_num = ipxfargc ( ) ! ! If at least one command line argument, it's the input file name. ! if ( arg_num < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter the input PLOT file name:' read ( *, '(a)', iostat = ios ) file_in_name if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_POINTS - Fatal error!' write ( *, '(a)' ) ' Unexpected read error!' stop end if else iarg = 1 ! ! Old style: ! call getarg ( iarg, file_in_name ) ! ! New style: ! ! call pxfgetarg ( iarg, file_in_name, ilen, ierror ) ! ! if ( ierror /= 0 ) then ! write ( *, '(a)' ) ' ' ! write ( *, '(a)' ) 'PLOT_POINTS - Fatal error!' ! write ( *, '(a)' ) ' Could not read command line argument.' ! stop ! end if end if ! ! If two command line arguments, the second one is the output file name. ! ! if ( arg_num < 2 ) then ! ! filedot_name = file_in_name ! call file_name_ext_swap ( filedot_name, 'ps' ) ! ! else ! ! iarg = 2 ! ! Old style: ! ! call getarg ( iarg, filedot_name ) ! ! New style: ! ! call pxfgetarg ( iarg, filedot_name, ilen, ierror ) ! ! if ( ierror /= 0 ) then ! write ( *, '(a)' ) ' ' ! write ( *, '(a)' ) 'PLOT_POINTS - Fatal error!' ! write ( *, '(a)' ) ' Could not read command line argument.' ! stop ! end if ! ! end if ! ! Now we know what to do. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_POINTS' write ( *, '(a)' ) ' Read point data from "' // trim ( file_in_name ) // '".' ! ! Get the "width" of the file. ! call file_column_count ( file_in_name, dim_num ) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' Points are assumed to have dimension ', dim_num ! ! Get the "length" of the file. ! call points_count ( file_in_name, dim_num, point_num ) ! ! Now allocate some data. ! allocate ( connect(point_num) ) allocate ( coord(dim_num,point_num) ) allocate ( coord_min(dim_num) ) allocate ( coord_max(dim_num) ) allocate ( tag(point_num) ) ! ! Read the coordinates into COORD. ! call points_read ( connect, file_in_name, dim_num, point_num, coord ) ! ! Determine the range. ! do i = 1, dim_num coord_min(i) = minval ( coord(i,1:point_num) ) coord_max(i) = maxval ( coord(i,1:point_num) ) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Coord Min Max' write ( *, '(a)' ) ' ' do i = 1, dim_num write ( *, '(i3,2f10.4)' ) i, coord_min(i), coord_max(i) end do ! ! Extend the range a bit. ! do i = 1, dim_num width = coord_max(i) - coord_min(i) coord_max(i) = coord_max(i) + grace * width coord_min(i) = coord_min(i) - grace * width end do x_index = 1 y_index = 2 ! ! Command loop begins here. ! do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter command:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DASH = dash plot of connected points.' write ( *, '(a)' ) 'DOT = dot plot of points.' write ( *, '(a)' ) 'DEL = Delaunay triangulation plot' write ( *, '(a)' ) 'KM = K-Means plot.' write ( *, '(a)' ) 'TH = thin the points.' write ( *, '(a)' ) 'TV = triangulated Voronoi diagram.' write ( *, '(a)' ) 'VOR = Voronoi diagram plot.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Q = quit.' write ( *, '(a,i6)' ) 'X = specify X index, currently ', x_index write ( *, '(a,i6)' ) 'Y = specify Y index, currently ', y_index write ( *, '(a)' ) ' ' read ( *, '(a)' ) command call s_cap ( command ) ! ! DASH ! Create a dash plot of connected points. ! if ( command(1:3) == 'DAS' ) then call file_name_inc ( file_out_name ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Write dash data to "' // & trim ( file_out_name ) // '".' call dash_plot ( connect, file_out_name, dim_num, point_num, coord, & coord_min, coord_max, x_index, y_index ) ! ! DELAUNAY ! Create a Delaunay plot. ! else if ( command(1:3) == 'DEL' ) then call file_name_inc ( file_out_name ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Create Delaunay plot "' // & trim ( file_out_name ) // '".' call delaunay_plot ( file_out_name, dim_num, point_num, coord, & coord_min, coord_max, x_index, y_index ) ! ! DOT ! Create a dot plot of points. ! else if ( command(1:3) == 'DOT' ) then call file_name_inc ( file_out_name ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Write dot data to "' // & trim ( file_out_name ) // '".' call dot_plot ( file_out_name, dim_num, point_num, & coord, coord_min, coord_max, x_index, y_index ) ! ! KMEANS ! Create a K-Means plot. ! We presume the cluster centers are each set off with blank records. ! We presume the data was tagged with a third cluster coordinate. ! else if ( command(1:1) == 'K' ) then tag(1:point_num) = nint ( coord(3,1:point_num) ) call file_name_inc ( file_out_name ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Create a K-Means plot "' // & trim ( file_out_name ) // '".' call kmeans_plot ( file_out_name, dim_num, point_num, coord, & coord_min, coord_max, x_index, y_index, connect, tag ) ! ! Q = Quit ! else if ( command(1:1) == 'Q' ) then exit ! ! TH = Thin ! else if ( command ( 1:2) == 'TH' ) then tol = 0.01E+00 * maxval ( coord_max(1:dim_num) - coord_min(1:dim_num) ) call points_thin ( connect, dim_num, point_num, coord, tol, point_num2 ) point_num = point_num2 ! ! TV ! Create a triangulated Voronoi plot. ! else if ( command(1:2) == 'TV' ) then call file_name_inc ( file_out_name ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Create triangulated Voronoi plot "' & // trim ( file_out_name ) // '".' call trivor_plot ( file_out_name, dim_num, point_num, coord, & coord_min, coord_max, x_index, y_index ) ! ! VORONOI ! Create a Voronoi plot. ! else if ( command(1:3) == 'VOR' ) then call file_name_inc ( file_out_name ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Create Voronoi plot "' // & trim ( file_out_name ) // '".' call voronoi_plot ( file_out_name, dim_num, point_num, coord, & coord_min, coord_max, x_index, y_index ) ! ! X_INDEX ! else if ( command(1:1) == 'X' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter new value for X index:' read ( *, * ) x_index write ( *, '(a,i6)' ) 'X index set to ', x_index ! ! Y_INDEX ! else if ( command(1:2) == 'Y' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter new value for Y index:' read ( *, * ) y_index write ( *, '(a,i6)' ) 'Y index set to ', y_index ! ! Unrecognized command. ! else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Unrecognized command: "'// trim ( command ) // '"' end if end do ! ! Deallocate arrays. ! deallocate ( connect ) deallocate ( coord ) deallocate ( coord_min ) deallocate ( coord_max ) deallocate ( tag ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_POINTS' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine angle_contains_ray_2d ( inside, x1, y1, x2, y2, x3, y3, x, y ) ! !******************************************************************************* ! !! ANGLE_CONTAINS_RAY_2D determines if an angle contains a ray, in 2D. ! ! ! Discussion: ! ! The angle is defined by the sequence of points (X1,Y1), (X2,Y2) ! and (X3,Y3). ! ! The ray is defined by the sequence of points (X2,Y2), (X,Y). ! ! Modified: ! ! 17 March 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X1, Y1, X2, Y2, X3, Y3, the X and Y coordinates of ! the angle. ! ! Input, real X, Y, the end point of the ray to be checked. ! The ray is assumed to have an origin at (X2,Y2). ! ! Output, logical INSIDE, is .TRUE. if the ray is inside ! the angle or on its boundary, and .FALSE. otherwise. ! implicit none ! real a1 real a2 real angle_deg_2d logical inside real x real x1 real x2 real x3 real y real y1 real y2 real y3 ! a1 = angle_deg_2d ( x1, y1, x2, y2, x, y ) a2 = angle_deg_2d ( x1, y1, x2, y2, x3, y3 ) if ( a1 <= a2 ) then inside = .true. else inside = .false. end if return end function angle_deg_2d ( x1, y1, x2, y2, x3, y3 ) ! !******************************************************************************* ! !! ANGLE_DEG_2D returns the angle swept out between two rays in 2D. ! ! ! Discussion: ! ! Except for the zero angle case, it should be true that ! ! ANGLE_DEG_2D(X1,Y1,X2,Y2,X3,Y3) ! + ANGLE_DEG_2D(X3,Y3,X2,Y2,X1,Y1) = 360.0 ! ! Modified: ! ! 14 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X1, Y1, X2, Y2, X3, Y3, define the rays ! ( X1-X2, Y1-Y2 ) and ( X3-X2, Y3-Y2 ) which in turn define the ! angle, counterclockwise from ( X1-X2, Y1-Y2 ). ! ! Output, real ANGLE_DEG_2D, the angle swept out by the rays, measured ! in degrees. 0 <= ANGLE_DEG_2D < 360. If either ray has zero length, ! then ANGLE_DEG_2D is set to 0. ! implicit none ! real angle_deg_2d real angle_rad_2d real pi real radians_to_degrees real x real x1 real x2 real x3 real y real y1 real y2 real y3 ! x = ( x1 - x2 ) * ( x3 - x2 ) + ( y1 - y2 ) * ( y3 - y2 ) y = ( x1 - x2 ) * ( y3 - y2 ) - ( y1 - y2 ) * ( x3 - x2 ) if ( x == 0.0E+00 .and. y == 0.0E+00 ) then angle_deg_2d = 0.0E+00 else angle_rad_2d = atan2 ( y, x ) if ( angle_rad_2d < 0.0E+00 ) then angle_rad_2d = angle_rad_2d + 2.0E+00 * pi ( ) end if angle_deg_2d = radians_to_degrees ( angle_rad_2d ) end if return end subroutine i_to_rgb ( i, ilo, ihi, r, g, b ) ! !******************************************************************************* ! !! I_TO_RGB returns a color on the perimeter of the RGB color hexagon. ! ! ! Modified: ! ! 12 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, indicates the color to select. ! ! Input, integer ILO, IHI, the lowest and highest values of I. ! ! Output, real R, G, B, RGB specifications for the color that lies ! at the angle specified by I, ILO and IHI, on the perimeter of the ! color hexagon. One value will be 1, and one value will be 0. ! implicit none ! real angle real b real g integer i integer ihi integer ilo real r ! angle = 360.0 * real ( i - ilo ) / real ( ihi + 1 - ilo ) angle = mod ( angle, 360.0E+00 ) if ( angle < 0.0E+00 ) then angle = angle + 360.0E+00 end if if ( angle <= 60.0E+00 ) then r = 1.0E+00 g = angle / 60.0E+00 b = 0.0E+00 else if ( angle <= 120.0E+00 ) then r = ( 120.0E+00 - angle ) / 60.0E+00 g = 1.0E+00 b = 0.0E+00 else if ( angle <= 180.0E+00 ) then r = 0.0E+00 g = 1.0E+00 b = ( angle - 120.0E+00 ) / 60.0E+00 else if ( angle <= 240.0E+00 ) then r = 0.0E+00 g = ( 240.0E+00 - angle ) / 60.0E+00 b = 1.0E+00 else if ( angle <= 300.0E+00 ) then r = ( angle - 240.0E+00 ) / 60.0E+00 g = 0.0E+00 b = 1.0E+00 else if ( angle <= 360.0E+00 ) then r = 1.0E+00 g = 0.0E+00 b = ( 360.0E+00 - angle ) / 60.0E+00 end if return end subroutine box_ray_int_2d ( xmin, ymin, xmax, ymax, xa, ya, xb, yb, xi, yi ) ! !******************************************************************************* ! !! BOX_RAY_INT_2D: intersection ( box, ray ) in 2D. ! ! ! Discussion: ! ! The box is assumed to be a rectangle with sides aligned on coordinate ! axes. ! ! The origin of the ray is assumed to be inside the box. This ! guarantees that the ray will intersect the box in exactly one point. ! ! Modified: ! ! 16 June 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real XMIN, YMIN, the lower left corner of the box. ! ! Input, real XMAX, YMAX, the upper right corner of the box. ! ! Input, real XA, YA, the origin of the ray, which should be ! inside the box. ! ! Input, real XB, YB, a second point on the ray. ! ! Output, real XI, YI, the point on the box intersected by the ray. ! implicit none ! logical inside integer ival integer side real xa real xb real xc real xd real xi real xmax real xmin real ya real yb real yc real yd real yi real ymax real ymin ! do side = 1, 4 if ( side == 1 ) then xc = xmin yc = ymin xd = xmax yd = ymin else if ( side == 2 ) then xc = xmax yc = ymin xd = xmax yd = ymax else if ( side == 3 ) then xc = xmax yc = ymax xd = xmin yd = ymax else if ( side == 4 ) then xc = xmin yc = ymax xd = xmin yd = ymin end if call angle_contains_ray_2d ( inside, xc, yc, xa, ya, xd, yd, xb, yb ) if ( inside ) then exit end if if ( side == 4 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BOX_RAY_INT_2D - Fatal error!' write ( *, '(a)' ) ' No intersection could be found.' stop end if end do call lines_exp_int_2d ( xa, ya, xb, yb, xc, yc, xd, yd, ival, xi, yi ) return end subroutine ch_cap ( c ) ! !******************************************************************************* ! !! CH_CAP capitalizes a single character. ! ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none ! character c integer itemp ! itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end function ch_eqi ( c1, c2 ) ! !******************************************************************************* ! !! CH_EQI is a case insensitive comparison of two characters for equality. ! ! ! Examples: ! ! CH_EQI ( 'A', 'a' ) is .TRUE. ! ! 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 function ch_is_digit ( c ) ! !******************************************************************************* ! !! CH_IS_DIGIT returns .TRUE. if a character is a decimal digit. ! ! ! Modified: ! ! 15 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, the character to be analyzed. ! ! Output, logical CH_IS_DIGIT, .TRUE. if C is a digit, .FALSE. otherwise. ! implicit none ! character c logical ch_is_digit ! if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then ch_is_digit = .true. else ch_is_digit = .false. end if return end subroutine ch_to_digit ( c, digit ) ! !******************************************************************************* ! !! CH_TO_DIGIT returns the integer value of a base 10 digit. ! ! ! Example: ! ! C DIGIT ! --- ----- ! '0' 0 ! '1' 1 ! ... ... ! '9' 9 ! ' ' 0 ! 'X' -1 ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, the decimal digit, '0' through '9' or blank ! are legal. ! ! Output, integer DIGIT, the corresponding integer value. If C was ! 'illegal', then DIGIT is -1. ! implicit none ! character c integer 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 subroutine dash_plot ( connect, file_name, dim_num, point_num, coord, & coord_min, coord_max, x, y ) ! !******************************************************************************* ! !! DASH_PLOT plots a set of points, and connects them. ! ! ! Modified: ! ! 06 March 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer CONNECT(POINT_NUM), connection indicator. ! If CONNECT(I) is: ! 0, it is an isolated point; ! 1, it is connected to the previous point only; ! 2, it is connected to the next point only; ! 3, it is connected to the previous and next points. ! ! Input, character ( len = * ) FILE_NAME, the name of the dot file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer POINT_NUM, the number of points. ! ! Input, real COORD(DIM_NUM,POINT_NUM), the point coordinates. ! ! Input, real COORD_MIN(DIM_NUM), COORD_MAX(DIM_NUM), the data ranges. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! implicit none ! integer dim_num integer point_num ! integer connect(point_num) real coord(dim_num,point_num) real coord_max(dim_num) real coord_min(dim_num) character ( len = * ) file_name integer file_unit integer i integer ierror integer x real xvec(4) integer y real yvec(4) ! call get_unit ( file_unit ) call ps_file_open ( file_name, file_unit, ierror ) call eps_file_head ( file_name ) call ps_page_head ( coord_min(x), coord_min(y), coord_max(x), coord_max(y) ) ! ! We define a PostScript clipping box. ! xvec(1:4) = (/ coord_min(x), coord_max(x), coord_max(x), coord_min(x) /) yvec(1:4) = (/ coord_min(y), coord_min(y), coord_max(y), coord_max(y) /) call ps_clip ( 4, xvec, yvec ) if ( point_num <= 350 ) then do i = 1, point_num call ps_mark_disk ( coord(x,i), coord(y,i) ) end do else do i = 1, point_num call ps_mark_point ( coord(x,i), coord(y,i) ) end do end if do i = 1, point_num-1 if ( connect(i) == 2 .or. connect(i) == 3 ) then call ps_line ( coord(x,i), coord(y,i), coord(x,i+1), coord(y,i+1) ) end if end do if ( connect(point_num) == 2 .or. connect(point_num) == 3 ) then call ps_line ( coord(x,point_num), coord(y,point_num), coord(x,1), & coord(y,1) ) end if call ps_page_tail call eps_file_tail call ps_file_close ( file_unit ) return end subroutine delaunay_plot ( file_name, dim_num, point_num, coord, coord_min, & coord_max, x, y ) ! !******************************************************************************* ! !! DELAUNAY_PLOT plots the Delaunay triangulation of a pointset. ! ! ! Modified: ! ! 26 April 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the output file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer POINT_NUM, the number of points. ! ! Input, real COORD(DIM_NUM,POINT_NUM), the point coordinates. ! ! Input, real COORD_MIN(DIM_NUM), COORD_MAX(DIM_NUM), the data ranges. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! implicit none ! integer dim_num integer point_num ! real coord(dim_num,point_num) real coord_max(dim_num) real coord_min(dim_num) real coord2(2,point_num) logical, parameter :: debug = .false. character ( len = * ) file_name integer file_unit integer i integer ierror integer ind(point_num) integer j1 integer j2 integer j3 integer nodtri(3,2*point_num) integer tri_num integer stack(point_num) integer tnbr(3,2*point_num) integer x real xvec(4) real xx integer y real yvec(4) real yy ! ! Compute the Delaunay triangulation. ! coord2(1,1:point_num) = coord(x,1:point_num) coord2(2,1:point_num) = coord(y,1:point_num) call ivec_identity ( point_num, ind ) call rtris2 ( point_num, point_num, coord2, ind, tri_num, nodtri, tnbr, stack, & ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DELAUNAY_PLOT - Fatal error!' write ( *, '(a,i6)' ) ' RTRIS2 returns IERROR = ', ierror return end if ! ! Print the Delaunay triangulation. ! if ( debug ) then call delaunay_print ( point_num, coord2, tri_num, nodtri, tnbr ) end if ! ! Plot the Delaunay triangulation. ! call get_unit ( file_unit ) call ps_file_open ( file_name, file_unit, ierror ) call eps_file_head ( file_name ) call ps_page_head ( coord_min(x), coord_min(y), coord_max(x), coord_max(y) ) ! ! We define a PostScript clipping box. ! xvec(1:4) = (/ coord_min(x), coord_max(x), coord_max(x), coord_min(x) /) yvec(1:4) = (/ coord_min(y), coord_min(y), coord_max(y), coord_max(y) /) call ps_clip ( 4, xvec, yvec ) xx = coord_min(x) yy = coord_min(y) call ps_moveto ( xx, yy ) ! call ps_label ( 'Delaunay triangulation' ) if ( point_num <= 350 ) then do i = 1, point_num call ps_mark_disk ( coord(x,i), coord(y,i) ) end do else do i = 1, point_num call ps_mark_point ( coord(x,i), coord(y,i) ) end do end if do i = 1, tri_num j1 = nodtri(1,i) j2 = nodtri(2,i) j3 = nodtri(3,i) call ps_triangle ( coord2(1,j1), coord2(2,j1), coord2(1,j2), & coord2(2,j2), coord2(1,j3), coord2(2,j3) ) end do call ps_page_tail call eps_file_tail call ps_file_close ( file_unit ) return end subroutine delaunay_print ( point_num, coord, tri_num, nodtri, tnbr ) ! !******************************************************************************* ! !! DELAUNAY_PRINT prints out information defining a Delaunay triangulation. ! ! ! Modified: ! ! 08 July 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer POINT_NUM, the number of points. ! ! Input, real COORD(2,POINT_NUM), the point coordinates. ! ! Input, integer TRI_NUM, the number of triangles. ! ! Input, integer NODTRI(3,TRI_NUM), the nodes that make up the triangles. ! ! Input, integer TNBR(3,TRI_NUM), the triangle neighbors on each side. ! implicit none ! integer tri_num integer point_num ! real coord(2,point_num) integer i integer i_wrap integer j integer k integer n1 integer n2 integer nodtri(3,tri_num) integer s integer t integer tnbr(3,tri_num) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DELAUNAY_PRINT' write ( *, '(a)' ) ' Information defining a Delaunay triangulation.' write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' The number of points is ', point_num call rmat_print ( point_num, point_num, 2, transpose ( coord ), & ' Point coordinates (transpose of internal array)' ) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' The number of triangles is ', tri_num write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Sets of three points are used as vertices of' write ( *, '(a)' ) ' the triangles. For each triangle, the points' write ( *, '(a)' ) ' are listed in counterclockwise order.' call imat_print ( tri_num, tri_num, 3, transpose ( nodtri ), & ' Nodes that make up triangles (transpose of internal array)' ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' On each side of a given triangle, there is either' write ( *, '(a)' ) ' another triangle, or a piece of the convex hull.' write ( *, '(a)' ) ' For each triangle, we list the indices of the three' write ( *, '(a)' ) ' neighbors, or (if negative) the codes of the' write ( *, '(a)' ) ' segments of the convex hull.' call imat_print ( tri_num, tri_num, 3, transpose ( tnbr ), & ' Indices of neighboring triangles (transpose of internal array)' ) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' The number of boundary points (and segments) is ', & 2 * point_num - tri_num - 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The segments that make up the convex hull can be' write ( *, '(a)' ) ' determined from the negative entries of the triangle' write ( *, '(a)' ) ' neighbor list.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' # Tri Side N1 N2' write ( *, '(a)' ) ' ' k = 0 do i = 1, tri_num do j = 1, 3 if ( tnbr(j,i) < 0 ) then s = - tnbr(j,i) t = s / 3 s = mod ( s, 3 ) + 1 k = k + 1 n1 = nodtri(s,t) n2 = nodtri(i_wrap(s+1,1,3),t) write ( *, '(5i4)' ) k, t, s, n1, n2 end if end do end do return end function diaedg ( x0, y0, x1, y1, x2, y2, x3, y3 ) ! !******************************************************************************* ! !! DIAEDG chooses a diagonal edge. ! ! ! 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. ! ! Modified: ! ! 19 February 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Parameters: ! ! Input, real X0, Y0, X1, Y1, X2, Y2, X3, Y3, the ! coordinates of the vertices of a quadrilateral, given in ! counter clockwise order. ! ! Output, integer DIAEDG, 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 ca real cb integer diaedg real dx10 real dx12 real dx30 real dx32 real dy10 real dy12 real dy30 real dy32 real s real tol real tola real tolb real x0 real x1 real x2 real x3 real y0 real y1 real y2 real y3 ! tol = 100.0E+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 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 ( ca > tola .and. cb > tolb ) then diaedg = -1 else if ( ca < -tola .and. cb < -tolb ) then diaedg = 1 else tola = max ( tola, tolb ) s = ( dx10 * dy30 - dx30 * dy10 ) * cb + ( dx32 * dy12 - dx12 * dy32 ) * ca if ( s > tola ) then diaedg = -1 else if ( s < -tola ) then diaedg = 1 else diaedg = 0 end if end if return end subroutine digit_inc ( c ) ! !******************************************************************************* ! !! DIGIT_INC increments a decimal digit. ! ! ! Example: ! ! Input Output ! ----- ------ ! '0' '1' ! '1' '2' ! ... ! '8' '9' ! '9' '0' ! 'A' 'A' ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, a digit to be incremented. ! implicit none ! character c integer digit call ch_to_digit ( c, digit ) if ( digit == -1 ) then return end if digit = digit + 1 if ( digit == 10 ) then digit = 0 end if call digit_to_ch ( digit, c ) return end subroutine digit_to_ch ( digit, c ) ! !******************************************************************************* ! !! DIGIT_TO_CH returns the character representation of a decimal digit. ! ! ! Example: ! ! DIGIT C ! ----- --- ! 0 '0' ! 1 '1' ! ... ... ! 9 '9' ! 17 '*' ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer DIGIT, the digit value between 0 and 9. ! ! Output, character C, the corresponding character, or '*' if DIGIT ! was illegal. ! implicit none ! character c integer digit ! if ( 0 <= digit .and. digit <= 9 ) then c = char ( digit + 48 ) else c = '*' end if return end subroutine dot_plot ( filedot_name, dim_num, point_num, coord, & coord_min, coord_max, x, y ) ! !******************************************************************************* ! !! DOT_PLOT plots the individual points. ! ! ! Modified: ! ! 24 April 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILEDOT_NAME, the name of the dot file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer POINT_NUM, the number of points. ! ! Input, real COORD(DIM_NUM,POINT_NUM), the point coordinates. ! ! Input, real COORD_MIN(DIM_NUM), COORD_MAX(DIM_NUM), the data ranges. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! implicit none ! integer dim_num integer point_num ! real coord(dim_num,point_num) real coord_max(dim_num) real coord_min(dim_num) character ( len = * ) filedot_name integer filedot_unit integer i integer ierror integer x real xvec(4) integer y real yvec(4) ! call get_unit ( filedot_unit ) call ps_file_open ( filedot_name, filedot_unit, ierror ) call eps_file_head ( filedot_name ) call ps_page_head ( coord_min(x), coord_min(y), coord_max(x), coord_max(y) ) ! ! We define a PostScript clipping box. ! xvec(1:4) = (/ coord_min(x), coord_max(x), coord_max(x), coord_min(x) /) yvec(1:4) = (/ coord_min(y), coord_min(y), coord_max(y), coord_max(y) /) call ps_clip ( 4, xvec, yvec ) if ( point_num <= 350 ) then do i = 1, point_num call ps_mark_disk ( coord(x,i), coord(y,i) ) end do else do i = 1, point_num call ps_mark_point ( coord(x,i), coord(y,i) ) end do end if call ps_page_tail call eps_file_tail call ps_file_close ( filedot_unit ) return end function enormsq0_2d ( x0, y0, x1, y1 ) ! !******************************************************************************* ! !! ENORMSQ0_2D computes the square of the Euclidean norm of (P1-P0) in 2D. ! ! ! Modified: ! ! 06 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X0, Y0, X1, Y1, the coordinates of the points ! P0 and P1. ! ! Output, real ENORMSQ0_2D, the square of the Euclidean norm of (P1-P0). ! implicit none ! real enormsq0_2d real x0 real x1 real y0 real y1 ! enormsq0_2d = ( x1 - x0 )**2 + ( y1 - y0 )**2 return end subroutine file_column_count ( file_name, ncolumn ) ! !******************************************************************************* ! !! 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 NCOLUMN 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. ! ! Modified: ! ! 21 June 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the file. ! ! Output, integer NCOLUMN, the number of columns assumed to be in the file. ! implicit none ! character ( len = * ) file_name logical got_one integer ios integer iunit character ( len = 256 ) line integer ncolumn ! ! Open the file. ! call get_unit ( iunit ) open ( unit = iunit, file = file_name, status = 'old', form = 'formatted', & access = 'sequential', iostat = ios ) if ( ios /= 0 ) then ncolumn = - 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Fatal error!' write ( *, '(a)' ) ' Could not open the file:' write ( *, '(a)' ) ' ' // trim ( file_name ) return end if ! ! Read one line, but skip blank lines and comment lines. ! got_one = .false. do read ( iunit, '(a)', iostat = ios ) line if ( ios /= 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 ( iunit ) do read ( iunit, '(a)', iostat = ios ) line if ( ios /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if got_one = .true. exit end do end if close ( unit = iunit ) if ( .not. got_one ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Warning!' write ( *, '(a)' ) ' The file does not seem to contain any data.' ncolumn = 0 return end if call word_count ( line, ncolumn ) return end subroutine file_column_range ( file_name, ncolumn, col_min, col_max ) ! !******************************************************************************* ! !! FILE_COLUMN_RANGE determines the minimum and maximum ranges of each column. ! ! ! Discussion: ! ! The file is assumed to be a simple text file. ! ! Each line of the file is presumed to consist of NCOLUMN words, separated ! by spaces. ! ! The routine computes the range of each column. ! ! Modified: ! ! 19 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the file. ! ! Input, integer NCOLUMN, the number of columns assumed to be in the file. ! ! Output, real COL_MIN(NCOLUM), COL_MAX(NCOLUMN), the minimum and maximum ! for each column. ! implicit none ! integer ncolumn ! real col_max(ncolumn) real col_min(ncolumn) character ( len = * ) file_name integer ierror integer ios integer iunit integer j character ( len = 256 ) line integer nrow real x(ncolumn) ! ! Open the file. ! call get_unit ( iunit ) open ( unit = iunit, file = file_name, status = 'old', form = 'formatted', & access = 'sequential', iostat = ios ) if ( ios /= 0 ) then ncolumn = - 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_RANGE - Fatal error!' write ( *, '(a)' ) ' Could not open the file:' write ( *, '(a)' ) ' ' // trim ( file_name ) return end if nrow = 0 col_min(1:ncolumn) = 0.0E+00 col_max(1:ncolumn) = 0.0E+00 do read ( iunit, '(a)', iostat = ios ) line if ( ios /= 0 ) then exit end if call s_to_rvec ( line, ncolumn, x, ierror ) if ( ierror /= 0 ) then exit end if nrow = nrow + 1 if ( nrow == 1 ) then col_min(1:ncolumn) = x(1:ncolumn) col_max(1:ncolumn) = x(1:ncolumn) else do j = 1, ncolumn col_min(j) = min ( col_min(j), x(j) ) col_max(j) = max ( col_max(j), x(j) ) end do end if end do close ( unit = iunit ) return end subroutine file_name_ext_get ( file_name, i, j ) ! !******************************************************************************* ! !! FILE_NAME_EXT_GET determines the "extension" of a file name. ! ! ! Definition: ! ! 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. ! ! Note: ! ! 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. ! ! Examples: ! ! FILE_NAME I J ! ! bob.for 4 7 ! N.B.C.D 6 7 ! Naomi. 6 6 ! Arthur 0 0 ! .com 1 1 ! ! Modified: ! ! 17 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, a file name to be examined. ! ! Output, integer 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 i integer j integer 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 ) ! !******************************************************************************* ! !! FILE_NAME_EXT_SWAP replaces the current "extension" of a file name. ! ! ! Definition: ! ! 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. ! ! Examples: ! ! Input Output ! ================ ========= ! FILE_NAME EXT FILE_NAME ! ! bob.for obj bob.obj ! bob.bob.bob txt bob.bob.txt ! bob yak bob.yak ! ! 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 i integer j integer len_max integer 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_name + 1 > len_max ) 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_name_inc ( file_name ) ! !******************************************************************************* ! !! FILE_NAME_INC generates the next filename in a series. ! ! ! Discussion: ! ! It is assumed that the digits in the name, whether scattered or ! connected, represent a number that is to be increased by 1 on ! each call. If this number is all 9's on input, the output number ! is all 0's. Non-numeric letters of the name are unaffected, and ! if the name contains no digits, then nothing is done. ! ! Examples: ! ! Input Output ! ----- ------ ! a7to11.txt a7to12.txt ! a7to99.txt a8to00.txt ! a9to99.txt a0to00.txt ! cat.txt cat.txt ! ! Modified: ! ! 09 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) FILE_NAME. ! On input, a character string to be incremented. ! On output, the incremented string. ! implicit none ! character c logical ch_is_digit character ( len = * ) file_name integer i integer lens ! lens = len_trim ( file_name ) do i = lens, 1, -1 c = file_name(i:i) if ( ch_is_digit ( c ) ) then call digit_inc ( c ) file_name(i:i) = c if ( c /= '0' ) then return end if end if end do return end subroutine get_unit ( iunit ) ! !******************************************************************************* ! !! GET_UNIT returns a free FORTRAN unit number. ! ! ! Discussion: ! ! A "free" FORTRAN unit number is an integer 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. ! ! Modified: ! ! 02 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IUNIT. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5 and 6). ! ! Otherwise, IUNIT is an integer 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. ! implicit none ! integer i integer ios integer iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 ) 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 i_modp ( i, j ) ! !******************************************************************************* ! !! I_MODP returns the nonnegative remainder of integer division. ! ! ! Formula: ! ! If ! NREM = I_MODP ( I, J ) ! NMULT = ( I - NREM ) / J ! then ! I = J * NMULT + NREM ! where NREM is always nonnegative. ! ! Comments: ! ! 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, I_MODP(A,360) is between 0 and 360, always. ! ! Examples: ! ! I J MOD I_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 ! ! Modified: ! ! 02 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the number to be divided. ! ! Input, integer J, the number that divides I. ! ! Output, integer I_MODP, the nonnegative remainder when I is ! divided by J. ! implicit none ! integer i integer i_modp integer j ! if ( j == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I_MODP - Fatal error!' write ( *, '(a,i6)' ) ' I_MODP ( I, J ) called with J = ', j stop end if i_modp = mod ( i, j ) if ( i_modp < 0 ) then i_modp = i_modp + abs ( j ) end if return end subroutine i_swap ( i, j ) ! !******************************************************************************* ! !! I_SWAP swaps two integer values. ! ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer I, J. On output, the values of I and ! J have been interchanged. ! implicit none ! integer i integer j integer k ! k = i i = j j = k return end function i_wrap ( ival, ilo, ihi ) ! !******************************************************************************* ! !! I_WRAP forces an integer to lie between given limits by wrapping. ! ! ! Example: ! ! ILO = 4, IHI = 8 ! ! I I_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 ! ! Modified: ! ! 15 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IVAL, an integer value. ! ! Input, integer ILO, IHI, the desired bounds for the integer value. ! ! Output, integer I_WRAP, a "wrapped" version of IVAL. ! implicit none ! integer i_modp integer i_wrap integer ihi integer ilo integer ival integer wide ! wide = ihi + 1 - ilo if ( wide == 0 ) then i_wrap = ilo else i_wrap = ilo + i_modp ( ival-ilo, wide ) end if return end subroutine imat_print ( lda, m, n, a, title ) ! !******************************************************************************* ! !! IMAT_PRINT prints an integer matrix. ! ! ! Modified: ! ! 04 September 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer M, the number of rows in A. ! ! Input, integer N, the number of columns in A. ! ! Input, integer A(LDA,N), the matrix to be printed. ! ! Input, character ( len = * ) TITLE, a title to be printed first. ! TITLE may be blank. ! implicit none ! integer lda integer n ! integer a(lda,n) integer i integer j integer jhi integer jlo integer m character ( len = * ) title ! if ( title /= ' ' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do jlo = 1, n, 10 jhi = min ( jlo + 9, n ) write ( *, '(a)' ) ' ' write ( *, '(6x,10(i7))' ) ( j, j = jlo, jhi ) write ( *, '(a)' ) ' ' do i = 1, m write ( *, '(i6,10i7)' ) i, a(i,jlo:jhi) end do end do return end subroutine ivec_identity ( n, a ) ! !******************************************************************************* ! !! IVEC_IDENTITY sets an integer vector to the identity vector A(I)=I. ! ! ! Modified: ! ! 09 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of elements of A. ! ! Output, integer A(N), the array to be initialized. ! implicit none ! integer n ! integer a(n) integer i ! do i = 1, n a(i) = i end do return end subroutine kmeans_plot ( file_name, dim_num, point_num, coord, coord_min, & coord_max, x, y, connect, tag ) ! !******************************************************************************* ! !! KMEANS_PLOT plots a K-Means diagram. ! ! ! Discussion: ! ! The "unconnected" points, which have CONNECT(I) == 0, are ! assumed to represent the centers. ! ! The "connected" points, which have CONNECT(I) /= 0, are assumed ! to be the data values to be clustered. ! ! A Voronoi diagram is made of the centers. ! ! Then the connected points are displayed, with colors corresponding ! to their assigned center. ! ! Modified: ! ! 08 March 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the output file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer POINT_NUM, the number of points. ! ! Input, real COORD(DIM_NUM,POINT_NUM), the point coordinates. ! ! Input, real COORD_MIN(DIM_NUM), COORD_MAX(DIM_NUM), the data ranges. ! These values may be set to the minimum and maximum values of the data, ! or can be used to clip the data, or to put some whitespace around it. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! This lets COOR_MIN and COORD_MAX have a dimension larger than 2. ! We simply specify which pair of indices will play the X and Y roles. ! ! Input, integer CONNECT(POINT_NUM), connection indicator. ! If CONNECT(I) is: ! 0, it is an isolated point; ! 1, it is connected to the previous point only; ! 2, it is connected to the next point only; ! 3, it is connected to the previous and next points. ! ! Input, integer TAG(POINT_NUM), a tag for each point, ! presumably the cluster index. ! implicit none ! integer dim_num integer point_num ! real b real center(dim_num,2*point_num) integer center_num integer connect(point_num) real coord(dim_num,point_num) real coord_max(dim_num) real coord_min(dim_num) real coord2(2,point_num) logical, parameter :: debug = .false. character ( len = * ) file_name integer file_unit real g integer i integer i_wrap integer i1 integer i2 integer i3 integer ierror integer ind(point_num) logical inside integer ix integer j integer jp1 integer jp2 integer k real n1 real n2 integer nodtri(3,2*point_num) real r integer tag(point_num) integer tag2(point_num) integer tag_lo integer tag_hi integer tri_num integer stack(point_num) integer tnbr(3,2*point_num) integer x real xvec(4) real x1 real x2 real x3 real x4 real x5 real x6 integer y real yvec(4) real y1 real y2 real y3 real y4 real y5 real y6 ! ! Determine which points represent the centers. ! center_num = 0 if ( debug ) then write ( *, * ) ' ' write ( *, * ) 'CENTERS:' write ( *, * ) ' ' end if do i = 1, point_num if ( connect(i) == 0 ) then center_num = center_num + 1 coord2(1,center_num) = coord(x,i) coord2(2,center_num) = coord(y,i) tag2(center_num) = tag(i) if ( debug ) then write ( *, '(3i4,2g14.6)' ) center_num, i, tag(i), coord(x,i), coord(y,i) end if end if end do ! ! Compute the Delaunay triangulation of the centers. ! call ivec_identity ( center_num, ind ) call rtris2 ( center_num, center_num, coord2, ind, tri_num, nodtri, tnbr, & stack, ierror ) ! ! Compute the intersection point of the perpendicular bisectors ! of each Delaunay triangle. ! do j = 1, tri_num i1 = nodtri(1,j) i2 = nodtri(2,j) i3 = nodtri(3,j) call triangle_circumcenter_2d ( coord2(1,i1), coord2(2,i1), coord2(1,i2), & coord2(2,i2), coord2(1,i3), coord2(2,i3), center(1,j), center(2,j) ) end do ! ! Plot the Voronoi tessellation. ! call get_unit ( file_unit ) call ps_file_open ( file_name, file_unit, ierror ) call eps_file_head ( file_name ) call ps_page_head ( coord_min(x), coord_min(y), coord_max(x), coord_max(y) ) ! ! We need to define a PostScript clipping box, so that we can let PostScript ! take care of clipping lines that lie outside the box. ! xvec(1:4) = (/ coord_min(x), coord_max(x), coord_max(x), coord_min(x) /) yvec(1:4) = (/ coord_min(y), coord_min(y), coord_max(y), coord_max(y) /) call ps_clip ( 4, xvec, yvec ) x1 = coord_min(x) y1 = coord_min(y) call ps_moveto ( x1, y1 ) ! call ps_label ( 'K-Means diagram' ) call ps_setting_int ( 'SET', 'MARKER_SIZE', 5 ) tag_lo = 1 tag_hi = center_num do i = 1, center_num call i_to_rgb ( tag2(i), tag_lo, tag_hi, r, g, b ) call ps_color_fill_set ( r, g, b ) call ps_mark_disk ( coord2(1,i), coord2(2,i) ) end do call ps_setting_int ( 'SET', 'MARKER_SIZE', 3 ) do i = 1, point_num if ( connect(i) /= 0 ) then call i_to_rgb ( tag(i), tag_lo, tag_hi, r, g, b ) call ps_color_line_set ( r, g, b ) call ps_mark_circle ( coord(x,i), coord(y,i) ) end if end do r = 0.1E+00 g = 1.0E+00 b = 0.1E+00 call ps_color_line_set ( r, g, b ) ! ! For each Delaunay triangle, I ! For each side J, ! do i = 1, tri_num do j = 1, 3 k = tnbr(j,i) ! ! If there is a neighboring triangle K on that side, ! connect the circumcenters. ! if ( k > 0 ) then if ( k > i ) then call ps_line ( center(1,i), center(2,i), center(1,k), center(2,k) ) end if ! ! If there is no neighboring triangle on that side, ! extend a line from the circumcenter of I in the direction of the ! outward normal to that side. ! else ix = nodtri(j,i) x1 = coord2(1,ix) y1 = coord2(2,ix) jp1 = i_wrap ( j+1, 1, 3 ) ix = nodtri(jp1,i) x2 = coord2(1,ix) y2 = coord2(2,ix) jp2 = i_wrap ( j+2, 1, 3 ) ix = nodtri(jp2,i) x3 = coord2(1,ix) y3 = coord2(2,ix) x4 = center(1,i) y4 = center(2,i) call line_exp_normal_2d ( x1, y1, x2, y2, n1, n2 ) x5 = x4 + n1 y5 = y4 + n2 call box_ray_int_2d ( coord_min(x), coord_min(y), coord_max(x), & coord_max(y), x4, y4, x5, y5, x6, y6 ) call ps_line ( x4, y4, x6, y6 ) end if end do end do call ps_color_line ( 'POP', r, g, b ) call ps_page_tail call eps_file_tail call ps_file_close ( file_unit ) return end subroutine line_exp2imp_2d ( x1, y1, x2, y2, a, b, c ) ! !******************************************************************************* ! !! LINE_EXP2IMP_2D converts an explicit line to implicit form in 2D. ! ! ! Formula: ! ! The explicit form of a line in 2D is: ! ! (X1,Y1), (X2,Y2). ! ! The implicit form of a line in 2D is: ! ! A * X + B * Y + C = 0 ! ! Modified: ! ! 24 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X1, Y1, X2, Y2. (X1,Y1) and (X2,Y2) are ! two points on the line. (X1,Y1) must be different ! from (X2,Y2). ! ! Output, real A, B, C, three coefficients which describe ! the line that passes through (X1,Y1) and (X2,Y2). ! implicit none ! real a real b real c real x1 real x2 real y1 real y2 ! ! Take care of degenerate cases. ! if ( x1 == x2 .and. y1 == y2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LINE_EXP2IMP_2D - Fatal error!' write ( *, '(a)' ) ' (X1,Y1) = (X2,Y2)' write ( *, '(a,2g14.6)' ) ' (X1,Y1) = ', x1, y1 write ( *, '(a,2g14.6)' ) ' (X2,Y2) = ', x2, y2 stop end if a = y2 - y1 b = x1 - x2 c = x2 * y1 - x1 * y2 return end subroutine line_exp_normal_2d ( x1, y1, x2, y2, n1, n2 ) ! !******************************************************************************* ! !! LINE_EXP_NORMAL_2D computes the normal to a line in 2D. ! ! ! Formula: ! ! The explicit form of a line in 2D is: ! ! (X1,Y1), (X2,Y2). ! ! Modified: ! ! 17 June 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X1, Y1, X2, Y2, two points on the line. ! ! Output, real N1, N2, the components of a unit normal vector to the line. ! If the two points are equal, then N1 = N2 = 0. ! implicit none ! real n1 real n2 real norm real x1 real x2 real y1 real y2 ! norm = sqrt ( ( x2 - x1 )**2 + ( y2 - y1 )**2 ) if ( norm == 0.0E+00 ) then n1 = 0.0E+00 n2 = 0.0E+00 return end if n1 = ( y2 - y1 ) / norm n2 = - ( x2 - x1 ) / norm return end subroutine lines_exp_int_2d ( x1, y1, x2, y2, x3, y3, x4, y4, ival, x, y ) ! !******************************************************************************* ! !! LINES_EXP_INT_2D determines where two explicit lines intersect in 2D. ! ! ! Formula: ! ! The explicit form of a line in 2D is: ! ! (X1,Y1), (X2,Y2). ! ! Modified: ! ! 16 June 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X1, Y1, X2, Y2, define the first line. ! ! Input, real X3, Y3, X4, Y4, define the second line. ! ! Output, integer IVAL, reports on the intersection: ! 0, no intersection, the lines may be parallel or degenerate. ! 1, one intersection point, returned in X, Y. ! 2, infinitely many intersections, the lines are identical. ! ! Output, real X, Y, if IVAl = 1, then X, Y contains ! the intersection point. Otherwise, X = 0, Y = 0. ! implicit none ! real a1 real a2 real b1 real b2 real c1 real c2 integer ival logical point_1 logical point_2 real x real x1 real x2 real x3 real x4 real y real y1 real y2 real y3 real y4 ! ival = 0 x = 0.0E+00 y = 0.0E+00 ! ! Check whether either line is a point. ! if ( x1 == x2 .and. y1 == y2 ) then point_1 = .true. else point_1 = .false. end if if ( x3 == x4 .and. y3 == y4 ) then point_2 = .true. else point_2 = .false. end if ! ! Convert the lines to ABC format. ! if ( .not. point_1 ) then call line_exp2imp_2d ( x1, y1, x2, y2, a1, b1, c1 ) end if if ( .not. point_2 ) then call line_exp2imp_2d ( x3, y3, x4, y4, a2, b2, c2 ) end if ! ! Search for intersection of the lines. ! if ( point_1 .and. point_2 ) then if ( x1 == x3 .and. y1 == y3 ) then ival = 1 x = x1 y = y1 end if else if ( point_1 ) then if ( a2 * x1 + b2 * y1 == c2 ) then ival = 1 x = x1 y = y1 end if else if ( point_2 ) then if ( a1 * x3 + b1 * y3 == c1 ) then ival = 1 x = x3 y = y3 end if else call lines_imp_int_2d ( a1, b1, c1, a2, b2, c2, ival, x, y ) end if return end subroutine lines_imp_int_2d ( a1, b1, c1, a2, b2, c2, ival, x, y ) ! !******************************************************************************* ! !! LINES_IMP_INT_2D determines where two implicit lines intersect in 2D. ! ! ! Formula: ! ! The implicit form of a line in 2D is: ! ! A * X + B * Y + C = 0 ! ! Modified: ! ! 16 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A1, B1, C1, define the first line. ! At least one of A1 and B1 must be nonzero. ! ! Input, real A2, B2, C2, define the second line. ! At least one of A2 and B2 must be nonzero. ! ! Output, integer IVAL, reports on the intersection. ! ! -1, both A1 and B1 were zero. ! -2, both A2 and B2 were zero. ! 0, no intersection, the lines are parallel. ! 1, one intersection point, returned in X, Y. ! 2, infinitely many intersections, the lines are identical. ! ! Output, real X, Y, if IVAL = 1, then X, Y contains ! the intersection point. Otherwise, X = 0, Y = 0. ! implicit none ! real a(2,2) real a1 real a2 real b(2,2) real b1 real b2 real c1 real c2 real det integer ival real x real y ! x = 0.0E+00 y = 0.0E+00 ! ! Refuse to handle degenerate lines. ! if ( a1 == 0.0E+00 .and. b1 == 0.0E+00 ) then ival = -1 return else if ( a2 == 0.0E+00 .and. b2 == 0.0E+00 ) then ival = -2 return end if ! ! Set up a linear system, and compute its inverse. ! a(1,1) = a1 a(1,2) = b1 a(2,1) = a2 a(2,2) = b2 call rmat2_inverse ( a, b, det ) ! ! If the inverse exists, then the lines intersect. ! Multiply the inverse times -C to get the intersection point. ! if ( det /= 0.0E+00 ) then ival = 1 x = - b(1,1) * c1 - b(1,2) * c2 y = - b(2,1) * c1 - b(2,2) * c2 ! ! If the inverse does not exist, then the lines are parallel ! or coincident. Check for parallelism by seeing if the ! C entries are in the same ratio as the A or B entries. ! else ival = 0 if ( a1 == 0.0E+00 ) then if ( b2 * c1 == c2 * b1 ) then ival = 2 end if else if ( a2 * c1 == c2 * a1 ) then ival = 2 end if end if end if return end function lrline ( xu, yu, xv1, yv1, xv2, yv2, dv ) ! !******************************************************************************* ! !! LRLINE determines where a point lies in relation to a directed line. ! ! ! Discussion: ! ! LRLINE determines whether a point is to the left of, right of, ! or on a directed line parallel to a line through given points. ! ! Modified: ! ! 14 July 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Parameters: ! ! Input, real XU, YU, XV1, YV1, XV2, YV2, are vertex coordinates; the ! directed line is parallel to and at signed distance DV to the left of ! the directed line from (XV1,YV1) to (XV2,YV2); (XU,YU) is the vertex for ! which the position relative to the directed line is to be determined. ! ! Input, real DV, the signed distance, positive for left. ! ! Output, integer LRLINE, is +1, 0, or -1 depending on whether (XU,YU) is ! to the right of, on, or left of the directed line. LRLINE is 0 if ! the line degenerates to a point. ! implicit none ! real dv real dx real dxu real dy real dyu integer lrline real t real, parameter :: tol = 0.0000001E+00 real tolabs real xu real xv1 real xv2 real yu real yv1 real yv2 ! 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**2 + dy**2 ) if ( tolabs < t ) then lrline = 1 else if ( -tolabs <= t ) then lrline = 0 else if ( t < -tolabs ) then lrline = -1 end if return end function points_avoid_point_nd ( ndim, n, xy_set, xy_test, tol ) ! !******************************************************************************* ! !! POINTS_AVOID_POINT_ND checks if a point is "far" from a set of points in ND. ! ! ! Discussion: ! ! The routine discards points that are too close to other points. ! The method used to check this is quadratic in the number of points, ! and may take an inordinate amount of time if there are a large ! number of points. ! ! The test point is "far enough" from an accepted point if ! the Euclidean distance is at least "TOL". In graphics, you probably ! can't distinguish items that are 1/100 of the graph size apart, ! and surely not 1/1000, so if WIDTH is the size of your image, you ! might try setting TOL to WIDTH/100. ! ! Modified: ! ! 24 April 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NDIM, the spatial dimension. ! ! Input, integer N, the number of points in the set. ! ! Input, real XY_SET(NDIM,N), the set of points to be avoided. ! ! Input, real XY_TEST(NDIM), a point to be tested. ! ! Input, real TOL, the tolerance to be used. ! ! Output, logical POINTS_AVOID_POINT_ND, is TRUE if XY_TEST is ! "far enough" from all the set of points. ! implicit none ! integer n integer ndim ! integer j logical points_avoid_point_nd real tol real xy_set(ndim,n) real xy_test(ndim) ! points_avoid_point_nd = .true. do j = 1, n if ( sum ( ( xy_set(1:ndim,j) - xy_test(1:ndim) )**2 ) < tol**2 ) then points_avoid_point_nd = .false. return end if end do return end subroutine points_count ( file_in_name, dim_num, point_num ) ! !******************************************************************************* ! !! POINTS_COUNT counts the valid point coordinates in a file. ! ! ! Discussion: ! ! The routine reads every line, and expects to find DIM_NUM ! real numbers on the line. ! ! It does not count lines that begin with a comment symbol '#'. ! ! Modified: ! ! 05 March 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_IN_NAME, the name of the input file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Output, integer POINT_NUM, the number of point coordinate records. ! implicit none ! integer dim_num ! integer bad_num integer comment_num character ( len = * ) file_in_name integer file_in_unit integer i integer ierror integer ios character ( len = 100 ) line integer point_num integer record_num real x(dim_num) ! call get_unit ( file_in_unit ) open ( unit = file_in_unit, file = file_in_name, status = 'old', & iostat = ios ) if ( ios /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'POINTS_COUNT - Fatal error!' write ( *, '(a)' ) ' Could not open the input file: ' // & trim ( file_in_name ) stop end if comment_num = 0 point_num = 0 record_num = 0 bad_num = 0 do read ( file_in_unit, '(a)', iostat = ios ) line if ( ios /= 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 call s_to_rvec ( line, dim_num, x, ierror ) if ( ierror == 0 ) then point_num = point_num + 1 else bad_num = bad_num + 1 end if end do close ( unit = file_in_unit ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'POINTS_COUNT:' write ( *, '(a,i6)' ) ' Number of records: ', record_num write ( *, '(a,i6)' ) ' Number of point records: ', point_num write ( *, '(a,i6)' ) ' Number of comment records: ', comment_num write ( *, '(a,i6)' ) ' Number of bad records: ', bad_num return end subroutine points_read ( connect, file_in_name, dim_num, point_num, coord ) ! !******************************************************************************* ! !! POINTS_READ reads point coordinates from a file. ! ! ! Modified: ! ! 06 March 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer CONNECT(POINT_NUM), connection indicator. ! If CONNECT(I) is: ! 0, it is an isolated point; ! 1, it is connected to the previous point only; ! 2, it is connected to the next point only; ! 3, it is connected to the previous and next points. ! ! Input, character ( len = * ) FILE_IN_NAME, the name of the input file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer POINT_NUM, the number of points. The program ! will stop reading data once POINT_NUM values have been read. ! ! Output, real COORD(DIM_NUM,POINT_NUM), the point coordinates. ! implicit none ! integer dim_num integer point_num ! integer connect(point_num) real coord(dim_num,point_num) character ( len = * ) file_in_name integer file_in_unit integer i integer ierror integer ios character ( len = 100 ) line integer p real x(dim_num) ! call get_unit ( file_in_unit ) open ( unit = file_in_unit, file = file_in_name, status = 'old', & iostat = ios ) if ( ios /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'POINTS_READ - Fatal error!' write ( *, '(a)' ) ' Could not open the input file: ' // & trim ( file_in_name ) stop end if connect(1:point_num) = 0 i = 0 p = 0 do while ( i < point_num ) read ( file_in_unit, '(a)', iostat = ios ) line if ( ios /= 0 ) then ierror = i exit end if if ( line(1:1) == '#' .or. len_trim ( line ) == 0 ) then p = 0 cycle end if call s_to_rvec ( line, dim_num, x, ierror ) if ( ierror /= 0 ) then p = 0 cycle end if i = i + 1 coord(1:dim_num,i) = x(1:dim_num) connect(i) = p if ( p == 1 .and. i > 1 ) then connect(i-1) = connect(i-1) + 2 end if p = 1 end do close ( unit = file_in_unit ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'POINTS_READ:' write ( *, '(a,i6)' ) ' Read coordinate data from file.' return end subroutine points_thin ( connect, dim_num, point_num, coord, tol, point_num2 ) ! !******************************************************************************* ! !! POINTS_THIN "thins" points that are too close to each other. ! ! ! Modified: ! ! 06 March 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer CONNECT(POINT_NUM), connection indicator. ! If CONNECT(I) is: ! 0, it is an isolated point; ! 1, it is connected to the previous point only; ! 2, it is connected to the next point only; ! 3, it is connected to the previous and next points. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer POINT_NUM, the number of points. ! ! Input/output, real COORD(DIM_NUM,POINT_NUM), the point coordinates. ! ! Input, real TOL, a tolerance to use for thinning. ! ! Output, integer POINT_NUM2, the number of points kept. ! implicit none ! integer dim_num integer point_num ! integer check integer connect(point_num) real coord(dim_num,point_num) integer point_num2 integer thin_num logical points_avoid_point_nd real tol ! point_num2 = 0 thin_num = 0 do check = 1, point_num ! ! Always keep the first point. ! if ( check == 1 ) then point_num2 = 1 ! ! Keep the next point if it's far enough away. ! else if ( points_avoid_point_nd & ( dim_num, point_num2, coord, coord(1,check), tol ) ) then point_num2 = point_num2 + 1 coord(1:dim_num,point_num2) = coord(1:dim_num,check) connect(point_num2) = connect(check) else thin_num = thin_num + 1 if ( connect(check) == 0 ) then else if ( connect(check) == 1 ) then connect(point_num2) = 1 else if ( connect(check) == 2 ) then connect(check+1) = 2 else if ( connect(check) == 3 ) then end if end if end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'POINTS_THIN:' write ( *, '(a,g14.6)' ) ' Thinning tolerance was ', tol write ( *, '(a,i6)' ) ' Number of points given was ', point_num write ( *, '(a,i6)' ) ' Number of points retained ', point_num2 write ( *, '(a,i6)' ) ' Number of points thinned out ', thin_num return end function radians_to_degrees ( angle ) ! !******************************************************************************* ! !! RADIANS_TO_DEGREES converts an angle from radians to degrees. ! ! ! Modified: ! ! 10 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ANGLE, an angle in radians. ! ! Output, real RADIANS_TO_DEGREES, the equivalent angle in degrees. ! implicit none ! real angle real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real radians_to_degrees ! radians_to_degrees = ( angle / pi ) * 180.0E+00 return end subroutine rhpsrt ( k, n, lda, a, map ) ! !******************************************************************************* ! !! RHPSRT sorts points into lexicographic order using heap sort. ! ! ! Discussion: ! ! This routine uses heapsort to obtain the permutation of N K-dimensional ! points so that the points are in lexicographic increasing order. ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Modified: ! ! 19 February 2001 ! ! Parameters: ! ! Input, integer K, the dimension of the points (for instance, 2 ! for points in the plane). ! ! Input, integer N, the number of points. ! ! Input, integer LDA, the leading dimension of array A in the calling ! routine; LDA should be at least K. ! ! Input, real A(LDA,N); A(I,J) contains the I-th coordinate ! of point J. ! ! Input/output, integer MAP(N). ! On input, the points of A with indices MAP(1), MAP(2), ..., ! MAP(N) are to be sorted ! ! On output, MAP contains a permutation of its input values, with the ! property that, lexicographically, ! A(*,MAP(1)) <= A(*,MAP(2)) <= ... <= A(*,MAP(N)) ! implicit none ! integer lda integer n ! real a(lda,n) integer i integer k integer map(n) ! do i = n/2, 1, -1 call rsftdw ( i, n, k, lda, a, map ) end do do i = n, 2, -1 call i_swap ( map(1), map(i) ) call rsftdw ( 1, i-1, k, lda, a, map ) end do return end function rless ( k, p, q ) ! !******************************************************************************* ! !! RLESS determines whether P is lexicographically less than Q. ! ! ! Discussion: ! ! P and Q are K-dimensional points. ! ! Modified: ! ! 19 February 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Parameters: ! ! Input, integer K, the spatial dimension of the points. ! ! Input, real P(K), Q(K), the points to be compared. ! ! Output, logical RLESS, is TRUE if P < Q, FALSE otherwise. ! implicit none ! integer k ! real cmax integer i real p(k) real q(k) logical rless real tol ! tol = 100.0E+00 * epsilon ( tol ) do i = 1, k cmax = max ( abs ( p(i) ), abs ( q(i) ) ) if ( abs ( p(i) - q(i) ) > tol * cmax .and. cmax > tol ) then if ( p(i) < q(i) ) then rless = .true. else rless = .false. end if return end if end do rless = .false. return end subroutine rmat_print ( lda, m, n, a, title ) ! !******************************************************************************* ! !! RMAT_PRINT prints a real matrix. ! ! ! Modified: ! ! 04 September 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer M, the number of rows in A. ! ! Input, integer N, the number of columns in A. ! ! Input, real A(LDA,N), the matrix to be printed. ! ! Input, character ( len = * ) TITLE, a title to be printed first. ! TITLE may be blank. ! implicit none ! integer lda integer n ! real a(lda,n) integer i integer j integer jhi integer jlo integer m character ( len = * ) title ! if ( title /= ' ' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do jlo = 1, n, 5 jhi = min ( jlo + 4, n ) write ( *, '(a)' ) ' ' write ( *, '(6x,5(i7,7x))' ) ( j, j = jlo, jhi ) write ( *, '(a)' ) ' ' do i = 1, m write ( *, '(i6,5g14.6)' ) i, a(i,jlo:jhi) end do end do return end subroutine rmat2_inverse ( a, b, det ) ! !******************************************************************************* ! !! RMAT2_INVERSE inverts a 2 by 2 real matrix using Cramer's rule. ! ! ! Modified: ! ! 16 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(2,2), the matrix to be inverted. ! ! Output, real B(2,2), the inverse of the matrix A. ! ! Output, real DET, the determinant of the matrix A. ! ! If DET is zero, then A is singular, and does not have an ! inverse. In that case, B is simply set to zero, and a ! message is printed. ! ! If DET is nonzero, then its value is roughly an estimate ! of how nonsingular the matrix A is. ! implicit none ! real a(2,2) real b(2,2) real det ! ! Compute the determinant. ! det = a(1,1) * a(2,2) - a(1,2) * a(2,1) ! ! If the determinant is zero, bail out. ! if ( det == 0.0E+00 ) then b(1:2,1:2) = 0.0E+00 return end if ! ! Compute the entries of the inverse matrix using an explicit formula. ! b(1,1) = + a(2,2) / det b(1,2) = - a(1,2) / det b(2,1) = - a(2,1) / det b(2,2) = + a(1,1) / det return end subroutine rmat_solve ( a, n, nrhs, info ) ! !******************************************************************************* ! !! RMAT_SOLVE uses Gauss-Jordan elimination to solve an N by N linear system. ! ! ! Modified: ! ! 08 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(N,N+NRHS), contains in rows and columns 1 ! to N the coefficient matrix, and in columns N+1 through ! N+NRHS, the right hand sides. On output, the coefficient matrix ! area has been destroyed, while the right hand sides have ! been overwritten with the corresponding solutions. ! ! Input, integer NRHS, the number of right hand sides. NRHS ! must be at least 0. ! ! Output, integer INFO, singularity flag. ! 0, the matrix was not singular, the solutions were computed; ! J, factorization failed on step J, and the solutions could not ! be computed. ! implicit none ! integer n integer nrhs ! real a(n,n+nrhs) real apivot real factor integer i integer info integer ipivot integer j integer k real temp ! info = 0 do j = 1, n ! ! Choose a pivot row. ! ipivot = j apivot = a(j,j) do i = j+1, n if ( abs ( a(i,j) ) > abs ( apivot ) ) then apivot = a(i,j) ipivot = i end if end do if ( apivot == 0.0E+00 ) then info = j return end if ! ! Interchange. ! do i = 1, n + nrhs call r_swap ( a(ipivot,i), a(j,i) ) end do ! ! A(J,J) becomes 1. ! a(j,j) = 1.0E+00 a(j,j+1:n+nrhs) = a(j,j+1:n+nrhs) / apivot ! ! A(I,J) becomes 0. ! do i = 1, n if ( i /= j ) then factor = a(i,j) a(i,j) = 0.0E+00 a(i,j+1:n+nrhs) = a(i,j+1:n+nrhs) - factor * a(j,j+1:n+nrhs) end if end do end do return end subroutine rsftdw ( l, u, k, lda, a, map ) ! !******************************************************************************* ! !! RSFTDW shifts A(*,MAP(L)) down a heap of size U. ! ! ! Modified: ! ! 19 February 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Parameters: ! ! Input, integer L, U, the lower and upper indices of part of the heap. ! ! Input, integer K, the spatial dimension of the points. ! ! Input, integer LDA, the leading dimension of A in the calling routine. ! ! Input, real A(LDA,N); A(I,J) contains the I-th coordinate ! of point J. ! ! Input/output, integer MAP(N). ! On input, the points of A with indices MAP(1), MAP(2), ..., ! MAP(N) are to be sorted ! ! On output, MAP contains a permutation of its input values, with the ! property that, lexicographically, ! A(*,MAP(1)) <= A(*,MAP(2)) <= ... <= A(*,MAP(N)) ! implicit none ! integer lda ! real a(lda,*) integer i integer j integer k integer l integer map(*) logical rless integer t integer u ! i = l j = 2 * i t = map(i) do while ( j <= u ) if ( j < u ) then if ( rless ( k, a(1,map(j)), a(1,map(j+1)) ) ) then j = j + 1 end if end if if ( rless ( k, a(1,map(j)), a(1,t) ) ) then exit end if map(i) = map(j) i = j j = 2 * i end do map(i) = t return end subroutine rtris2 ( npt, maxst, vcl, ind, ntri, til, tnbr, stack, ierr ) ! !******************************************************************************* ! !! RTRIS2 constructs a Delaunay triangulation of 2D 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. ! ! Modified: ! ! 19 February 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Parameters: ! ! Input, integer NPT, the number of vertices. ! ! Input, integer MAXST, the maximum size available for the STACK array; ! should be about NPT to be safe, but MAX(10,2*LOG2(NPT)) is usually enough. ! ! Input, real VCL(2,NPT), the coordinates of the vertices. ! ! Input/output, integer IND(NPT), the indices in VCL of the vertices ! to be triangulated. On output, IND has been permuted by the sort. ! ! Output, integer NTRI, the number of triangles in the triangulation; ! NTRI is equal to 2*NPT - NB - 2, where NB is the number of boundary ! vertices. ! ! Output, integer TIL(3,NTRI), the nodes that make up each triangle. ! The elements are indices of VCL. The vertices of the triangles are ! in counter clockwise order ! ! Output, integer TNBR(3,NTRI), 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; TNBR(J,I) refers to ! the neighbor along edge from vertex J to J+1 (mod 3) ! ! Workspace, integer STACK(MAXST), used for a stack of triangles for which ! the circumcircle test must be made. ! ! Output, integer IERR, an error flag, nonzero if an error occurred. ! implicit none ! integer maxst integer npt ! real cmax integer e integer i integer ierr integer ind(npt) integer j integer k integer l integer ledg integer lr integer lrline integer ltri integer m integer m1 integer m2 integer n integer ntri integer redg integer rtri integer stack(maxst) integer t real temp integer til(3,npt*2) integer tnbr(3,npt*2) real tol integer top real vcl(2,npt) ! ierr = 0 tol = 100.0E+00 * epsilon ( tol ) ierr = 0 ! ! Sort the vertices by increasing (x,y) and obtain initial triangle(s). ! call rhpsrt ( 2, npt, 2, vcl, ind ) m1 = ind(1) do i = 2, npt m = m1 m1 = ind(i) k = 0 do j = 1, 2 cmax = max ( abs ( vcl(j,m) ), abs ( vcl(j,m1) ) ) if ( abs ( vcl(j,m) - vcl(j,m1) ) > tol * ( cmax + 1.0E+00 ) ) then k = j exit end if end do if ( k == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RTRIS2 - Fatal error!' write ( *, '(a,i6)' ) ' Fails for point number I = ', i write ( *, '(a,i6)' ) ' M = ', m write ( *, '(a,i6)' ) ' M1 = ', m1 write ( *, '(a,2g14.6)' ) ' VCL(1,M), VCL(2,M) = ', vcl(1,m), vcl(2,m) write ( *, '(a,2g14.6)' ) ' VCL(1,M1), VCL(2,M1) = ', vcl(1,m1), vcl(2,m1) ierr = 224 return end if end do m1 = ind(1) m2 = ind(2) j = 3 do if ( j > npt ) then ierr = 225 return end if m = ind(j) lr = lrline ( vcl(1,m), vcl(2,m), vcl(1,m1), vcl(2,m1), vcl(1,m2), & vcl(2,m2), 0.0E+00 ) if ( lr /= 0 ) then exit end if j = j + 1 end do ntri = j - 2 if ( lr == -1 ) then til(1,1) = m1 til(2,1) = m2 til(3,1) = m tnbr(3,1) = -3 do i = 2, ntri m1 = m2 m2 = ind(i+1) til(1,i) = m1 til(2,i) = m2 til(3,i) = m tnbr(1,i-1) = -3 * i tnbr(2,i-1) = i tnbr(3,i) = i - 1 end do tnbr(1,ntri) = -3 * ntri - 1 tnbr(2,ntri) = -5 ledg = 2 ltri = ntri else til(1,1) = m2 til(2,1) = m1 til(3,1) = m tnbr(1,1) = -4 do i = 2, ntri m1 = m2 m2 = ind(i+1) til(1,i) = m2 til(2,i) = m1 til(3,i) = m tnbr(3,i-1) = i tnbr(1,i) = -3 * i - 3 tnbr(2,i) = i - 1 end do tnbr(3,ntri) = -3 * ntri tnbr(2,1) = -3 * ntri - 2 ledg = 2 ltri = 1 end if ! ! Insert vertices one at a time from outside 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, npt m = ind(i) m1 = til(ledg,ltri) if ( ledg <= 2 ) then m2 = til(ledg+1,ltri) else m2 = til(1,ltri) end if lr = lrline ( vcl(1,m), vcl(2,m), vcl(1,m1), vcl(2,m1), & vcl(1,m2), vcl(2,m2), 0.0E+00 ) if ( lr > 0 ) then rtri = ltri redg = ledg ltri = 0 else l = -tnbr(ledg,ltri) rtri = l / 3 redg = mod(l,3) + 1 end if call vbedg ( vcl(1,m), vcl(2,m), vcl, til, tnbr, ltri, ledg, rtri, redg ) n = ntri + 1 l = -tnbr(ledg,ltri) do t = l / 3 e = mod ( l, 3 ) + 1 l = -tnbr(e,t) m2 = til(e,t) if ( e <= 2 ) then m1 = til(e+1,t) else m1 = til(1,t) end if ntri = ntri + 1 tnbr(e,t) = ntri til(1,ntri) = m1 til(2,ntri) = m2 til(3,ntri) = m tnbr(1,ntri) = t tnbr(2,ntri) = ntri - 1 tnbr(3,ntri) = ntri + 1 top = top + 1 if ( top > maxst ) then ierr = 8 return end if stack(top) = ntri if ( t == rtri .and. e == redg ) then exit end if end do tnbr(ledg,ltri) = -3 * n - 1 tnbr(2,n) = -3 * ntri - 2 tnbr(3,ntri) = -l ltri = n ledg = 2 call swapec ( m, top, maxst, ltri, ledg, vcl, til, tnbr, stack, ierr ) if ( ierr /= 0 ) then return end if end do return end subroutine s_cap ( s ) ! !******************************************************************************* ! !! S_CAP replaces any lowercase letters by uppercase ones in a string. ! ! ! Modified: ! ! 16 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) S, the string to be transformed. ! implicit none ! character c integer i integer nchar character ( len = * ) s ! nchar = len ( s ) do i = 1, nchar c = s(i:i) call ch_cap ( c ) s(i:i) = c end do return end function s_index_last ( s, sub ) ! !******************************************************************************* ! !! 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. ! ! 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 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 i integer j integer llen1 integer llen2 character ( len = * ) s integer 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 ( llen2 > llen1 ) 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_r ( s, r, ierror, lchar ) ! !******************************************************************************* ! !! S_TO_R reads a real number from a string. ! ! ! Discussion: ! ! This 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 real number. ! ! Legal input is: ! ! 1 blanks, ! 2 '+' or '-' sign, ! 2.5 spaces ! 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. ! ! Examples: ! ! S R ! ! '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) ! ! Modified: ! ! 12 February 2001 ! ! 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 R, the real value that was read from the string. ! ! Output, integer 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 LCHAR, the number of characters read from ! the string to form the number, including any terminating ! characters such as a trailing comma or blanks. ! implicit none ! logical ch_eqi character c integer ierror integer ihave integer isgn integer iterm integer jbot integer jsgn integer jtop integer lchar integer nchar integer ndig real r real rbot real rexp real rtop character ( len = * ) s character, parameter :: TAB = char ( 9 ) ! nchar = len_trim ( s ) ierror = 0 r = 0.0E+00 lchar = - 1 isgn = 1 rtop = 0.0E+00 rbot = 1.0E+00 jsgn = 1 jtop = 0 jbot = 1 ihave = 1 iterm = 0 do lchar = lchar + 1 c = s(lchar+1:lchar+1) ! ! Blank or TAB character. ! if ( c == ' ' .or. c == TAB ) then if ( ihave == 2 ) then else if ( ihave == 6 .or. ihave == 7 ) then iterm = 1 else if ( ihave > 1 ) then ihave = 11 end if ! ! Comma. ! else if ( c == ',' .or. c == ';' ) then if ( ihave /= 1 ) then iterm = 1 ihave = 12 lchar = lchar + 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 ( ihave >= 6 .and. ihave <= 8 ) then ihave = 9 else iterm = 1 end if ! ! 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. lge ( c, '0' ) .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.0E+00 * rtop + real ( ndig ) else if ( ihave == 5 ) then rtop = 10.0E+00 * rtop + real ( ndig ) rbot = 10.0E+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 .or. lchar+1 >= nchar ) 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 LCHAR is equal to NCHAR. ! if ( iterm /= 1 .and. lchar+1 == nchar ) then lchar = 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 return end if ! ! Number seems OK. Form it. ! if ( jtop == 0 ) then rexp = 1.0E+00 else if ( jbot == 1 ) then rexp = 10.0E+00**( jsgn * jtop ) else rexp = jsgn * jtop rexp = rexp / jbot rexp = 10.0E+00**rexp end if end if r = isgn * rexp * rtop / rbot return end subroutine s_to_rvec ( s, n, rvec, ierror ) ! !******************************************************************************* ! !! S_TO_RVEC reads a real vector from a string. ! ! ! Modified: ! ! 19 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be read. ! ! Input, integer N, the number of values expected. ! ! Output, real RVEC(N), the values read from the string. ! ! Output, integer IERROR, error flag. ! 0, no errors occurred. ! -K, could not read data for entries -K through N. ! implicit none ! integer n ! integer i integer ierror integer ilo integer lchar real rvec(n) character ( len = * ) s ! i = 0 ilo = 1 do while ( i < n ) i = i + 1 call s_to_r ( s(ilo:), rvec(i), ierror, lchar ) if ( ierror /= 0 ) then ierror = -i exit end if ilo = ilo + lchar end do return end subroutine swapec ( i, top, maxst, btri, bedg, vcl, til, tnbr, stack, ierr ) ! !******************************************************************************* ! !! SWAPEC 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. ! ! Modified: ! ! 14 July 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Parameters: ! ! Input, integer I, the index in VCL of the new vertex. ! ! Input/output, integer TOP, the index of the top of the stack. ! On output, TOP is zero. ! ! Input, integer MAXST, the maximum size available for the STACK array. ! ! Input/output, integer 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, real VCL(2,*), the coordinates of the vertices. ! ! Input/output, integer TIL(3,*), the triangle incidence list. May be ! updated on output because of swaps. ! ! Input/output, integer TNBR(3,*), 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(1: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 IERR is set to 8 for abnormal return. ! implicit none ! integer maxst ! integer a integer b integer bedg integer btri integer c integer diaedg integer e integer ee integer em1 integer ep1 integer f integer fm1 integer fp1 integer i integer i_wrap integer ierr integer l integer r integer s integer stack(maxst) integer swap integer t integer til(3,*) integer tnbr(3,*) integer top integer tt integer u real vcl(2,*) real x real y ! ! Determine whether triangles in stack are Delaunay, and swap ! diagonal edge of convex quadrilateral if not. ! ierr = 0 x = vcl(1,i) y = vcl(2,i) do if ( top <= 0 ) then exit end if t = stack(top) top = top - 1 if ( til(1,t) == i ) then e = 2 b = til(3,t) else if ( til(2,t) == i ) then e = 3 b = til(1,t) else e = 1 b = til(2,t) end if a = til(e,t) u = tnbr(e,t) if ( tnbr(1,u) == t ) then f = 1 c = til(3,u) else if ( tnbr(2,u) == t ) then f = 2 c = til(1,u) else f = 3 c = til(2,u) end if swap = diaedg ( x, y, vcl(1,a), vcl(2,a), vcl(1,c), vcl(2,c), & vcl(1,b), vcl(2,b) ) if ( swap == 1 ) then em1 = i_wrap ( e - 1, 1, 3 ) ep1 = i_wrap ( e + 1, 1, 3 ) fm1 = i_wrap ( f - 1, 1, 3 ) fp1 = i_wrap ( f + 1, 1, 3 ) til(ep1,t) = c til(fp1,u) = i r = tnbr(ep1,t) s = tnbr(fp1,u) tnbr(ep1,t) = u tnbr(fp1,u) = t tnbr(e,t) = s tnbr(f,u) = r if ( tnbr(fm1,u) > 0 ) then top = top + 1 stack(top) = u end if if ( s > 0 ) then if ( tnbr(1,s) == u ) then tnbr(1,s) = t else if ( tnbr(2,s) == u ) then tnbr(2,s) = t else tnbr(3,s) = t end if top = top + 1 if ( top > maxst ) 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 ( tnbr(ee,tt) > 0 ) tt = tnbr(ee,tt) if ( til(1,tt) == a ) then ee = 3 else if ( til(2,tt) == a ) then ee = 1 else ee = 2 end if end do tnbr(ee,tt) = l end if if ( r > 0 ) then if ( tnbr(1,r) == t ) then tnbr(1,r) = u else if ( tnbr(2,r) == t ) then tnbr(2,r) = u else tnbr(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 ( tnbr(ee,tt) > 0 ) tt = tnbr(ee,tt) if ( til(1,tt) == b ) then ee = 3 else if ( til(2,tt) == b ) then ee = 1 else ee = 2 end if end do tnbr(ee,tt) = l end if end if end do return end subroutine timestamp ( ) ! !******************************************************************************* ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! character ( len = 8 ) ampm integer d character ( len = 8 ) date integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone ! call date_and_time ( date, time, zone, 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 ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine triangle_circumcenter_2d ( x1, y1, x2, y2, x3, y3, x, y ) ! !******************************************************************************* ! !! TRIANGLE_CIRCUMCENTER_2D computes the circumcenter of a triangle in 2D. ! ! ! Discussion: ! ! The exscribed circle of a triangle is the circle that passes through ! the three vertices of the triangle. The exscribed circle contains ! the triangle, but it is not necessarily the smallest triangle to do so. ! ! The circumcenter is the center of the exscribed circle. ! If all angles of the triangle are no greater than 90 degrees, then ! the center of the exscribed circle will lie inside the triangle. ! Otherwise, the center will lie outside the circle. ! ! Surprisingly, the diameter of the circle can be found by solving ! a 2 by 2 linear system. This is because the vectors P2 - P1 ! and P3 - P1 are secants of the circle, and each forms a right ! triangle with the diameter. Hence, the dot product of ! P2 - P1 with the diameter is equal to the square of the length ! of P2 - P1, and similarly for P3 - P1. This determines the ! diameter vector originating at P1. ! ! Reference: ! ! Adrian Bowyer and John Woodwark, ! A Programmer's Geometry, ! Butterworths, 1983. ! ! Modified: ! ! 30 April 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X1, Y1, X2, Y2, X3, Y3, the triangle corners. ! ! Output, real X, Y, the coordinates of the circumcenter of the triangle. ! If the linear system is singular, then X = Y = 0. ! implicit none ! integer, parameter :: n = 2 integer, parameter :: nrhs = 1 ! real a(n,n+nrhs) real enormsq0_2d integer info real x real x1 real x2 real x3 real y real y1 real y2 real y3 ! ! Set up the linear system. ! a(1,1) = x2 - x1 a(1,2) = y2 - y1 a(1,3) = enormsq0_2d ( x1, y1, x2, y2 ) a(2,1) = x3 - x1 a(2,2) = y3 - y1 a(2,3) = enormsq0_2d ( x1, y1, x3, y3 ) ! ! Solve the linear system. ! call rmat_solve ( a, n, nrhs, info ) ! ! Compute X, Y. ! if ( info /= 0 ) then x = 0.0E+00 y = 0.0E+00 else x = x1 + 0.5E+00 * a(1,n+1) y = y1 + 0.5E+00 * a(2,n+1) end if return end subroutine triangle_contains_point_2d ( x1, y1, x2, y2, x3, y3, x, y, inside ) ! !******************************************************************************* ! !! TRIANGLE_CONTAINS_POINT_2D finds if a point is inside a triangle in 2D. ! ! ! Modified: ! ! 12 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X1, Y1, X2, Y2, X3, Y3, the coordinates of ! the corners of the triangle. ! ! Input, real X, Y, the point to be checked. ! ! Output, logical INSIDE, is .TRUE. if (X,Y) is inside ! the triangle or on its boundary, and .FALSE. otherwise. ! implicit none ! integer, parameter :: n = 2 integer, parameter :: nrhs = 1 ! real a(n,n+nrhs) real c1 real c2 integer info logical inside real x real x1 real x2 real x3 real y real y1 real y2 real y3 ! ! Set up the linear system ! ! ( X2-X1 X3-X1 ) C1 = X-X1 ! ( Y2-Y1 Y3-Y1 ) C2 Y-Y1 ! ! which is satisfied by the barycentric coordinates of (X,Y). ! a(1,1) = x2 - x1 a(1,2) = x3 - x1 a(1,3) = x - x1 a(2,1) = y2 - y1 a(2,2) = y3 - y1 a(2,3) = y - y1 ! ! Solve the linear system. ! call rmat_solve ( a, n, nrhs, info ) if ( info /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TRIANGLE_CONTAINS_POINT_2D - Fatal error!' write ( *, '(a)' ) ' The linear system is singular.' write ( *, '(a)' ) ' The input data does not form a proper triangle.' stop end if c1 = a(1,3) c2 = a(2,3) ! ! If the point is in the triangle, its barycentric coordinates ! must both be nonnegative, and sum to no more than 1. ! if ( c1 < 0.0E+00 .or. c2 < 0.0E+00 ) then inside = .false. else if ( c1 + c2 > 1.0E+00 ) then inside = .false. else inside = .true. end if return end subroutine trivor_plot ( file_name, dim_num, point_num, coord, coord_min, & coord_max, x, y ) ! !******************************************************************************* ! !! TRIVOR_PLOT plots the triangulated Voronoi diagram of a set of points in 2D. ! ! ! Discussion: ! ! The routine first determines the Delaunay triangulation. ! The Voronoi diagram is then determined from this information. ! Finally, each Voronoi cell is triangulated, using the cell centroid ! and the vertices. ! ! Modified: ! ! 20 June 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the output file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer POINT_NUM, the number of points. ! ! Input, real COORD(DIM_NUM,POINT_NUM), the point coordinates. ! ! Input, real COORD_MIN(DIM_NUM), COORD_MAX(DIM_NUM), the data ranges. ! These values may be set to the minimum and maximum values of the data, ! or can be used to clip the data, or to put some whitespace around it. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! This lets COOR_MIN and COORD_MAX have a dimension larger than 2. ! We simply specify which pair of indices will play the X and Y roles. ! implicit none ! integer dim_num integer point_num ! real b real center(dim_num,2*point_num) real coord(dim_num,point_num) real coord_max(dim_num) real coord_min(dim_num) real coord2(2,point_num) logical, parameter :: debug = .false. character ( len = * ) file_name integer file_unit real g integer i integer i_wrap integer i1 integer i2 integer i3 integer ierror integer ind(point_num) logical inside integer ix integer j integer jp1 integer jp2 integer k real n1 real n2 integer nodtri(3,2*point_num) integer tri_num real r integer stack(point_num) integer tnbr(3,2*point_num) integer x real xvec(4) real x1 real x2 real x3 real x4 real x5 real x6 integer y real yvec(4) real y1 real y2 real y3 real y4 real y5 real y6 ! ! Compute the Delaunay triangulation. ! coord2(1,1:point_num) = coord(x,1:point_num) coord2(2,1:point_num) = coord(y,1:point_num) call ivec_identity ( point_num, ind ) call rtris2 ( point_num, point_num, coord2, ind, tri_num, nodtri, tnbr, & stack, ierror ) if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Data nodes:' write ( *, '(a)' ) ' ' do j = 1, point_num write ( *, '(i4,3g14.6)' ) j, coord2(1:2,j) end do end if if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Nodes in each Delaunay triangle:' write ( *, '(a)' ) ' ' do j = 1, tri_num write ( *, '(i4,3i4)' ) j, nodtri(1:3,j) end do end if if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Neighboring triangles of each Delaunay triangle:' write ( *, '(a)' ) ' ' do j = 1, tri_num write ( *, '(i4,3i4)' ) j, tnbr(1:3,j) end do end if ! ! Compute the intersection point of the perpendicular bisectors ! of each Delaunay triangle. ! do j = 1, tri_num i1 = nodtri(1,j) i2 = nodtri(2,j) i3 = nodtri(3,j) call triangle_circumcenter_2d ( coord2(1,i1), coord2(2,i1), coord2(1,i2), & coord2(2,i2), coord2(1,i3), coord2(2,i3), center(1,j), center(2,j) ) end do if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Circumcenters of Delaunay Triangles:' write ( *, '(a)' ) ' ' do j = 1, tri_num write ( *, '(i4,2g14.6)' ) j, center(1,j), center(2,j) end do end if call get_unit ( file_unit ) call ps_file_open ( file_name, file_unit, ierror ) call eps_file_head ( file_name ) call ps_page_head ( coord_min(x), coord_min(y), coord_max(x), coord_max(y) ) ! ! We need to define a PostScript clipping box, so that we can let PostScript ! take care of clipping lines that lie outside the box. ! xvec(1:4) = (/ coord_min(x), coord_max(x), coord_max(x), coord_min(x) /) yvec(1:4) = (/ coord_min(y), coord_min(y), coord_max(y), coord_max(y) /) call ps_clip ( 4, xvec, yvec ) x1 = coord_min(x) y1 = coord_min(y) call ps_moveto ( x1, y1 ) if ( point_num <= 350 ) then do i = 1, point_num call ps_mark_disk ( coord2(1,i), coord2(2,i) ) end do else do i = 1, point_num call ps_mark_point ( coord2(1,i), coord2(2,i) ) end do end if if ( point_num <= 350 ) then do i = 1, tri_num call ps_mark_circle ( center(1,i), center(2,i) ) end do end if ! ! For each Delaunay triangle, I ! For each side J, ! do i = 1, tri_num do j = 1, 3 k = tnbr(j,i) ! ! If there is a neighboring triangle K on that side, ! connect the circumcenters. ! if ( k > 0 ) then if ( k > i ) then call ps_line ( center(1,i), center(2,i), center(1,k), center(2,k) ) end if ! ! If there is no neighboring triangle on that side, ! extend a line from the circumcenter of I in the direction of the ! outward normal to that side. ! else ix = nodtri(j,i) x1 = coord(x,ix) y1 = coord(y,ix) jp1 = i_wrap ( j+1, 1, 3 ) ix = nodtri(jp1,i) x2 = coord(x,ix) y2 = coord(y,ix) jp2 = i_wrap ( j+2, 1, 3 ) ix = nodtri(jp2,i) x3 = coord(x,ix) y3 = coord(y,ix) x4 = center(1,i) y4 = center(2,i) call line_exp_normal_2d ( x1, y1, x2, y2, n1, n2 ) x5 = x4 + n1 y5 = y4 + n2 call box_ray_int_2d ( coord_min(x), coord_min(y), coord_max(x), & coord_max(y), x4, y4, x5, y5, x6, y6 ) call ps_line ( x4, y4, x6, y6 ) end if end do end do ! ! For each triangle, join its circumcenter to each of the three nodes, ! using a gray line. This will display a triangulation of the Voronoi ! polygons. ! r = 0.8E+00 g = 0.8E+00 b = 0.8E+00 call ps_color_line_set ( r, g, b ) do i = 1, tri_num x1 = center(1,i) y1 = center(2,i) do j = 1, 3 k = nodtri(j,i) x2 = coord(x,k) y2 = coord(y,k) call ps_line ( x1, y1, x2, y2 ) end do end do call ps_page_tail call eps_file_tail call ps_file_close ( file_unit ) return end subroutine vbedg ( x, y, vcl, til, tnbr, ltri, ledg, rtri, redg ) ! !******************************************************************************* ! !! 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. ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Modified: ! ! 18 June 2001 ! ! Parameters: ! ! Input, real X, Y, the coordinates of a point outside the convex hull ! of the current triangulation. ! ! Input, real VCL(1:2,1:*), the coordinates of the vertices. ! ! Input, integer TIL(1:3,1:*), the triangle incidence list. ! ! Input, integer TNBR(1:3,1:*), 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 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 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 REDG, the edge of triangle RTRI that is visible ! from (X,Y). 1 <= REDG <= 3. ! implicit none ! integer a integer b integer e integer i_wrap integer l logical ldone integer ledg integer lr integer lrline integer ltri integer redg integer rtri integer t real temp integer til(3,*) integer tnbr(3,*) real vcl(2,*) real x real 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 = -tnbr(redg,rtri) t = l / 3 e = mod ( l, 3 ) + 1 a = til(e,t) if ( e <= 2 ) then b = til(e+1,t) else b = til(1,t) end if temp = 0.0E+00 lr = lrline ( x, y, vcl(1,a), vcl(2,a), vcl(1,b), vcl(2,b), temp ) 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 = til(e,t) e = i_wrap ( e-1, 1, 3 ) do while ( tnbr(e,t) > 0 ) t = tnbr(e,t) if ( til(1,t) == b ) then e = 3 else if ( til(2,t) == b ) then e = 1 else e = 2 end if end do a = til(e,t) temp = 0.0E+00 lr = lrline ( x, y, vcl(1,a), vcl(2,a), vcl(1,b), vcl(2,b), temp ) if ( lr <= 0 ) then exit end if end do ltri = t ledg = e return end subroutine voronoi_plot ( file_name, dim_num, point_num, coord, coord_min, & coord_max, x, y ) ! !******************************************************************************* ! !! VORONOI_PLOT plots the Voronoi diagram of a set of points in 2D. ! ! ! Discussion: ! ! The routine first determines the Delaunay triangulation. ! The Voronoi diagram is then determined from this information. ! ! Modified: ! ! 08 March 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the output file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer POINT_NUM, the number of points. ! ! Input, real COORD(DIM_NUM,POINT_NUM), the point coordinates. ! ! Input, real COORD_MIN(DIM_NUM), COORD_MAX(DIM_NUM), the data ranges. ! These values may be set to the minimum and maximum values of the data, ! or can be used to clip the data, or to put some whitespace around it. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! This lets COOR_MIN and COORD_MAX have a dimension larger than 2. ! We simply specify which pair of indices will play the X and Y roles. ! implicit none ! integer dim_num integer point_num ! real center(dim_num,2*point_num) real coord(dim_num,point_num) real coord_max(dim_num) real coord_min(dim_num) real coord2(2,point_num) logical, parameter :: debug = .false. character ( len = * ) file_name integer file_unit integer i integer i_wrap integer i1 integer i2 integer i3 integer ierror integer ind(point_num) logical inside integer ix integer j integer jp1 integer jp2 integer k real n1 real n2 integer nodtri(3,2*point_num) integer tri_num integer stack(point_num) integer tnbr(3,2*point_num) integer x real xvec(4) real x1 real x2 real x3 real x4 real x5 real x6 integer y real yvec(4) real y1 real y2 real y3 real y4 real y5 real y6 ! ! Compute the Delaunay triangulation. ! coord2(1,1:point_num) = coord(x,1:point_num) coord2(2,1:point_num) = coord(y,1:point_num) call ivec_identity ( point_num, ind ) call rtris2 ( point_num, point_num, coord2, ind, tri_num, nodtri, tnbr, & stack, ierror ) if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Data nodes:' write ( *, '(a)' ) ' ' do j = 1, point_num write ( *, '(i4,3g14.6)' ) j, coord2(1:2,j) end do end if if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Nodes in each Delaunay triangle:' write ( *, '(a)' ) ' ' do j = 1, tri_num write ( *, '(i4,3i4)' ) j, nodtri(1:3,j) end do end if if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Neighboring triangles of each Delaunay triangle:' write ( *, '(a)' ) ' ' do j = 1, tri_num write ( *, '(i4,3i4)' ) j, tnbr(1:3,j) end do end if ! ! Compute the intersection point of the perpendicular bisectors ! of each Delaunay triangle. ! do j = 1, tri_num i1 = nodtri(1,j) i2 = nodtri(2,j) i3 = nodtri(3,j) call triangle_circumcenter_2d ( coord2(1,i1), coord2(2,i1), coord2(1,i2), & coord2(2,i2), coord2(1,i3), coord2(2,i3), center(1,j), center(2,j) ) end do if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Circumcenters of Delaunay Triangles:' write ( *, '(a)' ) ' ' do j = 1, tri_num write ( *, '(i4,2g14.6)' ) j, center(1,j), center(2,j) end do end if ! ! Plot the Voronoi tessellation. ! call get_unit ( file_unit ) call ps_file_open ( file_name, file_unit, ierror ) call eps_file_head ( file_name ) call ps_page_head ( coord_min(x), coord_min(y), coord_max(x), coord_max(y) ) ! ! We need to define a PostScript clipping box, so that we can let PostScript ! take care of clipping lines that lie outside the box. ! xvec(1:4) = (/ coord_min(x), coord_max(x), coord_max(x), coord_min(x) /) yvec(1:4) = (/ coord_min(y), coord_min(y), coord_max(y), coord_max(y) /) call ps_clip ( 4, xvec, yvec ) x1 = coord_min(x) y1 = coord_min(y) call ps_moveto ( x1, y1 ) ! call ps_label ( 'Voronoi tessellation' ) if ( point_num <= 350 ) then do i = 1, point_num call ps_mark_disk ( coord2(1,i), coord2(2,i) ) end do else do i = 1, point_num call ps_mark_point ( coord2(1,i), coord2(2,i) ) end do end if if ( point_num <= 350 ) then do i = 1, tri_num call ps_mark_circle ( center(1,i), center(2,i) ) end do end if ! ! For each Delaunay triangle, I ! For each side J, ! do i = 1, tri_num do j = 1, 3 k = tnbr(j,i) ! ! If there is a neighboring triangle K on that side, ! connect the circumcenters. ! if ( k > 0 ) then if ( k > i ) then call ps_line ( center(1,i), center(2,i), center(1,k), center(2,k) ) end if ! ! If there is no neighboring triangle on that side, ! extend a line from the circumcenter of I in the direction of the ! outward normal to that side. ! else ix = nodtri(j,i) x1 = coord(x,ix) y1 = coord(y,ix) jp1 = i_wrap ( j+1, 1, 3 ) ix = nodtri(jp1,i) x2 = coord(x,ix) y2 = coord(y,ix) jp2 = i_wrap ( j+2, 1, 3 ) ix = nodtri(jp2,i) x3 = coord(x,ix) y3 = coord(y,ix) x4 = center(1,i) y4 = center(2,i) call line_exp_normal_2d ( x1, y1, x2, y2, n1, n2 ) x5 = x4 + n1 y5 = y4 + n2 call box_ray_int_2d ( coord_min(x), coord_min(y), coord_max(x), & coord_max(y), x4, y4, x5, y5, x6, y6 ) call ps_line ( x4, y4, x6, y6 ) end if end do end do call ps_page_tail call eps_file_tail call ps_file_close ( file_unit ) return end subroutine word_count ( s, nword ) ! !******************************************************************************* ! !! WORD_COUNT counts the number of "words" in a string. ! ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be examined. ! ! Output, integer NWORD, the number of "words" in the string. ! Words are presumed to be separated by one or more blanks. ! implicit none ! logical blank integer i integer lens integer 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