* Date: Tue, 24 Nov 92 10:20:47 -0700
* From: seeger@gem.LANL.GOV

      SUBROUTINE PGPLOT10
C
C     Subroutines to implement high-level functions from the author's 
C       PLOT-10 library, and to interpret calls to PLOT-10 primatives 
C       as corresponding entries in the PGPLOT library.  When logarithmic
C       axes have been defined by calling GRAPH, the PLOT-10 move and draw
C       calls will automatically take logarithms before plotting.  Additional
C       high-level calls (CONTOUR, CURVE, ERRBAR) which duplicate existing 
C       PGPLOT functions (but which support logarithmic scaling) are included 
C       at the end of this file.
C
C     P. A. Seeger, Los Alamos National Laboratory, Oct. 6, 1992
C
C     Entries:
C     ANCHO     DASHA     DASHR     DRAWA     DRAWR     DRWABS    DRWREL    
C     DSHABS    DSHREL    DWINDO    ERASE     FINITT    GRAPH     HDCOPY    
C     HLABEL    INITT     MOVABS    MOVEA     MOVER     MOVREL    OSTRING   
C     PNTABS    PNTREL    POINTA    POINTR    RESET     SCURSR    SEEDW     
C     SEELOC    SEETW     SETCOLOR  SWINDO    SYMBOL    TERM      TWINDO    
C     VCURSR    VLABEL    VWINDO    
C
C     Auxiliary Externals:
C     CHECK_TEXT   LINAXIS   LOGAXIS   
C
C     PGPLOT Externals:
C     PGBEGIN   PGBOX     PGCURSE   PGDRAW    PGEND     PGETXT    PGMOVE    
C     PGMTEXT   PGNUMB    PGPOINT   PGQINF    PGQPOS    PGQVP     PGSCH     
C     PGSCI     PGSLS     PGSLW     PGVPORT   PGWINDOW  
C
      IMPLICIT NONE
C
C     Variables which may occur in calling sequences of entry points
      CHARACTER TITLE*(*),SUBTITLE*(*),XLABEL*(*),YLABEL*(*),CH*1
      REAL XMN,XMX,X,YMN,YMX,Y
      INTEGER LOGX,LOGY,L,IX,IY,LX,LY
C
C     Local temporary variables:
      CHARACTER XOPTION*7,YOPTION*8,XFORMAT*6,YFORMAT*6,STRING*12,
     1 NEWSTRING*255
      REAL XX,XMAJOR,XMINOR,YY,YMAJOR,YMINOR
      INTEGER I,J,IXMAJOR,IXTICK,IYMAJOR,IYTICK,NC,LL,IX2,IY2
      LOGICAL LPOINT
C
C     Local variables to be saved between entries:
      CHARACTER TERMTYPE*40,HDCPYTYPE*40
      REAL ASPECT,X1,X2,XMIN,XMAX,XPERPIX,X0,Y1,Y2,YMIN,YMAX,YPERPIX,Y0
      INTEGER LINE,ISYMB,ITERM,IHDCPY
      LOGICAL XLOG,YLOG,PLOT10
      SAVE PLOT10,TERMTYPE,HDCPYTYPE,ITERM,IHDCPY,ASPECT,LINE,ISYMB,
     1     X1,X2,XMIN,XMAX,XLOG,XPERPIX,X0,
     2     Y1,Y2,YMIN,YMAX,YLOG,YPERPIX,Y0
      DATA PLOT10, TERMTYPE,HDCPYTYPE,ITERM,IHDCPY,ASPECT,LINE,ISYMB
     1    /.FALSE.,'?',     '?',      1,    1,     0.75,  -1,  -1/
      DATA X1,X2,   XPERPIX,X0,Y1,Y2,  YPERPIX,Y0
     1    /0.,1023.,1.,     0.,0.,767.,1.,     0./
C
      ENTRY INITT(X)
C     Initialize graphics terminal device
      plot10 = .true.
      call pgbegin(0, termtype(1:iterm), 1, 1)
      if (termtype.eq.'?') then
         call pgqinf('dev/type',termtype,iterm)
         iterm = min0(iterm,index(termtype,'/')+3)
      end if
      return
C
      ENTRY HDCOPY
C     Initialize hardcopy output device
      call pgbegin(0, hdcpytype(1:ihdcpy), 1, 1)
      if (termtype.eq.'?') then
         call pgqinf('dev/type',hdcpytype,ihdcpy)
         ihdcpy = min0(ihdcpy,index(hdcpytype,'/')+3)
      end if
      return
C
      ENTRY ERASE
      ENTRY RESET
C     Erase the screen
      call pgpage
      ENTRY TERM(IX,IY)
      return
C
      ENTRY FINITT(X,Y)
C     Exit graphics
      call pgend
      return
C
C     Define graph area in world co-ordinates, to fill the available
C       viewport area.  Place two lines of titles at top, and axis labels
C       on bottom and left.  Draw box with tick marks and labels; logarithmic
C       if the corresponding LOGX or LOGY value is non-zero.  XMN, XMX, YMN,
C       and YMX will be changed to rounded values.
C
      ENTRY GRAPH(TITLE, SUBTITLE, XLABEL, YLABEL,
     1            XMN, XMX, LOGX,   YMN, YMX, LOGY)
      call pgetxt
      call pgqvp(3,x1,x2,y1,y2)
      aspect = (y2-y1+1.)/(x2-x1+1.)
      call pgvport(0.20*aspect, 1.0-0.05*aspect, 0.10, 0.85 )
