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 color_contours(array, val, nval, x, y, nx, ny, indx,
     1                  xlow, xhi, ylow, yhi,
     2                  max_cont_colors, do_y_flip,  y_sum)
      implicit none
c
c     Programer: John Faricelli 18-Apr-89
c
c     This FORTRAN program performs a color contour map plot of a matrix of
c     data values. Input to the program consists of a real matrix
c     of data values, a vector containing the contour values to be
c     plotted (in val). 
c
c     max_cont_colors are number of colors reserved for color contours
c     (numbered 2..max_cont_colors+1)
c
      logical do_y_flip
      integer indx, nval, nx, ny, max_cont_colors
      real array(indx,1), val(nval), y_sum, xlow, xhi, ylow, yhi
      real x(nx), y(ny)
c
      include 'gksdefs.cmn'
      include 'sysdep.cmn'
c
c     locals
c
      real FUZZ
      parameter (FUZZ = 1.0E-6)
c
      logical box_saved
      integer l, m, npts, ix1, ix2, iy1, iy2, ival, offset
      real  x_saved(4), y_saved(4),
     1      x_band(8), y_band(8)
      real cvalue1, cvalue2, zmin, zmax, x1, x2, y1, y2
c
c     Find subset of array to be plotted, so we don't do more work
c     than necessary.
      call find_greatest_subset_of_array(x, y, xlow, xhi, ylow, yhi,
     1        nx, ny, ix1, ix2, iy1, iy2)
c
c     Color map for contour map is assumed to go in order from 2 to n+1
c     Color #1 is reserved for axes and annotations.
      call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
c
c     do for each contour line
      do 3 ival = 1, nval+1
         if( ival .eq. 1 ) then
            cvalue1  = -PM_MAX_FLOAT
         else
            cvalue1  = val(ival-1)
         endif
         if( ival .eq. nval+1 ) then
            cvalue2  = PM_MAX_FLOAT
         else
            cvalue2 = val(ival)
         endif
c        Check for the case of a contour value of zero.
c        If so, boost it up a little, so the contour plotting
c        routine doesn't go crazy if the data has some
c        numerical "noise" around zero.
         if( cvalue1 .eq. 0.0 ) cvalue1 = FUZZ
         if( cvalue2 .eq. 0.0 ) cvalue2 = FUZZ
c
c        Colors are numbered 2..max_cont_colors+1 (reserve color 1 for axes)
         offset = MAX(max_cont_colors - nval,1) - 2
         call GKS$SET_FILL_COLOR_INDEX
     1           (MOD(offset+ival,max_cont_colors)+2)
c
c        Scan array, generating line segments
         do 1 l = ix1, ix2-1
            box_saved = .FALSE.
            x1 = x(l)
            x2 = x(l+1)
*** test. jvf
            if( x1 .ge. xhi .or. x2 .le. xlow ) goto 200
            do 2 m=iy1, iy2-1
c              Skip to bottom of loop if cvalue outside of data
               zmin = MIN(array(l,m),array(l,m+1),
     1                    array(l+1,m),array(l+1,m+1))
               zmax = MAX(array(l,m),array(l,m+1),
     1                    array(l+1,m),array(l+1,m+1))
               if( zmin .gt. cvalue2 .or.  cvalue1 .gt. zmax ) goto 2
c
c              If box is totally within a color band, then save
c              box for output with adjoining boxes.
               if( zmin .ge. cvalue1 .and. zmax .le. cvalue2 ) then
                  if( box_saved ) then
                     y_saved(3) = y(m+1)
                     y_saved(4) = y(m+1)
                  else
                     x_saved(1) = x1
                     x_saved(2) = x2
                     x_saved(3) = x2
                     x_saved(4) = x1
                     y_saved(1) = y(m)
                     y_saved(2) = y(m)
                     y_saved(3) = y(m+1)
                     y_saved(4) = y(m+1)
                     box_saved = .TRUE.
                  endif
                  goto 2
               endif
c              Else we have a partial box
c              Emit any full boxes we may have stored
               if( box_saved ) then
                  call fill_color_band(x_saved, y_saved, 4, y_sum,
     1                                 do_y_flip)
                  box_saved = .FALSE.
               endif
c              Construct partial box fill area
               npts = 0
               y1 = y(m)
               y2 = y(m+1)
c              Lower X
               call find_band(x1, x2, y1, array(l,m), array(l+1,m),
     1                        cvalue1, cvalue2, x_band, y_band, npts)
c              Right Y
               call find_band(y1, y2, x2, array(l+1,m), array(l+1,m+1),
     1                        cvalue1, cvalue2, y_band, x_band, npts)
c              Upper X
               call find_band(x2, x1, y2, array(l+1,m+1), array(l,m+1),
     1                        cvalue1, cvalue2, x_band, y_band, npts)
c              Left Y
               call find_band(y2, y1, x1, array(l,m+1), array(l,m),
     1                        cvalue1, cvalue2, y_band, x_band, npts)
c              This should "never" happend, but test for it anyway
               if( npts .le. 0 .or. npts .gt. 8 ) then
                  call plt$error_message('color_contours',
     1               'Internal error in color band generation')
               else
