program curvafea(book, curvafeap, curvaout, curvafeamarks, curvaxyin,
                 output);

(* curvafea: DNA curvature to lister features

  Dr. Thomas D. Schneider
  National Cancer Institute
  Laboratory of Experimental and Computational Biology
  Frederick, Maryland  21702-1201
  toms@ncifcrf.gov
  permanent email: toms@alum.mit.edu (use only if first address fails)
  http://www.lecb.ncifcrf.gov/~toms/

 *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 1.06; (* of curvafea.p 2003 Sep 2
2003 Sep  2, 1.06: reduce output
2003 Aug 31, 1.05: accept blank and * lines in curvaout
2003 Aug 27, 1.04: sequence orientation bug
2003 Aug 26, 1.03: functinal
2003 Aug 26, 1.02: create curvaxyin
2003 Aug 26, 1.01: functional display of curvature on lister map!
2003 Aug 25, 1.00: origin from live.p 1.53
*)
updateversion = 1.00; (* defines lowest acceptable current parameter file *)
(* end module version *)

(* begin module describe.curvafea *)
(*
name
   curvafea: DNA curvature to lister features

synopsis
   curvafea(book: in, curvafeap: in, curvaout: in,
            curvafeamarks: out, curvaxyin: out,
            output: out)

files

   book:  A book from the delila system.

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

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

      2. zerobase: character and one integer OR one integer.  The coordinate
         to set the phase of the wave by.  This base will be red.

         If the first character is 'f' then the first base is the zerobase.

         If the first character is 'l' then the last base is the zerobase.

         When the first character is 'f' or 'l', it is followed by the
         zerobase as an offset.  So f 0 is the first base, f 1 is the second
         base, etc.

         If the first character is not 'f' or 'l' then an integer is given
         that is the phase offset.

         When a set of coordinates are used (6th parameter and beyond), 'f'
         and 'l' apply to each coordinate pair.  As a result, the phases of
         the colors WILL NOT MATCH.  If neither 'f' nor 'l' are used, then
         there is only one zero base for all pairs, so THE PHASES WILL
         MATCH.

      3. wavelength: one real.  The wave length in bases.  Usually there are
         10.6 base pairs per turn of DNA, but using 10.0 allows checking the
         program.  It can have other useful value, for example, one can set
         it to 3 to see translation frames.

      4. wavevertical (real):  how much to move the colored bar up or down
         (in character heights = bits).  This is preadjusted (by a constant
         (shiftdown) to fit just between the lister map numbering and the
         sequence.

      5. waveheight (real):  the height of the colored bar (in character
         heights = bits).  0.25 is the normal value that fits the line
         between the lister map numbering and the sequence.

      6. and all following lines: firstbase lastbase: two integers per line.
         The first and last bases to create marks for.  There may be an
         indefinite number of lines.  You can also use the program with
         different parameters and concatenate the curvafeamarks files together.
         WARNING:  Coordinate systems can jump in the Delila system, so this
         method will fail across jumps.  You can set the ranges by hand.
         This method can be used to mark only portions of the sequence.
         The program bookends can be used to find the end numbers of book
         pieces.

         IF PARAMETERs 6 ARE MISSING then the book will be used to obtain the
         coordinate system.  This method should always work, but it gives
         curvafea marks for the entire book.

   curvaout:  One or more '.map' files from the curva program.  These
      give the coordinate followed by the curvature on each line.  To
      make it possible to concatenate .map files together with
      spacers, blank lines and '*' at the starts of line are accepted
      (the latter as comments).

   curvafeamarks:  marks for the lister program

   curvaxyin:  curve in xyin format for the xyplo program

   output: messages to the user

description

   The program makes waves for lister by using the marks definition.  The
   waves are colored strips on the number line (or wherever the lister line
   starts).  The strip shows the phase of the DNA, allowing the orientation
   of binding sites to be compared.

   Since lister marks must be in order of the coordinate system, integrating
   the curvafea marks with other marks would require a lot of hand editing.  The
   mergemarks program allows marks to be smoothly joined together.

examples

1.34      version of curvafea that this parameter file is designed for.
f 0       zerobase: to set the phase of the wave by.  This base will be red.
10.6      wavelength: The wave length in bases.  Usually 10.6.
0         wavevertical: amount to shift colored wave vertically
0.25      waveheight: height of colored wave 

Optionally, there can be:

-150  92  firstbase lastbase: range of sequence to mark
... (any number of repeats of the last line above may be made)

documentation

   The CURVATURE (ie curva) program is described at:

@article{Shpigelman.Bolshoy1993,
author = "E. S. Shpigelman
 and E. N. Trifonov
 and A. Bolshoy",
title = "{CURVATURE: software for the analysis of curved DNA}",
journal = "Comput Appl Biosci",
volume = "9",
pages = "435--440",
year = "1993"}

see also

   {example parameter file:} curvafeap

   {The curves program is available from:}
   http://ext.haifa.ac.il/~leokol/curvature/

   {A script that takes a book as input and uses this program  to
   generate various files that can be used with Delila programs is:}
   run.curva

   {the program that uses the curvafeamarks:} lister.p
   {another program that could use curvafea marks:}  makelogo.p
   {a program to merge curvafeamarks with others: } mergemarks.p
   {a program to obtain the range of a book: } bookends.p

   {a program for plotting the curvaxyin: } xyplo.p

   {Documentation on how the curvafea program produces the color wave:}
   http://www-lecb.ncifcrf.gov/~toms/livecolors/

   {This program was created by modifying: } live.p

author
   Thomas Dana Schneider

bugs

   NOTES from live.p:

   When one gives a pair of coordinates, curvafea creates marks running between
   the pair.  This will *fail* when there is a coordinate jump.  Also, it
   means that if the coordinate system is decreasing, the wave colors will be
   *reversed*.  When the book is used for coordinates, the wave is determined
   on internal (1 to n) coordinates and so the colors will never be
   reversed.  It is not clear yet whether the latter solution is a problem or
   not.  (Is it a bug or a feature?  :-)  Are there times when one would like
   the colors to reverse when the piece reverses?

   The answer to the above is YES: the colors must reverse for curva!
zzz check that this is true

   The code could be imbedded into lister.

technical notes

   The definition of the marks is in makelogo.p

   More information about the colors is in xyplo.p

*)
(* end module describe.curvafea *)

(* begin module curvafea.const *)
   shiftdown = -0.350; (* amount to shift the strip down, so that it fits
                          nicely between the numbering and DNA. *)
   saturation = 0.5; (* how saturated the colors are *)
(* end module curvafea.const *)

(* begin module book.const *)
(* constants needed for book manipulations *)
      dnamax = 1024;  (* length of dna arrays *)
      namelength = 100; (* maximum key name length *)
      linelength = 80; (* maximum line readable in book *)
(* end module book.const version = 7.56; {of delmod.p 2000 Nov 21} *)

type

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

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

      (* types defined in book definition *)

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

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

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

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

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

(* begin module base.type *)
   (* define the four nucleotide bases *)
   base = (a,c,g,t);
(* end module base.type version = 7.56; {of delmod.p 2000 Nov 21} *)

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


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

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

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


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

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

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

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

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

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

(* end module book.type version = 7.56; {of delmod.p 2000 Nov 21} *)

   parameters = record
     (* parameters to control the program *)
     parameterversion: real; (* parameter version number *)
     firstbase: integer; (* first base to mark *)
     lastbase: integer; (* last base to mark *)
     zerobase: integer; (* the base to be 'red' *)
     wavelength: real; (* wave length in bases *)
     wavevertical: real; (* wave vertical shift in bits *)
     waveheight: real; (* height of the bar in bits *)
     basephase: char; (* phase starts on the first (f) or last (l) base,
                         ' ' means specified zero *)
   end;

var
   book, (* file used by this program *)
   curvafeap, (* file used by this program *)
   curvaout, (* file used by this program *)
   curvafeamarks, (* file used by this program *)
   curvaxyin  (* file used by this program *)
      : text; (* file used by this program *)

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

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

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

(* ************************************************************************ *)
(* end module book.var version = 7.56; {of delmod.p 2000 Nov 21} *)

(* begin module halt *)
procedure halt;
(* stop the program.  the procedure performs a goto to the end of the
   program.  you must have a label:
      label 1;
   declared, and also the end of the program must have this label:
      1: end.
   examples are in the module libraries.
   this is the only goto in the delila system. *)
begin
      writeln(output,' program halt.');
      goto 1
end;
(* end module halt version = 7.56; {of delmod.p 2000 Nov 21} *)

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

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

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

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

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

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

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

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

procedure cleardna(var l: dnaptr);
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 *)
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 = 7.56; {of delmod.p 2000 Nov 21} *)

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

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

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

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

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

