program matrix(xbook,ybook,hlist,mlist,matrixp,output);
(* matrix: dot matricies of two books
      by thomas schneider
      copyright (c) 1987

modules from: delman, delmods, prgmods, auxmods
*)

label 1; (* end of program *)

const
(* begin module version *)
      version = 3.28; (* of matrix 1987 feb 13 *)
      (* original version 1981 *)
(* end module version *)

(* begin module describe.matrix *)
(*
name
      matrix: dot matrices for helices between two books

synopsis
      matrix(xbook: in, ybook: in, hlist: in, mlist: out,
             matrixp: in, output: out)

files
      xbook: a book from the delila system

      ybook: a book from the delila system.  If you want to look for
         structures in one sequence, then use the program copy to make
         a copy of xbook in ybook.

      hlist: the helix listing for xbook and ybook made by program helix

      mlist: the matrices listed out.  Sequences from the x book are printed
         vertically, while those from the y book are horizontal.
         helices are printed as a set of numbers:
            1 means gt base pair
            2 means at base pair
            3 means gc base pair
         if mlist is wider than your printer, use the split program.

      matrixp: parameters to control the mlist
         If matrixp is empty, default values are used. otherwise, the first
         line contains one number.  If this number is a positive integer,
         it specifies the minimum length helix in base pairs from hlist to
         record in mlist;
         if this number is a negative real number, it specifies the maximum
         energy in kcal of the helixes written in mlist.

      output: messages to the user

description
      Matrix produces a dot matrix for the two books.  Only helices
      of some length (or longer) or of some maximum energy (or less) are
      printed.  The helices are made using program helix.

documentation
      delman.use.comparison
      J. V. Maizel, Jr. and R. P. Lenk PNAS 78: 7665-7609 (1981)

see also
      helix, dotmat, split, keymat

author
      Thomas D. Schneider

bugs
      none known

technical notest
      The constant maxarray defines the maximum area that the program
      can handle.

*)
(* end module describe.matrix *)

