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 surf_plot_data_3d
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'surf.cmn'
      include 'defaults.cmn'
c
c     Defines for call to hidden line algorithm
      integer TOP, BOTTOM, BOTH
      parameter (TOP    = 1)
      parameter (BOTTOM = 2)
      parameter (BOTH   = 3)
c
      real, allocatable ::  pts_x(:),pts_y(:), pts_z(:)
      real  xloc(COMP_ARRAYSIZE), yloc(COMP_ARRAYSIZE)
      real  temp_x(5),temp_y(5), eye_equiv
      integer   i,j, k, ni, nj, dj, djp1, djm1
      integer   num_levels, istat, array_size
      real  level_values(10), min_value, max_value, dscale
      integer fill_color_index  
      logical do_separate, do_color_fill
      logical linvert_x, linvert_y, lswap_xy
c
      save array_size, pts_x, pts_y, pts_z
c
      if( .NOT. allocated(pts_x) ) then
         allocate(pts_x(1:surf_curve%npoints),stat=istat)
         allocate(pts_y(1:surf_curve%npoints),stat=istat)
         allocate(pts_z(1:surf_curve%npoints),stat=istat)
         array_size = surf_curve%npoints
         if( istat .ne. 0 ) then
            call display_message(
     1   'Postmini internal error: could not allocate temporary array')
            return
         endif
      else if( surf_curve%npoints .gt. array_size ) then
c        Reallocate larger arrays
         deallocate(pts_x,stat=istat)
         deallocate(pts_y,stat=istat)
         deallocate(pts_z,stat=istat)
         allocate(pts_x(1:surf_curve%npoints),stat=istat)
         allocate(pts_y(1:surf_curve%npoints),stat=istat)
         allocate(pts_z(1:surf_curve%npoints),stat=istat)
         array_size = surf_curve%npoints
         if( istat .ne. 0 ) then
            call display_message(
     1   'Postmini internal error: could not allocate temporary array')
            return
         endif
      endif
c
      call surf_tr_3d_2d(surf_curve%scaled_xposition,
     +                   surf_curve%scaled_yposition,
     +                   surf_curve%scaled_zposition,
     +                   pts_x, pts_y, pts_z, surf_curve%npoints)
c      
c     Perform different hidden line algorigthms, depending on what
c     user requested.
c      
c     AREA_FILL algorithm ("painter's algorithm"):
c      
c     Plot each mesh "square", first with area fill
c     to blot out anything behind it, then draw a polyline
c     to define the square.
c      
c     Advantages:
c        Fast (give fast area fill on the output device) and accurate.
c      
c     Disadvantages:
c        Doesn't work if device doesn't support area fill, or if
c        a filled primitive does not cover what is beneath it
c        (e.g. as on a pen plotter).
c      
c     HORIZON_FCN algorithm:
c      
c     Maintain a list of line segments which represents the upper and
c     lower horizons of the surface as currently plotted (line by line).
c     Horizons are updated as we go along.
c      
c     Advantages:
c       Reasonably fast, can plot bottom and top of surface
c       in different colors, works for all output devices.
c      
c     Disadvantages:
c        Can run out of storage space for very complicated shapes.
c        Can fail if a perpective transformation is used, because
c        horizon function can become multiple-valued.
c      
c     First, we inquire on the workstation.
c     If we are on a plotter (and we are using the horizon function
c     algorithm), then run through algorithm twice (once for top and
c     bottom) to avoid pen color changes.
c      
c     Color devices use color index 2 (usually red) for area fill
c     monochrome devices use the background color (index 0).
c
c     Test for pen plotters
      if( ws_type .ge. 51 .AND. ws_type .LE. 56 ) then
         do_separate = .TRUE.
      else
         do_separate = .FALSE.
      endif
c
c     Convert eye position to range of 0 to 360 degrees
      eye_equiv = AMOD(surf_eye%theta,360.0)
      if( eye_equiv  .LE. 0.0 ) eye_equiv = eye_equiv + 360.0
c      
      if( surf_hidden_line_method .eq. HIDDEN_LINES_WITH_HORIZON_FCN )
     +   then
c      
c        Shift data around so plotting always starts closest to eye position.
c        MUST be done for algorithm to work! Use pts_z array as temporary
c        storage.
c      
         if( eye_equiv .GT. 0 .AND. eye_equiv .LE. 90.0 ) then
               k = 1
               do i=surf_curve%ny,1,-1
                 do j=(surf_curve%nx-1),0,-1
                    pts_z(k) = pts_x(i + j*surf_curve%ny)
                    k = k + 1
                 end do
               end do
               do i=1,k
                  pts_x(i) = pts_z(i)
               end do
               k = 1
               do i=surf_curve%ny,1,-1
                 do j=(surf_curve%nx-1),0,-1
                    pts_z(k) = pts_y(i + j*surf_curve%ny)
                    k = k + 1
                 end do
               end do
               do i=1,k
                  pts_y(i) = pts_z(i)
               end do
               ni = surf_curve%nx
               nj = surf_curve%ny
         else if(eye_equiv .GT. 90.0 .AND. eye_equiv .LE. 180.0 ) then
               k = 1
               do j=0,(surf_curve%nx-1),1
                  do i=surf_curve%ny,1,-1
                    pts_z(k) = pts_x(i + j*surf_curve%ny)
                    k = k + 1
                 end do
               end do
               do i=1,k
                  pts_x(i) = pts_z(i)
               end do
               k = 1
               do j=0,(surf_curve%nx-1),1
                  do i=surf_curve%ny,1,-1
                    pts_z(k) = pts_y(i + j*surf_curve%ny)
                    k = k + 1
                 end do
               end do
               do i=1,k
                  pts_y(i) = pts_z(i)
               end do
               ni = surf_curve%ny
               nj = surf_curve%nx
         else if(eye_equiv .GT. 180.0 .AND. eye_equiv .LE. 270.0 ) then
               k = 1
               do i=1,surf_curve%ny
                 do j=0,(surf_curve%nx-1)
                    pts_z(k) = pts_x(i + j*surf_curve%ny)
                    k = k + 1
                 end do
               end do
               do i=1,k
                  pts_x(i) = pts_z(i)
               end do
               k = 1
               do i=1,surf_curve%ny
                 do j=0,(surf_curve%nx-1)
                    pts_z(k) = pts_y(i + j*surf_curve%ny)
                    k = k + 1
                 end do
               end do
               do i=1,k
                  pts_y(i) = pts_z(i)
               end do
               ni = surf_curve%nx
               nj = surf_curve%ny
         else 
               k = 1
               do j=(surf_curve%nx-1),0,-1
                  do i=1,surf_curve%ny
                    pts_z(k) = pts_x(i + j*surf_curve%ny)
                    k = k + 1
                 end do
               end do
               do i=1,k
                  pts_x(i) = pts_z(i)
               end do
               k = 1
               do j=(surf_curve%nx-1),0,-1
                  do i=1,surf_curve%ny
                    pts_z(k) = pts_y(i + j*surf_curve%ny)
                    k = k + 1
                 end do
               end do
               do i=1,k
                  pts_y(i) = pts_z(i)
               end do
               ni = surf_curve%ny
               nj = surf_curve%nx
         endif
         if( do_separate ) then
            call hidden(pts_x,pts_y,ni,nj,TOP)
            call hidden(pts_x,pts_y,ni,nj,BOTTOM)
         else
            call hidden(pts_x,pts_y,ni,nj,BOTH)
         endif
