!
! 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