!
! The grapheps Postscript functions
! See http://swiss.csail.mit.edu/~jaffer/Docupage/grapheps
!
module grapheps
  integer, parameter :: GSTRM = 22
contains
  subroutine pre_grapheps(outstr, xbound, ybound)
    integer, intent(in) :: outstr, xbound, ybound
    character (len=5) :: def0
#if IFORT || SUN
    character (len=24) :: fdate
#endif
    def0='0 def'
    write(outstr,'(a)')          '%!PS-Adobe-3.0 EPSF-3.0'
    write(outstr,'(a,2(1x,i5))') '%%BoundingBox: 0 0', xbound, ybound
    write(outstr,'(a)')          '%%Title: Sib-pair plot via grapheps'
    write(outstr,'(2a)')         '%%CreationDate: ', fdate()
    write(outstr,'(a)')          '%%EndComments'
    write(outstr,'(2i2,2(1x,i6))') 0, 0, xbound, ybound
    write(outstr,'(a)') '/plotdict 100 dict def'
    write(outstr,'(a)') 'plotdict begin'
    write(outstr,'(a)') '% Get dimensions the preamble left on the stack.'
    write(outstr,'(a)') '4 array astore /whole-page exch def'
    write(outstr,'(a)') '% Definitions so that internal assignments are bound before setting.'
    write(outstr,'(2a)') '/DATA ', def0
    write(outstr,'(2a)') '/DEN ', def0
    write(outstr,'(2a)') '/DIAG ', def0
    write(outstr,'(2a)') '/DIAG2 ', def0
    write(outstr,'(2a)') '/DLTA ', def0
    write(outstr,'(2a)') '/EXPSN ', def0
    write(outstr,'(2a)') '/GPROCS ', def0
    write(outstr,'(a)') '/GD 6 def'
    write(outstr,'(a)') '/GR 3 def'
    write(outstr,'(2a)') '/IDX ', def0
    write(outstr,'(2a)') '/ISIZ ', def0
    write(outstr,'(2a)') '/MAX ', def0
    write(outstr,'(2a)') '/MIN ', def0
    write(outstr,'(2a)') '/NUM ', def0
    write(outstr,'(2a)') '/PLOT-bmargin ', def0
    write(outstr,'(2a)') '/PLOT-lmargin ', def0
    write(outstr,'(2a)') '/PLOT-rmargin ', def0
    write(outstr,'(2a)') '/PLOT-tmargin ', def0
    write(outstr,'(2a)') '/PROC ', def0
    write(outstr,'(2a)') '/ROW ', def0
    write(outstr,'(2a)') '/TXT ', def0
    write(outstr,'(2a)') '/WPAGE ', def0
    write(outstr,'(2a)') '/X-COORD ', def0
    write(outstr,'(2a)') '/XDX ', def0
    write(outstr,'(2a)') '/XOFF ', def0
    write(outstr,'(2a)') '/XPARTS ', def0
    write(outstr,'(2a)') '/XRNG ', def0
    write(outstr,'(2a)') '/XSCL ', def0
    write(outstr,'(2a)') '/XSTEP ', def0
    write(outstr,'(2a)') '/XSTEPH ', def0
    write(outstr,'(2a)') '/XTSCL ', def0
    write(outstr,'(2a)') '/XWID ', def0
    write(outstr,'(2a)') '/Y-COORD ', def0
    write(outstr,'(2a)') '/YDX ', def0
    write(outstr,'(2a)') '/YHIT ', def0
    write(outstr,'(2a)') '/YOFF ', def0
    write(outstr,'(2a)') '/YPARTS ', def0
    write(outstr,'(2a)') '/YRNG ', def0
    write(outstr,'(2a)') '/YSCL ', def0
    write(outstr,'(2a)') '/YSTEP ', def0
    write(outstr,'(2a)') '/YSTEPH ', def0
    write(outstr,'(2a)') '/YTSCL ', def0
    write(outstr,'(2a)') '/graphrect ', def0
    write(outstr,'(2a)') '/plotrect ', def0
    write(outstr,'(a)') '% Here are the procedure-arrays for passing as the third argument to'
    write(outstr,'(a)') '% plot-column.  Plot-column moves to the first coordinate before'
    write(outstr,'(a)') '% calls to the first procedure.  Thus both line and scatter graphs are'
    write(outstr,'(a)') '% supported.  Many additional glyph types can be produced as'
    write(outstr,'(a)') '% combinations of these types.  This is best accomplished by calling'
    write(outstr,'(a)') '% plot-column with each component.'
    write(outstr,'(a)') '% GD and GR are the graphic-glyph diameter and radius.'
    write(outstr,'(a)') '% DIAG and DIAG2, used in /cross are diagonal and twice diagonal.'
    write(outstr,'(a)') '% gtrans maps x, y coordinates on the stack to 72dpi page coordinates.'
    write(outstr,'(a)') '% Render line connecting points'
    write(outstr,'(a)') '/line    [{} {lineto} {}] bind def'
    write(outstr,'(a)') '/mountain [{currentpoint 2 copy pop bottomedge moveto lineto}'
    write(outstr,'(a)') '	   {lineto}'
    write(outstr,'(a)') '	   {currentpoint pop bottomedge lineto closepath fill}] bind def'
    write(outstr,'(a)') '/cloud    [{currentpoint 2 copy pop topedge moveto lineto}'
    write(outstr,'(a)') '	   {lineto}'
    write(outstr,'(a)') '	   {currentpoint pop topedge lineto closepath fill}] bind def'
    write(outstr,'(a)') '% Render lines from x-axis to points'
    write(outstr,'(a)') '/impulse [{} {moveto XRNG 0 get 0 gtrans exch pop'
    write(outstr,'(a)') '	      currentpoint pop exch lineto} {}] bind def'
    write(outstr,'(a)') '/bargraph [{} {exch GR sub exch dup'
    write(outstr,'(a)') '	       XRNG 0 get 0 gtrans exch pop % y=0'
    write(outstr,'(a)') '	       exch sub GD exch rectstroke} {}] bind def'
    write(outstr,'(a)') '% Solid round dot.'
    write(outstr,'(a)') '/disc    [{GD setlinewidth 1 setlinecap}'
    write(outstr,'(a)') '			 {moveto 0 0 rlineto} {}] bind def'
    write(outstr,'(a)') '% Minimal point -- invisible if linewidth is 0.'
    write(outstr,'(a)') '/point   [{1 setlinecap} {moveto 0 0 rlineto} {}] bind def'
    write(outstr,'(a)') '% Square box.'
    write(outstr,'(a)') '/square  [{} {GR sub exch GR sub exch GD dup rectstroke} {}] bind def'
    write(outstr,'(a)') '% Square box at 45.o'
    write(outstr,'(a)') '/diamond [{}'
    write(outstr,'(a)') '	  {2 copy GR add moveto'
    write(outstr,'(a)') '	   GR neg GR neg rlineto GR GR neg rlineto'
    write(outstr,'(a)') '		   GR GR rlineto GR neg GR rlineto'
    write(outstr,'(a)') '	   closepath}'
    write(outstr,'(a)') '	  {}] bind def'
    write(outstr,'(a)') '% Plus Sign'
    write(outstr,'(a)') '/plus	 [{}'
    write(outstr,'(a)') '	  {       GR sub moveto  0 GD rlineto'
    write(outstr,'(a)') '	   GR neg GR neg rmoveto GD 0 rlineto}'
    write(outstr,'(a)') '	  {}] bind def'
    write(outstr,'(a)') '% X Sign'
    write(outstr,'(a)') '/cross   [{/DIAG GR .707 mul def /DIAG2 DIAG 2 mul def}'
    write(outstr,'(a)') '	  {exch DIAG sub exch DIAG add moveto DIAG2 dup neg rlineto'
    write(outstr,'(a)') '			  DIAG2 neg 0 rmoveto     DIAG2 dup rlineto}'
    write(outstr,'(a)') '	  {}] bind def'
    write(outstr,'(a)') '% Triangle pointing upward'
    write(outstr,'(a)') '/triup   [{}'
    write(outstr,'(a)') '	  {GR 1.12 mul add moveto GR neg GR -1.62 mul rlineto'
    write(outstr,'(a)') '	   GR 2 mul 0 rlineto GR neg GR  1.62 mul rlineto'
    write(outstr,'(a)') '	   closepath}'
    write(outstr,'(a)') '	  {}] bind def'
    write(outstr,'(a)') '% Triangle pointing downward'
    write(outstr,'(a)') '/tridown [{}'
    write(outstr,'(a)') '	  {GR 1.12 mul sub moveto GR neg GR  1.62 mul rlineto'
    write(outstr,'(a)') '	   GR 2 mul 0 rlineto GR neg GR -1.62 mul rlineto'
    write(outstr,'(a)') '	   closepath}'
    write(outstr,'(a)') '	  {}] bind def'
    write(outstr,'(a)') '/pentagon [{}'
    write(outstr,'(a)') '	   {gsave translate 0 GR moveto 4 {72 rotate 0 GR lineto} repeat'
    write(outstr,'(a)') '		  closepath stroke grestore}'
    write(outstr,'(a)') '	   {}] bind def'
    write(outstr,'(a)') '/circle   [{stroke} {GR 0 360 arc stroke} {}] bind def'
    write(outstr,'(a)') '% ( TITLE ) ( SUBTITLE )'
    write(outstr,'(a)') '/title-top'
    write(outstr,'(a)') '{ dup stringwidth pop -2 div plotrect 0 get plotrect 2 get 2 div add add'
    write(outstr,'(a)') '  plotrect 1 get plotrect 3 get add pointsize .4 mul add moveto show'
    write(outstr,'(a)') '  dup stringwidth pop -2 div plotrect 0 get plotrect 2 get 2 div add add'
    write(outstr,'(a)') '  plotrect 1 get plotrect 3 get add pointsize 1.4 mul add moveto show'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '% ( TITLE ) ( SUBTITLE )'
    write(outstr,'(a)') '/title-bottom'
    write(outstr,'(a)') '{ dup stringwidth pop -2 div plotrect 0 get plotrect 2 get 2 div add add'
    write(outstr,'(a)') '  plotrect 1 get pointsize -2 mul add moveto show'
    write(outstr,'(a)') '  dup stringwidth pop -2 div plotrect 0 get plotrect 2 get 2 div add add'
    write(outstr,'(a)') '  plotrect 1 get pointsize -1 mul add moveto show'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '% Plots column K against column J of given two-dimensional ARRAY.'
    write(outstr,'(a)') '% The arguments are:'
    write(outstr,'(a)') '%   [ ARRAY J K ] J and K are column-indexes into ARRAY'
    write(outstr,'(a)') '%   [ PREAMBLE RENDER POSTAMBLE ] Plotting procedures:'
    write(outstr,'(a)') '%       PREAMBLE  - Executed once before plotting row'
    write(outstr,'(a)') '%       RENDER    - Called with each pair of coordinates to plot'
    write(outstr,'(a)') '%       POSTAMBLE - Called once after plotting row (often does stroke)'
    write(outstr,'(a)') '/plot-column'
    write(outstr,'(a)') '{ /GPROCS exch def aload pop /YDX exch def /XDX exch def /DATA exch def'
    write(outstr,'(a)') '  /GD glyphsize def'
    write(outstr,'(a)') '  /GR GD .5 mul def'
    write(outstr,'(a)') '  gsave'
    write(outstr,'(a)') '    /ROW DATA 0 get def ROW XDX get ROW YDX get gtrans moveto'
    write(outstr,'(a)') '    GPROCS 0 get exec % preamble'
    write(outstr,'(a)') '    /PROC GPROCS 1 get def DATA {dup XDX get exch YDX get gtrans PROC} forall'
    write(outstr,'(a)') '    GPROCS 2 get exec stroke % postamble'
    write(outstr,'(a)') '  grestore'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '/partition-page'
    write(outstr,'(a)') '{ /YPARTS exch def /XPARTS exch def /WPAGE exch def'
    write(outstr,'(a)') '  /XWID WPAGE 2 get XPARTS div def /YHIT WPAGE 3 get YPARTS div def'
    write(outstr,'(a)') '  /Y-COORD WPAGE 1 get def'
    write(outstr,'(a)') '  YPARTS'
    write(outstr,'(a)') '  { /X-COORD WPAGE 0 get WPAGE 2 get add XWID sub def'
    write(outstr,'(a)') '    XPARTS {[X-COORD Y-COORD XWID YHIT]'
    write(outstr,'(a)') '	    /X-COORD X-COORD XWID sub def} repeat'
    write(outstr,'(a)') '    /Y-COORD Y-COORD YHIT add def'
    write(outstr,'(a)') '  } repeat'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '% The arguments are:'
    write(outstr,'(a)') '%   [ MIN-X MIN-Y DELTA-X DELTA-Y ] whole graph rectangle'
    write(outstr,'(a)') '%   [ MIN-COLJ MAX-COLJ ] Numerical range of plot data'
    write(outstr,'(a)') '%   [ MIN-COLK MAX-COLK ] Numerical range of plot data'
    write(outstr,'(a)') '% and the implicit current clippath'
    write(outstr,'(a)') '/setup-plot'
    write(outstr,'(a)') '{ /YRNG exch def /XRNG exch def /graphrect exch def'
    write(outstr,'(a)') '  /PLOT-bmargin pointsize 2.4 mul def'
    write(outstr,'(a)') '  /PLOT-tmargin pointsize 2.4 mul def'
    write(outstr,'(a)') '  /PLOT-lmargin lmargin-template stringwidth pop pointsize 1.2 mul add def'
    write(outstr,'(a)') '  /PLOT-rmargin rmargin-template stringwidth pop pointsize 1.2 mul add def'
    write(outstr,'(a)') '  /plotrect [ graphrect 0 get PLOT-lmargin add'
    write(outstr,'(a)') '	      graphrect 1 get PLOT-bmargin add'
    write(outstr,'(a)') '	      graphrect 2 get PLOT-lmargin sub PLOT-rmargin sub'
    write(outstr,'(a)') '	      graphrect 3 get PLOT-bmargin sub PLOT-tmargin sub ] def'
    write(outstr,'(a)') '  /XOFF XRNG 0 get def /YOFF YRNG 0 get def'
    write(outstr,'(a)') '  /XSCL plotrect 2 get XRNG aload pop exch sub div def'
    write(outstr,'(a)') '  /YSCL plotrect 3 get YRNG aload pop exch sub div def'
    write(outstr,'(a)') '  /XOFF XOFF plotrect 0 get XSCL div sub def'
    write(outstr,'(a)') '  /YOFF YOFF plotrect 1 get YSCL div sub def'
    write(outstr,'(a)') '  /YTSCL plotrect 3 get YRNG aload pop exch sub find-tick-scale def'
    write(outstr,'(a)') '  /YSTEP YTSCL 0 get 3 mod 0 eq {6} {8} ifelse 5 mul yuntrans def'
    write(outstr,'(a)') '  /XTSCL plotrect 2 get XRNG aload pop exch sub find-tick-scale def'
    write(outstr,'(a)') '  /XSTEP XTSCL 0 get 3 mod 0 eq {12} {10} ifelse 5 mul xuntrans def'
    write(outstr,'(a)') '  /YSTEPH YSTEP 2 div def'
    write(outstr,'(a)') '  /XSTEPH XSTEP 2 div def'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '% gtrans is the utility routine mapping data coordinates to view space.'
    write(outstr,'(a)') '% plot-column sets up XOFF, XSCL, and YSCL and uses it.'
    write(outstr,'(a)') '/gtrans {exch XOFF sub XSCL mul exch YOFF sub YSCL mul} bind def'
    write(outstr,'(a)') '%/guntrans {exch XSCL div XOFF add exch YSCL div YOFF add} bind def'
    write(outstr,'(a)') '% /ytrans {YTSCL aload pop div mul} bind def'
    write(outstr,'(a)') '% /xtrans {XTSCL aload pop div mul} bind def'
    write(outstr,'(a)') '/yuntrans {YTSCL aload pop exch div mul} bind def'
    write(outstr,'(a)') '/xuntrans {XTSCL aload pop exch div mul} bind def'
    write(outstr,'(a)') '/zero-in-range? {dup 0 get 0 le exch 1 get 0 ge and} bind def'
    write(outstr,'(a)') '/y-axis'
    write(outstr,'(a)') '{ XRNG zero-in-range?'
    write(outstr,'(a)') '    { 0 YRNG 0 get gtrans moveto 0 YRNG 1 get gtrans lineto stroke} if'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '/x-axis'
    write(outstr,'(a)') '{ YRNG zero-in-range?'
    write(outstr,'(a)') '    {XRNG 0 get 0 gtrans moveto XRNG 1 get 0 gtrans lineto stroke} if'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '% Find data range in column K of two-dimensional ARRAY.'
    write(outstr,'(a)') '%   ARRAY'
    write(outstr,'(a)') '%   K  is the column-index into ARRAY'
    write(outstr,'(a)') '/column-range'
    write(outstr,'(a)') '{ /IDX exch def dup /MIN exch 0 get IDX get def /MAX MIN def'
    write(outstr,'(a)') '  {IDX get dup dup MIN lt {/MIN exch def} {pop} ifelse'
    write(outstr,'(a)') '	       dup MAX gt {/MAX exch def} {pop} ifelse} forall'
    write(outstr,'(a)') '  [MIN MAX]'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '/min {2 copy lt {pop} {exch pop} ifelse} bind def'
    write(outstr,'(a)') '/max {2 copy gt {pop} {exch pop} ifelse} bind def'
    write(outstr,'(a)') '/combine-ranges'
    write(outstr,'(a)') '{ aload pop 3 2 roll aload pop exch 4 3 roll min 3 1 roll max 2 array astore}'
    write(outstr,'(a)') 'bind def'
    write(outstr,'(a)') '/pad-range'
    write(outstr,'(a)') '{ exch aload pop /MAX exch def /MIN exch def'
    write(outstr,'(a)') '  /EXPSN exch 100 div MAX MIN sub mul def'
    write(outstr,'(a)') '  [ MIN EXPSN sub  MAX EXPSN add ]'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '/snap-range'
    write(outstr,'(a)') '{dup aload pop exch sub 1 exch find-tick-scale aload pop'
    write(outstr,'(a)') '     /DEN exch def /NUM exch def 1 NUM div DEN mul /DLTA exch def'
    write(outstr,'(a)') '     aload pop /MAX exch def /MIN exch def'
    write(outstr,'(a)') '     [ DLTA MAX MIN sub sub 2 div dup MIN exch sub exch MAX add ]'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '% Given the width (or height) and the data-span, returns an array of'
    write(outstr,'(a)') '% numerator and denominator (NUM DEN)'
    write(outstr,'(a)') '%'
    write(outstr,'(a)') '% NUM will be 1, 2, 3, 4, 5, 6, or 8 times a power of ten.'
    write(outstr,'(a)') '% DEN will be a power of ten.'
    write(outstr,'(a)') '%'
    write(outstr,'(a)') '% NUM   ISIZ'
    write(outstr,'(a)') '% === < ===='
    write(outstr,'(a)') '% DEN   DLTA'
    write(outstr,'(a)') '/find-tick-scale'
    write(outstr,'(a)') '{/DLTA exch def /ISIZ exch def'
    write(outstr,'(a)') ' /DEN 1 def'
    write(outstr,'(a)') ' {DLTA ISIZ le {exit} if /DEN DEN 10 mul def /ISIZ ISIZ 10 mul def} loop'
    write(outstr,'(a)') ' /NUM 1 def'
    write(outstr,'(a)') ' {DLTA 10 mul ISIZ ge {exit} if /NUM NUM 10 mul def /DLTA DLTA 10 mul def} loop'
    write(outstr,'(a)') ' [[8 6 5 4 3 2 1] {/MAX exch def MAX DLTA mul ISIZ le {MAX exit} if} forall'
    write(outstr,'(a)') '  NUM mul DEN]'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '/rule-vertical'
    write(outstr,'(a)') '{ /XWID exch def'
    write(outstr,'(a)') '  /TXT exch def'
    write(outstr,'(a)') '  /X-COORD exch def'
    write(outstr,'(a)') '  X-COORD type [] type eq {/X-COORD X-COORD 0 get def} if'
    write(outstr,'(a)') '  gsave'
    write(outstr,'(a)') '    X-COORD plotrect 1 get plotrect 3 get 2 div add translate'
    write(outstr,'(a)') '    TXT stringwidth pop -2 div'
    write(outstr,'(a)') '    XWID 0 gt { 90 rotate PLOT-lmargin} {-90 rotate PLOT-rmargin} ifelse'
    write(outstr,'(a)') '    pointsize 1.2 mul sub moveto TXT show'
    write(outstr,'(a)') '  grestore'
    write(outstr,'(a)') '  YRNG 0 get YSTEP div ceiling YSTEP mul  YSTEP  YRNG 1 get'
    write(outstr,'(a)') '  { /YDX exch def 0 YDX gtrans /Y-COORD exch def pop'
    write(outstr,'(a)') '    X-COORD Y-COORD moveto XWID 0 rlineto stroke'
    write(outstr,'(a)') '    /TXT YDX 20 string cvs def'
    write(outstr,'(a)') '    X-COORD'
    write(outstr,'(a)') '    XWID 0 gt {TXT stringwidth pop sub ( ) stringwidth pop sub'
    write(outstr,'(a)') '	       Y-COORD pointsize .3 mul sub moveto}'
    write(outstr,'(a)') '	      {Y-COORD pointsize .3 mul sub moveto ( ) show} ifelse'
    write(outstr,'(a)') '    TXT show} for'
    write(outstr,'(a)') '  YRNG 0 get YSTEPH div ceiling YSTEPH mul  YSTEPH  YRNG 1 get'
    write(outstr,'(a)') '  { /YDX exch def 0 YDX gtrans /Y-COORD exch def pop'
    write(outstr,'(a)') '    X-COORD Y-COORD moveto XWID 2 div 0 rlineto stroke} for'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '/rule-horizontal'
    write(outstr,'(a)') '{ /YHIT exch def'
    write(outstr,'(a)') '  /TXT exch def'
    write(outstr,'(a)') '  /Y-COORD exch def'
    write(outstr,'(a)') '  Y-COORD type [] type eq {/Y-COORD Y-COORD 1 get def} if'
    write(outstr,'(a)') '  plotrect 0 get plotrect 2 get 2 div add TXT stringwidth pop -2 div add'
    write(outstr,'(a)') '  Y-COORD'
    write(outstr,'(a)') '  YHIT 0 gt {pointsize -2 mul} {pointsize 1.4 mul} ifelse add moveto TXT show'
    write(outstr,'(a)') '  XRNG 0 get XSTEP div ceiling XSTEP mul  XSTEP  XRNG 1 get'
    write(outstr,'(a)') '  { dup 0 gtrans pop /X-COORD exch def'
    write(outstr,'(a)') '    X-COORD Y-COORD moveto 0 YHIT rlineto stroke'
    write(outstr,'(a)') '    /TXT exch 10 string cvs def'
    write(outstr,'(a)') '    X-COORD TXT stringwidth pop 2.0 div sub'
    write(outstr,'(a)') '    Y-COORD YHIT 0 gt {pointsize sub} {pointsize .3 mul add} ifelse'
    write(outstr,'(a)') '    moveto TXT show'
    write(outstr,'(a)') '  } for'
    write(outstr,'(a)') '  XRNG 0 get XSTEPH div ceiling XSTEPH mul  XSTEPH  XRNG 1 get'
    write(outstr,'(a)') '  { 0 gtrans pop Y-COORD moveto 0 YHIT 2 div rlineto stroke} for'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '/grid-verticals'
    write(outstr,'(a)') '{ XRNG 0 get XSTEPH div ceiling XSTEPH mul  XSTEPH  XRNG 1 get'
    write(outstr,'(a)') '  { 0 gtrans pop /X-COORD exch def'
    write(outstr,'(a)') '    X-COORD plotrect 1 get moveto 0 plotrect 3 get rlineto} for'
    write(outstr,'(a)') '  stroke'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '/grid-horizontals'
    write(outstr,'(a)') '{ YRNG 0 get YSTEPH div ceiling YSTEPH mul  YSTEPH  YRNG 1 get'
    write(outstr,'(a)') '  { 0 exch gtrans /Y-COORD exch def pop'
    write(outstr,'(a)') '    plotrect 0 get Y-COORD moveto plotrect 2 get 0 rlineto} for'
    write(outstr,'(a)') '  stroke'
    write(outstr,'(a)') '} bind def'
    write(outstr,'(a)') '/leftedge {plotrect 0 get} bind def'
    write(outstr,'(a)') '/rightedge {plotrect dup 0 get exch 2 get add} bind def'
    write(outstr,'(a)') '/topedge {plotrect dup 1 get exch 3 get add} bind def'
    write(outstr,'(a)') '/bottomedge {plotrect 1 get} bind def'
    write(outstr,'(a)') '/outline-rect {aload pop rectstroke} bind def'
    write(outstr,'(a)') '/fill-rect {aload pop rectfill} bind def'
    write(outstr,'(a)') '/clip-to-rect {aload pop rectclip} bind def'
    write(outstr,'(a)') '/gstack [] def'
    write(outstr,'(a)') '/gpush {gsave /gstack [ gstack pointsize glyphsize ] def} bind def'
    write(outstr,'(a)') '/gpop {/gstack gstack aload pop /glyphsize exch def /pointsize exch def def grestore} bind def'
    write(outstr,'(a)') '% Default parameters'
    write(outstr,'(a)') '% The legend-templates are strings used to reserve horizontal space'
    write(outstr,'(a)') '/lmargin-template (-.0123456789) def'
    write(outstr,'(a)') '/rmargin-template (-.0123456789) def'
    write(outstr,'(a)') '% glyphsize is the graphic-glyph size; GR, graphic radius, is'
    write(outstr,'(a)') '% glyphsize/2.  Line width, set by "setlinewidth", must be much less'
    write(outstr,'(a)') '% than glyphsize for readable glyphs.'
    write(outstr,'(a)') '/glyphsize 6 def'
    write(outstr,'(a)') '% pointsize is the height of text characters in "points", 1/72 inch; 0.353.mm'
    write(outstr,'(a)') '/pointsize 12 def'
    write(outstr,'(a)') '% Set default font'
    write(outstr,'(a)') '/Helvetica pointsize selectfont'
    write(outstr,'(a)') 'gsave'
  end subroutine pre_grapheps