c
      else if(surf_hidden_line_method .eq. HIDDEN_LINES_WITH_FILL_AREA)
     +   then
c        Else do the area fill algorithm
c        If we have more that two colors, use multi-color contour fill.
c        Otherwise, use background color as fill pattern.
         if( surf_ws%ncolors .GT. 2 ) then
            do_color_fill = .TRUE.
c           Color map for contours set in surf subroutine, before we
c           write anything to screen. Avoids possible newpage when we
c           write annotations out.
            num_levels = 9 ! Fixed value for now
c           Compute default bucket values
            min_value =  PM_MAX_FLOAT
            max_value = -PM_MAX_FLOAT
            do i=1,surf_curve%npoints
              if( min_value .GT. surf_curve%zposition(i) )
     +            min_value = surf_curve%zposition(i)
              if( max_value .LT. surf_curve%zposition(i) )
     +            max_value = surf_curve%zposition(i)
            end do
c           User defined min/max color bands
            if( l_z_colormap_min ) then
               if( l_z_axis_log ) then
                  min_value = ALOG10(z_colormap_min)
               else
                  min_value = z_colormap_min
               endif
            endif
            if( l_z_colormap_max ) then
               if( l_z_axis_log ) then
                  max_value = ALOG10(z_colormap_max)
               else
                  max_value = z_colormap_max
               endif
            endif
c
c           Modification for Christian Schiebl, make color bars follow
c           Z axis min/max
            min_value = surf_zaxis%min
            max_value = surf_zaxis%max
c
            do i=1,num_levels
               level_values(i) = min_value + (max_value - min_value)*
     +                           (FLOAT(i)/(FLOAT(num_levels+1)))
            end do
         else
            num_levels = 1
            do_color_fill = .FALSE.
         endif
c
         fill_color_index = 0
         call gks$set_fill_color_index(fill_color_index)
         call gks$set_fill_int_style(GKS$K_INTSTYLE_SOLID)
