program frieze
!
!*******************************************************************************
!
!! FRIEZE demonstrates the use of blending to tile a region with a pattern.
!
!
!  Discussion:
!
!    We used this program to see what blending would look like.  We were
!    really interested in 3D models, but wanted to see what was involved
!    with a simple 2D problem first.
!
!    This program requires the PS_WRITE routines in order to create a
!    PostScript file with an image of the tiled region.
!
!  Reference:
!
!    W N Gordon and Charles A Hall,
!    Construction of Curvilinear Coordinate Systems and Application to
!      Mesh Generation,
!    International Journal of Numerical Methods in Engineering,
!    Volume 7, pages 461-477, 1973.
!
!    Joe Thompson, Bharat Soni, Nigel Weatherill,
!    Handbook of Grid Generation,
!    CRC Press,
!    1999.
!
!  Modified:
!
!    28 September 1999
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FRIEZE'
  write ( *, '(a)' ) '  Demonstrate the use of blending to "tile"'
  write ( *, '(a)' ) '  a region with a pattern.'
!
!  Draw the pattern.
!
  call draw_pattern
!
!  Draw the region
!
  call draw_region
!
!  Draw the region subdivided into cells
!
  call draw_cells
!
!  Draw the frieze.
!
  call draw_frieze

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FRIEZE'
  write ( *, '(a)' ) '  Normal end of the demonstration.'

  stop
end
subroutine draw_frieze
!
!*******************************************************************************
!
!! DRAW_FRIEZE draws the frieze.
!
!
!  Modified:
!
!    28 September 1999
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer, parameter :: max_end = 100
  integer, parameter :: max_pat = 100
  integer, parameter :: n = 4
  integer, parameter :: nfine = 21
!
  character ( len = 80 ) filename
  integer i
  integer icol
  integer iend(max_end)
  integer ierror
  integer iunit
  integer irow
  integer ncol
  integer num_end
  integer num_pat
  integer nrow
  real r
  real re
  real rw
  real s
  real sn
  real ss
  real x(nfine)
  real xcval(4)
  real xmax
  real xmin
  real xne
  real xnw
  real xpat(max_pat)
  real xse
  real xsw
  real y(nfine)
  real ycval(n)
  real ymax
  real ymin
  real yne
  real ynw
  real ypat(max_pat)
  real yse
  real ysw
!
!  Get the pattern to be used as a tile.
!
  call pattern ( iend, max_end, num_end, max_pat, num_pat, xpat, ypat )
!
!  Open the PostScript file.
!
  iunit = 1
  filename = 'frieze.ps'

  call ps_file_open ( filename, iunit, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FRIEZE'
    write ( *, '(a,i6)' ) '  PostScript file creation error ', ierror
    stop
  end if

  xmin = -3.0 E+00
  ymin =  2.0E+00
  xmax = +3.0E+00
  ymax =  4.5E+00

  call ps_file_head ( filename, xmax, xmin, ymax, ymin )

  call ps_page_head
!
!  Query the BOUNDARY routine for the (X,Y) locations of the corners.
!
  call boundary_2d ( 1.0E+00, 0.0E+00, xse, yse )
  call boundary_2d ( 1.0E+00, 1.0E+00, xne, yne )
  call boundary_2d ( 0.0E+00, 1.0E+00, xnw, ynw )
  call boundary_2d ( 0.0E+00, 0.0E+00, xsw, ysw )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '(X,Y) corners of the total region:'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) 'XSE, YSE = ', xse, yse
  write ( *, '(a,2g14.6)' ) 'XNE, YNE = ', xne, yne
  write ( *, '(a,2g14.6)' ) 'XNW, YNW = ', xnw, ynw
  write ( *, '(a,2g14.6)' ) 'XSW, YSW = ', xsw, ysw
