c
c     Copyright (c) 1999,2000, Compaq Computer Corporation
c
c     See the file "license.terms" for information on usage and redistribution
c     of this file, and for a DISCLAIMER OF ALL WARRANTIES.
c
      subroutine hidden(pts_x, pts_y, nx, ny, ikind)
      implicit none
c
      include 'threed.cmn'
c
c     For debugging
c      include 'sys$library:gksdefs.for'
c
      integer nx, ny, ikind
      real pts_x(nx*ny), pts_y(nx*ny)
c
      integer ABOVE, BELOW, INVISIBLE
      parameter (ABOVE=1, BELOW=2, INVISIBLE=3)
c
      integer TOP_SURFACE, BOTTOM_SURFACE, BOTH_SIDES
      parameter (TOP_SURFACE=1, BOTTOM_SURFACE=2, BOTH_SIDES=3)
c
      logical intersects, intersects2, do_top, do_bottom, folds_back
      integer i, j, k, offset, status
      integer n_add_min, n_add_max
      real ymin, ymax, xcur, ycur, xnext, ynext, x_inter, y_inter
      real add_to_max(2,4), add_to_min(2,4)
      real dx1, dx2, dy1, dy2, cp
c
c     Initialize:
c
      if( ikind .lt. 1 .or. ikind .gt. 3 ) ikind = BOTH_SIDES
c
      if( ikind .eq. TOP_SURFACE .or. ikind .eq. BOTH_SIDES ) then
         do_top = .TRUE.
      else
         do_top = .FALSE.
      endif
      if( ikind .eq. BOTTOM_SURFACE .or. ikind .eq. BOTH_SIDES ) then
         do_bottom = .TRUE.
      else
         do_bottom = .FALSE.
      endif
c
      nhorizon_max = nx + ny - 1
      nhorizon_min = nx + ny - 1
      do i=1,ny
         horz_max_x(i) = pts_x((ny-i)*nx + 1)
         horz_min_x(i) = pts_x((ny-i)*nx + 1)
         horz_max_y(i) = pts_y((ny-i)*nx + 1)
         horz_min_y(i) = pts_y((ny-i)*nx + 1)
      end do
      do i=2,nx
         horz_max_x(ny+i-1) = pts_x(i)
         horz_min_x(ny+i-1) = pts_x(i)
         horz_max_y(ny+i-1) = pts_y(i)
         horz_min_y(ny+i-1) = pts_y(i)
      end do
c
c     Test for horizon which "folds back" on itself
      folds_back = .FALSE.
      do i=1,nhorizon_max - 1
         if( horz_max_x(i+1) .lt. horz_max_x(i) ) folds_back = .TRUE.
      end do
      if( folds_back ) then
         call plt$error_message('hidden',
     1      'Initial horizon is not monotonically increasing')
         return
      endif
c
c     Draw the border
      call GKS$SET_PLINE_COLOR_INDEX(1)
      if( do_top ) call GKS$POLYLINE(nx+ny-1,horz_max_x,horz_max_y)
c
      do j=2,ny
c
         offset = (j-1)*nx
c
         xcur = pts_x(offset+1)
         ycur = pts_y(offset+1)
         xnext = pts_x(offset+2)
         ynext = pts_y(offset+2)
c
c        Determine "ABOVE" or "BELOW" from sign of angle of (cur->next) to
c        (cur->previous point on initial horizon).
c        Sign of cross product of vectors will give us this info
         dx1 = pts_x((j-2)*nx+1) - xcur
         dy1 = pts_y((j-2)*nx+1) - ycur
         dx2 = xnext - xcur
         dy2 = ynext - ycur
         cp = dx1*dy2 - dx2*dy1
      call draw_line(0.0, 0.0, 0.0, 0.0, .TRUE.)
         if(cp .ge. 0.0) then
            status = ABOVE
            call GKS$SET_PLINE_COLOR_INDEX(1)
         else
            status = BELOW
            call GKS$SET_PLINE_COLOR_INDEX(2)
         endif
c
         do i=2,nx
c
            n_add_min = 0
            n_add_max = 0
c
            do k=1,2
            if( k .eq. 1 ) then
               xnext = pts_x(offset+i)
               ynext = pts_y(offset+i)
            else
               xnext = pts_x(offset+i-nx)
               ynext = pts_y(offset+i-nx)
            endif
c
            call min_max_horizon(xnext, ymin, ymax)
