program xyplo(xyin, xyout, xyplop, xyplom, warnings, output);
(* xyplo: x, y data plotter (pronounced: "zyplo")

  Thomas D. Schneider, Ph.D.
  toms@alum.mit.edu
  https://alum.mit.edu/www/toms  

  modules: prgmods, matmods, delman, dops *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 9.51; (* of xyplo.p 2024 Jan 17
2024 Jan 17, 9.51: y axis 90 degree rotation into two 45 for LaTeX/skim?
                   no ... just needed more space on left side!!!
2016 Oct 05, 9.50: debug denplo output, Kevin's 'testdenplo'
2016 Oct 03, 9.49: clean up
2016 Oct 03, 9.48: debug ylabel - solved!
2016 Oct 03, 9.47: debug ylabel
2016 Oct 03, 9.46: move xlabel and ylabel with %x and %y
2016 Sep 29, 9.45: ylabel on side; use '^' to get old method
2016 Sep 29, 9.44: ylabel on side ...  I modified makeylabel it worked!!
2014 Oct 23, 9.43: bounding box on white background
2014 Oct 23, 9.42: backup
2014 Oct 23, 9.41: write a white background rectangle over the bounding box
                   This makes conversion to jpg give a white background.
2014 Mar 05, 9.40: y intercept error code from Greg Alvord - incomplete!
2014 Feb 06, 9.39: clean up
2014 Feb 06, 9.38: M mark can be defined in xyplom with '\name{'
2011 Oct 20, 9.37: clean up
2011 Oct 20, 9.36: no y range variation?
2011 Oct 20, 9.35: backup
2011 Aug 01, 9.34: don't put %% in xyout from xyplom
2011 Aug 01, 9.33: If x axis tic lengths changed with Xticscale or ticscale,
                   then move the label down out of the way! In makexlabel.
2011 Jul 31, 9.32: document M
2011 Jul 31, 9.31: turn debug off
2011 Jul 31, 9.30: debug what happens when using M symbol-to-plot
                   but there is nothing defined in xyplom.

                   When there are many symbols and M is used, the
                   symbols are named 0-9,a-z,A-Z.

                   The first ten cases put 0 to 9 on the stack.  These
                   prevented any graphic and so there was no symbol
                   written.  The numbers stayed on the stack.

                   When 'a' went on the stack, it executed the /a and
                   that was ok.

                   When 'b' went on the stack, it failed since /b was
                   not defined.  I located the problem with the new
                   pslocate script.

                   The lack of use of the numbers explains the large
                   stack at the time of the crash

                   We must look for the flagstring in xyplom!

2011 Jul 31, 9.29: warnedaboutusermarks
2011 Mar 03, 9.28: gettoken must parse on tabs too!
2010 Dec 07, 9.27: rgb switch!
2010 Dec 02, 9.26: document tic controls
2010 Dec 02, 9.25: more tic controls - X, Y independent!
2010 Nov 02, 9.23: y axis peak rounding to avoid missing highest number?
2010 Nov 02, 9.22: obtain new graphstring from dops.p
2010 Nov 02, 9.21: dang parenthesis 2008 Aug 5 was not right!
2010 Oct 26, 9.20: if setfont in xyplom, don't set font!
2010 Oct 21, 9.19: balance gsave and grestore.  Sheesh.
2010 Apr 09, 9.18: Plotdata issue
2010 Apr 09, 9.17: c and C are given as diameters in the xyin file
2008 Aug 05, 9.16: handle parenthesis in titles
2007 Aug 30, 9.15: upgrade axis - subtics to the end (finished)
2007 Aug 30, 9.14: upgrade axis - subtics to the end
2007 May 24, 9.13: count data points so warning is more informative
2007 May 05, 9.12: stop if xyin is empty
2006 Dec 06, 9.11: tag regression line in output
2006 Dec 06, 9.10: add r^2 to output
2006 Sep 07, 9.09: Fisher's z' transformation
2005 Aug 14, 9.08: remove duplicate warnings
2005 Jul 26, 9.07: upgrade from prgmod to get tabs ignored
2005 Jul 26, 9.06: user linetype can be any of l.-in
2004 Aug 25, 9.05: bug: in edgecontrol ury was a function of x
2004 Aug  2, 9.04: Accept '#' as the start of a comment line.
2004 Jan 19, 9.03: change color adjustment to be 0.84*hue+0.16
2003 Dec 21, 9.02: setlinewidth bug
2003 Aug 28, 9.01: better error message for missing parameterversion
2003 Aug 28, 9.00: substitute in corrected stringreal from dops.p
2003 Aug 28, 8.99: about to substitute in corrected stringreal from dops.p
2003 Aug 28, 8.98: bug: at 0.0, roundoff makes -0.00!
2003 Aug 28, 8.97: cleanup
2003 Aug 28, 8.96: bug1992: - captured and fixed
2003 Aug 28, 8.95: bug1992: - stringinteger cannot handle -0 of course!
2003 Aug 28, 8.94: bug1992: - sign missing on y axis
2003 Aug 28, 8.93: read and write edge and parameterversion params correctly
2003 Aug 27, 8.92: bug: userlines set nil in readparams
2003 Aug 22, 8.91: cleanup message
2003 Aug 22, 8.90: space in line section bombs.
2003 Aug 22, 8.89: upgrading xyplop to create edge control parameters
2003 Aug 22, 8.88: bug in boundary box fixed
2003 Aug 22, 8.87: fix bounding box values to be integer
2003 Aug 22, 8.86: clean up
2003 Aug 22, 8.84: user control of bounding box edges
2003 Aug 22, 8.83: clean up
2003 Aug 22, 8.82: graph controls bounding box
2003 Jul 31, 8.81: tabs accepted as blanks
2003 Mar 11, 8.80: xyplop.std in see also
2002 Nov  7, 8.79: fill bottom of skybox
2002 Nov  7, 8.78: skybox (s symbol) allows histogram
2002 May  5, 8.77: pure comment xyin crashes program; fix skipcopy
2002 May  5, 8.76: gpc ready; warnings file used in main now.
2001 Jun  1, 8.75: define postscript bounding box
2000 Jun 21, 8.74: no warning file unless there are warnings
2000 May 26, 8.73: introduce maxinterval to prevent too many intervals
2000 Mar 30, 8.72: allow removal of axis
2000 Jan 28, 8.71: L symbol for drawing thick lines
1999 Dec 4, 8.69: put warnings into warnings file to keep out of way.
1999 Dec 4, 8.68: Warning for neg log values gives both x and y
1999 Dec 4, 8.67: blank lines and lines that begin with "*" are ignored in
                  user line section.  Negative values ignored for log scale
                  rather than halt.
1999 Dec 3, 8.66: on a log plot, values <= zero are dropped rather than
                  messing up the graph.
1999 July 30, 8.65: C, B colored lines
1999 June 1, 8.64: correct documentation: colors range from 0 to 0.85.
1996 August 5 xyplo reads user defined symbols in xyplom
1996 August 5 xyplo now accepts blank lines (and skips them)
1996 March 21 conversion from inches to cm finally!
origin 1985 dec 24 from linreg *)
(* end module version *)

(* begin module describe.xyplo *)
(*
name
   xyplo: x, y data plotter (pronounced: "zyplo")

synopsis
   xyplo(xyin: in, xyout: out, xyplop: inout, xyplom: in,
         warnings: out, output: out)

files
   xyin: A set of header lines that begin with asterisk ('*') or pound
       sign ('#') are copied to output.  Remaining lines are the data
       in columns, ending with end of file.  Tabs may separate data. 
       Missing columns are not allowed.  See the demonstration file
       xyin.demo for an example.  Once the first data line has been
       read, lines that begin with an '*' or '#" or that are entirely
       blank will be ignored.  This allows one to place comments or
       other information deeper into the file without having xyplo
       object.  The # allows files to be read by gnuplot too.

   xyplop:  Parameters to control the plot, on lines as shown.  The major
       sections of the parameter file are separated by lines that are used by
       the program as separators.  A separator line may begin with blanks, and
       these must be followed by asterisks, as shown below.  These lines simply
       make the file easier to deal with, but you must have them in the file!
       The easiest way to create a xyplop file is to copy the demonstration
       file (xyplop.std or xyplop.demo) and modify that to suite your needs.

       xzero yzero         amounts to move the graph origin (cm)
       zx min max          (character, real, real) if zx='x' then set xaxis
       zy min max          (character, real, real) if zy='y' then set yaxis
                           These two lines set the minimum and maximum range of
                           the data to graph.  Other characters mean the
                           program automatically uses the range of the data.
       xinterval yinterval xsubintervals ysubintervals:
                           These 4 parameters, all on one line, define the
                           number of numbered intervals on each axes to plot
                           and the number of unnumbered subintervals.  Note
                           that on a regular scale, going from 10 to 20
                           requires 10 intervals (to get whole numbers), but on
                           a log scale, going from 10 to 100 requires 9
                           intervals (to get tics every 10), not 10 intervals!

       xwidth    ywidth    width of numbers in characters
       xdecimal  ydecimal  number of decimal places
       xsize     ysize     size of axes in cm
       xlabel              the x axis label
       ylabel              the y axis label
                           NEW: 2016 Sep 29 the y label is now rotated and
                           on the left side of the graph.  To restore
                           the original behavior of the y label above
                           the y axis, put '^' in front of the title.
       zc                  if zc='c' then a crosshairs put on zero of x and y
                                 'x' then only X axis is plotted
                                 'X' then only X axis and crosshairs
                                 'y' then only Y axis is plotted
                                 'Y' then only Y axis and crosshairs
                                 'n' then neither axis nor crosshairs
                                 'N' then neither axis with crosshairs
                                 'i' then numbering and tic marks but no line
                                     (invisible)
                           Otherwise, both axes are plotted without crosshairs.
       zxl base            if zxl='l' then convert the x axes to a log scale
                           using the indicated base.
                           zxl='L' is like 'l' but the numbers given on the
                           axis have not had the log taken of them.
       zyl base            if zyl='l' then convert the y axes to a log scale
                           using the indicated base.
                           zyl='L' is like 'l' but the numbers given on the
                           axis have not had the log taken of them.

       * define columns to read data from ***********************************
                           This section defines which column of xyin contains
                           what kind of data.  You can use a column only once.

       xcolumn   ycolumn   columns of xyin that determine the
                           location of the symbol
       symbol-column       the xyin column to read symbols from
                           if zero, then use the first symbol defined below
       xscolumn  yscolumn  columns of xyin that determine the size  of the
                           symbol.  If zero, then no data is expected.

                           NOTE: for most symbols this is the entire size of
                           the symbol.  For the I beam symbol, the yscolumn is
                           half of the total size plotted.  Thus one may use
                           standard deviations and obtain a symbol of 2
                           standard deviations high centered on the y
                           coordinate.

       colorkind hucolumn sacolumn brcolumn
                           define the hue saturation brightness columns.
                           These control the color of the rectangle symbol.
                           1 0 0 is black (assumed if columns are all zero)
                           1 0 1 is white
                           Hue runs from red (value 0) through the spectrum
                           to violet (value 1).  See the technical notes for
                           further details.

                           If the first character of the line is a
                           digit, blank or 'h' then the hsb (hue
                           saturation brigtness) color model is used.

                           If the first character of the line is 'r'
                           then the rgb (red green blue) color model
                           is used.

       * define one or more symbols *****************************************

                           Each of these sections defines one of the symbols by
                           specifying what to do for each symbol flag seen in
                           the symbol column.  There may be as many symbols as
                           will fit in memory.

                           The last of these sections must contain just a '.'
                           as the 'symbol-to-plot'.  This is required to end
                           the symbol definition section since there are an
                           indefinite number of symbols.

       symbol-to-plot      (character) Most symbols are plotted at the
                           coordinates given in xcolumn and ycolumn.

                           'c' plot a circle
                           'C' plot a circle; with connections in color
                           'b' plot a box
                           'B' plot a box; with connections in color
                           'x' plot an x
                           '+' plot a plus
                           'I' plot an I beam symbol
                           'd' plot a box with central dot
                           'p' point (or dot) alone.
                           'm' mark according to a xyplom definition
                           'M' mark according to a xyplom definition,
                               with connections in color

                           'R' plot a filled rectangle in color.  Unlike the
                           other symbols, which are centered on the data, the
                           lower right hand corner of this rectangle is placed
                           on the data.  This allows the user more control on
                           placement.

                           's' plot a histogram from zero to the height

                           'r' like 'R' but gray scale.  The brightness
                           column is used for controlling the brightness.

                           'f' Means to plot the symbol-flag (defined below).
                           The 'f' type allows several symbols to be made each
                           with its own regression and connection lines, but
                           plotted with the entire flag string in xyin.  The
                           symbols are distinguished by their first character.
                           The symbol-flag in xyplop should be set to the string
                           that one desires to be recognized.
                           'f' will center the string.
                           'F' will left justify the string.

                           'g' Means 'grab bag'.  The 'g' type has lower
                           priority than any other symbol.  Xyplo searches
                           through all the available symbols looking for a
                           match to the symbol-flag.  If a symbol-flag cannot
                           be found, then the data are assigned to the
                           'grab-bag'.  The program uses the symbol-flag on the
                           graph.
                           The symbol-flag in xyplop can be anything.
                           'g' will center the string.
                           'G' will left justify the string.

                           'L' will make just lines using the current colors.

                           The symbol underscore (_) in xyin is converted
                           to a blank to allow the appearance of separated
                           words.

                           One can do grab-bag connected curves without symbols
                           by setting g and the symbol-flag to ' '.  One can
                           also set the symbol-to-plot to blank (or other
                           unrecognized symbol) to get specific connected
                           curves.  In this case, the symbols MUST be connected
                           or the program will object (invisible symbol and
                           invisible connection means data loss).

                           'm' Means to look up the symbol-flag name
                            in the xyplom file and to use the PostScript
                           definition there to create the symbol.

                           'M' like 'm' but draw the connecting line
                           in color and if there is no symbol in the
                           xyplom, don't draw any symbol.  This allows
                           one to make curves without symbols.

       symbol-flag         The string of characters that indicates that this
                           symbol should be plotted.  Eg, if the
                           'symbol-to-plot' is I and the flag is x, then
                           whenever an x is seen in the symbol column, an I
                           beam will be plotted.  The flag can be more than one
                           character long, but (unfortunately) it cannot
                           contain blanks.

       symbol-sizex        Side in cm on the x axis of the symbol.  If this
                           value is negative, the data in xscolumn is used to
                           determine the size.  For circles, sizex determines
                           the diameter, sizey is ignored.

       symbol-sizey        Side in cm on the y axis of the symbol.  If this
                           value is negative, the data in yscolumn is used to
                           determine the size.  For circles, sizeX determines
                           the diameter but a positive number is still required
                           for sizey.

       connection linetype size   If the first character is 'c' then the
                           symbols will be connected by lines of linetype as
                           defined below.  (Linetype must follow the c
                           immediately, without blanks.)

       linetype  size      linetype is a character defining the kind of
                           regression line to plot for this symbol:
                           'l' means do regression line
                           'i' invisible,
                           '.' dotted
                           '-' dashed
                           'n' means no line.

                           '-' and '.' require a size in cm for the
                           spacing.  The others also require a number, but it
                           is ignored.

       * end the symbol definitions with a period (left justified!) *********
.
       * define zero or more user defined lines *****************************
       linetype m b size   One or more lines to be drawn on the plot, m and b
                           are slope and intercept.  Linetype and size are
                           define as for the symbol connection lines.  blank
                           lines and lines that begin with "*" are ignored.
                           in this section
                           linetype is defined as for the regression lines.

       * end of the line definitions and start of more parameters ***********

          edgecontrol edgeleft, edgeright, edgelow, edgehigh:
          edgecontrol is a single character that controls how the bounding
          box of the figure is handled.  If it is 'n' then the bounding box
          will be the page parameters defined in constants inside the program
          (llx, lly, urx, ury AND changes as set by the previous parameter
          line).  If the parameter is 'p', there are four real numbers that
          define the edges around the clist in cm.  To allow a map to be
          imbedded into another figure, its size must be defined in
          PostScript (with %%BoundingBox).  By setting these four numbers,
          the edges are defined.  Negative values are allowed, so one may
          move the edges as desired.

          (New as of 2003 Aug 22)

   xyout: regression results, ready for PostScript input.  (See technical
          notes.)

   xyplom: A file containing definitions of additional symbols,
      written in PostScript.  The symbols must be of the form:

             /mysymbol { % xsize ysize mysymbol -
                % xsize x size in points
                % ysize y size in points
                /xsize exch def
                /ysize exch def
                ...
             } bind def

      Xyplo translates to the proper location and sets the color before
      calling the routine.  This allows the xyin file to control the color of
      the symbols, but of course it can be overridden by the user's routine.
      It is even possible for the symbol to use both colors!

      Be sure to define the symbol in the xyplop.

      If the xyplom file is empty, there are no user defined symbols.

      FONT CONTROL (New as of 2010 Oct 26)
         %F  If the xyplom file contains any line starting with '%F'
         then the user defines the font.

         %T If the xyplom file contains any line starting with '%T'
         then the user can scale the tic marks on both axes.

         %X If the xyplom file contains any line starting with '%X'
         then the user can scale the tic marks on the X axis.

         %Y If the xyplom file contains any line starting with '%Y'
         then the user can scale the tic marks on the Y axis.

         %x If the xyplom file contains any line starting with '%x'
         then the user can move the x label up or down (cm).

         %y If the xyplom file contains any line starting with '%y'
         then the user can move the y label up or down (cm).

         %W If the xyplom file contains any line starting with '%W'
         then a white rectangle is written on the background.  If
         figure is converted to a jpg using the ImageMagick convert
         program, the background will not be black.
         The white rectangle is determined from the bounding box.
         You can control the bounding box using edgecontrol.

         Example:

            %F % font control from xyplom
            /Helvetica-Bold findfont
            16 scalefont
            setfont
            2 setlinewidth % sets linewidth
            %T 1.5 % tic mark scale factor

   warnings: extensive warning messages.  These get in the way if they go to
      output.  One line appears on the output file if there are any
      warnings.

   output: messages to the user

description
   The data in the xyin file are converted to graphics in the PostScript
   language on the xyout file, under control of the parameters set in xyplop.
   There are several distinct sections of the parameters:

   1. The first set of parameters determine the overall characteristics of the
      graph.
   2. The second set of parameters defines the columns of xyin to be read.
   3. The next section of the parameter file defines one or more symbols to be
      plotted on the graph.  If desired, a linear regression is performed
      between the data columns, and this may be graphed for each symbol.  The
      invisible option allows one to obtain the regression data without the
      graph.  Regression data include the correlation coefficient
      and Fisher's z'.
   4. A section with just a period ends the symbols section.
   5. The last section contains lines you define.

   Recommended procedure for using xyplo: obtain a copy of xyplop.demo and
   xyin.demo, set permission to read them for yourself (on a Unix system use
   chmod), and copy them to the names xyplop and xyin.  Try them out as is.  If
   you don't get a graph, doing your own data will not do any good!  Then
   convert the xyplop to your own use by changing the xyplop.demo file and
   substitute your xyin file.  This way the complexity of xyplop can be held at
   bay.

see also

   {Basic example files:}
   xyin, xyout, xyplop, xyplom,

   {Standard parameter file: } xyplop.std

   {Demonstration examples: }
   xyplop.demo, xyin.demo,
   xyplop.test, xyin.test,
   xyplop.mul, xyin.mul,
   xyin.logs
   xyplop.xn xyplop.xl xyplop.xL
   xyplop.yn xyplop.yl xyplop.yL
   xyin.genbank xyplop.genbank

   {Related program - generate histogram:}
   genhis.p

   {Related programs - graphics routines:}
   dops.p
   doodle.p

   {Technical note:
   To define the bounding box for the graph must be defined.
   The postscript program}
   https://alum.mit.edu/www/toms/ftp/printerarea.ps
   {will give the values for your printer.}
   {substitute in the values at %%BoundingBox}

   {Xyplo now accepts the '#' as the start of a comment
   so that data files can also be used for gnuplot.  See:}
   http://www.gnuplot.info/

   {Confidence Limits for the Correlation,
   Fisher's z':}
   http://sportsci.org/resource/stats/sscorr.html#fisherz

   {who says that to get the confidence limits of the correlation}
   {coefficient,}

     {"use the Fisher z transformation: z = 0.5log[(1 + r)/(1 - r)]. }
     {The transformed correlation (z) is normally distributed with}
     {variance 1/(n - 3), so the 95% confidence limits are given by z}
     {± 1.96/sqrt(n - 3).  You then have to back-transform these}
     {limits to correlation coefficients using the equation r =}
     {[(e^(2z) - 1)/(e^(2z) + 1)]."}

   {See also:}
   {Fisher RA (1921). On the probable error of a coefficient of
   correlation deduced from a small sample. Metron 1, 3-32.}
   http://davidmlane.com/hyperstat/A50760.html

author
   Thomas Schneider

technical notes

   The program originally generated output in the pic format.  One
   could then run this through pic and troff to produce a graph. 
   However, the program has been modified to eliminate the pic
   notation (by substituting modules from dops rather than domods). 
   All lines outside the graphics now are preceeded by a %, which is
   beginning of a comment in PostScript.  Thus the output of the
   program can be run directly into a PostScript interpreter.  This
   saves on both memory and speed of graphing since the intermediate
   file is no longer created.

   Colors in PostScript are defined with hue, saturation and
   brightness with the sethsbcolor function.  (Xyplo does not use red,
   green, blue model because it does not give the most useful
   continuous scale, though maybe someday I'll put in a switch.)  The
   standard hue runs from red at 0 to red at 1 with the Roy G. Biv in
   between.  (The famous physicist Roy G. Biv stands for: red, orange,
   yellow, green, blue, indigo, violet, the colors of the spectrum. 
   On 1992 September 16 I realized that, amazingly, the wavelength in
   nm match quite nicely!

      700 red
      750 orange
      600 yellow
      650 green
      500 blue
      550 indego
      400 violet

   Ok, Ok, back to xyplo!  Since the hue runs in a circle, the
   spectrum is not generated from the hue range 0 to 1.  Xyplo
   therefore converts the input numbers by multiplying by 0.84 and
   adding 0.16.  This gives the color range from yellow through red
   corresponding to values 0 to 1.  (Note:  adding 0.16 give the range
   from yellow through red, but that is not a spectrum which runs from
   red through violet.)

   As of version 8.50 (1996 March 21) the xyplo program now uses cm
   (YEA!) which means that ALL previous graphs are out of date.  If
   you happen to be stuck in a backwards, primitive country, tough
   luck.  Bite the bullet.

   Minor unobvious things have prevented people from getting graphs. 
   Most problems occur when badly formed xyplop files are used, and
   the program has no way to tell what the difficulty is.  More checks
   have been put it, so the program can detect most oddly formed xylop
   and xyin files.  Check your xyplop carefully.

   Setting the Default Printer Area:

   The size of the page in PostScript is determined with 4 constants
   llx, lly, urx and ury.  These must be set correctly for each
   printer.  These can be easily obtained by printing the file:

   https://alum.mit.edu/www/toms/ftp/printerarea.ps

   In this program these parameters are set as 4 constants defaultllx,
   defaultlly, defaulturx and defaultury.  The parameters for page
   edges define whether to use the defaults or to compute the size to
   show.

bugs

ENHANCEMENTS:
* xyplo should apply a PostScript to clip everything outside the drawing area 

* xyplo has a numerical drift problem.  For extremely large numbers of
data points, the computed position differes from the position rendered
by PostScript.  This should be corrected by using gsave and restore,
but it is not clear how to modify the draw routines to do this
cleanly.  As it is the program tracks the location of where it THINKS
it is on the graphical display.

ENHANCEMENTS:
* xyplo: cannot handle error bars in log mode

* xyplo: user defined lines have to be type l??

* xyplo: define symbols in postscript - should speed graphics a lot

* xyplo bug: when the last line of the file is empty, the program halts (that's
fine).  The count of the location of the error does NOT include the * lines
though!  Yet it says:  "at line 27 of data (INCLUDING * lines)" Correct this.

* xyplo bug from Mark (email 1992 May 27)

*)
(* end module describe.xyplo *)

(* solved bugs:

* MAJOR XYPLO BUG: xyplo bug1992
  negative 0.50 did not show the minus sign- demo
 is in ~/current/xyplo.bug [directory is lost?]
 also in ~/sites/yeast/gcn4/asym/affinity.
 problem is probably in xaxis calls picnumber, leading to stringreal.
 with width 6 decimal 2 , a negative number like -0.40 comes out 0.40,
 -0.80 as 0.80 but (!) -1.20 is -1.20!
2003 Aug 28
Extremely odd - bug found with same description as this one!
~/sites/yeast/gcn4/asym/affinity still exists.
Date: Sep  5 1992 ... but the graph there is ok.
SOLVED!!
Simply:
* set so that there are enough tics
* minus signs of tic marks between -1 and 0 are not shown.
* True for both axes.
* SOLVED 2003 Aug 28 - new stringreal routine in dops.p

* this was probably bug1992:
* xyplo bug:  when program choses range, negative values are lost?
******** Line 928: [Warning 371]  routine not referenced, code deleted: LINEA.
******** Line 546: [Information 279]  Value assigned to Y is never used;
         assignment is eliminated by optimizer.
******** Line 546: [Information 279]  Value assigned to X is never used;
         assignment is eliminated by optimizer.
******** Line 837: [Information 279]  Value assigned to Y is never used;
         assignment is eliminated by optimizer.
******** Line 1358: [Information 279]  Value assigned to DTHETA is never used;
         assignment is eliminated by optimizer.
******** Line 1699: [Information 279]  Value assigned to NOTEQUAL is never used;
         assignment is eliminated by optimizer.
******** Line 2688: [Warning 243]  Variable Y was not initialized before this
         use.
******** Line 2688: [Warning 243]  Variable X was not initialized before this
         use.

*)

(* begin module xyplo.const *)
   dsafety = 5; (* decimal place safety; assures at least this many
                   decimal places are printed for the line regression data *)

   (* distances in cm *)
   YaxisXshift = -0.127; (* for Y axis, shift numbers in X direction *)
   YaxisYshift = -0.127; (* for Y axis, shift numbers in Y direction *)
   Yticlength  = +0.127; (* Y axis tic mark length *)
   XaxisXshift = +0.000; (* for X axis, shift numbers in X direction *)
   XaxisYshift = -0.508; (* for X axis, shift numbers in Y direction *)
   Xticlength  = +0.127; (* X axis tic mark length *)
   labelabovey =  0.508; (* amount above the Y axis to put the label *)
   labelbelowx =  1.143; (* amount below the X axis to put the label *)

   (* default bounding box definitions.  See the technical notes for how to
      set this for a different printer. *)

   (* The following bounding box is for the Canon Color Laser Copier 1150. *)
   defaultllx =   7.10999; (* default for llx, lower left x *)
   defaultlly =   7.01995; (* default for lly, lower left y *)
   defaulturx = 588.15;    (* default for urx, upper right x *)
   defaultury = 784.98;    (* default for ury, upper right y *)

   (* PostScript constants *)
   pwid = 8; (* width in character places to print PostScript numbers *)
   pdec = 5; (* decimal places to print PostScript numbers *)
{2010 Aug 09: yyy this did not help the drift problem!}
{
   pwid = 10; (* width in character places to print PostScript numbers *)
   pdec = 7; (* decimal places to print PostScript numbers *)
}

   maxinterval = 1000; (* upper bound on xinterval yinterval, a safety net *)
(* end module xyplo.const *)

(* begin module xyplo.interact.const *)
      maxstring = 300; (* the maximum string *)
(* end module xyplo.interact.const *)

(* begin module pic.const *)
   pi = 3.14159265354; (* circumference divided by diameter of circle *)
   picwidth = 8; (* width of numbers printed to the file *)
   picdecim = 5; (* number of decimal places for numbers *)
   charwidth = 0.15875; (* the width of characters in cm (ie, cm/char)
                        this allows centering of strings. *)
   (* note: for the Times-Roman  font, 0.0625 is a good value.
            for the Courier-Bold font, 0.08   is a good value. *)
   dotfactor = 0.015875; (* the size of dots *)
{
   defscale = 72; (* default scale factor.  coordinate units per in *)
}
   defscale = 28.35; (* default scale factor.  coordinate units per cm *)
{ making this change would be a big shock to all the programs that
use it, unfortunately.  A major user is xyplo. }
(* end module pic.const version = 2.80; (@ of dops.p 2010 Nov 03 *)

type
(* begin module interact.type *)
(* begin module string.type *)
   stringptr = ^string; (* pointer to a string *)
   string = record (* a string of characters *)
      letters: array[1..maxstring] of char; (* the letters in the string *)
      length: integer; (* the number of characters in the string *)
      current: integer; (* the letter we are working on *)
      next: stringptr; (* the next string in a series *)
   end;
(* end module string.type version = 5.27; (@ of prgmod.p 2005 Aug 06 *)
(* end module interact.type version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module xyplo.type *)
      line = record (* define a line *)
         linetype: char;     (* i invisible, - dashed, . dotted, l means line,
                                n means no line, see drawa *)
         m, b: real;         (* slope and intercept as in y=m*x+b *)
         linesize: real;     (* size of dash or dots for line in cm *)

         (* internal records on the line. see procedure  regress for defs. *)
         sumx,sumy,sumxsqd,sumysqd,sumxy,
              ex,ey,varx,vary,covxy,r: real;
         n: integer;
      end;

      linesptr = ^lines; (* a pointer to a set of lines *)
      lines = record (* a set of line *)
         l: line; (* one line *)
         next: linesptr; (* the following lines *)
      end;

      symbol = record (* define a symbol *)
         symboltoplot: char;  (* symbol to plot: character c=circle, b=box
                                x, +, I=Ibeam, d=dotted box, p=point,
                                f=use flag, g = grab-bag from flag,
                                m = call postscript in xyplom,
                                M = call postscript in xyplom in color *)
         symbolflag: string;  (* what characters to recognize for this symbol *)
         symbolxsize: real;   (* x diameter or side in cm, if negative
                              convert the data in xscolumn to the size *)
         symbolysize: real;   (* y diameter or side in cm, as for x *)

         doconnection: boolean; (* connect the symbols with connecttype *)
         connecttype: char;     (* i invisible, - dashed, . dotted,
                                l means line, n means no line, see drawa *)
         connectsize: real;  (* spacing on the connected line dashes or dots *)
         didlastpoint: boolean; (* true if the previous point was plotted
                                   this implies that we can draw from that point
                                   to the present point *)
         oldx, oldy: real;   (* the previous coordinate for connecting *)

         doline: boolean;    (* do regression line *)
         l: line;            (* the data on the line for this symbol *)

         symbolhasxyplomdef: boolean; (* The symbolflag is found in xyplom *)
      end;

      symbolsptr = ^symbols; (* a pointer to a set of symbols *)
      symbols = record (* the defined  symbols *)
         s: symbol; (* one symbol *)
         next: symbolsptr; (* the following lines *)
      end;

      xyploparam = record (* parameters to control the plot *)

         parameterversion: real; (* parameter version number *)

         (* Definitions for the BoundingBox of the encapsulated
            PostScript (eps): *)
         llx, lly,         (* lower left x and y *)
         urx, ury:         (* upper right x and y *)
            real;

         edgecontrol: char; (* if 'p' then use page instead of edges *)
         edgeleft : real; (* left  edge margin *)
         edgelow  : real; (* lower edge margin *)
         edgeright: real; (* right edge margin *)
         edgehigh : real; (* high  edge margin *)

         setx, sety:         (* the user sets x and y max and min *)
            boolean;
         xmin, xmax,         (* minimum and maximum for x axis plot *)
         ymin, ymax,         (* minimum and maximum for y axis plot *)
         xscale, yscale:     (* the scale factors to multiply x and y by
                                to make them fit the graph *)
            real;
         xzero, yzero: real; (* the location of the zero coordinate of
                                the final graph (cm) *)
         xinterval,yinterval,(* number of intervals on x and y to plot *)
         xsubintervals,ysubintervals:(* number of subintervals on x and y to plot *)
            integer;
         xwidth,   ywidth,   (* width of numbers on graph in characters *)
         xdecimal, ydecimal: (* number of decimal places for numbers *)
            integer;
         xsize,    ysize:    (* size of axes in cm *)
            real;
         xlabel,             (* the x axis label *)
         ylabel:             (* the y axis label *)
            string;
         rotateylabel: boolean; (* rotate the y axis label instead
            of putting it on the top of the y axis *)

         crosshairs: boolean; (* if true then cross hairs put on zero
                                of x and y *)
         doXaxis: boolean;    (* true if the x axis should be plotted *)
         doYaxis: boolean;    (* true if the x axis should be plotted *)
         doaxisline: boolean; (* line on axis is plotted *)

         logxscale: boolean; (* if true then convert x scale to log scale *)
         logxnormal: boolean; (* if true then print x scale without taking log *)
         xbase: real;        (* the base to use for the log x scale *)
         logxbase: real;      (* the natural log of base, for converting x
                                values to log to the base by division *)
         logyscale: boolean; (* if true then convert y scale to log scale *)
         logynormal: boolean; (* if true then print y scale without taking log *)
         ybase: real;        (* the base to use for the log y scale *)
         logybase: real;      (* the natural log of base, for converting y
                                values to log to the base by division *)

         xcolumn,  ycolumn,  (* columns of input to chose *)
         scolumn,            (* the xyin column to read symbols from *)
         xscolumn,  yscolumn, (* columns of xyin that determine the size
                                 of the symbols *)
         hucolumn, sacolumn, brcolumn (* columns of xyin that determine the
                                 color of the symbols *)
            : integer;
         (* these are true when at least one variable needs
            the corresponding column data *)
         colorkind: char; (* kind - digit, ' ' or h = hsb; 'r' = rgb.
            These change the meaning of the color columns. *)
         needscolumn, needxscolumn, needyscolumn: boolean;
         needcocolumn: boolean;

         usersymbols: symbolsptr; (* all the symbols *)
         userlines: linesptr;(* a set of user defined lines *)

         xrect, yrect: real; (* current sizes of rectangles.  These are
            kept track of so that rectinit only will be called when
            the size changes *)
         startrectangles: boolean; (* When the parameters are read, this flag
            is set, it is true if a rectangle symbol was found.  The code for
            starting a rectangle is written later. *)

         usermarks: boolean; (* if true, user definitions were copied from
            xyplom to the xyout *)
         warnedaboutusermarks: boolean; (* we warned that there were no
            usermarks so don't do it again. *)

         fontcontrol: char; (* 't' = times, 'c'= courier, n= define in xyplom *)
         ticscale: real;  (* scale factor for both X and Y tic marks *)
         Xticscale: real; (* scale factor for X tic marks *)
         Yticscale: real; (* scale factor for Y tic marks *)
         Xlabelshift: real; (* X label shift vertically *)
         Ylabelshift: real; (* Y label shift horizontally *)
      end;
(* end module xyplo.type *)

var
(* begin module xyplo.var *)
   xyin, (* input data *)
   xyout, (* output regression graph *)
   xyplom, (* file from which to read the user defined marks *)
   xyplop, (* file from which to read the parameters *)
   warnings: text; (* warnings *)
   warningcount: integer; (* the number of warnings *)
   datacount: integer; (* count of data items to make warning more useful *)
   uniquestring: string; (* a unique string for symbolflag *)
       (* this variable was used 2011 Jul 31 for debugging; it may be
          removed now if desired.  However, it might still be useful
          to change the definitions given by the user in xyplom
          to be more unique in the xyout file.  That is, while the xyplom
          is being copied, this string would be added to the end of the
          definition and the calls would also be modified that way.
          So a symbolflag 'a' would become aUNIQUE which would not
          be linkely to cause a problem, while 'a' is already used. *)
(* end module xyplo.var *)

(* begin module pic.var *)
   inpicture: boolean; (* true if we are drawing the picture,
                          ie, startpic has been called *)
   picxglobal, picyglobal: real; (* absolute location in the graph *)
   pictolerance: real; (* 10 raised to the picwidth,
      to detect values close to zero *)

   scale: real; (* scale factor.  graphic coordinate units per cm *)

(* NONSTANDARD for efficient use of postscript, keep track of
whether there is a current path *)
   inpath: boolean;

(* NONSTANDARD keep track of number of segments drawn so that
they can be stroked.  This (probably) solves the problem of the
Apple printer dying because it can't handle the data. *)
   segments: integer;

   xsideold, ysideold: real; (* current size of a rectangle.  see rectsize *)
(* end module pic.var version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module halt *)
procedure halt;
(* stop the program.  the procedure performs a goto to the end of the
   program.  you must have a label:
      label 1;
   declared, and also the end of the program must have this label:
      1: end.
   examples are in the module libraries.
   this is the only goto in the delila system. *)
begin
      writeln(output,' program halt.');
      goto 1
end;
(* end module halt version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module skipblanks *)

(* 2003 July 31: tab is considered a blank character *)

function isblank(c: char): boolean;
(* is the character c blank or tab? *)
const
   tab = 9; (* tab character *)
begin
   isblank := (c = ' ') or (ord(c) = tab)
end;

procedure skipblanks(var thefile: text);
(* skip over blanks until a non-blank, or end of line, is found *)
begin
   while isblank(thefile^) and not eoln(thefile) do get(thefile);
end;

procedure skipnonblanks(var thefile: text);
(* skip over nonblanks until a blank, or end of line, is found *)
begin
   while (not isblank(thefile^)) and not eoln(thefile) do get(thefile);
end;

procedure skipcolumn(var thefile: text);
(* skip over a data column *)
begin
   skipblanks(thefile); skipnonblanks(thefile)
end;
(* end module skipblanks version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module copyaline *)
procedure copyaline(var fin, fout: text);
(* copy a line from file fin to file fout *)
begin (* copyaline *)
      while not eoln(fin) do begin
         fout^ := fin^;
         put(fout);
         get(fin)
      end;
      readln(fin);
      writeln(fout);
end; (* copyaline *)
(* end module copyaline version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module interact.clearstring *)
(* begin module clearstring *)
(* These modules clear strings in various ways *)
(* ---- *)
procedure emptystring(var ribbon: string);
(* empty the contents of the string but do NOT
remove the pointer.  This is useful for clearing
one string within a linked list of them. *)
var   index: integer; (* to the ribbon *)
begin (* clearstring *)
   with ribbon do begin
      for index := 1 to maxstring do letters[index] := ' ';
      length := 0;
      current := 0;
   end
end; (* emptystring *)
(* ---- *)
procedure clearstring(var ribbon: string);
(* empty the string and remove the pointer *)
begin (* clearstring *)
   with ribbon do begin
      emptystring(ribbon);
      next := nil;
   end
end; (* clearstring *)
(* ---- *)
procedure initializestring(var ribbon: string);
(* start the string with a nil pointer.   This routine
should be called before doing linked list work.  This allows
the standard string routines to clear the string without
killing the pointer.  This is now deprecated, do not use it
since clearstring still clears the next pointer. *)
begin (* initializestring *)
   writeln(output,'remove initializestring routine!');
   writeln(output,'replace it with clearstring routine!');
   halt; (* to force deprecation *)
   clearstring(ribbon);
   ribbon.next := nil;
end; (* initializestring *)
(* end module clearstring version = 5.27; (@ of prgmod.p 2005 Aug 06 *)
(* end module interact.clearstring version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module interact.getstring *)
procedure getstring(var afile: text; var buffer: string;
                    var gotten: boolean);
(* get a line (as a string) from a file not using string calls.  this lets
one obtain lines from a file without interactive prompts *)
var   index: integer; (* of buffer *)
begin (* getstring *)
      clearstring(buffer);
      if eof(afile)
      then gotten := false
      else begin
         index := 0;
         while (not eoln(afile)) and (index < maxstring) do begin
            index := succ(index);
            read(afile, buffer.letters[index])
         end;
         if not eoln(afile) then begin
            writeln(output, ' getstring: a line exceeds maximum string size (',
                            maxstring:1,')');
            halt
         end;
         buffer.length := index;
         buffer.current := 1;
         readln(afile);
         gotten := true
      end
end; (* getstring *)
(* end module interact.getstring version = 5.22; (@ of prgmod.p 2005 Jul 12 *)

(* begin module interact.writestring *)
(* begin module writestring *)
procedure writestring(var tofile: text; var s: string);
(* write the string s to file tofile, no writeln *)
var   i: integer; (* index to s *)
begin (* writestring *)
      with s do for i := 1 to length do write(tofile, letters[i])
end; (* writestring *)
(* end module writestring version = 5.27; (@ of prgmod.p 2005 Aug 06 *)
(* end module interact.writestring version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module linear.regression *)
(* linear regression procedure
origin: 1980 august 5
      by thomas schneider
      copyright 1986 *)

procedure regress ( control: char; (* described below *)
                    x, y: real;    (* the data pairs *)
                var sumx, sumy, sumxsqd, sumysqd, sumxy, (* internal records *)
(* mean, variance and covariance of x and y: *)
                    ex, ey, varx, vary, covxy,
(* correlation coefficient: *)
                    r,
(* m = slope, b = y intercept as in y = mx + b *)
                    m, b
                        : real;
(* number of data pairs entered: *)
                var n: integer);
(* linear regression for variables x and y.
   the control variable has only three acceptable states:
      c clear all variables to zero
        (note that r is set to 2.0, so that a program can check whether
         or not the results have been calculated yet. the recommended test
         is not for equality, since round-off could cause problems.  use
              if r > 1.5 then (no calculations yet)
                         else (calculations done)    )
      e enter the x and y values into the internal sums
      r results are calculated from the sums
*)
begin
      if not( control in ['c','e','r']) then begin
         writeln(output,' linear regression control variable value, "',control,
                        '" is not acceptable.');
         writeln(output,' it must be in ["c","e","r"]');
         halt
      end

      else case control of
         'c': begin (* clear *)
            x:=0.0; y:=0.0;
            sumx:=0.0; sumy:=0.0;
            sumxsqd:=0.0; sumysqd:=0.0;
            sumxy:=0.0;
            ex:=0.0; ey:=0.0;
            varx:=0.0; vary:=0.0;
            covxy:=0.0;
            r:=2.0; (* this (impossible) value can act as a flag that
                       shows that the calculations have not been made. *)

            m:=0.0; b:=0.0;
            n:=0
         end;

         'e': begin (* enter data *)
            sumx:=sumx + x;
            sumy:=sumy + y;
            sumxsqd:=sumxsqd + x*x;
            sumysqd:=sumysqd + y*y;
            sumxy:=sumxy + x*y;
            n:=succ(n)
         end;

         'r': begin (* calculate results *)
            (* check for conditions that would bomb the program *)
            if n = 0 then begin
               writeln(output,'regress: n is 0; no samples recorded');
               halt
            end;
            if n = 1 then begin
               writeln(output,'regress: n is 1; regression impossible');
               halt
            end;

            ex:=sumx/n;
            ey:=sumy/n;

            varx:=(sumxsqd/n) - ex*ex;
            vary:=(sumysqd/n) - ey*ey;
            if varx=0.0 then begin
               writeln(output,'regress: variance of x is zero;',
                              ' regression impossible');
               halt
            end;
            if vary=0.0 then begin
               writeln(output,'regress: variance of y is zero;',
                              ' regression impossible');
               halt
            end;

            covxy:=(sumxy/n) - ex*ey;

            r:=covxy/(sqrt(varx * vary));

            m:=covxy/varx;
            b:=ey - m*ex
         end
      end
end;
(* end module linear.regression version = 'matmod 1.98 86 nov 23 tds/gds'; *)

(* begin module pic.functions *)
(* ********************************************************************** *)
(* begin module pic.await *)
procedure await;
(* Wait for user to type a carriage return.  the routine assumes that
there is a global file called input. *)
begin
(* the infinite way:
    writeln(output);
    writeln(output,'*********************************');
    writeln(output,'* Use control-c to kill program *');
    writeln(output,'*********************************');
    while true do begin end;*)
end;
(* end module pic.await version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.startpic *)
procedure startpic(var afile:text; setscale,x,y: real; thefont: char);
(* open the graphics field, with the given scale, and at (x,y)
in that scale.  scale is in device coordinates per cm.
The font is chosen with thefont; t = Times-Roman, c = Courier-Bold *)
(* start pic output to file afile, set the globals *)
(* NONSTANDARD *)
(* this is the actual "world" coordinates used: *)
(*                        xmin,  xmax,  ymin, ymax *)
(*    ns; if (setwindow(-5.0/scale, +5.0/scale, -5.0/scale, +5.0/scale)*)
begin
   writeln(afile,'gsave'); (* save the current graphics state *)

{2005 Aug 6: get rid of these finally
   writeln(afile,'% initgraphics'); (* make sure the printer
           is ready to print, without this, sometimes an Apple laserwriter
           will print the graph upside down, tiny and backwards! *)
   writeln(afile,'% clear erasepage'); (* clean residue from before *)
}

   scale := setscale; (* set the global scale *)

   case thefont of
      'c': begin
           writeln(afile,'/Courier-Bold findfont'); (* locate the font *)
           writeln(afile,10:1,' scalefont'); (* set the font size in points*)
           writeln(afile,'setfont'); (* put the font into the current font *)
      end;
      't': begin
           writeln(afile,'/Times-Roman findfont'); (* locate the font *)
           writeln(afile,12:1,' scalefont'); (* set the font size in points*)
           writeln(afile,'setfont'); (* put the font into the current font *)
      end;
      'n': begin
           writeln(afile,'% /Times-Roman findfont'); (* locate the font *)
           writeln(afile,'%', 12:1,' scalefont'); (* set the font size in points*)
           writeln(afile,'% setfont'); (* put the font into the current font *)
      end;
   end;


(* If the following statement is done then it will work on the sun, but will
kill the applewriter!!!!

Sun's non-standard PostScript extension, setlinewidth has default 0,
as stated in the Read This First and the NeWS Manual.
This draws very quickly with 1 bit wide lines.
If '1 setlinequality' is not done, then one cannot set the width of lines.
So to use PostScript on the screen, I must first do
'1 setlinequality'.  However, if I send this code to the Applewriter,
it kills PostScript on the Applewriter and I get no output whatsoever!
(It took me several hours to figure this out, since once PostScript is killed
on the Applewriter, the NEXT output is also smashed and I had to figure that
out also...)

So a standard PostScript program will not work correctly with the
default.  "Correcting" the PostScript program so that it works on the
Sun means that it BOMBS on the Applewriter.

The default for setlinequality should be '1 setlinequality' so that
the same PostScript code can be used both on the Sun and with other devices.
If you want speed, use the nonstandard form.
An alternative is to redefine setlinequality so that '0 setlinequality'
does give correct results with standard PostScript.

Please review this, Randy.  I think that Sun should fix it.
   writeln(afile,'1 setlinequality');  makes lines at least 1 bit wide *)

   (* set the scale to cm
   writeln(afile,
                 scale:picwidth:picdecim,' ',
                 scale:picwidth:picdecim,' scale'); *)

   (* define some things in postscript *)
   (* doline allows less stuff to be put in the output file.
      it takes two numbers off the stack, copies them, draws a line
      to them as coordinates. *)
(* replaced by 'currentpoint translate'
   writeln(afile,'/doline { 2 copy lineto } def');
*)
   (* define a function that makes cm out of a number *)
(* do this all internally here, it's faster
   writeln(afile,'/i { ',scale:picwidth:picdecim,' mul} def');
*)

   (* move to the start point on the page *)
   writeln(afile,
                     (x*scale):picwidth:picdecim,
                 ' ',(y*scale):picwidth:picdecim,
                 ' translate');

    writeln(afile);
    writeln(afile,'% Define functions so the text produced is smaller');
    writeln(afile,'/a {stroke newpath 0 0} def     % special for arc');
    writeln(afile,'/c {stroke 0 0 moveto} def      % current point');
    writeln(afile,'/f {findfont 10 scalefont setfont} def');
    writeln(afile,'   % to set fonts simply use the f function.  Example:');
    writeln(afile,'   %/Symbol f (\142) /Courier-Bold f (-galactosidase');
    writeln(afile,'/l {lineto} def');
    writeln(afile,'/m {moveto} def');
    writeln(afile,'/n {stroke newpath 0 0 moveto} def'); (* new segment *)
    writeln(afile,'/rl {rlineto} def');
    writeln(afile,'/rm {rmoveto} def');
    writeln(afile,'/s {newpath 0 0 moveto} def     % Start path ');
    writeln(afile,'/t {currentpoint translate} def % translate ');
    writeln(afile,'/x {show} def                   % show teXt ');
    writeln(afile);

   (* start out the pathway *)
   inpath := false;
   (* start the number of segments written: *)
   segments := 0;

   (* now for the normal pic stuff: *)
   inpicture := true;
   picxglobal := 0.0;
   picyglobal := 0.0;
   pictolerance := trunc(exp(picwidth*ln(10))+0.5)
(*;writeln(output,'pictolerance = ',pictolerance:picwidth:picdecim);*)
end;
(* end module pic.startpic version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.stoppic *)
procedure stoppic(var afile:text);
(* stop pic output to file afile *)
(* NONSTANDARD *)
begin
    if inpath then begin
       writeln(afile,'stroke');
       inpath := false
    end;

    writeln(afile,'showpage');
    writeln(afile,'grestore'); (* restore the current graphics state
                                  to what it was before the startpic *)

    await;

    inpicture := false;
end;
(* end module pic.stoppic version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.drawr *)
procedure drawr(var afile: text; dx,dy: real; visibility: char;
               spacing: real);
(* make a line to file afile by relative draw of dx,dy with visibility
  i invisible
  - dashed
  . dotted
  l line
with the dashes or dots separated by the spacing given
(this has no effect with invisible and line). *)
(* NONSTANDARD *)
var
   ddx,ddy: real; (* changes in dx and dy for dots and dashes *)
   dr: real; (* the hypotenuse, the distance actually drawn *)
   on: boolean; (* draw linesegment if true *)
   y: real; (* the variable for tracking dots and dashes *)
   r: integer; (* number of times to cycle for dots and dashes *)
   ss: real; (* precalculated value to make things a bit faster *)
   theta: real; (* angle of the line *)
procedure checkseg(var afile: text);
(* NONSTANDARD checks how many segments have been written, if
more than 'buffer', stroke them to the postscript page *)
const buffer = 10;
begin
   if segments >= buffer
   then begin
   (* New segment: writeln(afile,'stroke newpath 0 0 moveto'); *)
      writeln(afile,'n');
      segments := 0
   end
   else segments := segments + 1;
end;
begin (* drawr *)
   if not inpath then begin
      (* starts from current coordinates *)
      (* Start path: writeln(afile,'newpath 0 0 moveto'); *)
      writeln(afile,'s');
      inpath := true
   end
   else checkseg(afile);
(* checks
if not (visibility in ['l','i','.','-'])
then writeln(afile,'%YELLLLLL!!!',visibility,'!');
writeln(afile,'% ',visibility,' line');*)

   (* put these on the stack, they will always be used *)
   write(afile,    (dx*scale):picwidth:picdecim,
               ' ',(dy*scale):picwidth:picdecim);
   case visibility of
  'l','i': begin
              case visibility of
                 'i': write(afile,' m');
                 'l': write(afile,' l');
              end
           end;
  '.','-': begin (* make up our own dots and dashes *)
              writeln(afile); (* move away from the (dx,dy) on the stack *)
              if spacing <= 0.0 then begin
                 writeln(output,'drawr: spacing zero with . or - line');
                 halt
              end;
              if dx = 0.0
              then begin
                 ddx := 0.0; (* avoid division by zero *)
                 ddy := scale*spacing;
                 if dy < 0 then ddy := - ddy; (* this makes sure that
                    we draw lines straight down if that was the request *)
              end
              else begin
                 (* find out the angle of the slope, intentionally
                    lose the sign *)
                 theta := arctan(abs(dy/dx));

                 ddx := scale*spacing*cos(theta);
                 ddy := scale*spacing*sin(theta);

                 (* return the sign to the little buggers *)
                 if dx < 0 then ddx := -ddx;
                 if dy < 0 then ddy := -ddy;
              end;
              y := 0;
              case visibility of
                 '.': ss := scale*dotfactor;
                 '-': on := true;
              end;

              dr := sqrt(dx*dx+dy*dy);
              for r := 1 to round(dr/spacing) do begin
                 case visibility of
                 '-': begin
                         write(afile,    (ddx):picwidth:picdecim,
                                     ' ',(ddy):picwidth:picdecim);
                         if on
                         then writeln(afile,' rl')
                         else writeln(afile,' rm');
                         on := not on
                      end;
                 '.': begin
                         (* put out a dot like in dotr *)
                         write(afile,
                                    +ss:picwidth:picdecim,' 0 rl');
                         write(afile,' ',
                                    -ss:picwidth:picdecim,' 0 rl');
                         write(afile,' ',(ddx):picwidth:picdecim,
                                     ' ',(ddy):picwidth:picdecim);
                         writeln(afile,' rm');
                      end;
                 end
              end;
              (* let's make really sure we got there!! *)
              writeln(afile,' m'); (* pulled from the stack *)
           end;
   end;

   (* an elegant way to make postscript keep a global record is
   to translate the coordinates! *)
   (* writeln(afile,' currentpoint translate'); *)
   writeln(afile,' t');

   picxglobal := picxglobal + dx;
   picyglobal := picyglobal + dy;
end;
(* end module pic.drawr version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.mover *)
procedure mover(var afile: text; dx,dy: real);
(* move relative the amount (dx, dy). *)
begin
   drawr(afile,dx,dy,'i',0.0);
end;
(* end module pic.mover version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.liner *)
procedure liner(var afile: text; dx,dy: real);
(* draw a line the relative amount (dx, dy). *)
begin
   drawr(afile,dx,dy,'l',0.0);
end;
(* end module pic.liner version = 2.80; (@ of dops.p 2010 Nov 03 *)


(* begin module pic.drawa *)
procedure drawa(var afile: text; x,y: real; visibility: char;
               spacing: real);
(* make a line to file afile to absolute coordinate x,y with visibility
  i invisible
  - dashed
  . dotted
  l line
with the dashes or dots separated by the spacing given
(this has no effect with invisible and line). *)
var
   dx, dy: real; (* differences between current and desired
                    locations *)
begin
   dx := x - picxglobal;
   dy := y - picyglobal;

   drawr(afile,dx,dy,visibility,spacing)
end;
(* end module pic.drawa version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.movea *)
procedure movea(var afile: text; x,y: real);
(* move to absolute x and y *)
begin
   drawa(afile,x,y,'i',0.0);
end;
(* end module pic.movea version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.linea *)
procedure linea(var afile: text; x,y: real);
(* draw a line from current position to absolute x and y *)
begin
   drawa(afile,x,y,'l',0.0);
end;
(* end module pic.linea version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.graphstring *)
procedure graphstring(var tofile: text; var s: string; justification: char);
(* graph the string s.  If it is recognized as a quoted string (surrounded
by double quotes), graph it without the quotes and center it.
Otherwise justify it based on the justification character:
'l' left, 'c' centered, 'r' right.
For right and centered justification, the drawing point is the same
as before the string was done.  For left justification it is at the
right of the string to allow more to be added on there.
If not in picture (global variable inpicture), there is no output.
2008 Aug 05.  '(' and ')' need to be protected by putting '\' in
front of them.
2010 Nov 02 ONLY IF THERE IS NO PROCEEDING '\'!
*)
(* NONSTANDARD: PostScript dependent code.  Since different fonts
have different sized characters, one must rely on the PostScript
to handle the justification of the string. *)
var i: integer; (* index to s, and temporary storage *)
    quoted: boolean; (* true if the string is quoted *)
    skipping: boolean; (* true if skipping leading blanks *)
procedure postscriptprotect(var afile: text; c: char);
(* 2008 Aug 05: protect parenthesis *)
begin
{
2010 Nov 02 DO NOT DO THIS!!!  It messes up beta
   if (c = '(') or
      (c = ')')
   then write(afile,'!\');
Instead, protect only if \ is given.
}
end;
begin
   if (inpicture and (s.length > 0))
   then with s do begin
      if length > 2
      then if (letters[1]='"') and (letters[length]='"')
           then quoted := true
           else quoted := false
      else quoted := false;

      (* override so quoted strings are always centered *)
      if quoted then justification := 'c';

      (* do the non-standard postscript: *)
      if justification <> 'l' then write(tofile,'gsave ');

      (* do postscript to complete previous path *)
      (* set current point: writeln(tofile,'stroke 0 0 moveto'); *)
      writeln(tofile,'c');

      if justification = 'c' then begin
         (* when centering, skip leading blanks *)
         if letters[1] = ' ' then skipping := true
                             else skipping := false;
      end
      else skipping := false;

      write(tofile,'('); (* begin postscript literal *)
      if quoted (* take it literally *)
      then for i := 2 to length-1 do begin
               if letters[i-1] <> '\'
               then postscriptprotect(tofile,letters[i]);
               write(tofile,letters[i]);
           end
      else for i := 1 to length
           do if skipping then begin (* skip leading blanks *)
                 if letters[i] <> ' '
                 then begin
                      skipping := false;
                      if i = 1 then postscriptprotect(tofile,letters[i])
                               else if letters[i-1] <> '\'
                                    then postscriptprotect(tofile,letters[i]);
                      write(tofile,letters[i]);
                 end
              (* else skip the blank by not writing it *)
              end
              else begin
                  if i = 1 then postscriptprotect(tofile,letters[i])
                           else if letters[i-1] <> '\'
                                then postscriptprotect(tofile,letters[i]);
                  write(tofile,letters[i]);
              end;
      write(tofile,')'); (* end postscript literal *)

      if justification = 'c' (* center the string *)
      then write(tofile,' dup stringwidth pop neg 2 div 0 rmoveto')
      else if justification = 'r' (* rigth justify the string *)
      then write(tofile,' dup stringwidth pop neg 0 rmoveto');

      writeln(tofile,' x'); (* show the literal *)
      inpath := false; (* force new path from here *)

      if justification <> 'l' then write(tofile,'grestore ');
   end
(* There is no output if not in picture
   else begin
      writestring(tofile,s);
      writeln(tofile)
   end
 *)
end;
(* end module pic.graphstring version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.stringinteger *)
procedure stringinteger(number: integer; var name: string;
                        width: integer; leadingzeros: boolean);
(* make the string from the number, start putting characters in
after the current length point. use width characters.
if leadingzeros is true, trail zeros before the number. *)
var
   bigdigit: integer; (* the location of the biggest digit *)
   dig: integer; (* number of digits in the number *)
   place: integer; (* place to write the next digit of the number *)
   sign: integer; (* the sign of the number *)
begin
   with name do begin
      if number < 0
      then begin
         sign := -1;
         length := length + 1; (* provide room for the sign!! *)
         number := -number;
         if leadingzeros then begin
            writeln(output,'WARNING: stringinteger: the sign of a negative',
                           ' number with leading zeros is lost');
         end
      end
      else sign := +1;

      (* log 10 of the number plus 1 is the number of digits in the number.
      On this sun computer ln(1000)/ln(10) is 2.9999, which when
      truncated gives 2, rather than the desired 3.  To avoid this
      kind of problem, 0.1 is added. *)
      if number > 9
      then dig := trunc(ln(number+0.1)/ln(10))+1
      else dig := 1;

      if dig > width then begin
         writeln(output,'stringinteger: number width too small');
         writeln(output,dig:1,' digit number (',number:1,')');
         writeln(output,'does not fit in ',width:1,' characters');
         halt
      end;
      if leadingzeros
      then bigdigit := length + 1 (* no sign if leading zeros *)
      else begin
         bigdigit := length + width - dig + 1;

         if (bigdigit <= length) and (sign < 0) then begin
            writeln(output,'stringinteger: no room for sign');
            halt
         end;
      end;
      if sign < 0 then letters[bigdigit-1] := '-';

      for place := length + width downto bigdigit do begin
         case (number mod 10) of
            0: letters[place] := '0';
            1: letters[place] := '1';
            2: letters[place] := '2';
            3: letters[place] := '3';
            4: letters[place] := '4';
            5: letters[place] := '5';
            6: letters[place] := '6';
            7: letters[place] := '7';
            8: letters[place] := '8';
            9: letters[place] := '9';
         end;
         number := number div 10;
      end;
      length := length + width;
   end
end;
(* end module pic.stringinteger version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.stringreal *)
procedure stringreal(number: real; var name: string;
                     width, decimal: integer);
(* make the string from the real number, start putting characters in
at the start point. use width characters and decimal characters
after the decimal place *)
   (* note that the rounding operation to get the digits below zero
      must be done first.  then the digits above zero can be lopped off.
      this makes 99.99 come out correctly to 100.0 (to 1 decimal place)
      otherwise, 99.99 -> 0.99 -> 1.0 (rounded) -> 10 (print with 1 decimal
      place), and stringinteger won't be happy about that.
    2003 Aug 28.  corrected missing minus sign for -1 < number <= 0.
*)
var
   abovezero: integer; (* the number shifted above the decimal place, to
      'decimal' positions (and rounded) *)
   shift: integer; (* power of ten used to shift a number around
      relative to the decimal point *)
   sign: integer; (* the sign of the number *)
   thedecimal: integer; (* integer version of the decimal part of the number *)
   theupper: integer; (* integer version of the upper part of the number *)
   signspot: integer; (* of the spot the sign will go. *)
begin
   (* sanity check: *)
   if name.length + width > maxstring then begin
      writeln(output,'real number =',number:1,
                     ' would exceed maxstring = ',maxstring:1);
      halt;
   end;
   if number < 0 then sign := -1
                 else sign := +1;

   number := abs(number); (* make positive *)

   (* the amount to shift the number above zero *)
   shift := round(exp(decimal*ln(10))); (* amount to move above zero *)
   abovezero := round(number*shift); (* move above zero, round off *)
   theupper := trunc(abovezero/shift);
   thedecimal := abovezero - shift*theupper;
{
writeln(output,'            stringreal:  number  = ',number:pwid:pdec);
writeln(output,'            stringreal:    sign  = ',sign:pwid);
writeln(output,'            stringreal: theupper = ',theupper:pwid);
}

   (* create the actual real number *)
   (* before decimal point *)
   (* provide a space for the sign in the resulting string: *)

   with name do begin (* put in the decimal point *)
      (* force a space for the sign by making the number negative *)
      signspot := length + 1; (* take note of the spot the sign will go. *)
      stringinteger(sign*theupper,name,width-decimal-1,false);

      (* 2003 Aug 28
         There is a very special case, known as bug1992.
         when the (number > -1) and (number < 0) the upper part
         of the number is zero (theupper = 0) BUT as an integer
         the sign cannot be passed to stringinteger, since
         -0 is of course 0 (usually, or sometimes).
         SO we have to handle that case and put a minus sign
         in 'by hand' here. *)
      if (sign < 1) and (theupper = 0)
         (* if number is tic = -0.000000000000000055511151231257827021181583405
            (a real example!!) then we would get -0.0 on rounding.
            SO round to the number of decimal places.
            The number of decimal places:
                  10^decimal = exp(ln(10^decimal)) = exp(decimal*ln(10))
         *)
         and (round(exp(decimal*ln(10))*number) <> 0)
      then begin
         (* ok, starting at signspot, move to the right until
            we are snug up against the number *)
         while letters[signspot+1]=' ' do signspot := succ(signspot);
         letters[signspot]:='-';
      end;

{
write(output,
   'stringinteger(',sign*theupper:1,',"');
writestring(output, name);
write(output, '",', width-decimal-1:1,',',false);
writeln(output,')');
}

      (* put in the decimal point *)
      length := length + 1;
      letters[length] := '.';
   end;
   stringinteger(thedecimal,name,decimal,true); (* after decimal point *)
end;
(* end module pic.stringreal version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.picnumber *)
procedure picnumber(var afile: text;
                    dx, dy, number: real; width, decimal: integer;
                    justification: char);
(* Supply graphic commands for a 'number' whose center is at the relative point
(dx, dy) from the current point, 'width' characters wide and 'decimal'
characters beyond the decimal point.
If the width is zero, no number is produced.
procedure stringnumber(number: integer; start: integer; var name: string);
the location after the call is the same as before the call.
The string is optionally justified: left, centered or right: lcr. *)
var
   name: string; (* the string to pack the number into for shipping out *)
begin
   if width > 0 then begin
      mover(afile,dx,dy);

      clearstring(name);

      if decimal>0
      then stringreal(number,name,width,decimal)
      else stringinteger(round(number),name,width,false);

      graphstring(afile, name, justification);
      mover(afile,-dx,-dy);
   end
end;
(* end module pic.picnumber version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.xtic *)
procedure xtic(var afile: text; length, dx, dy, number: real;
               width, decimal: integer;
               logxnormal: boolean;
               logxbase: real);
(* produce a tic mark for the x axis of "length" long.
Supply a number whose center is at the relative point (dx, dy)
from the end to the tick, 'width' characters wide and 'decimal'
characters beyond the decimal point.
If the width is zero, no number is produced.
the location after the call is the same as before the call.
If logxnormal is true, then raise the number
to logxbase. *)
begin
   liner(afile,0.0,-length);
   if logxnormal
   then picnumber(afile,dx,dy,exp(number*logxbase),width,decimal,'c')
   else picnumber(afile,dx,dy,number,width,decimal,'c');
   mover(afile,0.0,length);
end;
(* end module pic.xtic version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.ytic *)
procedure ytic(var afile: text; length, dx, dy: real;
               number: real;
               width, decimal: integer;
               logynormal: boolean;
               logybase: real);
(* produce a tic mark for the y axis of "length" long.
Supply a number whose right side is started at the relative point (dx, dy)
from the end to the tick, 'width' characters wide and 'decimal'
characters beyond the decimal point.
If the width is zero, no number is produced.
the location after the call is the same as before the call.
If logynormal is true, then raise the number
to logybase. *)
begin
   liner(afile,-length,0.0);

   (* convert the number if we are doing logynormal: *)
   if logynormal then number := exp(number*logybase);

   picnumber(afile,dx,dy,number,width,decimal,'r');
   mover(afile,length,0.0);
end;
(* end module pic.ytic version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.doaxis *)
procedure doaxis(var afile: text;
         theaxis: char;
         doaxisline: boolean;
         alength,fromtic,interval,totic: real;
         subintervals: real;
         length, dx, dy: real;
         width, decimal: integer;
         logscale, lognormal: boolean;
         logbase: real);
(* draw an axis starting from the current position.
Which axis it is is defined by theaxis, 'x' (horizontal) or 'y' (vertical).
Combining the code for both axes into one procedure is a little
slower, but drawing the axis does note ever take significant time,
and this allows improvements to be made on both axes simultaneously.
The length of the axis is alength.
If doaxisline is true then the axis line is drawn.
The axis is labeled with numbers starting with fromtic
at intervals given up to totic.
The remaining variables describe the form of the tic marks as in ytic.
If the width is zero, no number is produced.
the location after the call is the same as before the call.
If logscale and lognormal is true, then raise the tic numbers to logbase.
 *)
var
   half: real; (* half the jump interval.  By adding this to the while loops,
                  we assure that the very last tic gets done, and isn't lost
                  due to roundoff *)
   jump: real; (* the space to move on the graph between tic marks *)
   jumpdistance: real; (* the total jumps made.  this may not be
      a simple function of the input variables since they may
      not work out to an exact number of jumps *)
   tic: real; (* the numerical value of the tic label *)

   dosubtics: boolean; (* do sub tics *)
   subtic: real; (* the numerical value of the (unlabeled) subtic *)
   subinterval: real; (* the numerical interval between subtics *)
   subjump: real; (* the space to move on the graph between subtic marks *)
   halfsubinterval: real; (* half a subjump, see half *)

   currentspot: real; (* current graphing spot *)
   oldspot: real; (* previous graphing spot *)
   axisscale: real; (* axis scaling factor *)

begin
   writeln(afile,'gsave');

{
writeln(output,'In doaxis');
writeln(output,'interval=',interval:10:4);
writeln(output,'subintervals=',subintervals:10:4);
writeln(output,'logbase=',logbase:10:4);
}

   if theaxis = 'x' then begin
      if doaxisline then begin
         liner(afile,+alength,0.0);
         mover(afile,-alength,0.0);
      end;
   end
   else begin
      if doaxisline then begin
         liner(afile,0.0,+alength);
         mover(afile,0.0,-alength);
      end;
   end;

   if totic = fromtic then begin
      writeln(output,'doaxis: ',theaxis,' axis fromtic and totic',
                     ' cannot be equal');
      halt;
   end;

   if (alength = 0.0) or (interval = 0.0) then begin
      writeln(output,'doaxis: neither ',
                     theaxis,' axis length nor interval can be zero');

      halt;
   end;

   axisscale := alength / (totic - fromtic);
   jump := axisscale * interval;
   jumpdistance := 0;
   half := interval / 2.0;

   if subintervals > 1 then begin
      dosubtics := true;
      subinterval := interval/subintervals;
      halfsubinterval := subinterval / 2.0;
      subjump := jump/subintervals;
   end
   else begin
      dosubtics := false;
      subinterval := 0;
      halfsubinterval := 0;
      subjump := 0;
   end;
{
writeln(output,'fromtic  = ',fromtic:10:4);
writeln(output,'totic    = ',totic:10:4);
writeln(output,'interval = ',interval:10:4);
writeln(output,'half     = ',half:10:4);
}

   tic := fromtic;
   if interval > 0.0 then while tic <= totic+interval do begin
{
writeln(output,'* tic=',tic:10:4);
}
      (* adding interval/2 makes sure that the y axis does include
         the last point but that the x axis doesn't get an extra
         from (presumably) rounding effects *)
      if tic <= totic+(interval/2) then begin
         if theaxis = 'x'
         then xtic(afile,length,dx,dy,tic,width,decimal,lognormal,logbase)
         else ytic(afile,length,dx,dy,tic,width,decimal,lognormal,logbase);
      end;

      (* 2007 Aug 30 the extra interval makes subtics go to the end
         of the graph rather than ending at the last tic mark *)
      if tic <= totic + interval then begin
{
writeln(output,'totic+half + interval=',totic+half + interval:10:4);
writeln(output,'totic+interval=',totic+interval:10:4);
      if tic <= totic+half + 1 then begin
writeln(output,'TIC=',tic:10:4);
writeln(afile,'% tic=',tic:10:4);
mover(afile,0.05,0.0);
}
         if dosubtics then begin (* do subtic marks *)
            if logscale then begin (* do subtic marks on log scale *)
               (* subtic starts as a "normal" number (ie, no log taken) at tic: *)
{
writeln(output,'2^tic=',exp(tic*logbase):10:4);
writeln(output,'2^(tic+interval)=',exp((tic+interval)*logbase):10:4);
}
               subtic := exp(tic*logbase);
               (* subtic will proceed to the same but at tic+interval.
                  We divide that into the subintervals. *)
{
writeln(output,'halfsubinterval=',halfsubinterval:10:4,' original');
}
               subinterval := (exp((tic+interval)*logbase) - subtic)/subintervals;
               halfsubinterval := subinterval/2.0;

{
writeln(output,'subtic=         ',subtic:10:4);
writeln(output,'subinterval=    ',subinterval:10:4);
writeln(output,'halfsubinterval=',halfsubinterval:10:4);
}

               oldspot := axisscale * tic;
               while subtic < exp(logbase*(tic+interval)) - halfsubinterval
               do begin

                  (* although tic is on a log scale,
                     we have to have subtic on the regular scale
                     to alter the positions of the subtics *)

                  (* if subinterval is constant,
                     the following makes linearly spaced marks: *)
                  subtic := subtic + subinterval;

                  (* the actual jumps have to be in the log form: *)
                  currentspot := axisscale*ln(subtic)/logbase;
                  subjump := currentspot - oldspot;
{
writeln(output,'       SUBTIC=',subtic:10:4);
writeln(output,'       ln(SUBTIC)/logbase=',ln(subtic)/logbase:10:4);
writeln(output,'       currentspot=',currentspot:10:4);
writeln(output,'       subjump=',subjump:10:4);
writeln(output,'       oldspot=',oldspot:10:4);
writeln(afile,'% subtic=',subtic:10:4);
}
                  oldspot := currentspot;


                  if theaxis = 'x' then begin
                     xtic(afile,length/2,dx,dy,0,0,0,lognormal,logbase);
                     mover(afile,subjump,0.0);
                  end
                  else begin
                     ytic(afile,length/2,dx,dy,0,0,0,lognormal,logbase);
                     mover(afile,0.0,subjump);
                  end;

                  jumpdistance := jumpdistance + subjump;
               end
            end
            else begin (* do subtic marks on regular scale *)
               subtic := tic;
               while subtic < tic+interval-halfsubinterval do begin
                  subtic := subtic + subinterval;

                  if theaxis = 'x' then begin
                     mover(afile,subjump,0.0);
                     if subtic <= totic
                     then xtic(afile,length/2,dx,dy,0,0,0,lognormal,logbase);
                  end
                  else begin
                     mover(afile,0.0,subjump);
                     if subtic <= totic
                     then ytic(afile,length/2,dx,dy,0,0,0,lognormal,logbase);
                  end;

                  jumpdistance := jumpdistance + subjump;
               end
            end
         end
         else begin (* do regular tic marks *)
            if theaxis = 'x' then mover(afile,jump,0.0)
                             else mover(afile,0.0,jump);
            jumpdistance := jumpdistance + jump
         end
      end;
      tic := tic + interval;
   end
   else if interval < 0.0 then while tic >= totic-half do begin
      if dosubtics then writeln(output,
                        'Sorry, no subtics with negative scales');

      if theaxis = 'x'
      then xtic(afile,length,dx,dy,tic,width,decimal,lognormal,logbase)
      else ytic(afile,length,dx,dy,tic,width,decimal,lognormal,logbase);

      tic := tic + interval;
      if tic >= totic-half then begin
         if theaxis = 'x' then mover(afile,jump,0.0)
                          else mover(afile,0.0,jump);
         jumpdistance := jumpdistance + jump
      end
   end;

   if theaxis = 'x' then mover(afile,-jumpdistance,0.0)
                    else mover(afile,0.0,-jumpdistance);
   writeln(afile,'grestore');
end;
(* end module pic.doaxis version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.xaxis *)
procedure xaxis(var afile: text;
         doaxisline: boolean;  (* line on axis is plotted *)
         axlength,fromtic,interval,totic: real;
         xsubintervals: real;
         length, dx, dy: real;
         width, decimal: integer;
         logxscale, logxnormal: boolean;
         logxbase: real);
(* draw an x axis starting from the current position. *)
begin
   doaxis(afile,
         'x',
         doaxisline,
         axlength,fromtic,interval,totic,
         xsubintervals,
         length, dx, dy,
         width, decimal,
         logxscale, logxnormal,
         logxbase)
end;
(* end module pic.xaxis version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.yaxis *)
procedure yaxis(var afile: text;
         doaxisline: boolean;  (* line on axis is plotted *)
         aylength,fromtic,interval,totic: real;
         ysubintervals: real;
         length, dx, dy: real;
         width, decimal: integer;
         logyscale, logynormal: boolean;
         logybase: real);
(* draw an y axis starting from the current position. *)
begin
   doaxis(afile,
         'y',
         doaxisline,
         aylength,fromtic,interval,totic,
         ysubintervals,
         length, dx, dy,
         width, decimal,
         logyscale, logynormal,
         logybase)
end;
(* end module pic.yaxis version = 2.80; (@ of dops.p 2010 Nov 03 *)
(* ********************************************************************** *)
(* end module pic.functions version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.cboxr *)
procedure cboxr(var afile: text; width, height: real);
(* make a box to file afile with width in the x direction
and height in the y direction as given.
the box is centered at the current position.
the box is relative to the current position, so it
returns to original position afterwards *)
var h2,w2: real; (* height and width over 2 *)
begin
   h2 := height/2;
   w2 := width/2;
   mover(afile,-w2,-h2);
   liner(afile,0.0,height);
   liner(afile,width,0.0);
   liner(afile,0.0,-height);
   liner(afile,-width,0.0);
   mover(afile,w2,h2);
end;
(* end module pic.cboxr version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.skybox *)
procedure skybox(var afile: text; x,y, width, height: real);
(* make a box to file afile with width in the x direction.
The top center of the box is at x,y and the bottom is at x,0.
height is not used in this version.
*)
var h2,w2: real; (* height and width over 2 *)
begin
   h2 := height/2;
   w2 := width/2;
   movea(afile,x-h2,0.0);
   linea(afile,x+h2,0.0);
   linea(afile,x+h2,y);
   linea(afile,x-h2,y);
   linea(afile,x-h2,0.0);
end;
(* end module pic.skybox *)

(* begin module pic.dotr *)
procedure dotr(var afile: text);
(* draw a dot at the current position *)
begin
   (* a zero length line won't work here, since some systems
      really believe it and draw nothing.  (eg. PostScript on the
      Applewriter draws nothing, but PostScript on the Sun puts one pixle *)
   drawr(afile, dotfactor,0.0,'l',0.0);
end;
(* end module pic.dotr version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.degtorad *)
function degtorad(angle: real):real;
(* convert angle in degrees to radians *)
begin
   degtorad := (angle / 360) * 2 * pi
end;
(* end module pic.degtorad version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.polrec *)
procedure polrec(r,theta: real; var x,y: real);
(* convert polar to rectangular coordinates,
theta is in radians *)
begin
   x := r*cos(theta);
   y := r*sin(theta)
end;
(* end module pic.polrec version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.arc *)
procedure arc(var thefile: text; angle1, angle2, radius: real;
              steps: integer);
(* create an arc in thefile going from angle1 to angle2 (degrees) in the
positive direction of angle, with the given radius.
use the given number of steps to make it.  return to the same position
as before the arc was drawn. *)
var
   dtheta: real; (* change in theta *)
(* s: integer; (@ index to the steps *)
   theta: real; (* current angle *)
   x,y: real; (* coordinates around starting point *)
(* zerox,zeroy: real; (@ starting location, center of curve *)
begin
(*
   zerox := picxglobal;
   zeroy := picyglobal;
*)
   theta := degtorad(angle1);
   dtheta := degtorad( abs(angle2-angle1)/steps );
   polrec(radius,theta,x,y);
(* can't do this for postscript arc:
   movea(thefile,zerox+x,zeroy+y);
      ' ',scale*zerox:  picwidth:picdecim,
      ' ',scale*zeroy:  picwidth:picdecim,
*)

   (* NONSTANDARD postscript, much faster *)
   writeln(thefile,
   (* 'stroke newpath', *)(* force there to be no current point *)
   (* ' 0 0', *)
      'a',
      ' ',scale*radius: picwidth:picdecim,
      ' ',angle1:       picwidth:picdecim,
      ' ',angle2:       picwidth:picdecim);
   write(thefile, 'arc');
   if angle2 < angle1 then write(thefile,'n'); (* for negative draws *)
   (* origin move: writeln(thefile,' stroke newpath 0 0 moveto'); *)
   writeln(thefile,' n');
   (* the moveto puts us back to the origin *)

(*
   for s := 1 to steps do begin
      theta := theta + dtheta;
      polrec(radius,theta, x,y);
      linea(thefile,zerox+x,zeroy+y);
   end;

   movea(thefile,zerox,zeroy)
*)
end;
(* end module pic.arc version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.circler *)
procedure circler(var afile: text; radius: real);
(* make a circle at the current position of some radius. *)
var
   steps: integer; (* number of steps to make the circle *)
begin
   (* number of segments increases with diameter, but the constant
      still should be a function of how good it looks on a particular
      graphic system, I'm afraid.  However, there should be a lower
      bound on the number of steps, so even small circles look good *)
   if radius < 1.0 then steps := 25
                   else steps := round(radius*25);

   arc(afile,0.0,360.0,radius,steps);
end;
(* end module pic.circler version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.ibeam *)
procedure ibeam(var afile: text; width, height: real);
(* Make an ibeam shaped symbol to file afile with width in the x direction
and height in the y direction.  Center it at the current position.
Put a circle at the center, with radius 1/4th the width
Return to original position afterwards. *)
var h2,w2: real; (* height and width over 2 *)
    r: real; (* the radius of the circle *)
begin
   h2 := height/2;
   w2 := width/2;
   mover(afile,-w2,-h2);
   liner(afile,width,0.0);
   mover(afile,-width,height);
   liner(afile,width,0.0);
   mover(afile,-w2,0.0);
   liner(afile,0.0,-height);
   mover(afile,0.0,h2);
   r := width/8;
{ this is silly for the new printers
   if r < 0.025 then r := 0.025; (* small circles do not come out well *)
}
   circler(afile,r);
end;
(* end module pic.ibeam version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.xr *)
procedure xr(var afile: text; width, height: real);
(* make an x to file afile with width in the x direction
and height in the y direction as given.
the box is centered at the current position.
the box is relative to the current position, so it
returns to original position afterwards *)
var h2,w2: real; (* height and width over 2 *)
begin
   h2 := height/2;
   w2 := width/2;
   mover(afile,-w2,-h2);
   liner(afile,width,height);
        mover(afile,0,-height);
   liner(afile,-width,height);
   mover(afile,w2,-h2);
end;
(* end module pic.xr version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.plusr *)
procedure plusr(var afile: text; width, height: real);
(* make a plus sign to file afile with width in the x direction
and height in the y direction as given.
the box is centered at the current position.
the box is relative to the current position, so it
returns to original position afterwards *)
var h2,w2: real; (* height and width over 2 *)
begin
   h2 := height/2;
   w2 := width/2;
   mover(afile,-w2,0);
   liner(afile,width,0.0);
   mover(afile,-w2,h2);
   liner(afile,0.0,-height);
   mover(afile,0,h2);
end;
(* end module pic.plusr version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.rectinit *)
procedure rectinit(var outfile: text);
(* create the definition of a rectangle.  Rectsize must be called
to initialize and to change the size of the rectangle. *)
begin
   writeln(outfile,'/rct');
   writeln(outfile,'{gsave');
   writeln(outfile,'  newpath');
   writeln(outfile,'  0 0 moveto');
   writeln(outfile,'  xs 0 lineto');  (* xs is the x side length *)
   writeln(outfile,'  xs ys lineto'); (* ys is the x side length *)
   writeln(outfile,'  0 ys lineto');
   writeln(outfile,'  closepath fill');
   writeln(outfile,'grestore} def');
end;
(* end module pic.rectinit version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.rectsize *)
procedure rectsize(var afile: text;
                   var xsideold, ysideold, xside, yside: real);
(* determine if the values of xside and yside have changed from xsideold and
ysideold.  If either has changed, write the instructions to change the size of
the rectangle into afile. *)
begin
   if xside <> xsideold then begin
      writeln(afile,'/xs ',
                    (xside*scale):picwidth:picdecim,
                    ' def'); (* xs is the x side length used in rectinit *)
      xsideold := xside;
   end;

   if yside <> ysideold then begin
      writeln(afile,'/ys ',
                    (yside*scale):picwidth:picdecim,
                    ' def'); (* ys is the y side length used in rectinit *)
      ysideold := yside;
   end
end;
(* end module pic.rectsize version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.rectdo *)
procedure rectdo(var afile: text);
(* Make a rectangle with lower right hand corner at (x,y) and the
given side (position and color are determined by earlier calls
to setcolor and move).  Return to original position afterwards. *)
begin
   writeln(afile, ' rct');
end;
(* end module pic.rectdo version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* begin module pic.setgray *)
procedure setgray(var afile: text;
                      brightness: real);
(* set the gray scale to the requested one.
Range of the variables is 0 to 1. *)
const
   colfield = 8; (* width of numbers printed to the file *)
   colwidth = 4; (* number of decimal places for numbers *)
   (* PostScript on a Sun 4 cannot handle 5 decimal places!  Use 4 or less *)
begin
   write  (afile,'n'); (* be sure it's started cleanly *)
   write  (afile,' ',brightness:colfield:colwidth);
   writeln(afile,' setgray');
end;
(* end module pic.setgray version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* LOCK begin module pic.setcolor *)
procedure setcolor(var afile: text; colorkind: char;
                       hue, saturation, brightness: real);
(* Set the color to the requested one.  Range of the variables is 0 to 1.
   Colors in PostScript are defined with hue, saturation and brightness with
   the sethsbcolor function.  The standard hue runs from red at 0 to red
   at 1 with the Roy G. Biv in between.  The famous physicist Roy G. Biv
   stands for: red, orange, yellow, green, blue, indigo, violet, the colors
   of the spectrum.  On 1992 September 16 I realized that, amazingly, the
   wavelength in nm match quite nicely!

      700 red
      750 orange
      600 yellow
      650 green
      500 blue
      550 indego
      400 violet

   Since the hue runs in a circle, the spectrum is not generated from
   the hue range 0 to 1.  The routine setcolor therefore converts the
   input numbers by multiplying by 0.84 and adding 0.16.  This gives
   the color range from yellow through red corresponding to values 0
   to 1.  (Note: adding 0.16 give the range from yellow through red,
   but that is not a spectrum which runs from red through violet.)

   If colorkind is 'h' use hsb.
   If colorkind is 'r' use rgb.

*)
const
   colfield = 8; (* width of numbers printed to the file *)
   colwidth = 4; (* number of decimal places for numbers *)
   (* PostScript on a Sun 4 cannot handle 5 decimal places!  Use 4 or less *)
   huefactor = 0.84; (* number of decimal places for numbers *)
begin
   write  (afile,'n'); (* be sure it's started cleanly *)
{
   write  (afile,' ',hue*huefactor:       colfield:colwidth);
}
   if colorkind = 'h' then begin
      write  (afile,' ',hue*huefactor + (1-huefactor): colfield:colwidth);
      write  (afile,' ',saturation:colfield:colwidth);
      write  (afile,' ',brightness:colfield:colwidth);
      writeln(afile,' sethsbcolor');
   end
   else begin
      write  (afile,' ',hue: colfield:colwidth);  (* red! *)
      write  (afile,' ',saturation:colfield:colwidth);  (* green! *)
      write  (afile,' ',brightness:colfield:colwidth);  (* blue! *)
      writeln(afile,' setrgbcolor');
   end;
end;
(* LOCK end module pic.setcolor version = 2.73; (@ of dops.p 2003 Aug 28 *)

(* begin module pic.boxintercept *)
procedure boxintercept(xmin,ymin,xmax,ymax,m,b: real;
                       var intercept: boolean; var x1,y1,x2,y2: real);
(* does the line y=m*x+b intercept the box defined by the points
(xmin,ymin) and (xmax,ymax)?  if so, intercept is true and
the intercept points are given by (x1,y1) and (x2,y2) *)
var
   count: integer; (* number of sides the line intersects *)
      xlo,xhi,ylo,yhi: boolean; (* whether the line intersects the
         box at the low value of x, etc *)
function fny(x: real):real;
(* calculate the y value given the x *)
begin fny := m*x+b end;

function fnx(y: real):real;
(* calculate the x value given the y *)
begin fnx := (y-b)/m end;

function between(a,b,c: real):boolean;
(* is b between a and c? *)
begin between:=(a<=b) and (b<=c) end;

procedure normalcases;
(* analyze for the usual cases when the slope m is not zero *)
begin (* normalcases *)
(* writeln(output,'m=',m:5:4,' b=',b:5:4,' in normalcases');*)
(*  ymin, ymax, xmin and xmax are the coordinates of the box;
    xlo, yhi, xhi, and ylo are which sides are intersected

                             |
     ymax +----yhi----+      |
          |           |      |
         xlo         xhi     |
          |           |      |
     ymin +----ylo----+      |
        xmin        xmax     |
                             |

*)
      xlo := between(ymin,fny(xmin),ymax);
      xhi := between(ymin,fny(xmax),ymax);
      ylo := between(xmin,fnx(ymin),xmax);
      yhi := between(xmin,fnx(ymax),xmax);
(*
writeln(output,'xlo = ',xlo);
writeln(output,'xhi = ',xhi);
writeln(output,'ylo = ',ylo);
writeln(output,'yhi = ',yhi);
*)

      intercept := true; (* optimistic *)

      (* simplify cases which intersect corners.  These are the ones
         where more than two side intersections are true. *)
      count := 0;
      if xlo then count := succ(count);
      if xhi then count := succ(count);
      if ylo then count := succ(count);
      if yhi then count := succ(count);

      if count > 2 then begin
	 (* one corner must be intersected.  Simplify by preserving the
	    opposing intersections.  When there are 4 intersections, only one
	    pair is preserved.  The calculations will be correct either way. *)
         if xlo and xhi then begin
            yhi := false;
            ylo := false;
         end
         else if  ylo and yhi then begin
            xhi := false;
            xlo := false;
         end
         else begin
            writeln(output,'error in between count!');
            halt
         end
      end;

           if xlo and xhi then begin x1 := xmin;      x2 := xmax        end
      else if xlo and ylo then begin x1 := xmin;      x2 := fnx(ymin)   end
      else if xlo and yhi then begin x1 := xmin;      x2 := fnx(ymax)   end
      else if xhi and ylo then begin x1 := xmax;      x2 := fnx(ymin)   end
      else if xhi and yhi then begin x1 := xmax;      x2 := fnx(ymax)   end
      else if ylo and yhi then begin x1 := fnx(ymin); x2 := fnx(ymax)   end
      else intercept := false;

      if intercept then begin
         y1 := fny(x1);
         y2 := fny(x2)
      end

end; (* normalcases *)
begin (* boxintercept *)
      (* note: abs(m) is required to protect against negative zero... *)
      if abs(m) = 0.0
      then begin
         intercept := between(ymin,b,ymax);
         if intercept then begin
            x1 := xmin; y1 := b;
            x2 := xmax; y2 := b;
         end
      end
      else normalcases
end; (* boxintercept *)
(* end module pic.boxintercept version = 2.80; (@ of dops.p 2010 Nov 03 *)

(* ************************************************************************ *)
(* ************************************************************************ *)

(* begin module checknumber *)
function checknumber(var afile: text): boolean;
(* check that there is a number next in the file.  If not, return
false.   This is useful for protection when reading a parameter file. *)
var ok: boolean; (* result of this check *)
procedure conclude;
begin
   writeln(output,'Including this character, the rest of the data line is:');
   copyaline(afile,output);
   ok := false;
end;
begin
   ok := true; (* be optimistic *)
   if eof(afile)
   then begin
      ok := false;
      write  (output,'A number was expected on a data line, but');
      writeln(output,' the end of the file was found instead.');
   end
   else begin
      skipblanks(afile);
      if eoln(afile) then begin
         write  (output,'A number was expected on a data line, but');
         writeln(output,' the end of the line was found instead.');
         conclude;
      end;
      if not (afile^ in
         ['0','1','2','3','4','5','6','7','8','9','.','-','+']) then begin
         write  (output,'A number was expected on a data line, but');
         writeln(output,' the character "',afile^,'" was found instead.');
         conclude;
      end;
   end;
   checknumber := ok
end;
(* end module checknumber version = 5.22; (@ of prgmod.p 2005 Jul 12 *)

(* begin module xyplo.tellwarning *)
procedure tellwarning;
(* tell the user that there are warnings.  Unfortunately
these are globals: warnings: text, warningcount: integer. *)
begin
   warningcount := succ(warningcount);
   if warningcount = 1 then begin
      rewrite(warnings);
      writeln(warnings,'xyplo ',version:4:2);
   end;
   writeln(warnings);
   write(warnings,warningcount:1,' ');
end;
(* end module xyplo.tellwarning *)

(* begin module marksdo *)
procedure marksdo(var afile: text; width, height: real; flagstring: string);
(* make the postscript mark defined by the user *)
begin
   if flagstring.length = 0
   then begin
      writeln(output,'WARNING: user defined mark not used in xyin');
   end
   else begin
      writeln(afile);
      writeln(afile,'gsave % marksdo');
      write(afile,    scale*width:picwidth:picdecim,
                  ' ',scale*height: picwidth:picdecim,' ');
      writestring(afile,flagstring);
      writeln(afile);
      writeln(afile,'grestore % marksdo');
      writeln(afile);
   end
end;
(* end module marksdo *)

(* begin module xyplo.gettoken *)
procedure gettoken(var infile: text; var flag: char; var tokenstring: string;
                   var gots: boolean);
(* get the tokenstring from infile.  the first character of the string
is placed in flag.  gots is true if a non-blank token was found.
2011 Mar 03: tabs are also breaks! *)
var done: boolean; (* done reading characters *)
begin
   clearstring(tokenstring);
   done := false;
   while not done do with tokenstring do begin
      if eoln(infile)
      then done := true
      else if (infile^=' ') or (infile^=chr(9)) (* parse on space or tab *)
           then done := true
           else begin
              length := length + 1;
              read(infile,letters[length]);
           end
   end;
   gots := tokenstring.length <> 0;
   flag := tokenstring.letters[1];
end;
(* end module xyplo.gettoken *)

(* begin module xyplo.equalstring *)
function equalstring(a,b: string): boolean;
(* are the two strings equal up to their lengths? *)
var
   notequal: boolean; (* true if the strings are not equal *)
   j: integer; (* index to positions in the strings *)
begin
   if a.length <> b.length
   then equalstring := false
   else begin
     j:=1;
     notequal := false; (* until proven guilty *)

     repeat
        notequal := (a.letters[j] <> b.letters[j]);
        j := j + 1;
     until (j > a.length) or notequal;

     equalstring := (not notequal) and (j > a.length)
   end
end;
(* end module xyplo.equalstring *)

(* begin module xyplo.bar *)
procedure bar(var f: text; c: char);
(* produce a bar of the character c to separate parts of the file f *)
var i: integer; (* index *)
begin
   for i:=1 to 50 do write(f, c);
   writeln(f)
end;
(* end module xyplo.bar *)

(* begin module xyplo.readsymbol *)
procedure readsymbol(var f: text; symb: char; var p: xyploparam; var s: symbol);
(* read the symbol definition from f *)
var flag: char; (* first character of symbolflag string *)
    gotsymbol: boolean; (* the symbol string was read ok *)
procedure t;
(* test that the file is still there to read *)
begin
      if eof(f) then begin
         writeln(output,'missing symbol parameters');
         halt
      end
end;
begin
   s.symboltoplot := symb;
   with s,s.l do begin
      (* what to recognize for this symbol *)
      t; gettoken(f,flag,symbolflag,gotsymbol);
      if not gotsymbol then begin
         writeln(output,'Could not find symbol flag for symbol "',symb,'"');
         writeln(output,'while reading symbols in the xyplop file.');
         writeln(output,'Trying to read a symbol from this token:');
         writestring(output,symbolflag);
         writeln(output);
         writeln(output,'The rest of the file contains: ');
         while not eof(f) do copyaline(f,output);
         halt
      end;
      readln(f);

      t; readln(f,symbolxsize); (* its size on the x axis *)
      if (symbolxsize<=0) and (not p.needxscolumn) then begin
         writeln(output,'While reading symbol "',symb,
                        '" with flag "',flag,'" you asked that its');
         writeln(output,'x size be determined from a column.');
         writeln(output,'unfortunately you did not specify a positive column!');
         halt
      end;
(*writeln(output,'in readsymbol: needxscolumn = ',p.needxscolumn);*)
(*writeln(output,'symbolxsize = ',symbolxsize:10:2);*)

      t; readln(f,symbolysize); (* its size on the y axis *)
      if (symbolysize<=0) and (not p.needyscolumn) then begin
         writeln(output,'While reading symbol "',symb,
                        '" with flag "',flag,'" you asked that its');
         writeln(output,'y size be determined from a column.');
         writeln(output,'unfortunately you did not specify a positive column!');
         halt
      end;
(*writeln(output,'in readsymbol: needyscolumn = ',p.needyscolumn);*)
(*writeln(output,'symbolysize = ',symbolysize:10:2);*)

      if (symbolxsize = 0.0) or (symbolysize = 0.0) then begin
         writeln(output,'symbol size cannot be zero');
         halt
      end;

      (* decide on connection lines *)
      if f^='c' then begin
         get(f);
         if (f^ = 'n') or (f^='i')
         then doconnection := false (* oh well *)
         else begin
            doconnection := true;
            read(f,connecttype);
            if (connecttype = '.') or (connecttype = '-') then begin
               skipnonblanks(f);
               read(f,connectsize);
            end
            else connectsize := 0.05;
            didlastpoint := false;
         end
      end
      else doconnection := false;
      readln(f);

      (* read in the symbol line definition *)
      if not eoln(f) then begin
         read(f,linetype);
         if (linetype='l') or (linetype='.') or
            (linetype='-') or (linetype='i')
         then doline := true
         else if linetype='n' then doline := false
         else begin
            writeln(output,'for symbol ',symboltoplot,
                           ' line type must be one of "lin.-"');
               halt
         end;
         if (linetype='-') or (linetype='.') then begin
            t; skipnonblanks(f); read(f,linesize)
         end
         else linesize := 0.05
      end;
      readln(f);
      symbolhasxyplomdef := false;
   end;
   if (symb = 'r') or (symb = 'R')
   then p.startrectangles := true; (* start rectangles later *)
end;
(* end module xyplo.readsymbol *)

(* begin module copyfile *)
procedure copyfile(var fin, fout: text);
(* copy the rest of file fin to fout *)
begin
   while not eof(fin) do copyaline(fin, fout);
end;
(* end module copyfile version = 5.22; (@ of prgmod.p 2005 Jul 12 *)

(* begin module xyplo.upgrade1 *)
procedure upgrade1(var xyplop: text; p: xyploparam);
(* upgrade the parameter file xyplop *)
var
   internal: text; (* a place to hold the old xyplop *)
begin with p do begin
   reset(xyplop);
   rewrite(internal);

   (* copy xyplop to internal *)
   copyfile(xyplop, internal);

   (* copy internal to xyplop *)
   reset(internal);
   rewrite(xyplop);
   copyfile(internal, xyplop);

   (* add the new material to the end: *)
   writeln(xyplop,'          **** more parameters');
   write  (xyplop,'p ',
                  edgeleft:4:2,' ',edgeright:4:2,' ',edgelow:4:2,' ',edgehigh:4:2,
                  ' edgecontrol (p=page),');
   writeln(xyplop,' edgeleft, edgeright, edgelow, edgehigh in cm');
   writeln(xyplop, version:4:2,' version of',
         ' xyplo that this parameter file is designed for.');

   reset(xyplop); (* ready to start reading again *)
end end;
(* end module xyplo.upgrade1 *)

(* begin module xyplo.readparam *)
procedure readparam(var f: text; var p: xyploparam);
(* read the parameters from f *)
var
   checkout: boolean; (* if true, all variable values are ok *)
   moreparameters: boolean; (* more parameters will be following because
      a blank character is ending the line definitions *)
   gotten: boolean; (* for finding the label strings *)
   symb: char; (* the symbol to plot or '.' to end reading symbols *)
   symbolcount: integer; (* the number of symbols read so far,
      for informing the user about errors *)
   ul: linesptr; (* a pointer for constructing the user lines *)
   us: symbolsptr; (* a pointer for constructing the user symbols *)

procedure helpem;
(* help the poor user figure out where the problem is *)
begin
   writeln(output);
   writeln(output,'To help you locate the problem, here''s the rest');
   writeln(output,'of the xyplop file:');
   writeln(output,'"""""""""""""""""""""""""""""""""""""""""""""""""""""""');
   while not eof(xyplop) do copyaline(xyplop,output);
   writeln(output,'"""""""""""""""""""""""""""""""""""""""""""""""""""""""');
   halt
end;
procedure cn;
(* short version of call to check number *)
begin
  checkout := checknumber(f);
  if not checkout then helpem; (* avoid snowballing *)
end;
procedure t;
(* test that the file is still there to read *)
begin
      if eof(f) then begin
         writeln(output,'Missing parameters: unexpected end of xyplop');
         halt
      end
end;

procedure num(var f: text);
(* test that a number follows *)
begin
   skipblanks(f);
   if not (f^ in ['-','+','0','1','2','3','4','5','6','7','8','9']) then begin
      writeln(output,'found this character: "',f^,'" when expecting a number',
                     ' in the parameter file');
      helpem
   end
end;
function softnumbertest(var f: text): boolean;
(* test that a number follows, but don't die if there is none *)
begin
   skipblanks(f);
   if not (f^ in ['-','+','0','1','2','3','4','5','6','7','8','9']) then begin
      bar(output,'-');
      writeln(output,'NOTE: in xyplop, another number is now allowed for',
        ' this line:');
      copyaline(f,output);
      softnumbertest := false
   end
   else softnumbertest := true
end;

procedure testcolumns(a,b,c,d,e,f,g,h: integer);
(* test that the 8 columns do not conflict spatially *)
const
   columns = 8; (* number of columns.  Must correspond with the number of
      arguments to testcolumns *)
var
   all: array[1..8] of integer; (* the locations of the columns *)
   die: boolean; (* program will stop if any errors are found *)
   x,y: integer; (* index to all *)
function notok(a,b:integer):boolean;
(* check for the locations of two columns: if they are equal
and and positive, then the grabdata routine will not be able
to read them both *)
begin (* notok *)
   if (a=b) and (a>0) then notok:=true
                      else notok:=false
end; (* notok *)
begin (* testcolumns *)
   all[1] := a; all[2] := b; all[3] := c; all[4] := d;
   all[5] := e; all[6] := f; all[7] := g; all[8] := h;
   die := false ;

   for x := 1 to columns-1 do begin
      for y := x+1 to columns do begin
         if notok(all[x],all[y])
         then begin
            if not die then writeln(output,'no two columns can have the',
                                           ' same positive position');
            writeln(output,'columns in conflict are numbered ', all[x]:1);
            die := true
         end
      end
   end;
   if die then helpem
end; (* testcolumns *)
procedure nextsection(id: char);
(* make sure that there is a dashed line between sections,
id is the name of the section:
  c columns
  s symbols
  l lines
  m more parameters
  *)
procedure nosection;
(* object if we did not find the next section *)
begin (* nosection *)
   writeln(output,'ERROR: SECTION DIVIDER NOT FOUND');
   writeln(output,'Section divider lines in xyplop may (but need not)',
                  ' begin with blanks');
   writeln(output,'and these must be followed by at least one asterisk ("*").');
   writeln(output,'Check that your xyplop matches the current documentation.');
   writeln(output,'The error was found while');
   write(output,'Xyplo was looking for the section divider before ');
   case id of
      'c': writeln(output,'the COLUMN section.');
      's': writeln(output,'SYMBOL number ',symbolcount:1,'.');
      'l': writeln(output,'the LINE section.');
      'm': writeln(output,'the MORE PARAMETERS section.');
   end;
   helpem;

end; (* nosection *)
begin (* nextsection *)
   skipblanks(f);
   if eoln(f) then nosection;
   if f^<>'*' then nosection;
   readln(f);
end; (* nextsection *)

begin (* readparam *)
      checkout := true; (* be optimistic *)
      reset(f);
      with p do begin
         (* read in graph shift amounts *)
         t; num(f); read(f,xzero); num(f); readln(f,yzero);

         (* read in max and min for x *)
         setx := false;
         t; if not eoln(f)
         then begin
            if f^='x' then begin
               setx := true;
               skipnonblanks(f);
               t; num(f); read(f,xmin); num(f); read(f,xmax);
               if xmin >= xmax then begin
                  writeln(output,'xmin (',ymin:1:1,') cannot be >= ',
                                 'xmax (',xmax:1:1,') in xyplop');
                  halt
               end;
            end
         end;
         t; readln(f);
         if not setx then begin
            xmin := maxint; xmax := -maxint;
         end;

         (* read in max and min for y *)
         sety := false;
         t; if not eoln(f)
         then begin
            if f^='y' then begin
               sety := true;
               skipnonblanks(f);
               t; num(f); read(f,ymin); num(f); read(f,ymax);
               if ymin >= ymax then begin
                  writeln(output,'ymin (',ymin:1:1,') cannot be >= ',
                                 'ymax (',ymax:1:1,') in xyplop');
                  halt
               end;
            end
         end;
         t; readln(f);
         if not sety then begin
            ymin := maxint; ymax := -maxint;
         end;

         (* read in the xinterval and yinterval, and stay on the line *)
         t;
         num(f); read(f,xinterval); num(f); read(f,yinterval);
         if (xinterval<=0) or (yinterval<=0) then begin
            writeln(output,'x and y interval must be positive');
            halt
         end;

         if (xinterval>=maxinterval) or (yinterval>=maxinterval) then begin
            writeln(output,'x and y interval must be less than maxinterval',
                           ' (=',maxinterval:1,')');
            writeln(output,'You really can''t read a graph with that',
                           ' many intervals!!');
            halt
         end;

         (* continue the same line to read in the xintervals and yintervals *)

         if softnumbertest(f) then begin
            read(f, xsubintervals);
            if xsubintervals<=0 then begin
               writeln(output,'xsubintervals must be positive');
               halt
            end;

            if softnumbertest(f) then begin
               read(f, ysubintervals);
               if ysubintervals<=0 then begin
                  writeln(output,'ysubintervals must be positive');
                  halt
               end;
               readln(f);
            end
            else begin
               writeln(output,'You may now define the number of',
                              ' y subtic marks');
               bar(output,'-');
            end;
         end
         else begin
            writeln(output,'You may now define the number of',
                           ' x and y subtic marks');
            bar(output,'-');
         end;

         if (xsubintervals>=maxinterval) or (ysubintervals>=maxinterval) then begin
            writeln(output,'x and y subintervals must be less than maxinterval',
                           ' (=',maxinterval:1,')');
            writeln(output,'You really can''t read a graph with that',
                           ' many intervals!!');
            halt
         end;

         (* read in the xwidth and ywidth *)
         t; num(f); read(f,xwidth); num(f); readln(f,ywidth);
         if (xwidth<=0) or (ywidth<=0) then begin
            writeln(output,'x and y number widths must be positive');
            halt
         end;

         (* read in the xdecimal and ydecimal *)
         t; num(f); read(f,xdecimal); num(f); readln(f,ydecimal);
         if (xdecimal<0) or (ydecimal<0) then begin
            writeln(output,'x and y decimal places must be zero or positive');
            halt
         end;

         (* check that the x widths and decimal places are sensible *)
         if (xdecimal > 0) and (xwidth < xdecimal + 2) then begin
            writeln(output,'The width of numbers on the axis (xwidth) must',
                           ' be at least two larger');
            writeln(output,'than the number of decimal places (xdecimal),',
                           ' to provide for the sign and');
            writeln(output,'decimal place.');
            halt
         end;

         (* check that the y widths and decimal places are sensible *)
         if (ydecimal > 0) and (ywidth < ydecimal + 2) then begin
            writeln(output,'The width of numbers on the axis (ywidth) must',
                           ' be at least two larger');
            writeln(output,'than the number of decimal places (ydecimal),',
                           ' to provide for the sign and');
            writeln(output,'decimal place.');
            halt
         end;

         (* read in the xsize and ysize *)
         t; num(f); read(f,xsize); num(f); readln(f,ysize);
         if (xsize<0) or (ysize<=0) then begin
            writeln(output,'x and y size must be positive');
            halt
         end;

         (* get the labels for the axes *)
         getstring(f,xlabel,gotten);
         if not gotten then begin
            writeln(output,'cannot find xlabel');
            halt
         end;
         if f^ = '^'
         then begin
            rotateylabel := false;
            get(f); (* skip the '^' *)
         end
         else rotateylabel := true;
         getstring(f,ylabel,gotten);
         if not gotten then begin
            writeln(output,'cannot find ylabel');
            halt
         end;

         (* define cross hairs, and set which axis to plot *)
         crosshairs := false;
         t; if not eoln(f)
         then begin
            if (f^='c') or (f^='X') or (f^='Y') or (f^='N')
            then crosshairs := true;

            (* the default is true *)
            doXaxis := true;
            doYaxis := true;

            (* turn off the other axis *)
            if (f^='x') or (f^='X') then doYaxis := false;
            if (f^='y') or (f^='Y') then doXaxis := false;
            if (f^='n') or (f^='N') then begin
               doXaxis := false;
               doYaxis := false;
            end;
            if (f^='i')
            then doaxisline := false
            else doaxisline := true;
         end;
         readln(f);

         (* make type of x axis *)
         logxscale := false;
         logxnormal := false;
         t; if not eoln(f)
         then if (f^='l') or (f^='L') then begin
            logxscale := true;
            if f^='L' then logxnormal := true;
            skipnonblanks(f);
            num(f); read(f,xbase);
            if xbase <= 1.0 then begin
               writeln(output,'log x base must be > 1');
               halt
            end;
            logxbase := ln(xbase);
         end;
         readln(f);

         (* make type of y axis *)
         logyscale := false;
         logynormal := false;
         t; if not eoln(f)
         then if (f^='l') or (f^='L') then begin
            logyscale := true;
            if f^='L' then logynormal := true;
            skipnonblanks(f);
            num(f); read(f,ybase);
            if ybase <= 1.0 then begin
               writeln(output,'log y base must be > 1');
               halt
            end;
            logybase := ln(ybase);
         end;
         readln(f);

         (* skip the line between sections of the file *)
         nextsection('c');

         (* define the columns to read data from *)
         (* read the symbol columns *)
         t; num(f); read(f,xcolumn);
            num(f); read(f,ycolumn);
         readln(f);

         if (xcolumn<=0) or (ycolumn<=0) then begin
            writeln(output, ' in xyplop, the defined locations of',
                            ' both x and y columns must be positive');
            halt
         end;

         t; num(f); read(f,scolumn);
         readln(f);

         t; num(f); read(f,xscolumn);

            num(f); read(f,yscolumn);
         readln(f);

{qqq}
         t; if (f^ in [' ','0','1','2','3','4','5','6','7','8','9','h'])
         then colorkind := 'h'
         else
            if (f^ = 'r')
            then colorkind := 'r'
            else begin
               writeln(output,'colorkind must be: blank, digit, h or r');
               halt;
            end;
         t; if (f^ in ['h','r']) then get(f);
         t; num(f); read(f,hucolumn); (* hue column *)
         t; num(f); read(f,sacolumn); (* saturation column *)
         t; num(f); read(f,brcolumn); (* brightness column *)
         readln(f);

         testcolumns(xcolumn,ycolumn,scolumn,xscolumn,yscolumn,
                     hucolumn,sacolumn,brcolumn);

         needscolumn := (scolumn>0);
         needxscolumn := (xscolumn>0);
         needyscolumn := (yscolumn>0);

         needcocolumn := (hucolumn>0) or (sacolumn>0) or (brcolumn>0);

         (* Make sure that the rectangles have a size that is not likely
           to be requested, so that the proper size will be set on the
           first attempt to make a rectangle. *)
         xrect := -maxint;
         yrect := -maxint;

         symbolcount:=1;
         (* skip the line between sections of the file *)
         nextsection('s');

         (* define all of the symbols *)
         t; readln(f,symb);
         if symb = '.' then begin
            writeln(output,'you must define at least one symbol');
            halt;
         end;
         new(usersymbols);
         us := usersymbols;
         startrectangles := false;
         repeat
(*writeln(output,'reading symbols');*)
            readsymbol(f,symb,p,us^.s);

            symbolcount:=symbolcount+1;
            (* skip the line between sections of the file *)
            nextsection('s');

            t; readln(f,symb);
            if symb <> '.' then begin (* make the next one *)
               new(us^.next);
               us := us^.next
            end
            else us^.next := nil;
         until symb = '.';

         (* skip the line between sections of the file *)
         nextsection('l');

(*writeln(output,'user defined lines');*)
         moreparameters := false;
         (* read in user defined lines *)
         userlines := nil; (* set nil no matter what *)
         if not eof(f)
         then if f^ = ' ' then begin
            moreparameters := true;
         end
         else begin
(*writeln(output,'inside user defined lines');*)

            new(userlines);
            ul := userlines;

            while (not eof(f) and (not moreparameters)) do with ul^.l do begin
               if eoln(f) then readln(f)
               else if f^='*' then readln(f)
               else if f^=' ' then begin
                       moreparameters := true;
                    end
               else begin
                  t; read(f,linetype);
                  if not (linetype in ['l','.','-','i','n']) then begin
                     writeln(output,'user defined line type',
                                    ' must be one of "lin.-"');
                     writeln(output,'you had "',linetype,'"');
                     halt
                  end;
                  t; num(f); read(f,m);
                  t; num(f); read(f,b);
                  if (linetype='-') or (linetype='.') then begin
                     t; num(f); read(f,linesize);
                  end
                  else linesize := 0;
                  readln(f);

                  if (not eof(f)) and (f^<>' ') then begin
                     if not eoln(f) then begin
                        if f^<>'*' then begin
                           new(ul^.next);
                           ul := ul^.next
                        end
                     end
                  end
                  else ul^.next := nil
               end
            end
         end;

         (* read additional parameters ("more") *)

         (* first, set defaults *)
         llx := defaultllx;
         ury := defaultury;
         urx := defaulturx;
         lly := defaultlly;
         edgecontrol := 'p';
         edgeleft    := 1.5;
         edgeright   := 0.5;
         edgelow     := 1.5;
         edgehigh    := 1.0;
         parameterversion := 0.0;

         if moreparameters then begin
            (* skip the line between sections of the file *)
            nextsection('m');

            read(f, edgecontrol);
            if (edgecontrol = 'p')
            then begin
               cn; read(f, edgeleft);
               cn; read(f, edgeright);
               cn; read(f, edgelow);
               cn; read(f, edgehigh);

{
writeln(output,'xsize',xsize:11:5);
writeln(output,'ysize',ysize:11:5);
writeln(output,'defscale',defscale:11:5);
}

               llx := round((xzero         - edgeleft )* defscale);
               lly := round((yzero         - edgelow  )* defscale);
               urx := round((xzero + xsize + edgeright)* defscale);
               ury := round((yzero + ysize + edgehigh )* defscale);
            end;
            readln(f);

            cn; readln(f, parameterversion);
         end
         else begin
             writeln(output,'*******************************************************');
             writeln(output,'* more parameters are now available, see              *');
             writeln(output,'* ',
                 'https://alum.mit.edu/www/toms/delila/xyplo.html *');
             writeln(output,'*******************************************************');
             writeln(output,' They are being added to the parameter file ***********');

             upgrade1(xyplop, p);

         end

(*;writeln(output,'end of readparam');*)
      end;
end;
(* end module xyplo.readparam *)

(* begin module xyplo.loglabel *)
procedure loglabel(var f: text; s: string; base: real; justification: char);
(* graph the string s to file f as the log to the base given in base *)
var
   i: integer; (* index *)
   n: string; (* to modify the string *)
begin
(*
   write(f,'log\d',trunc(base):1'\u');
   writestring(f,s);
   write(f,')');
*)
   clearstring(n);
   with n do begin
      letters[1] := 'l';
      letters[2] := 'o';
      letters[3] := 'g';
      length := 3;
      if base < 10 then stringinteger(trunc(base),n,1,false)
                   else stringinteger(trunc(base),n,2,false);
      length := length + 1;
      letters[length] := '{';
      for i := 1 to s.length do n.letters[i+n.length] := s.letters[i];
      length := length + s.length + 1; (* for the above and the below *)
      letters[length] := '}';
      graphstring(f,n,justification);
   end
end;
(* end module xyplo.loglabel *)

(* begin module xyplo.logstring *)
procedure logstring(var f: text; s: string; base: real);
(* write the string s to file f as the log to the base given in base *)
begin
   write(f,'log',trunc(base):1,'(');
   writestring(f,s);
   write(f,')');
end;
(* end module xyplo.logstring *)

(* begin module xyplo.comment *)
procedure comment(var f: text);
(* put a PostScript comment start out to file f *)
begin write(f,'% ') end;
(* end module xyplo.comment *)

(* begin module xyplo.writeparam *)
procedure writeparam(var f: text; p: xyploparam);
(* write all the values in p out to file f *)
var
      us: symbolsptr; (* a pointer to the user symbols *)
      ul: linesptr; (* a pointer to the user lines *)
begin
with p do begin
(* writeln(f,'.sp 1'); (@ give room on output @) *)
comment(f);
writeln(f,'user specified parameters:');
comment(f);
bar(f,'*'); (* ************************** *)
comment(f);
bar(f,'*'); (* ************************** *)
comment(f);
if setx then write(f,'x')
        else write(f,'z');
if logxscale and (not logxnormal)
then write(f,' ',exp(logxbase*xmin):xwidth:(xdecimal+1), '  ',
        exp(logxbase*xmax):xwidth+dsafety:(xdecimal+1+dsafety),
        ' minimum and maximum for x axis (log scale)')
else begin
     write(f,' ',xmin:xwidth:(xdecimal+1), '  ',
        xmax:xwidth+dsafety:(xdecimal+1+dsafety),
        ' minimum and maximum for x axis');
     if logxnormal then write(f,' (log scale)');
end;
if setx then writeln(f,' set by user')
        else writeln(f,' determined from data');

comment(f);
if sety then write(f,'y')
        else write(f,'z');
if logyscale and (not logynormal)
then write(f,' ',exp(logybase*ymin):ywidth:(ydecimal+1), '  ',
        exp(logybase*ymax):ywidth+dsafety:(ydecimal+1+dsafety),
        ' minimum and maximum for y axis (log scale) ')
else begin
     write(f,' ',ymin:ywidth:(ydecimal+1), '  ',
        ymax:ywidth+dsafety:(ydecimal+1+dsafety),
        ' minimum and maximum for y axis ');
     if logynormal then write(f,'(log scale)');
end;
if sety then writeln(f,' set by user')
        else writeln(f,' determined from data');

comment(f);
writeln(f,xinterval:10,yinterval:10,' number of intervals on x and y to plot');

comment(f);
writeln(f,xwidth:10,   ywidth:10,' width of numbers on graph in characters');

comment(f);
writeln(f,xdecimal:10, ydecimal:10,  ' number of decimal places for numbers');

comment(f);
writeln(f,xsize:10:3,    ysize:10:3, ' size of axes in cm');

comment(f);
writestring(f,xlabel);writeln(f,'    the x axis label');

comment(f);
writestring(f,ylabel);writeln(f,'    the y axis label');

comment(f);
if not crosshairs then write(f,'no ');
writeln(f,'cross hairs put on zero of x and y');

comment(f);
if logxscale then begin
   if logxnormal then write(f,'L')
                 else write(f,'l');
   write(f,' ',xbase:5:3, ' log scale on x axis');
   if logxnormal then writeln(f,' (normal numbers on scale)')
                 else writeln(f,' (log of numbers on scale)')
end
else writeln(f,'no log scale on x axis');

comment(f);
if logyscale then begin
   if logynormal then write(f,'L')
                 else write(f,'l');
   write(f,' ',ybase:5:3, ' log scale on y axis');
   if logynormal then writeln(f,' (normal numbers on scale)')
                 else writeln(f,' (log of numbers on scale)')
end
else writeln(f,'no log scale on y axis');

(*
comment(f);
if logyscale then writeln(f,'l ',ybase:5:3, ' log scale on y axis')
             else writeln(f,'no log scale on y axis');
*)

comment(f);
write(f,'* column choices: ');
bar(f,'*'); (* ************************** *)
comment(f);
   writeln(f,xcolumn:5,  ycolumn:5,' columns of input chosen',
                                     ' for x and y respectively');
comment(f);
   writeln(f,scolumn:5,  ' ':5,    ' column that determines symbols');
comment(f);
   writeln(f,xscolumn:5, yscolumn:5,' columns that determine symbol size');

comment(f);
   writeln(f,hucolumn:5,sacolumn:5,brcolumn:5,' columns that determine color');

comment(f);
write(f,'* user defined symbols: ');
bar(f,'*'); (* ************************** *)

us := usersymbols;
while us<>nil do with us^.s,us^.s.l do begin
comment(f);
   writeln(f,symboltoplot,' ':12,' symbol-to-plot:',
                          'c(circle)bd(dotted box)x+Ifgpr(rectangle)');
comment(f);
   writestring(f,symbolflag);
   writeln(f,' ':12,' character string in xyin to indicate this symbol');
comment(f);
   writeln(f,symbolxsize:1:3,' ':8,' symbol x size in cm');
comment(f);
   writeln(f,symbolysize:1:3,' ':8,' symbol y size in cm');

comment(f);
   if doconnection
   then writeln(f, 'connected ',connecttype, ' ', connectsize:5:3,
                   ' the symbols are connected by lines')
   else writeln(f,'no connected lines');

comment(f);
   write(f,' ',linetype); (* the space protects the character from troff *)
   if (linetype='-') or (linetype='.') then write(f,' ',linesize:5:3);
   write(f,' l=line .=dotted -=dashed i=invisible n=no line');
   if (linetype='-') or (linetype='.') then write(f,'; spacing in cm');
   writeln(f);
comment(f);
bar(f,'*'); (* ************************** *)
   us := us^.next
end;

comment(f);
writeln(f,'.'); (* a period ends the section *)
comment(f);
write(f,'* user defined lines: ');
bar(f,'*'); (* ************************** *)

ul:=userlines;
while ul<>nil do with ul^.l do begin
comment(f);
   write(f,' ',linetype); (* the space protects the character from troff *)
   write(f,'  ',m:xwidth+dsafety:(xdecimal+1+dsafety));
   write(f,'  ',b:xwidth+dsafety:(xdecimal+1+dsafety));
   write(f,'  ',linesize:1:3);
   writeln(f,'   user defined line: linetype, m, b, linesize');
   ul:=ul^.next
end;
comment(f);
if userlines = nil then writeln(f,' (none)');
comment(f);
bar(f,'*'); (* ************************** *)
comment(f);
bar(f,'*'); (* ************************** *)

comment(f);
write(f,' ', edgecontrol);
write(f,' ',edgeleft :pwid:pdec);
write(f,' ',edgeright:pwid:pdec);
write(f,' ',edgelow  :pwid:pdec);
write(f,' ',edgehigh :pwid:pdec);
writeln(f,' edgecontrol (p=page), edgeleft, edgeright, edgelow, edgehigh in cm');

comment(f);
writeln(f,' ', parameterversion:4:2,
             ' version of xyplop');

comment(f);
writeln(f); (* 2010 oct 21 close the comment!! SHEESH!!!!!! *)

end (* of do with p *)
end; (* of writeparam *)
(* end module xyplo.writeparam *)

(* begin module xyplo.regressiondata *)
procedure regressiondata(var f: text; p: xyploparam);
(* give the data on the regression lines *)
var
   us: symbolsptr; (* pointer to the user symbols *)
   z: real; (* Fisher's z' *)
function inverseFisher(z: real): real;
(* compute the inverse Fisher function
http://sportsci.org/resource/stats/sscorr.html#fisherz
*)
begin
  inverseFisher := (exp(2*z) - 1)/(exp(2*z) + 1);
end; (* inverseFisher *)
begin
(* writeln(f,'.KS'); (@ start a keep, so that the lines stay
                        together if there is a page break *)
(* writeln(f,'.nf'); (@ stop filling *)
(*writeln(output,'in regressiondata');*)
   us := p.usersymbols;
   while us <> nil do with us^.s do begin
      if doline then with p,l do begin
      (* writeln(f,'.sp 1'); (@ give room on output *)
comment(f);
      writeln(f);
comment(f);
         write(f,'Data on the regression line for symbol "',symboltoplot,'" ');
         write(f,'with flag "');
         writestring(f,symbolflag);
         writeln(f,'":');
comment(f);
         writeln(f,n:1,' data points');
comment(f);
         writeln(f,'mean     x = ',ex:xwidth+dsafety:(xdecimal+1+dsafety),
                         ' and y = ',ey:xwidth+dsafety:(xdecimal+1+dsafety));
comment(f);
         writeln(f,'variance x = ',
                          varx:xwidth+dsafety:(xdecimal+1+dsafety),
                         ' and y = ',vary:xwidth+dsafety:(xdecimal+1+dsafety));
comment(f);
         writeln(f,'  r = ',r:xwidth+dsafety:(xdecimal+1+dsafety));
comment(f);
         writeln(f,'r^2 = ',(r*r):xwidth+dsafety:(xdecimal+1+dsafety));
         (* new as of 2006 Sep 7: *)
comment(f);
         write(f,'Fisher''s z'' = ');
         if (r > -1.0) and (r < 1.0) then begin
            z := 0.5* (ln(1+r) - ln(1-r));
            writeln(f, z :xwidth+dsafety:(xdecimal+1+dsafety));
comment(f);
            if n > 3 then begin
               write(f, '95% confidence limits on correlation r:');
               write(f,    ' ', inverseFisher(z-1.96/sqrt(n-3))
                        :xwidth+dsafety:(xdecimal+1+dsafety));
               write(f, ' to ', inverseFisher(z+1.96/sqrt(n-3))
                        :xwidth+dsafety:(xdecimal+1+dsafety));
            end
            else begin
               write(f, '(no confidence limits, n < 3)');
            end;
            writeln(f);
         end
         else begin
            writeln(f,' (undefined)');
         end;
comment(f);
         writeln(f,'covariance = ',
                          covxy:xwidth+dsafety:(xdecimal+1+dsafety));
comment(f);
         write(f,'regression: ');
         if p.logyscale then logstring(f,ylabel,ybase)
                        else writestring(f,p.ylabel);
         write  (f,' = ',m:xwidth+dsafety:(xdecimal+1+dsafety), ' * ');
         if p.logxscale then logstring(f,xlabel,xbase)
(* logstring(f,xlabel,xbase) *)
                        else writestring(f,p.xlabel);
         writeln(f,' + ',b:xwidth+dsafety:(xdecimal+1+dsafety));

         (* 2014 mar 05: rescue of Greg Alvord's method to compute
            the standard deviation of the y intercept.
            See the email at the end of the code
            This is incomplete - the parenthesis he gave don't match!

         writeln(f,' standard deviation of intercept:',
x:1
{
sqrt( (sum(Xsqr) / n*sum(X-Xbar)**2) * s).
}
s.d.(intercept)      =     sqrt( (sum(Xsqr) /
n*sum(X-Xbar)**2) * s.     ????)
                 );
*)


      end;
      us := us^.next
   end;
(*
   writeln(f,'.KE'); (@ end the keep @)
   writeln(f,'.fi'); (@ start filling again @)
   writeln(f,'.sp 1'); (@ give room on output @) *)
end;
(* end module xyplo.regressiondata *)

(* begin module xyplo.skipcopy *)
procedure skipcopy(var infile,outfile: text; copy: boolean;
                   var linenumber: integer);
(* skip and copy from infile to outfile until no more asterisks found,
if copy is true.  if copy is false, just skip.  return the current
line number of the file *)
var   copying: boolean; (* continue copying *)
begin
      linenumber := 0;
      reset(infile);
      if eof(infile) then begin
         writeln(output,'skipcopy: xyin file is empty');
         halt
      end;
      copying := true;
      while copying do begin
(*writeln(output,'skipcopy: while copying...');*)
        if eof(infile)
        then copying := false
        else if eoln(infile)
        then begin
           copying := false
        end
        else begin
           if (infile^='*') or (infile^='#')
           then begin
              if copy then begin
                         write(outfile,'% '); (* PostScript Comment *)
                         copyaline(infile,outfile)
                      end
                      else readln(infile);
              linenumber := linenumber + 1
           end
           else copying := false
        end
     end;
end;
(* end module xyplo.skipcopy *)

(* begin module xyplo.grabdata *)
procedure grabdata(var infile: text; p: xyploparam; var x,y, xs,ys: real;
                   var hue,saturation,brightness: real;
                   var flag: char;
                   var flagstring: string;
                   var linenumber: integer;
                   var gotdata: boolean);
(* obtain the x and y data points, the x size and y size (xs, ys)
of the symbol (if relevant) and the symbol s from one line of the infile.
the flag for the kind of the symbol to plot is one character in flag,
while the entire string is in flagstring.
The color parameters are returned in hue, saturation and brightness.
Keep track of the linenumber in case something goes wrong.
The routine skips blank lines and lines that begin with an "*".
If no data were found, gotdata is false. *)
var
   column: integer; (* the column of the datum *)
   done: boolean; (* done reading a data line (or eof) *)
   droppoint: boolean; (* mechanism for dropping non-positive points on log
                          scale *)
   gotx, goty, (* we have gotten the x and y coordinates *)
   gotxs, gotys, (* we have gotten the x and y symbol sizes *)
   gots: boolean; (* we got the symbol to plot *)
   gothu, gotsa,gotbr: boolean; (* got the color columns *)
   satisfied: boolean; (* true when we have data for all regressions
                          on the current line *)
   xrawdata, yrawdata: real; (* raw x and y values before conversion to logs *)
procedure die;
(* tell the user where we were and stop the program *)
begin (* to die *)
   writeln(output,'at line ',linenumber:1,' of data (INCLUDING * lines)');
   halt
end; (* die *)
procedure grab;
var
   c: integer; (* index to the flagstring for removing underscore characters *)

begin
(*writeln(output,'in grabdata');*)
   droppoint := false;
   hue := 0.0;
   saturation := 0.0;
   brightness := 0.0;
   column := 0;
   satisfied := false;
   gots := not p.needscolumn;
   gotx := false;
   goty := false;
   gotxs := not p.needxscolumn;
   gotys := not p.needyscolumn;
   gothu := not p.needcocolumn;
   gotsa := not p.needcocolumn;
   gotbr := not p.needcocolumn;
   flag := ' '; (* control the value of this *)
(*writeln(output,'needxscolumn=',p.needxscolumn);*)
(*writeln(output,'needyscolumn=',p.needyscolumn);*)
   clearstring(flagstring);

   (* read in the column data *)
   while not satisfied do begin
      column := column + 1;
      skipblanks(infile); (* always move to start of column *)
      if eoln(infile) then begin
         writeln(output,'found end of line before all data columns',
                        ' were found');
         if not gots  then writeln(output,'missing symbol flag column');
         if not gotx  then writeln(output,'missing x data column');
         if not goty  then writeln(output,'missing y data column');
         if not gotxs then writeln(output,'missing x symbol size column');
         if not gotys then writeln(output,'missing y symbol size column');
         if not gothu then writeln(output,'missing hue column');
         if not gotsa then writeln(output,'missing saturation column');
         if not gotbr then writeln(output,'missing brightness column');
         die
      end;
(*writeln(output,'not satisfied next column = ',column:1);*)
(*writeln(output,'gots,gotx,goty,gotxs,gotys',gots,gotx,goty,gotxs,gotys);*)
(*writeln(output,'next char: ',infile^);*)
      if column = p.xcolumn then begin
         read(infile, x);
         if p.logxscale then begin
            if x <= 0.0 then begin
               tellwarning;
               writeln(warnings,'WARNING: all x values must be positive',
                                ' for logxscale');
               writeln(warnings,'The value you have is ',x:10:5);
               writeln(warnings,'IT WAS DROPPED FROM THE DATA SET');
               droppoint := true;
            end
            else begin
               xrawdata := x;
               x := ln(x)/p.logxbase;
               droppoint := false;
            end
         end;
         gotx := true
      end
      else if column = p.ycolumn then begin
         read(infile, y);
         if p.logyscale then begin
            if y <= 0.0 then begin
               (* be more lenient about this problem!! *)
               tellwarning;
               writeln(warnings,'WARNING: all y values must be positive',
                                ' for logyscale');
               write  (warnings,'At x = ',x:10:5);
               writeln(warnings,' the value you have is y = ',y:10:5);
               writeln(warnings,'IT WAS DROPPED FROM THE DATA SET');
               droppoint := true;
            end
            else begin
               yrawdata := y;
               y := ln(y)/p.logybase;
               droppoint := false;
            end
         end;
         goty := true
      end
      else if column = p.xscolumn then begin
         if p.needxscolumn then begin
            read(infile, xs);
            if p.logxscale then begin
               if xs <= 0.0 then begin
                  writeln(output,'all x values for symbol size',
                                 ' must be positive for logxscale');
                  writeln(output,'The value you have is ',xs:10:5);
                  die;
               end;
{see notes below for ys
               xs := ln(xs)/p.logxbase;
}
            end;
            gotxs := true
         end
         else skipnonblanks(infile);
      end
      else if column = p.yscolumn then begin
         if p.needyscolumn then begin
            read(infile, ys);
            if p.logyscale then begin
               if ys <= 0.0 then begin
                  writeln(output,'all y values for symbol size',
                                 ' must be positive for logyscale');
                  writeln(output,'The value you have is ',ys:10:5);
                  die;
               end;

{according to calc.p procedure logrange,
   loga.estimate := ln(a.estimate)/ln(base);
but
   loga.error := a.error/ a.estimate;
So this should be ys := ys/y.  Unfortunately we need the original
value of y and we don't have that yet.  So we must
pospone the calculation until the entire data line has been
obtained.
It's wrong to do this here:
               ys := ln(ys)/p.logybase;
It's right to do this later:
               ys := ys/yrawdata;
}
            end;
            gotys := true
         end
         else skipnonblanks(infile);
      end

      (* color columns *)
      else if column = p.hucolumn then begin
         if p.needcocolumn then begin
            read(infile, hue);
            gothu := true
         end
         else begin
            hue := 0.0;
            skipnonblanks(infile);
         end
      end
      else if column = p.sacolumn then begin
         if p.needcocolumn then begin
            read(infile, saturation);
            gotsa := true
         end
         else begin
            saturation := 0.0;
            skipnonblanks(infile);
         end
      end
      else if column = p.brcolumn then begin
         if p.needcocolumn then begin
            read(infile, brightness);
            gotbr := true
         end
         else begin
            brightness := 0.0;
            skipnonblanks(infile);
         end
      end

      else if column = p.scolumn
           then if p.needscolumn then begin
{
write(output,'grabdata: BEFORE gettoken flagstring = "');
writestring(output,flagstring);
writeln(output,'"');
}
                                    gettoken(infile,flag,flagstring,gots);

{
write(output,'grabdata: AFTER gettoken flagstring = "');
writestring(output,flagstring);
writeln(output,'"');
}
                                    (* now convert the underscores into
                                       blank characters *)
                                    with flagstring
                                       do for c := 1 to length
                                       do if letters[c] = '_'
                                          then letters[c] := ' '
                                 end
                                 else skipnonblanks(infile)
           else skipnonblanks(infile); (* skip the column *)

      satisfied := gotx and goty and gotxs and gotys and gots
                   and gothu and gotsa and gotbr
   end;

   if  p.needxscolumn and p.logxscale then xs := xs/xrawdata;
   if  p.needyscolumn and p.logyscale then ys := ys/yrawdata;

   readln(infile);
(*write(output,'x,y=',x:4:2,'  ',y:4:2,'  flag=',flag,'  ');*)
(*writeln(output,'xs,ys=',xs:10:8,'  ',ys:10:8);*)
end;
begin
   gotdata := false;
   done := false;
   while not done do begin
      if eof(infile) (* close shop on end of file  *)
      then done := true
      else if eoln(infile) then begin (* skip blank lines *)
         readln(infile);
         linenumber := succ(linenumber)
      end
      else if (infile^ = '*') or (infile^ = '#')
      then begin (* skip comment lines *)
         readln(infile);
         linenumber := succ(linenumber)
      end
      else begin (* handle actual data *)
         grab;
         gotdata := true;
         done := true
      end
   end;
   if droppoint then gotdata := false;
{
write(output,'END grabdata: flagstring = "');
writestring(output,flagstring);
writeln(output,'"');
}
end;
(* end module xyplo.grabdata *)

(* begin module xyplo.bounds *)
procedure bounds(var infile, outfile: text; var p: xyploparam);
(* determine the bounds of the graph, and gather data and do the
liner regressions.  the procedure uses skipcopy
to get to the data and copy it to the outfile.
The procedure makes a pass through the data, and so does these
three things all at once. *)
var
   done: boolean; (* done with a loop *)

   linenumber: integer; (* the current line being read *)
   us: symbolsptr; (* a pointer to the user symbols *)

   noregressions: boolean; (* no regressions were done for the data
      read.  that is, there was no symbol identified for the data.
      this means that we must check to see if there is a 'g' symbol *)

   (* some variables for reading *)
   x,y: real; (* x and y coordinates *)
   xs,ys: real; (* x and y sizes of the symbol *)
   hu,sa,br: real; (* hue, saturation and brightness of the symbol *)
   flag: char; (* the flag for the symbol to use *)
   flagstring: string; (* the entire flag string for the symbol to use *)
   gotdata: boolean; (* data found *)
begin
(*writeln(output,'in bounds');*)
   skipcopy(infile,outfile,true,linenumber);

   with p do if setx then if logxscale and (not logxnormal) then begin
      if (xmin <= 0) or (xmax <= 0) then begin
         writeln(output,'xmin and xmax must be positive for log scale');
         halt
      end;
      xmin := ln(xmin)/logxbase;
      xmax := ln(xmax)/logxbase;
   end;

   with p do if sety then if logyscale and (not logynormal) then begin
      if (ymin <= 0) or (ymax <= 0) then begin
         writeln(output,'ymin and ymax must be positive for log scale');
         halt
      end;
      ymin := ln(ymin)/logybase;
      ymax := ln(ymax)/logybase;
   end;

   (* if there was no request for setting the symbols, regression
      is done only on the first column.  needscolumn determines whether
      or not there was a request, while doline of the first symbol
      tells if the line is to be regressed *)
   if not p.needscolumn and p.usersymbols^.s.doline
   then writeln(output,'regressing on first symbol only');

   (* clear the regression variables *)
   us := p.usersymbols;
   while us<>nil do begin
      with us^.s
      do if doline then with l
      do regress('c',x,y,sumx,sumy,sumxsqd,sumysqd,sumxy,
              ex,ey,varx,vary,covxy,r,m,b,n);
      us := us^.next
   end;

   (* read through the data to find bounds and tabulate regression *)
   while not eof(infile) do begin
      grabdata(infile,p,x,y,xs,ys,hu,sa,br,flag,flagstring,linenumber,gotdata);
      if gotdata then begin
(*writeln(output,'x,y=',x:4:2,'  ',y:4:2,'  flag=',flag);*)
{writeln(output,'output, x,y=',x:4:2,'  ',y:4:2,'  flag=',flag);}
         if not p.setx then begin
            if x<p.xmin then p.xmin:=x;
            if x>p.xmax then p.xmax:=x;
         end;
         if not p.sety then begin
            if y<p.ymin then p.ymin:=y;
            if y>p.ymax then p.ymax:=y;
         end;

         if not p.needscolumn then begin
            (* do the regression into the first symbol only *)
            with p.usersymbols^.s.l
            do regress('e',x,y,sumx,sumy,sumxsqd,sumysqd,sumxy,
                       ex,ey,varx,vary,covxy,r,m,b,n);
         end
         else begin (* regress into the appropriate symbol *)
            us := p.usersymbols;
            (* find the symbol *)
            while us<>nil do begin
               noregressions := true;
               with us^.s
               do if doline then if equalstring(flagstring , symbolflag)
                  then with l do begin
                     regress('e',x,y,sumx,sumy,sumxsqd,sumysqd,sumxy,
                                ex,ey,varx,vary,covxy,r,m,b,n);
                     noregressions := false;
                  end;
               us := us^.next
            end;
            if noregressions then begin (* regress into the first 'g' symbol *)
               (* look for a symbol that prints the flag itself as grabbag *)
               us := p.usersymbols;
               done := (us = nil);
               while not done do begin
                  done := (us^.s.symboltoplot = 'g') or
                          (us^.s.symboltoplot = 'G');
                  if not done then begin
                     us := us^.next;
                     done := (us = nil)
                  end
               end;
               if us <> nil then with us^.s
               do if doline
                  then with l
                     do regress('e',x,y,sumx,sumy,sumxsqd,sumysqd,sumxy,
                                ex,ey,varx,vary,covxy,r,m,b,n);
            end
         end
      end
   end;

   (* clear the regression variables *)
   us := p.usersymbols;
   while us<>nil do begin
      with us^.s
      do if doline then with l
      do begin
         if n = 0 then begin
            write(output,'no data found for symbol ',symboltoplot);

            write(output,' (which has flag ');
            writestring(output,symbolflag);
            write(output,'):');

            writeln(output,' no regression done');
            doline := false;
         end
         else if n = 1 then begin
            write(output,'only one data point found for symbol ',symboltoplot);

            write(output,' (which has flag ');
            writestring(output,symbolflag);
            write(output,'):');

            writeln(output,' no regression done');
            doline := false;
         end
         else regress('r',x,y,sumx,sumy,sumxsqd,sumysqd,sumxy,
              ex,ey,varx,vary,covxy,r,m,b,n);
      end;
      us := us^.next
   end;

   (* given the bounds, we can now determine the scale factor *)
   with p do begin
      if xmax=xmin then begin
         writeln(output,' no variation to x range of graph');
         xscale := 1
      end
      else xscale := xsize/(xmax-xmin);
      if ymax=ymin then begin
         writeln(output,' no variation to y range of graph');
         yscale := 1
      end
      else yscale := ysize/(ymax-ymin);
   end
end;
(* end module xyplo.bounds *)

(* begin module xyplo.locate *)
function locatex(p: xyploparam; x: real): real;
(* convert x to the correct plot location in cm *)
begin locatex := (x-p.xmin)*p.xscale (* +p.xzero *) end;

function locatey(p: xyploparam; y: real): real;
(* convert y to the correct plot location in cm *)
begin locatey := (y-p.ymin)*p.yscale (* +p.yzero *) end;
(* end module xyplo.locate *)

(* begin module xyplo.crosshairs *)
procedure crosshairs(var outfile: text; p: xyploparam);
(* make the crosshairs of the graph *)
begin
   with p do begin
      (* put on the cross hairs *)
      if crosshairs then begin
         (* vertical crosshairs: *)
         if (xmin < 0.0) and (xmax > 0.0) then begin
            movea(outfile,locatex(p,0.0),locatey(p,ymin));
            yaxis(outfile, doaxisline, locatey(p,ymax)-locatey(p,ymin),
              ymin,(ymax-ymin)/yinterval,ymax,ysubintervals,
              -0.05,-0.2,0.0,0,0,logynormal,logyscale,logybase);
            yaxis(outfile, doaxisline,
              ysize, ymin,(ymax-ymin)/yinterval,ymax,ysubintervals,
              +0.05,-0.2,0.0,0,0,logynormal,logyscale,logybase);
         end;
         (* horizontal crosshairs: *)
         if (ymin < 0.0) and (ymax > 0.0) then begin
            movea(outfile,locatex(p,xmin),locatey(p,0.0));
            xaxis(outfile, doaxisline, locatex(p,xmax)-locatex(p,xmin),
              xmin,(xmax-xmin)/xinterval,xmax,xsubintervals,
              -0.05,0.0,-0.12,0,0,logxnormal,logxscale,logxbase);
            xaxis(outfile, doaxisline, xsize, xmin,(xmax-xmin)/xinterval,
              xmax,xsubintervals,
              +0.05,0.0,-0.12,0,0,logxnormal,logxscale,logxbase);
         end
      end
   end
end;
(* end module xyplo.crosshairs *)

(* begin module xyplo.makexlabel *)
procedure makexlabel(var outfile: text; p: xyploparam);
(* make the x label centered just below the x axis *)
begin
   with p do begin
      (* put the labels inside the picture: *)
      movea(outfile, locatex(p,(xmax-xmin)/2.0 + xmin),
                     locatey(p,ymin)- labelbelowx
{zb}
                    (* move the x axis label down to account for the  
                    new tic lengths! *)
                    -ticscale*Xticscale*Xticlength
                    + Xlabelshift);
      if logxscale
      then if logxnormal then begin
              graphstring(outfile,xlabel,'l');
              writeln(outfile,'(, log scale) x');
           end
           else loglabel(outfile,xlabel,xbase,'c')
      else graphstring(outfile,xlabel,'c');
   end
end;
(* end module xyplo.makexlabel *)

(* begin module xyplo.makeylabel *)
procedure makeylabel(var outfile: text; p: xyploparam);
(* make the y label left justified above the y axis.
If 'rotate is true, rotate the label to the left side *)
const
   ydebug = false; (* set to true to debug this procedure *)
var
   shiftleft: real; (* amount to shift the y label left *)
   number: real; (* the number to use for shifing the y label left *)
   isminus: integer; (* 1: the number has a minus sign; 0: none *)
   digits: integer; (* digits to represent the number *)
begin
if ydebug then writeln(output,'********** debugging makeylabel');
   with p do begin if rotateylabel then begin
      (* rotate the label and put it on the left side of the graph *)
if ydebug then write  (output,'********** "');
if ydebug then writestring(output,ylabel);
if ydebug then writeln(output,'" = ',ylabel.length:1,' characters');

      (* compute how much to move left based on:
         characters in ymin or ymax
         number of characters
         width of characters
      *)
      shiftleft := 0.00;
if ydebug then writeln(output,'********** shiftleft = ',shiftleft:1:1);
      if abs(ymin) > abs(ymax)
      then begin
         number := abs(ymin);
if ydebug then writeln(output,'********** ymin = ',ymin:5:1);
         if ymin < 0
         then isminus := 1
         else isminus := 0
      end
      else begin
         number := abs(ymax);
if ydebug then writeln(output,'********** ymax = ',ymax:5:1);
         if ymax < 0
         then isminus := 1
         else isminus := 0
      end;
if ydebug then writeln(output,'********** number   = ',number:5:1);
if ydebug then writeln(output,'********** isminus  = ',isminus:1);

if ydebug then writeln(output,'********** ywidth   = ',ywidth:1);
if ydebug then writeln(output,'********** ydecimal = ',ydecimal:1);

      (* digits code from proc stringinteger: *)
      if number > 9
      then digits := trunc(ln(number+0.1)/ln(10))+1
      else digits := 1;
      if ydecimal > 0
      then digits := digits + 1 + ydecimal; (* period plus decimals *)
if ydebug then writeln(output,'********** digits = ',digits:1);

      (* shift like in makexlabel: *)
      shiftleft := shiftleft + ticscale*Yticscale*Yticlength;
if ydebug then writeln(output,'********** shiftleft = ',shiftleft:1:1);

      (* shift by approximate character width in cm: *)
      shiftleft := shiftleft + 0.3*(digits+1+isminus); (* 1 is the decimal *)
if ydebug then writeln(output,'********** shiftleft = ',shiftleft:1:1);

{
      if Ylabelshift <> 0.0 (* it is set to exactly 0.0 initially *)
      then shiftleft := Ylabelshift; (* fudge factor from person *)
This would make the shift ignore the number width.
}
      shiftleft := shiftleft + Ylabelshift; (* fudge factor from person *)
if ydebug then writeln(output,'********** shiftleft = ',shiftleft:1:1);

{zzz}

      movea(outfile, locatex(p,xmin)-shiftleft,
                     locatey(p,(ymax-ymin)/2.0 + ymin) );
{
bad:
                     locatey(p,(ymax)/2.0) + ymin);
based on:
      movea(outfile, locatex(p,(xmax-xmin)/2.0 + xmin),
                     locatey(p,ymin)- labelbelowx
}

      writeln(outfile, 'gsave');
(* 2024 Jan 17: LaTeX/skim fail to show rotated y axis.  Replace
      writeln(outfile, '90 rotate'); OLD original
      writeln(outfile, '45 rotate 45 rotate'); That failed.
*)
      writeln(outfile, '90 rotate');
      if logyscale
      then if logynormal then begin
              graphstring(outfile,ylabel,'c');
              writeln(outfile,'(, log scale) x');
           end
           else loglabel(outfile,ylabel,ybase,'c')
      else graphstring(outfile,ylabel,'c');
      writeln(outfile, 'grestore');
   end

   else begin (* original method: label above the y axis *)
      movea(outfile, locatex(p,xmin),
                     locatey(p,ymax)+labelabovey);
      if logyscale
      then if logynormal then begin
              graphstring(outfile,ylabel,'l');
              writeln(outfile,'(, log scale) x');
           end
           else loglabel(outfile,ylabel,ybase,'l')
      else graphstring(outfile,ylabel,'l');

   end;
   end;
end;
(* end module xyplo.makeylabel *)

(* begin module xyplo.makeaxes *)
procedure makeaxes(var outfile: text; p: xyploparam);
(* make the axes of the graph *)
begin
   with p do begin

      movea(outfile, locatex(p,xmin),locatey(p,ymin));
      if doXaxis
         then xaxis(outfile,                          (* names in xaxis: *)
                    doaxisline,                       (* " *)
                    locatex(p,xmax)-locatex(p,xmin),  (* axlength *)
                    xmin,                             (* fromtic *)
                    (xmax-xmin)/xinterval,            (* interval *)
                    xmax,                             (* totic *)
                    xsubintervals,                    (* " *)
               ticscale*Xticscale*Xticlength,         (* length %T *) 
                    XaxisXshift,                      (* dx *)
{zb}
                    XaxisYshift,                      (* dy *)
                    xwidth,                           (* width *)
                    xdecimal,                         (* decimal *)
                    logxscale,                        (* " *)
                    logxnormal,                       (* " *)
                    logxbase);                        (* " *)

      if doYaxis then yaxis(outfile, doaxisline,
         locatey(p,ymax)-locatey(p,ymin),
         ymin,(ymax-ymin)/yinterval,ymax,ysubintervals,
         ticscale*Yticscale*Yticlength,YaxisXshift,YaxisYshift,
         ywidth,ydecimal,logyscale,logynormal,logybase);
   end;
end;
(* end module xyplo.makeaxes *)

(* begin module xyplo.plotdata *)
procedure plotdata(var infile, outfile, warnings: text; p: xyploparam);
(* read from infile and plot the data on outfile using the
parameters in p *)
const
   maxbad = 10; (* maximum number of bad characters that can be handled *)
var
   b: integer; (* index to badlist *)
   badlist: array[1..maxbad] of char; (* a list of flags that have
     been rejected for plotting *)
   badnumber: integer; (* current count of the number of bad flags *)
   foundbad: boolean; (* found the bad flag in our list *)
   linenumber: integer; (* the current line being read *)
   us: symbolsptr; (* points to a user defined symbol *)

   (* variables for reading data from the file *)
   x,y: real; (* x and y coordinates *)
   xs,ys: real; (* x and y sizes of the symbol *)
   hu,sa,br: real; (* hue, saturation and brightness of the symbol *)
   flag: char; (* the flag found in reading the data *)
   flagstring: string; (* the entire flag string found in reading the data *)
   gotdata: boolean; (* data found *)
   s: char; (* the symbol to use, converted from the flag *)
procedure findsymbol;
(* figure out which symbol to plot *)
var done: boolean; (* done searching for grab bag *)
begin
(*writeln(output,'findsymbol');*)
   (* we now have the flag, now find out which symbol it corresponds to *)
   us := p.usersymbols;
   with us^.s, p do begin
      if p.needscolumn
      then begin

(* debug
us := usersymbols;
while us <> nil do begin
  write(output,'  symbolflag="');
  writestring(output,us^.s.symbolflag);
  write(output,'"');
  write(output,' | symbol="' ,us^.s.symboltoplot,'"');
  writeln(output,'"');
  us := us^.next;
end;
*)

         (* first look for a match between the flagstring
            and the symbolflags *)
         us := usersymbols;
         while (not equalstring(flagstring , us^.s.symbolflag))
               and (us^.next <> nil)
               do us := us^.next;
         (* Note: if there is only one symbol, we can't advance.
            We also can't advance past the last symbol *)

         (* now, if the last one was the wrong one, kill it: *)
         if (not equalstring(flagstring , us^.s.symbolflag))
         then us := nil;

         (* look for a grabbag symbol *)
         if us=nil then begin (* if could not find the symbol *)

(*
            writeln(output,'looking for the grabbag symbol');
*)
            us := usersymbols;

            done := false;
            while not done do begin

(*
               if us <> nil then writeln(output,'looking at symbol "',
                                 us^.s.symboltoplot,'"');
*)

               if us = nil
               then done := true
               else if (us^.s.symboltoplot = 'g') or
                       (us^.s.symboltoplot = 'G')
                    then done := true
                    else us := us^.next;
            end;

            if us = nil then begin
               writeln(output,'Could not find a symbol for the flagstring',
                              ' and there is no grabbag.');
               write(output,'The flagstring is: "');
               writestring(output,flagstring);
               writeln(output,'"');
               writeln(output,'The known symbols are: ');
               us := usersymbols;
               while us <> nil do begin
                 write(output,'  symbolflag="');
                 writestring(output,us^.s.symbolflag);
                 write(output,'"');

                 write(output,' | symboltoplot="' ,us^.s.symboltoplot,'"');
                 writeln(output,'"');
                 us := us^.next;
               end;

               halt
            end
         end
      end
   end
end; (* findsymbol *)
procedure makesymbol;
(* make the symbol *)
var
   docolor: boolean; (* do the connection in color *)
begin
   (* use the symbol defined by us *)
   s := us^.s.symboltoplot;
   docolor := (s = 'C') or (s = 'B') or (s = 'M');

   with us^.s, p do begin
      (* decide if the point is within the field of the plot *)
      if (x>=xmin) and (x<=xmax) and
         (y>=ymin) and (y<=ymax) then begin

{
writeln(outfile,'% x = ',x:5:2);
writeln(outfile,'% y = ',y:5:2);
writeln(outfile,'% locatex(p,x) = ',locatex(p,x):5:2);
writeln(outfile,'% locatey(p,y) = ',locatey(p,x):5:2);
}
         if doconnection then begin
{
writeln(outfile,'% doconnection');
}
            if didlastpoint then begin
               if docolor then setcolor(outfile,colorkind,hu,sa,br);
               movea(outfile,locatex(p,oldx),locatey(p,oldy));
               drawa(outfile,locatex(p,x),locatey(p,y),
                     connecttype,connectsize);
                if docolor then setgray(outfile,0.0);
            end
            else movea(outfile,locatex(p,x),locatey(p,y));

            (* prepare variables for next time *)
            didlastpoint := true;
            oldx := x;
            oldy := y;
         end
         else movea(outfile,locatex(p,x),locatey(p,y));
{
writeln(outfile,'% done locating ');
}

         if symbolxsize <= 0.0 then xs := p.xscale*xs
                               else xs := symbolxsize;
         if symbolysize <= 0.0 then ys := p.yscale*ys
                               else ys := symbolysize;

         if s in ['b','d'] then begin
            (* make lines match on boxes *)
            writeln(outfile,'2 setlinecap');
         end;

              if s='b' then cboxr(outfile,xs,ys)
         else if s='B' then begin
                               setcolor(outfile,colorkind,hu,sa,br);
                               cboxr(outfile,xs,ys);
                               setgray(outfile,0.0);
                            end

         else if s='c' then circler(outfile,xs/2)

         else if s='C' then begin
                               setcolor(outfile,colorkind,hu,sa,br);
                               circler(outfile,xs/2);
                               setgray(outfile,0.0);
                            end

         else if s='d' then begin
                               dotr(outfile);
                               cboxr(outfile,xs,ys);
                            end
         else if s='p' then dotr(outfile)
         else if s='I' then ibeam(outfile,xs,2*ys) (* see documentation *)
         else if s='x' then xr(outfile,xs,ys)
         else if s='+' then plusr(outfile,xs,ys)

         else if s='r' then begin
                               setgray(outfile,br);
                               rectsize(outfile,xrect,yrect,xs,ys);
                               rectdo(outfile);
                               setgray(outfile,0.0);
                            end

         else if s='R' then begin
                               setcolor(outfile,colorkind,hu,sa,br);
                               rectsize(outfile,xrect,yrect,xs,ys);
                               rectdo(outfile);
                               setgray(outfile,0.0);
                            end

{
         else if s='s' then skybox(outfile,x,y,xs,ys)
}
         else if s='s' then skybox(outfile,locatex(p,x),locatey(p,y),xs,ys)

         else if s='L' then begin
                               setcolor(outfile,colorkind,hu,sa,br);
                               {2.54 cm/inch / 72 points/inch = cm / point}
                               writeln(outfile, xs / (2.54 / 72) :7:5,
                                       ' setlinewidth % connectsize');
                            end

         else if (s='m') or (s='M') then begin
                               if usermarks then begin
                                  setcolor(outfile,colorkind,hu,sa,br);
                                  if not p.needscolumn then begin
                                     writeln(output,'user defined marks in',
                                    ' xyplom require a (non-zero) symbol',
                                    ' column to be defined in xyplop');
                                     halt
                                  end;

                                  (* we have already looked for this
                                  flagstring in xyplom! Result is
                                  symbolhasxyplomdef *)
                                  if symbolhasxyplomdef
                                  then begin
                                     marksdo(outfile,xs,ys,flagstring);
                                     setgray(outfile,0.0);
                                  end;
                               end
                               else begin
                                  if not p.warnedaboutusermarks then begin
                                     writeln(output,'WARNING: NO USER',
                                                    ' MARKS DEFINED');
                                     p.warnedaboutusermarks := true;
                                  end
                               end
                            end

         else if (s='f') or (s='g')
              or (s='F') or (s='G')
              then begin
                 if (s='f') or (s='g')
                 then graphstring(outfile,flagstring,'c')
                 else graphstring(outfile,flagstring,'l')
              end
         else if not doconnection then begin
            write(output,'unknown symbol "',s,'", with flag "');
            writestring(output,flagstring);
            writeln(output,'": no symbol plotted and points not connected');
         end;

         if s in ['b','d'] then begin
            (* make lines end flat again *)
            writeln(outfile,'1 setlinecap');
         end;

         datacount := datacount + 1;
      end
      else begin
         datacount := datacount + 1;
         tellwarning;
         write(warnings,'skipping point # ',datacount:1,' (');

         if xdecimal > 0
         then write(warnings,x:xwidth:xdecimal)
         else write(warnings,round(x):xwidth);

         write(warnings,',');

         if ydecimal > 0
         then write(warnings,y:ywidth:ydecimal)
         else write(warnings,round(y):ywidth);

         writeln(warnings,') because it is out of the plotting window');

         didlastpoint := false
      end;
   end;

end; (* makesymbol *)
begin (* plotdata *)
(*writeln(output,'in plotdata');*)
   linenumber := 0;
   badnumber := 0;
   while not eof(infile) do begin
      grabdata(infile,p,x,y,xs,ys,hu,sa,br,flag,flagstring,linenumber,gotdata);
      if gotdata then begin

(*
writeln(output,'data grabbed');
writeln(output,'   flag = "',flag,'"');
writeln(output,'make the symbol xs,xy=',xs:10:5,ys:10:5);
writeln(output,'hu,sa,br',hu:10:5,sa:10:5,br:10:5);
*)

         findsymbol;
         if us=nil then begin
            (* first see if it is in the bad list already *)
            foundbad := false;
            b := 1;
            while (b <= badnumber) and (not foundbad)
            do if badlist[b] = flag then foundbad := true
                                    else b := b + 1;

            if not foundbad then begin
               if badnumber < maxbad then badnumber := badnumber + 1;
               badlist[badnumber] := flag;
               write(output,'flag "',flag,'" has no symbol');
               write(output,' and there is no grab-bag:');
               writeln(output,' this datum will not be plotted');
            end
         end
         else makesymbol;
      end
   end;
   (* this used to be in makesymbol *)
   if p.fontcontrol <> 'n'
   then begin
      writeln(outfile, '1 setlinewidth % in makesymbol'); (* reset linewidth *)
   end;
end;
(* end module xyplo.plotdata *)

(* begin module xyplo.plotuserlines *)
procedure plotuserlines(var outfile: text; p: xyploparam);
(* plot the user defined lines.  'n' means don't plot. *)
var
   ul: linesptr; (* a pointer to the user lines *)
   intercept: boolean; (* whether or not the regression line intercepts
      the region displayed *)
   x1,y1,x2,y2: real; (* interception points of the regression line
      on the display window *)
begin
(*writeln(output); *)
   (* create the user defined lines *)
   ul := p.userlines;
   while ul <> nil do begin
      with p, ul^.l do begin

(*
writeln(output);
writeln(output,'user line:',
               ' m = ', m:10:5,
               ' b = ', b:10:5);
*)
         boxintercept(xmin,ymin,xmax,ymax,m,b,intercept,x1,y1,x2,y2);
(*
writeln(output, ' intercept = ', intercept);

         if intercept then writeln(output,'TRUE INTERCEPT')
                      else writeln(output,'FALSE INTERCEPT');
*)

         if intercept then begin

(*
writeln(output,'box:',
            ' xmin=', xmin:10:5,
            ' ymin=', ymin:10:5,
            ' xmax=', xmax:10:5,
            ' ymax=', ymax:10:5,
            ' x1=', x1:10:5,
            ' y1=', y1:10:5,
            ' x2=', x2:10:5,
            ' y2=', y2:10:5);
*)

            if not (linetype in ['n','i']) then begin
               if not (linetype in ['.','-','l','n']) then begin
                  writeln(output,'error in plotuserlines: ');
                  writeln(output,'linetype (=',linetype,
                                 ') is not one of "l.-in"');
                  halt
{ uuu }
               end;

               movea(outfile,locatex(p,x1),locatey(p,y1));
               drawa(outfile, locatex(p,x2),locatey(p,y2),
                              linetype, linesize);
            end;
         end
      end;
      ul := ul^.next
   end;
end;
(* end module xyplo.plotuserlines *)

(* begin module xyplo.plotregressionlines *)
procedure plotregressionlines(var outfile: text; p: xyploparam);
(* plot the regression line for each symbol *)
var
   us: symbolsptr; (* a pointer to the user symbols *)
   intercept: boolean; (* whether or not the regression line intercepts
      the region displayed *)
   x1,y1,x2,y2: real; (* interception points of the regression line
      on the display window *)
begin
   with p do begin
      us := usersymbols;
      while us <> nil do with us^.s do begin
         if doline then with l do begin
            boxintercept(xmin,ymin,xmax,ymax,m,b,intercept,x1,y1,x2,y2);
            if intercept then begin
               movea(outfile,locatex(p,x1),locatey(p,y1));
               drawa(outfile, locatex(p,x2),locatey(p,y2), linetype, linesize);
            end
         end;
         us := us^.next
      end;
   end;
end;
(* end module xyplo.plotregressionlines *)

(* begin module xyplo.telllines *)
procedure telllines(var outfile: text; p: xyploparam);
(* describe the line intercepts *)
var
   ul: linesptr; (* a pointer to the user lines *)
   us: symbolsptr; (* a pointer to the user symbols *)

   intercept: boolean; (* whether or not the regression line intercepts
      the region displayed *)
   x1,y1,x2,y2: real; (* interception points of the regression line
      on the display window *)
begin
(*
   writeln(outfile,'.fi'); (@ start filling text again @)
   writeln(outfile,'.sp 1'); (@ give room on output @)
   writeln(outfile,'.KS'); (@ start a keep, so that the lines stay
                              together if there is a page break *)

   with p do begin
      us := usersymbols;
      while us <> nil do with us^.s do begin
         if doline then with l do begin
         (* writeln(outfile,'.sp'); (@ make space on the display *)
            boxintercept(xmin,ymin,xmax,ymax,m,b,intercept,x1,y1,x2,y2);
comment(outfile);
            write(outfile,'The regression line for symbol "',symboltoplot);
            write(outfile,'" with flag "');
            writestring(outfile,symbolflag);
            writeln(outfile,'"');
comment(outfile);
            if intercept then begin
               if logxnormal then begin
                  x1 := exp(x1*logxbase);
                  x2 := exp(x2*logxbase);
               end;
               if logynormal then begin
                  y1 := exp(y1*logybase);
                  y2 := exp(y2*logybase);
               end;
               writeln(outfile,'intercepts the display at points (',
                 x1:xwidth:(xdecimal+1), ', ',y1:ywidth:(ydecimal+1),
                 ') and (',
                 x2:xwidth:(xdecimal+1), ', ',y2:ywidth:(ydecimal+1),')')
            end
            else begin
               writeln(outfile,' does not intercept the display');
comment(outfile);
               writeln(output,'the regression line does not',
                              ' intercept the display')
            end
         end;
         us := us^.next
      end
   end;

   (* now tell the user defined lines *)
   ul := p.userlines;
   while ul <> nil do begin
   (* writeln(outfile,'.sp'); (@ make space on the display *)
      with p, ul^.l do begin
         boxintercept(xmin,ymin,xmax,ymax,m,b,intercept,x1,y1,x2,y2);
comment(outfile);
         writeln(outfile,'the user defined line of slope = ',m:10:5,
                         ' and intercept = ',b:10:5);

comment(outfile);
         if intercept then begin
            if logxnormal then begin
               x1 := exp(x1*logxbase);
               x2 := exp(x2*logxbase);
            end;
            if logynormal then begin
               y1 := exp(y1*logybase);
               y2 := exp(y2*logybase);
            end;
            writeln(outfile,' intercepts the display',
                 ' at points (',
                 x1:xwidth:(xdecimal+1), ', ',y1:ywidth:(ydecimal+1),
                 ') and (',
                 x2:xwidth:(xdecimal+1), ', ',y2:ywidth:(ydecimal+1),')');
         end
         else writeln(outfile,' does not intercept the display');
      end;
      ul := ul^.next
   end;
(* writeln(outfile,'.KE'); (@ end the keep @)
   writeln(outfile,'.nf'); (@ return to no filling text @)
   writeln(outfile,'.sp 2'); (@ make space on the display *)
end;
(* end module xyplo.telllines *)

(* begin module xyplo.makeusermarks *)
procedure makeusermarks(var xyplom, outfile: text;
                        var usermarks: boolean;
                        var warnedaboutusermarks: boolean;
                        var fontcontrol: char;
                        var ticscale: real;
                        var Xticscale: real;
                        var Yticscale: real;
                        var Xlabelshift: real;
                        var Ylabelshift: real;
                        usersymbols: symbolsptr);
(* Make the user marks file by copying it to the outfile.
For symbols that are m or M check for definition. *)
var
   control: char; (* control character *)
   us: symbolsptr; (* a symbol *)
   psdef: string; (* possible symbol flag definition in xyplom *)
   debug: boolean; (* set to true to see debug data *)
begin
   debug := true;
   debug := false;
   reset(xyplom);
   warnedaboutusermarks := false;
   if eof(xyplom)
   then usermarks := false
   else begin
      (* Check here that there is at least one defined mark
         This requires reading the PostScript file and figuring
         out what is in it. *)
(* Check here that each symbol flag is defined or mark it
as undefined.  When M, just don't plot. *)
(* are symbols defined already at this point?
they are read in by readparam so, yes, they are called before
makeusermarks.
*)
(*
if debug then begin
us := usersymbols;
while us <> nil do begin
  write(output,'  symbolflag="');
  writestring(output,us^.s.symbolflag);
  write(output,'"');
  write(output,' | symbol="' ,us^.s.symboltoplot,'"');
  writeln(output,'"');
  us := us^.next;
end;
end;
*)

      reset(xyplom);
      while not eof(xyplom) do begin
         if eoln(xyplom) then begin
            readln(xyplom);
         end
         else begin
            read(xyplom, control);
            if control = '/' then begin
               clearstring(psdef);
(*
               while (not eoln(xyplom)) and (control <> ' ') do begin
               while (not eoln(xyplom)) and (not (control in [' ','}']))
               while (not eoln(xyplom)) and (control <> ' ') WORKS
                                        and (control <> '{') WORKS
*)
               while (not eoln(xyplom)) and (not (control in [' ','}','%']))
               do begin
                   read(xyplom, control);
(*
                   if control <> ' ' then begin
                   if (control <> ' ') and (control <> '{') WORKS
*)
                   if not (control in [' ','{','%'])
                   then begin
                      psdef.length := psdef.length + 1;
                      psdef.letters[psdef.length] := control;
                   end;
               end;
if debug then begin
write(output,'Checking for symbol defined in xyplom: "/');
writestring(output,psdef);
writeln(output,'"');
end;
               if psdef.length > 0 then begin
                  us := usersymbols;
                  while us <> nil do begin
if debug then begin
  write(output,'  symbolflag="');
  writestring(output,us^.s.symbolflag);
  write(output,'"');
  write(output,' | symbol="' ,us^.s.symboltoplot,'"');
  write(output,'"');
end;

if debug then begin
write(output,'symbol: "');
writestring(output,psdef); writeln(output,'= psdef BUBBA00');
writestring(output,us^.s.symbolflag); writeln(output,'= us^.s.symbolflag BUBBA01');
end;
                    if equalstring(us^.s.symbolflag,psdef) then begin
if debug then begin
write(output,' is it');
end;

                       usermarks := true;
                       us^.s.symbolhasxyplomdef := true;
                    end
                    else begin
if debug then begin
write(output,' is not it');
end;
                    end;
if debug then begin
writeln(output)
end;
                    us := us^.next;
                  end;
               end;
            end;
         end
      end;

      (* copy the xylopm and read additional graph controls *)
      reset(xyplom);
      while not eof(xyplom) do begin
         (* watch for font trigger %F *)
         if not eoln(xyplom) then begin
            read(xyplom, control);
            write(outfile,control); (* preserve the character *)
            if control = '%' then begin
               if not eoln(xyplom) then begin
                  read(xyplom, control);
                  write(outfile,control); (* preserve the character *)
                  if control = 'F' then begin (* %F *)
                     fontcontrol := 'n'; (* set the font control!!! *)
                     write(output,'Fonts are controlled from xyplom, ');
                     writeln(output,'user is responsible for correct form!');
                     writeln(output,'Example for xyplom -----');
                     writeln(output,'%F % font control from xyplom');
                     writeln(output,'/Helvetica-Bold findfont');
                     writeln(output,'16 scalefont');
                     writeln(output,'setfont');
                     writeln(output,'2 setlinewidth % sets linewidth');
                     writeln(output,'%T 1.5 % tic mark scale factor');
                     writeln(output,'------------------------');
                  end;

                  if control = 'T' then begin (* %T *)
                     read(xyplom, ticscale);
                     writeln(outfile, ticscale:5:1,
                           ' %T tic mark scale factor');
                     writeln(output, ticscale:5:1,
                           ' %T tic mark scale factor');
                  end;

                  if control = 'X' then begin (* %X *)
                     read(xyplom, Xticscale);
                     writeln(outfile, Xticscale:5:1,
                           ' %X Xtic mark scale factor');
                     writeln(output, Xticscale:5:1,
                           ' %X Xtic mark scale factor');
                  end;

                  if control = 'Y' then begin (* %Y *)
                     read(xyplom, Yticscale);
                     writeln(outfile, Yticscale:5:1,
                           ' %Y Ytic mark scale factor');
                     writeln(output, Yticscale:5:1,
                           ' %Y Ytic mark scale factor');
                  end;

                  if control = 'x' then begin (* %x *)
                     read(xyplom, Xlabelshift);
                     writeln(outfile, Xlabelshift:5:1,
                           ' %x label shift');
                     writeln(output, Xlabelshift:5:1,
                           ' %x label shift');
                  end;

                  if control = 'y' then begin (* %y *)
                     read(xyplom, Ylabelshift);
                     writeln(outfile, Ylabelshift:5:1,
                           ' %y label shift');
                     writeln(output, Ylabelshift:5:1,
                           ' %y label shift');
                  end;

                  if control = 'W' then begin (* %W *)
{maybe read a color later
                     read(xyplom, redcolor, greencolor, bluecolor);
ccc
}
                     (* 2014 Oct 23: write white rectangle to background
                        so that one gets a white background when
                        converting to jpg using imagemagick. *)
                     writeln(outfile,'% white background on bounding box:');
                     writeln(outfile,'gsave');
                     writeln(outfile,'llx lly moveto');
                     writeln(outfile,'urx lly lineto');
                     writeln(outfile,'urx ury lineto');
                     writeln(outfile,'llx ury lineto');
                     writeln(outfile,'llx lly lineto');
                     writeln(outfile,'1 setgray');
                     writeln(outfile,'%1 1 0 setrgbcolor % for testing');
                     writeln(outfile,'fill');
                     writeln(outfile,'grestore');
                     writeln(outfile);
                  end;

               end;
            end;
         end;
         copyaline(xyplom, outfile);
      end;
   end;
{ The original code here was simply:
      while not eof(xyplom) do copyaline(xyplom, outfile);
}
end;
(* end module xyplo.makeusermarks *)

(* begin module xyplo.generategraph *)
procedure generategraph(var infile, outfile, warnings: text; p: xyploparam;
                        fontcontrol:char);
(* generate the graph in the outfile using the data in infile
and the parameters given in p *)
var
   linenumber: integer; (* the current linenumber, to keep skipcopy happy *)
begin
   skipcopy(infile,outfile,false,linenumber);

(* writeln(outfile,'.bp'); *)
   startpic(outfile,defscale,p.xzero,p.yzero,fontcontrol);
   if p.startrectangles then rectinit(outfile);

   (* plot data first so that lines and so forth can go on top *)

{yyy}
   (* 2010 Apr 09.  The problem is that there is a slight drift
      in postscript presumably due to numerical inaccuracy (it probably
      only holds 5 decimal places).  So the axes are not placed
      exactly in the same location plotting before and afterwards.
      The way around this is proper use of gsave/grestore.
      But xyplo was written before I converted it to PostScript,
      So it will take a careful job to rebuild its basic
      drawing routines.  Putting gsave/grestore around the call
      to plotdata does not resolve the issue.
      Increasing pwid and pdec did not help!
      *)

{
setcolor(outfile,0.5,1.0,1.0); (* make axes blue  before drift *)
makeaxes(outfile,p);
}

{
with p do begin
movea(outfile, locatex(p,xmin),locatey(p,ymin));
end;
writeln(outfile,'gsave % [ plotdata start drift protection');
}
   plotdata(infile, outfile, warnings, p);

{
writeln(outfile,'grestore % ] plotdata start drift protection');
}

   setgray(outfile,0.0); (* make sure that the axes etc are black! *)
   makeaxes(outfile,p);
{zb}
   makexlabel(outfile,p);
   makeylabel(outfile,p);

   crosshairs(outfile,p);
   plotuserlines(outfile, p);
   plotregressionlines(outfile, p);

   stoppic(outfile);
end;
(* end module xyplo.generategraph *)

(* begin module xyplo.themain *)
procedure themain(var infile, outfile, xyplop, xyplom, warnings: text);
(* the main procedure of the program, reading from infile for data
and producing results to outfile.  parameters obtained from xyplop *)
var
      p: xyploparam; (* the parameters controlling the program *)
begin (* the main *)
   writeln(output,'xyplo ',version:4:2);

   reset(infile);
   readln(infile); (* required for gpcc compiler to detect eof! *)
   if eof(infile) then begin
      writeln(output, 'xyin file is empty');
      halt
   end;
   reset(infile);

   rewrite(outfile);
   warningcount := 0;
   datacount := 0;

   readparam(xyplop,p); (* obtain the parameters *)

   clearstring(uniquestring);
   uniquestring.length := 6;
   uniquestring.letters[1] := 'U';
   uniquestring.letters[2] := 'N';
   uniquestring.letters[3] := 'I';
   uniquestring.letters[4] := 'Q';
   uniquestring.letters[5] := 'U';
   uniquestring.letters[6] := 'E';
{
write(output,'"');
writestring(output,uniquestring);
writeln(output,'" = uniquestring');
}

   writeln(outfile,'%!PS-Adobe-2.0 EPSF-2.0');
   writeln(outfile,'%%Title: xyplo ',version:4:2);
   writeln(outfile,'%%Creator: Tom Schneider');
   (* define bounding box.  The postscript program
   https://alum.mit.edu/www/toms/ftp/printerarea.ps
   will give the values for your printer.
   *)

   with p do begin
      writeln(outfile,'%%BoundingBox:',
                ' ',round(llx):5,
                ' ',round(lly):5,
                ' ',round(urx):5,
                ' ',round(ury):5);
      writeln(outfile,'%%DocumentFonts:');
      writeln(outfile,'%%EndComments');
      writeln(outfile,'%%EndProlog');

      writeln(outfile); (* this blank line makes sure that stupid
         printers think that this IS a real PostScript file, even
         though there are many comment lines following! *)

      writeln(outfile,'/defaultllx ',defaultllx:5:1,' def');
      writeln(outfile,'/defaultlly ',defaultlly:5:1,' def');
      writeln(outfile,'/defaulturx ',defaulturx:5:1,' def');
      writeln(outfile,'/defaultury ',defaultury:5:1,' def');
      writeln(outfile,'/llx ',llx:5:1,' def');
      writeln(outfile,'/lly ',lly:5:1,' def');
      writeln(outfile,'/urx ',urx:5:1,' def');
      writeln(outfile,'/ury ',ury:5:1,' def');

      writeln(outfile,'/cmfactor 72 2.54 div def % defines points -> centimeters');
      writeln(outfile,'/cm { cmfactor mul} def % defines centimeters');
      writeln(outfile,'/edgecontrol (',edgecontrol,') def');
      writeln(outfile,'/edgeleft  ',edgeleft: pwid:pdec,' cm def');
      writeln(outfile,'/edgeright ',edgeright:pwid:pdec,' cm def');
      writeln(outfile,'/edgelow   ',edgelow:  pwid:pdec,' cm def');
      writeln(outfile,'/edgehigh  ',edgehigh: pwid:pdec,' cm def');

   end;

   (* determine the bounds of the plot, run the regressions and
      calculate a few more things to keep around in p *)
   bounds(infile,outfile,p);

   writeparam(outfile,p); (* print the parameters *)

   regressiondata(outfile, p); (* print the regression data *)
   regressiondata(output, p); (* print the regression data *)

   telllines(outfile, p); (* describe the line intercepts *)
   telllines(output, p); (* describe the line intercepts *)

   (* Defaults for font control and tic scale: *)
   p.fontcontrol := 't';
   p.ticscale := 1.0;
   p.Xticscale := 1.0;
   p.Yticscale := 1.0;
   p.Xlabelshift := 0.0;
   p.Ylabelshift := 0.0;
   makeusermarks(xyplom, outfile, p.usermarks, p.warnedaboutusermarks,
                 p.fontcontrol,
                 p.ticscale, p.Xticscale, p.Yticscale,
                 p.Xlabelshift, p.Ylabelshift,
                 p.usersymbols);

   (* The bounds calls create warnings because the xyin file
   is read with grabdata.  The same set of warnings is generated
   in generategraph, which means that there would be two for each
   warning.  To avoid this, reset warningcount to zero.  This forces
   the file to be cleared and only one warning given per case. *)
   warningcount := 0;
   datacount := 0;

   (* create the graph: *)
   generategraph(infile, outfile, warnings, p, p.fontcontrol);

   if warningcount > 0 then begin
      write  (output,'**********');
      write  (output,' There are ',warningcount:1,
                     ' WARNINGS in file warnings');
      writeln(output,' **********');
   end;

   writeln(outfile,'%%Trailer');
   writeln(outfile,'%%Pages: 1');

(* writeln(outfile,'.fi'); (% start filling text again *)
end;
(* end module xyplo.themain *)

begin
      themain(xyin, xyout, xyplop, xyplom, warnings);
1: end.

{
From gwa@css.ncifcrf.gov Thu Jan  7 09:12 EST 1999
Received: from mail.ncifcrf.gov (mail.NCIFCRF.GOV [129.43.100.100]) by ncisun1-nf0.ncifcrf.gov (8.8.7/8.7.3) with ESMTP id JAA09003 for <toms@ncisun1.ncifcrf.gov>; Thu, 7 Jan 1999 09:12:14 -0500 (EST)
Received: from css.ncifcrf.gov (css.NCIFCRF.GOV [129.43.31.200])
        by mail.ncifcrf.gov (8.8.8/8.8.8) with ESMTP id JAA05676
        for <toms@alum.mit.edu>; Thu, 7 Jan 1999 09:12:18 -0500 (EST)
Received: by css.NCIFCRF.GOV with Internet Mail Service (5.5.1960.3)
        id <CJKZ25QV>; Thu, 7 Jan 1999 09:12:13 -0500
Message-ID: <E185275AD4A7D1119A2500A024BF2A42131655@css.NCIFCRF.GOV>
From: "Alvord, Greg" <gwa@css.ncifcrf.gov>
To: toms@alum.mit.edu
Cc: gwa@css.ncifcrf.gov
Subject: Standard deviation (error) of the intercept
Date: Thu, 7 Jan 1999 09:12:03 -0500 
MIME-Version: 1.0
X-Mailer: Internet Mail Service (5.5.1960.3)
Content-Type: text/plain
Content-Length: 1159


        Tom -

        Yes, the answer is straightforward.  The y intercept is, of
course, one of the parameters that is estimated in the linear
regression.  Since it is a parameter, the standard deviation is usually
referred to as the standard error.  It should be reported in just about
any statistical package that you would use to do linear regression.  

        In this e-mail, I don't have a simple, quick way to provide
mathematical notation, but

        Let s =  sqrt(MSE) in the regression analysis - the square root
of the mean square error,

                Xsqr = X squared,

                sum(Xsqr) = the summation of the all the Xsqr values,

                Xbar = mean of X's,

                n = number of observations (number of Y's),

                * = multiplication,

                ** = raise to power (i.e., square it if you see ...**2).

        Then the standard deivation (error) of the intercept is


                s.d.(intercept)      =     sqrt( (sum(Xsqr) /
n*sum(X-Xbar)**2) * s.

        I'm in conference all day.  Hope this helps. 

                        Greg.

W. Gregory Alvord, Ph.D.
Director, Statistical Services
Computer and Statistical Services
NCI-FCRDC
Phone:  (301) 846-5101
Fax:  (301) 846-6196
E-mail:  gwa@css.ncifcrf.gov
}