!
!  Draw the boundary of the region.
!
  call ps_line_rgb ( 0.7E+00, 0.7E+00, 0.7E+00 )

  r = 1.0E+00
  do i = 1, nfine
    s = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  s = 1.0E+00
  do i = nfine, 1, -1
    r = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  r = 0.0E+00
  do i = nfine, 1, -1
    s = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  s = 0.0E+00
  do i = 1, nfine
    r = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )
!
!  Map the reference element into subregion IROW, JCOL.
!
  ncol = 6
  nrow = 2

  do irow = 1, nrow

    sn = real ( irow ) / real ( nrow )
    ss = real ( irow - 1 ) / real ( nrow )

    do icol = 1, ncol

      re = real ( icol ) / real ( ncol )
      rw = real ( icol - 1 ) / real ( ncol )

      call blend1d2d ( re, ss, xcval(1), ycval(1) )
      call blend1d2d ( re, sn, xcval(2), ycval(2) )
      call blend1d2d ( rw, sn, xcval(3), ycval(3) )
      call blend1d2d ( rw, ss, xcval(4), ycval(4) )

      do i = 1, 4
        write ( *, * ) xcval(i), ycval(i)
      end do

      call ps_line_rgb ( 0.0E+00, 0.0E+00, 0.4E+00 )

      call drawcell ( re, rw, sn, ss )

      call ps_line_rgb ( 1.0E+00, 0.0E+00, 0.0E+00 )

      call drawpat ( iend, num_end, num_pat, xpat, ypat, re, rw, sn, ss )

    end do
  end do
!
!  Finish up the PostScript file.
!
  call ps_page_tail

  call ps_file_tail

  call ps_file_close ( iunit )

  return
end
subroutine drawpat ( iend, num_end, num_pat, xpat, ypat, re, rw, sn, ss )
!
!*******************************************************************************
!
!! DRAWPAT draws a copy of the pattern.
!
!
!  Reference:
!
!    W N Gordon and Charles A Hall,
!    Construction of Curvilinear Coordinate Systems and Application to
!      Mesh Generation,
!    International Journal of Numerical Methods in Engineering,
!    Volume 7, pages 461-477, 1973.
!
!    Joe Thompson, Bharat Soni, Nigel Weatherill,
!    Handbook of Grid Generation,
!    CRC Press,
!    1999.
!
!  Modified:
!
!    02 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IEND(NUM_END), contains the indices of XPAT and YPAT
!    that represent ends of line segments.  The first line segment involves
!    the entries IEND(1) + 1 through IEND(2) of XPAT and YPAT, and the last
!    one involves IEND(NUM_END-1) + 1 through IEND(NUM_END).
!
!    Input, integer NUM_END, the number of entries used in IEND, and one
!    more than the number of line segments used to draw the pattern.
!
!    Input, integer NUM_PAT, the number of entries used in XPAT and YPAT.
!
!    Input, real XPAT(NUM_PAT), YPAT(MAX_PAT), the X and Y coordinates of
!    points that define the line segments used to draw the pattern.
!
!    Input, real RE, RW, SN, SS, the extreme (east and west) R coordinates
!    and extreme (north and south) S coordinates of the subcell.
!
  implicit none
!
  integer, parameter :: MAX_POINT = 20
  integer num_end
  integer num_pat
!
  integer i
  integer iend(num_end)
  integer ipat
  integer j
  integer npoint
  real r
  real re
  real rw
  real s
  real sn
  real ss
  real xpoint(MAX_POINT)
  real xpat(num_pat)
  real ypoint(MAX_POINT)
  real ypat(num_pat)

  ipat = 0
  do i = 1, num_end - 1

    npoint = 0
    do j = iend(i)+1, iend(i+1)
      npoint = npoint + 1
      ipat = ipat + 1
      r = ( xpat(ipat) * re + ( 1.0E+00 - xpat(ipat) ) * rw )
      s = ( ypat(ipat) * sn + ( 1.0E+00 - ypat(ipat) ) * ss )
      call blend1d2d ( r, s, xpoint(npoint), ypoint(npoint) )
    end do

    call ps_line ( npoint, xpoint, ypoint )

  end do

  return