c
            if( ynext .eq. ymax .and. ymax .eq. ymin ) then
               if( status .eq. ABOVE ) then
                  call intersect_max(xcur, ycur, xnext, ynext,
     1                               x_inter, y_inter, intersects)
                  if( intersects ) then
                     if( do_top )
     1                  call draw_line(xcur, ycur, x_inter, y_inter,
     1                                 .FALSE.)
                     n_add_max = n_add_max + 1
                     if( n_add_max .gt. 2 )
     1                  call plt$error_message('hidden',
     2                     'Overflowed point storage')
                     add_to_max(n_add_max,1) = xcur
                     add_to_max(n_add_max,2) = ycur
                     add_to_max(n_add_max,3) = x_inter
                     add_to_max(n_add_max,4) = y_inter
                  else
                     if( do_top )
     1                  call draw_line(xcur, ycur, xnext, ynext,
     1                                 .FALSE.)
                     n_add_max = n_add_max + 1
                     if( n_add_max .gt. 2 )
     1                  call plt$error_message('hidden',
     2                     'Overflowed point storage')
                     add_to_max(n_add_max,1) = xcur
                     add_to_max(n_add_max,2) = ycur
                     add_to_max(n_add_max,3) = xnext
                     add_to_max(n_add_max,4) = ynext
                  endif
               else if( status .eq. BELOW ) then
c                 Possibilities are: intersect with min, intersect with max
c                 or neither (totally below)
                  call intersect_min(xcur, ycur, xnext, ynext,
     1                               x_inter, y_inter, intersects)
                  if( intersects ) then
                     if( do_bottom ) then
                        call draw_line(xcur, ycur, x_inter, y_inter,
     1                                 .FALSE.)
                        if( n_add_min+1 .gt. 2 ) then
                            call plt$error_message('hidden',
     1                         'Overflowed point storage')
                        else
                           n_add_min = n_add_min + 1
                           add_to_min(n_add_min,1) = xcur
                           add_to_min(n_add_min,2) = ycur
                           add_to_min(n_add_min,3) = x_inter
                           add_to_min(n_add_min,4) = y_inter
                        endif
                     endif
                     if( k .eq. 1 ) status = INVISIBLE
                  endif
                  call intersect_max(xcur, ycur, xnext, ynext,
     1                               x_inter, y_inter, intersects2)
                  if( intersects2 ) then
      call draw_line(0.0, 0.0, 0.0, 0.0, .TRUE.)
                     call GKS$SET_PLINE_COLOR_INDEX(1)
                     if( do_top ) then
                        call draw_line(x_inter, y_inter, xnext, ynext,
     1                                 .FALSE.)
                        if( n_add_max+1 .gt. 2 ) then
                            call plt$error_message('hidden',
     1                         'Overflowed point storage')
                        else
                           n_add_max = n_add_max + 1
                           add_to_max(n_add_max,1) = x_inter
                           add_to_max(n_add_max,2) = y_inter
                           add_to_max(n_add_max,3) = xnext
                           add_to_max(n_add_max,4) = ynext
                        endif
                     endif
                     if( k .eq. 1 ) status = ABOVE
                  endif
c                 If point neither intersects min or max, it must
c                 be totally below.
                  if( .not. intersects .and. .not. intersects2 ) then
                     if( do_bottom ) then
                        call draw_line(xcur, ycur, xnext, ynext,
     1                                 .FALSE.)
                        if( n_add_min+1 .gt. 2 ) then
                            call plt$error_message('hidden',
     1                         'Overflowed point storage')
                        else
                           n_add_min = n_add_min + 1
                           add_to_min(n_add_min,1) = xcur
                           add_to_min(n_add_min,2) = ycur
                           add_to_min(n_add_min,3) = xnext
                           add_to_min(n_add_min,4) = ynext
                        endif
                     endif
                  endif
               else
