program scan(book, ribl, scanp, data, scanfeatures, scaninst, output);
(* scan: scan a book with an Ribl weight matrix and generate a vector

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

module libraries: delman, delmods, matmod
*)

label 1;  (* end of program *)

const
(* begin module version *)
   version = 3.67; (* of scan.p 2017 Aug 07
2017 Aug 07, 3.66: stop showfeature for every feature
2017 Jul 14, 3.65: namelength/linelength set to 200/210
2017 Jul 12, 3.64: namelength set to 200 to allow even longer names
2017 Mar 15, 3.63: modules upgraded from prgmod: no censor
2017 Feb 14, 3.62: modules upgraded from prgmod
2017 Feb 14, 3.61: backup of modules upgraded from prgmod
2017 Feb 13, 3.60: backup
2012 Oct 12, 3.59: update patent statement
2011 Dec 10, 3.58: identify problem with ardht.ribl
2011 Jun 15, 3.57: clean up
2011 Jun 15, 3.56: scanp name must match ribl name!
2007 Mar 26, 3.55: cleanup
2007 Mar 26, 3.54: bug: getocp requires piece to start with nil!
2007 Mar 21, 3.53: upgrade from matmod
2007 Mar 21, 3.52: clean up
2007 Mar 21, 3.51: try more sensible tolerance in simpson - that worked!
2007 Mar 21, 3.50: debug 64 bit compiler ...
2007 Mar 21, 3.49: debug 64 bit compiler ...
2007 Mar 21, 3.48: debug 64 bit compiler hangs at ... simpson integration!
                   POINT 5.5.3.1.6.2.5 simpson inner loop i=52055
2007 Mar 21, 3.47: debug 64 bit compiler hangs at 5.5.3.1
2007 Mar 21, 3.46: debug 64 bit compiler hangs at 5.5.3
2007 Mar 21, 3.45: debug 64 bit compiler hangs at 5.5
2007 Mar 21, 3.44: debug 64 bit compiler hangs at 5
2007 Mar 13, 3.43: stringforFullName from stringforName
2005 May 30, 3.42: numberrealstring fixed - no leading spaces unless < thewidth
2005 May 30, 3.41: numberrealstring fixed - no leading spaces
2005 May 28, 3.40: fancy digitcounter
2005 May 28, 3.39: cleanup, this time for real!!
2005 May 26, 3.38: sigh. invention of emptystring for numberrealstring
2005 May 26, 3.37: more cleanup
2005 May 26, 3.36: cleanup
2005 May 26, 3.35: loss of 'otherquotestring' ... (ie no 'bits')
2005 May 26, 3.34: use clearfeature instead of version 1.33's method
2005 May 26, 3.33: segmentation fault because of new feature parts.
                   The solution is to set:
      f^.others := nil;
      f^.mothers := nil;
2005 May 23, 3.32: add patent message
2004 Dec 17, 3.31: scanp example must have version number
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 April 26: fix bug so that scansequence can handle circular sequences!
1999 April 24: failed attempt to remove leading blanks on numbers
1999 Jan 18: add references: Schneider.Ri and Schneider.walker
1998 August 1: fix getocp to make scaninst results simpler
orgchange:=false; chrchange:=false; has to be called at the start.
1998 July 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 July 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 January 26: namelength set to 100 to allow long names
1997 April 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 April 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 March 28: deltaRi implemented
1997 March 28: parameter version control implemented (2.69)
1997 Jan 4: scan now works with single arrays of 10 million bases
1996 November 18: data file has "* bits" on top so alist can use it (2.68)
1996 June 15: moved matrix symmetry to ri.
1996 March 22: added scaninst and parameters to turn on and off output.
1996 March 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 *)
updateversion = 2.70; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.scan *)
(*
name
   scan: scan a book with an Ribl weight matrix and generate a vector

synopsis
   scan(book: in, ribl: in, scanp: in,
        data: out, scanfeatures: out, scaninst: out, output: out)

files

   book: a book from the delila system

   ribl: a weight matrix from sites or ri programs.
      Lines that start with * are notes.  the next line contains the matrix
      FROM-TO coordinates, this is followed by the matrix in the order A, C, G,
      T from FROM to TO.

   scanp: parameters to control the program.

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

      seqs: One integer on the first line is the number of sequences to scan
         to produce the vector.  0 = none, positive = that number; negative =
         all.

      Ri range : Two real numbers on the second line give the range of
         information content to report in the data file.

      Z score range: Two real numbers on the third line give the range of the
         Z score 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.  Note: normally
         one would want the lower number to be zero.

      Probability range: Two real numbers on the fourth line give the range
         of probability to report in the data file.  The probability of a
         site is determined from the mean and standard deviation of the Ri
         distribution.  Note: normally one would want the lower number to be
         zero.

      fromwanted towanted range: two integers that define the FROM-TO range
         of the ribl matrix to use for computations.  This is independent of
         the range displayed in the walker.

      ways:  One integer.  2 means scan both the sequence and its
         complement.  1 means simply scan the sequence.  0 means to let the
         program figure it out.  The Ri program determines the symmetry of
         the matrix.  If it is symmetrical, it will only scan one way.  If it
         is asymmetrical, both scans are done.

      sitedefinition:  If the first non-blank character on the line is 'd',
         then the rest of the line contains a definition of how to write out
         the sites.  If no site is defined, the scanfeatures file will not be
         written to.  See program lister.p for details.  The basic format for
         an ASCII definition looks like this:

         define "Fis" "-" "[0]" "[0]" -7  0 +7

         For a walker it looks like:

         define "Fis" " w" "  " "  " -7 +7

         NOTE: the range for walker display (given in this site definition)
         is independent of the range of the weight matrix used for
         computation (given in the fromwanted and towanted parameters).

      print definitions:  Any number of lines that define how to print the
         "other" feature string in each feature definition.  The data that may
         be printed are the same as those in the data file.  They are:

         #           width
         length      width
         name        width
         coordinate  width
         orientation width
         Ri          width decimal
         Z           width decimal
         probability width decimal
         string      "quote string"
         .           end of print definitions

         If the first character on a line is '#', the line defines the
         width for the coordinate of the number of the DNA piece from the book.

         If the first character on a line is 'l' or 'L', the line defines the
         width for the length of the DNA piece in the book.

         If the first character on a line is 'n' or 'N', the line defines the
         width for the name of the DNA piece in the book.

         If the first character on a line is 'f' or 'F', the line defines the
         width for the fullname of the DNA piece in the book.

         If the first character on a line is 'c' or 'C', the line defines the
         width for the coordinate of the zero base of the site.

         If the first character on a line is 'o' or 'O', the line defines the
         width for the orientation of the site.  If the width is 1, the
         orientation is given as + or -, if ithe width is larger the
         orientation is given as -1 or +1.

         If the first character on a line is 'r' or 'R', the line defines the
         width and decimal fields for the individual information in bits.  The
         word "bits" is attached to the end of the string.

         If the first character on a line is 'p' or 'P', the line defines the
         width and decimal fields for the probability of the site.

         If the first character on a line is 'z' or 'Z', the line defines the
         width and decimal fields for the Z score of the site.

         If the first character is 's' or 'S' then the line defines a string to
         insert.

         The end of the file or a period "." ends the print definitions.

         The lines may be put in any order and this defines the order that they
         will be printed to the "other" string.  If the first character is not
         found (as, for example by having a blank in front of it), the
         corresponding data will not be printed.  This gives the user full
         control of the "other" string contents.

         The only kind of definition that may be repeated is the "string".
         This allows the user to put whatever they desire between the data
         items.

      file output definitions:  The first three characters on the line define
         which files will be output.  Capital characters turn on the output.
         Small characters turn it off.  The files are data, (scan)features,
         and (scan)inst so the characters are d, f and i, respectively.  Thus
         DfI turns on the data and scaninst files and leaves the scanfeatures
         off.  (Unidentified characters default to upper case.)

      normalizeRi:  The first character is defines how to normalize
         the reported Ri values.  The Ri value at coordinate zero
         is called Ri0.

         n: normal: scan and report Ri

         s: subtract: compute Ri(l) - Ri(0) at each position l

         d: divide: compute Ri(l) / Ri(0) at each position l

         The s and d modes are usually to be used in conjunction with
         renumbering by Delila (the 'default coordinate zero' command).

      instfrom, instto: range of Delila instructions produced in scaninst
         if that file is created.

   data: The results.  Comments are lines that begin with '*'.  The columns are
      defined in comments in the file.  The matrix is searched over both the
      sequence and its complement.  Ri is reported, as is the Z and probability
      based on the mean and st.dev.

   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:

      @ K01789 229.0 -1 "dnaA" "+12.2 bits " 12.200338    -0.473212     0.318031

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

   output: messages to the user

description
   The Ri(b,l) weight matrix is scanned across the sequences in the book to
   produce a vector.

examples

   Example scanp files:

3.00    version of scan that this parameter file is designed for.
-1      number of seqs to scan 0 = none, positive = that number; negative = all
0       information content at or above which to report in the data file.
100     Z score below which to report in the data file
0       probability at or above which to report in the data file.
-10 +10 desired region of the ribl weight matrix to use
0       0: program figures it out; 1: one way scan; 2: two way scan.
define "Fis" "-" "[0]" "[0]" -7  0 +7
string "data at:" string: A string listed at the feature
coordinate 5
string " Ri = "   string: A string listed at the feature
Ri 5 1  Riwidth Ridecimal: character places for reporting bits to scanfeatures
string " Z = "    string: A string listed at the feature
Z 4 1   z score
string " p = "    string: A string listed at the feature
probability 5 2
.       end of print definitions
DFI     dfi: data, features, inst: files output
n       normalizeRi: n: normal, s: Ri(l)-Ri(0), d: Ri(l)/Ri(0)
-50 +50 instfrom, instto: range to make the scaninst file (if made)
   scanp: parameters to control the program.

3.00    version of scan that this parameter file is designed for.
-1      number of seqs to scan 0 = none, positive = that number; negative = all
0   100 information content at or above which to report in the data file.
0   100 Z score below which to report in the data file
0   1   probability at or above which to report in the data file.
-10 +10 desired region of the ribl weight matrix to use
1       0: program figures it out; 1: one way scan; 2: two way scan.
define "Fis" " w" "  " "  " -10 10
 string "@"
 coordinate 4
 string "|Ri="
Ri 4 1  Riwidth Ridecimal: character places for reporting bits to scanfeatures
string " bits"
 string "|Z="
 Z 7 4   z score
 string "|p="
 probability 5 3
.  end of print definitions
dFi     dfi: data, features, inst: files output
n       normalizeRi: n: normal, s: Ri(l)-Ri(0), d: Ri(l)/Ri(0)   
-50 +50 range to make the scaninst file (if made)

documentation

@article{Schneider.Ri,
author = "T. D. Schneider",
title = "Information Content of Individual Genetic Sequences",
journal = "J. Theor. Biol.",
volume = "189",
number = "4",
pages = "427-441",
note = "http://www-lecb.ncifcrf.gov/$\sim$toms/paper/ri/",
comment = "indiv.tex",
comment = "Submitted, April 1997",
year = "1997"}

@article{Schneider.walker,
author = "T. D. Schneider",
title = "Sequence Walkers:
a graphical method to display how binding proteins
interact with {DNA} or {RNA} sequences",
journal = "Nucl. Acids Res.",
volume = "25",
comment = "walker.tex, November 1, issue 21",
note = "http://www-lecb.ncifcrf.gov/$\sim$toms/paper/walker/,
erratum: NAR 26(4): 1135, 1998",
pages = "4408-4415",
year = "1997"}

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

author
   Thomas Dana Schneider

bugs
   * The quote strings in the parameter file are not recorded and so are not
   reproduced in the data file comments.
   * Blank characters are placed around the quote strings.

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.

   To provide upwards compatability, scanp files of version 2.90 or less will
   be interpreted by the old definitions for the bounds of Ri, Z and p:

      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.

   It is not advisable to rely on this feature, as it will go away at some
   point.

*)
(* end module describe.scan *)

