****************************************************************************
c
c      SHAREWARE NOTICE              Copyright (C) D.I. Hoyer, 1990/1991.
c     ==================             GRAPHLIB.for  v3.0
c
c     Please register by sending US$35-00 or Aus$40-00 to the address below
c     if you find these Fortran subroutines useful. Registered users will
c     receive the latest version, plus technical support on queries related
c     to these routines.
c
c                    David I Hoyer
c                    P.O. Box 1743
c                    Macquarie Centre NSW 2113
c                    AUSTRALIA
c
c     This subroutine library is a shareware product. Copies of the original
c     unmodified programs and manual on disk may be made and distributed as 
c     required, as long as they are not charged for. You may not modify the 
c     source code or manual except for your personal use. 
c
c     Requirements  :  ANSI Fortran 77, full language.
c
****************************************************************************
c
c     No responsibility is accepted for any errors in this software, 
c     or for any loss or damage resulting from using it.
c
****************************************************************************

      SUBROUTINE PRTGRF(IOFF,MMAXX,MMAXY,LOHI,IGRAPH)
*
* Print graph IGRAPH on the dot matrix printer.
*
* Set the following character constants after referring to the printer manual:
*
* LSP7 and LSP12 = Set line spacing to 7/72" and 12/72" (1/6") respectively.
* LSP1  = Set line spacing to 1/216"  [or  1/144"]   (for hi-res plotting)
* LSP20 = Set line spacing to 20/216" [or 13/144"]   ( "    "       "    )
* GFXON = Set graphics mode on, with number of dot columns across page
*             Epson, Star NX-10 : <ESC> 'K'...  =  60 dpi (lo-res)
*                                 <ESC> 'L'...  = 120 dpi (hi-res)
*
* LOHI : Two options are available, viz low and high density.
* 1= Low density = 72 dots per inch vertical and 60 horizontal. Each
*                  row of integers is printed as two lines of 7 rows
*                  each on the printer. 8"x10" = 41k. IGRAPH(480,43)
* 2= High   "    = 144 dpi vert and 120 horiz. This is printed in two
*                  interleaved rows, with 1/216" or 1/144" line feed.
*                  8"x10" requires about 164 kbytes : IGRAPH(960,86)
*
* BUT NOTE: Some printers have a basic vertical line spacing of 60 & 180 dpi
* --------  rather than the 72 discussed above. This is indicated if the
*    printer manual refers to vertical line settings as n/60 or n/180 dpi
*    instead of n/72 or n/216. If this is the case, try the following:
*    LSP7  - set line spacing to 7/60"
*    LSP1  - set line spacing to 1/180"
*    LSP20 - set line spacing to 20/180"
* Then make sure your main program calls PREP(... with DPIV = 60.
* (for 60 dots per inch vertical)
*
* If the printer cannot do line feeds of 1/216" to 1/144" then only lo-res
* graphs can be printed.
*
*
      INTEGER*2 IGRAPH
      CHARACTER*5 LSP12
      CHARACTER*4 GFXON
      CHARACTER*3 LSP1,LSP20
      CHARACTER*2 LSP7
      CHARACTER BLANK*100,LINE1*1200,LINE2*1200
      DIMENSION IGRAPH(MMAXX,MMAXY)
      DATA BLANK/'
     $                                              '/
*
* This subroutine is set up for IBM, Epson, Star type dot matrix printers
* which have a smallest line feed of 1/216". Some printers have a smallest
* line feed of 1/144", in which case you should re-define LSP20 to :
*      LSP20 = CHAR(27)//'3'//CHAR(13)
*
* If you have a printer which uses different graphics commands, or if the
* graph doesn't print correctly, consult the printer manual and change
* the appropriate parameters for LSP7, LSP12, LSP1, LSP20, GFXON (and also
* remember to set the CHARACTER*n statements above to the correct sizes):
* LSP7 sometimes gives trouble - if the lo-res plots appear "expanded"
* try   LSP7 = CHAR(27)//'A'//CHAR(7)
*  or   LSP7 = CHAR(27)//'A'//CHAR(7)//CHAR(27)//'2'
* and change CHARACTER*2 LSP7 to   CHARACTER*3 or 5 respectively.
*
      LSP7 = CHAR(27)//'1'
      LSP12 = CHAR(27)//'A'//CHAR(12)//CHAR(27)//'2'
      LSP1 = CHAR(27)//'3'//CHAR(1)
      LSP20 = CHAR(27)//'3'//CHAR(20)
*
      IOFF = MAX(1,IOFF)
      I2 = MMAXX/256
      I1 = MMAXX - 256*I2
      IF(LOHI.EQ.1) THEN
*       Set 60 dots per inch across the page (lo-res)
        GFXON = CHAR(27)//'K'//CHAR(I1)//CHAR(I2)
      ELSE
*       Set 120 dots per inch across the page (hi-res)
        GFXON = CHAR(27)//'L'//CHAR(I1)//CHAR(I2)
      ENDIF
      WRITE(6,303) LSP7
 303  FORMAT(1X,A5)
      DO 10 IROW=MMAXY, 1, -1
        DO 30 ICOL=1, MMAXX
          ICH1 = IGRAPH(ICOL,IROW)/128
          ICH2 = 2*(IGRAPH(ICOL,IROW) - 128*ICH1)
          ICH1 = 2*ICH1
          IF(ICH1.EQ.26) ICH1=18
          IF(ICH2.EQ.26) ICH2=18
          LINE1(ICOL:ICOL) = CHAR(ICH1)
          LINE2(ICOL:ICOL) = CHAR(ICH2)
  30    CONTINUE
        IF(LOHI.EQ.1) THEN
          WRITE(6,101) BLANK(1:IOFF),GFXON,LINE1(1:MMAXX),BLANK(1:IOFF),
     $     GFXON,LINE2(1:MMAXX)
 101      FORMAT(1X,A,A5,A/1X,A,A5,A)
        ELSE
          WRITE(6,404) BLANK(1:IOFF),GFXON,LINE1(1:MMAXX),LSP1
          WRITE(6,404) BLANK(1:IOFF),GFXON,LINE2(1:MMAXX),LSP20
 404      FORMAT(1X,A,A5,A,A3)
        ENDIF 
  10  CONTINUE
* return the printer to 1/6" spacing
      WRITE(6,303) LSP12   
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE SETBIT(I,J)
*
* Set bit j in integer i. This method is marginally faster than using
* the IOR function which is available with some Fortran compilers.
*
      INTEGER*2 I
      DIMENSION MASK(0:14)
      DATA MASK/16384,8192,4096,2048,1024,512,256,128,64,32,16,8,4,2,1/
      IF(MOD(I,MASK(J-1)).LT.MASK(J)) I = I + MASK(J)
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE PREP(XMIN,YMIN,XMAX,YMAX,LOGX,LOGY,LOHI,IVH,
     $ XL,YL,BWLEFT,BWRGHT,BWTOP,BWBOTT,DPIV,DPIH)
*
* Prepare the input data for plotting.
* Calc graph sizes and limits, log conversions etc.
* This subroutine MUST be called before any plotting is done.
*
*  XMIN, XMAX = Min & max x values to define edges of plotting area.
*  YMIN, YMAX = Min & max y values to define edges of plotting area.
*  LOGX, LOGY = 0 to specify log scales on x or y-axes,
*              1 to specify linear scales.
*  LOHI   = 1 for Low resolution plot (about 4 times quicker than Hi-Res),
*           2 for Hi-Res plot.
*  IVH    = 1 for vertical (portrait) graph format,
*           2 for horizontal (landscape) format.
*  XL     = Length of x-axis, cm.
*  YL     = Length of y-axis, cm.
*  BWLEFT = Border width to left of axes, cm.
*  BWRGHT = Border width to right of axes, cm.
*  BWTOP  = Border width above axes, cm.
*  BWBOTT = Border width below axes, cm.
*  DPIH   = Horizontal dots per inch for lo-res plotting (normally 60)
*  DPIV   = Vertical dots per inch for lo-res plotting (normally 72 or 60)
*
*   The border widths allow room outside the graph axes for plotting labels,
*   titles etc. These add to the overall plotting area. Note that the total
*   plotting height and width (including borders) may not exceed the capacity
*   of IGRAPH, otherwise a run-time error will occur. The total plotting
*   area available depends on the dimension of IGRAPH in the main program,
*   and the number of dots per inch for the printer being used.
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN2,YMIN2,XMAX2,YMAX2,
     $   MAXX,MAXY,MMAXX,MMAXY,LOHI2,IVH2,LOGX2,LOGY2,DPIV2,DPIH2
*
*  XMINA..YMAXA = Min/max coordinate values of the edges of the plotting area.
*  They are derived from XMIN..YMAX and axis lengths and border widths.
*
* If log scales are specified, then take logs of axis details
      IF(LOGX.EQ.0) THEN
        XMIN = ALOG10(XMIN)
        XMAX = ALOG10(XMAX)
      ENDIF
      IF(LOGY.EQ.0) THEN
        YMIN = ALOG10(YMIN)
        YMAX = ALOG10(YMAX)
      ENDIF
* Allow extra space above, below, left & right, for labels and headings.
* The numbers refers to cm here.
      XMINA = XMIN - BWLEFT/XL*(XMAX-XMIN)
      XMAXA = XMAX + BWRGHT/XL*(XMAX-XMIN)
      YMINA = YMIN - BWBOTT/YL*(YMAX-YMIN)
      YMAXA = YMAX + BWTOP/YL*(YMAX-YMIN)
      IF(IVH.EQ.1) THEN
*------ Vertical format (portrait, or normal)
        WIDTH = XL + BWLEFT + BWRGHT
        HEIGHT = YL + BWTOP + BWBOTT
      ELSE
*------ If horizontal (landscape) format selected, transform width/height
        WIDTH = YL + BWTOP + BWBOTT
        HEIGHT = XL + BWLEFT + BWRGHT
      ENDIF
* Set the number of dots vert. and horiz. for required graph size
      IF(LOHI.EQ.1) THEN
        MAXX = INT(WIDTH*DPIH/2.54+1.5)
        MAXY = INT(HEIGHT*DPIV/2.54+1.5)
      ELSE
        MAXX = INT(WIDTH*2*DPIH/2.54+1.5)
        MAXY = INT(HEIGHT*2*DPIV/2.54+1.5)
      ENDIF
      MMAXX = MAXX
      MMAXY = MAXY/14+1
      IF(IVH.EQ.2) THEN
        I = MAXX
        MAXX = MAXY
        MAXY = I
      ENDIF
      XMIN2 = XMIN
      YMIN2 = YMIN
      XMAX2 = XMAX
      YMAX2 = YMAX
      LOHI2 = LOHI
      IVH2 = IVH
      LOGX2 = LOGX
      LOGY2 = LOGY
      DPIV2 = DPIV
      DPIH2 = DPIH
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE CLRBOX(X1,Y1,X2,Y2,LTYP,INOUT,IGRAPH)
*
* Clear a rectangular area of the graph, and draw a box around
* the cleared area as an additional option.
* See also subroutine  LGDBOX(...
*
* NOTE: This sub-program contains the IAND(-,-) statement which is not
*       standard Fortran-77. However, most F77 compilers recognise it.
*       If yours does not recognise it you will have to delete this
*       subroutine, and the reference to it in the main program.
*
* (X1,Y1) = coordinates of top left corner of box to be cleared
* (X2,Y2) = coordinates of bottom right corner of box to be cleared
* LTYP    = Line type for border around box (0=none, 1=solid etc.)
* IGRAPH  = The graph memory array
*
* Useful for clearing an area which may contain grid lines, or
* other data, before drawing a legend table or other text.
*

      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      DIMENSION IGRAPH(MMAXX,MMAXY), MASK1(14), MASK2(14)
      DIMENSION MASKH1(14), MASKH2(14)
      DATA MASK2/16382,16380,16376,16368,16352,16320,16256,16128,
     $           15872,15360,14336,12288,8192,0/
      DATA MASK1/1,3,7,15,31,63,127,255,511,1023,2047,4095,8191,16383/
      DATA MASKH2/16382,16254,16252,15996,15992,15480,15472,14448,
     $            14432,12384,12352,8256,8192,0/
      DATA MASKH1/1,129,131,387,391,903,911,1935,1951,3999,4031,
     $            8127,8191,16383/

* Check boundary restrictions
      IF(INOUT.EQ.0) THEN
        X1 = AMAX0(XMIN,AMIN0(X1,XMAX))
        X2 = AMAX0(XMIN,AMIN0(X2,XMAX))
        Y1 = AMAX0(YMIN,AMIN0(Y1,YMAX))
        Y2 = AMAX0(YMIN,AMIN0(Y2,YMAX))
      ELSE
        X1 = AMAX0(XMINA,AMIN0(X1,XMAXA))
        X2 = AMAX0(XMINA,AMIN0(X2,XMAXA))
        Y1 = AMAX0(YMINA,AMIN0(Y1,YMAXA))
        Y2 = AMAX0(YMINA,AMIN0(Y2,YMAXA))
      ENDIF

* Calc integer values of the box coordinates.
      IF(IVH.EQ.1) THEN
        IX1 = INT((X1-XMINA)/(XMAXA-XMINA)*FLOAT(MAXX-1)+1.5)
        IY1 = INT((Y1-YMINA)/(YMAXA-YMINA)*FLOAT(MAXY-1)+0.5)
        IX2 = INT((X2-XMINA)/(XMAXA-XMINA)*FLOAT(MAXX-1)+1.5)
        IY2 = INT((Y2-YMINA)/(YMAXA-YMINA)*FLOAT(MAXY-1)+0.5)
      ELSE
        IX1 = MAXY - INT((Y1-YMINA)/(YMAXA-YMINA)*FLOAT(MAXY-1)+0.5)
        IY1 = INT((X1-XMINA)/(XMAXA-XMINA)*FLOAT(MAXX-1)+0.5)
        IX2 = MAXY - INT((Y2-YMINA)/(YMAXA-YMINA)*FLOAT(MAXY-1)+0.5)
        IY2 = INT((X2-XMINA)/(XMAXA-XMINA)*FLOAT(MAXX-1)+0.5)
      ENDIF
* Sort ix1, ix2 into order, and sort iy1, iy2 into order
      IF(IX2.LT.IX1) THEN
        I = IX1
        IX1 = IX2
        IX2 = I
      ENDIF
      IF(IY2.LT.IY1) THEN
        I = IY1
        IY1 = IY2
        IY2 = I
      ENDIF
* Prepare for clearing
      JY2 = IY2/14
      JY1 = (IY1-2)/14+2
      JJY2 = IY2 - 14*JY2
      JJY1 = IY1 - 14*(JY1-2)
* For ix1 to ix2 clear whole integers jy1 to jy2
      IF(JY1.LE.JY2) THEN
        DO 10 I=IX1,IX2
          DO 20 J=JY1,JY2
            IGRAPH(I,J) = 0
  20      CONTINUE
  10    CONTINUE
      ENDIF
* For ix1 to ix2 clear partial integer jy2+1
      IF(JY2.LT.MAXY.AND.JJY2.GT.0.AND.JJY2.LT.15) THEN
        IF(LOHI.EQ.1) J = MASK2(JJY2)
        IF(LOHI.EQ.2) J = MASKH2(JJY2)
        DO 30 I=IX1,IX2
          IGRAPH(I,JY2+1) = IAND(IGRAPH(I,JY2+1),J)
  30    CONTINUE
      ENDIF
* For ix1 to ix2 clear partial integer jy1-1
      IF(JY1.GT.1.AND.JJY1.LT.15.AND.JJY1.GT.0) THEN
        IF(LOHI.EQ.1) J = MASK1(JJY1)
        IF(LOHI.EQ.2) J = MASKH1(JJY1)
        DO 40 I=IX1,IX2
          IGRAPH(I,JY1-1) = IAND(IGRAPH(I,JY1-1),J)
  40    CONTINUE
      ENDIF
* Now draw the border
      IF(LTYP.GT.0) THEN
        CALL LINE(LTYP,X1,Y2,X2,Y2,IGRAPH,NPP,0,INOUT)
        CALL LINE(LTYP,X2,Y2,X2,Y1,IGRAPH,NPP,0,INOUT)
        CALL LINE(LTYP,X2,Y1,X1,Y1,IGRAPH,NPP,0,INOUT)
        CALL LINE(LTYP,X1,Y1,X1,Y2,IGRAPH,NPP,0,INOUT)
      ENDIF
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE LGDBOX(LGDPOS,LINES,ICHRS,LGDSIZ,LTYP,INOUT,IGRAPH)
*
* Clear a rectangular area of the graph for a legend table, and draw a
* box around the cleared area as an additional option.
* See also subroutine  CLRBOX(...
*
* LGDPOS = Legend table position (1 to 8)
* LINES  = Number of lines (entries) in the legend table
* ICHRS  = Max number of characters for any legend
*             (to determine the width of the legend box)
* LGDSIZ = Size of legend text
* LTYP   = Line type for border around box (0=none, 1=solid etc.)
* INOUT  = 0 to clip at the axis boundary, 1 to allow overlapping the axes
* IGRAPH = The graph memory array
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      DIMENSION IGRAPH(MMAXX,MMAXY)
      ISIZE = MAX0(1,(LGDSIZ*LOHI+1)/2)
      DX = (XMAXA-XMINA)/FLOAT(MAXX)*FLOAT(ISIZE)
      DY = (YMAXA-YMINA)/FLOAT(MAXY)*FLOAT(ISIZE)
      IF(LGDPOS.LT.1.OR.LGDPOS.GT.8) LGDPOS = 1
      YYY = FLOAT(LINES*13+6)*DY
      Y1 = YMAX - 5*DY
      Y2 = Y1 - YYY
      LGDP = LGDPOS
      IF(LGDPOS.GT.4) THEN
        Y1 = YMIN+5*DY
        Y2 = Y1 + YYY
        LGDP = LGDPOS - 4
      ENDIF
      X1 = XMIN + 3.*DX + FLOAT(LGDP-1)*(XMAX-XMIN)/4.
      X2 = X1 + DX*FLOAT(ICHRS)*6 + 33*DX*FLOAT(LOHI)/FLOAT(ISIZE)
      CALL CLRBOX(X1,Y1,X2,Y2,LTYP,INOUT,IGRAPH)
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE AXES(NDIVX,IGRDX1,NSDIVX,IGRDX2,NDPX,XLABEL,
     $                NDIVY,IGRDY1,NSDIVY,IGRDY2,NDPY,YLABEL,IORIYV,
     $           TITLE,ISZVAL,VALPOS,ISZXYL,XLBPOS,YLBPOS,ISZTTL,TTLPOS,
     $           JUSTTL,LGDSIZ,LGDPOS,LGDLNS,LGDCHS,LGDTYP,IGRAPH)
*
* Draw a set of axes on linear or logarithmic scales.
* Size and position of labels and axis values can be specified.
* A cleared box for the legend table can be specified.
*
*  NDIVX  = No. of major divisions along x-axis, with axis values printed
*             (set to 0 for log scales, and divisions are calculated)
*  IGRDX1 = Type of grid lines on major x-axis divisions
*             (0=none, 1=solid, 2=dotted etc)
*  NSDIVX = No. of secondary divisions between each major division on x-axis
*  IGRDX2 = Type of grid lines on minor x-axis divisions (0=none, 1..5)
*  NDPX   = No. of decimal places for x-axis values (ignored for log scales)
*  XLABEL = The label or title of the x-axis
*  NDIVY...YLABEL = As for NDIVX...XLABEL, for y-axis
*  IORIYV = Orientation of y-axis values.
*             (0=numbers parallel to y-axis, 1=perpendicular to y-axis)
*  TITLE  = The graph title
*  ISZVAL = The text size for printing axis values
*  VALPOS = To adjust the distance between axis and centre of axis value.
*             (1=default, <1=closer to axis, >1=further from axis)
*  ISZXYL = The text size for printing x and y-axis labels
*  XLBPOS = To adjust the distance between x-axis and x-axis label.
*             1=default, <1=closer to axis, >1=further from axis.
*  YLBPOS = As for XLBPOS, but for y-axis. eg. If IORIYV=1 and the y-axis
*           values are several characters long, set YLBPOS>1 to shift the
*           label further away from the axis.
*  ISZTTL = The text size for printing the graph title
*  TTLPOS = To adjust the distance between axis and centre of graph title.
*             (>0 = Above top axis, <0 = below bottom axis
*             magnitude: 1=default, <1=closer to axis, >1=further from axis)
*  JUSTTL = Justification of graph title. (-1=left, 0=centre, 1=right)
*  LGDSIZ = Size of text for legend table (for drawing the legend box)
*  LGDPOS = Position of the legend box.
*             (1..4 = top left to top right, 5..8 = bottom left to bottom rt)
*  LGDLNS = Number of lines of text to be allocated for legend box
*             (set to 0 for no legend box)
*  LGDCHS = Max number of characters in a legend text (for the legend box)
*  LGDTYP = Line type for surrounding the legend box (0=none, 1=solid etc)
*  IGRAPH = The graph memory array
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      CHARACTER*80 TITLE, XLABEL, YLABEL, AXVAL, FMT
      DIMENSION IGRAPH(MMAXX,MMAXY)
      CALL LINE(1,XMIN,YMIN,XMAX,YMIN,IGRAPH,NPP,0,1)
      CALL LINE(1,XMIN,YMIN,XMIN,YMAX,IGRAPH,NPP,0,1)
      CALL LINE(1,XMIN,YMAX,XMAX,YMAX,IGRAPH,NPP,0,1)
      CALL LINE(1,XMAX,YMAX,XMAX,YMIN,IGRAPH,NPP,0,1)
* These are for spacing the labels
      ISIZE1 = MAX0(1,(ISZTTL*LOHI+1)/2)
      ISIZE2 = MAX0(1,(ISZXYL*LOHI+1)/2)
      ISIZE3 = MAX0(1,(ISZVAL*LOHI+1)/2)
      YY1 = 13*ABS(TTLPOS)*(YMAXA-YMINA)/FLOAT(MAXY)*FLOAT(ISIZE1)
      YY2 = 13*ABS(XLBPOS)*(YMAXA-YMINA)/FLOAT(MAXY)*FLOAT(ISIZE2)
      YY3 = 7*ABS(VALPOS)*(YMAXA-YMINA)/FLOAT(MAXY)*FLOAT(ISIZE3)
      XX2 = 13*ABS(YLBPOS)*(XMAXA-XMINA)/FLOAT(MAXX)*FLOAT(ISIZE2)
      XX3 = 7*ABS(VALPOS)*(XMAXA-XMINA)/FLOAT(MAXX)*FLOAT(ISIZE3)
* Plot graph heading
      X = (XMAX+XMIN)/2.
      IF(JUSTTL.LT.0) X = XMIN
      IF(JUSTTL.GT.0) X = XMAX
      IF(TTLPOS.GE.0) Y = YMAX + YY1
      IF(TTLPOS.LT.0) Y = YMIN - YY1 - YY2 - YY3
      CALL TEXT(X,Y,TITLE,1,ISZTTL,JUSTTL,IGRAPH,1)
* Plot x-axis and y-axis labels
      DO 30 IAX=0, 1
        IF(IAX.EQ.0) THEN
          NDIVS = AMAX0(1,NDIVX)
          NSDIVS = AMAX0(1,NSDIVX)
          LOGS = LOGX
          IGRDS1 = IGRDX1
          IGRDS2 = IGRDX2
          AXMIN = XMIN
          AXMAX = XMAX
          Y = YMIN - YY2 - YY3
          X = (XMAX+XMIN)/2.
          CALL TEXT(X,Y,XLABEL,1,ISZXYL,0,IGRAPH,1)
        ELSE
          NDIVS = AMAX0(1,NDIVY)
          NSDIVS = AMAX0(1,NSDIVY)
          LOGS = LOGY
          IGRDS1 = IGRDY1
          IGRDS2 = IGRDY2
          AXMIN = YMIN
          AXMAX = YMAX
          X = XMIN - XX2 - XX3
          Y = (YMAX+YMIN)/2.
          CALL TEXT(X,Y,YLABEL,2,ISZXYL,0,IGRAPH,1)
        ENDIF
* Set size of tick marks along axis
        IF(IAX.EQ.0) AXLL = (YMAXA-YMINA)*3./FLOAT(MAXY-1)
        IF(IAX.EQ.1) AXLL = (XMAXA-XMINA)*3./FLOAT(MAXX-1)
        IF(LOHI.EQ.2) AXLL = 2.*AXLL
* Calc number of major/minor tick marks on axis
        JDIV = NSDIVS
        IF(LOGS.NE.0) THEN
*         Linear axis scale
          I1 = 0
          I2 = NDIVS
        ELSE
*         Log scale
          I1 = INT(AXMIN) - 1
          I2 = INT(AXMAX) + 1
        ENDIF
* Plot tick mark(s) on the axis
        DO 10 I=I1, I2
          DO 20 J = 1, JDIV
            AXL = AXLL
            IF(J.EQ.1) THEN
*             major division
              AXL = AXLL*1.5
              IGRID = IGRDS1
            ELSE
*             minor division
              IGRID = IGRDS2
            ENDIF
            IF(LOGS.NE.0) THEN
              XY = AXMIN+(AXMAX-AXMIN)/FLOAT(NDIVS)*(I+(J-1)/
     $             FLOAT(NSDIVS))
            ELSE
              XY = ALOG10(10.**FLOAT(I)*FLOAT(J))
            ENDIF
            IF(IAX.EQ.0) THEN
              CALL LINE(1,XY,YMAX,XY,YMAX-AXL,IGRAPH,NPP,0,0)
              CALL LINE(1,XY,YMIN,XY,YMIN+AXL,IGRAPH,NPP,0,0)
            ELSE
              CALL LINE(1,XMAX,XY,XMAX-AXL,XY,IGRAPH,NPP,0,0)
              CALL LINE(1,XMIN,XY,XMIN+AXL,XY,IGRAPH,NPP,0,0)
            ENDIF
            IF(IGRID.GT.0) THEN
              NPP = 0
              IF(IAX.EQ.0) CALL LINE(IGRID,XY,YMIN,XY,YMAX,IGRAPH,NPP,0,
     $                               0)
              IF(IAX.EQ.1) CALL LINE(IGRID,XMIN,XY,XMAX,XY,IGRAPH,NPP,0,
     $                               0)
            ENDIF
* Plot the axis value at this position
            IF(J.EQ.1.AND.IAX.EQ.0.AND.XY.GE.XMIN.AND.XY.LE.XMAX) THEN
              IF(LOGS.EQ.0) NDPX = MAX0(0,INT(0.1-XY))
              WRITE(FMT,"('(F20.',I1,')')") NDPX
              WRITE(AXVAL,FMT) XY
              IF(LOGS.EQ.0) WRITE(AXVAL,FMT) 10.**(XY)
              IB = 0
  40          IB = IB + 1
              IF(AXVAL(IB:IB).EQ.' '.OR.AXVAL(IB:IB).EQ.'0') GOTO 40
              IF(IB.GT.1) AXVAL = AXVAL(IB:LENG(AXVAL))
              IF(LENG(AXVAL).EQ.1.AND.AXVAL(1:1).EQ.'.') AXVAL = '0'
              IB = LENG(AXVAL)
              IF(AXVAL(IB:IB).EQ.'.') AXVAL = AXVAL(1:IB-1)
              Y = YMIN - YY3
              X = XY
              CALL TEXT(X,Y,AXVAL,1,ISZVAL,0,IGRAPH,1)
            ELSE IF(J.EQ.1.AND.IAX.EQ.1.AND.XY.GE.YMIN.AND.XY.LE.YMAX)
     $      THEN
              IF(LOGS.EQ.0) NDPY = MAX0(0,INT(0.1-XY))
              WRITE(FMT,"('(F20.',I1,')')") NDPY
              WRITE(AXVAL,FMT) XY
              IF(LOGS.EQ.0) WRITE(AXVAL,FMT) 10.**(XY)
              IB = 0
  50          IB = IB + 1
              IF(AXVAL(IB:IB).EQ.' '.OR.AXVAL(IB:IB).EQ.'0') GOTO 50
              IF(IB.GT.0) AXVAL = AXVAL(IB:LENG(AXVAL))
              IF(LENG(AXVAL).EQ.1.AND.AXVAL(1:1).EQ.'.') AXVAL = '0'
              IB = LENG(AXVAL)
              IF(AXVAL(IB:IB).EQ.'.') AXVAL = AXVAL(1:IB-1)
              X = XMIN - XX3
              Y = XY
              IF(IORIYV.EQ.0) CALL TEXT(X,Y,AXVAL,2,ISZVAL,0,IGRAPH,1)
              IF(IORIYV.NE.0) CALL TEXT(X,Y,AXVAL,1,ISZVAL,1,IGRAPH,1)
            ENDIF
  20      CONTINUE
  10    CONTINUE
  30  CONTINUE
      IF(LGDLNS.GT.0) CALL LGDBOX(LGDPOS,LGDLNS,LGDCHS,LGDSIZ,
     $                            LGDTYP,1,IGRAPH)
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE PLOTD(NPTS,LTYP,MARK,MSIZE,LEGEND,IORI,ITXSIZ,
     $  JUSTIF,INOUT,IGRAPH,STRNG,IFN,P,X,Y,ILEGND,LGDPOS,LGDSIZ)
*
* Plot the data for this curve, or this data set. This subroutine is not
* essential, but gives easy access to most of the plotting functions.
* Note that only some of the input parameters are used on each call,
* depending on the value of NPTS.
*
*   NPTS  = No of data points for current data set, or..
*           0 to plot a function curve,
*           -1 for plotting a text string,
*           -2 to clear a rectangular area of the graph, with optional border.
*   LTYP  = Line type for joining points.
*           0 = no line, 1 to 5 straight lines,
*           -5 to -1 for cubic spline fit (smooth curve) in line type 1..5.
*   MARK  = Symbol to be plotted at each point
*           (1 = dot, 2..8 = symbol, 32..126 = ASCII character)
*   MSIZE = Size of MARK (1..n) Try 3 for a start.
*   LEGEND= Text to describe each data set. Blank to suppress.
*   IORI  = Orientation for plotting a text string (if NPTS=-1).
*           1 = Normal (vertical), 2 = 90 deg anti-clockwise
*           3 = Upside down,       4 = 90 deg clockwise
*   ITXSIZ= Size of text to be plotted (1 to n).  2 = "normal".
*   JUSTIF= Text justification. -1=left, 0=centre, 1=right.
*   INOUT = Border restriction. 0=plot only inside axes, 1=anywhere on graph
*   IGRAPH= The graph memory array
*   STRNG = A string of text to be printed on the graph
*   IFN   = Function number (pre-compiled in FUNCT) to plot if NPTS=0.
*   P     = The array of parameters to be passed to function IFN.
*   X,Y   = The array of points to be plotted (when NPTS>0).
*           Also used for specifying coordinate values for plotting text,
*           box and function plots. ie. when NPTS =...
*                -2 : (x1,y1) and (x2,y2) are the box coordinates
*                -1 : place text at (x1,y1)
*                 0 : plot function from x1 to x2.
*   ILEGND= The legend number in the legend table. Updated automatically
*           whenever another entry is printed in the table.
*   LGDPOS= Position of the legend box.
*           1..4 = top left to top right, 5..8 = bottom left to bottom right.
*           Should be the same as the value used in subroutine AXES(..
*   LGDSIZ= Size of text for legend table (for drawing the legend box)
*
* First plot the legend for this data set, if required
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      CHARACTER*80 LEGEND, STRNG
      DIMENSION IGRAPH(MMAXX,MMAXY), X(*), Y(*), P(*)
* Write the legend into the legend table
      IF(NPTS.GE.0.AND.(LENG(LEGEND).GT.1.OR.LEGEND(1:1).NE.' ')) THEN
        ISIZE = MAX0(1,(LGDSIZ*LOHI+1)/2)
        DY = (YMAXA-YMINA)/FLOAT(MAXY)*FLOAT(ISIZE)
        ILEGND = ILEGND + 1
        YY = YMAX - 13*DY*FLOAT(ILEGND)
        IF(LGDPOS.LT.1.OR.LGDPOS.GT.8) LGDPOS = 1
        LGDP = LGDPOS
        IF(LGDPOS.GT.4) LGDP = LGDPOS - 4
        IF(LGDPOS.GT.4) YY = YMIN + 13*DY*FLOAT(ILEGND)
        DX = (XMAXA-XMINA)/FLOAT(MAXX)*FLOAT(LOHI)
        XX = XMIN + 18.*DX + FLOAT(LGDP-1)*(XMAX-XMIN)/4.
        CALL MARKPT(MARK,MSIZE,XX,YY,IGRAPH,1)
        XX1 = XX - 13.*DX
        XX2 = XX + 13.*DX
        NPP = 0
        IF(LTYP.NE.0) CALL LINE(IABS(LTYP),XX1,YY,XX2,YY,IGRAPH,NPP,0,1)
        XX = XX + 19.*DX
        CALL TEXT(XX,YY,LEGEND,1,LGDSIZ,-1,IGRAPH,1)
      ENDIF
* If NPTS=-2 then clear a rectangular graph area (x1,y1) to (x2,y2)
      IF(NPTS.EQ.-2) THEN
        IF(LOGX.EQ.0) THEN
          X(1) = ALOG10(X(1))
          X(2) = ALOG10(X(2))
        ENDIF
        IF(LOGY.EQ.0) THEN
          Y(1) = ALOG10(Y(1))
          Y(2) = ALOG10(Y(2))
        ENDIF
        CALL CLRBOX(X(1),Y(1),X(2),Y(2),LTYP,INOUT,IGRAPH)
* If NPTS=-1 then text is to be plotted at (x,y) coordinates.
      ELSEIF(NPTS.EQ.-1) THEN
        IF(LOGX.EQ.0) X(1) = ALOG10(X(1))
        IF(LOGY.EQ.0) Y(1) = ALOG10(Y(1))
        CALL TEXT(X(1),Y(1),STRNG,IORI,ITXSIZ,JUSTIF,IGRAPH,INOUT)
* If NPTS=0 then a function plot is required.
      ELSEIF(NPTS.EQ.0) THEN
        IF(LOGX.EQ.0) THEN
*-------- Log scale on x-axis, so take logs of X1, X2
          X(1) = ALOG10(X(1))
          X(2) = ALOG10(X(2))
        ENDIF
        CALL FUNCT(IFN,P,LTYP,X(1),X(2),IGRAPH,INOUT)
* Else plot the set of points
      ELSEIF(NPTS.GT.0) THEN
        IF(LOGX.EQ.0) THEN
*-------- Log scale for x-axis, so take logs of x values
          DO 20 J=1, NPTS
            X(J) = ALOG10(X(J))
  20      CONTINUE
        ENDIF
        IF(LOGY.EQ.0) THEN
*-------- Log scale for y-axis...
          DO 30 J=1, NPTS
            Y(J) = ALOG10(Y(J))
  30      CONTINUE
        ENDIF
*------ Plot the points
        CALL POINTS(MARK,MSIZE,LTYP,NPTS,X,Y,IGRAPH,INOUT)
      ENDIF
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE POINT(XX,YY,IGRAPH,INOUT)
*
* Plot a single point at the coordinate (XX,YY) on IGRAPH
*
* INOUT = 0 to plot the point only if it is inside the axes
*         1 to plot the point anywhere on the page
*
* Each integer in the array IGRAPH has 16 bits, and each bit which
* equals 1 is to be printed as a dot.  However, only 14 are used
* for various technical reasons. In other words the integer at
* IGRAPH(5,2) would contain 14 dot positions corresponding to column
* 5 and rows 17 to 32 of the graph to be printed on the dot matrix
* printer.  MASKHI converts for Hi-res 144 dpi vertical printing, in
* which the row is printed as two interleaved rows 1/144" apart.
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      DIMENSION IGRAPH(MMAXX,MMAXY), MASKHI(14)
      DATA MASKHI/1,8,2,9,3,10,4,11,5,12,6,13,7,14/
      X = XX
      Y = YY
      IF(INOUT.EQ.0) THEN
* clip at axis boundary
        X = AMAX1(XMIN,AMIN1(XMAX,X))
        Y = AMAX1(YMIN,AMIN1(YMAX,Y))
      ELSE
        X = AMAX1(XMINA,AMIN1(XMAXA,X))
        Y = AMAX1(YMINA,AMIN1(YMAXA,Y))
      ENDIF
      IF(IVH.EQ.1) THEN
*------ Vertical format (normal)
        ICOL = INT((X-XMINA)/(XMAXA-XMINA)*FLOAT(MAXX-1)+1.5)
        IY = INT((Y-YMINA)/(YMAXA-YMINA)*FLOAT(MAXY-1)+0.5)
      ELSE
*------ If horizontal format selected, then transform values
        ICOL = MAXY - INT((Y-YMINA)/(YMAXA-YMINA)*FLOAT(MAXY-1)+0.5)
        IY = INT((X-XMINA)/(XMAXA-XMINA)*FLOAT(MAXX-1)+0.5)
      ENDIF
      IROW = IY/14+1
      IBIT = IROW*14-IY
      IF(LOHI.EQ.2) IBIT = MASKHI(IBIT)
      CALL SETBIT(IGRAPH(ICOL,IROW),IBIT)
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE MARKPT(IPT,ISIZE,X,Y,IGRAPH,INOUT)
*
* Plot a mark (symbol) of size ISIZE, centred at the point (X,Y) on IGRAPH.
* The symbol shapes are stored as characters, and accessed from TEXT,
* so changes to TEXT may alter these definitions.
*
*  ISIZE = Size of the mark to be plotted (1..n)
*   IPT  =  1 : point
*           2 : open octagon
*           3 : filled  "
*           4 : open square
*           5 : filled  "
*           6 : open triangle
*           7 : filled  "
*           8 : cross
*           9 : plus
*          10 : star
*          11 : open diamond
*          12 : filled  "
*      13-31  : Might be used later. Blank for now.
*      32-126 : plot the corresponding ASCII character (Orientation = 1)
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      CHARACTER*80 CH, BLNK
      DIMENSION IGRAPH(MMAXX,MMAXY)
      DATA BLNK /' '/
      CH = BLNK
      CALL POINT(X,Y,IGRAPH,INOUT)
      IORI = 1
      IF(IPT.GT.1.AND.IPT.LE.126) THEN 
        CH(1:1) = CHAR(IPT)
        CALL TEXT(X,Y,CH,IORI,ISIZE,0,IGRAPH,INOUT)
      ENDIF
      IF(IPT.EQ.10) THEN
        CH(1:1) = CHAR(8)
        CALL TEXT(X,Y,CH,IORI,ISIZE,0,IGRAPH,INOUT)
      ENDIF
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE TEXT(X,Y,STRNG,IORI,ITXSIZ,JUSTIF,IGRAPH,INOUT)
*
* Plot a string of text. Also used to plot graph symbols, using character
* positions 1 to 31 in the ASCII sequence to define these symbols.
*
* (X,Y) = Graph coordinates of the centre of the justification character
*            (actually the centre of upper case characters)
* STRNG = The character string to be plotted
* IORI  = Orientation. 1 = Normal (vertical)
*                      2 = Rotated 90 deg anti-clockwise
*                      3 = Upside down
*                      4 = Rotated 90 deg clockwise
* ITXSIZ= Text size, 1 to n. (Steps of two for lo-res : 2, 4, 6, ...)
* JUSTIF= Justification. -1=Left justified, 0=centre, 1=right.
*                     Y---->   abcde          abcde    abcde
*                     X---->   ^                ^          ^
* IGRAPH= The graph memory array
* INOUT = Border restriction. 0=plot only inside axes, 1=anywhere on graph
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH, ICH, NB, NI, CHARS
      CHARACTER*80 STRNG
      CHARACTER CH
      DIMENSION IGRAPH(MMAXX,MMAXY), CHARS(3,126), MASK2(16),
     $ IDOTS(10,0:6)
      DATA MASK2/1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,
     $           32768/
      DATA CHARS/    0,    0,    0, 4208, 1090, 7185,28784, 1987, 7199,
     $            4344, 1090,15889,
     $           28920, 1987,15903, 8240, 1153, 3082,24624, 1921, 3086,
     $            8328,  257, 8714,16416, 1984, 2052,16416, 1984, 2052,
     $            8224, 1089, 2058,24608, 1985, 2062,    0,    0,    0,
     $               0,    0,    0,    0,    0,    0,    0,    0,    0,
     $               0,    0,    0,    0,    0,    0,    0,    0,    0,
     $               0,    0,    0,    0,    0,    0,    0,    0,    0,
     $               0,    0,    0,    0,    0,    0,    0,    0,    0,
     $               0,    0,    0,    0,    0,    0,    0,    0,    0,
     $               0,    0,    0,    0,    0,    0,    0,    0,    0,
     $               0,    0,    0,    0, 4000,    0,    0,    7,   56,
     $           30800,17031, 5183,20552, 4066, 9237, 8584,16646, 8969,
     $           18648, 2724, 1297,    0, 3077,    0,    0,18368,   32,
     $            2048, 1988,    0,24648, 2016, 4614,16416, 1984, 2052,
     $            6656,  112,    0,16416,  256, 2052, 6144,   96,    0,
     $            8200,  256, 8200, 2296,18468,15904, 2048,20450,    0,
     $            6276,18596,12580,18692,19236,17972, 8240,17537, 1087,
     $            2504,18981,20008,18552,18722, 1572,  256, 2532,24616,
     $           18648,18724,13860,18624, 2340,15397,22528,  865,    0,
     $           23040,  881,    0, 8224,17473,   32, 8272,  641, 5130,
     $            2048, 1092, 2058,  128, 2212,12324, 2296,19364,15658,
     $            8316, 2178, 7954,30980,18727,13860, 2296,18468, 8736,
     $           30980,18471,15904,18940,18724,16676,16892, 2308,16420,
     $            2296,18468,20004,16892,  256,32516, 2048,20452,   32,
     $            2056, 2080,16447,16892,  640,16657, 2556,16416,  256,
     $             508,  770,32528,  508,  514,32516, 2296,18468,15904,
     $           16892, 2308,12324, 2296,18596,16161,16892, 2436,12581,
     $           18632,18724, 9764,  256, 4068,16416, 2552,16416,32256,
     $            4592,   32,31745, 2552,16832,32256, 8588,  257,25354,
     $             384,  481,24584,10508,18724,24872,30720,18471,   32,
     $             128,  257,  514, 2048,18468,   63,   64, 2050, 4112,
     $             513, 4104,   64,    0, 3072,   40,10248,17057, 3850,
     $            2556,16929, 3592, 2104,16929, 4360, 2104,16929,32520,
     $           10296,17057, 3338,30752, 2307,   36, 2617,21033, 8072,
     $             508,  513, 3848, 2048,19425,    0,  513,25096,   47,
     $            8700,  384,  265, 2048,20452,    0,   60,  481, 3848,
     $             124,  513, 3848, 2104,16929, 3592, 2175,16929, 3592,
     $            2104,16929, 8136,16508,  512, 4104,10276,17057, 4618,
     $           28736,16935,  520, 2168,16416, 7936, 4208,   32, 7169,
     $            2168,16576, 7680,20548,  128, 4357, 2680,20520, 8064,
     $            6212,17057, 4364,16384,18112,16672,    0, 3808,    0,
     $            2308, 1732,    4,  128, 1028, 8200/
      ISIZE = MAX0(1,(ITXSIZ*LOHI+1)/2)
      IF(IORI.LT.1.OR.IORI.GT.4) IORI = 1
* Calc integer value of (X,Y), by analogy with POINT.
* First calc integer limits of plotting area
      IF(IVH.EQ.1) THEN
        IF(INOUT.EQ.0) THEN
          IMINX = INT((XMIN-XMINA)/(XMAXA-XMINA)*FLOAT(MAXX-1)+1.5)
          IMAXX = INT((XMAX-XMINA)/(XMAXA-XMINA)*FLOAT(MAXX-1)+1.5)
          IMINY = INT((YMIN-YMINA)/(YMAXA-YMINA)*FLOAT(MAXY-1)+0.5)
          IMAXY = INT((YMAX-YMINA)/(YMAXA-YMINA)*FLOAT(MAXY-1)+0.5)
        ELSE
          IMINX = 0
          IMAXX = MAXX
          IMINY = 0
          IMAXY = MAXY
        ENDIF
        IORIX = IORI
        IX = INT((X-XMINA)/(XMAXA-XMINA)*FLOAT(MAXX-1)+1.5)
        IY = INT((Y-YMINA)/(YMAXA-YMINA)*FLOAT(MAXY-1)+0.5)
      ELSE
        IF(INOUT.EQ.0) THEN
          IMINX = INT((XMIN-XMINA)/(XMAXA-XMINA)*FLOAT(MAXY-1)+0.5)
          IMAXX = INT((XMAX-XMINA)/(XMAXA-XMINA)*FLOAT(MAXY-1)+0.5)
          IMINY = INT((YMIN-YMINA)/(YMAXA-YMINA)*FLOAT(MAXX-1)+1.5)
          IMAXY = INT((YMAX-YMINA)/(YMAXA-YMINA)*FLOAT(MAXX-1)+1.5)
        ELSE
          IMINX = 0
          IMAXX = MAXY
          IMINY = 0
          IMAXY = MAXX
        ENDIF
        IORIX = IORI + 1
        IF(IORIX.GT.4) IORIX = 1
        IX = MAXY - INT((Y-YMINA)/(YMAXA-YMINA)*FLOAT(MAXY-1)+0.5)
        IY = INT((X-XMINA)/(XMAXA-XMINA)*FLOAT(MAXX-1)+0.5)
      ENDIF
* Calc offset for justification
      IF(JUSTIF.LT.0) IJUST = 0
      IF(JUSTIF.EQ.0) IJUST = (LENG(STRNG)-1)*3
      IF(JUSTIF.GT.0) IJUST = (LENG(STRNG)-1)*6
      DO 60 I=1, LENG(STRNG)
* Calc bottom left of this character
        GOTO(10,20,30,40) IORIX
  10    IC = IX + (I*6-9-IJUST)*ISIZE
        IR = IY - 6*ISIZE
        GOTO 50
  20    IC = IX + 6*ISIZE
        IR = IY + (I*6-9-IJUST)*ISIZE
        GOTO 50
  30    IC = IX + (9-I*6+IJUST)*ISIZE
        IR = IY + 6*ISIZE
        GOTO 50
  40    IC = IX - 6*ISIZE
        IR = IY + (9-I*6+IJUST)*ISIZE
* Plot the character
  50    IF(ISIZE.GT.1) THEN
          DO 200 J=1, 10
            DO 210 K=0, 6
              IDOTS(J,K) = 0
 210        CONTINUE
 200      CONTINUE
        ENDIF
        CH = STRNG(I:I)
        ICH = ICHAR(CH)
        DO 70 J=1, 5
          DO 80 K=1, 9
            NI = (J*9+K+5)/15
            NB = J*9+K+6-15*NI
            IF(MOD(CHARS(NI,ICH),MASK2(NB+1)).GE.MASK2(NB)) THEN
              CALL CHRDOT(IORIX,IC,IR,J*ISIZE,K*ISIZE,IMINX,IMAXX,
     $                    IMINY,IMAXY,IGRAPH)
              IDOTS(K,J) = 1
            ENDIF
  80      CONTINUE
  70    CONTINUE
        IF(ISIZE.GT.1) THEN
          DO 320 ISZ=1, ISIZE-1
            DO 300 J=1, 5
              DO 310 K=1, 9
                JJ = ISIZE*J
                KK = ISIZE*K
                IICH = 0
                IF(ICH.EQ.3.OR.ICH.EQ.5.OR.ICH.EQ.7.OR.ICH.EQ.12) IICH=1
                IF (IDOTS(K,J).EQ.1) THEN
                  IF(IDOTS(K,J+1).EQ.1) CALL CHRDOT(IORIX,IC,IR,JJ+ISZ,
     $             KK,IMINX,IMAXX,IMINY,IMAXY,IGRAPH)
                  IF(IDOTS(K+1,J).EQ.1) THEN
                    CALL CHRDOT(IORIX,IC,IR,JJ,KK+ISZ,IMINX,IMAXX,
     $                          IMINY,IMAXY,IGRAPH)
                  ELSE
                    IF(IDOTS(K,J+1).NE.1.AND.IDOTS(K+1,J+1).EQ.1) CALL
     $               CHRDOT(IORIX,IC,IR,JJ+ISZ,KK+ISZ,IMINX,IMAXX,
     $                      IMINY,IMAXY,IGRAPH)
                    IF(IDOTS(K,J-1).NE.1.AND.IDOTS(K+1,J-1).EQ.1) CALL
     $               CHRDOT(IORIX,IC,IR,JJ-ISZ,KK+ISZ,IMINX,IMAXX,
     $                      IMINY,IMAXY,IGRAPH)
                  ENDIF
                  IF(IICH.EQ.1) THEN
                    IF(IDOTS(K+1,J+1).EQ.1) CALL CHRDOT(IORIX,IC,IR,
     $           JJ+ISZ,KK+ISZ,IMINX,IMAXX,IMINY,IMAXY,IGRAPH)
                    IF(IDOTS(K+1,J-1).EQ.1) CALL CHRDOT(IORIX,IC,IR,
     $           JJ-ISZ,KK+ISZ,IMINX,IMAXX,IMINY,IMAXY,IGRAPH)
                  ENDIF
                ENDIF
 310          CONTINUE
 300        CONTINUE
 320      CONTINUE
        ENDIF
  60  CONTINUE
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE LINE(LTYP,X1,Y1,X2,Y2,IGRAPH,NPP,IFPBL,INOUT)
*
* Plot a line from (X1,Y1) to (X2,Y2) on IGRAPH.
* Calls POINT to plot the individual points.
*
* LTYP  = Line type : 1 =  continuous line
*                     2 =  .................
*                     3 =  . . . . . . . . .
*                     4 =  - - - - - - - - -
*                     5 =  -- . -- . -- . --
*
* IGRAPH= The graph memory array
* NPP   = the number of points plotted. Used for the line type patterns.
* IFPBL = 0 to print both end-points of the line,
*         else blank first end-point (for plotting chains of lines)
* INOUT = Border restriction. 0=plot only inside axes, 1=anywhere on graph
*
* DOT(ijk) = Array of patterns for dotted lines
*            i=0 or 1 for dot or no dot, j=LTYP, k=LOHI
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      DIMENSION IGRAPH(MMAXX,MMAXY), DOT(12,5,2)
      DATA DOT/1,1,1,1,1,1,1,1,1,1,1,1,
     $         1,0,1,0,1,0,1,0,1,0,1,0,
     $         1,0,0,1,0,0,1,0,0,1,0,0,
     $         1,1,1,1,0,0,1,1,1,1,0,0,
     $         1,1,1,0,1,0,1,1,1,0,1,0,
     $         1,1,1,1,1,1,1,1,1,1,1,1,
     $         1,0,0,1,0,0,1,0,0,1,0,0,
     $         1,0,0,0,0,0,1,0,0,0,0,0,
     $         1,1,1,1,1,1,1,1,0,0,0,0,
     $         1,1,1,1,1,1,1,0,0,1,0,0/
      IF(NPP.LT.1) NPP = 1
      IF(LTYP.LT.0.OR.LTYP.GT.5) LTYP = 1
* Calc min. no. of points for an unbroken line
      PNTS = ABS((X2-X1)/(XMAXA-XMINA)*FLOAT(MAXX))
      PNTS = AMAX1(1.,PNTS,ABS((Y2-Y1)/(YMAXA-YMINA)*FLOAT(MAXY)))
      IPNTS = INT(PNTS+0.5)
      DX = (X2-X1)/PNTS
      DY = (Y2-Y1)/PNTS
      ISTART = 1
      X = X1
      Y = Y1
      IF(IFPBL.EQ.0) THEN
        ISTART = 0
        X = X-DX
        Y = Y-DY
      ENDIF
      DO 10 I=ISTART, IPNTS
        IF(NPP.GT.12) NPP = 1
        X = X + DX
        Y = Y + DY
        IF(DOT(NPP,LTYP,LOHI).EQ.1) CALL POINT(X,Y,IGRAPH,INOUT)
        NPP = NPP + 1
  10  CONTINUE
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE SORT(NPTS,X,Y)
*
* Sort the data arrays X and Y into ascending order (of X values)
* NPTS = the number of values to be sorted.
*
      DIMENSION X(*),Y(*)
      M = NPTS
      L = M/2 + 1
  10  IF(L.GT.1) THEN
        L = L - 1
        XX = X(L)
        YY = Y(L)
      ELSE
        XX = X(M)
        YY = Y(M)
        X(M) = X(1)
        Y(M) = Y(1)
        M = M - 1
        IF(M.EQ.1) GOTO 30
      ENDIF
      I = L
      J = 2*L
      IF(J.LE.M) THEN
  20    IF(J.LT.M.AND.X(J).LT.X(J+1)) J = J+1
        IF(XX.LT.X(J)) THEN
          X(I) = X(J)
          Y(I) = Y(J)
          I = J
          J = 2*J
        ELSE
          J = M + 1
        ENDIF
        IF(J.LE.M) GO TO 20
      ENDIF
      X(I) = XX
      Y(I) = YY
      GO TO 10
  30  X(1) = XX
      Y(1) = YY
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE SPLINE(NPIN,NPOUT,X,Y,XX,YY)
*
* Fit a cubic spline to the NPIN (X,Y) pairs, and return NPOUT evenly
* spaced fitted (XX,YY) data pairs. Also returns sorted (X,Y).
*
      DIMENSION X(*), Y(*), XX(*), YY(*), S(500), C(500)
      CALL SORT(NPIN,X,Y)
      AI = X(2)-X(1)
      PI = (Y(2)-Y(1))/AI 
      S(1) = -1 
      C(1) = 0
      DI = -AI 
      CI = 0
      DO 10 I=2,NPIN-1
        A1 = X(I+1)-X(I) 
        Z = 2*(A1+AI)-DI
        P2 = (Y(I+1)-Y(I))/A1
        C(I) = (6*(P2-PI)-CI)/Z
        CI = C(I)*A1 
        PI = P2
        S(I) = A1/Z
        DI = S(I)*A1 
        AI = A1
  10  CONTINUE
      S(NPIN) = C(NPIN-1)/(1+S(NPIN-1))
      J = NPIN
      DO 20 I=1,NPIN-1
        J = J-1
        S(J) = C(J)-S(J)*S(J+1)
  20  CONTINUE
      IF(NPOUT.GT.500) NPOUT = 500
      DX = (X(NPIN)-X(1))/(NPOUT-1)
      XX(1) = X(1) 
      YY(1) = Y(1) 
      J = 1
      DO 30 I=1,NPIN-1 
        IF(XX(J)+DX.LT.X(I+1)) THEN
  40      J = J+1
          XX(J) = XX(J-1)+DX
          XV = XX(J)-X(I) 
          T = 2*S(I)+XV*(S(I)-S(I+1))/(X(I)-X(I+1))+S(I+1)
          YY(J)=Y(I)+XV*((Y(I)-Y(I+1))/(X(I)-X(I+1))+(XX(J)-X(I+1))*T/6)
          IF(XX(J)+DX.LT.X(I+1)) GOTO 40
        ENDIF
  30  CONTINUE
      XX(NPOUT) = X(NPIN)
      YY(NPOUT) = Y(NPIN)
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE POINTS(IPT,ISIZE,LTYP,NPTS,X,Y,IGRAPH,INOUT)
*
* Plot an array of points [X(i),Y(i)] on IGRAPH
*
* IPT   = Symbol to be plotted at each point
*         (1 = dot, 2..8 = symbol, 32..126 = ASCII character)
* ISIZE = Size of mark (1..n) Try 2 or 3 for a start.
* LTYP  = Line type for joining points.
*         0 = no line,
*         1 to 5 straight lines,
*         -5 to -1 for cubic spline fit (smooth curve).
* NPTS  = No of points to be plotted
* X,Y   = the array of points
* IGRAPH= The graph memory array
* INOUT = Border restriction. 0=plot only inside axes, 1=anywhere on graph
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      DIMENSION X(*), Y(*), IGRAPH(MMAXX,MMAXY), XX(500), YY(500)
      DO 10 I=1, NPTS
        CALL MARKPT(IPT,ISIZE,X(I),Y(I),IGRAPH,INOUT)
  10  CONTINUE
      NPP = 1
      IF(LTYP.GT.0) THEN
        DO 20 I=2,NPTS
          CALL LINE(LTYP,X(I-1),Y(I-1),X(I),Y(I),IGRAPH,NPP,1,INOUT)
  20    CONTINUE
      ELSEIF(LTYP.LT.0) THEN
*       Spline fit
*------ First sort into ascending x values, then calc and plot the spline
        NP = INT(ABS((X(NPTS)-X(1))/(XMAXA-XMINA)*(MAXX-1)/5.))
        CALL SPLINE(NPTS,NP,X,Y,XX,YY)
        DO 30 I=2, NP
        CALL LINE(-LTYP,XX(I-1),YY(I-1),XX(I),YY(I),IGRAPH,NPP,1,INOUT)
  30    CONTINUE
      ENDIF         
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE FUNCT(IFN,P,LTYP,X1,X2,IGRAPH,INOUT)
*
* Plot a continuous function from X1 to X2 on IGRAPH
*
* IFN   = Function number in subroutine GRAPHFNS
* P     = The array of parameters to be passed to function IFN.
* LTYP  = Line type for plotting the function.
*         0 = no line, 1 to 5 = solid, dotted etc..
* IGRAPH= The graph memory array
* INOUT = Border restriction. 0=plot only inside axes, 1=anywhere on graph
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      DIMENSION IGRAPH(MMAXX,MMAXY), P(*)
      IPNTS = INT(ABS((X2-X1)/(XMAXA-XMINA)*(MAXX-1)/4.))
      IF(IPNTS.LT.1) IPNTS = 1
      XB = X1
* Take logs as necessary if log scales specified
      IF(LOGX.GT.0) YB = FGRAPH(IFN,P,X1)
      IF(LOGX.EQ.0) YB = FGRAPH(IFN,P,10.**X1)
      IF(LOGY.EQ.0) YB = ALOG10(YB)
      DX = (X2-X1)/FLOAT(IPNTS)
      NPP = 1
      DO 10 I=1, IPNTS
        XA = XB
        YA = YB
        XB = X1 + I*DX
        IF(LOGX.GT.0) YB = FGRAPH(IFN,P,XB)
        IF(LOGX.EQ.0) YB = FGRAPH(IFN,P,10.**XB)
        IF(LOGY.EQ.0) YB = ALOG10(YB)
        CALL LINE(LTYP,XA,YA,XB,YB,IGRAPH,NPP,1,INOUT)
  10  CONTINUE
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE CLRGRF(IGRAPH)
*
* Set IGRAPH to zero - ie clear the graph
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      DIMENSION IGRAPH(MMAXX,MMAXY)
      DO 10 ICOL=1, MMAXX
        DO 20 IROW=1, MMAXY
          IGRAPH(ICOL,IROW) = 0
  20    CONTINUE
  10  CONTINUE
      RETURN
      END

*-----------------------------------------------------------------------

      SUBROUTINE CHRDOT(IORIX,IC,IR,J,K,IMINX,IMAXX,IMINY,IMAXY,
     $                  IGRAPH)
*
* For internal use by subroutine TEXT.
* Set this dot to on (for plotting characters).
* IC,IR = Col, row for bottom left of character
* J,K   = Offset in dots to current dot position
* IMINX..IMAXY = Integer limits of plotting area
*
      COMMON /GRF/XMINA,XMAXA,YMINA,YMAXA,XMIN,YMIN,XMAX,YMAX,MAXX,
     $            MAXY,MMAXX,MMAXY,LOHI,IVH,LOGX,LOGY,DPIV,DPIH
      INTEGER*2 IGRAPH
      DIMENSION IGRAPH(MMAXX,MMAXY), MASKHI(14)
      DATA MASKHI/1,8,2,9,3,10,4,11,5,12,6,13,7,14/
      GOTO(110,120,130,140) IORIX
 110  ICOL = IC + J
      IROW = IR + K
      GOTO 150
 120  ICOL = IC - K
      IROW = IR + J
      GOTO 150
 130  ICOL = IC - J
      IROW = IR - K
      GOTO 150
 140  ICOL = IC + K
      IROW = IR - J
 150  IF(ICOL.GT.IMINX.AND.ICOL.LT.IMAXX.AND.IROW.GT.IMINY.AND.IROW.
     $ LT.IMAXY) THEN
        IRR = IROW/14+1
        IBIT = IRR*14-IROW
        IF(LOHI.EQ.2) IBIT = MASKHI(IBIT)
        CALL SETBIT(IGRAPH(ICOL,IRR),IBIT)
      ENDIF
      RETURN
      END

*-----------------------------------------------------------------------

      FUNCTION LENG(STRING)
*
* Returns the length of STRING, excluding trailing blanks
*
      CHARACTER*80 STRING
      I = LEN(STRING)
      IF(I.LE.1) GOTO 20
  10  I = I - 1
      IF(STRING(I:I).EQ.' '.AND.I.GT.1) GOTO 10
  20  LENG = I
      RETURN
      END

*-----------------------------------------------------------------------