c        Do eight different orderings, depending on what quadrant eye
c        is in. Try to make ordering such that we progress towards the
c        eye position.
         if( eye_equiv .GE. 0.0 .AND. eye_equiv .LE. 45.0 ) then
            do j=1,(surf_curve%nx - 1)
               dj = (j-1)*surf_curve%ny
               djp1 = j*surf_curve%ny
               do i=1,(surf_curve%ny - 1)
                  if( do_color_fill ) then
                     temp_x(1) = surf_curve%zposition( i+1  + dj)
                     temp_x(2) = surf_curve%zposition( i+1  + djp1)
                     temp_x(3) = surf_curve%zposition( i    + djp1)
                     temp_x(4) = surf_curve%zposition( i    + dj)
                     call fill_3d_contours(temp_x,
     +                       surf_curve%xposition( i+1  + dj),
     +                       surf_curve%xposition( i+1  + djp1),
     +                       surf_curve%yposition( i+1  + dj),
     +                       surf_curve%yposition( i    + djp1),
     +                       num_levels, level_values)
                  endif
                  temp_x(1) = pts_x(i + dj)
                  temp_x(5) = temp_x(1)
                  temp_y(1) = pts_y(i + dj)
                  temp_y(5) = temp_y(1)
                  temp_x(2) = pts_x( i+1  + dj)
                  temp_y(2) = pts_y( i+1  + dj)
                  temp_x(3) = pts_x( i+1  + djp1)
                  temp_y(3) = pts_y( i+1  + djp1)
                  temp_x(4) = pts_x(i + djp1)
                  temp_y(4) = pts_y(i + djp1)
                  if( .NOT. do_color_fill )
     +               call gks$fill_area(4,temp_x,temp_y)
                  call gks$polyline(5,temp_x,temp_y)
               end do
            end do
         else if( eye_equiv .GT. 45.0  .AND. eye_equiv .LE. 90.0 ) then
            do i=1,(surf_curve%ny - 1)
               do j=1,(surf_curve%nx - 1)
                  dj = (j-1)*surf_curve%ny
                  djp1 = j*surf_curve%ny
                  if( do_color_fill ) then
                  temp_x(1) = surf_curve%zposition( i    + dj)
                  temp_x(2) = surf_curve%zposition( i    + djp1)
                  temp_x(3) = surf_curve%zposition( i+1  + djp1)
                  temp_x(4) = surf_curve%zposition( i+1  + dj)
                  call fill_3d_contours(temp_x,
     +                             surf_curve%xposition( i    + dj),
     +                             surf_curve%xposition( i    + djp1),
     +                             surf_curve%yposition( i    + dj),
     +                             surf_curve%yposition( i+1  + djp1),
     +                             num_levels, level_values)
                  endif
                  temp_x(1) = pts_x(i + dj)
                  temp_y(1) = pts_y(i + dj)
                  temp_x(2) = pts_x(i + djp1)
                  temp_y(2) = pts_y(i + djp1)
                  temp_x(3) = pts_x( i+1  + djp1)
                  temp_y(3) = pts_y( i+1  + djp1)
                  temp_x(4) = pts_x( i+1  + dj)
                  temp_y(4) = pts_y( i+1  + dj)
                  temp_x(5) = pts_x(i + dj)
                  temp_y(5) = pts_y(i + dj)
                  if( .NOT. do_color_fill )
     +               call gks$fill_area(4,temp_x,temp_y)
                  call gks$polyline(5,temp_x,temp_y)
               end do
            end do
         else if( eye_equiv .GT. 90.0 .AND. eye_equiv .LE. 135.0 ) then
            do i=1,(surf_curve%ny - 1)
               do j=(surf_curve%nx - 1),1,-1
                  dj = j*surf_curve%ny
                  djm1 = (j-1)*surf_curve%ny
                  if( do_color_fill ) then
                  temp_x(1) = surf_curve%zposition( i    + dj)
                  temp_x(2) = surf_curve%zposition( i    + djm1)
                  temp_x(3) = surf_curve%zposition( i+1  + djm1)
                  temp_x(4) = surf_curve%zposition( i+1  + dj)
                  call fill_3d_contours(temp_x,
     +                             surf_curve%xposition( i    + dj),
     +                             surf_curve%xposition( i    + djm1),
     +                             surf_curve%yposition( i    + dj),
     +                             surf_curve%yposition( i+1  + djm1),
     +                             num_levels, level_values)
                  endif
                  temp_x(1) = pts_x(i + dj)
                  temp_y(1) = pts_y(i + dj)
                  temp_x(2) = pts_x(i + djm1)
                  temp_y(2) = pts_y(i + djm1)
                  temp_x(3) = pts_x( i+1  + djm1)
                  temp_y(3) = pts_y( i+1  + djm1)
                  temp_x(4) = pts_x( i+1  + dj)
                  temp_y(4) = pts_y( i+1  + dj)
                  temp_x(5) = pts_x(i + dj)
                  temp_y(5) = pts_y(i + dj)
                  if( .NOT. do_color_fill )
     +               call gks$fill_area(4,temp_x,temp_y)
                  call gks$polyline(5,temp_x,temp_y)
               end do
            end do
         else if( eye_equiv .GT. 135.0 .AND. eye_equiv .LE. 180.0 ) then
            do j=(surf_curve%nx - 1),1,-1
               dj = j*surf_curve%ny
               djm1 = (j-1)*surf_curve%ny
               do i=1,(surf_curve%ny - 1)
                  if( do_color_fill ) then
                  temp_x(1) = surf_curve%zposition( i+1  + dj)
                  temp_x(2) = surf_curve%zposition( i+1  + djm1)
                  temp_x(3) = surf_curve%zposition( i    + djm1)
                  temp_x(4) = surf_curve%zposition( i    + dj)
                  call fill_3d_contours(temp_x,
     +                             surf_curve%xposition( i+1  +  dj),
     +                             surf_curve%xposition( i+1  + djm1),
     +                             surf_curve%yposition( i+1  +  dj),
     +                             surf_curve%yposition( i    + djm1),
     +                             num_levels, level_values)
                  endif
                  temp_x(1) = pts_x(i + dj)
                  temp_y(1) = pts_y(i + dj)
                  temp_x(2) = pts_x( i+1  + dj)
                  temp_y(2) = pts_y( i+1  + dj)
                  temp_x(3) = pts_x( i+1  + djm1)
                  temp_y(3) = pts_y( i+1  + djm1)
                  temp_x(4) = pts_x(i + djm1)
                  temp_y(4) = pts_y(i + djm1)
                  temp_x(5) = pts_x(i + dj)
                  temp_y(5) = pts_y(i + dj)
                  if( .NOT. do_color_fill )
     +               call gks$fill_area(4,temp_x,temp_y)
                  call gks$polyline(5,temp_x,temp_y)
               end do
            end do
         else if( eye_equiv .GT. 180.0 .AND. eye_equiv .LE. 225.0 ) then
            do j=(surf_curve%nx - 1),1,-1
               dj = j*surf_curve%ny
               djm1 = (j-1)*surf_curve%ny
               do i=surf_curve%ny,2,-1
                  if( do_color_fill ) then
                  temp_x(1) = surf_curve%zposition( i-1  + dj)
                  temp_x(2) = surf_curve%zposition( i-1  + djm1)
                  temp_x(3) = surf_curve%zposition( i    + djm1)
                  temp_x(4) = surf_curve%zposition( i    + dj)
                  call fill_3d_contours(temp_x,
     +                             surf_curve%xposition( i-1  + dj),
     +                             surf_curve%xposition( i-1  + djm1),
     +                             surf_curve%yposition( i-1  + dj),
     +                             surf_curve%yposition( i    + djm1),
     +                             num_levels, level_values)
                  endif
                  temp_x(1) = pts_x(i + dj)
                  temp_y(1) = pts_y(i + dj)
                  temp_x(2) = pts_x( i-1  + dj)
                  temp_y(2) = pts_y( i-1  + dj)
                  temp_x(3) = pts_x( i-1  + djm1)
                  temp_y(3) = pts_y( i-1  + djm1)
                  temp_x(4) = pts_x(i + djm1)
                  temp_y(4) = pts_y(i + djm1)
                  temp_x(5) = pts_x(i + dj)
                  temp_y(5) = pts_y(i + dj)
                  if( .NOT. do_color_fill )
     +               call gks$fill_area(4,temp_x,temp_y)
                  call gks$polyline(5,temp_x,temp_y)
               end do
            end do
         else if( eye_equiv .GT. 225.0 .AND. eye_equiv .LE. 270.0 ) then
            do i=surf_curve%ny,2,-1
               do j=(surf_curve%nx - 1),1,-1
                  dj = j*surf_curve%ny
                  djm1 = (j-1)*surf_curve%ny
                  if( do_color_fill ) then
                  temp_x(1) = surf_curve%zposition( i    + dj)
                  temp_x(2) = surf_curve%zposition( i    + djm1)
                  temp_x(3) = surf_curve%zposition( i-1  + djm1)
                  temp_x(4) = surf_curve%zposition( i-1  + dj)
                  call fill_3d_contours(temp_x,
     +                             surf_curve%xposition( i    + dj),
     +                             surf_curve%xposition( i    + djm1),
     +                             surf_curve%yposition( i    + dj),
     +                             surf_curve%yposition( i-1  + djm1),
     +                             num_levels, level_values)
                  endif
                  temp_x(1) = pts_x(i + dj)
                  temp_y(1) = pts_y(i + dj)
                  temp_x(2) = pts_x(i + djm1)
                  temp_y(2) = pts_y(i + djm1)
                  temp_x(3) = pts_x( i-1  + djm1)
                  temp_y(3) = pts_y( i-1  + djm1)
                  temp_x(4) = pts_x( i-1  + dj)
                  temp_y(4) = pts_y( i-1  + dj)
                  temp_x(5) = pts_x(i + dj)
                  temp_y(5) = pts_y(i + dj)
                  if( .NOT. do_color_fill )
     +               call gks$fill_area(4,temp_x,temp_y)
                  call gks$polyline(5,temp_x,temp_y)
               end do
            end do
         else if( eye_equiv .GT. 270.0 .AND. eye_equiv .LE. 315.0 ) then
            do i=surf_curve%ny,2,-1
               do j=1,(surf_curve%nx - 1)
                  dj = (j-1)*surf_curve%ny
                  djp1 = j*surf_curve%ny
                  if( do_color_fill ) then
                  temp_x(1) = surf_curve%zposition(i +     dj)
                  temp_x(2) = surf_curve%zposition(i +     djp1)
                  temp_x(3) = surf_curve%zposition( i-1  + djp1)
                  temp_x(4) = surf_curve%zposition( i-1  + dj)
                  call fill_3d_contours(temp_x,
     +                             surf_curve%xposition( i    + dj),
     +                             surf_curve%xposition( i    + djp1),
     +                             surf_curve%yposition( i    + dj),
     +                             surf_curve%yposition( i-1  + djp1),
     +                             num_levels, level_values)
                  endif
                  temp_x(1) = pts_x(i + dj)
                  temp_y(1) = pts_y(i + dj)
                  temp_x(2) = pts_x(i + djp1)
                  temp_y(2) = pts_y(i + djp1)
                  temp_x(3) = pts_x( i-1  + djp1)
                  temp_y(3) = pts_y( i-1  + djp1)
                  temp_x(4) = pts_x( i-1  + dj)
                  temp_y(4) = pts_y( i-1  + dj)
                  temp_x(5) = pts_x(i + dj)
                  temp_y(5) = pts_y(i + dj)
                  if( .NOT. do_color_fill )
     +               call gks$fill_area(4,temp_x,temp_y)
                  call gks$polyline(5,temp_x,temp_y)
               end do
            end do
         else    ! Eye position in 315 to 360
            do j=1,(surf_curve%nx - 1)
               dj = (j-1)*surf_curve%ny
               djp1 = j*surf_curve%ny
               do i=surf_curve%ny,2,-1
                  if( do_color_fill ) then
                  temp_x(1) = surf_curve%zposition( i-1  + dj)
                  temp_x(2) = surf_curve%zposition( i-1  + djp1)
                  temp_x(3) = surf_curve%zposition( i    + djp1)
                  temp_x(4) = surf_curve%zposition( i    + dj)
                  call fill_3d_contours(temp_x,
     +                             surf_curve%xposition( i-1  + dj),
     +                             surf_curve%xposition( i-1  + djp1),
     +                             surf_curve%yposition( i-1  + dj),
     +                             surf_curve%yposition( i    + djp1),
     +                             num_levels, level_values)
                  endif
                  temp_x(1) = pts_x(i + dj)
                  temp_y(1) = pts_y(i + dj)
                  temp_x(2) = pts_x( i-1 + dj)
                  temp_y(2) = pts_y( i-1  + dj)
                  temp_x(3) = pts_x( i-1  + djp1)
                  temp_y(3) = pts_y( i-1  + djp1)
                  temp_x(4) = pts_x(i + djp1)
                  temp_y(4) = pts_y(i + djp1)
                  temp_x(5) = pts_x(i + dj)
                  temp_y(5) = pts_y(i + dj)
                  if( .NOT. do_color_fill )
     +               call gks$fill_area(4,temp_x,temp_y)
                  call gks$polyline(5,temp_x,temp_y)
               end do
            end do
         endif
