program pbmlib_prb
!
!*******************************************************************************
!
!! PBMLIB_PRB calls the PBMLIB test routines.
!
  implicit none
!
  integer, parameter :: max_char = 100
  integer, parameter :: max_col = 5
  integer, parameter :: max_row = 7
!
  integer bits(max_row,max_col,max_char)
  integer ierror
  integer ipoint(0:255)
  integer nchar
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'PBMLIB_PRB'
  write ( *, '(a)' ) '  Tests the PBMLIB portable bit map routines.'

  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
  call test07
  call test08
  call test09
  call test10

  call test11
  call test12

  call test13 ( bits, ipoint, max_char, max_col, max_row, nchar, ierror )

  if ( ierror == 0 ) then
    call test14 ( bits, ipoint, max_char, max_col, max_row )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'Skipping TEST14 because of error in TEST13!'
  end if

  call test15 ( bits, ipoint, max_char, max_col, max_row, nchar )

  call test16

  call test17

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

  stop
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 tests PBMA_WRITE.
!
  implicit none
!
  integer, parameter :: ncol = 200
  integer, parameter :: nrow = 200
!
  integer b(nrow,ncol)
  character ( len = 80 ) filename
  integer i
  integer ierror
  integer j
  real r
  real test
  real x
  real xc
  real y
  real yc
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  PBMA_WRITE writes an ASCII PBM file.'

  filename = 'test01.pbm'

  xc = real ( nrow ) / 2.0E+00
  yc = real ( ncol ) / 2.0E+00
  r = real ( min ( nrow, ncol ) ) / 3.0E+00

  do i = 1, nrow
    do j = 1, ncol

      x = real ( i )
      y = real ( j ) 

      test = r - sqrt ( ( x - xc )**2 + ( y - yc )**2 )

      if ( abs ( test ) <= 3.0E+00 ) then
        b(i,j) = 1
      else
        b(i,j) = 0
      end if

    end do
  end do

  call pbma_write ( filename, ierror, nrow, ncol, b )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PBMA_WRITE returns IERROR = ', ierror
  end if

  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 tests PGMA_WRITE.
!
  implicit none
!
  integer, parameter :: ncol = 200
  integer, parameter :: nrow = 200
!
  character ( len = 80 ) filename
  integer g(nrow,ncol)
  integer i
  integer ierror
  integer j
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  PGMA_WRITE writes an ASCII PGM file.'

  filename = 'test02.pgm'

  do i = 1, nrow
    do j = 1, ncol

      g(i,j) = mod ( i, j )

    end do
  end do

  call pgma_write ( filename, ierror, nrow, ncol, g )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PGMA_WRITE returns IERROR = ', ierror
  end if

  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 tests PPMA_WRITE.
!
  implicit none
!
  integer, parameter :: ncol = 20
  integer, parameter :: nrow = 10
!
  integer b(nrow,ncol)
  character ( len = 80 ) filename
  integer g(nrow,ncol)
  integer i
  integer ierror
  integer j
  integer r(nrow,ncol)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  PPMA_WRITE writes an ASCII PPM file.'

  filename = 'test03.ppm'

  do i = 1, nrow
    do j = 1, ncol

      if ( i >= j ) then

        r(i,j) = 255
        g(i,j) = 0
        b(i,j) = 0

      else if ( j <= 10 ) then

        r(i,j) = 0
        g(i,j) = 255
        b(i,j) = 255

      else 

        r(i,j) = 0

        g(i,j) = ( i - 1 ) * 255 * ( ncol - j ) &
          / ( ( ncol - 11 ) * ( nrow - 1 ) )

        b(i,j) = ( i - 1 ) * 255 * ( j - 11 ) &
          / ( ( ncol - 11 ) * ( nrow - 1 ) )

      end if

    end do
  end do

  call ppma_write ( filename, ierror, nrow, ncol, r, g, b )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PPMA_WRITE returns IERROR = ', ierror
  end if

  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 tests PBMB_WRITE.
!
  implicit none
