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

  character c1
  character c1_cap
  character c2
  character c2_cap
  logical ch_eqi

  c1_cap = c1
  c2_cap = c2

  call ch_cap ( c1_cap )
  call ch_cap ( c2_cap )

  if ( c1_cap == c2_cap ) then
    ch_eqi = .true.
  else
    ch_eqi = .false.
  end if

  return
end
subroutine ch_to_digit ( c, digit )

!*****************************************************************************80
!
!! CH_TO_DIGIT returns the 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 dist_print ( n, dist, title )

!*****************************************************************************80
!
!! DIST_PRINT prints the distance matrix.
!
!  Modified:
!
!    13 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of data items.
!
!    Input, real ( kind = 8 ) DIST(N,N), the distance matrix.
!
!    Input, character ( len = * ) TITLE, an optional title.
!
  implicit none

  integer n

  real ( kind = 8 ) dist(n,n)
  character ( len = * ) title

  if ( len_trim ( title ) == 0 ) then

    call r8mat_print ( n, n, dist, '  The distance matrix:' )

  else

    call r8mat_print ( n, n, dist, title )

  end if

  return
end
subroutine dist_read ( file_name, n, dist )

!*****************************************************************************80
!
!! DIST_READ reads a distance matrix from a file.
!
!  Discussion:
!
!    The data is stored in a file, with each row of the distance
!    matrix written on one or more lines.
!
!    Blank lines and comment lines (beginning with '#') are ignored.
!
!    Individual data values should be separated by spaces or commas.
!
!  Modified:
!
!    08 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file to read.
!
!    Input, integer N, the number of data items.
!
!    Output, real ( kind = 8 ) DIST(N,N), the distance matrix values.
!
  implicit none

  integer n

  real ( kind = 8 ) dist(n,n)
  character ( len = * ) file_name
  integer i
  integer ierror
  integer input
  integer ios
  integer j
  integer last
  integer length
  character ( len = 80 ) line
  integer line_num
  real ( kind = 8 ) value

  call get_unit ( input )

  open ( unit = input, file = file_name, status = 'old' )

  dist(1:n,1:n) = huge ( dist(1,1) )

  i = 1
  j = 0
  line_num = 0

  do
!
!  Have we read enough data?
!
    if ( i == n .and. j == n ) then
      exit
    end if
!
!  Have we read too much data?
!
    if ( n < i .or. n < j ) then
      exit
    end if
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if

    line_num = line_num + 1
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else
!
!  LAST points to the last character associated with the previous 
!  data value read from the line.
!
      last = 0

      do
!
!  Try to read another value from the line.
!
        call s_to_r8 ( line(last+1:), value, ierror, length )
!
!  If we could not read a new value, it's time to read a new line.
!
        if ( ierror /= 0 ) then
          exit
        end if
!
!  Update the pointer.
!
        last = last + length
!
!  If we read a new value, where do we put it?
!
        j = j + 1

        if ( n < j ) then
          j = 1
          i = i + 1
          if ( n < i ) then
            exit
          end if
        end if

        dist(i,j) = value
!
!  If you reached the end of the row, it's time to read a new line.
!
        if ( j == n ) then
          exit
        end if

      end do

    end if

  end do

  close ( unit = input )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DIST_READ:'
  write ( *, '(a,i6)' ) '  Number of lines read was ', line_num

  return
end
subroutine dist_write ( file_name, n, dist )

!*****************************************************************************80
!
!! DIST_WRITE writes a distance matrix to a file.
!
!  Modified:
!
!    23 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file to write.
!
!    Input, integer N, the number of data items.
!
!    Input, real ( kind = 8 ) DIST(N,N), the data values.
!
  implicit none

  integer n

  real ( kind = 8 ) dist(n,n)
  character ( len = * ) file_name
  character ( len = 20 ) form
  integer i
  integer npline
  integer output

  call get_unit ( output )

  open ( unit = output, file = file_name, status = 'replace' )

  write ( output, '(a)' ) '# ' // trim ( file_name )
  write ( output, '(a)' ) '#'
  write ( output, '(a)' ) '#  Created by DIST_WRITE.'
  write ( output, '(a)' ) '#'
  write ( output, '(a)' ) '#  N by N distance matrix.'
  write ( output, '(a,i6)' ) '#  Number of points N is ', n
  write ( output, '(a)' ) '#'

  if ( all ( dist(1:n,1:n) == aint ( dist(1:n,1:n) ) ) ) then

    npline = min ( n, 10 )

    write ( form, '(a,i6,a)' ) '(', npline, 'i6)'

    do i = 1, n

      write ( output, form ) int ( dist(i,1:n) )

    end do

  else

    npline = min ( n, 10 )

    write ( form, '(a,i6,a)' ) '(', npline, 'g14.6)'

    do i = 1, n

      write ( output, form ) dist(i,1:n)

    end do

  end if

  close ( unit = output )

  return
end
subroutine dms_print ( n, lat_dms, long_dms, title )