c      
         if( surf_ws%ncolors .GT. 2 ) then
c           Plot color map legend (undo any data scaling)
            dscale = 10.0**surf_data_exp
            do i=1,num_levels
               level_values(i) = dscale*level_values(i)
            end do
            if( .NOT. use_stroke_prec ) then
c              Switch to hardware font for VTxxx to get right sizing
               call gks$set_text_fontprec(1,GKS$K_TEXT_PRECISION_STRING)
            endif
            call plot_color_map_legend(level_values, num_levels,
     1           max_value*dscale, .TRUE., .FALSE.)
         endif
      else if( surf_hidden_line_method .eq. HIDDEN_LINES_HISTOGRAM )
     +   then
c        Else do the histogram algorithm
c        If we have more that two colors, use multi-color contour fill.
c        Otherwise, use background color as fill pattern.
         if( surf_ws%ncolors .GT. 2 ) then
            do_color_fill = .TRUE.
c           Color map for contours set in surf subroutine, before we
c           write anything to screen. Avoids possible newpage when we
c           write annotations out.
            num_levels = 9 ! Fixed value for now
c           Compute default bucket values
            min_value =  PM_MAX_FLOAT
            max_value = -PM_MAX_FLOAT
            do i=1,surf_curve%npoints
              if( min_value .GT. surf_curve%zposition(i) )
     +            min_value = surf_curve%zposition(i)
               if( max_value .LT. surf_curve%zposition(i) )
     +            max_value = surf_curve%zposition(i)
            end do
c           User defined min/max color bands
            if( l_z_colormap_min ) then
               if( l_z_axis_log ) then
                  min_value = ALOG10(z_colormap_min)
               else
                  min_value = z_colormap_min
               endif
            endif
            if( l_z_colormap_max ) then
               if( l_z_axis_log ) then
                  max_value = ALOG10(z_colormap_max)
              else
                  max_value = z_colormap_max
               endif
            endif
c
c           Modification for Christian Schiebl, make color bars follow
c           Z axis min/max
            min_value = surf_zaxis%min
            max_value = surf_zaxis%max
c
            do i=1,num_levels
               level_values(i) = min_value + (max_value - min_value)*
     +                           (FLOAT(i)/(FLOAT(num_levels+1)))
            end do
         else
            num_levels = 1
            do_color_fill = .FALSE.
         endif
c
         fill_color_index = 0
         call gks$set_fill_color_index(fill_color_index)
         call gks$set_fill_int_style(GKS$K_INTSTYLE_SOLID)
c        Do four different orderings, depending on what quadrant eye
c        is in. Try to make ordering such that we progress from the furthest
c        position away from the eye position.
c
c        Get x and y grid locations (scaled coordinates)
c 
         do j=1,surf_curve%nx
            dj = 1 + (j-1)*surf_curve%ny
            xloc(j) = surf_curve%scaled_xposition(dj)
         end do
         do i=1,surf_curve%ny
            yloc(i) = surf_curve%scaled_yposition(i)
         end do
