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
c    Subroutine to get data value interactively from a 2D plot
c
      subroutine get_probe_location(window_id,
     1           h_win_min, h_win_max, v_win_min, v_win_max,
     2           x1_ndc, x2_ndc, y1_ndc, y2_ndc, linvert,
     3           x1, y1, do_probe)
      implicit none
c
      integer window_id
      real h_win_min, h_win_max, v_win_min, v_win_max
      real x1_ndc, x2_ndc, y1_ndc, y2_ndc
      real x1, y1
      logical linvert, do_probe
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
c
      integer lookup_ws_id, work_ws_type
      integer work_id,
     1        xform, pet_type, dev_num, data_record(2),
     2        record_buffer_length, input_status, err_stat
      real    echo_area(4), box_x(2), box_y(2)
c
      if( .NOT. gks_is_open ) return
c
      work_id = lookup_ws_id(window_id)
      if( work_id .eq. PLT$EMPTY_WINDOW_SLOT ) then
c        Bad window id
         call plt$error_message('get_probe_location','Bad window id')
         do_probe = .FALSE.
         return
      endif
c
      work_ws_type = ws_struct(work_id)%type
      if( work_ws_type .ne. GKS$K_VSII .and.
     1    work_ws_type .ne. GKS$K_VSII_GPX .and.
     2    work_ws_type .ne. GKS$K_DECWINDOWS .and.
     3    work_ws_type .ne. GKS$K_MOTIF ) then
         call plt$error_message('get_probe_location',
     1        'Not a workstation')
         do_probe = .FALSE.
         return
      endif
c
      call activate_device_for_zoom(window_id,
     1       h_win_min, h_win_max, v_win_min, v_win_max,
     2       x1_ndc, x2_ndc, y1_ndc, y2_ndc,
     3       echo_area, err_stat)
c
c     Put locator at initial point
c     Use echo type 2 to get cross hairs
      dev_num = 1
      xform = 1
      pet_type = 2
c     gks$init_locator modifies record_buffer_length.
c     Reset it to expected value.
      record_buffer_length = 0
      call gks$init_locator(work_id, dev_num, 0.5, 0.5,
     1     xform, pet_type, echo_area,
     2     data_record, record_buffer_length)
      call gks$request_locator(work_id, dev_num, input_status,
     1     xform, box_x(1), box_y(1))
      if( input_status .ne. 1 ) then
c        Input was cancelled
         err_stat = PLT$INPUT_CANCELLED
         goto 1000
      endif
c
1000  continue
c
c     Deactivate window
      call deactivate_plot_device
      if( err_stat .eq. PLT$INPUT_CANCELLED ) then
         do_probe = .FALSE.
      else
         x1 = box_x(1)
         if( linvert ) then
            y1 = (v_win_max+v_win_min) - box_y(1)
         else
            y1 = box_y(1)
         endif
c
c        Limit to current bounds
         if( x1 .lt. h_win_min ) x1 = h_win_min
         if( x1 .gt. h_win_max ) x1 = h_win_max
         if( y1 .lt. v_win_min ) y1 = v_win_min
         if( y1 .gt. v_win_max ) y1 = v_win_max
         do_probe = .TRUE.
      endif
c
      return
      end
c
      subroutine probe_interpolate(data_array,
     1           xp, yp, do_log, data_is_trigrid, probe_value, ierr)
      implicit none
c
      include 'const.h'
      include 'mesh.cmn'
c
      integer ierr
      real data_array(ni,nj)
      real xp, yp, probe_value
      logical do_log, data_is_trigrid
c
c     Find value in data_array interpolated to position (x1,y1)
c     and return. Use log interpolation if specified.
c
c     ierr returns 0 on success, 1 on failure (coordinates out of
c     bounds)
c
      integer i, j, getindex
      real    x1, x2, y1, y2, delx, dely, h1, h2
c
      if( data_is_trigrid ) then
         call probe_interp_tri(data_array, tri_npts, tri_ne,
     1                         cord, nop,
     1                         xp, yp, do_log, probe_value, ierr)
      else
c        Find box which surrounds value (xp,yp)
         i = getindex(xp, xyz(1,idisti), ni)
         if( i .eq. 0 ) then
            ierr = 1
            goto 1000
         endif
