program dnaplot(dnain, dnaout, dnaplotp, positions, dnasymbols, output);
(* dnaplot: plot Ri values for one or more large DNA sequences
   modules: prgmod, dops *)

label 1; (* end of program *)

const
(* begin module version *)
   version = 3.48; (* of dnaplot.p 2005 May 23
2005 May 23, 3.47: upgrade for gpc, patent message
2000 Jul 25, 3.46: remove pascru bomb
1998 June 8, 3.45: previous
1996 October 16: piece names go onto the plot
origin before 1993 August 11 *)
(* end module version *)

(* begin module describe.dnaplot *)
(*
name
   dnaplot: plot Ri values for one or more large DNA sequences

synopsis
   dnaplot(dnain: in, dnaout: out, dnaplotp: in,
           positions: in, dnasymbols: in, output: out)

files
   dnain:  A data input file created by scan.  It contains header lines that
       begin with asterisks ('*') that are copied to dnaout.  Remaining lines
       are the data in columns, ending with end of file.

   dnaout: output in PostScript format

   dnaplotp: Parameter file for dnaplot, which is configured as followed:

       frompos topos         These two integers are the positions on the
                             sequence which the graph will represent.  If,
                             instead, the first character on the line is 'r',
                             then these numbers are read from the positions
                             file for each sequence.

       sCol pCol cCol vCol   columns numbers for reading from the dnain file:
                             sCol = piece number
                             pCol = column number of piece name
                             cCol = piece coordinate
                             vCol = value (Ri evaluation)

       numperpg              number of graphs per page
       numperln              number of base pairs per line
       bitlower bitupper     lower and upper bounds of bits to display

       orix oriy             x, y origin of plot (in cm)
       xaxlength yaxlength   length of the x and y axes in cm
       showaxis              t=true means show coordinate axis to dnaout
       xinterval yinterval   size of intervals on axes to plot
       xsubint   ysubint     x and y sub intervals to mark
       xwidth    ywidth      width of numbers in characters
       xdecimal  ydecimal    number of decimal places
       xticlength xticdx xticdy  length of tic mark and shift of number (cm)
       yticlength yticdx yticdy  length of tic mark and shift of number (cm)
       sequencelabel         t=true means print sequence number on graphs
       xaxislabel            the label for the x axis
       yaxislabel            the label for the y axis
       plottype              How to draw the plot:
                             z = lines from zero to value
                             b = lines from bottom of graph to value
       dodash                Whether to put vertical dashed lines around
                             segments of continuous sequence.  This is
                             important for distinguishing between the absence
                             of sequence and low Ri values, but this often gets
                             in the way, so it can be turned off:
                             d = do dashes
                             n = no dashes
       thelinewidth          The width of the plot lines in bases.

   positions:  If the first character on the first line of the parameter file
      is an 'r' then this file will be read to determine the positions to graph
      for each sequence.  The file consists of pairs of integers, one pair per
      line, representing the first (frompos) and last (topos) coordinates to be
      plotted.

   dnasymbols: If the file is not empty, then it contains information on how
      and where to plot special symbols to make marks on the graph.  Each line
      has 5 values:
         symbol (first character): the type of symbol to draw:
            c = circle
            s = square
            t = triangle 
         symboltype (second character): the way to draw the symbol:
            s = stroked as a solid line
            f = filled
            d = dotted line 
         symbolplacement (third character): where to put the symbol on the graph:
            a = use absolute location (given by symbolbits) on graph
            r = use relative location (given by symbolbits) from current Ri value
         symbolbits (real): the distance in bits
         symbolsize (real): diameter of a circle or side of square and triangle
            relative to the spacing between graph lines.  A value of 1 fits
            between the lines.
         piece number (integer):  The number of the fragment to mark,
            as given in the dnain file
         piece coordinate (integer):  the coordinate on the piece to mark
            as given in the dnain file

        ********************************************************************
        * The symbols MUST be in increasing order of position in the plot! *
        ********************************************************************

        Lines that are empty or begin with "*" are ignored.

   output: messages to the user

description
   dnaplot creates a PostScript graph of information content (or
   other values) versus position on a DNA sequence.

examples

dnasymbols
The line:
csa 15 0.5 1 100
means place a circle, stroked, at absolute 15 bits, 0.5 size, 
for sequence 1 at coordinate 100.    

documentation

see also
   scan.p, xyplo.p, dbbk.p

author
   Stacy L. Bartram, modified by Tom Schneider

bugs

   The program cannot handle negatively numbered base systems because the axis
   cannot give decreasing numbers (yet or ever).  These simply come out as
   blank (or the program will halt, depending on the version).

technical notes

*)
(* end module describe.dnaplot *)

(* constants continued *)
(* begin module dnaplot.const *)
   infofield = 12; (* size of field for printing information in bits *)
   infodecim = 6; (* number of decimal places for printing information *)
       (* these are used for conlist only *)
   nfield = 4; (* size of field for printing n, the number of sites *)
(* end module dnaplot.const *)

(* begin module pic.const *)
   pi = 3.14159265354; (* circumference divided by diameter of circle *)
   picwidth = 8; (* width of numbers printed to the file *)
   picdecim = 5; (* number of decimal places for numbers *)
   charwidth = 0.0625; (* the width of characters in inches (ie, inches/char)
                        this allows centering of strings. *)
   (* note: for the Times-Roman  font, 0.0625 is a good value.
            for the Courier-Bold font, 0.08   is a good value. *)
   dotfactor = 0.00625; (* the size of dots *)
   defscale = 28.35; (* default scale factor.  coordinate units per cm *)