c                 Emit partial box fill area
                  call fill_color_band(x_band, y_band, npts, y_sum,
     1                                 do_y_flip)
               endif
c
2           continue
200         continue
c           Emit any full boxes we may have stored before we switch to
c           a new x value.
            if( box_saved ) then
               call fill_color_band(x_saved, y_saved, 4, y_sum,
     1                              do_y_flip)
               box_saved = .FALSE.
            endif
1        continue
c        Next cvalue
3     continue
c
      return
      end
c
      subroutine fill_color_band(x, y, npts, y_sum, do_y_flip)
      implicit none
c
      integer npts
      real x(npts), y(npts), y_sum
      logical do_y_flip
c
      integer i, n
      real xtemp(10), ytemp(10), yoff(10), t1, t2
c
      if( npts .gt. 10 ) return
c
c     Delete any duplicate points we might have created
      n = 1
      xtemp(1) = x(1)
      ytemp(1) = y(1)
      do i=2,npts
c        What this code is trying to do is figure out if two x,y pairs
c        are either identical, or close within one machine precision.
         t1 = nearest(x(i-1),x(i)-x(i-1))
         t2 = nearest(y(i-1),y(i)-y(i-1))
         if( (x(i-1) .ne. x(i) .and. t1 .ne. x(i)) .or.
     1       (y(i-1) .ne. y(i) .and. t2 .ne. y(i)) ) then
            n = n + 1
            xtemp(n) = x(i)
            ytemp(n) = y(i)
         endif
      end do
      if( n .le. 2 ) return
c
      if( do_y_flip ) then
         do i=1,n
            yoff(i) = (y_sum - ytemp(i))
         end do
         call GKS$FILL_AREA(n, xtemp, yoff)
      else
         call GKS$FILL_AREA(n, xtemp, ytemp)
      endif
c
      return
      end
c
      subroutine find_band(u1, u2, v, f1, f2,
     1                     cvalue1, cvalue2, coord1, coord2, npts)
      implicit none
c
      integer npts
      real u1, u2, v, f1, f2, cvalue1, cvalue2
      real coord1(8), coord2(8)
c
c     Given a linear function on the interval (u1,u2; v,v) with function values
c     (f1,f2), find the coordinates where the function is
c     cvalue1 <= f <= cvalue2. Place interpolated u coordinates in
c     array coord1, v in coord2.
c
c     u2 may be larger than u1.
c     cvalue2 must be greater than cvalue1.
c
      logical add_point
      real temp, du, df
c
c     Dispose of impossible cases first
      if( MAX(f1,f2) .lt. cvalue1 .or. MIN(f1,f2) .gt. cvalue2 ) return
c
      du = u2 - u1
      df = f2 - f1
c
      if( f1 .gt. cvalue2 ) then
c        We are going from a high to low area. Look for intersection
c        with cvalue2.
         temp = u1 + (cvalue2-f1)*du/df
         add_point = .FALSE.
         if( du .ge. 0.0 ) then
            if( u1 .le. temp .and. temp .le. u2 ) add_point = .TRUE.
         else
            if( u2 .le. temp .and. temp .le. u1 ) add_point = .TRUE.
         endif
         if( add_point ) then
            npts = npts + 1
            coord1(npts) = temp
            coord2(npts) = v
         endif
         if( add_point ) then
c           Look for end point in contour range
            if( f2 .le. cvalue2 .and. f2 .ge. cvalue1 ) then
               npts = npts + 1
               coord1(npts) = u2
               coord2(npts) = v
            else
c              Look for intersection with cvalue1
               temp = u1 + (cvalue1-f1)*du/df
               if( du .ge. 0.0 ) then
                  if( u1 .le. temp .and. temp .le. u2 )
     1               add_point = .TRUE.
               else
                  if( u2 .le. temp .and. temp .le. u1 )
     1               add_point = .TRUE.
               endif
               if( add_point ) then
                  npts = npts + 1
                  coord1(npts) = temp
                  coord2(npts) = v
               endif
            endif
         endif
      else if( f1 .lt. cvalue1 ) then
c        We are going from a low area to a high area. Look for intersection
c        with cvalue1.
         temp = u1 + (cvalue1-f1)*du/df
         add_point = .FALSE.
         if( du .ge. 0.0 ) then
            if( u1 .le. temp .and. temp .le. u2 ) add_point = .TRUE.
         else
            if( u2 .le. temp .and. temp .le. u1 ) add_point = .TRUE.
         endif
         if( add_point ) then
            npts = npts + 1
            coord1(npts) = temp
            coord2(npts) = v
         endif
         if( add_point ) then
c           Look for end point in contour range
            if( f2 .le. cvalue2 .and. f2 .ge. cvalue1 ) then
               npts = npts + 1
               coord1(npts) = u2
               coord2(npts) = v
            else
c              Look for intersection with cvalue2
               temp = u1 + (cvalue2-f1)*du/df
               if( du .ge. 0.0 ) then
                  if( u1 .le. temp .and. temp .le. u2 )
     1               add_point = .TRUE.
               else
                  if( u2 .le. temp .and. temp .le. u1 )
     1               add_point = .TRUE.
               endif
               if( add_point ) then
                  npts = npts + 1
                  coord1(npts) = temp
                  coord2(npts) = v
               endif
            endif
         endif
      else