c
         if( eye_equiv .GE. 0.0 .AND. eye_equiv .LE. 90.0 ) then
            linvert_x = .FALSE.
            linvert_y = .FALSE.
            lswap_xy = .FALSE.
            call surf_plot_3d_histo(xloc, yloc,
     1           surf_curve%scaled_zposition,
     1           surf_curve%nx, surf_curve%ny,
     1           linvert_x, linvert_y, lswap_xy,
     1           surf_zaxis%min,surf_zaxis%max,level_values,num_levels,
     1           do_color_fill)
         else if( eye_equiv .GT. 90.0 .AND. eye_equiv .LE. 180.0 ) then
            linvert_x = .TRUE.
            linvert_y = .FALSE.
            lswap_xy = .FALSE.
            call surf_plot_3d_histo(xloc, yloc,
     1           surf_curve%scaled_zposition,
     1           surf_curve%nx, surf_curve%ny,
     1           linvert_x, linvert_y, lswap_xy,
     1           surf_zaxis%min,surf_zaxis%max,level_values,num_levels,
     1           do_color_fill)
         else if( eye_equiv .GT. 180.0 .AND. eye_equiv .LE. 270.0 ) then
            linvert_x = .TRUE.
            linvert_y = .TRUE.
            lswap_xy = .TRUE.
            call surf_plot_3d_histo(xloc, yloc,
     1           surf_curve%scaled_zposition,
     1           surf_curve%nx, surf_curve%ny,
     1           linvert_x, linvert_y, lswap_xy,
     1           surf_zaxis%min,surf_zaxis%max,level_values,num_levels,
     1           do_color_fill)
         else
            linvert_x = .FALSE.
            linvert_y = .TRUE.
            lswap_xy = .TRUE.
            call surf_plot_3d_histo(xloc, yloc,
     1           surf_curve%scaled_zposition,
     1           surf_curve%nx, surf_curve%ny,
     1           linvert_x, linvert_y, lswap_xy,
     1           surf_zaxis%min,surf_zaxis%max,level_values,num_levels,
     1           do_color_fill)
         endif
c      
         if( surf_ws%ncolors .GT. 2 ) then
c           Plot color map legend (undo any data scaling)
            dscale = 10.0**surf_data_exp
            do i=1,num_levels
               level_values(i) = dscale*level_values(i)
            end do
            if( .NOT. use_stroke_prec ) then
c              Switch to hardware font for VTxxx to get right sizing
               call gks$set_text_fontprec(1,GKS$K_TEXT_PRECISION_STRING)
            endif
            call plot_color_map_legend(level_values, num_levels,
     1              max_value*dscale, .TRUE., .FALSE.)
         endif
c
      endif
      
c     Reset defaults
      call gks$set_fill_color_index(1)
      call gks$set_fill_int_style(GKS$K_INTSTYLE_HOLLOW)
c
      return
      end
c
      subroutine surf_plot_tridata_3d
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'surf.cmn'
      include 'defaults.cmn'
c
c     Defines for call to hidden line algorithm
      integer TOP, BOTTOM, BOTH
      parameter (TOP    = 1)
      parameter (BOTTOM = 2)
      parameter (BOTH   = 3)
c
      real, allocatable :: pts_x(:),pts_y(:), pts_z(:)
      real  temp_x(5),temp_y(5), temp_z(5), eye_equiv
      integer   i,j, k, ie
      integer   num_levels, istat, array_size
      real  level_values(10), min_value, max_value, dscale, y_sum
      integer fill_color_index
      logical do_separate, do_color_fill
      logical linvert_x, linvert_y, lswap_xy, do_y_flip
      type eye_element
         integer i
         real    dist
      end type
      type (eye_element) , allocatable :: eye_order(:)
      real pm_eye_dist
      external pm_eye_compare
c
      save array_size, pts_x, pts_y, pts_z, eye_order
c
      if( .NOT. allocated(pts_x) ) then
         allocate(pts_x(1:surf_curve%npoints),stat=istat)
         allocate(pts_y(1:surf_curve%npoints),stat=istat)
         allocate(pts_z(1:surf_curve%npoints),stat=istat)
         allocate(eye_order(1:2*surf_curve%npoints),stat=istat)
         array_size = surf_curve%npoints
         if( istat .ne. 0 ) then
            call display_message(
     1   'Postmini internal error: could not allocate temporary array')
            return
         endif
      else if( surf_curve%npoints .gt. array_size ) then
c        Reallocate larger arrays
         deallocate(pts_x,stat=istat)
         deallocate(pts_y,stat=istat)
         deallocate(pts_z,stat=istat)
         deallocate(eye_order,stat=istat)
         allocate(pts_x(1:surf_curve%npoints),stat=istat)
         allocate(pts_y(1:surf_curve%npoints),stat=istat)
         allocate(pts_z(1:surf_curve%npoints),stat=istat)
         allocate(eye_order(1:2*surf_curve%npoints),stat=istat)
         array_size = surf_curve%npoints
         if( istat .ne. 0 ) then
            call display_message(
     1   'Postmini internal error: could not allocate temporary array')
            return
         endif
      endif
c
      call surf_tr_3d_2d(surf_curve%scaled_xposition,
     +                   surf_curve%scaled_yposition,
     +                   surf_curve%scaled_zposition,
     +                   pts_x, pts_y, pts_z, surf_curve%npoints)
c      
c     Test for pen plotters
      if( ws_type .ge. 51 .AND. ws_type .LE. 56 ) then
         do_separate = .TRUE.
      else
         do_separate = .FALSE.
      endif
c
c     Convert eye position to range of 0 to 360 degrees
      eye_equiv = AMOD(surf_eye%theta,360.0)
      if( eye_equiv  .LE. 0.0 ) eye_equiv = eye_equiv + 360.0
c
      if(surf_hidden_line_method .eq. HIDDEN_LINES_WITH_FILL_AREA .or.
     +   surf_hidden_line_method .eq. HIDDEN_LINES_WITH_HORIZON_FCN )
     +   then
c        Do the area fill algorithm
c        If we have more that two colors, use multi-color contour fill.
c        Otherwise, use background color as fill pattern.
         if( surf_ws%ncolors .GT. 2 ) then
            do_color_fill = .TRUE.
