subroutine leaf_chaos ( n )

!*****************************************************************************80
!
!! leaf_chaos() draws a leaf using an iterated function system.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    23 August 2025
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Scott Bailey, Theodore Kim, Robert Strichartz,
!    Inside the Levy dragon,
!    American Mathematical Monthly,
!    Volume 109, Number 8, October 2002, pages 689-703.
!
!    Michael Barnsley, Alan Sloan,
!    A Better Way to Compress Images,
!    Byte Magazine,
!    Volume 13, Number 1, January 1988, pages 215-224.
!
!    Michael Barnsley,
!    Fractals Everywhere,
!    Academic Press, 1988,
!    ISBN: 0120790696,
!    LC: QA614.86.B37.
!
!    Michael Barnsley, Lyman Hurd,
!    Fractal Image Compression,
!    Peters, 1993,
!    ISBN: 1568810008,
!    LC: TA1632.B353
!
!    Alexander Dewdney,
!    Mathematical Recreations,
!    Scientific American,
!    Volume 262, Number 5, May 1990, pages 126-129.
!
!    Bernt Wahl, Peter VanRoy, Michael Larsen, Eric Kampman,
!    Exploring Fractals on the Mac,
!    Addison Wesley, 1995,
!    ISBN: 0201626306,
!    LC: QA614.86.W34.
!
!  Input:
!
!    integer n: the number of iterations.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ), save, dimension ( 2, 2 ) :: A0 = reshape ( (/ &
     0.800D+00,  0.000D+00, &
     0.000D+00,  0.800D+00 /), (/ 2, 2 /) )
  real ( kind = rk8 ), save, dimension ( 2, 2 ) :: A1 = reshape ( (/ &
     0.500D+00,  0.000D+00, &
     0.000D+00,  0.500D+00 /), (/ 2, 2 /) )
  real ( kind = rk8 ), save, dimension ( 2, 2 ) :: A2 = reshape ( (/ &
     0.355D+00,  0.355D+00, &
    -0.355D+00,  0.355D+00 /), (/ 2, 2 /) )
  real ( kind = rk8 ), save, dimension ( 2, 2 ) :: A3 = reshape ( (/ &
     0.355D+00, -0.355D+00, &
     0.355D+00,  0.355D+00 /), (/ 2, 2 /) )

  real ( kind = rk8 ), save, dimension ( 2, 1 ) :: b0 = reshape ( (/ &
    0.100D+00, 0.040D+00 /), (/ 2, 1 /) )
  real ( kind = rk8 ), save, dimension ( 2, 1 ) :: b1 = reshape ( (/ &
    0.250D+00, 0.400D+00 /), (/ 2, 1 /) )
  real ( kind = rk8 ), save, dimension ( 2, 1 ) :: b2 = reshape ( (/ &
    0.266D+00, 0.078D+00 /), (/ 2, 1 /) )
  real ( kind = rk8 ), save, dimension ( 2, 1 ) :: b3 = reshape ( (/ &
    0.378D+00, 0.434D+00 /), (/ 2, 1 /) )

  character ( len = 80 ) command_filename
  integer command_unit
  character ( len = 80 ) data_filename
  integer data_unit
  integer i
  integer, external :: i4_uniform_ab
  integer j
  integer n
  real ( kind = rk8 ), allocatable :: x(:,:)
!
!  Allocate space.
!
  allocate ( x(1:2,1:n+1) )
!
!  Random starting point in the unit square.
!
  call random_number ( harvest = x(1:2,1) )
!
!  Repeatedly compute and store A*x+b
!
  do i = 1, n
    j = i4_uniform_ab ( 0, 3 )
    if ( j == 0 ) then
      x(1:2,i+1) = matmul ( A0, x(1:2,i) ) + b0(1:2,1)
    elseif ( j == 1 ) then
      x(1:2,i+1) = matmul ( A1, x(1:2,i) ) + b1(1:2,1)
    elseif ( j == 2 ) then
      x(1:2,i+1) = matmul ( A2, x(1:2,i) ) + b2(1:2,1)
    elseif ( j == 3 ) then
      x(1:2,i+1) = matmul ( A3, x(1:2,i) ) + b3(1:2,1)
    end if
  end do
!
!  Create datafile.
!
  call get_unit ( data_unit )
  data_filename = 'leaf_chaos_data.txt'
  open ( unit = data_unit, file = data_filename, status = 'replace' )
  do i = 2, n + 1
    write ( data_unit, '(2x,g14.6,2x,g14.6)' ) x(1:2,i)
  end do
  close ( unit = data_unit )
  write ( *, '(a)' ) '  Created data file "' // trim ( data_filename ) // '".'
!
!  Create command file.
!
  call get_unit ( command_unit )
  command_filename = 'leaf_chaos_commands.txt'
  open ( unit = command_unit, file = command_filename, status = 'replace' )
  write ( command_unit, '(a)' ) '# ' // trim ( command_filename )
  write ( command_unit, '(a)' ) '#'
  write ( command_unit, '(a)' ) '# Usage:'
  write ( command_unit, '(a)' ) '#  gnuplot < ' // trim ( command_filename )
  write ( command_unit, '(a)' ) '#'
  write ( command_unit, '(a)' ) 'set term png'
  write ( command_unit, '(a)' ) 'set output "leaf_chaos.png"'
  write ( command_unit, '(a)' ) 'set title "leaf chaos"'
  write ( command_unit, '(a)' ) 'set grid'
  write ( command_unit, '(a)' ) 'set size square'
  write ( command_unit, '(a)' ) 'unset key'
  write ( command_unit, '(a)' ) 'set style data lines'
  write ( command_unit, '(a)' ) &
    'plot "leaf_chaos_data.txt" using 1:2 with points pt 7 ps 0.5'
  write ( command_unit, '(a)' ) 'quit'
  close ( unit = command_unit )
  write ( *, '(a)' ) &
    '  Created command file "' // trim ( command_filename ) // '".'

  return
end
subroutine get_unit ( iunit )

!*****************************************************************************80
!
!! get_unit() returns a free Fortran unit number.
!
!  Discussion:
!
!    A "free" Fortran unit number is a value between 1 and 99 which
!    is not currently associated with an I/O device.  A free Fortran unit
!    number is needed in order to open a file with the OPEN command.
!
!    If IUNIT = 0, then no free Fortran unit could be found, although
!    all 99 units were checked (except for units 5, 6 and 9, which
!    are commonly reserved for console I/O).
!
!    Otherwise, IUNIT is a value between 1 and 99, representing a
!    free Fortran unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    26 October 2008
!
!  Author:
!
!    John Burkardt
!
!  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_uniform_ab ( a, b )

!*****************************************************************************80
!
!! i4_uniform_ab() returns a scaled pseudorandom I4.
!
!  Discussion:
!
!    An I4 is an integer value.
!
!    The pseudorandom number will be scaled to be uniformly distributed
!    between A and B.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    27 August 2021
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer a, b: the limits of the interval.
!
!  Output:
!
!    integer i4_uniform_ab: a number between A and B.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer a
  integer b
  integer i4_uniform_ab
  real ( kind = rk8 ) r

  call random_number ( harvest = r )
  i4_uniform_ab = a + int ( ( b + 1 - a ) * r )

  return
end