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_box_gks
      implicit none
c
      include 'const.h'
      include 'gksdefs.cmn'
      include 'surf.cmn'
c
      integer xform_number
      integer n
      real x(5), y(5)
c
      integer gks$select_xform, istatus
      real THIRD
      parameter (THIRD = 1.0/3.0)
      real TWO_THIRD
      parameter (TWO_THIRD = 2.0/3.0)
c
      data n/5/
      data x/THIRD, 1., 1., THIRD, THIRD/
      data y/0., 0.,TWO_THIRD, TWO_THIRD, 0./
c      
c     Change transformation to the plot area.
c      
      xform_number = BOX_XFORM
c      
      istatus = gks$select_xform(xform_number)    
      call gks$set_clipping(GKS$K_NOCLIP)
      call gks$polyline(n,x,y)
      call gks$set_clipping(GKS$K_CLIP)
c
      return
      end
c      
      subroutine surf_set_plotwindow_gks
      implicit none
c
c     surf_set_plotwindow_gks performs the plot transformation
c
      include 'const.h'
      include 'surf.cmn'
c
      real    window(4)
      real    viewport(4)
      integer status, surf_xform_gks
      integer xform_number
c
      data window/0.,1.0,0.0,0.9/
      data xform_number/PLOT_XFORM/
c
c     First we set the window transformation to map to the plot area.
c      
      viewport(1) = window(1)*surf_global_xfactor + surf_global_xoffset
      viewport(2) = window(2)*surf_global_xfactor + surf_global_xoffset
      viewport(3) = window(3)*surf_ws%ratio*surf_global_yfactor +
     +                        surf_global_yoffset
      viewport(4) = window(4)*surf_ws%ratio*surf_global_yfactor +
     +                        surf_global_yoffset
      status = surf_xform_gks(xform_number,window,viewport)
      if(status .ne. SURF_NORMAL) 
     +   call surf_error('INVALID TRANSFORMATION FOR PLOT AREA',
     +                   status,1)
      return
      end      
c      
      subroutine surf_title_gks
      implicit none
c
      include 'const.h'
      include 'gksdefs.cmn'
      include 'surf.cmn'
c
      real   xposition, yposition(3)
      real   title_window(4), title_viewport(4)
      integer i, status, xform_number, surf_xform_gks
      real    char_height
      logical no_action, err
c
      data xposition/0.5/
      data yposition/0.965, 0.935, 0.905/
      data title_window/0.,1.0,0.9,1.0/
      data char_height/0.020/
      data xform_number/TITLE_XFORM/
c      
c     First we set the window transformation to map to the title area
c     of the plot.
c
      call gks$set_text_color_index(1)
      title_viewport(1) = title_window(1)*surf_global_xfactor +
     +                            surf_global_xoffset
      title_viewport(2) = title_window(2)*surf_global_xfactor +
     +                            surf_global_xoffset
      title_viewport(3) = title_window(3)*surf_ws%ratio*
     +                    surf_global_yfactor + surf_global_yoffset
      title_viewport(4) = title_window(4)*surf_ws%ratio*
     +                    surf_global_yfactor + surf_global_yoffset
      status = surf_xform_gks(xform_number,title_window,title_viewport)
      if (status .ne. SURF_NORMAL) 
     +    call surf_error('INVALID TRANSFORMATION FOR TITLE',status,1)
      
c          
      call gks$set_text_height(char_height)
      call gks$set_text_fontprec(surf_ws%font,
     +                           surf_ws%precision)
      call gks$set_text_align(GKS$K_TEXT_HALIGN_CENTER,
     +                        GKS$K_TEXT_VALIGN_NORMAL)
      no_action = .FALSE.  ! Tell mathmode to output string
      do i=1,2
     	  call mathmode(surf_ws%id, xposition, yposition(i),
     +                 surf_title(i), no_action, err)
      end do
c     Use GKS$TEXT for the third title, since postmini puts the
c     datafile in this line, and there may be backslashes in the
c     filename under Win32. This is really hacky, but...
      call GKS$TEXT(xposition, yposition(3), surf_title(3))
c      
      return
      end
c      
      subroutine surf_plot_axis_gks(axis, xstart, ystart, zstart,
     +                              xend, yend, zend, type)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'gksdefs.cmn'
      include 'surf.cmn'
      include 'io.cmn'
c
      real    xstart,ystart,zstart
      real    xend,yend,zend
      record  /axis_struct/ axis
      integer type