c        Point is inbetween contour values. Include it.
         npts = npts + 1
         coord1(npts) = u1
         coord2(npts) = v
c        Check for intersections. Three possibilities:
c        Intersect with cvalue1, cvalue2 or inbetween
         if( f2 .gt. cvalue2 ) then
            temp = u1 + (cvalue2-f1)*du/df
            add_point = .FALSE.
            if( du .ge. 0.0 ) then
               if( u1 .le. temp .and. temp .le. u2 ) add_point = .TRUE.
            else
               if( u2 .le. temp .and. temp .le. u1 ) add_point = .TRUE.
            endif
            if( add_point ) then
               npts = npts + 1
               coord1(npts) = temp
               coord2(npts) = v
            endif
         else if( f2 .lt. cvalue1 ) then
            temp = u1 + (cvalue1-f1)*du/df
            add_point = .FALSE.
            if( du .ge. 0.0 ) then
               if( u1 .le. temp .and. temp .le. u2 ) add_point = .TRUE.
            else
               if( u2 .le. temp .and. temp .le. u1 ) add_point = .TRUE.
            endif
            if( add_point ) then
               npts = npts + 1
               coord1(npts) = temp
               coord2(npts) = v
            endif
         else
            npts = npts + 1
            coord1(npts) = u2
            coord2(npts) = v
         endif
c        
      endif
c
      return
      end
c
      subroutine color_contours_tri(array, cord, nop, ne,
     1                  tri2reg, imat,
     2                  indx, val, nval,
     3                  hmin, hmax, vmin, vmax,
     4                  value_in_oxide,
     5                  max_cont_colors, do_y_flip,  y_sum)
      implicit none
c
c     Programer: John Faricelli 10-May-93
c
c     This FORTRAN program performs a color contour map plot of a matrix of
c     data values. Input to the program consists of a real matrix
c     of data values, a vector containing the contour values to be
c     plotted (in val). Mesh is assumed to be on a triangular grid.
c
c     max_cont_colors are number of colors reserved for color contours
c     (numbered 2..max_cont_colors+1)
c
      logical value_in_oxide, do_y_flip
      integer indx, nval, max_cont_colors, ne
      integer nop(3,ne), tri2reg(ne), imat(1)
      real array(indx), val(nval), y_sum, hmin, hmax, vmin, vmax,
     1     cord(2,indx)
c
      include 'gksdefs.cmn'
      include 'sysdep.cmn'
c
c     locals
c
      real FUZZ
      parameter (FUZZ = 1.0E-6)
c
      integer ie, npts, ival, offset, j1, j2, j3, j4
      real x1, x2, x3, x4, y1, y2, y3, y4, cvalue1, cvalue2,
     1     zmin, zmax, fz(3), xe(3), ye(3), poly(6,2)
c
c     Color map for contour map is assumed to go in order from 2 to n+1
c     Color #1 is reserved for axes and annotations.
      call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
CCC FOR DEBUG, SET NEXT STATEMENT CCC
CCC      call gks$set_defer_state(1, GKS$K_ASAP, GKS$K_IRG_ALLOWED)
c
c     do for each triangular element
      do 3 ie = 1,ne
c        Check for insulating material
c        Note: for now, this works only for PISCES, where material code
c        for insulator is < 0.
         if(imat(tri2reg(ie)) .lt. 0 .and. .NOT. value_in_oxide) goto 3
         xe(1) = cord(1,nop(1,ie))
         ye(1) = cord(2,nop(1,ie))
         xe(2) = cord(1,nop(2,ie))
         ye(2) = cord(2,nop(2,ie))
         xe(3) = cord(1,nop(3,ie))
         ye(3) = cord(2,nop(3,ie))
c        Check: is triangle outside bounding box?
**         if( MAX(xe(1),xe(2),xe(3)) .lt. hmin .or.
**     1       MIN(xe(1),xe(2),xe(3)) .gt. hmax .or.
**     2       MAX(ye(1),ye(2),ye(3)) .lt. vmin .or.
**     3       MIN(ye(1),ye(2),ye(3)) .gt. vmax ) goto 3
         if( MAX(xe(1),xe(2),xe(3)) .le. hmin .or.
     1       MIN(xe(1),xe(2),xe(3)) .ge. hmax .or.
     2       MAX(ye(1),ye(2),ye(3)) .le. vmin .or.
     3       MIN(ye(1),ye(2),ye(3)) .ge. vmax ) goto 3
c
c        Get values at triangle nodes
         fz(1) = array(nop(1,ie))
         fz(2) = array(nop(2,ie))
         fz(3) = array(nop(3,ie))
         call pm_color_fill_contours_tri(xe, ye, fz, val, nval,
     1           max_cont_colors, y_sum, do_y_flip)
c     Next element
3     continue
c
      return
      end