!
! Draw simple plot
!
! typ 1=scatterplot 2=jittered dotplot 3=mountain 4=bargraph
!     10=square Q-Q plot
!
  subroutine xy_grapheps(outstr, nvals, xvals, yvals,  &
                         xlab, ylab, title, gstyle, gratio, typ)
    integer, intent(in) :: outstr
    integer, intent(in) :: nvals
    double precision, dimension(nvals), intent(inout) :: xvals, yvals
    character (len=*), intent(in) :: xlab, ylab, title
    character (len=*), intent(in) :: gstyle
    double precision, intent(in) :: gratio
    integer, intent(in) :: typ
!
    integer :: i, linewidth
    character (len=len(gstyle)) :: style
    double precision :: rang, xmin, xmax, ymin, ymax
    integer :: xwid, ywid
! functions
    real :: random
    linewidth=1
    style=gstyle
    ywid=500
    xwid=int(gratio*dfloat(ywid))
    if (style==' ') style='point'
    xmax=-1.0d99
    xmin=+1.0d99
    ymax=-1.0d99
    ymin=+1.0d99
    do i=1, nvals
      xmin=min(xmin, xvals(i))
      xmax=max(xmax, xvals(i))
      ymin=min(ymin, yvals(i))
      ymax=max(ymax, yvals(i))
    end do
    if (typ == 10) then
      xmax=max(ymax, xmax)
      ymax=xmax
      xmin=min(ymin, xmin)
      ymin=xmin
      linewidth=2
      if (nvals < 50000) linewidth=5
    end if
    rang=xmax-xmin