c                 Possibilities are: intersect with min, intersect with max
                  call intersect_min(xcur, ycur, xnext, ynext,
     1                               x_inter, y_inter, intersects)
                  if( intersects ) then
      call draw_line(0.0, 0.0, 0.0, 0.0, .TRUE.)
                     call GKS$SET_PLINE_COLOR_INDEX(2)
                     if( do_bottom ) then
                        call draw_line(x_inter, y_inter, xnext, ynext,
     1                                 .FALSE.)
                        if( n_add_min+1 .gt. 2 ) then
                            call plt$error_message('hidden',
     1                         'Overflowed point storage')
                        else
                           n_add_min = n_add_min + 1
                           add_to_min(n_add_min,1) = x_inter
                           add_to_min(n_add_min,2) = y_inter
                           add_to_min(n_add_min,3) = xnext
                           add_to_min(n_add_min,4) = ynext
                        endif
                     endif
                     if( k .eq. 1 ) status = BELOW
                  endif
                  call intersect_max(xcur, ycur, xnext, ynext,
     1                               x_inter, y_inter, intersects)
                  if( intersects ) then
      call draw_line(0.0, 0.0, 0.0, 0.0, .TRUE.)
                     call GKS$SET_PLINE_COLOR_INDEX(1)
                     if( do_top ) then
                        call draw_line(x_inter, y_inter, xnext, ynext,
     1                                 .FALSE.)
                        if( n_add_max+1 .gt. 2 ) then
                            call plt$error_message('hidden',
     1                         'Overflowed point storage')
                        else
                           n_add_max = n_add_max + 1
                           add_to_max(n_add_max,1) = x_inter
                           add_to_max(n_add_max,2) = y_inter
                           add_to_max(n_add_max,3) = xnext
                           add_to_max(n_add_max,4) = ynext
                        endif
                     endif
                     if( k .eq. 1 ) status = ABOVE
                  endif
               endif
            else if( ynext .ge. ymax .and. status .eq. ABOVE ) then
               if( do_top ) call draw_line(xcur, ycur, xnext, ynext,
     1                                     .FALSE.)
               n_add_max = n_add_max + 1
               if( n_add_max .gt. 2 )
     1            call plt$error_message('hidden',
     2                                   'Overflowed point storage')
               add_to_max(n_add_max,1) = xcur
               add_to_max(n_add_max,2) = ycur
               add_to_max(n_add_max,3) = xnext
               add_to_max(n_add_max,4) = ynext
c
            else if( ynext .le. ymin .and. status .eq. BELOW ) then
               if( do_bottom ) call draw_line(xcur, ycur, xnext, ynext,
     1                                 .FALSE.)
               n_add_min = n_add_min + 1
               if( n_add_min .gt. 2 )
     1            call plt$error_message('hidden',
     2                                   'Overflowed point storage')
               add_to_min(n_add_min,1) = xcur
               add_to_min(n_add_min,2) = ycur
               add_to_min(n_add_min,3) = xnext
               add_to_min(n_add_min,4) = ynext
c
            else if( ynext .lt. ymax .and. status .eq. ABOVE ) then
               call intersect_max(xcur, ycur, xnext, ynext,
     1                            x_inter, y_inter, intersects)
               if( intersects ) then
                  if( do_top )
     1               call draw_line(xcur, ycur, x_inter, y_inter,
     1                              .FALSE.)
                  n_add_max = n_add_max + 1
                  if( n_add_max .gt. 2 )
     1               call plt$error_message('hidden',
     2                                      'Overflowed point storage')
                  add_to_max(n_add_max,1) = xcur
                  add_to_max(n_add_max,2) = ycur
                  add_to_max(n_add_max,3) = x_inter
                  add_to_max(n_add_max,4) = y_inter
               endif
               if( ynext .lt. ymin ) then
                  call intersect_min(xcur, ycur, xnext, ynext,
     1                 x_inter, y_inter, intersects)
                  if( .NOT. intersects ) then
                     call plt$error_message('hidden',
     1                  'Expected intersection (ynext <= ymin & ABOVE)')
                  else
                     if( k .eq. 1 ) status = BELOW
      call draw_line(0.0, 0.0, 0.0, 0.0, .TRUE.)
                     call GKS$SET_PLINE_COLOR_INDEX(2)
                     if( do_bottom )
     1                  call draw_line(x_inter, y_inter, xnext, ynext,
     1                                 .FALSE.)
                     n_add_min = n_add_min + 1
                     if( n_add_min .gt. 2 )
     1               call plt$error_message('hidden',
     2                                      'Overflowed point storage')
                     add_to_min(n_add_min,1) = x_inter
                     add_to_min(n_add_min,2) = y_inter
                     add_to_min(n_add_min,3) = xnext
                     add_to_min(n_add_min,4) = ynext
                  endif
               else
                  if( k .eq. 1 ) status = INVISIBLE
               endif