(* constants continued *)

(* 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 = 7.70; {of delmod.p 2005 Sep 15} *)

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.70; {of delmod.p 2005 Sep 15} *)

      (* 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.70; {of delmod.p 2005 Sep 15} *)

(* begin module datetime.type *)
   (* array for dates *)
   datetimearray = packed array[1..datetimearraylength] of char;
(* end module datetime.type version = 7.70; {of delmod.p 2005 Sep 15} *)

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

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.70; {of delmod.p 2005 Sep 15} *)

      book,  (* a book from the sequence library *)
      ribl,  (* weight matrix *)
      scanp,  (* program parameters *)
      data,  (* result of the scan *)
      scanfeatures, (* result of the scan as features *)
      scaninst: (* result of the scan as delila instructions *)
        text;

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

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* ************************************************************************ *)
(* end module package.primitive version = 7.70; {of delmod.p 2005 Sep 15} *)

(* 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).
As of 2000 February 18, the Sun Pascal compiler requires a formatting
statement.  This statement allows the date to be generated in this
standard Delila format in a single call.  Information about the
formatting statement is available on the manual page for date in Unix.
If a computer does not have this method, see the 'oldgetdatetime' routine
in delmod.p (https://alum.mit.edu/www/toms/delila/delmod.html)
for some conversion code.

GPC Functions:
function  GetUnixTime (var MicroSecond : Integer) : UnixTimeType;

http://agnes.dida.physik.uni-essen.de/~gnu-pascal/gpc_109.html#SEC109

7.10.8 Date And Time Routines 

procedure GetTimeStamp (var t : TimeStamp); 
function Date (t : TimeStamp) : packed array [1 .. DateLength] of Char; 
function Time (t : TimeStamp) : packed array [1 .. TimeLength] of Char; 

DateLength and TimeLength are implementation dependent constants. 

GetTimeStamp (t) fills the record `t' with values. If they are valid, the Boolean
flags are set to True. 

TimeStamp is a predefined type in the Extended Pascal standard. It may be
extended in an implementation, and is indeed extended in GPC. For the full
definition of `TimeStamp', see section 8.255 TimeStamp. 

*)
var
    t: TimeStamp; 