(* 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;
      while (not eoln(thefile)) do begin
         i:=succ(i);
         read(thefile,acharacter);
         l^.letters[i]:=acharacter
      end;
      l^.length:=i;
      l^.next:=nil;
      readln(thefile);
      theline := succ(theline)
end;
(* end module book.brline version = 7.56; {of delmod.p 2000 Nov 21} *)

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

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

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

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

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

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

         (* read note key *)
         if readnumber then brnotenumber(thefile,theline,note)
                       else brnote(thefile,theline,note)
      end
end;
(* end module book.brheader version = 7.56; {of delmod.p 2000 Nov 21} *)

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

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

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

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

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

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

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

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

(* ************************************************************************ *)
(* end module package.brpiece version = 7.56; {of delmod.p 2000 Nov 21} *)

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

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

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

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);
   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);
                     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);
                     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;
(* end module book.getocp version = 7.56; {of delmod.p 2000 Nov 21} *)

(* ************************************************************************ *)
(* end module package.getocp version = 7.56; {of delmod.p 2000 Nov 21} *)

procedure readparameters(var curvafeap: text; var p: parameters);
(* read the parameters p from curvafeap *)
begin
   with p do begin
      reset(curvafeap);
      readln(curvafeap, parameterversion);
      if parameterversion < updateversion then begin
         writeln(output, 'You have an old parameter file!');
         halt
      end;

      if (curvafeap^ = 'f') or (curvafeap^ = 'l')
      then begin
         readln(curvafeap,basephase, zerobase);
      end
      else begin
         basephase := ' ';
         readln(curvafeap,zerobase);
      end;

      readln(curvafeap,wavelength);
      readln(curvafeap,wavevertical);
      readln(curvafeap,waveheight);
      (* the remainder of the parameters are read later: *)
      (* readln(curvafeap,firstbase, lastbase); *)

      (* checks *)

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

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

   end;