c           Color map for contours set in surf subroutine, before we
c           write anything to screen. Avoids possible newpage when we
c           write annotations out.
            num_levels = 9 ! Fixed value for now
c           Compute default bucket values
            min_value =  PM_MAX_FLOAT
            max_value = -PM_MAX_FLOAT
            do i=1,surf_curve%npoints
              if( min_value .GT. surf_curve%zposition(i) )
     +            min_value = surf_curve%zposition(i)
              if( max_value .LT. surf_curve%zposition(i) )
     +            max_value = surf_curve%zposition(i)
            end do
c           User defined min/max color bands
            if( l_z_colormap_min ) then
               if( l_z_axis_log ) then
                  min_value = ALOG10(z_colormap_min)
               else
                  min_value = z_colormap_min
               endif
            endif
            if( l_z_colormap_max ) then
               if( l_z_axis_log ) then
                  max_value = ALOG10(z_colormap_max)
               else
                  max_value = z_colormap_max
               endif
            endif
c
c           Modification for Christian Schiebl, make color bars follow
c           Z axis min/max
            min_value = surf_zaxis%min
            max_value = surf_zaxis%max
c
            do i=1,num_levels
               level_values(i) = min_value + (max_value - min_value)*
     +                           (FLOAT(i)/(FLOAT(num_levels+1)))
            end do
         else
            num_levels = 1
            do_color_fill = .FALSE.
         endif
c
         fill_color_index = 0
         call gks$set_fill_color_index(fill_color_index)
         call gks$set_fill_int_style(GKS$K_INTSTYLE_SOLID)
c
c        Sort by eye order here...
c        Use scaled coordinates because it looks like the eye position
c        is in a normalized coordinate system
         do i=1,surf_curve%nelem
            eye_order(i)%i    = i
c           Compute distance of element to eye position
            do j=1,3
               temp_x(j) =
     1            surf_curve%scaled_xposition(surf_curve%nop(j,i))
               temp_y(j) =
     1            surf_curve%scaled_yposition(surf_curve%nop(j,i))
            end do
            eye_order(i)%dist = pm_eye_dist(temp_x, temp_y,
     1         surf_eye%xposition, surf_eye%yposition)
         end do
         call pm_qsort(eye_order, surf_curve%nelem, 8, pm_eye_compare)
c
         y_sum = 0.0
         do_y_flip = .FALSE.
         do i=1,surf_curve%nelem
c
            if( .NOT. surf_curve%tri_clip(eye_order(i)%i) ) then
               ie = eye_order(i)%i
               temp_x(1) = pts_x(surf_curve%nop(1,ie))
               temp_x(4) = temp_x(1)
               temp_y(1) = pts_y(surf_curve%nop(1,ie))
               temp_y(4) = temp_y(1)
               temp_x(2) = pts_x(surf_curve%nop(2,ie))
               temp_y(2) = pts_y(surf_curve%nop(2,ie))
               temp_x(3) = pts_x(surf_curve%nop(3,ie))
               temp_y(3) = pts_y(surf_curve%nop(3,ie))
               if( do_color_fill ) then
                  temp_z(1) = surf_curve%zposition(surf_curve%nop(1,ie))
                  temp_z(2) = surf_curve%zposition(surf_curve%nop(2,ie))
                  temp_z(3) = surf_curve%zposition(surf_curve%nop(3,ie))
                  call pm_color_fill_contours_tri(temp_x, temp_y,
     1                    temp_z, level_values, num_levels,
     1                    max_cont_colors, y_sum, do_y_flip)
               else
                  call gks$fill_area(3,temp_x,temp_y)
               endif
               call gks$polyline(4,temp_x,temp_y)
            endif
         end do
c      
         if( surf_ws%ncolors .GT. 2 ) then
c           Plot color map legend (undo any data scaling)
            dscale = 10.0**surf_data_exp
            do i=1,num_levels
               level_values(i) = dscale*level_values(i)
            end do
            if( .NOT. use_stroke_prec ) then
c              Switch to hardware font for VTxxx to get right sizing
               call gks$set_text_fontprec(1,GKS$K_TEXT_PRECISION_STRING)
            endif
            call plot_color_map_legend(level_values, num_levels,
     1           max_value*dscale, .TRUE., .FALSE.)
         endif
      else if( hidden_line_method .eq. HIDDEN_LINES_HISTOGRAM )
     +   then
c        Else do the histogram algorithm
c        If we have more that two colors, use multi-color contour fill.
c        Otherwise, use background color as fill pattern.
         if( surf_ws%ncolors .GT. 2 ) then
            do_color_fill = .TRUE.
c           Color map for contours set in surf subroutine, before we
c           write anything to screen. Avoids possible newpage when we
c           write annotations out.
            num_levels = 9 ! Fixed value for now
c           Compute default bucket values
            min_value =  PM_MAX_FLOAT
            max_value = -PM_MAX_FLOAT
            do i=1,surf_curve%npoints
              if( min_value .GT. surf_curve%zposition(i) )
     +            min_value = surf_curve%zposition(i)
               if( max_value .LT. surf_curve%zposition(i) )
     +            max_value = surf_curve%zposition(i)
            end do
c           User defined min/max color bands
            if( l_z_colormap_min ) then
               if( l_z_axis_log ) then
                  min_value = ALOG10(z_colormap_min)
               else
                  min_value = z_colormap_min
               endif
            endif
            if( l_z_colormap_max ) then
               if( l_z_axis_log ) then
                  max_value = ALOG10(z_colormap_max)
              else
                  max_value = z_colormap_max
               endif
            endif
c
c           Modification for Christian Schiebl, make color bars follow
c           Z axis min/max
            min_value = surf_zaxis%min
            max_value = surf_zaxis%max
c
            do i=1,num_levels
               level_values(i) = min_value + (max_value - min_value)*
     +                           (FLOAT(i)/(FLOAT(num_levels+1)))
            end do
         else
            num_levels = 1
            do_color_fill = .FALSE.
         endif
c
         fill_color_index = 0
         call gks$set_fill_color_index(fill_color_index)
         call gks$set_fill_int_style(GKS$K_INTSTYLE_SOLID)
c
c        ... implement plot here ...
c      
         if( surf_ws%ncolors .GT. 2 ) then
c           Plot color map legend (undo any data scaling)
            dscale = 10.0**surf_data_exp
            do i=1,num_levels
               level_values(i) = dscale*level_values(i)
            end do
            if( .NOT. use_stroke_prec ) then