c
            else if( ynext .le. ymin .and. status .eq. INVISIBLE ) then
               call intersect_min(xcur, ycur, xnext, ynext,
     1                 x_inter, y_inter, intersects)
               if( intersects ) then
                  if( k .eq. 1 ) status = BELOW
      call draw_line(0.0, 0.0, 0.0, 0.0, .TRUE.)
                  call GKS$SET_PLINE_COLOR_INDEX(2)
                  if( do_bottom )
     1               call draw_line(x_inter, y_inter, xnext, ynext,
     1                                 .FALSE.)
                  n_add_min = n_add_min + 1
                  if( n_add_min .gt. 2 )
     1               call plt$error_message('hidden',
     2                                      'Overflowed point storage')
                  add_to_min(n_add_min,1) = x_inter
                  add_to_min(n_add_min,2) = y_inter
                  add_to_min(n_add_min,3) = xnext
                  add_to_min(n_add_min,4) = ynext
               endif
c
            else if( ynext .ge. ymax .and. status .eq. INVISIBLE ) then
               call intersect_max(xcur, ycur, xnext, ynext,
     1                 x_inter, y_inter, intersects)
               if( intersects ) then
                  if( k .eq. 1 ) status = ABOVE
      call draw_line(0.0, 0.0, 0.0, 0.0, .TRUE.)
                  call GKS$SET_PLINE_COLOR_INDEX(1)
                  if( do_top )
     1               call draw_line(x_inter, y_inter, xnext, ynext,
     1                                 .FALSE.)
                  n_add_max = n_add_max + 1
                  if( n_add_max .gt. 2 )
     1               call plt$error_message('hidden',
     2                                      'Overflowed point storage')
                  add_to_max(n_add_max,1) = x_inter
                  add_to_max(n_add_max,2) = y_inter
                  add_to_max(n_add_max,3) = xnext
                  add_to_max(n_add_max,4) = ynext
               endif
c
            else if( ynext .ge. ymin .and. status .eq. BELOW ) then
               call intersect_min(xcur, ycur, xnext, ynext,
     1                 x_inter, y_inter, intersects)
               if( intersects ) then
                  if( do_bottom )
     1               call draw_line(xcur, ycur, x_inter, y_inter,
     1                              .FALSE.)
                  n_add_min = n_add_min + 1
                  if( n_add_min .gt. 2 )
     1               call plt$error_message('hidden',
     2                                      'Overflowed point storage')
                  add_to_min(n_add_min,1) = xcur
                  add_to_min(n_add_min,2) = ycur
                  add_to_min(n_add_min,3) = x_inter
                  add_to_min(n_add_min,4) = y_inter
               endif
               if( ynext .gt. ymax ) then
                  call intersect_max(xcur, ycur, xnext, ynext,
     1                               x_inter, y_inter, intersects)
                  if( .NOT. intersects ) then
                     call plt$error_message('hidden',
     1                  'Expected intersection (ynext >= ymax & BELOW)')
                  else
      call draw_line(0.0, 0.0, 0.0, 0.0, .TRUE.)
                     call GKS$SET_PLINE_COLOR_INDEX(1)
                     if( do_top )
     1                  call draw_line(x_inter, y_inter, xnext, ynext,
     1                                 .FALSE.)
                     n_add_max = n_add_max + 1
                     if( n_add_max .gt. 2 )
     1               call plt$error_message('hidden',
     2                                      'Overflowed point storage')
                     add_to_max(n_add_max,1) = x_inter
                     add_to_max(n_add_max,2) = y_inter
                     add_to_max(n_add_max,3) = xnext
                     add_to_max(n_add_max,4) = ynext
                     if( k .eq. 1 ) status = ABOVE
                  endif
               else
                  if( k .eq. 1 ) status = INVISIBLE
               endif
            endif
c
         if( k .eq. 1 ) then
            xcur = xnext
            ycur = ynext
         endif
         end do
c
         call update_max_horizon(n_add_max,add_to_max)
         call update_min_horizon(n_add_min,add_to_min)
c
c        Debugging code
c         call gks$deactivate_ws(1)
c         call gks$activate_ws(2)
c         call gks$clear_ws(2,GKS$K_CLEAR_CONDITIONALLY)
c         call gks$polyline(nhorizon_max,horz_max_x,horz_max_y)
c         call gks$deactivate_ws(2)
c         call gks$activate_ws(3)
c         call gks$clear_ws(3,GKS$K_CLEAR_CONDITIONALLY)
c         call gks$polyline(nhorizon_min,horz_min_x,horz_min_y)
c         call gks$deactivate_ws(3)
c         call gks$activate_ws(1)
c        END debug
c
         end do