c        Test for xp at edge of mesh
         if( i + 1 .gt. ni ) then
c           Back up one mesh increment
            i = i - 1
         endif
c
         j = getindex(yp, xyz(1,idistj), nj)
         if( j .eq. 0 ) then
            ierr = 1
            goto 1000
         endif
c        Test for yp at edge of mesh
         if( j + 1 .gt. nj ) then
c           Back up one mesh increment
            j = j - 1
         endif
c
c        Now interpolate value (use bilinear interpolate for now)
         x1 = xyz(i  ,idisti)
         x2 = xyz(i+1,idisti)
         delx = x2 - x1
         y1 = xyz(j  ,idistj)
         y2 = xyz(j+1,idistj)
         dely = y2 - y1
         h1 = (xp - x1)/(x2 - x1)
         h2 = (yp - y1)/(y2 - y1)
         probe_value = data_array(i,j) +
     1              h1*(data_array(i+1,j) - data_array(i,j)) +
     2              h2*(data_array(i,j+1) - data_array(i,j)) +
     3              h1*h2*(data_array(i+1,j+1) - data_array(i+1,j)
     4                     - data_array(i,j+1) + data_array(i,j) )
         if( do_log ) then
c           When logarithmic interpolation is specified, the data array
c           is assumed to have our 'special' log10 applied to it already.
c           Expand value back to true scale.
            probe_value = SIGN(10.0**MIN(38.0,ABS(probe_value)),
     1                         probe_value)
         endif
         ierr = 0
      endif
c
1000  continue
      end
c
      subroutine probe_interp_tri(data_array, npts, tri_ne,
     1           cord, nop, xp, yp, do_log, probe_value, ierr)
      implicit none
c
      integer ierr, npts, tri_ne, nop(3,tri_ne)
      real data_array(npts), cord(2,npts)
      real xp, yp, probe_value
      logical do_log
c
c     locals
      integer ie, j, inext(3)
      real x(3), y(3), z(3), d(3), s, t, dy1, dy2, x1, x2, z1, z2,
     1     a, b, c
      logical lfound
      data inext/2,3,1/
c
      lfound = .FALSE.
c     Find element which surrounds value (xp,yp)
      do 100 ie=1,tri_ne
         do j=1,3
            x(j) = cord(1,nop(j,ie))
            y(j) = cord(2,nop(j,ie))
         end do
c        Is point inside this triangle?
c        Do quick bounding box test first
         if( xp .gt. MAX(x(1),x(2),x(3)) .or.
     1       xp .lt. MIN(x(1),x(2),x(3)) .or.
     2       yp .gt. MAX(y(1),y(2),y(3)) .or.
     3       yp .lt. MIN(y(1),y(2),y(3)) ) goto 100
c        Full check: (xp,yp) inside triangle if it always lies
c        to the left (right) of each triangle segment.
         do j=1,3
            if( x(j) .eq. x(inext(j)) ) then
c              Vertical edge doesn't need to be checked,
c              since any point outside is caught by earlier bounding
c              box check. Nicely avoids zero divide.
               d(j) = 0.0
            else
c              Vertical distance between test point and edge
               t = (x(inext(j)) - x(j))
               s = (y(inext(j)) - y(j))/t
               d(j) = yp - (s*(xp-x(j)) + y(j))
c              Flip sign if x decreasing
               if( t .lt. 0.0 ) d(j) = -d(j)
            endif
         end do
c        Check for either all + or all -, so we can handle any
c        definition of triangle (clockwise or counter-clockwise)
         if((d(1).ge.0.0 .and. d(2).ge.0.0 .and. d(3).ge.0.0) .or.
     1      (d(1).le.0.0 .and. d(2).le.0.0 .and. d(3).le.0.0)) then
            lfound = .TRUE.
            goto 101
         endif
c         
100   continue
101   continue
c
      if( lfound ) then