end
subroutine draw_pattern
!
!*******************************************************************************
!
!! DRAW_PATTERN draws the pattern that will be repeated in the frieze.
!
!
!  Modified:
!
!    28 September 1999
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer, parameter :: max_end = 100
  integer, parameter :: max_pat = 100
  integer, parameter :: MAX_POINT = 20
  integer, parameter :: nfine = 21
!
  character ( len = 80 ) filename
  integer i
  integer iend(max_end)
  integer ierror
  integer iunit
  integer ipat
  integer j
  integer npoint
  integer num_end
  integer num_pat
  real x(nfine)
  real xmax
  real xmin
  real xpat(max_pat)
  real xpoint(MAX_POINT)
  real y(nfine)
  real ymax
  real ymin
  real ypat(max_pat)
  real ypoint(MAX_POINT)
!
!  Get the pattern to be used as a tile.
!
  call pattern ( iend, max_end, num_end, max_pat, num_pat, xpat, ypat )
!
!  Open the PostScript file.
!
  iunit = 1
  filename = 'pattern.ps'

  call ps_file_open ( filename, iunit, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FRIEZE'
    write ( *, '(a,i6)' ) '  PostScript file creation error ', ierror
    stop
  end if

  xmin =  0.0E+00
  ymin =  0.0E+00
  xmax =  1.0E+00
  ymax =  1.0E+00

  call ps_file_head ( filename, xmax, xmin, ymax, ymin )

  call ps_page_head
!
!  Draw the boundary of the region.
!
  call ps_line_rgb ( 0.7E+00, 0.7E+00, 0.7E+00 )

  x(1) = xmin
  y(1) = ymin
  x(2) = xmax
  y(2) = ymin
  x(3) = xmax
  y(3) = ymax
  x(4) = xmin
  y(4) = ymax
  x(5) = xmin
  y(5) = ymin

  call ps_line ( 5, x, y )

  call ps_line_rgb ( 1.0E+00, 0.0E+00, 0.0E+00 )

  ipat = 0
  do i = 1, num_end - 1

    npoint = 0
    do j = iend(i)+1, iend(i+1)
      npoint = npoint + 1
      ipat = ipat + 1
      xpoint(npoint) = xpat(ipat)
      ypoint(npoint) = ypat(ipat)
    end do

    call ps_line ( npoint, xpoint, ypoint )

  end do
!
!  Finish up the PostScript file.
!
  call ps_page_tail

  call ps_file_tail

  call ps_file_close ( iunit )

  return
end
subroutine drawcell ( re, rw, sn, ss )
!
!*******************************************************************************
!
!! DRAWCELL draws the borders of a given cell.
!
!
!  Discussion:
!
!    The internal parameter NPOINT determines how many points are
!    used to draw each of the four lines that define the border.
!    Increase it to draw a wiggly boundary more accurately.
!
!  Reference:
!
!    W N Gordon and Charles A Hall,
!    Construction of Curvilinear Coordinate Systems and Application to
!      Mesh Generation,
!    International Journal of Numerical Methods in Engineering,
!    Volume 7, pages 461-477, 1973.
!
!    Joe Thompson, Bharat Soni, Nigel Weatherill,
!    Handbook of Grid Generation,
!    CRC Press,
!    1999.
!
!  Modified:
!
!    02 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real RE, RW, SN, SS, the extreme (east and west) R coordinates
!    and extreme (north and south) S coordinates of the subcell.
!
  implicit none
!
  integer, parameter :: npoint = 5
!
  integer i
  real r
  real re
  real rw
  real s
  real sn
  real ss
  real xpoint(npoint)
  real ypoint(npoint)
!
!  East cell border.
!
  r = re
  do i = 1, npoint
    s = ( real ( npoint - i ) * ss + real ( i - 1 ) * sn ) / real ( npoint - 1 )
    call blend1d2d ( r, s, xpoint(i), ypoint(i) )
  end do

  call ps_line ( npoint, xpoint, ypoint )
!
!  North cell border.
!
  s = sn
  do i = 1, npoint
    r = ( real ( npoint - i ) * re + real ( i - 1 ) * rw ) / real ( npoint - 1 )
    call blend1d2d ( r, s, xpoint(i), ypoint(i) )
  end do

  call ps_line ( npoint, xpoint, ypoint )
!
!  West cell border.
!
  r = rw
  do i = 1, npoint
    s = ( real ( npoint - i ) * sn + real ( i - 1 ) * ss ) / real ( npoint - 1 )
    call blend1d2d ( r, s, xpoint(i), ypoint(i) )
  end do

  call ps_line ( npoint, xpoint, ypoint )
!
!  South cell border.
!
  s = ss
  do i = 1, npoint
    r = ( real ( npoint - i ) * rw + real ( i - 1 ) * re ) / real ( npoint - 1 )
    call blend1d2d ( r, s, xpoint(i), ypoint(i) )
  end do

  call ps_line ( npoint, xpoint, ypoint )

  return
end
subroutine drawdiag 
!
!*******************************************************************************
!
!! DRAWDIAG draws the diagonals of the entire region.
!
!
!  Reference:
!
!    W N Gordon and Charles A Hall,
!    Construction of Curvilinear Coordinate Systems and Application to
!      Mesh Generation,
!    International Journal of Numerical Methods in Engineering,
!    Volume 7, pages 461-477, 1973.
!
!    Joe Thompson, Bharat Soni, Nigel Weatherill,
!    Handbook of Grid Generation,
!    CRC Press,
!    1999.
!
!  Modified:
!
!    02 July 1999
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer, parameter :: nfine = 21
!
  integer i
  real r
  real re
  real rw
  real s
  real sn
  real ss
  real x(nfine)
  real y(nfine)
!
  re = 1.0E+00
  rw = 0.0E+00
  sn = 1.0E+00
  ss = 0.0E+00

  do i = 1, nfine
    r = ( real ( nfine - i ) * rw + real ( i - 1 ) * re ) / real ( nfine - 1 )
    s = ( real ( nfine - i ) * sn + real ( i - 1 ) * ss ) / real ( nfine - 1 )
    call blend1d2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  do i = 1, nfine
    r = ( real ( nfine - i ) * rw + real ( i - 1 ) * re ) / real ( nfine - 1 )
    s = ( real ( nfine - i ) * ss + real ( i - 1 ) * sn ) / real ( nfine - 1 )
    call blend1d2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  return
end
subroutine draw_cells
!
!*******************************************************************************
!
!! DRAW_CELLS makes an image of the region subdivided into cells.
!
!
!  Modified:
!
!    28 September 1999
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer, parameter :: n = 4
  integer, parameter :: nfine = 21
!
  character ( len = 80 ) filename
  integer i
  integer icol
  integer ierror
  integer iunit
  integer irow
  integer ncol
  integer nrow
  real r
  real re
  real rw
  real s
  real sn
  real ss
  real x(nfine)
  real xcval(4)
  real xmax
  real xmin
  real xne
  real xnw
  real xse
  real xsw
  real y(nfine)
  real ycval(n)
  real ymax
  real ymin
  real yne
  real ynw
  real yse
  real ysw
!
!  Open the PostScript file.
!
  iunit = 1
  filename = 'cells.ps'

  call ps_file_open ( filename, iunit, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FRIEZE'
    write ( *, '(a,i6)' ) '  PostScript file creation error ', ierror
    stop
  end if

  xmin = -3.0E+00
  ymin =  2.0E+00
  xmax = +3.0E+00
  ymax =  4.5E+00

  call ps_file_head ( filename, xmax, xmin, ymax, ymin )

  call ps_page_head
!
!  Query the BOUNDARY routine for the (X,Y) locations of the corners.
!
  call boundary_2d ( 1.0E+00, 0.0E+00, xse, yse )
  call boundary_2d ( 1.0E+00, 1.0E+00, xne, yne )
  call boundary_2d ( 0.0E+00, 1.0E+00, xnw, ynw )
  call boundary_2d ( 0.0E+00, 0.0E+00, xsw, ysw )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '(X,Y) corners of the total region:'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) 'XSE, YSE = ', xse, yse
  write ( *, '(a,2g14.6)' ) 'XNE, YNE = ', xne, yne
  write ( *, '(a,2g14.6)' ) 'XNW, YNW = ', xnw, ynw
  write ( *, '(a,2g14.6)' ) 'XSW, YSW = ', xsw, ysw
!
!  Draw the boundary of the region.
!
  call ps_line_rgb ( 0.7E+00, 0.7E+00, 0.7E+00 )

  r = 1.0E+00
  do i = 1, nfine
    s = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  s = 1.0E+00
  do i = nfine, 1, -1
    r = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  r = 0.0E+00
  do i = nfine, 1, -1
    s = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  s = 0.0E+00
  do i = 1, nfine
    r = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )
