module AA_point_data_module !*****************************************************************************80 ! !! AA_POINT_DATA_MODULE keeps track of some point data. ! ! Discussion: ! ! The name of this module has 'AA' preprended to it so that ! when this file is split up, and the pieces compiled, this ! piece will be compiled before the others. That's necessary ! so that routines that use this module will not be compiled ! first! Ah, the joys of learning modules! ! ! The arrays stored here are allocatable. The documentation ! lists their intended dimensions, which depend on NODE_NUM ! and DIM_NUM. ! ! Modified: ! ! 15 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Module data, integer CONNECT(NODE_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. ! ! Module data, real ( kind = 8 ) COORD(DIM_NUM,NODE_NUM), the point ! coordinates. ! ! Module data, real ( kind = 8 ) COORD_MIN(DIM_NUM), COORD_MAX(DIM_NUM), ! the data ranges. ! ! Module data, real ( kind = 8 ) RADIUS(NODE_NUM), the radius of a circle ! around each point. ! ! Module data, integer TAG(NODE_NUM), a tag for each point, ! presumably the cluster index. ! implicit none integer, allocatable, dimension ( : ) :: connect real ( kind = 8 ), allocatable, dimension ( :, : ) :: coord real ( kind = 8 ), allocatable, dimension ( : ) :: coord_max real ( kind = 8 ), allocatable, dimension ( : ) :: coord_min real ( kind = 8 ), allocatable, dimension ( : ) :: radius integer, allocatable, dimension ( : ) :: tag end program main !*****************************************************************************80 ! !! MAIN is the main program for PLOT_POINTS. ! ! Discussion: ! ! PLOT_POINTS plots the points in a file. ! ! Most recent change was that SHADOW is no longer a "toggle" ! command. Rather, SHADOW turns shadowing on, and NOSHADOW turns ! it off. ! ! Also, I started to set up a LINE command, but never got done with it. ! Any shoemaker's elves with nothing else to do, please come fix this! ! ! Usage: ! ! The program can be invoked by: ! ! plot_points plot_file_name ! ! or: ! ! plot_points ! ! in which case the user will be asked to supply the input file name. ! ! Modified: ! ! 19 November 2003 ! ! Author: ! ! John Burkardt ! use AA_point_data_module implicit none integer, parameter :: line_max = 10 integer arg_num logical :: box_requested = .false. real ( kind = 8 ), dimension ( 2, 2 ) :: box_xy = reshape ( & (/ 0.0D+00, 0.0D+00, 1.0D+00, 1.0D+00 /), (/ 2, 2 /) ) logical :: circle_requested = .false. real ( kind = 8 ) :: circle_r = 0.0D+00 real ( kind = 8 ) :: circle_x = 0.0D+00 real ( kind = 8 ) :: circle_y = 0.0D+00 real ( kind = 8 ), allocatable, dimension ( :, : ) :: coord2 character ( len = 100 ) command logical, parameter :: delaunay = .true. integer dim_num integer dim_num2 character ( len = 256 ) :: file_in_name = ' ' integer file_in_unit character ( len = 256 ) :: file_out_name = ' ' character ( len = 256 ) :: file_out_name_user = ' ' character ( len = 256 ) :: file_out_name_default = 'plot000.eps' integer first integer i integer iarg integer iargc integer :: ierror = 0 integer ilen integer ios integer ipxfargc integer last integer length integer lens integer :: line_num = 0 real ( kind = 8 ), dimension ( 2, 2, line_max ) :: line_xy integer :: marker_size = 4 integer next real ( kind = 8 ) plot_max(2) real ( kind = 8 ) plot_min(2) character ( len = 80 ) point_file integer node_num integer node_num2 logical s_eqi logical :: shadow = .false. character ( len = 80 ) :: title = ' ' real ( kind = 8 ) tol logical, parameter :: voronoi = .true. logical, parameter :: walls = .true. integer :: x_index = 1 integer :: y_index = 2 call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_POINTS' write ( *, '(a)' ) ' FORTRAN90 version' 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 ( 1 <= arg_num ) then 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 call data_read ( dim_num, file_in_name, node_num ) end if ! ! Command loop begins here. ! do if ( command(1:1) /= '#' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter command (H for help):' write ( *, '(a)' ) ' ' end if read ( *, '(a)' ) command ! ! # is a comment. ! Echo it in case we're saving the output. ! if ( command(1:1) == '#' ) then write ( *, '(a)' ) trim ( command ) ! ! BALLOON plot ! Plot each point with its maximum circle. ! else if ( s_eqi ( command(1:3), 'BAL' ) ) then if ( len_trim ( file_out_name_user ) == 0 ) then call file_name_inc ( file_out_name_default ) file_out_name = file_out_name_default else file_out_name = file_out_name_user file_out_name_user = ' ' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Write balloon plot data to "' // & trim ( file_out_name ) // '".' dim_num2 = 2 allocate ( coord2(1:2,1:node_num) ) coord2(1,1:node_num) = coord(x_index,1:node_num) coord2(2,1:node_num) = coord(y_index,1:node_num) call radius_maximus ( dim_num2, node_num, coord2, walls, radius ) deallocate ( coord2 ) call plot_size ( dim_num, box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, x_index, y_index, plot_min, plot_max ) call balloon_plot ( box_requested, box_xy, circle_requested, circle_x, & circle_y, circle_r, file_out_name, dim_num, & node_num, coord, plot_min, plot_max, x_index, y_index, shadow, & title, radius ) ! ! BOX x1, y1, x2, y2 ! specify a box to be drawn. ! else if ( s_eqi ( command(1:3), 'BOX' ) ) then box_requested = .true. command(1:3) = ' ' command = adjustl ( command ) if ( command(4:4) == '=' ) then command(4:4) = ' ' end if if ( len_trim ( command ) /= 0 ) then call s_to_r8vec ( command, 4, box_xy(1:2,1:2), ierror ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter lower left (x,y), then upper right (x,y):' read ( *, *, iostat = ios ) box_xy(1:2,1:2) if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_POINTS - Fatal error!' write ( *, '(a)' ) ' I/O error reading box coordinates.' stop end if end if write ( *, '(a)' ) ' ' write ( *, '(a,2f9.4)' ) & ' Lower left box corner = ', box_xy(1,1), box_xy(2,1) write ( *, '(a,2f9.4)' ) & ' Upper right box corner = ', box_xy(1,2), box_xy(2,2) ! ! CIRCLE x y r ! Plot a single circle of given size. ! else if ( s_eqi ( command(1:6), 'CIRCLE' ) ) then circle_requested = .true. command(1:6) = ' ' command = adjustl ( command ) if ( command(7:7) == '=' ) then command(7:7) = ' ' command = adjustl ( command ) end if if ( len_trim ( command ) /= 0 ) then first = 1 call s_to_r8 ( command(first:), circle_x, ierror, length ) first = first + length call s_to_r8 ( command(first:), circle_y, ierror, length ) first = first + length call s_to_r8 ( command(first:), circle_r, ierror, length ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter center (x,y), then radius r:' read ( *, *, iostat = ios ) circle_x, circle_y, circle_r if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_POINTS - Fatal error!' write ( *, '(a)' ) ' I/O error reading circle coordinates.' stop end if end if write ( *, '(a)' ) ' ' write ( *, '(a,2f9.4)' ) ' Circle center = ', circle_x, circle_y write ( *, '(a,2f9.4)' ) ' Circle radius = ', circle_r ! ! DASH ! Create a dash plot of connected points. ! else if ( s_eqi ( command(1:3), 'DAS' ) ) then if ( len_trim ( file_out_name_user ) == 0 ) then call file_name_inc ( file_out_name_default ) file_out_name = file_out_name_default else file_out_name = file_out_name_user file_out_name_user = ' ' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Write dash data to "' // & trim ( file_out_name ) // '".' call plot_size ( dim_num, box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, x_index, y_index, plot_min, plot_max ) call dash_plot ( box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, connect, file_out_name, & dim_num, node_num, coord, plot_min, plot_max, x_index, y_index, & shadow, title ) ! ! DELAUNAY ! Create a Delaunay plot. ! else if ( s_eqi ( command(1:3), 'DEL' ) ) then if ( len_trim ( file_out_name_user ) == 0 ) then call file_name_inc ( file_out_name_default ) file_out_name = file_out_name_default else file_out_name = file_out_name_user file_out_name_user = ' ' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Create Delaunay plot "' // & trim ( file_out_name ) // '".' call plot_size ( dim_num, box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, x_index, y_index, plot_min, plot_max ) call delaunay_plot ( box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, file_out_name, & dim_num, node_num, coord, plot_min, plot_max, x_index, y_index, & shadow, title ) ! ! DOT ! Create a dot plot of points. ! else if ( s_eqi ( command(1:3), 'DOT' ) ) then if ( len_trim ( file_out_name_user ) == 0 ) then call file_name_inc ( file_out_name_default ) file_out_name = file_out_name_default else file_out_name = file_out_name_user file_out_name_user = ' ' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Write dot data to "' // & trim ( file_out_name ) // '".' call plot_size ( dim_num, box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, x_index, y_index, plot_min, plot_max ) call dot_plot ( box_requested, box_xy, circle_requested, circle_x, & circle_y, circle_r, file_out_name, dim_num, & node_num, coord, plot_min, plot_max, x_index, y_index, & shadow, title, marker_size ) ! ! H ! else if ( s_eqi ( command(1:1), 'H' ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BALLOON a "balloon" plot of points.' 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)' ) 'READ filename read another input file.' write ( *, '(a)' ) 'OUTPUT filename name the next output file.' write ( *, '(a)' ) 'BOX x1 y1 x2 y2 draw a single box in the plot.' write ( *, '(a)' ) 'NOBOX do not draw a box in the plot.' write ( *, '(a)' ) 'CIRCLE x y r draw a single circle in the plot.' write ( *, '(a)' ) 'NOCIRCLE do not draw a circle in the plot.' write ( *, '(a)' ) 'LINE x1 y1 x2 y2 draw a line in the plot.' write ( *, '(a)' ) 'NOLINE do not draw lines in the plot.' write ( *, '(a)' ) 'MARKER_SIZE = n Specify marker size.' write ( *, '(a,i6)' ) ' Current value = ', marker_size write ( *, '(a)' ) 'SHADOW mark X and Y axes for each point.' write ( *, '(a)' ) 'NOSHADOW cancel shadow command.' write ( *, '(a)' ) 'TITLE enter a title for the next plot.' write ( *, '(a)' ) 'X = n specify data index to use as X.' write ( *, '(a,i6)' ) ' Current value = ', x_index write ( *, '(a)' ) 'Y = n specify data index to use as Y.' write ( *, '(a,i6)' ) ' Current value = ', y_index write ( *, '(a)' ) '# comment make a comment.' ! ! 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 ( s_eqi ( command(1:1), 'K' ) ) then if ( len_trim ( file_out_name_user ) == 0 ) then call file_name_inc ( file_out_name_default ) file_out_name = file_out_name_default else file_out_name = file_out_name_user file_out_name_user = ' ' end if tag(1:node_num) = nint ( coord(3,1:node_num) ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Create a K-Means plot "' // & trim ( file_out_name ) // '".' call plot_size ( dim_num, box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, x_index, y_index, plot_min, plot_max ) call kmeans_plot ( box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, file_out_name, dim_num, & node_num, coord, plot_min, plot_max, x_index, y_index, connect, & tag, shadow, title ) ! ! LINE x1, y1, x2, y2 ! specify a line to be drawn. ! else if ( s_eqi ( command(1:4), 'LINE' ) ) then line_num = line_num + 1 command(1:4) = ' ' command = adjustl ( command ) if ( command(5:5) == '=' ) then command(5:5) = ' ' command = adjustl ( command ) end if if ( len_trim ( command ) /= 0 ) then first = 1 call s_to_r8 ( command(first:), line_xy(1,1,line_num), ierror, length ) first = first + length call s_to_r8 ( command(first:), line_xy(2,1,line_num), ierror, length ) first = first + length call s_to_r8 ( command(first:), line_xy(1,2,line_num), ierror, length ) first = first + length call s_to_r8 ( command(first:), line_xy(2,2,line_num), ierror, length ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter first (x,y), then second (x,y):' read ( *, *, iostat = ios ) line_xy(1:2,1:2,line_num) if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_POINTS - Fatal error!' write ( *, '(a)' ) ' I/O error reading line coordinates.' stop end if end if write ( *, '(a)' ) ' ' write ( *, '(a,2f9.4)' ) & ' Line begins at ', line_xy(1,1,line_num), line_xy(2,1,line_num) write ( *, '(a,2f9.4)' ) & ' Line ends at ', line_xy(1,2,line_num), line_xy(2,2,line_num) ! ! MARKER_SIZE = n ! else if ( s_eqi ( command(1:11), 'MARKER_SIZE' ) ) then call s_blank_delete ( command ) if ( command(12:12) == '=' ) then call s_to_i4 ( command(13:), marker_size, ierror, length ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter new value for marker size:' read ( *, * ) marker_size end if write ( *, '(a,i6)' ) 'Marker size set to ', marker_size ! ! NOBOX ! specify no box to be drawn. ! else if ( s_eqi ( command(1:5), 'NOBOX' ) ) then box_requested = .false. ! ! NOCIRCLE ! specify no circle to be drawn. ! else if ( s_eqi ( command(1:8), 'NOCIRCLE' ) ) then circle_requested = .false. ! ! NOLINE ! specify no lines to be drawn. ! else if ( s_eqi ( command(1:6), 'NOLINE' ) ) then line_num = 0 ! ! NOSHADOW = Turn "shadow" OFF ! else if ( s_eqi ( command(1:4), 'NOSH' ) ) then shadow = .false. write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SHADOWING is turned OFF.' ! ! OUTPUT filename ! else if ( s_eqi ( command(1:6), 'OUTPUT' ) ) then call s_blank_delete ( command ) if ( command(7:7) == '=' ) then file_out_name_user = command(8:) else file_out_name_user = command(7:) end if if ( len_trim ( file_out_name_user ) == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter the name of the next output file.' read ( *, '(a)' ) file_out_name_user end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' Next output file is ' // trim ( file_out_name_user ) ! ! Q = Quit ! else if ( s_eqi ( command(1:1), 'Q' ) ) then exit ! ! READ ! else if ( s_eqi ( command(1:4), 'READ' ) ) then call s_blank_delete ( command ) if ( command(5:5) == '=' ) then file_in_name = command(6:) else file_in_name = command(5:) end if if ( len_trim ( file_in_name ) == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter the name of the new input file.' read ( *, '(a)' ) file_in_name end if if ( allocated ( connect ) ) then deallocate ( connect ) end if if ( allocated ( coord ) ) then deallocate ( coord ) end if if ( allocated ( coord_max ) ) then deallocate ( coord_max ) end if if ( allocated ( coord_min ) ) then deallocate ( coord_min ) end if if ( allocated ( radius ) ) then deallocate ( radius ) end if if ( allocated ( tag ) ) then deallocate ( tag ) end if dim_num = 0 node_num = 0 call data_read ( dim_num, file_in_name, node_num ) ! ! SH = Shadow ! else if ( s_eqi ( command(1:2), 'SH' ) ) then shadow = .true. write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SHADOWING is turned ON.' write ( *, '(a)' ) ' Each point makes a small mark on X and Y axis.' ! ! TH = Thin ! else if ( s_eqi ( command(1:2), 'TH' ) ) then tol = 0.01D+00 * maxval ( coord_max(1:dim_num) - coord_min(1:dim_num) ) call points_thin ( connect, dim_num, node_num, coord, tol, node_num2 ) node_num = node_num2 ! ! TITLE ! else if ( s_eqi ( command(1:5), 'TITLE' ) ) then next = 6 length = len_trim ( command ) do next = next + 1 if ( length < next ) then exit end if if ( command(next:next) /= ' ' .and. command(next:next) /= '=' ) then exit end if end do if ( next <= length ) then title = command(next:) else title = ' ' end if if ( len_trim ( title ) == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter the title:' read ( *, '(a)' ) title end if title = adjustl ( title ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TITLE:' write ( *, '(a)' ) ' "' // trim ( title ) // '"' ! ! TV ! Create a triangulated Voronoi plot. ! else if ( s_eqi ( command(1:2), 'TV' ) ) then if ( len_trim ( file_out_name_user ) == 0 ) then call file_name_inc ( file_out_name_default ) file_out_name = file_out_name_default else file_out_name = file_out_name_user file_out_name_user = ' ' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Create triangulated Voronoi plot "' & // trim ( file_out_name ) // '".' call plot_size ( dim_num, box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, x_index, y_index, plot_min, plot_max ) call trivor_plot ( box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, file_out_name, dim_num, & node_num, coord, plot_min, plot_max, x_index, y_index, & shadow, title ) ! ! VORONOI ! Create a Voronoi plot. ! else if ( s_eqi ( command(1:3), 'VOR' ) ) then if ( len_trim ( file_out_name_user ) == 0 ) then call file_name_inc ( file_out_name_default ) file_out_name = file_out_name_default else file_out_name = file_out_name_user file_out_name_user = ' ' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Create Voronoi plot "' // & trim ( file_out_name ) // '".' call plot_size ( dim_num, box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, x_index, y_index, plot_min, plot_max ) call voronoi_plot ( box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, file_out_name, dim_num, & node_num, coord, plot_min, plot_max, x_index, y_index, & shadow, title ) ! ! X_INDEX ! else if ( s_eqi ( command(1:1), 'X' ) ) then call s_blank_delete ( command ) if ( command(2:2) == '=' ) then call s_to_i4 ( command(3:), x_index, ierror, length ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter new value for X index:' read ( *, * ) x_index end if write ( *, '(a,i6)' ) 'X index set to ', x_index ! ! Y_INDEX ! else if ( s_eqi ( command(1:1), 'Y' ) ) then call s_blank_delete ( command ) if ( command(2:2) == '=' ) then call s_to_i4 ( command(3:), y_index, ierror, length ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter new value for Y index:' read ( *, * ) y_index end if write ( *, '(a,i6)' ) 'Y index set to ', y_index ! ! Unrecognized command. ! else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_POINTS: Unrecognized command:' write ( *, '(a)' ) ' "'// trim ( command ) // '"' end if end do ! ! Deallocate arrays. ! if ( allocated ( connect ) ) then deallocate ( connect ) end if if ( allocated ( coord ) ) then deallocate ( coord ) end if if ( allocated ( coord_max ) ) then deallocate ( coord_max ) end if if ( allocated ( coord_min ) ) then deallocate ( coord_min ) end if if ( allocated ( radius ) ) then deallocate ( radius ) end if if ( allocated ( tag ) ) then deallocate ( tag ) end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_POINTS' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine angle_contains_ray_2d ( inside, x1, y1, x2, y2, x3, y3, x, y ) !*****************************************************************************80 ! !! 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 ( kind = 8 ) X1, Y1, X2, Y2, X3, Y3, the X and Y coordinates of ! the angle. ! ! Input, real ( kind = 8 ) 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 ( kind = 8 ) a1 real ( kind = 8 ) a2 real ( kind = 8 ) angle_deg_2d logical inside real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) y real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) 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 ) !*****************************************************************************80 ! !! 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: ! ! 10 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) 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 ( kind = 8 ) 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 ( kind = 8 ) angle_deg_2d real ( kind = 8 ) angle_rad_2d real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) radians_to_degrees real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) y real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 x = ( x1 - x2 ) * ( x3 - x2 ) + ( y1 - y2 ) * ( y3 - y2 ) y = ( x1 - x2 ) * ( y3 - y2 ) - ( y1 - y2 ) * ( x3 - x2 ) if ( x == 0.0D+00 .and. y == 0.0D+00 ) then angle_deg_2d = 0.0D+00 else angle_rad_2d = atan2 ( y, x ) if ( angle_rad_2d < 0.0D+00 ) then angle_rad_2d = angle_rad_2d + 2.0D+00 * pi end if angle_deg_2d = radians_to_degrees ( angle_rad_2d ) end if return end subroutine angle_to_rgb ( angle, r, g, b ) !*****************************************************************************80 ! !! ANGLE_TO_RGB returns a color on the perimeter of the color hexagon. ! ! Modified: ! ! 13 December 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) ANGLE, the angle in the color hexagon. ! The sextants are defined by the following points: ! 0 degrees, 1, 0, 0, red; ! 60 degrees, 1, 1, 0, yellow; ! 120 degrees, 0, 1, 0, green; ! 180 degrees, 0, 1, 1, cyan; ! 240 degrees, 0, 0, 1, blue; ! 300 degrees, 1, 0, 1, magenta. ! ! Output, real ( kind = 8 ) R, G, B, RGB specifications for the color ! that lies at the given angle, on the perimeter of the color hexagon. One ! value will be 1, and one value will be 0. ! implicit none real ( kind = 8 ) angle real ( kind = 8 ) angle2 real ( kind = 8 ) b real ( kind = 8 ) g real ( kind = 8 ), parameter :: degrees_to_radians = & 3.14159265D+00 / 180.0D+00 real ( kind = 8 ) r angle = mod ( angle, 360.0D+00 ) if ( angle < 0.0D+00 ) then angle = angle + 360.0D+00 end if if ( angle <= 60.0D+00 ) then angle2 = degrees_to_radians * 3.0D+00 * angle / 4.0D+00 r = 1.0D+00 g = tan ( angle2 ) b = 0.0D+00 else if ( angle <= 120.0D+00 ) then angle2 = degrees_to_radians * 3.0D+00 * angle / 4.0D+00 r = cos ( angle2 ) / sin ( angle2 ) g = 1.0D+00 b = 0.0D+00 else if ( angle <= 180.0D+00 ) then angle2 = degrees_to_radians * 3.0D+00 * ( angle - 120.0D+00 ) / 4.0D+00 r = 0.0D+00 g = 1.0D+00 b = tan ( angle2 ) else if ( angle <= 240.0D+00 ) then angle2 = degrees_to_radians * 3.0D+00 * ( angle - 120.0D+00 ) / 4.0D+00 r = 0.0D+00 g = cos ( angle2 ) / sin ( angle2 ) b = 1.0D+00 else if ( angle <= 300.0D+00 ) then angle2 = degrees_to_radians * 3.0D+00 * ( angle - 240.0D+00 ) / 4.0D+00 r = tan ( angle2 ) g = 0.0D+00 b = 1.0D+00 else if ( angle <= 360.0D+00 ) then angle2 = degrees_to_radians * 3.0D+00 * ( angle - 240.0D+00 ) / 4.0D+00 r = 1.0D+00 g = 0.0D+00 b = cos ( angle2 ) / sin ( angle2 ) end if return end subroutine balloon_plot ( box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, filedot_name, dim_num, & node_num, coord, plot_min, plot_max, x, y, shadow, title, radius ) !*****************************************************************************80 ! !! BALLOON_PLOT plots points with maximal nonintersecting circles. ! ! Modified: ! ! 28 March 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical BOX_REQUESTED, is true if the user has specified ! a box to be drawn. ! ! Input, real ( kind = 8 ) BOX_XY(2,2), the coordinates of the lower left and ! upper right corners of the requested box. ! ! Input, character ( len = * ) FILEDOT_NAME, the name of the dot file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer NODE_NUM, the number of points. ! ! Input, real ( kind = 8 ) COORD(DIM_NUM,NODE_NUM), the point coordinates. ! ! Input, real ( kind = 8 ) PLOT_MIN(2), PLOT_MAX(2), the minimum and maximum ! X and Y values to plot. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! ! Input, logical SHADOW, is true if the user would like the ! X and Y axes to be marked at the shadow of each point. ! ! Input, character ( len = * ) TITLE, a title for the plot. ! ! Input, real ( kind = 8 ) RADIUS(NODE_NUM), the radius of the circle ! around each point. ! implicit none integer dim_num integer node_num real ( kind = 8 ) b real ( kind = 8 ) bitx real ( kind = 8 ) bity logical box_requested real ( kind = 8 ) box_xy(2,2) logical circle_requested real ( kind = 8 ) circle_r real ( kind = 8 ) circle_x real ( kind = 8 ) circle_y real ( kind = 8 ) coord(dim_num,node_num) character ( len = * ) filedot_name integer filedot_unit real ( kind = 8 ) g integer i integer ierror real ( kind = 8 ) plot_max(2) real ( kind = 8 ) plot_min(2) real ( kind = 8 ) r real ( kind = 8 ) radius(node_num) logical shadow character ( len = * ) title integer x real ( kind = 8 ) xvec(4) integer y real ( kind = 8 ) yvec(4) call get_unit ( filedot_unit ) call ps_file_open ( filedot_name, filedot_unit, ierror ) call eps_file_head ( filedot_name, 36, 36, 576, 756 ) call ps_page_head ( plot_min(1), plot_min(2), plot_max(1), plot_max(2) ) ! ! Initialize the line width to 1. ! call ps_line_width ( 1 ) ! ! Print title, if requested. ! if ( 0 < len_trim ( title ) ) then call ps_font_size ( 0.30D+00 ) call ps_moveto ( plot_min(1), plot_max(2) ) call ps_label ( title ) end if ! ! Define a PostScript clipping box. ! xvec(1:4) = (/ plot_min(1), plot_max(1), plot_max(1), plot_min(1) /) yvec(1:4) = (/ plot_min(2), plot_min(2), plot_max(2), plot_max(2) /) call ps_clip ( 4, xvec, yvec ) ! ! Because filled objects can obscure lines, draw them first. ! r = 0.75D+00 g = 0.75D+00 b = 1.00D+00 call ps_color_fill_set ( r, g, b ) do i = 1, node_num call ps_circle_fill ( coord(x,i), coord(y,i), radius(i) ) end do r = 0.0D+00 g = 0.0D+00 b = 0.0D+00 call ps_color_fill_set ( r, g, b ) ! ! Shadow the points, if requested. ! if ( shadow ) then bitx = 0.025D+00 * ( plot_max ( 1 ) - plot_min ( 1 ) ) bity = 0.025D+00 * ( plot_max ( 2 ) - plot_min ( 2 ) ) do i = 1, node_num call ps_line ( plot_min(1), coord(y,i), plot_min(1) + bitx, coord(y,i) ) call ps_line ( coord(x,i), plot_min(2), coord(x,i), plot_min(2) + bity ) end do end if ! ! Draw a user-specified box, if requested. ! if ( box_requested ) then call ps_comment ( 'User-requested box:' ) call ps_line ( box_xy(1,1), box_xy(2,1), box_xy(1,2), box_xy(2,1) ) call ps_line ( box_xy(1,2), box_xy(2,1), box_xy(1,2), box_xy(2,2) ) call ps_line ( box_xy(1,2), box_xy(2,2), box_xy(1,1), box_xy(2,2) ) call ps_line ( box_xy(1,1), box_xy(2,2), box_xy(1,1), box_xy(2,1) ) end if ! ! Draw a user-specified circle, if requested. ! if ( circle_requested ) then call ps_comment ( 'User-requested circle:' ) call ps_circle ( circle_x, circle_y, circle_r ) end if r = 0.0D+00 g = 0.0D+00 b = 0.0D+00 call ps_color_fill_set ( r, g, b ) if ( node_num <= 350 ) then do i = 1, node_num call ps_mark_disk ( coord(x,i), coord(y,i) ) end do else do i = 1, node_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 subroutine box_clip_line_2d ( xmin, ymin, xmax, ymax, x1, y1, x2, y2, x3, y3, & x4, y4, ival ) !*****************************************************************************80 ! !! BOX_CLIP_LINE_2D uses a box to clip a line segment in 2D. ! ! Discussion: ! ! The box is assumed to be a rectangle with sides aligned on coordinate ! axes. ! ! Modified: ! ! 18 June 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) XMIN, YMIN, XMAX, YMAX, the minimum and maximum ! X and Y values, which define the box. ! ! Input, real ( kind = 8 ) X1, Y1, X2, Y2, the coordinates of the ! endpoints of the line segment. ! ! Output, real ( kind = 8 ) X3, Y3, X4, Y4, the clipped coordinates. ! ! Output, integer IVAL: ! -1, no part of the line segment is within the box. ! 0, no clipping was necessary. The line segment is entirely within ! the box. ! 1, (X1,Y1) was clipped. ! 2, (X2,Y2) was clipped. ! 3, (X1,Y1) and (X2,Y2) were clipped. ! implicit none integer ival logical l1 logical l2 real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) x4 real ( kind = 8 ) xmax real ( kind = 8 ) xmin real ( kind = 8 ) y real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) y4 real ( kind = 8 ) ymax real ( kind = 8 ) ymin l1 = .false. l2 = .false. x3 = x1 y3 = y1 x4 = x2 y4 = y2 ! ! Require that XMIN <= X. ! if ( x3 < xmin .and. x4 < xmin ) then ival = -1 return end if if ( x3 < xmin .and. xmin <= x4 ) then x = xmin y = y3 + ( y4 - y3 ) * ( x - x3 ) / ( x4 - x3 ) x3 = x y3 = y l1 = .true. else if ( xmin <= x3 .and. x4 < xmin ) then x = xmin y = y3 + ( y4 - y3 ) * ( x - x3 ) / ( x4 - x3 ) x4 = x y4 = y l2 = .true. end if ! ! Require that X <= XMAX. ! if ( xmax < x3 .and. xmax < x4 ) then ival = -1 return end if if ( xmax < x3 .and. x4 <= xmax ) then x = xmax y = y3 + ( y4 - y3 ) * ( x - x3 ) / ( x4 - x3 ) x3 = x y3 = y l1 = .true. else if ( x3 <= xmax .and. xmax < x4 ) then x = xmax y = y3 + ( y4 - y3 ) * ( x - x3 ) / ( x4 - x3 ) x4 = x y4 = y l2 = .true. end if ! ! Require that YMIN <= Y. ! if ( y3 < ymin .and. y4 < ymin ) then ival = -1 return end if if ( y3 < ymin .and. ymin <= y4 ) then y = ymin x = x3 + ( x4 - x3 ) * ( y - y3 ) / ( y4 - y3 ) y3 = y x3 = x l1 = .true. else if ( ymin <= y3 .and. y4 < ymin ) then y = ymin x = x3 + ( x4 - x3 ) * ( y - y3 ) / ( y4 - y3 ) y4 = y x4 = x l2 = .true. end if ! ! Require that Y <= YMAX. ! if ( ymax < y3 .and. ymax < y4 ) then ival = -1 return end if if ( ymax < y3 .and. y4 <= ymax ) then y = ymax x = x3 + ( x4 - x3 ) * ( y - y3 ) / ( y4 - y3 ) y3 = y x3 = x l1 = .true. else if ( y3 <= ymax .and. ymax < y4 ) then y = ymax x = x3 + ( x4 - x3 ) * ( y - y3 ) / ( y4 - y3 ) y4 = y x4 = x l2 = .true. end if ival = 0 if ( l1 ) then ival = ival + 1 end if if ( l2 ) then ival = ival + 2 end if return end subroutine box_ray_int_2d ( xmin, ymin, xmax, ymax, xa, ya, xb, yb, xi, yi ) !*****************************************************************************80 ! !! 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 ( kind = 8 ) XMIN, YMIN, the lower left corner of the box. ! ! Input, real ( kind = 8 ) XMAX, YMAX, the upper right corner of the box. ! ! Input, real ( kind = 8 ) XA, YA, the origin of the ray, which should be ! inside the box. ! ! Input, real ( kind = 8 ) XB, YB, a second point on the ray. ! ! Output, real ( kind = 8 ) XI, YI, the point on the box intersected ! by the ray. ! implicit none logical inside integer ival integer side real ( kind = 8 ) xa real ( kind = 8 ) xb real ( kind = 8 ) xc real ( kind = 8 ) xd real ( kind = 8 ) xi real ( kind = 8 ) xmax real ( kind = 8 ) xmin real ( kind = 8 ) ya real ( kind = 8 ) yb real ( kind = 8 ) yc real ( kind = 8 ) yd real ( kind = 8 ) yi real ( kind = 8 ) ymax real ( kind = 8 ) 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 ) !*****************************************************************************80 ! !! 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 ) !*****************************************************************************80 ! !! 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 ) !*****************************************************************************80 ! !! 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 ) !*****************************************************************************80 ! !! 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 circle_points ( x0, y0, r, n, x, y ) !*****************************************************************************80 ! !! CIRCLE_POINTS returns N equally spaced points on a circle in 2D. ! ! Discussion: ! ! The first point is always ( X0 + R, Y0 ), and subsequent points ! proceed counterclockwise around the circle. ! ! Modified: ! ! 28 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X0, Y0, the coordinates of the center of ! the circle. ! ! Input, real ( kind = 8 ) R, the radius of the circle. ! ! Input, integer N, the number of points desired. N must be at least 1. ! ! Output, real ( kind = 8 ) X(N), Y(N), the coordinates of points ! on the circle. ! implicit none integer n real ( kind = 8 ) angle integer i real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) r real ( kind = 8 ) x0 real ( kind = 8 ) x(n) real ( kind = 8 ) y0 real ( kind = 8 ) y(n) do i = 1, n angle = ( 2.0D+00 * pi * real ( i - 1, kind = 8 ) ) & / real ( n, kind = 8 ) x(i) = x0 + r * cos ( angle ) y(i) = y0 + r * sin ( angle ) end do return end subroutine dash_plot ( box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, connect, file_name, dim_num, & node_num, coord, plot_min, plot_max, x, y, shadow, title ) !*****************************************************************************80 ! !! DASH_PLOT plots a set of points, and connects them. ! ! Modified: ! ! 28 March 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer CONNECT(NODE_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 NODE_NUM, the number of points. ! ! Input, real ( kind = 8 ) COORD(DIM_NUM,NODE_NUM), the point coordinates. ! ! Input, real ( kind = 8 ) PLOT_MIN(2), PLOT_MAX(2), the minimum and maximum ! X and Y values to plot. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! implicit none integer dim_num integer node_num ! real ( kind = 8 ) bitx real ( kind = 8 ) bity logical box_requested real ( kind = 8 ) box_xy(2,2) logical circle_requested real ( kind = 8 ) circle_r real ( kind = 8 ) circle_x real ( kind = 8 ) circle_y integer connect(node_num) real ( kind = 8 ) coord(dim_num,node_num) character ( len = * ) file_name integer file_unit integer i integer ierror real ( kind = 8 ) plot_max(2) real ( kind = 8 ) plot_min(2) logical shadow character ( len = * ) title integer x real ( kind = 8 ) xvec(4) integer y real ( kind = 8 ) yvec(4) call get_unit ( file_unit ) call ps_file_open ( file_name, file_unit, ierror ) call eps_file_head ( file_name, 36, 36, 576, 756 ) call ps_page_head ( plot_min(1), plot_min(2), plot_max(1), plot_max(2) ) ! ! Initialize the line width to 1. ! call ps_line_width ( 1 ) ! ! Print title, if requested. ! if ( 0 < len_trim ( title ) ) then call ps_font_size ( 0.30D+00 ) call ps_moveto ( plot_min(1), plot_max(2) ) call ps_label ( title ) end if ! ! Define a PostScript clipping box. ! xvec(1:4) = (/ plot_min(1), plot_max(1), plot_max(1), plot_min(1) /) yvec(1:4) = (/ plot_min(2), plot_min(2), plot_max(2), plot_max(2) /) call ps_clip ( 4, xvec, yvec ) ! ! Shadow the points, if requested. ! if ( shadow ) then bitx = 0.025D+00 * ( plot_max ( 1 ) - plot_min ( 1 ) ) bity = 0.025D+00 * ( plot_max ( 2 ) - plot_min ( 2 ) ) do i = 1, node_num call ps_line ( plot_min(1), coord(y,i), plot_min(1) + bitx, coord(y,i) ) call ps_line ( coord(x,i), plot_min(2), coord(x,i), plot_min(2) + bity ) end do end if ! ! Draw a user-specified box, if requested. ! if ( box_requested ) then call ps_comment ( 'User-requested box:' ) call ps_line ( box_xy(1,1), box_xy(2,1), box_xy(1,2), box_xy(2,1) ) call ps_line ( box_xy(1,2), box_xy(2,1), box_xy(1,2), box_xy(2,2) ) call ps_line ( box_xy(1,2), box_xy(2,2), box_xy(1,1), box_xy(2,2) ) call ps_line ( box_xy(1,1), box_xy(2,2), box_xy(1,1), box_xy(2,1) ) end if ! ! Draw a user-specified circle, if requested. ! if ( circle_requested ) then call ps_comment ( 'User-requested circle:' ) call ps_circle ( circle_x, circle_y, circle_r ) end if ! ! Draw the points, with big or little marks as appropriate. ! if ( node_num <= 350 ) then do i = 1, node_num call ps_mark_disk ( coord(x,i), coord(y,i) ) end do else do i = 1, node_num call ps_mark_point ( coord(x,i), coord(y,i) ) end do end if do i = 1, node_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(node_num) == 2 .or. connect(node_num) == 3 ) then call ps_line ( coord(x,node_num), coord(y,node_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 data_read ( dim_num, file_in_name, node_num ) !*****************************************************************************80 ! !! DATA_READ reads the point data from a file. ! ! Modified: ! ! 16 February 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer DIM_NUM, the number of spatial dimensions. ! ! Input, character ( len = * ) FILE_IN_NAME, the file to be read. ! ! Output, integer NODE_NUM, the number of points. ! use AA_point_data_module implicit none integer dim_num character ( len = * ) file_in_name integer i integer node_num write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DATA_READ' 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, node_num ) ! ! Now allocate some data. ! allocate ( connect(node_num) ) allocate ( coord(dim_num,node_num) ) allocate ( coord_min(dim_num) ) allocate ( coord_max(dim_num) ) allocate ( radius(node_num) ) allocate ( tag(node_num) ) ! ! Read the coordinates into COORD. ! call points_read ( connect, file_in_name, dim_num, node_num, coord ) ! ! Determine the range. ! do i = 1, dim_num coord_min(i) = minval ( coord(i,1:node_num) ) coord_max(i) = maxval ( coord(i,1:node_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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DATA_READ' write ( *, '(a)' ) ' The data has been read.' return end subroutine delaunay_plot ( box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, file_name, dim_num, & node_num, coord, plot_min, plot_max, x, y, shadow, title ) !*****************************************************************************80 ! !! DELAUNAY_PLOT plots the Delaunay triangulation of a pointset. ! ! Modified: ! ! 03 January 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical BOX_REQUESTED, is true if the user has specified ! a box to be drawn. ! ! Input, real ( kind = 8 ) BOX_XY(2,2), the coordinates of the lower left and ! upper right corners of the requested box. ! ! Input, character ( len = * ) FILE_NAME, the name of the output file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer NODE_NUM, the number of points. ! ! Input, real ( kind = 8 ) COORD(DIM_NUM,NODE_NUM), the point coordinates. ! ! Input, real ( kind = 8 ) PLOT_MIN(2), PLOT_MAX(2), the minimum and maximum ! X and Y values to plot. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! ! Input, logical SHADOW, is true if the user would like the ! X and Y axes to be marked at the shadow of each point. ! ! Input, character ( len = * ) TITLE, a title for the plot. ! implicit none integer dim_num integer node_num real ( kind = 8 ) bitx real ( kind = 8 ) bity logical box_requested real ( kind = 8 ) box_xy(2,2) logical circle_requested real ( kind = 8 ) circle_r real ( kind = 8 ) circle_x real ( kind = 8 ) circle_y real ( kind = 8 ) coord(dim_num,node_num) real ( kind = 8 ) coord2(2,node_num) logical, parameter :: debug = .false. character ( len = * ) file_name integer file_unit integer i integer ierror integer j1 integer j2 integer j3 integer triangle_node(3,2*node_num) real ( kind = 8 ) plot_max(2) real ( kind = 8 ) plot_min(2) logical shadow character ( len = * ) title integer triangle_neighbor(3,2*node_num) integer triangle_num integer x real ( kind = 8 ) xvec(4) real ( kind = 8 ) xx integer y real ( kind = 8 ) yvec(4) real ( kind = 8 ) yy ! ! Compute the Delaunay triangulation. ! coord2(1,1:node_num) = coord(x,1:node_num) coord2(2,1:node_num) = coord(y,1:node_num) call dtris2 ( node_num, coord2, triangle_num, triangle_node, triangle_neighbor ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DELAUNAY_PLOT:' write ( *, '(a,i6)' ) ' The number of triangles is ', triangle_num ! ! Print the triangulation. ! if ( debug ) then call triangulation_order3_print ( node_num, triangle_num, coord2, & triangle_node, triangle_neighbor ) 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, 36, 36, 576, 756 ) call ps_page_head ( plot_min(1), plot_min(2), plot_max(1), plot_max(2) ) ! ! Initialize the line width to 1. ! call ps_line_width ( 1 ) ! ! Print title, if requested. ! if ( 0 < len_trim ( title ) ) then call ps_font_size ( 0.30D+00 ) call ps_moveto ( plot_min(1), plot_max(2) ) call ps_label ( title ) end if ! ! Define a PostScript clipping box. ! xvec(1:4) = (/ plot_min(1), plot_max(1), plot_max(1), plot_min(1) /) yvec(1:4) = (/ plot_min(2), plot_min(2), plot_max(2), plot_max(2) /) call ps_clip ( 4, xvec, yvec ) ! ! Shadow the points, if requested. ! if ( shadow ) then bitx = 0.025D+00 * ( plot_max ( 1 ) - plot_min ( 1 ) ) bity = 0.025D+00 * ( plot_max ( 2 ) - plot_min ( 2 ) ) do i = 1, node_num call ps_line ( plot_min(1), coord(y,i), plot_min(1) + bitx, coord(y,i) ) call ps_line ( coord(x,i), plot_min(2), coord(x,i), plot_min(2) + bity ) end do end if ! ! Draw a user-specified box, if requested. ! if ( box_requested ) then call ps_comment ( 'User-requested box:' ) call ps_line ( box_xy(1,1), box_xy(2,1), box_xy(1,2), box_xy(2,1) ) call ps_line ( box_xy(1,2), box_xy(2,1), box_xy(1,2), box_xy(2,2) ) call ps_line ( box_xy(1,2), box_xy(2,2), box_xy(1,1), box_xy(2,2) ) call ps_line ( box_xy(1,1), box_xy(2,2), box_xy(1,1), box_xy(2,1) ) end if ! ! Draw a user-specified circle, if requested. ! if ( circle_requested ) then call ps_comment ( 'User-requested circle:' ) call ps_circle ( circle_x, circle_y, circle_r ) end if xx = plot_min(1) yy = plot_min(2) call ps_moveto ( xx, yy ) ! call ps_label ( 'Delaunay triangulation' ) if ( node_num <= 350 ) then do i = 1, node_num call ps_mark_disk ( coord(x,i), coord(y,i) ) end do else do i = 1, node_num call ps_mark_point ( coord(x,i), coord(y,i) ) end do end if ! ! Increase the line width for small problems. ! if ( triangle_num < 50 ) then call ps_line_width ( 3 ) else if ( triangle_num < 250 ) then call ps_line_width ( 2 ) else call ps_line_width ( 1 ) end if do i = 1, triangle_num j1 = triangle_node(1,i) j2 = triangle_node(2,i) j3 = triangle_node(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 function diaedg ( x0, y0, x1, y1, x2, y2, x3, y3 ) !*****************************************************************************80 ! !! 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 ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) X0, Y0, X1, Y1, X2, Y2, X3, Y3, the ! coordinates of the vertices of a quadrilateral, given in ! counter clockwise order. ! ! 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 ( kind = 8 ) ca real ( kind = 8 ) cb integer diaedg real ( kind = 8 ) dx10 real ( kind = 8 ) dx12 real ( kind = 8 ) dx30 real ( kind = 8 ) dx32 real ( kind = 8 ) dy10 real ( kind = 8 ) dy12 real ( kind = 8 ) dy30 real ( kind = 8 ) dy32 real ( kind = 8 ) s real ( kind = 8 ) tol real ( kind = 8 ) tola real ( kind = 8 ) tolb real ( kind = 8 ) x0 real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) y0 real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 tol = 100.0D+00 * epsilon ( tol ) dx10 = x1 - x0 dy10 = y1 - y0 dx12 = x1 - x2 dy12 = y1 - y2 dx30 = x3 - x0 dy30 = y3 - y0 dx32 = x3 - x2 dy32 = y3 - y2 tola = tol * max ( abs ( dx10 ), abs ( dy10 ), abs ( dx30 ), abs ( dy30 ) ) tolb = tol * max ( abs ( dx12 ), abs ( dy12 ), abs ( dx32 ), abs ( dy32 ) ) ca = dx10 * dx30 + dy10 * dy30 cb = dx12 * dx32 + dy12 * dy32 if ( tola < ca .and. tolb < cb ) then diaedg = -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 ( tola < s ) then diaedg = -1 else if ( s < -tola ) then diaedg = 1 else diaedg = 0 end if end if return end subroutine digit_inc ( c ) !*****************************************************************************80 ! !! 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 ) !*****************************************************************************80 ! !! 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 ( box_requested, box_xy, circle_requested, circle_x, & circle_y, circle_r, file_name, dim_num, & node_num, node_xy, plot_min, plot_max, x, y, shadow, title, marker_size ) !*****************************************************************************80 ! !! DOT_PLOT makes a plot of a set of points. ! ! Modified: ! ! 06 July 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical BOX_REQUESTED, is true if the user has specified ! a box to be drawn. ! ! Input, real ( kind = 8 ) BOX_XY(2,2), the coordinates of the lower left and ! upper right corners of the requested box. ! ! Input, character ( len = * ) FILE_NAME, the name of the dot file. ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, integer NODE_NUM, the number of points. ! ! Input, real ( kind = 8 ) NODE_XY(DIM_NUM,NODE_NUM), the point coordinates. ! ! Input, real ( kind = 8 ) PLOT_MIN(2), PLOT_MAX(2), the minimum and maximum ! X and Y values to plot. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! ! Input, logical SHADOW, is true if the user would like the ! X and Y axes to be marked at the shadow of each point. ! ! Input, character ( len = * ) TITLE, a title for the plot. ! ! Input, integer MARKER_SIZE, controls the size of the dots. ! A value of 3, 4, or 5 is typical. ! implicit none integer dim_num integer node_num real ( kind = 8 ) b real ( kind = 8 ) bitx real ( kind = 8 ) bity logical box_requested real ( kind = 8 ) box_xy(2,2) logical circle_requested real ( kind = 8 ) circle_r real ( kind = 8 ) circle_x real ( kind = 8 ) circle_y character ( len = 40 ) date_time integer delta character ( len = * ) file_name integer file_unit real ( kind = 8 ) g integer i integer ierror integer marker_size real ( kind = 8 ) node_xy(dim_num,node_num) real ( kind = 8 ) plot_max(2) real ( kind = 8 ) plot_min(2) real ( kind = 8 ) r logical shadow character ( len = * ) title integer x real ( kind = 8 ) x_max real ( kind = 8 ) x_min integer x_ps integer :: x_ps_max = 576 integer :: x_ps_max_clip = 594 integer :: x_ps_min = 36 integer :: x_ps_min_clip = 18 real ( kind = 8 ) x_scale real ( kind = 8 ) xvec(4) integer y real ( kind = 8 ) y_max real ( kind = 8 ) y_min integer y_ps integer :: y_ps_max = 666 integer :: y_ps_max_clip = 684 integer :: y_ps_min = 126 integer :: y_ps_min_clip = 108 real ( kind = 8 ) y_scale real ( kind = 8 ) yvec(4) call timestring ( date_time ) ! ! We need to do some figuring here, so that we can determine ! the range of the data, and hence the height and width ! of the piece of paper. ! x_max = maxval ( node_xy(x,1:node_num) ) x_min = minval ( node_xy(x,1:node_num) ) if ( box_requested ) then x_max = max ( x_max, box_xy(1,1) ) x_max = max ( x_max, box_xy(1,2) ) x_min = min ( x_min, box_xy(1,1) ) x_min = min ( x_min, box_xy(1,2) ) end if x_scale = x_max - x_min x_max = x_max + 0.05D+00 * x_scale x_min = x_min - 0.05D+00 * x_scale x_scale = x_max - x_min y_max = maxval ( node_xy(y,1:node_num) ) y_min = minval ( node_xy(y,1:node_num) ) if ( box_requested ) then y_max = max ( y_max, box_xy(2,1) ) y_max = max ( y_max, box_xy(2,2) ) y_min = min ( y_min, box_xy(2,1) ) y_min = min ( y_min, box_xy(2,2) ) end if y_scale = y_max - y_min y_max = y_max + 0.05D+00 * y_scale y_min = y_min - 0.05D+00 * y_scale y_scale = y_max - y_min if ( x_scale < y_scale ) then delta = nint ( real ( x_ps_max - x_ps_min, kind = 8 ) & * ( y_scale - x_scale ) / ( 2.0D+00 * y_scale ) ) x_ps_max = x_ps_max - delta x_ps_min = x_ps_min + delta x_ps_max_clip = x_ps_max_clip - delta x_ps_min_clip = x_ps_min_clip + delta x_scale = y_scale else if ( y_scale < x_scale ) then delta = nint ( real ( y_ps_max - y_ps_min, kind = 8 ) & * ( x_scale - y_scale ) / ( 2.0D+00 * x_scale ) ) y_ps_max = y_ps_max - delta y_ps_min = y_ps_min + delta y_ps_max_clip = y_ps_max_clip - delta y_ps_min_clip = y_ps_min_clip + delta y_scale = x_scale end if call get_unit ( file_unit ) call ps_file_open ( file_name, file_unit, ierror ) call eps_file_head ( file_name, x_ps_min, y_ps_min, x_ps_max, y_ps_max ) call ps_page_head ( x_min, y_min, x_max, y_max ) ! ! Initialize the line width to 1. ! call ps_line_width ( 1 ) ! ! Print title, if requested. ! if ( 0 < len_trim ( title ) ) then call ps_font_size ( 0.30D+00 ) call ps_moveto ( plot_min(1), plot_max(2) ) call ps_label ( title ) end if ! ! Define a PostScript clipping box. ! xvec(1) = x_min xvec(2) = x_max xvec(3) = x_max xvec(4) = x_min yvec(1) = y_min yvec(2) = y_min yvec(3) = y_max yvec(4) = y_max call ps_clip ( 4, xvec, yvec ) ! ! Shadow the points, if requested. ! if ( shadow ) then bitx = 0.025D+00 * ( plot_max ( 1 ) - plot_min ( 1 ) ) bity = 0.025D+00 * ( plot_max ( 2 ) - plot_min ( 2 ) ) do i = 1, node_num call ps_line ( plot_min(1), node_xy(y,i), plot_min(1) + bitx, & node_xy(y,i) ) call ps_line ( node_xy(x,i), plot_min(2), node_xy(x,i), & plot_min(2) + bity ) end do end if ! ! Draw a user-specified box, if requested. ! if ( box_requested ) then call ps_comment ( 'User-requested box:' ) call ps_line ( box_xy(1,1), box_xy(2,1), box_xy(1,2), box_xy(2,1) ) call ps_line ( box_xy(1,2), box_xy(2,1), box_xy(1,2), box_xy(2,2) ) call ps_line ( box_xy(1,2), box_xy(2,2), box_xy(1,1), box_xy(2,2) ) call ps_line ( box_xy(1,1), box_xy(2,2), box_xy(1,1), box_xy(2,1) ) end if ! ! Draw a user-specified circle, if requested. ! if ( circle_requested ) then call ps_comment ( 'User-requested circle:' ) call ps_circle ( circle_x, circle_y, circle_r ) end if call ps_comment ( 'Set color to bluish green:' ) call ps_setting_int ( 'SET', 'MARKER_SIZE', marker_size ) r = 0.000D+00 g = 0.750D+00 b = 0.150D+00 call ps_color_fill_set ( r, g, b ) if ( node_num <= 1100 ) then do i = 1, node_num call ps_mark_disk ( node_xy(x,i), node_xy(y,i) ) end do else do i = 1, node_num call ps_mark_point ( node_xy(x,i), node_xy(y,i) ) end do end if call ps_page_tail ( ) call eps_file_tail ( ) call ps_file_close ( file_unit ) return end subroutine dtris2 ( node_num, node_xy, triangle_num, triangle_node, & triangle_neighbor ) !*****************************************************************************80 ! !! DTRIS2 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: ! ! 25 August 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input/output, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates ! of the nodes. On output, the vertices have been sorted into ! dictionary order. ! ! Output, integer TRIANGLE_NUM, the number of triangles in the triangulation; ! TRIANGLE_NUM is equal to 2*NODE_NUM - NB - 2, where NB is the number ! of boundary vertices. ! ! Output, integer TRIANGLE_NODE(3,TRIANGLE_NUM), the nodes that make up each ! triangle. The elements are indices of P. The vertices of the triangles ! are in counter clockwise order. ! ! Output, integer TRIANGLE_NEIGHBOR(3,TRIANGLE_NUM), the triangle neighbor ! list. Positive elements are indices of TIL; negative elements are used ! for links of a counter clockwise linked list of boundary edges; ! LINK = -(3*I + J-1) where I, J = triangle, edge index; ! TRIANGLE_NEIGHBOR(J,I) refers to the neighbor along edge from vertex J ! to J+1 (mod 3). ! implicit none integer, parameter :: dim_num = 2 integer node_num real ( kind = 8 ) cmax integer e integer i integer ierr integer indx(node_num) integer j integer k integer l integer ledg integer lr integer lrline integer ltri integer m integer m1 integer m2 integer n real ( kind = 8 ) node_xy(dim_num,node_num) integer redg integer rtri integer stack(node_num) integer t real ( kind = 8 ) tol integer top integer triangle_neighbor(3,node_num*2) integer triangle_num integer triangle_node(3,node_num*2) tol = 100.0D+00 * epsilon ( tol ) ierr = 0 ! ! Sort the vertices by increasing (x,y). ! call r82vec_sort_heap_index_a ( node_num, node_xy, indx ) call r82vec_permute ( node_num, node_xy, indx ) ! ! Make sure that the data nodes are "reasonably" distinct. ! m1 = 1 do i = 2, node_num m = m1 m1 = i k = 0 do j = 1, dim_num cmax = max ( abs ( node_xy(j,m) ), abs ( node_xy(j,m1) ) ) if ( tol * ( cmax + 1.0D+00 ) & < abs ( node_xy(j,m) - node_xy(j,m1) ) ) then k = j exit end if end do if ( k == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - 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)' ) ' NODE_XY(M) = ', node_xy(1:dim_num,m) write ( *, '(a,2g14.6)' ) ' NODE_XY(M1) = ', node_xy(1:dim_num,m1) ierr = 224 stop end if end do ! ! Starting from nodes M1 and M2, search for a third point M that ! makes a "healthy" triangle (M1,M2,M) ! m1 = 1 m2 = 2 j = 3 do if ( node_num < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' ierr = 225 stop end if m = j lr = lrline ( node_xy(1,m), node_xy(2,m), node_xy(1,m1), & node_xy(2,m1), node_xy(1,m2), node_xy(2,m2), 0.0D+00 ) if ( lr /= 0 ) then exit end if j = j + 1 end do ! ! Set up the triangle information for (M1,M2,M), and for any other ! triangles you created because points were collinear with M1, M2. ! triangle_num = j - 2 if ( lr == -1 ) then triangle_node(1,1) = m1 triangle_node(2,1) = m2 triangle_node(3,1) = m triangle_neighbor(3,1) = -3 do i = 2, triangle_num m1 = m2 m2 = i+1 triangle_node(1,i) = m1 triangle_node(2,i) = m2 triangle_node(3,i) = m triangle_neighbor(1,i-1) = -3 * i triangle_neighbor(2,i-1) = i triangle_neighbor(3,i) = i - 1 end do triangle_neighbor(1,triangle_num) = -3 * triangle_num - 1 triangle_neighbor(2,triangle_num) = -5 ledg = 2 ltri = triangle_num else triangle_node(1,1) = m2 triangle_node(2,1) = m1 triangle_node(3,1) = m triangle_neighbor(1,1) = -4 do i = 2, triangle_num m1 = m2 m2 = i+1 triangle_node(1,i) = m2 triangle_node(2,i) = m1 triangle_node(3,i) = m triangle_neighbor(3,i-1) = i triangle_neighbor(1,i) = -3 * i - 3 triangle_neighbor(2,i) = i - 1 end do triangle_neighbor(3,triangle_num) = -3 * triangle_num triangle_neighbor(2,1) = -3 * triangle_num - 2 ledg = 2 ltri = 1 end if ! ! Insert the vertices one at a time from outside the convex hull, ! determine visible boundary edges, and apply diagonal edge swaps until ! Delaunay triangulation of vertices (so far) is obtained. ! top = 0 do i = j+1, node_num m = i m1 = triangle_node(ledg,ltri) if ( ledg <= 2 ) then m2 = triangle_node(ledg+1,ltri) else m2 = triangle_node(1,ltri) end if lr = lrline ( node_xy(1,m), node_xy(2,m), node_xy(1,m1), & node_xy(2,m1), node_xy(1,m2), node_xy(2,m2), 0.0D+00 ) if ( 0 < lr ) then rtri = ltri redg = ledg ltri = 0 else l = -triangle_neighbor(ledg,ltri) rtri = l / 3 redg = mod(l,3) + 1 end if call vbedg ( node_xy(1,m), node_xy(2,m), node_num, node_xy, triangle_num, & triangle_node, triangle_neighbor, ltri, ledg, rtri, redg ) n = triangle_num + 1 l = -triangle_neighbor(ledg,ltri) do t = l / 3 e = mod ( l, 3 ) + 1 l = -triangle_neighbor(e,t) m2 = triangle_node(e,t) if ( e <= 2 ) then m1 = triangle_node(e+1,t) else m1 = triangle_node(1,t) end if triangle_num = triangle_num + 1 triangle_neighbor(e,t) = triangle_num triangle_node(1,triangle_num) = m1 triangle_node(2,triangle_num) = m2 triangle_node(3,triangle_num) = m triangle_neighbor(1,triangle_num) = t triangle_neighbor(2,triangle_num) = triangle_num - 1 triangle_neighbor(3,triangle_num) = triangle_num + 1 top = top + 1 if ( node_num < top ) then ierr = 8 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' write ( *, '(a)' ) ' Stack overflow.' stop end if stack(top) = triangle_num if ( t == rtri .and. e == redg ) then exit end if end do triangle_neighbor(ledg,ltri) = -3 * n - 1 triangle_neighbor(2,n) = -3 * triangle_num - 2 triangle_neighbor(3,triangle_num) = -l ltri = n ledg = 2 call swapec ( m, top, ltri, ledg, node_num, node_xy, triangle_num, & triangle_node, triangle_neighbor, stack, ierr ) if ( ierr /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DTRIS2 - Fatal error!' write ( *, '(a)' ) ' Error return from SWAPEC.' stop end if end do ! ! Now account for the sorting that we did. ! do i = 1, 3 do j = 1, triangle_num triangle_node(i,j) = indx ( triangle_node(i,j) ) end do end do call perm_inv ( node_num, indx ) call r82vec_permute ( node_num, node_xy, indx ) return end subroutine eps_file_head ( file_name, x_ps_min, y_ps_min, x_ps_max, & y_ps_max ) !*****************************************************************************80 ! !! EPS_FILE_HEAD writes header information to an encapsulated PostScript file. ! ! Discussion: ! ! The file should contain the description of only one page, but this ! is not currently checked. ! ! Modified: ! ! 12 April 2005 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Henry McGilton, Mary Campione, ! PostScript by Example, ! Addison-Wesley, ! ISBN: 0-201-63228-4 ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the output file. ! ! Input, integer X_PS_MIN, Y_PS_MIN, X_PS_MAX, Y_PS_MAX, the minimum ! and maximum X and Y values of the data, in PostScript units. Any data ! that lies outside this range will not show up properly. A reasonable ! set of values might be 0, 0, 612, 792, or, for a half inch margin, ! 36, 36, 576, 756. ! implicit none character ( len = 8 ) date character ( len = * ) file_name real ( kind = 8 ) line_blue real ( kind = 8 ) line_green real ( kind = 8 ) line_red integer state integer unit integer x_ps_max integer x_ps_min integer y_ps_max integer y_ps_min ! ! Determine if the PostScript state is acceptable. ! call ps_setting_int ( 'GET', 'STATE', state ) if ( state /= 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'EPS_FILE_HEAD - Fatal error!' write ( *, '(a,i9)' ) ' PostScript state is ', state write ( *, '(a)' ) ' PostScript state 1 is required.' return end if ! ! Initialization ! call ps_default ( ) ! ! Get the unit number. ! call ps_setting_int ( 'GET', 'UNIT', unit ) call date_and_time ( date ) ! ! Write the prolog. ! write ( unit, '(a)' ) '%!PS-Adobe-3.0 EPSF-3.0' write ( unit, '(a)' ) '%%Creator: ps_write.f90' write ( unit, '(a)' ) '%%Title: ' // trim ( file_name ) write ( unit, '(a)' ) '%%CreationDate: '// trim ( date ) write ( unit, '(a)' ) '%%Pages: 1' write ( unit, '(a,4i6)' ) '%%BoundingBox:', & x_ps_min, y_ps_min, x_ps_max, y_ps_max write ( unit, '(a)' ) '%%Document-Fonts: Times-Roman' write ( unit, '(a)' ) '%%LanguageLevel: 1' write ( unit, '(a)' ) '%%EndComments' write ( unit, '(a)' ) '%%BeginProlog' write ( unit, '(a)' ) '/inch {72 mul} def' write ( unit, '(a)' ) '%%EndProlog' ! ! Set the font. ! write ( unit, '(a)' ) '/Times-Roman findfont' write ( unit, '(a)' ) '1.00 inch scalefont' write ( unit, '(a)' ) 'setfont' ! ! Set the line color. ! line_red = 0.0D+00 line_green = 0.0D+00 line_blue = 0.0D+00 call ps_color_line ( 'SET', line_red, line_green, line_blue ) ! ! Reset the state. ! state = 2 call ps_setting_int ( 'SET', 'STATE', state ) return end subroutine eps_file_tail ( ) !*****************************************************************************80 ! !! EPS_FILE_TAIL writes trailer information to an encapsulated PostScript file. ! ! Discussion: ! ! Looks like that penultimate 'end' line is not wanted, so I commented ! it out. ! ! Modified: ! ! 05 March 2002 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Henry McGilton, Mary Campione, ! PostScript by Example, ! Addison-Wesley, ! ISBN: 0-201-63228-4 ! ! Parameters: ! ! None ! implicit none integer num_pages integer state integer unit ! ! Determine if the PostScript state is acceptable. ! call ps_setting_int ( 'GET', 'STATE', state ) if ( state == 3 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'EPS_FILE_TAIL - Warning!' write ( *, '(a)' ) ' A page was open. It is being forced closed.' state = 2 call ps_setting_int ( 'SET', 'STATE', state ) end if if ( state /= 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'EPS_FILE_TAIL - Fatal error!' write ( *, '(a,i9)' ) ' PostScript state is ', state write ( *, '(a)' ) ' PostScript state 2 is required.' return end if ! ! Get the unit number. ! call ps_setting_int ( 'GET', 'UNIT', unit ) ! ! Retrieve the number of pages. ! call ps_setting_int ( 'GET', 'NUM_PAGES', num_pages ) if ( 1 < num_pages ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'EPS_FILE_TAIL - Warning!' write ( *, '(a)' ) ' An encapsulated PostScript file describes ONE page.' write ( *, '(a,i9,a)' ) ' This file describes ', num_pages, ' pages.' write ( *, '(a)' ) ' It is not a legal EPS file.' end if ! ! Write the epilog. ! write ( unit, '(a)' ) '%%Trailer' ! write ( unit, '(a)' ) 'end' write ( unit, '(a)' ) '%%EOF' ! ! Zero out the number of pages. ! num_pages = 0 call ps_setting_int ( 'SET', 'NUM_PAGES', num_pages ) ! ! Reset the state. ! state = 4 call ps_setting_int ( 'SET', 'STATE', state ) return end subroutine file_column_count ( file_name, ncolumn ) !*****************************************************************************80 ! !! FILE_COLUMN_COUNT counts the number of columns in the first line of a file. ! ! Discussion: ! ! The file is assumed to be a simple text file. ! ! Most lines of the file is presumed to consist of 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 s_word_count ( line, ncolumn ) return end subroutine file_column_range ( file_name, ncolumn, col_min, col_max ) !*****************************************************************************80 ! !! 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 ( kind = 8 ) COL_MIN(NCOLUM), COL_MAX(NCOLUMN), the ! minimum and maximum for each column. ! implicit none integer ncolumn real ( kind = 8 ) col_max(ncolumn) real ( kind = 8 ) col_min(ncolumn) character ( len = * ) file_name integer ierror integer ios integer iunit integer j character ( len = 256 ) line integer nrow real ( kind = 8 ) 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.0D+00 col_max(1:ncolumn) = 0.0D+00 do read ( iunit, '(a)', iostat = ios ) line if ( ios /= 0 ) then exit end if call s_to_r8vec ( 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 ) !*****************************************************************************80 ! !! 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 ) !*****************************************************************************80 ! !! 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_max < len_name + 1 ) then return end if len_name = len_name + 1 file_name(len_name:len_name) = '.' i = len_name + 1 else i = i + 1 file_name(i:j) = ' ' end if file_name(i:) = ext return end subroutine file_name_inc ( file_name ) !*****************************************************************************80 ! !! 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 ) !*****************************************************************************80 ! !! 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. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is 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. ! ! Modified: ! ! 18 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IUNIT, the free unit number. ! implicit none integer i integer ios integer iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end function i4_log_2 ( i ) !*****************************************************************************80 ! !! I4_LOG_2 returns the integer part of the logarithm base 2 of |I|. ! ! Example: ! ! I I4_LOG_2 ! ! 0 -1 ! 1, 0 ! 2, 1 ! 3, 1 ! 4, 2 ! 5, 2 ! 6, 2 ! 7, 2 ! 8, 3 ! 9, 3 ! 10, 3 ! 127, 6 ! 128, 7 ! 129, 7 ! ! Modified: ! ! 13 January 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the number whose logarithm base 2 is desired. ! ! Output, integer I4_LOG_2, the integer part of the logarithm base 2 of ! the absolute value of I. ! For positive I4_LOG_2(I), it should be true that ! 2**I4_LOG_2(X) <= |I| < 2**(I4_LOG_2(I)+1). ! The special case of I4_LOG_2(0) returns -HUGE(). ! implicit none integer i4_log_2 integer i integer i_abs if ( i == 0 ) then i4_log_2 = - huge ( i4_log_2 ) else i4_log_2 = 0 i_abs = abs ( i ) do while ( 2 <= i_abs ) i_abs = i_abs / 2 i4_log_2 = i4_log_2 + 1 end do end if return end function i4_modp ( i, j ) !*****************************************************************************80 ! !! I4_MODP returns the nonnegative remainder of integer division. ! ! Formula: ! ! If ! NREM = I4_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, I4_MODP(A,360) is between 0 and 360, always. ! ! Examples: ! ! I J MOD I4_MODP Factorization ! ! 107 50 7 7 107 = 2 * 50 + 7 ! 107 -50 7 7 107 = -2 * -50 + 7 ! -107 50 -7 43 -107 = -3 * 50 + 43 ! -107 -50 -7 43 -107 = 3 * -50 + 43 ! ! 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 I4_MODP, the nonnegative remainder when I is ! divided by J. ! implicit none integer i integer i4_modp integer j if ( j == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4_MODP - Fatal error!' write ( *, '(a,i6)' ) ' I4_MODP ( I, J ) called with J = ', j stop end if i4_modp = mod ( i, j ) if ( i4_modp < 0 ) then i4_modp = i4_modp + abs ( j ) end if return end subroutine i4_swap ( i, j ) !*****************************************************************************80 ! !! I4_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 subroutine i4_to_angle ( i, angle ) !*****************************************************************************80 ! !! I4_TO_ANGLE maps integers to points on a circle. ! ! Discussion: ! ! The angles are intended to be used to select colors on a color ! hexagon whose 6 vertices are red, yellow, green, cyan, blue, ! magenta. ! ! Example: ! ! I X ANGLE ! ! 0 0/3 0 ! 1 1/3 120 ! 2 2/3 240 ! ! 3 1/6 60 ! 4 3/6 180 ! 5 5/6 300 ! ! 6 1/12 30 ! 7 3/12 90 ! 8 5/12 150 ! 9 7/12 210 ! 10 9/12 270 ! 11 11/12 330 ! ! 12 1/24 15 ! 13 3/24 45 ! 14 5/24 75 ! etc ! ! Modified: ! ! 14 January 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the index of the desired color. ! ! Output, real ( kind = 8 ) ANGLE, an angle, measured in degrees, ! between 0 and 360. ! implicit none real ( kind = 8 ) angle integer i integer i4_log_2 integer i1 integer i2 integer i3 integer i4 if ( 0 <= abs ( i ) .and. abs ( i ) <= 2 ) then angle = 120.0D+00 * real ( abs ( i ), kind = 8 ) else i1 = i4_log_2 ( abs ( i ) / 3 ) i2 = abs ( i ) + 1 - 3 * 2**i1 i3 = 2 * ( i2 - 1 ) + 1 i4 = 3 * 2**( i1 + 1 ) angle = 360.0D+00 * real ( i3, kind = 8 ) / real ( i4, kind = 8 ) end if return end subroutine i4_to_rgb ( i, r, g, b ) !*****************************************************************************80 ! !! I4_TO_RGB maps integers to RGB colors. ! ! Modified: ! ! 21 January 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the index of the desired color. ! ! Output, integer R, G, B, the RGB specifications for a color. ! implicit none real ( kind = 8 ) angle integer b integer g integer i integer r real ( kind = 8 ) rb real ( kind = 8 ) rg real ( kind = 8 ) rr ! ! Red ! if ( i == 1 ) then r = 255 g = 0 b = 0 ! ! Green ! else if ( i == 2 ) then r = 0 g = 255 b = 0 ! ! Blue ! else if ( i == 3 ) then r = 0 g = 0 b = 255 ! ! Cyan ! else if ( i == 4 ) then r = 0 g = 255 b = 255 ! ! Magenta ! else if ( i == 5 ) then r = 255 g = 0 b = 255 ! ! Yellow ! else if ( i == 6 ) then r = 255 g = 255 b = 0 ! ! Brown5 ! else if ( i == 7 ) then r = 139 g = 35 b = 35 ! ! Orange ! else if ( i == 8 ) then r = 255 g = 165 b = 0 ! ! Goldenrod5 ! else if ( i == 9 ) then r = 139 g = 105 b = 20 ! ! Medium Purple ! else if ( i == 10 ) then r = 147 g = 112 b = 219 ! ! Coral ! else if ( i == 11 ) then r = 255 g = 127 b = 80 ! ! Pink5 ! else if ( i == 12 ) then r = 139 g = 99 b = 108 ! ! GreenYellow ! else if ( i == 13 ) then r = 173 g = 255 b = 47 ! ! Aquamarine ! else if ( i == 14 ) then r = 127 g = 255 b = 212 ! ! Pale Green3 ! else if ( i == 15 ) then r = 124 g = 205 b = 124 ! ! Burlywood ! else if ( i == 16 ) then r = 222 g = 184 b = 135 ! ! Cornsilk3 ! else if ( i == 17 ) then r = 205 g = 200 b = 177 ! ! Lemon_Chiffon3 ! else if ( i == 18 ) then r = 205 g = 201 b = 165 ! ! Maroon ! else if ( i == 19 ) then r = 176 g = 48 b = 96 ! ! Slate_Blue2 ! else if ( i == 20 ) then r = 131 g = 111 b = 255 else call i4_to_angle ( i, angle ) call angle_to_rgb ( angle, rr, rg, rb ) r = min ( int ( rr * 255 ), 255 ) g = min ( int ( rg * 255 ), 255 ) b = min ( int ( rb * 255 ), 255 ) end if return end function i4_wrap ( ival, ilo, ihi ) !*****************************************************************************80 ! !! I4_WRAP forces an integer to lie between given limits by wrapping. ! ! Example: ! ! ILO = 4, IHI = 8 ! ! I I4_WRAP ! ! -2 8 ! -1 4 ! 0 5 ! 1 6 ! 2 7 ! 3 8 ! 4 4 ! 5 5 ! 6 6 ! 7 7 ! 8 8 ! 9 4 ! 10 5 ! 11 6 ! 12 7 ! 13 8 ! 14 4 ! ! 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 I4_WRAP, a "wrapped" version of IVAL. ! implicit none integer i4_modp integer i4_wrap integer ihi integer ilo integer ival integer wide wide = ihi + 1 - ilo if ( wide == 0 ) then i4_wrap = ilo else i4_wrap = ilo + i4_modp ( ival-ilo, wide ) end if return end subroutine i4mat_transpose_print ( m, n, a, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT prints an I4MAT, transposed. ! ! Modified: ! ! 28 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), an M by N matrix to be printed. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer m integer n integer a(m,n) character ( len = * ) title call i4mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT_SOME prints some of the transpose of an I4MAT. ! ! Modified: ! ! 09 February 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: incx = 10 integer m integer n integer a(m,n) character ( len = 7 ) ctemp(incx) integer i integer i2 integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7)') i end do write ( *, '('' Row '',10a7)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' write ( *, '(a)' ) ' ' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(i7)' ) a(i,j) end do write ( *, '(i5,1x,10a7)' ) j, ( ctemp(i), i = 1, inc ) end do end do write ( *, '(a)' ) ' ' return end subroutine i4vec_heap_d ( n, a ) !*****************************************************************************80 ! !! I4VEC_HEAP_D reorders an array of integers into a descending heap. ! ! Discussion: ! ! A descending heap is an array A with the property that, for every index J, ! A(2*J) <= A(J) and A(2*J+1) <= A(J), (as long as the indices ! 2*J and 2*J+1 are legal). ! ! Diagram: ! ! A(1) ! / \ ! A(2) A(3) ! / \ / \ ! A(4) A(5) A(6) A(7) ! / \ / \ ! A(8) A(9) A(10) A(11) ! ! Modified: ! ! 15 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! A Nijenhuis and H Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Parameters: ! ! Input, integer N, the size of the input array. ! ! Input/output, integer A(N). ! On input, an unsorted array. ! On output, the array has been reordered into a heap. ! implicit none integer n integer a(n) integer i integer ifree integer key integer m ! ! Only nodes N/2 down to 1 can be "parent" nodes. ! do i = n/2, 1, -1 ! ! Copy the value out of the parent node. ! Position IFREE is now "open". ! key = a(i) ifree = i do ! ! Positions 2*IFREE and 2*IFREE + 1 are the descendants of position ! IFREE. (One or both may not exist because they exceed N.) ! m = 2 * ifree ! ! Does the first position exist? ! if ( n < m ) then exit end if ! ! Does the second position exist? ! if ( m + 1 <= n ) then ! ! If both positions exist, take the larger of the two values, ! and update M if necessary. ! if ( a(m) < a(m+1) ) then m = m + 1 end if end if ! ! If the large descendant is larger than KEY, move it up, ! and update IFREE, the location of the free position, and ! consider the descendants of THIS position. ! if ( a(m) <= key ) then exit end if a(ifree) = a(m) ifree = m end do ! ! Once there is no more shifting to do, KEY moves into the free spot IFREE. ! a(ifree) = key end do return end subroutine i4vec_indicator ( n, a ) !*****************************************************************************80 ! !! I4VEC_INDICATOR sets an integer vector to the indicator vector. ! ! 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 i4vec_sort_heap_a ( n, a ) !*****************************************************************************80 ! !! I4VEC_SORT_HEAP_A ascending sorts an integer array using heap sort. ! ! Modified: ! ! 15 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! A Nijenhuis and H Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Parameters: ! ! Input, integer N, the number of entries in the array. ! ! Input/output, integer A(N). ! On input, the array to be sorted; ! On output, the array has been sorted. ! implicit none integer n integer a(n) integer n1 if ( n <= 1 ) then return end if ! ! 1: Put A into descending heap form. ! call i4vec_heap_d ( n, a ) ! ! 2: Sort A. ! ! The largest object in the heap is in A(1). ! Move it to position A(N). ! call i4_swap ( a(1), a(n) ) ! ! Consider the diminished heap of size N1. ! do n1 = n-1, 2, -1 ! ! Restore the heap structure of A(1) through A(N1). ! call i4vec_heap_d ( n1, a ) ! ! Take the largest object from A(1) and move it to A(N1). ! call i4_swap ( a(1), a(n1) ) end do return end subroutine i4vec_sorted_unique ( n, a, unique_num ) !*****************************************************************************80 ! !! I4VEC_SORTED_UNIQUE gets the unique elements in a sorted integer array. ! ! Modified: ! ! 09 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of elements in A. ! ! Input/output, integer A(N). On input, the sorted ! integer array. On output, the unique elements in A. ! ! Output, integer UNIQUE_NUM, the number of unique elements in A. ! implicit none integer n integer a(n) integer itest integer unique_num unique_num = 0 if ( n <= 0 ) then return end if unique_num = 1 do itest = 2, n if ( a(itest) /= a(unique_num) ) then unique_num = unique_num + 1 a(unique_num) = a(itest) end if end do return end subroutine kmeans_plot ( box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, file_name, dim_num, & node_num, coord, plot_min, plot_max, x, y, connect, tag, shadow, title ) !*****************************************************************************80 ! !! 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. The third coordinate of ! such points is presumed to be an integer that indicates which ! cluster they belong to. ! ! A Voronoi diagram is made of the centers. ! ! Then the connected points are displayed, with colors corresponding ! to their assigned center. ! ! Modified: ! ! 22 July 2003 ! ! 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 NODE_NUM, the number of points. ! ! Input, real ( kind = 8 ) COORD(DIM_NUM,NODE_NUM), the point coordinates. ! ! Input, real ( kind = 8 ) PLOT_MIN(2), PLOT_MAX(2), the minimum and maximum ! X and Y values to plot. ! ! Input, integer X, Y, the indices of the two dimensions to plot. ! ! Input, integer CONNECT(NODE_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(NODE_NUM), a tag for each point, ! presumably the cluster index. ! implicit none integer dim_num integer node_num real ( kind = 8 ) b integer b_int real ( kind = 8 ) bitx real ( kind = 8 ) bity logical box_requested real ( kind = 8 ) box_xy(2,2) integer center_num logical circle_requested real ( kind = 8 ) circle_r real ( kind = 8 ) circle_x real ( kind = 8 ) circle_y integer connect(node_num) real ( kind = 8 ) coord(dim_num,node_num) real ( kind = 8 ) coord2(2,node_num) logical, parameter :: debug = .false. character ( len = * ) file_name integer file_unit real ( kind = 8 ) g integer g_int integer i integer i4_wrap integer i1 integer i2 integer i3 integer ierror integer ind(node_num) logical inside integer ix integer j integer jp1 integer jp2 integer k real ( kind = 8 ) n1 real ( kind = 8 ) n2 real ( kind = 8 ) plot_max(2) real ( kind = 8 ) plot_min(2) real ( kind = 8 ) r integer r_int logical shadow real ( kind = 8 ) t(2,3) integer tag(node_num) integer tag2(node_num) integer tag_lo integer tag_hi character ( len = * ) title real ( kind = 8 ) triangle_center(dim_num,2*node_num) integer triangle_neighbor(3,2*node_num) integer triangle_node(3,2*node_num) integer triangle_num integer x real ( kind = 8 ) xvec(4) real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) x4 real ( kind = 8 ) x5 real ( kind = 8 ) x6 integer y real ( kind = 8 ) yvec(4) real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) y4 real ( kind = 8 ) y5 real ( kind = 8 ) y6 ! ! Determine which points represent the centers. ! center_num = 0 if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CENTERS:' write ( *, '(a)' ) ' ' end if do i = 1, node_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 if ( center_num < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'KMEANS_PLOT - Fatal error!' write ( *, '(a)' ) ' There are no centers.' return end if if ( center_num == node_num ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'KMEANS_PLOT - Fatal error!' write ( *, '(a)' ) ' There are only centers, and no data points.' return end if ! ! Compute the Delaunay triangulation of the centers. ! call dtris2 ( center_num, coord2, triangle_num, triangle_node, & triangle_neighbor ) ! ! Compute the intersection point of the perpendicular bisectors ! of each Delaunay triangle. ! do j = 1, triangle_num i1 = triangle_node(1,j) i2 = triangle_node(2,j) i3 = triangle_node(3,j) t(1:2,1) = coord2(1:2,i1) t(1:2,2) = coord2(1:2,i2) t(1:2,3) = coord2(1:2,i3) call triangle_circumcenter_2d ( t, triangle_center(1:2,j) ) end do if ( debug ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'KMEANS_PLOT:' write ( *, '(a)' ) ' The data has been calculated.' 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, 36, 36, 576, 756 ) call ps_page_head ( plot_min(1), plot_min(2), plot_max(1), plot_max(2) ) ! ! Initialize the line width to 1. ! call ps_line_width ( 1 ) ! ! Print title, if requested. ! if ( 0 < len_trim ( title ) ) then call ps_font_size ( 0.30D+00 ) call ps_moveto ( plot_min(1), plot_max(2) ) call ps_label ( title ) end if ! ! Define a PostScript clipping box. ! xvec(1:4) = (/ plot_min(1), plot_max(1), plot_max(1), plot_min(1) /) yvec(1:4) = (/ plot_min(2), plot_min(2), plot_max(2), plot_max(2) /) call ps_clip ( 4, xvec, yvec ) ! ! Shadow the points, if requested. ! if ( shadow ) then bitx = 0.025D+00 * ( plot_max ( 1 ) - plot_min ( 1 ) ) bity = 0.025D+00 * ( plot_max ( 2 ) - plot_min ( 2 ) ) do i = 1, node_num call ps_line ( plot_min(1), coord(y,i), plot_min(1) + bitx, coord(y,i) ) call ps_line ( coord(x,i), plot_min(2), coord(x,i), plot_min(2) + bity ) end do end if ! ! Draw a user-specified box, if requested. ! if ( box_requested ) then call ps_comment ( 'User-requested box:' ) call ps_line ( box_xy(1,1), box_xy(2,1), box_xy(1,2), box_xy(2,1) ) call ps_line ( box_xy(1,2), box_xy(2,1), box_xy(1,2), box_xy(2,2) ) call ps_line ( box_xy(1,2), box_xy(2,2), box_xy(1,1), box_xy(2,2) ) call ps_line ( box_xy(1,1), box_xy(2,2), box_xy(1,1), box_xy(2,1) ) end if ! ! Draw a user-specified circle, if requested. ! if ( circle_requested ) then call ps_comment ( 'User-requested circle:' ) call ps_circle ( circle_x, circle_y, circle_r ) end if x1 = plot_min(1) y1 = plot_min(2) call ps_moveto ( x1, y1 ) ! call ps_label ( 'K-Means diagram' ) call ps_setting_int ( 'SET', 'MARKER_SIZE', 8 ) tag_lo = 1 tag_hi = center_num do i = 1, center_num call i4_to_rgb ( tag2(i), r_int, g_int, b_int ) r = real ( r_int, kind = 8 ) / 255.0D+00 g = real ( g_int, kind = 8 ) / 255.0D+00 b = real ( b_int, kind = 8 ) / 255.0D+00 call ps_color_fill_set ( r, g, b ) call ps_mark_disk ( coord2(1,i), coord2(2,i) ) end do ! ! Mark the cluster generators. ! call ps_setting_int ( 'SET', 'MARKER_SIZE', 4 ) do i = 1, node_num if ( connect(i) /= 0 ) then call i4_to_rgb ( tag(i), r_int, g_int, b_int ) r = real ( r_int, kind = 8 ) / 255.0D+00 g = real ( g_int, kind = 8 ) / 255.0D+00 b = real ( b_int, kind = 8 ) / 255.0D+00 call ps_color_line_set ( r, g, b ) call ps_mark_circle ( coord(x,i), coord(y,i) ) end if end do ! ! Construct the Voronoi region boundaries. ! r = 0.0D+00 g = 0.0D+00 b = 0.0D+00 call ps_color_fill_set ( r, g, b ) ! ! For each Delaunay triangle, I ! For each side J, ! do i = 1, triangle_num do j = 1, 3 k = triangle_neighbor(j,i) ! ! If there is a neighboring triangle K on that side, ! connect the circumcenters. ! if ( 0 < k ) then if ( i < k ) then call ps_line ( triangle_center(1,i), triangle_center(2,i), & triangle_center(1,k), triangle_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 = triangle_node(j,i) x1 = coord2(1,ix) y1 = coord2(2,ix) jp1 = i4_wrap ( j+1, 1, 3 ) ix = triangle_node(jp1,i) x2 = coord2(1,ix) y2 = coord2(2,ix) jp2 = i4_wrap ( j+2, 1, 3 ) ix = triangle_node(jp2,i) x3 = coord2(1,ix) y3 = coord2(2,ix) x4 = triangle_center(1,i) y4 = triangle_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 ( plot_min(1), plot_min(2), plot_max(1), & plot_max(2), 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 ) !*****************************************************************************80 ! !! 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 ( kind = 8 ) 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 ( kind = 8 ) A, B, C, three coefficients which describe ! the line that passes through (X1,Y1) and (X2,Y2). ! implicit none real ( kind = 8 ) a real ( kind = 8 ) b real ( kind = 8 ) c real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) y1 real ( kind = 8 ) 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 ) !*****************************************************************************80 ! !! 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 ( kind = 8 ) X1, Y1, X2, Y2, two points on the line. ! ! Output, real ( kind = 8 ) 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 ( kind = 8 ) n1 real ( kind = 8 ) n2 real ( kind = 8 ) norm real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) y1 real ( kind = 8 ) y2 norm = sqrt ( ( x2 - x1 )**2 + ( y2 - y1 )**2 ) if ( norm == 0.0D+00 ) then n1 = 0.0D+00 n2 = 0.0D+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 ) !*****************************************************************************80 ! !! 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 ( kind = 8 ) X1, Y1, X2, Y2, define the first line. ! ! Input, real ( kind = 8 ) 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 ( kind = 8 ) X, Y, if IVAl = 1, then X, Y contains ! the intersection point. Otherwise, X = 0, Y = 0. ! implicit none real ( kind = 8 ) a1 real ( kind = 8 ) a2 real ( kind = 8 ) b1 real ( kind = 8 ) b2 real ( kind = 8 ) c1 real ( kind = 8 ) c2 integer ival logical point_1 logical point_2 real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) x4 real ( kind = 8 ) y real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) y4 ival = 0 x = 0.0D+00 y = 0.0D+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 ) !*****************************************************************************80 ! !! 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 ( kind = 8 ) A1, B1, C1, define the first line. ! At least one of A1 and B1 must be nonzero. ! ! Input, real ( kind = 8 ) 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 ( kind = 8 ) X, Y, if IVAL = 1, then X, Y contains ! the intersection point. Otherwise, X = 0, Y = 0. ! implicit none real ( kind = 8 ) a(2,2) real ( kind = 8 ) a1 real ( kind = 8 ) a2 real ( kind = 8 ) b(2,2) real ( kind = 8 ) b1 real ( kind = 8 ) b2 real ( kind = 8 ) c1 real ( kind = 8 ) c2 real ( kind = 8 ) det integer ival real ( kind = 8 ) x real ( kind = 8 ) y x = 0.0D+00 y = 0.0D+00 ! ! Refuse to handle degenerate lines. ! if ( a1 == 0.0D+00 .and. b1 == 0.0D+00 ) then ival = -1 return else if ( a2 == 0.0D+00 .and. b2 == 0.0D+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 r8mat2_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.0D+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.0D+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 ) !*****************************************************************************80 ! !! LRLINE determines if a point is left of, right or, or on a directed line. ! ! Discussion: ! ! The directed line is parallel to, and at a signed distance DV from ! a directed base line from (XV1,YV1) to (XV2,YV2). ! ! Modified: ! ! 14 July 2001 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) XU, YU, the coordinates of the point whose ! position relative to the directed line is to be determined. ! ! Input, real ( kind = 8 ) XV1, YV1, XV2, YV2, the coordinates of two points ! that determine the directed base line. ! ! Input, real ( kind = 8 ) DV, the signed distance of the directed line ! from the directed base line through the points (XV1,YV1) and (XV2,YV2). ! DV is positive for a line to the left of the base line. ! ! Output, integer LRLINE, the result: ! +1, the point is to the right of the directed line; ! 0, the point is on the directed line; ! -1, the point is to the left of the directed line. ! implicit none real ( kind = 8 ) dv real ( kind = 8 ) dx real ( kind = 8 ) dxu real ( kind = 8 ) dy real ( kind = 8 ) dyu integer lrline real ( kind = 8 ) t real ( kind = 8 ) tol real ( kind = 8 ) tolabs real ( kind = 8 ) xu real ( kind = 8 ) xv1 real ( kind = 8 ) xv2 real ( kind = 8 ) yu real ( kind = 8 ) yv1 real ( kind = 8 ) yv2 tol = 100.0D+00 * epsilon ( tol ) dx = xv2 - xv1 dy = yv2 - yv1 dxu = xu - xv1 dyu = yu - yv1 tolabs = tol * max ( abs ( dx ), abs ( dy ), abs ( dxu ), & abs ( dyu ), abs ( dv ) ) t = dy * dxu - dx * dyu + dv * sqrt ( dx * dx + dy * dy ) if ( tolabs < t ) then lrline = 1 else if ( -tolabs <= t ) then lrline = 0 else lrline = -1 end if return end subroutine perm_inv ( n, p ) !*****************************************************************************80 ! !! PERM_INV inverts a permutation "in place". ! ! Modified: ! ! 25 July 2000 ! ! Parameters: ! ! Input, integer N, the number of objects being permuted. ! ! Input/output, integer P(N), the permutation, in standard index form. ! On output, P describes the inverse permutation ! implicit none integer n integer i integer i0 integer i1 integer i2 integer ierror integer is integer p(n) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PERM_INV - Fatal error!' write ( *, '(a,i6)' ) ' Input value of N = ', n stop end if is = 1 do i = 1, n i1 = p(i) do while ( i < i1 ) i2 = p(i1) p(i1) = -i2 i1 = i2 end do is = -sign ( 1, p(i) ) p(i) = sign ( p(i), is ) end do do i = 1, n i1 = -p(i) if ( 0 <= i1 ) then i0 = i do i2 = p(i1) p(i1) = i0 if ( i2 < 0 ) then exit end if i0 = i1 i1 = i2 end do end if end do return end subroutine plot_size ( dim_num, box_requested, box_xy, circle_requested, & circle_x, circle_y, circle_r, x, y, plot_min, plot_max ) !*****************************************************************************80 ! !! PLOT_SIZE determines a plot range. ! ! Modified: ! ! 28 March 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer DIM_NUM, the number of spatial dimensions. ! ! Input, logical BOX_REQUESTED, is true if the user has specified ! a box to be drawn. ! ! Input, real ( kind = 8 ) BOX_XY(2,2), the coordinates of the lower left and ! upper right corners of the requested box. ! ! Input, integer X, Y, the dimensions to be used for X and Y. ! ! Output, real ( kind = 8 ) plot_min(2), plot_max(2), the ! minimum and maximum coordinates to use in the plot. ! use AA_point_data_module implicit none logical box_requested real ( kind = 8 ) box_xy(2,2) logical circle_requested real ( kind = 8 ) circle_r real ( kind = 8 ) circle_x real ( kind = 8 ) circle_y integer dim_num real ( kind = 8 ), parameter :: grace = 0.05D+00 integer i real ( kind = 8 ) plot_max(2) real ( kind = 8 ) plot_min(2) real ( kind = 8 ) width integer x integer y integer z do i = 1, 2 if ( i == 1 ) then z = x else if ( i == 2 ) then z = y end if plot_max(i) = coord_max(z) plot_min(i) = coord_min(z) if ( box_requested ) then plot_max(i) = max ( plot_max(i), box_xy(i,2) ) plot_min(i) = min ( plot_min(i), box_xy(i,1) ) end if if ( circle_requested ) then if ( i == 1 ) then plot_max(i) = max ( plot_max(i), circle_x + circle_r ) plot_min(i) = min ( plot_min(i), circle_x - circle_r ) else if ( i == 2 ) then plot_max(i) = max ( plot_max(i), circle_y + circle_r ) plot_min(i) = min ( plot_min(i), circle_y - circle_r ) end if end if width = plot_max(i) - plot_min(i) plot_max(i) = plot_max(i) + grace * width plot_min(i) = plot_min(i) - grace * width end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Plot Plot' write ( *, '(a)' ) ' Coord Min Max' write ( *, '(a)' ) ' ' do i = 1, 2 write ( *, '(i3,2f10.4)' ) i, plot_min(i), plot_max(i) end do return end function point_inside_box_2d ( x1, y1, x2, y2, x, y ) !*****************************************************************************80 ! !! POINT_INSIDE_BOX_2D determines if a point is inside a box in 2D. ! ! Definition: ! ! A "box" is defined by its "left down" corner and its ! "right up" corner, and all the points between. It is ! assumed that the sides of the box align with coordinate directions. ! ! Modified: ! ! 01 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) X1, Y1, X2, Y2, the two corners of the box. ! ! Input, real ( kind = 8 ) X, Y, the point to be checked. ! ! Output, logical POINT_INSIDE_BOX_2D, is .TRUE. if (X,Y) is inside the ! box, or on its boundary, and .FALSE. otherwise. ! implicit none logical point_inside_box_2d real ( kind = 8 ) x real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) y real ( kind = 8 ) y1 real ( kind = 8 ) y2 if ( x1 <= x .and. x <= x2 .and. & y1 <= y .and. y <= y2 ) then point_inside_box_2d = .true. else point_inside_box_2d = .false. end if return end function points_avoid_point_nd ( dim_num, n, xy_set, xy_test, tol ) !*****************************************************************************80 ! !! 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 DIM_NUM, the spatial dimension. ! ! Input, integer N, the number of points in the set. ! ! Input, real ( kind = 8 ) XY_SET(DIM_NUM,N), the set of points ! to be avoided. ! ! Input, real ( kind = 8 ) XY_TEST(DIM_NUM), a point to be tested. ! ! Input, real ( kind = 8 ) 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 dim_num integer j logical points_avoid_point_nd real ( kind = 8 ) tol real ( kind = 8 ) xy_set(dim_num,n) real ( kind = 8 ) xy_test(dim_num) points_avoid_point_nd = .true. do j = 1, n if ( sum ( ( xy_set(1:dim_num,j) - xy_test(1:dim_num) )**2 ) < tol**2 ) then points_avoid_point_nd = .false. return end if end do return end subroutine points_count ( file_in_name, dim_num, node_num ) !*****************************************************************************80 ! !! 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 NODE_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 node_num integer record_num real ( kind = 8 ) 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 node_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_r8vec ( line, dim_num, x, ierror ) if ( ierror == 0 ) then node_num = node_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: ', node_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, node_num, coord ) !*****************************************************************************80 ! !! POINTS_READ reads point coordinates from a file. ! ! Modified: ! ! 06 March 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer CONNECT(NODE_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 NODE_NUM, the number of points. The program ! will stop reading data once NODE_NUM values have been read. ! ! Output, real ( kind = 8 ) COORD(DIM_NUM,NODE_NUM), the point coordinates. ! implicit none integer dim_num integer node_num integer connect(node_num) real ( kind = 8 ) coord(dim_num,node_num) character ( len = * ) file_in_name integer file_in_unit integer i integer ierror integer ios character ( len = 100 ) line integer p real ( kind = 8 ) 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:node_num) = 0 i = 0 p = 0 do while ( i < node_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 e