c              Switch to hardware font for VTxxx to get right sizing
               call gks$set_text_fontprec(1,GKS$K_TEXT_PRECISION_STRING)
            endif
            call plot_color_map_legend(level_values, num_levels,
     1              max_value*dscale, .TRUE., .FALSE.)
         endif
c
      endif
      
c     Reset defaults
      call gks$set_fill_color_index(1)
      call gks$set_fill_int_style(GKS$K_INTSTYLE_HOLLOW)
c
      return
      end
c      
      subroutine surf_scale_coords(x_temp, y_temp, z_temp, npts)
      implicit none
c
      include 'const.h'
      include 'surf.cmn'
c
      integer  npts
      real x_temp(npts), y_temp(npts), z_temp(npts)
c
      call surf_scale_array(0.0, 1.0, surf_xaxis%min, surf_xaxis%max, 
     +                      x_temp, x_temp, npts, surf_xaxis%log)
      call surf_scale_array(0.0, 1.0, surf_yaxis%min, surf_yaxis%max,
     +                      y_temp, y_temp, npts, surf_yaxis%log)
      call surf_scale_array(0.0, 1.0, surf_zaxis%min, surf_zaxis%max,
     +                      z_temp, z_temp, npts, surf_zaxis%log)
      return
      end

c
      subroutine surf_plot_3d_histo(xloc, yloc, zloc,
     1           nx, ny, linvert_x, linvert_y, lswap_xy,
     1           plot_min, plot_max, level_values, num_levels,
     1           do_color_fill)
      implicit none
c
      integer nx, ny, num_levels
      real xloc(nx), yloc(ny), zloc(nx*ny), level_values(num_levels),
     1     plot_min, plot_max
      logical linvert_x, linvert_y, lswap_xy, do_color_fill
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'io.cmn'
      include 'plot.cmn'
c
c     locals
c
      integer i, j, k, dj
      integer istart, iend, istep, iside, jstart, jend, jstep, jside
      integer ilev, offset
      real  xbox_3d(5), ybox_3d(5), zbox_3d(5),
     1      xtemp_3d(5), ytemp_3d(5), ztemp_3d(5),
     1      xbox_2d(5), ybox_2d(5), zbox_2d(5),
     1      levels(25)
c
c     Code ASSUMES that fill color has been set to background (index=0)
c     before code has been entered, so it does not have to set the fill
c     color index if color fill is not to be used.
c
      if( num_levels .gt. 25 ) then
         write(luttyo,100)
100      format(' ',
     1          '(surf_plot_3d_histo) internal error: too many levels')
         return
      endif
c
c     Define scaled level values to match scaled z value (0..1)
      do k=1,num_levels
         levels(k) = (level_values(k) - plot_min)/(plot_max - plot_min)
      end do
c
      if( linvert_x ) then
         jstart = nx
         jend   = 1
         jstep  = -1
      else
         jstart = 1
         jend   = nx
         jstep  = 1
      endif
c
      if( linvert_y ) then
         istart = ny
         iend   = 1
         istep  = -1
      else
         istart = 1
         iend   = ny
         istep  = 1
      endif
c
      if ( lswap_xy ) then
         do i=istart,iend,istep
c           Get box coordinates in scaled 3D space
            if( i .eq. 1 ) then
               ybox_3d(1) = yloc(i)
            else
               ybox_3d(1) = 0.5*(yloc(i-1) + yloc(i))
            endif
            if( i .eq. ny ) then
               ybox_3d(2) = yloc(i)
            else
               ybox_3d(2) = 0.5*(yloc(i) + yloc(i+1))
            endif
            ybox_3d(3) = ybox_3d(2)
            ybox_3d(4) = ybox_3d(1)
            ybox_3d(5) = ybox_3d(1)
c
            do j=jstart,jend,jstep
               dj = (j-1)*ny
               if( j .eq. 1 ) then
                  xbox_3d(1) = xloc(j)
               else
                  xbox_3d(1) = 0.5*(xloc(j) + xloc(j-1))
               endif
               xbox_3d(2) = xbox_3d(1)
               if( j .eq. nx ) then
                  xbox_3d(3) = xloc(j)
               else
                  xbox_3d(3) = 0.5*(xloc(j) + xloc(j+1))
               endif
               xbox_3d(4) = xbox_3d(3)
               xbox_3d(5) = xbox_3d(1)
c
               do k=1,5
                  zbox_3d(k) = zloc(i + dj)
               end do
c
               if( do_color_fill ) then
c                 Select fill color
                  if( zbox_3d(1) .le. levels(1) ) then
                     ilev = 1
                  else if( zbox_3d(1) .gt. levels(num_levels)) then
                     ilev = num_levels+1
                  else
                     ilev = 2
                     do k=1,num_levels-1
                        if( zbox_3d(1) .gt. levels(k) .and.
     1                      zbox_3d(1) .le. levels(k+1) ) then
                           goto 1000
                        endif
                       ilev = ilev + 1
                     end do
1000                 continue
                 endif
c                 Colors are numbered 2..max_cont_colors+1
c                 (reserve color 1 for axes)
                  offset = MAX(max_cont_colors - num_levels,1) - 2
                  call GKS$SET_FILL_COLOR_INDEX
     1                 (MOD(offset+ilev,max_cont_colors)+2)
               endif
c
c              Draw in the sides of the rectangular
c              prism (with apropriate fill to get hidden surfaces)..
c
               if(linvert_x) then
                  jside = 1
               else
                  jside = 3
               endif
               do k=1,5
                  xtemp_3d(k) = xbox_3d(jside)
               end do
c
               ytemp_3d(1) = ybox_3d(1)
               ytemp_3d(2) = ybox_3d(2)
               ytemp_3d(3) = ybox_3d(2)
               ytemp_3d(4) = ybox_3d(1)
               ytemp_3d(5) = ybox_3d(1)
c
               ztemp_3d(1) = zbox_3d(1)
               ztemp_3d(2) = zbox_3d(1)
               ztemp_3d(3) = 0.0
               ztemp_3d(4) = 0.0
               ztemp_3d(5) = zbox_3d(1)
c
               call surf_tr_3d_2d(xtemp_3d, ytemp_3d, ztemp_3d,
     1                xbox_2d, ybox_2d, zbox_2d, 5)
               call gks$fill_area(5,xbox_2d,ybox_2d)
               call gks$polyline(5,xbox_2d,ybox_2d)