c
      subroutine find_band_tri(x1, x2, y1, y2, z1, z2,
     1                     cvalue1, cvalue2, coord1, coord2, npts)
      implicit none
c
      integer npts
      real x1, x2, y1, y2, z1, z2, cvalue1, cvalue2
      real coord1(8), coord2(8)
c
c     Given a linear function on the interval (u1,u2; v1,v2)
c     with function values (f1,f2), find the coordinates where the function is
c     cvalue1 <= f <= cvalue2. Place interpolated u coordinates in
c     array coord1, interpolated v in coord2.
c
c     u2 may be larger than u1
c     v2 may be larger than v1
c     cvalue2 must be greater than cvalue1
c
      real tempu, tempv, du, dv, df, u1, u2, v1, v2, f1, f2
      logical lrev
c
c     Dispose of impossible cases first
      if( MAX(z1,z2) .lt. cvalue1 .or. MIN(z1,z2) .gt. cvalue2 ) return
c
      if( z1 .lt. z2 ) then
         u1 = x1
         u2 = x2
         v1 = y1
         v2 = y2
         f1 = z1
         f2 = z2
         lrev = .FALSE.
      else
         u1 = x2
         u2 = x1
         v1 = y2
         v2 = y1
         f1 = z2
         f2 = z1
         lrev = .TRUE.
      endif
      du = u2 - u1
      dv = v2 - v1
      df = f2 - f1
c
      if(df .eq. 0.0 .or.
     1  (cvalue1 .le. f1 .and. f2 .le. cvalue2) ) then
c        Degenerate case: use entire edge
         npts = npts + 1
         coord1(npts) = u1
         coord2(npts) = v1
         npts = npts + 1
         coord1(npts) = u2
         coord2(npts) = v2
      else if( cvalue1 .le. f1 .and. f2 .gt. cvalue2 ) then
c        One intersection point at cvalue2
         npts = npts + 1
         coord1(npts) = u1
         coord2(npts) = v1
         tempu = u1 + (f2-cvalue2)*du/df
         tempv = v1 + (f2-cvalue2)*dv/df
         npts = npts + 1
         coord1(npts) = tempu
         coord2(npts) = tempv
      else if( cvalue1 .gt. f1 .and. f2 .le. cvalue2 ) then
c        One intersection point at cvalue1
         tempu = u1 + (cvalue1-f1)*du/df
         tempv = v1 + (cvalue1-f1)*dv/df
         npts = npts + 1
         coord1(npts) = tempu
         coord2(npts) = tempv
         npts = npts + 1
         coord1(npts) = u2
         coord2(npts) = v2
      else
c        Two intersection points at cvalue1 and cvalue2
         tempu = u1 + (cvalue1-f1)*du/df
         tempv = v1 + (cvalue1-f1)*dv/df
         npts = npts + 1
         coord1(npts) = tempu
         coord2(npts) = tempv
         tempu = u1 + (f2-cvalue2)*du/df
         tempv = v1 + (f2-cvalue2)*dv/df
         npts = npts + 1
         coord1(npts) = tempu
         coord2(npts) = tempv
      endif
      if( lrev ) then
         tempu = coord1(npts-1)
         tempv = coord2(npts-1)
         coord1(npts-1) = coord1(npts)
         coord2(npts-1) = coord2(npts)
         coord1(npts) = tempu
         coord2(npts) = tempv
      endif
c
      return
      end
c
      subroutine contur_tri(array, cord, nop, ne,
     1                  tri2reg, imat,
     2                  indx, val, nval,
     3                  hmin, hmax, vmin, vmax, value_in_oxide,
     4                  make_linetype_dashed, do_y_flip, y_sum,
     5                  do_label, do_solid, do_colormat)
      implicit none
c
c     Programer: John Faricelli 10-May-93
c
c     This FORTRAN program performs a line contour map plot of a matrix of
c     data values. Input to the program consists of a real matrix
c     of data values, a vector containing the contour values to be
c     plotted (in val). Mesh is assumed to be on a triangular grid.
c
c
      logical value_in_oxide, make_linetype_dashed, do_y_flip,
     1        do_label, do_solid, do_colormat
      integer indx, nval, ne
      integer nop(3,ne), tri2reg(ne), imat(1)
      real array(indx), val(nval), y_sum, hmin, hmax, vmin, vmax,
     1     cord(2,indx)
c
      include 'gksdefs.cmn'
      include 'sysdep.cmn'
c
c     locals
c
      real FUZZ
      parameter (FUZZ = 1.0E-6)
c
c     Define parametes which allocate local storage.
c     NPOLYMAX is the longest polyline that will be emitted.
c     This parameter may have to be adjusted depending on the limitations
c     of the underlying graphics system.
      integer NSEGMAX
      parameter (NSEGMAX = 300)
      integer NPOLYMAX
      parameter (NPOLYMAX = 255)
c
      integer i, npts, lseg, ival, idum, ie, node(3),
     1        icolor, mat_color(20), num_colors_assigned
      real seg(NSEGMAX,4), xpts(5), ypts(5),
     1     xtemp(NPOLYMAX), ytemp(NPOLYMAX), x(3), y(3), z(3), yoff(3)
      real cvalue, x1, x2, y1, y2, zmin, zmax
      logical used(NSEGMAX), lout