!
  integer, parameter :: ncol = 200
  integer, parameter :: nrow = 200
!
  integer b(nrow,ncol)
  character ( len = 80 ) filename
  integer i
  integer ierror
  integer j
  real r
  real test
  real x
  real xc
  real y
  real yc
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  PBMB_WRITE writes a binary PBM file.'

  filename = 'test04.pbm'

  xc = real ( nrow ) / 2.0E+00
  yc = real ( ncol ) / 2.0E+00
  r = real ( min ( nrow, ncol ) ) / 3.0E+00

  do i = 1, nrow
    do j = 1, ncol

      x = real ( i )
      y = real ( j ) 

      test = r - sqrt ( ( x - xc )**2 + ( y - yc )**2 )

      if ( abs ( test ) <= 3.0E+00 ) then
        b(i,j) = 1
      else
        b(i,j) = 0
      end if

    end do
  end do

  call pbmb_write ( filename, ierror, nrow, ncol, b )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PBMB_WRITE returns IERROR = ', ierror
  end if

  return
end
subroutine test05
!
!*******************************************************************************
!
!! TEST05 tests PGMB_WRITE.
!
  implicit none
!
  integer, parameter :: ncol = 200
  integer, parameter :: nrow = 200
!
  character ( len = 80 ) filename
  integer g(nrow,ncol)
  integer i
  integer ierror
  integer j
  integer maxg
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  PGMB_WRITE writes a binary PGM file.'

  filename = 'test05.pgm'

  maxg = 0

  do i = 1, nrow
    do j = 1, ncol

      g(i,j) = mod ( i, j )
      maxg = max ( maxg, g(i,j) )

    end do
  end do

  do i = 1, nrow
    do j = 1, ncol

      g(i,j) = nint ( real ( 255 * g(i,j) ) / real ( maxg ) )

    end do
  end do

  call pgmb_write ( filename, ierror, nrow, ncol, g )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PGMB_WRITE returns IERROR = ', ierror
  end if

  return
end
subroutine test06
!
!*******************************************************************************
!
!! TEST06 tests PPMB_WRITE.
!
  implicit none
!
  integer, parameter :: ncol = 20
  integer, parameter :: nrow = 10
!
  integer b(nrow,ncol)
  character ( len = 80 ) filename
  integer g(nrow,ncol)
  integer i
  integer ierror
  integer j
  integer r(nrow,ncol)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  PPMB_WRITE writes a binary PPM file.'

  filename = 'test06.ppm'

  do i = 1, nrow
    do j = 1, ncol

      if ( i >= j ) then

        r(i,j) = 255
        g(i,j) = 0
        b(i,j) = 0

      else if ( j <= 10 ) then

        r(i,j) = 0
        g(i,j) = 255
        b(i,j) = 255

      else 

        r(i,j) = 0

        g(i,j) = ( i - 1 ) * 255 * ( ncol - j ) &
          / ( ( ncol - 11 ) * ( nrow - 1 ) )

        b(i,j) = ( i - 1 ) * 255 * ( j - 11 ) &
          / ( ( ncol - 11 ) * ( nrow - 1 ) )

      end if

    end do
  end do

  call ppmb_write ( filename, ierror, nrow, ncol, r, g, b )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PPMB_WRITE returns IERROR = ', ierror
  end if

  return
end
subroutine test07
!
!*******************************************************************************
!
!! TEST07 tests PBMA_READ.
!
  implicit none
!
  integer, parameter :: max_b = 200*200
!
  integer b(max_b)
  character ( len = 80 ) filename
  integer ierror
  integer ncol
  integer nrow
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  PBMA_READ reads an ASCII PBM file.'

  filename = 'test01.pbm'

  call pbma_read ( filename, ierror, max_b, nrow, ncol, b  )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PBMA_READ returns IERROR = ', ierror
  end if

  return
end
subroutine test08
!
!*******************************************************************************
!
!! TEST08 tests PGMA_READ.
!
  implicit none
!
  integer, parameter :: max_g = 200*200