c
               xtemp_3d(1) = xbox_3d(1)
               xtemp_3d(2) = xbox_3d(3)
               xtemp_3d(3) = xbox_3d(3)
               xtemp_3d(4) = xbox_3d(1)
               xtemp_3d(5) = xbox_3d(1)
c
               if( linvert_y ) then
                  iside =  1
               else
                  iside =  2
               endif
               do k=1,5
                  ytemp_3d(k) = ybox_3d(iside)
               end do
c
               call surf_tr_3d_2d(xtemp_3d, ytemp_3d, ztemp_3d,
     1                xbox_2d, ybox_2d, zbox_2d, 5)
               call gks$fill_area(5,xbox_2d,ybox_2d)
               call gks$polyline(5,xbox_2d,ybox_2d)
c
c              Draw top of rectangular prism
               call surf_tr_3d_2d(xbox_3d, ybox_3d, zbox_3d,
     1                xbox_2d, ybox_2d, zbox_2d, 5)

               call gks$fill_area(5,xbox_2d,ybox_2d)
               call gks$polyline(5,xbox_2d,ybox_2d)
            end do
         end do
      else
         do j=jstart,jend,jstep
            dj = (j-1)*ny
            if( j .eq. 1 ) then
               xbox_3d(1) = xloc(j)
            else
               xbox_3d(1) = 0.5*(xloc(j) + xloc(j-1))
            endif
            xbox_3d(2) = xbox_3d(1)
            if( j .eq. nx ) then
               xbox_3d(3) = xloc(j)
            else
               xbox_3d(3) = 0.5*(xloc(j) + xloc(j+1))
            endif
            xbox_3d(4) = xbox_3d(3)
            xbox_3d(5) = xbox_3d(1)
c
            do i=istart,iend,istep
c              Get box coordinates in scaled 3D space
               if( i .eq. 1 ) then
                  ybox_3d(1) = yloc(i)
               else
                  ybox_3d(1) = 0.5*(yloc(i-1) + yloc(i))
               endif
               if( i .eq. ny ) then
                  ybox_3d(2) = yloc(i)
               else
                  ybox_3d(2) = 0.5*(yloc(i) + yloc(i+1))
               endif
               ybox_3d(3) = ybox_3d(2)
               ybox_3d(4) = ybox_3d(1)
               ybox_3d(5) = ybox_3d(1)
c
               do k=1,5
                  zbox_3d(k) = zloc(i + dj)
               end do
c
               if( do_color_fill ) then
c                 Select fill color
                  if( zbox_3d(1) .le. levels(1) ) then
                     ilev = 1
                  else if( zbox_3d(1) .gt. levels(num_levels)) then
                     ilev = num_levels+1
                  else
                     ilev = 2
                     do k=1,num_levels-1
                        if( zbox_3d(1) .gt. levels(k) .and.
     1                      zbox_3d(1) .le. levels(k+1) ) then
                           goto 2000
                       endif
                        ilev = ilev + 1
                     end do
2000                 continue
                  endif
               endif
c
c              Colors are numbered 2..max_cont_colors+1
c              (reserve color 1 for axes)
               offset = MAX(max_cont_colors - num_levels,1) - 2
               call GKS$SET_FILL_COLOR_INDEX
     1              (MOD(offset+ilev,max_cont_colors)+2)
c
c              Draw in the sides of the rectangular
c              prism (with apropriate fill to get hidden surfaces)..
c
               if(linvert_x) then
                  jside = 1
               else
                  jside = 3
               endif
               do k=1,5
                  xtemp_3d(k) = xbox_3d(jside)
               end do
c
               ytemp_3d(1) = ybox_3d(1)
               ytemp_3d(2) = ybox_3d(2)
               ytemp_3d(3) = ybox_3d(2)
               ytemp_3d(4) = ybox_3d(1)
               ytemp_3d(5) = ybox_3d(1)
c
               ztemp_3d(1) = zbox_3d(1)
               ztemp_3d(2) = zbox_3d(1)
               ztemp_3d(3) = 0.0
               ztemp_3d(4) = 0.0
               ztemp_3d(5) = zbox_3d(1)
c
               call surf_tr_3d_2d(xtemp_3d, ytemp_3d, ztemp_3d,
     1                xbox_2d, ybox_2d, zbox_2d, 5)
               call gks$fill_area(5,xbox_2d,ybox_2d)
               call gks$polyline(5,xbox_2d,ybox_2d)
c
               xtemp_3d(1) = xbox_3d(1)
               xtemp_3d(2) = xbox_3d(3)
               xtemp_3d(3) = xbox_3d(3)
               xtemp_3d(4) = xbox_3d(1)
               xtemp_3d(5) = xbox_3d(1)
c
               if( linvert_y ) then
                  iside =  1
               else
                  iside =  2
               endif
               do k=1,5
                  ytemp_3d(k) = ybox_3d(iside)
               end do
c
               call surf_tr_3d_2d(xtemp_3d, ytemp_3d, ztemp_3d,
     1                xbox_2d, ybox_2d, zbox_2d, 5)
               call gks$fill_area(5,xbox_2d,ybox_2d)
               call gks$polyline(5,xbox_2d,ybox_2d)
c
c              Draw top of rectangular prism
               call surf_tr_3d_2d(xbox_3d, ybox_3d, zbox_3d,
     1                xbox_2d, ybox_2d, zbox_2d, 5)

               call gks$fill_area(5,xbox_2d,ybox_2d)
               call gks$polyline(5,xbox_2d,ybox_2d)
            end do
         end do
      endif
c
      end
c
      integer function pm_eye_compare(el1, el2)
      implicit none
c
      type eye_element
         integer i
         real    dist
      end type
      type (eye_element) el1, el2
c
      if( el1%dist .gt. el2%dist ) then
         pm_eye_compare = -1
      else if( el1%dist .lt. el2%dist ) then
         pm_eye_compare = 1
      else
         pm_eye_compare = 0
      endif
c
      return
      end
c
      real function pm_eye_dist(x, y, eye_xpos, eye_ypos)
      implicit none
c
      real x(3), y(3), eye_xpos, eye_ypos
c
      integer i
      real tmp(3)
c
      do i=1,3
        tmp(i) = sqrt( (eye_xpos - x(i))**2 + (eye_ypos - y(i))**2 )
      end do
c
      pm_eye_dist = MIN(tmp(1), tmp(2), tmp(3))
c
      return
      end