c
      if( do_colormat ) then
         call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
         num_colors_assigned = 0
         do ie = 1,ne
            do i=1,3
               node(i) = nop(i,ie)
               x(i)    = cord(1,node(i))
               y(i)    = cord(2,node(i))
            end do
c
c           Check for element outside window
            x1 =  PM_MAX_FLOAT
            x2 = -PM_MAX_FLOAT
            y1 =  PM_MAX_FLOAT
            y2 = -PM_MAX_FLOAT
            do i=1,3
               x1 = MIN(x1,x(i))
               x2 = MAX(x2,x(i))
               y1 = MIN(y1,y(i))
               y2 = MAX(y2,y(i))
            end do
            if( hmin .gt. x2 .or. hmax .lt. x1 .or.
     1          vmin .gt. y2 .or. vmax .lt. y1 ) then
               lout = .TRUE.
            else
               lout = .FALSE.
            endif
            if( .NOT. lout ) then
               icolor = -1
               do i=1,num_colors_assigned
                  if( mat_color(i) .eq. imat(tri2reg(ie)) ) then
                     icolor = i+1
                  endif
               end do
               if( icolor .eq. -1 ) then
                  num_colors_assigned = MIN(num_colors_assigned+1,20)
                  icolor = num_colors_assigned + 1
                  mat_color(num_colors_assigned) = imat(tri2reg(ie))
               endif
               call GKS$SET_FILL_COLOR_INDEX(icolor)
               if( do_y_flip ) then
                  do i=1,3
                     yoff(i) = (y_sum - y(i))
                  end do
                  call GKS$FILL_AREA(3, x, yoff)
               else
                  call GKS$FILL_AREA(3, x, y)
               endif
            end if
         end do
      endif
c
c     do for each contour line
      do 3 ival = 1, nval
         cvalue = val(ival)
c        Check for the case of a contour value of zero.
c        If so, boost it up a little, so the contour plotting
c        routine doesn't go crazy if the data has some
c        numerical "noise" around zero.
         if( cvalue .eq. 0.0 ) cvalue = FUZZ
c
         if( make_linetype_dashed ) then
            call set_line_index(1)
            call GKS$SET_PLINE_LINETYPE(GKS$K_LINETYPE_DASHED)
         else
            if( do_solid ) then
               call set_line_index(1)
            else
               call set_line_index(ival)
            endif
         endif
         call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
c
         lseg = 0
c        Scan array, generating line segments
         do 1 ie = 1,ne
c
c           Check for insulating material with no data value
c           Note: for now, this works only for PISCES, where material code
c           for insulator is < 0.
            if(imat(tri2reg(ie)) .lt. 0 .and. .NOT. value_in_oxide)
     1         goto 1
c
            do 10 i=1,3
               node(i) = nop(i,ie)
               x(i)    = cord(1,node(i))
               y(i)    = cord(2,node(i))
               z(i)    = array(node(i))
10          continue
c
c           Check for element outside window
c
c           Check for element outside window
            x1 =  PM_MAX_FLOAT
            x2 = -PM_MAX_FLOAT
            y1 =  PM_MAX_FLOAT
            y2 = -PM_MAX_FLOAT
            do i=1,3
               x1 = MIN(x1,x(i))
               x2 = MAX(x2,x(i))
               y1 = MIN(y1,y(i))
               y2 = MAX(y2,y(i))
            end do
            if( hmin .gt. x2 .or. hmax .lt. x1 .or.
     1          vmin .gt. y2 .or. vmax .lt. y1 ) then
               lout = .TRUE.
            else
               lout = .FALSE.
            endif
            if( lout ) goto 1
c
c           Check for data values outside contour range
            zmin = MIN(z(1),z(2),z(3))
            zmax = MAX(z(1),z(2),z(3))
            if( zmin .gt. cvalue .or. zmax .lt. cvalue ) goto 1
c
            npts = 0
c           Check for case where function value is constant
            if( z(1) .eq. z(2) .and. z(2) .eq. z(3) ) then
               npts = 4
               xpts(1) = x(1)
               xpts(2) = x(2)
               xpts(3) = x(3)
               xpts(4) = x(1)
               ypts(1) = y(1)
               ypts(2) = y(2)
               ypts(3) = y(3)
               ypts(4) = y(1)
            else
               call pm_find_contour_in_tri(cvalue,z,x,y,zmax,zmin,
     1                                     x1,y1,x2,y2,idum,idum)
               npts = 2
               xpts(1) = x1
               xpts(2) = x2
               ypts(1) = y1
               ypts(2) = y2
            endif
c           Create line segments from points
            if( npts .eq. 2 ) then
               lseg = lseg + 1
               seg(lseg,1) = xpts(1)
               seg(lseg,2) = xpts(2)
               seg(lseg,3) = ypts(1)
               seg(lseg,4) = ypts(2)
            else if( npts .eq. 3 ) then