(* begin module pluckdigit *)
function pluckdigit(number, logplace:integer): char;
(* return the digit at the place value ('logplace') position of number.
example:
      pluckdigit(13625, 3) = 3
      pluckdigit(13625, 4) = 1
This routine was taken from module numberdigit in prgmod.p, but is
modified so as not to give the sign.  Instead it gives zeros
above the digits.  'myabsolute' replaced 'absolute', which
is apparently a keyword for GPC.  The name is kept for
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 *)
begin (* pluckdigit *)
      place:=1;
      for count:=1 to logplace do place:=10*place;
      if number=0 then begin
         acharacter:='0'
      end
      else begin
         myabsolute:=number;
         if myabsolute >= place
         then digit
         else acharacter := '0'
      end;
      pluckdigit:=acharacter
end; (* pluckdigit *)
(* end module pluckdigit *)

begin

 (* according to:
http://agnes.dida.physik.uni-essen.de/~gnu-pascal/gpc_109.html#SEC109
*)
   GetTimeStamp(t); 

(* Predefined time stamp:
http://agnes.dida.physik.uni-essen.de/~gnu-pascal/gpc_389.html#SEC389
  TimeStamp = {@@packed} record
      DateValid,
      TimeValid   : Boolean;
      Year        : Integer;
      Month       : 1 .. 12;
      Day         : 1 .. 31;
      DayOfWeek   : 0 .. 6;  { 0 means Sunday }
      Hour        : 0 .. 23;
      Minute      : 0 .. 59;
      Second      : 0 .. 61; { to allow for leap seconds }
      MicroSecond : 0 .. 999999
    end;
*)

  with t do begin
    if TimeValid then begin
{
       writeln(output,'valid time');
       writeln(output,'year   =',year:4);
       writeln(output,'month  =',month:2);
       writeln(output,'day    =',day:2);
       writeln(output,'hour   =',hour:2);
       writeln(output,'minute =',minute:2);
       writeln(output,'second =',second:2);
}
       adatetime := 'year/mm/dd hh:mm:ss';
       adatetime[ 1] := pluckdigit(Year,3);
       adatetime[ 2] := pluckdigit(Year,2);
       adatetime[ 3] := pluckdigit(Year,1);
       adatetime[ 4] := pluckdigit(Year,0);
       adatetime[ 6] := pluckdigit(Month,1);
       adatetime[ 7] := pluckdigit(Month,0);
       adatetime[ 9] := pluckdigit(Day,1);
       adatetime[10] := pluckdigit(Day,0);
       adatetime[12] := pluckdigit(Hour,1);
       adatetime[13] := pluckdigit(Hour,0);
       adatetime[15] := pluckdigit(Minute,1);
       adatetime[16] := pluckdigit(Minute,0);
       adatetime[18] := pluckdigit(Second,1);
       adatetime[19] := pluckdigit(Second,0);
    end
    else begin
       writeln(output,'getdatetime: invalid time!');
       halt;
    end
  end;
{
   Sun compiler method:
   date(adatetime,'%Y/%m/%d %H:%M:%S');
}
end;
(* end module getdatetime version = 1.17; (@ of timegpc.p 2002 Oct 9 *)

(* 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.17; (@ of timegpc.p 2002 Oct 9 *)

(* 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.17; (@ of timegpc.p 2002 Oct 9 *)

(* 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.17; (@ of timegpc.p 2002 Oct 9 *)
 
(*[[*)
(* 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.17; (@ of timegpc.p 2002 Oct 9 *)
(*]]*)

(* ************************************************************************ *)
(* end module package.datetime version = 1.17; (@ of timegpc.p 2002 Oct 9 *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* ************************************************************************ *)
(* end module package.brpiece version = 7.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* ************************************************************************ *)
(* end module package.getocp version = 7.70; {of delmod.p 2005 Sep 15} *)

(* 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.70; {of delmod.p 2005 Sep 15} *)

(* 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 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 clearfeature *)
procedure clearfeature(var f: featuretypeptr);
(* clear out the feature f *)
begin
   with f^ do begin
      clearstring(id);
      coordinate := 0;
      orientation := +1;
      clearstring(nametag);
      clearstring(othertag);
      definition := nil;
      duplicateprevious := nil;
      Ri := 0.0;
      Z := 0.0;
      probability := 0.0;
      Tparam := ' ';
      Aparam := 0.0;
      Bparam := 0.0;
      Cparam := 0.0;
      Dparam := 0.0;
      relcoo := 0.0;
      evencoordinate := 0.0;
      unsatisfied := true;
      fromrange := 0;
      torange := 0;
      featurenumber := 0;
      desiredline := 0;
      others := nil;
      mothers := nil;
      usedInGroup := false; (* used in group *)
      next := nil;
   end
end;
(* end module clearfeature 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 showafeature *)
procedure showafeature(var afile: text; f: featuretypeptr);
(* show the feature f to file afile *)
var
   o: featurecarrierptr; (* other features that this one points to *)
   totals: integer; (* totals pointing to f *)
   gaps: integer; (* gaps pointing to f *)
begin
   totals := 0;
   gaps := 0;
   with f^ do begin
      writeln(afile,'FEATURE #', featurenumber:1,':');
      write(afile,'id: ');
      writequotestring(afile,id);
      writeln(afile);

      writeln(afile,'coordinate: ',coordinate:1:1);

      write(afile,'orientation: ');
      if orientation > 0 then writeln(afile,'+',orientation:1)
                         else writeln(afile,    orientation:1);

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

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

      write(afile,'unsatisfied: ');
      if unsatisfied then writeln(afile,'true')
                     else writeln(afile,'false');

      writeln(afile,'fromrange: ',fromrange:5:1);
      writeln(afile,'  torange: ',  torange:5:1);
      writeln(afile,'desiredline: ',  desiredline:5);
      if duplicateprevious = nil
      then writeln(afile,'NOT a duplicate feature')
      else writeln(afile,'DUPLICATE of feature #',
                   duplicateprevious^.featurenumber:1);

      if others <> nil then begin
         writeln(afile,'This feature points at these OTHERS:');
         write(afile,'feature -> ');
         o := others;
         while o <> nil do begin
{
write(afile,'-LOOP-');
if o = others then write(afile, '|o = others');
if o <> others then write(afile, '|o <> others');
if not (o = others) then write(afile, '|not (o = others)');
write(afile, '|');
}
{KKKKK}
            if o <> others then write(afile,', ');
            if o^.afeature <> nil then begin
               write(afile,'[');
               writeNOTABstring(afile,o^.afeature^.nametag);
               write(afile,']');
            end
            else write(afile,'[EMPTY]');
            o := o^.next;
         end;
         writeln(afile);
         o := others;
         while o <> nil do begin
            if o^.afeature <> nil then begin
               write(afile,'= ');
               writeafeature(afile,o^.afeature);
               writeln(afile);
            end;
            o := o^.next;
         end;
      end
      else begin
         writeln(afile,'This feature does not point to any OTHERS.');
      end;

      if mothers <> nil then begin
         writeln(afile,'The MOTHERS of this feature', 
                       ' (that point to this feature) are:');
         o := mothers;
         while o <> nil do begin
            if o^.afeature <> nil then begin
               write(afile,'feature <- "');
               writeNOTABstring(afile,o^.afeature^.nametag);
               writeln(afile,'" #',o^.afeature^.featurenumber:1);
{
               writeln(afile,'"');
}
               if o^.afeature^.definition^.istotal
               then totals := succ(totals);
               if o^.afeature^.definition^.isgap
               then gaps := succ(gaps);
            end;
{writeln(output,'BUBBA: mothers part 1 about to "o := o^.next;" ');}
            o := o^.next;
         end;
         o := mothers;
{writeln(output,'BUBBA: mothers part 1 is done');}
         while o <> nil do begin
            if o^.afeature <> nil then begin
               write(afile,'= ');
               writeafeature(afile,o^.afeature);
               writeln(afile);
            end;
            o := o^.next;
         end;
{writeln(output,'BUBBA: mothers part 2 is done');}
      end
      else begin
         writeln(afile,'This feature has no MOTHERS.');
      end;
{writeln(output,'BUBBA: showafeature is done');}

   end;
   write(afile,gaps:1,' gap');
   if gaps <> 1 then write(afile,'s'); (* stolen from plural *)
   writeln(afile);
   write(afile,totals:1,' total');
   if totals <> 1 then write(afile,'s'); (* stolen from plural *)
   writeln(afile);
   if (totals > 1) then begin
      writeln(output,'showafeature: more than one total!');
      halt;
   end;
   if (gaps >1) and (totals < 1) then begin
      writeln(output,'showafeature: ',gaps:1,' gaps but no total!');
{NNN}
      halt;
   end;
end;
(* end module showafeature 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,'* Paramters 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 scan.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);
         clearstring(prebuiltstring^);
         pbs := prebuiltstring;
      end
      else begin
         new(pbs^.next);
         pbs := pbs^.next;
         clearstring(pbs^);
      end;
      pointer := pbs
   end
end;

begin with parameters do begin
   reset(scanp);

   readln(scanp, parameterversion);
   if parameterversion < updateversion 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
{zzzmmm}
   end
   else begin
      sitedefined := false;
      sitedefinition := nil;
   end;
   readln(scanp);

   (* read the controls for the strings to be written *)
   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;

   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 = '.');

{uuu}

{The following code inserts things into the prebuilt string
to prove that one can do that
stringforRi^.letters[1] := 'R';
stringforRi^.length := 1;
stringforCoordinate^.letters[1] := 'C';
stringforCoordinate^.length := 1;
writeln(output);
write  (output,'prebuiltstring = "');
writestringlines(output, prebuiltstring);
writeln(output,'"');
writeln(output,'JOINED:');
clearstring(otherquotestring);
joinstring(prebuiltstring,otherquotestring);
write(output,'prebuiltstring ');
writequotestring(output,otherquotestring);
clearstring(otherquotestring);
writeln(output);
}

   (* 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 *)
   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 scan.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 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'; *)

(* begin module scan.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);
(* 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. *)
var
   b: base; (* the base at position l *)
   i: integer; (* internal coodinates for a piece *)
   f: featuretypeptr; (* a feature found *)
   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 *)
   Ri0: real; (* Ri normalized by Ri at zero coordinate *)
   tol: real; (* tolerance in the resulting p *)
   thelength: integer; (* the length of apiece *)
   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 evalutation at coordinate i out if it passes the criterion.
specify the orientation *)
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 coordinate *)
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;

   (* 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
      sites := sites + 1;

      (* coordinate of the evaluation *)
      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;

      if sitedefined then with parameters do begin

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

         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);

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

         f^.definition^.matrix := riblmatrix;

         if f^.definition^.matrix^.symmetry = 'e'
         then begin
            f^.evencoordinate := f^.coordinate + orientation*0.5;
         end
         else f^.evencoordinate := f^.coordinate;

         if outfeat then begin
            writeafeature(features, f);
            writeln(features);
(*
showafeature(output, f);{zzzqqq}
*)
         end;

      end;

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

         write(inst,'name ');
         writequotestring(inst,nametag);
         write(inst,';');

         write(inst,' piece ');
         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;

   (* build the feature *)
   if sitedefined then begin
      (* build as much of the feature to be written out as we can *)
      new(f);
      clearfeature(f);
      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
      *)
      if (0 < inttopie(-fromwanted+1,apiece)) or
         (0 > inttopie(piecelength(apiece), apiece)-towanted)
      then begin
         writeln(output, 'piece without zero base found; cannot normalize');
         halt
      end;

      i := pietoint(0, apiece);
      Ri := 0.0;
      for l := fromwanted to towanted do begin
         b := getbase(i, apiece);
         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 *)
      startpoint := -fromwanted + 1;
      if startpoint < 1 then startpoint := 1;
      endpoint := thelength-towanted;
      if endpoint > thelength then endpoint := thelength;
      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;
         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 *)
         startpoint := towanted + 1;
         if startpoint < 1 then startpoint := 1;
         endpoint := thelength + fromwanted;
         if endpoint > thelength then endpoint := thelength;
         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

   else begin (* handle CIRCULAR case *)
      startpoint := 1;
      endpoint := thelength;
      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 scan.scansequence *)