!
  character ( len = 80 ) filename
  integer g(max_g)
  integer ierror
  integer maxcol
  integer ncol
  integer nrow
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  PGMA_READ reads an ASCII PGM file.'

  filename = 'test02.pgm'

  call pgma_read ( filename, ierror, maxcol, max_g, nrow, ncol, g  )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PGMA_READ returns IERROR = ', ierror
  end if

  return
end
subroutine test09
!
!*******************************************************************************
!
!! TEST09 tests PPMA_READ.
!
  implicit none
!
  integer, parameter :: max_p = 20*10
!
  integer b(max_p)
  character ( len = 80 ) filename
  integer g(max_p)
  integer ierror
  integer maxcol
  integer ncol
  integer nrow
  integer r(max_p)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  PPMA_READ reads an ASCII PPM file.'

  filename = 'test03.ppm'

  call ppma_read ( filename, ierror, maxcol, max_p, nrow, ncol, r, g, b  )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PPMA_READ returns IERROR = ', ierror
  end if

  return
end
subroutine test10
!
!*******************************************************************************
!
!! TEST10 tests PBMB_READ.
!
  implicit none
!
  integer, parameter :: max_g = 200*200
!
  integer b(max_g)
  character ( len = 80 ) filename
  integer ierror
  integer ncol
  integer nrow
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  PBMB_READ reads a binary PBM file.'

  filename = 'test04.pbm'

  call pbmb_read ( filename, ierror, max_g, nrow, ncol, b )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PBMB_READ returns IERROR = ', ierror
  end if

  return
end
subroutine test11
!
!*******************************************************************************
!
!! TEST11 tests PGMB_READ.
!
  implicit none
!
  integer, parameter :: max_g = 200*200
!
  character ( len = 80 ) filename
  integer g(max_g)
  integer ierror
  integer maxcol
  integer ncol
  integer nrow
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  PGMB_READ reads a binary PGM file.'

  filename = 'test05.pgm'

  call pgmb_read ( filename, ierror, maxcol, max_g, nrow, ncol, g )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PGMB_READ returns IERROR = ', ierror
  end if

  return
end
subroutine test12
!
!*******************************************************************************
!
!! TEST12 tests PPMB_READ.
!
  implicit none
!
  integer, parameter :: max_g = 20*10
!
  integer b(max_g)
  character ( len = 80 ) filename
  integer g(max_g)
  integer ierror
  integer maxcol
  integer ncol
  integer nrow
  integer r(max_g)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  PPMB_READ reads a binary PPM file.'

  filename = 'test06.ppm'

  call ppmb_read ( filename, ierror, maxcol, max_g, nrow, ncol, r, g, b )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PPMB_READ returns IERROR = ', ierror
  end if

  return
end
subroutine test13 ( bits, ipoint, maxchar, maxcol, maxrow, nchar, ierror )
!
!*******************************************************************************
!
!! TEST13 tests FONT_READ.
!
  implicit none
!
  integer maxchar
  integer maxcol
  integer maxrow
!
  integer bits(maxrow,maxcol,maxchar)
  integer ierror
  integer inunit
  integer ios
  integer ipoint(0:255)
  integer nchar
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST13'
  write ( *, '(a)' ) '  FONT_READ reads in a simple bit map font.'

  ierror = 0
  inunit = 1

  open ( unit = inunit, file = 'alphabits.txt', status = 'old', &
    iostat = ios )

  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TEST13'
    write ( *, '(a)' ) '  Could not open the font file!'
    return
  end if

  call font_read ( bits, ierror, inunit, ipoint, maxchar, maxcol, maxrow, &
    nchar )

  close ( unit = inunit )

  return
end
subroutine test14 ( bits, ipoint, maxchar, maxcol, maxrow )
!
!*******************************************************************************
!
!! TEST14 tests FONT_PRINT.
!
  implicit none
!
  integer maxchar
  integer maxcol
  integer maxrow
