subroutine hide (x,y,xg,g,xh,h,ng,maxdim,n1,nfns,title,xlnth,ylnth
     *   ,xmin,deltax,ymin,deltay,iform)
c
c
c @(#)hide.f	1.1  4/3/89 
c This subroutine produces a 2-dimensional representation of
c a 3-dimensional figure or surface.
c The first call to hide is for initialization and plotting
c the curve farthest in the foreground.  on each subsequent
c call, a curve farther in the backround is plotted.
c x is the abcissa array for the curve to be plotted by
c hide on this call.  the x values must be increasing.
c if x(i) ge x(i+1) for some i, maxdim will be set to zero,
c and a return will be executed.
c y is the ordinate array.
c g vs xg is the current visual maximum function on each
c return from hide.
c xh and h are working arrays.
c on each return from hide, ng is the number of points in
c the current maximum function.
c on the first call, ng is a nonpositive integer which
c specifies certain options.
c -1 do not draw the 8 1/2 by 11 inch border.
c -2 plot unhidden minimum rather than maximum.  in this
c    case, g vs xg will be the negative of the visual
c    minimum function.
c -3 do not plot border, plot minimum rather than maximum.
c  0 plot border, plot maximum.
c if the border is drawn, its left, bottom corner will be
c where the plotting reference point was just before the
c first call to hide, and the reference point will be moved
c 1 inch right and 2 inches up.  if the border is not drawn,
c the reference point will not be moved by hide.
c maxdim is the dimension in the calling program of the
c arrays xg, g, xh, and h.  if one of these arrays would
c have been overflowed, maxdim is set equal to its negative,
c and a return is executed.
c n1 is the number of points (x(i),y(i)) to be plotted in
c a given call to hide.
c if n1 is less than 0, y vs x will not be plotted, but on
c subsequent calls, plotting will be done as if
c ((x(i),y(i)),i=1,-n1) had been plotted (where unhidden).
c n1 will be returned as its absolute value.
c nfns is the total no. of curves to be plotted for this
c graph if translating the arrays to simulate stepping in
c the depth dimension is desired.  if no translation is
c desired, nfns should be negative.  if the same translation
c as in the previous call to hide is desired, nfns should be
c zero.  the nfns = 0 option may be specified for individual
c curves after the first for a given graph.  all
c translations which are performed will have equal step size
c determined by the values in the initial call for xlnth,
c ylnth, and nfns.
c title is an 80-character title.
c if title(1) = 6hnottle, the title will not be plotted.
c title(1) and xh(1) or h(1) may be the same location if the
c title is not needed after it is plotted.
c xlnth is the length in inches of the horizontal axis.
c if xlnth is less than 0, the x-axis and the depth axis
c will not be drawn.  in any case, unless this option is
c suppressed through nfns, the ith curve will be translated
c (i - 1) * (9.0 - abs(xlnth)) / (nfns - 1) inches to the left.
c this plus a similar vertical translation is done to simulate
c stepping in the depth dimension.
c xmin - (9.0 - abs(xlnth)) * deltax will be the abcissa
c value at the plotting reference point (which is where the
c horizontal and vertical axes would intersect if drawn*.
c ylnth is the length of the vertical axis in inches.
c if ylnth is less than 0, the vertical and depth axes will
c not be drawn.  but in any case, unless this option is
c suppressed through nfns, the ith curve will be translated
c (i - 1) * (6.0 - abs(ylnth)) / (nfns - 1) inches up
c to simulate stepping in the depth dimension.  ymin -
c (6.0 - abs(ylnth)) * deltay will be the ordinate value
c at the plotting reference point.
c if translations are performed, x and y will be restored to
c their original values before the return to the calling
c program.
c note that if abs(xlnth) = 9, and abs(ylnth) = 6, there will be
c no translation, and, if border and axes are not drawn, the
c dimensions of the plot are unspecified.
c if the axes and border are drawn, the top of the vertical
c axis and the right end of the horizontal axis are fixed
c relative to the border, and the depth axis joins the left
c end of the horizontal axis and the bottom of the vertical axis.
c xmin is a lower bound for x.
c deltax is the x data increment per inch for the plot.
c xmin and deltax determine the plotting scale for x.
c (see above.)
c ymin and deltay, similarly, determine the scale for y.
c if an error return is made from hide, all further calls
c will result only in the execution of a return unless
c maxdim is reset to a positive value.
c
c conversion to gould: 6/27/85 (jfb)
c
c
      dimension x(1),y(1),xg(1),g(1),h(1),xh(1)
      integer iform
      character*(*) title
      character*6 nottle
c
c
c the only purpose of the following equivalence statement
c is to save storage.
c
c
      equivalence (k1,iwhich) , (k2,slope) , (fnsm1,z1) , (iggp1,k1) , (
     *   k1,n2)
c
c
c eps1 is the relative abcissa increment used to simulate
c discontinuities in the maximum function.
c
c
      data eps1 / 1.0e-6 /
      data nottle / 'nottle'/
c
c
c the following statement function computes the ordinate on
c the line joining (xi,yi) and (xip1,yip1) corresponding to
c the abcissa xx.
c
c
      f(xx,xi,yi,xip1,yip1) = yi+(xx-xi)*(yip1-yi)/(xip1-xi)
c
      if (maxdim.le.0) return
      ifplot = 1
      if (n1.le.0) then
         n1 = -n1
         ifplot = 0
      endif
      do 10 i = 2, n1
         if (x(i-1).lt.x(i)) go to 10
         maxdim = 0
         go to 20
   10 continue
      if (ng.gt.0) go to 60
      if (n1+4.0.le.maxdim) go to 30
      maxdim = -maxdim
   20 return
c
c we want sign = 1 if we are plotting maximum, = -1 if
c minimum
c
   30 sign = 1.0e+00
      if (ng.lt.-1) sign = -1.0e+00
c
c the kth curve to be plotted will (optionally) be
c translated by the vector (-dxin,dyin) * (k - 1) to
c simulate stepping in the depth dimension.
c
      fnsm1 = 0.0e+00
      if (nfns.gt.0) then
         fnsm1 = nfns-1
         dxin = (8.0e+00-abs(xlnth))*deltax/fnsm1
         dyin = (6.0e+00-abs(ylnth))*deltay/fnsm1
      endif
c
c systems routine plot moves the pen to a point whose
c coordinates are specified in inches by the first two
c parameters.  the pen is picked up if the absolute value of
c the third parameter is 3, is put down if 2, and is left as
c after last call if 1.  if the third parameter is negative,
c a new reference point will be established.
c
      if (ng.ne.-1.and.ng.ne.-3) then
c
c draw 8 1/2 by 11 inch border.
c
         call plot (0.0e+00,0.0e+00,3)
         call plot (8.00e+00,0.0e+00,2)
         call plot (8.00e+00,6.0e+00,1)
         call plot (0.0e+00,6.0e+00,1)
         call plot (0.0e+00,0.0e+00,1)
         call plot (1.0e+00,2.0e+00,-3)
      endif
c
c call systems routine to plot the 80-character title.
c the first two arguments are the coordinates in inches
c relative to the reference point of the lower left-hand
c corner of the first character.  the third argument
c determines the height in inches of the characters.  the
c fifth argument gives the angle relative to horizontal of
c the plotted characters.
c
      ltit = len(title)
      if (title.ne.nottle) call symbol (-0.28e+00,-1.0e+00,0.14e+00,    
     *   title,0.0e+00,ltit)
      if (xlnth.ge.0.0e+00) then
c
c call systems routine to draw the horizontal axis.  the
c left end is specified in inches relative to the reference
c point by the first two arguments.
c
         call axis (8.0e+00-xlnth,0.0e+00,' ',-1,xlnth,0.e+00,xmin,     
     *      deltax,0.e+00,iform)
         if (ylnth.lt.0.0e+00) go to 40
c
c draw the depth axis.
c
         call plot (8.0e+00-xlnth,0.0e+00,3)
         call plot (0.0e+00,6.0e+00-ylnth,2)
      endif
      if (ylnth.ge.0.0e+00) call axis (0.0e+00,6.0e+00-ylnth,' ',1,ylnth
     *   ,90.0e+00,ymin,deltay,-1,iform)
c
c curves successively farther in the background will be
c plotted where they are not hidden by g vs xg.  g vs xg
c will be updated each time a new curve is drawn and will be
c the visual maximum (or minimum) function of the curves
c already plotted.
c
   40 indext = 3
      do 50 j = 1, n1
         xg(indext) = x(j)
         g(indext) = sign*y(j)
         indext = indext+1
   50 continue
c
c the following precautionary step is used in place of a
c test in subroutine lookup to see if the value for which we
c want an index is outside the table.
c the last xg value will be set equal to the last abcissa
c of the curve to be plotted in the next call to hide.
c
      eps = eps1*(abs(xmin)+abs(deltax))
      ng = n1+4
      xg(1) = -fnsm1*dxin+xmin-abs(xmin)-abs(xg(3))-1.0e+00
      xg(2) = xg(3)-eps
      xg(n1+3) = xg(n1+2)+eps
      zz = ymin
      if (sign.lt.0.0e+00) zz = -ymin-50.0e+00*deltay
      g(1) = zz
      g(2) = zz
      g(n1+3) = zz
      g(ng) = zz
c
c call systems routine to produce a line plot of
c (x(i), y(i), i=1,n1) - this is the curve farthest in the
c foreground.
c xstart is the x value at the reference point.
c
      xstart = xmin-(9.0e+00-abs(xlnth))*deltax
c
      if (ifplot.eq.1) call pdatax (x,y,n1,xstart,deltax,ymin,deltay)
      dxkk = 0.0e+00
      dykk = 0.0e+00
      relinc = deltax/deltay
      xg(ng) = sign
      return
c
c statement 5000 is reached if any except the curve farthest
c in the foreground is to be plotted.
c
   60 sign = xg(ng)
      xg(ng) = x(n1)
c
c translate the arrays before plotting to simulate stepping
c in the depth dimension.
c
      if (nfns) 100, 80, 70
   70 dxkk = dxkk+dxin
      dykk = dykk+dyin
   80 do 90 j = 1, n1
         y(j) = sign*(y(j)+dykk)
         x(j) = x(j)-dxkk
   90 continue
  100 call lookup (x(1),xg(1),jj)
      if (jj.ge.maxdim) go to 220
      do 110 j = 1, jj
         xh(j) = xg(j)
         h(j) = g(j)
  110 continue
      ig = jj+1
      xh(ig) = x(1)
      h(ig) = f(x(1),xg(jj),g(jj),xg(ig),g(ig))
c
c we will be making table lookups for an increasing sequence
c of numbers - therefore, we do not have to search from the
c first of the (xg and x) tables each time.  hence indexg
c and indext.
c
      indexg = jj
      indext = 1
      z1 = x(1)
      f1 = h(ig)-y(1)
      it = 2
      jj = ig
      if (h(ig).lt.y(1)) then
         if (jj.ge.maxdim) go to 220
         jj = ig+1
         h(jj) = y(1)
         xh(jj) = z1+eps
      endif
      last = 0
      x1 = z1
c
c find the first zero, z2, of the function g-y to the right
c of z1.
c
  120 if (xg(ig).ge.x(it)) then
c
c do not jump if we are to look for a zero between x1 and
c x(i).
c
         iwhich = 0
         x2 = x(it)
         f2 = f(x2,xg(ig-1),g(ig-1),xg(ig),g(ig))-y(it)
         it = it+1
      else
c
c come to 1001 if we are to look for a zero between x1 and
c xg(ig).
c
         x2 = xg(ig)
         iwhich = 1
         f2 = g(ig)-f(x2,x(it-1),y(it-1),x(it),y(it))
         ig = ig+1
      endif
c
c the function (g - y) has a zero z2 such that x1 le z2 le x2
c if and only if (g - y at x1) * (g - y at x2) le 0.
c (g - y is assumed, for plotting purposes, to be linear on
c each interval (x1, x2).)
c
      if (f1*f2.le.0.e+00) then
         if (f1.ne.f2) then
            slope = (f2-f1)/(x2-x1)
            igg = ig-1-iwhich
            itt = it-2+iwhich
            if (abs(slope*relinc).gt.eps1) go to 130
c
c if g and y differ imperceptibly (for plotting purposes)
c on the interval (x1, x2), set z2 = x2.  this step prevents
c division by zero.
c
            z2 = x2
            go to 150
c
c otherwise, compute the zero z2.
c
  130       z2 = x1-f1/slope
            go to 150
         endif
      endif
c
c if no zero was found between x1 and x2, continue the
c search for zeroes.
c
      x1 = x2
      f1 = f2
      if (it.le.n1) go to 120
c
c if the end of the x table has been reached, consider the
c interval from the last zero found to the end of the x
c table (plot, update maximum function as indicated).
c
  140 last = 1
      z2 = x(n1)
      call lookup (z2,xg(indexg),igg)
      igg = indexg+igg-1
      itt = n1-1
c
c it is necessary to plot y vs x on the interval (z1, z2)
c only if y is unhidden at each zz such that z1 lt zz lt z2.
c we choose zz near the left end of the interval for
c efficiency in the table lookup.
c note that it is more efficient to choose this value for zz
c than, say, 0.99 * x(indext) + 0.01 * x(indext + 1), which
c would eliminate one of the two table lookups, but would
c necessitate a test to determine if zz was between z1 and z2.
c
  150 zz = 0.99e+00*z1+0.01e+00*z2
      call lookup (zz,x(indext),k1)
      call lookup (zz,xg(indexg),k2)
      k1 = k1+indext-1
      k2 = k2+indexg-1
      if (f(zz,x(k1),y(k1),x(k1+1),y(k1+1)).le.f(zz,xg(k2),g(k2),xg(k2+1
     *   ),g(k2+1))) then
c
c
c if y is hidden between z1 and z2, update the maximum
c function.
c for generality, the maximum function is updated even if
c this is the (nfns)th curve.
c
c
         if (jj+igg-indexg.ge.maxdim) go to 220
         if (indexg.ne.igg) then
            j1 = indexg+1
            do 160 i = j1, igg
               jj = jj+1
               xh(jj) = xg(i)
               h(jj) = g(i)
  160       continue
         endif
         jj = jj+1
         xh(jj) = z2
         h(jj) = f(z2,xg(igg),g(igg),xg(igg+1),g(igg+1))
         indexg = igg
         indext = itt
      else
c
c if y is not hidden between z1 and z2, update the maximum
c function and plot.
c
         ngraph = itt-indext+2
         if (jj+ngraph-1.gt.maxdim) go to 220
         n2 = jj
         if (ngraph.ne.2) then
            j1 = indext+1
            do 170 i = j1, itt
               jj = jj+1
               xh(jj) = x(i)
               h(jj) = y(i)
  170       continue
         endif
         jj = jj+1
         xh(jj) = z2
         h(jj) = f(z2,x(itt),y(itt),x(itt+1),y(itt+1))
c
c call systems routine to produce line plot of
c (xh(i), h(i), i=n2, n2 + ngraph - 1).
c
         if (ifplot.eq.1) call pdatax (xh(n2),h(n2),ngraph,xstart,deltax
     *      ,sign*ymin,sign*deltay)
c
c
         indext = itt
         indexg = igg
      endif
      if (last.ne.1) then
         x1 = x2
         f1 = f2
         z1 = z2
c
c after plotting and/or updating the maximum function on the
c interval (z1, z2), search for the next zero if the end of
c the abcissa table xt has not been reached.
c
         if (it.le.n1) go to 120
         go to 140
      endif
c
c after y vs x has been plotted, finish updating and store
c the new maximum function.
c allow for the possibility that the previous maximum
c function extends to the right of the function just
c plotted.
c
      if (xg(ng).le.xg(ng-1)) ng = ng-1
      if (xg(ng).gt.x(n1)) then
         if (jj+3+ng-igg.gt.maxdim) go to 220
         xh(jj+1) = xh(jj)+eps
         jj = jj+1
         h(jj) = f(x(n1),xg(igg),g(igg),xg(igg+1),g(igg+1))
         iggp1 = igg+1
         do 180 j = iggp1, ng
            jj = jj+1
            xh(jj) = xg(j)
            h(jj) = g(j)
  180    continue
      endif
      ng = jj+2
      if (ng.gt.maxdim) go to 220
      do 190 i = 1, jj
         g(i) = h(i)
         xg(i) = xh(i)
  190 continue
c
c the following precautionary step is used in place of a
c test in subroutine lookup to see if the value for which we
c want an index is outside the table.
c the last xg value will be set equal to the last abcissa
c of the next curve to be plotted.
c
      xg(jj+1) = xg(jj)+eps
      g(jj+1) = ymin+dykk
      if (sign.lt.0.e+00) g(jj+1) = -ymin-50.0e+00*deltay+dykk
      g(ng) = g(jj+1)
c
c restore arrays x and y before returning.
c
  200 if (nfns.ge.0) then
         do 210 i = 1, n1
            x(i) = x(i)+dxkk
            y(i) = sign*y(i)-dykk
  210    continue
      endif
      xg(ng) = sign
      return
c
c if statement 700 is reached, dimensions would have been
c exceeded.  see comments on calling sequence for hide.
c
  220 maxdim = -maxdim
      go to 200
      end
c
c
      subroutine lookup (x,xtbl,j)
c
c this subroutine is called by hide to perform a table
c lookup.  because of precautions taken in hide, a test to
c see if x is outside the table is unnecessary.
c
      dimension xtbl(1)
      j = 2
   10 if (xtbl(j)-x) 20, 30, 40
   20 j = j+1
      go to 10
   30 return
   40 j = j-1
      return
      end
c
c
      subroutine pdatax (x,y,n,xm,dx,ym,dy)
c
c     calcomp/dipl compatable data ploting routine
c
      dimension x(n),y(n)
      data cx,cy / 2*0.0e+00 /
      px(i) = (x(i)-xm)/dx
      py(i) = (y(i)-ym)/dy
c
c
      i1 = 1
      i2 = 1
      test1 = amax1(abs(cx-px(1)),abs(cy-py(1)))
      test2 = amax1(abs(cx-px(n)),abs(cy-py(n)))
      if (test1.ge.test2) then
c
c     if (amax1(abs(cx-px(1)),abs(cy-py(1))) .lt. amax1(abs(cx-px(n)),
c    x     abs(cy-py(n)))) go to 10
c
         i1 = n
         i2 = -i2
      endif
      call plot (px(i1),py(i1),3)
      do 10 i3 = 2, n
         i1 = i1+i2
         call plot (px(i1),py(i1),2)
   10 continue
      cx = px(i1)
      cy = py(i1)
      return
      end