program ranmap
!
!*******************************************************************************
!
!! RANMAP carries out a weighted affine mapping iteration.
!
!
!  Discussion:
!
!    The map coefficients for the dragon are "slightly" wrong, and
!    those for the tree seem very wrong.
!
!    The random map coefficients must be processed to guarantee that
!    they map [0,1]x[0,1] into itself.
!
!  Modified:
!
!    08 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Michael Barnsley,
!    Fractals Everywhere,
!    Academic Press, 1988.
!
!    Michael Barnsley and Lyman Hurd,
!    Fractal Image Compression,
!    Peters, 1993.
!
!    Michael Barnsley and Alan Sloan,
!    A Better Way to Compress Images,
!    Byte Magazine, January 1988.
!
!    Bernt Wahl,
!    Exploring Fractals on the Mac,
!    Addison Wesley, 1995.
!
  implicit none
!
  real a
  real b
  real, allocatable, dimension (:,:) :: coord
  real, save, dimension ( 2 ) :: coord_max = (/ 1.0E+00, 1.0E+00 /)
  real, save, dimension ( 2 ) :: coord_min = (/ 0.0E+00, 0.0E+00 /)
  character ( len = 80 ) :: filedot_name = 'ranmap.ps'
  character ( len = 80 ) :: filepoint_name = 'ranmap.dat'
  integer i
  integer iunit
  integer k
  integer map
  real, allocatable, dimension ( : ) :: map_bin
  real, allocatable, dimension (:,:,:) :: map_matrix
  real, allocatable, dimension ( : ) :: map_weight
  integer n
  integer nmap
  logical rvec_ascends
  real s
  real weight_sum
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'RANMAP'
  write ( *, '(a)' ) '  Demonstrate the use of weighted affine'
  write ( *, '(a)' ) '  mappings to draw fractals.'
!
!  Get the number of points to plot.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Choose the number of points to plot'
  read ( *, * ) n
!
!  Get the map option.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Choose a pattern:'
  write ( *, '(a)' ) '  1 = cross'
  write ( *, '(a)' ) '  2 = leaf'
  write ( *, '(a)' ) '  3 = dragon'
  write ( *, '(a)' ) '  4 = triangle'
  write ( *, '(a)' ) '  5 = fern'
  write ( *, '(a)' ) '  6 = tree'
  write ( *, '(a)' ) '  7 = a random map (not debugged)'
  write ( *, '(a)' ) '  8 = specify your own.'

  read ( *, * ) map
