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 contur(array, val, nval, x, y, nx, ny, indx,
     1                  xlow, xhi, ylow, yhi,
     2                  make_linetype_dashed, do_y_flip,  y_sum,
     3                  do_label, do_solid)
      implicit none
c
c     Programer: John Faricelli 30-Sep-88
c
c     This FORTRAN program performs a contour 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.
c
c     This subroutine replaces the "ancient" contour plotting routine
c     with a different algorithm. The new algorithm creates all the
c     contour segments in temporary storage and then plots them in
c     connected order. This greatly improves pen plotting speed.
c     The interpolation code is more transparent now, if somewhat more
c     lengthy. The new code appears to be reasonably efficient compared
c     to the old algorithm, even when there are hundreds of contour segments
c     to be plotted. In tests, the new algorithm appears to be at least
c     as fast as the old, and sometimes noticeably faster.
c
c     Logical do_label controls whether labels are placed on contour
c     lines.
c
c     Logical do_solid controls whether contours are done in solid lines
c     (line color 1). Otherwise we use a different line style to
c     distinguish them.
c
      logical make_linetype_dashed, do_y_flip, do_label, do_solid
      integer indx, nval, nx, ny
      real array(indx,1), val(1), y_sum, xlow, xhi, ylow, yhi
      real x(1), y(1)
c
      include 'gksdefs.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 l, m, npts, lseg, ix1, ix2, iy1, iy2, ival
      real seg(NSEGMAX,4), xpts(5), ypts(5),
     1     xtemp(NPOLYMAX), ytemp(NPOLYMAX)
      real cvalue, zmin, zmax, x1, x2, y1, y2, dx, dy,
     1     xpos, ypos, temp
      logical used(NSEGMAX)
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     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
c
         lseg = 0
c        Scan array, generating line segments
         do 1 l = ix1, ix2-1
            x1 = x(l)
            x2 = x(l+1)
            dx = (x2 - x1)
            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. cvalue .or.  cvalue .gt. zmax ) goto 2
c              Generate points
               npts = 0
               y1 = y(m)
               y2 = y(m+1)
               dy = y2 - y1
c              Lower x
               temp = array(l+1,m) - array(l,m)
               if( temp .eq. 0.0 ) then
                  if( array(l,m) .eq. cvalue ) then
                     npts = npts + 1
                     xpts(npts) = x2
                     ypts(npts) = y1
                  endif
               else
                  if( cvalue .gt. MIN(array(l,m),array(l+1,m)) .and.
     1                cvalue .le. MAX(array(l,m),array(l+1,m)) ) then
                     xpos = x1 + dx*(cvalue - array(l,m))/temp
                     npts = npts + 1
                     xpts(npts) = xpos
                     ypts(npts) = y1
                  endif
               endif
c              Right Y
               temp = array(l+1,m+1) - array(l+1,m)
               if( temp .eq. 0.0 ) then
                  if( array(l+1,m) .eq. cvalue ) then
                     npts = npts + 1
                     xpts(npts) = x2
                     ypts(npts) = y2
                  endif
               else
                  if( cvalue .gt. MIN(array(l+1,m),array(l+1,m+1)) .and.
     1                cvalue .le. MAX(array(l+1,m),array(l+1,m+1))) then
                     ypos = y1 + dy*(cvalue - array(l+1,m))/temp
                     npts = npts + 1
                     xpts(npts) = x2
                     ypts(npts) = ypos
                  endif
               endif
c              Upper x
               temp = array(l+1,m+1) - array(l,m+1)
               if( temp .eq. 0.0 ) then
                  if( array(l,m+1) .eq. cvalue ) then
                     npts = npts + 1
                     xpts(npts) = x1
                     ypts(npts) = y2
                  endif
               else
                  if( cvalue .gt. MIN(array(l,m+1),array(l+1,m+1)) .and.
     1                cvalue .le. MAX(array(l,m+1),array(l+1,m+1))) then
                     xpos = x1 + dx*(cvalue - array(l,m+1))/temp
                     npts = npts + 1
                     xpts(npts) = xpos
                     ypts(npts) = y2
                  endif
               endif
c              Left Y
               temp = array(l,m+1) - array(l,m)
               if( temp .eq. 0.0 ) then
                  if( array(l,m) .eq. cvalue ) then
                     npts = npts + 1
                     xpts(npts) = x1
                     ypts(npts) = y1
                  endif
               else
                  if( cvalue .gt. MIN(array(l,m),array(l,m+1)) .and.
     1                cvalue .le. MAX(array(l,m),array(l,m+1))) then
                     ypos = y1 + dy*(cvalue - array(l,m))/temp
                     npts = npts + 1
                     xpts(npts) = x1
                     ypts(npts) = ypos
                  endif
               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