c
      end do
c     Flush any remaining line segments
      call draw_line(0.0, 0.0, 0.0, 0.0, .TRUE.)
c
      return
      end
c
      subroutine update_max_horizon(n, points)
      implicit none
c
      include 'threed.cmn'
c
      integer n
      real points(2,4)
c
      real xl, xr, yl, yr
      logical joined, exact_left, exact_right, found_left, found_right
      integer i, j, k, nseg, istart, iend, ishift
c
      if( n .eq. 0 .or. n .gt. 2 ) return
c
      joined = .false.
      if( n .eq. 2 ) then
         joined = points(1,3) .eq. points(2,1) .and.
     1            points(1,4) .eq. points(2,2)
         if( joined ) then
c           Treat insertion as a single segment
            nseg = 1
         else
            nseg = 2
         endif
      else
         nseg = 1
      endif
c
      do j=1,nseg
c
         found_left  = .false.
         found_right = .false.
         exact_left  = .false.
         exact_right = .false.
c
         xl = points(j,1)
         yl = points(j,2)
         xr = points(j,3)
         yr = points(j,4)
         if( joined ) then
            xr = points(2,3)
            yr = points(2,4)
         endif
c
         do i=1,nhorizon_max
c
c           Test for common case of an adjoining, unclipped region
            if( joined .and. i .lt. (nhorizon_max-1) ) then
               if( xl .eq. horz_max_x(i) .and.
     1             yl .eq. horz_max_y(i) .and.
     2             xr .eq. horz_max_x(i+2) .and.
     3             yr .eq. horz_max_y(i+2) ) then
c                 Simply replace next point by common point of two segs
                  horz_max_x(i+1) = points(1,3)
                  horz_max_y(i+1) = points(1,4)
                  return
               endif
            endif
c
c           Look for place to insert line segment
            if( horz_max_x(i) .eq. xl .and.
     1          horz_max_y(i) .eq. yl .and. .not. found_left) then
               found_left = .true.
               istart = i
               exact_left = .true.
            else if( horz_max_x(i) .gt. xl .and. .not. found_left) then
               found_left = .true.
               istart = i-1
            else if( horz_max_x(i) .eq. xr .and.
     1               horz_max_y(i) .eq. yr ) then
               found_right = .true.
               iend = i
               exact_right = .true.
            else if( horz_max_x(i) .gt. xr ) then
               found_right = .true.
               iend = i
            endif
            if( found_right ) goto 100
c
         end do
c
c        Now insert segment, shifting horizon points as needed
c
100      continue
         ishift = 3 - (iend - istart)
         if(exact_left) ishift = ishift - 1
         if(exact_right) ishift = ishift - 1
         if( joined ) ishift = ishift + 1
c
         if( (nhorizon_max + ishift) .gt. MAXHORZPTS ) then
            call plt$error_message('update_max_horizon',
     1                    'Overflowed horizon array')
            return
         endif
c
         if( ishift .gt. 0 ) then
c           Number of inserted points is greater than number we
c           are deleting. Shift values down to make room.
            do k=nhorizon_max,iend,-1
               horz_max_x(k+ishift) = horz_max_x(k)
               horz_max_y(k+ishift) = horz_max_y(k)
            end do
         else if( ishift .lt. 0 ) then
c           Number of inserted points is less than number we
c           are deleting. Shift values up to free up space.
            do k=iend,nhorizon_max
               horz_max_x(k+ishift) = horz_max_x(k)
               horz_max_y(k+ishift) = horz_max_y(k)
            end do
         endif
c
         k = istart + 1
         if( .not. exact_left ) then
            horz_max_x(k) = xl
            horz_max_y(k) = yl
            k = k + 1
         endif
         if( joined ) then
            horz_max_x(k) = points(1,3)
            horz_max_y(k) = points(1,4)
            k = k + 1
         endif
         if( .not. exact_right ) then
            horz_max_x(k) = xr
            horz_max_y(k) = yr
         endif
c
         nhorizon_max = nhorizon_max + ishift
c
      end do
c
      return
      end
c
      subroutine update_min_horizon(n, points)
      implicit none
c
      include 'threed.cmn'
c
      integer n
      real points(2,4)
c
      real xl, xr, yl, yr
      logical joined, exact_left, exact_right, found_left, found_right
      integer i, j, k, nseg, istart, iend, ishift
c
      if( n .eq. 0 .or. n .gt. 2 ) return