!
!  Set the number of maps.
!
  if ( map == 1 ) then
    nmap = 5
  else if ( map == 2 ) then
    nmap = 4
  else if ( map == 3 ) then
    nmap = 2
  else if ( map == 4 ) then
    nmap = 3
  else if ( map == 5 ) then
    nmap = 4
  else if ( map == 6 ) then
    nmap = 4
  else if ( map == 7 ) then
    call i_random ( 2, 6, nmap )
  else if ( map == 8 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'How many maps will be used?'
    read ( *, * ) nmap
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of maps to use is ', nmap
!
!  Set the map bins.
!
  allocate ( map_bin(1:nmap-1) )

  if ( map == 1 ) then
    map_bin(1:nmap-1) = (/ 0.20, 0.40, 0.60, 0.80 /)
  else if ( map == 2 ) then
    map_bin(1:nmap-1) = (/ 0.25, 0.50, 0.75 /)
  else if ( map == 3 ) then
    map_bin(1:nmap-1) = (/ 0.50 /)
  else if ( map == 4 ) then
    map_bin(1:nmap-1) = (/ 0.33, 0.66 /)
  else if ( map == 5 ) then
    map_bin(1:nmap-1) = (/ 0.75, 0.85, 0.95 /)
  else if ( map == 6 ) then
    map_bin(1:nmap-1) = (/ 0.05, 0.20, 0.60 /)

  else if ( map == 7 ) then

    allocate ( map_weight(1:nmap) )
    call random_number ( harvest = map_weight(1:nmap) )
    weight_sum = sum ( map_weight(1:nmap) )
    map_weight(1:nmap) = map_weight(1:nmap) / weight_sum
    s = 0.0
    do i = 1, nmap-1
      s = s + map_weight(i)
      map_bin(i) = s
    end do
    deallocate ( map_weight )

  else if ( map == 8 ) then

    do

      write ( *, '(a)' ) ' '
      write ( *, '(a,i6,a)' ) '  Enter an increasing set of ', nmap-1, ' values.'
      write ( *, '(a)' ) '  These, along with 0 and 1, define the intervals'
      write ( *, '(a)' ) '  associated with each map.'
      write ( *, '(a)' ) ' '
      do i = 1, nmap-1
        read ( *, * ) map_bin(i)
      end do

      if ( .not. rvec_ascends ( nmap-1, map_bin ) ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Your bin values are not acceptable because'
        write ( *, '(a)' ) '  they are not ascending.'
        cycle
      end if

      if ( map_bin(1) < 0.0E+00 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Your bin values are not acceptable because'
        write ( *, '(a)' ) '  MAP_BIN(1) < 0.'
        cycle
      end if

      if ( map_bin(nmap-1) > 1.0E+00 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Your bin values are not acceptable because'
        write ( *, '(a)' ) '  MAP_BIN(NMAP-1) > 1.'
        cycle
      end if

      exit

    end do

  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The map bin values:'
  write ( *, '(a)' ) ' '
  do i = 1, nmap
    if ( i == 1 ) then
      a = 0.0E+00
    else
      a = map_bin(i-1)
    end if
    if ( i < nmap ) then
      b = map_bin(i)
    else
      b = 1.0E+00
    end if
    write ( *, '(3f8.4)' ) a, b, b - a
  end do
!
!  Set the map coefficients.
!
  allocate ( map_matrix(2,3,nmap) )

  if ( map == 1 ) then
    call cross_map ( map_matrix )
  else if ( map == 2 ) then
    call leaf_map ( map_matrix )
  else if ( map == 3 ) then
    call dragon_map ( map_matrix )
  else if ( map == 4 ) then
    call triangle_map ( map_matrix )
  else if ( map == 5 ) then
    call fern_map ( map_matrix )
  else if ( map == 6 ) then
    call tree_map ( map_matrix )

  else if ( map == 7 ) then

    call random_number ( harvest = map_matrix(1:2,1:3,1:nmap) )

  else if ( map == 8 ) then

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  For each map, enter 2 rows of 3 values.'
    write ( *, '(a)' ) ' '
    do k = 1, nmap
      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) '  For map K = ', k
      write ( *, '(a)' ) ' '
      do i = 1, 2
        read ( *, * ) map_matrix(i,1:3,k)
      end do
    end do
  end if

  do k = 1, nmap
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) ' Map # ', k
    write ( *, '(a)' ) ' '
    do i = 1, 2
      write ( *, '(3f8.4)' ) map_matrix(i,1:3,k)
    end do
  end do

  allocate ( coord(2,n) )
!
!  Open a file, carry out the mapping, and close the file.
!
  call get_unit ( iunit )

  open ( unit = iunit, file = filepoint_name, status = 'replace' )

  call random_map ( iunit, n, nmap, map_bin, map_matrix, coord )

  close ( unit = iunit )

  deallocate ( map_bin )

  call dot_plot ( filedot_name, n, coord, coord_min, coord_max )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  X coordinate range:'
  write ( *, '(2g14.6)' ) minval ( coord(1,1:n) ), maxval ( coord(1,1:n) )
  write ( *, '(a)' )  ' '
  write ( *, '(a)' ) '  Y coordinate range:'
  write ( *, '(2g14.6)' ) minval ( coord(2,1:n) ), maxval ( coord(2,1:n) )

  deallocate ( coord )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The point data has been saved in the file ' // &
    trim ( filepoint_name )
  write ( *, '(a)' ) '  A point plot has been saved in the file ' // &
    trim ( filedot_name )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'RANMAP'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine cross_map ( map_matrix )
!
!*******************************************************************************
!
!! CROSS_MAP sets the cross map matrix.
!
!
!  Modified:
!
!    08 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real MAP_MATRIX(2,3,NMAP), the set of NMAP affine maps.
!    Each map has the form of a 2 by 3 matrix, so that
!      Xnew = A11 * X + A12 * Y + A13
!      Ynew = A21 * X + A22 * Y + A23
!
  implicit none
!
  real map_matrix(2,3,5)
!
  map_matrix(1:2,1:3,1) = reshape ( source = &
    (/ 0.333, 0.000, &
       0.000, 0.333, &
       0.333, 0.000 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,2) = reshape ( source = &
    (/ 0.333, 0.000, &
       0.000, 0.333, &
       0.000, 0.333 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,3) = reshape ( source = &
    (/ 0.333, 0.000, &
       0.000, 0.333, &
       0.333, 0.333 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,4) = reshape ( source = &
    (/ 0.333, 0.000, &
       0.000, 0.333, &
       0.666, 0.333 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,5) = reshape ( source = &
    (/ 0.333, 0.000, &
       0.000, 0.333, &
       0.333, 0.666 /), shape = (/ 2, 3 /) )

  return
end
subroutine dot_plot ( filedot_name, num_pts, coord, coord_min, coord_max )
!
!*******************************************************************************
!
!! DOT_PLOT plots the individual points.
!
!
!  Modified:
!
!    24 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILEDOT_NAME, the name of the dot file.
!
!    Input, integer NUM_PTS, the number of points.
!
!    Input, real COORD(2,NUM_PTS), the point coordinates.
!
!    Input, real COORD_MIN(2), COORD_MAX(2), the data ranges.
!
  implicit none
!
  integer num_pts
!
  real coord(2,num_pts)
  real coord_max(2)
  real coord_min(2)
  character ( len = * ) filedot_name
  integer filedot_unit
  integer i
  integer ierror
!
  call get_unit ( filedot_unit )

  call ps_file_open ( filedot_name, filedot_unit, ierror )

  call eps_file_head ( filedot_name )

  call ps_page_head ( coord_min(1), coord_min(2), coord_max(1), coord_max(2) )

  do i = 1, num_pts
    call ps_mark_point ( coord(1,i), coord(2,i) )
  end do

  call ps_page_tail

  call eps_file_tail

  call ps_file_close ( filedot_unit )

  return
end
subroutine dragon_map ( map_matrix )
!
!*******************************************************************************
!
!! DRAGON_MAP sets the dragon map matrix.
!
!
!  Modified:
!
!    08 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real MAP_MATRIX(2,3,NMAP), the set of NMAP affine maps.
!    Each map has the form of a 2 by 3 matrix, so that
!      Xnew = A11 * X + A12 * Y + A13
!      Ynew = A21 * X + A22 * Y + A23
!
  implicit none
!
  real map_matrix(2,3,2)
!
  map_matrix(1:2,1:3,1) = reshape ( source = &
    (/ 0.500, -0.500, &
       0.500,  0.500, &
       0.125,  0.675 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,2) = reshape ( source = &
    (/ 0.500, -0.500, &
       0.500,  0.500, &
      -0.125,  0.325 /), shape = (/ 2, 3 /) )

  return
end
subroutine fern_map ( map_matrix )
!
!*******************************************************************************
!
!! FERN_MAP sets the fern map matrix.
!
!
!  Modified:
!
!    08 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real MAP_MATRIX(2,3,NMAP), the set of NMAP affine maps.
!    Each map has the form of a 2 by 3 matrix, so that
!      Xnew = A11 * X + A12 * Y + A13
!      Ynew = A21 * X + A22 * Y + A23
!
  implicit none
!
  real map_matrix(2,3,4)
!
  map_matrix(1:2,1:3,1) = reshape ( source = &
    (/ 0.8560, -0.0205, &
       0.0414,  0.8580, &
       0.0700,  0.1470 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,2) = reshape ( source = &
    (/ 0.2440, 0.1760, &
      -0.3850, 0.2240, &
       0.3930, 0.1020 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,3) = reshape ( source = &
    (/ -0.1440,  0.1810, &
        0.3900,  0.2590, &
        0.5270, -0.0140 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,4) = reshape ( source = &
    (/ 0.0000, 0.0310, &
       0.0000, 0.2160, &
       0.4860, 0.0500 /)  , shape = (/ 2, 3 /) )

  return
end
subroutine get_unit ( iunit )
!
!*******************************************************************************
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!
!  Discussion:
!
!    A "free" FORTRAN unit number is an integer between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!  Modified:
!
!    02 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5 and 6).
!
!    Otherwise, IUNIT is an integer between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
  implicit none
!
  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )

      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if

  end do

  return
end
subroutine i_random ( ilo, ihi, i )
!
!*******************************************************************************
!
!! I_RANDOM returns a random integer in a given range.
!
!
!  Modified:
!
!    23 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ILO, IHI, the minimum and maximum acceptable values.
!
!    Output, integer I, the randomly chosen integer.
!
  implicit none
!
  integer i
  integer ihi
  integer ilo
  real r
!
  call random_number ( harvest = r )

  i = ilo + int ( r * real ( ihi + 1 - ilo ) )
!
!  In case of oddball events at the boundary, enforce the limits.
!
  i = max ( i, ilo )
  i = min ( i, ihi )

  return
end
subroutine leaf_map ( map_matrix )
!
!*******************************************************************************
!
!! LEAF_MAP sets the leaf map matrix.
!
!
!  Modified:
!
!    08 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real MAP_MATRIX(2,3,NMAP), the set of NMAP affine maps.
!    Each map has the form of a 2 by 3 matrix, so that
!      Xnew = A11 * X + A12 * Y + A13
!      Ynew = A21 * X + A22 * Y + A23
!
  implicit none
!
  real map_matrix(2,3,4)
!
  map_matrix(1:2,1:3,1) = reshape ( source = &
    (/ 0.800, 0.000, &
       0.000, 0.800, &
       0.100, 0.040 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,2) = reshape ( source = &
    (/ 0.500, 0.000, &
       0.000, 0.500, &
       0.250, 0.400 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,3) = reshape ( source = &
    (/  0.355, 0.355, &
       -0.355, 0.355, &
        0.266, 0.078 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,4) = reshape ( source = &
    (/ 0.355, -0.355, &
       0.355,  0.355, &
       0.378,  0.434 /), shape = (/ 2, 3 /) )

  return
end
subroutine r_to_bin_uneven ( nbin, xbin, x, bin )
!
!*******************************************************************************
!
!! R_TO_BIN_UNEVEN places X in one of several unevenly spaced bins.
!
!
!  Discussion:
!
!    The XBIN array is assumed to be sorted.
!
!  Example:
!
!    NBIN = 5
!    XBIN(1:4) = (/ 0.0, 2.0, 8.0, 9.0 /)
!
!    so bins are
!
!    1  ( -Inf,   0 )
!    2  (    0,   2 )
!    3  (    2,   8 )
!    4  (    8,   9 )
!    5  (    9, Inf )
!
!    X   BIN
!
!   -7    1
!   -3    1
!    0    1
!    0.1  2
!    1    2
!    3    3
!    8    3
!    9.5  5
!   13    5
!
!  Modified:
!
!    04 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NBIN, the number of bins.
!
!    Input, real XBIN(NBIN-1), the dividing values for the bins.
!
!    Input, real X, a value to be placed in a bin.
!
!    Output, integer BIN, the index of the bin to which X is assigned.
!
  implicit none
!
  integer nbin
!
  integer bin
  real x
  real xbin(nbin-1)
!
  bin = 1

  do while ( bin < nbin )

    if ( x <= xbin(bin) ) then
      return
    end if

    bin = bin + 1

  end do

  return
end
subroutine random_map ( iunit, n, nmap, map_bin, map_matrix, coord )
!
!*******************************************************************************
!
!! RANDOM_MAP randomly applies one of several affine maps to a point.
!
!
!  Discussion:
!
!    The user specifies a set of affine maps on [0,1]x[0,1], and their
!    associated weights or selection probabilities.
!
!    This routine then repeatedly selects a map at random, and applies
!    it to the current point to get the next point.
!
!    The computed points are written to an open I/O unit.
!
!  Modified:
!
!    03 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IUNIT, the FORTRAN unit to which the points should
!    be written.
!
!    Input, integer N, the number of points to draw, or times to apply
!    the transformation.
!
!    Input, integer NMAP, the number of affine maps.
!
!    Input, real MAP_BIN(NMAP-1), the cumulative probabilities of the maps.
!    The entries of MAP_BIN should be nonnegative, increasing, and no greater
!    than 1.
!    MAP_BIN(1) is the probability that map 1 will be chosen.
!    MAP_BIN(2) is the probability that map 1 or 2 will be chosen, and so on.
!    MAP_BIN(NMAP) is implicitly 1.
!
!    Input, real MAP_MATRIX(2,3,NMAP), the set of NMAP affine maps.
!    Each map has the form of a 2 by 3 matrix, so that
!      Xnew = A11 * X + A12 * Y + A13
!      Ynew = A21 * X + A22 * Y + A23
!
!    Output, real COORD(2,N), the coordinates of the points.
!
  implicit none
!
  integer n
  integer nmap
!
  real coord(2,n)
  integer i
  integer iunit
  integer map
  real map_matrix(2,3,nmap)
  real map_bin(nmap-1)
  real p
  real x(3)
!
!  Pick a random starting point.
!
!  The third component is 1, to allow us to implement affine mapping
!  using matrix multiplication.
!
  call random_number ( harvest = x(1:2) )
  x(3) = 1.0E+00
!
!  Choose a random affine map and apply it to the current point.
!
  do i = 1, n

    call random_number ( harvest = p )

    call r_to_bin_uneven ( nmap, map_bin, p, map )

    x(1:2) = matmul ( map_matrix(1:2,1:3,map), x(1:3) )

    coord(1:2,i) = x(1:2)

    write ( iunit, '(2g14.6)' ) x(1:2)

  end do

  return
end
function rvec_ascends ( n, x )
!
!*******************************************************************************
!
!! RVEC_ASCENDS determines if a real vector is (weakly) ascending.
!
!
!  Example:
!
!    X = ( -8.1, 1.3, 2.2, 3.4, 7.5, 7.5, 9.8 )
!
!    RVEC_ASCENDS = TRUE
!
!    The sequence is not required to be strictly ascending.
!
!  Modified:
!
!    07 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the size of the array.
!
!    Input, real X(N), the array to be examined.
!
!    Output, logical RVEC_ASCENDS, is TRUE if the entries of X ascend.
!
  implicit none
!
  integer n
!
  integer i
  logical rvec_ascends
  real x(n)
!
  rvec_ascends = .false.

  do i = 1, n-1
    if ( x(i) > x(i+1) ) then
      return
    end if
  end do

  rvec_ascends = .true.

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

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

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

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

  return
end
subroutine tree_map ( map_matrix )
!
!*******************************************************************************
!
!! TREE_MAP sets the tree map matrix.
!
!
!  Modified:
!
!    08 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real MAP_MATRIX(2,3,NMAP), the set of NMAP affine maps.
!    Each map has the form of a 2 by 3 matrix, so that
!      Xnew = A11 * X + A12 * Y + A13
!      Ynew = A21 * X + A22 * Y + A23
!
  implicit none
!
  real map_matrix(2,3,4)
!
  map_matrix(1:2,1:3,1) = reshape ( source = &
    (/ 0.000, 0.000, &
       0.000, 0.500, &
       0.500, 0.000 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,2) = reshape ( source = &
    (/ 0.100, 0.000, &
       0.000, 0.100, &
       0.450, 0.150 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,3) = reshape ( source = &
    (/  0.420, 0.420, &
       -0.420, 0.420, &
        0.290,-0.010 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,4) = reshape ( source = &
    (/  0.420, -0.420, &
        0.420,  0.420, &
        0.290,  0.410 /), shape = (/ 2, 3 /) )

  return
end
subroutine triangle_map ( map_matrix )
!
!*******************************************************************************
!
!! TRIANGLE_MAP sets the triangle map matrix.
!
!
!  Modified:
!
!    08 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real MAP_MATRIX(2,3,NMAP), the set of NMAP affine maps.
!    Each map has the form of a 2 by 3 matrix, so that
!      Xnew = A11 * X + A12 * Y + A13
!      Ynew = A21 * X + A22 * Y + A23
!
  implicit none
!
  real map_matrix(2,3,3)
!
  map_matrix(1:2,1:3,1) = reshape ( source = &
    (/ 0.500, 0.000, &
       0.000, 0.500, &
       0.000, 0.000 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,2) = reshape ( source = &
    (/ 0.500, 0.000, &
       0.000, 0.500, &
       0.500, 0.000 /), shape = (/ 2, 3 /) )

  map_matrix(1:2,1:3,3) = reshape ( source = &
    (/ 0.500, 0.000, &
       0.000, 0.500, &
       0.250, 0.500 /), shape = (/ 2, 3 /) )

  return
end