C
      call pgsch(2.0)
      call pgslw(3)
      call check_text(title,newstring,j,plot10)
      call pgmtext('T', 2.0, 0.5, 0.5, newstring(1:j))
C
      call pgsch(1.5)
      call pgslw(2)
      call check_text(subtitle,newstring,j,plot10)
      call pgmtext('T', 1.2, 0.5, 0.5, newstring(1:j))
C
      call check_text(xlabel,newstring,j,plot10)
      call pgmtext('B', 2.16, 0.5, 0.5, newstring(1:j))
C
      call check_text(ylabel,newstring,j,plot10)
      call pgmtext('L', 4.0, 0.5, 0.5, newstring(1:j))
      call pgsch(1.0)
      call pgslw(1)
C
      xlog = logx.gt.0
      ylog = logy.gt.0
      if (xlog) then
         call logaxis(xmn, xmx, ixmajor, ixtick, xminor, xformat)
         xmin = alog10(xmn)
         xmax = alog10(xmx)
         xmajor = 1.
         ixtick = 9
         xoption = 'BCLNTS'
         if (.not.ylog) xoption = 'ABCLNTS'
      else
         call linaxis(xmn, xmx, ixmajor, ixtick, xformat)
         xmin = xmn
         xmax = xmx
         xmajor = (xmax-xmin)/float(ixmajor)
         xoption = 'BCNTS'
         if (.not.ylog) xoption = 'ABCNTS'
      end if
C
      if (ylog) then
         call logaxis(ymn, ymx, iymajor, iytick, yminor, yformat)
         ymin = alog10(ymn)
         ymax = alog10(ymx)
         ymajor = 1.
         iytick = 9
         yoption = 'BCLNTSV'
         if (.not.xlog) yoption = 'ABCLNTSV'
      else
         call linaxis(ymn, ymx, iymajor, iytick, yformat)
         ymin = ymn
         ymax = ymx
         ymajor = (ymax-ymin)/float(iymajor)
         yoption = 'BCNTSV'
         if (.not.xlog) yoption = 'ABCNTSV'
      end if
C
      call pgwindow(xmin, xmax, ymin, ymax)
      call pgbox(xoption, xmajor, ixtick, yoption, ymajor, iytick)
      xperpix = (xmax-xmin)/(x2-x1)/(1.-0.25*aspect)
      x0 = x1 + 0.20*aspect*(x2-x1)
      yperpix = (ymax-ymin)/(y2-y1)/0.75
      y0 = y1 + 0.10*(y2-y1)
C
      if (xlog) then
C        Make sure ends of X-axis are labeled
         if (amod(abs(xmin)+0.01,1.) .gt. 0.02) then
            i = nint(xmin-0.5)
            j = nint(xmn/10.**i)
            call pgnumb(j, i, 0, string, nc)
            call pgmtext('B', 1.25, 0., 0.5, string(1:nc))
         end if
         if (amod(abs(xmax)+0.01,1.) .gt. 0.02) then
            i = nint(xmax-0.5)
            j = nint(xmx/10.**i)
            call pgnumb(j, i, 0, string, nc)
            call pgmtext('B', 1.25, 1., 0.5, string(1:nc))
         end if
      end if
C
      if (ylog) then
C        Make sure ends of Y-axis are labeled
         if (amod(abs(ymin)+0.01,1.) .gt. 0.02) then
            i = nint(ymin-0.5)
            j = nint(ymn/10.**i)
            call pgnumb(j, i, 0, string, nc)
            call pgmtext('LV', aspect, 0., 1., string(1:nc))
         end if
         if (amod(abs(ymax)+0.01,1.) .gt. 0.02) then
            i = nint(ymax-0.5)
            j = nint(ymx/10.**i)
            call pgnumb(j, i, 0, string, nc)
            call pgmtext('LV', aspect, 1., 1., string(1:nc))
         end if
      end if
      return
C
C     Following entries either move to a new point without drawing, or
C       move to a new point and draw a single dot or a symbol, accounting 
C       for possibility of log scales.
C
      ENTRY MOVEA(X,Y)
      lpoint = .false.
      xx = x
      yy = y
      go to 100
C
      ENTRY MOVABS(IX,IY)
      lpoint = .false.
      xx = xmin + xperpix*(float(ix)-x0)
      yy = ymin + yperpix*(float(iy)-y0)
      go to 110
C
      ENTRY POINTA(X,Y)
      lpoint = .true.
      ll = -1
      xx = x
      yy = y
      go to 100
C
      ENTRY SYMBOL(X,Y,L)
      lpoint = .true.
      ll = l
      xx = x
      yy = y
      go to 100
C
      ENTRY ANCHO(L)
      lpoint = .true.
      ll = l
      call pgqpos(xx,yy)
      go to 110
C
      ENTRY PNTABS(IX,IY)
      lpoint = .true.
      ll = -1
      xx = xmin + xperpix*(float(ix)-x0)
      yy = ymin + yperpix*(float(iy)-y0)
      go to 110
C
      ENTRY MOVREL(IX,IY)
      lpoint = .false.
      go to 80
C
      ENTRY PNTREL(IX,IY)
      lpoint = .true.
      ll = -1
80    call pgqpos(xx,yy)
      xx = xx + xperpix*float(ix)
      yy = yy + yperpix*float(iy)
      go to 110
C
      ENTRY MOVER(X,Y)
      lpoint = .false.
      go to 90
C
      ENTRY POINTR(X,Y)
      lpoint = .true.
      ll = -1