(* 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 digitcounter *)
procedure digitcounter(var afile: text; var counter, previous: integer);
(* Make a digit counter on a single line.  The digits are written and
overwritten the next time by using control H.  To initialize, set the
'previous' variable to zero.  To terminate, put a
"writeln(afile,'whateveryouwant');".  If you put
"write(afile,'babble-counter: ');" before the call, the numbers will
show up on the same line and your message can indicate what the
counter is about. *)
const
   ordControlH = 8; (* the ordinal of control H so that we do not
      have to put an actual control H into the source code.
      The chr() of this will give the control H character *)
var
   digits: integer; (* number of digits in in previous number *)
   digitcount: integer; (* index to digits *)
begin
   if previous > 0 then begin
      (* backspace the PREVIOUS number of digits *)
      (* Compute the number of digits.  For sequencenumber = 1000, we
      add an extra 0.5 to avoid rounding that would give 3 instead of
      4 as the number of digits.  Other cases work fine with this. *)
      digits :=  trunc(ln(previous+0.5)/ln(10)+1);
{
writeln(output, digits:1, ' ', counter:1);
      for digitcount := 1 to previous do write(output, 'H');
}
      for digitcount := 1 to digits do write(output, chr(ordControlH));
{
writeln(output, digits:1, ' ', counter:1, ' ',previous:1);
}
   end;
   write(output, counter:1);
   previous := counter;