c
      joined = .false.
      if( n .eq. 2 ) then
         joined = points(1,3) .eq. points(2,1) .and.
     1            points(1,4) .eq. points(2,2)
         if( joined ) then
c           Treat insertion as a single segment
            nseg = 1
         else
            nseg = 2
         endif
      else
         nseg = 1
      endif
c
      do j=1,nseg
c
         found_left  = .false.
         found_right = .false.
         exact_left = .false.
         exact_right = .false.
c
         xl = points(j,1)
         yl = points(j,2)
         xr = points(j,3)
         yr = points(j,4)
         if( joined ) then
            xr = points(2,3)
            yr = points(2,4)
         endif
c
         do i=1,nhorizon_min
c
c           Test for common case of an adjoining, unclipped region
            if( joined .and. i .lt. (nhorizon_min-1) ) then
               if( xl .eq. horz_min_x(i) .and.
     1             yl .eq. horz_min_y(i) .and.
     2             xr .eq. horz_min_x(i+2) .and.
     3             yr .eq. horz_min_y(i+2) ) then
c                 Simply replace next point by common point of two segs
                  horz_min_x(i+1) = points(1,3)
                  horz_min_y(i+1) = points(1,4)
                  return
               endif
            endif
c
c           Look for place to insert line segment
            if( horz_min_x(i) .eq. xl .and.
     1          horz_min_y(i) .eq. yl .and. .not. found_left) then
               found_left = .true.
               istart = i
               exact_left = .true.
            else if( horz_min_x(i) .gt. xl .and. .not. found_left) then
               found_left = .true.
               istart = i-1
            else if( horz_min_x(i) .eq. xr .and.
     1               horz_min_y(i) .eq. yr) then
               found_right = .true.
               iend = i
               exact_right = .true.
            else if( horz_min_x(i) .gt. xr ) then
               found_right = .true.
               iend = i
            endif
            if( found_right ) goto 100
c
         end do
c
c        Now insert segment, shifting horizon points as needed
c
100      continue
         ishift = 3 - (iend - istart)
         if(exact_left) ishift = ishift - 1
         if(exact_right) ishift = ishift - 1
         if( joined ) ishift = ishift + 1
c
         if( (nhorizon_min + ishift) .gt. MAXHORZPTS ) then
            call plt$error_message('update_min_horizon',
     1                    'Overflowed horizon array')
            return
         endif
c
         if( ishift .gt. 0 ) then
c           Number of inserted points is greater than number we
c           are deleting. Shift values down to make room.
            do k=nhorizon_min,iend,-1
               horz_min_x(k+ishift) = horz_min_x(k)
               horz_min_y(k+ishift) = horz_min_y(k)
            end do
         else if( ishift .lt. 0 ) then
c           Number of inserted points is less than number we
c           are deleting. Shift values up to free up space.
            do k=iend,nhorizon_min
               horz_min_x(k+ishift) = horz_min_x(k)
               horz_min_y(k+ishift) = horz_min_y(k)
            end do
         endif
c
         k = istart + 1
         if( .not. exact_left ) then
            horz_min_x(k) = xl
            horz_min_y(k) = yl
            k = k + 1
         endif
         if( joined ) then
            horz_min_x(k) = points(1,3)
            horz_min_y(k) = points(1,4)
            k = k + 1
         endif
         if( .not. exact_right ) then
            horz_min_x(k) = xr
            horz_min_y(k) = yr
         endif
c
         nhorizon_min = nhorizon_min + ishift
c
      end do
c
      return
      end
c
      subroutine min_max_horizon(x, ymin, ymax)
      implicit none
c
      include 'threed.cmn'
c
      real x, ymin, ymax
c
      integer i
c
c      if( x .le. horz_max_x(1) ) then
c         ymax = horz_max_y(1)
c      else if( x .ge. horz_max_x(nhorizon_max) ) then
c         ymax = horz_max_y(nhorizon_max)
c      else
         do i=2,nhorizon_max
           if( x .lt. horz_max_x(i) ) then
c             interpolate
              ymax = horz_max_y(i-1) + 
     1           (horz_max_y(i) - horz_max_y(i-1))*
     2           (x - horz_max_x(i-1))/(horz_max_x(i) - horz_max_x(i-1))
              goto 100
           endif
         end do
       ymax = horz_max_y(nhorizon_max)