90    continue
      call pgqpos(xx,yy)
      if (xlog) then
         xx = (10.**xx) + x
      else
         xx = xx + x
      end if
      if (ylog) then
         yy = (10.**yy) + y
      else
         yy = yy + y
      end if
100   continue
      if (xlog) then
         if (xx.le.0.) then
            xx = -38.
         else
            xx = alog10(xx)
         end if
      end if
      if (ylog) then
         if (yy.le.0.) then
            yy = -38.
         else
            yy = alog10(yy)
         end if
      end if
110   if (lpoint) then
         call pgpoint(1, xx, yy, ll)
      else
         call pgmove(xx,yy)
      end if
      return
C
C     Following entries draw a solid or dashed line from the current
C       location to a new point, accounting for possibility of log scales.
C
      ENTRY DRAWA(X,Y)
      ll = 0
      xx = x
      yy = y
      go to 200
C
      ENTRY DRWABS(IX,IY)
      ll = 0
      xx = xmin + xperpix*(float(ix)-x0)
      yy = ymin + yperpix*(float(iy)-y0)
      go to 210
C
      ENTRY DASHA(X,Y,L)
      ll = l
      xx = x
      yy = y
      go to 200
C
      ENTRY DSHABS(IX,IY,L)
      ll = l
      xx = xmin + xperpix*(float(ix)-x0)
      yy = ymin + yperpix*(float(iy)-y0)
      go to 210
C
      ENTRY DRWREL(IX,IY)
      ll = 0
      go to 180
C
      ENTRY DSHREL(IX,IY,L)
      ll = l
180   call pgqpos(xx,yy)
      xx = xx + xperpix*float(ix)
      yy = yy + yperpix*float(iy)
      go to 210
C
      ENTRY DRAWR(X,Y)
      ll = 0
      go to 190
C
      ENTRY DASHR(X,Y,L)
      ll = l
190   continue
      call pgqpos(xx,yy)
      if (xlog) then
         xx = (10.**xx) + x
      else
         xx = xx + x
      end if
      if (ylog) then
         yy = (10.**yy) + y
      else
         yy = yy + y
      end if
200   continue
      if (xlog) then
         if (xx.le.0.) then
            xx = -38.
         else
            xx = alog10(xx)
         end if
      end if
      if (ylog) then
         if (yy.le.0.) then
            yy = -38.
         else
            yy = alog10(yy)
         end if
      end if
210   if (ll.lt.0) then
         call pgmove(xx,yy)
      else
         if (ll.ne.line) then
            line = mod(ll,5)
            call pgsls(line+1)
         end if
         call pgdraw(xx,yy)
      end if
      return
C
C     Set boundaries of viewport or of world co-ordinate system
C
      ENTRY VWINDO(XMN, X, YMN, Y)
      xmax = xmn+x
      ymax = ymn+y
      go to 250
C
      ENTRY DWINDO(XMN, XMX, YMN, YMX)
      xmax = xmx
      ymax = ymx
250   xmin = xmn
      ymin = ymn
      call pgwindow(xmin, xmax, ymin, ymax)
      xperpix = (xmax-xmin)/(x2-x1+1.)
      x0 = x1
      yperpix = (ymax-ymin)/(y2-y1+1.)
      y0 = y1
      return
C
      ENTRY SWINDOW(IX, LX, IY, LY)
      ix2 = x+lx
      iy2 = y+ly
      go to 260
C
      ENTRY TWINDO(IX, LX, IY, LY)
      ix2 = lx
      iy2 = ly
260   call pgqvp(3, x1, x2, y1, y2)
      call pgvport(float(ix)/(x2-x1), float(ix2)/(x2-x1),
     1             float(iy)/(y2-y1), float(iy2)/(y2-y1))
      call pgqvp(3, x1, x2, y1, y2)
      xperpix = (xmax-xmin)/(x2-x1+1.)
      x0 = x1
      yperpix = (ymax-ymin)/(y2-y1+1.)
      y0 = y1
      return
C
C     Write horizontal or vertical text strings 
C
      ENTRY HLABEL(TITLE, IX, IY, LX, LY)
      string = 'B'
      xx = xperpix*(float(ix)-x0)/(xmax-xmin)
      ly = nint((y2-y1)/40.)
      lx = (3*ly)/4
      yy = -(float(iy)-y0)/((y2-y1)/40.)
      go to 300
C
      ENTRY OSTRING(TITLE,L)
      call pgqpos(xx,yy)
      if (l.eq.0) then
         string = 'B'
         xx = (xx-xmin)/(xmax-xmin)
         yy = (yy-ymin)/yperpix/(y2-y1)*40.
      else
         string = 'L'
         xx = (yy-ymin)/(ymax-ymin)
         yy = (xx-xmin)/yperpix/(y2-y1)*40.
      end if
      go to 300
C
      ENTRY VLABEL(TITLE, IX, IY, LX, LY)
      string = 'L'
      xx = yperpix*(float(iy)-y0)/(ymax-ymin)
      ly = nint((y2-y1)/40.)
      lx = (3*ly)/4
      yy = -(float(iy)-y0)/((y2-y1)/40.)
300   continue
      call check_text(title,newstring,j,plot10)
      call pgmtext(string, yy, xx, 0.5, newstring(1:j))
      return
C
C     Find ("see") various window parameters
C
      ENTRY SEETW(IX, LX, IY, LY)
      ix = x1
      lx = x2
      iy = y1
      ly = y2
      return