c      
      real   xstart_2d,ystart_2d,xend_2d,yend_2d
      real   tic_x(100),tic_y(100),tic_z(100)
      real   tic_2d_x(100),tic_2d_y(100),tic_2d_z(100)
      real   scaled_tic(100)
      real   xs(100),ys(100),xe(100),ye(100)
      real   longtic_x(100), longtic_y(100), longtic_z(100)
      real   longtic_2d_x(100), longtic_2d_y(100), longtic_2d_z(100)
      real   shorttic_values(100)
      real   scaled_longtic(100),longtic_values(100)
      integer num_of_longtics, num_of_tics
      real   xupvec,yupvec
      real   xmid,ymid
      real   xlong,ylong
      real   xlabel, ylabel
      real   char_height
      real   dist
      integer ndecades
      integer i,ii
      real    major_tic_value, dlab
      integer exp_flag, ilen, string_len
      character*20 value
      logical merr
c
      char_height = 0.015
c      
      call surf_xform_line_3d(xstart,ystart,zstart,
     +                        xend,yend,zend,xstart_2d,ystart_2d,
     +                        xend_2d,yend_2d)
      call surf_draw_line_gks(xstart_2d,ystart_2d,xend_2d,yend_2d,
     +                        axis.color)
c      
c     Set the text attributes to write the label at mid point.
c
      xupvec = - (yend_2d - ystart_2d)
      yupvec =   (xend_2d - xstart_2d)
      if (yupvec .le. 1.0E-4 ) then
          yupvec = -yupvec
          xupvec = -xupvec
      endif
      dist = SQRT(yupvec*yupvec + xupvec * xupvec)
c
c     Calculate the mid point of the axis in 2D and go down more in the
c     direction opposite to the upvec
      xmid = (xstart_2d + xend_2d)/2.0
      ymid = (ystart_2d + yend_2d)/2.0
      if( type .NE. SURF_ZAXIS_TYPE ) then
         dlab = 0.05
      else
c        Z axis label is offset more so there is room for tic mark labels
         if(axis.nonumbering) then
            dlab = 0.05
         else
            dlab = 0.07
         endif
      endif
      xlabel = xmid - (xupvec*dlab)/dist
      ylabel = ymid - (yupvec*dlab)/dist
c
      call gks$set_text_height(char_height)
      call gks$set_text_color_index(axis.color)
      call gks$set_text_fontprec(surf_ws%font,
     +                           surf_ws%precision)
      if( type .NE. SURF_ZAXIS_TYPE ) then
         call gks$set_text_align(GKS$K_TEXT_HALIGN_CENTER,
     +                           GKS$K_TEXT_VALIGN_TOP)
**         call gks$set_text_align(GKS$K_TEXT_HALIGN_CENTER,
**     +                           GKS$K_TEXT_VALIGN_HALF)
         call gks$set_text_upvec(xupvec,yupvec)
         ilen = string_len(axis.label)
**         call gks$text (xlabel,ylabel,axis.label(1:ilen))
         call mathmode(surf_ws%id, xlabel, ylabel, axis.label,
     +                 .FALSE., merr)
      else
         call gks$set_text_align(GKS$K_TEXT_HALIGN_CENTER,
     +                           GKS$K_TEXT_VALIGN_BOTTOM)
c        Z axis gets a vertical label so we can use mathmode
         call gks$set_text_upvec(-1.0 , 0.0)
         call mathmode(surf_ws%id, xlabel, ylabel, axis.label,
     +                 .FALSE., merr)
      endif
c      
      if(axis.tic_spacing .eq. 0.0 ) 
     +    axis.tic_spacing = (axis.max - axis.min) / 10.0
c      
      if(.NOT. axis.nonumbering) then
c        Now place the minimum and maximum.
CC         xmin = xstart_2d - (xupvec*0.03)/dist
CC         ymin = ystart_2d - (yupvec*0.03)/dist
CC         xmax = xend_2d - (xupvec*0.03)/dist
CC         ymax = yend_2d - (yupvec*0.03)/dist
c      
         char_height = 0.75 * char_height
         if( type .eq. SURF_ZAXIS_TYPE ) then
c           Make labels go perpendicular to Z axis for better readability
            call gks$set_text_upvec((xend_2d - xstart_2d),
     +                              (yend_2d - ystart_2d))
         endif
         call gks$set_text_height(char_height)
c      
         major_tic_value = axis.tic_spacing * axis.longtic_frequency