!
!  Map the reference element into subregion IROW, JCOL.
!
  ncol = 6
  nrow = 2

  do irow = 1, nrow

    sn = real ( irow ) / real ( nrow )
    ss = real ( irow - 1 ) / real ( nrow )

    do icol = 1, ncol

      re = real ( icol ) / real ( ncol )
      rw = real ( icol - 1 ) / real ( ncol )

      call blend1d2d ( re, ss, xcval(1), ycval(1) )
      call blend1d2d ( re, sn, xcval(2), ycval(2) )
      call blend1d2d ( rw, sn, xcval(3), ycval(3) )
      call blend1d2d ( rw, ss, xcval(4), ycval(4) )

      do i = 1, 4
        write ( *, * ) xcval(i), ycval(i)
      end do

      call ps_line_rgb ( 0.0E+00, 0.0E+00, 0.4E+00 )

      call drawcell ( re, rw, sn, ss )

    end do
  end do
!
!  Finish up the PostScript file.
!
  call ps_page_tail

  call ps_file_tail

  call ps_file_close ( iunit )

  return
end
subroutine draw_region
!
!*******************************************************************************
!
!! DRAW_REGION makes an image of the region.
!
!
!  Modified:
!
!    29 September 1999
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer, parameter :: n = 4
  integer, parameter :: nfine = 21