c              Bifurcation point
               lseg = lseg + 1
               seg(lseg,1) = xpts(1)
               seg(lseg,2) = xpts(2)
               seg(lseg,3) = ypts(1)
               seg(lseg,4) = ypts(2)
               lseg = lseg + 1
               seg(lseg,1) = xpts(2)
               seg(lseg,2) = xpts(3)
               seg(lseg,3) = ypts(2)
               seg(lseg,4) = ypts(3)
            else if( npts .eq. 4 ) then
c              Saddle point
               lseg = lseg + 1
               seg(lseg,1) = xpts(1)
               seg(lseg,2) = xpts(2)
               seg(lseg,3) = ypts(1)
               seg(lseg,4) = ypts(2)
               lseg = lseg + 1
               seg(lseg,1) = xpts(3)
               seg(lseg,2) = xpts(4)
               seg(lseg,3) = ypts(3)
               seg(lseg,4) = ypts(4)
            endif
c           Safety valve. Make sure we don't overrun line segment storage
c           We can add at most 2 line segments per loop trip
            if( lseg .ge. NSEGMAX-2) then
                call plot_contour_segments(seg, used, lseg, NSEGMAX,
     1                  xtemp, ytemp, NPOLYMAX, do_y_flip, y_sum,
     2                  cvalue, do_label)
                lseg = 0
            endif
1        continue
c        Now plot segments in order
         call plot_contour_segments(seg, used, lseg, NSEGMAX, 
     1           xtemp, ytemp, NPOLYMAX, do_y_flip, y_sum, cvalue,
     2           do_label)
c        Next cvalue
3     continue
c
      return
      end
c
      subroutine pm_find_contour_in_tri(conval,funct,xe,ye,fmax,fmin,
     1                                  x1,y1,x2,y2,j1,j2)
      implicit none
c
c     Find the intersections of a contour line with a triangle
c
c     conval  - the contour value
c     funct   - 3 triangle function values
c     xe      - 3 xcoords of triangle
c     ye      - 3 ycoords of triangle
c     fmax,fmin - max/min of funct
c     x1,y1   - first intersection found (in order of side number)
c     x2,y2   - second
c     j1,j2   - intersected sides.
c
      real conval, funct(3), xe(3), ye(3), fmax, fmin, x1, y1, x2, y2
      integer j1, j2
c
c     locals
c
      real x, y
      integer inext(3), iprev(3), is, isn, isp, count
      double precision tau, funfr, funto, u1, u2
      data inext/2,3,1/,iprev/3,1,2/
c
      j1=0
      j2=0
      count=0
      if( (conval.lt.fmin) .or. (conval .gt. fmax) ) return
c
      do 100 is=1,3
         isn=inext(is)
         isp=iprev(is)

c        See if this side intersects contour
         funfr=funct(isp)
         funto=funct(isn)
         if((funfr.gt.conval .and. funto.gt.conval) .or.
     1      (funfr.lt.conval .and. funto.lt.conval)) goto 100

         if(funfr.eq.funto) then
c           Oops, function equals conval on entire side.
            j1=isp
            j2=isn
            x1=xe(isp)
            y1=ye(isp)
            x2=xe(isn)
            y2=ye(isn)
            return
         endif

c        Normal case - function values span conval
         tau=(DBLE(conval) - funfr) / (funto - funfr)
         x=(xe(isn) - xe(isp))*tau + xe(isp)
         u1 = ye(isn)
         u2 = ye(isp)
c        Perform y computations in double precision to avoid small round
c        off problems that confuse the routine which assembles contour line
c        segments into matching pairs. X values do not seem to have this
c        problem.
         y=(u1 - u2)*tau + u2
         count=count+1
         if(count.eq.1) then
            x1=x
            y1=y
            j1=is
         else
            x2=x
            y2=y
            j2=is
         endif
100   continue
c
      return
      end
c
      subroutine color_contours_histo(array, val, nval, x, y,
     1                  nx, ny, indx,
     1                  xlow, xhi, ylow, yhi,
     2                  max_cont_colors, do_y_flip,  y_sum)
      implicit none
c
c     Programer: John Faricelli 2-Apr-96
c
c     This FORTRAN program performs a color contour map plot of a matrix of
c     data values. Input to the program consists of a real matrix
c     of data values, a vector containing the contour values to be
c     plotted (in val). 
c
c     max_cont_colors are number of colors reserved for color contours
c     (numbered 2..max_cont_colors+1)
c
c     It is similar to color_contours(), but plots a 'histogram' style
c     plot, with one color filling the entire box.
c
      logical do_y_flip
      integer indx, nval, nx, ny, max_cont_colors
      real array(indx,1), val(nval), y_sum, xlow, xhi, ylow, yhi
      real x(nx), y(ny)
c
      include 'gksdefs.cmn'
      include 'sysdep.cmn'
c
c     locals
c
      real FUZZ
      parameter (FUZZ = 1.0E-6)
c
      logical box_saved
      integer l, m, ix1, ix2, iy1, iy2, ival, offset
      real  x_saved(4), y_saved(4)
      real cvalue1, cvalue2, zval, x1, x2, y1, y2