!
  integer bits(maxrow,maxcol,maxchar)
  integer ipoint(0:255)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST14'
  write ( *, '(a)' ) '  FONT_PRINT prints the font we just read in.'

  call font_print ( bits, ipoint, maxchar, maxcol, maxrow )

  return
end
subroutine test15 ( bits, ipoint, maxchar, maxcol, maxrow, nchar )
!
!*******************************************************************************
!
!! TEST15 tests FONT_DATA.
!
  implicit none
!
  integer maxchar
  integer maxcol
  integer maxrow
!
  integer bits(maxrow,maxcol,maxchar)
  integer ipoint(0:255)
  integer nchar
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST15'
  write ( *, '(a)' ) '  FONT_DATA prints out a font data statement.'

  call font_data ( bits, ipoint, maxchar, maxcol, maxrow, nchar )

  return
end
subroutine test16
!
!*******************************************************************************
!
!! TEST16 tests BITCHR75.
!
  implicit none
!
  character, parameter, dimension ( 0:1 ) :: bit = (/ ' ', '*' /)
  character c
  integer i
  integer j
  integer k
  integer pattern(7,5)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST16'
  write ( *, '(a)' ) '  BITCHR75 returns a 7 x 5 representation of'
  write ( *, '(a)' ) '    SOME characters.'

  do k = 32, 96

    c = char ( k )

    call bitchr75 ( c, pattern )

    write ( *, '(a)' ) ' '
    write ( *, '(a,i6,a)' ) 'Character ', k, ' is "' // c // '"'
    write ( *, '(a)' ) ' '

    do i = 1, 7
      write ( *, '(5a1)' ) ( bit ( pattern(i,j) ), j = 1, 5 )
    end do

  end do

  return
end
subroutine test17
!
!*******************************************************************************
!
!! TEST17 tests PPMA_WRITE.
!
!
!  Discussion:
!
!    This example uses color to record which starting values converge to
!    which roots in a Newton iteration in the complex plane.
!
!  Modified:
!
!    22 April 2002
!
  implicit none
!
  integer, parameter :: ncol = 1152
  integer, parameter :: nrow = 864
!
  integer b(nrow,ncol)
  character ( len = 80 ) filename
  integer g(nrow,ncol)
  integer i
  integer ierror
  integer it
  integer, parameter :: it_max = 20
  integer j
  integer r(nrow,ncol)
  real r_pi
  real radius
  complex root1
  complex root2
  complex root3
  real theta
  real, parameter :: tol = 0.01E+00
  real x
  real, parameter :: xmax = 1.4E+00
  real, parameter :: xmin = -1.4E+00
  real y
  real, parameter :: ymax = 1.2E+00
  real, parameter :: ymin = -1.2E+00
  complex z
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST17'
  write ( *, '(a)' ) '  PPMA_WRITE writes an ASCII PPM file.'
  write ( *, '(a)' ) '  This displays the basins of attraction for'
  write ( *, '(a)' ) '  Newton''s method applied to a particular'
  write ( *, '(a)' ) '  nonlinear equation in the complex plane.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The equation is Z^3 = -1, which has three roots.'
  write ( *, '(a)' ) ' '

  radius = 1.0E+00
  theta = r_pi ( ) / 3.0E+00
  call rt_to_xy ( radius, theta, x, y )
  root1 = cmplx ( x, y )
  write ( *, '(a,2g14.6)' ) '  ROOT1 = ', root1

  theta = r_pi ( )
  call rt_to_xy ( radius, theta, x, y )
  root2 = cmplx ( x, y )
  write ( *, '(a,2g14.6)' ) '  ROOT2 = ', root2

  theta = 5.0E+00 * r_pi ( ) / 3.0E+00
  call rt_to_xy ( radius, theta, x, y )
  root3 = cmplx ( x, y )
  write ( *, '(a,2g14.6)' ) '  ROOT3 = ', root3

  filename = 'test17.ppm'

  do i = 1, nrow

    y = ( ymax * real ( nrow - i ) + ymin * real ( i - 1 ) ) &
      / real ( nrow - 1 )

    do j = 1, ncol

      x = ( xmin * real ( ncol - j ) + xmax * real ( j - 1 ) ) &
        / real ( ncol - 1 )

      z = cmplx ( x, y )

      r(i,j) = 100
      g(i,j) = 100
      b(i,j) = 100

      do it = 1, it_max