!*****************************************************************************80
!
!! DMS_PRINT prints the latitude and longitude in degrees/minutes/seconds.
!
!  Modified:
!
!    23 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of data items.
!
!    Input, integer LAT_DMS(4,N), LONG_DMS(4,N), the latitudes
!    and longitudes, in degrees, minutes and seconds.  The fourth
!    argument is +1/-1 for North/South latitude or East/West longitude.
!
!    Input, character ( len = * ) TITLE, an optional title.
!
  implicit none

  integer n

  integer i
  character lat_char
  integer lat_dms(4,n)
  character long_char
  integer long_dms(4,n)
  character ( len = * ) title

  if ( len_trim ( title ) /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  #    Latitude            Longitude'
  write ( *, '(a)' ) '       (Deg/Min/Sec)       (Deg/Min/Sec)'
  write ( *, '(a)' ) '---  ---------------      ---------------'
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(i3,2x,3i4,2x,a1,5x,3i4,2x,a1)' ) i, &
      lat_dms(1:3,i), lat_char ( lat_dms(4,i) ), &
      long_dms(1:3,i), long_char ( long_dms(4,i) )
  end do

  return
end
subroutine dms_read ( file_name, n, lat_dms, long_dms )

!*****************************************************************************80
!
!! DMS_READ reads DMS data from a file.
!
!  Discussion:
!
!    Blank lines and comment lines (beginning with '#') are ignored.
!
!    Individual data values should be separated by spaces or commas.
!
!  Modified:
!
!    15 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file to read.
!
!    Input, integer N, the number of data items.
!
!    Output, integer LAT_DMS(4,N), LONG_DMS(4,N), the latitude and
!    longitudes, in degrees, minutes and seconds.  The fourth
!    argument is +1/-1 for North/South latitude or East/West longitude. 
!
  implicit none

  integer n
  integer, parameter :: ncol = 8

  character ( len = * ) file_name
  integer i
  integer ierror
  integer input
  integer ios
  integer j
  integer last
  integer lat_dms(4,n)
  integer length
  character ( len = 80 ) line
  integer line_num
  integer long_dms(4,n)
  logical s_eqi
  integer value
  integer vector(ncol)
  character ( len = 80 ) word

  call get_unit ( input )

  open ( unit = input, file = file_name, status = 'old' )

  lat_dms(1:4,1:n) = huge ( lat_dms(1,1) )
  long_dms(1:4,1:n) = huge ( long_dms(1,1) )

  i = 1
  j = 0
  line_num = 0

  do
!
!  Have we read enough data?
!
    if ( i == n .and. j == ncol ) then
      exit
    end if
!
!  Have we read too much data?
!
    if ( n < i .or. ncol < j ) then
      exit
    end if
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if

    line_num = line_num + 1
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else
!
!  LAST points to the last character associated with the previous 
!  data value read from the line.
!
      last = 0

      do
!
!  Try to read another value from the line.
!  Note that, confusingly, J right now indicates the column
!  of the PREVIOUS data item that was read.
!
        if ( j == 3 ) then

          call s_to_w ( line(last+1:), word, ierror, length )

          if ( ierror /= 0 ) then
            exit
          end if

          if ( s_eqi ( word(1:1), 'N' ) ) then
            value = +1
          else if ( s_eqi ( word(1:1), 'S' ) ) then
            value = -1
          else
            value = 0
          end if

        else if ( j == 7 ) then

          call s_to_w ( line(last+1:), word, ierror, length )

          if ( ierror /= 0 ) then
            exit
          end if

          if ( s_eqi ( word(1:1), 'E' ) ) then
            value = +1
          else if ( s_eqi ( word(1:1), 'W' ) ) then
            value = -1
          else
            value = 0
          end if

        else

          call s_to_i4 ( line(last+1:), value, ierror, length )

          if ( ierror /= 0 ) then
            exit
          end if

        end if
!
!  Update the pointer.
!
        last = last + length
!
!  If we read a new value, where do we put it?
!
        j = j + 1

        if ( ncol < j ) then
          j = 1
          i = i + 1
          if ( n < i ) then
            exit
          end if
        end if

        vector(j) = value
!
!  If you reached the end of the row, it's time to read a new line.
!
        if ( j == ncol ) then
          lat_dms(1:4,i) = vector(1:4)
          long_dms(1:4,i) = vector(5:8)
          exit
        end if

      end do

    end if

  end do

  close ( unit = input )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DMS_READ:'
  write ( *, '(a,i6)' ) '  Number of lines read was ', line_num

  return
end
subroutine dms_to_dist ( n, lat_dms, long_dms, dist )

!*****************************************************************************80
!
!! DMS_TO_DIST creates a distance matrix from latitudes and longitudes.
!
!  Discussion:
!
!    A distance function is used which is appropriate for the earth.
!
!  Modified:
!
!    15 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of data items.
!
!    Input, integer LAT_DMS(4,N), LONG_DMS(4,N), the latitude and
!    longitude, in degrees, minutes, and seconds, for each point.
!    The fourth argument is +1/-1 for North/South latitude or 
!    East/West longitude.
!
!    Output, real ( kind = 8 ) DIST(N,N), the distance matrix.  Distances are
!    measured in miles.
!
  implicit none

  integer n

  real ( kind = 8 ) dist(n,n)
  integer lat_dms(4,n)
  integer long_dms(4,n)
  integer i
  integer j

  do i = 1, n
    dist(i,i) = 0.0D+0
    do j = i+1, n
      call dms_to_distance_earth ( lat_dms(1,i), long_dms(1,i), &
        lat_dms(1,j), long_dms(1,j), dist(i,j) )
      dist(j,i) = dist(i,j)
    end do
  end do

  return
end
subroutine dms_to_distance_earth ( lat1_dms, long1_dms, lat2_dms, &
  long2_dms, dist )

!*****************************************************************************80
!
!! DMS_TO_DISTANCE_EARTH finds the distance between two points on the earth.
!
!  Discussion:
!
!    The positions of the the points are given as longitude and
!    latitude, measured in degrees, minutes, and seconds.
!
!    The distance is measured on the surface of the earth, which
!    is approximated by a sphere.
!
!  Modified:
!
!    15 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer LAT1_DMS(4), LONG1_DMS(4), the latitude and 
!    longitude of the first point.  The fourth
!    argument is +1/-1 for North/South latitude or East/West longitude.
!
!    Input, integer LAT2_DMS(4), LONG2_DMS(4), the latitude and 
!    longitude of the second point.  The fourth
!    argument is +1/-1 for North/South latitude or East/West longitude.
!
!    Output, real ( kind = 8 ) DIST, the distance between the points, in miles.
!
  implicit none

  real ( kind = 8 ) dist
  real ( kind = 8 ) dms_to_radians
  integer lat1_dms(4)
  real ( kind = 8 ) lat1_rad
  integer lat2_dms(4)
  real ( kind = 8 ) lat2_rad
  integer long1_dms(4)
  real ( kind = 8 ) long1_rad
  integer long2_dms(4)
  real ( kind = 8 ) long2_rad
  real ( kind = 8 ), parameter :: radius = 3958.89D+00
  real ( kind = 8 ) theta

  lat1_rad = dms_to_radians ( lat1_dms )
  long1_rad = dms_to_radians ( long1_dms )

  lat2_rad = dms_to_radians ( lat2_dms )
  long2_rad = dms_to_radians ( long2_dms )

  theta = acos ( sin ( lat1_rad ) * sin ( lat2_rad ) &
               + cos ( lat1_rad ) * cos ( lat2_rad ) &
               * cos ( long1_rad - long2_rad ) )

  dist = radius * theta

  return
end
function dms_to_radians ( dms )

!*****************************************************************************80
!
!! DMS_TO_RADIANS converts degrees, minutes, seconds to radians.
!
!  Modified:
!
!    15 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DMS(4), the measurement of an angle in
!    degrees, minutes and seconds.  The fourth
!    argument is +1/-1 for North/South latitude or East/West longitude.
!
!    Output, real ( kind = 8 ) DMS_TO_RADIANS, the measurement of the same
!    angle in radians.
!
  implicit none

  real ( kind = 8 ) d_real
  integer dms(4)
  real ( kind = 8 ) dms_to_radians
  integer i4_sign
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00

  d_real = real ( i4_sign ( dms(4) ), kind = 8 ) * &
    real ( dms(3) + 60 * dms(2) + 3600 * dms(1), kind = 8 ) / 3600.0D+00

  dms_to_radians = pi * d_real / 180.0D+00

  return
end
subroutine dms_to_xy ( n, lat_dms, long_dms, xy )

!*****************************************************************************80
!
!! DMS_TO_XY creates pseudo X,Y coordinates from latitudes and longitudes.
!
!  Discussion:
!
!    Essentially, the latitude and longitude information is treated
!    as though the Earth were a cylinder.  As long as the the
!    data is relatively close on the sphere (and far from either
!    pole!) the distortion will not be too severe.  If the data
!    is closely clustered, and also near the equator, the
!    positions will be relatively accurate.
!
!  Modified:
!
!    13 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of data items.
!
!    Input, integer LAT_DMS(4,N), LONG_DMS(4,N), the latitude and
!    longitude, in degrees, minutes, and seconds, for each point.
!    The fourth argument is +1/-1 for North/South latitude or 
!    East/West longitude.
!
!    Output, real ( kind = 8 ) XY(N,2), the XY coordinates, in miles.
!
  implicit none

  integer n

  real ( kind = 8 ) dms_to_radians
  integer lat_dms(4,n)
  integer long_dms(4,n)
  integer i
  real ( kind = 8 ) phi
  real ( kind = 8 ), parameter :: radius = 3958.89D+00
  real ( kind = 8 ) theta
  real ( kind = 8 ) xy(n,2)

  do i = 1, n
    theta = dms_to_radians ( long_dms(1,i) )
    phi = dms_to_radians ( lat_dms(1,i) )
    xy(i,1) = radius * theta
    xy(i,2) = radius * phi
  end do

  return
end
subroutine dms_write ( file_name, n, lat_dms, long_dms )

!*****************************************************************************80
!
!! DMS_WRITE writes a latitude, longitude file.
!
!  Modified:
!
!    15 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file to write.
!
!    Input, integer N, the number of data items.
!
!    Input, integer LAT_DMS(4,N), LONG_DMS(4,N), the data values.
!
  implicit none

  integer n

  character ( len = * ) file_name
  integer i
  character lat_char
  integer lat_dms(4,n)
  character long_char
  integer long_dms(4,n)
  integer output

  call get_unit ( output )

  open ( unit = output, file = file_name, status = 'replace' )

  write ( output, '(a)' ) '# ' // trim ( file_name )
  write ( output, '(a)' ) '#'
  write ( output, '(a)' ) '#  Created by DMS_WRITE.'
  write ( output, '(a)' ) '#'
  write ( output, '(a)' ) '#  Latitude, Longitude in degrees, minutes, seconds'
  write ( output, '(a,i6)' ) '#  Number of points N is ', n
  write ( output, '(a)' ) '#'

  do i = 1,  n

    write ( output, '(3i5,2x,3i5)' ) &
      lat_dms(1:3,i), lat_char ( lat_dms(4,i) ), &
      long_dms(1:3,i), long_char ( long_dms(4,i) )

  end do


  close ( unit = output )

  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_sign ( x )

!*****************************************************************************80
!
!! I4_SIGN evaluates the sign of an I4.
!
!  Discussion:
!
!    This function differs from the intrinsic SIGN function, because
!    it returns a value of 0 if the input argument is 0.
!
!  Modified:
!
!    09 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer X, the number whose sign is desired.
!
!    Output, integer I4_SIGN, a result based on the sign of X:
!
!    -1, if X < 0.
!     0, if X = 0.
!    +1, if X > 0.
!
  implicit none

  integer i4_sign
  integer x

  if ( x < 0 ) then
    i4_sign = -1
  else if ( x == 0 ) then
    i4_sign = 0
  else if ( 0 < x ) then
    i4_sign = +1
  end if

  return
end
function i4_to_a ( i )

!*****************************************************************************80
!
!! I4_TO_A returns the I-th alphabetic character.
!
!  Examples:
!
!    I  I4_TO_A
!
!   -8  ' '
!    0  ' '
!    1  'A'
!    2  'B'
!   ..
!   26  'Z'
!   27  'a'
!   52  'z'
!   53  ' '
!   99  ' '
!
!  Modified:
!
!    23 February 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, the index of the letter to be returned.
!    0 is a space;
!    1 through 26 requests 'A' through 'Z', (ASCII 65:90);
!    27 through 52 requests 'a' through 'z', (ASCII 97:122);
!
!    Output, character I4_TO_A, the requested alphabetic letter.
!
  implicit none

  integer, parameter :: cap_shift = 64
  integer i
  character i4_to_a
  integer, parameter :: low_shift = 96

  if ( i <= 0 ) then
    i4_to_a = ' '
  else if ( 1 <= i .and. i <= 26 ) then
    i4_to_a = char ( cap_shift + i )
  else if ( 27 <= i .and. i <= 52 ) then
    i4_to_a = char ( low_shift + i - 26 )
  else if ( 53 <= i ) then
    i4_to_a = ' '
  end if

  return
end
function lat_char ( i )

!*****************************************************************************80
!
!! LAT_CHAR returns a character for negative or positive latitude.
!
!  Modified:
!
!    15 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, is negative for negative latitude, and positive
!    for positive latitude.
!
!    Output, character LAT_CHAR, is 'S' for negative latitude, and
!    'N' for positive latitude.
!
  implicit none

  integer i
  character lat_char

  if ( i < 0 ) then
    lat_char = 'S'
  else if ( 0 < i ) then
    lat_char = 'N'
  else
    lat_char = '?'
  end if

  return
end
function long_char ( i )

!*****************************************************************************80
!
!! LONG_CHAR returns a character for negative or positive longitude.
!
!  Modified:
!
!    15 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer I, is negative for negative longitude, and positive
!    for positive longitude.
!
!    Output, character LONG_CHAR, is 'W' for negative longitude, and
!    'E' for positive longitude.
!
  implicit none

  integer i
  character long_char

  if ( i < 0 ) then
    long_char = 'W'
  else if ( 0 < i ) then
    long_char = 'E'
  else
    long_char = '?'
  end if

  return
end
subroutine main_read_code ( file_main, file_code )

!*****************************************************************************80
!
!! MAIN_READ_CODE reads the name of the code file from the main file.
!
!  Discussion:
!
!    FILE_CODE is the name of a file containing short codes for the 
!    cities.
!
!    There MAY be a record in the main file of the form
!
!    "code  key_code.txt"
!
!  Modified:
!
!    07 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_MAIN, the name of the file to read.
!
!    Output, character ( len = * ) FILE_CODE, the name of the code file,
!    or ' ' if no information was found.
!
  implicit none

  logical done
  character ( len = * ) file_code
  character ( len = * ) file_main
  integer input
  integer ios
  character ( len = 80 ) line
  logical s_eqi
  character ( len = 80 ) word

  call get_unit ( input )

  open ( unit = input, file = file_main, status = 'old' )

  file_code = ' '

  do
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else

      done = .true.

      call word_next_read ( line, word, done )

      if ( done ) then
        cycle
      end if

      if ( .not. s_eqi ( word, 'code' ) ) then
        cycle
      end if

      call word_next_read ( line, file_code, done )

      exit

    end if

  end do

  close ( unit = input )

  return
end
subroutine main_read_dist ( file_main, file_dist )

!*****************************************************************************80
!
!! MAIN_READ_DIST reads the name of the distance file from the main file.
!
!  Discussion:
!
!    FILE_DIST is the name of a file containing the city-to-city
!    distance matrix.
!
!    There MAY be a record in the main file of the form
!
!    "dist  key_dist.txt"
!
!  Modified:
!
!    07 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_MAIN, the name of the file to read.
!
!    Output, character ( len = * ) FILE_DIST, the name of the distance file,
!    or ' ' if no information was found.
!
  implicit none

  logical done
  character ( len = * ) file_dist
  character ( len = * ) file_main
  integer input
  integer ios
  character ( len = 80 ) line
  logical s_eqi
  character ( len = 80 ) word

  call get_unit ( input )

  open ( unit = input, file = file_main, status = 'old' )

  file_dist = ' '

  do
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else

      done = .true.

      call word_next_read ( line, word, done )

      if ( done ) then
        cycle
      end if

      if ( .not. s_eqi ( word, 'dist' ) ) then
        cycle
      end if

      call word_next_read ( line, file_dist, done )

      exit

    end if

  end do

  close ( unit = input )

  return
end
subroutine main_read_dms ( file_main, file_dms )

!*****************************************************************************80
!
!! MAIN_READ_DMS reads the name of the DMS file from the main file.
!
!  Discussion:
!
!    FILE_DMS is the name of a file containing the longitude and latitude
!    of each city in degrees/minutes/seconds.
!
!    There MAY be a record in the main file of the form
!
!    "dms  key_dms.txt"
!
!  Modified:
!
!    09 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_MAIN, the name of the file to read.
!
!    Output, character ( len = * ) FILE_DMS, the name of the DMS file,
!    or ' ' if no information was found.
!
  implicit none

  logical done
  character ( len = * ) file_dms
  character ( len = * ) file_main
  integer input
  integer ios
  character ( len = 80 ) line
  logical s_eqi
  character ( len = 80 ) word

  call get_unit ( input )

  open ( unit = input, file = file_main, status = 'old' )

  file_dms = ' '

  do
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else

      done = .true.

      call word_next_read ( line, word, done )

      if ( done ) then
        cycle
      end if

      if ( .not. s_eqi ( word, 'dms' ) ) then
        cycle
      end if

      call word_next_read ( line, file_dms, done )

      exit

    end if

  end do

  close ( unit = input )

  return
end
subroutine main_read_geom ( file_main, geom )

!*****************************************************************************80
!
!! MAIN_READ_GEOM reads the name of the geometry from the main file.
!
!  Discussion:
!
!    GEOM is the name of the geometry of the city data.
!    Typical values include:
!    none - no special geometry
!    plane - the points lie in a plane
!    sphere - the points lie on a sphere
!
!    There MAY be a record in the main file of the form
!
!    "geom  geom_value"
!
!  Modified:
!
!    07 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_MAIN, the name of the file to read.
!
!    Output, character ( len = * ) GEOM, the name of the geometry,
!    or ' ' if no information was found.
!
  implicit none

  logical done
  character ( len = * ) file_main
  character ( len = * ) geom
  integer input
  integer ios
  character ( len = 80 ) line
  logical s_eqi
  character ( len = 80 ) word

  call get_unit ( input )

  open ( unit = input, file = file_main, status = 'old' )

  geom = ' '

  do
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else

      done = .true.

      call word_next_read ( line, word, done )

      if ( done ) then
        cycle
      end if

      if ( .not. s_eqi ( word, 'geom' ) ) then
        cycle
      end if

      call word_next_read ( line, geom, done )

      exit

    end if

  end do

  close ( unit = input )

  return
end
subroutine main_read_name ( file_main, file_name )

!*****************************************************************************80
!
!! MAIN_READ_NAME reads the name of the name file from the main file.
!
!  Discussion:
!
!    FILE_NAME is the name of a file containing the city names.
!
!    There MAY be a record in the main file of the form
!
!    "name  key_name.txt"
!
!  Modified:
!
!    07 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_MAIN, the name of the file to read.
!
!    Output, character ( len = * ) FILE_NAME, the name of the name file,
!    or ' ' if no information was found.
!
  implicit none

  logical done
  character ( len = * ) file_main
  character ( len = * ) file_name
  integer input
  integer ios
  character ( len = 80 ) line
  logical s_eqi
  character ( len = 80 ) word

  call get_unit ( input )

  open ( unit = input, file = file_main, status = 'old' )

  file_name = ' '

  do
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else

      done = .true.

      call word_next_read ( line, word, done )

      if ( done ) then
        cycle
      end if

      if ( .not. s_eqi ( word, 'name' ) ) then
        cycle
      end if

      call word_next_read ( line, file_name, done )

      exit

    end if

  end do

  close ( unit = input )

  return
end
subroutine main_read_size ( file_main, n )

!*****************************************************************************80
!
!! MAIN_READ_SIZE reads the problem size N from the main file.
!
!  Discussion:
!
!    The problem size is N, the number of cities.
!
!    There should always be a record in the main file of the form
!
!    "size  7"
!
!  Modified:
!
!    07 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_MAIN, the name of the file to read.
!
!    Output, integer N, the problem size.
!
  implicit none

  logical done
  character ( len = * ) file_main
  integer ierror
  integer input
  integer ios
  integer length
  character ( len = 80 ) line
  integer n
  logical s_eqi
  character ( len = 80 ) word

  call get_unit ( input )

  open ( unit = input, file = file_main, status = 'old' )

  n = 0

  do
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else

      done = .true.

      call word_next_read ( line, word, done )

      if ( done ) then
        cycle
      end if

      if ( .not. s_eqi ( word, 'size' ) ) then
        cycle
      end if

      call word_next_read ( line, word, done )

      call s_to_i4 ( word, n, ierror, length )

      exit

    end if

  end do

  close ( unit = input )

  return
end
subroutine main_read_weight ( file_main, file_weight )

!*****************************************************************************80
!
!! MAIN_READ_WEIGHT reads the name of the weight file from the main file.
!
!  Discussion:
!
!    FILE_WEIGHT is the name of a file containing the weight
!    vector.
!
!    There MAY be a record in the main file of the form
!
!    "weight  key_weight.txt"
!
!  Modified:
!
!    09 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_MAIN, the name of the file to read.
!
!    Output, character ( len = * ) FILE_WEIGHT, the name of the weight file,
!    or ' ' if no information was found.
!
  implicit none

  logical done
  character ( len = * ) file_weight
  character ( len = * ) file_main
  integer input
  integer ios
  character ( len = 80 ) line
  logical s_eqi
  character ( len = 80 ) word

  call get_unit ( input )

  open ( unit = input, file = file_main, status = 'old' )

  file_weight = ' '

  do
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else

      done = .true.

      call word_next_read ( line, word, done )

      if ( done ) then
        cycle
      end if

      if ( .not. s_eqi ( word, 'weight' ) ) then
        cycle
      end if

      call word_next_read ( line, file_weight, done )

      exit

    end if

  end do

  close ( unit = input )

  return
end
subroutine main_read_xy ( file_main, file_xy )

!*****************************************************************************80
!
!! MAIN_READ_XY reads the name of the XY file from the main file.
!
!  Discussion:
!
!    FILE_XY is the name of a file containing (X,Y) coordinate data.
!
!    There MAY be a record in the main file of the form
!
!    "xy  key_xy.txt"
!
!  Modified:
!
!    07 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_MAIN, the name of the file to read.
!
!    Output, character ( len = * ) FILE_XY, the name of the XY file,
!    or ' ' if no information was found.
!
  implicit none

  logical done
  character ( len = * ) file_main
  character ( len = * ) file_xy
  integer input
  integer ios
  character ( len = 80 ) line
  logical s_eqi
  character ( len = 80 ) word

  call get_unit ( input )

  open ( unit = input, file = file_main, status = 'old' )

  file_xy = ' '

  do
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else

      done = .true.

      call word_next_read ( line, word, done )

      if ( done ) then
        cycle
      end if

      if ( .not. s_eqi ( word, 'xy' ) ) then
        cycle
      end if

      call word_next_read ( line, file_xy, done )

      exit

    end if

  end do

  close ( unit = input )

  return
end
subroutine points_dist_sphere ( lat1, long1, lat2, long2, radius, dist )

!*****************************************************************************80
!
!! POINTS_DIST_SPHERE finds the distance between two points on a sphere.
!
!  Discussion:
!
!    The positions of the the points are given as longitude and
!    latitude, in radians.
!
!    The distance is measured on the surface of the sphere.
!
!  Modified:
!
!    05 June 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) LAT1, LONG1, LAT2, LONG2, the latitude and
!    longitude of the two points, in radians.
!
!    Input, real ( kind = 8 ) RADIUS, the radius of the sphere.
!
!    Output, real ( kind = 8 ) DIST, the distance between the points.
!
  implicit none

  real ( kind = 8 ) dist
  real ( kind = 8 ) lat1
  real ( kind = 8 ) lat2
  real ( kind = 8 ) long1
  real ( kind = 8 ) long2
  real ( kind = 8 ) radius
  real ( kind = 8 ) theta

  theta = acos ( sin ( lat1 ) * sin ( lat2 ) &
               + cos ( lat1 ) * cos ( lat2 ) * cos ( long1 - long2 ) )

  dist = radius * theta

  return