C
      ENTRY SEEDW(XMN,XMX,YMN,YMX)
      if (xlog) then
         xmn = 10.**xmin
         xmx = 10.**xmax
      else
         xmn = xmin
         xmx = xmax
      end if
      if (ylog) then
         ymn = 10.**ymin
         ymx = 10.**ymax
      else
         ymn = ymin
         ymx = ymax
      end if
      return
C
      ENTRY SEELOC(IX,IY)
      call pgqpos(xx,yy)
      ix = x0+(xx-xmin)/xperpix
      iy = y0+(yy-ymin)/yperpix
      return
C
      ENTRY SETCOLOR(L)
      call pgsci(l)
      return
C
      ENTRY SCURSR(L, IX, IY)
      xx = xmin+xperpix*(float(ix)-x0)
      yy = ymin+yperpix*(float(iy)-y0)
      call pgcurse(xx, yy, ch)
      l = ichar(ch)
      if (l.gt.0) then
         ix = x0+(xx-xmin)/xperpix
         iy = y0+(yy-ymin)/yperpix
      end if
      return
C
      ENTRY VCURSR(L, X, Y)
      if (xlog) then
         if (x.le.0.) then
            xx = xmin
         else
            xx = alog10(x)
         end if
      else
         xx = x
      end if
      if (ylog) then
         if (y.le.0.) then
            yy = ymin
         else
            yy = alog10(y)
         end if
      else
         yy = y
      end if
      call pgcurse(xx,yy,ch)
      l = ichar(ch)
      if (l.ne.0) then
         if (xlog) then
            x = 10.**xx
         else
            x = xx
         end if
         if (ylog) then
            y = 10.**yy
         else
            y = yy
         end if
      end if
      return
C
      END   
C
      SUBROUTINE CHECK_TEXT(IN,OUT,KOUT,PLOT10_FLAG)
C
C     Decode special characters in PLOT10 string to PGPLOT
C
      IMPLICIT NONE
      CHARACTER IN*(*),OUT*(*)
      INTEGER KOUT
      LOGICAL PLOT10_FLAG
C
      CHARACTER SPECIAL(5)*1,CH2*2,SYMBOL(8)*6
      INTEGER I,J,K,IIN,JIN,IOUT,MODE,KMAX
      DATA SPECIAL/'<', '>', '?', '#', '&'/
      DATA SYMBOL/'\(845)', '\(847)', '\(840)', '\(846)', '\(841)',
     1            '\(842)', '\(843)', '\(852)'/
C
C     Omit trailing blanks and nulls
      do 2 jin=len(in),3,-1
         if (in(jin:jin).ne.' ' .and. in(jin:jin).ne.char(0)) go to 3
2     continue
C     Omit terminal '$'
3     if (in(jin:jin).eq.'$') jin = jin-1
C
      kmax = len(out)
      j = 0
      k = 0
      mode = 2
10    continue
         j = j+1
         if (plot10_flag) then
C           Look for PLOT10 special characters
            do 80 i=1,5
               if (in(j:j).eq.special(i)) then
                  if ((i.eq.1.or.i.eq.2).and.(i.ne.mode)) then
                     k = k+3
                     out(k-2:k) = '\f1'
                  else if (i.eq.5) then
C                    Convert next 2 characters to lower case for testing
                     ch2(1:1) = char(ior(ichar(in(j+1:j+1)),32))
                     ch2(2:2) = char(ior(ichar(in(j+2:j+2)),32))
                     if (ch2.eq.'ex') then
C                       End of superscript
                        j = j+2
                        k = k+2
                        out(k-1:k) = '\d'
                     else if (ch2(1:1).eq.'e') then
C                       Beginning of superscript
                        j = j+1
                        k = k+2
                        out(k-1:k) = '\u'
                     else if (ch2.eq.'lx') then
C                       End of subscript
                        j = j+2
                        k = k+2
                        out(k-1:k) = '\u'
                     else if (ch2(1:1).eq.'l') then
C                       Beginning of subscript
                        j = j+1
                        k = k+2
                        out(k-1:k) = '\d'
                     end if
                  end if
                  mode = i
                  go to 100
               end if
80          continue
C           Not a special case, just keep the character
            if (mode.eq.3 .or. mode.eq.4) then
C              Character is Greek, must be preceded with flag
               k = k+2
               out(k-1:k) = '\g'
            end if
            k = k+1
            if (mode.eq.1 .or. mode.eq.4) then
C              Character must be lower case
               out(k:k) = char(ior(ichar(in(j:j)),32))
            else
C              No case modification
               out(k:k) = in(j:j)
            end if
         else
