program multiscan(book, ribl, scanp, histog, multiscanp, data,
                  scanfeatures, scaninst, output);
(* multiscan: scan a Delila system book for many multipart binding sites

  modified scan written by:
  Ryan Shultzaberger

  modified by:

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

module libraries: delman, delmods
*)

label 1;  (* end of program *)

(* NOTE: code marked 'TODO' should be replaced with appropriate calls !!! *)

const
(* begin module version *)
   version = 4.56; (* of multiscan.p 2022 Oct 08
2022 Oct 08, 4.55: data file should be named multiscan, not scan!
2017 Jul 14, 4.54: namelength/linelength set to 200/210
2017 Jul 12, 4.53: namelength set to 200 to allow even longer names
2017 Jun 01, 4.52: report coodir in data file '*' lines
2017 Mar 15, 4.51: modules upgraded from prgmod: no censor
2017 Feb 22, 4.50: timep2c upgrade
2017 Feb 14, 4.49: backup
2017 Feb 14, 4.48: upgrade featurenumber
2017 Feb 14, 4.47: modules writeafeature and definition.feature.type 
                   moved to prgmod.p and returned here.
2017 Feb 05, 4.46: upgrade PatentMessage
2016 Oct 08, 4.44: clean up
2016 Oct 05, 4.43: add feature # display??? (not done so far)
2016 Oct 05, 4.42: add group # display in multiscanp
2016 Oct 05, 4.41: remove 10halt
2016 Oct 05, 4.40: tighten code slightly; add more clearmultifeatureptr calls
2016 Oct 04, 4.39: fig1Bfull.sequ.all: Segmentation fault (continue)
                   SUCCESS: new function: clearmultifeatureptr did it ...
2016 Oct 03, 4.38: fig1Bfull.sequ.all: Segmentation fault
2016 Oct 03, 4.37: clean up then revert to 4.36 for further debugging
2016 Oct 03, 4.36: fig1Bfull.sequ.nosite: Segmentation fault
                   multifeatwrite should not write when flist^.mstart=nil
2016 Sep 28, 4.35: failed for nar cover ... reversed coordinates have a bug
2016 Sep 20, 4.34: failed for nar cover ... reversed coordinates have a bug
                   if true then begin TTT nBUBBA THIS WORKED!
2016 Sep 19, 4.33: failed for nar cover ... reversed coordinates have a bug
2016 Aug 29, 4.32: why gap is slid?  Independently of coordinates, use:
                   write(scanfeatures,' ',    flist^.m^.parentevencoord:3:1);
2016 Aug 29, 4.31: fix bar ascii
2016 Aug 25, 4.30: found Gap draw bug: gap is now defined on right side like total
2016 Aug 19, 4.29: two factor group parts.
2016 Aug 18, 4.28: new feature write
2016 Aug 15, 4.26: time to thread the group together
2016 Aug 11, 4.25: look at why several gap lines are written
                   if groupnumber > 2 then halt; zzzooo
2016 Aug 10, 4.24: gap proper but several gaps written
2016 Aug 09, 4.23: gap reversal bug
2016 Aug 09, 4.22: write group mark
2016 Aug 04, 4.21: add groups
2016 Jul 26, 4.20: debug fail to make map for unshifted coo
2016 Jul 26, 4.19: remove small warnings from fpc
2016 Jul 05, 4.18: p2c version
2016 Jul 03, 4.17: 10.6.8 work
2016 Jul 03, 4.16: (work in progress)
2016 Jul 03, 4.15: clean up
2016 Jul 03, 4.14: more gpc objections Ryan found
2016 Jun 30, 4.13: gpc objections Ryan found on pietoint arg 1 type
                   make Pietointr for rounding before calling Pietoint
2016 Jun 30, 4.12: more debugging coordinates; remove non 1-n code
2016 Jun 27, 4.11: clean up
2016 Jun 27, 4.10: still excessive sites found.
2016 Jun 27, 4.09: work on segmentation fault; WEIRD setting => solved
2016 Jun 24, 4.08: work on segmentation fault; remove blanks with rb;
                   the segmentation fault is in multiloop
                   it crashes in the second time through the loop:
                   'while not totalloop do begin'
2016 Jun 24, 4.07: work on all cases of ^.coordinate - many should be
                   calling pietoint!
                   {zzzqqqkkk another coordinate problem location}
2016 Jun 24, 4.06: switch coordinate comparison logic to internal coordinates
                   {zzzqqqkkk another coordinate problem location}
2016 Jun 24, 4.05: switch coordinate comparison logic to internal coordinates
                   This fixed yjtD_thrL-Codirectional-dir-.book !!!
2016 Jun 24, 4.04: debug coordinate system problems: qqqlll
                   zzzqqqlll wrong logic! never compare piece coords!
2016 Jun 23, 4.03: debug coordinate system problems: Ri0 is not used
2016 Jun 23, 4.02: debug coordinate system problems: normalboolean is not used
2016 Jun 23, 4.01: debug coordinate system problems
2016 Jan 12, 4.00: one orderseed in time module
2015 Mar 23, 3.99: backup
2015 Mar 23, 3.98: more debug CASES
2015 Mar 23, 3.97: debug CASES
2013 Feb 07, 3.96: report actual number of 'sites found in piece' to data file
2012 Sep 21, 3.95: Correction of total bar location
@ NC_000913 6459.0 -1 "ir" " 7.3 bits"   7.347458   0.276926   0.390919
@ NC_000913 6467.0 -1 "sd" " 5.1 bits"   5.138375  -0.260022   0.397423
define "sd-(8)-ir 6459 Gap" "-" "}{" "}{" 0.0 -8
@ NC_000913 6459.0 -1 "sd-(8)-ir 6459 Gap" "2.4 bits" 2.446019 2.276413 9.344296
define "sd-ir 6459 total 10.0" "-" "||" "||"          0          8
@ NC_000913 6467 -1 "sd-ir 6459 total 10.0" "bits" 10.039814 0.000000 0.000000
For a complementary sequence the total should STILL be at the zero base!
In this case it should report the total at 6459 NOT 6467!
2007 Nov 29, 3.94: gap location was 3.5 but should have been 4.0 (TDS)
                   The locations of all ribl scans are the locations
                   of the zero base and so are always integers.  So
                   gap bars always run from integer to integer.  I
                   'defanged' the location that made non-integer
                   values for the locations.  (Note: OTHER kinds of
                   features can have non-integer locations.)  aaa
2007 Nov 28, 3.93: orgchange/chrchange fixed for scaninst (TDS)
2007 Nov 28, 3.92: fixed scaninst so it only reports for the parent site (RS)
2007 Nov 27, 3.91: scaninst name is long name from book instead of ribl
2007 Apr 25, 3.90: Fixed filter bug for sites with same -10 but
                   opposite orientation.
2007 Apr 07, 3.89: call orderseed
2007 Mar 21, 3.88: update simpson
2007 Mar 14, 3.87: improve warning for
                   RYE DID NOT CONSIDER ALL CASES
                   piecelength(apiece)-towanted = thelength
2007 Mar 14, 3.86: clean up
2007 Mar 14, 3.85: the define string also needs stringforFullName code
2007 Mar 14, 3.84: clean up stringforFullName code
2007 Mar 13, 3.83: tighten stringforFullName code
2007 Mar 13, 3.82: stringforFullName from stringforName
2007 Feb  8, 3.81: update reference
2005 Sep  1, 3.80: KNOWN BUG NOT SOLVED: in scansequence
                   endpoint is not set when
                   piecelength(apiece)-towanted = thelength
2005 Sep  1, 3.79: others and mothers must be copied in copyafeature
                   and initialized in scansequence.
2005 Sep  1, 3.78: Gap uncertainty was given in P position.
                   Change definition to be:
                      gapsurp - gap surprisal (in Ri position)
                      MinGapSurprisal - minimum gap (in P position)
                      MaxGapSurprisal - maximum gap (in Z position)
2005 Sep  1, 3.77: Ri, P and Z values were missing or wrong! TDS
                   definition needed to be copied in copyafeature!
2005 Aug  3, 3.76: cleanup
2005 Aug  3, 3.75: pick up new writeafeature from prgmod
2005 Aug  3, 3.74: why is a definition missing?  It doesn't matter:
                   allow empty definition in writeafeature!
2005 Aug  3, 3.73: remove initializestring routine
2005 Aug  3, 3.72: upgrade to new prgmod:
                   protect writeafeature from nil definition
2005 Aug  3, 3.71: Segmentation fault crash at count of seqs
2005 May 23, 3.70: add patent message
2004 Aug 12, 3.69: blank space between sets of (walker, Gap and total)
2004 Aug 11, 3.68: extra space needed in features before @ first quote
2004 Aug  4, 3.67: delmod upgrade.
2004 Aug  4, 3.66: cleanup.
2004 Aug  4, 3.65: TDS: apiece := nil was a memory leak, remove dispose!
2004 Aug  4, 3.64: RS: Fixed bug that prevented multiple pieces from being used
2004 Jan 27, 3.63: TDS: multiscan outputs Gap min, value, max for simplicity
2004 Jan 27, 3.62: RS: multiscan outputs Gap and total worst, value, best
2004 Jan 26, 3.61: RS: multiscan outputs Gap and total worst, value, best
2004 Jan 21, 3.60: TDS: upgrade modules, set getocp variables correctly
2003 Jul 30, 3.59: functional version from Ryan Shultzaberger
2003 Jul 24, 3.56  Reformatted the documentation for the program.
2003 May 20, 3.55  The program would write out the total information line
                    wrong in the scanfeatures output if the site had a negative
                    orientation.   This was fixxed.
2003 Feb  6, 3.50:  Documentation brush-up by TDS
2003 Jan 30, 3.47  New parameter file so that the matrix names are stated.
     Several checks have been added to minimize user confusion.  Changed
     scanfeatures output.
2002 Sep 23, 3.44: fixed a bug that wouldn't allow for program to complete
     if there were multiple non-essential sites.
2002 Jan 22, 3.31: scan becomes multiscan
2001 Aug 25, 3.30: spelling correction
2000 Mar 22, 3.29: fix synopsis line, add "," before 'output'. ...
2000 Jan  8, 3.28: set dnamax large so as to scan Genomes quickly
1999 Oct  6, 3.27: fix bug: scaninst names had two double quotes.
1999 Jul  9, 3.22: numbers compress together but will expand as needed.
    This means that higher digits never will be cut off!
1999 Jul  6, 3.14: upgrade to new ribl definition
1999 Jun  9: produce evencoordinate in the feature for even symmetry sites.
1999 Apr 26: fix bug so that scansequence can handle circular sequences!
1999 Apr 24: failed attempt to remove leading blanks on numbers
1999 Jan 18: add references: Schneider.Ri and Schneider.walker
1998 Aug 1: fix getocp to make scaninst results simpler
orgchange:=false; chrchange:=false; has to be called at the start.
1998 Jul 27: ranges for cutoffs implemented for scanp versions
    higher than 2.90.   For scanp using 2.90 or less, the old definitions
    will still work:

      Ri cutoff : One real on the second line is the information content at or
      above which to report in the data file.

      Z score cutoff: One real on the third line is the Z score at or below
      which to report in the data file.  A negative sign will be converted to a
      positive sign so that this parameter limits the range of acceptable sites
      to an interval on the real line.

      Probability cutoff: One real on the fourth line is the lowest probability
      which to report in the data file.  The probability of a site is
      determined from the mean and standard deviation of the Ri distribution.

1998 Jul 27: version 2.99 was only putting out one orientation on
   the scaninst; corrected in routine scansequence.
1998 May 30: names written to data file are only their length
1998 Jan 26: namelength set to 100 to allow long names
1997 Apr 13: 2.82:  Redefinition of orientation to be from coodir
   (direction (+/-) relative to genetic map) rather than piedir (direction
   (+/-) relative to coordinates) so that scans of two pieces in opposite
   orientation give the SAME results and so that lister can identify and
   remove the duplication.
1997 Apr 12: 2.75: corrected fromwanted/towanted range.  The ENTIRE Ribl is
   read in by getriblmatrix but only part is used.  The matrix is now defined
   -minribl to +maxribl.  The simplification of code will prevent bugs.
1997 Mar 28: deltaRi implemented
1997 Mar 28: parameter version control implemented (2.69)
1997 Jan  4: scan now works with single arrays of 10 million bases
1996 Nov 18: data file has "* bits" on top so alist can use it (2.68)
1996 Jun 15: moved matrix symmetry to ri.
1996 Mar 22: added scaninst and parameters to turn on and off output.
1996 Mar 22: stdev = 0 forced to give z score of Rsequence
  so the program does not blow up!
create features for input to the lister program: 1995 Dec 20
reading ribl instead of histog for mean and st.dev: 1992 sep 8
reading mean and st.dev directly from ri 1992 Sep 6
reading histog for mean and st.dev: 1992 sep 3
Limit portion of the ribl desired for the scan: 1991 May 31
Scan both strands: 1991 March 22
generalize to calculate Berg and von Hippel measure: 1990 Nov 19
generalize to accept weight matrix from Ri: 1990 Sep 26
search T7: 1988 feb 24
origin: 1988 february 24 from parse *)

(* defines lowest acceptable current multiscan parameter file *)
updateversionm = 1.03;

(* defines lowest acceptable current scan parameter file *)
updateversions = 2.70;

(* end module version *)