!
!  This is one step of Newton's iteration for f(z) = z^3 + 1.
!
        z = z - ( z * z * z + 1.0E+00 ) / ( 3.0E+00 * z * z )
!
!  This is one step of Sterlings iteration for f(z) = z^3 + 1.
!  We have commented this one out.
!
!       z = z - ( z * z * z + 1.0E+00 ) / ( 3.0E+00 * ( z - z**3 - 1.0E+00 )**2 )
!
!  It the current iterate is close enough to a root, exit.
!  But mark the starting point with a shade of red, green or blue
!  depending on which root it converged to.
!
!  To make visualization easier, we gray out the red, green or blue
!  depending on the number of iterations taken.
!
!  Also, every third iteration we set the color to white, to make the
!  bands easier to see.
!
        if ( abs ( z - root1 ) < tol ) then

          if ( mod ( it, 3 ) /= 2 ) then
            r(i,j) = ( ( it_max + 1 - it ) * 255 + ( it - 1 ) * 127 ) / it_max
            g(i,j) = ( ( it_max + 1 - it ) *   0 + ( it - 1 ) * 127 ) / it_max
            b(i,j) = ( ( it_max + 1 - it ) *   0 + ( it - 1 ) * 127 ) / it_max
          else
            r(i,j) = 255
            g(i,j) = 255
            b(i,j) = 255
          end if

          exit

        else if ( abs ( z - root2 ) < tol ) then

          if ( mod ( it, 3 ) /= 2 ) then
            r(i,j) = ( ( it_max + 1 - it ) *   0 + ( it - 1 ) * 127 ) / it_max
            g(i,j) = ( ( it_max + 1 - it ) * 255 + ( it - 1 ) * 127 ) / it_max
            b(i,j) = ( ( it_max + 1 - it ) *   0 + ( it - 1 ) * 127 ) / it_max
          else
            r(i,j) = 255
            g(i,j) = 255
            b(i,j) = 255
          end if

          exit

        else if ( abs ( z - root3 ) < tol ) then

          if ( mod ( it, 3 ) /= 2 ) then
            r(i,j) = ( ( it_max + 1 - it ) *   0 + ( it - 1 ) * 127 ) / it_max
            g(i,j) = ( ( it_max + 1 - it ) *   0 + ( it - 1 ) * 127 ) / it_max
            b(i,j) = ( ( it_max + 1 - it ) * 255 + ( it - 1 ) * 127 ) / it_max
          else
            r(i,j) = 255
            g(i,j) = 255
            b(i,j) = 255
          end if

          exit

        end if

      end do

    end do

  end do
!
!  Now fill with gray the three roots.
!
  do i = 1, nrow

    y = ( ymax * real ( nrow - i ) + ymin * real ( i - 1 ) ) &
      / real ( nrow - 1 )

    do j = 1, ncol

      x = ( xmin * real ( ncol - j ) + xmax * real ( j - 1 ) ) &
        / real ( ncol - 1 )

      z = cmplx ( x, y )

      if ( abs ( z - root1 ) < 4.0E+00 * tol ) then

        r(i,j) = 100
        g(i,j) = 100
        b(i,j) = 100

      else if ( abs ( z - root2 ) < 4.0E+00 * tol ) then

        r(i,j) = 100
        g(i,j) = 100
        b(i,j) = 100

      else if ( abs ( z - root3 ) < 4.0E+00 * tol ) then

        r(i,j) = 150
        g(i,j) = 150
        b(i,j) = 150

      end if

    end do

  end do
!
!  Now write all this data to a file.
!
  call ppma_write ( filename, ierror, nrow, ncol, r, g, b )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) 'PPMA_WRITE returns IERROR = ', ierror
  end if

  return
end