program denplo(data, denplop, denploxyin, denploxyplom, denploxyplop,
               output);

(* denplo: density plot in color

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

 *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 2.23; (* of denplo.p 2022 May 04
2022 May 04, 2.23: add column to indicate diagonal.  yaverage column constant.
2017 Jul 10, 2.22: use bolder axes
2014 Oct 23, 2.21: allow shrinkfactorX shrinkfactorY to be > 1
2014 Mar 28, 2.20: key number rounding?
2014 Feb 26, 2.19: not upgrading correctly: two shrinkYav's
2013 Nov 25, 2.18: improve documentation for smoothkind
2013 Nov 16, 2.17: smoothing function - implement
2013 Nov 16, 2.16: smoothing function - almost ready to implement
2013 Nov 15, 2.15: binsarray under pointer as preparation for smoothing
2013 Nov 15, 2.14: correct shrinkfactorX message
2013 Nov 15, 2.13: axis control default is not to be log.
2013 Nov 15, 2.12: log axis works awesomely!!!!!!!! clean up
2013 Nov 15, 2.11: continue to implement log for axes
2013 Nov 15, 2.10: implement log for axes
2012 Mar 11, 2.09: clean up
2012 Mar 11, 2.08: put yaverage exactly at right place.
2012 Mar 07, 2.07: debug: don't compute y average if X is outside the range!
2012 Mar 04, 2.06: y average plot (clean up)
2012 Mar 04, 2.05: y average plot (functional)
2011 Oct 06, 2.04: clean up
2011 Oct 06, 2.03: debug crash: keyintervals is set by BinMax - if
                   flat distribution BinMax = 0 but keyintervals must
                   be > 0.
2010 Nov 02, 2.02: key width for fraction too wide
2009 Apr 13, 2.01: set BinMax in parameters
2009 Apr 11, 2.00: extend cutoff parameter too range
2009 Apr 02, 1.99: '('and')' in denploxyplop should NOT be protected with '\'!
2007 Jan 10, 1.98: fix putreal - taken from makelogo
2005 Dec 29, 1.97: spelling brightness
2005 Jun  1, 1.96: If the peak search starts outside the range of the
                   data, just give a warning and then let it do the
                   best it can.
2005 May 31, 1.95: more information when peak search fails
2005 May 29, 1.94: correct denplop upgrade to include SearchRange
2005 May 29, 1.93: make output file not wrap lines
2005 May 29, 1.92: make xyplop file not wrap lines
2005 Mar 30, 1.91: allow thelinewidth to be < 0 (no line around key boxes)
2004 Aug  2, 1.90: prepare xyin output for use by gnuplot: blanks, #
2004 Aug  1, 1.89: cleanup
2004 Aug  1, 1.88: cleanup
2004 Aug  1, 1.87: success!
2004 Aug  1, 1.86: SearchRange made into a parameter
2004 Aug  1, 1.85: peak finding algorithm is functional, black if exceeded!
2004 Aug  1, 1.84: core algorithm finds a near by peak!
2004 Aug  1, 1.83: test core algorithm! (stupid bug fixed)
2004 Aug  1, 1.82: test core algorithm!
2004 Aug  1, 1.81: Begin core algorithm! (report global max is ok)
2004 Aug  1, 1.80: allow switch to local max, set up initial Xbin,Ybin.
2004 Aug  1, 1.79: start to set up local maximum computation
2004 Jul 31, 1.78: set up local maxima variables
2004 Jul 31, 1.77: goal: find local maxima
2004 Feb 11, 1.76: cleanup
2004 Feb 11, 1.75: debugging ...
2004 Feb 11, 1.74: black strip bug: truncation problem in computing the bin
2004 Jan 20, 1.73: denplop.ri-example
2004 Jan 19, 1.72: complete keyintervals
2004 Jan 19, 1.71: make keyintervals logical
2004 Jan 19, 1.70: edge adjustment!!
2004 Jan 19, 1.69: fraction mode bug; using hue*0.84+0.16
2004 Jan 19, 1.68: brush up example documentation
2004 Jan 19, 1.67: clean up
2004 Jan 19, 1.66: implement reading color key parameters
2004 Jan 19, 1.65: color key parameters documented
2004 Jan 19, 1.64: color key matrix fully matches counts
2004 Jan 19, 1.63: color key max solved: must account for xyplo colors!
2004 Jan 19, 1.62: color key max is off, key numbers given.
2004 Jan 19, 1.61: color key max is off.
2004 Jan 19, 1.60: color key zero is correct.
2004 Jan 19, 1.59: color key partially functional
2004 Jan 19, 1.58: create color key in xyplom
2004 Jan 19, 1.57: change from rounding to truncating to drop into bins
2003 Aug 31, 1.56: Yintervals or Xintervals may be 1
2003 Aug 31, 1.53: MaxDisplayIntervals -> XDisplayIntervals YDisplayIntervals
2003 Aug 31, 1.52: greyscale or color parameter: cog
2003 Aug 31, 1.51: fix module structure at postscriptstring
2003 Aug 30, 1.50: need enough decimals in denploxyin to accurately plot!
2003 Aug 30, 1.49: link for movie 
2003 Aug 30, 1.48: refine documentation - denplo.data example
2003 Aug 30, 1.47: catch unused variable
2003 Aug 30, 1.46: Upgrade documentation, rename xyin denploxyin
2003 Aug 30, 1.45: add denplo.gif
2003 Aug 28, 1.44: compiling ...
2003 Aug 28, 1.43: necessary prgmod functions in ...
2003 Aug 28, 1.42: begin to create \t and \b functions ...
2003 Aug 28, 1.41: cleanup
2003 Aug 28, 1.40: rearrange xyin columns (and denxyplop)
2003 Aug 28, 1.39: cutoff is either counts or fraction of max
2003 Aug 28, 1.38: controls for X and Y axis labels (import strings)
2003 Aug 28, 1.37: controls for X and Y axis labels (documentation)
2003 Aug 28, 1.36: cleanup
2003 Aug 28, 1.35: background coloring functions.
2003 Aug 28, 1.34: correct background coloring
2003 Aug 28, 1.33: tie up loose ends; bins not filled correctly
2003 Aug 28, 1.32: correct placement of shrunken squares!
2003 Aug 28, 1.31: fine adjusting xyplop and other parameters
2003 Aug 28, 1.30: make sure intervals is as requested
2003 Aug 28, 1.29: lovely!! clean up
2003 Aug 28, 1.28: positioning with shrinkfactorX and shrinkfactorY
2003 Aug 28, 1.27: shrinkfactorX and shrinkfactorY
2003 Aug 28, 1.26: bug: SymbolSizeX and SymbolSizeY is jumpX, jumpY ...
2003 Aug 28, 1.25: bug: SymbolSizeX and SymbolSizeY are negative!
2003 Aug 28, 1.24: bug: location of squares is off ...
2003 Aug 28, 1.23: smart crosshairs
2003 Aug 28, 1.22: xyplop output is now called denploxyplop, crosshairs
2003 Aug 27, 1.21: MaxIntervals renamed MaxDisplayIntervals
2003 Aug 27, 1.20: cleanup
2003 Aug 27, 1.19: document fillonedatum error condition situation better
2003 Aug 24, 1.18: cutoffvalue
2003 Aug 24, 1.17: control over background where there are no data.
2003 Aug 24, 1.16: bug: always squares: SymbolSizeY not set correctly
2003 Aug 24, 1.15: further testing, MaxDisplayIntervals reinstated, documented.
2003 Aug 24, 1.14: cleanup; MaxDisplayIntervals variable DROPPED from parmeters.
2003 Aug 24, 1.13: correct location of squares (continued).
2003 Aug 24, 1.12: correct location of squares (continued).
2003 Aug 24, 1.11: correct location of squares again.
2003 Aug 24, 1.10: extend outer edge of graph to cover requested field
2003 Aug 24, 1.09: graph grey over entire requested field
2003 Aug 24, 1.08: report original coordinates AND square corner for plotting
2003 Aug 24, 1.07: streamline the code, center the squares
2001 Aug  4, 1.06: documentation upgrade
2001 Aug  4, 1.05: graph symbolsize under proper control
2001 Aug  4, 1.04: size and corner control
2001 Aug  4, 1.03: the plot is beautiful!
2001 Aug  4, 1.02: beginning to plot ...
2001 Aug  4, 1.01: able to read data ...
2001 Aug  3, 1.00: origin from denri
*)
updateversion = 2.16; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.denplo *)
(*
name
   denplo: density plot in color

synopsis
   denplo(data: in, denplop: inout,
          denploxyin: out, denploxyplom: out, denploxyplop: out,
          output: out)

files

   data:  data pairs

   denplop:  parameters to control the program.  The file must contain the
      following parameters, one per line:

      0. parameterversion: The version number of the program.  This allows
         the user to be warned if an old parameter file is used.

      1. numbertoprocess: number of data items to process (negative = all)
         This allows the user to explore a part of a large data set
         before committing to a long process time.  I suppose it would
         also allow you to make a movie of the data building up!

      2. XminValue, XmaxValue, Xintervals: (real, real, integer) range
         for x, and the number of intervals

      3. YminValue, YmaxValue, Yintervals: (real, real, integer) range
         for Y, and the number of intervals

      4. NuWi, NuDe: (integer) define the width and
         decimal places of output numbers.  NuWi = 1, Pascal will give
         the number out anyway!  This is a big space saver.

      5. Xcolumn, Ycolumn: (integer) column numbers for X and Y data.
         They can be in any order.

      6. Xsize, Ysize: (real) size of the plotting area, cm.

      7. Xcorner, Ycorner: (real) lower left corner of the plotting area, cm.

      8. XDisplayIntervals,    YDisplayIntervals
         XDisplaySubIntervals, YDisplaySubIntervals
         : (integer)

         Number of numeric intervals for xyplo to display on the axis
         This keeps the number of numbers on the side of the graph
         reasonable but does not affect the graph itself.  It allows
         one to have a highly dense plot without labeling every data
         interval.

         If one of the DisplayInterval values is less than or equal to
         zero, then the corresponding bin interval will be used.  For
         example, if XDisplayIntervals = 0 then the display interval
         used will be given by the value of Xintervals.

         Sub intervals apply all the time.

      9. xwidth, ywidth: (integer) axis number width  

     10. xdecimal, ydecimal: (integer) axis number decimals 

     11. background, grey (character, real) color to give when there
         are no data:

         b = black
         g = grey; the grey value is given (between 0.0 and 1.0)
         w = white
         c = colorful. For debugging:
             a complex formula is used to change the colors of squares
             over the background so one can debug the program.
             Data are plotted on top of this.  So most of the time
             this option is not useful to people.

         For white no background square needs to be written, so it
         isn't.  When there are sparse data above the cutoff, this
         option can be extremely efficient because very little data
         needs to be written to the denploxyin file.  For black, grey
         or colorful every square is written all the time.

     12. cutoffvalue (character, real, [real]):

         The first character on the line defines the kind of cutoff:

          'c': use the number of counts in a bin (one value used)
          'f': use the fraction of this bin to the fullest bin (one value used)
          'C': range of counts to plot (two values needed)
          'F': range of fractions to plot (two values needed)

         For 'c' and 'f', only one number is needed.  The number is
         the cutoff value.  Counts are reported for values higher than
         or equal to the cutoff.  For example, in c mode a zero cutoff
         will have denplo display cases with 0 or more counts.

         For 'C' and 'F', two numbers are need, this is the range of
         counts or fractions to use.

         2009 Apr 11: cutoff values are now inclusive, so if a count
         is 5 and the range is 5 to 10 it will be plotted.

     13. shrinkfactorX, shrinkfactorY (real): numbers above zero

         A value of 1 will make the display be completely covered
         while one less than 1 will leave gaps that allow one to
         distinguish between the squares.  Note: values larger than 1
         will cause overlapping squares, so something like 1.01 will
         slightly overlap the squares, guaranteeing that there will be
         no space at all between them.

     14. XaxisLabel: (string, entire line): the label for the x axis.

     15. YaxisLabel: (string, entire line): the label for the y axis.

     16. cog: (character): color or grey

         The first character on the line defines the value:

          'c': use colors in the spectrum
          'g': use grey scale

     COLOR KEY DISPLAY CONTROLS

     The color key is a set of rectangles surrounded by a black edge. 
     To the right are numbers that report the density.  The numbers
     can either report the number that the color corresponds to or the
     density relative to the highest density in the plot.  The key is
     placed into the denploxyplom file where it can then be passed to
     xyplo as a xyplom file.

     17. keyintervals: The number of boxes in the key will be
         keyintervals+1 (to allow for a zero).  If keyintervals is
         zero, then no color key is shown.  Example:  To get
         boxes numbered 0.0, 0.0 ... 10.0, use 10 for keyintervals.
     18. keytype (character): f: fraction of maximum, n: raw number.
         The way to display the key numbers.
     19. keyX: X coordinate of the lower left corner of the key (cm)
     20. keyY: Y coordinate of the lower left corner of the key (cm)
     21. keyXsize: X size of key boxes (cm)
     22. keyysize: Y size of key boxes (cm)
     23. keyshrinkfactor: a number between 0 and 1 which defines
         how much the key boxes should shrink.
     24. thelinewidth: the width of a line around each key box (cm)
         If thelinewidth is zero than the line is the default width.
         If less than zero, no line is drawn.
     25. keydecimals: number of decimal places for the
         numbers placed to the right of each key box.
     26. keyfontsize: the font size for the key numbers. 10 is a good
         place to stat.

     27. 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.

     FINDING LOCAL PEAKS

     28. findpeak, startX, startY, SearchRange:

         l: findpeak is a single character for which 'l' means to
         search for a local peak starting at coordinate (startX,
         startY) instead of using the global maximum.  If the local
         peak found is not the global maximum, then positions that
         have more counts than the local peak are shown by black.  The
         search is made in a series of squares, starting at (startX,
         startY).  The side of each square is determined by
         SearchRange: side = 2*Searchrange+1.  That is, if SearchRange
         = 0, the side is 1 and the peak is set at (startX, startY). 
         When SearchRange = 3, the side is 3 and 9 positions are
         searched.  The square is moved to the position that has the
         highest local density.  This is repeated until no higher
         position is found.  The algorithm is deterministic (no random
         numbers are used) and the result may depend on the order of
         scanning the square.

         maxpeak, BinMax:

         m: If the first character is 'm' then instead of finding the
         peak the second parameter is used for the maximum BinMax. 
         This allows one to force the maximum color to be a particular
         value and so allows comparison of plots on the same scale.

     29. yaverage, shrinkYav: (character, positive real) if 'y' then
         compute the average of y values in the bins and display.  The
         average is shown by a black bar with width the width of the
         displayed rectangles.  The height of the bar is determined by
         shrinkYav.  shrinkYav is the factor by which to scale the bar
         vertically.  If shrinkYav = 1, then the bar will be the same
         size as the display rectangles.  (Note that its position will
         vary up and down depending on where the average is so it will
         not necessarily match the rectangles.)  If shrinkYav = 2 then
         the bar will be half the size of the rectangles (vertically).

     30.  zxl xbase  if zxl='l' then convert the x axes to a log scale
                     using the indicated base.

     31.  zyl ybase  if zyl='l' then convert the y axes to a log scale
                     using the indicated base.
          The resulting log plot is very pretty!

     31.  smoothkind (char), smoothsize (real)
          If smoothkind is '-' then no smoothing is done.

          If smoothkind is 's' then smoothing is in a square.  A
          square of side smoothsize cells is averaged and placed into
          the center (if smootsize is odd) or to the positive side (if
          smoothsize is even).

          Once the data are collected into the first array, they are
          averaged into a second array of the same size as the first
          array.  The edges are lost (set to zero) depending on the
          size.  The second array is then printed.

   denploxyin: control file for xyplo: density of data pairs
      represented in color.  Copy or link this file to xyin.

   denploxyplom: control file for xyplo: a postscript file
      that contains marks.  It also contains the color key if
      that is being used.  Copy or link this file to xyplom.

   denploxyplop: control file for xyplo: parameters
      Copy or link this file to xyplop.
      Documentation is in the xyplo.p program.

   output: messages to the user

description

   Denplo takes a file containing pairs of numbers and (using the
   xyplo program) plots the density of the numbers in color.

   The user specifies a rectangular region on the X-Y plane.  Then
   each axis is divided into the number of intervals requested by the
   user.  These define a set of data bins onthe X-Y plane.  The data
   are dropped into the bins and counted.  The output of denplo
   consists of the three files that control the xyplo plotting
   program.

   For example, one output file is called a denploxypop = a denplo
   xyplop = density plot, x-y plotter parameter file.  This is
   pronounced 'den-plo-zye-plop'!  The denploxyin is 'den-plo-zin' and
   the denploxyplom is 'den-plo-zye-plom'.  :-)

   Since these are not the names used by xyplo (so that xyplo files
   are not overwritten), you need to copy (Unix cp) or move (Unix mv)
   them into the names that xyplo uses:

   cp denploxyin   xyin
   cp denploxyplom xyplom
   cp denploxyplop xyplop

   Then you can run xyplo to get the xyout (xy-out = 'zyout'), which
   is a eps PostScript file.