end;

procedure writeparameters(var f: text; var p: parameters);
(* read the parameters p from curvafeap *)
const
   wid = 10; (* width for putting numbers *)
   dec =  5; (* decimal for putting numbers *)
begin
   with p do begin
     writeln(f, '* ', parameterversion:wid:2, ' parameterverion');
     writeln(f, '* ', firstbase:wid, ' firstbase');
     writeln(f, '* ', lastbase:wid, ' lastbase');
{
writeln(output,'basephase =',basephase);
}
     case basephase of
        'f': writeln(f, '* first base + ',zerobase:wid,' is zerobase');
        'l': writeln(f,  '* last base + ',zerobase:wid,' is zerobase');
        ' ': writeln(f, '* ', zerobase:wid, ' zerobase');
     end;
     writeln(f, '* ', wavelength:wid:dec, ' wavelength');
     writeln(f, '* ', wavevertical:wid:dec, ' wavevertical');
     writeln(f, '* ', waveheight:wid:dec, ' waveheight');
   end;
end;

procedure makecurvafeamarks(var curvaout, f: text; p: parameters);
(* using the curve data in curvaout, make curvafea marks into file f*)
begin

rewrite(f);
writeln(f,'* curvafea ',version:4:2,', create curve features');
writeln(f);
writeparameters(f, p);