! jitter
    if (typ==2) then
      do i=1, nvals
        xvals(i)=xvals(i)+dble(0.1*(random()-0.5))
      end do
      xmax=xmax+1.0d0
      xmin=xmin-1.0d0
    else if (typ == 4) then
      xmax=xmax+0.1d0*rang
      xmin=xmin-0.1d0*rang
    else
      xmax=xmax+0.03d0*rang
      xmin=xmin-0.03d0*rang
    end if
    rang=ymax-ymin
    ymax=ymax+0.03d0*rang
    ymin=ymin-0.03d0*rang
    if (typ == 3) then
      ymin=0.0d0
      ymax=max(4.0d0, ymax)
    else if (typ == 4) then
      ymin=0.0d0
    end if
    call pre_grapheps(outstr, xwid, ywid)
    if (typ == 2) then
      write(outstr,'(a)')  '/vlines', '['
      write(outstr,*)  '[', 1, 2, ymin, ']'
      write(outstr,*)  '[', 1, 2, ymax, ']'
      write(outstr,'(a)')  '] def'
    else if (typ == 4) then 
      write(outstr,*)  '/glyphsize', 320/nvals, ' def'
    else if (typ == 10) then
      write(outstr,'(a)')  '/identity', '['
      write(outstr,*)  '[', xmin, ymin, ']'
      write(outstr,*)  '[', xmax, ymax, ']'
      write(outstr,'(a)')  '] def'
    end if
    write(outstr,'(a)')  '/Data', '['
    do i=1, nvals
      write(outstr,*) '[', xvals(i), yvals(i), ']'
    end do
    write(outstr,'(a)')  '] def'
    write(outstr,*)      'whole-page [ ', xmin, xmax, ' ] [', ymin, ymax, '] setup-plot'
    if (typ /= 4) then
      write(outstr,*)  &
        '(', trim(title), ') (', trim(ylab) ,' versus ', trim(xlab),') title-top'
    else
      write(outstr,*)  &
        '(Frequency bargraph for "', trim(xlab),'") (', trim(title), ') title-top'
    end if
    write(outstr,'(a)')  'plotrect outline-rect'
    write(outstr,*)      'leftedge (', ylab, ') 5 rule-vertical'
    write(outstr,'(a)')  'gpush'
    write(outstr,'(a)')  'plotrect clip-to-rect'
    if (trim(style) == 'mountain') then
      write(outstr,'(a)') '.9 .9 .9 setrgbcolor'
    end if
    write(outstr,'(i0,a)') linewidth, ' setlinewidth'
    write(outstr,'(3a)') '[ Data 0 1 ] ', trim(style), ' plot-column'
    write(outstr,'(a)')  '1 setlinewidth'
    write(outstr,'(a)')  'gpop'
    write(outstr,*)      'bottomedge (', trim(xlab), ') 5 rule-horizontal'