c      endif
100   continue
c
c      if( x .le. horz_min_x(1) ) then
c         ymin = horz_min_y(1)
c      else if( x .ge. horz_min_x(nhorizon_min) ) then
c         ymin = horz_min_y(nhorizon_min)
c      else
         do i=2,nhorizon_min
           if( x .lt. horz_min_x(i) ) then
c             interpolate
              ymin = horz_min_y(i-1) + 
     1           (horz_min_y(i) - horz_min_y(i-1))*
     2           (x - horz_min_x(i-1))/(horz_min_x(i) - horz_min_x(i-1))
              goto 200
           endif
         end do
       ymin = horz_min_y(nhorizon_min)
c      endif

200   continue
      return
      end
c
      subroutine draw_line(x1, y1, x2, y2, lflush)
      implicit none
c
      real x1, y1, x2, y2
      logical lflush
c
c     New version batches calls to GKS$POLYLINE
c
c     If lflush = .TRUE., then ignore x1,x2,y1,y2 and emit any remaining
c     points.
c
      integer NMAX, i
      parameter (NMAX = 20)
      real xtemp(NMAX+1), ytemp(NMAX+1)
c
      data i/0/
      save i, xtemp, ytemp
c
      if( lflush ) then
         if( i .gt. 0 ) then
            call GKS$POLYLINE(i, xtemp, ytemp)
            i = 0
         endif
         return
      endif
c
      if( i .eq. 0 ) then
         xtemp(1) = x1
         xtemp(2) = x2
         ytemp(1) = y1
         ytemp(2) = y2
         i = 2
      else
         if( i .ge. NMAX ) then
            call GKS$POLYLINE(i, xtemp, ytemp)
            xtemp(1) = x1
            xtemp(2) = x2
            ytemp(1) = y1
            ytemp(2) = y2
            i = 2
         else if( xtemp(i) .eq. x1 .and. ytemp(i) .eq. y1 ) then
            i = i + 1
            xtemp(i) = x2
            ytemp(i) = y2
         else if( xtemp(i-1) .eq. x1 .and. ytemp(i-1) .eq. y1 ) then
c           Note, "i" will always be at least two, so we
c           can safely use xtemp(i-1) in above if test.
c           Also, the size of x/ytemp was increased by one, so we
c           can always safely add two points.
c           This section was added because of the peculiar way
c           the hidden line drawing routine emits line segments.
            xtemp(i+1) = xtemp(i-1)
            ytemp(i+1) = ytemp(i-1)
            i = i + 2
            xtemp(i) = x2
            ytemp(i) = y2
         else
            call GKS$POLYLINE(i, xtemp, ytemp)
            xtemp(1) = x1
            xtemp(2) = x2
            ytemp(1) = y1
            ytemp(2) = y2
            i = 2
         endif
      endif
c
      return
      end
c
      subroutine intersect_max(x1, y1, x2, y2, x_i, y_i, intersects)
      implicit none
c
      include 'threed.cmn'
c
      real x1, y1, x2, y2, x_i, y_i
      logical intersects
c
      integer i, istart, iend
c
      intersects = .false.
c
c     Find range of horizon maximum which might intersect segment.
c     NOTE WELL. If horizon function "folds back" on itself,
c     this function should still work, but I haven't gone through all the
c     ramifications.
c
c     NOTE WELL: It is possible for the segment to intesect with the
c     horizon function in multiple places. We only return the
c     first intersection, assuming this will be the most common case.
c     THIS IS AN APPROXIMATION!
c
c     Left range
      istart = 0
      do i = 1,nhorizon_max
         if( horz_max_x(i) .ge. x1 ) then
            istart = i
            goto 100
         endif
      end do
100   continue
c     Test for segment to right of horizon function
      if( istart .eq. 0 ) return
      if( istart .ne. 1 ) istart = istart - 1
c
c     Right range
      iend = 0
      do i = istart,nhorizon_max
         if( horz_max_x(i) .ge. x2 ) then
            iend = i
            goto 200
         endif
      end do
