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 fill_3d_contours(z, x1, x2, y1, y2, nval, val)
      implicit none
c
      integer nval
      real z(4), x1, x2, y1, y2, val(nval)
c
c     Plots color contours in a mesh square, using area fill.
c     Coordinates are tranformed from 3D to 2D space before plotting.
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
c
      integer ival, npts, offset
      real zmin, zmax, cvalue1, cvalue2
      real x_temp(8), y_temp(8), z_temp(8), xtr(8), ytr(8), ztr(8)
      real FUZZ
      parameter (FUZZ = 1.0E-6)

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        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
c        Skip to bottom of loop if cvalue outside of data
         zmin = MIN(z(1),z(2),z(3),z(4))
         zmax = MAX(z(1),z(2),z(3),z(4))
         if( zmin .gt. cvalue2 .or.  cvalue1 .gt. zmax ) goto 3
c
c        If box is totally within a color band, then output and exit
         if( zmin .ge. cvalue1 .and. zmax .le. cvalue2 ) then
            x_temp(1) = x1
            x_temp(2) = x2
            x_temp(3) = x2
            x_temp(4) = x1
            y_temp(1) = y1
            y_temp(2) = y1
            y_temp(3) = y2
            y_temp(4) = y2
            z_temp(1) = z(1)
            z_temp(2) = z(2)
            z_temp(3) = z(3)
            z_temp(4) = z(4)
c           Scale to 0.0,1.0
            call surf_scale_coords(x_temp, y_temp, z_temp, 4)
c           Transform 3d to 2d
            call surf_tr_3d_2d(x_temp, y_temp, z_temp, xtr, ytr, ztr, 4)
            call GKS$FILL_AREA(4, xtr, ytr)
            return
         endif
c        Construct partial box fill area
         npts = 0
c        Lower X
         call find_band_3d(x1, x2, y1, z(1), z(2),
     1                  cvalue1, cvalue2, x_temp, y_temp, z_temp, npts)
c        Right Y
         call find_band_3d(y1, y2, x2, z(2), z(3),
     1                  cvalue1, cvalue2, y_temp, x_temp, z_temp, npts)
c        Upper X
         call find_band_3d(x2, x1, y2, z(3), z(4),
     1                  cvalue1, cvalue2, x_temp, y_temp, z_temp, npts)
c        Left Y
         call find_band_3d(y2, y1, x1, z(4), z(1),
     1                  cvalue1, cvalue2, y_temp, x_temp, z_temp, npts)
c        This should "never" happend, but test for it anyway
         if( npts .le. 0 .or. npts .gt. 8 ) then
            call plt$error_message('fill_3d_contours',
     1         'Internal error in color band generation')
         else
c           Emit partial box fill area
c           Scale to 0.0,1.0
            call surf_scale_coords(x_temp, y_temp, z_temp, npts)
c           Transform 3d to 2d
            call surf_tr_3d_2d(x_temp, y_temp, z_temp,
     1                         xtr, ytr, ztr, npts)
            call GKS$FILL_AREA(npts, xtr, ytr)
         endif
c        Next cvalue
3     continue
c
      return
      end
c
      subroutine find_band_3d(u1, u2, v, f1, f2,
     1           cvalue1, cvalue2, coord1, coord2, coordz, npts)
      implicit none
c
      integer npts
      real u1, u2, v, f1, f2, cvalue1, cvalue2
      real coord1(8), coord2(8), coordz(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, z coordinate in coordz
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 no intersection 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
            coordz(npts) = cvalue2
         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
               coordz(npts) = f2
            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
                  coordz(npts) = cvalue1
               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
            coordz(npts) = cvalue1
         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
               coordz(npts) = f2
            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
                  coordz(npts) = cvalue2
               endif
            endif
         endif
      else
c        Point is inbetween contour values. Include it.
         npts = npts + 1
         coord1(npts) = u1
         coord2(npts) = v
         coordz(npts) = f1
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
               coordz(npts) = cvalue2
            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
               coordz(npts) = cvalue1
            endif
         else
            npts = npts + 1
            coord1(npts) = u2
            coord2(npts) = v
            coordz(npts) = f2
         endif
c        
      endif
c
      return
      end
c
      subroutine plot_color_map_legend(contour_values, num_contours,
     1           rmax, do_fill, do_solid)
      implicit none
c
      integer num_contours
      real    contour_values(num_contours), rmax
      logical do_fill, do_solid
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
c
      integer i, offset, ilen, ierr
      real xbox(5), ybox(5)
      character*50 num_label
c
c     draw legend on side...
      call GKS$SET_WINDOW(1, 0.0, 0.15, 0.1, 0.8)
      call GKS$SET_VIEWPORT(1, 0.825*ratio_x, 1.0*ratio_x,
     1                         0.1*ratio_y, 0.80*ratio_y)
      call GKS$SET_CLIPPING(GKS$K_NOCLIP)
      call GKS$SELECT_XFORM(1)
      call tune_charsize(0.01, 0.04, -0.1, 0.015, 0.025, 0.0)
      call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_NORMAL,
     1                        GKS$K_TEXT_VALIGN_HALF)