! Vertical guidelines at "n" and "y"
    if (typ == 2) then
      write(outstr,'(a)')  'gpush', '[ 5 2 ] 0 setdash', '.9 .9 .9 setrgbcolor'
      write(outstr,'(a)')  '[ vlines 0 2 ] line plot-column'
      write(outstr,'(a)')  '[ vlines 1 2 ] line plot-column', 'gpop'
! Line of identity
    else if (typ == 10) then
      write(outstr,'(a)')  'gpush', '[ 5 2 ] 0 setdash', '.9 .9 .9 setrgbcolor'
      write(outstr,'(a)')  '[ identity 0 1 ] line plot-column', 'gpop'
    end if
! 
    write(outstr,'(a)')  'grestore', 'end', 'showpage'
  end subroutine xy_grapheps
!
! Draw scatterplot with different glyphs for each category
!
  subroutine scatter_grapheps(outstr, nvals, symbols, xvals, yvals,  &
                              slab, xlab, ylab, title, gratio)
    integer, intent(in) :: outstr
    integer, intent(in) :: nvals
    integer, dimension(nvals), intent(inout) :: symbols
    double precision, dimension(nvals), intent(inout) :: xvals, yvals
    character (len=*), intent(in) :: slab, xlab, ylab, title
    double precision, intent(in) :: gratio