200   continue
c     Test for segment to left of horizon function (can't intersect)
      if( istart .eq. 1  .and. iend .eq. 1 ) return
      if( iend .eq. 0 ) iend = nhorizon_max
c
c     Loop to determine intersection. Segment may intersect with more
c     than one part of the horizon function, but we only handle first
c     case (this is an approximation that can be improved later).
      do i = istart,iend-1
         call try_intersect(x1,y1, x2,y2, horz_max_x(i),horz_max_y(i),
     1        horz_max_x(i+1),horz_max_y(i+1), x_i,y_i,
     2        intersects)
         if(intersects) goto 300
      end do
300   continue
c
      return
      end
c
      subroutine intersect_min(x1, y1, x2, y2, x_i, y_i, intersects)
      implicit none
c
      include 'threed.cmn'
c
      real x1, y1, x2, y2, x_i, y_i
      logical intersects
c
      integer i, istart, iend
c
      intersects = .false.
c
c     Find range of horizon minimum which might intersect segment.
c     NOTE WELL. If horizon function "folds back" on itself,
c     this function should still work, but I haven't gone through all the
c     ramifications.
c
c     NOTE WELL: It is possible for the segment to intesect with the
c     horizon function in multiple places. We only return the
c     first intersection, assuming this will be the most common case.
c     THIS IS AN APPROXIMATION!
c
c     Left range
      istart = 0
      do i = 1,nhorizon_min
         if( horz_min_x(i) .ge. x1 ) then
            istart = i
            goto 100
         endif
      end do
100   continue
c     Test for segment to right of horizon function
      if( istart .eq. 0 ) return
      if( istart .ne. 1 ) istart = istart - 1
c
c     Right range
      iend = 0
      do i = istart,nhorizon_min
         if( horz_min_x(i) .ge. x2 ) then
            iend = i
            goto 200
         endif
      end do
200   continue
c     Test for segment to left of horizon function (can't intersect)
      if( istart .eq. 1  .and. iend .eq. 1 ) return
      if( iend .eq. 0 ) iend = nhorizon_min
c
c     Loop to determine intersection. Segment may intersect with more
c     than one part of the horizon function, but we only handle first
c     case (this is an approximation that can be improved later).
      do i = istart,iend-1
         call try_intersect(x1,y1, x2,y2, horz_min_x(i),horz_min_y(i),
     1        horz_min_x(i+1),horz_min_y(i+1), x_i,y_i,
     2        intersects)
         if(intersects) goto 300
      end do
300   continue
c
      return
      end
c
      subroutine try_intersect(x1, y1, x2, y2, u1, v1, u2, v2,
     1                         xi, yi, intersects)
      implicit none
c
      include 'sysdep.cmn'
c
c     Given segments (x1,y1) to (x2,y2) and (u1,v1) to (u2,v2),
c     find if they intersect. If so, then return intersection point in
c     xi, yi. It is assumed that x2 > x1 and u2 > u1.
c
      real x1, y1, x2, y2, u1, v1, u2, v2, xi, yi
      logical intersects
c
      real m1, m2, b1, b2, miny, maxy, minv, maxv
c
      intersects = .FALSE.
c
c     Check bounding box around segments first, to see if
c     there is any chance of an intersection.
c
      miny = MIN(y1,y2)
      maxy = MAX(y1,y2)
      minv = MIN(v1,v2)
      maxv = MAX(v1,v2)
c
      if( x2 .le. u1 .or. x1 .ge. u2 .or.
     1    maxy .le. minv .or. miny .ge. maxv ) return
c
c     Test for segment endpoint match (no intersection)
c     Cases of x1 .eq. u2 and x2 .eq. u1 were take care of by above
c     check.
      if( x1 .eq. u1 .and. y1 .eq. v1 ) return
      if( x2 .eq. u2 .and. y2 .eq. v2 ) return
c
c     Compute intersection point
      m1 = PM_MAX_FLOAT
      m2 = PM_MAX_FLOAT
      if( x1 .ne. x2 ) then
         m1 = (y2 - y1)/(x2 - x1)
         b1 = y1 - m1*x1
      endif
      if( u1 .ne. u2 ) then
         m2 = (v2 - v1)/(u2 - u1)
         b2 = v1 - m2*u1
      endif
      if( m1 .eq. m2 ) return  ! Slopes equal
      if( m1 .eq. PM_MAX_FLOAT ) then
         xi = x1
         yi = x1*m2 + b2
      else if( m2 .eq. PM_MAX_FLOAT ) then
         xi = u1
         yi = u1*m1 + b1
      else
         xi = (b1 - b2)/(m2 - m1)
         yi = (m2*b1 - m1*b2)/(m2 - m1)
      endif
c
c     Finally, check to see if intersection point is within segment
c     bounding box (check on x is sufficient)
      if( x1   .lt. xi .and. xi .lt. x2   .and.
     1    u1   .lt. xi .and. xi .lt. u2   ) intersects = .TRUE.
c
      return
      end