!
  character ( len = 80 ) filename
  integer i
  integer icol
  integer ierror
  integer iunit
  integer irow
  integer ncol
  integer nrow
  real r
  real re
  real rw
  real s
  real sn
  real ss
  real x(nfine)
  real xcval(4)
  real xmax
  real xmin
  real xne
  real xnw
  real xse
  real xsw
  real y(nfine)
  real ycval(n)
  real ymax
  real ymin
  real yne
  real ynw
  real yse
  real ysw
!
!  Open the PostScript file.
!
  iunit = 1
  filename = 'region.ps'

  call ps_file_open ( filename, iunit, ierror )

  if ( ierror /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FRIEZE'
    write ( *, '(a,i6)' ) '  PostScript file creation error ', ierror
    stop
  end if

  xmin = -3.0E+00
  ymin =  2.0E+00
  xmax = +3.0E+00
  ymax =  4.5E+00

  call ps_file_head ( filename, xmax, xmin, ymax, ymin )

  call ps_page_head
!
!  Query the BOUNDARY routine for the (X,Y) locations of the corners.
!
  call boundary_2d ( 1.0E+00, 0.0E+00, xse, yse )
  call boundary_2d ( 1.0E+00, 1.0E+00, xne, yne )
  call boundary_2d ( 0.0E+00, 1.0E+00, xnw, ynw )
  call boundary_2d ( 0.0E+00, 0.0E+00, xsw, ysw )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '(X,Y) corners of the total region:'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) 'XSE, YSE = ', xse, yse
  write ( *, '(a,2g14.6)' ) 'XNE, YNE = ', xne, yne
  write ( *, '(a,2g14.6)' ) 'XNW, YNW = ', xnw, ynw
  write ( *, '(a,2g14.6)' ) 'XSW, YSW = ', xsw, ysw