c      
CC         call format_axis_label(axis.min, major_tic_value, value,
CC     +                          exp_flag)
CC         if (xstart_2d .gt. xend_2d ) then
CC            call gks$set_text_align(GKS$K_TEXT_HALIGN_RIGHT,
CC     +                              GKS$K_TEXT_VALIGN_TOP)
CC         else
CC            call gks$set_text_align(GKS$K_TEXT_HALIGN_LEFT,
CC     +                              GKS$K_TEXT_VALIGN_TOP)
CC         endif
CC         if( type .eq. SURF_ZAXIS_TYPE ) then
CC           call gks$set_text_align(GKS$K_TEXT_HALIGN_RIGHT,
CC     +                             GKS$K_TEXT_VALIGN_BOTTOM)
CC         endif
CC         ilen = string_len(value)
CC         call gks$text(xmin,ymin,value(1:ilen))
c      
CC         call format_axis_label(axis.max, major_tic_value, value,
CC     +                          exp_flag)
CC         if(xstart_2d .gt. xend_2d) then
CC           call gks$set_text_align(GKS$K_TEXT_HALIGN_LEFT,
CC     +                             GKS$K_TEXT_VALIGN_TOP)
CC         else
CC           call gks$set_text_align(GKS$K_TEXT_HALIGN_RIGHT,
CC     +                             GKS$K_TEXT_VALIGN_TOP)
CC         endif
CC         if( type .eq. SURF_ZAXIS_TYPE ) then
CC           call gks$set_text_align(GKS$K_TEXT_HALIGN_RIGHT,
CC     +                             GKS$K_TEXT_VALIGN_HALF)
CC         endif
CC         ilen = string_len(value)
CC         call gks$text(xmax, ymax, value(1:ilen))
      endif
c      
      call surf_tic_calculation(axis.log, axis.max, axis.min,
     +        axis.tic_spacing, axis.longtic_frequency,
     +        num_of_tics, num_of_longtics, shorttic_values,
     +        longtic_values, 100) 
c      
      if (num_of_tics .gt. 100) then
          write(luttyo,100)
100       format(' ','*** WARNING : MORE SHORT TICS REQUIRED THAN ',
     +               'ALLOWED ***'/
     +           ' ','*** PROGRAM WILL PROCEED WITH THE LIMIT NUMBER ',
     +               '     ***')
          num_of_tics = 100
      endif
c
      if( type .eq. SURF_XAXIS_TYPE ) then
c
        do i=1,num_of_longtics
           longtic_y(i) = ystart
           longtic_z(i) = zstart
           longtic_x(i) = longtic_values(i)
         end do
         if( num_of_longtics .gt. 0) then
         call surf_scale_array(0.0,1.0,axis.min,axis.max,longtic_x,
     +                         scaled_longtic,num_of_longtics,axis.log)
         call surf_tr_3d_2d (scaled_longtic,longtic_y,longtic_z,
     +             longtic_2d_x,longtic_2d_y,longtic_2d_z,
     +             num_of_longtics)
         endif
      
         do i=1,num_of_tics
            tic_y(i) = ystart
            tic_z(i) = zstart
            if (axis.log) then
               ndecades = INT(FLOAT(i/9))
               tic_x(i) = axis.min * 10.0**ndecades*
     +                    FLOAT(i+1-ndecades*9)
            else
               tic_x(i) = shorttic_values(i)
            endif
         end do
         call surf_scale_array(0.0,1.0,axis.min,axis.max,tic_x,
     +                         scaled_tic,num_of_tics,axis.log)
         call surf_tr_3d_2d (scaled_tic,tic_y,tic_z,tic_2d_x,tic_2d_y,
     +                  tic_2d_z,num_of_tics)
c
      else if(type .eq. SURF_YAXIS_TYPE) then
         do i=1,num_of_longtics
            longtic_x(i) = xstart
            longtic_z(i) = zstart
            longtic_y(i) = longtic_values(i)
         end do
         if( num_of_longtics .gt. 0) then
         call surf_scale_array(0.0,1.0,axis.min,axis.max,longtic_y,
     +                         scaled_longtic,num_of_longtics,axis.log)
         call surf_tr_3d_2d (longtic_x,scaled_longtic,longtic_z,
     +             longtic_2d_x,longtic_2d_y,longtic_2d_z,
     +             num_of_longtics)
         endif
      
         do i=1,num_of_tics
            tic_x(i) = xstart
            tic_z(i) = zstart
            if (axis.log) then
               ndecades = INT(FLOAT(i/9))
               tic_y(i) =  axis.min * 10.0**ndecades*
     +                     FLOAT(i+1-ndecades*9)
            else
               tic_y(i) = shorttic_values(i)
            endif
         end do
         call surf_scale_array(0.0,1.0,axis.min,axis.max,tic_y,
     +                         scaled_tic,num_of_tics,axis.log)
         call surf_tr_3d_2d (tic_x,scaled_tic,tic_z,tic_2d_x,tic_2d_y,
     +          tic_2d_z,num_of_tics)