(* end module pic.const version = 2.66; (@ of dops.p 1994 Oct 6 *)

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

(* begin module datetime.const *)

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

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

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

(* begin module dnaplot.type *)
   params = record
     frompos, topos: integer; (* positions on sequence graph will represent *)
     readpositions: boolean; (* if true, obtain frompos and topos from
                                the positions file *)

     (* columns of the dnain file:
       sequence number (integer), piece name (string),
       coordinate number (integer), value (real) *)
     sCol: integer; (* sequence number *)
     pCol: integer; (* piece name *)
     cCol: integer; (* coordinate number *)
     vCol: integer; (* value *)

     numperpg: integer; (* number of horizontal lines per page *)
     numperln: integer; (* number of base pairs per line *)
     bitspercm: real;   (* height of vertical lines for information content *)
     bitlower, bitupper: real;  (* lower and upper bounds of bits to display *)

     orix, oriy: real;     (* x, y origin of plot (in cm) *)
     xaxlength, yaxlength: real; (* length of the x and y axes in cm *)
     showaxis: boolean;    (* show coordinate axis *)
     xinterval, yinterval, (* number of intervals on axes to plot *)
     xsubint, ysubint,     (* number of sub intervals on axes to mark *)
     xwidth, ywidth,       (* width of numbers in characters *)
     xdecimal, ydecimal: integer; (* number of decimal places *)
     xticlength, xticdx, xticdy, (* length of tic mark, shift of number (cm) *)
     yticlength, yticdx, yticdy: (* length of tic mark, shift of number (cm) *)
        real;
     sequencelabel: boolean; (* true means print sequence number on graphs *)
     xaxislabel, (* label for x axis *)
     yaxislabel: (* label for y axis *)
        string;
     plottype: char; (* type of plot to produce *)
     dodash: char; (* do dashes or not *)
     thelinewidth: real; (* the width of plot lines in bases *)
   end;
(* end module dnaplot.type *)

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

{ junk this silly code eventually:
     infonum: real;        (* information above which to report to dnaout *)
     MagicNumber: real;    (* print Ri even if it is below infonum *)
from parameters:
* the next two are not valid in this version and must not be included
       infonum               information above which to report to dnaout
       MagicNumber           print Ri even if it is below infonum
from technical notes
   In Kenn Rudd's database, there are unknown sections of the E. coli genome.
   Delila was designed with the idea that unknown sequence would always be
   determined.  So the program dbbk which converts to Delila format substitutes
   an A at every unknown base.  When a string of A's is evaluated by an Ribl
   matrix, it gives a unique value, called the "MagicNumber" that can therefore
   be used to detect the regions of unknown sequence.  This program is set up
   so that when it comes across the MagicNumber, it assumes the sequence is
   missing.  If the user requests that values below zero not be plotted, the
   positions with the MagicNumber will be plotted anyway, to indicate the
   regions where there is no sequence.
in writing loop:
          if (vVal = MagicNumber) then vVal := 0;
}

var
(* begin module dnaplot.var *)
   dnain, (* input data from scan *)
   dnaout, (* output - PostScript graph instructions *)
   dnaplotp, (* file from which to read the parameters *)
   positions, (* file from which to read the positions to plot *)
   dnasymbols: text; (* file from which to read the plot symbols *)
(* end module dnaplot.var *)

(* begin module pic.var *)
   inpicture: boolean; (* true if we are drawing the picture,
                          ie, startpic has been called *)
   picxglobal, picyglobal: real; (* absolute location in the graph *)
   pictolerance: real; (* 10 raised to the picwidth,
      to detect values close to zero *)

   scale: real; (* scale factor.  graphic coordinate units per inch *)

(* NONSTANDARD for efficient use of postscript, keep track of
whether there is a current path *)
   inpath: boolean;

(* NONSTANDARD keep track of number of segments drawn so that
they can be stroked.  This (probably) solves the problem of the
Apple printer dying because it can't handle the data. *)
   segments: integer;

   xsideold, ysideold: real; (* current size of a rectangle.  see rectsize *)
(* end module pic.var version = 2.66; (@ of dops.p 1994 Oct 6 *)

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

(* begin module package.datetime *) 
(* ************************************************************************ *)
(* begin module getdatetime *)
procedure getdatetime(var adatetime: datetimearray);
(* get the date and time into a single array from the system clock.
   adatetime contains the date:
       1980/06/09 18:49:11
         ye mo da ho mi se
(year, month, day, hour, minute, second).
As of 2000 February 18, the Sun Pascal compiler requires a formatting
statement.  This statement allows the date to be generated in this
standard Delila format in a single call.  Information about the
formatting statement is available on the manual page for date in Unix.
If a computer does not have this method, see the 'oldgetdatetime' routine
in delmod.p (https://alum.mit.edu/www/toms/delila/delmod.html)
for some conversion code.

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

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

7.10.8 Date And Time Routines 

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

DateLength and TimeLength are implementation dependent constants. 

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

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

*)
var
    t: TimeStamp; 

(* begin module pluckdigit *)
function pluckdigit(number, logplace:integer): char;
(* return the digit at the place value ('logplace') position of number.
example:
      pluckdigit(13625, 3) = 3
      pluckdigit(13625, 4) = 1
This routine was taken from module numberdigit in prgmod.p, but is
modified so as not to give the sign.  Instead it gives zeros
above the digits.  'myabsolute' replaced 'absolute', which
is apparently a keyword for GPC.  The name is kept for
to keep the code looking similar to its origin.
*)
var
      place: integer; (* the exponent of logplace *)
      count: integer; (* used to make place *)
      myabsolute: integer; (* the absolute value of number *)
      acharacter: char; (* the character to be returned *)
procedure digit;
(* extract a digit at the place position *)
var
      tenplace: integer; (* ten times place *)
      z: integer; (* an intermediate value *)
      d: integer; (* the digit extracted *)
begin (* digit *)
      tenplace:=10*place;
      z:=myabsolute-((myabsolute div tenplace)*tenplace);
      if place = 1
         then d:=z
         else d:= z div place;
      case d of
         0: acharacter:='0';
         1: acharacter:='1';
         2: acharacter:='2';
         3: acharacter:='3';
         4: acharacter:='4';
         5: acharacter:='5';
         6: acharacter:='6';
         7: acharacter:='7';
         8: acharacter:='8';
         9: acharacter:='9';
      end
end; (* digit *)
begin (* pluckdigit *)
      place:=1;
      for count:=1 to logplace do place:=10*place;
      if number=0 then begin
         acharacter:='0'
      end
      else begin
         myabsolute:=number;
         if myabsolute >= place
         then digit
         else acharacter := '0'
      end;
      pluckdigit:=acharacter
end; (* pluckdigit *)
(* end module pluckdigit *)

begin

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

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

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

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

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

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

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

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

procedure timeseed(var seed: real);
(* read the computer date and time.  reverse the order of the digits
and put a decimal point in front.  this gives a fraction between
zero and one that varies quite quickly, and is always unique (if the
computer has sufficient accuracy).  it is to be used as a seed to
a random number generator. *)
var
   adatetime: datetimearray; (* a date and time *)
begin (* timeseed *)
   getdatetime(adatetime);
{
writeln(output,'timeseed: adatetime: ',adatetime);
}
   makeseed(adatetime, seed);
end; (* timeseed *)
(* end module timeseed version = 1.17; (@ of timegpc.p 2002 Oct 9 *)
 
(*[[*)
(* begin module limitdate *)
procedure limitdate(a,b,c,d: char; limitdatetime: datetimearray);
(* test whether the current time is before the limit.
If it is later, halt the program *)
var
   adatetime: datetimearray; (* a date and time *)
   Dday: real; (* the critical day *)
   now: real;  (* this very moment *)
begin
   getdatetime(adatetime);
{
writeln(output,'adatetime:',adatetime);
writeln(output,'adatetime[1]: ', adatetime[1]);
writeln(output,'adatetime[2]: ', adatetime[2]);
writeln(output,'adatetime[3]: ', adatetime[3]);
writeln(output,'adatetime[4]: ', adatetime[4]);
writeln(output,'adatetime[5]: ', adatetime[5]);
}

   orderseed(adatetime, now);

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

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

   orderseed(limitdatetime, Dday);

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

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

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

(* begin module interact.getstring *)
procedure getstring(var afile: text; var buffer: string;
                    var gotten: boolean);
(* get a line (as a string) from a file not using string calls.  this lets
one obtain lines from a file without interactive prompts *)
var   index: integer; (* of buffer *)
begin (* getstring *)
      clearstring(buffer);
      if eof(afile)
      then gotten := false
      else begin
         index := 0;
         while (not eoln(afile)) and (index < maxstring) do begin
            index := succ(index);
            read(afile, buffer.letters[index])
         end;
         if not eoln(afile) then begin
            writeln(output, ' getstring: a line exceeds maximum string size (',
                            maxstring:1,')');
            halt
         end;
         buffer.length := index;
         buffer.current := 1;
         readln(afile);
         gotten := true
      end
end; (* getstring *)
(* end module interact.getstring version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

(* begin module skipblanks *)

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

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

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

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

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

(* begin module onetoken *)
procedure onetoken(var afile: text; var buffer: string;
                  var gotten: boolean);
(* Get a string from a file not using string calls.  This lets one
obtain lines from a file without interactive prompts.  if end of file
is found, gotten is false. *)
const
   tabcharacter = 9; (* the ordinal of the tab character *)
var   index: integer; (* of buffer *)
      done: boolean;
begin (* onetoken *)
   skipblanks(afile);
   emptystring(buffer);
   done := false;
   if eof(afile)
   then gotten := false
   else begin
      index := 0;
      while (not eoln(afile)) and (index < maxstring)
         and not done
      do begin
         index := succ(index);
         read(afile, buffer.letters[index]);
         if    (buffer.letters[index] = ' ')
            or (ord(buffer.letters[index]) = tabcharacter)
         then begin
            done := true;
            index := pred(index);
         end
      end;
      buffer.length := index;
      buffer.current := 1;
      gotten := true
   end
end;
(* end module onetoken version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

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

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

(* ****************************************************************** *)
(* graphics for axes *)

(* begin module dnaplot.startpic *)
procedure startpic(var afile:text; setscale,x,y: real; thefont: char);
(* open the graphics field, with the given scale, and at (x,y)
in that scale.  scale is in device coordinates per inch.
The font is chosen with thefont; t = Times-Roman, c = Courier-Bold *)
(* start pic output to file afile, set the globals *)
(* NONSTANDARD *)
(* this is the actual "world" coordinates used: *)
(*                        xmin,  xmax,  ymin, ymax *)
(*    ns; if (setwindow(-5.0/scale, +5.0/scale, -5.0/scale, +5.0/scale)*)
begin
   scale := setscale; (* set the global scale *)

   case thefont of
      'c': begin
           writeln(afile,'/Courier-Bold findfont'); (* locate the font *)
           writeln(afile,10:1,' scalefont'); (* set the font size in points*)
      end;
      't': begin
           writeln(afile,'/Times-Roman findfont'); (* locate the font *)
           writeln(afile,12:1,' scalefont'); (* set the font size in points*)
      end;
   end;

   writeln(afile,'setfont'); (* put the font into the current font *)

   (* set the scale to inches
   writeln(afile,
                 scale:picwidth:picdecim,' ',
                 scale:picwidth:picdecim,' scale'); *)

   (* define some things in postscript *)
   (* doline allows less stuff to be put in the output file.
      it takes two numbers off the stack, copies them, draws a line
      to them as coordinates. *)
(* replaced by 'currentpoint translate'
   writeln(afile,'/doline { 2 copy lineto } def');
*)
   (* define a function that makes inches out of a number *)
(* do this all internally here, it's faster
   writeln(afile,'/i { ',scale:picwidth:picdecim,' mul} def');
*)

   (* move to the start point on the page *)
   writeln(afile,
                     (x*scale):picwidth:picdecim,
                 ' ',(y*scale):picwidth:picdecim,
                 ' translate');

    writeln(afile);
    writeln(afile,'% Define functions so the text produced is smaller');
    writeln(afile,'/a {stroke newpath 0 0} def     % special for arc');
    writeln(afile,'/c {stroke 0 0 moveto} def      % current point');
    writeln(afile,'/f {findfont 10 scalefont setfont} def');
    writeln(afile,'   % to set fonts simply use the f function.  Example:');
    writeln(afile,'   %/Symbol f (\142) /Courier-Bold f (-galactosidase');
    writeln(afile,'/l {lineto} def');
    writeln(afile,'/m {moveto} def');
    writeln(afile,'/n {stroke newpath 0 0 moveto} def'); (* new segment *)
    writeln(afile,'/rl {rlineto} def');
    writeln(afile,'/rm {rmoveto} def');
    writeln(afile,'/s {newpath 0 0 moveto} def     % Start path ');
    writeln(afile,'/t {currentpoint translate} def % translate ');
    writeln(afile,'/x {show} def                   % show teXt ');
    writeln(afile);

   (* start out the pathway *)
   inpath := false;
   (* start the number of segments written: *)
   segments := 0;

   (* now for the normal pic stuff: *)
   inpicture := true;
   picxglobal := 0.0;
   picyglobal := 0.0;
   pictolerance := trunc(exp(picwidth*ln(10))+0.5)
(*;writeln(output,'pictolerance = ',pictolerance:picwidth:picdecim);*)
end;

(* MODIFIED from pic.startpic version = 2.66; (@ of dops.p 1994 Oct 6 *)
(* end module dnaplot.startpic *)

(* begin module pic.drawr *)
procedure drawr(var afile: text; dx,dy: real; visibility: char;
               spacing: real);
(* make a line to file afile by relative draw of dx,dy with visibility
  i invisible
  - dashed
  . dotted
  l line
with the dashes or dots separated by the spacing given
(this has no effect with invisible and line). *)
(* NONSTANDARD *)
var
   ddx,ddy: real; (* changes in dx and dy for dots and dashes *)
   dr: real; (* the hypotenuse, the distance actually drawn *)
   on: boolean; (* draw linesegment if true *)
   y: real; (* the variable for tracking dots and dashes *)
   r: integer; (* number of times to cycle for dots and dashes *)
   ss: real; (* precalculated value to make things a bit faster *)
   theta: real; (* angle of the line *)
procedure checkseg(var afile: text);
(* NONSTANDARD checks how many segments have been written, if
more than 'buffer', stroke them to the postscript page *)
const buffer = 10;
begin
   if segments >= buffer
   then begin
   (* New segment: writeln(afile,'stroke newpath 0 0 moveto'); *)
      writeln(afile,'n');
      segments := 0
   end
   else segments := segments + 1;
end;
begin (* drawr *)
   if not inpath then begin
      (* starts from current coordinates *)
      (* Start path: writeln(afile,'newpath 0 0 moveto'); *)
      writeln(afile,'s');
      inpath := true
   end
   else checkseg(afile);
(* checks
if not (visibility in ['l','i','.','-'])
then writeln(afile,'%YELLLLLL!!!',visibility,'!');
writeln(afile,'% ',visibility,' line');*)

   (* put these on the stack, they will always be used *)
   write(afile,    (dx*scale):picwidth:picdecim,
               ' ',(dy*scale):picwidth:picdecim);
   case visibility of
  'l','i': begin
              case visibility of
                 'i': write(afile,' m');
                 'l': write(afile,' l');
              end
           end;
  '.','-': begin (* make up our own dots and dashes *)
              writeln(afile); (* move away from the (dx,dy) on the stack *)
              if spacing <= 0.0 then begin
                 writeln(output,'drawr: spacing zero with . or - line');
                 halt
              end;
              if dx = 0.0
              then begin
                 ddx := 0.0; (* avoid division by zero *)
                 ddy := scale*spacing;
                 if dy < 0 then ddy := - ddy; (* this makes sure that
                    we draw lines straight down if that was the request *)
              end
              else begin
                 (* find out the angle of the slope, intentionally
                    lose the sign *)
                 theta := arctan(abs(dy/dx));

                 ddx := scale*spacing*cos(theta);
                 ddy := scale*spacing*sin(theta);

                 (* return the sign to the little buggers *)
                 if dx < 0 then ddx := -ddx;
                 if dy < 0 then ddy := -ddy;
              end;
              y := 0;
              case visibility of
                 '.': ss := scale*dotfactor;
                 '-': on := true;
              end;

              dr := sqrt(dx*dx+dy*dy);
              for r := 1 to round(dr/spacing) do begin
                 case visibility of
                 '-': begin
                         write(afile,    (ddx):picwidth:picdecim,
                                     ' ',(ddy):picwidth:picdecim);
                         if on
                         then writeln(afile,' rl')
                         else writeln(afile,' rm');
                         on := not on
                      end;
                 '.': begin
                         (* put out a dot like in dotr *)
                         write(afile,
                                    +ss:picwidth:picdecim,' 0 rl');
                         write(afile,' ',
                                    -ss:picwidth:picdecim,' 0 rl');
                         write(afile,' ',(ddx):picwidth:picdecim,
                                     ' ',(ddy):picwidth:picdecim);
                         writeln(afile,' rm');
                      end;
                 end
              end;
              (* let's make really sure we got there!! *)
              writeln(afile,' m'); (* pulled from the stack *)
           end;
   end;

   (* an elegant way to make postscript keep a global record is
   to translate the coordinates! *)
   (* writeln(afile,' currentpoint translate'); *)
   writeln(afile,' t');

   picxglobal := picxglobal + dx;
   picyglobal := picyglobal + dy;
end;
(* end module pic.drawr version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.mover *)
procedure mover(var afile: text; dx,dy: real);
(* move relative the amount (dx, dy). *)
begin
   drawr(afile,dx,dy,'i',0.0);
end;
(* end module pic.mover version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.liner *)
procedure liner(var afile: text; dx,dy: real);
(* draw a line the relative amount (dx, dy). *)
begin
   drawr(afile,dx,dy,'l',0.0);
end;
(* end module pic.liner version = 2.66; (@ of dops.p 1994 Oct 6 *)


(* begin module pic.drawa *)
procedure drawa(var afile: text; x,y: real; visibility: char;
               spacing: real);
(* make a line to file afile to absolute coordinate x,y with visibility
  i invisible
  - dashed
  . dotted
  l line
with the dashes or dots separated by the spacing given
(this has no effect with invisible and line). *)
var
   dx, dy: real; (* differences between current and desired
                    locations *)
begin
   dx := x - picxglobal;
   dy := y - picyglobal;

   drawr(afile,dx,dy,visibility,spacing)
end;
(* end module pic.drawa version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.movea *)
procedure movea(var afile: text; x,y: real);
(* move to absolute x and y *)
begin
   drawa(afile,x,y,'i',0.0);
end;
(* end module pic.movea version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.linea *)
procedure linea(var afile: text; x,y: real);
(* draw a line from current position to absolute x and y *)
begin
   drawa(afile,x,y,'l',0.0);
end;
(* end module pic.linea version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.graphstring *)
procedure graphstring(var tofile: text; var s: string; justification: char);
(* graph the string s.  If it is recognized as a quoted string (surrounded
by double quotes), graph it without the quotes and center it.
Otherwise justify it based on the justification character:
'l' left, 'c' centered, 'r' right.
For right and centered justification, the drawing point is the same
as before the string was done.  For left justification it is at the
right of the string to allow more to be added on there.
If not in picture (global variable inpicture), there is no output *)
(* NONSTANDARD: PostScript dependent code.  Since different fonts
have different sized characters, one must rely on the PostScript
to handle the justification of the string. *)
var i: integer; (* index to s, and temporary storage *)
    quoted: boolean; (* true if the string is quoted *)
    skipping: boolean; (* true if skipping leading blanks *)
begin
   if (inpicture and (s.length > 0))
   then with s do begin
      if length > 2
      then if (letters[1]='"') and (letters[length]='"')
           then quoted := true
           else quoted := false
      else quoted := false;

      (* override so quoted strings are always centered *)
      if quoted then justification := 'c';

      (* do the non-standard postscript: *)
      if justification <> 'l' then write(tofile,'gsave ');

      (* do postscript to complete pervious path *)
      (* set current point: writeln(tofile,'stroke 0 0 moveto'); *)
      writeln(tofile,'c');

      if justification = 'c' then begin
         (* when centering, skip leading blanks *)
         if letters[1] = ' ' then skipping := true
                             else skipping := false;
      end
      else skipping := false;

      write(tofile,'('); (* begin postscript literal *)
      if quoted (* take it literally *)
      then for i := 2 to length-1 do write(tofile,letters[i])
      else for i := 1 to length
           do if skipping then begin (* skip leading blanks *)
                 if letters[i] <> ' '
                 then begin
                      skipping := false;
                      write(tofile,letters[i])
                 end
              (* else skip the blank by not writing it *)
              end
              else write(tofile,letters[i]);
      write(tofile,')'); (* end postscript literal *)

      if justification = 'c' (* center the string *)
      then write(tofile,' dup stringwidth pop neg 2 div 0 rmoveto')
      else if justification = 'r' (* rigth justify the string *)
      then write(tofile,' dup stringwidth pop neg 0 rmoveto');

      writeln(tofile,' x'); (* show the literal *)
      inpath := false; (* force new path from here *)

      if justification <> 'l' then write(tofile,'grestore ');
   end
(* There is no output if not in picture
   else begin
      writestring(tofile,s);
      writeln(tofile)
   end
 *)
end;
(* end module pic.graphstring version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.stringinteger *)
procedure stringinteger(number: integer; var name: string;
                        width: integer; leadingzeros: boolean);
(* make the string from the number, start putting characters in
after the current length point. use width characters.
if leadingzeros is true, trail zeros before the number. *)
var
   bigdigit: integer; (* the location of the biggest digit *)
   dig: integer; (* number of digits in the number *)
   place: integer; (* place to write the next digit of the number *)
   sign: integer; (* the sign of the number *)
begin
   with name do begin
      if number < 0
      then begin
         sign := -1;
         length := length + 1; (* provide room for the sign!! *)
         number := -number;
         if leadingzeros then begin
            writeln(output,'WARNING: stringinteger: the sign of a negative',
                           ' number with leading zeros is lost');
         end
      end
      else sign := +1;

      (* log 10 of the number plus 1 is the number of digits in the number.
      On this sun computer ln(1000)/ln(10) is 2.9999, which when
      truncated gives 2, rather than the desired 3.  To avoid this
      kind of problem, 0.1 is added. *)
      if number > 9
      then dig := trunc(ln(number+0.1)/ln(10))+1
      else dig := 1;

      if dig > width then begin
         writeln(output,'stringinteger: number width too small');
         writeln(output,dig:1,' digit number (',number:1,')');
         writeln(output,'does not fit in ',width:1,' characters');
         halt
      end;
      if leadingzeros
      then bigdigit := length + 1 (* no sign if leading zeros *)
      else begin
         bigdigit := length + width - dig + 1;

         if (bigdigit <= length) and (sign < 0) then begin
            writeln(output,'stringinteger: no room for sign');
            halt
         end;
      end;
      if sign < 0 then letters[bigdigit-1] := '-';

      for place := length + width downto bigdigit do begin
         case (number mod 10) of
            0: letters[place] := '0';
            1: letters[place] := '1';
            2: letters[place] := '2';
            3: letters[place] := '3';
            4: letters[place] := '4';
            5: letters[place] := '5';
            6: letters[place] := '6';
            7: letters[place] := '7';
            8: letters[place] := '8';
            9: letters[place] := '9';
         end;
         number := number div 10;
      end;
      length := length + width;
   end
end;
(* end module pic.stringinteger version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.stringreal *)
procedure stringreal(number: real; var name: string;
                     width, decimal: integer);
(* make the string from the real number, start putting characters in
at the start point. use width characters and decimal characters
after the decimal place *)
   (* note that the rounding operation to get the digits below zero
      must be done first.  then the digits above zero can be lopped off.
      this makes 99.99 come out correctly to 100.0 (to 1 decimal place)
      otherwise, 99.99 -> 0.99 -> 1.0 (rounded) -> 10 (print with 1 decimal
      place), and stringinteger won't be happy about that. *)
var
   abovezero: integer; (* the number shifted above the decimal place, to
      'decimal' positions (and rounded) *)
   shift: integer; (* power of ten used to shift a number around
      relative to the decimal point *)
   sign: integer; (* the sign of the number *)
   thedecimal: integer; (* integer version of the decimal part of the number *)
   theupper: integer; (* integer version of the upper part of the number *)
begin
   if number < 0 then sign := -1
                 else sign := +1;

   number := abs(number); (* make positive *)

   (* the amount to shift the number above zero *)
   shift := round(exp(decimal*ln(10))); (* amount to move above zero *)
   abovezero := round(number*shift); (* move above zero, round off *)
   theupper := trunc(abovezero/shift);
   thedecimal := abovezero - shift*theupper;

   (* create the actual real number *)
   (* before decimal point *)
   stringinteger(sign*theupper,name,width-decimal-1,false);
   with name do begin (* put in the decimal point *)
      length := length + 1;
      letters[length] := '.';
   end;
   stringinteger(thedecimal,name,decimal,true); (* after decimal point *)
end;
(* end module pic.stringreal version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.picnumber *)
procedure picnumber(var afile: text;
                    dx, dy, number: real; width, decimal: integer;
                    justification: char);
(* Supply graphic commands for a 'number' whose center is at the relative point
(dx, dy) from the current point, 'width' characters wide and 'decimal'
characters beyond the decimal point.
If the width is zero, no number is produced.
procedure stringnumber(number: integer; start: integer; var name: string);
the location after the call is the same as before the call.
The string is optionally justified: left, centered or right: lcr. *)
var
   name: string; (* the string to pack the number into for shipping out *)
begin
   if width > 0 then begin
      mover(afile,dx,dy);

      clearstring(name);

      if decimal>0
      then stringreal(number,name,width,decimal)
      else stringinteger(round(number),name,width,false);

      graphstring(afile, name, justification);
      mover(afile,-dx,-dy);
   end
end;
(* end module pic.picnumber version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.xtic *)
procedure xtic(var afile: text; length, dx, dy, number: real;
               width, decimal: integer;
               logxnormal: boolean;
               logxbase: real);
(* produce a tic mark for the x axis of "length" long.
Supply a number whose center is at the relative point (dx, dy)
from the end to the tick, 'width' characters wide and 'decimal'
characters beyond the decimal point.
If the width is zero, no number is produced.
the location after the call is the same as before the call.
If logxnormal is true, then raise the number
to logxbase. *)
begin
   liner(afile,0.0,-length);
   if logxnormal
   then picnumber(afile,dx,dy,exp(number*logxbase),width,decimal,'c')
   else picnumber(afile,dx,dy,number,width,decimal,'c');
   mover(afile,0.0,length);
end;
(* end module pic.xtic version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.ytic *)
procedure ytic(var afile: text; length, dx, dy: real;
               number: real;
               width, decimal: integer;
               logynormal: boolean;
               logybase: real);
(* produce a tic mark for the y axis of "length" long.
Supply a number whose right side is started at the relative point (dx, dy)
from the end to the tick, 'width' characters wide and 'decimal'
characters beyond the decimal point.
If the width is zero, no number is produced.
the location after the call is the same as before the call.
If logynormal is true, then raise the number
to logybase. *)
begin
   liner(afile,-length,0.0);

   (* convert the number if we are doing logynormal: *)
   if logynormal then number := exp(number*logybase);

   picnumber(afile,dx,dy,number,width,decimal,'r');
   mover(afile,length,0.0);
end;
(* end module pic.ytic version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.doaxis *)
procedure doaxis(var afile: text;
         theaxis: char;
         alength,fromtic,interval,totic: real;
         subintervals: real;
         length, dx, dy: real;
         width, decimal: integer;
         logscale, lognormal: boolean;
         logbase: real);
(* draw an axis starting from the current position.
Which axis it is is defined by theaxis, 'x' (horizontal) or 'y' (vertical).
Combining the code for both axes into one procedure is a little
slower, but drawing the axis does note ever take significant time,
and this allows improvements to be made on both axes simultaneously.
The length of the axis is alength.
The axis is labeled with numbers starting with fromtic
at intervals given up to totic.
The remaining variables describe the form of the tic marks as in ytic.
If the width is zero, no number is produced.
the location after the call is the same as before the call.
If logscale and lognormal is true, then raise the tic numbers to logbase.
 *)
var
   half: real; (* half the jump interval.  By adding this to the while loops,
                  we assure that the very last tic gets done, and isn't lost
                  due to roundoff *)
   jump: real; (* the space to move on the graph between tic marks *)
   jumpdistance: real; (* the total jumps made.  this may not be
      a simple function of the input variables since they may
      not work out to an exact number of jumps *)
   tic: real; (* the numerical value of the tic label *)

   dosubtics: boolean; (* do sub tics *)
   subtic: real; (* the numerical value of the (unlabeled) subtic *)
   subinterval: real; (* the numerical interval between subtics *)
   subjump: real; (* the space to move on the graph between subtic marks *)
   halfsubinterval: real; (* half a subjump, see half *)

   currentspot: real; (* current graphing spot *)
   oldspot: real; (* previous graphing spot *)
   axisscale: real; (* axis scaling factor *)

begin

{
writeln(output,'In doaxis');
writeln(output,'interval=',interval:10:4);
writeln(output,'subintervals=',subintervals:10:4);
writeln(output,'logbase=',logbase:10:4);
}

   if theaxis = 'x' then begin
      liner(afile,+alength,0.0);
      mover(afile,-alength,0.0);
   end
   else begin
      liner(afile,0.0,+alength);
      mover(afile,0.0,-alength);
   end;

   if totic = fromtic then begin
      writeln(output,'doaxis: ',theaxis,' axis fromtic and totic',
                     ' cannot be equal');
      halt;
   end;

   if (alength = 0.0) or (interval = 0.0) then begin
      writeln(output,'doaxis: neither ',
                     theaxis,' axis length nor interval can be zero');

      halt;
   end;

   axisscale := alength / (totic - fromtic);
   jump := axisscale * interval;
   jumpdistance := 0;
   half := interval / 2.0;

   if subintervals > 1 then begin
      dosubtics := true;
      subinterval := interval/subintervals;
      halfsubinterval := subinterval / 2.0;
      subjump := jump/subintervals;
   end
   else begin
      dosubtics := false;
      subinterval := 0;
      halfsubinterval := 0;
      subjump := 0;
   end;

   tic := fromtic;
   if interval > 0.0 then while tic <= totic+half do begin
      if theaxis = 'x'
      then xtic(afile,length,dx,dy,tic,width,decimal,lognormal,logbase)
      else ytic(afile,length,dx,dy,tic,width,decimal,lognormal,logbase);

      tic := tic + interval;
      if tic <= totic+half then begin
{
writeln(output,'TIC=',tic:10:4);
writeln(afile,'% tic=',tic:10:4);
mover(afile,0.05,0.0);
}
         if dosubtics then begin (* do subtic marks *)
            if logscale then begin (* do subtic marks on log scale *)
            (* subtic starts as a "normal" number (ie, no log taken) at tic: *)
{
writeln(output,'2^tic=',exp(tic*logbase):10:4);
writeln(output,'2^(tic+interval)=',exp((tic+interval)*logbase):10:4);
}
               subtic := exp(tic*logbase);
               (* subtic will proceed to the same but at tic+interval.
                  We divide that into the subintervals. *)
{
writeln(output,'halfsubinterval=',halfsubinterval:10:4,' original');
}
               subinterval := (exp((tic+interval)*logbase) - subtic)/subintervals;
               halfsubinterval := subinterval/2.0;

{
writeln(output,'subtic=         ',subtic:10:4);
writeln(output,'subinterval=    ',subinterval:10:4);
writeln(output,'halfsubinterval=',halfsubinterval:10:4);
}

               oldspot := axisscale * tic;
               while subtic < exp(logbase*(tic+interval)) - halfsubinterval
               do begin

                  (* although tic is on a log scale,
                     we have to have subtic on the regular scale
                     to alter the positions of the subtics *)

                  (* if subinterval is constant,
                     the following makes linearly spaced marks: *)
                  subtic := subtic + subinterval;

                  (* the actual jumps have to be in the log form: *)
                  currentspot := axisscale*ln(subtic)/logbase;
                  subjump := currentspot - oldspot;
{
writeln(output,'       SUBTIC=',subtic:10:4);
writeln(output,'       ln(SUBTIC)/logbase=',ln(subtic)/logbase:10:4);
writeln(output,'       currentspot=',currentspot:10:4);
writeln(output,'       subjump=',subjump:10:4);
writeln(output,'       oldspot=',oldspot:10:4);
writeln(afile,'% subtic=',subtic:10:4);
}
                  oldspot := currentspot;


                  if theaxis = 'x' then begin
                     xtic(afile,length/2,dx,dy,0,0,0,lognormal,logbase);
                     mover(afile,subjump,0.0);
                  end
                  else begin
                     ytic(afile,length/2,dx,dy,0,0,0,lognormal,logbase);
                     mover(afile,0.0,subjump);
                  end;

                  jumpdistance := jumpdistance + subjump;
               end
            end
            else begin (* do subtic marks on regular scale *)
               subtic := tic;
               while subtic < tic+interval-halfsubinterval do begin
                  subtic := subtic + subinterval;

                  if theaxis = 'x' then begin
                     mover(afile,subjump,0.0);
                     xtic(afile,length/2,dx,dy,0,0,0,lognormal,logbase);
                  end
                  else begin
                     mover(afile,0.0,subjump);
                     ytic(afile,length/2,dx,dy,0,0,0,lognormal,logbase);
                  end;

                  jumpdistance := jumpdistance + subjump;
               end
            end
         end
         else begin (* do regular tic marks *)
            if theaxis = 'x' then mover(afile,jump,0.0)
                             else mover(afile,0.0,jump);
            jumpdistance := jumpdistance + jump
         end
      end
   end
   else if interval < 0.0 then while tic >= totic-half do begin
      if dosubtics then writeln(output,'Sorry, no subtics with negative scales');

      if theaxis = 'x'
      then xtic(afile,length,dx,dy,tic,width,decimal,lognormal,logbase)
      else ytic(afile,length,dx,dy,tic,width,decimal,lognormal,logbase);

      tic := tic + interval;
      if tic >= totic-half then begin
         if theaxis = 'x' then mover(afile,jump,0.0)
                          else mover(afile,0.0,jump);
         jumpdistance := jumpdistance + jump
      end
   end;

   if theaxis = 'x' then mover(afile,-jumpdistance,0.0)
                    else mover(afile,0.0,-jumpdistance)
end;
(* end module pic.doaxis version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.xaxis *)
procedure xaxis(var afile: text;
         axlength,fromtic,interval,totic: real;
         xsubintervals: real;
         length, dx, dy: real;
         width, decimal: integer;
         logxscale, logxnormal: boolean;
         logxbase: real);
(* draw an x axis starting from the current position. *)
begin
   doaxis(afile,
         'x',
         axlength,fromtic,interval,totic,
         xsubintervals,
         length, dx, dy,
         width, decimal,
         logxscale, logxnormal,
         logxbase)
end;
(* end module pic.xaxis version = 2.66; (@ of dops.p 1994 Oct 6 *)

(* begin module pic.yaxis *)
procedure yaxis(var afile: text;
         aylength,fromtic,interval,totic: real;
         ysubintervals: real;
         length, dx, dy: real;
         width, decimal: integer;
         logyscale, logynormal: boolean;
         logybase: real);
(* draw an y axis starting from the current position. *)
begin
   doaxis(afile,
         'y',
         aylength,fromtic,interval,totic,
         ysubintervals,
         length, dx, dy,
         width, decimal,
         logyscale, logynormal,
         logybase)
end;
(* end module pic.yaxis version = 2.66; (@ of dops.p 1994 Oct 6 *)

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

(* begin module dnaplot.ReadColumns *)
(* This routine is no longer used. *)
procedure ReadColumns(var fin: text;
                      xCol, yCol: integer;
                      var xVal: integer;
                      var yVal: real);
(* reads data xVal (integer) and yVal (real)
from columns xCol and yCol in fin *)
var
  col: integer; (* current column being read *)
  done: boolean; (* done skipping lines *)
begin (* ReadColumns *)
  (* skip data lines *)
  done := false;
  while not done do begin
    if eof(fin)
    then done := true
    else if fin^ = '*'
         then readln(fin)
         else done := true
  end;
  if not eof(fin) then begin
    if xCol < yCol then begin
      for col := 1 to xCol-1 do skipcolumn(fin);
      read(fin, xVal);
      for col := xCol+1 to yCol-1 do skipcolumn(fin);
      read(fin, yVal);
    end
    else if yCol < xCol then begin
      for col := 1 to yCol-1 do skipcolumn(fin);
      read(fin, yVal);
      for col := yCol+1 to xCol-1 do skipcolumn(fin);
      read(fin, xVal);
    end;
    readln(fin);
  end;
end;
(* end module dnaplot.ReadColumns *)

(* begin module dnaplot.skipdata *)
procedure skipdata(var fin: text);
(* skip data lines in fin that begin with asterisk *)
var
  done: boolean; (* done skipping lines *)
begin
  (* skip data lines *)
  done := false;
  while not done do begin
    if eof(fin)
    then done := true
    else if fin^ = '*'
         then readln(fin)
         else done := true
  end;
end;
(* end module dnaplot.skipdata *)

(* begin module dnaplot.grabcolumns *)
procedure grabcolumns(var fin: text;
                      xCol, pCol, yCol, zCol: integer;
                      var xVal: integer;
                      var pVal: string;
                      var yVal: integer;
                      var zVal: real);
(* read data fromfin:
    xVal (integer) yVal (integer) and zVal (real)
from columns
    xCol yCol zCol
The columns must be in order.
The procedure skips blank data columns after reading
so that the end of file is reached after reading the
last piece of data. *)
var
  col: integer; (* current column being read *)
  gotten: boolean; (* piece name was gotten *)
begin
  if not eof(fin) then begin
    for col := 1 to xCol-1 do skipcolumn(fin);
    read(fin, xVal);

    for col := xCol+1 to pCol-1 do skipcolumn(fin);
    onetoken(fin,pVal,gotten);
    if not gotten then begin
       writeln(output,'failed to find piece name');
       halt
    end;
{
write(output,'grabcolumns piece name string: "');
writestring(output,pVal);
writeln(output,'"');
}

    for col := pCol+1 to yCol-1 do skipcolumn(fin);
    read(fin, yVal);

    for col := yCol+1 to zCol-1 do skipcolumn(fin);
    read(fin, zVal);
    readln(fin);
  end;
  skipdata(fin);
end;
(* end module dnaplot.grabcolumns *)

(* begin module dnaplot.readparam *)
procedure readparam(var finp: text; var parameters: params);
(* reads the parameters from finp *)
var
   gotten: boolean; (* true if a label string was read *)
begin (* readparam *)
  with parameters do begin
    reset(finp);

    if finp^='r'
    then begin
       readpositions := true;
       readln(finp);
    end
    else begin
       readpositions := false;
       readln(finp, frompos, topos);
    end;

(*
    readln(finp, infonum);
    readln(finp, MagicNumber);
*)
    readln(finp, sCol, pCol, cCol, vCol);
{
writeln(output, ' s=',sCol:1, ' p=',pCol:1, ' c=',cCol:1,' v=', vCol:1);
}

    readln(finp, numperpg);
    readln(finp, numperln);

    readln(finp, bitlower, bitupper);

    readln(finp, orix, oriy);
    readln(finp, xaxlength, yaxlength);
    if finp^ = 't' then showaxis := true
                   else showaxis := false;
    readln(finp);

    readln(finp, xinterval, yinterval);
    readln(finp, xsubint, ysubint);
    readln(finp, xwidth, ywidth);
    readln(finp, xdecimal, ydecimal);
    readln(finp, xticlength, xticdx, xticdy);
    readln(finp, yticlength, yticdx, yticdy);

    if yaxlength <= 0.0 then begin
      writeln(output, 'ERROR: yaxislength cannot be less than 0');
      halt;
    end;

    bitspercm := (bitupper - bitlower)/yaxlength;

    if (sCol = cCol) or
       (cCol = vCol) or
       (vCol = sCol)
    then begin
      writeln(output, 'ERROR: sCol cCol and vCol columns read',
                      ' cannot be equal');
      halt;
    end;

    if (sCol > cCol) or
       (cCol > vCol) or
       (vCol < sCol)
    then begin
      writeln(output, 'ERROR: sCol cCol and vCol columns read',
                      ' must be in increasing order');
      halt;
    end;

    if sCol <= 0 then begin
      writeln(output, 'ERROR: sCol cannot be less than 0');
      halt;
    end;

    if cCol <= 0 then begin
      writeln(output, 'ERROR: cCol cannot be less than 0');
      halt;
    end;

    if vCol <= 0 then begin
      writeln(output, 'ERROR: vCol cannot be less than 0');
      halt;
    end;

    if finp^ = 't' then sequencelabel := true
                   else sequencelabel := false;
    readln(finp);

    getstring(finp,xaxislabel,gotten);
    if not gotten then writeln(output,'no xaxislabel');
    getstring(finp,yaxislabel,gotten);
    if not gotten then writeln(output,'no yaxislabel');

    readln(finp,plottype);
    if not (plottype in ['z','b']) then begin
      writeln(output, 'plottype must be z or b');
      halt;
    end;

    if eof(finp) then begin
       writeln(output, 'missing dodash parameter');
       halt
    end;

    readln(finp,dodash);
    if not (dodash in ['d','n']) then begin
      writeln(output, 'dodash must be d or n');
      halt;
    end;

    if eof(finp) then begin
       writeln(output, 'missing thelinewidth parameter');
       halt
    end;

    readln(finp,thelinewidth);
    if thelinewidth <=0.0 then begin
      writeln(output, 'thelinewidth must be positive');
      halt;
    end;

  end;
end;
(* end module dnaplot.readparam *)

(* begin module dnaplot.writeparam *)
(* write the finp parameters to output *)
procedure writeparam(var fout: text; parameters: params);
begin (* writeparam *)
with parameters do begin
  writeln(fout, '% **************************************************');
  writeln(fout, '% User Specified Parameters:');

  if readpositions then begin
  writeln(fout, '% reading frompos and topos from the positions file');
  end
  else begin
  writeln(fout, '     /frompos ', frompos:infofield,' def');
  writeln(fout, '       /topos ', topos:infofield,' def');
  end;

(*
  writeln(fout, '     /infonum ', infonum:infofield:infodecim,' def');
  writeln(fout, ' /MagicNumber ', MagicNumber:infofield:infodecim,' def');
*)
  writeln(fout, '        /sCol ', sCol:infofield,' def');
  writeln(fout, '        /pCol ', pCol:infofield,' def');
  writeln(fout, '        /cCol ', cCol:infofield,' def');
  writeln(fout, '        /vCol ', vCol:infofield,' def');

  writeln(fout, '    /numperpg ', numperpg:infofield,' def');
  writeln(fout, '    /numperln ', numperln:infofield,' def');
  writeln(fout, '    /bitlower ', bitlower:infofield:infodecim,' def');
  writeln(fout, '    /bitupper ', bitupper:infofield:infodecim,' def');

  writeln(fout, '        /orix ', orix:infofield:infodecim,' def');
  writeln(fout, '        /oriy ', oriy:infofield:infodecim,' def');
  writeln(fout, '   /xaxlength ', xaxlength:infofield:infodecim,' def');
  writeln(fout, '   /yaxlength ', yaxlength:infofield:infodecim,' def');

  write  (fout, '    /showaxis ');
  if showaxis then write(fout,'true') else write(fout,'false');
  writeln(fout, ' def');

  writeln(fout, '   /xinterval ', xinterval:infofield,' def');
  writeln(fout, '   /yinterval ', yinterval:infofield,' def');
  writeln(fout, '   /xsubint   ', xsubint:infofield,' def');
  writeln(fout, '   /ysubint   ', ysubint:infofield,' def');
  writeln(fout, '      /xwidth ', xwidth:infofield,' def');
  writeln(fout, '      /ywidth ', ywidth:infofield,' def');
  writeln(fout, '    /xdecimal ', xdecimal:infofield,' def');
  writeln(fout, '    /ydecimal ', ydecimal:infofield,' def');

  writeln(fout, '  /xticlength ', xticlength:infofield:infodecim,' def');
  writeln(fout, '      /xticdx ', xticdx:infofield:infodecim,' def');
  writeln(fout, '      /xticdy ', xticdy:infofield:infodecim,' def');

  writeln(fout, '  /yticlength ', yticlength:infofield:infodecim,' def');
  writeln(fout, '      /yticdx ', yticdx:infofield:infodecim,' def');
  writeln(fout, '      /yticdy ', yticdy:infofield:infodecim,' def');
  writeln(fout, '   /bitspercm ', bitspercm:infofield:infodecim,' def');
write(fout, '/xaxislabel: ('); writestring(fout,xaxislabel); writeln(fout,') def');
write(fout, '/yaxislabel: ('); writestring(fout,yaxislabel); writeln(fout,') def');
  writeln(fout, '    /plottype (', plottype, ') def');

  write(fout, '/dodash ');
  if dodash = 'd' then write(fout,'true')
                  else write(fout,'false');
  writeln(fout, ' def');

  writeln(fout, '/thelinewidth ', thelinewidth:infofield:infodecim, ' def');

  writeln(fout, '% **************************************************');
  writeln(fout, '%');
end;
end;
(* end module dnaplot.writeparam *)

(* begin module dnaplot.rightjustifynumber *)
procedure rightjustifynumber(var fout: text; var anumber: integer);
(* number the axis by right justifying the number *)
begin
  (* right justify the string, with a one space gap on the right: *)
  writeln(fout, '(',anumber:1,' )');
  writeln(fout,' dup stringwidth pop neg 0 rmoveto');
  writeln(fout, '(',anumber:1,') show');
  writeln(fout, '0 0 moveto');
end;
(* end module dnaplot.rightjustifynumber *)

(* begin module makelogo.marktype *)
procedure makemarktype(var f: text; marktype: char);
(* make the mark type to the file f *)
begin
   case marktype of
      (* for filling, we stroke first!  This ensures that
         the symbol edges are identical whether stroked or not! *)
      'f': writeln(f,' gsave stroke grestore fill');
      's': writeln(f,' stroke');
      'd': writeln(f,' [3] 0 setdash stroke');
   end;
end;
(* end module makelogo.marktype version = 5.67; (@ of prgmod.p 2017 Sep 29 *)

(* begin module dnaplot.getsymbol *)
procedure getsymbol(var dnasymbols: text;
                    var symbol, symboltype, symbolplacement: char;
                    var symbolbits, symbolsize: real;
                    var symbolpiece, symbolcoordinate: integer);
(* get the symbol information from dnasymbols.  Skip comments
that begin with "*". *)
var
   clearofcomments: boolean; (* are we past all empty lines and comments? *)
begin
   clearofcomments := false;
   repeat
      if eof(dnasymbols)
      then clearofcomments := true
      else if eoln(dnasymbols)
           then readln(dnasymbols)
           else if dnasymbols^ = '*'
                then readln(dnasymbols)
                else clearofcomments := true
   until clearofcomments;

   if eof(dnasymbols) then begin
      symbol := ' ';
   end
   else begin
      read(dnasymbols, symbol);
      read(dnasymbols, symboltype);
      read(dnasymbols, symbolplacement);
      read(dnasymbols, symbolbits);
      read(dnasymbols, symbolsize);
      read(dnasymbols, symbolpiece);
      read(dnasymbols, symbolcoordinate);
      readln(dnasymbols);

      (* make it so the user gives the diameter, and this
         becomes the radius internally. *)
      if symbol = 'c' then symbolsize := symbolsize/2.0;
   end;
end;
(* end module dnaplot.getsymbol *)

(* begin module dnaplot.writepostscript *)
procedure writepostscript(var fin, fout: text; parameters: params;
                          var positions, dnasymbols: text);
(* read data from the fin file and write PostScript to the fout file.
Use positions to determine the range to plot if parameter readpositions
is true.  Use the dnasymbols to put marks on the graph *)
var
  currentpos: integer; (* the currently plotted position *)
  cVal: integer; (* the coordinate column in fin *)
  graphsperpage: integer; (* current count of graphs on this page *)
  inout: char; (* i = inside data, o = outside.  This allows the program
     to mark the edges of where the data lies. *)
  piecename: string; (* name of the piece *)
  newsequence: boolean; (* true just after we start a new sequence.
     This is used to trigger making a new set of axes *)
  pagenumber: integer; (* the current page *)
  sVal: integer; (* the sequence number column in fin *)
  vVal: real;    (* the value column in fin *)
  seqnum: integer; (* the current sequence number *)
  symbol: char; (* the current symbol: c(ircle), s(quare), t(riangle)
     and blank meaning that there are no further symbols *)
  symboltype: char; (* how to mark the current symbol, see makemarktype *)
  symbolplacement: char; (* how to place the current symbol, a(bsolute)
     on the graph or r(elative) to the current Ri value *)
  symbolbits: real; (* where to place the current symbol on the absolute
     or relative scale *)
  symbolsize: real; (* size of symbols relative to graph spacing *)
  symbolpiece: integer; (* the piece number on which to plot the next symbol *)
  symbolcoordinate: integer; (* the coordinate to plot the next symbol *)
  ygraphsize: real; (* the size into which the graphs fit vertically *)

procedure drawxaxis;
(* draw the x axis *)
begin with parameters do begin
  writeln(fout);
  writeln(fout,'% draw x axis');
  writeln(fout,'gsave % [ xaxis');
  writeln(fout,'0 ',bitlower/bitspercm *defscale:infofield:infodecim,
               ' translate');
  writeln(fout,'0 0 moveto');
  xaxis(fout,
        xaxlength*defscale,currentpos,xinterval,currentpos+numperln,
        xsubint,
        xticlength*defscale, xticdx*defscale, xticdy*defscale,
        xwidth, xdecimal,
        false, false, 2.0);
  (* now put the label out *)
  writeln(fout,(xaxlength /2.0 * defscale):infofield:infodecim,
               ((xticlength+3*xticdy)*defscale):infofield:infodecim,
               ' moveto');
  write(fout,'(');
  if sequencelabel then begin
     write(fout,'sequence: ',seqnum:1);
{zzz}
     write(fout,', ');
     writestring(fout,piecename);
     write(fout,', ');
  end;
  writestring(fout,xaxislabel); writeln(fout,')');
  writeln(fout,' dup stringwidth pop neg 2 div 0 rmoveto % center');
  writeln(fout,'show');
  writeln(fout,'grestore % ] xaxis');
  if plottype = 'b' then begin
     writeln(fout, '0 0 moveto');
     writeln(fout, xaxlength:infofield:infodecim, ' cm  0 cm lineto');
     writeln(fout, 'stroke');
  end;
end; end;

procedure drawyaxis;
(* since tic marks are drawn individually, the total
length is larger than the stack can handle and one gets
a limit check in ghostscript.  So put them out each time... *)
begin with parameters do begin
  (* begin plotting the axes *)
{
  writeln(fout, '0 0 moveto');
  writeln(fout, '0 cm yaxlength cm lineto');
  writeln(fout, 'stroke');
}
  (* move to zero of y axis then draw *)
  writeln(fout);
  writeln(fout,'% draw y axis');
  writeln(fout,'gsave % [ yaxis');
  writeln(fout,'0 ',bitlower/bitspercm *defscale:infofield:infodecim,
               ' translate');
  writeln(fout,'0 0 moveto');
  yaxis(fout,
        yaxlength*defscale,bitlower,yinterval,bitupper,
        ysubint,
        yticlength*defscale, yticdx*defscale, yticdy*defscale,
        ywidth, ydecimal,
        false, false, 2.0);
  (* now put the label out *)
  writeln(fout,((0)*defscale):infofield:infodecim,
               ((yaxlength+xticlength) * defscale):infofield:infodecim,
               ' moveto');
  write  (fout, '('); writestring(fout,yaxislabel); writeln(fout,') show');
  writeln(fout,'grestore % ] yaxis');
  writeln(fout);
(* end plotting the axes *)
end; end; (* drawyaxis *)

procedure stepcurrent;
(* step one position and
do a page or axes depending on the value of currentpos and graphsperpage *)
var
  dopage: boolean; (* if true, start a page *)
  doaxes: boolean; (* if true, draw the axes *)
begin with parameters do begin

  currentpos := succ(currentpos);
{
writeln(output,'stepcurrent currentpos=',currentpos:1);
}

  (* The last part of the logic below is that the very last base was done on
     the previous graph or page, so don't it do again *)

  doaxes := ((((currentpos-frompos) mod numperln) = 0)
             and (currentpos <> topos))
            or newsequence;

  (* only when we have completed a graph do we count it! *)
  if doaxes then graphsperpage := succ(graphsperpage);

  dopage := (graphsperpage >= numperpg) or (pagenumber = 0);

{
writeln(output,'graphsperpage = ',graphsperpage:1);
writeln(output,'pagenumber = ',pagenumber:1);
writeln(output,'currentpos=frompos = ',(currentpos=frompos):1);
writeln(output,'dopage = ',dopage);
writeln(output,'doaxes = ',doaxes);
writeln(output,'newsequence = ',newsequence);
}

  (* note: one cannot put the grestore and gsave inside the startpage and
     endpage functions:  for some reason they are ignored! *)
  if dopage then begin

     (* finish the last data segment *)
     if inout = 'i' then begin
        writeln(fout,'dodash {segmentmark} if');
        inout := 'o';
     end;

     pagenumber := succ(pagenumber);
     graphsperpage := 0;
     writeln(output, 'starting page ',pagenumber:1);
     if pagenumber > 1 then begin
        writeln(fout, 'showpage thepage restore');
     end;
     writeln(fout);
     writeln(fout, '%%Page: ',pagenumber:1,' ',pagenumber:1);
     writeln(fout, '/thepage save def startpage % ',(pagenumber):1);
     doaxes := true; (* new axes are needed on the new page *)
  end;

  if doaxes then begin

     (* finish the last data segment *)
     if inout = 'i' then begin
        writeln(fout,'dodash {segmentmark} if');
        inout := 'o';
     end;

     writeln(output, 'drawing axis @ ', currentpos:1);
     writeln(fout,' 0 ',-ygraphsize:infofield:infodecim,' cm translate');
     if showaxis then begin
        drawxaxis;
        drawyaxis;
     end;

     case plottype of (* start out at the right place *)
       'b': writeln(fout,'0 ',bitlower:infofield:infodecim,' bits moveto');
       'z': writeln(fout,'0 0 moveto');
     end;
     newsequence := false;
  end;

end; end; (* stepcurrent *)

procedure drawsymbols;
(* select the symbol, the way to draw it, its size, its position *)
begin
  if symboltype <> ' ' then begin
    if (currentpos = symbolcoordinate) and
       (sVal = symbolpiece) then begin
      writeln(fout);
      writeln(fout,'gsave % symbol @ ',currentpos:nfield);
      writeln(fout,'  currentpoint pop % x position');
(* This is a method to put the symbol BETWEEN lines by modifying x:
  writeln(fout,'spaceperbase 2 div neg add');
*)

      case symbolplacement of
        'a': write(fout, symbolbits:infofield:infodecim,' bits');
        'r': write(fout, (vVal+symbolbits):infofield:infodecim,' bits');
      end;
      writeln(fout, ' % y position');

      (* size is based on spacing between lines *)
      writeln(fout,symbolsize:infofield:infodecim,
                   ' spaceperbase mul % size');

      case symbol of
        'c': write(fout,'  circlesymbol');
        't': write(fout,'  trianglesymbol');
        's': write(fout,'  squaresymbol');
      end;

      makemarktype(fout, symboltype);

      writeln(fout,'grestore');

      (* move to next symbol *)
      getsymbol(dnasymbols, symbol, symboltype, symbolplacement,
                symbolbits, symbolsize, symbolpiece, symbolcoordinate);
    end;
  end;
end;

begin (* writepostscript *)
with parameters do begin

  writeln(fout, '% { % writepostscript');
  writeln(fout, '/Times-Roman findfont');
  writeln(fout, '12 scalefont');
  writeln(fout, 'setfont');
  writeln(fout);
  startpic(fout, 1.0, 0.0, 0.0,'t');

  (* begin setting scale factors *)
  writeln(fout, '/cmfactor 72 2.54 div def % defines points -> centimeters');
  writeln(fout, '/cm {cmfactor mul} def % defines centimeters');
  writeln(fout, '/bits {',  bitspercm:infofield:infodecim, ' div cm } def',
                '  % defines bits');
  writeln(fout, '/totalseq ', topos:infofield,  frompos:infofield, ' sub def',
                '  % the length of the total sequence');
  writeln(fout, '/spacing ', yaxlength:infofield:infodecim, numperpg:infofield,
                ' div cm def  % the space between lines');
  writeln(fout, '/numberperln ', numperln:infofield, ' def % bases per line');
  writeln(fout, '/spaceperbase ', xaxlength:infofield:infodecim, ' cm ',
                ' numperln div def % space between base pairs');
  writeln(fout, '/linewidthfactor spaceperbase thelinewidth mul def');

  writeln(fout);

  (* set up a page *)
  ygraphsize := yaxlength;
  (* + xticlength + 3*xticdy + 3*yticdy;*)
  if xticlength < 0 then ygraphsize := ygraphsize - 3*xticlength;
  if xticdy < 0 then ygraphsize := ygraphsize - 3*xticdy;
  if yticdy < 0 then ygraphsize := ygraphsize - 3*yticdy;

  writeln(fout, '/startpage {');
  writeln(fout,    orix:infofield:infodecim, ' cm',
               ' ',(oriy + numperpg*ygraphsize):infofield:infodecim,
               ' cm translate');
  writeln(fout, 'erasepage');
  writeln(fout,'0 0 moveto');
  writeln(fout, '} def');
  writeln(fout);

  (* The function d is named with a single character
     to save space and increase speed in the final PostScript file *)
  writeln(fout);
{writeln(fout,'/debug 1 def');}
  writeln(fout, '/d {% drawmarks');
{writeln(fout, 'debug 0 0 setrgbcolor');}
  writeln(fout, '/thebits exch def');
  writeln(fout, '/theheight thebits bits def');
  writeln(fout, 'currentpoint');
  writeln(fout, 'pop % dump y');
  writeln(fout, '/xcurrent exch def');
  writeln(fout, 'thebits bitlower ne {');
{thebits converted but bitlower not!! }
{writeln(fout, '0 debug 0 setrgbcolor');}
  writeln(fout, '  xcurrent theheight lineto');
  writeln(fout, '  linewidthfactor setlinewidth stroke');
  writeln(fout, '  1 setlinewidth');
  writeln(fout, '} if');

  case plottype of
     'z': begin
        writeln(fout, 'xcurrent spaceperbase add 0 cm moveto',
                      ' % drawmarks from zero');
     end;
     'b': begin
        writeln(fout, 'xcurrent spaceperbase add ',
                      bitlower/bitspercm:infofield:infodecim,
                      ' cm moveto % drawmarks from bottom');
     end;
  end;
{writeln(fout, '0 0 0 setrgbcolor');(* debug *)}
  writeln(fout, '} def');
  writeln(fout);

  writeln(fout,'/segmentmark { % mark at the contiguous segment of plot');
  writeln(fout,'gsave % [');
{  writeln(fout,'1 0 0 setrgbcolor '); mark red for testing}
  writeln(fout,'currentpoint pop 0 translate % use x position but not y');
  writeln(fout,'0 ',bitlower/bitspercm:infofield:infodecim, ' cm moveto');
  writeln(fout,'0 ',bitupper/bitspercm:infofield:infodecim, ' cm lineto');
  writeln(fout,'[2 4] 0 setdash % DO DASH'); (* turn dash on *)
  writeln(fout,'stroke % DO DASH');
  writeln(fout,'[] 0 setdash % DO DASH'); (* turn dash off *)
  writeln(fout,'grestore % ]');
  writeln(fout,'} bind def');
  writeln(fout);

  writeln(fout,'/presegmentmark { % mark before the contiguous segment of plot');
  writeln(fout,'gsave % [');
  writeln(fout,'  currentpoint pop % x position');
  writeln(fout,'  spaceperbase sub % new x position');
  writeln(fout,'  dup 0 gt { % don''t mark left of axis');
  writeln(fout,'    0 translate');
  writeln(fout,'    0 0 moveto');
  writeln(fout,'    segmentmark');
  writeln(fout,'  } if');
  writeln(fout,'grestore % ]');
  writeln(fout,'} bind def');
  writeln(fout);

  (* define symbols and their use.  Code taken from makelogo 8.07 *)
  writeln(fout,'/circlesymbol { % x y radius circlesymbol - (path)');
  writeln(fout,'newpath 0 360 arc');
  writeln(fout,'closepath} bind def');
  writeln(fout);

  writeln(fout,'/sqrt3 3 sqrt def');
  writeln(fout,'/trianglesymbol { % x y radius trianglesymbol - (path)');
  writeln(fout,'/r exch def');
  writeln(fout,'/sqrt3r sqrt3 r mul def');
  writeln(fout,'translate');
  writeln(fout,'% firstaxis');
  writeln(fout,'120 rotate');
  writeln(fout,'% secondaxis');
  writeln(fout,'0 r translate');
  writeln(fout,'-120 rotate');
  writeln(fout,'newpath');
  writeln(fout,'0 0 moveto');
  writeln(fout,'sqrt3r 0 lineto');
  writeln(fout,'-300 rotate');
  writeln(fout,'% thirdaxis');
  writeln(fout,'sqrt3r 0 lineto');
  writeln(fout,'closepath} bind def');
  writeln(fout);
 
  writeln(fout,'/squaresymbol { % x y side squaresymbol - (path)');
  writeln(fout,'/side exch def');
  writeln(fout,'translate');
  writeln(fout,'side 2 div neg dup translate');
  writeln(fout,'newpath');
  writeln(fout,'0 0 moveto');
  writeln(fout,'0 side lineto');
  writeln(fout,'side side lineto');
  writeln(fout,'side 0 lineto');
  writeln(fout,'closepath} bind def');
  writeln(fout);

  writeln(fout, '%%EndProlog');

  writeln(fout);
  (**** Start the main loop through the data, making the graph ************)

  reset(positions);

  reset(dnasymbols);
  getsymbol(dnasymbols, symbol, symboltype, symbolplacement, symbolbits,
            symbolsize, symbolpiece, symbolcoordinate);

  pagenumber := 0;
  graphsperpage := 0;
  seqnum := -maxint; (* force new sequence function in main loop *)
  currentpos := frompos-1;
  newsequence := false;

  skipdata(fin);
  repeat
    grabcolumns(fin, sCol, pCol, cCol, vCol, sVal, piecename, cVal, vVal);
{
writeln(output,'grabbed:',
' sVal=',sVal:1,
' cVal=',cVal:1,
' vVal=',vVal:1);
}
 
     if seqnum <> sVal then begin
        newsequence := true;
        if readpositions then readln(positions, frompos, topos);
        seqnum := sVal;
        writeln(output);
        write(output,'Sequence ',seqnum:1,', ');
{zzz}
        writestring(output,piecename);
        writeln(output,', positions: ', frompos:1, ' to ',topos:1);
        currentpos := frompos-1;
        stepcurrent;
        inout := 'o';  (* we start outside the data *)
     end;

    if (cVal >= frompos) and (cVal <= topos) then begin
       (* skip over blanks and catch up to cVal *)
       (* If curpos is GREATER than cVal, then the sequence
          numbering is DECREASING and the direction of plotting
          geez... can this be handled?  Well, the graph coordinates
          cannot go negative!  Looks like this is a bug for now. *)
       if currentpos > cVal then begin
          writeln(output);
          writeln(output,'The current position',
                         ' currentpos = ', currentpos:1,
                         ' exceeds',
                         ' cVal = ', cVal:1,
                         ' at sequence ',seqnum:1,'.');
          write(output,'This program cannot handle negative');
          writeln(output,' coordinate systems because');
          writeln(output,'the axes can only be drawn positively.');
          write  (output,'This error may also have occured because');
          writeln(output,' the program cannot handle scans');
          writeln(output,'in both directions.');
          halt
       end;

       while currentpos < cVal do begin
          (* if inside the data, finish it with out mark *)
          if inout = 'i' then begin
             writeln(fout,'dodash {segmentmark} if');
             inout := 'o';
          end;

          drawsymbols;
          (* draw the data line here *)
          case plottype of
            'b': write(fout,bitlower:infofield:infodecim,' d'); (* drawmarks *)
            'z': write(fout,'0 d'); (* drawmarks *)
          end;
          writeln(fout, ' % @ ', currentpos:nfield);
          stepcurrent;
       end;

       if (currentpos <> cVal) then begin
          if not eof(fin) then begin
             writeln(output,'ERROR:',
                            ' currentpos = ', currentpos:1,
                            ' <> ',
                            ' cVal =', cVal:1,
                            ' at sequence ',seqnum:1);
             halt
          end
          else begin (* else we just hit end of data *)
             (* if inside the data, finish it with out mark *)
             if inout = 'i' then begin
                writeln(fout,'dodash {segmentmark} if');
                inout := 'o';
             end;
          end
       end
       else begin
          if (vVal < bitlower) then vVal := bitlower;

          (* if outside the data, start it with going in mark *)
          if inout = 'o' then begin
             writeln(fout,'dodash {presegmentmark} if');
             inout := 'i';
          end;

          drawsymbols;
          write(fout, vVal:infofield:infodecim,' d'); (* drawmarks *)
          writeln(fout, ' % @ ', cVal:nfield);
          stepcurrent;
       end;
    end;
  until eof(fin);

  (* final mark of all if needed *)
  if inout = 'i' then begin
     writeln(fout,'dodash {segmentmark} if');
     inout := 'o';
  end;

  writeln(fout, 'showpage');
  writeln(fout, '%%Trailer');

  if symbol <> ' ' then begin
     writeln(output,'WARNING: There are unused dnasymbols:');
     while symbol <> ' ' do begin
        writeln(output);
        writeln(output,'symbol:          ',symbol:infofield);
        writeln(output,'symboltype:      ',symboltype:infofield);
        writeln(output,'symbolplacement: ',symbolplacement:infofield);
        writeln(output,'symbolbits:      ',symbolbits:infofield:infodecim);
        writeln(output,'symbolsize:      ',symbolsize:infofield:infodecim);
        writeln(output,'symbolpiece:     ',symbolpiece:infofield);
        writeln(output,'symbolcoordinate:',symbolcoordinate:infofield);
        getsymbol(dnasymbols, symbol, symboltype, symbolplacement,
                  symbolbits, symbolsize, symbolpiece, symbolcoordinate);
     end;
  end;

end;
end;
(* end module dnaplot.writepostscript *)

(* begin module dnaplot.makeheader *)
procedure makeheader(var fin: text; parameters: params; var fout: text);
(* Reads the header lines from fin and writes them to fout. *)
var
  totalnumpages: integer; (* total number of pages that will be produced *)
begin (* makeheader *)

  with parameters do begin
    if not readpositions then begin
       (* determine how many pages the process will produce *)
       totalnumpages := round((((topos-frompos)/numperln)/numperpg));
       write(output, 'For each sequence, dnaplot will produce ',
                      totalnumpages:1, ' page');
       if totalnumpages <> 1 then write(output,'s');
       writeln(output,' of graphs.');
    end;

    (* do our PostScript duty *)
    rewrite(fout);
    writeln(fout,'%!PS-Adobe-2.0 EPSF-2.0');
    writeln(fout,'%%Title: dnaplot ',version:4:2);
    writeln(fout,'%%Creator: Thomas D. Schneider');
    writeln(fout,'%%EndComments');
    writeln(fout);
  end;

  reset(fin);
  repeat
     write(fout, '% ');
     copyaline(fin, fout); (* copy commented lines of fin to fout *)
   until fin^ <> '*';

  writeln(fout);
end;
(* end module dnaplot.makeheader *)

(* begin module dnaplot.themain *)
procedure themain(var fin, fout, finp, positions, dnasymbols: text);
(* the main procedure of the program *)
var
   parameters: params; (* parameters read from finp *)

begin (* themain *)
  writeln(output, 'dnaplot ', version:4:2);

  PatentMessage(output);

  limitdate ('9','9','9','9','    /99/99 99:99:99');
  readparam(finp, parameters);
  makeheader(fin, parameters, fout);
  writeparam(fout, parameters);
  writepostscript(fin, fout, parameters, positions, dnasymbols);
end;
(* end module dnaplot.themain *)

begin
   themain(dnain, dnaout, dnaplotp, positions, dnasymbols);
1: end.