end;
(* end module digitcounter version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module scan.themain *)
procedure themain(var book, ribl, scanp, data, scanfeatures, scaninst: text);
(* the main procedure *)
var
   count: integer; (* number of sequences done so far *)
   previouscount: integer; (* previous count value *)
   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 *)
   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 *)

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
   writeln(output, 'scan ',version:4:2);
   PatentMessage(output);
   limitdate ('9','9','9','9','    /99/99 99:99:99');

   (* initialize the variables *)
   reset(ribl);
   getriblmatrix(ribl,riblmatrix);

   readparameters(scanp, riblmatrix, parameters);

   (* compare the riblmatrix name to the sitedefinition name *)
   if not equalstring(parameters.sitedefinition^.nametag,
              riblmatrix^.riblname)
   then begin
      writeln(output,'====================================================');
{qqq}
      write(output,'site name in scanp: ');
      writequotestring(output,parameters.sitedefinition^.nametag);
      writeln(output);
      write(output,'site name in  ribl: ');
      writequotestring(output,riblmatrix^.riblname);
      writeln(output);
      writeln(output,'These names should match.');
      writeln(output,'====================================================');
      halt;
   end;

   if parameters.outdata then begin
      rewrite(data);
      (* by making the first line as follows, alist can use the file directly *)
      writeln(data,'* bits');
      writeln(data, '* scan ',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);
   end
   else writeparameters(output, parameters);

   if parameters.outinst then begin
      rewrite(scaninst);
      writeln(scaninst, 'title "scan ',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;

   if parameters.outfeat then begin
      with parameters do if sitedefined then begin
         rewrite(scanfeatures);
         writeln(scanfeatures, '* scan ',version:4:2);
         reset(book);
         copyaline(book, scanfeatures);
         writeln(scanfeatures); (* blank line *)
         writeadefinition(scanfeatures,sitedefinition);
         writeln(scanfeatures); (* blank line *)
{
zzzmmm this is done in writeadefinition! NO
but sitedefinition needs to link to riblmatrix!
}
         putriblmatrix(scanfeatures,riblmatrix)
      end;
   end;

   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 piece number');
         writeln(data,'* 2 piece length');
         writeln(data,'* 3 piece name');
         writeln(data,'* 4 piece coordinate');
         writeln(data,'* 5 matrix orientation',
                      ' (+1 = as in book, -1 = complement)');
         writeln(data,'* 6 Ri evaluation (bits per site)');
         writeln(data,'* 7 Z');
         writeln(data,'* 8 probability');

         with parameters do
         if normalboolean then begin
            write(data,'* 9 normalized Ri:');
            case normalizeRi of
               's': writeln(data,' ',' delta Ri');
               'd': writeln(data,' ',' div Ri');
            end;
         end
         else writeln(data,'* 9 0, meaning: not normalizing Ri');

         writeln(data,'*');
      end;

      tellsymmetry(output,oneway,ways);

      totalsites := 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;

      apiece := nil; (* getocp requires the piece to start at nil *)
      orgchange := false; orgopen := false;
      chrchange := false; chropen := false;
      piechange := false; pieopen := false;

      while (not eof(book)) and (count < todo) do begin
         getocp(book, theline,
                org,    orgchange, orgopen,
                chr,    chrchange, chropen,
                apiece, piechange, pieopen);

         if not eof(book) then begin
            count := count + 1;
            if (count mod countmark) = 0
            then digitcounter(output, count, previouscount);

            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 writeln(scaninst,' ');
                  write(scaninst,'chromosome ');
                  writename(scaninst,chr.hea.keynam);
                  write(scaninst,';');
               end;
               if chrchange or chrchange
               then writeln(scaninst);
            end;

            scansequence(data, scanfeatures, scaninst,
                         outdata, {outfeat,} outinst,
                         riblmatrix, fromwanted, towanted,
                         apiece,
                           cutDNRi,
                           cutUPRi,
                           cutDNZ,
                           cutUPZ,
                           cutDNP,
                           cutUPP,
                         sites,positionsevaluated,
                         not oneway, riblmatrix^.mean, riblmatrix^.stdev,
                         sitedefined, sitedefinition,
                         parameters);

            totalsites := totalsites + sites;
            totalpositionsevaluated := totalpositionsevaluated
                                     + positionsevaluated;
            if outdata then begin
               writeln(data,'* ',sites:1,' sites found in piece',
                          ' ',apiece^.key.hea.keynam.letters);
               writeln(data,'* ',positionsevaluated:1,' positions evaluated');
            end;
            clearpiece(apiece); (* clear the piece for reuse *)
         end
      end;

      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 scan.themain *)
begin
   themain(book, ribl, scanp, data, scanfeatures, scaninst);
1: end.