!
    integer :: i
    double precision :: rang, xmin, xmax, ymin, ymax
    integer :: xwid, ywid
    integer, dimension(10) :: slevels
    character (len=8), dimension(10) :: style = (/  &
      'circle  ', 'disc    ', 'square  ', 'triup   ', 'plus    ',  &
      'diamond ', 'cross   ', 'pentagon', 'tridown ', 'point   ' /)

    ywid=500
    xwid=int(gratio*dfloat(ywid))
    xmax=-1.0d99
    xmin=+1.0d99
    ymax=-1.0d99
    ymin=+1.0d99
    do i=1, nvals
      xmin=min(xmin, xvals(i))
      xmax=max(xmax, xvals(i))
      ymin=min(ymin, yvals(i))
      ymax=max(ymax, yvals(i))
    end do
    rang=xmax-xmin
    xmax=xmax+0.03d0*rang
    xmin=xmin-0.03d0*rang
    rang=ymax-ymin
    ymax=ymax+0.03d0*rang
    ymin=ymin-0.03d0*rang
    call pre_grapheps(outstr, xwid, ywid)
    do j=1, 10
      slevels(j)=0
      write(outstr,'(a,i0/a)')  '/Data', j, '['
      do i=1, nvals
      if (symbols(i) == j) then
        slevels(j)=slevels(j)+1
        write(outstr,*) '[', xvals(i), yvals(i), ']'
      end if
      end do
      write(outstr,'(a)')  '] def'
    end do
    write(outstr,*)  &
      'whole-page [ ', xmin, xmax, ' ] [', ymin, ymax, '] setup-plot'
    write(outstr,*)  &
      '(', trim(title), ') ("', trim(ylab) ,'" versus "', trim(xlab),  &
      '" \(symbol type represent "', trim(slab), '"\) ) title-top'
    write(outstr,'(a)')  'plotrect outline-rect'
    write(outstr,*)      'leftedge (', ylab, ') 5 rule-vertical'
    write(outstr,'(a)')  'gpush'
    write(outstr,'(a)')  'plotrect clip-to-rect'
    do j=1, 10
    if (slevels(j) > 0) then
      write(outstr,'(a,i0,3a)')  &
        '[ Data', j, ' 0 1 ] ', trim(style(j)), ' plot-column'
    end if
    end do
    write(outstr,'(a)')  'gpop'
    write(outstr,*)      'bottomedge (', trim(xlab), ') 5 rule-horizontal'
    write(outstr,'(a)')  'grestore', 'end', 'showpage'
  end subroutine scatter_grapheps
end module grapheps
!
! Example call
!  
program doplot
  use grapheps 
  implicit none
  integer, parameter :: GSTRM=8
  character (len=256) :: fil = 'example.ps'
  character (len=10) :: style
  character (len=80) :: xlab, ylab, title
  integer :: i, ios, nobs
  integer, dimension(:), allocatable :: symbols
  double precision, dimension(:), allocatable :: xvals, yvals

  style='circle'
  title='Example plot'
  xlab='X variable'
  ylab='Y variable'
  read(*,*) nobs
  if (nobs > 0) then
    open(GSTRM, file=trim(fil), iostat=ios)
    if (ios /= 0) then
      write(outstr, '(a)') 'ERROR:  Unable to open "', trim(fil), '".'
      stop
    end if
    allocate(xvals(nobs), yvals(nobs))
    do i=1, nobs
      read(*,*) xvals(i), yvals(i)
    end do
    call xy_grapheps(GSTRM, nobs, xvals, yvals, xlab, ylab, title, style, 1.2d0, 1)
    close(GSTRM, status='keep')
  end if
end program doplot