c
c     Find subset of array to be plotted, so we don't do more work
c     than necessary.
      call find_greatest_subset_of_array(x, y, xlow, xhi, ylow, yhi,
     1        nx, ny, ix1, ix2, iy1, iy2)
c
c     Color map for contour map is assumed to go in order from 2 to n+1
c     Color #1 is reserved for axes and annotations.
      call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
c
c     do for each contour value
      do 3 ival = 1, nval+1
         if( ival .eq. 1 ) then
            cvalue1  = -PM_MAX_FLOAT
         else
            cvalue1  = val(ival-1)
         endif
         if( ival .eq. nval+1 ) then
            cvalue2  = PM_MAX_FLOAT
         else
            cvalue2 = val(ival)
         endif
c        Check for the case of a contour value of zero.
c        If so, boost it up a little, so the contour plotting
c        routine doesn't go crazy if the data has some
c        numerical "noise" around zero.
         if( cvalue1 .eq. 0.0 ) cvalue1 = FUZZ
         if( cvalue2 .eq. 0.0 ) cvalue2 = FUZZ
c
c        Colors are numbered 2..max_cont_colors+1 (reserve color 1 for axes)
         offset = MAX(max_cont_colors - nval,1) - 2
         call GKS$SET_FILL_COLOR_INDEX
     1           (MOD(offset+ival,max_cont_colors)+2)
c
c        Scan array, generating boxes
         do 1 l = ix1, ix2
            box_saved = .FALSE.
            if( l .eq. ix1 ) then
               x1 = x(ix1)
            else
               x1 = 0.5*(x(l) + x(l-1))
            endif
            if( l .eq. ix2 ) then
               x2 = x(ix2)
            else
               x2 = 0.5*(x(l) + x(l+1))
            endif
            do 2 m=iy1, iy2
c              Skip to bottom of loop if cvalue outside of data
               zval = array(l,m)
               if( zval .gt. cvalue2 .or.  cvalue1 .gt. zval ) then
c                 Emit any full boxes we may have stored
                  if( box_saved ) then
                     call fill_color_band(x_saved, y_saved, 4, y_sum,
     1                                    do_y_flip)
                     box_saved = .FALSE.
                  endif
                  goto 2
               endif
c              Save box for (possible) output with adjoining boxes.
               if( m .eq. iy2 ) then
                  y2 = y(m)
               else
                  y2 = 0.5*(y(m)+y(m+1))
               endif
               if( box_saved ) then
                  y_saved(3) = y2
                  y_saved(4) = y2
               else
                  if( m .eq. iy1 ) then
                     y1 = y(m)
                  else
                     y1 = 0.5*(y(m)+y(m-1))
                  endif
                  x_saved(1) = x1
                  x_saved(2) = x2
                  x_saved(3) = x2
                  x_saved(4) = x1
                  y_saved(1) = y1
                  y_saved(2) = y1
                  y_saved(3) = y2
                  y_saved(4) = y2
                  box_saved = .TRUE.
               endif
               goto 2
2           continue
c           Emit any full boxes we may have stored before we switch to
c           a new x value.
            if( box_saved ) then
               call fill_color_band(x_saved, y_saved, 4, y_sum,
     1                              do_y_flip)
               box_saved = .FALSE.
            endif
1        continue
c        Next cvalue
3     continue
c
      return
      end
c
      subroutine pm_color_fill_contours_tri(xe, ye, fz, val, nval,
     1                           max_cont_colors, y_sum, do_y_flip)
      implicit none
c
      integer nval, max_cont_colors
      real xe(3), ye(3), fz(3), val(nval), y_sum
      logical do_y_flip
c
      include 'const.h'
      include 'sysdep.cmn'
c
c     Locals
      integer npts, ival, offset, j1, j2, j3, j4
      real x1, x2, x3, x4, y1, y2, y3, y4, cvalue1, cvalue2,
     1     zmin, zmax, poly(6,2)
c
      real FUZZ
      parameter (FUZZ = 1.0E-6)
c
      zmin = MIN(fz(1),fz(2),fz(3))
      zmax = MAX(fz(1),fz(2),fz(3))
c     Do for each contour value
      do 2 ival = 1, nval+1
         if( ival .eq. 1 ) then
            cvalue1  = -PM_MAX_FLOAT
         else
            cvalue1  = val(ival-1)
         endif
         if( ival .eq. nval+1 ) then
            cvalue2  = PM_MAX_FLOAT
         else
            cvalue2 = val(ival)
         endif
c        Check for the case of a contour value of zero.
c        If so, boost it up a little, so the contour plotting
c        routine doesn't go crazy if the data has some
c        numerical "noise" around zero.
         if( cvalue1 .eq. 0.0 ) cvalue1 = FUZZ
         if( cvalue2 .eq. 0.0 ) cvalue2 = FUZZ
c
c         Scan element, generating line segments
c         Skip to bottom of loop if cvalue outside of data
         if( zmin .gt. cvalue2 .or.  cvalue1 .gt. zmax ) goto 2
c
c        Colors are numbered 2..max_cont_colors+1 (reserve color 1 for axes)
         offset = MAX(max_cont_colors - nval,1) - 2
         call GKS$SET_FILL_COLOR_INDEX
     1              (MOD(offset+ival,max_cont_colors)+2)