!
!  Draw the boundary of the region.
!
  call ps_line_rgb ( 0.1E+00, 0.5E+00, 1.0E+00 )

  r = 1.0E+00
  do i = 1, nfine
    s = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  s = 1.0E+00
  do i = nfine, 1, -1
    r = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  r = 0.0E+00
  do i = nfine, 1, -1
    s = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )

  s = 0.0E+00
  do i = 1, nfine
    r = real ( i - 1 ) / real ( nfine - 1 )
    call boundary_2d ( r, s, x(i), y(i) )
  end do

  call ps_line ( nfine, x, y )
!
!  Finish up the PostScript file.
!
  call ps_page_tail

  call ps_file_tail

  call ps_file_close ( iunit )

  return
end
subroutine boundary_2d ( r, s, x, y )
!
!*******************************************************************************
!
!! BOUNDARY_2D returns (X,Y) points on a side of the boundary.
!
!
!  Discussion:
!
!    The boundary is divided into four segments:
!
!      BOTTOM: ( 3 * cos ( ( 3-2*r)*pi/4 ),   3 * sin ( ( 3-2*r)*pi/4 ) );
!      RIGHT:  ( (3+s) * sqrt(2)/2,          (3+s)*sqrt(2)/2            );
!      TOP:    ( 4 * cos ( ( 2*r+1)*pi/4 ),   4 * sin ( (2*r+1)*pi/4 )  );
!      LEFT:   ( -(4-s)*sqrt(2)/2,           (4-s)*sqrt(2)/2            );
!
!    A
!    |
!    1  *-----------*
!    |  |           |
!    S  |           |
!    |  |           |
!    0  *-----------*
!    |
!    +--0-----R-----1--->
!
!    I'm assuming that R and S both go from 0 to 1 exactly.  Since writing
!    this code, I've come to prefer to allow the more general case where
!    the ranges of R and S are allowed to be other values.
!
!  Reference:
!
!    W N Gordon and Charles A Hall,
!    Construction of Curvilinear Coordinate Systems and Application to
!      Mesh Generation,
!    International Journal of Numerical Methods in Engineering,
!    Volume 7, pages 461-477, 1973.
!
!    Joe Thompson, Bharat Soni, Nigel Weatherill,
!    Handbook of Grid Generation,
!    CRC Press,
!    1999.
!
!  Modified:
!
!    10 December 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real R, S, the (R,S) coordinates of a point on the boundary.
!
!    Output, real X, Y, the (X,Y) coordinates of the point.
!
  implicit none
!
  real angle
  real, parameter :: pi = 3.14159265E+00
  real r
  real radius
  real s
  real x
  real y
!
  if ( s == 0.0E+00 ) then
    angle = ( 3.0E+00 - 2.0E+00 * r ) * pi / 4.0E+00
    x = 3.0E+00 * cos ( angle )
    y = 3.0E+00 * sin ( angle )
  else if ( r == 1.0E+00 ) then
    x = ( 3.0E+00 + s ) * sqrt ( 2.0E+00 ) / 2.0E+00
    y = ( 3.0E+00 + s ) * sqrt ( 2.0E+00 ) / 2.0E+00
  else if ( s == 1.0E+00 ) then
    angle = ( 3.0E+00 - 2.0E+00 * r ) * pi / 4.0E+00
    radius = 4.0E+00 + 0.2E+00 * sin ( 4.0E+00 * pi * r )
    x = radius * cos ( angle )
    y = radius * sin ( angle )
  else if ( r == 0.0E+00 ) then
    x = - ( 3.0E+00 + s ) * sqrt ( 2.0E+00 ) / 2.0E+00
    y =   ( 3.0E+00 + s ) * sqrt ( 2.0E+00 ) / 2.0E+00
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'BOUNDARY:'
    write ( *, '(a)' ) '  Illegal side coordinates!'
    write ( *, '(a,2g14.6)' ) '  (R,S) = ', r, s
    stop
  end if

  return