c        Interpolate on plane (coded linear for now)
c        Code assumes that the triangle is 'well behaved' -- that is --
c        each point has unique coordinates. Should be no problem.
         z(1) = data_array(nop(1,ie))
         z(2) = data_array(nop(2,ie))
         z(3) = data_array(nop(3,ie))
         dy1 = y(3) - y(2)
         dy2 = y(2) - y(1)
         if( dy1 .eq. 0.0 ) then
            a = (z(3) - z(2))/(x(3) - x(2))
            b = ((z(2) - z(1)) - a*(x(2) - x(1)) )/dy2
         elseif ( dy2 .eq. 0.0 ) then
            a = (z(2) - z(1))/(x(2) - x(1))
            b = ((z(3) - z(2)) - a*(x(3) - x(2)) )/dy1
         else
            x1 = (x(3) - x(2))/dy1
            x2 = (x(2) - x(1))/dy2
            z1 = (z(3) - z(2))/dy1
            z2 = (z(2) - z(1))/dy2
            if( (x1 - x2) .eq. 0.0 ) then
               a  = 0.0
            else
               a  = (z1 - z2)/(x1 - x2)
            endif
            b  = z1 - a*x1
         endif
         c  = z(1) - a*x(1) - b*y(1)
         probe_value = a*xp + b*yp + c
         if( do_log ) then
c           When logarithmic interpolation is specified, the data array
c           is assumed to have our 'special' log10 applied to it already.
c           Expand value back to true scale.
            probe_value = SIGN(10.0**MIN(38.0,ABS(probe_value)),
     1                         probe_value)
         endif
         ierr = 0
      else
         ierr = 1
         probe_value = 0.0
      endif
c
      return
      end
c
      subroutine probe_multi_interpolate(dataset, xp, yp,
     1                                   probe_value, ierr)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'compare_2d.cmn'
c
      integer dataset, ierr
      real xp, yp, probe_value
c
c     Find value in dataset "dataset" interpolated to position (x1,y1)
c     and return.
c
c     ierr returns 0 on success, 1 on coordinates out of bounds,
c     -1 on other failure
c
      integer slot, i, j, nx, ny, joff, joff1, getindex
      real    x1, x2, y1, y2, delx, dely, h1, h2
      logical do_log
c
      slot = slot_index_2d(dataset)
      if( slot .eq. -1 ) then
         ierr = -1
         return
      endif
c
      if( multi_2d(slot)%data_is_trigrid ) then
         do_log = .FALSE.
         call probe_interp_tri(multi_2d(slot)%data,
     1        multi_2d(slot)%tri_npts, multi_2d(slot)%tri_ne,
     1        multi_2d(slot)%cord, multi_2d(slot)%nop,
     1        xp, yp, do_log, probe_value, ierr)
      else
         nx = multi_2d(slot)%nx
         ny = multi_2d(slot)%ny
c        Find box which surrounds value (xp,yp)
         i = getindex(xp, multi_2d(slot)%x, nx)
         if( i .eq. 0 ) then
            ierr = 1
            goto 1000
         endif
c        Test for xp at edge of mesh
         if( i + 1 .gt. nx ) then
c           Back up one mesh increment
            i = i - 1
         endif
c
         j = getindex(yp, multi_2d(slot)%y, ny)
         if( j .eq. 0 ) then
            ierr = 1
            goto 1000
         endif
c        Test for yp at edge of mesh
         if( j + 1 .gt. ny ) then
c           Back up one mesh increment
            j = j - 1
         endif
c
c        Now interpolate value (use bilinear interpolate for now)
         x1 = multi_2d(slot)%x(i)
         x2 = multi_2d(slot)%x(i+1)
         delx = x2 - x1
         y1 = multi_2d(slot)%y(j)
         y2 = multi_2d(slot)%y(j+1)
         dely = y2 - y1
         h1 = (xp - x1)/(x2 - x1)
         h2 = (yp - y1)/(y2 - y1)
         joff  = nx*(j-1)
         joff1 = nx*j
         probe_value =  multi_2d(slot)%data(i+joff) +
     1              h1*(multi_2d(slot)%data(i+1+joff) -
     1                  multi_2d(slot)%data(i+joff)) +
     2              h2*(multi_2d(slot)%data(i+joff1) -
     2                  multi_2d(slot)%data(i+joff)) +
     3              h1*h2*(multi_2d(slot)%data(i+1+joff1) -
     3                     multi_2d(slot)%data(i+1+joff) -
     4                     multi_2d(slot)%data(i+joff1) +
     4                     multi_2d(slot)%data(i+joff) )
         ierr = 0
      endif
c
1000  continue
      end