writeln(f);
writeln(f,'u'); (* user definitions *)
writeln(f,'% fs is the font size vertically');
writeln(f,'% charwidth is assumed to exist (it does in makelogo and lister)');
writeln(f,'% as the width of the characters horizontally.');
writeln(f,'% define font size as a unit of measure:');
writeln(f,'/fs {fontsize mul} def');

(*
      Definitions of marks according to makelogo:

           U Call the user defined symbol.  The U must be followed by 4
             coordinates numbers: x1 y1 x2 y2.  The x1 and x2 are in bases,
             while y1 and y2 are in bits.  The remainder of the line is copied
             to the logo file, so you can have more parameters there.  End the
             line with the name of one of your defined symbols.
           * a comment line
           % a comment line

        If marksymbol is b, l or U, 4 more parameters are required:
           base coordinate: a real number that determines end 1
           bits coordinate: a real number that determines end 1
           base coordinate: a real number that determines end 2
           bits coordinate: a real number that determines end 2
        The line is drawn from end 1 to end 2 while these ends define box
        diagonal.  Note that the center of a base is defined as an integer, so
        one must add 0.5 to base coordinates to put a boxes around a base.  You
        may make the user symbol use these coordinates however you want.

*)

writeln(f);
writeln(f,'/wavebox { % TailX TailY HeadX HeadY wbcolor wavebox -');
writeln(f,'% make a box for the wave n charcters vertically from zero');
writeln(f,'/wbcolor exch def');
writeln(f,'/HeadY exch def');
writeln(f,'/HeadX exch def');
writeln(f,'/TailY exch def');
writeln(f,'/TailX exch def');
writeln(f,'gsave');
writeln(f,'% shift left 1 character:');
writeln(f,'charwidth neg 0 translate');
(* writeln(f,'charwidth neg wbcolor 0.5 sub cm translate'); for testing*)
writeln(f,'TailX TailY moveto');
writeln(f,'HeadX TailY lineto');
writeln(f,'HeadX HeadY lineto');
writeln(f,'TailX HeadY lineto');
writeln(f,'TailX TailY lineto');
writeln(f,'wbcolor ',saturation:5:3,' 1 sethsbcolor');
writeln(f,'fill');
writeln(f,'grestore');
(*
%    TailX     TailY     HeadX     HeadY     Shift
%    bases      bits     bases      bits     bases
*)
writeln(f,'} def');
writeln(f,'!'); (* end of user definitions *)
writeln(f);

end;

procedure markat(var curvafeamarks: text; x, y, waveheight, color: real);
(* mark a box at base x and y characters vertically.
The procedure is a wrapper to make generation easy.
The lister program needs two user defined coordinates
("head" and "tail" of an arrow or corners of a rectangle).
These are not convenient in this program.  Here it is easier
to think about the base location we want to box.
*)
begin
(*
writeln(curvafeamarks, '*    TailX     TailY     HeadX     HeadY     color');
writeln(curvafeamarks, '*    bases      bits     bases      bits          ');
*)
writeln(curvafeamarks, 'U ',x:10:5,
                    ' ',y:10:5,
                    ' ',(x+1):10:5,
                    ' ',(y+waveheight):10:5,
                    ' ',color:10:5,
                    ' wavebox');
end;

function setcolor(thezero: integer;
                  internal: integer;
                  wavelength: real;
                  huge: integer
                  ): real;
(* Set the color to cycle over the spectrum, using internal coordinates.
thezero is where the base should be tagged red, internal is the base being
worked on, wavelength is the wave length of the color cycle change.

This is a particularly difficult routine because one needs to map the
difference between internal and thezero into a number between 0 and 1.  The
hard part is that functions like trunc act *differently* on the two sides of
the origin, which would make an unacceptable jump in the final color.
Methods that try to account for the trunc function near zero tend to fail.
The simplest solution is to get away from zero by adding a huge number,
taking the trunc, and then subtracting away the huge number.  This huge
number doesn't really need to be very large, it only needs to be beyond the
range of the coordinate system.

*)

var
   thecolor: real; (* color computed from position *)
   x: real;  (* the distance from the zero coordinate, in wavelengths *)
   y: real;  (* a number between 0 and 1, cycling according to x *)