end
subroutine r8mat_print ( m, n, a, title )

!*****************************************************************************80
!
!! R8MAT_PRINT prints an R8MAT, with an optional title.
!
!  Modified:
!
!    02 October 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer M, the number of rows in A.
!
!    Input, integer N, the number of columns in A.
!
!    Input, real ( kind = 8 ) A(M,N), the matrix.
!
!    Input, character ( len = * ) TITLE, a title to be printed.
!
  implicit none

  integer m
  integer n

  real ( kind = 8 ) a(m,n)
  integer i
  integer j
  integer jhi
  integer jlo
  character ( len = * ) title

  if ( title /= ' ' ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  if ( maxval ( abs ( a(1:m,1:n) ) ) < 1000000.0D+00 ) then

    do jlo = 1, n, 5
      jhi = min ( jlo + 4, n )
      write ( *, '(a)' ) ' '
      write ( *, '(6x,5(i7,7x))' ) ( j, j = jlo, jhi )
      write ( *, '(a)' ) ' '
      do i = 1, m
        write ( *, '(i6,5f14.6)' ) i, a(i,jlo:jhi)
      end do
    end do

  else

    do jlo = 1, n, 5
      jhi = min ( jlo + 4, n )
      write ( *, '(a)' ) ' '
      write ( *, '(6x,5(i7,7x))' ) ( j, j = jlo, jhi )
      write ( *, '(a)' ) ' '
      do i = 1, m
        write ( *, '(i6,5g14.6)' ) i, a(i,jlo:jhi)
      end do
    end do

  end if

  return
end
subroutine r8vec_print ( n, a, title )

!*****************************************************************************80
!
!! R8VEC_PRINT prints a real vector.
!
!  Discussion:
!
!    If all the entries are integers, the data if printed
!    in integer format.
!
!  Modified:
!
!    19 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components of the vector.
!
!    Input, real ( kind = 8 ) A(N), the vector to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  implicit none

  integer n

  real ( kind = 8 ) a(n)
  integer i
  character ( len = * ) title

  if ( title /= ' ' ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, '(a)' ) ' '

  if ( all ( a(1:n) == aint ( a(1:n) ) ) ) then
    do i = 1, n
      write ( *, '(i6,i6)' ) i, int ( a(i) )
    end do
  else if ( all ( abs ( a(1:n) ) < 1000000.0D+00 ) ) then
    do i = 1, n
      write ( *, '(i6,f14.6)' ) i, a(i)
    end do
  else
    do i = 1, n
      write ( *, '(i6,g14.6)' ) i, a(i)
    end do
  end if

  return
end
function s_eqi ( s1, s2 )

!*****************************************************************************80
!
!! S_EQI is a case insensitive comparison of two strings for equality.
!
!  Examples:
!
!    S_EQI ( 'Anjana', 'ANJANA' ) is .TRUE.
!
!  Modified:
!
!    14 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S1, S2, the strings to compare.
!
!    Output, logical S_EQI, the result of the comparison.
!
  implicit none

  character c1
  character c2
  integer i
  integer len1
  integer len2
  integer lenc
  logical s_eqi
  character ( len = * ) s1
  character ( len = * ) s2

  len1 = len ( s1 )
  len2 = len ( s2 )
  lenc = min ( len1, len2 )

  s_eqi = .false.

  do i = 1, lenc

    c1 = s1(i:i)
    c2 = s2(i:i)
    call ch_cap ( c1 )
    call ch_cap ( c2 )

    if ( c1 /= c2 ) then
      return
    end if

  end do

  do i = lenc + 1, len1
    if ( s1(i:i) /= ' ' ) then
      return
    end if
  end do

  do i = lenc + 1, len2
    if ( s2(i:i) /= ' ' ) then
      return
    end if
  end do

  s_eqi = .true.

  return
end
subroutine s_rep_one ( s1, sub1, sub2, s2 )

!*****************************************************************************80
!
!! S_REP_ONE replaces the first occurrence of SUB1 with SUB2.
!
!  Discussion:
!
!    The input and output strings may coincide.
!
!  Modified:
!
!    07 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S1, the initial string.
!
!    Input, character ( len = * ) SUB1, the string to be replaced.
!
!    Input, character ( len = * ) SUB2, the replacement string.
!
!    Output, character ( len = * ) S2, the final string.
!
  implicit none

  integer i1
  integer i2
  integer i3
  integer i4
  character ( len = * ) s1
  character ( len = * ) s2
  character ( len = 255 ) s3
  character ( len = * ) sub1
  character ( len = * ) sub2

  s3 = ' '

  i1 = index ( s1, sub1 )

  if ( i1 == 0 ) then

    s3 = s1

  else

    s3(1:i1-1) = s1(1:i1-1)

    i2 = len_trim ( sub2 )
    s3(i1:i1+i2-1) = sub2(1:i2)

    i3 = i1 + len_trim ( sub1 )
    i4 = len_trim ( s1 )

    s3(i1+i2:i1+i2+1+i4-i3) = s1(i3:i4)

  end if

  s2 = s3

  return
end
subroutine s_to_i4 ( s, ival, ierror, last )

!*****************************************************************************80
!
!! S_TO_I4 reads an I4 from a string.
!
!  Modified:
!
!    28 June 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, a string to be examined.
!
!    Output, integer IVAL, the integer value read from the string.
!    If the string is blank, then IVAL will be returned 0.
!
!    Output, integer IERROR, an error flag.
!    0, no error.
!    1, an error occurred.
!
!    Output, integer LAST, the last character of S used to make IVAL.
!
  implicit none

  character c
  integer i
  integer ierror
  integer isgn
  integer istate
  integer ival
  integer last
  character ( len = * ) s

  ierror = 0
  istate = 0
  isgn = 1
  ival = 0

  do i = 1, len_trim ( s )

    c = s(i:i)
!
!  Haven't read anything.
!
    if ( istate == 0 ) then

      if ( c == ' ' ) then

      else if ( c == '-' ) then
        istate = 1
        isgn = -1
      else if ( c == '+' ) then
        istate = 1
        isgn = + 1
      else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
        istate = 2
        ival = ichar ( c ) - ichar ( '0' )
      else
        ierror = 1
        return
      end if
!
!  Have read the sign, expecting digits.
!
    else if ( istate == 1 ) then

      if ( c == ' ' ) then

      else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
        istate = 2
        ival = ichar ( c ) - ichar ( '0' )
      else
        ierror = 1
        return
      end if
!
!  Have read at least one digit, expecting more.
!
    else if ( istate == 2 ) then

      if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then
        ival = 10 * ival + ichar ( c ) - ichar ( '0' )
      else
        ival = isgn * ival
        last = i - 1
        return
      end if

    end if

  end do
!
!  If we read all the characters in the string, see if we're OK.
!
  if ( istate == 2 ) then
    ival = isgn * ival
    last = len_trim ( s )
  else
    ierror = 1
    last = 0
  end if

  return
end
subroutine s_to_r8 ( s, r, ierror, lchar )

!*****************************************************************************80
!
!! S_TO_R8 reads an R8 from a string.
!
!  Discussion:
!
!    This routine will read as many characters as possible until it reaches
!    the end of the string, or encounters a character which cannot be
!    part of the real number.
!
!    Legal input is:
!
!       1 blanks,
!       2 '+' or '-' sign,
!       2.5 spaces
!       3 integer part,
!       4 decimal point,
!       5 fraction part,
!       6 'E' or 'e' or 'D' or 'd', exponent marker,
!       7 exponent sign,
!       8 exponent integer part,
!       9 exponent decimal point,
!      10 exponent fraction part,
!      11 blanks,
!      12 final comma or semicolon.
!
!    with most quantities optional.
!
!  Examples:
!
!    S                 R
!
!    '1'               1.0
!    '     1   '       1.0
!    '1A'              1.0
!    '12,34,56'        12.0
!    '  34 7'          34.0
!    '-1E2ABCD'        -100.0
!    '-1X2ABCD'        -1.0
!    ' 2E-1'           0.2
!    '23.45'           23.45
!    '-4.2E+2'         -420.0
!    '17d2'            1700.0
!    '-14e-2'         -0.14
!    'e2'              100.0
!    '-12.73e-9.23'   -12.73 * 10.0**(-9.23)
!
!  Modified:
!
!    12 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, the string containing the
!    data to be read.  Reading will begin at position 1 and
!    terminate at the end of the string, or when no more
!    characters can be read to form a legal real.  Blanks,
!    commas, or other nonnumeric data will, in particular,
!    cause the conversion to halt.
!
!    Output, real ( kind = 8 ) R, the real value that was read from the string.
!
!    Output, integer IERROR, error flag.
!
!    0, no errors occurred.
!
!    1, 2, 6 or 7, the input number was garbled.  The
!    value of IERROR is the last type of input successfully
!    read.  For instance, 1 means initial blanks, 2 means
!    a plus or minus sign, and so on.
!
!    Output, integer LCHAR, the number of characters read from
!    the string to form the number, including any terminating
!    characters such as a trailing comma or blanks.
!
  implicit none

  character c
  logical ch_eqi
  integer ierror
  integer ihave
  integer isgn
  integer iterm
  integer jbot
  integer jsgn
  integer jtop
  integer lchar
  integer nchar
  integer ndig
  real ( kind = 8 ) r
  real ( kind = 8 ) rbot
  real ( kind = 8 ) rexp
  real ( kind = 8 ) rtop
  character ( len = * ) s
  character, parameter :: TAB = char ( 9 )

  nchar = len_trim ( s )
  ierror = 0
  r = 0.0D+00
  lchar = - 1
  isgn = 1
  rtop = 0.0D+00
  rbot = 1.0D+00
  jsgn = 1
  jtop = 0
  jbot = 1
  ihave = 1
  iterm = 0

  do

    lchar = lchar + 1
    c = s(lchar+1:lchar+1)
!
!  Blank or TAB character.
!
    if ( c == ' ' .or. c == TAB ) then

      if ( ihave == 2 ) then

      else if ( ihave == 6 .or. ihave == 7 ) then
        iterm = 1
      else if ( 1 < ihave ) then
        ihave = 11
      end if
!
!  Comma.
!
    else if ( c == ',' .or. c == ';' ) then

      if ( ihave /= 1 ) then
        iterm = 1
        ihave = 12
        lchar = lchar + 1
      end if
!
!  Minus sign.
!
    else if ( c == '-' ) then

      if ( ihave == 1 ) then
        ihave = 2
        isgn = - 1
      else if ( ihave == 6 ) then
        ihave = 7
        jsgn = - 1
      else
        iterm = 1
      end if
!
!  Plus sign.
!
    else if ( c == '+' ) then

      if ( ihave == 1 ) then
        ihave = 2
      else if ( ihave == 6 ) then
        ihave = 7
      else
        iterm = 1
      end if
!
!  Decimal point.
!
    else if ( c == '.' ) then

      if ( ihave < 4 ) then
        ihave = 4
      else if ( 6 <= ihave .and. ihave <= 8 ) then
        ihave = 9
      else
        iterm = 1
      end if
!
!  Exponent marker.
!
    else if ( ch_eqi ( c, 'E' ) .or. ch_eqi ( c, 'D' ) ) then

      if ( ihave < 6 ) then
        ihave = 6
      else
        iterm = 1
      end if
!
!  Digit.
!
    else if ( ihave < 11 .and. lge ( c, '0' ) .and. lle ( c, '9' ) ) then

      if ( ihave <= 2 ) then
        ihave = 3
      else if ( ihave == 4 ) then
        ihave = 5
      else if ( ihave == 6 .or. ihave == 7 ) then
        ihave = 8
      else if ( ihave == 9 ) then
        ihave = 10
      end if

      call ch_to_digit ( c, ndig )

      if ( ihave == 3 ) then
        rtop = 10.0D+00 * rtop + real ( ndig, kind = 8 )
      else if ( ihave == 5 ) then
        rtop = 10.0D+00 * rtop + real ( ndig, kind = 8 )
        rbot = 10.0D+00 * rbot
      else if ( ihave == 8 ) then
        jtop = 10 * jtop + ndig
      else if ( ihave == 10 ) then
        jtop = 10 * jtop + ndig
        jbot = 10 * jbot
      end if
!
!  Anything else is regarded as a terminator.
!
    else
      iterm = 1
    end if
!
!  If we haven't seen a terminator, and we haven't examined the
!  entire string, go get the next character.
!
    if ( iterm == 1 .or. nchar <= lchar+1 ) then
      exit
    end if

  end do
!
!  If we haven't seen a terminator, and we have examined the
!  entire string, then we're done, and LCHAR is equal to NCHAR.
!
  if ( iterm /= 1 .and. lchar+1 == nchar ) then
    lchar = nchar
  end if
!
!  Number seems to have terminated.  Have we got a legal number?
!  Not if we terminated in states 1, 2, 6 or 7!
!
  if ( ihave == 1 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then

    ierror = ihave

    return
  end if
!
!  Number seems OK.  Form it.
!
  if ( jtop == 0 ) then
    rexp = 1.0D+00
  else

    if ( jbot == 1 ) then
      rexp = 10.0D+00**( jsgn * jtop )
    else
      rexp = jsgn * jtop
      rexp = rexp / jbot
      rexp = 10.0D+00**rexp
    end if

  end if

  r = isgn * rexp * rtop / rbot

  return
end
subroutine s_to_w ( s, w, ierror, last )

!*****************************************************************************80
!
!! S_TO_W reads the next blank-delimited word from a string.
!
!  Modified:
!
!    15 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, a string to be examined.
!
!    Output, character ( len = * ) W, the word that was read.
!
!    Output, integer IERROR, an error flag.
!    0, no error.
!    1, an error occurred.
!
!    Output, integer LAST, the last character of S used to make W.
!
  implicit none

  character c
  integer first
  integer i
  integer ierror
  integer istate
  integer last
  character ( len = * ) s
  character ( len = * ) w

  w = ' '
  ierror = 0
  istate = 0
  first = 0
  last = 0
  i = 0

  do

    i = i + 1

    if ( len_trim ( s ) < i ) then

      if ( istate == 0 ) then
        ierror = 1
        last = 0
      else
        last = i-1
        w = s(first:last)
      end if

      exit

    end if

    c = s(i:i)

    if ( istate == 0 ) then

      if ( c /= ' ' ) then
        first = i
        istate = 1
      end if

    else if ( istate == 1 ) then

      if ( c == ' ' ) then
        last = i - 1
        w = s(first:last)
        exit
      end if

    end if

  end do

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end
subroutine word_count ( s, nword )

!*****************************************************************************80
!
!! WORD_COUNT counts the number of "words" in a string.
!
!  Modified:
!
!    14 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, the string to be examined.
!
!    Output, integer NWORD, the number of "words" in the string.
!    Words are presumed to be separated by one or more blanks.
!
  implicit none

  logical blank
  integer i
  integer lens
  integer nword
  character ( len = * ) s
!
  nword = 0
  lens = len ( s )

  if ( lens <= 0 ) then
    return
  end if

  blank = .true.

  do i = 1, lens

    if ( s(i:i) == ' ' ) then
      blank = .true.
    else if ( blank ) then
      nword = nword + 1
      blank = .false.
    end if

  end do

  return
end
subroutine word_next_read ( s, word, done )

!*****************************************************************************80
!
!! WORD_NEXT_READ "reads" words from a string, one at a time.
!
!  Special cases:
!
!    The following characters are considered to be a single word,
!    whether surrounded by spaces or not:
!
!      " ( ) { } [ ]
!
!    Also, if there is a trailing comma on the word, it is stripped off.
!    This is to facilitate the reading of lists.
!
!  Modified:
!
!    23 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, a string, presumably containing words
!    separated by spaces.
!
!    Output, character ( len = * ) WORD.
!
!    If DONE is FALSE, then WORD contains the "next" word read.
!    If DONE is TRUE, then WORD is blank, because there was no more to read.
!
!    Input/output, logical DONE.
!
!    On input with a fresh string, set DONE to TRUE.
!
!    On output, the routine sets DONE:
!      FALSE if another word was read,
!      TRUE if no more words could be read.
!
  implicit none

  logical done
  integer ilo
  integer, save :: lenc = 0
  integer, save :: next = 1
  character ( len = * ) s
  character, parameter :: TAB = char ( 9 )
  character ( len = * ) word
!
!  We "remember" LENC and NEXT from the previous call.
!
!  An input value of DONE = TRUE signals a new line of text to examine.
!
  if ( done ) then

    next = 1
    done = .false.
    lenc = len_trim ( s )

    if ( lenc <= 0 ) then
      done = .true.
      word = ' '
      return
    end if

  end if
!
!  Beginning at index NEXT, search the string for the next nonblank,
!  which signals the beginning of a word.
!
  ilo = next
!
!  ...S(NEXT:) is blank.  Return with WORD = ' ' and DONE = TRUE.
!
  do

    if ( lenc < ilo ) then
      word = ' '
      done = .true.
      next = lenc + 1
      return
    end if
!
!  If the current character is blank, skip to the next one.
!
    if ( s(ilo:ilo) /= ' ' .and. s(ilo:ilo) /= TAB ) then
      exit
    end if

    ilo = ilo + 1

  end do
!
!  ILO is the index of the next nonblank character in the string.
!
!  If this initial nonblank is a special character,
!  then that's the whole word as far as we're concerned,
!  so return immediately.
!
  if ( s(ilo:ilo) == '"' .or. &
       s(ilo:ilo) == '(' .or. &
       s(ilo:ilo) == ')' .or. &
       s(ilo:ilo) == '{' .or. &
       s(ilo:ilo) == '}' .or. &
       s(ilo:ilo) == '[' .or. &
       s(ilo:ilo) == ']' ) then

    word = s(ilo:ilo)
    next = ilo + 1
    return

  end if
!
!  Now search for the last contiguous character that is not a
!  blank, TAB, or special character.
!
  next = ilo + 1

  do while ( next <= lenc )

    if ( s(next:next) == ' ' ) then
      exit
    else if ( s(next:next) == TAB ) then
      exit
    else if ( s(next:next) == '"' ) then
      exit
    else if ( s(next:next) == '(' ) then
      exit
    else if ( s(next:next) == ')' ) then
      exit
    else if ( s(next:next) == '{' ) then
      exit
    else if ( s(next:next) == '}' ) then
      exit
    else if ( s(next:next) == '[' ) then
      exit
    else if ( s(next:next) == ']' ) then
      exit
    end if

    next = next + 1

  end do
!
!  Ignore a trailing comma.
!
  if ( s(next-1:next-1) == ',' ) then
    word = s(ilo:next-2)
  else
    word = s(ilo:next-1)
  end if

  return
end
subroutine weight_read ( file_name, n, weight )

!*****************************************************************************80
!
!! WEIGHT_READ reads weights from a file.
!
!  Discussion:
!
!    The data is stored in a file, with each weight on one line.
!
!    Blank lines and comment lines (beginning with '#') are ignored.
!
!    Individual data values should be separated by spaces or commas.
!
!  Modified:
!
!    08 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file to read.
!
!    Input, integer N, the number of data items.
!
!    Output, real ( kind = 8 ) WEIGHT(N), the weights.
!
  implicit none

  integer n
  integer, parameter :: ncol = 1

  character ( len = * ) file_name
  integer i
  integer ierror
  integer input
  integer ios
  integer j
  integer last
  integer length
  character ( len = 80 ) line
  integer line_num
  real ( kind = 8 ) value
  real ( kind = 8 ) weight(n)

  call get_unit ( input )

  open ( unit = input, file = file_name, status = 'old' )

  weight(1:n) = huge ( weight(1) )

  i = 1
  j = 0
  line_num = 0

  do
!
!  Have we read enough data?
!
    if ( i == n .and. j == ncol ) then
      exit
    end if
!
!  Have we read too much data?
!
    if ( n < i .or. ncol < j ) then
      exit
    end if
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if

    line_num = line_num + 1
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else
!
!  LAST points to the last character associated with the previous 
!  data value read from the line.
!
      last = 0

      do
!
!  Try to read another value from the line.
!
        call s_to_r8 ( line(last+1:), value, ierror, length )
!
!  If we could not read a new value, it's time to read a new line.
!
        if ( ierror /= 0 ) then
          exit
        end if
!
!  Update the pointer.
!
        last = last + length
!
!  If we read a new value, where do we put it?
!
        j = j + 1

        if ( ncol < j ) then
          j = 1
          i = i + 1
          if ( n < i ) then
            exit
          end if
        end if

        weight(i) = value
!
!  If you reached the end of the row, it's time to read a new line.
!
        if ( j == ncol ) then
          exit
        end if

      end do

    end if

  end do

  close ( unit = input )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'WEIGHT_READ:'
  write ( *, '(a,i6)' ) '  Number of lines read was ', line_num

  return
end
subroutine xy_print ( n, xy, title )

!*****************************************************************************80
!
!! XY_PRINT prints XY planar locations.
!
!  Modified:
!
!    07 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of data items.
!
!    Input, real ( kind = 8 ) XY(N,2), the data values.
!
!    Input, character ( len = * ) TITLE, an optional title.
!
  implicit none

  integer n

  integer i
  character ( len = * ) title
  real ( kind = 8 ) xy(n,2)

  if ( len_trim ( title ) /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(i3,2x,2g14.6)' ) i, xy(i,1:2)
  end do

  return
end
subroutine xy_read ( file_name, n, xy )

!*****************************************************************************80
!
!! XY_READ reads XY planar locations from a file.
!
!  Discussion:
!
!    The data is stored in a file, with each (X,Y) pair on one line.
!
!    Blank lines and comment lines (beginning with '#') are ignored.
!
!    Individual data values should be separated by spaces or commas.
!
!  Modified:
!
!    07 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file to read.
!
!    Input, integer N, the number of data items.
!
!    Output, real ( kind = 8 ) XY(N,2), the data values.
!
  implicit none

  integer n
  integer, parameter :: ncol = 2

  character ( len = * ) file_name
  integer i
  integer ierror
  integer input
  integer ios
  integer j
  integer last
  integer length
  character ( len = 80 ) line
  integer line_num
  real ( kind = 8 ) value
  real ( kind = 8 ) xy(n,ncol)

  call get_unit ( input )

  open ( unit = input, file = file_name, status = 'old' )

  xy(1:n,1:ncol) = huge ( xy(1,1) )

  i = 1
  j = 0
  line_num = 0

  do
!
!  Have we read enough data?
!
    if ( i == n .and. j == ncol ) then
      exit
    end if
!
!  Have we read too much data?
!
    if ( n < i .or. ncol < j ) then
      exit
    end if
!
!  Read the next line from the file.
!
    read ( input, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if

    line_num = line_num + 1
!
!  Skip blank lines and comment lines.
!
    if ( len_trim ( line ) == 0 ) then

    else if ( line(1:1) == '#' ) then

    else
!
!  LAST points to the last character associated with the previous 
!  data value read from the line.
!
      last = 0

      do
!
!  Try to read another value from the line.
!
        call s_to_r8 ( line(last+1:), value, ierror, length )
!
!  If we could not read a new value, it's time to read a new line.
!
        if ( ierror /= 0 ) then
          exit
        end if
!
!  Update the pointer.
!
        last = last + length
!
!  If we read a new value, where do we put it?
!
        j = j + 1

        if ( ncol < j ) then
          j = 1
          i = i + 1
          if ( n < i ) then
            exit
          end if
        end if

        xy(i,j) = value
!
!  If you reached the end of the row, it's time to read a new line.
!
        if ( j == ncol ) then
          exit
        end if

      end do

    end if

  end do

  close ( unit = input )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'XY_READ:'
  write ( *, '(a,i6)' ) '  Number of lines read was ', line_num

  return
end
subroutine xy_to_dist ( n, xy, dist )

!*****************************************************************************80
!
!! XY_TO_DIST creates a distance matrix from a list of (X,Y) plane positions.
!
!  Discussion:
!
!    The euclidean planar distance is used.
!
!  Modified:
!
!    07 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of data items.
!
!    Input, real ( kind = 8 ) XY(N,2), the X,Y locations of the data.
!
!    Output, real ( kind = 8 ) DIST(N,N), the distance matrix.
!
  implicit none

  integer n

  real ( kind = 8 ) dist(n,n)
  integer i
  integer j
  real ( kind = 8 ) xy(n,2)

  do i = 1, n
    dist(i,i) = 0.0D+0
    do j = i+1, n
      dist(i,j) = sqrt ( ( xy(i,1) - xy(j,1) )**2 &
                       + ( xy(i,2) - xy(j,2) )**2 )
      dist(j,i) = dist(i,j)
    end do
  end do

  return
end
subroutine xy_write ( file_name, n, xy )

!*****************************************************************************80
!
!! XY_WRITE writes XY coordinates to a file.
!
!  Modified:
!
!    13 November 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file to write.
!
!    Input, integer N, the number of data items.
!
!    Input, real ( kind = 8 ) XY(N,2), the X,Y coordinates.
!
  implicit none

  integer n

  character ( len = * ) file_name
  integer i
  integer output
  real ( kind = 8 ) xy(n,2)

  call get_unit ( output )

  open ( unit = output, file = file_name, status = 'replace' )

  write ( output, '(a)' ) '# ' // trim ( file_name )
  write ( output, '(a)' ) '#'
  write ( output, '(a)' ) '#  Created by XY_WRITE.'
  write ( output, '(a)' ) '#'
  write ( output, '(a)' ) '#  N by 2 XY coordinate list.'
  write ( output, '(a,i6)' ) '#  Number of points N is ', n
  write ( output, '(a)' ) '#'

  if ( any ( xy(1:n,1:2) /= aint ( xy(1:n,1:2) ) ) ) then

    do i = 1, n

      write ( output, '(2g14.6)' ) xy(i,1:2)

    end do

  else

    do i = 1, n

      write ( output, '(2i6)' ) int ( xy(i,1:2) )

    end do

  end if

  close ( unit = output )

  return
end