C
C           Look for PGPLOT symbol numbers
            if (in(j:j).eq.'\') then
               i = ichar(in(j+1:j+1)) - ichar('0')
               if (i.ge.1 .and. i.le.8) then
                  j = j+1
                  k = k+6
                  out(k-5:k) = symbol(i)
                  go to 100
               end if
            end if
C           Not symbol number, copy unmodified character to output
            k = k+1
            out(k:k) = in(j:j)
         end if
C
100   if (j.lt.jin .and. k.lt.kmax) go to 10
      kout = min0(k,kmax)
      return
      END
C
      SUBROUTINE LINAXIS(XONE,XTWO,MAJOR,MINOR,FORMAT)
C
C  SCALE ENDS OF LINEAR AXIS TO ROUNDED NUMBERS
C
C    An axis is defined to include XONE and XTWO, with each end rounded to   
C  be an integer number of units of the form (1, 2, or 5) x 10**n.  The
C  resulting number MAJOR of major divisions will be between 4 and 10, with
C  MINOR minor divisions in each step.  A character string (FORMAT = 'Fww.dd')
C  is generated to use as the format for labeling the axis.
C
C  P. A. Seeger, Los Alamos National Laboratory, May 24, 1986
C    Modified to prevent log of X=0., Aug. 31, 1986
C    Use E instead of F format if width > 10, Nov. 21, 1987
C    Guard against axis ends > 10**38, Feb. 5, 1990
C
C    No Externals
C
      IMPLICIT NONE
      REAL*4 XONE,XTWO,XMIN,XMAX,AXLEN,UNIT
      INTEGER MAJOR,MINOR,IPOWR,NN,IW,ID
      CHARACTER*6 FORMAT
C
      XMIN = AMAX1(-0.7E38,AMIN1(XONE,XTWO))
      XMAX = AMIN1( 0.7E38,AMAX1(XONE,XTWO))
      IF (XMAX.EQ.XMIN) THEN
         IF (XMAX.LT.0.) THEN
            XMAX = 0.
         ELSE IF (XMIN.GT.0.) THEN
            XMIN = 0.
         ELSE 
            XMAX = 0.001
         END IF
      END IF
      AXLEN = XMAX-XMIN
      IPOWR = NINT(ALOG10(AXLEN)-1.05)
      UNIT = 10.**IPOWR
      NN = AXLEN/UNIT
      IF (NN.GE.15) THEN
         UNIT=UNIT*5.
         MINOR=1
         IF (NN.LE.30) MINOR=5
      ELSE IF (NN.GE.7) THEN
         UNIT=UNIT*2.
         MINOR=2
         IF (NN.LE.10) MINOR=4
      ELSE
         MINOR=2
         IF (NN.LE.4) MINOR=5
      END IF
C
      IF (XMIN.GE.0.) THEN
         XMIN = UNIT*AINT((XMIN+0.01*AXLEN)/UNIT)
      ELSE
         XMIN = UNIT*AINT((XMIN+0.01*AXLEN)/UNIT-1.)
      END IF
      IF (XMAX.GE.0.) THEN
         XMAX = UNIT*AINT((XMAX-0.01*AXLEN)/UNIT+1.)
      ELSE
         XMAX = UNIT*AINT((XMAX-0.01*AXLEN)/UNIT)
      END IF
C
      MAJOR = NINT((XMAX-XMIN)/UNIT)
      IF (XTWO.GT.XONE) THEN
         XONE = XMIN
         XTWO = XMAX
      ELSE
         XONE = XMAX
         XTWO = XMIN
      END IF
C
      ID = MAX0(0,-IPOWR)
      IF (XMIN.LT.0.) XMIN=-10.*XMIN
      IF (XMAX.LT.0.) XMAX=-10.*XMAX
      IPOWR = ALOG10(AMAX1(XMIN,XMAX,1.))+0.001
      IW = IPOWR+ID+2
      IF (IW.LT.10) THEN
         WRITE (FORMAT,100) IW,ID
100      FORMAT ('F',I2,'.',I2)
      ELSE
         FORMAT = 'E 9. 2'
      END IF
C
      RETURN
      END
C
      SUBROUTINE LOGAXIS(XONE,XTWO,MAJOR,MINOR,UNIT,FORMAT)
C
C  SCALE ENDS OF LOGARITHMIC AXIS TO ROUNDED NUMBERS
C    
C    An axis is defined to include XONE and XTWO, with each end rounded to   
C  be (1, 2, 3, or 5) x some power of 10.  The resulting axis spans MAJOR
C  decades, including partial decades at one or both ends.  The lowest decade
C  on the axis has MINOR units of size UNIT; if it is a full decade, MINOR = 9
C  and UNIT = min(X).  A character string (FORMAT = 'Fww.dd') is generated to
C  use as the format for labeling the axis.
C
C  P. A. Seeger, Los Alamos National Laboratory, May 20, 1986
C      Modified to avoid log(X) when X.LE.0., June 17, 1986
C      Corrected FORMAT when XMAX<1., June 19, 1987
C      Use E instead of F format when width > 10, Nov. 21, 1987
C      Guard against axis ends > 10**38, Feb. 5, 1990
C
C    No Externals
C
      IMPLICIT NONE
      REAL*4 XONE,XTWO,UNIT,XMIN,XMAX,XLOG,LOG2,LOG3,LOG5,DELTA
      INTEGER MAJOR,MINOR,IMIN,IMAX,IW,ID
      CHARACTER*6 FORMAT
      PARAMETER (LOG2=0.3010300,LOG3=0.4771213,LOG5=0.6989700)
C
      XMIN = AMIN1(XONE,XTWO)
      XMAX = AMAX1(XONE,XTWO)
      IF (XMAX.LE.0.) XMAX = 1.
      IF (XMIN.LE.0.) THEN
         XMIN = XMAX/1000.
      ELSE IF (XMAX.GE.1.E38) THEN
         XMAX = XMIN*1.E6
      ELSE
         XMIN = AMAX1(XMIN,XMAX/1.E24)
      END IF
      IF (XMIN.EQ.XMAX) THEN
         XMAX = XMAX*2.
         XMIN = XMIN/2.
      END IF
      DELTA = 0.01*ALOG10(XMAX/XMIN)
C
C     Look at small end of axis first
      XLOG = ALOG10(XMIN)+DELTA
      IMIN = NINT(XLOG-0.5)
      XMIN = 10.**IMIN
      UNIT = XMIN
      XLOG = XLOG-FLOAT(IMIN)
      IF (XLOG.GT.LOG5) THEN
         XMIN = 5.*XMIN
         MINOR = 5
      ELSE IF (XLOG.GT.LOG3) THEN
         XMIN = 3.*XMIN
         MINOR = 7
      ELSE IF (XLOG.GT.LOG2) THEN
         XMIN = 2.*XMIN
         MINOR = 8
      ELSE
         MINOR = 9
      END IF
      ID = MAX0(0,-IMIN)
C
C     Now do "same" thing at larger end
      XLOG = ALOG10(XMAX)-DELTA
      IMAX = NINT(XLOG-0.5)
C     Compute total number of decades included
      MAJOR = IMAX-IMIN+1
      XMAX = 10.**IMAX
      XLOG = XLOG-FLOAT(IMAX)
      IF (XLOG.LT.LOG2) THEN
         XMAX = 2.*XMAX
      ELSE IF (XLOG.LT.LOG3) THEN
         XMAX = 3.*XMAX
      ELSE IF (XLOG.LT.LOG5) THEN
         XMAX = 5.*XMAX
      ELSE
         XMAX = 10.*XMAX
         IMAX = IMAX+1
      END IF
      IF (MAJOR.EQ.1) MINOR=NINT((XMAX-XMIN)/UNIT)
C
      IW = MAX0(0,IMAX)+ID+2
      IF (IW.LT.10) THEN
         WRITE (FORMAT,100) IW,ID
100      FORMAT ('F',I2,'.',I2)
      ELSE
         FORMAT = 'E 9. 2'
      END IF
      IF (XONE.LT.XTWO) THEN
         XONE = XMIN
         XTWO = XMAX
      ELSE
         XONE = XMAX
         XTWO = XMIN
      END IF
C
      RETURN
      END
C
C***********************************************************************
C
      SUBROUTINE CONTOUR(Z,NRZ,X,NX,Y,NY,CV,NCV,LINE,ZMAX,BITMAP)
C
C  DRAW CONTOURS THROUGH EQUAL VALUES IN AN ARRAY
C
C  From Collected Algorithms from ACM #531 "Contour Plotting [J6]"
C   by: William V. Snyder, 1978
C  Copied from GSAS; transferred to IBM-PC, Dec. 8, 1987 (P.A.Seeger)
C  Added X, Y, and LINE to calling sequence, Dec. 9, 1987 (PAS) 
C  Restructured, closer to "standard", Dec. 11, 1987 (PAS)
C  Changed I and J to II and JJ in innermost loop, Dec. 13, 1987 (PAS)
C  Revised scan pattern from spiral to linear, Aug. 13, 1988 (PAS)
C
C  Externals
C    DASHA     FILL0     LGETMARK  MOVEA
C
C  Arguments in calling sequence:
C    Z         R(NRZ,*) Input    Array of values to be contoured; nodes must
C                                lie on a topologically rectangular grid.
C    NRZ       I        Input    Number of rows declared for array Z
C    X         R(*)     Input    Values of X at grid points of array Z
C    NX        I        Input    Limit for 1st subscript of Z and X
C    Y         R(*)     Input    Values of Y at grid points of array Z
C    NY        I        Input    Limit for 2nd subscript of Z and Y
C    CV        R(*)     Input    Values of contour levels
C    NCV       I        Input    No. of contour levels
C    LINE      I(*)     Input    Line style for each contour level
C    ZMAX      R        Input    Maximum Z to be considered; grid lines at
C                                a node with value above ZMAX are excluded
C    BITMAP    I(*)     Scratch  Work area of size (2*NX*NY*NCV+7)/8 bytes
C
      IMPLICIT NONE
      INTEGER NRZ,NX,NY,NCV,LINE(*)
      REAL*4 Z(NRZ,*),X(*),Y(*),CV(*),ZMAX
      INTEGER*4 BITMAP(*) 
C
      LOGICAL LGETMARK
C
C  Local variables in CONTOUR:
C    CVAL      R        Contour-line value being traced
C    DELZ      R        Change in Z when moving 1 cell right or up
C    IADD      I        Bit address in BITMAP, starting at zero
C    IBORDER   I        Edge of current cell which is on a border
C    ICV       I        Index of contour line being traced
C    IEDGE     I        Edge of new cell where contour line enters
C    IFLAG     I        1=continue, 2=start at boundary, 3=start in interior,
C                       4=end at boundary, 5=close contour, 6=none found yet
C    I,J       I        Subscripts for search
C    IJ        I(2)     Equivalent to I,J
C    IJMAX     I(2)     Local copies of NX and NY
C    I1(2),I2(2),I3(6)  Used for subscript computations
C    II,JJ     I        Cell with continuation of contour line being traced
C    INDEX     I        L-1 + 2*(I-1 + NX*(J-1))
C    K         I        Index of cell edges: 1=bottom, 2=left, 3=top, 4=right
C    KS        I        Next cell boundary to cross
C    L,LL      I        Orientation flags: 1 is horizontal line, 2 vertical
C    LBORDER   L        Flag to show that only border lines are to be done
C    NI        I        Number of edges of cell which contour crosses
C    XINT      R(4)     Intersections of contour with edges of cell,
C                          in order bottom, left, top, right
C    XX,YY     R        Plotting coordinates
C    Z1,Z2     R        Smaller and larger values at segment ends
C    Z3,Z4     R        Values at ends of segment for continuation of contour
C
      REAL*4 XINT(4),Z1,Z2,Z3,Z4,DELZ,CVAL,XX,YY
      INTEGER I1(2),I2(2),I3(6),I,J,IJ(2),IJMAX(2),K,L,LL,II,JJ,ICV,
     , IBORDER,IEDGE,IFLAG,NI,KS
      INTEGER*4 INDEX,MAXINDX,NBITS,IADD
      LOGICAL LBORDER
      EQUIVALENCE (IJ(1),I),(IJ(2),J)
      DATA I1/1,0/,  I2/1,-1/,  I3/1,0,0,1,1,0/
C
      IJMAX(1) = NX
      IJMAX(2) = NY
C  Clear bit map
      NBITS = 2*NX*NY*NCV
      CALL FILL0(BITMAP,NBITS)
      IFLAG = 6
C
C  Search every cell in rectangular array for a line segment such that:
C    1. the end points are not excluded because Z > ZMAX
C    2. current contour <= Z at one end and > Z at other end
C    3. no mark has been recorded for this contour on this segment
C
C  Set start at lower left corner of array; do borders first, then interior
C
      LBORDER = .TRUE.
      MAXINDX = 2*NX*NY-2
100   CONTINUE
      DO 500 INDEX=0,MAXINDX
         L = MOD(INDEX,2)+1
         I = MOD(INDEX/2,NX)+1
         J = (INDEX/2)/NX+1
         IF (Z(I,J).LE.ZMAX) THEN
C           Node itself is within non-excluded range
            II = I+I1(L)
            JJ = J+I1(3-L)
            IF (II.LE.NX .AND. JJ.LE.NY .AND. Z(II,JJ).LE.ZMAX) THEN
C              Both ends of grid line within range; test if "border" segment
               IBORDER = 0
               IF (IJ(3-L).EQ.1 .OR.
     1             Z(I-I1(3-L),J-I1(L)).GT.ZMAX .OR.
     2             Z(I+I2(L),J+I2(3-L)).GT.ZMAX)   IBORDER = 1
               IF (IJ(3-L).GE.IJMAX(3-L) .OR.
     1             Z(I+I1(3-L),J+I1(L)).GT.ZMAX .OR.
     2             Z(I+1,J+1).GT.ZMAX)     IBORDER = IBORDER+2 
C              1st time, do ONLY borders (including edges of exclusions)
               IF (IBORDER.NE.3 .AND. (LBORDER.EQV.(IBORDER.NE.0))) THEN
C                 Examine this line segment this pass
                  Z1 = AMIN1(Z(I,J),Z(II,JJ))
                  Z2 = AMAX1(Z(I,J),Z(II,JJ))
                  DELZ = Z(II,JJ)-Z(I,J)
                  DO 400 ICV=1,NCV
C                    Test for all possible contours crossing this grid line;
C                     first check if already done, and set bit to show done
                     IADD = ICV-1 + NCV*INDEX
                     IF (.NOT.LGETMARK(BITMAP,IADD) .AND. 
     1                   CV(ICV).GT.Z1 .AND. CV(ICV).LE.Z2) THEN
C                       Found one we haven't done yet!!!  Interpolate.
                        CVAL = CV(ICV)
C                       Decide which edge we are entering cell from
                        IEDGE = L
                        IF (IBORDER.EQ.2) IEDGE = IEDGE+2
                        XINT(IEDGE) = (CVAL-Z(I,J))/DELZ
C                       Move "pen" to starting point of contour
                        XX = X(I)+XINT(IEDGE)*(X(II)-X(I))
                        YY = Y(J)+XINT(IEDGE)*(Y(JJ)-Y(J))
                        IFLAG = 3
                        IF (LBORDER) IFLAG = 2
                        CALL MOVEA(XX,YY)
C
C  Follow this contour until it hits boundary or closes on itself
                        II = I
                        JJ = J
                        IFLAG = 1
C                       "DO WHILE (IFLAG.LT.4)"
200                     CONTINUE
C                          If haven't moved to next cell yet, do so
                           IF (IEDGE.EQ.3) JJ = JJ-1
                           IF (IEDGE.EQ.4) II = II-1
                           NI = 1
                           DO 300 K = 1,4
C                             Test interpolation on other 3 edges
                              IF (K.NE.IEDGE) THEN
                                 Z3 = Z(II+I3(K),JJ+I3(K+1))
                                 Z4 = Z(II+I3(K+1),JJ+I3(K+2))
                                 IF (CVAL.GT.AMIN1(Z3,Z4) .AND. 
     .                               CVAL.LE.AMAX1(Z3,Z4)) THEN
C                                   The contour also crosses this edge
                                    IF (K.EQ.1 .OR. K.EQ.4) THEN
                                       XINT(K) = (CVAL-Z4)/(Z3-Z4)
                                    ELSE
                                       XINT(K) = (CVAL-Z3)/(Z4-Z3)
                                    END IF
C                                   Count how many crossings
                                    NI = NI+1
                                    KS = K
                                 END IF
                              END IF
300                        CONTINUE
C
                           IF (NI.NE.2) THEN
C  The contour crosses all four edges of the cell being examined. Choose the 
C  lines top-to-left and bottom-to-right if the interpolation point on the top 
C  edge is less than the interpolation point on the bottom edge. Otherwise, 
C  choose the other pair. This method produces the same result if the axes are 
C  reversed. The contour may close at any edge, but must not cross itself 
C  inside any cell.
                              IF (XINT(3).GE.XINT(1)) THEN
                                 KS = 3-IEDGE
                                 IF (KS.LE.0) KS = KS+4
                              ELSE
                                 KS = 5-IEDGE
                              END IF
                           END IF
C
C  Determine if the contour will close or run into a boundary at edge KS of the 
C  current cell.
                           IF (KS.LE.2) THEN
                              LL = KS
                              IEDGE = KS+2
                           ELSE
C                             Must move to adjacent cell before test for closure
                              II = II+I3(KS)
                              JJ = JJ+I3(KS+2)
                              LL = KS-2
                              IEDGE = KS-2
                           END IF
                           IADD = ICV-1+NCV*(LL-1+2*(II-1+NX*(JJ-1)))
                           IF (LGETMARK(BITMAP,IADD)) THEN
C                             We've already been here; contour has closed 
                              IFLAG = 5
                           ELSE IF (LL.EQ.2.AND.(II.EQ.1.OR.II.GE.NX)
     1                     .OR. LL.EQ.1.AND.(JJ.EQ.1.OR.JJ.GE.NY)) THEN
C                             Segment is actual boundary of plot
                              IFLAG = 4
                           ELSE IF (Z(II-I1(3-LL),JJ-I1(LL)).GT.ZMAX 
     1                     .OR. Z(II+I2(LL),JJ+I2(3-LL)).GT.ZMAX
     2                     .OR. Z(II+I1(3-LL),JJ+I1(LL)).GT.ZMAX
     3                     .OR. Z(II+1,JJ+1).GT.ZMAX) THEN
C                             Segment is boundary of an excluded cell
                              IFLAG = 4
                           END IF
C
C                          Draw piece of contour
                           XINT(IEDGE) = XINT(KS)
                           IF (LL.EQ.1) THEN
                              XX = X(II)+XINT(IEDGE)*(X(II+1)-X(II))
                              YY = Y(JJ)
                           ELSE
                              XX = X(II)
                              YY = Y(JJ)+XINT(IEDGE)*(Y(JJ+1)-Y(JJ))
                           END IF
                           CALL DASHA(XX,YY,LINE(ICV))
C
                        IF (IFLAG.LT.4) GO TO 200
C                       Reset II and JJ before looking for next contour
                        II = I+I1(L)
                        JJ = J+I1(3-L)
                     END IF
400               CONTINUE
               END IF
            END IF
         END IF
500   CONTINUE
C
      IF (.NOT.LBORDER) RETURN
      LBORDER = .FALSE.
      GO TO 100
      END
C  
C  CONTOUR SUBROUTINES FOR BIT MANIPULATION
C
      SUBROUTINE FILL0(BITMAP,N)
C  Fills entire BITMAP with zeros
      IMPLICIT NONE
      INTEGER N,I,LOOP
      INTEGER*4 BITMAP(*)
C
      LOOP = (N-1)/32+1
      DO 10 I=1,LOOP
         BITMAP(I) = 0
10    CONTINUE
      RETURN
      END
C
      LOGICAL FUNCTION LGETMARK(BITMAP,N)
C  Tests bit in BITMAP, and then sets it to one
      IMPLICIT NONE
      INTEGER N,NWORD,NBIT
      INTEGER*4 BITMAP(*)
      LOGICAL BTEST
C                                        
      NWORD = N/32+1
      NBIT = MOD(N,32)
      LGETMARK = BTEST(BITMAP(NWORD),NBIT)
      IF (.NOT.LGETMARK) BITMAP(NWORD) = IBSET(BITMAP(NWORD),NBIT)
      RETURN
      END
C
C******************************************************************************
C
      SUBROUTINE CURVE(X,Y,N,LDASH,ISYM)
C
C  PLOT N POINTS AT (X(I),Y(I)), USING SYMBOL CORRESPONDING TO ISYM, AND
C    CONNECTING WITH LINE DESCRIBED BY LDASH.
C
C  P. A. Seeger, Los Alamos National Laboratory, May 24, 1986
C
      INTEGER N,LDASH,ISYM,I
      REAL*4 X(*),Y(*)
C
C  Externals
C     DASHA     MOVEA     SYMBOL
C
      CALL MOVEA(X(1),Y(1))
      DO 100 I=1,N
         CALL DASHA(X(I),Y(I),LDASH)
         CALL SYMBOL(X(I),Y(I),ISYM)
100   CONTINUE
C
      RETURN
      END
C
C******************************************************************************
C
      SUBROUTINE ERRBARS(X,Y,DY,N,INC,LOGY)
C
C  PLOT ERROR BARS ON EVERY INCth POINT OF AN ARRAY OF N POINTS 
C    (X(I),Y(I)+-DY(I)).  LOGY=1 IF ORDINATE IS LOGARITHMIC.
C
C  P. A. Seeger, Los Alamos National Laboratory, Nov. 24, 1987
C
      IMPLICIT NONE
      INTEGER N,INC,LOGY,I
      REAL*4 X(*),Y(*),DY(*),D,R,Y1,Y2
C
C  Externals
C     DRAWA     MOVEA     
C
      DO 100 I=1,N,INC
         D = ABS(DY(I))
         IF (LOGY.EQ.1) THEN
            IF (Y(I).LE.0.) THEN
               Y1 = 2.E-38
               Y2 = 2.E-38
            ELSE IF (D.LT.0.05*Y(I)) THEN
               Y1 = Y(I)-D
               Y2 = Y(I)+D
            ELSE IF (D.LT.88.*Y(I)) THEN
               R = EXP(D/Y(I))
               Y1 = Y(I)/R
               Y2 = Y(I)*R
            ELSE
               Y1 = 2.E-38
               Y2 = 1.E+38
            END IF
         ELSE
            Y1 = Y(I)-D
            Y2 = Y(I)+D
         END IF
         CALL MOVEA(X(I),Y1)
         CALL DRAWA(X(I),Y2)
100   CONTINUE
C
      RETURN
      END