begin

   x := (internal - thezero)/wavelength;
   y := x + huge - trunc(x+huge);

   thecolor := 0.85 * y; (* rescale to make spectrum *)
   setcolor := thecolor;

{
write(output,'thezero = ',thezero:6:3,' ');
write(output,'wavelength = ',wavelength:6:3,' ');
}
{
write(output,'internal = ',internal:6:3,' ');
write(output,'huge = ',huge:6:3,' ');
writeln(output,'x = ',x:6:3, ', y=',y:6:3);
}
{
write(output,'internal-thezero=',(internal-thezero):5);
write(output,'thecolor = ',thecolor:5:3);
write(output,', thezero  = ',thezero:5);
write(output,', internal = ',internal:5);
writeln(output);
}

end;

procedure nobookloop(thezero: integer; (* zero point *)
                     cooexternal: integer; (* external coordinate *)
                     p: parameters; 
                     huge: integer); 
(* do the inner part of a loop *)
var
   thecolor: real; (* color computed from position *)
begin
   with p do begin
      thecolor := setcolor(thezero, cooexternal, wavelength, huge);
{writeln(output,'thecolor = ',thecolor:5:3);}
      markat(curvafeamarks,cooexternal,shiftdown+wavevertical, waveheight, thecolor);
   end;
end;

procedure bookloop(thezero: integer; (* zero point *)
                   internal: integer; (* internal coordinate *)
                   pie: pieceptr; (* piece of DNA *)
                   p: parameters; 
                   huge: integer;
                   var cinternal: integer;
                   var curvaout, curvaxyin: text
                   ); 
(* do the inner part of a loop, with zero base at thezero and
internal as the *internal* coordinate.
Read the curve data from curvaout. *)
(* thezero: the computed zero base *)
const
   wid = 10; (* width    of numerical output *)
   dec =  5; (* decimals of numericaloutput *)
   saturation = 1.0; (* color definition *)
   brightness = 1.0; (* color definition *)
var
   hue: real; (* color computed from position *)
   curvature: real; (* the curvature from curvaout *)
   piecoordinate: integer; (* piece coordinates *)
   sign: integer; (* orientation of the piece *)
begin
   with p do begin
      with pie^.key do begin
         if piedir = coodir
         then sign := +1
         else sign := -1
      end;
{
This is the part of the live.p program that is replaced:
      hue := setcolor(thezero, internal, wavelength, huge);
}
{writeln(output,'hue = ',hue:5:3);}

{
write(output,    ' internal: ',  internal:3);
}

      if (cinternal <= internal) and (not eof(curvaout))
      then begin
         if cinternal < internal
         then begin
            if eoln(curvaout)
            then readln(curvaout) (* skip blank lines *)
            else if curvaout^ = '*'
            then readln(curvaout) (* skip comment lines *)
            else readln(curvaout, cinternal, curvature)
         end
         else begin (* use previous value *) end;

{
write(output, ' cinternal: ', cinternal:3);
write(output, ' curvature: ',curvature:5:3);
}

{
writeln(curvaxyin, ' curvature: ',curvature:5:3);
if internal > 15 then begin
writeln(output); halt;
end;
}

         if internal = cinternal then begin
{
            write(output, ' internal');
writeln(output); writeln(output,' HOLD ON NELLY'); halt;
}

            if curvature <= 1
            then hue := 0.84*curvature + 0.16
            else hue := 1.0; (* safety feature *)

            piecoordinate := inttopie(internal, pie);
            markat(curvafeamarks, piecoordinate, {inttopie(internal, pie),}
                   shiftdown+wavevertical, waveheight, hue);

            writeln(curvaxyin,      internal:wid,
                               ' ', curvature:wid:dec,
                               ' ', hue:wid:dec,
                               ' ', saturation:wid:dec,
                               ' ', brightness:wid:dec,
                               ' ', sign*piecoordinate:wid
                               );

{zzz here also write out the xyin data so it can be
plotted on a logo ... oh!  do like a sine wave?}

         end;

      end
      else begin
{
            write(output, ' -------');
}
      end;
{
writeln(output,' %');
}
   end;