c
      else if(type .eq. SURF_ZAXIS_TYPE) then
        do i=1,num_of_longtics
           longtic_x(i) = xstart
           longtic_y(i) = ystart
           longtic_z(i) = longtic_values(i)
         end do
         if( num_of_longtics .gt. 0) then
         call surf_scale_array(0.0,1.0,axis.min,axis.max,longtic_z,
     +                        scaled_longtic,num_of_longtics,axis.log)
         call surf_tr_3d_2d (longtic_x,longtic_y,scaled_longtic,
     +           longtic_2d_x,longtic_2d_y,longtic_2d_z,
     +           num_of_longtics)
         endif
      
         do i=1,num_of_tics
            tic_x(i) = xstart
            tic_y(i) = ystart
            if(axis.log) then
               ndecades = INT(FLOAT(i/9))
               tic_z(i) = axis.min*10.0**ndecades*
     +                    FLOAT(i+1-ndecades*9)
            else
               tic_z(i) = shorttic_values(i)
            endif
         end do
         call surf_scale_array(0.0,1.0,axis.min,axis.max,tic_z,
     +                         scaled_tic,num_of_tics,axis.log)
         call surf_tr_3d_2d (tic_x,tic_y,scaled_tic,tic_2d_x,tic_2d_y,
     +           tic_2d_z,num_of_tics)
      endif
c
      if (.NOT. axis.nonumbering) then
         call gks$set_text_align(GKS$K_TEXT_HALIGN_CENTER,
     +                           GKS$K_TEXT_VALIGN_TOP)
         if( type .eq. SURF_ZAXIS_TYPE ) then
           call gks$set_text_align(GKS$K_TEXT_HALIGN_RIGHT,
     +                             GKS$K_TEXT_VALIGN_HALF)
         endif
         do i=1,num_of_longtics
           xlong = longtic_2d_x(i) - (xupvec*0.03)/dist
           ylong = longtic_2d_y(i) - (yupvec*0.03)/dist
           call format_axis_label(longtic_values(i), major_tic_value,
     +                            value, exp_flag)
           ilen = string_len(value)
           call gks$text(xlong, ylong, value(1:ilen))
         end do
      endif
      do i=1,num_of_tics
          xs(i) = tic_2d_x(i)
          ys(i) = tic_2d_y(i)
          xe(i) = tic_2d_x(i) - (xupvec*0.01)/dist
          ye(i) = tic_2d_y(i) - (yupvec*0.01)/dist
      end do
      if(axis.log) then
          do i=1, num_of_tics/9
             ii = 9*i
             xe(ii) = tic_2d_x(ii) - (xupvec*0.02)/dist
             ye(ii) = tic_2d_y(ii) - (yupvec*0.02)/dist
          end do
      else
          ii = 1
          do i=1, num_of_tics
             if( shorttic_values(i) .eq. longtic_values(ii) ) then
                xe(i) = tic_2d_x(i) - (xupvec*0.02)/dist
                ye(i) = tic_2d_y(i) - (yupvec*0.02)/dist
                ii = ii + 1
             endif
          end do
      endif
      do i=1, num_of_tics
          call surf_draw_line_gks(xs(i),ys(i),xe(i),ye(i),axis.color)
      end do
c
c     Reset the gks attributes
c
      xupvec = 0.0
      yupvec = 1.0
      call gks$set_text_upvec (xupvec,yupvec)
      call gks$set_text_align (GKS$K_TEXT_HALIGN_NORMAL,
     +                         GKS$K_TEXT_VALIGN_NORMAL)
      call gks$set_text_color_index (1)
c
      return      
      end      
c      
      subroutine surf_draw_line_gks(sx,sy,ex,ey,color)
      implicit none
c
c     draw_line_gks() draws a line in the 2D plane
c      
      real   sx,sy,ex,ey
      integer    color
c      
      real   x(2),y(2)
c      
      x(1) = sx
      x(2) = ex
      y(1) = sy
      y(2) = ey
      call gks$set_pline_color_index(color)
      call gks$polyline(2,x,y)
c     Reset
      color = 1
      call gks$set_pline_color_index(color)
c      
      return
      end
c      
      subroutine surf_plot_gks(x,y,n)
      implicit none