c
      if( do_fill ) then
c        Plot color map legend.
c        Colors are numbered 2..max_cont_colors+1 (reserve color 1 for axes)
         offset = MAX(max_cont_colors - num_contours,1) - 2
         do 50 i=1,num_contours+1
            call GKS$SET_FILL_COLOR_INDEX(
     1              MOD(offset+i,max_cont_colors)+2)
            xbox(1) = 0.02
            xbox(2) = 0.04
            xbox(3) = 0.04
            xbox(4) = 0.02
            ybox(1) = 0.8  - 0.06*FLOAT(i)
            ybox(2) = ybox(1)
            ybox(3) = 0.86 - 0.06*FLOAT(i)
            ybox(4) = ybox(3)
            call GKS$FILL_AREA(4, xbox, ybox)
c
            if( i .le. num_contours ) then
CC               write(num_label,2000) contour_values(i)
CC               call GKS$TEXT(0.05,  0.8 - 0.06*FLOAT(i),num_label)
               call format_mathmode_sci(contour_values(i),num_label,
     1                                  ilen)
               call mathmode(ws_id, 0.05,  0.8 - 0.06*FLOAT(i),
     1                       num_label(1:ilen), .FALSE., ierr)
            else
c              Check for case where user gave maximum contour
c              value larger than any data value (rmax)
               if( rmax .gt. contour_values(num_contours) ) then
CC                  write(num_label,2000) rmax
CC                  call GKS$TEXT(0.05,  0.8 - 0.06*FLOAT(i),num_label)
                  call format_mathmode_sci(rmax,num_label,ilen)
                  call mathmode(ws_id,0.05,  0.8 - 0.06*FLOAT(i),
     1                          num_label(1:ilen), .FALSE., ierr)
               endif
            endif
CC2000        format(' ',1pe9.2)
50       continue
c        Draw box around color legend
         call GKS$SET_PLINE_COLOR_INDEX(1)
         xbox(1) = 0.02
         xbox(2) = 0.04
         xbox(3) = 0.04
         xbox(4) = 0.02
         xbox(5) = 0.02
c        Lowest Y coordinate of legend is in ybox(1),ybox(2)
         ybox(3) = 0.80
         ybox(4) = ybox(3)
         ybox(5) = ybox(1)
         call GKS$POLYLINE(5, xbox, ybox)
         call GKS$SET_TEXT_ALIGN(GKS$K_TEXT_HALIGN_NORMAL,
     1                           GKS$K_TEXT_VALIGN_NORMAL)
c
      else
c        Do simple poly line legend
         do 60 i=1,num_contours
            if( do_solid ) then
               call set_line_index(1)
            else
               call set_line_index(i)
            endif
            call move_abs(0.02,  0.8 - 0.06*FLOAT(i))
            call draw_abs(0.04,  0.8 - 0.06*FLOAT(i))
60       continue
         do 61 i=1,num_contours
CC            write(num_label,2003) contour_values(i)
CC2003        format(1pe9.2)
CC            call GKS$TEXT(0.05,  0.8 - 0.06*FLOAT(i),num_label)
             call format_mathmode_sci(contour_values(i),num_label,
     1                                ilen)
             call mathmode(ws_id, 0.05,  0.8 - 0.06*FLOAT(i),
     1                     num_label(1:ilen), .FALSE., ierr)
61       continue
      endif
c
      return
      end