c
2           continue
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 plot_contour_segments(seg, used, lseg, maxseg,
     1              x, y, maxpoly, do_y_flip, ysum, cvalue, do_label)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
c
c     Plot contour segments in order
c     If do_y_flip is .TRUE., then invert y coordinates using
c     ysum as maximum y
c
      integer lseg, maxseg, maxpoly
      real seg(maxseg,4), x(maxpoly), y(maxpoly), ysum, cvalue
      logical do_y_flip, used(maxseg), do_label
c
c     locals
c
      integer i, j, l, imid, ierr, ilen, string_len
      logical more_segs, found_matching_seg
      character*8 label
      real delx, dely, denom, xv, yv, world_win(4), world_vp(4),
     1     xmid, ymid, x1, x2, y1, y2, concat_x(4), concat_y(4),
     2     extent_x(4), extent_y(4)
c
      do i=1,lseg
         used(i) = .FALSE.
      end do
c
1000  continue
c        Look for a segment not plotted yet
         more_segs = .FALSE.
         do i=1,lseg
            if( .NOT. used(i) ) then
c              Found a segment not plotted
               more_segs = .TRUE.
               x(1) = seg(i,1)
               x(2) = seg(i,2)
               y(1) = seg(i,3)
               y(2) = seg(i,4)
               used(i) = .TRUE.
               l = 2
               goto 100
            endif
         end do
100      if( .NOT. more_segs ) goto 2000  ! exit
c        Now look for neighboring segments
         found_matching_seg = .TRUE.
         do while( found_matching_seg )
            found_matching_seg = .FALSE.
            do i=1,lseg
               if( .NOT. used(i) .AND. l .lt. maxpoly ) then
                  if(x(l) .EQ. seg(i,1) .AND. y(l) .EQ. seg(i,3)) then
c                    Append segment on end
                     l = l + 1
                     x(l) = seg(i,2)
                     y(l) = seg(i,4)
                     found_matching_seg = .TRUE.
                     used(i) = .TRUE.
                  else if(x(l) .EQ. seg(i,2) .AND.
     1                    y(l) .EQ. seg(i,4)) then
c                    Append segment on end
                     l = l + 1
                     x(l) = seg(i,1)
                     y(l) = seg(i,3)
                     found_matching_seg = .TRUE.
                     used(i) = .TRUE.
                  else if(x(1) .EQ. seg(i,2) .AND.
     1                    y(1) .EQ. seg(i,4)) then
c                    Prepend segment at top of list
                     do j=l,1,-1
                        x(j+1) = x(j)
                        y(j+1) = y(j)
                     end do
                     l = l+1
                     x(1) = seg(i,1)
                     y(1) = seg(i,3)
                     found_matching_seg = .TRUE.
                     used(i) = .TRUE.
                  else if(x(1) .EQ. seg(i,1) .AND.
     1                    y(1) .EQ. seg(i,3)) then
c                    Prepend segment at top of list
                     do j=l,1,-1
                        x(j+1) = x(j)
                        y(j+1) = y(j)
                     end do
                     l = l+1
                     x(1) = seg(i,2)
                     y(1) = seg(i,4)
                     found_matching_seg = .TRUE.
                     used(i) = .TRUE.
                  endif
               endif
            end do
         end do
c        FINALLY!!!!! Plot something!!!!
         if( do_y_flip ) then
            do i=1,l
               y(i) = (ysum - y(i))
            end do
         endif
         call gks$polyline(l,x,y)
c
c        Plot label on contour line if asked...
c
         if( do_label ) then
c           Select world coordinate xformation similar to NDC space
c           in plot2d to avoid text distortion when rotating text.
            call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_CENTER,
     1                              GKS$K_TEXT_VALIGN_HALF)
            call gks$inq_xform(1,ierr,world_win,world_vp)
            call gks$set_window(2, world_vp(1), world_vp(2),
     1                             world_vp(3), world_vp(4))
            call gks$set_viewport(2, world_vp(1), world_vp(2),
     1                               world_vp(3), world_vp(4))
            call gks$select_xform(2)
c
c           Set up character size for label and format
            call tune_charsize(0.01, 0.01, 0.0, 0.01, 0.01, 0.0)
            call format_contour_label(cvalue, label)