c
c        Check for triangle completely within contour values
c        (a common case)
         if( cvalue1 .le. zmin .and. zmax .lt. cvalue2 ) then
            npts=3
            call fill_color_band(xe, ye,
     1           npts, y_sum, do_y_flip)
c           Break out of loop and return
            goto 3
         endif
c
         npts = 0
c
c        Watch out for the degenerate case.
         if(fz(1).eq.fz(2) .and. fz(2).eq.fz(3)) then
            npts=3
            poly(1,1)=xe(1)
            poly(1,2)=ye(1)
            poly(2,1)=xe(2)
            poly(2,2)=ye(2)
            poly(3,1)=xe(3)
            poly(3,2)=ye(3)
         else
c           Interpolate edges.
            call pm_find_contour_in_tri(cvalue1, fz,xe,ye,zmax,zmin,
     1                                  x1,y1, x2,y2, j1,j2)
            call pm_find_contour_in_tri(cvalue2, fz,xe,ye,zmax,zmin,
     1                                  x3,y3, x4,y4, j3,j4)
c           Determine what polygon to draw.
c           Can be off either end, intersect two of the sides,
c           or intersect all three sides, hence the multiple cases.
         if(j1.eq.0 .and. j3.eq.0) then
            npts=3
            poly(1,1)=xe(1)
            poly(1,2)=ye(1)
            poly(2,1)=xe(2)
            poly(2,2)=ye(2)
            poly(3,1)=xe(3)
            poly(3,2)=ye(3)
         else if(j1 .eq. 0) then
            if(fz(6-j3-j4).lt.cvalue2) then
               npts=3
               poly(1,1)=xe(6-j3-j4)
               poly(1,2)=ye(6-j3-j4)
               poly(2,1)=x3
               poly(2,2)=y3
               poly(3,1)=x4
               poly(3,2)=y4
            else
               npts=4
               poly(1,1)=xe(j3)
               poly(1,2)=ye(j3)
               poly(2,1)=xe(j4)
               poly(2,2)=ye(j4)
               poly(3,1)=x3
               poly(3,2)=y3
               poly(4,1)=x4
               poly(4,2)=y4
            endif
         else if(j3 .eq. 0) then
            if(fz(6-j1-j2).gt.cvalue1) then
               npts=3
               poly(1,1)=xe(6-j1-j2)
               poly(1,2)=ye(6-j1-j2)
               poly(2,1)=x1
               poly(2,2)=y1
               poly(3,1)=x2
               poly(3,2)=y2
            else
               npts=4
               poly(1,1)=xe(j1)
               poly(1,2)=ye(j1)
               poly(2,1)=xe(j2)
               poly(2,2)=ye(j2)
               poly(3,1)=x1
               poly(3,2)=y1
               poly(4,1)=x2
               poly(4,2)=y2
            endif
         else if( (j1.eq.j3) .and. (j2.eq.j4) ) then
            npts=4
            poly(1,1)=x1
            poly(1,2)=y1
            poly(2,1)=x2
            poly(2,2)=y2
            poly(3,1)=x4
            poly(3,2)=y4
            poly(4,1)=x3
            poly(4,2)=y3
         else if(j1.eq.j3) then
            npts=5
            poly(1,1)=x1
            poly(1,2)=y1
            poly(2,1)=x2
            poly(2,2)=y2
            poly(3,1)=xe(j1)
            poly(3,2)=ye(j1)
            poly(4,1)=x4
            poly(4,2)=y4
            poly(5,1)=x3
            poly(5,2)=y3
         else if(j1.eq.j4) then
            npts=5
            poly(1,1)=x1
            poly(1,2)=y1
            poly(2,1)=x2
            poly(2,2)=y2
            poly(3,1)=xe(j1)
            poly(3,2)=ye(j1)
            poly(4,1)=x3
            poly(4,2)=y3
            poly(5,1)=x4
            poly(5,2)=y4
         else if(j2.eq.j4) then
            npts=5
            poly(1,1)=x2
            poly(1,2)=y2
            poly(2,1)=x1
            poly(2,2)=y1
            poly(3,1)=xe(j2)
            poly(3,2)=ye(j2)
            poly(4,1)=x3
            poly(4,2)=y3
            poly(5,1)=x4
            poly(5,2)=y4
         else if(j2.eq.j3) then
            npts=5
            poly(1,1)=x2
            poly(1,2)=y2
            poly(2,1)=x1
            poly(2,2)=y1
            poly(3,1)=xe(j2)
            poly(3,2)=ye(j2)
            poly(4,1)=x4
            poly(4,2)=y4
            poly(5,1)=x3
            poly(5,2)=y3
         else
c           Can't happen...
         endif
         endif
c
         if( npts .gt. 1 )
c           Do we need to sort points to be contiguous?
c           Don't think so, as long as we do points in order
c           around the triangle.
     1      call fill_color_band(poly(1,1), poly(1,2),
     2           npts, y_sum, do_y_flip)
c     Next cvalue
2     continue
c
3     continue
c
      return
      end