(* more constants: *)
      headerlines = 7; (* the number of lines in the header of hlist *)
      nodot = ':'; (* the character to print for background,
                      note: it must not blank, so that the user can
                      align the helixes by the numbers. *)
      pagelength = 48; (* lines printed of matrix on one page in mlist *)
      maxarray = 150000; (* size of the bit matrix *)

      prime = ''''; (* single quote *)

(* begin module book.const *)
(* constants needed for book manipulations *)

      dnamax = 3000;  (* length of dna arrays *)
      namelength = 20; (* maximum key name length *)
      linelength = 80; (* maximum line readable in book *)

(* end module book.const version = 'delmod 6.51 85 apr 17 tds/gds' *)

type
      bit = (zero,one);

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


      (* base types *)
      base = (a,c,g,t);
      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 = 'delmod 6.51 85 apr 17 tds/gds' *)

var
      xbook, (* one book *)
      ybook, (* the other book *)
      hlist, (* list of helixes *)
      mlist, (* the dot matricies *)
      matrixp: (* parameters *)
         text;

      xseq, yseq: pieceptr; (* the sequences from the books *)
      xnumber, ynumber: integer; (* the numbers of the sequences *)
      minlenmaxene: real; (* a parameter read from matrixp;
                          the minimum size helix to print
                          or the maximum energy of the helix *)

      (* variables for the bit array *)
(* the program simulates a rectangular bit array (1,1) to (xlength,ylength)
in the linear bitarray. *)
      bitarray: packed array[1..maxarray] of bit;
      index: integer; (* direct index into bitarray *)

      (* globals that determine the dimensions used to access bitarray *)
      xlength, ylength: integer; (* of the sequences *)

      length: integer; (* current helix length *)

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

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

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

(* ************************************************************************ *)
(* end module book.var version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module package.primitive *)
(* ************************************************************************ *)
(* begin module halt *)
procedure halt;
(* stop the program.  the procedure performs a goto to the end of the
   program.  you must have a label:
      label 1;
   declared, and also the end of the program must have this label:
      1: end.
   examples are in the module libraries.
   this is the only goto in the delila system. *)
begin
      writeln(output,' program halt.');
      goto 1
end;
(* end module halt version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module unlimitln *)
procedure unlimitln(var afile: text);
(* this procedure removes a stupid system dependent limit on the number of
lines that one can write to a file.  you may remove it from the code if
your system does not want or need this.  suggested method: place comments
around the contents of the procedure. *)
begin
      linelimit(afile, maxint); (* set 'infinite' lines allowed for afile *)
end;
(* end module unlimitln version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* 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 = 'delmod 6.51 85 apr 17 tds/gds' *)

(* 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 = 'delmod 6.51 85 apr 17 tds/gds' *)

(* ************************************************************************ *)
(* end module package.primitive version = 'delmod 6.51 85 apr 17 tds/gds' *)

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

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

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

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

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

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

procedure cleardna(var l: dnaptr);
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 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 *)
var   i: integer; (* an intermediate value *)
begin
      with pie^.key do begin
         case piedir of
            plus:  if p>=piebeg
                   then i:=p-piebeg+1
                   else i:=(p-coobeg)+(cooend-piebeg)+2;
            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 *)
var   p: integer; (* an intermediate value *)
begin
      with pie^.key do begin
         case piedir of
            plus: begin
               p:=piebeg+(i-1);
               if p>cooend
                  then if coocon=circular
                     then p:=p-(cooend-coobeg+1)
            end;
            minus: begin
               p:=piebeg-(i-1);
               if p<coobeg
                  then if coocon=circular
                     then p:=p+(cooend-coobeg+1)
            end
         end;
         inttopie:=p
      end
end;

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

(* end module book.basis version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.getto *)
function getto(var thefile: text; ch: chset): char;
(* search the file for a character in the first line which is a
   member of the set ch. *)
var achar: char;
begin
      achar:=' ';
      while (not(achar in ch)) and (not eof(thefile))
         do readln(thefile,achar);
      if (achar in ch) then getto:=achar
                       else getto:=' '
end;
(* end module book.getto version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.skipstar *)
procedure skipstar(var thefile:text);
(* skip start of line (or star = '*'). *)
begin (* skipstar *)
      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; (* skipstar *)
(* end module book.skipstar version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brreanum *)
procedure brreanum(var thefile: text; var reanum: real);
(* read a real number from the file *)
begin
      skipstar(thefile);
      readln(thefile,reanum);
end;
(* end module book.brreanum version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brnumber *)
procedure brnumber(var thefile: text; var num: integer);
(* read a number from the file *)
begin
      skipstar(thefile);
      readln(thefile,num)
end;
(* end module book.brnumber version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brname *)
procedure brname(var thefile: text; 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)
end; (* brname *)
(* end module book.brname version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brline *)
procedure brline(var thefile: text; var l: lineptr);
(* read a line from the file *)
var
      i,j: integer;
      acharacter: char;
begin
      skipstar(thefile);
      i:=0;
      while (not eoln(thefile)) do begin
         i:=succ(i);
         read(thefile,acharacter);
         l^.letters[i]:=acharacter
      end;
      if i<l^.length then for j:=i+1 to l^.length do l^.letters[j]:=' ';
      l^.length:=i;
      l^.next:=nil;
      readln(thefile)
end;
(* end module book.brline version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brdirect *)
procedure brdirect(var thefile: text; var direct: direction);
(* read a direction *)
var   ch: char;
begin
      skipstar(thefile);
      readln(thefile,ch);
      if ch='+' then direct:=plus
                else direct:=minus
end;
(* end module book.brdirect version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brconfig *)
procedure brconfig(var thefile: text; var config: configuration);
(* read a configuration *)
var   ch: char;
begin
      skipstar(thefile);
      readln(thefile,ch);
      if ch='l' then config:=linear
                else config:=circular
end;
(* end module book.brconfig version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brnotenumber *)
procedure brnotenumber(var thefile: text; 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);
         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)
            until thefile^ = 'n';
            readln(thefile)
         end
         else readln(thefile)
      end
end; (* brnotenumber *)
(* end module book.brnotenumber version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brnote *)
procedure brnote(var thefile: text; 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);
         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,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)
         end
         else readln(thefile)
      end
end; (* brnote *)
(* end module book.brnote version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brheader *)
procedure brheader(var thefile: text; var hea: header);
(* read the header of a key. *)
begin
      with hea do begin
         (* read key name *)
         brname(thefile,keynam);

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

         (* read note key *)
         if readnumber then brnotenumber(thefile,note)
                       else brnote(thefile,note)
      end
end;
(* end module book.brheader version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brpiekey *)
procedure brpiekey(var thefile: text; var pie: piekey);
(* read piece key *)
begin
      with pie do begin
         brheader(thefile,hea);
         brreanum(thefile,mapbeg);
         brconfig(thefile,coocon);
         brdirect(thefile,coodir);
         brnumber(thefile,coobeg);
         brnumber(thefile,cooend);
         brconfig(thefile,piecon);
         brdirect(thefile,piedir);
         brnumber(thefile,piebeg);
         brnumber(thefile,pieend);
      end
end;
(* end module book.brpiekey version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brdna *)
procedure brdna(var thefile: text; var dna: dnaptr);
(* read in dna from thefile *)
(* 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,['d']);
      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 *)
         read(thefile,ch); (* ch is either '*' or 'd' *)
      end;
      readln(thefile)
end;
(* end module book.brdna version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brpiece *)
procedure brpiece(var thefile: text; var pie: pieceptr);
(* read in a piece *)
begin
      brpiekey(thefile,pie^.key);
      if numbered or (not skipunnum)
         then brdna(thefile,pie^.dna)
end;
(* end module book.brpiece version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.brinit *)
procedure brinit(var book: text);
(* 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;
end; (* brinit *)
(* end module book.brinit version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* ************************************************************************ *)
(* end module package.brpiece version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.getpiece *)
procedure getpiece(var thefile: text; var pie: pieceptr);
(* move to and read in the next piece in the book *)
var   ch: char;
begin
      ch:=getto(thefile,['p']); (* get to the next p(iece) in the book *)
      if ch<>' ' then begin
         brpiece(thefile,pie);
         ch:=getto(thefile,['p']); (* read past closing p *)
      end
end;
(* end module book.getpiece version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* ************************************************************************ *)
(* end module package.getpiece version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module book.getbase *)
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;
(* end module book.getbase version = 'delmod 6.51 85 apr 17 tds/gds' *)

(* begin module package.numbar *)
(* ************************************************************************ *)
(* begin module numberdigit *)
function numberdigit(number, logplace:integer): char;
(* return the digit at the place value ('logplace') position of number.
example:
      numberdigit(13625, 3) = 3
      numberdigit(13625, 4) = 1
*)
var
      place: integer; (* the exponent of logplace *)
      count: integer; (* used to make place *)
      absolute: 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:=absolute-((absolute div tenplace)*tenplace);
      if place = 1
         then d:=z
         else d:= z div place;
      case d of
         0: acharacter:='0';
         1: acharacter:='1';
         2: acharacter:='2';
         3: acharacter:='3';
         4: acharacter:='4';
         5: acharacter:='5';
         6: acharacter:='6';
         7: acharacter:='7';
         8: acharacter:='8';
         9: acharacter:='9';
      end
end; (* digit *)
procedure sign;
(* put a negative sign out or a positive sign *)
begin (* sign *)
      if number <0 then acharacter:='-'
                   else acharacter:='+'
end; (* sign *)
begin (* numberdigit *)
      place:=1;
      for count:=1 to logplace do place:=10*place;

      if number=0 then begin
         if place=1 then acharacter:='0'
                    else acharacter:=' '
      end
      else begin
         absolute:=abs(number);
         if absolute < (place div 10)
            then acharacter:=' '
            else if absolute >= place
                    then digit
                    else sign
      end;
      numberdigit:=acharacter
end; (* numberdigit *)
(* end module numberdigit version = 'prgmod 3.96  85 mar 18 tds'; *)

(* begin module numbersize *)
function numbersize(n: integer):integer;
(* calculate amount of space to be reserved for the integer n *)
const ln10 = 2.30259; (* natural log of 10 - for conversion to log base 10 *)
      epsilon = 0.00001; (* a small number to correct log base 10 errors *)
begin (* numbersize *)
      if n = 0
         then numbersize:=1
         else numbersize:=trunc(ln(abs(n))/ln10 + epsilon) + 2;
(* the epsilon assures that we do not lose a place due to
   roundoff.  eg, sometimes log base 10 of 10 would be
   0.9999 instead of 1, and we would not do it right...
   note: this will fail for very large numbers on the order of 1/epsilon. *)
(* the 2 is for the sign and last digit *)
end; (* numbersize *)
(* end module numbersize version = 'prgmod 3.96  85 mar 18 tds'; *)

(* begin module numberbar *)
procedure numberbar(var afile: text; spaces, firstnumber, lastnumber: integer;
                    var linesused: integer);
(* write a bar of numbers to a file, with several spaces before.
   the number of lines used is returned *)
var
      logplace: integer; (* the log of the digit being looked at *)
      spacecount: integer; (* count of spaces *)
      number: integer; (* the current number being written *)
begin
      if abs(firstnumber) > abs(lastnumber)
         then linesused:= numbersize(firstnumber)
         else linesused:= numbersize(lastnumber);

      for logplace:=linesused-1 downto 0 do begin
         for spacecount:=1 to spaces do write(afile,' ');
         for number:=firstnumber to lastnumber
            do write(afile,numberdigit(number,logplace));
         writeln(afile)
      end
end;
(* end module numberbar version = 'prgmod 3.96  85 mar 18 tds'; *)
(* ************************************************************************ *)
(* end module package.numbar version = 'prgmod 3.96  85 mar 18 tds'; *)

(* begin module complines *)
function complines(var a,b: text): boolean;
(* compare one line from each of files a and b starting at the
   current file positions.  both lines are compared up to and including
   the end-of-line mark (i.e., they are both readln-ed).  return true if
   the lines are identical and false if they are not.  ignore trailing
   blanks.  this is done by letting chara or charb assume the
   value of a blank when eoln is true.  when both lines are eoln,
   the next line is inspected.  if one line is at end-of-file and
   the other is not, 'false' is returned.  if both lines are at
   eof, 'true' is returned. *)
var
      chara,charb: char; (*  characters in files a and b *)
      done: boolean; (* stop looking for identity *)
begin
      if eof(a) and eof(b)
         then complines := true
         else if eof(a) or eof(b)
            then complines := false
            else begin (* neither file is eof. *)
               done := false;
               while not done do begin
                  if eoln(a) then chara := ' '
                             else read(a,chara);

                  if eoln(b) then charb := ' '
                             else read(b,charb);

                  if eoln(a) and eoln(b)
                     then begin
                        done := true;
                        complines := true
                     end
                     else if chara <> charb then begin
                        done := true;
                        complines := false
                     end;
               end;
               readln(a);
               readln(b);
            end
end; (* complines *)
(* end module complines version = 'auxmod 1.37 85 apr 4 gds/tds'; *)

(* begin module hlistvsbooks *)
procedure hlistvsbooks (var hlist : text; (* list of helices from helix *)
                     var xbook,ybook : text (* the delila books *) );
(* checks the title line in xbook with the corresponding line in hlist,
   ditto ybook.  if they are not identical, error messages are sent to
   output and the program halts. written by britta singer *)
var   xasx, (* true if xbook input file corresponds to xbook in hlist *)
      yasy, (* true if ybook input file corresponds to ybook in hlist *)
      xasy, (* true if xbook input file corresponds to ybook in hlist *)
      yasx: (* true if ybook input file corresponds to xbook in hlist *)
                     boolean;

procedure gettitleline(var xbook,ybook,hlist : text);
      (* gets to the lines in the xbook, ybook, and hlist input files
         that contain the titles and thus are to be compared *)
begin (* gettitleline *)
      reset(xbook);
      reset(ybook);
      reset(hlist);
      readln(hlist)
end; (* gettitleline *)

procedure findtitle(var hlist : text);
      (* places the cursor at the beginning of the title in hlist *)
var    n : integer; (* counter for skipping characters *)
begin (* findtitle *)
      for n := 1 to 4 do
         get(hlist)
end; (* findtitle *)

begin (* hlistvsbooks *)
      gettitleline(xbook,ybook,hlist);
      findtitle(hlist);
      xasx := complines(xbook,hlist);
      findtitle(hlist);
      yasy := complines(ybook,hlist);
      gettitleline(xbook,ybook,hlist);
      findtitle(hlist);
      yasx := complines(ybook,hlist);
      findtitle(hlist);
      xasy := complines(xbook,hlist);
      if not xasx
         then writeln(output,
            ' xbook input file does not correspond to xbook in hlist.');
      if not yasy
         then writeln(output,
            ' ybook input file does not correspond to ybook in hlist.');
      if xasy and (not xasx)
         then writeln(output,
            ' xbook input file corresponds to ybook in hlist.');
      if yasx and (not yasy)
         then writeln(output,
            ' ybook input file corresponds to xbook in hlist.');
      if xasy and yasx and (not xasx) and (not yasy)
         then writeln(output,
            ' xbook and ybook input files are reversed. try again.');
      if not (xasx and yasy)
         then halt
end; (* hlistvsbooks *)
(* end module hlistvsbooks version = 'auxmod 1.37 85 apr 4 gds/tds'; *)

procedure initialize;
(* set up all the variables *)
var
      i: integer; (* index for counting lines in hlist *)
      getenergy: boolean; (* do we want to get energies?  this is used only
         to check that minlenmaxene is consistent with the hlist *)
begin (* initialize *)
      writeln(output,' matrix ',version:4:2);
      brinit(xbook);
      brinit(ybook);

      hlistvsbooks(hlist,xbook,ybook);
      reset(hlist);
      rewrite(mlist);
      unlimitln(mlist);

      new(xseq);
      new(yseq);

      (* put header in mlist *)
      writeln(mlist,' matrix ',version:4:2,' from:');
      if copylines(hlist,mlist,headerlines) <> headerlines then begin
         writeln(output,' hlist is too short to even have a correct header');
         halt
      end;

      (* read minlenmaxene from hlist *)
      reset(hlist);
      for i := 1 to 3 do readln(hlist);
      readln(hlist,minlenmaxene);
      readln(hlist);
      get(hlist); (* skip the blank *)
      if hlist^ = 'e' then getenergy := true
                      else getenergy := false;
      for i := 6 to headerlines do readln(hlist);

      (* read parameters *)
      reset(matrixp);
      if not eof(matrixp) then readln(matrixp,minlenmaxene);

      (* check consistency of this minlenmaxene with the hlist *)
      if not getenergy
         then if minlenmaxene <= 0.0
            then begin
               write  (output,' no energies listed in this hlist,');
               writeln(output,' change parameter');
               writeln(output,' to request positive lengths');
               halt
            end;


      if minlenmaxene > 0
      then writeln(mlist,' only helixes longer than or equal to ',
                   trunc(minlenmaxene):1,' base pairs are printed')
      else writeln(mlist,' only helixes stronger than or equal to ',
                   minlenmaxene:5:2,' kcal are printed');
      writeln(mlist);
      writeln(mlist,' symbols used:');
      write(mlist,' ',nodot,' every helix through this point is ');
      if minlenmaxene > 0
      then writeln(mlist, 'less than ', trunc(minlenmaxene):1,
                          ' base pairs long.')
      else writeln(mlist,'weaker than ', minlenmaxene:5:2,' kcals.');
      writeln(mlist,' 1 gt base pair');
      writeln(mlist,' 2 at base pair');
      writeln(mlist,' 3 gc base pair');
end; (* initialize *)

(*****************************************************************************)
(* procedures to access the bit array as if it were dimensioned
   [1..xlength], 1..ylength]
      clearbitarray, writebitarray, readbitarray *)

procedure setindex(x,y: integer);
(* set the value of index based on (x,y) and (xlength, ylength) *)
begin
      index:=x+(y-1) * xlength
end;

procedure clearbitarray;
(* clear bitarray using xlength and ylength globals *)
var
      index: integer; (* index to the array *)
      maximum: integer; (* the amount of bitarray used *)
begin
      maximum:=xlength * ylength;
      if maximum > maxarray then begin
         writeln(output,' array size (maxarray) exceeded');
         halt
      end;

      for index:=1 to maximum do bitarray[index]:=zero
end;

procedure writebitarray(x,y: integer);
(* turn on bitarray at (x,y) *)
begin
(*write(output,'writebitarray',x,y);debug*)
      setindex(x,y);
      bitarray[index]:=one
(*;if bitarray[index]=one then writeln(output,'wrote a one')
                        else writeln(output,'wrote a zero');*)
end;

function readbitarray(x, y: integer): bit;
(* read the bit at (x,y) in bitarray *)
begin
      setindex(x,y);
(*;if bitarray[index]=one then writeln(output,'read a one====at ',index)
                        else writeln(output,'read a zero at ',index);*)
      readbitarray:=bitarray[index]
end;

(* begin module findcolon *)
procedure findcolon(var thefile: text);
(* move the file to the characters just past ': ' (colon, space) *)
var
      found: boolean;
procedure die;
begin writeln(output,' no helix list data'); halt end;
begin
      found:=false;
      while not found do begin
         if eof(thefile)
            then die
            else if thefile^ = ':'
               then begin
                  get(thefile);
                  if eof(thefile)
                     then die
                     else if thefile^ = ' '
                        then begin
                           get(thefile);
                           found:=true
                        end
               end
               else get(thefile)
      end
end;
(* end module findcolon version = 'auxmod 1.37 85 apr 4 gds/tds'; *)

procedure checknames;
(* check that the names in the hlist correspond
to those in the books.  hlist is ready to read the names. *)
begin
      readln(hlist) (* ignore for now *)
end;

(* begin module gethelix *)
procedure gethelix(var hlist: text;
                   var x,y: integer;
                   var length: integer;
                       getenergy: boolean;
                   var energy: real;
                   var done: boolean);
(* get the (x,y,length) info on some helix listed in hlist.
if the x,y pair changes, then done is true and no values are
returned.  checknames may be used at this point to start the
next pair.  if getenergy is true, then the routine expects that there
will be an energy on the helix line *)
begin (* gethelix *)
      if eof(hlist)
         then done:=true
         else begin
            get(hlist); (* skip the blank *)
            if hlist^ <> ' '
               then done:=true
               else begin
                  done:=false;
                  findcolon(hlist); read(hlist,x);
                  findcolon(hlist); read(hlist,y);
                  findcolon(hlist); read(hlist,length);
                  if getenergy then read(hlist, energy)
                               else energy := 0;
                  readln(hlist)
(*;writeln(output,'got helix x,y,length,energy=',x,y,length,energy);*)
               end
         end
end; (* gethelix *)
(* end module gethelix version = 'auxmod 1.37 85 apr 4 gds/tds'; *)

procedure fillbitarray;
(* fill bitarray with the hlist helixes *)
var
      done: boolean; (* done with the fill *)
      x0,y0: integer; (* 5 prime ends of the helix *)
      l: integer; (* index to the helix *)
      getenergy: boolean; (* are there energies to get? *)
      energy: real; (* the energy of a helix *)
      x,y: integer; (* location in the array *)
begin (* fillbitarray *)
      done:=false;
      getenergy := (minlenmaxene <= 0);
      while not done do begin
         gethelix(hlist,x0,y0,length,getenergy,energy,done);
(*writeln(output,'(',x0:1,',',y0:1,')0 ',length:1); *)
         if (not done) and
            (((not getenergy) and (length >= minlenmaxene)) or
             (getenergy and (energy <= minlenmaxene))         ) then begin
            x:=pietoint(x0,xseq);
            y:=pietoint(y0,yseq) + (length-1); (* to 3 prime end *)
(*writeln(output,'(',x:1,',',y:1,') ',length:1); *)
            for l:=1 to length do begin
               writebitarray(x,y);
               x:=x+1;
               y:=y-1
            end
         end
      end
end; (* fillbitarray *)

(* begin module dotmat.piecebar *)
procedure piecebar(var afile: text;
                       spaces: integer;
                       pie: pieceptr;
                       firstnumber,lastnumber: integer;
                   var linesused: integer);
(* write a bar of numbers of the piece pie to a file,
with several spaces before.  the number of lines used is returned.
based on procedure numbar.  external procedures required:
  inttopie
  numbersize
  numberdigit
*)
var
      logpower: integer; (* the log of the digit being looked at *)
      spacecount: integer; (* count of spaces *)
      number: integer; (* the current number being written *)
      numbers: file of integer; (* a file for storing all the piece numbers *)
      width: integer; (* the lines used by one number *)
      index: integer; (* to pie *)
begin
      (* calculate lines used and fill up the numbers file *)
      linesused:=0;
      rewrite(numbers);
      for index:=firstnumber to lastnumber do begin
         number:=inttopie(index,pie);
         write(numbers,number);
         width:=numbersize(number);
         if width > linesused then linesused:=width
      end;

      for logpower:=linesused-1 downto 0 do begin
         for spacecount:=1 to spaces do write(afile,' ');
         reset(numbers);
         while not eof(numbers) do begin
            read(numbers,number);
            write(afile,numberdigit(number,logpower))
         end;
         writeln(afile)
      end
end;
(* end module dotmat.piecebar *)

(* begin module dotmat.makedna *)
procedure makedna(var afile: text;
                      spaces: integer;
                      pie: pieceptr;
                      firstnumber, lastnumber: integer);
(* write a line of dna to a file, preceeded by several spaces *)
var
      spacecount: integer; (* count spaces *)
      index: integer; (* to pie *)
begin
      for spacecount:=1 to spaces - 3 do write(afile,' ');
      write(afile,'5',prime,' ');
      for index:=firstnumber to lastnumber
         do write(afile,basetochar(getbase(index,pie)));
      writeln(afile,' 3',prime)
end;
(* end module dotmat.makedna *)

(* begin module dotmat.makepageheader *)
procedure makepageheader(var pageheader: text);
var
      lines: integer; (* dummy until lines are counted by matrix *)
begin
      rewrite(pageheader);
      page(pageheader);
      writeln(pageheader,' matrix, ',
         ' x down: '  ,xnumber:1,' ',xseq^.key.hea.keynam.letters,
         ' y across: ',ynumber:1,' ',yseq^.key.hea.keynam.letters);
      writeln(pageheader);
      piecebar(pageheader,10,yseq,1,ylength,lines);
      makedna(pageheader,10,yseq,1,ylength);
      writeln(pageheader,'        5',prime);   (* 8 spaces *)
end;
(* end module dotmat.makepageheader *)

(* begin module dotmat.printpageheader *)
procedure printpageheader(var pageheader: text);
(* print the page header *)
begin
      reset(pageheader);
      while not eof(pageheader) do begin
         while not eoln(pageheader) do begin
            write(mlist,pageheader^);
            get(pageheader)
         end;
         readln(pageheader);
         writeln(mlist)
      end
end;
(* end module dotmat.printpageheader *)

procedure printbitarray;
(* print out the bit array to mlist *)
var
      pageheader: text; (* the header for pages of mlist *)
      bx, by: base; (* bases at x and y of xseq and yseq respectively *)
      x,y: integer; (* indices *)
begin (* print bit array *)
      makepageheader(pageheader);
      for x:=1 to xlength do begin
         if ((x mod pagelength) = 0) or (x = 1)
            then printpageheader(pageheader);
         bx := getbase(x, xseq);
         write(mlist,' ',inttopie(x,xseq):6,' ',basetochar(bx),' ');
         for y:=1 to ylength
            do if readbitarray(x,y) = one
               then begin (* make a number at the place *)
                  by := getbase(y, yseq);
                  if bx = complement(by)
                  then case bx of
                     a: write(mlist,'2');
                     c: write(mlist,'3');
                     g: write(mlist,'3');
                     t: write(mlist,'2');
                  end
                  else write(mlist,'1') (* for gt pairs *)
               end
               else write(mlist,nodot);
         writeln(mlist)
      end;
      writeln(mlist,'        3',prime);   (* 8 spaces *)
end;

procedure printmatrix;
(* print the matrix for xseq and yseq *)
begin
      checknames;
      clearbitarray;
      fillbitarray;
      printbitarray
end;

begin
      initialize;
      while not eof(xbook) do begin
         getpiece(xbook,xseq);
         xnumber:=number;
         xlength:=piecelength(xseq);
         if not eof(xbook) then begin
            reset(ybook);
            while not eof(ybook) do begin
               getpiece(ybook,yseq);
               ynumber:=number;
               ylength:=piecelength(yseq);
               if not eof(ybook) then begin
                  printmatrix;
                  clearpiece(yseq)
               end
            end;
            clearpiece(xseq)
         end
      end;
1: end. (* of matrix *)