end
subroutine blend1d2d ( r, s, x, y )
!
!*******************************************************************************
!
!! BLEND1D2D uses transfinite interpolation on a 2D cell.
!
!
!  Reference:
!
!    W N Gordon and Charles A Hall,
!    Construction of Curvilinear Coordinate Systems and Application to
!      Mesh Generation,
!    International Journal of Numerical Methods in Engineering,
!    Volume 7, pages 461-477, 1973.
!
!    Joe Thompson, Bharat Soni, Nigel Weatherill,
!    Handbook of Grid Generation,
!    CRC Press,
!    1999.
!
!  Modified:
!
!    02 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real R, S, the (R,S) coordinates of a point.
!
!    Output, real X, Y, the (X,Y) coordinates of the point.
!
  implicit none
!
  real r
  real s
  real x
  real xe
  real xn
  real xne
  real xnw
  real xs
  real xse
  real xsw
  real xw
  real y
  real ye
  real yn
  real yne
  real ynw
  real ys
  real yse
  real ysw
  real yw
!
!  Find the (X,Y) coordinates of the corners.
!
  call boundary_2d ( 1.0E+00, 0.0E+00, xse, yse )
  call boundary_2d ( 1.0E+00, 1.0E+00, xne, yne )
  call boundary_2d ( 0.0E+00, 1.0E+00, xnw, ynw )
  call boundary_2d ( 0.0E+00, 0.0E+00, xsw, ysw )
!
!  Find the (X,Y) coordinates of corresponding points on the sides.
!
  call boundary_2d ( r, 1.0E+00, xn, yn )
  call boundary_2d ( r, 0.0E+00, xs, ys )
  call boundary_2d ( 0.0E+00, s, xw, yw )
  call boundary_2d ( 1.0E+00, s, xe, ye )
!
!  Now interpolate the (X,Y) coordinates of the point in the interior.
!
  x =                          ( 1.0E+00 - s ) * xs &
           -             r   * ( 1.0E+00 - s ) * xse &
           +             r                     * xe &
           -             r   *             s   * xne &
           +                               s   * xn &
           - ( 1.0E+00 - r ) *             s   * xnw &
           + ( 1.0E+00 - r )                   * xw &
           - ( 1.0E+00 - r ) * ( 1.0E+00 - s ) * xsw

  y =                          ( 1.0E+00 - s ) * ys &
           -             r   * ( 1.0E+00 - s ) * yse &
           +             r                     * ye &
           -             r   *             s   * yne &
           +                               s   * yn &
           - ( 1.0E+00 - r ) *             s   * ynw &
           + ( 1.0E+00 - r )                   * yw &
           - ( 1.0E+00 - r ) * ( 1.0E+00 - s ) * ysw

  return
end
subroutine pattern ( iend, max_end, num_end, max_pat, num_pat, xpat, ypat )
!
!*******************************************************************************
!
!! PATTERN defines the tiling pattern.
!
! 
!  Reference:
!
!    W N Gordon and Charles A Hall,
!    Construction of Curvilinear Coordinate Systems and Application to
!      Mesh Generation,
!    International Journal of Numerical Methods in Engineering,
!    Volume 7, pages 461-477, 1973.
!
!    Joe Thompson, Bharat Soni, Nigel Weatherill,
!    Handbook of Grid Generation,
!    CRC Press,
!    1999.
!   
!  Modified:
!
!    02 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IEND(MAX_END), contains the indices of XPAT and YPAT
!    that represent ends of line segments.  The first line segment involves
!    the entries IEND(1) + 1 through IEND(2) of XPAT and YPAT, and the last
!    one involves IEND(NUM_END-1) + 1 through IEND(NUM_END).
!
!    Input, integer MAX_END, the maximum number of entries in IEND.
!
!    Output, integer NUM_END, the number of entries used in IEND, and one
!    more than the number of line segments used to draw the pattern.
!
!    Input, integer MAX_PAT, the maximum number of entries in XPAT and YPAT.
!
!    Output, integer NUM_PAT, the number of entries used in XPAT and YPAT.
!
!    Output, real XPAT(MAX_PAT), YPAT(MAX_PAT), the X and Y coordinates of
!    points that define the line segments used to draw the pattern.
!
  implicit none