end;

(* begin module curvafea.themain *)

procedure themain(var book, curvafeap, curvaout, curvafeamarks,
                  curvaxyin: text);
(* the main procedure of the program *)
var
   internal: integer; (* internal coordinate corresponding to l *)
   inthezero: integer; (* the computed zero base in internal coordinates *)
   l: integer; (* position on sequence *)
   p: parameters; 
   thezero: integer; (* the computed zero base in external coordinates *)
   huge: integer; (* a large positive integer *)

   (* parameters for reading the book: *)
   theline: integer;
   org: orgkey;   orgchange, orgopen: boolean;
   chr: chrkey;   chrchange, chropen: boolean;
   pie: pieceptr; piechange, pieopen:  boolean;

   cinternal: integer; (* the internal coordinate according
                          to curvaout *)

begin
   writeln(output,'curvafea ',version:4:2);
   readparameters(curvafeap, p);
   makecurvafeamarks(curvaout, curvafeamarks, p);

   reset(curvaout);
   cinternal := 0;

   rewrite(curvaxyin);
   writeln(curvaxyin,'* curvafea ',version:4:2);
   writeln(curvaxyin,'* columns:');
   writeln(curvaxyin,'* 1 internal coordinate');
   writeln(curvaxyin,'* 2 curvature units (relative to nucleosome bending)');
   writeln(curvaxyin,'* 3 hue for color');
   writeln(curvaxyin,'* 4 saturation for color');
   writeln(curvaxyin,'* 5 brightness for color');
   writeln(curvaxyin,'* 6 piecoordinate: coordinate on the book');

   if not eof(curvafeap) then begin
      writeln(output,'Warning: not using book; coordinates may not be correct:');
      writeln(output,'assuming 1 to n coordinates');
      huge := trunc(maxint/2);
      while not eof(curvafeap) do with p do begin
         readln(curvafeap,firstbase, lastbase);
         case basephase of
            'f': thezero := firstbase + zerobase;
            'l': thezero := firstbase + zerobase;
            ' ': thezero :=             zerobase;
         end;
         inthezero := thezero; (* assuming 1 to n coordinates *)
         if firstbase <= lastbase
         then for l := firstbase to     lastbase do nobookloop(inthezero,l,p, huge)
         else for l := firstbase downto lastbase do nobookloop(inthezero,l,p, huge)
      end;
   end
   else begin
      (* read from the book *)
      brinit(book,theline);
      orgchange := false; orgopen := false;
      chrchange := false; chropen := false;
      piechange := false; pieopen := false;
      pie := nil;
      with p do repeat
         getocp(book, theline,
                org, orgchange, orgopen,
                chr, chrchange, chropen,
                pie, piechange, pieopen);
         if pieopen then begin
            if number > 1 then begin
               writeln(curvafeamarks,'p - curvafea skip to piece ',number:1);
               writeln(curvafeamarks);
            end;
            writeln(curvafeamarks,'* piece ',number:1);
{
writeln(output, zerobase:5,' is zerobase');
writeln(output, 'basephase = "',basephase,'"');
}
            case basephase of
               'f': begin
                    thezero := pietoint(pie^.key.piebeg, pie) + zerobase;
{ this can't be right! double call to the pietoint!
                    inthezero := pietoint(thezero, pie);
}
                    inthezero := thezero;
               end;
               'l': begin
                    thezero := pietoint(pie^.key.pieend, pie) + zerobase;
{ see above
                    inthezero := pietoint(thezero, pie);
}
                    inthezero := thezero;
               end;
               ' ': begin
                    thezero :=                                  zerobase;
                    inthezero := thezero;
               end;
            end;
{
writeln(output, inthezero:5,' is inthezero');
}
            huge := 2* piecelength(pie);
            for internal := 1 to piecelength(pie)
            do bookloop(inthezero,internal,pie,p, huge,
                        cinternal, curvaout, curvaxyin);
         end
      until not pieopen
   end
end;
(* end module curvafea.themain *)

begin
   themain(book, curvafeap, curvaout, curvafeamarks, curvaxyin);
1: end.