examples

   example denplop file:

1.88     0. parameterversion: denplo version this parameter file is designed for
-1       1. numbertoprocess: number of data items to process (negative = all)
0 10 10  2. XminValue,XmaxValue,Xintervals: (real,real,integer) X range,ints
0 10 10  3. YminValue,YmaxValue,Yintervals: (real,real,integer) Y range,ints
8 5      4. NumberWidth, NumberDecimals: (integer)
1 2      5. Xcolumn, Ycolumn: (integer) column numbers for X and Y data.
10 10    6. Xsize, Ysize: (real) plotting area, cm
3  7     7. Xcorner, Ycorner: (real) lower left corner of the plotting area, cm.
10 10 2 2 8. X/YDisplayIntervals X/YDisplaySubIntervals
8 8      9. xwidth, ywidth: (integer) axis number width
0 0     10. xdecimal, ydecimal: (integer) axis number decimals
g 0.9   11. background: b = black, g = grey, w = white, c = colorful (debug),
c 0 100 12. cutofftype [c,f,C,F] cutoffvalue(min and max)
0.8 0.8 13. shrinkfactorX, shrinkfactorY (real) 
X axis label
Y axis label
c       16. cog: color or grey plot
5       17. keyintervals: # of boxes in the key (not including the zero)
n       18. keytype (character): f: fraction of maximum, n: raw number.
10.2    19. keyX: X coordinate of the lower left corner of the key (cm)
0.1     20. keyY: Y coordinate of the lower left corner of the key (cm)
0.5     21. keyXsize: X size of key boxes (cm)
1.0     22. keyYsize: Y size of key boxes (cm)
0.8     23. keyshrinkfactor: between 0 and 1: box shrink.
0.005   24. thelinewidth: the width of a line around each key box (cm)
1       25. keydecimals: number of decimal places for the key
12      26. keyfontsize: the font size for the key numbers
p 1.50 2.50 1.50 1.50 edgecontrol, edgeleft, edgeright, edgelow, edgehigh cm 
- 0 0   28. findpeak(l), startX, startY, SearchRange | maxpeak(m), BinMax

see also

   {Example data file (copy to 'data' to use):} denplo.data
   {Example Parameter file:} denplop

   {Given these two files, you can run denplo.  To do this:
      cp denplo.data  data
      denplo
      cp denploxyin   xyin
      cp denploxyplom xyplom
      cp denploxyplop xyplop
    followed by running} xyplo.p {:}
      xyplo}

   {The final results use the general plotting program:} xyplo.p

   {The resulting example density plot is:}
   https://alum.mit.edu/www/toms/icons/denplo.gif

   {--------------------------------------}

   {a larger example, with large key scale:} denplop.ri-example

   {--------------------------------------}

   {Information about PostScript:}
   https://alum.mit.edu/www/toms/postscript.html
   {This describes how to convert from PostScript to PDF.}

   {Program that describes label controls:} makelogo.p
   {See bugs for further information on this partially implemented
   feature.}

   {How to make a movie:}
   https://alum.mit.edu/www/toms/paper/ev/movie/

   {Related density program:} denri.p

   {Information about the color hue formula:} diana.p xyplo.p

author

   Thomas Dana Schneider

bugs

   Smarter crosshairs:  xyplo needs an upgrade so that the crosshairs
   can be put in for x and y and independently the x and y axes.
   Then this program could control them.

   One could implement strings like '\c' within the axes labels to
   produce data counts into the labels.  (A mechanism like that in
   makelogo could be used.)

      \t total counts
      \b maximum bin count
      \i start italics and end italics

   The function in makelogo is putreal, inside module
   makelogo.postscriptstring in makelogo.p.  This routine is now in
   the program but not fully integrated so \t and \b are not available
   yet.  You can, however, do italics and special symbols;
   instructions are in the makelogo.p program.

technical notes

   Constants set in the program:
   maxValue is the largest number of allowed bins allowed on each axis.

   Postscript hue colors range from red (at 0) to red (at 1).  To
   avoid this ambiguity, only part of the hue range is used.  I have
   found that this formula is useful for 0 <= r <= 1: hue := 0.84*r +
   0.16. See diana.p and xyplo.p module pic.setcolor.  The
   denploxyplom has this correction built in since the xyplom is not
   color corrected.  For this reason, the values in the color scale in
   denploxyplom will not match those in the denploxyin.  However, when
   plotted with xyplo they are identical.

   The edge control method is from xyplo.p.

*)
(* end module describe.denplo *)

(* begin module denplo.const *)
   minValue = 0; (* lowest  Value allowed (range of data storage) *)
   maxValue = 1000; (* highest Value allowed (range of data storage) *)

   (* 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 *)

   komment = '#'; (* the first character of comments in denploxyin *)

(* end module denplo.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.73; (@ of dops.p 2003 Aug 28 *)

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 = 4.86; (@ of prgmod.p 2004 Sep 8 *)
(* end module interact.type version = 7.67; {of delmod.p 2004 Sep 8} *)

   (* bins to store the data.  By making it a pointer,
      a second array can be created only as needed.*)
   binsarrayptr = ^binsarray;
   binsarray = array[minValue..maxValue, minValue..maxValue] of integer;

var
   data, (* file used by this program *)
   denplop, (* file used by this program *)
   denploxyin, (* file used by this program *)
   denploxyplom, (* file used by this program *)
   denploxyplop: text; (* file used by this program *)

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

(* 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 = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module copyaline *) (* for transfering header info *)
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 = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module skipblanks *)

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

function isblank(c: char): boolean;
(* is the character c blank or tab? *)
begin
   isblank := (c = ' ') or (c = chr(9))
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 = 7.67; {of delmod.p 2004 Sep 8} *)

(* begin module clearstring *)
procedure clearstring(var ribbon: string);
(* empty the string *)
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; (* 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. *)
begin (* initializestring *)
   clearstring(ribbon);
   ribbon.next := nil;