c      
      integer  n
      real x(n),y(n)
c      
      call gks$polyline (n,x,y)
c
      return
      end
c
       
      integer function surf_xform_gks(xform_number, xform_window,
     1                 xform_viewport)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'surf.cmn'
c
      integer xform_number
      integer gks$set_window, gks$set_viewport, gks$select_xform
      real    xform_window(4)
      real    xform_viewport(4)
c      
      integer status
c      
      status = gks$set_window(xform_number, xform_window(1),
     +         xform_window(2), xform_window(3), xform_window(4))
      if ( status .ne. GKS$_SUCCESS ) then
         surf_xform_gks = -1
         return
      endif
      status = gks$set_viewport (xform_number, xform_viewport(1), 
     +         xform_viewport(2),xform_viewport(3), xform_viewport(4))
      if ( status .ne. GKS$_SUCCESS ) then
         surf_xform_gks = -1
         return
      endif
      status =  gks$select_xform (xform_number)
      if ( status .ne. GKS$_SUCCESS ) then
         surf_xform_gks = -1
         return
      endif
c
      surf_xform_gks = SURF_NORMAL
      return
      end
c      
      subroutine surf_error(message,status,severity)
      implicit none
c
      include 'sysdep.cmn'
      include 'io.cmn'
c      
      character*(*) message
      integer status,severity
c
      write(luttyo,100) message, status
100   format(//
     +       ' ','***************************************************'/
     +       ' ','****  ',a,'  ****'/
     +       ' ','****  Status = ',i5,'  ****'/
     +       ' ','***************************************************')
      if(severity .ne. 0) then
         stop
      endif
c
      return
      end
c      
      subroutine surf_tic_calculation(log, max, min, spacing,
     +                                frequency, i_minor, i_major,
     +                                minor_values, major_values, isize)
      implicit none
c
      logical log
      real    max, min, spacing
      integer frequency, i_minor, i_major, isize
      real    major_values(isize), minor_values(isize)
c      
      integer  i
      real ll, delta, cur, tic_major
c      
      if (.NOT. log) then
c        Find first major tic mark below min
         tic_major = spacing*FLOAT(frequency)
         if( min .lt. 0.0 ) then
            cur = FLOAT(INT(min/tic_major) - 1)*tic_major
         else
            cur = FLOAT(INT(min/tic_major))*tic_major
         endif
c        Normalize to 0.0 to 1.0 scale
         cur = (cur - min)/(max - min)
         delta = tic_major/
     1           (FLOAT(frequency)*(max - min))
         i = 0
         i_minor = 0
         i_major = 0
c        We use 1.0001 instead of 1.0 as a safeguard agains fp rounding
         do while(cur .le. 1.0001)
c           Plot tic marks only if within axis bounds
c           We use -0.0001 instead of 0.0 as a safeguard agains fp rounding
            if(cur .ge. -0.0001) then
c              Determine where a big tic mark should go
               if( MOD(i,frequency) .eq. 0 ) then
                  i_major = i_major + 1
                  if( i_major .le. isize )
     +               major_values(i_major) = cur*(max-min) + min
               endif
               i_minor = i_minor + 1
               if( i_minor .le. isize )
     +            minor_values(i_minor) = cur*(max-min) + min
            endif
            cur = cur + delta
            i = i + 1
         enddo
      else
         i_major = INT(ALOG10 (max) - ALOG10(min) - 1.0)
         ll = min * 10.0**(i_major + 1)
         if (ll .LT. 0.99999*max ) i_major = i_major + 1
         i_minor = 9 * ALOG10(max/min) - 1
         ll = FLOAT(i_minor)/10.
         ll = min*10.0**ll
         if (ll .GT. max) i_minor = i_minor - 1
         if (i_major .GT. 10)  then
c           Multiply by 1.0000001 to avoid rounding
            i_major = i_major / 10
            if (i_major .GT. isize) then
               i_major = isize
            endif
            do i=1,i_major
               major_values(i) = min * 1.0000001*(10.0**(10.0*i))
            end do
         else
            if (i_major .GT. isize) then
               i_major = isize
            endif
            do i=1,i_major
                 major_values(i) = min * 1.0000001*(10.0**i)
            end do
         endif
      endif
c
      end
c      
      subroutine surf_clear_ws_gks
      implicit none
c
c     Clear ws gks --- clears the workstation
c
      include 'const.h'
      include 'gksdefs.cmn'
      include 'surf.cmn'
c
      call gks$clear_ws(surf_ws%id, GKS$K_CLEAR_ALWAYS)
c
      return
      end
