program delmod(book, output);
(* delila module library

  Dr. Thomas D. Schneider
  National Institutes of Health
  National Cancer Institute
  Center for Cancer Research Nanobiology Program
  Molecular Information Theory Group
  Frederick, Maryland  21702-1201
  toms@ncifcrf.gov
  permanent email: toms@alum.mit.edu (use only if first address fails)
  http://www.ccrnp.ncifcrf.gov/~toms/

  Thomas Schneider and Gary Stormo *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 7.72; {of delmod.p 2007 Jul 23}
(*
2007 Jul 23: 7.72: cleanup
2007 Jun 22: 7.71: "The piece name in the book" - ignore dotted names
2005 Sep 15: 7.70: linelength increased, brline protected
2005 Sep 13: 7.69: spelling correction: asterisk
2005 Jun 23: 7.68: orderseed is now outside timeseed package
2004 Sep  8: 7.67: upgrade skipblanks to match prgmod
2004 Aug  4: 7.66: upgrade cleardna documentation
2004 Jul 27: 7.65: bug: "The piece name in the book: ... does not match the
                   inst file name ..."  The inst file name was all dots.
                   used variable p instead of p1 of for loop.  fixed.
2004 Jul 18: 7.64: dircomplement/dirhomologous handled
2003 Jan 13: 7.63: getocp must handle mapunit
2003 Jan 13: 7.61: bworg and bwchr wrote to 'book' instead of thefile
2003 May  3: 7.60: ex0bk description improved.
2002 Sep  5: 7.59: cleanup
2001 Mar 16: 7.58: clearname has a var variable!
2001 Mar 15: 7.57: improve getocp documentation about how to initialize
2000 Nov 21: 7.56: improve getocp documentation about how to initialize
2000 Nov 16: 7.55: make it skip quotes in proc.align.
2000 Nov 16: 7.54: fix reading of name and set instructions in proc. align
2000 Jul 30: 7.53: use gpc time modules!
2000 Jul 30: 7.52: no basic longer source of time functions.
2000 Jul 23: 7.51: update documentation
2000 Jul 23: 7.50: gut showfreedna - incompatable with gpcc
2000 Jun 26: 7.49: introduce iwgetsimple
2000 Jun 26: 7.48: writeline to write lineptr type
2000 Jun 26: 7.47: upgrade iwget to control return
2000 Jun 26: 7.46: upgrade iworgchr to include open booleans
2000 Apr 19: 7.43: duplicate parts in bworg removed (it's in bworgkey)
2000 Mar 22: 7.42: bug in getbase: now does circular grabs for safety.
2000 Feb 18: 7.40: final solution to y2k problem: date call with formatting
1999 Dec 13: 7.38: Solve y2k bug for 10 years (see getdatetime)
1999 July 24: 7.37: adjust maxminalignment for first base alignment
1999 July 13: 7.33: revamp getbase and putbase to work from length
     of dna parts rather than dnamax.  This solves bugs in delila.
1999 May 5: 7.14: added crash function
1999 April 28: 7.11: iwget same function taken from search.
1999 April 27: 7.09: fixed equalname: need to set i initially
1999 Mar 13: 6.99: brpiece assumes that one starts at the word 'piece'.
     This means that it now reads the complete object.
1999 Mar 9: 6.95: brpiece and other read routines give change in line
     number in the file that they read.  This allows delila to
     use these routines.  MANY programs will have to be changed...
     but it's not a hard change.
1999 Mar 4: strings were added so that the program compiles
1999 Mar 4: timeseed dates were corrected
1998 Aug 1: getocp has to set 
1998 Jan 26: namelength set to 100 to allow long names
1998 Jan 4: dnamax set to 10 million for faster scans
1997 April 22: Year 2000 date solved; program
   accepts old format AND new format *)
(* end module version *)

(* begin module describe.delmod *)
(*
name

   delmod: delila module library

synopsis

   delmod(book: in, output: out)

files

   book: any book from the delila system, or an empty file.

   output: the version of delmod is printed along with test results
      if the book is not empty.  Successful compilation and running of the
      program indicates that the modules are correct.

description

   Delmod is a collection of modules used by delila system programs.  The
   easiest way to obtain a list of the modules is to run the module program
   using delmod for both sin and modlib (with dummy files for the other
   input).  There are a number of information modules, indicated by names
   beginning with 'info.'.  There are also a number of packages of modules
   that pickup other modules.  These begin with 'package.'.  You should note
   that some modules are constants, others types, etc.  These must remain in
   their proper location to allow compilation.

   The delmod program will report the current date and time in the standard
   Delila format.

   The delmod program will read a delila book if the 'book' file is not empty
   and report some information about it.

examples

   You can use an empty file for book or a good book to use to test
   delmod is ex0bk.  To do this, get the ex0bk file (see link below)
   and then rename it 'book'.

see also

   module.p {The module program is use to transfer modules from this library
   to other programs.}

   {Example book for testing:} ex0bk

   {Other module libraries:} prgmod.p matmod.p

   {General discussion on compiling Delila programs:} 
   http://www.lecb.ncifcrf.gov/~toms/delila.html#How.To.Compile 

   {Time modules may be moved into this module library, but
   the ultimate source is now in three places, not here: }
   timesun.p
   timep2c.p
   timegpc.p

author

   Thomas D. Schneider and Gary D. Stormo

bugs

technical notes

*)
(* end module describe.delmod *)

(* packages of modules ******************************************************)
(* begin module info.package *)
(* information on packages of modules

package.getpiece
      picks up procedure getpiece and the procedures it requires.
      this allows one to scan a book and detect pieces.
      the user must still specify const, type, var and procedures
      to pickup bases from the pieces.

package.getocp
      picks up procedure getocp and the procedures it requires.
      this allows one to scan a book and detect organisms, chromosomes and
      pieces.  the user must still specify const, type, var and procedures
      to pickup bases from the pieces.

package.nextbase
      this package allows a programmer to pull out bases from a book with
      little programming effort (at the expense of flexability).
      see documentation in module info.nextbase.

package.bwrite
      these routines allow one to write out parts of a book, without worry
      (for the most part) on the structure of the book.

package.iwrite
      picks up procedures that allow one to easily write delila instructions.
      the typical use is as output from a book searching program.
      note: module copylines is required above this package.

package.align
      picks up procedure align and the procedures it uses.  it does not pick
      up the other align procedures.  align (and associated procedures, if
      pickedup) allow one to read a book and instructions for the book
      concurrently.  the pieces of the book are returned along with the
      aligning base (in internal coordinates).

package.datetime
      interfaces to the system clock.  these make date and time calls
      system independent.
      THIS PACKAGE IS NOW TO BE TAKEN FROM ONE OF THE
      TIME PROGRAMS:  timesun.p, timegpc.p timep2c.p

package.primitive
      this is a package of primitive programming functions that many programs
      will use.  it includes a halt function and line copying functions.
*)
(* end module info.package *)

(* begin module package.brpiece *)
(* ************************************************************************ *)
(* begin module book.basis *)
(* end module book.basis *)

(* begin module book.getto *)
(* end module book.getto *)

(* begin module book.skipstar *)
(* end module book.skipstar *)

(* begin module book.brreanum *)
(* end module book.brreanum *)

(* begin module book.brnumber *)
(* end module book.brnumber *)

(* begin module book.brname *)
(* end module book.brname *)

(* begin module book.brline *)
(* end module book.brline *)

(* begin module book.brdirect *)
(* end module book.brdirect *)

(* begin module book.brconfig *)
(* end module book.brconfig *)

(* begin module book.brnotenumber *)
(* end module book.brnotenumber *)

(* begin module book.brnote *)
(* end module book.brnote *)

(* begin module book.brheader *)
(* end module book.brheader *)

(* begin module book.copyheader *)
(* end module book.copyheader *)

(* begin module book.brpiekey *)
(* end module book.brpiekey *)

(* begin module book.brdna *)
(* end module book.brdna *)

(* begin module book.brpiece *)
(* end module book.brpiece *)

(* begin module book.brinit *)
(* end module book.brinit *)

(* ************************************************************************ *)
(* end module package.brpiece *)

(* begin module package.getpiece *)
(* ************************************************************************ *)
(* begin module package.brpiece *)
(* end module package.brpiece *)

(* begin module book.getpiece *)
(* end module book.getpiece *)

(* ************************************************************************ *)
(* end module package.getpiece *)

(* begin module package.getocp *)
(* ************************************************************************ *)
(* begin module package.brpiece *)
(* end module package.brpiece *)

(* begin module book.brorgkey *)
(* end module book.brorgkey *)

(* begin module book.brchrkey *)
(* end module book.brchrkey *)

(* begin module book.getocp *)
(* end module book.getocp *)

(* ************************************************************************ *)
(* end module package.getocp *)

(* begin module package.nextbase *)
(* ************************************************************************ *)
(* begin module package.getpiece *)
(* end module package.getpiece *)

(* begin module book.stepbase *)
(* end module book.stepbase *)

(* begin module nextbase *)
(* end module nextbase *)

(* ************************************************************************ *)
(* end module package.nextbase *)

(* begin module package.bwrite *)
(****************************************************************************)
(* this is a package of procedures for writing books,
   by gary stormo,  aug 17, 1982 *)

(* begin module book.bwbasics *)
(* end module book.bwbasics *)

(* begin module book.bworg *)
(* end module book.bworg *)

(* begin module book.bwchr *)
(* end module book.bwchr *)

(* begin module book.bwdna *)
(* end module book.bwdna *)

(* begin module book.bwpie *)
(* end module book.bwpie *)

(* begin module book.bwref *)
(* end module book.bwref *)

(* begin module book.bwgen *)
(* end module book.bwgen *)

(* begin module book.bwtra *)
(* end module book.bwtra *)

(* begin module book.bwmar *)
(* end module book.bwmar *)

(****************************************************************************)
(* end module package.bwrite *)

(* begin module package.iwrite *)
(* ************************************************************************ *)
(* begin module book.iwcombk *)
(* end module book.iwcombk *)

(* begin module book.iwname *)
(* end module book.iwname *)

(* begin module book.iworg *)
(* end module book.iworg *)

(* begin module book.iwchr *)
(* end module book.iwchr *)

(* begin module book.iwpie *)
(* end module book.iwpie *)

(* begin module book.iworgchr *)
(* end module book.iworgchr *)

(* begin module book.iwget *)
(* end module book.iwget *)

(* begin module book.iwget2 *)
(* end module book.iwget2 *)

(* begin module book.iwgetsimple *)
(* end module book.iwgetsimple *)

(* ************************************************************************ *)
(* end module package.iwrite *)

(* begin module LOCKED.package.datetime *)
(* ************************************************************************ *)
(* as of 2000 July 30,
the source of package.datetime is now the timeXXX.p programs *)
(* begin module getdatetime *)
(* end module getdatetime *)

(* begin module readdatetime *)
(* end module readdatetime *)

(* begin module writedatetime *)
(* end module writedatetime *)

(* begin module timeseed *)
(* end module timeseed *)
 
(*[[*)
(* begin module limitdate *) 
(* end module limitdate *)

(*]]*)
(* ************************************************************************ *)
(* end module LOCKED.package.datetime *)

(* begin module package.primitive *)
(* ************************************************************************ *)
(* begin module halt *)
(* end module halt *)

(* begin module copyaline *)
(* end module copyaline *)

(* begin module copylines *)
(* end module copylines *)

(* ************************************************************************ *)
(* end module package.primitive *)

(* begin module package.align *)
(* ************************************************************************ *)
(* begin module package.getpiece *)
(* end module package.getpiece *)

(* begin module findblank *)
(* end module findblank *)

(* begin module findnonblank *)
(* end module findnonblank *)

(* begin module align.align *)
(* end module align.align *)

(* ************************************************************************ *)
(* end module package.align *)


(* begin module book.const ***************************************************)
(* constants needed for book manipulations *)
      dnamax = 1024;  (* length of dna arrays *)
      namelength = 100; (* maximum key name length *)
      linelength = 200; (* maximum line readable in book *)
(* end module book.const *)

(* begin module datetime.const *)

   datetimearraylength = 19; (* length of dataarray for dates,
       It is just long enough to include the 4 digit year - solving the
       year 2000 problem:
1980/06/09 18:49:11
123456789 123456789 
         1         2
*)

(* end module datetime.const version = 1.15; (@ of timegpc.p 2000 Oct 11 *)

(* trigger definitions follow to make program compilable.
The prgmod.p program has the triggers *)
(* module filler.const *)
      fillermax = 50; (* the size of the filler array for a string *)
(* module filler.const from prgmod.p 4.20 *)

(* begin module interact.const *)
(* begin module string.const *)
      maxstring = 2000; (* the maximum string *)
(* end module string.const version = 4.86; (@ of prgmod.p 2004 Sep 8 *)
(* end module interact.const version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

type

(* begin module datetime.type *)
   (* array for dates *)
   datetimearray = packed array[1..datetimearraylength] of char;
(* end module datetime.type version = 1.15; (@ of timegpc.p 2000 Oct 11 *)

(* begin module base.type *)
   (* define the four nucleotide bases *)
   base = (a,c,g,t);
(* end module base.type version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* 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 *)
(* end module base.type *)

      (* 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 = 1.01; (@ of testtime.p 1997 Jan 11 *)

(* begin module amino.type *)
      aminoacid = packed array[1..3] of char;
(* end module amino.type *)

(* trigger definitions follow to make program compilable.
The prgmod.p program has the triggers *)

(* begin module interact.type *)
(* begin module string.type *)
   stringptr = ^string; (* pointer to a string *)
   string = record (* a string of characters *)
      letters: array[1..maxstring] of char; (* the letters in the string *)
      length: integer; (* the number of characters in the string *)
      current: integer; (* the letter we are working on *)
      next: stringptr; (* the next string in a series *)
   end;
(* end module string.type version = 4.86; (@ of prgmod.p 2004 Sep 8 *)
(* end module interact.type version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* module filler.type *)
      (* the following is an array used to fill a string.
      it is convenient to have it much shorter than the maxstring, so that
      it is easy to fill the string using procedure fillstring.
      the user must declare the value of constant fillermax. *)
      filler = packed array[1..fillermax] of char;
(* module filler.type from prgmod.p 4.20 *)
 
(* module trigger.type *)
      trigger = record (* an object to be searched for *)
         seek: string; (* the characters looked for *)
         state: integer; (* how close to triggering we are *)
         skip: boolean; (* trigger not found- skip the line *)
         found: boolean (* the trigger was found *)
      end;
(* module trigger.type from prgmod.p 4.20 *)

var
      book: text; (* for testing delmods *)

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


(* ************************************************************************ *)
(* primitive procedures and functions for everyday use *)

(* begin module crash *)
procedure crash;
(* Crash the program by trying to open a nonexistant file.  This allows
tracing by the dbx program.  To use:  insert call into the halt program
or whereever a traceable stop is desired. *)
var
   bogus: text; (* boghous internal file *)
begin
  writeln(output,' program crash.');
  reset(bogus);
end;
(* end module crash *)

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

(* begin module skipblanks *)

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

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

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

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

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

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

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

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

{
This module is commented out for the moment, since strings
are not (yet) defined in delmod.  But the module
program can pick it up perfectly well, since it doesn't
know about this kind of comment!!
(* begin module isnamestring *)
function isnamestring(n: name; s: string): boolean;
(* is the delila-type name n the same as the string s? *)
var
   c: integer; (* index to n and s *)
   done: boolean; (* done checking *)
begin
(*
write(output,'isnamestring s(',s.length:1,'):');
writestring(output,s);
write(output,', n(',n.length:1,'):');
for c := 1 to n.length do write(output,n.letters[c]);
writeln(output);
*)
   if n.length = s.length then begin
      c := 1;
      done := false;
      while not done do begin
(*
writeln(output,'   ',s.letters[c],' ' ,n.letters[c]);
*)
         if s.letters[c] = n.letters[c]
         then begin
           c := succ(c);
           if c > s.length then begin
              isnamestring := true;
              done := true
           end
         end
         else begin
            isnamestring := false;
            done := true
         end
      end
   end
   else isnamestring := false;
end;  
(* end module isnamestring *)
}


(*****************************************************************************)
(* dating procedures *********************************************************)
(*****************************************************************************)

(* 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 (http://www.lecb.ncifcrf.gov/~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.15; (@ of timegpc.p 2000 Oct 11 *)

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

(* begin module oldgetdatetime *)
{ This old module contains methods for converting dates, but it
is messy compared to the new one installed on 2000 Feb 18.
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
X THE ROUTINE IS NEUTRALIZED BY A COMMENT TO PREVENT THE COMPLER FROM X
X SEEING THE NOW "UNSAFE" datetime(adate) CALL.                       X
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
procedure oldgetdatetime(var adatetime: datetimearray);
(* get the date and time into a single array from the system clock *)
(* this procedure works on the pascal compiler of a Pyramid 90x computer *)
(* adatetime contains the date:
       1980/06/09 18:49:11
         ye mo da ho mi se
(year, month, day, hour, minute, second) *)
var
      adate, atime: alfa; (* ie, packed array[1..10] of char; *)
(*
      (* the month array is only used on computers where month is expreessed
         in characters *)
      month: packed array[1..3] of char;
*)
      index: integer; (* index for times *)
begin
      (* clear *)
      for index:=1 to 20 do adatetime[index]:='?';

      date(adate);
      time(atime);
(* use the following line to determine the characters returned by the
date and functions.  Comment out the line when this routine is in
normal use *)
(*
writeln(output,'getdatetime: adate[',adate,'] atime[',atime,']');
*)

(* Unix was not (as of 1997 Mar 16) supplying the
"19", so we put it in ourselves:
      adatetime[1]:='1';
      adatetime[2]:='9';

However, even by 1999 Dec 13 Sun Pascal is *still* not Y2K compliant, as
shown by this code:

writeln(output,'adate: "',adate,'"');
gives:
adate: "12/13/99  "

adate: "12/13/99  "
atime: "16:57:46  "

Here is a workaround that will last 10 years:
*)

      (* Have we passed y2k? *)
      if ((adate[7]='9') and (adate[8]='9'))
      then begin
         adatetime[1]:='1';
         adatetime[2]:='9';
      end
      else begin
         adatetime[1]:='2';
         adatetime[2]:='0';
      end;

      (* year *)
      for index:=1 to 2 do adatetime[index+2]:=adate[index+6];
      adatetime[3+2]:='/';

      (* month *)
      for index:=4 to 5 do adatetime[index+2]:=adate[index-3];
      adatetime[6+2]:='/';

      (* day *)
      for index:=7 to 8 do adatetime[index+2]:=adate[index-6+3];
      if adatetime[7+2] = ' ' then adatetime[7+2] := '0'; (* safety *)
      adatetime[9+2]:=' ';

      (* time *)
      for index:=10 to 17 do adatetime[index+2]:=atime[index-9];

(*
writeln(output,'final result: ',adatetime);
*)

(*
Suddenly on 1997 January 11 this no longer worked under UNIX!!!
      (* year *)
      for index:=1 to 2 do adatetime[index]:=adate[index+7];
      adatetime[3]:='/';
      (* month *)
      for index:=4 to 6 do month[index-3]:=adate[index];
           if month='Jan' then begin adatetime[4]:='0'; adatetime[5]:='1' end
      else if month='Feb' then begin adatetime[4]:='0'; adatetime[5]:='2' end
      else if month='Mar' then begin adatetime[4]:='0'; adatetime[5]:='3' end
      else if month='Apr' then begin adatetime[4]:='0'; adatetime[5]:='4' end
      else if month='May' then begin adatetime[4]:='0'; adatetime[5]:='5' end
      else if month='Jun' then begin adatetime[4]:='0'; adatetime[5]:='6' end
      else if month='Jul' then begin adatetime[4]:='0'; adatetime[5]:='7' end
      else if month='Aug' then begin adatetime[4]:='0'; adatetime[5]:='8' end
      else if month='Sep' then begin adatetime[4]:='0'; adatetime[5]:='9' end
      else if month='Oct' then begin adatetime[4]:='1'; adatetime[5]:='0' end
      else if month='Nov' then begin adatetime[4]:='1'; adatetime[5]:='1' end
      else if month='Dec' then begin adatetime[4]:='1'; adatetime[5]:='2' end;
      adatetime[6]:='/';
      (* day *)
      for index:=7 to 8 do adatetime[index]:=adate[index-6];
      if adatetime[7] = ' ' then adatetime[7] := '0'; (* safety *)
      adatetime[9]:=' ';
*)
end;
}
(* end module oldgetdatetime version = 1.01; (@ of testtime.p 1997 Jan 11 *)

(* 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
 *)
(* 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,' old datetime (only 2 year digits) read: ');
      for index:=1 to datetimearraylength do write(thefile,adatetime[index]);
      writeln(output);
   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(thefile,adatetime[index]);
      writeln(output);
      halt
   end;
end;
(* end module readdatetime version = 1.15; (@ of timegpc.p 2000 Oct 11 *)

(* 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.15; (@ of timegpc.p 2000 Oct 11 *)

(* 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 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.15; (@ of timegpc.p 2000 Oct 11 *)

(* begin module orderseed *)
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;
(* end module orderseed version = 1.15; (@ of timegpc.p 2000 Oct 11 *)

(*[[*)
(* 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: http://www.lecb.ncifcrf.gov/~toms/walker/contacts.html');
       halt
   end
end;
(* end module limitdate version = 1.15; (@ of timegpc.p 2000 Oct 11 *)
(*]]*)

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

(* begin module book.name *)
procedure clearname(var n: name);
(* clear the name n *)
var
  i: integer; (* index to piece name *)
begin
  n.length := 0;
  for i := 1 to namelength do n.letters[i] := ' ';
end;

procedure writename(var f: text; n: name);
(* write the name n to file f *)
var
  i: integer; (* index to piece name *)
begin
   for i := 1 to n.length do
       write(f,n.letters[i]);
end;

procedure copyname(a: name; var b: name);
(* copy name a to name b *)
var
  i: integer; (* index to piece name *)
begin
   for i := 1 to a.length
   do b.letters[i] := a.letters[i];
   b.length := a.length;
end;

function equalname(a,b: name): boolean;
(* is name a equal to name b? *)
var
  i: integer; (* index to piece name *)
  same: boolean; (* temporary variable to hold the answer *)
begin
   same := true; (* what optimism!! *)
   if b.length = a.length
   then begin
      i := 1;
      while same and (i <= a.length) do begin
         same := (b.letters[i] = a.letters[i]);
         i := succ(i);
      end
   end
   else same := false;
   equalname := same;
end;
(* end module book.name *)

(* begin module book.other *)
(* other useful dna, line and name manipulating functions *)
(* begin module book.name *)
(* end module book.name *)

procedure emptydna(var l: dnaptr);
(* empty all of the dna in l onto the freedna list *)
begin
   while l<>nil do cleardna(l);
end;

procedure emptyline(var l: lineptr);
(* empty all of the line in l onto the freeline list *)
begin
   while l<>nil do clearline(l);
end;

procedure copyline(fromline: lineptr; var toline: lineptr);
(* copy a set of lines *)
var
   f,t : lineptr; (* internal pointers *)
begin
   (* destroy previous to line *)
   emptyline(toline);
   if fromline <> nil then begin
      getline(toline);
      f:=fromline;
      t:=  toline;
      while f<>nil do begin
         t^.letters:=f^.letters;
         t^.length:=f^.length;
         f:=f^.next;
         if f<>nil
         then begin
            getline(t^.next);
            t:=t^.next
         end
         else t^.next:=nil
      end
   end
end;

procedure copydna(var fromdna,todna: dnaptr);
(* copy a set of dna *)
var
   adna, memdna: dnaptr;
begin
   (* destroy previous to dna *)
   while todna<>nil do cleardna(todna);
   adna:=fromdna;
   if adna<>nil then begin
      getdna(memdna);
      todna:=memdna;
      while adna<>nil do begin
         memdna^.part:=adna^.part;
         memdna^.length:=adna^.length;
         adna:=adna^.next;
         if adna<>nil
         then begin
            getdna(memdna^.next);
            memdna:=memdna^.next
         end
         else memdna^.next:=nil
      end
   end;
end;

procedure copypiece(Apiece: pieceptr; var Bpiece: pieceptr);
(* copy the pice a to piece b *)
begin
   Bpiece^.key.hea.keynam := Apiece^.key.hea.keynam;
   copyline(Apiece^.key.hea.fulnam, Bpiece^.key.hea.fulnam);
   copyline(Apiece^.key.hea.note,   Bpiece^.key.hea.note);
   Bpiece^.key.mapbeg := Apiece^.key.mapbeg;
   Bpiece^.key.coocon := Apiece^.key.coocon;
   Bpiece^.key.coodir := Apiece^.key.coodir;
   Bpiece^.key.coobeg := Apiece^.key.coobeg;
   Bpiece^.key.cooend := Apiece^.key.cooend;
   Bpiece^.key.piecon := Apiece^.key.piecon;
   Bpiece^.key.piedir := Apiece^.key.piedir;
   Bpiece^.key.piebeg := Apiece^.key.piebeg;
   Bpiece^.key.pieend := Apiece^.key.pieend;
   copydna(Apiece^.dna,Bpiece^.dna);
end;

function between(a,b,c: integer): boolean;
(* is b between a and c? *)
(* this is an inclusive between *)
begin
      between:=((a<=b) and (b<=c))
                       or
               ((c<=b) and (b<=a))
end;

function within(pie: pieceptr; p: integer): boolean;
(* is p (external coordinates) within the piece pie? *)
(* note 1: if coocon is linear then piecon must be linear.
   note 2: does the piece not go over the coordinate boundaries?
   note 3: if coocon is circular and piecon is circular,
           then one has the entire piece, so we can ask if p is within
           the coordinate system.
*)
begin
      with pie^.key do begin
         case coocon of
            linear: within:=between(piebeg,p,pieend); (* note 1 *)
            circular: case piecon of
               linear: case piedir of
                  dirhomologous, (* handle case, may not be right *)
                  plus:  if pieend>=piebeg (* note 2 *)
                         then within:=between(piebeg,p,pieend)
                         else within:=between(piebeg,p,cooend) or
                                      between(coobeg,p,pieend);
                  dircomplement, (* handle case, may not be right *)
                  minus: if pieend<=piebeg (* note 2 *)
                         then within:=between(piebeg,p,pieend)
                         else within:=between(piebeg,p,coobeg) or
                                      between(cooend,p,pieend)
               end;
               circular:  within:=between(coobeg,p,cooend) (* note 3 *)
            end
         end
      end
end;

function withininternal( pie: pieceptr; p: integer): boolean;
(* Is the internal position p inside the piece pie? *)
begin
   withininternal := (p >= 1) and
                     (p <= piecelength(pie))
end;

(* end module book.other *)

{
********************************************************************************

In this section are old versions of getbase.

(* OLDbegin module book.getbase *)
This version was from before 1999 June
function getbase(position: integer; pie: pieceptr):base;
(* get a base from the nth position (internal coordinates) of the
   piece.  no protection is made against positions outside the piece *)
var
      workdna: dnaptr;
      p: integer; (* the last base of the dna part *)
begin
      workdna:=pie^.dna;
      p:=dnamax;
      while position>p do begin
         p:=p+dnamax;
         workdna:=workdna^.next
      end;
      getbase:=workdna^.part[position-(p-dnamax)]
end;
(* OLDend module book.getbase *)

(* ANOTHER OLD begin module book.getbase *)
This version was until march 22, when the lister program
showed a bomb for grabbing a circle in libdef.test:

title "version = 1.08 of libdef.inst 1999 June 11  Testing Delila";
organism E.xamples;
chromosome E.xamples;
piece CIRCULAR;
name "get all;"; get all piece;

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. *)
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 *)
begin
(*
writeln(output,'NEW getbase: position=',position:1,'^^^^^^^^^^^^^^^^^^^^');
*)
   if position < 0 then begin
      writeln(output,'error in getbase: request position (= ', position:1,
                     ') before end of piece');
      halt
   end;
   workdna:=pie^.dna;
   p:=workdna^.length;
   while position > p do begin
(*
writeln(output,' workdna^.length=',workdna^.length:1);
*)
      workdna := workdna^.next;
      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;
(* ANOTHER OLD end module book.getbase *)
********************************************************************************
}

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

(* begin module fixpiececoordinate *)
procedure fixpiececoordinate(var pie: pieceptr; excess: integer;
                        coordinateside: direction);
(* Fix the piece coordinates for insertions or deletions.  Coordinateside is
the end of the coordinate system that gets changed *)
begin
{
writeln(output,'fixpiececoordinate: excess = ',excess:1);
}
   with pie^.key do begin
      case coordinateside of
         minus: begin
            case piedir of
               minus: begin
                  (* piedir minus coordinateside minus *)
                  pieend := pieend - excess;
                  coobeg := coobeg - excess;
               end;
               plus: begin
                  (* piedir plus coordinateside minus *)
                  piebeg := piebeg - excess;
                  coobeg := coobeg - excess;
               end;
            end;
         end;
         plus: begin
            case piedir of
               minus: begin
                  (* piedir minus coordinateside plus *)
                  piebeg := piebeg + excess;
                  cooend := cooend + excess;
               end;
               plus: begin
                  (* piedir plus coordinateside plus *)
                  pieend := pieend + excess;
                  cooend := cooend + excess;
               end;
            end;
         end;
      end;
   end;
end;
(* end module fixpiececoordinate *)

(* begin module book.putbase *)
procedure putbase(b: base; position: integer; var pie: pieceptr;
                  coordinateside: direction);
(* put a base b into the nth position (internal coordinates) of the
   piece.  Protection is made against positions outside the piece.
NEW:
* If the base is before coordinate 1, the program halts.
* If the base is after the end of the sequence, extra space is made, and
  the coordinate system is changed on the coordinate side given.
*)
var
   excess: integer; (* the implied insertion size *)
   workdna: dnaptr; (* working dna segment *)
   p: integer; (* the last base of the dna part or current length *)
   pielength: integer; (* current length of the piece *)
{
z: integer; (* index to the piece for debuging stuff *)
}
begin
   pielength := piecelength(pie);
{
writeln(output,'putbase ==================================');
write  (output,'putbase: b = ',basetochar(b),', pielength=',pielength:1);
writeln(output,', position=',position:1);
}
   if position < 1 then begin
      writeln(output,'putbase: can not put bases',
                     ' before the start of the piece (position < 1)');
      writeln(output,'Program error!  Please report it to',
                     ' toms@ncifcrf.gov.');
      halt;
   end;
   workdna:=pie^.dna;
   if position > pielength then begin (* add to end of piece *)
{
writeln(output,'putbase: INCREMENT WORKDNA');
}
      (* since position > pielength, excess is always positive *)
      excess := position - pielength;
      fixpiececoordinate(pie,excess,coordinateside);
{
writeln(output,'putbase: excess=',excess:1);
}

      (* find the last segment *)
      p := workdna^.length;
      while workdna^.next <> nil do begin
{
writeln(output,'looping');
}
         workdna := workdna^.next;
         p := p + workdna^.length;
      end;
{
writeln(output,'USED SPACE is p=',p:1);
}

      if workdna^.length + excess <= dnamax then begin
         (* fill into the available segment *)
{
writeln(output,'NO NEW SEGMENTS NEEDED');
}
         workdna^.length := workdna^.length + excess;
         workdna^.part[workdna^.length] := b
      end

      else begin (* add new segments *)
{
write(output,'NEW SEGMENTS NEEDED:');
write(output,' p =',p:1);
write(output,' excess =',excess:1);
writeln(output);
}

            (* make the rest of the last segment useable by increasing
               its length up to dnamax, first increment p: *)
            p := p + (dnamax - workdna^.length);
            workdna^.length := dnamax;

         (* build enough segments to accommodate the new sequence *)
         while p < position do begin
{
writeln(output,' start adding segment, p = ',p:1,'--------------++++++++++');
}
            (* now make a new segment: *)
            p := p + dnamax;
            getdna(workdna^.next);
            workdna := workdna^.next;
            workdna^.length := dnamax; (* make it all available to fill *)
{
for z := 1 to dnamax do workdna^.part[z] := x;
writeln(output,'----------add segment, p = ',p:1);
}
         end;
         (* put the base in at the end *)
         workdna^.length := dnamax -(p-position);
         workdna^.part[workdna^.length] := b

{
;writeln(output,'placement, workdna^.length = ',workdna^.length:1);
showsegments(output,pie^.dna);
}

      end;
   end


{
;writeln(output,'putbase: final workdna^.length=',workdna^.length:1);
}
   else begin (* insert into the current piece *)
{
writeln(output,'putbase: NO EXCESS');
}
      (* locate segment in which to put the base *)
      p := workdna^.length;
      while position > p do begin
{
writeln(output,'THE PLACE =============!! p = ',p:1);
}
         if workdna^.next <> nil
         then workdna := workdna^.next;
         p := p + workdna^.length; (* only count the filled part *)
      end;
{
writeln(output,'FINAL     =============== p = ',p:1);
writeln(output,'FINAL     =============== position=',position:1);
writeln(output,'FINAL     ===============  (position > p)=',
                                           (position > p));
writeln(output,'FINAL     ===============  (workdna^.next <> nil)=',
                                           (workdna^.next <> nil));
;writeln(output,'putbase: about to smash: workdna^.length=',workdna^.length:1);
}
      workdna^.part[ workdna^.length - (p-position) ] := b;
{
;writeln(output,'putbase: after smash: workdna^.length=',workdna^.length:1);
}
{
writeln(output,'base ',b,' placed into ',
      workdna^.length - (p-position):1);
showdnasegment(output,workdna, workdna^.length);
writeln(output);
}
   end;

{
showsegments(output,pie^.dna);
;writeln(output,'putbase: p=',p:1);
;writeln(output,'putbase: dnamax=',dnamax:1);
;writeln(output,'putbase: position=',position:1);
;writeln(output,'==========================================');
}

end;
(* end module book.putbase *)

(* begin module book.stepbase *)
function stepbase(startdna: dnaptr; var dna: dnaptr; var d: dnarange): base;
(* advance d by one base in dna and then return the base at the new d.
   (this means that one should initialize d to zero)
   if we go past the last base, we restart at startdna.
   note: d is not the number of the base... it is used as a
   record for stepbase.  do not mess with it, and do not use it to find
   out what base you are on.  use a separate counter.  *)
begin
      if (d=dnamax) or (d=dna^.length) then begin
         d:=1;
         dna:=dna^.next;
         if dna=nil then dna:=startdna
      end
      else d:=succ(d);
      stepbase:=dna^.part[d]
end;
(* end module book.stepbase *)

(* procedures needed to read books *******************************************)

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

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

(* read book variables ****************************************************)
(* these procedure read attributes from the book.  they are all
   prefixed by b to indicate this.  *)

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

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

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

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

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

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

(* begin module book.brstate *)
procedure brstate(var thefile: text; var theline: integer; var sta: state);
(* read a state *)
var   ch: char;
begin
      skipstar(thefile);
      read(thefile,ch);
      readln(thefile,ch);
      theline := succ(theline);
      if ch='n' then sta:=on
                else sta:=off
end;
(* end module book.brstate *)

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


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

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

(* 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 = 5.62; (@ of search.p 1993 January 9 *)

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

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

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

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

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

(* begin module book.getpiece *)
procedure getpiece(var thefile: text; var theline: integer; var pie: pieceptr);
(* move to and read in the next piece in the book *)
var   ch: char;
begin
   ch:=getto(thefile,theline,['p']); (* get to the next p(iece) in the book *)
   if ch<>' ' then begin
      brpiece(thefile,theline,pie);
{
1999 june 2: removed this:
      ch:=getto(thefile,theline,['p']); (* read to end of p *)
}
{
bbb - now done in brpiece
      readln(thefile); (* read past piece *)
      theline := succ(theline);
}
   end
   else clearpiece(pie);
end;
(* end module book.getpiece *)

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

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

(* the procedures needed to do book writing *********************************)
(* original version 81/11/17 by gary stormo *)

(* begin module book.bwbasics *)

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

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

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

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

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

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

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

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

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

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

procedure bwheader(var book: text; hea: header);
(* write a key header to the book *)
begin
   with hea do begin
      bwname(book,keynam);
      bwline(book,fulnam);
      bwnote(book,note);
   end
end;
 
procedure bworgkey(var book: text; org: orgkey);
(* write the organism key *)
begin
   with org do begin
      bwheader(book,hea);
      bwline  (book,mapunit)
   end
end;
 
procedure bwchrkey(var book: text; chr: chrkey);
(* write the chromosome key *)
begin
   with chr do begin
      bwheader(book,hea);
      bwreanum(book,mapbeg);
      bwreanum(book,mapend)
   end
end;

(* end module book.bwbasics *)

(* begin module book.bworg *)
procedure bworg(var thefile: text; org: orgkey;
                var chropen, orgopen: boolean);
(* this writes the organism key 'org' to 'thefile', and returns 'orgopen'
   as true.  if there is already a chromosome or organism open they
   are closed. *)
begin
      if chropen then begin
         writeln(thefile,'chromosome');
         chropen := false;
      end;
      if orgopen then writeln(thefile,'organism');
      writeln(thefile,'organism');
      bworgkey(thefile,org);
      orgopen := true;
end;
(* end module book.bworg *)

(* begin module book.bwchr *)
procedure bwchr(var thefile: text; chr: chrkey; var chropen: boolean);
(* write the chromosome key 'chr' to 'thefile' and return chropen as
   true.  if a chromosome is already open it is closed. *)
begin
      if chropen then writeln(thefile,'chromosome');
      writeln(thefile,'chromosome');
      bwchrkey(thefile,chr);
      chropen := true;
end;
(* end module book.bwchr *)

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

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

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

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

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

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

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

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

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

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

      bwstartline(thefile);
      writeln(thefile,key.pieend:1);
   end;
   bwdna(thefile,pie^.dna);
   writeln(thefile,'piece');
end;
(* end module book.bwpie *)

(* begin module book.bwref *)
procedure bwref(var book: text; ref: reference);
(* write a key reference to the book *)
begin
   with ref do begin
      bwname  (book,pienam);
      bwreanum(book,mapbeg);
      bwdirect(book,refdir);
      bwnumber(book,refbeg);
      bwnumber(book,refend)
   end
end;
(* end module book.bwref *)

(* begin module book.bwgen *)
procedure bwgen(var thefile: text; gene: genkey);
(* this proecdure writes to 'thefile' the information in 'gene',
   properly formatted; *)
begin
      writeln(thefile,'gene');
      with gene do
      begin
         bwheader(thefile,hea);
         bwref(thefile,ref);
      end;
      writeln(thefile,'gene');
end;
(* end module book.bwgen *)

(* begin module book.bwtra *)
procedure bwtra(var thefile: text; trans: trakey);
(* this proecdure writes to 'thefile' the information in 'trans',
   properly formatted; *)
begin
      writeln(thefile,'transcript');
      with trans do
      begin
         bwheader(thefile,hea);
         bwref(thefile,ref);
      end;
      writeln(thefile,'transcript');
end;
(* end module book.bwtra *)

(* begin module book.bwmar *)
procedure bwmar(var thefile: text; mark: marker);
(* this proecdure writes to 'thefile' the information in 'mark',
   properly formatted; *)
var i : integer;
begin
      writeln(thefile,'marker');
      with mark.key do
      begin
         bwheader(thefile,hea);
         bwref(thefile,ref);
         bwstate(thefile,sta);
         bwstartline(thefile);
         for i := 1 to phenotype^.length do
            write(thefile,phenotype^.letters[i]);
         writeln(thefile);
      end;
      bwdna(thefile,mark.dna);
      writeln(thefile,'marker');
end;
(* end module book.bwmar *)


(* delila instruction writing routines ***************************************)
(* this set of routines allows one to write delila instructions.
   the types are found in book.type *)

(* begin module book.iwcombk *)
procedure iwcombk(var book, afile: text);
(* make a comment in the file that says the name of the book *)
begin
      write(afile,'(* ');
      reset(book);
      if copylines(book, afile, 1) = 0 then begin
         writeln(output,
            ' book is empty, can not write comment for instructions');
         halt
      end;
      writeln(afile,'*)');
end;
(* end module book.iwcombk *)

(* begin module book.iwname *)
procedure iwname(var thefile: text; thename: name);
(* write the name to the file *)
var   c: integer;
begin
      for c:=1 to thename.length
         do write(thefile, thename.letters[c])
end;
(* end module book.iwname *)

(* begin module book.iworg *)
procedure iworg(var afile: text; org: orgkey);
(* write an organism specification.  no writeln is done to allow
   write orgchr to do this. *)
begin
      write(afile,'organism ');
      iwname(afile,org.hea.keynam);
      write(afile,';');
end;
(* end module book.iworg *)

(* begin module book.iwchr *)
procedure iwchr(var afile: text; chr: chrkey);
(* write an chromosome specification.  no writeln is done to allow
   write orgchr to do this. *)
begin
      write(afile,'chromosome ');
      iwname(afile,chr.hea.keynam);
      write(afile,';');
end;
(* end module book.iwchr *)

(* begin module book.iwpie *)
procedure iwpie(var afile: text; pie: piekey);
(* write a piece specification *)
begin
      write(afile,'piece ');
      iwname(afile,pie.hea.keynam);
      writeln(afile,';');
end;
(* end module book.iwpie *)

(* begin module book.iworgchr *)
procedure iworgchr(var afile: text;
                       org: orgkey; orgchange, orgopen: boolean;
                       chr: chrkey; chrchange, chropen: boolean);
(* write both organism and chromosome specifications, based
on whether the organism or chromosome changed (orgchange and chrchange)
and whether they are currently open (orgopen, chropen).
See getocp in the br routines. *)
begin
   if orgchange and orgopen then iworg(afile,org);
   if orgchange and chrchange and orgopen and chropen then write(afile,' ');
   if chrchange and chropen then iwchr(afile,chr);
   if (orgchange and orgopen) or (chrchange and chropen) then writeln(afile)
end;
(* end module book.iworgchr *)

(* begin module book.iwget *)
procedure iwget(var afile: text; pie: pieceptr;
                    fromplace, pieceplace, toplace: integer;
                    flip: boolean; insttype: integer;
                    carriagereturn: boolean);
(* print a get delila instruction in the orientation of pie,
   from fromplace to toplace pieceplace. +/- direction;
If flip is false, the piece direction is as on the piece,
if it is true, the it is the opposite direction.
insttype: instruction type.
          insttype=1 means the form   get from p  -/+f to same +/-t dir +/-;
          insttype=2 means the form   get from p1 -/+f to p2   +/-t dir +/-;
          where p, p1 and p2 are locations
carriagereturn: if true, add a carriage return to the end of the line.
 *)
procedure iwposition(relative: integer; sameallowed: boolean);
(* write the <POSITION> *)
procedure iwrelative(relative: integer);
begin
           if relative>=0 then write(afile,' +', relative:1)
      else if relative< 0 then write(afile,' ',relative:1)
end;
begin (* iwposition *)
      if (insttype = 1) and sameallowed
      then write(afile,' same')
      else write(afile,' ',pieceplace:1);
      case pie^.key.piedir of
         plus:  iwrelative(+relative);
         minus: iwrelative(-relative);
      end
end;
begin (* iwget *)
      write(afile,'get from');
      iwposition(fromplace, false);
      write(afile,' to');
      iwposition(toplace, true);
      write(afile,' direction');
      case pie^.key.piedir of
                  dirhomologous, (* handle case, may not be right *)
         plus:  case flip of
                  false: write(afile,' +');
                   true: write(afile,' -');
                end;
                  dircomplement, (* handle case, may not be right *)
         minus: case flip of
                  false: write(afile,' -');
                   true: write(afile,' +');
                end
      end;         
      write(afile,';');
      if carriagereturn then writeln(afile);
end;
(* end module book.iwget *)

(* begin module book.iwget2 *)
procedure iwget2(var afile: text; pie: pieceptr;
                    fromplace, place1, toplace, place2: integer;
                    flip: boolean;
                    carriagereturn: boolean);
(* print a get Delila instruction in the orientation of pie,
The form of the instructions is:
  get from place1 +/-fromplace
        to place2 +/-toplace direction +/-;
If flip is false, the piece direction is as on the piece,
if it is true, the it is the opposite direction.
carriagereturn: if true, add a carriage return to the end of the line.
*)
procedure iwposition(place, relative: integer);
   procedure iwrelative(relative: integer);
      begin
              if relative>=0 then write(afile,' +', relative:1)
         else if relative< 0 then write(afile,' ',relative:1)
      end;
   begin (* iwposition *)
         write(afile,' ',place:1);
         case pie^.key.piedir of
            plus:  iwrelative(+relative);
            minus: iwrelative(-relative);
         end
   end;
begin (* iwget2 *)
      write(afile,'get from');
      iwposition(place1,fromplace);
      write(afile,' to');
      iwposition(place2,toplace);
      write(afile,' direction');
      case pie^.key.piedir of
                  dirhomologous, (* handle case, may not be right *)
         plus:  case flip of
                  false: write(afile,' +');
                   true: write(afile,' -');
                end;
                  dircomplement, (* handle case, may not be right *)
         minus: case flip of
                  false: write(afile,' -');
                   true: write(afile,' +');
                end
      end;
      write(afile,';');
      if carriagereturn then writeln(afile);
end;
(* end module book.iwget2 *)

(* begin module book.iwgetsimple *)
procedure iwgetsimple(var afile: text; pie: pieceptr;
                      fromplace, toplace: integer;
                      flip: boolean;
                      carriagereturn: boolean);
(* print a get delila instruction in the orientation of pie,
   from fromplace to toplace pieceplace. +/- direction;
If flip is false, the piece direction is as on the piece,
if it is true, the it is the opposite direction.
Format:
   get from [pieceplace+fromplace]
         to [pieceplace+  toplace] +/- direction;
where [] means to precompute the value.
carriagereturn: if true, add a carriage return to the end of the line.
 *)
begin (* iwget *)
      write(afile,'get from');
      write(afile,' ',fromplace:1);
      write(afile,' to');
      write(afile,' ',  toplace:1);
      write(afile,' direction');
      case pie^.key.piedir of
                  dirhomologous, (* handle case, may not be right *)
         plus:  case flip of
                  false: write(afile,' +');
                   true: write(afile,' -');
                end;
                  dircomplement, (* handle case, may not be right *)
         minus: case flip of
                  false: write(afile,' -');
                   true: write(afile,' +');
                end
      end;         
      write(afile,';');
      if carriagereturn then writeln(afile);
end;
(* end module book.iwgetsimple *)

(* modules for book alignment ************************************************)

(* begin module findblank *)
procedure findblank(var afile: text);
(* read a file to find the next blank character *)
var   ch: char;
begin
      repeat read(afile,ch) until ch = ' '
end;
(* end module findblank *)

(* begin module findnonblank *)
procedure findnonblank(var afile: text; var ch: char);
(* find the next non blank character in a file, return it in ch. *)
begin
      ch:=' ';
      while (not eof(afile)) and
            (ch = ' ')
      do begin
         read(afile,ch);
         if eoln(afile) then readln(afile)
      end
end;
(* end module findnonblank *)

(* begin module oldalign.align *)
procedure oldalign(var inst, book: text; var theline: integer;
                 var pie:  pieceptr;
                 var length, alignedbase: integer);
(* documentation on align is in module info.align and
delman.use.aligned.books *)
const maximumrange = 500; (* if the alignment point is more than this
   distance from the piece ends, the program halts in an attempt to catch
   the alignment bug... 1991 Jan 11 It appears that the rewrite of the
   code has removed the bug, but the check will be kept. *)
var
   ch: char; (* a character in inst *)
   comment: boolean; (* true means we are inside a comment *)
   done: boolean; (* done finding an aligning get *)
   thebase: integer; (* the base read in *)
begin
   if not eof(book) then begin (* if there is still more to the book ... *)
      getpiece(book,theline,pie); (* read in the piece *)
      if not eof(book) then begin (* if we found a piece ... *)
         length:=pietoint(pie^.key.pieend,pie); (* calculate piece length *)

         (* now find inst the next occurance of 'get' *)
         done := false;
         while not done do begin
            if eof(inst) then begin (* no instructions? *)
               alignedbase := 1; (* simply align by the first base *)
               done := true
            end
            else begin
               if inst^ = '(' then begin (* skip comment *)
                  get(inst);
                  if inst^ = '*' then comment := true;
    {if comment then write(output,'COMMENT: (');}
                  while comment do begin
                     if eof(inst) then begin
                        writeln(output,' in procedure align:');
                        writeln(output,' an instruction comment does not end!');
                        halt
                     end;
    {write(output,inst^);}
                     get(inst);
                     if inst^ = '*' then begin
                        get(inst);
    {if inst^ = ')' then writeln(output,'*)');}
                        if inst^ = ')' then comment := false
                     end;
                  end;
               end;
               if inst^ = 'g' then begin
                  get(inst);
                  if inst^ = 'e' then begin
                  get(inst);
                     if inst^ = 't' then begin
                        get(inst);
                        if inst^ = ' ' then begin
                           findnonblank(inst,ch); (* get to "from" *)
                           findblank(inst); (* get past "from" *)
                           read(inst,thebase); (* read in the alignedbase *)
    {writeln(output,'thebase=',thebase:1);}
                           alignedbase:=pietoint(thebase,pie);
                           done := true
                        end;
                     end;
                  end;
               end;
               get(inst); (* move along now *)
            end;
         end;

         if (alignedbase <= -maximumrange) or
            (alignedbase > length + maximumrange) then begin
            writeln(output,' in procedure align:');
            writeln(output,' read in base was ',thebase:1);
            writeln(output,' in internal coordinates: ',alignedbase:1);
            writeln(output,' maximum range was ',maximumrange:1);
            writeln(output,' piece length was ',length:1);
            with pie^.key.hea.keynam
               do writeln(output,' piece name: ',letters:length);
            writeln(output,' piece number: ',number:1);
            writeln(output,' aligned base is too far away... see the code');
            halt
         end
      end
   end
end;
(* end module oldalign.align *)

(* begin module interact.clearstring *)
(* begin module clearstring *)
procedure clearstring(var ribbon: string);
(* empty the string *)
var   index: integer; (* to the ribbon *)
begin (* clearstring *)
      with ribbon do begin
         for index := 1 to maxstring do letters[index] := ' ';
         length := 0;
         current := 0;
      end
end; (* clearstring *)
procedure initializestring(var ribbon: string);
(* start the string with a nil pointer.   This routine
should be called before doing linked list work.  This allows
the standard string routines to clear the string without
killing the pointer. *)
begin (* initializestring *)
   clearstring(ribbon);
   ribbon.next := nil;
end; (* initializestring *)
(* end module clearstring version = 4.86; (@ of prgmod.p 2004 Sep 8 *)
(* end module interact.clearstring version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* begin module filler.fillstring *)
procedure fillstring(var s: string; a: filler);
(* this procedure makes it reasonably easy to fill the string s with
characters.  one calls the procedure as: *)
(*                           1         2         3         4         5 *)
(*                  12345678901234567890123456789012345678901234567890 *)
(*   fillstring(s, 'this-is-the-string                                ');
the two comments make it easy to line the characters up. also, for this
example, it was assumed that the length of filler as defined by the
constant fillermax was 50. *)
var
      length: integer; (* of the string without trailing blanks *)
      index: integer; (* of s *)
begin (* fillstring *)
      clearstring(s);
      length := fillermax;
      while (length > 1) and (a[length] = ' ') do length := pred(length);
      if (length = 1) and (a[length] = ' ') then begin
         writeln(output, 'fillstring: the string is empty');
         halt
      end;

      for index := 1 to length do s.letters[index] := a[index];
      s.length := length;
      s.current := 1
end; (* fillstring *)
(* end module filler.fillstring version = 4.86; (@ of prgmod.p 2004 Sep 8 *)
 
(* begin module filler.filltrigger *)
procedure filltrigger(var t: trigger;
                      a: filler);
(* fill the trigger t *)
begin (* filltrigger *)
         fillstring(t.seek,a)
end; (* fillstring *)
(* end module filler.filltrigger version = 4.86; (@ of prgmod.p 2004 Sep 8 *)

(* begin module trigger.proc *)
(* this module allows one to scan a series of characters, as from
an array or a file, and to "trigger" or detect a simple string
in the series.  the advantage of the trigger is that several triggers
can "observe" a stream of characters at once, each looking for a
different thing.
some other modules required: interact.const, interact.type *)

procedure resettrigger(var t: trigger);
(* reset the trigger to ground state *)
begin (* resettrigger *)
      with t do begin
         state := 0;
         skip := false;
         found := false
      end
end; (* resettrigger *)

procedure testfortrigger(ch: char; var t: trigger);
(* look at the character ch.
   if it is part of the trigger (at the current trigger state),
       then the trigger state goes higher.
   if it is not part of the trigger then the trigger state is reset,
      skip is true and one should skip onward to find the trigger.
   if the trigger is found, found is true.
1996 Sep 12: Bug found!  In the case of a trigger "ab", the program
used to miss it for situations like "aab".  This was because at the
first a it would step up.  Then it would see the second a and recognize
that was not part of ab.  It would fail to realize that it could be
the start of a new one.  The code now accounts for that possibility. *)
begin (* testfortrigger *)
      with t do begin
         state := succ(state);
{
writestring(list,seek);
writeln(list,'testfortrigger seek.letters[',state:1,']:',
seek.letters[state],' ch:',ch);
}
         if seek.letters[state] = ch
         then begin
            skip := false;
            if state = seek.length then found := true
                                   else found := false
         end
         else begin
            (* it failed.  But wait!  It could be the beginning
            of a NEW trigger string! *)
            if seek.letters[1] = ch then begin
               state := 1;
               skip := false;
               found := false
            end
            else begin (* reset trigger *)
               state := 0;
               skip := true;
               found := false
            end
         end
      end
end; (* testfortrigger *)
(* end module trigger.proc version = 4.86; (@ of prgmod.p 2004 Sep 8 *)


(* begin module align.align *)
procedure align(var inst, book: text; var theline: integer;
                var pie:  pieceptr;
                var length, alignedbase: integer);
(* documentation on align is in module info.align and
delman.use.aligned.books.
1996 Sep 12: The routine now uses the trigger functions found
in prgmod.  The bug in the oldalign routine (that it misses the
end of comments that end in a series of asterisks) has been fixed.
It now checks that the piece corresponds to the book. *)
const maximumrange = 10000; (* if the alignment point is more than this
   distance from the piece ends, the program halts in an attempt to catch
   the alignment bug... 1991 Jan 11 It appears that the rewrite of the
   code has removed the bug, but the check will be kept. *)
   semicolon = ';'; (* end of delila instruction *)
var
   ch: char; (* a character in inst *)
   p: integer; (* index to a piece name *)
   p1: integer; (* another index to a piece name *)
   done: boolean; (* done finding an aligning get *)
   thebase: integer; (* the base read in *)
   indefault: boolean; (* true when within a default statement.  These
      can contain the word 'piece', which must be ignored. *)

   gettrigger: trigger; (* trigger to find 'get' *)
   defaulttrigger: trigger; (* trigger to find 'default' *)
   nametrigger: trigger; (* trigger to find 'name' *)
   piecetrigger: trigger; (* trigger to find 'piece' *)
   settrigger: trigger; (* trigger to find 'set' *)

   begincomment: trigger; (* trigger to find '(-*' (ignore the dash!) *)
   endcomment: trigger; (* trigger to find '*-)'  (ignore the dash!) *)
   begincurly: trigger; (* trigger to find comments: '{' *)
   endcurly: trigger; (* trigger to find comments: '}' *)

   quote1trigger: trigger; (* trigger to find single quote ' *)
   quote2trigger: trigger; (* trigger to find double quote " *)

   dotteddone: boolean; (* a dot '.' has been found in the name - ignore
      the rest of the name - for comparisons with mutations. *)

{
procedure rd(var f: text; var ch: char);
(* read ch from f allowing inspection of the result *)
begin
   read(f,ch);
   write(output,ch);
   write(list,ch);
   write(output,'<',ch,'>');
end;

procedure rdln(var f: text);
(* readln f allowing inspection of the result *)
begin
   readln(f);
   writeln(output);
   writeln(list);
end;
}

procedure skipcomment(var f: text);
(* skip an entire comment *)
var
   comment: boolean; (* true means we are inside a comment *)
begin
   (* skip to end of comment *)
   resettrigger(endcomment);
   comment := true;
   while comment do begin
      if eof(f) then begin
         writeln(output,'A comment does not end!');
         halt
      end;
      if eoln(f) 
      then readln(f)
{ rdln(f) }
      else begin
{write(output,'<');
rd(f,ch);
write(output,'>');}
         read(f,ch);
         testfortrigger(ch, endcomment);
         if endcomment.found
         then comment := false;
      end
   end
end;

procedure skipcurly(var f: text);
(* skip an entire comment made by {}*)
var
   comment: boolean; (* true means we are inside a comment *)
begin
   (* skip to end of comment *)
   resettrigger(endcurly);
   comment := true;
   while comment do begin
      if eof(f) then begin
         writeln(output,'A comment does not end!');
         halt
      end;
      if eoln(f) 
      then readln(f)
{ rdln(f) }
      else begin
{write(output,'<');
rd(f,ch);
write(output,'>');}
         read(f,ch);
         testfortrigger(ch, endcurly);
         if endcurly.found
         then comment := false;
      end
   end
end;

procedure skipquote(quote: trigger);
(* skip an entire quote of either the ' or " persuasion *)
var
   kind: char; (* the kind of quote, ' or " *)
begin
   kind := quote.seek.letters[1];
{writeln(output,'skipquote ',kind);}
   repeat
      findnonblank(inst,ch); (* get to the quote *)
   until (ch = kind) or eof(inst);
   if ch <> kind then begin
      writeln(output,'end of quote starting with ',kind,' not found');
      halt;
   end;
end;

begin
 filltrigger(defaulttrigger,'default');
     filltrigger(gettrigger,'get    ');
    filltrigger(nametrigger,'name   ');
   filltrigger(piecetrigger,'piece  ');
     filltrigger(settrigger,'set    ');
   filltrigger(begincomment,'(*     ');
     filltrigger(endcomment,'*)     ');
     filltrigger(begincurly,'{      ');
       filltrigger(endcurly,'}      ');
  filltrigger(quote1trigger,'''      ');
  filltrigger(quote2trigger,'"      ');

  resettrigger(defaulttrigger);
  resettrigger(gettrigger);
  resettrigger(nametrigger);
  resettrigger(piecetrigger);
  resettrigger(settrigger);

  resettrigger(begincomment);
  resettrigger(begincurly);
  resettrigger(quote1trigger);
  resettrigger(quote2trigger);

   indefault := false;

   if not eof(book) then begin (* if there is still more to the book ... *)
      getpiece(book,theline,pie); (* read in the piece *)
      if not eof(book) then begin (* if we found a piece ... *)
         length:=pietoint(pie^.key.pieend,pie); (* calculate piece length *)

         (* now find in inst the next occurance of 'get' *)
         done := false;
         while not done do begin
            if eof(inst) then begin (* no instructions? *)
               alignedbase := 1; (* simply align by the first base *)
               done := true
            end
            else begin
               if eoln(inst)
               then readln(inst)
{then rdln(inst)}
               else begin
{rd(inst,ch);}
                  read(inst,ch);

                  testfortrigger(ch, begincomment);
                  testfortrigger(ch, begincurly);
                  if begincomment.found or begincurly.found
                  then begin
                     if ch = '*' then begin
                        skipcomment(inst);
                        resettrigger(begincomment);
                     end
                     else begin
                        resettrigger(begincurly);
                        skipcurly(inst);
                     end
                  end
                  else begin (* we are not inside a comment *)

                     testfortrigger(ch, gettrigger);
                     if gettrigger.found then begin
                        findnonblank(inst,ch); (* get to "from" *)
                        findblank(inst); (* get past "from" *)
                        read(inst,thebase); (* read in the alignedbase *)
{writeln(output);writeln(output,'thebase = ',thebase:1);}
                        alignedbase:=pietoint(thebase,pie);
{writeln(output,'alignedbase=',alignedbase:1);}
                        done := true
                     end;

                     testfortrigger(ch, quote1trigger);
                     if quote1trigger.found then begin
                        skipquote(quote1trigger);
                     end;

                     testfortrigger(ch, quote2trigger);
                     if quote2trigger.found then begin
                        skipquote(quote2trigger);
                     end;

                     testfortrigger(ch, defaulttrigger);
                     if defaulttrigger.found then begin
                        indefault := true;
                        resettrigger(defaulttrigger)
                     end;
                     if ch = semicolon then indefault := false;

                     testfortrigger(ch, settrigger);
                     if settrigger.found then begin
                        indefault := true;
                        resettrigger(settrigger)
                     end;
                     if ch = semicolon then indefault := false;

                     (* check that piece names are correct *)
                     testfortrigger(ch, piecetrigger);
                     if not indefault then if piecetrigger.found then begin
                        skipblanks(inst); (* get to name *)
                        with pie^.key.hea.keynam do begin
{
                           for p := 1 to length do begin
}
                           (* 2007 Jun 22: replace loop with while
                              so that we can drop out when dotted
                              names are detected. *)
                           p := 1;
                           dotteddone := false;
                           while not dotteddone do begin
                              if eoln(inst)
                              then dotteddone := true
                              else begin
                                 read(inst,ch);
                                 (* ignore names after a dot *)
{
if ch = '.' then writeln(output,'inst dotteddone');
}
                                 if ch         = '.' then dotteddone := true;
                                 if letters[p] = '.' then dotteddone := true;
{
if ch = '.' then writeln(output,'book dotteddone');
writeln(output,'BUBBa ch = ',ch,'     ',p:1);
}
{zzz}
                                 if (letters[p] <> ch) and (not dotteddone) and
                                    (ch <> ';')
                                 then begin
                                    writeln(output,
                                           'The piece name in the book: ');
                                    writeln(output,letters:length);
                                    writeln(output,'does not match',
                                                 ' the inst file piece name:');
                                    (* write the letters that matched: *)
                                    for p1 := 1 to p-1 
                                        do write(output,letters[p1]);
                                    (* write the offending letter: *)
                                    write(output, ch);
                                    (* get the rest of the name and show it: *)
                                    done := eoln(inst);
                                    while not done do begin
                                      done := eoln(inst);
                                      if not done then begin
                                        read(inst,ch); 
                                        if (ch = ' ') or (ch = ';') 
                                        then done := true;
                                        if not done then write(output,ch);
                                      end; 
                                    end;
                                    writeln(output);
                                    (* mark the first letter
                                       that does not match: *)
                                    for p1 := 1 to p-1 do write(output,' ');
                                    write(output,'^'); 
                                    writeln(output);
                                    halt
                                 end;
                                 p := p + 1;
                                 if p > length then dotteddone := true;
                              end;
                           end
                        end;
                     end;

                  end
               end
            end
         end;

         if (alignedbase <= -maximumrange) or
            (alignedbase > length + maximumrange) then begin
            writeln(output,' In procedure align:');
            writeln(output,' read in base was ',thebase:1);
            writeln(output,' in internal coordinates: ',alignedbase:1);
            writeln(output,' maximum range was ',maximumrange:1);
            writeln(output,' piece length was ',length:1);
            with pie^.key.hea.keynam
               do writeln(output,' piece name: ',letters:length);
            writeln(output,' piece number: ',number:1);
            writeln(output,' aligned base is too far away... see the code');
            halt
         end
      end
   end
end;
(* end module align.align *)

(* begin module align.maxminalignment *)
procedure maxminalignment(var inst, book: text; var theline: integer;
                          var fromparam, toparam: integer;
                          alignmenttype: char);
(* prescan the book to find the range over which the pieces of the
   book are spread, relative to the aligned base.  the procedure uses
   the same variables that align does (so it can call align itself), and
   it returns the range in fromparam and toparam.
   alignmenttype: 'f' means alignment by First internal coordinate
      base, 'b' means alignment by Book, 'i' means alignment by Instructions.
 *)
const
      maximumrange = 500; (* the maximum size aligned piece;
                             this will presumably catch the alignment bug *)
var
      distance: integer; (* a distance to the aligned base *)
      pie:  pieceptr;
      length, alignedbase: integer;
begin
      new(pie);
      (* set an initial range for the two bounds *)
      fromparam:=+maxint;
        toparam:=-maxint;

      reset(book);
      reset(inst);
      while not eof(book) do begin
         case alignmenttype of
            'i': align(inst,book,theline,pie,length,alignedbase);
            'b','f':  begin
                 getpiece(book,theline,pie); (* read in the piece *)
                 length := piecelength(pie);
            end;
         end;

         if not eof(book) then begin
            case alignmenttype of
               'f': begin (* force alignment on first base *)
                    alignedbase := 0;
                    fromparam := 1;
                    distance:=length-alignedbase;
                    if toparam < distance then toparam:=distance;
               end;
               'i': begin (* use the alignedbase from the book *)
                    distance:=1-alignedbase;
                    if fromparam > distance then fromparam:=distance;

                    distance:=length-alignedbase;
                    if toparam < distance then toparam:=distance;
               end;
               'b': begin (* use the internal book *)
                    alignedbase := pietoint(0, pie);
                    distance:=1-alignedbase;
                    if fromparam > distance then fromparam:=distance;

                    distance:=length-alignedbase;
                    if toparam < distance then toparam:=distance;
               end;
            end;

            clearpiece(pie)
         end
      end;

      if toparam - fromparam > maximumrange then begin
         writeln(output,' in procedure maxminalignment:');
         writeln(output,' alignedbase = ',alignedbase:1);
         writeln(output,' fromparameter = ',fromparam:1);
         writeln(output,' toparameter = ',toparam:1);
         writeln(output,' this exceeds the maximum range allowed (',
                        maximumrange:1,')');
         writeln(output,' see notes in the procedure. ');
         halt
         (* notes:  if you desired this range, increase 'maximumrange'.
            otherwise, this may indicate a bug - either:
            1) locate the bug (and tell tom schneider, please...)
            2) reduce the size of the fragments, from one or the other end
               until the bombing is stopped. *)
      end;

      (* make the book readable again *)
      reset(book);
      reset(inst);
      dispose(pie)
end;
(* end module align.maxminalignment *)

(* begin module align.withinalignment *)
function withinalignment(alignedposition, alignedbase, length: integer)
                         :boolean;
(* this function tells one if an aligned position, relative to an aligned
   base in a piece of some length is within the piece. *)
var   p: integer;(* the position on the piece *)
begin
      p := alignedposition + alignedbase;
      withinalignment := (p>0) and (p<=length)
end;
(* end module align.withinalignment *)

(* amino acid procedures  ****************************************************)

(* this is a set of procedures that deal with amino acids.
   they require the book reading routines. *)

(* begin module amino.getaminoacid *)
procedure getaminoacid(position: integer; pie: pieceptr;
                       var aa: aminoacid);
(* get the amino acid corresponding to a codon starting in pie at position.
   the amino acid is blank ('   ') if it is outside a linear piece.
   amb = amber nonsense (uag)
   och = ochre nonsense (uaa)
   opa = opal  nonsense (uga)
*)
var
      b0, b1, b2: base; (* bases that correspond to the amino acid *)
      length: integer; (* length of the piece *)
      done: boolean; (* flag for out of bounds test *)
begin
      length:=pietoint(pie^.key.pieend,pie);
      done:=false;

      (* deal with position out of the piece *)
      if (position < 1) or (position > length - 2)
      then case pie^.key.piecon of
         linear: begin (* ignore out of bounds requests *)
            aa:='   ';
            done:=true
         end;
         circular: begin (* modify the local copy of position
                            to be within the circle. *)
            while position < 1 do position:=position + length;
            while position > length do position:=position - length;

            b0:=getbase(position,pie);

            (* step through the codon, making sure to stay on the circle. *)
            position:=succ(position);
            if position > length then position:=1;
            b1:=getbase(position,pie);

            position:=succ(position);
            if position > length then position:=1;
            b2:=getbase(position,pie);
         end
      end
      else begin (* pick up the bases at the position *)
         b0:=getbase(position+0,pie);
         b1:=getbase(position+1,pie);
         b2:=getbase(position+2,pie);
      end;

      if not done then case b0 of
         t: case b1 of
            t: case b2 of
               t: aa:='phe';
               c: aa:='phe';
               a: aa:='leu';
               g: aa:='leu';
            end;
            c: case b2 of
               t: aa:='ser';
               c: aa:='ser';
               a: aa:='ser';
               g: aa:='ser';
            end;
            a: case b2 of
               t: aa:='tyr';
               c: aa:='tyr';
               a: aa:='och';
               g: aa:='amb';
            end;
            g: case b2 of
               t: aa:='cys';
               c: aa:='cys';
               a: aa:='opa';
               g: aa:='trp';
            end;
         end;
         c: case b1 of
            t: case b2 of
               t: aa:='leu';
               c: aa:='leu';
               a: aa:='leu';
               g: aa:='leu';
            end;
            c: case b2 of
               t: aa:='pro';
               c: aa:='pro';
               a: aa:='pro';
               g: aa:='pro';
            end;
            a: case b2 of
               t: aa:='his';
               c: aa:='his';
               a: aa:='gln';
               g: aa:='gln';
            end;
            g: case b2 of
               t: aa:='arg';
               c: aa:='arg';
               a: aa:='arg';
               g: aa:='arg';
            end;
         end;
         a: case b1 of
            t: case b2 of
               t: aa:='ile';
               c: aa:='ile';
               a: aa:='ile';
               g: aa:='met';
            end;
            c: case b2 of
               t: aa:='thr';
               c: aa:='thr';
               a: aa:='thr';
               g: aa:='thr';
            end;
            a: case b2 of
               t: aa:='asn';
               c: aa:='asn';
               a: aa:='lys';
               g: aa:='lys';
            end;
            g: case b2 of
               t: aa:='ser';
               c: aa:='ser';
               a: aa:='arg';
               g: aa:='arg';
            end;
         end;
         g: case b1 of
            t: case b2 of
               t: aa:='val';
               c: aa:='val';
               a: aa:='val';
               g: aa:='val';
            end;
            c: case b2 of
               t: aa:='ala';
               c: aa:='ala';
               a: aa:='ala';
               g: aa:='ala';
            end;
            a: case b2 of
               t: aa:='asp';
               c: aa:='asp';
               a: aa:='glu';
               g: aa:='glu';
            end;
            g: case b2 of
               t: aa:='gly';
               c: aa:='gly';
               a: aa:='gly';
               g: aa:='gly';
            end;
         end;
      end;
end;
(* end module amino.getaminoacid *)

(* begin module basepair *)
function basepair(xseq,yseq: pieceptr; x,y: integer; guallowed: boolean)
                 :boolean;
(* does xseq basepair to yseq at x to y?
allow gu pair if guallowed is true.
a series of if-then statements are used for speed. *)
var
      bx,by: base; (* bases in the sequences *)
begin
      bx := getbase(x,xseq);
      by := getbase(y,yseq);

      if bx = complement(by)
      then basepair := true
      else if guallowed
         then if (bx=g) and (by=t)
            then basepair := true
            else if (bx=t) and (by=g)
               then basepair := true
               else basepair := false
         else basepair := false
end;
(* end module basepair *)

(* begin module freeenergy *)

(* this module contains:
      moveendsin: a procedure to find the core of a helix, where
         the core is the part that contributes the energy.
         gu pairs and single base pairs are removed from both ends.
      freeenergy: a function that uses moveendsin and then calculates the free
         energy of the core in kcal.
*)

procedure moveendsin(var x5bound, x3bound: integer; xpiece: pieceptr;
                     var y5bound, y3bound: integer; ypiece: pieceptr);
(* move the ends of the helix in.  gu pairs contribute no energy,
nor do single base pairs surrounded by gu pairs.  these are removed. *)
begin (* moveendsin *)
      (* move the bounds inward to the first occurance of base pairing.
         since a single base pair provides no energy, the successive position
         is also required to pair.  the energy beyond these bounds is zero. *)
      (* move one end in, decrease x3bound, increase y5bound
         note 1: x3bound cannot go below succ(x5bound) because the call to
         getbase would object to looking before the x5" end.  the x5 bound
         is brought up in the second while loop. *)

      while ((getbase(x3bound,xpiece) <> complement(getbase(y5bound,ypiece)))
             or
             (getbase(pred(x3bound),xpiece) <>
              complement(getbase(succ(y5bound),ypiece)))
            ) and (x3bound > succ(x5bound)) do begin
               x3bound := pred(x3bound);
               y5bound := succ(y5bound)
      end;

      (* move the other end in, increase x5bound, decrease y3bound
         note 2: for efficiency,x5bound can stop incrementation
         when it passes x3bound *)
      while ((getbase(x5bound,xpiece) <> complement(getbase(y3bound,ypiece)))
             or
             (getbase(succ(x5bound),xpiece) <>
              complement(getbase(pred(y3bound),ypiece)))
            ) and (x5bound < x3bound) do begin
               x5bound := succ(x5bound);
               y3bound := pred(y3bound)
      end
end; (* moveendsin *)

function freeenergy(x5start, x3start: integer; xpiece: pieceptr;
                    y5start, y3start: integer; ypiece: pieceptr): real;
(* evaluate the free energy of the pairing
      (x3start, y5start) to (x5start, y3start).
interior loops are evaluated only for the case where both sides are
equal in length.  the units of the energy returned are kcalorie.
based on tinoco et al nature new biology vol 246 pp 40-41, 1973.
      internal coordinates (1 to n) are used.  gu pairs are eliminated
from the ends by calling moveendsin.
*)
var
      x5, x3, y5, y3: integer; (* 5 and 3 prime ends on x and y pieces *)
      x5b, x3b, y5b, y3b: base; (* the bases corresponding to the above *)
      x5bound, x3bound, y5bound, y3bound: integer; (* the bounds of energy
                                                      calculation *)
      energy: real; (* the current total energy *)
      bulgebases: integer; (* the number of bases involved in a bulge *)
      bulging: boolean; (* flag for continuing a bulge after a single match *)
procedure readofftable;
(* create a 4 by 4 table based on table 1 in tinoco et al.
we know at this point that both base pairs complement, so only one
pair is used. *)
var
      deltag: real; (* temporary curve *)
begin
      case x5b of
         a: case x3b of
            a: deltag := -1.2;
            c: deltag := -2.2;
            g: deltag := -2.2;
            t: deltag := -1.8
         end;
         c: case x3b of
            a: deltag := -2.2;
            c: deltag := -5.0;
            g: deltag := -3.2;
            t: deltag := -2.2
         end;
         g: case x3b of
            a: deltag := -2.2;
            c: deltag := -5.0;
            g: deltag := -5.0;
            t: deltag := -2.2
         end;
         t: case x3b of
            a: deltag := -1.8;
            c: deltag := -2.2;
            g: deltag := -2.2;
            t: deltag := -1.2
         end;
      end;
      energy := energy + deltag
end;

procedure gupair;
(* in the tinoco paper, one pairing has exceptional energy: gu/ug *)
begin
      energy := energy - 0.3
end;

procedure continuebulge;
(* the interior loop bulge continues *)
begin
      bulgebases := bulgebases + 2
end;

procedure beginbulge;
(* begin a bulge, of the interior loop kind *)
begin
      if not bulging then begin
         bulgebases := 2;
         bulging := true
      end
      else continuebulge (* only a single base matches,
                            so it is to be ignored *)
end;

procedure endbulge;
(* close the bulge by recording its energy *)
begin
      (* perhaps it is only a single match *)
      if getbase(pred(x5),xpiece) = complement(getbase(succ(y3),ypiece))
      then begin (* it is really the end of a bulge *)
              if (2 <= bulgebases) and (bulgebases <= 6)
              then energy := energy + 2.0
         else if (7 <= bulgebases) and (bulgebases <= 20)
              then energy := energy + 3.0
         else (* bulgebases > 20 *)
              energy := energy + 1.0  + 2.0*ln(bulgebases)/ln(10);
         bulging := false
      end
      else continuebulge (* the single match is to be ignored *)
end;

begin (* freeenergy *)
      (* test consistency of input *)
      if (x5start - x3start) <> (y5start - y3start) then begin
         writeln(output,' function freeenergy:');
         writeln(output,' the supplied ends are not consistent with a helix.');
         halt
      end;

      x5bound := x5start;
      x3bound := x3start;
      y5bound := y5start;
      y3bound := y3start;

      moveendsin(x5bound,x3bound,xpiece,
                 y5bound,y3bound,ypiece);

      (* set external energy to zero *)
      energy := 0.0;

      (* no bulges yet *)
      bulging := false;

      if x5bound < x3bound then begin
         (* within this if, the range is bounded by base pairs *)
         x3 := x3bound;
         y5 := y5bound;
         x5 := pred(x3);
         y3 := succ(y5);

         (* these points (x,y/5,3) now define a dinucleotide pair
         between x and y.  now we must scan the
         the dinucleotide across to the other bound. *)
         repeat
            (* convert positons to bases *)
            x5b := getbase(x5,xpiece);
            x3b := getbase(x3,xpiece);
            y5b := getbase(y5,ypiece);
            y3b := getbase(y3,ypiece);

                 if (x3b = complement(y5b)) and
                    (x5b = complement(y3b))
                 then readofftable         (* everybody pairs *)
            else if (y5b = g) and (y3b = t) and
                    (x3b = t) and (x5b = g)
                 then gupair               (* a special case *)
            else if (x5b <> complement(y3b)) and
                    (x3b <> complement(y5b))
                 then continuebulge        (* nobody pairs *)
            else if x3b = complement(y5b)
                 then beginbulge           (* the first two pair *)
            else if x5b = complement(y3b)
                 then endbulge             (* the last two pair *)
            else halt;

            x3 := pred(x3);
            y5 := succ(y5);
            x5 := pred(x5);
            y3 := succ(y3);
(*if debugging then writeln(output,energy:4:1); *)
         until x5 < x5bound
      end;
      freeenergy := energy
end; (* freeenergy *)
(* end module freeenergy version = 'delmod 6.17 84 apr 12 tds/gds'; *)

(* begin module lowestenergy *)
function lowestenergy (s: pieceptr): real;
(* find the lowest possible energy of pairing to sequence s
by calculating the energy for pairing of s to its complement *)
var
      cs: pieceptr; (* complement of s *)
      slength: integer; (* lengths of s and cs *)
      p: integer; (* position *)
begin
      slength:=piecelength(s);

      (* see if the assumption noted below is true: *)
      if slength > dnamax then begin
         writeln(output,' in function lowestenergy:',
                        ' length of s > dnamax (',
                        slength:1,' > ',dnamax:1,').');
         writeln(output,' shorten s or increase dnamax.');
         halt
      end;

      (* make complement of s *)
      new(cs);
      clearpiece(cs);
      cs^.key:=s^.key;
      new(cs^.dna); (* assumption: that this will hold the entire piece... *)
      cs^.dna^.next:=nil;
      for p:=1 to slength
         do cs^.dna^.part[slength-p+1]:=complement(getbase(p,s));

      (* find the energy of pairing *)
      lowestenergy:=freeenergy(1,slength,s,1,slength,cs);
(*      if debugging then writeln(output,' lowest energy:',
                        freeenergy(1,slength,s,1,slength,cs):10:3);*)

      (* clean up *)
      clearpiece(cs);
         dispose(cs)
end;
(* end module lowestenergy *)

(* begin module info.nextbase *)
(* nextbase: book reading routines that return one base at a time

      the nextbase routines allow a programmer to access a delila library
without having to know any details of the library interface routines.
the routines are most efficient for reading a book of sequences
sequentially, that is, one base at a time.  for random access, the
user should learn how to use the other routines in the book modules.

      your program must contain several things:
1) the book reading routines must be inserted into your program.
   they are the modules book.const, book.type, book.var, and package.nextbase.
2) you must declare the several variables used by the two routines of nextbase.
   fortunately, you need not know how they are used inside the routines.
3) you must use the routine initnextbase to initilize variables
   used by the routine nextbase.
4) once these things are set up, all you need to do is call function nextbase,
   which will returns the next base.
5) if you want characters returned call
         acharacter:= basetochar(nextbase(....));
   other useful functions are contained in the book reading package.

*)
(* end module info.nextbase *)

(* begin module nextbase *)
procedure initnextbase(var book: text; var pie: pieceptr;
                       var lastbase, endofbook: boolean;
                       var theline: integer);
(* initialize variables for function nextbase.
   book is the book to be read,
   pie is the piece,
   lastbase is the flag that is true when we are at the
   last base of a piece,
   and endofbook is true if we are at the end of the book (see nextbase).
   theline records the current line number in the book. *)
begin
      brinit(book, theline);

      new(pie);
      with pie^ do begin
         with key.hea do begin
            fulnam:=nil;
            note:=nil
         end;
         dna:=nil
      end;

      lastbase:=true; (* this will trigger reading of the next piece *)
      if eof(book) then endofbook:=true
                   else endofbook:=false
end;

function nextbase(var book: text;  (* the book being read *)
                  var theline: integer; (* current line number of the book *)
(* the next variables can be ignored *)
                  var pie: pieceptr; (* the piece *)
                  var dnalink: dnaptr; (* the current link we are on *)
                  var dnalinkspot: dnarange; (* the spot in the dnalink *)
(* these are useful to the general user: *)
                  var dnaspot, (* integer in 1 to length, which base this is *)
                      length: integer; (* length of this piece *)
                  var lastbase, (* true if the base was the last one on
                                   the piece. *)
                      endofbook: boolean) (* true when we have reached
                                             the end of the book *)
                  : base;
(* the user simply declares variables for
   book, pie, dnalink, dnalinkspot (of the appropriate type)
note: you can convert the valuespot to published coordinates by
      pubcoords:= inttopie(dnaspot,pie);
warning: if the end of the book has been reached, then endofbook
is true, but the value returned by the function has no meaning.
*)
begin
   if not endofbook then begin
      if not lastbase then begin
         nextbase:=stepbase(pie^.dna,dnalink,dnalinkspot);
         dnaspot:=succ(dnaspot);
         if dnaspot = length then lastbase:=true
      end
      else begin (* we are at the last base of the previous piece *)
         clearpiece(pie);
         getpiece(book,theline,pie);
         if not eof(book) then begin
            dnalink:=pie^.dna;
            dnalinkspot:=0;
            dnaspot:=0;
            length:=piecelength(pie);
            lastbase:=false;
            endofbook:=false;
            nextbase:=nextbase(book,theline,pie,dnalink,dnalinkspot,
                               dnaspot,length,lastbase,endofbook)
         end
         else begin
            endofbook:=true; (* we have reached the end of the book *)
            lastbase:=false; (* we are no longer at the last base *)
            nextbase:=a (* a fake value *)
         end
      end
   end
end;
(* end module nextbase *)

(* begin module info.book.br.routines *)
(*
      information about the book reading (br) routines

these types, constants, variables, procedures, and functions make it
      relatively easy to read the books as defined in:
      'organism and recognition class library definition:
       a dna sequence data base'

the routines use several global variables. (we plan to eliminate this.)
readnumber: when an item (eg a piece) is read its notes are investigated.
      if readnumber is false, then the notes are read in line by line.
      if readnumber is true, then the notes are not available, because
      they are investigated for a number of the item:

number: if there is a number in the notes, this variable contains the value,
      and the variable:

numbered: is true (otherwise it is false).

skipunnum: if set to false, the next item is read in.  if skipunnum is
      true and an item is not numbered then the read routines will
      continue looking for a numbered item (skip unnumbered items).
      in this way one can select and read only the numbered items.

other global variables are:
      freeline, freedna
since these are used by the br routines, the user should not
use these names.


procedures:

procedure getdatetime(var adatetime: datetimearray);
      get the date and time into a single array from the system clock

procedure readdatetime (var thefile: text; var adatetime: datetimearray);
      read the date and time from the file

procedure writedatetime(var thefile: text; adatetime: datetimearray);
      expand the date and time out and print in the file

procedure brinit(var book: text);
      initializes a book file.  it is required before any book reads.

procedure getpiece(var book: text; var theline: integer; var pie: pieceptr);
      returns a pointer to the next piece in the book,
      or eof(book) is true when there is none.
      theline tracks the line number.

procedure clearpiece(var pie: pieceptr);
      returns the dna of a piece to free storage.
      it should be used when a piece is no longer needed.

procedure getocp(var thefile: text;  var theline: integer;
                 var org: orgkey;  var orgchange: boolean;
                                   var orgopen:   boolean;
                 var chr: chrkey;  var chrchange: boolean;
                                   var chropen:   boolean;
                 var pie: pieceptr; var pieopen:   boolean);
      get the next piece and its organism and chromosome keys.
      orgchange and chrchange indicate whether or not a new organism
      or chromosome was found.
      if a piece is not found, then eof(book) will be true.
      orgopen, chropen and pieopen are used by getocp to tell when it has
      entered an organism, chromosome or piece.
      they should be set to false initially.
      there should be one triplet for each book read.
      theline tracks the line number.

functions:

function pietoint(p: integer; pie: pieceptr): integer;
      converts a number in piece coordinates to a convenient internal
      coordinate system:
      piebeg becomes 1 and
      pieend becomes the length of the piece.

function inttopie(i: integer; pie: pieceptr):integer;
      the inverse of pietoint.

function getbase(i: integer; pie: pieceptr):base;
      returns the base at a given internal coordinate.

function stepbase(startdna: dnaptr; var dna: dnaptr; var d: integer): base;
      advance d by one base in dna and then return the base at the new d.
      (this means that one should initialize d to zero)
      if we go past the last base, we restart at startdna.

function basetochar(ba:base):char;
      converts a type base to a character.

function chartobase(ch:char):base;
      the inverse of basetochar.

function complement(ba:base):base;
      returns the complement of a base.

function piecelength(pie: pieceptr): integer;
      return the length of the dna in pie

dnamax is a constant which determines the size of each dnastring.
for most uses, 3000 works fine.  under some circumstances, other values
can be used, for example when a program will usually use far less than
3000 bases.  to calculate the number of bases, one can start with how
many machine words one wants to use per dna string.  for 100 words on a
machine with 60 bits per word:
       (60 bits/word)*(1base/2bits)*100 words=3000 bases
*)

(* end module info.book.br.routines *)

(* begin module info.book.br.dependencies *)
(*

 dependencies between book reading (br) routines

      the first word is the name of a procedure.  the other words on the line
are the procedures, functions or (package)s that the procedure calls or uses.
this set is useful for designing packages of procedures.
(note:  this is a condensed set.  not all actual calls are noted, if one
of the procedures called uses a lower procedure.  most of these are in
module basis.  the purpose of this module is to aid in construction, not to
document the structure.)

 basis (const) (type) (var) halt
 getbase basis
 stepbase basis
 getto  basis
 skipstar basis
 brreanum skipstar
 brnumber skipstar
 brname skipstar
 brline skipstar
 brdirect skipstar
 brconfig skipstar
 brstate skipstar
 brnotenumber
 brnote brline
 brheader brname brline brnotenumber brnote
 brorgkey brheader brline
 brchrkey brheader brreanum
 brpiekey brheader brreanum brconfig brdirect brnumber
 brdna getto
 brpiece brpiekey brdna
 getpiece getto brpiece
 getocp getto brorgkey brchrkey brpiece
 brinit
 copyheader

*)
(* end module info.book.br.dependencies *)

(* begin module info.align *)
(* sequence alignment routines

   align reads a set of delila instructions (inst), and the book
   created by delila from those instructions (book).  it returns
   the next piece in the book (pie), the length of the piece (length)
   in bases, and a number somewhere between 1 and length which can be
   used to align all the pieces of the book (alignbedbase).
   if no piece is found, then eof(book) is true, otherwise the piece
   is returned.
   if the instruction file has fewer get instructions than pieces in
   the book, then align will halt.

   the instructions are assumed to have a rigid format:
1. each get starts on a new line.
2. a "g" as the first character of a line is always a get.
3. all gets are in the form:
   get from # ...
   (where # is an integer)
4. the # is picked up as the number for aligning, and is converted
   to internal coordinates (see br routines), in the range 1 to
   length, inclusive.


   there are two other procedures in this package.  one finds the
   minimum and maximum range of the pieces in the book relative to
   the aligned base.  the other tests a position in 'aligned space'
   for inclusion in a particular piece.

   these procedures make it possible to scan in the space of aligned
   pieces rather than in the space of each piece.  for example, it
   allows one to list aligned pieces by running a for loop:
      for index:=fromparam to toparam
         do if withinalignment(index,alignedbase,length)
            then write(afile,basetochar(getbase(index+alignedbase,thepiece)))
            else write(afile,' ')
*)
(* end module info.align *)

(* begin module demo.nextbase *)
procedure demonextbase(var book, fout: text);
(* demonstration of the use of the nextbase routines *)
var
   (* these variables are all defined in nextbase *)
   pie: pieceptr;
   dnalink: dnaptr;
   dnalinkspot: dnarange;
   dnaspot, length: integer;
   lastbase, endofbook: boolean;
   theline: integer;
   (* this variable is needed to catch the value of nextbase,
      so that we can check that we have not reached the end of the book. *)
   character: char;
begin (* demonextbase *)
   initnextbase(book,pie,lastbase,endofbook,theline);
   theline := 0;
   writeln(fout,' demonstration of package.nextbase');
   writeln(fout,' the bases of the book are given,',
                  ' followed by two numbers:');
   writeln(fout,' the first is the internal coordinate of the base.');
   writeln(fout,' the second is the published coordinate of the base.');
   while not endofbook do begin
      character:=basetochar(
         nextbase(book,theline,pie,dnalink,dnalinkspot,
                  dnaspot,length,lastbase,endofbook));
      if not endofbook
         then writeln(fout,' ',character,' ',dnaspot:3,
                           ' ',inttopie(dnaspot,pie):3);
      if lastbase
         then writeln(fout,' end of a piece ',length:1,' bp long',
                           ' at line ',(theline):1, ' of the book');
   end
end; (* demonextbase *)
(* end module demo.nextbase *)

(* begin module demo.time *)
procedure demotime(var fout: text);
(* write the time to file fout *)
var   dateandtime: datetimearray; (* the date and time *)
      seed: real; (* a seed for a random number generater, made
         from the date and time written backwards *)
begin (* demotime *)
   getdatetime(dateandtime);
   writeln(fout);
   writeln(fout,'The date and time is: ');
   writedatetime(fout, dateandtime);
   writeln(fout,' <- This should be the current time and date.');
   write  (fout,'1980/06/09 18:49:11');
   writeln(fout,' <- Times and dates should look like this');
   writeln(fout,'year mo da ho mi se <- with parts in these positions');
   writeln(fout);

   timeseed(seed);
   writeln(fout,'A timeseed is ',seed:16:14);
   writeln(fout,'Timeseeds can be used to start random number generators.')
end; (* demotime *)
(* end module demo.time version = 1.15; (@ of timegpc.p 2000 Oct 11 *)

(* begin module delilamodules.themain *)
procedure themain;
(* the main procedure of the program *)
begin
   writeln(output, ' delila modules ', version:4:2);

   reset(book);
   if eof(book) then begin
      writeln(output,' if you provide delmod with a book,',
      ' then the nextbase routines will be tested.');
   end
   else demonextbase(book, output);

   demotime(output);

end;
(* end module delilamodules.themain *)

begin (* delmod *)
  themain;
1: end. (* delmod *)