end; (* initializestring *)
(* end module clearstring version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* 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 = 4.86; (@ of prgmod.p 2004 Sep 8 *)
 
(* 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 = 4.86; (@ of prgmod.p 2004 Sep 8 *)
(* end module interact.writestring version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* begin module copystring *)
procedure copystring(a: string; var b: string);
(* copy string a to b *)
var l: integer; (* index to the string *)
begin
   b.length := a.length;
   for l := 1 to a.length do b.letters[l] := a.letters[l]
end;
(* end module copystring version = 4.86; (@ of prgmod.p 2004 Sep 8 *)
(* begin module makelogo.protectcharacter *)
procedure protectcharacter(c: char;
                       var protectioncharacter: char;
                       var needed: boolean);
(* In PostScript, special characters must be protected against.
This routine looks at a character c and returns a protection character
if it is needed.
The parenthesis is used in PostScript to indicate the bounds of a
string, while the percent is the comment character.  The backslash also needs
protection, since it is the escape to indicate that the next character is part
of the string. *)
begin
   if c in ['(',')','%','\']
   then begin
      protectioncharacter := '\';
      needed := true
   end
   else begin
      protectioncharacter := ' ';
      needed := false
   end
end;
(* end module makelogo.protectcharacter version = 4.86; (@ of prgmod.p 2004 Sep 8 *)
(* begin module numberdigit *)
function numberdigit(number, logplace:integer): char;
(* return the digit at the place value ('logplace') position of number.
example:
      numberdigit(13625, 3) = 3
      numberdigit(13625, 4) = 1
2000 July 30 'myabsolute' replaced 'absolute', which
is apparently a keyword for GPC.  The name is kept
to keep the code looking similar to its origin.
*)
var
      place: integer; (* the exponent of logplace *)
      count: integer; (* used to make place *)
      myabsolute: integer; (* the absolute value of number *)
      acharacter: char; (* the character to be returned *)
procedure digit;
(* extract a digit at the place position *)
var
      tenplace: integer; (* ten times place *)
      z: integer; (* an intermediate value *)
      d: integer; (* the digit extracted *)
begin (* digit *)
      tenplace:=10*place;
      z:=myabsolute-((myabsolute div tenplace)*tenplace);
      if place = 1
         then d:=z
         else d:= z div place;
      case d of
         0: acharacter:='0';
         1: acharacter:='1';
         2: acharacter:='2';
         3: acharacter:='3';
         4: acharacter:='4';
         5: acharacter:='5';
         6: acharacter:='6';
         7: acharacter:='7';
         8: acharacter:='8';
         9: acharacter:='9';
      end
end; (* digit *)
procedure sign;
(* put a negative sign out or a positive sign *)
begin (* sign *)
      if number <0 then acharacter:='-'
                   else acharacter:='+'
end; (* sign *)
begin (* numberdigit *)
      place:=1;
      for count:=1 to logplace do place:=10*place;

      if number=0 then begin
         if place=1 then acharacter:='0'
                    else acharacter:=' '
      end
      else begin
         myabsolute:=abs(number);
         if myabsolute < (place div 10)
            then acharacter:=' '
            else if myabsolute >= place
                    then digit
                    else sign
      end;
      numberdigit:=acharacter
end; (* numberdigit *)
(* end module numberdigit version = 4.86; (@ of prgmod.p 2004 Sep 8 *)
(* begin module numbersize *)
function numbersize(n: integer):integer;
(* calculate amount of space to be reserved for the integer n *)
const ln10 = 2.30259; (* natural log of 10 - for conversion to log base 10 *)
      epsilon = 0.00001; (* a small number to correct log base 10 errors *)
var
   size: integer; (* intermediate result *)
begin (* numbersize *)
   if n = 0
   then numbersize:=1
   else begin
      size:=trunc(ln(abs(n))/ln10 + epsilon) + 1;
      (* the 1 is for the last digit *)
      (* the epsilon assures that we do not lose a place due to
         roundoff.  eg, sometimes log base 10 of 10 would be 0.9999
         instead of 1, and we would not do it right...  note: this will
         fail for very large numbers on the order of 1/epsilon. *)
      if n < 0 then size := succ(size);  (* account for minus sign *)
      numbersize := size;
   end
end; (* numbersize *)
(* end module numbersize version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

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

(* begin module postscriptstring *)
procedure postscriptstring(var instring, outstring, rawstring: string;
{
                           var symvec: text;
                           lowest, highest: integer;
}
                           rs, sd: real;
                           n: integer{;
                           havers, havesd: boolean}
);
(*
version = 2.00 of postscriptstring  2003 Aug 29
2003 Aug 29, 2.00: origin from makelogo.postscriptstring

Copy the instring to the outstring while protecting the
string by blocking postscript specific characters.  Example:
If the user typed    the result should be     so they get
      (                      \(                    (
     \)                       )                    \) (ie real postscript)
In the first case the user just said "I want a '('.  To do this,
we have to put protection on it.
In the second case the user said "I want a '(' as a postscript
command.  To do this, we have to remove the protection!  So this
procedure reverses the protection.  How curious.
The variable outlength reports the output length of the string
not counting real postscript commands so that the string
can be centered.  Of course since the user gives postscript
commands, they could do anything and that could mess up the
centering.  C'est la vie.
The rawstring is the string as the user would see it,
without control stuff or fancy postscript.

New commands:
   \i toggle italics on and off
   \n 5 produce number of sequences at coordinate 5
   \160 produce the Greek letter pi

New commands implemented 1999 Feb 20:
   \r produce Rsequence for the current range
   \s produce standard deviation for the current range

The symvec file is used to find the number of sequences
for the \n, \r and \d commands.  The lowest and highest values
are the range of the logo to use to compute Rs and sd.

This procedure requires these modules from prgmod.p:

   interact.const
   string.type
   interact.clearstring
   clearstring
   interact.getstring
   interact.writestring
   writestring
   copystring
   makelogo.protectcharacter
   numberdigit
   numbersize

 *)
const
   escape = '\'; (* the escape character is backslash *)
   infofield = 8; (* size of field for printing information in bits *)
   infodecim = 5; (* number of decimal places for printing information *)
   protecting = true; (* protect the user against PostScript requirements
      in user defined strings.  *)

var
   curr: char; (* the current character on the line *)
   decimals: integer; (* number of decimal places for next real number *)
   i: integer; (* index to s *)
   italic: boolean; (* if true, do italics *)
   insidestring: boolean; (* are we inside a postscript string? *)
   needed: boolean; (* is protection needed? *)
   prev: char; (* the previous character on the line *)
   protectionchar: char; (* is protection needed? *)

procedure getchar;
(* step forward one character *)
begin
   i := succ(i);
   prev := curr;
   curr := instring.letters[i];
end; (* getchar *)

procedure putchar(c: char);
(* put the character c onto the end of outstring *)
begin
   outstring.length := succ(outstring.length);
   if outstring.length > maxstring then begin
      writeln(output,'postscriptstring: label too long, increase maxstring');
      halt
   end;
   outstring.letters[outstring.length] := c;

   (* Count characters that are inside unprotected postscript parens.
      That is, count the characters that will really be printed *)
   if insidestring then begin
       rawstring.length := succ(rawstring.length);
       rawstring.letters[rawstring.length] := c;
{
writeln(output,'length = ',length:2,
               ' char = "',letters[length],'"',
               insidestring);
}
   end
end; (* putchar *)
procedure badstring;
begin
   writeln(output,'This string:');
   writestring(output,instring);
   writeln(output);
   writeln(output,'is bad because parenthesis must be paired like this:');
   writeln(output,'"\)"  stuff "\("');
{
   writeln(output,'I am emptying the logo file.');
   rewrite(logo);
}
   halt
end; (* badstring *)

procedure doitalics;
(* generate:
38\) \( E. coli \) IT \(LexA binding sites
38) (E. coli ) IT (LexA binding sites
*)
begin
   insidestring := false;
   if italic then begin (* complete italics *)
      putchar(')');
      putchar(' ');
      putchar('I'); putchar('T');
      putchar(' ');
      putchar('(');
   end
   else begin (* start italics *)
      putchar(')');
      putchar(' ');
      putchar('(');
   end;
   insidestring := true;
   italic := not italic
end; (* doitalics *)

procedure dosymbol;
(* from:
38 \160 SY LexA binding sites
generate:
38) (160) SY (LexA binding sites
*)
begin
   insidestring := false;
   putchar(')');

   (* the symbol counts as a single space inside the string *)
   insidestring := true;
   putchar(' ');
   insidestring := false;

   putchar('(');
   putchar('\');
   putchar(curr);
   (*get two and put them *)
   getchar; putchar(curr);
   getchar; putchar(curr);
   putchar(')');
   putchar(' ');
   putchar('S'); putchar('Y');
   putchar(' ');
   putchar('(');
   insidestring := true;
end; (* dosymbol *)

procedure getnum(var num: integer);
(* pick up a number from the string *)
var
   done: boolean; (* done finding number *)
   numberstarted: boolean; (* have we started reading the number? *)
   sign: integer; (* -1 or +1 *)
   firsti: integer; (* start point for reading *)
begin
   sign := +1;
   done := false;
   numberstarted := false;
{
writeln(output,'This string:');
writestring(output,instring);
writeln(output);
}
   firsti := i;
   num := 0;
   while not done do begin
      i := succ(i);
      if i > instring.length then begin
         done := true (* number defaults to zero *)
      end
      else begin
         curr := instring.letters[i];
{
writeln(output,'curr="',curr,'"');
}
         if curr = '-' then begin
            if sign = -1 then begin
               writeln(output,'strings can have only one "-" in \n numbers');
               halt
            end;
            sign := -1;
            numberstarted := true
         end
         else if curr = '+' then begin
            if sign = -1 then begin
               writeln(output,'you cannot have both + and - signs',
                              ' in \n numbers');
               halt
            end;
         end
         else if curr in ['0','1','2','3','4','5','6','7','8','9'] then begin
            num := 10*num + (ord(curr) - ord('0'));
{
writeln(output,'    current num:',num:1);
}
            numberstarted := true
         end
         else if (curr <> ' ') or numberstarted then begin
            (* allow leading blanks *)
            done := true
         end;
      end;
   end;
   num := sign * num;
{
writeln(output,'num=',num:1);
}

   (* set string reading variable for further analysis *)
   i := pred(i);
   prev := ' ';
   curr := instring.letters[i];

   if not numberstarted then begin
      writeln(output,'WARNING: a number was not found in this string:');
      writestring(output,instring);
      writeln(output);
      writeln(output,'zero is being used');
      num := 0;
      i := firsti;
   end;

end;

procedure don;
(* insert the number of sequences into the string *)
var
   c: char; (* a character perhaps to be part of the output number *)
{
   coo: integer; (* current coordinate in symvec *)
}
   num: integer; (* current number in symvec *)
   maxcoo: integer; (* maximum coordinate in symvec *)
   maxnum: integer; (* maximum number in symvec *)
{
   sd: real; (* current sd in symvec *)
   rs: real; (* current rs in symvec *)
}
   desiredcoo: integer; (* the coordinate to get the number of sequences *)
{
   symbols: integer; (* number of symbols *)
}
   s: integer; (* index to symbols *)
begin

   getnum(desiredcoo);

   (* locate the coordinate in the symvec *)
   maxnum := -maxint;
   maxcoo := -maxint;

{
   (* skip header *)
   reset(symvec);
   if not eof(symvec) then begin
      processsymvec(symvec);
      readln(symvec, symbols);
   end;

   coo := -maxint;
   while (not eof(symvec)) and (coo <> desiredcoo) do begin
      processsymvec(symvec);
      readln(symvec,coo,num,rs,sd);
      if num >= maxnum then begin
         maxnum := num;
         maxcoo := coo;
      end;
      for s := 1 to symbols do readln(symvec);
   end;

   if eof(symvec) then begin
      (* we could be eof(symvec) for two reasons.  One is that
         we read to the end of the file.  The other is that
         the file is empty. *)
      reset(symvec);
      if eof(symvec)
      then begin (* symvec is empty *)
         coo := 0;
         num := 0;
         writeln(output,'* The symvec is empty.');
         writeln(output,'  Using ',coo:1,' for the coordinate and ',
                        num:1,' for the number.');
      end
      else begin (* we read to the end *)
         writeln(output,'* Coordinate ',desiredcoo:1, ' is not in the symvec:');
         writeln(output,'  the maximum value (',maxnum:1,
                        ') at coordinate ',maxcoo:1, ' will be used instead.');
         coo := maxcoo;
         num := maxnum;
      end
   end;
}

   num := n;

   (* now put the number into the output strings *)
   for s := numbersize(num) downto 0 do begin
      c := numberdigit(num,s);
      if (c <> ' ') and (c <> '+') then putchar(c);
   end;

   n := num;
{
   writeln(output,'* ',n:1,' sequences found at coordinate ',coo:1);
}
   writeln(output,'* ',n:1,' sequences found at coordinate ,coo:1');
end; (* don *)

{
procedure getrssd;
(* get Rsequence and its sd *)
var
   b: integer; (* index to a symbol *)
   nl: integer; (* number of symbols at position l *)
   position: integer; (* a location in the aligned sequence (true coordinate) *)
   rsl: real; (* information at position l *)
   rsvar: real; (* variance of information at position l *)
   rstotal: real; (* sum of the rsl for the whole logo *)
   symbols: integer; (* number of symbols possible *)
   varhnb: real; (* variance of rstotal for the whole logo *)

begin
   reset(symvec);
   if eof(symvec)
   then begin (* symvec is empty *)
      writeln(output, 'symvec is empty: can''t do \r');
      halt
   end;

   processsymvec(symvec);
   readln(symvec, symbols);

   rstotal := 0.0;
   varhnb := 0.0;
   position := lowest;
   while (not eof(symvec) and (position <= highest)) do begin
      processsymvec(symvec);
      readln(symvec,position,nl,rsl,rsvar);

      if (position >= lowest) and (position <= highest) then begin
         rstotal := rstotal + rsl;
         varhnb := varhnb + rsvar;
      end;
 
      (* skip the symbols *)
      for b := 1 to symbols do begin
         processsymvec(symvec);
         readln(symvec); (* skip *)
      end;
   end;

   rs := rstotal;
   if varhnb >= 0
   then sd := sqrt(varhnb)
   else sd := 0.0;

   havers := true;
   havesd := true;

end;
}

procedure putreal(r: real; {wid,} dec: integer);
(* put a real number r into the string with dec decimal
places.  The total width of the number is not controlled;
there are no leading blanks.  If here are no decimal places,
the number is rounded.  If there are decimal places, the
number is rounded at the end.  (This is as one expects,
but the code is interesting.) *)
const
   debug = false; (* set true to debug *)
var
   c: char; (* part of the real number as a character *)
   m: integer; (* multiplier index *)
   num: integer; (* part of the real number *)
   s: integer; (* position in part of the number *)
begin
   if r < 0 then begin (* handle negative values *)
      c := '-';
      putchar(c);
      r := abs(r) (* set it positive for later work *)
   end;
   if dec > 0
   then begin
      (* 2007 Jan 10: BUG! *)
{
      (* move part below decimal above it: *)
      for m := 1 to dec do r := r * 10;
      (* round at that point only! *)
      r := round(r);
      (* put r back as it was, but rounded: *)
      for m := 1 to dec do r := r / 10;
      (* get part above decimal: *)
      num := trunc(r)
}
if debug then writeln(output, 'BB r initial   = ',r:10:5);

      (* move part below decimal above it: *)
      for m := 1 to dec do r := r * 10;
      (* round at that point only! *)
      r := round(r);
if debug then writeln(output, 'BB r rounded   = ',r:10:5);
      (* put r back as it was, but rounded: *)
      for m := 1 to dec do r := r / 10;
if debug then writeln(output, 'BB r divided   = ',r:10:5);
      (* get part above decimal: *)
      num := trunc(r);
if debug then writeln(output, 'BB r final     = ',r:10:5);
if debug then writeln(output, 'BB num = ',num:1);

   end
   else num := round(r);

   (* do part above decimal: *)
   for s := numbersize(num) downto 0 do begin
      c := numberdigit(num,s);
      if (c <> ' ') and (c <> '+') then putchar(c);
   end;

   (* do part below decimal: *)
   if dec > 0 then begin
      (* decimal point *)
      c := '.';
      putchar(c);
      (* digits after decimal point *)
      r := r - num; (* remove part above decimal point *)
if debug then writeln(output, 'BL r rm num    = ',r:10:5);

      (* move part below decimal above it: *)
      for m := 1 to dec do r := 10 * r;
if debug then writeln(output, 'BL r above     = ',r:10:5);
      num := round(r); (* round gives same result as Pascal *)
if debug then writeln(output, 'BL num         = ',num:10);
{
      (* 2007 Jan 10 bug here! - should be dec not numbersize! *)
      for s := numbersize(num) downto 0 do begin
}
if debug then writeln(output, 'BL numbersize(num) = ',numbersize(num):10);
if debug then writeln(output, 'BL dec             = ',dec            :10);
      for s := dec downto 0 do begin
         c := numberdigit(num,s);
{
2007 Jan 10: original code dropped spaces that should be zeros:
         if (c <> ' ') and (c <> '+') then putchar(c);
}
         (* 2007 Jan 10 Bug solution:
            spaces are zeros below the decimal! *)
         if c =' ' then c := '0';
         if (c <> '+') then putchar(c);
      end;
   end;
end;

procedure dors;
(* insert Rsequence into the string *)
begin
{
   if not havers then getrssd;
}
   putreal(rs, {infofield,} decimals);
   writeln(output,'* Rsequence for logo:',rs:infofield:decimals);
end;

procedure dosd;
(* insert standard deviation of Rsequence into the string *)
begin
{
   if not havesd then getrssd;
}
   putreal(sd, {infofield,} decimals);
   writeln(output,'* SD of Rsequence for logo:',sd:infofield:decimals);
end;

procedure dodecimal;
(* determine number of decimal places *)
begin
   getnum(decimals);
   writeln(output,'* decimal places set to: ',decimals:1);
end;

begin (* postscriptstring *)
   if protecting then begin
      clearstring(outstring);
      clearstring(rawstring);
      curr := ' ';
      decimals := infodecim; (* default value *)
      italic := false;
{in makelogo:
      insidestring := false; (* we start outside the string *)
EXCEPTION: here we STAY insidestring:
      insidestring := true; (* we start outside the string *)
      putchar('(');
}
      insidestring := true; (* we are now inside the string *)

      i := 0;
      while i < instring.length do begin;
         getchar;

         (* since we just got to a \) we are no longer inside the string *)
         if (prev = escape) and (curr = ')') then begin
            if not insidestring then badstring;
            insidestring := false;
         end;

         (* \\ gives \\ *)
         if (curr = escape) then begin
            if (prev = escape)
            then begin
               putchar(prev);
               putchar(curr);
            end;
         end
         else begin
            if prev <> escape 
            then begin
               protectcharacter(curr, protectionchar, needed);
               if needed then putchar(protectionchar);
            end;

            if (prev = escape) and (curr = 'i')
            then doitalics

            else if (prev = escape) and (curr = 'n')
            then don

            else if (prev = escape) and
                 (curr in ['0','1','2','3','4','5','6','7'])
            then dosymbol

            else if (prev = escape) and (curr = 'r')
            then dors

            else if (prev = escape) and (curr = 's')
            then dosd

            else if (prev = escape) and (curr = 'd')
            then dodecimal

            else putchar(curr);

         end;

         (* now that we are past the \( we are truely inside the string *)
         if (prev = escape) and (curr = '(') then begin
            if insidestring then badstring;
            insidestring := true;
         end;
      end;

      if insidestring
      then insidestring := false
      else begin
         writeln(output,'This string:');
         writestring(output,instring);
         writeln(output);
         writeln(output,'is bad because it needs a final "\(".');
{
         writeln(output,'I am emptying the logo file.');
         rewrite(logo);
}
         halt
      end;

      (* close italics if necessary *)
      if italic then doitalics;

{ DO NOT CLOSE STRING
      putchar(')');
}
   end
   else begin
      copystring(instring, outstring);
      copystring(instring, rawstring);
   end

{
;
write(output,' instring = "');
writestring(output,instring);
writeln(output,'"');

write(output,'outstring = "');
writestring(output,outstring);
writeln(output,'"');

write(output,'rawstring = "');
writestring(output,rawstring);
writeln(output,'"');
}

end; (* postscriptstring *)
(* end module postscriptstring *)

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

(* 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 = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* 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 = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* begin module denplo.themain *)
procedure themain(var data, denplop,
                  denploxyin, denploxyplom, denploxyplop: text);
(* the main procedure of the program *)
var
   background: char; (* control the background color *)
   grey: real; (* definition of grey *)

   cutofftype: char; (* the type of cutoff count to use,
                        c for counts or f for fraction. *)
   cutoffvaluemin: real; (* lowest  value to report *)
   cutoffvaluemax: real; (* highest value to report *)

{ original definition
   bins: array[minValue..maxValue, minValue..maxValue] of integer;
      (* bins to store the data *)
}
   bins: binsarrayptr; (* bins to store the data *)
   smoothingbins: binsarrayptr; (* bins for smoothing the data *)
   swapbins: binsarrayptr; (* a pointer to switch the bins *)
   ysums: array[minValue..maxValue] of real; (* store the sum of Y for the mean *)
   ycounts: array[minValue..maxValue] of integer; (* store Y counts for the mean *)
{zzz}
   numbertoprocess: integer; (* number of data items to process *)

   parameterversion: real; (* parameter version number *)

   NuWi: integer; (* width of output numbers *)
   NuDe: integer; (* decimals of output numbers *)

   XminValue, XmaxValue: real; (* range of the X values *)
   Xrange: real; (* Xrange := XmaxValue - XminValue; *)
   Xintervals: integer; (* number of intervals for bins of X values *)

   YminValue, YmaxValue: real; (* range of the Y values *)
   Yrange: real; (* Yrange := YmaxValue - YminValue; *)
   Yintervals: integer; (* number of intervals for bins of Y values *)

   Xcolumn, Ycolumn: integer; (* column numbers for X and Y data. *)
   xlessthany: boolean; (* Xcolumn is less than Ycolumn *)

   Xsize, Ysize: real; (* size of the plotting area, cm. *)

   Xcorner, Ycorner: real; (* lower left corner of the plotting area, cm. *)

   XDisplayIntervals, YDisplayIntervals:
       integer; (* maximum number of numeric
                   intervals for xyplo to display *)

   XDisplaySubIntervals, YDisplaySubIntervals:
       integer; (* maximum number of numeric
                   sub intervals for xyplo to display *)

   xwidth, ywidth: integer; (* axis number width *)
   xdecimal, ydecimal: integer; (* axis number decimals *)

   Xvalue, Yvalue: real; (* a pair of data items *)

   dataCount: integer; (* total data count *)

   inCount: integer; (* number of pairs in the X and Y range *)
   outCount: integer; (* number of pairs in the X and Y range *)

   XinCount: integer; (* number of pairs in the X range *)
   XoutCount: integer; (* number of pairs in the X range *)

   YinCount: integer; (* number of pairs in the Y range *)
   YoutCount: integer; (* number of pairs in the Y range *)

   Xin: boolean; (* the X value is inside *)
   Yin: boolean; (* the Y value is inside *)

   Xbin, Ybin: integer; (* current bin location for the data *)

   Xlo: integer; (* the lowest value of Xbin *)
   Ylo: integer; (* the lowest value of Ybin *)

   Xhi: integer; (* the highest value of Xbin *)
   Yhi: integer; (* the highest value of Ybin *)

   SymbolSizeX, SymbolSizeY: real; (* size of X and Y symbols *)
   shrinkfactorX, shrinkfactorY: real; (* shrink of Symbols relative to jumps *)

   jumpX: real; (* jump between X values *)
   jumpY: real; (* jump between Y values *)

   shiftX, shiftY: real; (* amount to shift  the square X and Y so that
       it is centered on the REGION of the bin *)

   XaxisLabel: string; (* a the X axis label *)
   YaxisLabel: string; (* a the Y axis label *)

   cog: char; (* color or grey plot *)

   BinMax: integer; (* the maximumly filled bin *)

   (* controls for the color key *)
   keytype: char; (* the type of number to report: fraction or number *)
   keyX, keyY: real; (* The coordinate to place the key, in cm. *)
   keyXsize, keyYsize: real; (* the size of the key in cm. *)
   keydecimals: integer; (* number of decimal places in the key *)
   keyintervals: integer; (* intervals on the color key *)
   keyshrinkfactor: real; (* factor by which to shrink the key boxes *)
   keyfontsize: real; (* the font size for the key numbers *)
   thelinewidth: real; (* the line width around the key boxes *)

   (* 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 *)

   findpeak: char; (* if 'l' then find peak starting at startX, startY *)
   startX: real; (* x coordinate of local peak finding *)
   startY: real; (* Y coordinate of local peak finding *)
   SearchRange: integer; (* 2*SearchRange+1 is the size of the side of
      a square within which the search for the maximum is done. *)

{zzz}
   yaverage: char; (* parameter to determin whether to compute y value averages *)
   shrinkYav: real; (* amount to shrink the yaverage mark by relative
     to the yjump *)
{ppp}

   (* 2013 Nov 15: log axis controls, variables from xyplo.p zzzlogaxis *)
         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 *)

   smoothkind: char; (* kind of smoothing to do.  - none; s square *)
   smoothsize: real; (* size of region to smooth *)
   smoothrange: integer; (* size of region to smooth, rounding of
                            smoothsize *)

(* begin module denplo.upgradefrom166to188 *)
procedure upgradefrom166to188(var denplop: text);
(* upgrade the denplop file from 1.66 to version 1.88.
Introduce artmode variable.
*)
const
   copylines = 27; (* the number of copied lines *)
var
   internal: text; (* a place to hold the old denplop *)
   line: integer; (* a line to be worked with *)
   parameterversion: real; (* parameter version number *)
begin
   parameterversion := 1.88;
   writeln(output, 'upgrading to version ',parameterversion:4:2,' ...');

   (* copy denplop to internal *)
   reset(denplop);
   readln(denplop); (* skip old parameter version line *)
   rewrite(internal);
   for line := 1 to copylines
   do copyaline(denplop, internal);

   (* write the NEW PARAMETER LINE *)
   writeln(internal,
      '- 0 0 1 28. findpeak(l), startX, startY, SearchRange');
{ppp}

   (* finish the copy of denplop to internal *)
   copyfile(denplop, internal);

   (* copy internal to denplop *)
   reset(internal);
   rewrite(denplop);
   writeln(denplop,parameterversion:4:2,'        version of',
         ' denplop that this parameter file is designed for.');
   copyfile(internal, denplop);

   (* add the new material at the end: *)
   (* none to add *)

   reset(denplop); (* ready to start reading again *)

end;
(* end module denplo.upgradefrom166to188 *)

(* begin module denplo.upgradefrom188to201 *)
procedure upgradefrom188to201(var denplop: text);
(* upgrade the denplop file from 1.88 to version 2.01.
Introduce artmode variable.
*)
const
   copylines = 28; (* the number of copied lines *)
var
   internal: text; (* a place to hold the old denplop *)
   line: integer; (* a line to be worked with *)
   parameterversion: real; (* parameter version number *)
begin
   parameterversion := 2.01;
   writeln(output, 'upgrading to version ',parameterversion:4:2,' ...');

   (* copy denplop to internal *)
   reset(denplop);
   readln(denplop); (* skip old parameter version line *)
   rewrite(internal);
   for line := 1 to copylines
   do copyaline(denplop, internal);

   (* write the NEW PARAMETER LINE *)
{
   writeln(internal, 'n ');
}
{ppp}

   (* finish the copy of denplop to internal *)
   copyfile(denplop, internal);

   (* copy internal to denplop *)
   reset(internal);
   rewrite(denplop);
   writeln(denplop,parameterversion:4:2,'        version of',
         ' denplop that this parameter file is designed for.');
   copyfile(internal, denplop);

{
2014 Feb 26 WRONG - second copy!
   (* add the new material at the end: *)
   writeln(denplop,'n 4     29. yaverage, shrinkYav:',
                   ' (char, real) "y": compute average of y values');
}
{zzz}

   reset(denplop); (* ready to start reading again *)

end;
(* end module denplo.upgradefrom188to201 *)


(* begin module denplo.upgradefrom201to210 *)
procedure upgradefrom201to210(var denplop: text);
(* upgrade the denplop file from 2.01 to version 2.10.
Introduce log scales.
*)
const
   copylines = 29; (* the number of copied lines *)
var
   internal: text; (* a place to hold the old denplop *)
   line: integer; (* a line to be worked with *)
   parameterversion: real; (* parameter version number *)
begin
   parameterversion := 2.10;
   writeln(output, 'upgrading to version ',parameterversion:4:2,' ...');

   (* copy denplop to internal *)
   reset(denplop);
   readln(denplop); (* skip old parameter version line *)
   rewrite(internal);
   for line := 1 to copylines
   do copyaline(denplop, internal);

   (* write the NEW PARAMETER LINE *)
{
   writeln(internal, 'n ');
}
{ppp}

   (* finish the copy of denplop to internal *)
   copyfile(denplop, internal);

   (* copy internal to denplop *)
   reset(internal);
   rewrite(denplop);
   writeln(denplop,parameterversion:4:2,'        version of',
         ' denplop that this parameter file is designed for.');
   copyfile(internal, denplop);

   (* add the new material at the end: zzzlogaxis *)
   writeln(denplop,'n 2     30. zxl (char), xbase (integer)');
   writeln(denplop,'n 2     31. zyl (char), ybase (integer)');

   reset(denplop); (* ready to start reading again *)

end;
(* end module denplo.upgradefrom201to210 *)

(* begin module denplo.upgradefrom210to216 *)
procedure upgradefrom210to216(var denplop: text);
(* upgrade the denplop file from 2.10 to version 2.16.
Introduce log scales.
*)
const
   copylines = 31; (* the number of copied lines *)
var
   internal: text; (* a place to hold the old denplop *)
   line: integer; (* a line to be worked with *)
   parameterversion: real; (* parameter version number *)
begin
   parameterversion := 2.16;
   writeln(output, 'upgrading to version ',parameterversion:4:2,' ...');

   (* copy denplop to internal *)
   reset(denplop);
   readln(denplop); (* skip old parameter version line *)
   rewrite(internal);
   for line := 1 to copylines
   do copyaline(denplop, internal);

   (* write the NEW PARAMETER LINE *)
{
   writeln(internal, 'n ');
}
{ppp}

   (* finish the copy of denplop to internal *)
   copyfile(denplop, internal);

   (* copy internal to denplop *)
   reset(internal);
   rewrite(denplop);
   writeln(denplop,parameterversion:4:2,'        version of',
         ' denplop that this parameter file is designed for.');
   copyfile(internal, denplop);

   (* add the new material at the end: smooth *)
   writeln(denplop,'- 2     32. smoothkind (char) "-s", smoothsize (integer)');

   reset(denplop); (* ready to start reading again *)

end;
(* end module denplo.upgradefrom201to210 *)

(* begin module denplo.upgradeparameters *)
procedure upgradeparameters(var denplop: text);
(* make sure that the parameters are the latest spiffy version *)
var
   parameterversion: real; (* parameter version number *)

(* Since parameters are real numbers, there can be round off errors in
the 15th decimal place that prevent comparisons.  However, I only use
parameters to 2 decimal places, so we can multipy by 100 and round to
do reliable comparisons. *)

function r100(a: real): integer;
(* simple multiply and round function *)
begin
 r100 := round(100*a);
end;

function equalversion(a,b: real): boolean;
(* Are versions a and b equal? *)
begin
   if r100(a) = r100(b)
   then equalversion := true
   else equalversion := false
end; (* equalversion *)

function lessthan(a,b: real): boolean;
(* Is version a less than b? *)
begin
   if r100(a) < r100(b)
   then lessthan := true
   else lessthan := false
end; (* lessthan *)

function greaterthan(a,b: real): boolean;
(* Is version a greater than b? *)
begin
   if r100(a) > r100(b)
   then greaterthan := true
   else greaterthan := false
end; (* greaterthan *)

begin
   readln(denplop, parameterversion);
   if    lessthan(parameterversion, updateversion) {or}
{
      greaterthan(parameterversion, versionupperbound)
}
   then begin
     writeln(output, '^GYou have an old parameter file!, version ',
              parameterversion:4:2,'!');
      writeln(output, '          version = ',          version:4:2);
      writeln(output, '    updateversion = ',    updateversion:4:2);
      writeln(output, ' parameterversion = ', parameterversion:4:2);
{
      writeln(output, 'versionupperbound = ',versionupperbound:4:2);
}

      if equalversion(parameterversion,1.66) then upgradefrom166to188(denplop);
      if lessthan(parameterversion,2.01) then upgradefrom188to201(denplop);
      (* implement log axis zzzlogaxis: *)
      if lessthan(parameterversion,2.10) then upgradefrom201to210(denplop);
      (* implement smoothing functions: *)
      if lessthan(parameterversion,2.11) then upgradefrom210to216(denplop);
{ppp}

      reset(denplop);
      readln(denplop, parameterversion);
      if lessthan(parameterversion , updateversion) then begin
         writeln(output, 'Sorry! I am unable to fully upgrade',
                         ' your parameter file');
         writeln(output, 'from version ', parameterversion:4:2,
                          ' to version ', updateversion:4:2,'!');
         writeln(output, 'Start from a fresh copy or edit this one.');
         halt;
      end
      else writeln(output, '... upgrade successful!');

      writeln(output, 'See this page for the new documentation:');
      writeln(output,
         'https://alum.mit.edu/www/toms/delila/denplo.html');
   end;

end;
(* end module denplo.upgradeparameters *)

procedure readparameters(var denplop: text);
(* read user defined parameters *)
var
   checkout: boolean; (* if true, all variable values are ok *)
   gotten: boolean; (* did we get the buffer? *)

procedure cn;
(* short version of call to check number *)
begin
  checkout := checknumber(denplop);
  if not checkout then begin
     writeln(output,'denplo readparameters: halt to avoid snoballing');
     writeln(output,'A number was expected but not found.');
     if eof(denplop)
     then writeln(output,'Unexpected End of File')
     else if eoln(denplop)
          then writeln(output,'Unexpected End Of Line')
          else writeln(output,'The bad character in denplop is: "',
                    denplop^,'"');
     halt; (* avoid snowballing *)
  end;
end;

procedure t;
(* test that the file is still there to read *)
begin
      if eof(denplop) then begin
         writeln(output,'Missing parameters: unexpected end of denplop');
         halt
      end
end;

begin
   reset(denplop);
   cn; upgradeparameters(denplop);

   cn;readln(denplop, numbertoprocess); (* 1 *)
   cn;readln(denplop, XminValue, XmaxValue, Xintervals); (* 2 *)
   cn;readln(denplop, YminValue, YmaxValue, Yintervals); (* 3 *)
   cn;readln(denplop, NuWi, NuDe); (* 4 *)
   cn;readln(denplop, Xcolumn, Ycolumn); (* 5 *)
   cn;readln(denplop, Xsize, Ysize); (* 6 *)
   cn;readln(denplop, Xcorner, Ycorner); (* 7 *)
   cn;readln(denplop, XDisplayIntervals,    YDisplayIntervals,
                   XDisplaySubIntervals, YDisplaySubIntervals ); (* 8 *)
   cn;readln(denplop, xwidth, ywidth); (* 9 *)
   cn;readln(denplop, xdecimal, ydecimal); (* 10 *)

   read(denplop, background);
   if background = 'g'
   then read(denplop, grey)
   else grey := 0.9;
   readln(denplop); (* 11 *)

   read(denplop, cutofftype, cutoffvaluemin); (* 12 *)
   if (cutofftype = 'C') or (cutofftype = 'F')
   then read(denplop, cutoffvaluemax)
   else cutoffvaluemax := maxint;
   readln(denplop);

   cn;readln(denplop, shrinkfactorX, shrinkfactorY); (* 13 *)

   getstring(denplop, XaxisLabel, gotten);
   if not gotten then begin
      writeln(output, 'end of file found when reading X axis label');
      halt
   end;

   getstring(denplop, YaxisLabel, gotten);
   if not gotten then begin
      writeln(output, 'end of file found when reading X axis label');
      halt
   end;

   readln(denplop, cog);

   cn;readln(denplop, keyintervals);
      readln(denplop, keytype);
   cn;readln(denplop, keyX);
   cn;readln(denplop, keyY);
   cn;readln(denplop, keyXsize);
   cn;readln(denplop, keyYsize);
   cn;readln(denplop, keyshrinkfactor);
   cn;readln(denplop, thelinewidth);
   cn;readln(denplop, keydecimals);
   cn;readln(denplop, keyfontsize);

   (* first, set defaults *)
   llx := defaultllx;
   ury := defaultury;
   urx := defaulturx;
   lly := defaultlly;

   read(denplop, edgecontrol);
   if (edgecontrol = 'p')
   then begin
      cn; read(denplop, edgeleft);
      cn; read(denplop, edgeright);
      cn; read(denplop, edgelow);
      cn; read(denplop, edgehigh);
{note: cn comes from xyplo.p}

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

      (* in xyplo.p xzero -> Xcorner; yzero -> Ycorner; *)
      llx := round((Xcorner         - edgeleft )* defscale);
      lly := round((Ycorner         - edgelow  )* defscale);
      urx := round((Xcorner + Xsize + edgeright)* defscale);
      ury := round((Ycorner + Xsize + edgehigh )* defscale);
   end;
   readln(denplop);

   (* read the start point for local search *)
   read(denplop, findpeak);
   if (findpeak = 'l')
   then begin
      cn; read(denplop, startX);
      cn; read(denplop, startY);
      cn; read(denplop, SearchRange);
   end
   else if (findpeak = 'm')
   then begin
      cn; read(denplop, BinMax);
      startX := 0;
      startY := 0;
      SearchRange := 0;
   end
   else begin
      startX := 0;
      startY := 0;
      SearchRange := 0;
   end;
   readln(denplop);
{ppp}

{zzz}
   readln(denplop, yaverage, shrinkYav);

   (* Read log axis parameters (code modified from xyplo) zzzlogaxis *)
         (* make type of x axis *)
         logxscale := false;
         logxnormal := false;
         xbase := 2;
         logxbase := ln(xbase);
         t; if not eoln(denplop)
         then if (denplop^='l') or (denplop^='L') then begin
            logxscale := true;
            if denplop^='L' then logxnormal := true;
            skipnonblanks(denplop);
            cn; read(denplop,xbase);
            if xbase <= 1.0 then begin
               writeln(output,'log x base must be > 1');
               halt
            end;
            logxbase := ln(xbase);
         end;
         readln(denplop);

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

   read(denplop, smoothkind);
   cn; readln(denplop, smoothsize); (* size of region to smooth *)

   (* check the input parameters ************************* *)

   if XmaxValue <= XminValue then begin
      writeln(output,'XmaxValue (',XmaxValue:NuWi:NuDe,')',
                     ' cannot be less than ',
                     'XminValue (',XminValue:NuWi:NuDe,')');
      halt;
   end;

   if YmaxValue <= YminValue then begin
      writeln(output,'YmaxValue (',YmaxValue:NuWi:NuDe,')',
                     ' cannot be less than ',
                     'YminValue (',YminValue:NuWi:NuDe,')');
      halt;
   end;

   if Xintervals < 1 then begin
      writeln(output,'Xintervals (=',Xintervals:1,')',
                     ' cannot be less than 2');
      halt;
   end;

   if Yintervals < 1 then begin
      writeln(output,'Yintervals (=',Yintervals:1,')',
                     ' cannot be less than 2');
      halt;
   end;

   if NuWi < 1 then begin
      writeln(output,'NuWi must be positive');
      halt;
   end;

   if NuDe < 1 then begin
      writeln(output,'NuDe must be positive');
      halt;
   end;

   (***** No!  If NuDe = 1, Pascal will fill the number out anyway!
          This is a big space saver.  2004 July 31.
   if NuDe > NuWi then begin
      writeln(output,'NuDe must be less than NuWi');
      halt;
   end;
   *******)

   if NuDe < 5 then begin
      writeln(output,'NuDe must be at least 5 decimal places');
      halt;
   end;

   if Xcolumn < 1 then begin
      writeln(output,'Xcolumn must be positive');
      halt;
   end;

   if Ycolumn < 1 then begin
      writeln(output,'Ycolumn must be positive');
      halt;
   end;

   if Xcolumn = Ycolumn then begin
      writeln(output,'X column cannot equal Ycolumn');
      halt;
   end;

   if Xsize <= 0.0 then begin
      writeln(output,'Xsize must be positive');
      halt;
   end;

   if Ysize <= 0.0 then begin
      writeln(output,'Ysize must be positive');
      halt;
   end;

   if XDisplayIntervals < 0 then begin
      writeln(output,'XDisplayIntervals must be zero or positive',
                     ' it is ', XDisplayIntervals:1);
      halt;
   end;

   if YDisplayIntervals < 0 then begin
      writeln(output,'YDisplayIntervals must be zero or positive',
                     ' it is ', YDisplayIntervals:1);
      halt;
   end;

   if XDisplaySubIntervals <= 0 then begin
      writeln(output,'XDisplaySubIntervals must be positive');
      halt;
   end;

   if YDisplaySubIntervals <= 0 then begin
      writeln(output,'YDisplaySubIntervals must be positive');
      halt;
   end;

   if xwidth <= 1 then begin
      writeln(output,'xwidth must be positive');
      halt;
   end;

   if ywidth <= 1 then begin
      writeln(output,'ywidth must be positive');
      halt;
   end;

   if xdecimal < 0 then begin
      writeln(output,'xdecimal must be positive or zero');
      halt;
   end;

   if ydecimal < 0 then begin
      writeln(output,'ydecimal must be positive or zero');
      halt;
   end;

   if not (background in ['b','g','w','c']) then begin
      writeln(output,'background must be one of bgwc');
      halt;
   end;

   if (grey < 0.0) or (grey > 1.0) then begin
      writeln(output,'grey must be between 0 and 1');
      halt;
   end;

   if not (cutofftype in ['c','f','C','F']) then begin
      writeln(output,'cutofftype must be one of "cfCF"');
      halt;
   end;

   if cutoffvaluemin < 0 then begin
      writeln(output,'cutoffvaluemin (',cutoffvaluemin:1:1,
                     ') must be positive or zero');
      halt;
   end;

   if cutoffvaluemax < 0 then begin
      writeln(output,'cutoffvaluemax (',cutoffvaluemax:1:1,
                     ') must be positive or zero');
      halt;
   end;

   if (cutofftype = 'f') or (cutofftype = 'F') then begin
      if (cutoffvaluemin < 0.0) or (cutoffvaluemin > 1.00) then begin
         writeln(output,'for cutofftype = "f",',
                        ' cutoffvaluemin must be between 0 and 1');
         halt;
      end;
   end;

   if (cutofftype = 'f') or (cutofftype = 'F') then begin
      if (cutoffvaluemax < 0.0) or (cutoffvaluemax > 1.00) then begin
         writeln(output,'for cutofftype = "f",',
                        ' cutoffvaluemax must be between 0 and 1');
         halt;
      end;
   end;

   if (cutofftype = 'C') or (cutofftype = 'F') then begin
      if (cutoffvaluemin > cutoffvaluemax) then begin
         writeln(output,'for cutoffvaluemin must be less than or equal to',
                        ' cutoffvaluemax ');
         halt;
      end;
   end;

{2014 Oct 23 Original restrictive code did not allow values > 1
   if (shrinkfactorX <= 0) or (shrinkfactorX > 1.0) then begin
      writeln(output,'shrinkfactorX must be > 0 and <= 1');
      halt;
   end;
   if (shrinkfactorY <= 0) or (shrinkfactorY > 1.0) then begin
      writeln(output,'shrinkfactorY must be > 0 and <= 1');
      halt;
   end;
}

   (* allow values > 1 to guarantee full overlap *)
   if shrinkfactorX <= 0 then begin
      writeln(output,'shrinkfactorX must be > 0');
      halt;
   end;

   if shrinkfactorY <= 0 then begin
      writeln(output,'shrinkfactorY must be > 0');
      halt;
   end;

   if not (cog in ['c','g']) then begin
      writeln(output,'cog (color or grey) must be one of "cg"');
      halt;
   end;

   if keyintervals < 0 then begin
      writeln(output,'No color key will be displayed');
   end;
      
   if not (keytype in ['f','n']) then begin
      writeln(output,'keytype must be one of "fn"');
      halt;
   end;

{
   if keyX then begin
      halt;
   end;

   if KeyY then begin
      halt;
   end;
}

   if keyXsize <= 0.0 then begin
      writeln(output,'keyXsize must be positive, it was ',
                      keyXsize :1);
      halt;
   end;

   if keyYsize <= 0.0 then begin
      writeln(output,'keyYsize must be positive, it was ',
                      keyYsize :1);
      halt;
   end;

   (* 2005 Mar 30: linewidth can now be negative - ie no line.
   if thelinewidth <= 0.0 then begin
      writeln(output,'thelinewidth must be positive, it was ',
                      thelinewidth :1);
      halt;
   end;
   *)

   if (keytype = 'f') and (keydecimals < 1) then begin
      writeln(output,'keydecimals must be positive for fractions.',
                     ' It was ', keydecimals:1);
   end;

{
not needed:
   else begin
      keydecimals := 0; (* Force zero decimal places for counts! *)
   end;
}
{yyy stupid, no decimals needed for counts!:
   if (keytype <> 'f') and (keydecimals < 0) then begin
      writeln(output,'keydecimals must be positive or zero for counts.',
                     ' It was ', keydecimals:1);
   end;
}

   if shrinkYav <= 0 then begin
      writeln(output,'shrinkYav must be positive');
      halt;
   end;

   if not (smoothkind in ['-','s']) then begin
      writeln(output,'smoothkind is "',smoothkind,'" BUT');
      writeln(output,'smoothkind must be one of "-s"');
      halt;
   end;
   if smoothsize < 0 then begin
      writeln(output,'smoothsize (=',smoothsize:5:2,')must be one positive');
      halt;
   end;

end;

procedure writeparameters(var out: text; k: char);
(* write the parameters to file out with character k at the
start of the line *)
begin
   writeln(out,k,' denplo ',version:4:2);
   writeln(out,k,' ', parameterversion:9:2,' parameterversion for denplo');
   writeln(out,k,' ', numbertoprocess:9,' numbertoprocess');

   writeln(out,k,' ', XminValue:NuWi:NuDe,
                ' ', XmaxValue:NuWi:NuDe,
                ' ', Xintervals:NuWi,
                ' Xminvalue, XmaxValue, Xintervals');

   writeln(out,k,' ', YminValue:NuWi:NuDe,
                ' ', YmaxValue:NuWi:NuDe,
                ' ', Yintervals:NuWi,
                ' Yminvalue, YmaxValue, Yintervals');

   writeln(out,k,' ', NuWi:1,
                ' ', NuDe:1,
                ' NuWi, NuDe');

   writeln(out,k,' ', Xcolumn:1, ' ', Ycolumn:1, ' Xcolumn, Ycolumn');
   writeln(out,k,' ', Xsize:NuWi:NuDe, ' ', Ysize:NuWi:NuDe,
                  ' Xsize, Ysize (cm)');
   writeln(out,k,' ', Xcorner:NuWi:NuDe, ' ', Ycorner:NuWi:NuDe,
                   ' Xcorner, Ycorner (cm)');
{
   writeln(out,k,' ', XDisplayIntervals:1, ' ', XDisplayIntervals:1,
                ' ', XDisplaySubIntervals:1, ' ', YDisplaySubIntervals:1,
                   ' XDisplayIntervals, YDisplayIntervals,',
                   ' XDisplaySubIntervals, YDisplaySubIntervals');
}
   writeln(out,k,' ', XDisplayIntervals:1, ' ', XDisplayIntervals:1,
                ' ', XDisplaySubIntervals:1, ' ', YDisplaySubIntervals:1,
                   ' (X,Y)DisplayIntervals,',
                   ' (X,Y)DisplaySubIntervals');
   writeln(out,k,' ', xwidth:1, ' ', ywidth:1, ' xwidth, ywidth');
   writeln(out,k,' ', xdecimal:1, ' ', ydecimal:1, ' xdecimal, ydecimal');
   writeln(out,k,' ', background,' ',grey:NuWi:NuDe,' background');
   writeln(out,k,' ', cutofftype,
                ' ', cutoffvaluemin:NuWi:NuDe,
                ' ', cutoffvaluemax:NuWi:NuDe,
                ' cutofftype cutoffvaluemin cutoffvaluemax');
   writeln(out,k,' ', shrinkfactorX:NuWi:NuDe,
                ' ', shrinkfactorY:NuWi:NuDe,
                ' shrinkfactorX, shrinkfactorY');

   write(out,k,' ');
   writestring(out, XaxisLabel);
   writeln(out);
   write(out,k,' ');
   writestring(out, YaxisLabel);
   writeln(out);

   writeln(out,k,' ', cog,' cog: color or grey plotting');
   writeln(out,k,' ',keyintervals:1,' keyintervals');
   writeln(out,k,' ',keytype,     ' keytype');
   writeln(out,k,' ',keyX:NuWi:NuDe,' keyX');
   writeln(out,k,' ',keyY:NuWi:NuDe,' keyY');
   writeln(out,k,' ',keyXsize:NuWi:NuDe,' keyXsize');
   writeln(out,k,' ',keyYsize:NuWi:NuDe,' keyYsize');
   writeln(out,k,' ',keyshrinkfactor:NuWi:NuDe,' keyshrinkfactor');
   writeln(out,k,' ',thelinewidth:NuWi:NuDe,' thelinewidth');
   writeln(out,k,' ',keydecimals:NuWi,' keydecimals');
   writeln(out,k,' ',keyfontsize:NuWi:NuDe,' keyfontsize');

   write  (out,k,' ',edgecontrol:1);
   write  (out,' ', edgeleft  :1:2);
   write  (out,' ', edgeright :1:2);
   write  (out,' ', edgelow   :1:2);
   write  (out,' ', edgehigh  :1:2);
   write  (out,' edgecontrol, edge(left, right, low, high)');
   writeln(out);

   write  (out,k,' ',findpeak:1);
   write  (out,' ', startX:1:NuDe);
   write  (out,' ', startY:1:NuDe);
   write  (out,' ', SearchRange:NuWi);
   writeln(out,' findpeak, startX, startY, SearchRange');
   writeln(out);

end; (* writeparameters *)

procedure buildxyplom(var l: text);
(* build the xyplom into file l *)
begin
rewrite(l);

writeln(l);
writeln(l,'/charwidth 12 def');
writeln(l,'/fontsize charwidth def');

writeln(l);
writeln(l,'% define fonts');
writeln(l,'/ffss {findfont fontsize scalefont setfont} def');
writeln(l,'/FontForStringRegular {/Times-Bold       ffss} def');
writeln(l,'/FontForStringItalic  {/Times-BoldItalic ffss} def');
writeln(l,'/FontForLogo          {/Helvetica-Bold   ffss} def');
writeln(l,'/FontForPrime         {/Symbol           ffss} def');
writeln(l,'/FontForSymbol        {/Symbol           ffss} def');

writeln(l);
writeln(l,'% make italics possible in titles');
writeln(l,'/IT {% TRstring ITstring IT -');
writeln(l,'  exch show');
writeln(l,'  FontForStringItalic');
writeln(l,'  show');
writeln(l,'  FontForStringRegular');
writeln(l,'} def');
writeln(l);

writeln(l);
writeln(l,'% make symbols possible in titles');
writeln(l,'/SY {% TRstring SYstring SY -');
writeln(l,'  exch show');
writeln(l,'  FontForSymbol');
writeln(l,'  show');
writeln(l,'  FontForStringRegular');
writeln(l,'} def');
writeln(l);

writeln(l);
writeln(l,'% make a bolder display');
writeln(l,'2 setlinewidth');
writeln(l,'%F');
writeln(l,'/Helvetica-Bold findfont');
writeln(l,'16 scalefont');
writeln(l,'setfont');
writeln(l,'%T 1.5');
writeln(l);

end;

{ccc}
(* begin module colormodule *)
procedure makecolor(cog: char; r: real;
                    var hue, saturation, brightness: real);
(* determine the color (hue, saturation, brightness)
based on cog (c, color; g, grey) and r for which
0 <= r <= 1. *)
begin
   case cog of
      'c': begin (* do color graphics *)
         (* the formula avoids red at hue = 0
            being the same as the red at hue = 1.
            see diana.p *)
{
         hue := 0.84*r + 0.16;
}
         hue := r; (* no correction, done by xyplo *)
{hhh}
         saturation := 1.0;
         brightness := 1.0;
      end;
      'g': begin (* do grey scale graphics *)
         hue := 0.0;
         saturation := 0.0;
         brightness := 1-r; (* white is low density *)
      end;
   end;
end;

procedure dobackground(
   background: char; (* background definition *)
   Xbin, Ybin: integer; (* bin location *)
   var r: real; (* b/BinMax if used *)
   var hue, saturation, brightness: real (* color definition *)
   );
(* set up the background coloring *)
(* 11. background: c = colorful (debug), g = grey, b = black, w = white *)
begin
   r := 0.0;
   case background of

      'b': begin (* black *)
         hue := 0.0;
         brightness := 0.0; (* 0 = black *)
         saturation := 0.0;
      end;

      'g': begin (* grey *)
         hue := 0.0;
         brightness := grey; (* 0.9 = grey *)
         saturation := 0.0;
      end;

      'w': begin (* white *)
         hue := 0.0;
         brightness := 1.0; (* 1.0 = white *)
         saturation := 0.0;
      end;

      'c': begin (* start with something *)

         (* just a function to distinguish various squares *)
         r := (Xbin/Xintervals + Ybin/Yintervals)/2;
         if Xbin > Ybin
         then r := 1 - r;
{hhh}
{
         hue := 0.84*r + 0.16;
}
         hue := r; (* no correction, done by xyplo *)
         saturation := 0.5;
         (* checkerboard for further distinction: *)
         if Xbin+Ybin mod 2 = 0
         then brightness := 1.0
         else brightness := 0.9;

      end;

   end;
end;

procedure determinecolors(
   cog: char; (* color or grey *)
   b: integer; (* number of data points *)
   BinMax: integer; (* BinMax number of data points *)
   var r: real; (* b/BinMax if used *)
   cutofftype: char; (* cutoff type *)
   cutoffvaluemin: real; (* cutoff value *)
   cutoffvaluemax: real; (* cutoff value *)
   Xbin, Ybin: integer; (* bin location *)
   var hue, saturation, brightness: real (* color definition *)
   );
(* determine the colors and return them as
hue, saturation, brightness.  If the number of data points
in the bin at (Xbin,Ybin) exceeds BinMax, it means that
a local peak was found that was not the global peak; this point
is higher than the recorded local peak.  So display this point
as black.
2009 Apr 11:  Note that the cutoff range changed
(from => to >) and (from =< to <)
so that background is NOT included when there is equality
to the cutoff value.  Cutoffs are inclusive. *)
procedure normal;
(* proceed with normal color determination *)
begin
   if ((cutofftype = 'f') or (cutofftype = 'F'))
      and ((b/BinMax < cutoffvaluemin) or (b/BinMax > cutoffvaluemax))
   then dobackground
        (background, Xbin, Ybin, r, hue, saturation, brightness)
   else begin
      makecolor(cog, r, hue, saturation, brightness);
   end;
end; (* normal *)
begin (* determinecolors *)
   if b = 0    (* there are no data *)
   then dobackground
        (background, Xbin, Ybin, r, hue, saturation, brightness)
   else if ((cutofftype = 'c') or (cutofftype = 'C'))
           and ((b < cutoffvaluemin) or (b > cutoffvaluemax))
   then dobackground
        (background, Xbin, Ybin, r, hue, saturation, brightness)
   else begin
      r := b/BinMax; (* density in this bin *)
      if findpeak = 'l' then begin
         if r > 1.0 then begin (* exceeded local peak, make black! *)
{
writeln(output,'     r = ',r:1:8);
writeln(output,'     b = ',b:8);
writeln(output,'BinMax = ',BinMax:8);
}
            hue := 0.0;
            brightness := 0.0;
            saturation := 0.0;
         end
         else normal
      end
      else normal
   end;
   (* White can be done faster by skipping the denploxyin
   printing: If it is a white background and b=0, don't plot. *)
end; (* determinecolors *)

procedure colorkey(var l: text;
      keytype: char;
      BinMax: integer;
      cog: char;
      keyX, keyY, keyXsize, keyYsize: real;
      keydecimals: integer;
      keyintervals: integer;
      keyshrinkfactor: real;
      keyfontsize: real;
      cutofftype: char; cutoffvaluemin, cutoffvaluemax: real;
      thelinewidth: real);
(* Make a color key in postscript to the l file.
BinMax: maximum data counts
cog: color or grey to plot.
keyX, keyY: The coordinate to place the key, in cm.
keyXsize, keyYsize: the size of the key in cm.
keydecimals: the number of decimal places to report in the key
keyintervals: number of intervals to show on the key.
keyshrinkfactor: factor by which to shrink the key boxes.
thelinewidth: width of lines around key boxes
*)
var
   hue: real;  (* conversion of r to a color *)
   saturation: real;  (* conversion of r to a color *)
   brightness: real;  (* conversion of r to a color *)
   i: integer; (* integer counter for r *)
   r: real; (* a saturation to plot *)
   Xbin, Ybin: integer; (* bin location *)
   BinValue: integer; (* bin value to report *)
   keywidth: integer; (* width of numbers based on BinMax *)
begin

   writeln(l,'% Color Key');
   writeln(l);
   Xbin := 0;
   Ybin := 0;


writeln(l,'/cmfactor 72 2.54 div def % defines points -> centimeters');
writeln(l,'/cm { cmfactor mul} def % defines centimeters');

writeln(l, '/keyX ', keyX:NuWi:NuDe, ' cm def');
writeln(l, '/keyY ', keyY:NuWi:NuDe, ' cm def');
writeln(l, '/keyXsize ', keyXsize:NuWi:NuDe, ' cm def');
writeln(l, '/keyYsize ', keyYsize:NuWi:NuDe, ' cm def');
writeln(l, '/keyshrinkfactor ', keyshrinkfactor:NuWi:NuDe, ' def');
writeln(l, '/keyfontsize ', keyfontsize:NuWi, ' def');
writeln(l, '/thelinewidth ',thelinewidth:NuWi:NuDe,' cm def');
writeln(l,'/Courier-Bold findfont keyfontsize scalefont setfont');

writeln(l, 'gsave');
writeln(l, 'keyX keyY translate');

writeln(l, '/keybox { % make the key box ');
writeln(l, 'sethsbcolor');
writeln(l, 'gsave');
writeln(l, 'keyXsize keyYsize scale');
writeln(l, '/s keyshrinkfactor def');

writeln(l, '0 0 moveto');
writeln(l, '0 s lineto');
writeln(l, 's s lineto');
writeln(l, 's 0 lineto');
writeln(l, 'closepath');
writeln(l, 'thelinewidth 0 ge {');
writeln(l, '   gsave');
writeln(l, '     thelinewidth setlinewidth');
writeln(l, '      0 setgray');
writeln(l, '      stroke');
writeln(l, '   grestore');
writeln(l, '} if');
writeln(l, 'fill');
writeln(l, 'grestore');
writeln(l, '} def');

writeln(l, '% columns: hue saturation  brightness sethsbcolor % r');

   if BinMax <= 0
   then keywidth := 1
   else begin
      if keytype = 'f' then begin (* add space for the decimals *)
         keywidth := 2 + keydecimals;
           (* The 2 covers the 0 and the decimal place. *)
      end
      else begin
         keywidth := trunc(ln(BinMax)/ln(10)+1+0.2);
           (* The 1 makes 9 be 1 decimal and 10 be two.
              The 0.2 insures no rounding effects. *)
      end;
   end;
{yyy}
   if keyintervals > BinMax then begin
      (* BinMax is zero on a flat distribution! *)
      if BinMax > 0 then keyintervals := BinMax
                    else keyintervals := 1;
   end;
   if keyintervals <= 0 then keyintervals := 1;
   for i := 0 to keyintervals do begin
      {2014 Mar 28: trunc gave 14 and 14 in hawk
      BinValue := trunc(BinMax*(i/(keyintervals)));
      }
      BinValue := trunc(BinMax*(i/(keyintervals))+0.5);

      determinecolors(
         cog,
         BinValue,
         BinMax,
         r,
         cutofftype, cutoffvaluemin, cutoffvaluemax,
         Xbin, Ybin,
         hue, saturation, brightness
         );

      (* xyplo color conversion is NOT applied to colors in xyplom, of
      course!  So it must be done NOW!  See: xyplo.p module
      pic.setcolor.  *)
      hue := 0.84 * hue + 0.16;
                       

      writeln(l, ' ',        hue:NuWi:NuDe,
                 ' ', saturation:NuWi:NuDe,
                 ' ', brightness:NuWi:NuDe,
                 ' keybox % ',i:1);

      write  (l, 'gsave');
      write  (l, ' keyXsize 0 moveto');
      write  (l, ' 0 setgray');
      case keytype of
         'n': write  (l, ' (',BinValue:keywidth,') show');
         'f': write  (l, ' (',i/(keyintervals):keywidth:keydecimals,') show');
      end;
      writeln(l, ' grestore');
      
      writeln(l, '0 keyYsize translate'); (* shift up *)

   end;

writeln(l, 'grestore');

end;


(* end module colormodule *)
{ccc}

procedure setvariables;
(* clear the variables and set up some general values *)
const
{
   debugging = false; (* whether to report debugging information *)
}
   debugging = true; (* whether to report debugging information *)
var
   X, Y: integer; (* indicies to the bin *)
begin
   writeln(output, 'clearing bins ...');
      for X := minValue to maxValue
   do for Y := minValue to maxValue
   do bins^[X,Y] := 0;

   if smoothkind <> '-' then begin (* clear smoothing bins *)
{zqsmooth}
      new(smoothingbins);
      writeln(output, 'clearing smoothing bins ...');
         for X := minValue to maxValue
      do for Y := minValue to maxValue
      do smoothingbins^[X,Y] := 0;
   end;

   for X := minValue to maxValue do begin (* clear for averaging *)
      ysums[X] := 0.0;
      ycounts[X] := 0;
   end;
   writeln(output, 'cleared');

   dataCount := 0;
   inCount := 0;
   outCount := 0;

   XinCount := 0;
   XoutCount := 0;

   YinCount := 0;
   YoutCount := 0;

   Xlo := +maxint;
   Xhi := -maxint;
   Ylo := +maxint;
   Yhi := -maxint;

   (* set variables *)

   xlessthany := (Xcolumn < Ycolumn);

   Xrange := XmaxValue - XminValue;
   Yrange := YmaxValue - YminValue;

   jumpX := Xrange / Xintervals;
   jumpY := Yrange / Yintervals;

   SymbolSizeX := shrinkfactorX*(Xrange / Xintervals);
   SymbolSizeY := shrinkfactorY*(Yrange / Yintervals);

   shiftX := +(1-shrinkfactorX)*jumpX/2;
   shiftY := +(1-shrinkfactorY)*jumpY/2;

if debugging then begin
   writeln(output, '==================================================');

   writeln(output, '  XminValue: ', XminValue:10:5, ' ');
   writeln(output, '  XmaxValue: ', XmaxValue:10:5, ' ');
   writeln(output, '     Xrange: ', Xrange:10:5, ' ');

   writeln(output, '  YminValue: ', YminValue:10:5, ' ');
   writeln(output, '  YmaxValue: ', YmaxValue:10:5, ' ');
   writeln(output, '     Yrange: ', Yrange:10:5, ' ');

   writeln(output, ' Xintervals: ', Xintervals:10, ' ');
   writeln(output, ' Yintervals: ', Yintervals:10, ' ');

   writeln(output, '      jumpX: ', jumpX:10:5, ' ');
   writeln(output, '      jumpY: ', jumpY:10:5, ' ');

   writeln(output, 'SymbolSizeX: ', SymbolSizeX:10:5, ' ');
   writeln(output, 'SymbolSizeY: ', SymbolSizeY:10:5, ' ');

   writeln(output, 'shiftX:      ', shiftX:10:5, ' ');
   writeln(output, 'shiftY:      ', shiftY:10:5, ' ');

   writeln(output, '==================================================');
end;

if SymbolSizeX <= 0 then begin
   writeln(output, 'SymbolSizeX: ', SymbolSizeX:10:5, ' is negative!');
   writeln(output, 'Program error!');
   halt
end;

if SymbolSizeY <= 0 then begin
   writeln(output, 'SymbolSizeY: ', SymbolSizeY:10:5, ' is negative!');
   writeln(output, 'Program error!');
   halt
end;
{
writeln(output, '==================================================');
}

end;

procedure grab(var Xvalue: real; Xcolumn: integer;
               var Yvalue: real; Ycolumn: integer);
(* grab the X value from column X and the Y value from column Y *)
var
   column: integer; (* count of column number *)
begin
{
copyaline(data,output);
halt;
}
   column := 1;
   while column < Xcolumn do begin
{
writeln(output,'skip x');
}
      skipcolumn(data);
      column := succ(column);
   end;
   read(data, Xvalue);
{
write  (output,'Xvalue: ',Xvalue:10:5);
writeln(output,' column: ',column:5);
}
   column := succ(column);
   while column < Ycolumn do begin
      skipcolumn(data);
      column := succ(column);
   end;
   read(data, Yvalue);
   readln(data);
{
write  (output,'Yvalue: ',Yvalue:10:5);
writeln(output,' column: ',column:5);
halt;
}
end;

procedure fillonedatum;
(* read and store one data item *)
procedure fail;
(* explain the failure *)
begin
  writeln(output,'PROGRAM ERROR:',
                 ' This datum is outside the acceptable range.');
  halt;
end; (* fail *)
begin
(*
writeln(output,'fillonedatum, dataCount =', dataCount:1);
*)
   if (data^ = '*') or (data^ = '#') or eoln(data)
   then readln(data)
   else begin
      dataCount := succ(dataCount);

      (* obtain data in either order ... *)
      if xlessthany
      then grab(Xvalue, Xcolumn, Yvalue, Ycolumn)
      else grab(Yvalue, Ycolumn, Xvalue, Xcolumn);

      if (XminValue <= Xvalue ) and
                      (Xvalue < XmaxValue)
      then begin

         (* 2004 Feb 11: the original computation was:
         Xbin := trunc( Xintervals * ((Xvalue - XminValue)/ Xrange));
         This computation would give 14.9999 that would trunc to 14,
         when the EXACT value was 15.  To solve (?) this problem
         The integer computations are now done first: *)
         Xbin := trunc((Xintervals * (Xvalue - XminValue))/ Xrange);

         (* protect the program *)
         if Xbin < minValue then begin
            writeln(output,'fillonedatum: X (= ', Xbin:1,
                           ') < minValue = ',minValue:1);
{
writeln(output,', Xbin: ',    Xbin:NuWi);
writeln(output,', minValue: ',minValue:NuWi);
writeln(output,', maxValue: ',maxValue:NuWi);
writeln(output,', Xintervals: ',Xintervals:NuWi);
writeln(output,', Xvalue: ',Xvalue:NuWi:NuDe);
writeln(output,', XminValue: ',XminValue:NuWi:NuDe);
writeln(output,', Xrange: ',Xrange:NuWi:NuDe);
}
            fail;
         end;
         if Xbin > maxValue then begin
            writeln(output,'fillonedatum: Xbin (= ',Xbin:1,
                           ') > maxValue = ',maxValue:1);
            fail;
         end;
         XinCount := succ(XinCount);
         Xin := true;
         if Xbin < Xlo then Xlo := Xbin;
         if Xbin > Xhi then Xhi := Xbin;
      end
      else begin
         Xbin := maxint;
         XoutCount := succ(XoutCount);
         Xin := false;
      end;

      if Xin then begin
         ysums[Xbin] := ysums[Xbin] + Yvalue;   (* sum for Y average *)
         ycounts[Xbin] := ycounts[Xbin] + 1; (* counts for Y average *)
      end;

      if (YminValue <= Yvalue ) and
                      (Yvalue < YmaxValue)
      then begin

         (* Original:
         Ybin := trunc( Yintervals * ((Yvalue - YminValue)/ Yrange));
         see notes for Xbin above *)
         Ybin := trunc((Yintervals * (Yvalue - YminValue))/ Yrange);

         (* protect the program *)
         if Ybin < minValue then begin
            writeln(output,'fillonedatum: Ybin (= ',Ybin:1,
                           ') < minValue = ',minValue:1);
            fail;
         end;
         if Ybin > maxValue then begin
            writeln(output,'fillonedatum: Ybin (= ',Ybin:1,
                           ') > maxValue = ',maxValue:1);
            fail;
         end;
         YinCount := succ(YinCount);
         Yin := true;
         if Ybin < Ylo then Ylo := Ybin;
         if Ybin > Yhi then Yhi := Ybin;
      end
      else begin
         Ybin := maxint;
         YoutCount := succ(YoutCount);
         Yin := false;
      end;

      if (Yin and Xin)
      then begin
         bins^[Xbin, Ybin] := succ(bins^[Xbin, Ybin]);
         inCount := succ(inCount);
      end
      else outCount := succ(outCount);
   end;
end;

procedure fillbins;
(* fill the bins with data *)
begin
   writeln(output, 'filling bins ...');
   reset(data);
   if numbertoprocess < 0
   then while not eof(data) do fillonedatum
   else while (not eof(data)) and (dataCount < numbertoprocess)
        do fillonedatum;
   writeln(output, 'filled, dataCount = ', dataCount:1);
end;

procedure findbinmax(var denploxyin: text;
                     var BinMax: integer);
(* Find the bin with the maximum number in it.
   BinMax is the maximumly filled bin.
   If findpeak is 'l' then a local maximum starting from
   (startX, startY) is searched for with the given SearchRange.
*)
const
   debugging = false; (* set to true to debug the peak finding algorithm *)
var
   b: integer; (* bin value *)
   done: boolean; (* done with local search algorithm *)
   nextX, nextY: integer; (* a possible next location *)
   newBinMax: integer; (* a possible new value for BinMax *)
   tryX, tryY: integer; (* a possible location to try *)
   step: integer; (* count of number of steps *)
   loops: integer; (* count of number of loops *)
   Xbin, Ybin: integer; (* current bin location for the data *)

procedure showlocation(X, Y: integer);
(* show the location at (X, Y) in user coordinates *)
begin
   Xvalue := Xrange*(X/Xintervals) + XminValue;
   Yvalue := Yrange*(Y/Yintervals) + YminValue;
   (* show internal coordinates only if debugging *)
   if debugging then write(output, '[',X:1,', ',Y:NuWi,']=>');
   writeln(output, 'at (',Xvalue:1:NuWi,', ',Yvalue:1:NuWi,')');
{
   writeln(output, 'at (',Xvalue:1:NuWi,', ',Yvalue:1:NuWi,') ',
           bins[X,Y]:1);
}
end; (* showlocation *)

function IsOutsideData(X, Y: integer): boolean;
(* Is the point (X, Y) outside the plot range? *)
begin
   if (X < Xlo) or
      (X > Xhi) or
      (Y < Ylo) or
      (Y > Yhi)
   then IsOutsideData := true
   else IsOutsideData := false;
end; (* IsOutsideData *)

procedure around(SearchRange: integer);
(* debug routine, show the stuff around (Xbin, Ybin) *)
var
   tryX, tryY: integer; (* a possible location to try *)
   loops: integer; (* count of number of loops *)
begin
   loops := 0;
{
   for tryX := Xbin - SearchRange to Xbin + SearchRange do begin
     for tryY := Ybin - SearchRange to Ybin + SearchRange do begin
}
   (* make the loop go so it prints x horizontal and y vertical: *)
   for tryY := Ybin + SearchRange downto Ybin - SearchRange do begin
     for tryX := Xbin - SearchRange to Xbin + SearchRange do begin
       if IsOutsideData(tryX, tryY)
       then begin
         write(output,' +');
       end
       else begin
         b := bins^[tryX, tryY];
         if (tryX = Xbin) and (tryY = Ybin)
         then write('(')
         else if (tryX-1 = Xbin) and (tryY = Ybin)
              then write(')')
              else write(' ');
         write(output,b:1);
       end;
       loops := succ(loops);
     end;
     writeln(output,' ',tryY:2);
{
     writeln(output,' ',loops:1);
}
   end;
end; (* range *)
begin (* findbinmax *)
{
findpeak := ' '; (* debug *)
findpeak := 'l'; (* debug *)
}

   if findpeak <> 'l' then begin
      writeln(output, 'Finding GLOBAL MAXIMUM for the entire plot');
      BinMax := 0;
      for Xbin := Xlo to Xhi do begin
         for Ybin := Ylo to Yhi do begin
            b := bins^[Xbin, Ybin];
            if b > BinMax then begin
               BinMax := b;
               nextX := Xbin;
               nextY := Ybin;
            end;
         end;
      end;
      Xbin := nextX;
      Ybin := nextY;
if debugging then around(2);
if debugging then showlocation(Xbin,Ybin);
   end
   else begin
      write(output, 'Finding LOCAL MAXIMUM starting from');
      writeln(output, ' (',startX:1:NuWi,', ',startY:1:NuWi,')');
      writeln(output, 'with a SearchRange of ',SearchRange :1);

{
      (* see notes in writedenploxyin on the packing function *)
      Xbin := trunc( Xintervals * ((startX - XminValue)/ Xrange));
      Ybin := trunc( Yintervals * ((startY - YminValue)/ Yrange));
}

      (* 2005 May 31: I had to add the 0.5 to make sure these matched
      the start point: *)
      Xbin := trunc( 0.5 + Xintervals * ((startX - XminValue)/ Xrange));
      Ybin := trunc( 0.5 + Yintervals * ((startY - YminValue)/ Yrange));

{
halt;
}
if debugging then showlocation(Xbin,Ybin);

      (* This should match the user startX and startY. *)
      if IsOutsideData(Xbin, Ybin) then begin
         writeln(output,'WARNING: Point (',
                        startX:1:NuWi,', ',startY:1:NuWi,
                        ') is outside the data!');
         showlocation(Xbin, Ybin);
      end
      else begin
         writeln(output,'Good, that point is inside the graph!');
      end;

      (* Core local maximum algorithm *)

      (* Let's start with something very simple *)
      (* yes, it looks at its old spot stupidly *)

if debugging then SearchRange := 3;
if debugging then around(SearchRange);

{
SearchRange := 1;
}
      BinMax := 0;
      newBinMax := 0;
      step := 0;
      loops := 0;
      done := false;
      nextX := Xbin;
      nextY := Ybin;
      showlocation(Xbin,Ybin); (* show the start again *)
      while not done do begin
         (* look at the 8 positions around the current (Xbin,Ybin)
            stupidly including (Xbin,Ybin) to find a better one *)
if debugging then writeln(output,'BinMax=',BinMax:1);
         for tryX := Xbin - SearchRange to Xbin + SearchRange do begin
            for tryY := Ybin - SearchRange to Ybin + SearchRange do begin
               if not IsOutsideData(tryX, tryY)
               then begin
                  b := bins^[tryX, tryY];
if debugging then write(output,'b=',b:1, ' ');
if debugging then showlocation(tryX,tryY);
                  if b > newBinMax then begin
                     newBinMax := b;
if debugging then writeln(output,'grab newBinMax=',newBinMax:1);
                     nextX := tryX;
                     nextY := tryY;
                     step := succ(step);
                  end;
               end
               else begin
if debugging then writeln(output,'outside of data');
if debugging then halt;
               end;
            end;
         end;

if debugging then writeln(output,'after square: newBinMax=',newBinMax:1);
if debugging then writeln(output,'after square:    BinMax=',   BinMax:1);
         if newBinMax = BinMax then begin
if debugging then writeln(output,'NO MOVE');
            (* there was no better one! *)
            done := true;
         end
         else begin
if debugging then writeln(output,'MOVE');
            (* move to the new location *)
            Xbin := nextX;
            Ybin := nextY;
            BinMax := newBinMax;
if debugging then writeln(output,'NEW location: ');
            write(output,bins^[Xbin,Ybin]:1,' ');
            showlocation(Xbin,Ybin);
         end;

if debugging then around(SearchRange);
         loops := succ(loops);
if debugging then writeln(output,'loops: ',loops:1,'------------------------------------');
if debugging then if loops > 4 then halt;

      end;
   end;

   write  (output, 'The highest number of points in a bin is: ',
                   BinMax:1,' ');
   showlocation(Xbin,Ybin);
end; (* findbinmax *)

procedure smoothsquare;
(* Smooth the data using smoothsize and smoothrange := trunc(smoothsize)
*)
const
   debugging = true;
var
   Xbin, Ybin: integer; (* current bin location for the data *)
{
   b: integer; (* bin value *)
}
   sum: integer; (* sum of surrounding bins *)
   low, high: integer; (* low and high of the range *)
   dx, dy: integer; (* change in x and change in y *)
begin
   writeln(output,'Smoothing in square with side ',smoothrange:1);
{zqsmooth}

   (* determine actual range. *)
   if trunc(smoothsize/2)*2 = round(smoothsize)  then begin
      writeln(output,'smoothsize = ',round(smoothsize):1,' is even');
      low  := -trunc(smoothsize/2) + 1;
      high := +trunc(smoothsize/2);
   end
   else begin
      writeln(output,'smoothsize = ',round(smoothsize):1,' is odd');
      low  := -trunc(smoothsize/2);
      high := +trunc(smoothsize/2);
   end;
   writeln(output, 'low =  ',low:1);
   writeln(output, 'high = ',high:1);
   (* test the range: *)
   if high - low + 1 <> round(smoothsize) then begin
      writeln(output, 'high - low + 1 =', (high - low + 1):1,
                      ' <>',
                      ' round(smoothsize) = ',
                      round(smoothsize):1);
      halt;
   end;

   (* The full array printed in the writedenplopxyin proceedure is:
      for Xbin := 0 to Xintervals - 1 do begin
         for Ybin := 0 to Yintervals - 1 do begin
      Here we reduce that by the low and high values: *)
   for Xbin := 0-low to (Xintervals - 1)-high do begin
      for Ybin := 0-low to (Yintervals - 1)-high  do begin
{
if debugging then begin
b := bins^[Xbin, Ybin];
write  (output,' Xbin: ',Xbin:2);
write  (output,' Ybin: ',Ybin:2);
write  (output,' b: ',b:NuWi);
writeln(output);
halt;
end;
}
         sum := 0;
         for dx := low to high do begin
         for dy := low to high do begin
           sum := sum + bins^[Xbin + dx, Ybin + dy];
         end;
         end;
         smoothingbins^[Xbin, Ybin] := sum;
      end;
   end;
   (* switch the bins!! *)
   swapbins := smoothingbins;
   bins := smoothingbins;
   smoothingbins := swapbins;
end;

procedure writedenploxyin(var denploxyin: text;
                          BinMax: integer; k: char;
                          yaverage: char);
(* write denploxyin and output results.
   BinMax is the maximumly filled bin.
   k is the first character to write on the line as a comment.
   If yaverage = 'y' then plot the average too.
*)
var
   b: integer; (* bin value *)
   r: real;  (* bin ratio = b/sum *)
   hue: real;  (* conversion of r to a color *)
   saturation: real;  (* conversion of r to a color *)
   brightness: real;  (* conversion of r to a color *)
   Xbin, Ybin: integer; (* current bin location for the data *)

begin
   writeln(output, 'writing denplopxyin ...');
   writeln(denploxyin);
   writeln(denploxyin,k,' definition of data columns:');

   writeln(denploxyin,k,'  1: X coordinate');
   writeln(denploxyin,k,'  2: Y coordinate');
   writeln(denploxyin,k,'  3: total counts at this position');
   writeln(denploxyin,k,'  4: (counts at this bit-position) /',
                           ' (maximum bin count)');
   writeln(denploxyin,k,'  5: X plotting coordinate');
   writeln(denploxyin,k,'  6: Y plotting coordinate');
   writeln(denploxyin,k,'  7: SymbolSizeX: X column size control');
   writeln(denploxyin,k,'  8: SymbolSizeY: Y column size control');
   writeln(denploxyin,k,'  9: hue');
   writeln(denploxyin,k,' 10: saturation');
   writeln(denploxyin,k,' 11: brightness');
   (* new as of 2022 May 04 *)
   writeln(denploxyin,k,' 12: - means y mean value, n = none');
   writeln(denploxyin,k,' 13: = means X = Y coordinate, n = not equal');

   writeln(denploxyin);
   writeln(denploxyin,k,' maximum bin count, Binmax: ',BinMax:1);
   writeln(denploxyin);

(* the packing function is:
   Xbin := trunc( Xintervals * ((Xvalue - XminValue)/ Xrange));

So for min values, it goes to zero.

   Xbinmax := trunc( Xintervals * ((XmaxValue - XminValue)/ Xrange));
   Ybinmax := trunc( Yintervals * ((YmaxValue - YminValue)/ Yrange));

but XmaxValue - XminValue is the Xrange, so this all drops out!

   Xbinmin := 0;
   Ybinmin := 0;

   Xbinmax := Xintervals;
   Ybinmax := Xintervals;
*)

   for Xbin := 0 to Xintervals - 1 do begin
      writeln(denploxyin); (* separate y lines for gnuplot.
         see gnuplot> help contour *)
      for Ybin := 0 to Yintervals - 1 do begin
         b := bins^[Xbin, Ybin];
{
write  (output,' Xbin: ',Xbin:2);
write  (output,' Ybin: ',Ybin:2);
write  (output,' b: ',b:NuWi);
writeln(output,' background: ',background:NuWi);
}

         determinecolors(
            cog,
            b,
            BinMax,
            r,
            cutofftype, cutoffvaluemin, cutoffvaluemax,
            Xbin, Ybin,
            hue, saturation, brightness
            );

         if not ( (background = 'w') and (b = 0)) then begin

{
write  (output,' Xbin: ',Xbin:2);
write  (output,' Ybin: ',Ybin:2);
write  (output,' b: ',b:NuWi);
writeln(output,' background: ',background:NuWi);
}

            (* the packing function is:
            Xbin := trunc( Xintervals * ((Xvalue - XminValue)/ Xrange));
               so unpacking should be the inverse:
            *)
            Xvalue := Xrange*(Xbin/Xintervals) + XminValue;
            Yvalue := Yrange*(Ybin/Yintervals) + YminValue;

{ttt}
{
if Xvalue > -5 then begin
   writeln(output, 'Xvalue = ', Xvalue:1);
   write  (output,' Xbin: ',Xbin:2);
   write  (output,' Ybin: ',Ybin:2);
   write  (output,' b: ',b:NuWi);
   writeln(output,' background: ',background:NuWi);
end;

if Xvalue > -4 then begin
   halt
end;
}

            write(denploxyin,
                             (Xvalue):NuWi:NuDe,        (*  1 *)
                         ' ',(Yvalue):NuWi:NuDe,        (*  2 *)
                         ' ',b:NuWi,                    (*  3 *)
                         ' ',r:NuWi:NuDe,               (*  4 *)
                         ' ',(Xvalue+shiftX):NuWi:NuDe, (*  5 *)
                         ' ',(Yvalue+shiftY):NuWi:NuDe, (*  6 *)
                         ' ',SymbolSizeX:NuWi:NuDe,     (*  7 *)
                         ' ',SymbolSizeY:NuWi:NuDe,     (*  8 *)
                         ' ',hue:NuWi:NuDe,             (*  9 *)
                         ' ',saturation:NuWi:NuDe,      (* 10 *)
                         ' ',brightness:NuWi:NuDe);     (* 11 *)
           if yaverage = 'y'
           then write(denploxyin, ' ','-')              (* 12 *)
           else write(denploxyin, ' ','n');             (* 12 new 2022 May 04*)
           if Xvalue = Yvalue
           then write(denploxyin, ' ','=')              (* 13 new 2022 May 04*)
           else write(denploxyin, ' ','n');             (* 13 new 2022 May 04*)
           writeln(denploxyin);

         end;
      end; (* end Ybin loop *)

      if yaverage = 'y' then begin (* average plot *)
{zzz}
         if ycounts[Xbin] > 0 then begin (* if there is data here ... *)
{
      writeln(output,
      XbinTEST:5,
      ' ',
      ysums[XbinTEST]:5:2,
      ' ',
      ycounts[XbinTEST]:5,
      ' ',
      ysums[XbinTEST]/ycounts[XbinTEST]:5:2
)
}
               (* write average line *)
               Yvalue := ysums[Xbin]/ycounts[Xbin];
{
writeln(output,'Yvalue=',Yvalue:5:2);
shrinkYav := 4;
shrinkYav := 8;
shrinkYav := 8;
}
               write(denploxyin,
                             (Xvalue):NuWi:NuDe,        (*  1 *)
                         ' ',(Yvalue):NuWi:NuDe,        (*  2 *)
                         ' ',b:NuWi,                    (*  3 *)
                         ' ',r:NuWi:NuDe,               (*  4 *)
                         ' ',(Xvalue+shiftX):NuWi:NuDe, (*  5 *)
               (* shift for the center of the symbol *)
                         ' ',(Yvalue -(jumpY/(2*shrinkYav))):NuWi:NuDe, (*  6 *)
                         ' ',SymbolSizeX:NuWi:NuDe,     (*  7 *)
               (* make the symbol smaller *)
                         ' ',jumpY/shrinkYav:NuWi:NuDe, (*  8 *)
                         ' ',hue:NuWi:NuDe,             (*  9 *)
                         ' ',saturation:NuWi:NuDe,      (* 10 *)
                         ' ',0.0{brightness}:NuWi:NuDe, (* 11 *)
                         ' ','+');                      (* 12 *)
         end;
      end; (* end average plot *)
 
   end; (* end Xbin loop *)
   writeln(denploxyin);
   writeln(output, 'finished writing denplopxyin');
end; (* writedenploxyin *)

(***********************************************************************)
procedure mkdenploxyplop(var denploxyplop: text);
(* make the denploxyplop file *)
var
   graphXintervals: integer; (* the X intervals for the graph *)
   graphYintervals: integer; (* the Y intervals for the graph *)
   crosshairsymbol: char; (* the symbol that defines the crosshair to use *)
   crosshairs: boolean; (* if true then cross hairs put on zero
                           of x and y *)

{
   (* strings for calling postscriptstring: *)
   instring, outstring, rawstring: string;
}

(*
from xyplo 9.01:
       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)
*)
begin

(* smart crosshairs: determine crosshairs automatically *)

if ((XminValue <= 0.0) and (XmaxValue >= 0.0)) or
   ((YminValue <= 0.0) and (YmaxValue >= 0.0))
then begin
   crosshairs := true;
   crosshairsymbol := 'c';
end
else begin
   crosshairs := false;
   crosshairsymbol := 'a';
end;

{
graphXintervals := Xintervals;
graphYintervals := Yintervals;
ddd
if   graphXintervals > MaxDisplayIntervals
then graphXintervals := MaxDisplayIntervals;
if   graphYintervals > MaxDisplayIntervals
then graphYintervals := MaxDisplayIntervals;
}

if   XDisplayIntervals > 0
then graphXintervals := XDisplayIntervals
else graphXintervals := Xintervals;

if   YDisplayIntervals > 0
then graphYintervals := YDisplayIntervals
else graphYintervals := Yintervals;

rewrite(denploxyplop);
writeln(denploxyplop,Xcorner:NuWi:NuDe,' ',Ycorner:NuWi:NuDe,' zerox zeroy');

writeln(denploxyplop,'x ',(XminValue):NuWi:NuDe,
                      ' ',(XmaxValue):NuWi:NuDe,
                        ' zx min max');
writeln(denploxyplop,'y ',(YminValue):NuWi:NuDe,
                      ' ',(YmaxValue):NuWi:NuDe,
                        ' zy min max');

writeln(denploxyplop,
                    graphXintervals:1,      ' ',graphYintervals:1,
               ' ', XDisplaySubIntervals:1, ' ',YDisplaySubIntervals:1,
               ' Display Intervals (X,Y)',
               ' Display SubIntervals (X,Y)');
(* 2005 May 29: remove the word 'Display' so that these
lines are not wraping:
               ' XDisplayIntervals YDisplayIntervals',
               ' XDisplaySubIntervals YDisplaySubIntervals ');
*)

writeln(denploxyplop,xwidth:1, ' ', ywidth:1, ' xwidth    ywidth    ');
writeln(denploxyplop,xdecimal:1, ' ', ydecimal:1, ' xdecimal  ydecimal  ');
writeln(denploxyplop,Xsize:NuWi:NuDe,' ', Ysize:NuWi:NuDe,
                     '     xsize     ysize     ');

{
writeln(denploxyplop,', data counts: ',inCount:1);
}

(* maybe implement \c to produce data counts into the strings someday ... *)

{
}
writestring(denploxyplop,XaxisLabel); writeln(denploxyplop);
writestring(denploxyplop,YaxisLabel); writeln(denploxyplop);

(* Do NOT use postscriptstring!  That takes ( and ) and puts \ infront!
2009 Apr 02
copystring(XaxisLabel, instring);
postscriptstring(instring, outstring, rawstring,
                           0.0, 0.0,
                           52);
writestring(denploxyplop, outstring); writeln(denploxyplop);

copystring(YaxisLabel, instring);
postscriptstring(instring, outstring, rawstring,
                           0.0, 0.0,
                           58);
writestring(denploxyplop, outstring); writeln(denploxyplop);
*)

{
postscriptstring(var instring, outstring, rawstring: string;
NOT USED                   var symvec: text;
NOT USED                   lowest, highest: integer;
                           rs, sd: real;
                           n: integer;
NOT USED                   havers, havesd: boolean
);
}



writeln(denploxyplop, crosshairsymbol,
                     '         zc                  ');

(* 2013 Nov 15 implement log scales
   zzzlogaxis make the n's here variable
   original hardwired code:
   writeln(denploxyplop,'n 2       zxl base            ');
   writeln(denploxyplop,'n 2       zyl base            ');
   *)
if logxscale then write(denploxyplop, 'l')
             else write(denploxyplop, 'n');
writeln(denploxyplop,' ',xbase:8:5,'      zxl base            ');
if logyscale then write(denploxyplop, 'l')
             else write(denploxyplop, 'n');
writeln(denploxyplop,' ',ybase:8:5,'      zyl base            ');

write  (denploxyplop,'          *******');
writeln(denploxyplop,          ' This xyplop was produced by denplo ',
                                version:4:2);
writeln(denploxyplop,'5 6       xcolumn   ycolumn   ');
writeln(denploxyplop,'0         symbol column       ');
writeln(denploxyplop,'7  8      xscolumn  yscolumn  ');
writeln(denploxyplop,'9 10 11   hue saturation brightness');
writeln(denploxyplop,'          ********************');
writeln(denploxyplop,'R         symbol-to-plot      ');
writeln(denploxyplop,'i         symbol-flag         ');
writeln(denploxyplop,'-1.0      symbol sizex        ');
writeln(denploxyplop,'-1.0      symbol sizey        ');
writeln(denploxyplop,'n         connection size     ');
writeln(denploxyplop,'n  0.125  linetype  size      ');
writeln(denploxyplop,'          ********************');
writeln(denploxyplop,'.');
writeln(denploxyplop,'          ********************');
{ writeln(denploxyplop,'l 0 0 0.125  User defined line'); }
writeln(denploxyplop,'          **** more parameters');

writeln(denploxyplop,
               edgecontrol,
               ' ',edgeleft:4:2,
               ' ',edgeright:4:2,
               ' ',edgelow:4:2,
               ' ',edgehigh:4:2,
               ' edgecontrol (p=page),',
               ' edge: left,right,low,high in cm');

{
writeln(denploxyplop,'p 1.50 1.50 1.50 1.50 edgecontrol (p=page),',
               ' edgeleft, edgeright, edgelow, edgehigh in cm');
}

(* this is nice, but it is not part of xyplo!!
writeln(denploxyplop,
               findpeak,
               ' ',startX:4:2,
               ' ',startY:4:2,
               ' findpeak (l=find locally),',
               ' start at (startX, startY)');
*)
{ppp}

writeln(denploxyplop,'8.91 version of xyplo that this parameter file is',
               ' designed for.');

end; (* mkdenploxyplop *)
(***********************************************************************)

procedure report(var f: text; k: char);
(* report values to file f, with k as the first character *)
begin
   writeln(f,k,' ', inCount:NuWi,' data points were INSIDE the plot');
   writeln(f,k,' ',outCount:NuWi,' data points were OUTSIDE the plot');
   writeln(f,k,' ',(outCount+inCount):NuWi, ' total data points');
{
   writeln(f,k,' Xlo: ',Xlo:NuWi);
   writeln(f,k,' Xhi: ',Xhi:NuWi);
   writeln(f,k,' Ylo: ',Ylo:NuWi);
   writeln(f,k,' Yhi: ',Yhi:NuWi);
}
end;

begin (* themain *)

   writeln(output,'denplo ',version:4:2);

   readparameters(denplop);
   writeparameters(output,' ');

   buildxyplom(denploxyplom);

   new(bins);

   (* set up smoothing *)
   if smoothkind <> '-' then begin
      smoothrange := trunc(smoothsize);
   end;

   setvariables;

   reset(data);
   fillbins;

   if smoothkind <> '-' then begin (* smoothing *)
      if smoothkind = 's' then begin (* square smoothing *)
         smoothsquare;
      end;
   end;

   rewrite(denploxyin);
   writeparameters(denploxyin, komment);
{yyy}
   (* when findpeak is 'b', BinMax is given by the user *)
   if findpeak <> 'm'
   then findbinmax(denploxyin, BinMax);
   writedenploxyin(denploxyin, BinMax, komment, yaverage);

   mkdenploxyplop(denploxyplop);

   report(output,' ');
   report(denploxyin, komment);

{
(* for testing: *)
   keyintervals := 5;
   keytype := 'f';
   keytype := 'n';
   keyX := 1.5;
   keyY := 0.0;
   keyXsize := 2.0;
   keyysize := 1.0;
   keydecimals := 1;
   keyshrinkfactor := 0.8;
   thelinewidth := 0.005;
}

   if keyintervals > 0
   then colorkey(denploxyplom,
      keytype,
      BinMax, cog,
      keyX+Xcorner, keyY+Ycorner, keyXsize, keyYsize,
      keydecimals,
      keyintervals,
      keyshrinkfactor,
      keyfontsize,
      cutofftype, cutoffvaluemin, cutoffvaluemax,
      thelinewidth);

   writeln(output,'denplo ',version:4:2, ' is DONE');

end;
(* end module denplo.themain *)

begin
   themain(data, denplop, denploxyin, denploxyplom, denploxyplop);
1: end.