!
  integer max_end
  integer max_pat
!
  integer iend(max_end)
  integer num_end
  integer num_pat
  real xpat(max_pat)
  real ypat(max_pat)

  iend(1) = 0

  xpat(1) = 0.0E+00
  ypat(1) = 0.8E+00
  xpat(2) = 0.1E+00
  ypat(2) = 0.83E+00
  xpat(3) = 0.26E+00
  ypat(3) = 0.9E+00
  xpat(4) = 0.3E+00
  ypat(4) = 1.0E+00

  iend(2) = 4

  xpat(5) = 0.1E+00
  ypat(5) = 0.83E+00
  xpat(6) = 0.2E+00
  ypat(6) = 0.6E+00

  iend(3) = 6

  xpat(7) = 0.0E+00
  ypat(7) = 0.4E+00
  xpat(8) = 0.1E+00
  ypat(8) = 0.3E+00
  xpat(9) = 0.2E+00
  ypat(9) = 0.4E+00

  iend(4) = 9

  xpat(10) = 0.1E+00
  ypat(10) = 0.3E+00
  xpat(11) = 0.2E+00
  ypat(11) = 0.3E+00
  xpat(12) = 0.3E+00
  ypat(12) = 0.0E+00

  iend(5) = 12

  xpat(13) = 0.3E+00
  ypat(13) = 0.5E+00
  xpat(14) = 0.4E+00
  ypat(14) = 0.3E+00
  xpat(15) = 0.5E+00
  ypat(15) = 0.3E+00
  xpat(16) = 0.5E+00
  ypat(16) = 0.1E+00
  xpat(17) = 0.6E+00
  ypat(17) = 0.0E+00

  iend(6) = 17

  xpat(18) = 0.4E+00
  ypat(18) = 0.9E+00
  xpat(19) = 0.5E+00
  ypat(19) = 0.8E+00
  xpat(20) = 0.5E+00
  ypat(20) = 0.6E+00
  xpat(21) = 0.6E+00
  ypat(21) = 0.5E+00
  xpat(22) = 0.5E+00
  ypat(22) = 0.3E+00
  xpat(23) = 0.6E+00
  ypat(23) = 0.2E+00
  xpat(24) = 0.7E+00
  ypat(24) = 0.0E+00

  iend(7) = 24

  xpat(25) = 0.7E+00
  ypat(25) = 1.0E+00
  xpat(26) = 0.8E+00
  ypat(26) = 0.8E+00
  xpat(27) = 0.9E+00
  ypat(27) = 0.9E+00
  xpat(28) = 1.0E+00
  ypat(28) = 0.8E+00
  xpat(29) = 0.9E+00
  ypat(29) = 0.7E+00
  xpat(30) = 0.9E+00
  ypat(30) = 0.5E+00
  xpat(31) = 0.9E+00
  ypat(31) = 0.2E+00
  xpat(32) = 1.0E+00
  ypat(32) = 0.4E+00

  iend(8) = 32

  xpat(33) = 0.6E+00
  ypat(33) = 0.5E+00
  xpat(34) = 0.8E+00
  ypat(34) = 0.6E+00
  xpat(35) = 0.9E+00
  ypat(35) = 0.5E+00

  iend(9) = 35

  xpat(36) = 0.5E+00
  ypat(36) = 0.8E+00
  xpat(37) = 0.6E+00
  ypat(37) = 1.0E+00

  iend(10) = 37

  num_pat = 37
  num_end = 10

  return
end