(* begin module describe.multiscan *)
(*

name
   multiscan: scan a Delila system book for many multipart binding sites

synopsis
   multiscan(book: in, ribl: in, scanp: in, histog: in, multiscanp: in,
   data: out, scanfeatures: out, scaninst: out, output: out);

files

   book: a book from the delila system

   ribl: a concatenated file containing all the ribls of all elements used in
         the scan.  The order of these ribls needs to match the order of
         elements in the scanp file file, and the order of elements in the
         multiscanp file.

   scanp: a concatenated file containing scanp files for all of the
          elements.  The scanp files are the same as with the scan program,
          and a detailed description of these files is in scan.p.

   multiscanp: the parameter file for multiscan.

      Parameterversion: Real;   This is the version of the paramter file, if
      it is less than the version used by the program, the program will not
      work.

      Maximum range:  Integer;   This is the range with which the sites have
      to be inside for them to be able to work together.  So if a minus
      10, a minus 35, and an Activator do not fall within this range then
      they will not be considered as a feasible combination of sites.  This
      is used so that the user can limit the spacing histograms, you may want
      to do a scan that has a smaller range than what your maximum histogram
      spacings allow, without having to change your histograms, so you change
      your maximum range to limit it.

      Total Ri cutoff: Real;  The total information for a site has to be
      above this value.  The equation of total information is given in
      greater detail in the program description.

      Filter toggle: char; This tells whether or not to do a filter based on
      the strongest parent site.  This is the first site in the list of
      sites.  So if your first site was a minus 10 site, followed by a -35,
      and activator, it would find the best combination of -35 and activators
      for the -10 and only report that combination.

   'f' means filter.  'n' means don't filter.

      Essential number cutoff: Real; So if a number is non-essential we only
      include it if it contributes positive information.  If you want your
      cutoff to be higher, you raise this value and the non-essential (info -
      GS) has to be greater than this to be included.  If it is negative it
      uses the scanp cutoff value only and doesn't subtract the gap, so if
      the scanp allows for it to be a site, then it will be included, no
      matter if it is positive or not, or if its GS is larger than it.

      Individual model string: "string", (string), char, char, char,
      char+integer; This is the string that defines how you treat each
      model.  The paramters are as follows:

   '#' if the character following the filter character (above)
      is '#', then the groupnumber is reported after the bits on each
      total line, e.g. in the form "bits-#16".  This is useful for
      debugging extremely complex maps.

 Model name: string; This is the name of the model.  This name needs
 to be the same as that in the Ribl file.  The model name needs to be
 in " ".

 Parent name: string;  The name that this model points to.  This needs
 to be the exact same name as the parent model.  The parent name needs
 to be in ( ).

 Essentiality: char;  Is it essential 'e' or non-essential 'n'.  See
 above for a better explanation of these terms.

 Orientation toggle: char;  Does the site have to have the same 's',
 different 'd', or any 'n' orientation relative to its parent.  So the
 minus 35 has to have the same orientation as the minus 10 so its
 toggle would be 's'.  But, an activator could have either orientation
 so its toggle would be 'n'.

 GS toggle: char;  It tells whether or not to include the gap surprisal.
            'g' you include the gap surprisal in the total information.
            'n' you do not include the gap surprisal.

     If the gap surprisal is turned off, then the program uses the
     furthest range in the histog that has some value greater than
     zero, and all bases in between, regardless of whether they are
     above zero.  Also it will write out the GS as 'NA' in the lister
     map, which means the GS is not counted (ie set to zero).

 Bar toggle: char; If you want to write out the gap surprsial bar then
 use the symbol 'b'.  Any other symbol will keep the GS bar from
 getting made.


 Overlap toggle and overlap difference: char and integer;  This
 controls the amount of overlapping that the site will allow.  You can
 have some overlap, no overlap, or complete overlap.
        'o' means allow overlap, no matter how much it overlaps.
        'n' means allow no overlap at all.
 'd' means allow overlap as long as the zero coordinates are further
  than the overlap distance apart.

   histog: concatenated histog files for each of the elements.  These files
      need to be in the same order as the multiscanp, scanp, and ribl files.
      Histog files need to be whole numbers and can be generated using
      genhis.

   data: The results.  Comments are lines that begin with '*'.  This file
      contains information containing all of the components in the scan.  A
      line is given for each element that gives the model name, coordinate,
      matrix orienation, Ri evaluation, gap size, and gap surprisal.  For
      each grouping the total information and the order of elements, and all
      elements used are reported.

   scanfeatures: The results in the "features" format for input to the lister
      program.  This consists of comment lines (beginning with "*"),

      definition lines (as shown above), and features of the form:
      @ U00096 3646159.0 -1 "p10" " 4.0 bits"     3.97    -0.0     0.0

      See program lister.p for details.

   scaninst:  The results are given in the form of delila instructions:
      name "dnaA"; piece K01789; get from 229 -100 to 229 +100 direction -;
      The name is the long name in the book.
      This is turned on by the last scanp parameter set in the scanp file.


   output: messages to the user

description

Multiscan searches a Delila book with an information-theory based
flexible binding site model and reports the located features.

Multiscan is the 3rd upgrade of scan. The first upgrade was to be able
to handle two part models, and was called Discan.  Discan worked by
comparing scanfeature files of two different models with some spacing
histogram to compute the total information of the site based on the
information of the two parts and the Gap surprisal of the
spacing between the two parts.  Biscan did the same computation but
did not use the scanfeatures output, instead it scanned the sequence
directly.

The biscan algorithm was to find the first site, and then based on the
spacing histogram provided, looked to see if a second site was near
by, using the same code that Scan uses.  If it found a suitable
site it would then report the site and all of the relevant information
about it.   Biscan is much faster and has more flexibility in control
than discan did, and therefore maded discan obsolete.

Discan and Biscan could only handle two part models, Multiscan can
handle an infinite number of models binding together with a large
number of controls over the binding characteristics of each model.
Multiscan can replace Biscan, but cannot replace scan (yet).

Here is a brief description of multiscan and the way that it can be
manipulated:

The models that are used by multiscan are either "essential" or
"non-essential".  Essential sites have to be present for a site to be
found (for example, in a promoter the minus 10 has to bind in order
for the polymerase to work).  Non-essential sites, are sites that do
not have be present for the system to work (for example and activator
does not have to be present for the RNAp to start initiation).  This
distinction between essential and non-essential sites is important in
understandind the multiscan algorithm.

Non-essential sites are only included if the site contributes to the
total information of the system.  That is, if some activator had an
individual information of 3 bits, but a gap surprisal value of 4 bits,
it would have an overall net nformation contribution of -1 bits to the
system, not improving the total strength of the site.  Therefore it
will not be included.

Anoter important concept is the idea of parent sites.  Since we based
the gap surprisal value on the distance between two sites, we have to
identify which model another model will bind relative to.  For
example, the minus 35 will bind relative to the minus 10 so its parent
is the minus 10.  Some activator will bind relative to the minus 35 so
its parent is the minus 35 and its grandparent is the minus 10.
Another activator may bind relative to the minus 35 so its parent is
also the minus 35 and can be thought of as the brother of the first
activator.  But like a parent, that molecule has to be bound in order
to have a child bind relative to it.  So if a parent does not find a
binding site, no children can bind.

The math for the calculating the total information of a site is the
same as presented in biscan and in the Shultaberger, Schneider
ribosome paper.  The only difference is that instead of two parts, we
now have unlimited parts.  So the math is:

Total info = Ri(UP) - GS(UP-35) + Ri(35) + Ri(10) - GS(10-35)  if
Ri(UP) - GS(UP-35) > 0 bits

else

Total info = Ri(10) + Ri(35) - GS(10-35).

The lister maps are fairly complicated now.  Since there are so many
components which could be part of a model, we used the following
techniques to help group them.  First:  All connecting bars have the
coordinate of the parent molecule.  That is if you had a p10, p35, and
Fis and p10 was the first site in your model (in the paramter file,
not spacially on the DNA), the connector bars for all the other
sites would contain its coordinate in their feature name tags.

Second: The total information line has the names of all the components
used in the model, and the order that they fall on the DNA.  So if Fis
was the most upstream site, followed by p35, and then p10 the line
would read Fis-p35-p10.

You use catted scanp files, catted ribl files, and catted histog files
as well as a multiscanp file.

examples

   an example multiscanp file:

1.03    version of multiscan that this parameter file is designed for.
200     maximum range in bases
5.0     total Ri cutoff
f       f means filter out duplicates of the features
0       essential number cutoff (Ri - gap), if negative use scanp cutoff, no gap.
"p10" (null) e s n n o3  parent, essential or non, orient tog, use GS, bar, overlap tog
"p35" (p10)  e s g n o3  parent, essential or non, orient tog, use GS, bar, overlap tog
"Fis" (p35)  e s n b o3  parent, essential or non, orient tog, use GS, bar, overlap tog
.

The order of gap surprisal numbers is now (2005 Sep 1) defined to be:
  gapsurp - gap surprisal (in Ri position)
  MinGapSurprisal - minimum gap (in P position)
  MaxGapSurprisal - maximum gap (in Z position)
This makes the computation of the sum easier.

documentation

@article{Schneider-Ri1997,
author = "T. D. Schneider
 and J. Spouge",
title = "{Information content of individual genetic sequences}",
journal = "J. Theor. Biol.",
volume = "189",
pages = "427--441",
pmid = "9446751",
note = "\htmladdnormallink
{http://dx.doi.org/10.1006/jtbi.1997.0540}
{http://dx.doi.org/10.1006/jtbi.1997.0540},
\htmladdnormallink
{https://alum.mit.edu/www/toms/papers/ri/}
{https://alum.mit.edu/www/toms/papers/ri/}",
year = "1997"}

@article{Schneider-walker1997,
author = "T. D. Schneider",
title = "Sequence Walkers:
a graphical method to display how binding proteins
interact with {DNA} or {RNA} sequences",
journal = "Nucleic Acids Res.",
volume = "25",
comment = "walker.tex, November 1, issue 21",
note = "\htmladdnormallink
{https://alum.mit.edu/www/toms/papers/walker/}
{https://alum.mit.edu/www/toms/papers/walker/},
erratum: NAR 26(4): 1135, 1998",
pages = "4408--4415",
pmid = "9336476",
pmcid = "PMC147041",
year = "1997"}

@article{Shultzaberger.Schneider2001,
author = "R. K. Shultzaberger
 and R. E. Bucheimer
 and K. E. Rudd
 and T. D. Schneider",
title = "{Anatomy of \emph{Escherichia coli}
Ribosome Binding Sites}",
journal = "J. Mol. Biol.",
volume = "313",
pages = "215--228",
pmid = "11601857",
note = "\htmladdnormallink
{https://alum.mit.edu/www/toms/papers/flexrbs/}
{https://alum.mit.edu/www/toms/papers/flexrbs/},
\htmladdnormallink
{http://dx.doi.org/10.1006/jmbi.2001.5040}
{http://dx.doi.org/10.1006/jmbi.2001.5040}",
year = "2001"}

@article{Shultzaberger.Schneider-flexprom2007,
author = "R. K. Shultzaberger
 and Zehua Chen
 and Karen A. Lewis
 and T. D. Schneider",
title = "{Anatomy of \emph{Escherichia coli} $\sigma^{70}$ promoters}",
journal = "Nucleic Acids Res.",
volume = "35",
pages = "771--788",
pmid = "17189297",
pmcid = "PMC1807945",
note = "\htmladdnormallink
{https://alum.mit.edu/www/toms/papers/flexprom/}
{https://alum.mit.edu/www/toms/papers/flexprom/}",
year = "2007"}

see also
   sites.p, ri.p, genhis.p, lister.p, dnaplot.p

   {The first published paper using the program:}
   https://alum.mit.edu/www/toms/paper/flexprom/

   {This example lister map using multiscan output
    shows multiple sequence walkers:}
https://alum.mit.edu/www/toms/papers/flexprom/icons/flexprom-fig5.jpg

author
   Ryan Kent Shultzaberger

bugs
   * Multiscan cannot handle negative coordinates!!
     (2016 Aug 30 fix in progress)

   * The large connecting bar may not be right if the series of features do
     not have the most important site as the furthest downstream, and a positive
     orientation.  You need to check this.

   * Sites with even symmetry have coordinates at position x.5, so the
     distribution has to have half base values (-0.5 0.5, 1.5, 2.5).
     I don't know if you can use them.

   * There is a problem with coupling non-essential sites.  Let me give an
     example:  You want to look at where the sigma 70 model (composed of a
     minus 10 and minus 35) binds relative to a RBS, but in your scan you do
     not want the promoter to be essential.  That is, you can identify RBSs
     and only include the promoter if it has positive information.  Since the
     promoter is a flexible thing, and you need both a p10 and p35 for the
     site to be complete, then the algorithm should only include the promoter
     if it contains both a p10, and p35.  The current algorithm does not know
     that both things need to be present and would include the p10 if it was
     good, even if the p35 wasn't.   This is a tricky concept.  My
     recommendation is to make all flexible elements essential, and that way
     you will not have any problems.

   * The rho normalization factor parameter has not been added, so the
     information contribution of any regulator is the Ri of the site.

   * It would be nice if the program could determine the range size
     based on the sizes of the individual histograms.

   * If some weird segmentation fault comes up, it may be related to freef
     being a multifeatureptr, instead of a featuretypeptr.  For some reason
     the program worked with freef as a featuretypeptr.  Try changing it back
     in themain and see what happens.

   * 2005 Sep 1: The Ri value is now properly reported for the features,
     but the P and Z values are set to zero.  This should be provided
     at some point.

   * 2016 Sep 20: multitype records use parameters: allparameters
     but should use the parameterlist instead - this would save memory!
     Those parameterlist is not even used!

technical notes
   The mean and standard deviation of the Ri distribution are stored just
   after the Ri(b,l) table in the ribl file.  They are produced automatically
   by the ri program.
*)
(* end module describe.multiscan *)

(* constants continued *)

   gapkill = 10000; (* if there are zero occurrences of a site, this will
                       set the gap so high that there is no way it could
                       be considered as a site *)

{
   WEIRD = maxint; (* weird coordinate for initiating. tds *)
}
   WEIRD = 10000000; (* weird coordinate for initiating. tds *)

(* begin module ribl.const *)
   minribl = -10001; (* lowest  ribl matrix from allowed *)
   maxribl = +10001; (* highest ribl matrix to   allowed *)
   defnegativeinfinity = -1000; (* default for negative infinity
                                   for the Ri(b,l) table *)
(* end module ribl.const version = 5.67; (@ of prgmod.p 2017 Sep 29 *)
(*]]*)

(* begin module info.const *)
   infofield = 10; (* size of field for printing information in bits *)
   infodecim = 6; (* number of decimal places for printing information *)
       (* these are used for conlist only *)
   nfield = 8; (* size of field for printing n, the number of sites *)
(* end module info.const version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module scan.const *)
   countmark = 50; (* how often to report the number of sequences scanned *)
   bvshow = false; (* if bvshow is true, then the Berg and von Hippel measure is
     calculated and reported.  It's not as useful as the individual
     information, so I will keep it off *)
(* end module scan.const *)

(* begin module interact.const *)
(* begin module string.const *)
      maxstring = 2000; (* the maximum string *)
(* end module string.const version = 5.67; (@ of prgmod.p 2017 Sep 29 *)
(* end module interact.const version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* LOCKbegin module book.const *)
(* constants needed for book manipulations *)
      dnamax = 5000000 {1024};  (* length of dna arrays *)
      namelength = 200; (* maximum key name length *)
      linelength = 210; (* maximum line readable in book *)
(* LOCKend module book.const version = 7.39; {of delmod.p 1999 Dec 13} *)

(* begin module datetime.const *)
   datetimearraylength = 19; (* length of dataarray for dates,
       It is just long enough to include the 4 digit year - solving the
       year 2000 problem:
1980/06/09 18:49:11
123456789 123456789 
         1         2
*)
(* end module datetime.const version = 1.33; (@ of timep2c.p 2017 Feb 21 *)

   maxhist = 1000; (* max size of the histogram array *)

type

(* begin module book.type *)
(* types needed for book manipulations *)

      chset = set of 'a'..'z';

      (* types defined in book definition *)

      alpha = packed array[1..namelength] of char; (* this is not alfa *)

      (* name is a left justified string with blanks following the
         characters *)
      name = record
         letters: alpha;
         length: 0..namelength (* zero means an unspecified structure *)
      end;

      lineptr = ^line;
      line = record (* a line of characters *)
         letters: packed array [1..linelength] of char;
         length: 0..linelength;
         next: lineptr
      end;

      direction = (plus, minus, dircomplement, dirhomologous);
      configuration = (linear, circular);
      state = (on, off);

      header = record (* header of key *)
         keynam: name; (* key name of structure *)
         fulnam: lineptr; (* full name of structure *)
         note:   lineptr (* note key *)
      end;

(* begin module base.type *)
   (* define the four nucleotide bases *)
   base = (a,c,g,t);
(* end module base.type version = 7.75; {of delmod.p 2013 Jun 06} *)

      (* sequence types *)
      dnaptr = ^dnastring;
      dnarange = 0..dnamax;
      seq = packed array[1..dnamax] of base;
      dnastring = record
         part: seq;
         length: dnarange;
         next: dnaptr
      end;


      orgkey = record (* organism key *)
         hea: header;
         mapunit: lineptr (* genetic map units *)
      end;

      chrkey = record (* chromosome key *)
         hea: header;
         mapbeg: real; (* number of genetic map beginning *)
         mapend: real (* number of genetic map ending *)
      end;

      pieceptr = ^piece;
      piekey = record (* piece key *)
         hea: header;
         mapbeg: real; (* genetic map beginning *)
         coocon: configuration; (* configruation (circular/linear) *)
         coodir: direction; (* direction (+/-) relative to genetic map *)
         coobeg: integer; (* beginning nucleotide *)
         cooend: integer; (* ending nucleotide *)
         piecon: configuration; (* configruation (circular/linear) *)
         piedir: direction; (* direction (+/-) relative to coordinates *)
         piebeg: integer; (* beginning nucleotide *)
         pieend: integer; (* ending nucleotide *)
      end;


      piece = record
         key: piekey;
         dna: dnaptr
      end;

      reference = record
         pienam : name;  (* name of piece referred to *)
         mapbeg : real;  (* genetic map beginning *)
         refdir : direction;  (* direction relative to coordinates *)
         refbeg : integer;  (* beginning nucleotide *)
         refend : integer;  (* ending nucleotide *)
      end;

      genkey = record  (* gene key *)
         hea : header;
         ref : reference;
      end;

      trakey = record  (* transcript key *)
         hea : header;
         ref : reference;
      end;

      markerptr = ^marker;
      markey = record (* marker key *)
         hea : header;
         ref : reference;
         sta : state;
         phenotype : lineptr;
         next : markerptr;
      end;

      marker = record
         key : markey;
         dna : dnaptr;
      end;

(* end module book.type version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module datetime.type *)
   (* array for dates *)
   datetimearray = packed array[1..datetimearraylength] of char;
(* end module datetime.type version = 1.33; (@ of timep2c.p 2017 Feb 21 *)

(* 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 *)
      previous: stringptr; (* the previous string in a series *)
      next: stringptr; (* the next string in a series *)
   end;
(* end module string.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)
(* end module interact.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module rsgra.type *)
(* begin module wave.type *)
      waveptr = ^waveparam; (* to link several wave definitions together *)
      waveparam = record (* parameters to define a cosine wave *)
         (* define a cosine wave: *)
         extreme: char; (* h or l, the high or low extreme to be defined *)
         wavelocation: real; (* the location in bases of the extreme *)
         wavebit: real; (* the location in bits of the extreme *)
         waveamplitude: real; (* the amplitude of the wave in bits *)
         wavelength: real; (* the wave length of the wave in bases *)
{ 2000 Nov 1.  This is difficult because the length of a curve is not
simply multiplied by a scale factor.
         BASEdashon: real; (* dashon interval in bases *)
         BASEdashoff: real; (* dashon interval in bases *)
         BASEdashoffset: real; (* dashon interval in bases *)
}
         dashon: real; (* the size of on dashes in cm.
                          dashon <= 0 means no dashes *)
         dashoff: real; (* the size of on dashes in cm *)
         dashoffset: real; (* the size of off dashes in cm *)
         thickness: real; (* thickness of the wave in cm.  <= 0 means default *)
         next: waveptr; (* the next wave *)
      end;
(* end module wave.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)
(* end module rsgra.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module scan.type *)
   (* define the ribl matrix.  Source library is prgmod.p
      REQUIRES ribl.const module
      REQUIRES base.type or book.type module
      REQUIRES wave.type module for wave definition! *)
   rblarray = array[a..t, minribl..maxribl] of real; (* real(B,L) *)
   ribltype = record
      riblname: string; (* name of the weight matrix *)
      riblheader: stringptr; (* misc header information, starting with '*' *)
      data: rblarray; (* real(B,L) *)
      numbers: array[a..t, minribl..maxribl] of integer; (* n(b,l) *)
      frombase, tobase: integer; (* range of the matrix *)
      mean, stdev: real; (* mean and standard deviation of distribution *)
      consensus: real; (* Ri value of consensus sequence *)
      anticonsensus: real; (* Ri value of anticonsensus sequence *)
      averageRi: real; (* average Ri for random sequence *)
      n: integer; (* number of sequences used to create matrix *)
      symmetry: char; (* the symmetry of the Riblmatrix *)
      waves: waveptr; (* cosine wave definitions for the matrix *)
      cmperbase, cmperbit: real; (* required for wave dash definition *)
      Ribound: real; (* the Ri boundary for this definition *) 
      Zbound: real; (* the Z boundary for this definition *) 
      Pbound: real; (* the probability boundary for this definition *) 
   end;
   ribltypeptr = ^ribltype; (* pointer to ribltype so it does
      not need to be created unless actually used *)
(* end module scan.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module definition.feature.type *)

   (* Features and their definitions for the lister program *)
   petaltype = record
      (* define parts of a colored rectangle around a walker *)
      el: integer; (* edgelinewidth: edge linewidth (integer) *)
      c : char;    (* color kind: how color is defined: r for RGB, h for HSB *)
      (* ---------------------------------------------- *)
      eh: real;    (* edgeh: edge hue        OR red     *)
      es: real;    (* edges: edge saturation OR blue    *)
      eb: real;    (* edgeb: edge brightness OR green   *)
      (* ---------------------------------------------- *)
      fh: real;    (* fillh: fill hue        OR red     *)
      fs: real;    (* fills: fill saturation OR blue    *)
      fb: real;    (* fillb: fill brightness OR green   *)
      (* ---------------------------------------------- *)
      spare: real;    (* spare variable *)
   end;

   definitiontypeptr = ^definetype; (* pointer to definetype *)
   definetype = record (* store a definition *)
(*
         define "Site" "-" "<v0*]" "[*0^>" -7 -3.5 0 +3.0 +5
*)
      nametag: string; (* name of the feature *)
      background: string; (* background padding characters *)
      negparts: string; (* parts of the feature display *)
      posparts: string; (* parts of the feature display *)
      locations: array [1..maxstring] of real; (* locations of the marks
         of the feature *)
      (* ---------------------------- *)
      (* computed values: *)
      marks: integer; (* the number of marks used in this feature *)
      min: real; (* one end of the feature (minimum location) *)
      max: real; (* the other end of the feature (maximum location) *)
      definitionnumber: integer; (* the number of the definition read in *)
      matrix: ribltypeptr; (* an ribl matrix if available *)
      (* The following boundaries are DIFFERENT from the originally defined
         boundaries of the matrix.  This allows the matrix to keep its bounds
         (eg for default values for scanning AND it allows the lister program
         to define separate boundaries for displaying green versus red on the
         zero bar.  *)
      RiboundFeature: real; (* the Ri boundary for this feature *)
      ZboundFeature: real; (* the Z boundary for this feature *)
      PboundFeature: real; (* the probability boundary for this feature *)
      petal: petaltype; (* color rectangle background *)
      isgap: boolean; (* the definition represents a gap between two
                         sequence walkers *)
      istotal: boolean; (* the definition represents a total feature for
                           a set of sequence walkers *)
      ismulti: boolean; (* the definition represents a sequence walker 
                           that is part of a set of sequence walkers *)
      next: definitiontypeptr; (* the next definition recorded *)
   end;

(* ---------------------------- *)
   featurecarrierptr = ^featurecarrier; (* pointer to featurecarrier *)
   featuretypeptr = ^featuretype; (* pointer to featuretype *)
   featurecarrier = record (* carries a list of features *)
      (* to make a list of features, we need a carrier that
         is not the original feature list. *)
      afeature: featuretypeptr; (* the feature being carried *)
      next: featurecarrierptr; (* the next carriage *)
   end;
   featuretype = record (* store a feature *)
(*
         @ K01789 270 +1 "Site" "Positive 270"
         @ K01789 270 -1 "Site" "Negative 270"
         @ K01789 280 +1 "Site" "Positive 280"
         @ K01789 229.0 -1 "dnaA" "+12.2 bits " 12.200338 -0.473212 0.318031
*)
      id: string; (* name of the sequence piece the feature is on *)
      featurenumber: integer; (* the number of the feature read in *)
      coordinate: real; (* location of zero base *)
      orientation: integer; (* -1 or +1 is in direction of piece *)
      nametag: string; (* the name tag - name of the feature *)
      othertag: string; (* the other tag - another tag for misc use *)
      (* ---------------------------- *)
      definition: definitiontypeptr; (* the definition of the feature *)
      (* ---------------------------- *)
      (* the following three values are only defined when there
         is a definition for the feature matrix *)
      Ri: real; (* individual information of this site *)
      Z: real; (* Z score for this site *)
      probability: real; (* probability for this site *)

      (* Three parameters replace Ri, Z and p when there is no matrix,
         Aparam, Bparam and Cparam.  If there is a matrix, Tparam
         will be blank.  Otherwise: *)
      Tparam: char; (* a parameter that defines the kind of parameter
                       that the remainder parameters are. *)
      Aparam: real; (* a parameter *)
      Bparam: real; (* a parameter *)
      Cparam: real; (* a parameter *)
      Dparam: real; (* a parameter *)
      relcoo: real; (* relative coordinate *)
         (* Since these parameters
         are now (2004 July 17) able to substitute for Ri, Z and
         probability, the Tparam cannot be placed in front of the
         three parameters in feature files.  Tparamvalues:
            h = hsb feature colors: hue, saturation, brightness
            r = rgb feature colors: red, green, blue
            b = bit range for the feature
            ' ' = defaults to sequence walker
         Summary:
         Tparam  Aparam     Bparam     Cparam     Dparam     For
          ' '    Ri         Z          P          (none)     walker
          'h'    hue        saturation brightness (none)     ascii lines
          'r'    red        green      blue       (none)     ascii lines
          'H'    hue        saturation brightness thickness  color bar
          'R'    red        green      blue       thickness  color bar
         In this scheme one can control the color of the bar as an
         ascii string OR one can control the colored rectangle behind
         the bar but NOT BOTH.
         *)
      (* ---------------------------- *)
      evencoordinate: real; (* location of zero base for asymmetric and odd
         symmetry sites, with 0.0 as the decimal.  For even symmetry sites
         it's the point of symmetry, with 0.5 as the decimal. *)
      (* ---------------------------- *)
      (* computed values: *)
      unsatisfied: boolean; (* true when the feature has not yet been used
                               during printing a line *)
      fromrange: real; (* the lower part of the feature *)
        torange: real; (* the higher part of the feature *)
      desiredline: integer; (* The line the feature wants to be printed on.
                            Zero means any value.  Line increases down the
                            page. *)
      (* There can be several features that this feature points to,
         From multiscan, Gap features have two and total features have
         at least two.  So others needs to be a list of features.
         This cannot be done on the feature list itself because that
         would wreck the feature list structure.  So the feature
         carriers are used.  Likewise, a single feature, such as
         an E. coli -35 promoter walker, can have several total features
         that point to it, representing several -10s.  So there can be
         several mother features. *)
      others:  featurecarrierptr; (* other features this feature points to! *)
      mothers: featurecarrierptr; (* The features that point to this feature *)
      (* ---------------------------- *)
      (* If this feature is duplicated, the next two pointers
      indicate the previous and next feature *)
      duplicateprevious: featuretypeptr;
{
      duplicatenext: featuretypeptr; (* making the chain go forward
        means we must search the entire list in featureinlist.
        This does not seem worthwhile yet. *)
}
{ppp}
      usedInGroup: boolean; (* feature has been used in the group *) 
      (* ---------------------------- *)
      next: featuretypeptr; (* the next feature recorded *)
   end;
(* end module definition.feature.type version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module scan.parameters *)
   allparameters = record (* all parameters *)
      parameterversion: real; (* parameter version number *)
      todo: integer; (* number of sequences to do *)
      cutDNRi: real; (* lowest  Ri value to print to data file *)
      cutUPRi: real; (* highest Ri value to print to data file *)
      cutDNZ: real; (* lowest  Z value to print to data file *)
      cutUPZ: real; (* highest Z value to print to data file *)
      cutDNP: real; (* lowest  probability print to data file *)
      cutUPP: real; (* highest probability print to data file *)
      fromwanted, towanted: integer; (* region of a ribl weight matrix
         the user requested to use for the scan *)
      ways: integer; (* 0 = program figures it out, 1 = scan sequence as is, 2 =
                        scan sequence and complement *)
      sitedefined: boolean; (* true if the site was defined *)
      sitedefinition: definitiontypeptr; (* how to write a binding site *)

      othercontrol: string; (* control data in otherquotestring *)
      othertypes: integer; (* number of things put into otherquotestring,
         as defined below; the parameter is set in readparameters *)

      numberwidth: integer; (* width of piece number in characters *)
      lengthwidth: integer; (* width of piece length in characters *)
      namewidth:   integer; (* width of piece name in characters *)
      fullnamewidth:   integer; (* width of piece fullname in characters *)
      coordinatewidth: integer; (* width of coordinate in characters *)
      orientationwidth: integer; (* width of orientation in characters *)

      Riwidth: integer; (* width of Ri in characters *)
      Ridecimal: integer; (* decimal places of Ri in characters *)

      Zwidth: integer; (* width of Z in characters *)
      Zdecimal: integer; (* decimal places of Z in characters *)

      Pdecimal: integer; (* decimal places of probability in characters *)
      Pwidth: integer; (* width of probability in characters *)

      otherquotestring: string; (* string in the feature *)

      prebuiltstring: stringptr; (* string containing pre-written information
         about the features.  These eventually get concatenated together and
         put into the f^.othertag.  Following are the locations to write
         various data into prebuiltstring *)
      stringforNumber: stringptr; (* string to write number *)
      stringforLength: stringptr; (* string to write Length *)
      stringforName: stringptr; (* string to write Name *)
      stringforFullName: stringptr; (* string to write FullName *)
      stringforCoordinate: stringptr; (* string to write coordinate *)
      stringforOrientation: stringptr; (* string to write Orientation *)
      stringforRi: stringptr; (* string to write Ri*)
      stringforZ: stringptr; (* string to write z *)
      stringforProbability: stringptr; (* string to write probability *)
      stringforOther: stringptr; (* string to write otherquotestring *)
      pbs: stringptr; (* pointer to the current prebuiltstring part *)

      (* dfi: data feature inst file output controls *)
      outdata: boolean; (* output data file *)
      outfeat: boolean; (* output scanfeatures file *)
      outinst: boolean; (* output scaninst file *)

      (* normalization controls *)
      normalizeRi: char; (* s = subtract, d = divide, other = ignore *)
      normalboolean: boolean; (* true when normalizeRi is s or d *)

      instfrom, instto: integer; (* from and to of scaninst Delila inst *)
   end;
(* end module scan.parameters *)

(* begin module discan.types *)
   box = record (* values related to distribution *)
      range: real;    (* the numbers that make up the range of the dist *)
      occur: real;    (* the number of occurrences at each position,
                         This can be a real if one is making up a
                         distribution. *)
      probdist: real; (* the log2 (prob for each postion) *)
      isitzero: boolean;    (* checks to see if the occurrence is 0 *)
   end;

   map = array[-maxhist..maxhist] of box; (* a set of boxes *)
(* end module discan.types *)

(* begin new types made by RYAN for multiscan *)
    ribllistptr = ^ribllist; (* pointer to a ribllist *)
    ribllist = record   (* a record of all the ribls *)
        matrix: ribltypeptr; (* an ribl matrix *)
        next: ribllistptr;
    end;

    parameterlistptr = ^parameterlist; (* pointer to a the parameterlist *)
    parameterlist = record   (* a record of all the parameters *)
        parameters: allparameters; (* the parameters *)
        next: parameterlistptr;
    end;

    multitypeptr = ^multitype; (* a pointer to a multitype *)
    multitype = record   (* a record of all the features for a site *)
        riblname: string; (* name of the weight matrix *)
        parentname: string; (* name of the parent weight matrix *)
        matrix: ribltypeptr; (* an ribl matrix *)
        parameters: allparameters; (* the parameters *)

        distsize: integer; (* number of positions for distribution *)
        correction: real; (* sample correction value *)
        store: map; (* disribution array *)
        lowhistog: integer;  (* low value of distribution *)
        highhistog: integer;  (* high value of distribution *)
        MinGapSurprisal: real; (* minimum gap surprisal - best spacing *)
        MaxGapSurprisal: real; (* maximum gap surprisal - worst spacing *)

        f: featuretypeptr; (* a list of features for one part of the model *)
        fstart: featuretypeptr; (* the first feature in the list of features *)
        fend: featuretypeptr; (* the last feature in the list of features *)
        fcurrent: featuretypeptr; (* the current feature in the list
                                     of features *)
        position: char; (* is it the top, mid, or end of the list *)
        posint: integer; (* the integer position of the model on the list *)
        coorddist: real; (* distance between two coords *)
        totalinfo: real; (* the total information for a site *)

        parent: integer; (* which model is this model dependent on? *)
        parentevencoord: real; (* the coordinate of the parent for the
                                  bar feature *)
        parentorient: integer; (* the orientation of the parent for the
                                  bar feature *)
        usetog: char; (* If "e" it is essential, if n you only
                         include if helpful *)
        orienttog: char; (* If "s" the models have to have the same orientation,
                            if "d" they have to have different orientations,
                            if "n" it doesn't matter. *)
        gaptog: char; (* If "g" the program considers the gap surprisal,
                         if "n" it doesn't *)

        bartog: char; (* If "b" then do write out the gap surprisal bar.
                         If "n" then do not write out the gap surprisal bar. *)

        gapsurp: real; (* the surprisal for that site *)

        oltog: char; (* overlap toggle,
                                     if 'n' no overlap is allowed,
                                     if 'd' some overlap is allowed at oldiff,
                                     if 'o' any overlap is fine.*)

        oldiff: integer; (* The space that you will allow between zero coords *)

        maxtotal: real; (* the information of the best site for colored
                           rectangles in lister *)

        foundsite: boolean; (* did you find a site in this list *)
        writeitout: boolean; (* write out this part of the model *)
        doneloop: boolean; (* are you done the loop? *)
        esscut: real; (* the essential value cutoff *)
        parentptr: multitypeptr; (* the pointer to the parent of this model *)
        location: integer; (* the location of the feature in the
                              string of features *)
        prev: multitypeptr; (* the previous item *)
        next: multitypeptr; (* the next item *)
    end;

   multifeatureptr = ^multifeature; (* a pointer to a multifeature *)
   multifeature = record   (* a record of all the features for a boxsize *)
      m: multitypeptr; (* the mlist containing all of the info for that site *)
      mstart: multitypeptr; (* the first node in the mlist *)
      totalinfo: real; (* the total information for a site *)
      writeitout: boolean; (* should you write it out *)
      next: multifeatureptr; (* a pointer to the next multifeature *)
   end;


{
   multifeatureptr = ^multifeature; (* a pointer to a multifeature *)
   multifeature = record   (* a record of all the features for a boxsize *)
      f: featuretypeptr; (* a list of features for one part of the model *)
      totalinfo: real; (* the total information for a site *)
      writeitout: boolean; (* should you write it out *)
      next: multifeatureptr; (* a pointer to the next multifeature *)
   end;
}

(* end new types made by RYAN for multiscan *)

var
(* begin module book.var *)
(* ************************************************************************ *)
(* global variables needed for book manipulations *)

      (* free storage: *)
      freeline: lineptr; (* unused lines *)
      freedna: dnaptr; (* unused dnas *)

      readnumber: boolean; (* whether to read a number from the notes, or
                              to read in the notes *)
      number: integer; (* the number of the item just read *)
      numbered: boolean; (* true when the item just read is numbered *)
      skipunnum: boolean; (* a control variable to allow skipping of
                             un-numbered items in the book *)

(* ************************************************************************ *)
(* end module book.var version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module feature.globals *)
(* globals needed by lister.p: *)
   groupnumber: integer; (* counted group number *)
{the following are not used for now:
   globaldefinitionnumber: integer; (* global counter for definitionnumber *)
   globalfeaturenumber: integer; (* global counter for featurenumber *)
}
(* end module feature.globals version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

      book,  (* a book from the sequence library *)
      ribl,  (* weight matrix *)
      scanp,  (* scanp parameters for each of the models used *)
      histog,  (* concatenated histog files *)
      multiscanp, (* overall parameters for this program *)
      data,  (* result of the scan *)
      scanfeatures, (* result of the scan as features *)
      scaninst: (* result of the scan as delila instructions *)
        text;

(* begin module package.primitive *)
(* ************************************************************************ *)
(* 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.75; {of delmod.p 2013 Jun 06} *)

(* begin module copyaline *)
procedure copyaline(var fin, fout: text);
(* copy a line from file fin to file fout *)
begin (* copyaline *)
      while not eoln(fin) do begin
         fout^ := fin^;
         put(fout);
         get(fin)
      end;
      readln(fin);
      writeln(fout);
end; (* copyaline *)
(* end module copyaline version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module copylines *)
function copylines(var fin, fout: text; n: integer): integer;
(* copy n lines of file fin to file fout.
the actual number of lines copied is returned. *)
var
      index: integer; (* the current line number *)
begin (* copylines *)
      index := 0;
      while (not eof(fin)) and (index < n) do begin
         copyaline(fin, fout);
         index := succ(index)
      end;

      copylines := index
end; (* copylines *)
(* end module copylines version = 7.75; {of delmod.p 2013 Jun 06} *)

(* ************************************************************************ *)
(* end module package.primitive version = 7.75; {of delmod.p 2013 Jun 06} *)

(* UNBUNDLE begin module package.datetime *)
(* ************************************************************************ *)
(* begin module getdatetime *)
procedure getdatetime(var adatetime: datetimearray);
(* Get the date and time into a single array from the system clock.
   adatetime contains the date: 
       1980/06/09 18:49:11
         ye mo da ho mi se
(year, month, day, hour, minute, second)
This version works after translation of the pascal by p2c
to C and then compiling with gcc. *) 
var 
   adate, atime: packed array[1..datetimearraylength] of char;
{
   adate, atime: alfa; (* ie, packed array[1..10] of char; *)
This old method won't work, since the last digit gets cut off!
}
   month: packed array[1..3] of char;
   index: integer; (* index for times *) 
begin 
(*
                1 
      12345678901
adate[13-DEC-1999]
atime[17:39:44.00]
*)

   date(adate);
   time(atime);

{
writeln(output,'br: adate[',adate,'] atime[',atime,']');
}

   (* transfer the year *)
   for index:=1 to 4 do adatetime[index]:=adate[index+7];
   adatetime[5]:='/';
   for index:=4 to 6 do month[index-3]:=adate[index];
        if month='JAN' then begin adatetime[6]:='0'; adatetime[7]:='1' end 
   else if month='FEB' then begin adatetime[6]:='0'; adatetime[7]:='2' end 
   else if month='MAR' then begin adatetime[6]:='0'; adatetime[7]:='3' end 
   else if month='APR' then begin adatetime[6]:='0'; adatetime[7]:='4' end 
   else if month='MAY' then begin adatetime[6]:='0'; adatetime[7]:='5' end 
   else if month='JUN' then begin adatetime[6]:='0'; adatetime[7]:='6' end 
   else if month='JUL' then begin adatetime[6]:='0'; adatetime[7]:='7' end 
   else if month='AUG' then begin adatetime[6]:='0'; adatetime[7]:='8' end 
   else if month='SEP' then begin adatetime[6]:='0'; adatetime[7]:='9' end 
   else if month='OCT' then begin adatetime[6]:='1'; adatetime[7]:='0' end 
   else if month='NOV' then begin adatetime[6]:='1'; adatetime[7]:='1' end 
   else if month='DEC' then begin adatetime[6]:='1'; adatetime[7]:='2' end;
   adatetime[8]:='/';
   for index:=7 to 8 do adatetime[index+2]:=adate[index-6];

   (* replace blanks with spaces in dates *)
   if adatetime[6]=' ' then adatetime[6]:='0';
   if adatetime[9]=' ' then adatetime[9]:='0';

   adatetime[9+2]:=' ';
   for index:=10 to 17 do adatetime[index+2]:=atime[index-9];
end;
(* end module getdatetime version = 1.33; (@ of timep2c.p 2017 Feb 21 *)

(* begin module readdatetime *)
procedure readdatetime (var thefile: text; var adatetime: datetimearray);
(* read the date and time from the file.  It must have this format:
123456789 123456789
         1         
1980/06/09 18:49:11
 *)
(*
2001 Jan 2: thefile is not written to!  Force 4 digit years.
2000 Oct 11: upgraded so that the p2c compiler does not object
to writing out the adatetime; added checks for the date.
*)
var
   index: integer; (* to the udatetime *)
   (* the following is an unpacked date time array, to avoid reading
   into a packed array.  reading into a packed array is not transportable *)
   udatetime: array[1..datetimearraylength] of char;
begin
   for index:=1 to datetimearraylength do read(thefile,udatetime[index]);
   pack(udatetime, 1, adatetime);
   if (adatetime[3]='/') and (adatetime[12]=':') then begin
      writeln(output,'You have an old datetime (only 2 year digits): ');
      for index:=1 to datetimearraylength do write(output,adatetime[index]);
      writeln(output);
      writeln(output,'Convert your database to 4 digit years.');
      halt;
   end;
   (* check the adatetime format.  Note that further checks for the
      other positions in the array could be done to be sure that
      they are numbers.  But this should be pretty good. *)
   if (adatetime[ 5]<>'/') or
      (adatetime[ 8]<>'/') or
      (adatetime[14]<>':') or
      (adatetime[17]<>':')
   then begin
      writeln(output,'readdatetime: bad date time read:');
      for index:=1 to datetimearraylength do write(output,adatetime[index]);
      writeln(output);
      halt
   end;
end;
(* end module readdatetime version = 1.33; (@ of timep2c.p 2017 Feb 21 *)

(* begin module writedatetime *)
procedure writedatetime(var thefile: text; adatetime: datetimearray);
(* expand the date and time out and print in the file *)
var
      index: integer; (* index of datetime *)
begin
      for index:=1 to datetimearraylength do write(thefile,adatetime[index])
end;
(* end module writedatetime version = 1.33; (@ of timep2c.p 2017 Feb 21 *)

(* begin module timeseed *)
(* Read the computer date and time.  Reverse the order of the digits and put
a decimal point in front.  This gives a fraction between zero and one that
varies quite quickly, and is always unique (if the computer has sufficient
accuracy).  It is to be used as a seed to a random number generator.  This
has the nice property that the seed changes every second and does not repeat
for thousands of years!  *)
procedure addtoseed(var seed, power: real; c: char);
(* add the digit represented by c to the seed at the power position *)
var
      n: integer; (* the character represented by c *)
begin (* addtoseed *)
      power := power/10;
{
writeln(output,'addtoseed, c = ',c);
writeln(output,'addtoseed, ord(c) = ',ord(c));
}
   n := ord(c) - ord('0');
   if (n < 0) or (n > 9) then begin
      writeln(output,'timeseed: error in datetime');
      writeln(output,'it contains "',c,'" which is not a number.');
      writeln(output,'The getdatetime routine must be fixed.');
      halt;
   end;
   seed := seed + power*n
end; (* addtoseed *)

procedure makeseed(adatetime: datetimearray; var seed: real);
(* convert adatetime to a real number in seed, reversed order
Here is the standard adatetime format:
123456789 123456789
         1         
1980/06/09 18:49:11
*)
var
   power: real; (* a digit of the seed such as 0.01 *)
begin
   seed := 0.0;
   power := 1.0;
   addtoseed(seed, power, adatetime[19]);
   addtoseed(seed, power, adatetime[18]);
   (* : *)
   addtoseed(seed, power, adatetime[16]);
   addtoseed(seed, power, adatetime[15]);
   (* : *)
   addtoseed(seed, power, adatetime[13]);
   addtoseed(seed, power, adatetime[12]);
   (*   *)
   addtoseed(seed, power, adatetime[10]);
   addtoseed(seed, power, adatetime[ 9]);
   (* / *)
   addtoseed(seed, power, adatetime[ 7]);
   addtoseed(seed, power, adatetime[ 6]);
   (* / *)
   addtoseed(seed, power, adatetime[ 4]);
   addtoseed(seed, power, adatetime[ 3]);
   addtoseed(seed, power, adatetime[ 2]);
   addtoseed(seed, power, adatetime[ 1]);
end;

procedure orderseed(adatetime: datetimearray; var seed: real);
(* convert adatetime to a real number in seed, normal order *)
var
   power: real; (* a digit of the seed such as 0.01 *)
begin
   seed := 0.0;
   power := 1.0;
   addtoseed(seed, power, adatetime[ 3]);
   addtoseed(seed, power, adatetime[ 4]);
   addtoseed(seed, power, adatetime[ 6]);
   addtoseed(seed, power, adatetime[ 7]);
   (* / *)
   addtoseed(seed, power, adatetime[ 9]);
   addtoseed(seed, power, adatetime[10]);
   (* / *)
   addtoseed(seed, power, adatetime[12]);
   addtoseed(seed, power, adatetime[13]);
   (*   *)
   addtoseed(seed, power, adatetime[15]);
   addtoseed(seed, power, adatetime[16]);
   (* : *)
   addtoseed(seed, power, adatetime[18]);
   addtoseed(seed, power, adatetime[19]);
end;

procedure timeseed(var seed: real);
(* read the computer date and time.  reverse the order of the digits
and put a decimal point in front.  this gives a fraction between
zero and one that varies quite quickly, and is always unique (if the
computer has sufficient accuracy).  it is to be used as a seed to
a random number generator. *)
var
   adatetime: datetimearray; (* a date and time *)
begin (* timeseed *)
   getdatetime(adatetime);
{
writeln(output,'timeseed: adatetime: ',adatetime);
}
   makeseed(adatetime, seed);
end; (* timeseed *)
(* end module timeseed version = 1.33; (@ of timep2c.p 2017 Feb 21 *)

(*[[*)
(* begin module limitdate *)
procedure limitdate(a,b,c,d: char; limitdatetime: datetimearray);
(* test whether the current time is before the limit.
If it is later, halt the program *)
var
   adatetime: datetimearray; (* a date and time *)
   Dday: real; (* the critical day *)
   now: real;  (* this very moment *)
begin
   getdatetime(adatetime);
{
writeln(output,'adatetime:',adatetime);
writeln(output,'adatetime[1]: ', adatetime[1]);
writeln(output,'adatetime[2]: ', adatetime[2]);
writeln(output,'adatetime[3]: ', adatetime[3]);
writeln(output,'adatetime[4]: ', adatetime[4]);
writeln(output,'adatetime[5]: ', adatetime[5]);
}

   orderseed(adatetime, now);

   if
      (limitdatetime[1] <> ' ') or
      (limitdatetime[2] <> ' ') or
      (limitdatetime[3] <> ' ') or
      (limitdatetime[4] <> ' ') then halt;

   limitdatetime[1] := a;
   limitdatetime[2] := b;
   limitdatetime[3] := c;
   limitdatetime[4] := d;

   orderseed(limitdatetime, Dday);

{
   writeln(output,'now:  ',now:20:10);
   writeln(output,'Dday: ',Dday:20:10);
}
   if now > Dday then begin
       writeln(output,'This program expired on ',limitdatetime);
       writeln(output,'See: https://alum.mit.edu/www/toms/walker/contacts.html');
       halt
   end
end;
(* end module limitdate version = 1.33; (@ of timep2c.p 2017 Feb 21 *)

(*]]*)
(* ************************************************************************ *)
(* UNBUNDLE end module package.datetime version = 7.40; {of delmod.p 2000 Feb 18} *)

(* begin module missparam *)
procedure missparam(var param: text);
(* look at param to see if the next parameter is missing
this is useful when reading in a series of parameters.  use it
just before readln of each parameter.*)
begin (* missparam *)
      if eof(param) then begin
         writeln(output,' missing parameter');
         halt
      end
end; (* missparam *)
(* end module missparam version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module package.getocp *)
(* ************************************************************************ *)
(* begin module package.brpiece *)
(* ************************************************************************ *)
(* begin module book.basis *)
(* procedures needed for book manipulations *)

(* get procedures should be used for all linked lists of records *)

procedure getline(var l: lineptr);
(* obtain a line from the free line list or by making a new one *)
begin
      if freeline<>nil
         then begin
            l:=freeline;
            freeline:=freeline^.next
         end
         else new(l);
      l^.length:=0;
      l^.next:=nil
end;

procedure getdna(var l: dnaptr);
begin
      if freedna<>nil
         then begin
            l:=freedna;
            freedna:=freedna^.next
         end
         else new(l);
      l^.length:=0;
      l^.next:=nil
end;

(* clear procedures should be called each time the records are no longer needed
         failure to do this may result in a stack overflow. *)

procedure clearline(var l: lineptr);
(* return a line to the free line list *)
var   lptr: lineptr;
begin
      if l<>nil then begin
         lptr:=l;
         l:=l^.next;
         lptr^.next:=freeline;
         freeline:=lptr
      end
end;

procedure writeline(var afile: text; l: lineptr;
                        carriagereturn: boolean);
(* write a line to a file, with carriage return if
carriagereturn is true. *)
var
   index: integer; (* index to characters in l *)
begin
   with l^ do begin
      for index := 1 to length
      do write(afile, letters[index]);
   end;
   if carriagereturn then writeln(afile);
end;

procedure showfreedna;
(* show the freedna list *)
var
   counter: integer; (* count of freedna list *)
   l: dnaptr; (* pointer into freedna list *)
begin
   l := freedna;
   counter := 0;
   while l <> nil do begin
      counter := succ(counter);
      write(output,counter:1);
      write(output, ', length = ',l^.length:1);
{
This is illegal according to gpc because one cannot
write a pointer to a text file.  It can be unearthed
for debugging.
      write(output, ', pointer id: ',l:1);
}
      writeln(output);
      l := l^.next
   end;
end;

procedure cleardna(var l: dnaptr);
(* clear the dna strutures to the free list *)
var   lptr: dnaptr;
begin
      if l<>nil then begin
         lptr:=l;
         l:=l^.next;
         lptr^.next:=freedna;
         freedna:=lptr
      end
end;

procedure clearheader(var h: header);
(* clear the header h (remove lines to free storage) *)
begin
      with h do begin
         clearline(fulnam);
         while note<>nil do clearline(note)
      end
end;

procedure clearpiece(var p: pieceptr);
(* clear the dna of the piece *)
begin
      while p^.dna<>nil do cleardna(p^.dna);
      clearheader(p^.key.hea)
end;

function chartobase(ch:char):base;
(* convert a character into a base *)
begin
      case ch of
         'a': chartobase:=a;
         'c': chartobase:=c;
         'g': chartobase:=g;
         't': chartobase:=t
      end
end;

function basetochar(ba:base):char;
(* convert a base into a character *)
begin
      case ba of
         a: basetochar:='a';
         c: basetochar:='c';
         g: basetochar:='g';
         t: basetochar:='t';
      end
end;

function complement(ba:base):base;
(* take the complement of ba *)
begin
      case ba of
         a: complement:=t;
         c: complement:=g;
         g: complement:=c;
         t: complement:=a;
      end
end;

function chomplement(b: char): char;
(* create the character complement of base b.  I must be getting hungry! *)
begin
  chomplement := basetochar(complement(chartobase(b)));
end;

function pietoint(p: integer; pie: pieceptr): integer;
(* p is a coordinate on the piece.
   we want to transform p into a number
   from 1 to n: an internal coordinate system for
   easy manipulation of piece coordinates *)
(* Note: the dirhomologous and dircomplement are treated as
plus and minus directions, which MIGHT NOT BE RIGHT! *)
var   i: integer; (* an intermediate value *)
begin
      with pie^.key do begin
         case piedir of
            dirhomologous,
            plus:  if p>=piebeg
                   then i:=p-piebeg+1
                   else i:=(p-coobeg)+(cooend-piebeg)+2;
            dircomplement,
            minus: if p<=piebeg
                   then i:=piebeg-p+1
                   else i:=(cooend-p)+(piebeg-coobeg)+2
         end;
         pietoint:=i
      end
end;

function inttopie(i: integer; pie: pieceptr):integer;
(* i is in the range 1 to some maximum.  it is an internal coordinate
   system for the program.  we want to do a
   coordinate transformation to obtain
   a value in the range of the piece called pie:
   i=1 corresponds to piebeg and
   i=its maximum corresponds to pieend *)
(* Note: the dirhomologous and dircomplement are treated as
plus and minus directions, which MIGHT NOT BE RIGHT! *)
var   p: integer; (* an intermediate value *)
begin
      with pie^.key do begin
         case piedir of
            dirhomologous,
            plus: begin
               p:=piebeg+(i-1);
               if p>cooend
                  then if coocon=circular
                     then p:=p-(cooend-coobeg+1)
            end;
            dircomplement,
            minus: begin
               p:=piebeg-(i-1);
               if p<coobeg
                  then if coocon=circular
                     then p:=p+(cooend-coobeg+1)
            end
         end;
         inttopie:=p
      end
end;

function piecelength(pie: pieceptr): integer;
(* return the length of the dna in pie *)
begin
      piecelength:=pietoint(pie^.key.pieend,pie)
end;

(* end module book.basis version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.getto *)
function getto(var thefile: text; var theline: integer; ch: chset): char;
(* search the file for a character in the first line which is a
   member of the set ch.
Note:  on 1999 March 10 the definition of this function was cleaned
up.  Instead of putting thefile on the line AFTER the charcter ch
has been found, it puts thefile ON the line.  Other routines like
brdna and brpiece have to move to the next line themselves.  This makes
getto give the OBJECT.  *)
var
   achar: char; (* a character in thefile *)
   done: boolean; (* done finding achar *)
begin
   achar:=' ';
   done := false;
   while not done do begin
      if eof(thefile)
      then done := true
      else begin
         achar := thefile^;
         if achar in ch
         then done := true
         else begin
            readln(thefile);
            theline := succ(theline)
         end
      end
   end;
   if (achar in ch) then getto:=achar
                    else getto:=' '
{ The old method -
      while (not(achar in ch)) and (not eof(thefile))
         do begin
            readln(thefile,achar);
            theline := succ(theline)
         end;
   if (achar in ch) then getto:=achar
                    else getto:=' '
}
end;
(* end module book.getto version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.skipstar *)
procedure skipstar(var thefile: text);
(* skip start of line (or star = '*'). *)
begin (* skipstar *)
   if eof(thefile) then begin
      writeln(output,' procedure skipstar: end of book found');
      halt
   end
   else begin
      if thefile^ <> '*' then begin
         writeln(output,' procedure skipstar: bad book');
         writeln(output,' "*" expected as first character on the line, but "',
                        thefile^,'" was found');
         halt
      end;
      get(thefile); (* skip the star *)
      if thefile^ <> ' ' then begin
         writeln(output,' procedure skipstar: bad book');
         writeln(output,' "* " expected on a line but "*',
                        thefile^,'" was found');
         halt
      end;
      get(thefile) (* skip the blank *)
   end
end; (* skipstar *)
(* end module book.skipstar version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brreanum *)
procedure brreanum(var thefile: text; var theline: integer; var reanum: real);
(* read a real number from the file *)
begin
      skipstar(thefile);
      readln(thefile,reanum);
      theline := succ(theline)
end;
(* end module book.brreanum version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brnumber *)
procedure brnumber(var thefile: text; var theline: integer; var num: integer);
(* read a number from the file *)
begin
      skipstar(thefile);
      readln(thefile,num);
      theline := succ(theline)
end;
(* end module book.brnumber version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brname *)
procedure brname(var thefile: text; var theline: integer; var nam: name);
(* read a name from the file *)
var   i: integer; (* an index to the name *)
      c: char; (* a character read *)
begin (* brname *)
      skipstar(thefile);
      with nam do begin
         length:=0;
         repeat
            length:=succ(length);
            read(thefile,c);
            letters[length] := c
         until (eoln(thefile)) or
               (length>=namelength) or
               (letters[length]=' ');
         if letters[length]=' ' then length:=length-1;
         if length<namelength
            then for i:=length+1 to namelength do letters[i]:=' '
      end;
      readln(thefile);
      theline := succ(theline)
end; (* brname *)
(* end module book.brname version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brline *)
procedure brline(var thefile: text; var theline: integer; var l: lineptr);
(* read a line from the file *)
var
      i: integer;
      acharacter: char;
begin
      skipstar(thefile);
      i:=0;
      (* protection added 2005 Sep 15 *)
      while (not eoln(thefile)) and (i < linelength)  do begin
         i:=succ(i);
         read(thefile,acharacter);
         l^.letters[i]:=acharacter
      end;
      (* protection added 2005 Sep 15 *)
      if not eoln(thefile) then begin
         writeln(output,'***********************************************');
         writeln(output,'* WARNING: brline: book line length exceeded');
         writeln(output,'* linelength > ',linelength:1,' characters');
         writeln(output,'* Only ',linelength:1,' characters read from book');
         writeln(output,'***********************************************');
      end;
      l^.length:=i;
      l^.next:=nil;
      readln(thefile);
      theline := succ(theline)
end;
(* end module book.brline version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brdirect *)
procedure brdirect(var thefile: text; var theline: integer; var direct: direction);
(* read a direction *)
var   ch: char;
begin
      skipstar(thefile);
      readln(thefile,ch);
      theline := succ(theline);
      if ch='+' then direct:=plus
                else direct:=minus
end;
(* end module book.brdirect version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brconfig *)
procedure brconfig(var thefile: text; var theline: integer; var config: configuration);
(* read a configuration *)
var   ch: char;
begin
      skipstar(thefile);
      readln(thefile,ch);
      theline := succ(theline);
      if ch='l' then config:=linear
                else config:=circular
end;
(* end module book.brconfig version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brnotenumber *)
procedure brnotenumber(var thefile: text; var theline: integer; var note: lineptr);
(* book note reading to obtain the number of the object.
the procedure returns the value of the number as a global.
(this is not such a good practice, but we are stuck with it for now.) *)
begin (* brnotenumber *)
      note:=nil;
      numbered := false;
      number := 0; (* force number to zero if there
         is no number at all *)
      (* the next character is n or * depending on whether there are notes *)
      if thefile^ = 'n' then begin
         readln(thefile);
         theline := succ(theline);
         if thefile^ <> 'n' then begin
            skipstar(thefile);
            if not eoln(thefile) then begin
               if thefile^ = '#' then begin
                  numbered := true;
                  get(thefile); (* move past the number symbol *)
                  read(thefile,number);
               end
            end;
            repeat
               readln(thefile);
               theline := succ(theline)
            until thefile^ = 'n';
            readln(thefile);
            theline := succ(theline)
         end
         else begin
            readln(thefile);
            theline := succ(theline)
         end
      end
end; (* brnotenumber *)
(* end module book.brnotenumber version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brnote *)
procedure brnote(var thefile: text; var theline: integer; var note: lineptr);
(* read note key *)
var
      newnote: lineptr; (* the new note *)
      previousnote: lineptr; (* the last line of the notes *)
begin (* brnote *)
      note:=nil;
      if thefile^ = 'n' then begin (* enter note *)
         readln(thefile);
         theline := succ(theline);
         if thefile^ <> 'n' then begin (* abort null note (n/n) *)
            getline(note);
            newnote:=note;
            while thefile^ <> 'n' do begin (* wait until end of note *)
               brline(thefile,theline,newnote);
               previousnote:=newnote;
               (* get next note *)
               getline(newnote^.next);
               newnote:=newnote^.next;
            end;
            (* last note was not used, so: *)
            clearline(newnote);
            previousnote^.next:=nil;
            readln(thefile);
            theline := succ(theline);
         end
         else begin
            readln(thefile);
            theline := succ(theline);
         end;
      end
end; (* brnote *)
(* end module book.brnote version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brheader *)
procedure brheader(var thefile: text; var theline: integer; var hea: header);
(* read the header of a key. *)
begin
      with hea do begin
         readln(thefile); (* move past the object name - new definition 1999 Mar 13 *)
         theline := succ(theline);
{bbb}
         (* read key name *)
         brname(thefile,theline,keynam);

         (* read full name *)
         getline(fulnam);
         brline(thefile,theline,fulnam);

         (* read note key *)
         if readnumber then brnotenumber(thefile,theline,note)
                       else brnote(thefile,theline,note)
      end
end;
(* end module book.brheader version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.copyheader *)
procedure copyheader(fromhea: header; var tohea: header);
(* copy the header fromhea into tohea.  Note that the
linked objects are NOT copied, but merely pointed to. *)
begin
   tohea.keynam.letters := fromhea.keynam.letters;
   tohea.keynam.length  := fromhea.keynam.length;
   tohea.note           := fromhea.note;
   tohea.fulnam         := fromhea.fulnam;
end;
(* end module book.copyheader version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brpiekey *)
procedure brpiekey(var thefile: text; var theline: integer; var pie: piekey);
(* read piece key, track the line number *)
begin
      with pie do begin
         brheader(thefile,theline,hea);
         brreanum(thefile,theline,mapbeg);
         brconfig(thefile,theline,coocon);
         brdirect(thefile,theline,coodir);
         brnumber(thefile,theline,coobeg);
         brnumber(thefile,theline,cooend);
         brconfig(thefile,theline,piecon);
         brdirect(thefile,theline,piedir);
         brnumber(thefile,theline,piebeg);
         brnumber(thefile,theline,pieend);
      end
end;
(* end module book.brpiekey version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brdna *)
procedure brdna(var thefile: text; var theline: integer; var dna: dnaptr);
(* read in dna from thefile, track the line *)
(* note: if the dna were circularized, by linking the last dnastring
to the first, then the cleardna routine could not clear properly,
and would loop forever... there is no reason to do that, since a simple
mod function will allow one to access the circle. *)
var
      ch: char;
      workdna: dnaptr;
begin
      getdna(dna);
      workdna:=dna;
      ch:=getto(thefile,theline,['d']);
      readln(thefile);
      theline := succ(theline);
      read(thefile,ch); (* skipstar *)
      while (ch = '*') do
      begin
         read(thefile,ch); (* skip blank *)
         repeat
            read(thefile,ch);
            if ch in ['a','c','g','t'] then begin
               if workdna^.length=dnamax then begin
                  getdna(workdna^.next);
                  workdna:=workdna^.next
               end;
               workdna^.length:=succ(workdna^.length);
               workdna^.part[workdna^.length]:=chartobase(ch)
            end
         until eoln(thefile);
         readln(thefile); (* go to next line *)
         theline := succ(theline);
         read(thefile,ch); (* ch is either '*' or 'd' *)
      end;
      readln(thefile); (* read past the d *)
      theline := succ(theline);
end;
(* end module book.brdna version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brpiece *)
procedure brpiece(var thefile: text; var theline: integer; var pie: pieceptr);
(* read in a piece, change theline to reflect the lines traversed *)
begin
{
   readln(thefile); (* move past the word 'piece' - new definition 1999 Mar 13 *)
   theline := succ(theline); (* BUG: was below! *)
bbb}
   brpiekey(thefile,theline,pie^.key);
   if numbered or (not skipunnum)
   then brdna(thefile,theline,pie^.dna);
   readln(thefile); (* move past the word 'piece' - new definition 1999 Mar 13 *)
   theline := succ(theline);
end;
(* end module book.brpiece version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brinit *)
procedure brinit(var book: text; var theline: integer);
(* check that the book is ok to read, and
   set up the global variables for br routines *)
begin (* brinit *)
      (* halt if the book is bad (first word is 'halt') or the first
         character is not * *)
      reset(book);
      if not eof(book) then begin
         (* check for the date line *)
         if book^ <> '*' then begin
            if book^ <> 'h'
            then writeln(output, ' this is not the first line of a book:')
            else writeln(output, ' bad book:');
            write(output, ' ');

            while not (eoln(book) or eof(book)) do begin
               write(output, book^);
               get(book)
            end;
            writeln(output);
            halt
         end
      end
      else begin
         writeln(output, ' book is empty');
         halt
      end;

      (* initialize free storage *)
      freeline:=nil;
      freedna:=nil;

      readnumber:=true; (* usually we read in numbers for items *)
      number:=0; (* arbitrary value *)
      numbered:=false; (* the piece has no number (none yet read in) *)
      skipunnum:=false;
   theline := 1;
end; (* brinit *)
(* end module book.brinit version = 7.75; {of delmod.p 2013 Jun 06} *)

(* ************************************************************************ *)
(* end module package.brpiece version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brorgkey *)
procedure brorgkey(var thefile: text; var theline: integer; var org: orgkey);
(* read organism key *)
begin
      with org do begin
{bbb}
         brheader(thefile,theline,hea);
         getline(mapunit);
         brline(thefile,theline,mapunit);
      end
end;
(* end module book.brorgkey version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.brchrkey *)
procedure brchrkey(var thefile: text; var theline: integer; var chr: chrkey);
(* read chromosome key *)
begin
      with chr do begin
{bbb}
         brheader(thefile,theline,hea);
         brreanum(thefile,theline,mapbeg);
         brreanum(thefile,theline,mapend);
      end
end;
(* end module book.brchrkey version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.getocp *)
procedure getocp(var thefile: text; var theline: integer;
                 var org: orgkey;   var orgchange, orgopen: boolean;
                 var chr: chrkey;   var chrchange, chropen: boolean;
                 var pie: pieceptr; var piechange, pieopen: boolean);
(* Get the next piece and its organism and chromosome keys.  The three change
variables indicate whether or not a new organism, chromosome or piece name
was found.  If a piece is not found, then pieopen will be false.  orgopen,
chropen and pieopen are used by getocp to tell when it has entered an
organism, chromosome or piece.  All booleans should be set to false
initially.  There should be one triplet for each book read.

It is important to initialize ALL variables, including pie:

   orgchange := false; orgopen := false;
   chrchange := false; chropen := false;
   piechange := false; pieopen := false;
   pie := nil;
   theline := 0;

1999 June 2  The book reading routines now treat data objects more
precisely.  Rather than test for eof, the endo of book occurs when pieopen is
returned as false.  A book reading loop now looks like this:

   repeat
      getocp(book, theline,
                   org, orgchange, orgopen,
                   chr, chrchange, chropen,
                   pie, piechange, pieopen);
      writeln(output,'pieopen: ',pieopen);
      if pieopen then begin
         writeln(output,'piece at line: ',theline:1);
      end;
   until not pieopen;
*)

var
      ch: char;
      newchr: chrkey;
      neworg: orgkey;
      newpie: pieceptr;
begin
      ch:='a';
      while not (ch in [' ','p']) do begin
         ch:=getto(thefile,theline,['o','c','p']);
         if ch <> ' ' then begin
            case ch of
               'o': if orgopen
                  then begin
                     readln(thefile);
                     (* move past the word 'organism'
                        - new definition 1999 Mar 13 *)
                     orgopen:=false (* close organism *)
                  end
                  else begin
                     brorgkey(thefile,theline,neworg);
                     if (neworg.hea.keynam.letters <> org.hea.keynam.letters) and
                        (neworg.hea.keynam.length <> org.hea.keynam.length)
                     then begin
{
writeln(output,'--------orgchanged!');
write  (output,'--------old org:"',   org.hea.keynam.letters);
writeln(output,              '" ',    org.hea.keynam.length:1);
write  (output,'--------new org:"',neworg.hea.keynam.letters);
writeln(output,              '" ',neworg.hea.keynam.length:1);
}
(*ccc*)
                        orgchange:=true;
                        copyheader(neworg.hea,org.hea);
                        (* move the mapunit over to the org! *)
                        org.mapunit := neworg.mapunit;
                        clearline(neworg.mapunit);

                     end
                     else orgchange:=false;
                     orgopen:=true;
                  end;
               'c': if chropen
                  then begin
                     readln(thefile);
                     (* move past the word 'chromosome'
                        - new definition 1999 Mar 13 *)
                     chropen:=false (* close chromosome *)
                  end
                  else begin
                     brchrkey(thefile,theline,newchr);
                     if (newchr.hea.keynam.letters <> chr.hea.keynam.letters) and
                        (newchr.hea.keynam.length <> chr.hea.keynam.length)
                     then begin
{
writeln(output,'--------chrchanged!');
write  (output,'--------old chr:"',   chr.hea.keynam.letters);
writeln(output,              '" ',    chr.hea.keynam.length:1);
write  (output,'--------new chr:"',newchr.hea.keynam.letters);
writeln(output,              '" ',newchr.hea.keynam.length:1);
}
                        chrchange:=true;
                        copyheader(newchr.hea,chr.hea);
                        (* move the map range over to the chr! *)
                        chr.mapbeg := newchr.mapbeg;
                        chr.mapend := newchr.mapend;
                     end
                     else chrchange:=false;
                     chropen:=true;
                  end;
               'p': if pieopen
                  then begin
                     pieopen:=false; (* close last piece *)
                     ch:='a' (* prevent falling out of the loop *)
                  end
                  else begin
                     new(newpie);
                     brpiece(thefile,theline,newpie);
                     if pie = nil
                     then piechange := true
                     else begin
                        if (newpie^.key.hea.keynam.letters
                         <> pie^.key.hea.keynam.letters) and
                           (newpie^.key.hea.keynam.length
                         <> pie^.key.hea.keynam.length)
                        then begin
                           piechange:=true;
                        end
                        else piechange:=false;
                     end;
                     pieopen:=true;
                     (* we always have to switch over to the new piece,
                        because although the name may be the same,
                        the DNA sequence could be different.  That is,
                        the book may contain two pieces with the same
                        name, and we want to be sure to search the
                        new one, not the old one. *)
                     if pie <> nil then begin
                        clearpiece(pie); (* save the links *)
                        dispose(pie); (* close up shop *)
                     end;
                     pie := newpie;
                  end
            end
         end
         else begin
            pieopen := false
         end
      end
end;
(* origin: search version = 6.39 *)
(* end module book.getocp version = 7.75; {of delmod.p 2013 Jun 06} *)

(* ************************************************************************ *)
(* end module package.getocp version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.bwbasics *)

procedure bwstartline(var book: text);
(* start a line of output to the book *)
begin
      write(book,'* ')
end;

procedure bwline(var book: text; l: lineptr);
(* write a line to the book *)
var   i: integer;
begin
   if l<>nil then begin
      bwstartline(book);
      if l^.length<>0
      then for i:=1 to l^.length do write(book,l^.letters[i]);
      writeln(book);
   end
end;

procedure bwtext(var book: text; lines: lineptr);
(* write a set of lines to the book *)
var
   l: lineptr;
begin
   l := lines;
   while l<>nil do begin
      bwline(book,l);
      l := l^.next;
   end;
end;

procedure bwnote(var book: text; note: lineptr);
(* writes the notes pointed to by 'note' to 'book' *)
begin
   if note<>nil then begin
      writeln(book,'note');
      bwtext(book,note);
      writeln(book,'note');
   end
end;

procedure bwnumber(var book: text; num: integer);
(* write a number to the book *)
begin
   bwstartline(book);
   writeln(book,num:1) (* pascal will expand the field as far as needed *)
end;

procedure bwreanum(var book: text; reanum: real);
(* write a real number to the book *)
begin
   bwstartline(book);
   writeln(book,reanum:1:2) (* pascal will expand the field *)
end;

procedure bwstate(var book: text; sta: state);
(* write a state to the book *)
begin
   bwstartline(book);
   case sta of
      on:  writeln(book,'on');
      off: writeln(book,'off')
   end
end;

procedure bwname(var book: text; nam: name);
(* write a name to the book *)
var   i: integer;
begin
   with nam do begin
      bwstartline(book);
write(book,'THINGY-"');
      for i:=1 to length do write(book,letters[i]);
write(book,'"-BANGY');
      writeln(book)
   end
end;

procedure bwdirect(var book: text; direct: direction);
(* write a direction to the book *)
begin
   bwstartline(book);
   case direct of
                  dirhomologous, (* handle case, may not be right *)
      plus:  writeln(book,'+');
                  dircomplement, (* handle case, may not be right *)
      minus: writeln(book,'-')
   end
end;

procedure bwconfig(var book: text; config: configuration);
(* write a configuration to the book *)
begin
   bwstartline(book);
   case config of
      linear:   writeln(book,'linear');
      circular: writeln(book,'circular')
   end
end;

procedure bwheader(var book: text; hea: header);
(* write a key header to the book *)
begin
   with hea do begin
      bwname(book,keynam);
      bwline(book,fulnam);
      bwnote(book,note);
   end
end;

procedure bworgkey(var book: text; org: orgkey);
(* write the organism key *)
begin
   with org do begin
      bwheader(book,hea);
      bwline  (book,mapunit)
   end
end;

procedure bwchrkey(var book: text; chr: chrkey);
(* write the chromosome key *)
begin
   with chr do begin
      bwheader(book,hea);
      bwreanum(book,mapbeg);
      bwreanum(book,mapend)
   end
end;

(* end module book.bwbasics version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.bwdna *)
procedure bwdna(var thefile: text; d: dnaptr);
(* write the dna pointed to by 'd' to 'thefile' *)
var i: integer; (* index to the sequence *)
    l: integer; (* index to the number of bases on the line *)
    newline: boolean; (* true when a new line should be started *)
begin
   writeln(thefile,'dna');
   newline := true;
   while (d <> nil) do begin
      for i := 1 to d^.length do begin
         if newline then begin
            bwstartline(thefile);
            l := 0;
            newline := false
         end;
         write(thefile,basetochar(d^.part[i]));
         l := succ(l);
         if (l mod 60 = 0) or  (* end of line *)
            ((i = d^.length) and (d^.next = nil)) (* last base *)
         then begin
            writeln(thefile);
            newline := true
         end;
      end;
      d := d^.next;
   end;
   if not newline then writeln(thefile); (* last base *)
   writeln(thefile,'dna');
end;
(* end module book.bwdna version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.bwpie *)
procedure bwpie(var thefile: text; pie: pieceptr);
(* writes the information pointed to by 'pie' to 'thefile' *)
begin
   writeln(thefile,'piece');
   with pie^ do begin
      bwheader(thefile,key.hea);

      bwstartline(thefile);
      writeln(thefile,key.mapbeg:1:2);

      bwstartline(thefile);
      if (key.coocon = circular) then writeln(thefile,'circular')
            else writeln(thefile,'linear');

      bwstartline(thefile);
      if (key.coodir = plus) then writeln(thefile,'+')
            else writeln(thefile,'-');

      bwstartline(thefile);
      writeln(thefile,key.coobeg:1);

      bwstartline(thefile);
      writeln(thefile,key.cooend:1);

      bwstartline(thefile);
      if (key.piecon = circular) then writeln(thefile,'circular')
            else writeln(thefile,'linear');

      bwstartline(thefile);
      if (key.piedir = plus) then writeln(thefile,'+')
            else writeln(thefile,'-');

      bwstartline(thefile);
      writeln(thefile,key.piebeg:1);

      bwstartline(thefile);
      writeln(thefile,key.pieend:1);
   end;
   bwdna(thefile,pie^.dna);
   writeln(thefile,'piece');
end;
(* end module book.bwpie version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module book.getbase *)
function getbase(position: integer; pie: pieceptr):base;
(* Get a base from the position (internal coordinates) of the
   piece.  Protection is made against positions outside the piece.
In the case of circles it would be convenient to wrap around when
requests are off the end.  So the routine will do a modular
wrap for positions outside the range 1 to the length.
This is a new feature as of 2000 March 22.
*)
var
   workdna: dnaptr; (* pointer to the dna part of pie *)
   p: integer; (* current count of bases into the workdna *)
   spot: integer; (* the last base of the dna part *)
   thelength: integer; (* the length of the piece *)
begin
{
writeln(output,'NEW getbase: position=',position:1,'^^^^^^^^^^^^^^^^^^^^');
}
   (* handle cases of position out of range by circular wrapping *)
   thelength := piecelength(pie);
   while position < 1         do position := position + thelength;
   while position > thelength do position := position - thelength;

   workdna:=pie^.dna;
   p:=workdna^.length;
   while position > p do begin
{
writeln(output,' workdna^.length=',workdna^.length:1);
}
      workdna := workdna^.next;
      if workdna = nil then begin
        writeln(output,'error in function getbase!');
        halt
      end;
      p := p + workdna^.length;
   end;
{
writeln(output,'p=',p:1);
}
   if workdna = nil then begin
      writeln(output,'error in getbase: request off end of piece');
      halt
   end
   else begin
      spot := workdna^.length - (p-position);
{
writeln(output,'spot=',spot:1);
showdnasegment(output,workdna, spot);
}
      if (spot <= 0) then begin
         writeln(output,'error in getbase, spot (= ',spot:1, ') must be positive');
         halt
      end;
      if (spot > workdna^.length) then begin
         writeln(output,'error in getbase, spot (=',spot:1,
                        ') must be less than length (=',workdna^.length:1,')');
         halt
      end;
{
writeln(output,'base = ', workdna^.part[spot]);
}
      getbase:=workdna^.part[spot]
   end
end;
(* end module book.getbase version = 7.75; {of delmod.p 2013 Jun 06} *)

(* begin module skipblanks *)

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

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

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

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

procedure skipcolumn(var thefile: text);
(* skip over a data column *)
begin
   skipblanks(thefile); skipnonblanks(thefile)
end;
(* end module skipblanks version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module clearstring *)
(* These modules clear strings in various ways *)
(* ---- *)
procedure emptystring(var ribbon: string);
(* empty the contents of the string but do NOT
remove the pointer.  This is useful for clearing
one string within a linked list of them. *)
var   index: integer; (* to the ribbon *)
begin (* clearstring *)
   with ribbon do begin
      for index := 1 to maxstring do letters[index] := ' ';
      length := 0;
      current := 0;
   end
end; (* emptystring *)
(* ---- *)
procedure clearstring(var ribbon: string);
(* empty the string and remove the pointer.
As of 2010 Dec 04 there is a previous pointer that
allows proper removal of the linked list using disposestring.
This clearstring routine is for the initial clearing
of the variables. *)
begin (* clearstring *)
   with ribbon do begin
      emptystring(ribbon);
      previous := nil;
      next := nil;
   end
end; (* clearstring *)
(* ---- *)
procedure disposestring(var ribbon: stringptr);
(* remove the linked list from the string in revere order *)
var
   t: stringptr; (* a link in the string *)
begin
   t := ribbon;
   (* find the end *)
   while t^.next <> nil do t := t^.next;
   (* move backwards and dispose links *)
   while t^.previous <> nil do begin
      t := t^.previous;
      dispose(t^.next);
   end;
   dispose(t);
end; (* disposestring *)
(* ---- *)
procedure initializestring(var ribbon: string);
(* start the string with a nil pointer.   This routine
should be called before doing linked list work.  This allows
the standard string routines to clear the string without
killing the pointer.  This is now deprecated, do not use it
since clearstring still clears the next pointer. *)
begin (* initializestring *)
   writeln(output,'remove initializestring routine!');
   writeln(output,'replace it with clearstring routine!');
   halt; (* to force deprecation *)
   clearstring(ribbon);
   ribbon.next := nil;
end; (* initializestring *)
(* end module clearstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module detabstring *)
procedure detabstring(var s: string);
(* Tabs are converted to blanks.  This allows the tab to be a special "solid"
kind of blank, but the string written won't vary with position on the line
depending on the tabstops. *)
const
   tabcharacter = 9; (* the ordinal of the tab character *)
var
   c: char; (* a character write out *)
   i: integer; (* index to the string s *)
begin
   with s do for i := 1 to length do begin
      c := letters[i];
{
for testing:
      if ord(c) = tabcharacter then c := 'T'; (* tab converted to blank *)
}
      if ord(c) = tabcharacter then c := ' '; (* tab converted to blank *)
      letters[i] := c
   end;
end;
(* end module detabstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module interact.writestring *)
(* begin module writestring *)
procedure writestring(var tofile: text; var s: string);
(* write the string s to file tofile, no writeln *)
var   i: integer; (* index to s *)
begin (* writestring *)
      with s do for i := 1 to length do write(tofile, letters[i])
end; (* writestring *)
(* end module writestring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)
(* end module interact.writestring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module readquotestring *)
procedure readquotestring(var afile: text; var s: string);
(* read the quoted string from afile, remove the quotes.  Trailing
blanks are converted to tab characters. *)
const
   tabcharacter = 9; (* the ordinal of the tab character *)
var
   c: char; (* a character read *)
   i: integer; (* index to the string s *)
begin
   emptystring(s);
   skipblanks(afile);
   with s do begin
      read(afile,c);
      if c <> '"' then begin
         writeln(output, 'quote string expected but " missing');
         writeln(output, 'character found instead:"',c,'"');
         halt
      end;
      if eoln(afile) then begin
         writeln(output, 'missing end of quote string " string');
         halt
      end;
      repeat
         read(afile,c);
         if c <> '"' then begin
            length := succ(length);
            letters[length] := c;
         end;
      until (c = '"')
            or eoln(afile) or (length = maxstring);
      if length = maxstring then begin
         writeln(output, 'A quote string exceeded ',maxstring:1,
                         ' characters.');
         writeln(output, 'Use a shorter string or increase constant',
                         ' maxstring.');
         halt
      end;
      if c <> '"' then begin
         writeln(output, 'The second " mark of a string is missing.');
         halt
      end;
      (* convert all blanks to tab characters *)
      for i := 1 to length do begin
         if letters[i] = ' '
         then letters[i] := chr(tabcharacter); (* set to tab character *)
      end;
   end;
end;
(* end module readquotestring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module readparanstring *)
procedure readparanstring(var afile: text; var s: string);
(* read the quoted string from afile, remove the quotes.  Trailing
blanks are converted to tab characters. *)
const
   tabcharacter = 9; (* the ordinal of the tab character *)
var
   c: char; (* a character read *)
   i: integer; (* index to the string s *)
begin
   clearstring(s);
   skipblanks(afile);
   with s do begin
      read(afile,c);
      if c <> '(' then begin
         writeln(output, 'paran string expected but ( missing');
         writeln(output, 'character found instead:"',c,'"');
         halt
      end;
      if eoln(afile) then begin
         writeln(output, 'missing end of paran string ) string');
         halt
      end;
      repeat
         read(afile,c);
         if c <> ')' then begin
            length := succ(length);
            letters[length] := c;
         end;
      until (c = ')')
            or eoln(afile) or (length = maxstring);
      if length = maxstring then begin
         writeln(output, 'A paran string exceeded ',maxstring:1,
                         ' characters.');
         writeln(output, 'Use a shorter string or increase constant',
                         ' maxstring.');
         halt
      end;
      if c <> ')' then begin
         writeln(output, 'The ) mark of a string is missing.');
         halt
      end;
      (* convert all blanks to tab characters *)
      for i := 1 to length do begin
         if letters[i] = ' '
         then letters[i] := chr(tabcharacter); (* set to tab character *)
      end;
   end;
end;
(* end module readparanstring version = 4.31; (@ of prgmod.p 1999 July 13 *)

(* begin module writequotestringB *)
procedure writequotestringB(var afile: text; s: string);
(* Write the string s to afile, surrounded by double quote marks.  Tabs are
converted to blanks.  This allows the tab to be a special "solid" kind of
blank, but the string written won't vary with position on the line depending
on the tabstops. *)
begin
   detabstring(s);
   writestring(afile, s);
end;
(* end module writequotestringB *)

(* begin module writequotestring *)
procedure writequotestring(var afile: text; s: string);
(* Write the string s to afile, surrounded by double quote marks.  Tabs are
converted to blanks.  This allows the tab to be a special "solid" kind of
blank, but the string written won't vary with position on the line depending
on the tabstops. *)
begin
   detabstring(s);
   write(afile, '"');
   writestring(afile, s);
   write(afile, '"');
end;
(* end module writequotestring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module grabtoken *)
procedure grabtoken(var thefile: text; var thestring: string);
(* skip any blanks and then grab the next token from the file *)
var
   c: char; (* a character in thefile *)
   done: boolean; (* done finding the name *)
begin
   skipblanks(thefile);
   done := false;
   with thestring do begin
      length := 0;
      while not done do begin
         if eoln(thefile)
         then done := true
         else begin
            read(thefile,c);
            if c = ' '
            then done := true
            else begin
               length := succ(length);
               letters[length] := c;
            end
         end
      end
   end
end;
(* end module grabtoken version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module equalstring *)
function equalstring(a, b: string): boolean;
(* Test for equality between two strings at current positions.
NOTE:  A compiler bug results if one directly tests this way:
   if thedefinition^.nametag = aname
The reason is completely not clear!  I showed that the parts of the
strings were identical, but the whole was not by this test.
For this reason it is *critical to test strings with equalstring.
*)
var   index: integer; (* index to both strings *)
      equal: boolean; (* are letters in a and b the same? *)
begin (* equalstring *)
      if a.length = b.length
      then begin
         index := 1;
         repeat
            equal := (a.letters[index] = b.letters[index]);
            index := succ(index)
         until (not equal) or (index > a.length);
         equalstring := equal
      end
      else equalstring := false
end; (* equalstring *)
(* end module equalstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module readstringline *)
procedure readstringline(var afile: text; var buffer: string);
(* get a line (as a string) from a file. *)
var   index: integer; (* of buffer *)
begin (* getstring *)
   clearstring(buffer);
   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, 'readstringline: a line exceeds maximum string size (',
                      maxstring:1,')');
      halt
   end;
   buffer.length := index;
   buffer.current := 1;
   readln(afile);
end; (* readstringline *)
(* end module readstringline version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module readheader *)
procedure readheader(var afile: text; ch: char; var theheader: stringptr);
(* Read the header from a file.  Only read lines that begin with ch.
2004 Dec 15: Handle end of file (eof) gracefully. *)
var
   h: stringptr; (* pointer to the riblheader *)
begin
   if eof(afile) then begin
      clearstring(theheader^);
   end
   else begin
      if afile^ = ch then begin
         new(theheader);
         h := theheader;
         while afile^='*' do begin
            readstringline(afile, h^);
            if not eof(afile) then begin
               if afile^='*' then begin
                  new(h^.next);
                  h := h^.next;
               end
               else h^.next := nil;
            end
         end
      end
      else begin
         theheader := nil
      end;
   end;
end;
(* end module readheader version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module writestringlines *)
procedure writestringlines(var afile: text; l: stringptr);
(* write lines of strings to a file. *)
var
   h: stringptr; (* pointer to l *)
begin (* getstring *)
   if l <> nil then begin
      h := l;
      while h <> nil do begin
         writestring(afile, h^);
         writeln(afile);
         h := h^.next;
      end;
   end
end;
(* end module writestringlines version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module copynamestring *)
procedure copynamestring(n: name; var s: string);
(* copy the delila type name n to the string s *)
var
   i: integer; (* index to s *)
begin
   clearstring(s);
   s.length := n.length;
   for i := 1 to s.length do begin
      s.letters[i] := n.letters[i];
   end
end;
(* end module copynamestring *)

(* begin module copylinestring *)
procedure copylinestring(n: lineptr; var s: string);
(* copy the delila type line n to the string s *)
var
   i: integer; (* index to s *)
begin
   clearstring(s);
   s.length := n^.length;
   for i := 1 to s.length do begin
      s.letters[i] := n^.letters[i];
   end
end;
(* end module copylinestring *)

(* begin module copystring *)
procedure copystring(a: string; var b: string);
(* copy string a to b.  The length, characters and
where the string points to are all copied. *)
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];
   b.current := a.current;
   b.next := a.next;
end;
(* end module copystring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module concatstring *)
procedure concatstring(var a, b: string; var c: string);
(* concatenate a and b and put in c *)
var   abindex: integer; (* index to both strings *)
      cindex: integer; (* index to both strings *)
begin (* concatstring *)
   emptystring(c);
   c.length := a.length + b.length;
   if c.length > maxstring then begin
      writeln(output,'string too long in concatstring');
      halt
   end;

   cindex := 0;
   for abindex := 1 to a.length do begin
      cindex := succ(cindex);
      c.letters[abindex] := a.letters[abindex]
   end;

   for abindex := a.length+1 to a.length+b.length do begin
      cindex := succ(cindex);
      c.letters[abindex] := b.letters[abindex-a.length]
   end;
end; (* concatstring *)
(* end module concatstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module attachstring *)
procedure attachstring(var a: string; b: string);
(* attach the string b onto the end of a *)
var
   hold: string; (* intermediate for doing the attachment *)
begin
   copystring(a, hold);
   concatstring(hold,b,a)
end;
(* end module attachstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module joinstring *)
procedure joinstring(parts: stringptr; var target: string);
(* the parts string contains one or more strings linked
together.  They are to be concatenated together and put on the end
of the target string *)
var
   p: stringptr; (* pointer into the parts *)
begin
   p := parts;
   while p <> nil do begin
      attachstring(target, p^);
      p := p^.next;
   end;
end;
(* end module joinstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module attachcharacterstring *)
procedure attachcharacterstring(c: char; var a: string);
(* attach the character c to the end of a *)
begin
   if a.length = maxstring then begin
      writeln(output,'attachcharacterstring: maxstring exceeded');
      halt
   end;
   with a do begin
      length := succ(length);
      letters[length] := c
   end
end;
(* end module attachcharacterstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

(* 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 = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* 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 = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module numberrealstring *)
procedure numberrealstring(var s: string;
                            n: real;
                            thewidth: integer;
                            thedecimal: integer);
(* Put the real number n into the string s taking up thewidth characters and
using thedecimal decimal places.  Note that the procedure multiplies n by a
perfectly integer power (which is quite clean) and then rounds that.  This
ensures that the last decimal place is rounded but no others.  Then the
program picks apart that clean number using numberdigit.  NOTE: If thewidth
is too small to fit the number then the size is increased.
Note:  The general rule is to throw away leading zeros EXCEPT in the
special case like 0.1!  This code does not account for that case.
2005 May 30: That's not acceptable.  No leading spaces are allowed
unless the string is shorter than thewidth.  Add a bit of code to do that.
*)
{
const dndebugpoint = -20; updebugpoint = -0;
const dndebugpoint = -1.1; updebugpoint = 100;
}
var
   decimalplace: integer; (* amount of space for a decimal place *)
   d: integer; (* one of the digit places *)
   digits: integer; (* the number n is multiplied by powers of 10 until
      all relevant digits are above zero.  Then it is rounded to produce
      'digits'.  This guarantees that the rounding is correct over the
      entire number. *)
   i: integer; (* index to the digits of the number *)
   power: integer; (* powers of 10 *)
   minus: boolean; (* whether we made a - *)
   newwidth: integer; (* the larger width to fit the number into *)
   sign: integer; (* amount of space for a sign *)
   symbol: char; (* a symbol to put into s *)
begin
   (* numberrealstring MUST NOT CALL clearstring(s); because
      clearstring WIPES OUT the rest of a linked list *)
   emptystring(s);
   if thewidth > 0 then with s do begin
      power := 1;
      for i := 1 to thedecimal do begin
         power := power * 10;
      end;
      digits := round(power*abs(n));
      if thedecimal > 0 then decimalplace := 1
                        else decimalplace := 0;
      if n < 0 then sign := 1
               else sign := 0;
      newwidth := numbersize(digits)+thedecimal+sign;

      (* add space for leading zero for numbers below abs(1).
      2005 May 30: This puts an extra space
      at the start of 0.97 when using only one decimal place.
      This must be removed later in some cases *)
      if abs(n) < 1.0 then newwidth := succ(newwidth);

      (* use the requested width if it is bigger *)
      if thewidth > newwidth
      then newwidth := thewidth;

{
if (n > dndebugpoint) and (n < updebugpoint) then begin
writeln(output);
write(output,'n = ',n:1:1);
write(output,'  n = ',n:thewidth:(thedecimal+3));
write(output,'  digits = ',digits:1);
write(output,'  trunc(n) = ',trunc(abs(n)):1);
write(output,'  numbersize(trunc(abs(n))) = ',numbersize(trunc(abs(n))):1);
writeln(output);
write(output,'  newwidth = ',newwidth:1);
write(output,'  thewidth = ',thewidth:1);
write(output,'  thedecimal = ',thedecimal:1);
writeln(output);
end;
}

      for d := newwidth-1-decimalplace downto 0 do begin
          symbol := numberdigit(digits,d);
{
if (n > dndebugpoint) and (n < updebugpoint) then begin
write(output,'        d = ',d:1:1);
write(output,'        symbol = ',symbol:1);
writeln(output);
end;
}
          if symbol = '-' then symbol := '0';
          if symbol = '+' then symbol := '0';
          if symbol = ' ' then symbol := '0';
          length := succ(length);
          letters[length] := symbol;
          if d = thedecimal then if thedecimal > 0 then begin
{
if (n > dndebugpoint) and (n < updebugpoint) then begin
write(output,'            at the decimal    = .');
writeln(output);
end;
}
             length := succ(length);
             letters[length] := '.';
          end;
      end;
{
if (n > dndebugpoint) and (n < updebugpoint) then begin
writequotestring(output,s); writeln(output);
end;
}

      (* convert leading zero's to blanks *)
      current := 1;
      while (letters[current] = '0') and (current < length) do begin
         if letters[current+1] <> '.' then letters[current] := ' ';
         current := succ(current);
{
if (n > dndebugpoint) and (n < updebugpoint) then begin
writequotestring(output,s); writeln(output);
end;
}

      end;
      (* now put the minus sign back *)
      if n < 0 then begin
         current := 1;
         minus := false;
         while letters[current] = ' ' do begin
            if letters[current+1] <> ' ' then begin
               letters[current] := '-';
               minus := true;
            end;
            current := succ(current);
            if current >= maxstring then begin
               write  (output,'numberrealstring:');
               writeln(output,' cannot place minus sign in string');
               writeln(output,'placing number at ',current:1);
               writeln(output,'the number is n =',n:newwidth:thedecimal);
               writestring(output,s);
               halt;
            end;
         end;
         if not minus then begin
            writeln(output,'numberrealstring ERROR:',
                           ' missing minus sign for number ',
                           n:newwidth:thedecimal);
            writeln(output,'Increase one of the number widths.');
            halt
         end;
      end;
   end;

   (* 2005 May 30: At this point the number 0.97 with one decimal
   place comes out as: " 1.0".  This is not acceptable, so remove
   any leading blanks. *)
   if s.length > thewidth
   then with s do while letters[1] = ' ' do begin
      for i := 1 to length - 1
      do letters[i] := letters[i + 1];
      length := length - 1;
   end;

{
if (n > dndebugpoint) and (n < updebugpoint) then begin
writeln(output,'s.length = ',s.length:1);
writeln(output,'n=',n:thewidth:thedecimal);
writequotestring(output,s); writeln(output);
end;
}

end;
(* end module numberrealstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module readname *)
procedure readname(var afile: text; var id: name);
(* read a name from afile, code is similar to brname *)
var
   c: char; (* a character read *)
   i: integer; (* an index *)
begin
   skipblanks(afile);
   with id do begin
      length:=0;
      repeat
         length:=succ(length);
         read(afile,c);
         letters[length] := c
      until (eoln(afile)) or
            (length>=namelength) or
            (c=' ');
      if letters[length]=' ' then length:=length-1;
      if length<namelength
         then for i:=length+1 to namelength do letters[i]:=' '
   end;

end;
(* end module readname version = 9.25; (@ of lister.p 1999 June 9 *)

(* begin module writename *)
procedure writename(var afile: text; s: name);
(* write the string s to afile *)
var i: integer; (* index to the string s *)
begin
   with s do for i := 1 to length do write(afile,letters[i]);
end;
(* end module writename version = 9.25; (@ of lister.p 1999 June 9 *)

(* begin module readwaveparameters *)
procedure readawaveparameter(var afile: text; var wp: waveparam);
{
                                 cmperbase, cmperbit: real);
}
(* Read one wave parameter from a file. *)
{
2000 Nov 1:  The following experimental code is not implemented:
See the coscurve program
(https://alum.mit.edu/www/toms/delila/coscurve.html) for details
about the conversion from bases to cm.  The code now requires
cm per base and cmperbit.
const
   lover2pi = 1.2160067233; (* L/2pi for L being the length of
       the cosine function from 0 to 2pi *)
var
   multiplier: real; (* factor to multiply by to convert from bases to cm. *)
}
begin
   readln(afile,wp.extreme);
   readln(afile,wp.wavelocation);
   readln(afile,wp.wavebit);
   readln(afile,wp.waveamplitude);
   readln(afile,wp.wavelength);
   if wp.wavelength <= 0.0 then begin
      writeln(output,'wave parameters: wavelength must be positive');
      halt
   end;
   (* one may skip the dash and thickness parameters
      of the last wave definition.  This allows
      backwards compatability, but should not be
      used in general. *)
   if eof(afile)
   then begin
      wp.dashon := 0;
      wp.dashoff := 0;
      wp.dashoffset := 0;
      wp.thickness := 0
   end
   else begin
      (* implement new method to read dash information
         without breaking previous wave parameter files *)
      if afile^ <> 'd'
      then begin (* default to old method *)
         readln(afile,wp.dashon);
         wp.dashoff := wp.dashon;
         wp.dashoffset := 0;
      end
      else begin
         get(afile); (* skip the 'd' *)
         (* reading the parameters directly in cm would require
            that the user compute or guess the values.
            So read in as bases and compute cm.
            Unfortunately the length of a curve is not the simple
            scale factor.  Make the user guess for now.
         *)
         readln(afile,wp.dashon, wp.dashoff, wp.dashoffset)

{
experimental code to account for actual wave length:
         readln(afile, wp.BASEdashon, wp.BASEdashoff, wp.BASEdashoffset);
         multiplier    := 1;
         multiplier    := lover2pi
 * wp.waveamplitude (* in bits *)
 * cmperbit (* converts to cm *)
 / wp.wavelength (* 10.6 bp/turn *)
* cmperbase
;
multiplier    := lover2pi;
multiplier    := 1.0;
multiplier    := 1.0/lover2pi;
         multiplier := lover2pi
                       * wp.waveamplitude (* in bits *)
                       * cmperbit (* converts to cm *)
                       * wp.wavelength (* 10.6 bp/turn *)
                       * cmperbase;
         multiplier := 1.0
         wp.dashon     := multiplier * wp.BASEdashon;
         wp.dashoff    := multiplier * wp.BASEdashoff;
         wp.dashoffset := multiplier * wp.BASEdashoffset;
with wp do begin
   writeln(output, 'cmperbase  = ',cmperbase:5:2);
   writeln(output, 'cmperbit   = ',cmperbit:5:2);
   writeln(output, 'dashon     = ',dashon:5:2);
   writeln(output, 'dashoff    = ',dashoff:5:2);
   writeln(output, 'dashoffset = ',dashoffset:5:2);
   writeln(output, 'waveamplitude   = ',wp.waveamplitude:5:2);
   writeln(output, 'wavelength      = ',wp.wavelength:5:2);
end;
}

      end;
      if eof(afile)
      then wp.thickness := 0
      else readln(afile,wp.thickness);
   end;
end;

procedure readwaveparameters(var afile: text; var w: waveptr);
{
                                 cmperbase, cmperbit: real);
}
(* read from afile the wave parameters w.  The routine is done either at end
of file, when a period ('.') appears or at a blank line.  This makes it
compatable with different reading mechanisms. *)
var
   done: boolean; (* done reading *)
   morethanone: boolean; (* more than one wave? *)
   p: waveptr; (* one of the definitions of a cosine wave *)
procedure waystoend;
(* if end of file, a blank line or a period, consider the job done *)
var
   clear: boolean; (* not eof, so clear to test for comment *)
begin
   if eof(afile) then done := true;
   clear := false;
   if not done (* skip any number of comments *)
   then while (not clear) and (not done) do begin
      if eof(afile) then begin
         done := true;
         clear := true;
      end
      else if afile^ = '*'
           then readln(afile)
           else clear := true
   end;
   if not done then if eoln(afile) then done := true;
   if not done then if afile^='.' then begin
      readln(afile);
      done := true;
   end;
end;
begin
   done := false;
   morethanone := false;
   w := nil;
   waystoend;
   if not done then begin
      new(w);
      p := w;
      while not done do begin
         waystoend;
         if done
         then p^.next := nil
         else begin
            if afile^='*'
            then readln(afile) (* skip any number of comments *)
            else begin
               if morethanone then begin
                  new(p^.next);
                  p := p^.next
               end;
{
               readawaveparameter(afile, p^, cmperbase, cmperbit);
}
               readawaveparameter(afile, p^);
               (* if we ever come back again there will be more than one: *)
               morethanone := true;
            end
         end
      end
   end
   else w := nil;
end;
(* end module readwaveparameters version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module writewaveparameters *)
procedure writewaveparameters(var afile: text; w: waveptr);
(* write wave parameters to a file *)
const
   wid = 4; (* width of numbers *)
var
   i: waveptr; (* pointer to the waves in w *)
begin
   i := w;
   while i <> nil do with i^ do begin
      writeln(afile,extreme,
                    '  extreme: char; h or l, the high or low extreme to be defined');
      writeln(afile,wavelocation:wid:1,
                    '  wavelocation: real; the location in bases of the extreme');
      writeln(afile,wavebit:wid:1,
                    '  wavebit: real; the location in bits of the extreme');
      writeln(afile,waveamplitude:wid:1,
                    '  waveamplitude: real; the amplitude of the wave in bits');
      writeln(afile,wavelength:wid:1,
                    '  wavelength: real; the wave length of the wave in bases');
{
      writeln(afile,dash:wid:1,
                    '  dash: real; the size of dashes in cm.',
                    '  dash <= 0 means no dashes');
      writeln(afile,    BASEdashon:wid:1,
                   ' ', BASEdashoff:wid:1,
                   ' ', BASEdashoffset:wid:1,
                   ' dashon, dashoff, dashoffset (in bases)');
}
      writeln(afile,    dashon:wid:1,
                   ' ', dashoff:wid:1,
                   ' ', dashoffset:wid:1,
                   ' dashon, dashoff, dashoffset (in bases)');
      writeln(afile,thickness:wid:1,
                    '  thickness: real; the thickness of the cosine wave.',
                    '  <=0 means default');
      i := i^.next
   end;
   writeln(afile,'.');
end;
(* end module writewaveparameters version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module getriblmatrix *)
procedure getriblmatrix(var afile: text; var matrix: ribltypeptr);
(* Get the matrix from a file, with the defining coordinate limits,
followed by the mean and standard deviation and other data.
The entire matrix is read in. *)
const
   clearvalue = 5.0; (* If a large bit position is ever seen (for DNA), we
      KNOW there is a bug in the code! *)
var
   b: base; (* a base in the matrix *)
   l: integer; (* a coordinate in the matrix *)
   thel: integer; (* the l value in the file (should match l) *)

   teststring: string; (* string to test for 'ribl ' at matrix start *) 
   riblstart: string; (* word at the start of the matrix *)

procedure testforend;
(* test for end of file *)
begin
   if eof(afile) then begin
      writeln(output,'part or all of the ribl matrix is missing');
      halt;
   end
end;
procedure skip;
(* skip comment lines *)
begin
   testforend;
   while afile^='*' do begin
     readln(afile);
     testforend
   end
end;
begin
   new(matrix);
   with matrix^ do begin

      (* read name of the matrix *)
      grabtoken(afile, riblstart);
      teststring.letters[1] := 'r'; 
      teststring.letters[2] := 'i';
      teststring.letters[3] := 'b';
      teststring.letters[4] := 'l';
      teststring.length := 4; 
 
      if not equalstring(riblstart, teststring) then begin
         writeln(output,'The ribl matrix must begin with ''ribl "name"''');   
         writequotestring(output,riblstart);
         writeln(output, ' seen instead');
         halt
      end;
 
      skipblanks(afile); (* skip any spaces *)
      readquotestring(afile, matrix^.riblname);
      readln(afile);

      readheader(afile,'*',riblheader);

      readln(afile,frombase,tobase);

      if frombase < minribl then begin
         writeln(output,'Matrix is too large for program.');
         writeln(output,'frombase (=',frombase:1,')');
         writeln(output,'minribl (=',minribl:1,')');
         writeln(output,'Decrease minribl or increase frombase');
         halt
      end;

      if tobase > maxribl then begin
         writeln(output,'Matrix is too large for program.');
         writeln(output,'tobase (=',tobase:1,')');
         writeln(output,'maxribl (=',maxribl:1,')');
         writeln(output,'Increase maxribl or decrease tobase');
         halt
      end;

      (* clear unused matrix for safety *)
      for l := minribl to frombase-1
      do for b := a to t do matrix^.data[b,l] := clearvalue;
      for l := tobase+1 to maxribl
      do for b := a to t do matrix^.data[b,l] := clearvalue;

      for l := frombase to tobase do begin
         testforend;
         for b := a to t do read(afile,matrix^.data[b,l]);
         read(afile,thel);
         if thel <> l then begin
            writeln(output, 'positions in ribl do not match from-to range');
            halt;
         end;
         for b := a to t do read(afile,matrix^.numbers[b,l]);
         readln(afile)
      end;

      skip; readln(afile,mean);
      skip; readln(afile,stdev);
      skip; readln(afile,consensus);
      skip; readln(afile,anticonsensus);
      skip; readln(afile,averageRi);
      skip; readln(afile,n);
      skip; readln(afile,symmetry);
      skip; readln(afile,Ribound);
      skip; readln(afile,Zbound);
      skip; readln(afile,Pbound);
{this is a guess - set the converters to neutral for now
      cmperbase := 1.0;
      cmperbit := 1.0;
      readwaveparameters(afile, waves, cmperbase, cmperbit);
}
      readwaveparameters(afile, waves);
   end;
end;
(* end module getriblmatrix version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module putriblmatrix *)
procedure putriblmatrix(var afile: text;
                    matrix: ribltypeptr);
(* write the ribl matrix to afile. *)
var
   lindex: integer; (* index to Riblmatrix, equivalent to l *)
begin
   if matrix <> nil then begin
      with matrix^ do begin
         write(afile,'ribl');
         write(afile,' ');
         writequotestring(afile, riblname);
         writeln(afile);
         (* write the header *)
         writestringlines(afile, riblheader);
         writeln(afile, frombase: nfield,' ', tobase: nfield,' frombase, tobase');
         for lindex := frombase to tobase do begin
            write(afile,' ',
                       ' ',data[a, lindex]:infofield:infodecim,
                       ' ',data[c, lindex]:infofield:infodecim,
                       ' ',data[g, lindex]:infofield:infodecim,
                       ' ',data[t, lindex]:infofield:infodecim,
                       ' ',lindex:nfield);
            write(afile,' ',numbers[a,lindex]:nfield,
                        ' ',numbers[c,lindex]:nfield,
                        ' ',numbers[g,lindex]:nfield,
                        ' ',numbers[t,lindex]:nfield);
            writeln(afile);
         end;
         writeln(afile,'*');
         writeln(afile,mean:infofield:infodecim, 
                       ' bits = mean (Rsequence of selected region)');
         writeln(afile,stdev:infofield:infodecim,
                       ' bits = standard deviation');
         writeln(afile,'*');
         writeln(afile,consensus:infofield:infodecim,
                       ' bits = Ri of consensus sequence',
                       ' from ',frombase:1,' to ',tobase:1);
         writeln(afile,anticonsensus:infofield:infodecim,
                       ' bits = Ri of anticonsensus sequence',
                       ' from ',frombase:1,' to ',tobase:1);
         writeln(afile,'*');
         writeln(afile,averageRi:infofield:infodecim,
                       ' bits = average Ri for random sequence',
                       ' from ',frombase:1,' to ',tobase:1);
         writeln(afile,'*');
         writeln(afile,n:1,
                       ' n, number of sequences used to create the matrix');
         writeln(afile,'*');
         case symmetry of
           'a': write(afile, 'asymmetric');
           'e': write(afile, 'even');
           'o': write(afile, 'odd');
         end;
         writeln(afile, ' symmetry of the matrix');
         writeln(afile,'*');
         writeln(afile,Ribound:infofield:infodecim,
                       ' Ri bound: lower bound on Ri');
         writeln(afile,Zbound:infofield:infodecim,
                       ' Z bound: lower bound on Z');
         writeln(afile,Pbound:infofield:infodecim,
                       ' P bound: upper probability');
         writeln(afile,'*');
         writewaveparameters(afile, waves);
         writeln(afile);
      end
   end
end;
(* end module putriblmatrix version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module writeadefinition *)
procedure writeadefinition(var afile: text; d: definitiontypeptr);
(* write the definition d to file afile *)
var i: integer; (* index to the string s *)
begin
   with d^ do begin
      write(afile,'define');
      write(afile,' '); writequotestring(afile,nametag);
      write(afile,' '); writequotestring(afile,background);
      write(afile,' '); writequotestring(afile,negparts);
      write(afile,' '); writequotestring(afile,posparts);
      for i := 1 to marks do write(afile,' ',locations[i]:1:1);
      writeln(afile);
{
      (* don't do this here - it's awkward if one wants
         to just write the parts above.  Just rember to write the
      matrix after the definition! *)
      writeln(afile);
      putriblmatrix(afile, matrix);
}
   end;
end;
(* end module writeadefinition version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module showadefinition *)
procedure showadefinition(var afile: text; d: definitiontypeptr);
(* show the definition d to file afile *)
const
   pwid = 8; (* width in character places to print PostScript numbers *)
   pdec = 5; (* decimal places to print PostScript numbers *)
var i: integer; (* index to the string s *)
begin
   writeln(afile,'DEFINITION: ');
   with d^ do begin
      write(afile,'nametag: ');
      writequotestring(afile,nametag);
      writeln(afile);

      write(afile,'background: ');
      writequotestring(afile,background);
      writeln(afile);

      write(afile,'negative parts: ');
      writequotestring(afile,negparts);
      writeln(afile);

      write(afile,'positive parts: ');
      writequotestring(afile,posparts);
      writeln(afile);

      write(afile,'locations: ');
      for i := 1 to marks do write(afile,' ',locations[i]:1:1);
      writeln(afile);

      writeln(afile,'marks: ',marks:1);
      writeln(afile,'min: ',min:5:1);
      writeln(afile,'max: ',max:5:1);

      write(afile,'petal: ');
      with petal do begin
         write(afile,el:1);
         write(afile,' ', c:1);
         write(afile,' ',eh:pwid:pdec);
         write(afile,' ',es:pwid:pdec);
         write(afile,' ',eb:pwid:pdec);
         write(afile,' ',fh:pwid:pdec);
         write(afile,' ',fs:pwid:pdec);
         write(afile,' ',fb:pwid:pdec);
         write(afile,' ',spare:pwid:pdec);
      end;
      writeln(afile);

{
      write(output,'THIS IS');
      if d^.istotal
      then writeln(output,      ' a total.')
      else if d^.isgap
           then writeln(output, ' a Gap.')
           else writeln(output, ' NEITHER a Gap nor a total.');
}
{sss}

   end;
end;
(* end module showadefinition version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module readadefinition *)
procedure readadefinition(var afile: text; var d: definitiontypeptr);
(* read a single definition from afile *)
var
   i: integer; (* an index *)
   location: real; (* a location on the sequence *)

procedure IsGapOrIsTotal(var d: definitiontypeptr);
(* Is this a definition which is for a Gap or a total connecting
multi-walkers?  If so, set isgap or istotal true (respectively).  *)
const
   tab = 9; (* tab characters separate name parts *)
   debugging = false; (* use to debug the code *)
var
   totstring: string; (* the string that indicates a IsGapOrIsTotal total *)
   gapstring: string; (* the string that indicates a IsGapOrIsTotal gap *)
   deflocation: integer; (* location in the nametag *)
   gaplocation: integer; (* location in the gapstring *)
   totlocation: integer; (* location in the total *)
   done: boolean; (* done comparing strings *)
begin

if debugging then begin
  writeln(output,'IsGapOrIsTotal test of this definition:'); {sss}
  writeadefinition(output,d); {sss}
end;

   clearstring(totstring);
   with totstring do begin
      letters[1] := chr(tab);
      letters[2] := 't';
      letters[3] := 'o';
      letters[4] := 't';
      letters[5] := 'a';
      letters[6] := 'l';
      length := 6;
   end;

   clearstring(gapstring);
   with gapstring do begin
      letters[1] := chr(tab);
      letters[2] := 'G';
      letters[3] := 'a';
      letters[4] := 'p';
      length := 4;
   end;

{
writeln(output, 'IsGapOrIsTotal: begin &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
}

   (* search for the Gap string *)
   done := false;
   d^.isgap := true;
   deflocation := d^.nametag.length;
   gaplocation := gapstring.length;
   while not done do begin
{
write(output,
'====>',
d^.nametag.letters[deflocation],
gapstring.letters[gaplocation]);
}
      if d^.nametag.letters[deflocation] <> gapstring.letters[gaplocation]
      then begin
         d^.isgap := false;
         done := true
      end;
      deflocation := deflocation - 1;
      gaplocation := gaplocation - 1;
{
write  (output,' deflocation = ',deflocation:1);
writeln(output,' gaplocation = ',gaplocation:1);
}
      if deflocation < 1 then done := true;
      if gaplocation < 1 then done := true;
   end;

(* just for debugging:
   if d^.isgap then begin
      if debugging then write(output, 'gap:   ');
      write(output, '====================');
      writestring(output, d^.nametag);
      writeln(output, '=======');
      writeadefinition(output, d);
      writeln(output);
      writeln(output, '====================');
{ok pointing to the features of this
has to happen with each feature of this type }
{      d^.isgap := true why set this???}
   end
   else begin
{
      if debugging then write(output, '       ');
      writestring(output, d^.nametag);
      writeln(output);
}
{      d^.isgap := false why set this???}
   end;
*)

   (* search for the total string *)
   d^.istotal := true;
   done := false;
   (* skip to the tab before the bits *)
   deflocation := d^.nametag.length;
   while not done do begin
{
write(output,d^.nametag.letters[deflocation]);
}
      if d^.nametag.letters[deflocation] = chr(tab)
      then done := true;
      deflocation := deflocation - 1;
      if deflocation < 1 then begin
         done := true;
         d^.istotal := false; (* fell off begin of string *)
      end;
   end;
{
writeln(output,'<-PRETOTAL');
}
   if d^.istotal = true then begin (* it still might be true *)
{
writeln(output,'...');
}
      totlocation := totstring.length;
      done := false;
{
this is defined above!
      deflocation := d^.nametag.length;
}
      while not done do begin
{
write(output,
'====>',
d^.nametag.letters[deflocation],
totstring.letters[totlocation]);
}
         if  d^.nametag.letters[deflocation]
         <> totstring.letters[totlocation]
         then begin
            d^.istotal := false;
            done := true
         end;
         deflocation := deflocation - 1;
         totlocation := totlocation - 1;
{
write  (output,' deflocation = ',deflocation:1);
writeln(output,' totlocation = ',totlocation:1);
}
         if deflocation < 1 then done := true;
         if totlocation < 1 then done := true;
      end;
   end;

   (* mark that this definition is part of a multiple set *)
   if d^.isgap or d^.istotal
   then d^.ismulti := true
   else d^.ismulti := false;

   if debugging then begin
{sss}
{
      writeadefinition(output, d);
}
      showadefinition(output, d);
{
if d^.istotal then begin
writeln(output, 'IsGapOrIsTotal: is total');
halt;
end;
}
      if d^.isgap
      then writeln(output,'Is a gap')
      else writeln(output,'Is NOT a gap');
      if d^.istotal
      then writeln(output,'Is a total')
      else writeln(output,'Is NOT a total');
      writeln(output, 'IsGapOrIsTotal: end &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');
   end;

end; (* IsGapOrIsTotal *)

begin (* readadefinition *)
{writeln(output,'readadefinition');}
   skipnonblanks(afile);
   skipblanks(afile);
   with d^ do begin
      readquotestring(afile, nametag);

      (* read background character from a quotestring *)
      readquotestring(afile,background);
      if (background.length < 1) or
         (background.length > 2) then begin
         writeln(output,'Background string must contain only 1 or 2',
                        ' characters.');
         write(output,'Instead, it contained: ');
         writequotestring(output,background);
         writeln(output);
         halt
      end;

      readquotestring(afile,negparts);
      readquotestring(afile,posparts);
      if negparts.length <> posparts.length then begin
         write(output,'While reading definition of ');
         writequotestring(output,nametag);
         writeln(output);
         write(output,'the negative part ');
         writequotestring(output,negparts);
         writeln(output,' and');
         write(output,'the positive part ');
         writequotestring(output,posparts);
         writeln(output,' are not the same length.');
         writeln(output,'(',negparts.length:1,'<>',posparts.length:1,')');
         halt
      end;

      marks := negparts.length;
      if marks = 0 then begin
         (* multiscan uses empty strings sometimes ... *)
         marks := 2
      end;

      for i := 1 to marks do read(afile,locations[i]);

      matrix := nil;
      RiboundFeature := -maxint;
      ZboundFeature := +maxint;
      PboundFeature := 0.0;

      ismulti := false;

      next := nil;

      (* checks *)

      (* now compute the absolute coordinates that the feature resides in *)
      min := +maxint;
      max := -maxint;
      for i := 1 to marks do begin
         location := locations[i];
         if location < min then min := location;
         if location > max then max := location;
      end;
      if min = +maxint then begin
         writeln(output,'readadefinition: program error at min');
         halt
      end;
      if max = -maxint then begin
         writeln(output,'readadefinition: program error at max');
         halt
      end;

      (* determine if this is a Gap or total type feature *)
      IsGapOrIsTotal(d);

   end;
{
writeln(output,'-----');
showadefinition(list,d);
writeln(output,'-----');
}
end;
(* end module readadefinition version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module writeNOTABstring *)
procedure writeNOTABstring(var tofile: text; var s: string);
(* write the string s to file tofile, no writeln, remove
all tab characters first. *)
var   notabs: string; (* no tab version of s *)
begin (* writeNOTABstring *)
   copystring(s,notabs);
   detabstring(notabs);
   writestring(tofile,notabs);
end; (* writeNOTABstring *)
(* end module writeNOTABstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module writeafeature *)
procedure writeafeature(var afile: text; f: featuretypeptr);
(* write the feature f to file afile in the format:
@ K00000 36 +1 "nametag" "othertag" 12.200338 -0.473212 0.318031 -26
where the last three numbers are the Ri, Z and p
followed by the relative coordinate [2016 aug 18].
Do not carriage return.
*)
var
   o: featurecarrierptr; (* other or mother features that this one
      points to *)
begin
   if f = nil then begin
      writeln(output, 'writeafeature: nil feature pointer!');
      halt;
   end;
   with f^ do begin
      write(afile,'@ ');
      writeNOTABstring(afile,id);
      write(afile,' ',coordinate:1:1);
      write(afile,' ');
      if orientation > 0 then write(afile,'+',orientation:1)
                         else write(afile,    orientation:1);
      write(afile,' ');
      writequotestring(afile,nametag);
      write(afile,' ');
      writequotestring(afile,othertag);
      if definition = nil then begin
{
         writeln(output); (* break away from the previous line *)
         writeln(output, 'writeafeature: nil definition pointer');
         writeln(output, 'for feature:');
         writeNOTABstring(output,id);
         write(output,' ');
         writequotestring(output,nametag);
         write(output,' ');
         writequotestring(output,othertag);
         writeln(output);
}
         write(afile,' ',0.0:infofield:infodecim);
         write(afile,' ',0.0:infofield:infodecim);
         write(afile,' ',0.0:infofield:infodecim);
         write(afile,' ',0.0:3:1); (* relative coordinate *)
      end
      else begin
         if definition^.matrix <> nil then begin
{
            write(afile,' ',Ri:infofield:infodecim);
            write(afile,' ',Z:infofield:infodecim);
            write(afile,' ',probability:infofield:infodecim);
}
            write(afile,' ',Ri:1:infodecim);
            write(afile,' ',Z:1:infodecim);
            write(afile,' ',probability:1:infodecim);
         end
         else begin (* parameters, new as of 2004 July 17 *)
            write(afile,' ',Tparam);
{
            write(afile,' ',Aparam:infofield:infodecim);
            write(afile,' ',Bparam:infofield:infodecim);
            write(afile,' ',Cparam:infofield:infodecim);
            write(afile,' ',Dparam:infofield:infodecim);
}
            write(afile,' ',Aparam:1:infodecim);
            write(afile,' ',Bparam:1:infodecim);
            write(afile,' ',Cparam:1:infodecim);
            write(afile,' ',Dparam:1:infodecim);
         end;
         write(afile,' ',relcoo:1:1); (* relative coordinate *)
         write(afile,' #',featurenumber:1); (* feature number *)
      end;
      if others <> nil then begin
         if others^.afeature <> nil then write(afile,' > ');
         o := others;
         while o <> nil do begin
            if o <> others then write(afile,', ');
            if o^.afeature <> nil then begin
               write(afile,'[');
               writeNOTABstring(afile,o^.afeature^.nametag);
               write(afile,']');
            end;
            o := o^.next;
         end;
      end;
   end
end;
(* originally from 9.29; (@ of lister.p 1999 July 5 *)
(* end module writeafeature version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module scan.writeparameters *)
procedure writeparameters(var afile: text; parameters: allparameters);
(* write the parameters to a file *)
var
   i: integer; (* index *)
begin with parameters do begin
   writeln(afile,'****************************************************');
   writeln(afile,'* Parameters used by scan: ');
   writeln(afile,'* todo:             ',todo:1);
   writeln(afile,'* cutDNRi:         ',cutDNRi:1:1);
   writeln(afile,'* cutUPRi:         ',cutUPRi:1:1);
   writeln(afile,'* cutDNZ:          ',cutDNZ:1:1);
   writeln(afile,'* cutUPZ:          ',cutUPZ:1:1);
   writeln(afile,'* cutDNP:          ',cutDNP:1:1);
   writeln(afile,'* cutUPP:          ',cutUPP:1:1);
   writeln(afile,'* fromwanted:       ',fromwanted:1);
   writeln(afile,'* towanted:         ',towanted:1);
   writeln(afile,'* ways:             ',ways:1);
   write  (afile,'* sitedefinition:   ');
   if sitedefined
   then writeadefinition(afile,sitedefinition)
   else writeln(afile,'(none)');

   for i := 1 to othertypes do begin
      case othercontrol.letters[i] of
         ' ': begin
            writeln(afile,'* (no control)');
         end;
         '#': begin
            writeln(afile,'* # number:',numberwidth:1);
         end;
         'l': begin
            writeln(afile,'* length:',lengthwidth:1);
         end;
         'n': begin
            writeln(afile,'* name:',namewidth:1);
         end;
         'f': begin
            writeln(afile,'* fullname:',fullnamewidth:1);
         end;
         'c': begin
            writeln(afile,'* coordinate:',coordinatewidth:1);
         end;
         'o': begin
            writeln(afile,'* orientation:',orientationwidth:1);
         end;
         'r': begin
            writeln(afile,'* Ri:',Riwidth:1,':',Ridecimal:1);
         end;
         'z': begin
            writeln(afile,'* Z:',Zwidth:1,':',Zdecimal:1);
         end;
         'p': begin
            writeln(afile,'* probability:',Pwidth:1,':',Pdecimal:1);
         end;
         's': begin
            write  (afile,'* string: "');
                   writequotestring(afile,otherquotestring);
            writeln(afile,'"');
         end;
         '.': begin
            writeln(afile,'* .: end of print definitions');
         end;
      end;
   end;

   write(afile,'* ');
   if outdata then write(afile,'D') else write(afile,'d');
   if outfeat then write(afile,'F') else write(afile,'f');
   if outinst then write(afile,'I') else write(afile,'i');
   writeln(afile,' DFI: capitals mean data, scanfeatures, scaninst written');

   writeln(afile,'* ',normalizeRi,
                 '  normalizeRi: n: normal, s: Ri(l)-Ri(0), d: Ri(l)/Ri(0)');

   writeln(afile,'* ',instfrom:1,' ',instto:1,
                 '  range to make the scaninst file (if made)');

   writeln(afile,'****************************************************');
end; end;
(* end module scan.writeparameters *)

(* begin module decapitalize *)
function decapitalize(c: char): char; 
(* convert the character c to lower case *) 
var 
   n: integer; (* c is the n'th letter of the alphabet *)
begin 
   n := ord(c); 
   if (n >= ord('A')) and 
      (n <= ord('Z'))   
   then c := chr( n - ord('A') + ord('a'))
   else c := chr(n);
   decapitalize := c
end; 
(* end module decapitalize version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

(* begin module multiscan.readparameters *)
procedure readparameters(var scanp: text;
                         riblmatrix: ribltypeptr;
                         var parameters: allparameters);
(* read from the scanp file the parameters *)
var
   checkout: boolean; (* if true, all variable values are ok *)
   l: integer; (* integer location of a mark *)
   i: integer; (* index to marks *)
   ot: char; (* other type control character *)

procedure cn;
begin
  checkout := checknumber(scanp);
  if not checkout then halt; (* avoid snowballing *)
end;

procedure nextlink(var pointer: stringptr);
(* Either start the prebuiltstring and its index pbs or make the next link.
Then make the pointer point to the new link*)
begin
   with parameters do begin
      if prebuiltstring = nil then begin
         new(prebuiltstring);
{
         initializestring(prebuiltstring^);
}
         clearstring(prebuiltstring^);
         pbs := prebuiltstring;
      end
      else begin
         new(pbs^.next);
         pbs := pbs^.next;
{
         initializestring(pbs^);
}
         clearstring(pbs^);
      end;
      pointer := pbs
   end
end;

begin with parameters do begin

   (* RYAN commented out the following so that it could read multiple paramter
      sets from 1 scanp file *)
   (* reset(scanp); *)

   readln(scanp, parameterversion);
   if parameterversion < updateversions then begin
      writeln(output, 'You have an old parameter file!');
      halt
   end;


   todo := maxint; (* do all sequences *)

   fromwanted := -maxint;
   towanted := +maxint;
   ways := 2;

   readln(scanp,todo);
   if todo < 0 then todo := maxint;

   if parameterversion <= 2.90 then begin
      cn; readln(scanp, cutDNRi);
      cutUPRi := +maxint;

      cutDNZ := 0;
      cn; readln(scanp, cutUPZ);

      cn; readln(scanp, cutDNP);
      cutUPP := 1;
   end
   else begin
      cn; readln(scanp, cutDNRi, cutUPRi);
      cn; readln(scanp, cutDNZ,  cutUPZ);
      cn; readln(scanp, cutDNP,  cutUPP);
   end;

   cn; read  (scanp, fromwanted);
   cn; readln(scanp, towanted);
   cn; readln(scanp, ways);
   skipblanks(scanp);
   if scanp^ = 'd'
   then begin
      sitedefined := true;
      new(sitedefinition);
      readadefinition(scanp, sitedefinition);
      (* check that the locations are inside the matrix *)
      with sitedefinition^, riblmatrix^ do begin
         for i := 1 to marks do begin
            (* convert the real number of the mark locations to an integer by
               rounding towards zero.  Thus a mark at -2.5 is allowed if the
               frombase is only -2 *)
            if locations[i] < 0
            then l := -trunc(abs(locations[i]))
            else l := +trunc(locations[i]);
            if (l < frombase) or (l > tobase)
            then begin
               writeln(output,'in scanp: mark definition (',l:1,
                              ') is out of matrix range ',
                              frombase:1,' to ',tobase:1);
               halt;
            end;
         end;
      end;
      (* insert the riblmatrix into the definition *)
      sitedefinition^.matrix := riblmatrix;
   end
   else begin
      sitedefined := false;
      sitedefinition := nil;
   end;
   readln(scanp);

   (* read the controls for the strings to be written *)
{
   initializestring(othercontrol);
}
   clearstring(othercontrol);
   othertypes := 0;
   prebuiltstring := nil;
   pbs := nil;
   stringforNumber := nil;
   stringforLength := nil;
   stringforName := nil;
   stringforFullName := nil;
   stringforCoordinate := nil;
   stringforOrientation := nil;
   stringforRi := nil;
   stringforZ := nil;
   stringforProbability := nil;
   stringforOther := nil;

   numberwidth := 0;
   lengthwidth := 0;
   namewidth := 0;
   fullnamewidth := 0;
   coordinatewidth := 0;
   orientationwidth := 0;
   Riwidth := 0;
   Zwidth := 0;
   Pwidth := 0;

{
   initializestring(otherquotestring);
}
   clearstring(otherquotestring);
   if sitedefined then repeat
      while eoln(scanp) do readln(scanp);
      if eof(scanp)
      then ot := ' '
      else read(scanp, ot);
      ot := decapitalize(ot);
      othercontrol.length := succ(othercontrol.length);
      if ot in ['#','l','n','f','c','o','r','z','p','s','.']
      then begin
         othercontrol.letters[othercontrol.length] := ot;
{
writeln(output);
writequotestring(output,othercontrol);
write(output,' ot = "',ot,'"');
writeln(output);
}
{zzz}

         case ot of

            '#': begin (* number *)
               if number > 0 then begin
                  writeln(output,'duplicate number definition');
                  halt
               end;
               skipnonblanks(scanp);
               readln(scanp, numberwidth);
               if numberwidth < 0 then begin
                  writeln(output,'numberwidth must be >= 0');
                  halt
               end;
               nextlink(stringforNumber);
{
               with pbs^ do begin
                  length := succ(length);
                  spotforNumber := length;
                  length := spotforNumber + numberwidth - 1;
zzz -1 ?
                  stringforNumber := pbs;
                  length := numberwidth;
               end;
}
            end;

            'l': begin (* length *)
               if lengthwidth > 0 then begin
                  writeln(output,'duplicate length definition');
                  halt
               end;
               skipnonblanks(scanp);
               readln(scanp, lengthwidth);
               if lengthwidth < 0 then begin
                  writeln(output,'lengthwidth must be >= 0');
                  halt
               end;
               nextlink(stringforLength);
{
               with pbs^ do begin
                  length := succ(length);
                  spotforLength := length;
                  length := spotforLength + lengthwidth - 1;
               end;
}
            end;

            'n': begin (* name *)
               if namewidth > 0 then begin
                  writeln(output,'duplicate name definition');
                  halt
               end;
               skipnonblanks(scanp);
               readln(scanp, namewidth);
               if namewidth < 0 then begin
                  writeln(output,'namewidth must be >= 0');
                  halt
               end;
               nextlink(stringforName);
{
               with pbs^ do begin
                  length := succ(length);
                  spotforName := length;
                  length := spotforName + namewidth - 1;
               end;
}
            end;

            'f': begin (* fullname *)
               if fullnamewidth > 0 then begin
                  writeln(output,'duplicate fullname definition');
                  halt
               end;
               skipnonblanks(scanp);
               readln(scanp, fullnamewidth);
               if fullnamewidth < 0 then begin
                  writeln(output,'fullnamewidth must be >= 0');
                  halt
               end;
               nextlink(stringforFullName);
{
               with pbs^ do begin
                  length := succ(length);
                  spotforName := length;
                  length := spotforName + fullnamewidth - 1;
               end;
}
            end;

            'c': begin (* coordinate *)
               if coordinatewidth > 0 then begin
                  writeln(output,'duplicate coordinate definition');
                  halt
               end;
               skipnonblanks(scanp);
               readln(scanp, coordinatewidth);
               if coordinatewidth < 0 then begin
                  writeln(output,'coordinatewidth must be >= 0');
                  halt
               end;
               nextlink(stringforCoordinate);
{
               with pbs^ do begin
                  length := succ(length);
                  spotforCoordinate := length;
                  length := spotforCoordinate + coordinatewidth - 1;
               end;
}
            end;

            'o': begin (* orientation *)
               if orientationwidth > 0 then begin
                  writeln(output,'duplicate orientation definition');
                  halt
               end;
               skipnonblanks(scanp);
               readln(scanp, orientationwidth);
               if orientationwidth < 0 then begin
                  writeln(output,'orientationwidth must be >= 0');
                  halt
               end;
               nextlink(stringforOrientation);
{zzzyyy}
            end;

            'r': begin (* Ri, individual information content, bits *)
               if Riwidth > 0 then begin
                  writeln(output,'duplicate Ri definition');
                  halt
               end;
               skipnonblanks(scanp);
               readln(scanp, Riwidth, Ridecimal);
               if Riwidth < 0 then begin
                  writeln(output,'Riwidth must be >= 0');
                  halt
               end;

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

               nextlink(stringforRi);
{
               with pbs^ do begin
                  length := succ(length);
                  spotforRi := length;
                  length := spotforRi + Riwidth - 1;
                  length := Riwidth;
                  stringforRi := pbs;
               end;
}
            end;

            'z': begin (* Z score *)
               if Zwidth > 0 then begin
                  writeln(output,'duplicate Z definition');
                  halt
               end;
               skipnonblanks(scanp);
               readln(scanp, Zwidth, Zdecimal);
               if Zwidth < 0 then begin
                  writeln(output,'Zwidth must be >= 0');
                  halt
               end;
               if Zdecimal < 1 then begin
                  writeln(output,'Zdecimal must be >= 1');
                  halt
               end;
               nextlink(stringforZ);
{
               with pbs^ do begin
                  length := succ(length);
                  spotforZ := length;
                  length := spotforZ + Zwidth - 1;
               end;
}
            end;

            'p': begin (* probability *)
               if Pwidth > 0 then begin
                  writeln(output,'duplicate probability definition');
                  halt
               end;
               skipnonblanks(scanp);
               readln(scanp, Pwidth, Pdecimal);
               if Pwidth < 0 then begin
                  writeln(output,'Pwidth must be >= 0');
                  halt
               end;
               if Pdecimal < 1 then begin
                  writeln(output,'Pdecimal must be >= 1');
                  halt
               end;
               nextlink(stringforProbability);
{
               with pbs^ do begin
                  length := succ(length);
                  spotforProbability := length;
                  length := spotforProbability + Pwidth - 1;
               end;
}
            end;

            's': begin (* string *)
               skipnonblanks(scanp);
               readquotestring(scanp,otherquotestring);
               readln(scanp);
               nextlink(stringforOther);
               copystring(otherquotestring,stringforOther^);
{
write(output,'stringforOther ');
writequotestring(output,stringforOther^);
writeln(output);
}

{
zzz
               with pbs^ do begin
                  length := succ(length);
                  spotforOther := length;
                  length := spotforOther + otherquotestring.length - 1;
                  for j := 1 to otherquotestring.length
                  do letters[j+spotforOther-1] := otherquotestring.letters[j];
               end;
               clearstring(otherquotestring);
zzz
}
            end;

            '.': begin (* end of print definitions *)
               readln(scanp);
            end;

         end;

{
writeln(output);
write(output,'pbs ');
writequotestring(output,pbs^);
writeln(output);

clearstring(otherquotestring);
joinstring(prebuiltstring,otherquotestring);
write(output,'prebuiltstring ');
writequotestring(output,otherquotestring);
clearstring(otherquotestring);
writeln(output);
}
      end
      else begin
         if not eof(scanp) then readln(scanp);
         othercontrol.letters[othercontrol.length] := ' ';
         if ot <> ' '
         then writeln(output,'unrecognized command in scanp: ',ot);
      end;
      othertypes := succ(othertypes);
   until eof(scanp) or (ot = '.');

   (* dfi *)
   if not eof(scanp) then begin
      if scanp^ = 'd' then outdata := false
                      else outdata := true;
      get(scanp);
      if scanp^ = 'f' then outfeat := false
                      else outfeat := true;
      get(scanp);
      if scanp^ = 'i' then outinst := false
                      else outinst := true;
      readln(scanp);
   end
   else begin
      outdata := true;
      outfeat := true;
      outinst := true;
   end;

   (* normalization *)
{zzzqqq this is not implemented}
   if not eof(scanp) then begin
      readln(scanp,normalizeRi);
      if (normalizeRi = 's') or
         (normalizeRi = 'd')
      then normalboolean := true
      else normalboolean := false
   end
   else begin
      normalizeRi := 'n';
      normalboolean := false
   end;

(*
if normalboolean
then writeln(output, 'normalboolean: TRUE')
else writeln(output, 'normalboolean: FALSE');
writeln(output, 'normalizeRi: ', normalizeRi);
*)

   (* normalization *)
   if not eof(scanp) then begin
      readln(scanp,instfrom, instto);
   end
   else begin
      instfrom := -100;
      instto   := +100;
   end;

   with parameters, riblmatrix^ do begin
      if (fromwanted < frombase) or (fromwanted > tobase)
      then begin
         writeln(output,'in scanp: range FROM wanted (',fromwanted:1,
                        ') is out of matrix range ',
                        frombase:1,' to ',tobase:1);
         halt;
      end;

      if (towanted < frombase) or (towanted > tobase)
      then begin
         writeln(output,'in scanp: range TO wanted (',towanted:1,
                        ') is out of matrix range ',
                        frombase:1,' to ',tobase:1);
         halt;
      end;
   end;

end; end;
(* end module multiscan.readparameters *)

(* begin module simpson *)
procedure simpson(upper : real; var answer,tol : real);
(* Perform a numerical integration of the normal distribution using
Simpson's rule.  The variable "lower" is the lower bound of the
integration.  The variable "upper" is the upper bound of the
integration.  "answer" is the value of the integration, and "tol"
is the tolerance of the integration procedure, and thus the error
of the calculation.  Written by Mark Shaner, 1992.  Taken from
Miller, Alan R. PASCAL PROGRAMS FOR SCIENTISTS AND ENGINEERS.
SYBEX Inc., 1981. pgs 260-291

2007 Mar 21:  Tom Schneider:  Using a tolerance of 1/maxint caused an
infinite loop on a 64 bit machine.  Presumably it was attempting to
compute the value to 20 or so decimal places.  To avoid this, the
tolerance is set to a more reasonable value.

*)

var
  i  (* a counter *): integer;

  x,x1,x2, (* the independent variables for calculating the value of
              functions *)
  pi, (* the value of pi *)
  val, (* the value of 1/sqrt(2*pi), it is used in calculating the value of the
          gaussian for any x value, and is defined as a variable in order to
          speed calculations. *)
  deltax, (* the distance between every x value and the subsequent one *)
  evensum, (* the sum of the area under each of the even numbered parabolas *)
  oddsum, (* the sum of the area under each of the odd numbered parabolas *)
  endsum, (* the sum of the area under the first and last parabolas *)
  endcor, (*  the value of the end correction, it is determined using dgauss *)
  sum1 : real; (* a place to store the previous sum so that it can be compared
                with the subsequent sum to determine if the tolerance level has
                been reached *)
  pieces : integer; (* the number of parabolas under the curve *)
  lower, (* the lower bound of the integration *)
  sum :real; (* the value of the area under the curve *)

begin
  pi := 4.0*arctan(1.0);
  val := 1/sqrt (2*pi);
  (* activate both files to be used *)
  lower := 0;
  upper := abs(upper);
{
  tol := 1/maxint;
}
  tol := 1.0E-6;
  pieces := 2;
  deltax := (upper-lower)/pieces;
  x := lower + deltax;
  oddsum := val*exp(-0.5*x*x);
  evensum := 0.0;
  x1 := lower;
  x2 := upper;
  endsum := val*exp(-0.5*x1*x1) + val*exp(-0.5*x2*x2);
  endcor := -(x1*val)*exp(-0.5*x1*x1) - (-(x2*val)*exp(-0.5*x2*x2));
  sum := (endsum + 4.0*oddsum)*deltax/3.0;

  repeat
    pieces := pieces*2;
    sum1 := sum;
    deltax := (upper - lower)/pieces;
    evensum := evensum + oddsum;
    oddsum := 0.0;
    for i := 1 to (pieces div 2) do
    begin
      x := lower + deltax*(2.0*i-1.0);
      oddsum := oddsum + val*exp(-0.5*x*x);
    end;
    sum := (7.0*endsum + 14.0*evensum + 16.0*oddsum + endcor*deltax)
          *deltax/15.0;
  until abs(sum-sum1) <= abs(tol*sum);
  answer := (0.5 - sum);
  if answer < 0.0 then answer := 0.0; (* safety for roundoff errors *)
(*
  writeln (upper:7:5,'      ',answer:7:5,'+/-',tol)
*)
end;
(* end module simpson version = 'matmod 2.07 2007mar21 tds/gds'; *)

function pietointR(p: real; pie: pieceptr): integer;
(* Convert the real number p from external piece coordinates
into internal 1 to n coordinates.  This is done using:
function pietoint(p: integer; pie: pieceptr): integer;
This routine just makes the many calls simpler.
*)
begin
{zzzppp}
   pietointr := pietoint(round(p), pie);
end;

(* begin module multiscan.scansequence *)
procedure scansequence(var {data, features,} inst: text;
                           {outdata, outfeat,} outinst: boolean;
                           riblmatrix: ribltypeptr;
                           fromwanted,towanted: integer;
                           apiece: pieceptr;

                           cutDNRi: real;
                           cutUPRi: real;
                           cutDNZ: real;
                           cutUPZ: real;
                           cutDNP: real;
                           cutUPP: real;

                       var sites, positionsevaluated: integer;
                           bothways: boolean;
                           mean, stdev: real;
                           sitedefined: boolean;
                           {sitedefinition: definitiontypeptr;}
                           parameters: allparameters;
                           thelength: integer;
                       var f: featuretypeptr;
                          initialover: boolean;
                       var freelist: featuretypeptr;
                       var freeliststart: featuretypeptr);
(* Scan the sequence apiece with the riblmatrix (using only the portion
fromwanted to towanted) and put the results to the three output files data,
features and inst if the corresponding booleans outdata outfeat and outinst
are turned on.

Old method:
Don't report sites below the cutoff Ri, above the absolute of
the cutoff Z score or below the cutoff probability P.

New method:  Don't report sites outside the given Ri, Z score or
probability P.

Report the number of
sites and the number of positions evaluated.  Scan the complementary strand
at the same time if bothways is true.  The routine is written to be
reasonably fast, so it uses a 'parameter' which is precalculated. *)
(* Note:  The start and endpoints of the scans are different for
the two orientations of asymmetric recognizers on linear pieces. *)

(* RYAN added 5 variables to scansequence, thelength, f, initialover, freelist,
   freeliststart. *)
var
   b: base; (* the base at position l *)
{zzzqqq i is needed}
   i: integer; (* internal coordinates for a piece, 1 to piecelength *)
   l: integer; (* standard coordinate of the base around aligned point *)
   lrequest: integer; (* requested position that needs to be wrapped
      back onto the piece in circular cases *)
   p: real; (* probability of result *)
   Ri: real; (* the value of the riblmatrix applied to the sequence at i *)
   Ric: real; (* Same as Ri, but for the complementary orientation *)
{not implemented:
   Ri0: real; (* Ri normalized by Ri at zero coordinate *)
}
   tol: real; (* tolerance in the resulting p *)
   startpoint: integer; (* the point on apiece to start scanning *)
   endpoint: integer; (* the point on apiece to end scanning *)
   Z: real; (* Z value of Ri with respect to mean and st. dev. *)

procedure writeitout(Ri: real; i: integer; orientation: integer);
(* write an Ri evaluation at coordinate i out if it passes the criterion.
specify the orientation *)
{2016 Sep 20:  TDS: apparently in this version it does NOT write
to a file, it puts the features that pass selections onto the linked
list.  So apparently this proceedure is now misnamed.}

const
    tolerance = 1e-6; (* how close Ri can be to the mean to
       consider it to be within 0 standard deviations of the mean *)
var thecoordinate: integer; (* the current external coordinate in the book *)
begin
   (* NB: if stdev is zero (that is there was no variation in the input
      sequences) then Z is 0.  The if avoids division by zero. *)
   if stdev > 0
   then Z := (Ri-mean)/stdev
   else begin
      if abs(Ri-mean) < tolerance
      then Z := 100.0
      else Z := 0.0;
   end;

(* XXXXXXXX *)

   (* compute the probability from the Z score *)
   if (abs(Z) > 0) and (abs(Z) < 9) then simpson(abs(Z), p, tol)
                                    else p := 0.0;
   if      (Ri >= cutDNRi) and (Ri     <= cutUPRi)
   and (abs(Z) >= cutDNZ)  and (abs(Z) <= cutUPZ)
   and      (p >= cutDNP)  and (p      <= cutUPP)
   then begin
{
writeln(output,'nBUBBA ++++++++++++++++ ',i:1,' inside +++++++ writeitout');
}
      sites := sites + 1;

      (* coordinate of the evaluation in book coordinates *)
      thecoordinate := inttopie(i, apiece);

{
      if outdata then begin
         if numbered
         then write(data, ' ', number:nfield, ' ')
         else write(data, ' (no.#) ');
         (* length of the sequence *)
         write(data, ' ', piecelength(apiece):nfield);
         (* name of the piece *)
         write(data, ' ');
         writename(data, apiece^.key.hea.keynam);
         write(data,' ',thecoordinate: nfield);
         (* orientation of the matrix *)
         write(data,' ',orientation:nfield);
         (* the evaluation *)
         write(data,' ',Ri:infofield:infodecim);
         write(data,' ',Z:infofield:infodecim);
         write(data,' ',p:infofield:infodecim);

         (* write normalized data column *)
         with parameters do begin
            if normalboolean
            then case normalizeRi of
               's': write(data,' ',(Ri-Ri0):infofield:infodecim);
               'd': write(data,' ',(Ri/Ri0):infofield:infodecim);
            end
            else write(data,' ',0:1);
         end;

         writeln(data);
      end;
}{zzzqqq NOT USED!!!}

      if sitedefined then with parameters do begin
(* Ryan added here, trying to get the features to add to a list *)

         if not initialover then begin
            new(f^.next);
            f:= f^.next;
            f^.next:= nil;
         end;

         if initialover then begin
            freelist:= freeliststart;
            if (freelist^.next^.next = nil) then begin
               new(f^.next);
               f:= f^.next;
               f^.next:= nil;
            end
            else begin
               f^.next:= freelist;
               f:= f^.next;

               freeliststart:= freeliststart^.next^.next;
               f^.next:= nil;

               freelist:= freeliststart;
               if freelist^.next = nil then begin
                  new(freelist^.next);
                  freelist:= freelist^.next;
                  freelist^.next:= nil;
               end;

            end;
         end;

         copynamestring(apiece^.key.hea.keynam, f^.id);

         f^.definition := sitedefinition;
         f^.nametag := sitedefinition^.nametag;

         (* fill up the feature and then write the feature out. *)
         f^.orientation := orientation;
         f^.coordinate := inttopie(i, apiece);

         (* construct the othertag string from its parts *)

         (* each of the stringfors is a pointer into the prebuiltstring *)
         if stringforNumber <> nil
         then numberrealstring(stringforNumber^, number, numberwidth, 0);
         if stringforLength <> nil
         then numberrealstring(stringforLength^, piecelength(apiece),
                           lengthwidth, 0);

         if stringforName <> nil
         then copynamestring(apiece^.key.hea.keynam,stringforName^);

         if stringforFullName <> nil
         then copylinestring(apiece^.key.hea.fulnam,stringforFullName^);

         if stringforCoordinate <> nil
         then numberrealstring(stringforCoordinate^, thecoordinate,
                           coordinatewidth, 0);
         if stringforOrientation <> nil then begin
            if orientationwidth = 1 then begin
               if orientation < 0 then stringforOrientation^.letters[1] := '-'
                                  else stringforOrientation^.letters[1] := '+';
               stringforOrientation^.length := 1;
            end
            else numberrealstring(stringforOrientation^, orientation,
                             orientationwidth, 0);
         end;

         if stringforRi <> nil
         then numberrealstring(stringforRi^, Ri, Riwidth, Ridecimal);
{zzzyyy}

         if stringforZ <> nil
         then numberrealstring(stringforZ^, Z, Zwidth, Zdecimal);
         if stringforProbability <> nil
         then numberrealstring(stringforProbability^, p, Pwidth, Pdecimal);

         clearstring(f^.othertag);
         joinstring(prebuiltstring, f^.othertag);

{
writequotestring(output,f^.othertag); writeln(output,' othertag');
}

{zzz}

         f^.Ri := Ri;
         f^.Z  := Z ;
         f^.probability  := p;

         f^.definition^.matrix := riblmatrix;

         (* 2005 Sep 1 TDS *)
         f^.others := nil;
         f^.mothers := nil;

(* 2007 Nov 30: gap distances are always integer because the zero
coordinates of ribl models are alway an integer.  So defang the
following code: *)
         if f^.definition^.matrix^.symmetry = 'e'
         then begin
{2016jun24 actually remove the old code!
            f^.evencoordinate := f^.coordinate + orientation*0.5;
}
{aaa}
{try to fix it:}
{removing the '+ orientation*0.5' does "solve" it: }
            f^.evencoordinate := f^.coordinate;
         end
         else f^.evencoordinate := f^.coordinate;

{
         if outfeat then begin
            writeafeature(features, f);
            writeln(features);
         end;
}
{TTT it's actually finding the sites but NOT writing them out}
(*
if outfeat then begin {nBUBBA}
writeln(output,'nBUBBA outfeat is true');
writeafeature(output, f); {nBUBBA}
writeln(output); {nBUBBA}
end; {nBUBBA}
*)

      end;

      (* if outinst then with parameters, f^ do begin *)
      if 1 = 2 then with parameters, f^ do begin
      (*
      name "dnaA"; piece K01789; get from 229 -100 to 229 +100 direction -;
      *)

         write(inst,'name ');
         (* 2007 Nov 27: Formerly it was the name of the ribl.
            That is not as useful as the name of the piece.  TDS.
            writequotestring(inst,nametag); *)
         write(inst,'"');
{qqq}
         writeline(inst, apiece^.key.hea.fulnam,false);
         write(inst,'"');
         write(inst,';');

         write(inst,' piece ');
{
         writename(inst,id);
}
         writestring(inst,id);
         write(inst,';');

         write(inst,' get from ',trunc(coordinate):1,' ');

         (* orientation relative to the original piece coordinates
            was determined outside this procedure *)
         if -orientation > 0
         then write(inst,'+',abs(instfrom) *(-orientation):1)
         else write(inst,    abs(instfrom) *(-orientation):1);

         write(inst,      ' to same '); (* simplify life! *)

         if orientation > 0 then write(inst,'+',abs(instto)   *orientation:1)
                            else write(inst,    abs(instto)   *orientation:1);

         write(inst,' direction ');
         if orientation > 0 then write(inst,'+;')
                         else write(inst,'-;');

         writeln(inst);
      end;

   end;
end;

begin (* scansequence *)
   sites := 0;
{writeln(output,'nBUBBA ++++++++++++++++++++++++++++++++++++ in scansequence');}

   (* build the feature *)
{
   if sitedefined then begin
      (* build as much of the feature to be written out as we can *)

      new(f);
}
{
      f^.id := apiece^.key.hea.keynam;
}
{
      copynamestring(apiece^.key.hea.keynam, f^.id);

      f^.definition := sitedefinition;
      f^.nametag := sitedefinition^.nametag;

   end;
}

  (* prepare for normalization *)
  with parameters do if normalboolean then begin
      if bothways then begin
         writeln(output, 'Scan cannot compute normalization when there are');
         writeln(output, 'two different orientations');
         halt
      end;

      (* check here that the piece has a zero.  It must be in this range:
         -frombase+1
         piecelength(apiece)-tobase
      *)
(* 2016 Jun 30 code removed - no longer makes sense
      if (0 < inttopie(-fromwanted+1,apiece)) or
         (0 > inttopie(piecelength(apiece), apiece)-towanted)
      then begin
         writeln(output, 'piece without zero base found; cannot normalize');
{zzzqqq but pieces don't have to have a zero!  this is the wrong language}
         halt
      end;
*)

{zzzqq apparently not used, can delete! unless variables are used later:
i - used elsewhere, must keep
Ri - next use is set to zero, not needed here
b - next use is set, not needed here
}
{ Can remove this code: it does nothing since Ri0 is not used!
      i := pietoint(0, apiece); (* zzzqqq this should not ever be needed!*)
      Ri := 0.0;
      for l := fromwanted to towanted do begin
         b := getbase(i, apiece); (*zzzqqq move up !*)
         Ri  := Ri  + riblmatrix^.data[b, l];
      end;
      Ri0 := Ri; (* capture zero base value *)
}

   end;

   (* scan the piece *)
   (* thelength := piecelength(apiece); *)
   if apiece^.key.piecon = linear then begin  (* LINEAR case *)
      if not initialover then begin
         startpoint := -fromwanted + 1;
         if startpoint < 1 then startpoint := 1;

         if piecelength(apiece)-towanted >= thelength
         then endpoint := thelength;
         if piecelength(apiece)-towanted < thelength
         then endpoint := piecelength(apiece)-towanted;

         (* 2005 Sep 1: bug - endpoint was not set in this case!
            writeln (output,'RYE DID NOT CONSIDER ALL CASES');
            2015 Mar 23.  Resolve the issue.
            if piecelength(apiece)-towanted = thelength then begin
            If towanted = 0, then
               piecelength(apiece) = thelength
            Both of the above cases give the same result for endpoint:
            endpoint := thelength;
            endpoint := piecelength(apiece);
            SO I changed:
            if piecelength(apiece)-towanted > thelength
            to
            if piecelength(apiece)-towanted >= thelength
         *)

{ Testing code for the above bug:
            writeln (output);
            writeln (output,'RYE DID NOT CONSIDER ALL CASES');
(*zzzggg*)
            write(output, 'FOR piece "');
            writename(output, apiece^.key.hea.keynam);
            writeln(output, '"');
            writeln(output, 'piecelength(apiece)-towanted = thelength');
            writeln (output,
              '[',piecelength(apiece):1,
             ' - ',towanted:1,
             ' = ',thelength:1,']');
            writeln(output,'thelength = ',
                            thelength:1);
            writeln(output,'piecelength(apiece) = ',
                            piecelength(apiece):1);
            writeln(output, 'fromwanted = ',fromwanted:1);
            writeln(output, '  towanted = ',  towanted:1);
(* SUGGESTED NEW CODE: *)
            endpoint := piecelength(apiece)-towanted-1;
(* or try this one: *)
            endpoint := thelength;
(* Which is right? ANSWER:  THE SECOND ONE! *)
         end;
}

      end;

      (* calculate the new startpoint and endpoint here *)
      if initialover then begin
         startpoint := thelength;
         endpoint := thelength;
         if endpoint + towanted > piecelength(apiece) then begin
            startpoint := piecelength(apiece)-towanted;
            endpoint := piecelength(apiece)-towanted;
         end;
      end;

      for i := startpoint to endpoint do begin
         Ri := 0.0;
         for l := fromwanted to towanted do begin
            b := getbase(i + l, apiece);
            Ri  := Ri  + riblmatrix^.data[b, l];
         end;
{
writeln(output,'nBUBBA +++++++++++++++++++++ inside +++++++ scansequence');
TTT
}
         if apiece^.key.piedir = apiece^.key.coodir
         then writeitout(Ri,i,+1)
         else writeitout(Ri,i,-1);
      end;

      if bothways then begin (* do complement of LINEAR *)
         if not initialover then begin
            startpoint := towanted + 1;
            if startpoint < 1 then startpoint := 1;

            if piecelength(apiece)+fromwanted > thelength
            then endpoint := thelength;
            if piecelength(apiece)+fromwanted < thelength
            then endpoint := piecelength(apiece)+fromwanted;

         end;

         (* calculate your new endpoint and startpoint here *)
         if initialover then begin
            startpoint := thelength;
            endpoint   := thelength;
            if endpoint - fromwanted > piecelength(apiece) then begin
               startpoint := piecelength(apiece) + fromwanted;
               endpoint   := piecelength(apiece) + fromwanted;
            end;
         end;

         for i := startpoint to endpoint do begin
            Ric := 0.0;
            for l := fromwanted to towanted do begin
               b := getbase(i - l, apiece);
               Ric := Ric + riblmatrix^.data[complement(b),  l]
            end;
         if apiece^.key.piedir = apiece^.key.coodir
         then writeitout(Ric,i,-1)
         else writeitout(Ric,i,+1);
         end;
      end;

      (* To calculate the positions evaluated, we find the region of the
         sequence which is used:
         positionsevaluated :=
            (piecelength(apiece)-towanted) - (-fromwanted+1) + 1;
         but this can be simplified: *)
      positionsevaluated := (piecelength(apiece) - towanted) + fromwanted;
      if positionsevaluated < 0 then positionsevaluated := 0;
   end  (* LINEAR case *)

   else begin (* handle CIRCULAR case *)
      if not initialover then begin
         startpoint := 1;
         endpoint := thelength;
      end;

      if initialover then begin
         startpoint:= thelength;
         endpoint:= thelength;
      end;

      for i := startpoint to endpoint do begin
         Ri := 0.0;
         Ric := 0.0;
         for l := fromwanted to towanted do begin
            b := getbase(i + l, apiece);
            Ri  := Ri  + riblmatrix^.data[b, l];
            if bothways then begin (* do complement *)
               lrequest := i - l;
               while lrequest > thelength do lrequest := lrequest - thelength;
               while lrequest < 1 do lrequest := lrequest + thelength;
               b := getbase(lrequest, apiece);
               Ric := Ric + riblmatrix^.data[complement(b),  l]
            end;
         end;
         if apiece^.key.piedir = apiece^.key.coodir
         then writeitout(Ri,i,+1)
         else writeitout(Ri,i,-1);
         if bothways then begin (* do complement *)
            if apiece^.key.piedir = apiece^.key.coodir
            then writeitout(Ric,i,-1)
            else writeitout(Ric,i,+1);
         end;
      end;
      positionsevaluated := piecelength(apiece);
   end;

   (* clear the feature *)


{
   if sitedefined then begin
      dispose(f);
   end;
}

end; (* scansequence *)
(* end module multiscan.scansequence *)

(* begin module multiscan.readdiparameters *)
procedure readdiparameters(var thefile: text;
                       var boxsize: integer;
                       var ricutoff: real;
                       var singleA: char;
                       var givegroupnumber: boolean;
                       var esscut: real);
(* Read the multiscan parameters from the parameter file *)
var
   parameterversion: real; (* parameter version number *)
   achar: char; (* a character read from thefile to determine givegroupnumber *)
begin
   reset(thefile);
   readln(thefile, parameterversion);

   if parameterversion < updateversionm then begin
      writeln(output, 'You have an old parameter file!');
      halt
   end;

   readln(thefile, boxsize); (* The range in which sites can be matched *)
   readln(thefile, ricutoff);  (* read ri cutoff value *)

   read(thefile, singleA);
   if not eoln(thefile) then begin
      read(thefile, achar);
      if achar = '#'
      then givegroupnumber := true
      else givegroupnumber := false
   end;
   readln(thefile);

   readln(thefile, esscut);
end;
(* end module multiscan.readdiparameters *)

(* begin module multiscan.distribution *)
procedure distribution(var thefile : text;
                       var distsize: integer;
                       var correction: real;
                       var store : map;
                       var lowhistog: integer;
                       var highhistog: integer);
(* read in distribution information from the file and store it in
   an array *)
var
   x: integer; (* general integer *)
   z: real; (* general real *)
   normvalue: real; (* the occurrence / sum *)
   sum: real; (* number of sites in distribution,
      generalized to any distribution. *)
   ln2: real; (* the natural log of 2 *)
   lowbound: real; (* the bottom of the histogram *)
   highbound: real; (* the top of the histogram *)
   numpos: integer; (* the number of positions *)
   c: char; (* first character on the line *)

begin
   read(thefile,c);
   while c = '*' do begin
      readln(thefile);
      read(thefile,c);
   end;

   lowbound:= 1000000;
   highbound:= -100000;

   numpos:= 0;
   while not eof(thefile) and (c <> '*') do begin
      (* read in the positions that make up range *)
      read(thefile,z);
      numpos:= numpos+1;

      distsize:= round(z);
      store[distsize].range:= distsize;
      readln(thefile,store[distsize].occur);

      if (store[distsize].range <= lowbound) and (store[distsize].occur <> 0)
         then lowbound := store[distsize].range;

      if (store[distsize].range >= highbound) and (store[distsize].occur <> 0)
         then highbound := store[distsize].range;

      (* read in the # of occurrences at each position *)
      if round(store[distsize].range) <> store[distsize].range then begin
         writeln(output, 'distribution values must be increments of 1');
         halt;
      end;

(* Now we can add a line between histog files, to make moving through files
   easier.*)

      if not eof (thefile) and not eoln(thefile) then read(thefile,c);
      if not eof(thefile) and eoln(thefile) then c := '*';
   end;

   sum := 0.0;
   for x:= round(lowbound) to round(highbound) do begin
      sum := sum + store[x].occur; (* the sum of occurrences at all postions *)
   end;

   if sum <= 0.0 then begin
      writeln(output,'the input histog distribution has only zero values!');
      halt;
   end;

   ln2 := ln(2);
   for x := round(lowbound) to round(highbound) do begin
      store[x].isitzero := false;
      normvalue := store[x].occur / sum; (* normalize distribution values *)
      if (normvalue = 0) then begin
          store[x].isitzero := true;
          store[x].probdist := gapkill;
      end
      else store[x].probdist := -1 * ln(normvalue) / ln2;
   end;

   (* calculate sample correction *)
   lowhistog:= round(lowbound);
   highhistog:= round(highbound);

   numpos:= highhistog - lowhistog + 1;

   correction := (numpos - 1 ) / ( 2 * ln( 2 ) * sum );

end;

(* end module multiscan.distribution *)

(* begin module multiscan.copyafeature *)
procedure copyafeature(var fone: featuretypeptr; ftwo: featuretypeptr);
(* Copy feature two over to feature one *)
const
  debugging = false; (* set true to show the copy results *)
begin
      fone^.id:= ftwo^.id;
      fone^.coordinate:= ftwo^.coordinate;
      fone^.orientation:= ftwo^.orientation;
      fone^.nametag:= ftwo^.nametag;
      fone^.othertag:= ftwo^.othertag;
(*       fone^.definition:= ftwo^.definition; *)
{
      fone^.definition^.nametag:= ftwo^.definition^.nametag;
      fone^.definition^.background:= ftwo^.definition^.background;
      fone^.definition^.negparts:= ftwo^.definition^.negparts;
      fone^.definition^.posparts:= ftwo^.definition^.posparts;
      fone^.definition^.locations:= ftwo^.definition^.locations;
      fone^.definition^.marks:= ftwo^.definition^.marks;
      fone^.definition^.min:= ftwo^.definition^.min;
      fone^.definition^.max:= ftwo^.definition^.max;
      fone^.definition^.number:= ftwo^.definition^.number;
      fone^.definition^.matrix:= ftwo^.definition^.matrix;
      fone^.definition^.RiboundFeature:= ftwo^.definition^.RiboundFeature;
      fone^.definition^.ZboundFeature:= ftwo^.definition^.ZboundFeature;
      fone^.definition^.PboundFeature:= ftwo^.definition^.PboundFeature;
}

      (* 2005 Sep 1: TDS asks: Why wasn't the definition copied? *)
      fone^.definition:= ftwo^.definition;

      fone^.Ri:= ftwo^.Ri;
      fone^.Z:= ftwo^.Z;
      fone^.probability:= ftwo^.probability;
      fone^.evencoordinate:= ftwo^.evencoordinate;
      fone^.unsatisfied:= ftwo^.unsatisfied;
      fone^.fromrange:= ftwo^.fromrange;
      fone^.torange:= ftwo^.torange;
      fone^.featurenumber:= ftwo^.featurenumber;
      fone^.desiredline:= ftwo^.desiredline;
{
      fone^.next:= ftwo^.next;
}

      fone^.others:= ftwo^.others; (* 2005 Sep 01 *)
      fone^.mothers:= ftwo^.mothers; (* 2005 Sep 01 *)

if debugging then begin
writeln(output,'copyafeature fone: ------');
writeafeature(output, fone); (* *)
writeln(output,'------'); (* *)
writeln(output,'copyafeature ftwo: ------');
writeafeature(output, ftwo); (* *)
writeln(output,'------'); (* *)
end;

end;
(* end module multiscan.copyafeature *)

(* begin module multiscan.clearmultifeatureptr *)
procedure clearmultifeatureptr(var mfp: multifeatureptr);
(* Clear variables of a multifeatureptr.  For safety's sake
and to prevent late nights for Tom, this must be done
after each new(mfp)
2016 Oct 05 *)
begin
   if mfp = nil then begin
      writeln(output,'clearmultifeatureptr: mfp is nil, cannot clear');
      halt;
   end;
   mfp^.m := nil;
   mfp^.mstart := nil;
   mfp^.totalinfo := 0;
   mfp^.writeitout := false;
   mfp^.next:= nil;
end;
(* end module multiscan.clearmultifeatureptr *)

(* begin module checksite *)
procedure checksite(var mlist: multitypeptr; var goodsite: boolean;
                    apiece: pieceptr);
(* Using the mlist's parameters determine whether you have a match *)

var
  coorddist: real; (* the distance between two coordinates *)

begin

   goodsite:= false;

(* check to see if spacing is good *)
   if mlist^.parentptr^.f <> nil then begin
{
      if mlist^.f^.orientation > 0 then begin
}
      if true then begin {TTT nBUBBA THIS WORKED!}
{TTT nBUBBA THIS WORKED! based on the idea that orientation should not matter}
{writeln(output,'nBUBBA checksite mlist^.f^.orientation > 0');}

{zzzqqqkkk
old:
          coorddist:=
          mlist^.f^.coordinate - mlist^.parentptr^.f^.coordinate;
}
          coorddist :=
                    pietointR(mlist^.f^.coordinate, apiece)
                    -
                    pietointR(mlist^.parentptr^.f^.coordinate, apiece);

       end
       else begin
{zzzqqqkkk here too:}
{old:
          coorddist :=
          mlist^.parentptr^.f^.coordinate - mlist^.f^.coordinate;
}
{writeln(output,'nBUBBA checksite mlist^.f^.orientation <= 0');}
          coorddist :=
                    pietointR(mlist^.parentptr^.f^.coordinate, apiece)
                    -
                    pietointR(mlist^.f^.coordinate, apiece);
       end;
   end
{
   else coorddist:= 10000000;
}
   else coorddist:= WEIRD; (* originally 10000000; *)
{writeln(output,'nBUBBA checksite coorddist=',coorddist:1:1);}

   mlist^.coorddist:= coorddist;

   if (coorddist >= mlist^.lowhistog) and (coorddist <= mlist^.highhistog) and
      (mlist^.parentptr^.f <> nil) then goodsite:= true;


(* check the orientation *)

   if mlist^.orienttog = 's' then begin
      if mlist^.f^.orientation <> mlist^.parentptr^.f^.orientation then
      goodsite:= false;
   end;

   if mlist^.orienttog = 'd' then begin
      if mlist^.f^.orientation = mlist^.parentptr^.f^.orientation then
      goodsite:= false;
   end;

end;
(* end module checksite *)

(* begin module copymlist *)
procedure copymlist (mlist, mliststart: multitypeptr;
                     var flist: multifeatureptr; (* features list *)
                     totalinfo: real);
(* copy the mlist info into the flist^.m *)

begin

   mlist:= mliststart;
   flist^.m:= flist^.mstart;
   flist^.totalinfo:= totalinfo;
   while mlist <> nil do begin

      if mlist^.writeitout then begin
         copyafeature(flist^.m^.f, mlist^.f);
      end;

      flist^.m^.writeitout:= mlist^.writeitout;
      if (flist^.m^.writeitout) then begin
         flist^.m^.gapsurp := mlist^.gapsurp;
         flist^.m^.coorddist := mlist^.coorddist;
         flist^.m^.posint:= mlist^.posint;
         if mlist^.parentptr <> nil then begin
            flist^.m^.parentevencoord:= mlist^.parentptr^.f^.coordinate;
            flist^.m^.parentorient:= mlist^.parentptr^.f^.orientation;
         end;
      end;

      mlist:= mlist^.next;
      flist^.m:= flist^.m^.next;
   end;
end;
(* end module copymlist *)

(* begin module comparelist *)
procedure comparelist (var mlist, mliststart: multitypeptr;
                       var flist: multifeatureptr;
                       var fliststart: multifeatureptr;
                           totalinfo: real; singleA: char;
                       var freef: multifeatureptr;
                       var freefstart: multifeatureptr;
                           apiece: pieceptr);
(* Compare the mlist to the flist, to see if you should add as a new feature *)
var
   matchbool: boolean; (* was there a match in the list? *)
begin


(* to initialize the list *)
{zzzqqqkkk
   if (flist^.totalinfo = 0) and (flist^.m^.f^.coordinate = -1000) then begin
}
{writeln(output,'BUBBA comparelist: 1');}
{writeln(output,'BUBBA comparelist: 1 flist^.totalinfo=', flist^.totalinfo:1:1);}
{writeln(output,'BUBBA comparelist: 1 flist^.m^.f^.coordinate =', flist^.m^.f^.coordinate:1:1);}

{
Originally -1000
      and (pietointR(flist^.m^.f^.coordinate,apiece) = WEIRD) then begin
This is not the right thing to do since WEIRD is a test for being
way out of the piece - to initialize things.  So test directly:
}
   if (flist^.totalinfo = 0)
      and (flist^.m^.f^.coordinate = WEIRD) then begin (* Originally -1000 *)
{HERE}

{writeln(output,'BUBBA comparelist: 2');}
      flist^.totalinfo:= 1;

      mlist:= mliststart;
      mlist:= mlist^.next;
      flist^.mstart:= flist^.m;
      while mlist <> nil do begin
         new(flist^.m^.next);
         new(flist^.m^.next^.f);
         flist^.m:= flist^.m^.next;
         flist^.m^.next:= nil;
         mlist:= mlist^.next;
      end;
{writeln(output,'BUBBA comparelist: 3');}

      copymlist(mlist, mliststart, flist, totalinfo);

      flist^.m:= flist^.mstart;
   end

   else begin

      matchbool:= false;

      flist:= fliststart;
      mlist:= mliststart;

      if singleA = 'f' then begin
         while (flist <> nil) and not matchbool do begin
            flist^.m:= flist^.mstart;
            mlist:= mliststart;
            matchbool:= true;
            while (flist^.m <> nil) and not (flist^.m^.writeitout) do begin
               flist^.m:= flist^.m^.next;
            end;
            while (mlist <> nil) and not (mlist^.writeitout) do begin
               mlist:= mlist^.next;
            end;
            if mlist^.writeitout <> flist^.m^.writeitout then matchbool:= false;
            (* if flist^.m^.f^.id = mlist^.f^.id then begin*)
            if equalstring(flist^.m^.f^.id, mlist^.f^.id) then begin
               if mlist^.f^.coordinate <> flist^.m^.f^.coordinate then
                  matchbool:= false;
               if (mlist^.f^.coordinate = flist^.m^.f^.coordinate) and
                  (mlist^.f^.orientation <> flist^.m^.f^.orientation) then
                  matchbool:= false;
               if (mlist^.f^.coordinate = flist^.m^.f^.coordinate) and
                  (mlist^.f^.orientation = flist^.m^.f^.orientation) then begin
                  if totalinfo > flist^.totalinfo then begin
                     copymlist(mlist, mliststart, flist, totalinfo);
                     matchbool:= true;
                  end
                  else matchbool:= true;
               end;
            end
            else matchbool:= false;
            flist:= flist^.next;
         end;
      end;

      flist:= fliststart;
      mlist:= mliststart;

(* check to see if this feature is already on the list *)
      while (flist <> nil) and not matchbool do begin
         flist^.m:= flist^.mstart;
         mlist:= mliststart;
         matchbool:= true;

         while flist^.m <> nil do begin
            if mlist^.writeitout <> flist^.m^.writeitout then matchbool:= false;
            if (mlist^.writeitout) and (flist^.m^.writeitout) then begin
               if mlist^.f^.coordinate <> flist^.m^.f^.coordinate then
               matchbool:= false;
            end;

            flist^.m:= flist^.m^.next;
            mlist:= mlist^.next;
         end;

         flist:= flist^.next;
      end;

(* if it is not, you need to add a new node to the list *)


      if not matchbool then begin
         flist:= fliststart;
         while flist^.next <> nil do flist:= flist^.next;

         freef:= freefstart;
         if freef^.next <> nil then while freef^.next^.next <> nil do
         freef:= freef^.next;

(* if there is a free node to use, then use it *)
         if freef^.next <> nil then begin
            flist^.next:= freef^.next;
            flist^.next^.next:= nil;
            flist:= flist^.next;
            freef^.next:= nil;
         end
(* if there is only one node on the freef list, you have to make a new list*)
         else begin
            new(flist^.next);
            flist:= flist^.next;
            flist^.next:= nil;

            new(flist^.m);
            flist^.mstart:= flist^.m;
            flist^.m^.next:= nil;

            new(flist^.m^.f);
            flist^.m^.f^.next:= nil;

            mlist:= mliststart;
            mlist:= mlist^.next;

            while mlist <> nil do begin
               new(flist^.m^.next);
               new(flist^.m^.next^.f);
               flist^.m:= flist^.m^.next;
               flist^.m^.next:= nil;
               mlist:= mlist^.next;
            end;
         end;

         copymlist(mlist, mliststart, flist, totalinfo);

         flist^.m:= flist^.mstart;
      end;

   end;

   flist:= fliststart;
   flist^.m:= flist^.mstart;
   mlist:= mliststart;

end;
(* end module comparelist *)

(* begin module overlapcheck *)
procedure overlapcheck(var mlist, mliststart: multitypeptr;
                       var killit: boolean;
                           apiece: pieceptr);
(* This procedure checks for overlap in sites and removes sites based
   on the parameters specified in multiscanp.  If there is a filter based
   on overlap, essential sites get priority over non-essential sites.
   If 2 essential sites overlap, and overlapping is not allowed, the site
   is killed.  If 2 non-essential sites overlap, the strongest is kept. *)

var
   checklist: multitypeptr; (* a pointer into the mlist to check overlaps *)
   overbool: boolean; (* does one site overlap another? true means yes *)
   overdist: real; (* distance between two coordinates *)
   chigh, clow: real; (* the coordinates of the first and last base in the
                         checklist site *)
   mhigh, mlow: real; (* the coordinates of the first and last base in the
                         mlist site *)
   internalcoo: real; (* internal coordinate computed once *)

begin
   checklist:= mliststart;
   while checklist <> nil do begin
      if (checklist^.oltog = 'n') or (checklist^.oltog = 'd') then begin
         if checklist^.writeitout then begin
            mlist:= mliststart;
            while mlist <> nil do begin
               if mlist^.writeitout then begin
                  if mlist^.posint <> checklist^.posint then begin

                     if checklist^.oltog = 'd' then begin
                        overbool:= false;
                        overdist:=
{zzzqqqkkk
                           abs(checklist^.f^.coordinate - mlist^.f^.coordinate);
}
                           abs(pietointR(checklist^.f^.coordinate,apiece)
                               -
                               pietointR(mlist^.f^.coordinate,apiece));
                        if overdist < checklist^.oldiff then overbool:= true;
                        if overbool then begin
                           if (mlist^.usetog = 'n') and
                           (checklist^.usetog = 'e') then
                           mlist^.writeitout:= false;

                           if (mlist^.usetog = 'e') and
                           (checklist^.usetog = 'e') then begin
                              mlist^.writeitout:= false;
                              killit:= true;
                           end;
                           if (mlist^.usetog = 'n') and
                           (checklist^.usetog = 'n') then begin
                              if (mlist^.f^.Ri - mlist^.gapsurp) <=
                              (checklist^.f^.Ri - checklist^.gapsurp) then
                              mlist^.writeitout:= false;
                           end;
                        end;
                     end;

                     if checklist^.oltog = 'n' then begin
                        overbool:= false;
{zzzqqqkkk
                        chigh:= checklist^.f^.coordinate +
                        (checklist^.f^.orientation *
                         checklist^.parameters.towanted);

                        clow:= checklist^.f^.coordinate +
                        (checklist^.f^.orientation *
                         checklist^.parameters.fromwanted);

                        mhigh:= mlist^.f^.coordinate +
                        (mlist^.f^.orientation * mlist^.parameters.towanted);

                        mlow:= mlist^.f^.coordinate +
                        (mlist^.f^.orientation * mlist^.parameters.fromwanted);
}

{
                        chigh:= pietointR(checklist^.f^.coordinate, apiece) +
                        (checklist^.f^.orientation *
                         checklist^.parameters.towanted);

                        clow:= pietointR(checklist^.f^.coordinate, apiece) +
                        (checklist^.f^.orientation *
                         checklist^.parameters.fromwanted);

                        mhigh:= pietointR(mlist^.f^.coordinate, apiece) +
                        (mlist^.f^.orientation * mlist^.parameters.towanted);

                        mlow:= pietointR(mlist^.f^.coordinate, apiece) +
                        (mlist^.f^.orientation * mlist^.parameters.fromwanted);
}

{zzzqqqnnn}
                        internalcoo := pietointR(checklist^.f^.coordinate, apiece);
                        chigh:= internalcoo +
                              (checklist^.f^.orientation *
                               checklist^.parameters.towanted);
                        clow:= internalcoo +
                              (checklist^.f^.orientation *
                               checklist^.parameters.fromwanted);

                        internalcoo := pietointR(mlist^.f^.coordinate, apiece);
                        mhigh:= internalcoo +
                              (mlist^.f^.orientation*mlist^.parameters.towanted);
                        mlow:= internalcoo +
                             (mlist^.f^.orientation*mlist^.parameters.fromwanted);

                        if mlist^.f^.orientation > 0 then begin
                           if (chigh <= mhigh) and (chigh >= mlow) then
                           overbool:= true;
                           if (clow <= mhigh) and (clow >= mlow) then
                           overbool:= true;
                        end;

                        if mlist^.f^.orientation < 0 then begin
                           if (chigh >= mhigh) and (chigh <= mlow) then
                           overbool:= true;
                           if (clow >= mhigh) and (clow <= mlow) then
                           overbool:= true;
                        end;

                        if overbool then begin
                           if (mlist^.usetog = 'n') and
                           (checklist^.usetog = 'e') then
                           mlist^.writeitout:= false;

                           if (mlist^.usetog = 'e') and
                           (checklist^.usetog = 'e') then begin
                              mlist^.writeitout:= false;
                              killit:= true;
                           end;
                           if (mlist^.usetog = 'n') and
                           (checklist^.usetog = 'n') then begin
                              if (mlist^.f^.Ri - mlist^.gapsurp) <=
                              (checklist^.f^.Ri - checklist^.gapsurp) then
                              mlist^.writeitout:= false;
                           end;
                        end;
                     end;
                  end;
               end;
               mlist:= mlist^.next;
            end;
         end;
      end;
      checklist:= checklist^.next;
   end;
end;

(* end module overlapcheck *)



(* begin module calcinfo *)
procedure calcinfo(mlist, mliststart: multitypeptr;
                   var flist: multifeatureptr;
                   var fliststart: multifeatureptr;
                       ricutoff: real;
                       singleA: char;
                   var freef: multifeatureptr;
                   var freefstart: multifeatureptr;
                       apiece: pieceptr);
(* This programs checks other parameters, calculates the total information
   and puts the features in an feature list *)

var
   killit: boolean; (* If true, you do not keep the site *)
   totalinfo: real; (* the total information for a site *)
   gapsurp: real; (* the gap surprisal of a site *)
   info: real; (* just a generic number to use to check things *)

begin
{writeln(output,'BUBBA calcinfo: 1');}
   killit:= false;
   mlist:= mliststart;
   while (mlist <> nil) and not killit do begin
{writeln(output,'BUBBA calcinfo: 2');}

(* check to see if the essentials have been met *)
      if (mlist^.f = nil) and (mlist^.usetog = 'e') then killit:= true;
{writeln(output,'BUBBA calcinfo: 3');}
      if (mlist^.f = nil) then mlist^.writeitout:= false;
{writeln(output,'BUBBA calcinfo: 4');}

      if mlist^.f <> nil then begin
{zzzqqqkkk this makes no sense!???
          if mlist^.f^.coordinate < 0 then mlist^.writeitout:= false;
}
{
          if pietointR(mlist^.f^.coordinate,apiece) < 0
}
          if mlist^.f^.coordinate = WEIRD (* I'M NOT SURE THIS IS RIGHT! *)
          then mlist^.writeitout:= false;
      end;

{writeln(output,'BUBBA calcinfo: 5');}
      if (mlist^.f <> nil) and not killit then begin
         if mlist^.gaptog = 'g' then begin
            gapsurp:= mlist^.store[round(mlist^.coorddist)].probdist
            + mlist^.correction;
            mlist^.gapsurp:= gapsurp;
         end
         else begin
            gapsurp:= 0;
            mlist^.gapsurp:= gapsurp;
         end;

{writeln(output,'BUBBA calcinfo: 6');}
         info:= mlist^.f^.Ri - gapsurp;
{writeln(output,'BUBBA calcinfo: 7');}

         if (mlist^.esscut >= 0) and (mlist^.usetog = 'n') then begin
{writeln(output,'BUBBA calcinfo: 8');}

            if (info <= mlist^.esscut) and (mlist^.usetog = 'n')
            then mlist^.writeitout:= false;

            if (info >= mlist^.esscut) then mlist^.writeitout:= true;
         end
         else begin
{writeln(output,'BUBBA calcinfo: 9');}
            if (info <= 0) and (mlist^.usetog = 'n')
            then mlist^.writeitout:= false
            else mlist^.writeitout:= true;
{writeln(output,'BUBBA calcinfo: 10');}
         end;
      end;

{writeln(output,'BUBBA calcinfo: 11');}
      mlist:= mlist^.next;
{writeln(output,'BUBBA calcinfo: 12');}
   end;

   mlist:= mliststart;
   while mlist <> nil do begin
      if mlist^.f <> nil then begin
{zzzqqqkkk
          if mlist^.f^.coordinate < 0 then mlist^.writeitout:= false;
}
{
          if pietointR(mlist^.f^.coordinate, apiece) < 0
NOWWWWW
}
          if mlist^.f^.coordinate = WEIRD
          then mlist^.writeitout:= false;
      end;
      mlist:= mlist^.next;
   end;
{writeln(output,'BUBBA calcinfo: 13');}

(* filter based on overlap parameters *)
   overlapcheck(mlist, mliststart, killit, apiece);
{writeln(output,'BUBBA calcinfo: 14');}


(* calc the total information *)
   mlist:= mliststart;
   totalinfo:= 0;
   while mlist <> nil do begin
      if mlist^.writeitout then begin
         totalinfo:= totalinfo + mlist^.f^.Ri;
         totalinfo:= totalinfo - mlist^.gapsurp;
      end;
      mlist:= mlist^.next;
   end;
{writeln(output,'BUBBA calcinfo: 15');}

   if totalinfo < ricutoff then killit:= true;
   mlist:= mliststart;
   mlist^.totalinfo:= totalinfo;
{writeln(output,'BUBBA calcinfo: 16');}


   if not killit then begin
      mlist:= mliststart;
      flist:= fliststart;
{{writeln(output,'BUBBA calcinfo: 17');}
      comparelist(mlist, mliststart, flist, fliststart, totalinfo,
                  singleA, freef, freefstart, apiece);
{writeln(output,'BUBBA calcinfo: 18');}
   end;
{writeln(output,'BUBBA calcinfo: 19');}

end;
(* end module calcinfo *)

(* begin module multiscan.paramcheck *)
procedure paramcheck(var mlist, mliststart: multitypeptr;
                     var sitefound: boolean;
                     var dogo: boolean; var flist: multifeatureptr;
                     var fliststart: multifeatureptr;
                     ricutoff: real; singleA: char;
                     var freef: multifeatureptr;
                     var freefstart: multifeatureptr;
                     apiece: pieceptr);

(* This procedure uses the paramters from multiscanp contained in the mlist
   to determine whether there is a site in the current loop and whether the
   program should move into the next loop *)

var
   goodsite: boolean; (* is it a goodsite? *)

begin

{writeln(output,'BUBBA paramcheck: 1');}

   dogo:= false;
   goodsite:= false;

   checksite(mlist,goodsite,apiece);
{writeln(output,'BUBBA paramcheck: 2');}

   if goodsite then begin

      sitefound:= true;
      if mlist^.position = 'm' then begin
         mlist^.fcurrent:= mlist^.f;

         mlist^.next^.doneloop:= false;
         mlist^.next^.f:= mlist^.fstart;
         dogo:= true;
      end;
{writeln(output,'BUBBA paramcheck: 3');}

      if mlist^.position = 'e' then begin
{writeln(output,'BUBBA paramcheck: 4');}
         mlist^.fcurrent:= mlist^.f;

{writeln(output,'BUBBA paramcheck: 5');}
         calcinfo(mlist, mliststart, flist, fliststart, ricutoff, singleA,
                  freef, freefstart, apiece);
{writeln(output,'BUBBA paramcheck: 6');}

         if mlist^.f = nil then mlist^.doneloop:= true;
      end;
   end;
end;
(* end module multiscan.paramcheck *)

(* begin module barname *)
procedure barname(var flist: multifeatureptr; apiece: pieceptr);
(* We want to determine the order of the sites, and write out the name
string based on the order *)
var
   counter: multifeatureptr;
   (* pointer into loop that determines order of sites *)
   x: integer; (* index to the numbar loop *)
   countbool: boolean; (* is the counter^.list = nil *)

   firstcounter: boolean; (* is it the first site in the counter *)

begin

(* Ok, you need some sort of boolean that will check if two numbers are the
same.  If they are they same you turn the boolean on, which will say, write
this number out twice.  No, maybe I need to count the occurrences of each
coordinate and then write that number out that many times *)

   new(counter);
   clearmultifeatureptr(counter); { FFF999 }
   firstcounter:= true;

   flist^.m:= flist^.mstart;
   if flist^.m^.f^.orientation > 0 then begin
      while flist^.m <> nil do begin
         if flist^.m^.writeitout then begin
            counter^.m:= flist^.mstart;
            flist^.m^.location:= 0;
            while counter^.m <> nil do begin
               if counter^.m^.writeitout then begin
{zzzqqqkkkk
                  if (flist^.m^.f^.coordinate > counter^.m^.f^.coordinate) then
                  flist^.m^.location:= flist^.m^.location + 1;
                  if (flist^.m^.f^.coordinate = counter^.m^.f^.coordinate)
}
                  if   pietointR(flist^.m^.f^.coordinate,apiece)
                     > pietointR(counter^.m^.f^.coordinate,apiece)
                  then flist^.m^.location:= flist^.m^.location + 1;

                  (* and (flist^.m^.f^.nametag <> counter^.m^.f^.nametag) and*)
                  if (  pietointR(flist^.m^.f^.coordinate,apiece)
                      = pietointR(counter^.m^.f^.coordinate,apiece))
                     and (not equalstring(  flist^.m^.f^.nametag,
                                          counter^.m^.f^.nametag))
                     and (flist^.m^.posint > counter^.m^.posint)
                  then flist^.m^.location:= flist^.m^.location + 1;

{
                  if (flist^.m^.f^.coordinate = counter^.m^.f^.coordinate)
}
                  (* and (flist^.m^.f^.nametag = counter^.m^.f^.nametag) then *)
                  if (  pietointR(flist^.m^.f^.coordinate,apiece)
                      = pietointR(counter^.m^.f^.coordinate,apiece))
                     and equalstring(flist^.m^.f^.nametag,
                                   counter^.m^.f^.nametag)
                  then flist^.m^.location:= flist^.m^.location + 1;

               end;
               counter^.m:= counter^.m^.next;
            end;
         end;
         if flist^.m^.writeitout = false then flist^.m^.location:= -10;
         flist^.m:= flist^.m^.next;
      end;
   end

   else begin
      while flist^.m <> nil do begin
         if flist^.m^.writeitout then begin
            counter^.m:= flist^.mstart;
            flist^.m^.location:= 0;
            while counter^.m <> nil do begin
               if counter^.m^.writeitout then begin
{zzzqqqkkk
                  if (flist^.m^.f^.coordinate < counter^.m^.f^.coordinate) then
}
                  if   pietointR(flist^.m^.f^.coordinate,apiece)
                     < pietointR(counter^.m^.f^.coordinate,apiece)
                  then flist^.m^.location:= flist^.m^.location + 1;

                  (* and (flist^.m^.f^.nametag <> counter^.m^.f^.nametag) and*)
                  if  (pietointR(flist^.m^.f^.coordinate,apiece)
                     = pietointR(counter^.m^.f^.coordinate, apiece))
                  and (not equalstring(flist^.m^.f^.nametag,
                           counter^.m^.f^.nametag))
                  and (flist^.m^.posint > counter^.m^.posint)
                  then flist^.m^.location:= flist^.m^.location + 1;

                  (* and (flist^.m^.f^.nametag = counter^.m^.f^.nametag) then *)
                  if (flist^.m^.f^.coordinate = counter^.m^.f^.coordinate)
                  and equalstring(flist^.m^.f^.nametag, counter^.m^.f^.nametag)
                  then flist^.m^.location:= flist^.m^.location + 1;

               end;
               counter^.m:= counter^.m^.next;
            end;
         end;
         if flist^.m^.writeitout = false then flist^.m^.location:= -10;
         flist^.m:= flist^.m^.next;
      end;
   end;

   flist^.m:= flist^.mstart;
   x:= 0;
   while flist^.m <> nil do begin
      if (flist^.m^.writeitout) then begin
         x:= x+1;
         counter^.m:= flist^.mstart;
         countbool:= false;
         while not countbool and (counter^.m^.location <> x) do begin
            if counter^.m^.next <> nil then begin
               counter^.m:= counter^.m^.next;
            end
            else countbool:= true;
         end;
         if not firstcounter then write(scanfeatures,'-');
         if firstcounter then firstcounter:= false;
         writequotestringB(scanfeatures,counter^.m^.f^.nametag);
      end;
      flist^.m:= flist^.m^.next;
   end;

   dispose(counter);
end;
(* end module barname *)

(* begin module barnumb *)
procedure barnumb(var scanfeatures: text;
                  var flist: multifeatureptr; lastcoord: real;
                  apiece: pieceptr);
(* Ryan's description: "We want to determine the order of the sites."
TDS: this is the site list at the end of a define feature.
Write the result into scanfeatures. *)
var
   counter: multifeatureptr; (* pointer used to determine order of sites *)
   x: integer; (* index to the numbar loop *)
   countbool: boolean; (* is the counter^.list = nil *)
begin

(* Ok, you need some sort of boolean that will check if two numbers are the
same.  If they are they same you turn the boolean on, which will say, write
this number out twice.  No, maybe I need to count the occurrences of each
coordinate and then write that number out that many times *)

   new(counter);
   clearmultifeatureptr(counter); { FFF999 }

   flist^.m:= flist^.mstart;
   if flist^.m^.f^.orientation > 0 then begin
      while flist^.m <> nil do begin
         if flist^.m^.writeitout then begin
            counter^.m:= flist^.mstart;
            flist^.m^.location:= 0;
            while counter^.m <> nil do begin
               if counter^.m^.writeitout then begin
{zzzqqqkkk
                  if (flist^.m^.f^.coordinate > counter^.m^.f^.coordinate) then
                  flist^.m^.location:= flist^.m^.location + 1;
}
                  if   pietointR(flist^.m^.f^.coordinate,apiece)
                     > pietointR(counter^.m^.f^.coordinate,apiece)
                  then flist^.m^.location:= flist^.m^.location + 1;

               (* and (flist^.m^.f^.nametag <> counter^.m^.f^.nametag) and *)
                  if (  pietointR(flist^.m^.f^.coordinate,apiece)
                      = pietointR(counter^.m^.f^.coordinate,apiece))
                     and (not equalstring (flist^.m^.f^.nametag,
                                         counter^.m^.f^.nametag))
                     and (flist^.m^.posint > counter^.m^.posint)
                  then flist^.m^.location:= flist^.m^.location + 1;

               (* and (flist^.m^.f^.nametag = counter^.m^.f^.nametag) then *)
                  if (  pietointR(flist^.m^.f^.coordinate,apiece)
                      = pietointR(counter^.m^.f^.coordinate,apiece))
                  and equalstring (flist^.m^.f^.nametag,
                                 counter^.m^.f^.nametag)
                  then flist^.m^.location:= flist^.m^.location + 1;

               end;
               counter^.m:= counter^.m^.next;
            end;
         end;
         if flist^.m^.writeitout = false then flist^.m^.location:= -10;
         flist^.m:= flist^.m^.next;
      end;
   end

   else begin
      while flist^.m <> nil do begin
         if flist^.m^.writeitout then begin
            counter^.m:= flist^.mstart;
            flist^.m^.location:= 0;
            while counter^.m <> nil do begin
               if counter^.m^.writeitout then begin
{zzzqqqkkk
                  if (flist^.m^.f^.coordinate < counter^.m^.f^.coordinate) then
                  flist^.m^.location:= flist^.m^.location + 1;
}
                  if   pietointR(flist^.m^.f^.coordinate,apiece)
                     < pietointR(counter^.m^.f^.coordinate,apiece)
                  then flist^.m^.location:= flist^.m^.location + 1;

               (* and (flist^.m^.f^.nametag <> counter^.m^.f^.nametag) and *)
                  if (flist^.m^.f^.coordinate = counter^.m^.f^.coordinate)
                     and not equalstring(flist^.m^.f^.nametag,
                                       counter^.m^.f^.nametag)
                     and (flist^.m^.posint > counter^.m^.posint)
                  then flist^.m^.location:= flist^.m^.location + 1;

               (* and (flist^.m^.f^.nametag = counter^.m^.f^.nametag) then *)
                  if (flist^.m^.f^.coordinate = counter^.m^.f^.coordinate)
                     and equalstring (flist^.m^.f^.nametag,
                                    counter^.m^.f^.nametag)
                  then flist^.m^.location:= flist^.m^.location + 1;

               end;
               counter^.m:= counter^.m^.next;
            end;
         end;
         if flist^.m^.writeitout = false then flist^.m^.location:= -10;
         flist^.m:= flist^.m^.next;
      end;
   end;

   flist^.m:= flist^.mstart;
   x:= 0;
   while flist^.m <> nil do begin
      if (flist^.m^.writeitout) then begin
         x:= x+1;
         counter^.m:= flist^.mstart;
         countbool:= false;
         while (not countbool) and (counter^.m^.location <> x) do begin
            if counter^.m^.next <> nil then begin
               counter^.m:= counter^.m^.next;
            end
            else countbool:= true;
         end;
         if counter^.m^.f^.orientation > 0 then
         write(scanfeatures,
{
               counter^.m^.f^.coordinate - lastcoord:10:0, ' ');
2012 Sep 21 LEAVE ALONE the number for total:
}
{zzzqqqkkk
               +(counter^.m^.f^.coordinate - lastcoord):10:0, ' ');
}
               +(pietointR(counter^.m^.f^.coordinate,apiece)
                          - lastcoord):1:1, ' ');

{2016 aug 04: this produced the E format:
define "p35-p10 9 total 2.2" "-" "||" "||" -2.200E+01  0.000E+00 
                          - lastcoord):10, ' ');
}

         if counter^.m^.f^.orientation < 0 then
         write(scanfeatures,
{
               -(counter^.m^.f^.coordinate - lastcoord):10, ' ');
2012 Sep 21 flip the number for total:
}
{zzzqqqkkk
               +(counter^.m^.f^.coordinate - lastcoord):10, ' ');
}
               +(pietointR(counter^.m^.f^.coordinate,apiece)
                - lastcoord):1:1, ' ');

{2016 aug 04: this produced the E format:
define "p35-p10 9 total 2.2" "-" "||" "||" -2.200E+01  0.000E+00 
                - lastcoord):10, ' ');
}

{ 2012 Sep 21:
conclusion: extra code here can be removed - not orientation dependent}

      end;
      flist^.m:= flist^.m^.next;
   end;

   dispose(counter);
end;
(* end module barnumb *)

(* begin module writebar *)
procedure writebar(mlist, mliststart: multitypeptr;
                   flist: multifeatureptr;
                   apiece: pieceptr;
                   givegroupnumber: boolean; groupnumber: integer);
{xxx}
(* Write out the tick bar that indicates the location of every site
and the total information for that site *)

var
   firstcoord: real; (* the first coordinate in the series of features,
        TDS: INTERNAL COORDINATES *)
   lastcoord: real; (* the last coordinate in the series of features,
        TDS: INTERNAL COORDINATES *)
   parcoord: real; (* the coordinate of the parent model
        TDS: INTERNAL COORDINATES *)

   done: boolean; (* finished searching the mliststart *)
   mlistFullName: multitypeptr; (* the location in mliststart of the
      full name *)

begin

(* If greatest coord of all, then set as firstcoord.  Then subtract everything
else from it *)

{define "total 6.04336" "-" "||" "||" -23.0 0.0
@ AL009126 203561.0 1 "total info = 6.04336" " "}

   firstcoord:= 1000000000;
   lastcoord:= -1000000000;

   flist^.m:= flist^.mstart;
{zzzqqqkkk
   if flist^.m <> nil then parcoord:= flist^.m^.f^.coordinate;
}
   if flist^.m <> nil then parcoord:= pietointR(flist^.m^.f^.coordinate,apiece);
   while flist^.m <> nil do begin
      if flist^.m^.writeitout then begin

{zzzqqqkkk
         if flist^.m^.f^.coordinate <= firstcoord then
            firstcoord:= flist^.m^.f^.coordinate;
         if flist^.m^.f^.coordinate >= lastcoord then
            lastcoord:= flist^.m^.f^.coordinate;
}
         if pietointR(flist^.m^.f^.coordinate,apiece) <= firstcoord
         then firstcoord:= pietointR(flist^.m^.f^.coordinate,apiece);
         if pietointR(flist^.m^.f^.coordinate,apiece) >= lastcoord
         then lastcoord:= pietointR(flist^.m^.f^.coordinate,apiece);

      end;

      flist^.m:= flist^.m^.next;
   end;


(* RYAN HERE *)

   (* write out the definition of the total for this multifeature*)
   flist^.m:= flist^.mstart;
   mlist:= mliststart;
   write(scanfeatures, 'define "');

   flist^.m:= flist^.mstart;
   barname(flist, apiece);
{
   flist^.m:= flist^.mstart;
   while (flist^.m <> nil) do begin
      if flist^.m^.writeitout then begin

         if flist^.m = flist^.mstart then begin
            writequotestringB(scanfeatures,flist^.m^.f^.nametag);
         end
         else begin
            write(scanfeatures,'-');
            writequotestringB(scanfeatures,flist^.m^.f^.nametag);
         end;

      end;

      flist^.m:= flist^.m^.next;
   end;
}

   write(scanfeatures,' ',inttopie(round(parcoord),apiece):1);

   (* 2007 Mar 13: If there is a fullname available, put it
      after the parcoord *)
{ original code:
   writeln(scanfeatures,' ',inttopie(parcoord,apiece):1,' total',
                        ' ',flist^.totalinfo:1:1,
                        '" "bits"',
   write  (scanfeatures, ' ',inttopie(parcoord,apiece):1);
}

   (* Search through all of the features in this mlist.
      If one of them has a stringforFullName, print it to the total
      'bar'.  Note that there can only be one such stringforFullName
      since the feature must reside on a particular piece.
      So we terminate the search after finding the first one.
      Note that mlistFullName is used later to write the Fullname
      to the @ line. *)
   mlistFullName:= mliststart;
   done := false;
   while (mlistFullName <> nil) and (not done) do begin
      if mlistFullName^.parameters.stringforFullName <> nil then begin
         write(scanfeatures, ' ');
         writestring(scanfeatures,
                     mlistFullName^.parameters.stringforFullName^);
         done := true
      end
      else mlistFullName := mlistFullName^.next;
   end;

   write(scanfeatures,' total ', flist^.totalinfo:1:1,'" "-" "');

   flist^.m:= flist^.mstart;
   while (flist^.m <> nil) do begin
      if flist^.m^.writeitout then begin
         write(scanfeatures,'|');
      end;

      flist^.m:= flist^.m^.next;
   end;

   write(scanfeatures,'" "');

   flist^.m:= flist^.mstart;
   while flist^.m <> nil do begin
      if flist^.m^.writeitout then begin
         write(scanfeatures,'|');
      end;
      flist^.m:= flist^.m^.next;
   end;
   write(scanfeatures,'" ');

   flist^.m:= flist^.mstart;
{writeln(scanfeatures,'ZZZQQQ-ALMOSTfinalpoint');}
   barnumb(scanfeatures, flist, lastcoord, apiece);
{QQQ}
{writeln(scanfeatures,'ZZZQQQyyy-finalpoint');}
   writeln(scanfeatures); (* end of the define --------------------------- *)

   (* Locate the last feature, the "total" feature *)
   flist^.m:= flist^.mstart;
   while (flist^.m^.writeitout <> true) and (flist^.m <> nil) do begin
      flist^.m:= flist^.m^.next;
   end;

   (* write the "total" feature *)
(* TODO NOTE:  this code should be replaced with a call to writeafeature !!! *)
{write(scanfeatures,'THIS_IS_IT_ZZZQQQ:');}
   write(scanfeatures,'@ ');
   writequotestringB(scanfeatures,flist^.m^.f^.id);
   (* 2012 Sep 21: bug was to use last coord here instead of parcoord:
   write(scanfeatures,' ',round(lastcoord):1,' ',
                    flist^.m^.f^.orientation:1, ' "');  {2004 aug 11}
   Orientation has to be set as the gap is set
ZZZQQQ
                  write(scanfeatures,' ',inttopie(parcoord,apiece):1,' Gap"');
                  if flist^.m^.f^.orientation > 0 then begin
                     write(scanfeatures,  ' "-" "{}" "{}" 0.0');
                     writeln(scanfeatures,' ',abs(flist^.m^.coorddist):1);
                  end;
                  if flist^.m^.f^.orientation < 0 then begin
                     write(scanfeatures,  ' "-" "}{" "}{" 0.0');
                     writeln(scanfeatures,' ',-abs(flist^.m^.coorddist):1);
                  end;
   *)
   write(scanfeatures,' ',inttopie(round(parcoord),apiece):1,' ');
   write(scanfeatures, flist^.m^.f^.orientation:1, ' "');

   flist^.m:= flist^.mstart;
   barname(flist, apiece);
{ old code:
   flist^.m:= flist^.mstart;
   while (flist^.m <> nil) do begin
      if flist^.m^.writeitout then begin
         if flist^.m = flist^.mstart then begin
            writequotestringB(scanfeatures,flist^.m^.f^.nametag);
         end
         else begin
            write(scanfeatures,'-');
            writequotestringB(scanfeatures,flist^.m^.f^.nametag);
         end;
      end;
      flist^.m:= flist^.m^.next;
   end;
}

   write  (scanfeatures, ' ',inttopie(round(parcoord),apiece):1);

   (* 2007 Mar 13: If there is a fullname available, put it out *)
   if mlistFullName <> nil then begin
      write(scanfeatures, ' ');
      writestring(scanfeatures, mlistFullName^.parameters.stringforFullName^);
   end;

{before #
   writeln(scanfeatures, ' total',
                        ' ',flist^.totalinfo:1:1,
                        '" "bits"',
                        ' ',flist^.totalinfo:1:infodecim,
                        ' ',0.0:3:1, (* 2016 aug 04 tighten the output *)
                        ' ',0.0:3:1, (* 2016 aug 04 tighten the output *)
                        ' ',0.0:3:1); (* 2016 aug 18 relcoo *)
}

   write(scanfeatures, ' total',
                       ' ',flist^.totalinfo:1:1,
                       '" "bits');
   if givegroupnumber
   then write(scanfeatures, '-#',groupnumber:1,'"')
   else write(scanfeatures, '"');
   write(scanfeatures, ' ',flist^.totalinfo:1:infodecim,
                       ' ',0.0:3:1, (* 2016 aug 04 tighten the output *)
                       ' ',0.0:3:1, (* 2016 aug 04 tighten the output *)
                       ' ',0.0:3:1); (* 2016 aug 18 relcoo *)
                       (* Note that relcoo = relative coordinate is zero
                          for the total by definition. *)
   writeln(scanfeatures);

{gggg}
{
                        ' ',0.0:1:infodecim,
                        ' ',0.0:1:infodecim);
}
{
and one should not be there:
                        ' ',flist^.totalinfo:1:infodecim,
                        ' ',mliststart^.maxtotal:1:infodecim);
who knows what the last one was!!!
}

end;
(* end module writebar *)

(* begin module multifeatwrite *)
procedure multifeatwrite(mlist, mliststart: multitypeptr;
                   var  flist, fliststart: multifeatureptr;
                   var totalsites: integer;  apiece: pieceptr;
                   var givegroupnumber: boolean;
                   var groupnumber: integer);
(* Sort features in order to filter, and write them out *)

var
   parcoord: real; (* the coordinate of the parent model INTERNAL coordinates *)

procedure writescaninst;
(* write out scaninst *)
begin
   (* qqq *)
   if (mliststart^.parameters.outinst) and (mlist = mliststart)
   then begin

      write(scaninst,'name ');
      (* 2007 Nov 27: Formerly it was the name of the ribl.
         That is not as useful as the name of the piece.  TDS.
         writequotestring(scaninst,nametag); *)
      write(scaninst,'"');
      writeline(scaninst, apiece^.key.hea.fulnam,false);
      write(scaninst,'"');
      write(scaninst,';');
      write(scaninst,' piece ');

      writestring(scaninst,flist^.m^.f^.id);
      write(scaninst,';');

      write(scaninst,' get from ',trunc(flist^.m^.f^.coordinate):1,' ');

      (* orientation relative to the original piece coordinates
         was determined outside this procedure *)
      if -flist^.m^.f^.orientation > 0
      then write(scaninst,'+',abs(mliststart^.parameters.instfrom)
                 *(-flist^.m^.f^.orientation):1)
      else write(scaninst,    abs(mliststart^.parameters.instfrom)
                 *(-flist^.m^.f^.orientation):1);

      write(scaninst,      ' to same '); (* simplify life! *)

      if flist^.m^.f^.orientation > 0 then write(scaninst,'+',
           abs(mliststart^.parameters.instto)   *flist^.m^.f^.orientation:1)
      else write(scaninst,    abs(mliststart^.parameters.instto)
           *flist^.m^.f^.orientation:1);

      write(scaninst,' direction ');
      if flist^.m^.f^.orientation > 0 then write(scaninst,'+;')
                      else write(scaninst,'-;');
      writeln(scaninst);

   end;
end; (* writescaninst *)

begin (* multifeatwrite *)
{writeln(output,'nBUBBA ---- multifeatwrite');}
{writeln(output,'BUBBA 10---- multifeatwrite');}
   flist:= fliststart;
   (* original code from Ryan does not handle missing sites:
   while flist <> nil do begin
wrong:
   while (flist <> nil) and (flist^.mstart <> nil) do begin
   *)
   while flist <> nil do begin (* 2016 oct 03 - I'm not going to
                                  add the extra indent for the next part *)
   if flist^.mstart <> nil then begin
{writeln(output,'BUBBA 11---- multifeatwrite');}
      totalsites:= totalsites + 1;
{if flist^.mstart = nil then writeln(output,'BUBBA 11.5 ---- multifeatwrite');}
      flist^.m:= flist^.mstart;
      mlist:= mliststart;
      if flist^.m <> nil then begin
{writeln(output,'BUBBA 12---- multifeatwrite sequence number = ', number:1);}
{zzzqqqkkk
         parcoord:= flist^.m^.f^.coordinate;
}
{writeln(output,'BUBBA 12.1 ---- multifeatwrite CRASH');}
{zzzqqqkkk}
{if flist = nil then writeln(output,'BUBBA 12.15 ---- multifeatwrite flist = nil')}
{               else writeln(output,'BUBBA 12.15 ---- multifeatwrite flist != nil');}

{if flist^.m = nil then writeln(output,'BUBBA 12.2 ---- flist^.m = nil')}
{                  else writeln(output,'BUBBA 12.2 ---- flist^.m <> nil');}
{writeln(output,'BUBBA 12.25 CRASH IS NEXT, number = ',number:1);}

(*
{CRASH ON THIS LINE: flist^.m^.f is pointing somewhere weird}
writeafeature(output,flist^.m^.f); writeln(output);
writeln(output,'BUBBA 12.30, past writeafeature, number = ',number:1); {BUBBA}
if flist^.m^.f = nil then writeln(output,'BUBBA 12.4 ---- flist^.m^.f nil');
if flist^.m^.f <>nil then writeln(output,'BUBBA 12.3 ---- flist^.m^.f not nil');
writeln(output,'BUBBA 12.42, number = ',number:1); { BUBBA }
{if number = 10 then halt { BUBBA } else writeln(output,'pass 10 halt');}
writeln(output,'BUBBA 12.45'); { BUBBA }
*)

         parcoord:= pietointR(flist^.m^.f^.coordinate,apiece);
{writeln(output,'BUBBA 12.5 ---- multifeatwrite');}
      end;
      writeln(data); (* 2013 Feb 07 blank START before each set of sites *)
      writeln(scanfeatures); (* 2004 aug 12: separate the parts *)
      groupnumber := groupnumber + 1;
{if groupnumber > 1 then begin writeln(output,'zzzooo 2 halt');halt;end;(*zzzooo*)}
      writeln(scanfeatures, 'group ', groupnumber:1);
{writeln(output,'BUBBA 13---- multifeatwrite');}
      while flist^.m <> nil do begin
{writeln(output,'BUBBA 14---- multifeatwrite');}
         if (flist^.m^.writeitout) then begin
            flist^.m^.f^.relcoo := flist^.m^.coorddist; (* 2016aug18 *) {gggg}
            writeafeature(scanfeatures, flist^.m^.f);
            writeln(scanfeatures);
            writequotestringB(data,flist^.m^.f^.nametag);
            writeln(data,round(flist^.m^.f^.coordinate):10,
                    flist^.m^.f^.orientation:5,
                    flist^.m^.f^.Ri:10:5, flist^.m^.coorddist:10:1,
                    flist^.m^.gapsurp:10:5);

            writescaninst;

(*XXXXXXXXXX*)

            if mlist^.parent <> 0 then begin
            if mlist^.bartog = 'b' then begin
               (* define the gap surprisal bar out *)
               if flist^.m^.f <> nil  then begin


{aaa}
(* write the Gap definition out **********************)
(* this is ok
write(output,'define "');
writequotestringB(output,flist^.m^.f^.nametag);
write(output,'-(',abs(flist^.m^.coorddist):1:0,')-');
writequotestringB(output,
                  mlist^.parentptr^.matrix^.riblname);
{zzzqqqkkk
write(output,' ',parcoord:1:0,' Gap"');
}
write(output,' ',inttopie(parcoord,apiece):1:0,' Gap"');
if flist^.m^.f^.orientation > 0 then begin
   write(output,  ' "-" "{}" "{}" 0.0');
   writeln(output,' ',abs(flist^.m^.coorddist):1:0);
end;
if flist^.m^.f^.orientation < 0 then begin
   write(output,  ' "-" "}{" "}{" 0.0');
   writeln(output,' ',-abs(flist^.m^.coorddist):1:0);
end;
*)

                  (* write the Gap definition out **********************)
(* TODO NOTE:  this code should be replaced with a call to writeadefinition !!! *)
                  write(scanfeatures,'define "');
                  writequotestringB(scanfeatures,flist^.m^.f^.nametag);
                  write(scanfeatures,'-(',round(abs(flist^.m^.coorddist)):1,')-');
                  writequotestringB(scanfeatures,
                                    mlist^.parentptr^.matrix^.riblname);
(* ZZZQQQ parcoord is the value to use for total: *)
{zzzqqqkkk
                  write(scanfeatures,' ',round(parcoord):1,' Gap"');
}
                  write(scanfeatures,' ',inttopie(round(parcoord),apiece):1,
                                       ' Gap"');

                  (* 2016 Aug 29: make clean standard gap bars: *)
                  write(scanfeatures,  ' "-" "||" "||"');
                  write(scanfeatures,' ',-abs(round(flist^.m^.coorddist)):1);
                  write(scanfeatures,' ',0:1);
                  writeln(scanfeatures);

{aaa}
(* write the Gap surprisal bar feature out ***********)
(*
write(output,'@ ');
writequotestringB(output,flist^.m^.f^.id);
if flist^.m^.parentevencoord <= flist^.m^.f^.evencoordinate
then write(output,' ',flist^.m^.parentevencoord:3:1);
if flist^.m^.parentevencoord > flist^.m^.f^.evencoordinate
then write(output,' ',flist^.m^.f^.evencoordinate:3:2); {aaa this one writes}
write(output,' ',flist^.m^.f^.orientation:1,' "');
writequotestringB(output,flist^.m^.f^.nametag);
write(output,'-(',abs(round(flist^.m^.coorddist)):1,')-');
writequotestringB(output,
                                    mlist^.parentptr^.matrix^.riblname);
write(output,' ',inttopie(parcoord,apiece):1,' Gap"');
if flist^.m^.gapsurp > 0
then write(output,' "',flist^.m^.gapsurp:1:1,' bits"')
else write(output,' "NA" ');
write(output, ' ', flist^.m^.gapsurp:1:infodecim);
write(output, ' ', mlist^.MinGapSurprisal:1:infodecim);
write(output, ' ', mlist^.MaxGapSurprisal:1:infodecim);
writeln(output);
*)

                  (* write the Gap surprisal bar feature out ***********)
(* TODO NOTE:  this code should be replaced with a call to writeafeature !!! *)
{write(scanfeatures,'ZZZQQQ-the@ is:');}
                  write(scanfeatures,'@ ');
                  writequotestringB(scanfeatures,flist^.m^.f^.id);

                  (* 2016 Aug 25: this is the original Gap logic for
                  which end of the gap to write out for its
                  coordinate:
                  if flist^.m^.parentevencoord <= flist^.m^.f^.evencoordinate
                  then write(scanfeatures,' ',    flist^.m^.parentevencoord:3:1);
                  if flist^.m^.parentevencoord >  flist^.m^.f^.evencoordinate
                  then write(scanfeatures,' ',    flist^.m^.f^.evencoordinate:3:1);
eg:
flist^.m^.parentevencoord=9.0
flist^.m^.f^.evencoordinate=4639662.0
                  ----
                  REVERSE THE LOGIC so that gaps are defined on the right side:
                  *)
{2016aug30 NO COORDINATE TESTING!
                  if flist^.m^.parentevencoord <= flist^.m^.f^.evencoordinate
                  then write(scanfeatures,' ',    flist^.m^.f^.evencoordinate:3:1)
                  else write(scanfeatures,' ',    flist^.m^.parentevencoord:3:1);
testing:
write(scanfeatures,' ',    9.0:3:1);
}
                  write(scanfeatures,' ',    flist^.m^.parentevencoord:3:1);
{GGGG}

{
writeln(output,
'flist^.m^.parentevencoord=',
 flist^.m^.parentevencoord:1:1);
writeln(output,
'flist^.m^.f^.evencoordinate=',
 flist^.m^.f^.evencoordinate:1:1);
}

(* 2007 Nov 20: TDS Ribls always have a zero coordinate and so always
have integer locations.  So gaps are always at integer locations.
The entire evencoordinate mechanism seems to be inappropriate.
I removed it by doing this:
write(scanfeatures,' ',flist^.m^.f^.coordinate:3:1);
Oddly, that made the gap bar run the wrong direction!
*)

                  write(scanfeatures,' ',flist^.m^.f^.orientation:1,' "');
                  writequotestringB(scanfeatures,flist^.m^.f^.nametag);
                  write(scanfeatures,'-(',abs(round(flist^.m^.coorddist)):1,')-');
                  writequotestringB(scanfeatures,
                                    mlist^.parentptr^.matrix^.riblname);
{zzzqqqkkk
                  write(scanfeatures,' ',parcoord:1:0,' Gap"');
}
                  write(scanfeatures,' ',inttopie(round(parcoord),apiece):1,
                                     ' Gap"');
                  if flist^.m^.gapsurp > 0
                  then write(scanfeatures,' "',flist^.m^.gapsurp:1:1,' bits"')
                  else write(scanfeatures,' "NA" ');
                  write(scanfeatures, ' ', flist^.m^.gapsurp:1:infodecim);
                  write(scanfeatures, ' ', mlist^.MinGapSurprisal:1:infodecim);
                  write(scanfeatures, ' ', mlist^.MaxGapSurprisal:1:infodecim);
         write(scanfeatures,' ',0.0:1:1); (* relative coordinate *)
{write(scanfeatures,'ZZZQQQBUBBA');}
                  writeln(scanfeatures);

{ old code:
                  if flist^.m^.gapsurp > 0
                  then write(scanfeatures,'',flist^.m^.gapsurp:1:1,' bits" ',
                                    flist^.m^.gapsurp:1:infodecim)
                  else write(scanfeatures,'NA" ',
                                    flist^.m^.gapsurp:1:infodecim);
                  writeln(scanfeatures,' 0.0000  0.0000');
}

{zzzddd}

               end;
            end;
            end

         end;
         flist^.m:= flist^.m^.next;
         mlist:= mlist^.next;
      end;
{writeln(output,'BUBBA 15---- multifeatwrite');}
{write(scanfeatures,'ZZZQQQ-LOOP_END-BUBBA');}

      (* 2013 Feb 07: tigher code here; report coordinate with total !!! *)
      write(data, 'totalinfo = ',flist^.totalinfo:infofield:infodecim,'   ');
{writeln(output,'BUBBA 15.0 ---- multifeatwrite');}
      flist^.m:= flist^.mstart;
{writeln(output,'BUBBA 15.1 ---- multifeatwrite');}
      while flist^.m <> nil do begin
{writeln(output,'BUBBA 15.2 ---- multifeatwrite');}
         if flist^.m^.writeitout then begin
            writequotestringB(data,flist^.m^.f^.nametag);
         end;
         if flist^.m^.next <> nil then begin
            if flist^.m^.next^.writeitout then write(data,'-');
         end;
         flist^.m:= flist^.m^.next;
      end;

(*
writeln(output,'BUBBA 15.3 ---- multifeatwrite WHAT???');
if flist  = nil then writeln(output,'BUBBA 15.35 ---- nil flist');
if flist <> nil then writeln(output,'BUBBA 15.35 ---- <> nil flist');
writeln(output,'BUBBA 15.4 ---- multifeatwrite');
if flist^.mstart = nil then writeln(output,'BUBBA 15.35 ---- nil flist^.mstart');
if flist^.mstart <> nil then writeln(output,'BUBBA 15.35 ---- <> nil flist^.mstart');
*)

      if flist^.mstart <> nil then begin

(*
if flist^.mstart^.f = nil then writeln(output,'BUBBA 15.35 ---- nil');
if flist^.mstart^.f <> nil then writeln(output,'BUBBA 15.35 ---- not nil flist^.mstart^.f');
if flist^.mstart^.f = nil then writeln(output,'BUBBA 15.35 ---- nil flist^.mstart^.f');
if flist^.mstart^.f <> nil then writeln(output,'BUBBA 15.35 ---- not nil flist^.mstart^.f');
*)

         write(data, ' ',round(flist^.mstart^.f^.coordinate):1);
{writeln(output,'BUBBA 15.9 ---- multifeatwrite');}
         writeln(data);
         writeln(data);
         writebar(mlist, mliststart, flist, apiece,givegroupnumber,groupnumber);
      end;

{writeln(output,'BUBBA 16---- multifeatwrite');}
{writeln(output,'BUBBA 17---- multifeatwrite');}
   end; (* 2016 oct 03 note: all the above code is not indented *)
   flist:= flist^.next;
   end;

{writeln(output,'BUBBA 21 ZZZQQQ here end of multifeatwrite');}
{writeln(scanfeatures,'ZZZQQQ here end of multifeatwrite');}
end;
(* end module multifeatwrite *)

(* begin module freeflist *)
procedure freeflist(var flist: multifeatureptr;
                    var fliststart: multifeatureptr;
                        mlist, mliststart: multitypeptr;
                    var totalsites: integer;
                        coord: integer;
                        boxsize: integer;
                    var freef: multifeatureptr;
                    var freefstart: multifeatureptr;
                         apiece: pieceptr;
                    var givegroupnumber: boolean;
                    var groupnumber: integer);

var
   inputlist: multifeatureptr; (* list to be written out *)
   doneloop: boolean; (* since the list is in order when you get to
                          the coordinate above the check kick out *)

begin

   flist:= fliststart;
   doneloop:= false;
   if flist <> nil then begin
      while (flist^.next <> nil) and not doneloop do begin
{zzzqqqlll}
{zzzqqqkkk another coordinate problem location:}
{
         if flist^.mstart^.f^.coordinate < (coord - 2 * boxsize) then begin
}
         if pietointR(flist^.mstart^.f^.coordinate, apiece)
            < (coord - 2 * boxsize) then begin

            fliststart:= flist^.next;
            flist^.next:= nil;
            inputlist:= flist;

            freef:= freefstart;
            while freef^.next <> nil do freef:= freef^.next;
            freef^.next:= flist;
            freef^.next^.next:= nil;

{writeln(output,'BUBBA multifeatwrite(mlist, mliststart, flist, inputlist, ...');}
            multifeatwrite(mlist, mliststart, flist, inputlist, totalsites,
                           apiece, givegroupnumber, groupnumber);

            flist:= fliststart;
            while flist <> nil do begin
               flist^.m:= flist^.mstart;
               flist:= flist^.next;
            end;
            flist:= fliststart;

            freef:= freefstart;
            freef:= freef^.next;
             while freef <> nil do begin
                freef^.m:= freef^.mstart;
                freef:= freef^.next;
             end;
             freef:= freefstart;

         end
         else begin
            doneloop:= true;
         end;


         flist:= fliststart;
      end;
   end;
   flist:= fliststart;
end;
(* end module freeflist *)

(* begin module multiscan.multiloop *)
procedure multiloop(var mlist, mliststart: multitypeptr;
                    var flist: multifeatureptr;
                    var fliststart: multifeatureptr;
                        ricutoff: real; singleA: char;
                    var freef: multifeatureptr;
                    var freefstart: multifeatureptr;
                    apiece: pieceptr);
(* Identify sites based on multiscan parameters *)

var
   sitefound: boolean; (* did you find a site? *)
   totalloop: boolean; (* is the entire loop finished *)
   dogo: boolean; (* A site has been found do go on to the next list *)

   {loopcount: integer; (* for debugging the loop *) (* BUBBA *)}
loopcount: integer; (* for debugging the loop *) (* nBUBBA *)
begin
{writeln(output,'BUBBA multiloop begin');}
{writeln(output,'nBUBBA multiloop begin TTT multiloop is the guilty party somehow');}
   mlist:= mliststart;
   mlist^.f:= mlist^.fstart;
   mlist^.f:= mlist^.f^.next;
   mlist^.fcurrent:= mlist^.f;
   dogo:= false;

   totalloop:= false;
   if mlist^.f = nil then totalloop:= true;
(*I'm guessing I put this in for a reason, but it is limiting my scan, if you
end up getting some sort of segmentation fault later, you may need to do
something where after you set mlist^.f:= mlist^.f^.next, you check if nil.
If nil, then doneloop, and kick out of something *)
   (* if not totalloop then if mlist^.f^.next = nil then totalloop:= true; *)

{loopcount := 1; (* BUBBA *)}
loopcount := 1; (* nBUBBA *)
   while not totalloop do begin
{writeln(output,'nBUBBA multiloop: totalloop loopcount = ', loopcount:1);}
{writeln(output,'BUBBA multiloop: totalloop loopcount = ', loopcount:1);}
{if loopcount = 2 then begin writeln(output,'BUBBA=',loopcount:1);halt;end;}

      if (mlist^.position = 'm') or (mlist^.position = 'e') then begin
         if mlist^.doneloop then mlist:= mlist^.prev;
         while mlist^.doneloop and (mlist^.position <> 't')
         do mlist:= mlist^.prev;
      end;
{writeln(output,'BUBBA multiloop: totalloop 2');}

      if mlist^.position = 'e' then begin
         mlist^.f:= mlist^.f^.next;
         mlist^.fcurrent:= mlist^.f;
         if mlist^.f = nil then mlist^.doneloop:= true;
      end;
{writeln(output,'BUBBA multiloop: totalloop 3');}

      if mlist^.position = 'm' then begin
        if mlist^.next^.doneloop then begin
if mlist^.f <> nil then begin
            mlist^.f:= mlist^.f^.next;
end;
            mlist^.fcurrent:= mlist^.f;
         end;
         if mlist^.f = nil then mlist^.doneloop:= true;

         mlist^.next^.f := mlist^.next^.fstart;
         mlist^.next^.doneloop:= false;
         if not (mlist^.doneloop) and (dogo) then mlist:= mlist^.next;
      end;
{writeln(output,'BUBBA multiloop: totalloop 4');}

      if mlist^.position = 't' then begin
         if mlist^.next^.doneloop then begin
            mlist^.f:= mlist^.f^.next;
            mlist^.fcurrent:= mlist^.f;
         end;
{writeln(output,'BUBBA multiloop: totalloop 5');}

         mlist^.next^.f := mlist^.next^.fstart;
         mlist^.next^.doneloop:= false;
         if mlist^.f = nil then mlist^.doneloop:= true;
         if not mlist^.doneloop then mlist:= mlist^.next;
      end;
{writeln(output,'BUBBA multiloop: totalloop 6');}

      sitefound:= false;
      if not mlist^.doneloop then begin
{writeln(output,'BUBBA multiloop: totalloop 6.1');}
         while not sitefound and not mlist^.doneloop do begin
{writeln(output,'BUBBA multiloop: totalloop 6.2');}
            paramcheck(mlist, mliststart, sitefound, dogo, flist, fliststart,
                       ricutoff, singleA, freef, freefstart,apiece);
{writeln(output,'BUBBA multiloop: totalloop 6.3');}
            if not sitefound then mlist^.f:= mlist^.f^.next;
{writeln(output,'BUBBA multiloop: totalloop 6.4');}
            if mlist^.f = nil then mlist^.doneloop:= true;
{writeln(output,'BUBBA multiloop: totalloop 6.5');}
         end;
      end;
{writeln(output,'BUBBA multiloop: totalloop 7');}

      if (mlist^.position = 'm') then begin
         if (mlist^.doneloop) and (mlist^.usetog = 'e') then
            mlist:= mlist^.prev;
         if (mlist^.doneloop) and (mlist^.usetog = 'n') then begin
            mlist^.next^.f := mlist^.next^.fstart;
            mlist^.next^.doneloop:= false;
            mlist:= mlist^.next;
         end;
      end;
{writeln(output,'BUBBA multiloop: totalloop 8');}

      if (mlist^.position = 'e') then begin
         if (not mlist^.foundsite) and (mlist^.usetog = 'n') then
         calcinfo(mlist, mliststart, flist, fliststart, ricutoff, singleA,
                  freef, freefstart, apiece);
         if (mlist^.doneloop) then mlist:= mlist^.prev;
      end;

      if (mlist^.position = 't') and (mlist^.doneloop) then totalloop:= true;

{writeln(output,'BUBBA multiloop: totalloop middle, loopcount =', loopcount:1);}
{if loopcount = 1 then halt (* BUBBA *)}
{loopcount := succ(loopcount); (* BUBBA *)}
loopcount := succ(loopcount); (* nBUBBA *)
   end;
{writeln(output,'BUBBA multiloop: totalloop END'); }

end;
(* end module multiscan.multiloop *)

(* begin module calcmaxtotal *)
procedure calcmaxtotal(var mlist: multitypeptr;
                       var mliststart: multitypeptr);
(* Calculate the maximum total information to use for range of colored
rectangles in lister. *)

var
   Ri: real; (* the ri of the consensus *)
   b: base; (* a base *)
   l: integer; (* the position in the matrix *)
   bestb: real; (* the information of the best base at that position *)
   total: real; (* the total highest Ri *)
   bestGap: real;  (* the gap surprisal of the best spacing in the histog *)
   worstGap: real; (* the gap surprisal of the worst spacing in the histog *)

begin

    mlist:= mliststart;
    b:= a;

    total:= 0;
    while mlist <> nil do begin
      Ri := 0;
      bestGap:= 1000;
      worstGap:= -1000;
      for l := mlist^.parameters.fromwanted
            to mlist^.parameters.towanted do begin
         bestb:= -10000;
         (* Ri  := Ri  + riblmatrix^.data[b, l]; *)
         b:= a;
         if mlist^.matrix^.data[b, l] > bestb then
            bestb:= mlist^.matrix^.data[b, l];
         b:= c;
         if mlist^.matrix^.data[b, l] > bestb then
            bestb:= mlist^.matrix^.data[b, l];
         b:= g;
         if mlist^.matrix^.data[b, l] > bestb then
            bestb:= mlist^.matrix^.data[b, l];
         b:= t;
         if mlist^.matrix^.data[b, l] > bestb then
            bestb:= mlist^.matrix^.data[b, l];
         Ri:= Ri + bestb;
      end;
      total:= total + Ri;

       if mlist^.gaptog = 'g' then begin
          for l:= mlist^.lowhistog to mlist^.highhistog do begin
             if mlist^.store[l].probdist < bestGap then
                bestGap:= mlist^.store[l].probdist;
             if mlist^.store[l].probdist > worstGap then
                worstGap:= mlist^.store[l].probdist;
          end;
       end
       else begin
          bestGap:= 0;
          worstGap:= 0;
       end;

       total:= total - bestGap;
       mlist^.MinGapSurprisal:= bestGap;
       mlist^.MaxGapSurprisal:= worstGap;
       mlist:= mlist^.next;
    end;

    mlist:= mliststart;
    mlist^.maxtotal:= total;

end;

(* end module calcmaxtotal *)

(* begin module PatentMessage *)
procedure PatentMessage(var f: text);
(* Write the patent message to file f
version = 3.02 of PatentMessage_module 2017 Jul 31
*)
procedure makebar;
begin
writeln(f,
'**********************************************************************');
end;
begin
makebar;
writeln(f,
'* By downloading this code you agree to the Source Code Use License: *');
writeln(f,
'* https://alum.mit.edu/www/toms/Source_Code_Use_License_Schneider.txt *');
writeln(f,
'* Contact: https://alum.mit.edu/www/toms/contacts.html                *');
makebar;
end;
(* end module PatentMessage version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module multiscan.themain *)
procedure themain(var book, ribl, scanp, data, scanfeatures, scaninst: text);
(* the main procedure *)
var
   count: integer; (* number of sequences done so far *)
   oneway: boolean; (* scan the sequences and not their complements *)
   parameters: allparameters; (* all of the user controlled parameters *)
   positionsevaluated: integer; (* positions evaluated in the piece *)
   sites: integer; (* number of sites found in the piece *)
   totalpositionsevaluated: integer; (* positions evaluated in the book *)
   totalsites: integer; (* number of sites found in the book *)
   previoussites: integer; (* previous number of sites found in the book *)
   riblmatrix: ribltypeptr; (* the weight matrix, Ri(b,l) *)

   (* book reading variables: *)
   org: orgkey;   orgchange, orgopen: boolean;
   chr: chrkey;   chrchange, chropen: boolean;
   apiece: pieceptr; piechange, pieopen: boolean;
   theline: integer; (* the line of the book *)

(* new multiscan variables for the main *)
   mlist: multitypeptr; (* the list of all the features for a site *)
   mliststart: multitypeptr; (* the first multifeature in the list *)
   holdmlist: multitypeptr; (* points to a certain mlist node *)
   numberofparts: integer; (* number of parts in the model, based on the
                              number of ribls *)
   boxsize: integer;  (* The size of the region that components can be found *)
   ricutoff: real;  (* total ri cutoff value *)

   singleA: char; (* a parameter that says to keep only the single best
      feature from the a features *)

   givegroupnumber: boolean; (* a parameter that says to add '#' and
      the group number to total lines at output on the line with totalinfo *)

   esscut: real; (* the number which the essential values have to be over, if
      negative it uses the scap cutoff and doesn't care about
      gap surprisal *)

   thelength: integer; (* the length of a piece *)
   newlength: integer; (* the new endpoint for the scan *)
   initialover: boolean; (* is the initial setup and scanning over *)

   freelist: featuretypeptr;  (* a list of not-used already made nodes for
                                 the mlist *)
   freeliststart: featuretypeptr;  (* The first object in the freelist *)

   (* coord: integer;  Original comment by Ryan:
      "the coordinate equivalent of the base in the scrolling loop"
      2016 Jul 26 I replaced this variable by internal coordinate calls so
      it is no longer needed.
   *)

   flist: multifeatureptr; (* a list of all the features *)
   fliststart: multifeatureptr; (* the first flist *)

   freef: multifeatureptr;  (* a list of not-used already made nodes for
                                 the flist *)
   freefstart: multifeatureptr;  (* The first object in the freef *)

   groupnumber: integer;  (* number of the group *)

procedure tellsymmetry(var f: text; oneway: boolean; ways: integer);
(* tell the directions of the scan to file f *)
begin
   if oneway
   then begin
       writeln(f,'* The sequences in the book will be scanned one way');
       if ways = 0
       then writeln(f,'* since the matrix is symmetric')
       else writeln(f,'* because you asked me to!');
   end
   else begin
       writeln(f,'* Scanning Both the sequences and their',
                 ' complements in the book');
       if ways = 0
       then writeln(f,'* since the matrix is asymmetric')
       else writeln(f,'* because you asked me to!');
   end;
   writeln(f,'*');
end;
begin (* themain *)
   writeln(output, 'multiscan ',version:4:2);

   PatentMessage(output);

   limitdate ('9','9','9','9','    /99/99 99:99:99');

   (* new multiscan reading routines to get parameters, ribls and histogs *)
   (* initialize the variables *)

   reset(ribl);
   reset(scanp);

   reset(multiscanp);
   reset(histog);

   readdiparameters(multiscanp, boxsize, ricutoff,
                    singleA, givegroupnumber,
                    esscut);

   new(mlist);
   new(mlist^.f);
{
   mlist^.f^.coordinate := -WEIRD; (* initializing originally -10000 *)
   mlist^.f^.coordinate := -10000; (* BUBBA TEST WEIRD *)
}
   mlist^.f^.coordinate := WEIRD; (* initializing originally -10000 *)
   mlist^.f^.next := nil;
   mlist^.fstart := mlist^.f;
   mlist^.prev := nil;
   mlist^.position := 't';
   mlist^.doneloop := false;
   mlist^.esscut := esscut;

   mliststart:= mlist;
   numberofparts:= 0;
   groupnumber := 0;

(* make the features list *)
{
   new(flist);
   new(flist^.m);
   new(flist^.m^.f);
   flist^.m^.f^.coordinate:= -1000;
   flist^.totalinfo:= 0;
   fliststart:= flist;
   flist^.next:= nil;
   flist^.m^.next:= nil;
   flist^.m^.f^.next:= nil;
}

   getriblmatrix(ribl,mlist^.matrix);
   readln(ribl);
   numberofparts:= numberofparts + 1;

   readparameters(scanp, mlist^.matrix, mlist^.parameters);
   readquotestring(multiscanp, mlist^.riblname);
   readparanstring(multiscanp, mlist^.parentname);

   if not equalstring(mlist^.riblname, mlist^.matrix^.riblname) then begin
      write(output, 'Parameter name: ');
      writequotestring(output, mlist^.riblname);
      write(output,' does not match ribl name: ');
      writequotestring(output,mlist^.matrix^.riblname);
      writeln(output,'.');
      halt;
   end;

   if not
   equalstring(mlist^.riblname, mlist^.parameters.sitedefinition^.nametag)
   then begin
      write(output, 'Parameter name: ');
      writequotestring(output, mlist^.riblname);
      write(output,' does not match scanp name: ');
      writequotestring(output,mlist^.parameters.sitedefinition^.nametag);
      writeln(output,'.');
      halt;
   end;

   skipblanks(multiscanp);
   read(multiscanp,mlist^.usetog);

   if mlist^.usetog <> 'e' then begin
      writeln(output,'The first model has to be essential.');
      halt;
   end;

   skipblanks(multiscanp);
   read(multiscanp,mlist^.orienttog);
   skipblanks(multiscanp);
   read(multiscanp,mlist^.gaptog);
   skipblanks(multiscanp);
   read(multiscanp, mlist^.bartog);
   skipblanks(multiscanp);
   read(multiscanp,mlist^.oltog);
   readln(multiscanp,mlist^.oldiff);

   mlist^.posint:= numberofparts;

   while not eof(ribl) do begin
      new(mlist^.next);
      new(mlist^.next^.f);
      mlist^.next^.esscut:= esscut;
{
      mlist^.next^.f^.coordinate:= -WEIRD; (* tds: initializing originally -10000 *)
}
      mlist^.next^.f^.coordinate:= WEIRD; (* tds: initializing originally -10000 *)
      mlist^.next^.doneloop:= false;
      mlist^.next^.fstart:= mlist^.next^.f;
      mlist^.next^.prev:= mlist;
      mlist:= mlist^.next;
      mlist^.next:= nil;

      getriblmatrix(ribl,mlist^.matrix);
      readln(ribl);
      numberofparts:= numberofparts + 1;

      readparameters(scanp, mlist^.matrix, mlist^.parameters);
      readquotestring(multiscanp, mlist^.riblname);
      readparanstring(multiscanp, mlist^.parentname);

      if not equalstring(mlist^.riblname, mlist^.matrix^.riblname) then begin
         write(output, 'Parameter name: ');
         writequotestring(output, mlist^.riblname);
         write(output,' does not match ribl name: ');
         writequotestring(output,mlist^.matrix^.riblname);
         writeln(output,'.');
         halt;
      end;

      if not
      equalstring(mlist^.riblname, mlist^.parameters.sitedefinition^.nametag)
      then begin
         write(output, 'Parameter name: ');
         writequotestring(output, mlist^.riblname);
         write(output,' does not match scanp name: ');
         writequotestring(output,mlist^.parameters.sitedefinition^.nametag);
         writeln(output,'.');
         halt;
      end;

      skipblanks(multiscanp);
      read(multiscanp,mlist^.usetog);
      skipblanks(multiscanp);
      read(multiscanp,mlist^.orienttog);
      skipblanks(multiscanp);
      read(multiscanp,mlist^.gaptog);
      skipblanks(multiscanp);
      read(multiscanp, mlist^.bartog);
      skipblanks(multiscanp);
      read(multiscanp,mlist^.oltog);
      readln(multiscanp,mlist^.oldiff);

      mlist^.posint:= numberofparts;

      if eof(histog) then begin
         writeln(output,'You do not have enough histogs',
                        ' for the number of components that you have.');
         writeln(output,'You need n-1 histogs.');
         halt;

      end;
      readln(histog);
      distribution(histog, mlist^.distsize, mlist^.correction, mlist^.store,
                   mlist^.lowhistog, mlist^.highhistog);

   end;

   (* end of multiscan reading routines to get parameters, ribls and histogs *)

   (* Set the parent link in each of the models to the right position *)
   mlist:= mliststart;
   if numberofparts > 1 then begin
      mlist:= mlist^.next;
      while mlist <> nil do begin
         holdmlist:= mliststart;
         while (holdmlist^.next <> nil) and
         not equalstring(mlist^.parentname, holdmlist^.riblname) do begin
            holdmlist:= holdmlist^.next;
         end;
         if holdmlist^.next = nil then begin
            write(output,'The parent name you gave: ');
            writequotestring(output,mlist^.parentname);
            writeln(output,' does not match any ribls.');
            halt;
         end;
         mlist^.parentptr:= holdmlist;
         mlist^.parent:= holdmlist^.posint;
         if mlist^.parent = mlist^.posint then begin
            writeln(output,'A model cannot point to itself!!!');
            halt;
         end;
         if holdmlist^.usetog <> 'e' then begin
   writeln(output,'For something to be a parent it has',
                  ' to be essential.... for now.');
            halt;
         end;

         mlist:= mlist^.next;
      end;
   end;

   mlist:= mliststart;

   parameters:= mlist^.parameters;
   riblmatrix:= mlist^.matrix;

   (* always write out data *)
   rewrite(data);
   (* by making the first line as follows, alist can use the file directly *)
   writeln(data,'* bits');
   writeln(data, '* multiscan ',version:4:2);
   writeln(data, '* with matrix from:');
   writeln(data, '***********************************************************');
   reset(ribl);
   while ribl^='*' do copyaline(ribl,data); (* copy the header *)
   writeln(data, '***********************************************************');
   writeln(data,'* scan matrix is FROM = ',
                riblmatrix^.frombase:1,', TO = ',riblmatrix^.tobase:1);
   with parameters do writeln(data,'* region used for the scan is FROM = ',
                fromwanted:1,', TO = ',towanted:1);
   writeparameters(data, parameters);

   writeln(data, '* Scan of book:');
   reset(book);
   copyaline(book,data);

   writeln(data,'*  mean = ',riblmatrix^.mean:infofield:infodecim);
   writeln(data,'* stdev = ',riblmatrix^.stdev:infofield:infodecim);
   writeln(output,'*  mean = ',riblmatrix^.mean:infofield:infodecim);
   writeln(output,'* stdev = ',riblmatrix^.stdev:infofield:infodecim);

   if parameters.outinst then begin
      rewrite(scaninst);
      writeln(scaninst, 'title "multiscan ',version:4:2,'";');
      writeln(scaninst, '(*');
      writeln(scaninst, '* Scan of book:');
      reset(book);
      copyaline(book,scaninst);
      writeparameters(scaninst, parameters);
      writeln(scaninst, '*)');
      writeln(scaninst);
      writeln(scaninst,'default out-of-range reduce-range;');
      writeln(scaninst,'default numbering piece;');
      writeln(scaninst);
   end;

   mlist:= mliststart;
   rewrite(scanfeatures);
   while mlist <> nil do begin
      if mlist^.parameters.outfeat then begin
         with mlist^.parameters do if mlist^.parameters.sitedefined then begin
            if mlist = mliststart then rewrite(scanfeatures);
            writeln(scanfeatures, '* multiscan ',version:4:2);
            reset(book);
            copyaline(book, scanfeatures);
            writeln(scanfeatures); (* blank line *)
            writeadefinition(scanfeatures,mlist^.parameters.sitedefinition);
            writeln(scanfeatures); (* blank line *)
   {
   zzzmmm this is done in writeadefinition! NO
   but sitedefinition needs to link to riblmatrix!
   }
            putriblmatrix(scanfeatures,mlist^.matrix);
         end;
      end;
      mlist:= mlist^.next;
   end;


   (* set position char in mlist  *)
   mlist:= mliststart;
   mlist:= mlist^.next;
   while mlist <> nil do begin
      if (mlist^.next = nil) then mlist^.position:= 'e';
      if (mlist^.next <> nil) then mlist^.position:= 'm';
      mlist:= mlist^.next;
   end;

   calcmaxtotal(mlist, mliststart);

   mlist:= mliststart; (* reset the mlist *)

   brinit(book,theline);

   count := 0;

   with parameters do begin

      case ways of
         0: oneway := riblmatrix^.symmetry <> 'a';
         1: oneway := true;
         2: oneway := false;
      end;

      if outdata then begin
         writeln(data,'*');

         writeln(data,'* lowest   Ri cutoff = ', cutDNRi:infofield:infodecim);
         writeln(data,'* highest  Ri cutoff = ', cutUPRi:infofield:infodecim);
         writeln(data,'* lowest    Z cutoff = ', cutDNZ :infofield:infodecim);
         writeln(data,'* highest   Z cutoff = ', cutUPZ :infofield:infodecim);
         writeln(data,'* lowest  probability = ', cutDNP :infofield:infodecim);
         writeln(data,'* highest probability = ', cutUPP :infofield:infodecim);

         tellsymmetry(data,oneway,ways);

         writeln(data,'* DEFINITION OF THE DATA COLUMNS:');
         writeln(data,'* 1 model name');
         writeln(data,'* 2 model coordinate');
         writeln(data,'* 3 matrix orientation',
                      ' (+1 = as in book, -1 = complement)');
         writeln(data,'* 4 Ri evaluation (bits per site)');
         writeln(data,'* 5 gap size');
         writeln(data,'* 6 gap surprisal');

         with parameters do

         writeln(data,'*');
      end;

      tellsymmetry(output,oneway,ways);

      totalsites := 0;
      previoussites := 0;
      totalpositionsevaluated := 0;
      writeln(output,'* count of the sequences scanned (every ',
                      countmark:1,'):');

      if riblmatrix^.stdev <= 0 then begin
         writeln(output,'WARNING: The standard deviation is zero.');
      end;
      orgchange := false; orgopen := false;
      chrchange := false; chropen := false;
      piechange := false; pieopen := false;
      apiece := nil;
      theline := 0;

      repeat
         getocp(book, theline,
                org,    orgchange, orgopen,
                chr,    chrchange, chropen,
                apiece, piechange, pieopen);
{writeln(output,'BUBBA 1');}

         if pieopen then begin

{zzzqqq remove when properly debugged - this should not matter!}
(* 2016 Jun 30: removed
         if inttopie(1,apiece) < 1 then begin
writeln(output,'Multiscan cannot handle negative coordinates.');
writeln(output,'Change your inst file so that it does not',
               ' fall below coordinate 1.');
            halt;
         end;
*)

            (* make the features list *)
            new(flist);
            (* 2016 oct 05: CALLING clearmultifeatureptr
               PREVENTED THE 12.25 crash! *)
            clearmultifeatureptr(flist); {FFF}
            new(flist^.m);
            new(flist^.m^.f);
            flist^.m^.f^.coordinate:= WEIRD; (* tds: init. originally -1000 *)

            (* None of clear variables below prevented the 12.25 crash
               but they are good to do! *)
            clearstring(flist^.m^.f^.id);
            clearstring(flist^.m^.f^.nametag);
            clearstring(flist^.m^.f^.othertag);
            flist^.m^.f^.others := nil;
            flist^.m^.f^.mothers := nil;
            flist^.m^.f^.duplicateprevious := nil;
            flist^.m^.f^.next := nil;
            flist^.totalinfo:= 0;

{writeln(output, 'BUBBA FFF 999 set fliststart:= flist; *****');}
            fliststart:= flist;
            flist^.next:= nil;
            flist^.m^.next:= nil;
            flist^.m^.f^.next:= nil;

            count := count + 1;
            if (count mod countmark) = 0 then writeln(output,count:1);

            if outinst then begin
               (* write organism and chromosome instructions *)
               if orgchange then begin
                  write(scaninst,'organism ');
                  writename(scaninst,org.hea.keynam);
                  write(scaninst,';');
                  (* if orgchange is true, chrchange MUST be true! *)
                  chrchange := true;
               end;
               if chrchange then begin
                  if orgchange then write(scaninst,' ');
                  write(scaninst,'chromosome ');
                  writename(scaninst,chr.hea.keynam);
                  write(scaninst,';');
               end;
               if chrchange or orgchange
               then writeln(scaninst);
               (* reset so these instructions are only written as needed: *)
               orgchange := false;
               chrchange := false;
            end;

            (* Below is the core of the new code, which makes this
               multiscan, instead of just scan. *)

{This seems to be in the wrong place?  it is done for every piece.}
            (* Initialize the freelist for featuretypeptr types, this
               is the list of nodes which were kicked off of mlist^.f
               lists *)
            new(freelist);
            freeliststart:= freelist;
            new(freelist^.next);
            freelist:= freelist^.next;
            freelist^.next:= nil;

            (* Initialize the freelist for multifeatureptr types, this
               is the list of nodes which were kicked off of mlist^.f
               lists *)
            new(freef);
            clearmultifeatureptr(freef);
            freefstart:= freef;
{
            freef^.next:= nil;
this is defined above - could remove from here.
procedure clearmultifeatureptr(var mfp: multifeatureptr);
(* clear a multifeatureptr *)
begin
   freef^.m := nil;
   freef^.mstart := nil;
   freef^.totalinfo := 0;
   freef^.writeitout := false;
   freef^.next:= nil;
end;
}

            (* First, we need to scan the first range of bases
               up to boxsize or the piece length (the smaller of these) *)
            thelength := piecelength(apiece);
            newlength:= boxsize;
{writeln(output,'nBUBBA boxsize = ', boxsize:1,'=================================');}
{writeln(output,'nBUBBA before if newlength = ', newlength:1,' =========');}
            if newlength > thelength then newlength:= thelength;
{writeln(output,'nBUBBA after  if newlength = ', newlength:1,' =========');}
            initialover:= false;

            while mlist <> nil do begin
{writeln(output,'BUBBA 2');}
{writeln(output,'nBUBBA => while mlist <> nil do begin');}
               mlist^.f:= mlist^.fstart;
               mlist^.f^.next:= nil;

{TTT Tom working point}
               scansequence({data, scanfeatures,} scaninst,
                            {mlist^.parameters.outdata,}
                            mlist^.parameters.outinst,
                            mlist^.matrix,
                            mlist^.parameters.fromwanted,
                            mlist^.parameters.towanted,
                            apiece,
                              mlist^.parameters.cutDNRi,
                              mlist^.parameters.cutUPRi,
                              mlist^.parameters.cutDNZ,
                              mlist^.parameters.cutUPZ,
                              mlist^.parameters.cutDNP,
                              mlist^.parameters.cutUPP,
                            sites, positionsevaluated,
                            not oneway, mlist^.matrix^.mean,
                            mlist^.matrix^.stdev,
                            mlist^.parameters.sitedefined,
                            {mlist^.parameters.sitedefinition,}
                            mlist^.parameters, newlength, mlist^.f, initialover,
                            freelist, freeliststart);
                            (* totalsites:= totalsites + sites; *)
{writeln(output,'nBUBBA sites = ',sites:1);}

               mlist^.fend:= mlist^.f;
               mlist:= mlist^.next;
            end;
{writeln(output,'BUBBA 3');}

(* We do the first analysis of the sites in the lists, make the features out
   list, and prepare to move on *)

            mlist:= mliststart;
            while mlist <> nil do begin
               mlist^.doneloop:= false;
               mlist:= mlist^.next;
            end;
{writeln(output,'BUBBA 4');}
(*
writeln(output,'nBUBBA before multiloop THAT FAILS!'); {TTT}
writeafeature(output,fliststart^.m^.f);writeln(output,'<-fliststart^.m^.f nBUBBA');
*)
            multiloop(mlist, mliststart, flist, fliststart, ricutoff, singleA,
                      freef, freefstart, apiece);
{writeln(output,'BUBBA 5');}
(*
writeafeature(output,fliststart^.m^.f);writeln(output,'<-fliststart^.m^.f nBUBBA');
writeln(output,'nBUBBA after multiloop');
*)

            (* Now that we have created the lists and have the first
               set of coords, we move only 1 base at a time *)
            initialover:= true;
            newlength:= newlength + 1;
(*
writeln(output,'nBUBBA 1 newlength = ',newlength:1);
writeln(output,'nBUBBA 2 thelength = ',thelength:1);
{the book piece lengths are 71 but:
BUBBA 1 newlength = 72
BUBBA 2 thelength = 71
so it never enters the loop!
newlength:= newlength - 1; does not help
}
writeln(output,'nBUBBA ++++++++++++++++++++++++++++++++++++ BEFORE that');
*)
(*
writeln(output,'BUBBA 5.0');
writeln(output,'BUBBA 5.01 newlength=', newlength:1);
writeln(output,'BUBBA 5.01 thelength=', thelength:1);
*)
            while newlength <= thelength do begin
{writeln(output,'BUBBA 5.1'); halt;}
{writeln(output,'nBUBBA 3 inside forbidden city - never gets here!');}
               mlist:= mliststart;
               while mlist <> nil do begin

(* right here is where you do the check to see if you need to kick it off:
   Inside the writeitout procedure, you need to put an if statement that
   says "if not initial don't make a new unless there are no extras on the
   free list" *)

(* ALSO: You could have duplicates at the very end of the file, because the
ends cycle if too long.  They will all be next to each other so if
f^.coordinate = f^.next^.coordinate then put the node on the free list, but
don't forget that based on orientation, your features could be shuffled by
one *)

{
                  RYANcoord:= inttopie(newlength, apiece) - boxsize;
}
{zzzqqq wrong logic! never compare piece coords!}
{where else is coord used? nowhere!}

                  mlist^.f:= mlist^.fstart;
                  mlist^.f:= mlist^.f^.next;
                  if mlist^.f <> nil then begin
(* 2016 Jun 24: This was Ryan's original logic:
                     if mlist^.f^.coordinate < RYANcoord then begin
zzzqqqlll wrong logic! never compare piece coords!
The problem is that the RYANcoord (originally named 'coord') is
compared using external coordinates.  ALL comparisons need to be
done using internal 1 to n coordinates!!!
*)
{
write(output,'BUBBA, boxsize: ', boxsize:1);
write(output,' mlist^.f^.coordinate: ', mlist^.f^.coordinate:5:1); (* BUBBA *)
write(output,' newlength: ', newlength:1); (* BUBBA *)
writeln(output); (* BUBBA *)
}
                     if pietointR(mlist^.f^.coordinate, apiece)
{
                        < newlength then begin
This works but is missing boxsize.
But adding boxsize is still wrong, so there must be more
locations where there are problems: zzzqqqkkk!}
                        < newlength - boxsize then begin
{writeln(output,'nBUBBA ++++++++++++++++++++++++++++++++++++ coo tst=never here');}
                        freelist^.next:= mlist^.f;
                        mlist^.f:= mlist^.fstart;
                        mlist^.f^.next:= mlist^.f^.next^.next;
                        freelist:= freelist^.next;
                        freelist^.next:= nil;
                     end;
                  end;

(* THIS IS SO YOU DON'T ADD YOUR NEW SITE TO THE END OF THE FREELIST IF
IT IS THE LAST SITE ON THE NODE !!!!!  I HATE THIS LINE OF CODE *)
                  if mlist^.fstart^.next = nil then mlist^.fend:= mlist^.fstart;

                  mlist^.f:= mlist^.fend;

{writeln(output,'nBUBBA ++++++++++++++++++++++++++++++++++++ scansequence');}
                  scansequence({data, scanfeatures,} scaninst,
                               {mlist^.parameters.outdata,}
                               mlist^.parameters.outinst,
                               mlist^.matrix,
                               mlist^.parameters.fromwanted,
                               mlist^.parameters.towanted,
                               apiece,
                                 mlist^.parameters.cutDNRi,
                                 mlist^.parameters.cutUPRi,
                                 mlist^.parameters.cutDNZ,
                                 mlist^.parameters.cutUPZ,
                                 mlist^.parameters.cutDNP,
                                 mlist^.parameters.cutUPP,
                               sites, positionsevaluated,
                               not oneway, mlist^.matrix^.mean,
                               mlist^.matrix^.stdev,
                               mlist^.parameters.sitedefined,
                               {mlist^.parameters.sitedefinition,}
                               mlist^.parameters, newlength,
                               mlist^.f, initialover,
                               freelist, freeliststart);
                  (* totalsites:= totalsites + sites; *)

                  mlist^.fend:= mlist^.f;
                  mlist:= mlist^.next;
               end;

               mlist:= mliststart;
               while mlist <> nil do begin
                  mlist^.doneloop:= false;
                  mlist:= mlist^.next;
               end;

{
               freeflist(flist, fliststart, mlist, mliststart, totalsites,
                 inttopie(newlength, apiece), boxsize, freef, freefstart, apiece);
}
               (* 2016 Jun 30: newlength should be internal coordinates
                  so don't convert to external ! *)
{TTT}
               freeflist(flist, fliststart, mlist, mliststart, totalsites,
                 newlength, boxsize, freef, freefstart, apiece,
                 givegroupnumber, groupnumber);

               multiloop(mlist, mliststart, flist, fliststart, ricutoff,
                         singleA, freef, freefstart, apiece);

               newlength:= newlength +1;
if (newlength mod 50000) = 0 then writeln(output,newlength:1);
{writeln(output,'nBUBBA ++++++++++++++++++++++++++++++++++++ never inside that');}
            end;
{writeln(output,'BUBBA 6');}
{writeln(output,'nBUBBA ++++++++++++++++++++++++++++++++++++ PAST that');}

            totalpositionsevaluated := totalpositionsevaluated
                                     + positionsevaluated;

(* Final feature sort and write the features out *)
            flist:= fliststart;
{zzzqqqkkk
            if (flist^.m^.f^.coordinate > 0) then begin
}

(*
{2016sep20: The problem is before this point since when it fails it gives WEIRD.
However the test is bad because it disallows negative legitimate coordinates.
Must devise a new test!  TTT}
write(output,'nBUBBA maxint = ',maxint:1,' ');
writeln(output,'nBUBBA WEIRD = ',WEIRD:1);
writeln(output,'nBUBBA flist^.m^.f^.coordinate = ',
                       flist^.m^.f^.coordinate:1:1); {nBUBBA}
if flist^.m^.f^.coordinate = WEIRD then writeln(output,'nBUBBA =WEIRD');
writeln(output,'nBUBBA pietointR(flist^.m^.f^.coordinate,apiece) = ',
                       pietointR(flist^.m^.f^.coordinate,apiece):1); {nBUBBA}
{
writeafeature(output,flist^.m^.f); writeln(output,'<- flist^.m^.f nBUBBA');
}
writeafeature(output,fliststart^.m^.f); writeln(output,'<- fliststart^.m^.f nBUBBA');
*)

{writeln(output,'BUBBA 7');}
            if pietointR(flist^.m^.f^.coordinate,apiece) > 0 then begin
{writeln(output,'nBUBBA ------------------------- multifeatwrite call');}
{writeln(output,'BUBBA 8 multifeatwrite(mlist, mliststart, flist, fliststart, ...');}
{writeln(output,'BUBBA 8.5 fliststart is bad ... ************');}
               multifeatwrite(mlist, mliststart, flist, fliststart, totalsites,
                              apiece, givegroupnumber, groupnumber);
{writeln(output,'BUBBA 8.5 ************** number = ', number:1); }
            end;
{writeln(output,'BUBBA 9');}

            if outdata then begin
{writeln(output,'nBUBBA ------------------------- doing outdata');}
{
               writeln(data,'* ',totalsites:1,' sites found in piece',
zq
}
               (* 2013 Feb 07 revised: give sites on this piece and coordinates *)
               write(data,'* ',totalsites-previoussites:1,
                          ' sites found in piece',
                          ' ',apiece^.key.hea.keynam.letters:
                              apiece^.key.hea.keynam.length,
                          ' from ',apiece^.key.piebeg:1,
                          ' to ',apiece^.key.pieend:1,
                          ' dir ');
               (* 2017 Jun 01: add direction data: *)
               if    (apiece^.key.piedir = dirhomologous)
                  or (apiece^.key.piedir = plus) 
               then write(data,'+')
               else write(data,'-');
               writeln(data, ' ', positionsevaluated:1,' bp evaluated');
{
zzzbbb
               writeln(data,'* ',positionsevaluated:1,' positions evaluated');
}
               previoussites := totalsites;
            end;
{writeln(output,'nBUBBA HEEERRRRRRRRREEEEEEEEEEEE clearpiece(apiece);');}
            clearpiece(apiece); (* clear the piece for reuse *)

{
writeln(data,'here');
freeflist(flist, fliststart, mlist, mliststart, totalsites,
inttopie(newlength, apiece), boxsize, freef, freefstart, apiece);
writeln(data,'here');
multifeatwrite(mlist, mliststart, flist, fliststart, totalsites, apiece);
}

            flist:= fliststart;
            dispose(flist);

         end;
      until (not pieopen) or (count >= todo);

      if outdata then begin
         writeln(data,'*');
         writeln(data,'* ',totalsites:1,' sites found in this book in');
         writeln(data,'* ',totalpositionsevaluated:1,
                         ' total positions evaluated');
         if totalsites > 0
         then writeln(data,'* ',
              ln(totalpositionsevaluated/totalsites)/ln(2):infodecim:infofield,
              ' effective Rfrequency (bits per site)');
      end;
      writeln(output,'*');
      writeln(output,'* ',totalsites:1,' sites found in this book');

   end;
end;
(* end module multiscan.themain *)
begin
   themain(book, ribl, scanp, data, scanfeatures, scaninst);
1: end.