c
c           Select midpoint of contour line segments as place to
c           put label. NOTE that since line segments are not equal in
c           length that the label does NOT fall in the middle of the
c           contour line. This is good from an esthetic point of view.
            imid = l/2
            x1 = x(imid)
            x2 = x(imid+1)
            y1 = y(imid)
            y2 = y(imid+1)
            xmid = 0.5*(x1+x2)
            ymid = 0.5*(y1+y2)
c           Rotate text to be parallel to contour line segment
c           Xform to current wc space so angles comes out right
            call xform_wc_wc(x1, y1, world_win, world_vp)
            call xform_wc_wc(x2, y2, world_win, world_vp)
            if( x2 .gt. x1 ) then
               delx = x2 - x1
               dely = y2 - y1
            else
               delx = x1 - x2
               dely = y1 - y2
            endif
            denom = SQRT(delx**2 + dely**2)
            if( denom .eq. 0.0 ) then
               xv = 0.0
               yv = 1.0
            else
               xv = -dely/denom
               yv = delx/denom
            endif
            call GKS$SET_TEXT_UPVEC(xv,yv)
c
            call xform_wc_wc(xmid, ymid, world_win, world_vp)
c           Draw background box for label (only works on workstations
c           that can draw in the background color).
            ilen = string_len(label)
            call GKS$INQ_TEXT_EXTENT(ws_id, xmid, ymid, label(1:ilen),
     1           ierr, concat_x, concat_y, extent_x, extent_y)
            call GKS$SET_FILL_COLOR_INDEX(0)
            call GKS$SET_FILL_INT_STYLE(GKS$K_INTSTYLE_SOLID)
            call GKS$FILL_AREA(4, extent_x, extent_y)
            call GKS$TEXT(xmid,ymid,label(1:ilen))
c           Restore transformation
            call gks$select_xform(1)
         end if
      goto 1000
c
2000  continue
      call GKS$SET_TEXT_UPVEC(0.0,1.0)
      return
      end
c
      subroutine xform_wc_wc(x, y, wc1, wc2)
      implicit none
c
      real x, y, wc1(4), wc2(4)
c
c     Tranform coodinate (x,y) from world coordinate system specified
c     by wc1 to wc2. Wc1/2 are vectors of length four which
c     specify the world coordiante XMIN, XMAX, YMIN, YMAX (in that order)
c
c     Xform to 0..1 space
      x = (x-wc1(1))/(wc1(2)-wc1(1))
      y = (y-wc1(3))/(wc1(4)-wc1(3))
c     Xform to wc2 space
      x = wc2(1) + x*(wc2(2)-wc2(1))
      y = wc2(3) + y*(wc2(4)-wc2(3))
c
      return
      end
c
      subroutine format_contour_label(value, label)
      implicit none
c
      real value
      character*(*) label
c
      integer i,j
      character*8 temp
      real abs_value
c
      if( LEN(label) .lt. 8 ) then
         call plt$error_message('format_contour_label',
     1      'Internal error: Label string too small (< 8 chars)')
         label = '*'
         return
      endif
c
      abs_value = ABS(value)
      if( abs_value .le. 1.0E-6 ) then
c        1.0E-6 is FUZZ value used in contour plotting routine
c        (value used if contour to be plotted is exactly 0.0)
         temp = '0.0'
      else if( abs_value .ge. 10000.0 .or. abs_value .lt. 0.001 ) then
         write(temp,200) value
200      format(1pe8.0)
c        Remove redundant characters in exponent and mantissa
c        to make string as short as possible.
         if( temp(7:7) .eq. '0' ) then
            temp(7:7) = temp(8:8)
            temp(8:8) = ' '
         endif
         if( temp(6:6) .eq. '+' ) then
            temp(6:7) = temp(7:8)
            temp(8:8) = ' '
         endif
         if( temp(4:4) .eq. '.' ) then
            temp(4:7) = temp(5:8)
            temp(8:8) = ' '
         endif
      else if( abs_value .ge. 10.0 ) then
         write(temp,300) value
300      format(f8.0)
      else if( abs_value .ge. 0.1 ) then
         write(temp,400) value
400      format(f8.1)
      else if( abs_value .ge. 0.01 ) then
         write(temp,500) value
500      format(f8.2)
      else if( abs_value .ge. 0.001 ) then
         write(temp,600) value
600      format(f8.3)
      endif
c     Return formatted string minus any initial blanks
      j = 0
      do i=1,8
         j = j + 1
         if( temp(i:i) .ne. ' ' ) goto 700
      end do
700   continue
      label = temp(j:8)
c
      return
      end
