program rsgra(rsdata, picture, rsgrap, wave, marks, output);
(* rsequence graph

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

module libraries required: delman, prgmods, dops *)

label 1; (* end of program *) 
  
const 
(* begin module version *)
   version = 5.06; (* of rsgra.p 2007 Feb 9
2007 Feb  9, 5.06: prgmod upgrade
2007 Feb  9, 5.05: avoid rotation by 90 degrees
                   crash because eof fails on empty wave file.
2003 Nov 18, 5.04: adjust location of figure on page.
2003 Nov 18, 5.01: label position axis
1998 Jun 19: point to rseq in see also
1993 Aug 19: minor changes
1991 May 31: rsgrap file becomes wave file and rsgrap takes
range of the graph to show.
origin 1986 feb 19 *)
(* end module version *)
  
(* begin module describe.rsgra *)
(*
name
   rsgra: rsequence graph
  
synopsis
   rsgra(rsdata: in, picture: out, rsgrap, wave: in, marks: in,
         output: out) 

files 
   rsdata: data file from rseq program

   picture: graph of rsequence in PostScript

   rsgrap: parameters to control the program.
      first line: two integers that define the from-to range to display

   wave: parameters to control the program.  empty file means no cosine
      wave, otherwise the parameters of the wave are given one per line:
      extreme: char;  h or l, the extreme high or low point on the curve
         defined by the wavelocation and wavebit
      wavelocation: real; the location in bases of the extreme
      wavebit: real; the location in bits of the extreme
      waveamplitude: real; the amplitude of the wave in bits
      wavelength: real; the wave length of the wave in bases

   marks: an empty file or a set of integers, one per line that are the
      locations of bases that should be specially marked on the graph.
      If the first line of the file begins with the letter 'b' and is
      followed by a real number, then this number defines the location of
      a bar to be placed on the graph immediately after the position given.

   output: messages to the user
  
description 
   Rsgra generates a graph of Rsequence versus position l.
   The user may overlay a cosine wave on the graph, to signify
   the orientation of the major grove, for example.

   See the discussion about the REMOVE feature in makelogo.p.

see also
   rseq.p, makelogo.p
  
author
   Thomas D. Schneider 
  
bugs

   2007 Feb 9: Adobe Readers look at PostScript files and if the file
   contains a 90 rotate, the reader TURNS THE WHOLE PAGE.  Stupid.  To
   avoid this, break 90 degree rotations into two 45 degree rotations.
   Amazingly simple way to fool an idiotic program.

   But it doesn't work.  So the rotation is not done now.

   For the same reason, the marks file cannot be empty.
   Use a large negative value to avoid having marks.

*)
(* end module describe.rsgra *)

(* begin module rsgra.const *)
      nfield = 4; (* size of field for printing n, the number of sites *)
(* end module rsgra.const *)
  
(* begin module pic.const *)
   pi = 3.14159265354; (* circumference divided by diameter of circle *)
   picfield = 8; (* width of numbers printed to the file *)
   picwidth = 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 = 81; (* default scale factor.  coordinate units per inch *)
(* end module pic.const version = 2.59; (@ of dops.p 1990 December 5 *)
  
(* begin module interact.const *) 
(* begin module string.const *)
      maxstring = 2000; (* the maximum string *)
(* end module string.const version = 5.31; (@ of prgmod.p 2006 Oct 10 *)
(* end module interact.const version = 5.31; (@ of prgmod.p 2006 Oct 10 *)
  
(* begin module rsgra.filler.const *) 
      fillermax = 21; (* the size of the filler array for a string *) 
(* end module rsgra.filler.const *)

type  
(* begin module rsgra.type *)
(* begin module wave.type *)
      waveptr = ^waveparam; (* to link several wave definitions together *)
      waveparam = record (* parameters to define a cosine wave *)
         (* define a cosine wave: *)
         extreme: char; (* h or l, the high or low extreme to be defined *)
         wavelocation: real; (* the location in bases of the extreme *)
         wavebit: real; (* the location in bits of the extreme *)
         waveamplitude: real; (* the amplitude of the wave in bits *)
         wavelength: real; (* the wave length of the wave in bases *)
{ 2000 Nov 1.  This is difficult because the length of a curve is not
simply multiplied by a scale factor.
         BASEdashon: real; (* dashon interval in bases *)
         BASEdashoff: real; (* dashon interval in bases *)
         BASEdashoffset: real; (* dashon interval in bases *)
}
         dashon: real; (* the size of on dashes in cm.
                          dashon <= 0 means no dashes *)
         dashoff: real; (* the size of on dashes in cm *)
         dashoffset: real; (* the size of off dashes in cm *)
         thickness: real; (* thickness of the wave in cm.  <= 0 means default *)
         next: waveptr; (* the next wave *)
      end;
(* end module wave.type version = 5.31; (@ of prgmod.p 2006 Oct 10 *)
(* end module rsgra.type version = 5.31; (@ of prgmod.p 2006 Oct 10 *)

(* begin module rs.type *)
   rstype = record (* types of data in the data table from rseq *)
      rstart, rstop: integer; (* range of the data *)
      l, (* position *)
      nal,ncl,ngl,ntl: integer; (* numbers of each base *)
      rsl, (* rsequence(l) *)
      rs,  (* running sum of rs *)
      varhnb, (* variance estimate for small sample correction *)
      sumvar: (* sum of varhnb *)
         real;
      nl: integer; (* = a+c+g+t *)
      ehnb: real; (* hg-e(n) *)
      flag: char; (* e = exact, a = approximate sampling statistics *)
   end;
(* end module rs.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 *)
      next: stringptr; (* the next string in a series *)
   end;
(* end module string.type version = 5.31; (@ of prgmod.p 2006 Oct 10 *)
(* end module interact.type version = 5.31; (@ of prgmod.p 2006 Oct 10 *)

(* begin module trigger.type *) 
      trigger = record (* an object to be searched for *)
         seek: string; (* the characters looked for *)
         state: integer; (* how close to triggering we are *)
         skip: boolean; (* trigger not found- skip the line *)
         found: boolean (* the trigger was found *)
      end;
(* end module trigger.type version = 5.31; (@ of prgmod.p 2006 Oct 10 *)

(* begin module filler.type *)
      (* the following is an array used to fill a string.
      it is convenient to have it much shorter than the maxstring, so that
      it is easy to fill the string using procedure fillstring.
      the user must declare the value of constant fillermax. *)
      filler = packed array[1..fillermax] of char;
(* end module filler.type version = 5.31; (@ of prgmod.p 2006 Oct 10 *)
  
var 
(* begin module rsgra.var *)
   rsdata, (* output of rseq program *)
   picture, (* output of this program, graph in PostScript *)
   rsgrap, (* parameters to control the program *)
   wave, (* parameters to control the drawing of the cosine wave *)
   marks: (* locations of marks *)
      text;
(* end module rsgra.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.59; (@ of dops.p 1990 December 5 *)

(* 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.31; (@ of prgmod.p 2006 Oct 10 *)
  
(* 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 *)
begin (* clearstring *)
   with ribbon do begin
      emptystring(ribbon);
      next := nil;
   end
end; (* clearstring *)
(* ---- *)
procedure initializestring(var ribbon: string);
(* start the string with a nil pointer.   This routine
should be called before doing linked list work.  This allows
the standard string routines to clear the string without
killing the pointer.  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.31; (@ of prgmod.p 2006 Oct 10 *)
(* end module interact.clearstring version = 5.31; (@ of prgmod.p 2006 Oct 10 *)

(* 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.31; (@ of prgmod.p 2006 Oct 10 *)
(* end module interact.writestring version = 5.31; (@ of prgmod.p 2006 Oct 10 *)
  
(* begin module trigger.proc *) 
(* this module allows one to scan a series of characters, as from
an array or a file, and to "trigger" or detect a simple string
in the series.  the advantage of the trigger is that several triggers
can "observe" a stream of characters at once, each looking for a
different thing.
some other modules required: interact.const, interact.type *)

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

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

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

      for index := 1 to length do s.letters[index] := a[index];
      s.length := length;
      s.current := 1
end; (* fillstring *)
(* end module filler.fillstring version = 5.31; (@ of prgmod.p 2006 Oct 10 *)

(* begin module filler.filltrigger *) 
procedure filltrigger(var t: trigger;
                      a: filler);
(* fill the trigger t *)
begin (* filltrigger *)
         fillstring(t.seek,a)
end; (* fillstring *)
(* end module filler.filltrigger version = 5.31; (@ of prgmod.p 2006 Oct 10 *)

(* 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.31; (@ of prgmod.p 2006 Oct 10 *)

(* begin module pic.functions *)
(* ********************************************************************** *)
(* begin module pic.await *)
procedure await;
(* Wait for user to type a carriage return.  the routine assumes that
there is a global file called input. *)
begin
(* the infinite way: 
    writeln(output);
    writeln(output,'*********************************');
    writeln(output,'* Use control-c to kill program *');
    writeln(output,'*********************************');
    while true do begin end;*)
end;
(* end module pic.await version = 2.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.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
   writeln(afile,'gsave'); (* save the current graphics state *)

{
Only use initgraphics if really necessary.  Since it restarts
the graphics, one can't put the figure inside another
PostScript figure.
   writeln(afile,'initgraphics'); (* make sure the printer
           is ready to print, without this, sometimes an Apple laserwriter
           will print the graph upside down, tiny and backwards! *)
}
   scale := setscale; (* set the global scale *)
   (* clean residue from before *)
   writeln(afile,'clear erasepage % REMOVE FOR PACKAGING INTO ANOTHER FIGURE');

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

(* If the following statement is done then it will work on the sun, but will
kill the applewriter!!!! 

Sun's non-standard PostScript extension, setlinewidth has default 0,
as stated in the Read This First and the NeWS Manual.
This draws very quickly with 1 bit wide lines.
If '1 setlinequality' is not done, then one cannot set the width of lines.
So to use PostScript on the screen, I must first do
'1 setlinequality'.  However, if I send this code to the Applewriter,
it kills PostScript on the Applewriter and I get no output whatsoever!
(It took me several hours to figure this out, since once PostScript is killed
on the Applewriter, the NEXT output is also smashed and I had to figure that
out also...)

So a standard PostScript program will not work correctly with the
default.  "Correcting" the PostScript program so that it works on the
Sun means that it BOMBS on the Applewriter.

The default for setlinequality should be '1 setlinequality' so that
the same PostScript code can be used both on the Sun and with other devices.
If you want speed, use the nonstandard form.
An alternative is to redefine setlinequality so that '0 setlinequality'
does give correct results with standard PostScript.

Please review this, Randy.  I think that Sun should fix it.
   writeln(afile,'1 setlinequality');  makes lines at least 1 bit wide *)

   (* set the scale to inches
   writeln(afile,
                 scale:picfield:picwidth,' ',
                 scale:picfield:picwidth,' 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:picfield:picwidth,' mul} def');
*)

   (* move to the start point on the page *)
   writeln(afile,
                     (x*scale):picfield:picwidth,
                 ' ',(y*scale):picfield:picwidth,
                 ' 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:picfield:picwidth);*)
end;
(* end module pic.startpic version = 2.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.stoppic *)
procedure stoppic(var afile:text);
(* stop pic output to file afile *)
(* NONSTANDARD *)
begin
    if inpath then begin
       writeln(afile,'stroke');
       inpath := false
    end;

    writeln(afile,'showpage % REMOVE FOR PACKAGING INTO ANOTHER FIGURE');
    writeln(afile,'grestore'); (* restore the current graphics state
                                  to what it was before the startpic *)

    await;

    inpicture := false;
end;
(* end module pic.stoppic version = 2.59; (@ of dops.p 1990 December 5 *)

(* 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:picfield,
               ' ',(dy*scale):picwidth:picfield);
   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:picfield,
                                     ' ',(ddy):picwidth:picfield);
                         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:picfield,' 0 rl');
                         write(afile,' ',
                                    -ss:picwidth:picfield,' 0 rl');
                         write(afile,' ',(ddx):picwidth:picfield,
                                     ' ',(ddy):picwidth:picfield);
                         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.59; (@ of dops.p 1990 December 5 *)

(* 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.59; (@ of dops.p 1990 December 5 *)

(* 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.59; (@ of dops.p 1990 December 5 *)


(* 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.59; (@ of dops.p 1990 December 5 *)

(* 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.59; (@ of dops.p 1990 December 5 *)

(* 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.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.graphstring *)
procedure graphstring(var tofile: text; var s: string; centered: boolean); 
(* graph the string s.  If it is recognized as a quoted string (surrounded
by double quotes), graph it without the quotes and center it.
Always center if centered is true.
Otherwise simply graph it.  if not in picture, just write it to output *)
(* NONSTANDARD *)
var i: integer; (* index to s, and temporary storage *)
    mv: real; (* holds amount to move, in plotting coordinates *)
    quoted: boolean; (* true if the string is quoted *)
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 centered := true;

      if centered then begin
         (* generate the calls to center the string.  Note: this must not
            be done for the pic program, which already centers *)
         if quoted then i := length-2
                   else i := length;
         mv := i*charwidth/2.0;
         mover(tofile,-mv,0.0);
      end;

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

      write(tofile,'('); (* begin postscript literal *)
      if quoted
      then for i := 2 to length-1 do write(tofile,letters[i])
      else for i := 1 to length   do write(tofile,letters[i]);
      write(tofile,')'); (* end postscript literal *)
      writeln(tofile,' x'); (* show the literal *)
      inpath := false; (* force new path from here *)

      if centered then begin
         (* restore to previous location *)
         mover(tofile,mv,0.0);
      end;
   end
   else begin
(* no output if not in picture 
      writestring(tofile,s);
      writeln(tofile) *)
   end
end;
(* end module pic.graphstring version = 2.59; (@ of dops.p 1990 December 5 *)

(* 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.59; (@ of dops.p 1990 December 5 *)

(* 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.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.picnumber *)
procedure picnumber(var afile: text;
                    dx, dy, number: real; width, decimal: integer;
                    centered: boolean);
(* 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 centered *)
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, centered);
      mover(afile,-dx,-dy);
   end
end;
(* end module pic.picnumber version = 2.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.xtic *)
procedure xtic(var afile: text; length, dx, dy, number: real;
               width, decimal: integer);
(* 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. *)
begin
   liner(afile,0.0,-length);
   picnumber(afile,dx,dy,number,width,decimal,true);
   mover(afile,0.0,length);
end;
(* end module pic.xtic version = 2.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.ytic *)
procedure ytic(var afile: text; length, dx, dy, number: real;
               width, decimal: integer);
(* produce a tic mark for the y 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. *)
begin
   liner(afile,-length,0.0);
   picnumber(afile,dx,dy,number,width,decimal,true);
   mover(afile,length,0.0);
end;
(* end module pic.ytic version = 2.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.xaxis *)
procedure xaxis(var afile: text;
         axlength,fromtic,interval,totic: real;
         length, dx, dy: real;
                width, decimal: integer);
(* draw an x axis starting from the current position.
the length of the xaxis is axlength.
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 xtic.
If the width is zero, no number is produced.
the location after the call is the same as before the call. *)
var
   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 *)
   totici2: real; (* an extra half interval guarantees the last tic is placed *)
begin
   liner(afile,axlength,0.0);
   mover(afile,-axlength,0.0);
   if totic = fromtic then begin
      writeln(output,'xaxis: fromtic and totic cannot be equal');
      halt;
   end;
   if (axlength = 0.0) or (interval = 0.0) then begin
      writeln(output,'xaxis: neither axlength nor interval can be zero');
      halt;
   end;
   jump := axlength * interval / (totic - fromtic);
   jumpdistance := 0;

   (* note!  the extra amount of interval/2 makes sure that the last tic
      is placed, if the computer does not add quite right and overshoots
      by a tich on the value of tic.  This much makes it surely work! *)
   tic := fromtic;
   if interval > 0.0 then begin
      totici2 := totic + interval/2.0;
      while tic <= totici2 do begin
         xtic(afile,length,dx,dy,tic,width,decimal);
         tic := tic + interval;
         if tic <= totici2 then begin
            mover(afile,jump,0.0);
            jumpdistance := jumpdistance + jump;
         end
      end
   end
   else if interval < 0.0 then begin
      totici2 := totic - interval/2.0;
      while tic >= totici2 do begin
         xtic(afile,length,dx,dy,tic,width,decimal);
         tic := tic + interval;
         if tic >= totici2 then begin
            mover(afile,jump,0.0);
            jumpdistance := jumpdistance + jump
         end
      end
   end;

   mover(afile,-jumpdistance,0.0)
end;
(* end module pic.xaxis version = 2.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.yaxis *)
procedure yaxis(var afile: text;
         aylength,fromtic,interval,totic: real;
         length, dx, dy: real;
                width, decimal: integer);
(* draw a y axis starting from the current position.
the length of the yaxis is aylength.
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. *)
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 *)
begin
   liner(afile,0.0,aylength);
   mover(afile,0.0,-aylength);
   if totic = fromtic then begin
      writeln(output,'yaxis: fromtic and totic cannot be equal');
      halt;
   end;
   if (aylength = 0.0) or (interval = 0.0) then begin
      writeln(output,'yaxis: neither aylength nor interval can be zero');
      halt;
   end;
   jump := aylength * interval / (totic - fromtic);
   jumpdistance := 0;
   half := interval / 2.0;

   tic := fromtic;
   if interval > 0.0 then while tic <= totic+half do begin
      ytic(afile,length,dx,dy,tic,width,decimal);
      tic := tic + interval;
      if tic <= totic+half then begin
         mover(afile,0.0,jump);
         jumpdistance := jumpdistance + jump
      end
   end
   else if interval < 0.0 then while tic >= totic-half do begin
      ytic(afile,length,dx,dy,tic,width,decimal);
      tic := tic + interval;
      if tic >= totic-half then begin
         mover(afile,0.0,jump);
         jumpdistance := jumpdistance + jump
      end
   end;

   mover(afile,0.0,-jumpdistance)
end;
(* end module pic.yaxis version = 2.59; (@ of dops.p 1990 December 5 *)
(* ********************************************************************** *)
(* end module pic.functions version = 2.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.degtorad *)
function degtorad(angle: real):real;
(* convert angle in degrees to radians *)
begin
   degtorad := (angle / 360) * 2 * pi 
end;
(* end module pic.degtorad version = 2.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.polrec *)
procedure polrec(r,theta: real; var x,y: real);
(* convert polar to rectangular coordinates,
theta is in radians *)
begin
   x := r*cos(theta);
   y := r*sin(theta)
end;
(* end module pic.polrec version = 2.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.arc *)
procedure arc(var thefile: text; angle1, angle2, radius: real;
              steps: integer);
(* create an arc in thefile going from angle1 to angle2 (degrees) in the
positive direction of angle, with the given radius.
use the given number of steps to make it.  return to the same position
as before the arc was drawn. *)
var
   dtheta: real; (* change in theta *)
(* s: integer; (@ index to the steps *)
   theta: real; (* current angle *)
   x,y: real; (* coordinates around starting point *)
(* zerox,zeroy: real; (@ starting location, center of curve *)
begin
(*
   zerox := picxglobal;
   zeroy := picyglobal;
*)
   theta := degtorad(angle1);
   dtheta := degtorad( abs(angle2-angle1)/steps );
   polrec(radius,theta,x,y);
(* can't do this for postscript arc:
   movea(thefile,zerox+x,zeroy+y);
      ' ',scale*zerox:  picfield:picwidth,
      ' ',scale*zeroy:  picfield:picwidth,
*)

   (* NONSTANDARD postscript, much faster *)
   writeln(thefile,
   (* 'stroke newpath', *)(* force there to be no current point *)
   (* ' 0 0', *)
      'a',
      ' ',scale*radius: picfield:picwidth,
      ' ',angle1:       picfield:picwidth,
      ' ',angle2:       picfield:picwidth);
   write(thefile, 'arc');
   if angle2 < angle1 then write(thefile,'n'); (* for negative draws *)
   (* origin move: writeln(thefile,' stroke newpath 0 0 moveto'); *)
   writeln(thefile,' n');
   (* the moveto puts us back to the origin *)

(*
   for s := 1 to steps do begin
      theta := theta + dtheta;
      polrec(radius,theta, x,y);
      linea(thefile,zerox+x,zeroy+y);
   end;

   movea(thefile,zerox,zeroy)
*)
end;
(* end module pic.arc version = 2.59; (@ of dops.p 1990 December 5 *)

(* begin module pic.circler *)
procedure circler(var afile: text; radius: real);
(* make a circle at the current position of some radius. *)
var
   steps: integer; (* number of steps to make the circle *)
begin
   (* number of segments increases with diameter, but the constant
      still should be a function of how good it looks on a particular
      graphic system, I'm afraid.  However, there should be a lower
      bound on the number of steps, so even small circles look good *)
   if radius < 1.0 then steps := 25
                   else steps := round(radius*25);

   arc(afile,0.0,360.0,radius,steps);
end;
(* end module pic.circler version = 2.59; (@ of dops.p 1990 December 5 *)

(* ********************************************************************** *)
  
(* 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.31; (@ of prgmod.p 2006 Oct 10 *)
(* ********************************************************************** *)

{ 2007 Feb 9 Block original routine which fails because the eof test
fails (an empty afile is not detected as empty, then the read fails
because it is empty).
(* begin module rsgra.readwaveparameters *)
procedure readwaveparameters(var afile: text; var p: waveparam);
(* read from afile the parameters p, as defined by type param *)
begin
writeln(output,'BUBBA 0');
reset(afile);
      with p do begin
         if eof(afile)
         then wave := false
         else begin
            wave := true;
writeln(output,'BUBBA 1');
            readln(afile,extreme);
            readln(afile,wavelocation);
            readln(afile,wavebit);
            readln(afile,waveamplitude);
            readln(afile,wavelength);
         end
      end
end;
(* end module rsgra.readwaveparameters *)
}

(* begin module readwaveparameters *)
procedure readawaveparameter(var afile: text; var wp: waveparam);
{
                                 cmperbase, cmperbit: real);
}
(* Read one wave parameter from a file. *)
{
2000 Nov 1:  The following experimental code is not implemented:
See the coscurve program
(http://www.ccrnp.ncifcrf.gov/~toms/delila/coscurve.html) for details
about the conversion from bases to cm.  The code now requires
cm per base and cmperbit.
const
   lover2pi = 1.2160067233; (* L/2pi for L being the length of
       the cosine function from 0 to 2pi *)
var
   multiplier: real; (* factor to multiply by to convert from bases to cm. *)
}
begin
   readln(afile,wp.extreme);
   readln(afile,wp.wavelocation);
   readln(afile,wp.wavebit);
   readln(afile,wp.waveamplitude);
   readln(afile,wp.wavelength);
   if wp.wavelength <= 0.0 then begin
      writeln(output,'wave parameters: wavelength must be positive');
      halt
   end;
   (* one may skip the dash and thickness parameters
      of the last wave definition.  This allows
      backwards compatability, but should not be
      used in general. *)
   if eof(afile)
   then begin
      wp.dashon := 0;
      wp.dashoff := 0;
      wp.dashoffset := 0;
      wp.thickness := 0
   end
   else begin
      (* implement new method to read dash information
         without breaking previous wave parameter files *)
      if afile^ <> 'd'
      then begin (* default to old method *)
         readln(afile,wp.dashon);
         wp.dashoff := wp.dashon;
         wp.dashoffset := 0;
      end
      else begin
         get(afile); (* skip the 'd' *)
         (* reading the parameters directly in cm would require
            that the user compute or guess the values.
            So read in as bases and compute cm.
            Unfortunately the length of a curve is not the simple
            scale factor.  Make the user guess for now.
         *)
         readln(afile,wp.dashon, wp.dashoff, wp.dashoffset)

{
experimental code to account for actual wave length:
         readln(afile, wp.BASEdashon, wp.BASEdashoff, wp.BASEdashoffset);
         multiplier    := 1;
         multiplier    := lover2pi
 * wp.waveamplitude (* in bits *)
 * cmperbit (* converts to cm *)
 / wp.wavelength (* 10.6 bp/turn *)
* cmperbase
;
multiplier    := lover2pi;
multiplier    := 1.0;
multiplier    := 1.0/lover2pi;
         multiplier := lover2pi
                       * wp.waveamplitude (* in bits *)
                       * cmperbit (* converts to cm *)
                       * wp.wavelength (* 10.6 bp/turn *)
                       * cmperbase;
         multiplier := 1.0
         wp.dashon     := multiplier * wp.BASEdashon;
         wp.dashoff    := multiplier * wp.BASEdashoff;
         wp.dashoffset := multiplier * wp.BASEdashoffset;
with wp do begin
   writeln(output, 'cmperbase  = ',cmperbase:5:2);
   writeln(output, 'cmperbit   = ',cmperbit:5:2);
   writeln(output, 'dashon     = ',dashon:5:2);
   writeln(output, 'dashoff    = ',dashoff:5:2);
   writeln(output, 'dashoffset = ',dashoffset:5:2);
   writeln(output, 'waveamplitude   = ',wp.waveamplitude:5:2);
   writeln(output, 'wavelength      = ',wp.wavelength:5:2);
end;
}

      end;
      if eof(afile)
      then wp.thickness := 0
      else readln(afile,wp.thickness);
   end;
end;

procedure readwaveparameters(var afile: text; var w: waveptr);
{
                                 cmperbase, cmperbit: real);
}
(* read from afile the wave parameters w.  The routine is done either at end
of file, when a period ('.') appears or at a blank line.  This makes it
compatable with different reading mechanisms. *)
var
   done: boolean; (* done reading *)
   morethanone: boolean; (* more than one wave? *)
   p: waveptr; (* one of the definitions of a cosine wave *)
procedure waystoend;
(* if end of file, a blank line or a period, consider the job done *)
var
   clear: boolean; (* not eof, so clear to test for comment *)
begin
   if eof(afile) then done := true;
   clear := false;
   if not done (* skip any number of comments *)
   then while (not clear) and (not done) do begin
      if eof(afile) then begin
         done := true;
         clear := true;
      end
      else if afile^ = '*'
           then readln(afile)
           else clear := true
   end;
   if not done then if eoln(afile) then done := true;
   if not done then if afile^='.' then begin
      readln(afile);
      done := true;
   end;
end;
begin
   done := false;
   morethanone := false;
   w := nil;
   waystoend;
   if not done then begin
      new(w);
      p := w;
      while not done do begin
         waystoend;
         if done
         then p^.next := nil
         else begin
            if afile^='*'
            then readln(afile) (* skip any number of comments *)
            else begin
               if morethanone then begin
                  new(p^.next);
                  p := p^.next
               end;
{
               readawaveparameter(afile, p^, cmperbase, cmperbit);
}
               readawaveparameter(afile, p^);
               (* if we ever come back again there will be more than one: *)
               morethanone := true;
            end
         end
      end
   end
   else w := nil;
end;
(* end module readwaveparameters version = 5.31; (@ of prgmod.p 2006 Oct 10 *)

(* begin module rsgra.header *)
procedure header(var infile, marks, outfile: text);
(* do the header of the plot *)
var
      index: integer; (* to lines in the infile *) 
begin
      reset(infile);
      reset(marks);
      rewrite(outfile);

(*
      writeln(outfile,'.nf');
      writeln(outfile,'% rsgra ',version:4:2);
*)

      writeln(outfile, '%!PS-Adobe-2.0 EPSF-2.0');
      writeln(outfile, '%%Title: rsgra ',version:4:2);
      writeln(outfile, '%%Creator: Tom Schneider');
      writeln(outfile, '%%BoundingBox:    0   0   630 900');
      writeln(outfile, '%%DocumentFonts:');
      writeln(outfile, '%%EndComments');
      writeln(outfile, '%%EndProlog');

      (* copy header lines to outfile *)
      for index := 1 to 3 do begin
         write(outfile,'% ');
         copyaline(infile,outfile);
      end;

      writeln(outfile);

      (* identify the marks *)
      writeln(outfile,'% MARKS: ');
      while not eof(marks) do begin
         write(outfile,'% ');
         copyaline(marks,outfile);
      end;

      reset(marks);
end;
(* end module rsgra.header *)

(* begin module rsgra.cosine *)
procedure cosine(var afile: text;
                 amplitude, phase, wavelength, base,
                 xmin,ymin,xmax,ymax: real);
(* draw a cosine function down the page:
         x := amplitude*cos( ((-y-phase)/wavelength)*2*pi) + base,
where the given parameters are in inches.
parts outside the box defined by xmin,ymin,xmax,ymax are not drawn.
pi is defined in pic.const *)
const
      step = 0.01; (* step size for the lines drawn *)
var
      inside1: boolean; (* point 1 was inside the bounds *)
      inside2: boolean; (* point 2 was inside the bounds *)
      x1,y1,x2,y2: real; (* internal coordinates *)
function c(y: real): real;
(* do the actual calculation of the cosine given y *)
begin
      c := amplitude*cos( ((-y-phase)/wavelength)*2*pi) + base
end;
function between(a,b,c: real):boolean;
(* is b between a and c? *)
begin between:=(a<=b) and (b<=c) end;
begin
      (* start graphing at the top of the plot *)
      y1 := ymax;
      x1 := c(y1);
      inside1 := between(ymin,y1,ymax) and between(xmin,x1,xmax);

      while y1 >= ymin do begin
         x2 := x1;
         y2 := y1;
         inside2 := inside1;

         y1 := y1 - step;
         x1 := c(y1);
         inside1 := between(ymin,y1,ymax) and between(xmin,x1,xmax);

         if inside1 and inside2 then begin
               movea(afile,x1,y1);
               linea(afile,x2,y2);
         end
      end
end;
(* end module rsgra.cosine *)

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

(* begin module rs.readrsrange *)
procedure readrsrange(var rsdata: text; var r: rstype);
(* read the range data from rsdata to r.  data is assumed to
be the rsdata file from program rseq. *)
var
   index: integer; (* for counting lines of rsdata *)
   skip: char; (* a character to skip the '*' on the line *)
begin
      for index := 1 to 11 do readln(rsdata);
      readln(rsdata, skip, r.rstart, r.rstop);
    (*  writeln(output, 'range: ',r.rstart:1,' ',r.rstop:1);*)
end;
(* end module rs.readrsrange *)
  
(* begin module rs.getrsbegin *)
procedure getrsbegin(var infile: text);
(* skip to the beginning of the data in a data file from
the rseq program *)
var
      ch: char; (* a character read from infile *)  
      begindata: trigger; (* a trigger to locate the beginning of the data *)
begin
(*                                       1         2  *) 
(*                              123456789012345678901 *)
filltrigger(begindata         ,'l    a    c    g    t');   
      resettrigger(begindata);  

      reset(infile);
      while not begindata.found do begin
         if eoln(infile) then readln(infile);
         if eof(infile) then begin
            writeln(output,'beginning of data not found');
            halt;
         end;
         read(infile,ch);
         testfortrigger(ch,begindata); 
      end;
      readln(infile);
end;
(* end module rs.getrsbegin *)

(* begin module rs.readrsdata *)
procedure readrsdata(var rsdata: text; var rdata: rstype);
(* read data from the data file of program rseq into the datatype *)
begin
   with rdata do begin
      read(rsdata,l,nal,ncl,ngl,ntl,rsl,rs,varhnb,sumvar,nl,ehnb);
      (* skip spaces to find the flag: *)
      while rsdata^=' ' do get(rsdata);
      readln(rsdata,flag);
(* writeln(output,'readrsdata: l a c g t flag = ',
l:1,' ',nal:1,' ',ncl:1,' ',ngl:1,' ',ntl:1,' ',flag); *)
   end
end;
(* end module rs.readrsdata *)
  
(* ********************************************************************** *)
(* ********************************************************************** *)

(* begin module themain *)
procedure themain(var infile, rsgrap, wave, marks, outfile: text); 
(* the main procedure of the program *) 
const
   duplicates = 2;  (* number of times to draw the main line *)
   halfdatawidth = 2.4; (* inches to shift the data to the left
      from the left edge of the graph.
      the lowpoint is added to this to find the dataposition *)
   leftbound = 3.0; (* distance of the left bound of the graph
      from the left side of the paper, (not including the data) *)
   spacing = 0.01;  (* spacing between the multiple drawings of
                       the main line *)
   xscale = 2.00;   (* inches per bit across the plot *)
   yscale = 0.1595;   (* vertical spacing between graph lines,
      ie, inches per base *)
      (* 0.15 inches worked nicely for a while, until I wanted the graph
         to match the standard printer output (unix pr).
         There were 39 bases in the range, and they should have been
         2.3 bases longer, so the new yscale is 0.15 * (39+2.3)/39. *)
var 
   bar: real; (* the size of an error bar *)
   barposition: integer; (* position after which to write a bar *)
   basetomark: integer; (* the base to be marked according to the marks file *)
   dataposition: real; (* amount to shift the data to the left from the
      zero of the graph *)
   dx,dy: real; (* changes of x and y to make multiple lines *)
   frombase, tobase: integer; (* range of the graph to display *)
   line: integer; (* an index to the multiple lines drawn *)
   linedx,linedy: real; (* line * dx or dy, the actual jump made for 
      multiple lines (done for speed) *)
   lowpoint: real; (* the lowest point in the entire plot *)
   notfirstpoint: boolean; (* true for all points except
      the first point plotted *)
{
   params: waveparam; (* parameters to control the program *)
}
   params: waveptr; (* wave parameters to control the program *)
   point: real; (* a low point in the graph *)
   position: integer; (* a location in the aligned sequence *)
   rdata: rstype; (* data from rseq *)
   xmin,ymin,xmax,ymax: real; (* the bounds of the graph, in inches *)
   x1,y1: real; (* the current end of the curve *)
   x2,y2: real; (* the previous end of the curve *)
   ybar: real; (* location of the bar in graphics coordinates *)
begin
      writeln(output,'rsgra ',version:4:2);

      (* copy header stuff *)
      header(infile,marks,outfile);

      reset(wave);
      readwaveparameters(wave,params);

      reset(rsgrap);
      if eof(rsgrap) then begin
         frombase := -maxint;
         tobase := maxint;
      end
      else readln(rsgrap,frombase,tobase);

      (* now read the rest of the parameters *)
      if eof(marks)
      then begin
         barposition := -maxint (* don't do the bar *)
      end
      else begin
         if marks^ = 'b' then begin
            get(marks); (* skip the 'b' *)
            readln(marks,barposition);
            writeln(output,'bar placed at position ',barposition:1);
         end;
      end;

      (* get the first base to mark *)
      if not eof(marks) then readln(marks,basetomark)
                        else basetomark := -maxint;

      (* find the range of the graph in bases *)
      reset(infile);
      readrsrange(infile,rdata);
{
      if params.wave
}
      if params <> nil
      then writeln(outfile,'% cosine wavelength is ',
                            params^.wavelength:4:2, ' bases');
      with rdata do begin
         if frombase < rstart then frombase := rstart;
         if tobase > rstop then tobase := rstop;
      end;

      (* find lowest bound of the graph, in bits *)
      getrsbegin(infile);
      lowpoint := maxint;
      with rdata do for position := rstart to rstop do begin
         (* skip lines with an '*' *)
         if infile^ <> '*' then begin
            readrsdata(infile,rdata);
            point := rsl - sqrt(varhnb);
            if point < lowpoint then lowpoint := point
         end
         else readln(infile)
      end;
      (* make lowpoint be just below zero, or the next tenth of
         a bit lower: *)
      if lowpoint >= 0.0 then lowpoint := -0.1
                         else lowpoint := -trunc( (-10*lowpoint) + 1)/10;
      dataposition := -lowpoint*xscale + halfdatawidth;

      (* define the bounds of the graph region, in inches *)
      with rdata do begin
         xmin := lowpoint * xscale;
         xmax := 2.0 * xscale;
         ymin := -(tobase-frombase+2)*yscale;
         ymax := 0.0;
      end;

{
      startpic(outfile,defscale,3.3,8.5,'c');
}

      startpic(outfile,defscale,1.0,5.0,'c');

writeln(outfile,'45 rotate');
writeln(outfile,'45 rotate');
      (* the coordinate (0,0) is the upper left hand corner of the graph,
         not including the data.  First create the location of the data: *)
      movea(outfile,-leftbound,0.0);
      movea(outfile,-dataposition,0.0);

      writeln(outfile,'gsave 1 2 scale');
      writeln(outfile,'(   l    a    c    g    t  Rs(l)) x');
      writeln(outfile,'grestore');

      (* placement of the label in bits along the Y axis: *)
      movea(outfile,0.25*xscale,3*yscale);
      writeln(outfile,'gsave 1 2 scale');
      writeln(outfile,'(Rs(l) = Rsequence(l), Information in bits) x');
      writeln(outfile,'grestore');

      (* placement of the label in bits along the X axis: *)
(*
      movea(outfile,-1.49*xscale,-8.0*yscale);
zzz
*)
      movea(outfile,
            -1.60*xscale,
            (
            (frombase-tobase)/2 (* the center *)
            + 8.5 (* adjust for size of the words, in cm *)
            )*yscale
            );
      writeln(outfile,'gsave 2 2 scale');

(*
   2007 Feb 9: Adobe Readers look at PostScript files and if the file
   contains a 90 rotate, the reader TURNS THE WHOLE PAGE.  Stupid.  To
   avoid this, break 90 degree rotations into two 45 degree rotations.
   Amazingly simple way to fool an idiotic program.
      writeln(outfile,'-90 rotate');
*)
(* that didn't work.  Don't rotate figure
      writeln(outfile,'-45 rotate');
      writeln(outfile,'-45 rotate');
*)
      writeln(outfile,'-45 rotate');
      writeln(outfile,'-45 rotate');

      writeln(outfile,'(Position L (in bases)) x');
      writeln(outfile,'grestore');

      (* now put the book id onto the graph *)
      (* x coordinate was xscale-dataposition,
         but this way more title will fit on: *)
      movea(outfile,(-2*dataposition)/2,4*yscale);
      write(outfile,'(');

      write(outfile,'rsgra ',version:4:2,' ');
      reset(infile); readln(infile); readln(infile);
      while not eoln(infile) do begin (* like copyaline but no writeln *)
         outfile^ := infile^;   put(outfile);   get(infile)
      end;
      writeln(outfile,') x % NOHEADER FOR PACKAGING INTO ANOTHER FIGURE');

      (* create upper x axis *)
      movea(outfile,0.0,0.0);
      xaxis(outfile, xmax, 0.0,1.0,2.0, -0.1,0.0,+0.1, 3,1);
      (* more ticks: *)
      mover(outfile,lowpoint*xscale,0.0);
      xaxis(outfile, (2.0-lowpoint)*xscale,
            lowpoint,0.1,2.0, -0.1,0.0,+0.0, 0,1);

      (* dashed line down the plot: *)
      movea(outfile,0.0,0.0);
      drawa(outfile,0.0,ymin,'-',0.03);

      (* create lower x axis *)
      xaxis(outfile, xmax, 0.0,1.0,2.0, +0.1,0.0,-0.1, 0,1);
      (* more ticks: *)
      mover(outfile,xmin,0.0);
      xaxis(outfile, (2.0-lowpoint)*xscale,
            lowpoint,0.1,2.0, +0.1,0.0,-0.0, 0,1);
      (* xaxis(outfile, axlength,fromtic,interval,totic: real;
              length, dx, dy: real;   width, decimal: integer); *)

      notfirstpoint := false; (* about to draw the first point, so there is
         no previous point *)

      (* prepare for reading the data for graphing *)
      getrsbegin(infile);

      (* create the curve of rsequence(l) *)
      for position := rdata.rstart to rdata.rstop do begin
         (* skip lines with an '*' *)
         if infile^ <> '*' then begin
            readrsdata(infile,rdata);

            if (position >= frombase) and (position <= tobase)
            then with rdata do begin

               (* move to the next data line and draw a bar *)
               if position = barposition then begin
                  (* the factor multiplying dataposition determines the
                     left end of the bar.  It is made to fit nicely. *)
                  (* the extra 0.5 makes sure that the position is between
                     bases, so it won't overwrite the marks *)
                  ybar := (frombase-position-1.5)*yscale;
(*
                  movea(outfile,-1.60*dataposition,ybar);
*)
                  movea(outfile,-dataposition,ybar);
                  linea(outfile, 2*xscale, ybar)
               end;

               (* move to the next data line and give the data *)
	       (* the extra amount added to the position places the numbers in
	          a better position *)
               movea(outfile,-dataposition,-(position-frombase+1.2)*yscale);
               writeln(outfile,'(',l:nfield,
                                   ' ',nal:nfield,' ',ncl:nfield,
                                   ' ',ngl:nfield,' ',ntl:nfield,
                                   rsl:7:2,') x');

               if notfirstpoint then begin
                  (* multiple draw line from previous point *)
                  x2 := x1;                (* set (x2,y2) to previous point *)
                  y2 := y1;
                  x1 := rsl*xscale;        (* update (x1,y1) *)
                  y1 := (frombase-position-1)*yscale;
   
                  (* the slope of the line from (x2,y2) to (x1,y1) is
                     m := (y1-y2)/(x1-x2);
                     the multiple lines must be drawn perpendicular to this,
                     m' := -1/m;
                     the arctan of this gives the angle to move the lines in.
                     conversion from polar to rectangular coordinates
                     completes the job.  note that y2-y1 will never be zero. *)
                  polrec(spacing,arctan((x1-x2)/(y2-y1)),dx,dy);
   
                  for line := -duplicates to duplicates do begin
                     linedx := line*dx;
                     linedy := line*dy;
                     movea(outfile, x2+linedx, y2+linedy);
                     linea(outfile, x1+linedx, y1+linedy)
                  end
               end
               else begin (* the first point does not have a line drawn to it *)
                  x1 := rsl*xscale;
                  y1 := (frombase-position-1)*yscale;
                  notfirstpoint := true;
               end;
   
               (* create error bar and circle *)
               (* move back to the location of the new data point: *)
               movea(outfile, x1, y1);
               bar := sqrt(varhnb)*xscale;
               mover(outfile,+bar,0.0);
               liner(outfile,-bar,0.0);
               circler(outfile,yscale/3);
               liner(outfile,-bar,0.0);

               (* create specially marked bases *)
               if basetomark = l then begin
                  (*movea(outfile,lowpoint+xscale*0.025,*)
                  movea(outfile,lowpoint,
                                -(position-frombase+1)*yscale);
                  circler(outfile,yscale/3);
                  writeln(output,'marking base ',l:1);
                  if not eof(marks) then readln(marks,basetomark)
                                    else basetomark := -maxint;
               end
            end
         end
         else readln(infile)
      end;
   
      (* overlay a cosine wave as specified by the parameters *)
      with rdata do begin
{
         if params.wave then with params do begin
}
         if params <> nil then with params^ do begin
            if extreme = 'h'
            then cosine(outfile,           (* high center point cosine wave *)
                    +xscale*waveamplitude,
                    yscale*(wavelocation-frombase+1), yscale*wavelength,
                    xscale*(wavebit-waveamplitude),
                    xmin,ymin,xmax,ymax)
            else cosine(outfile,           (* low center point cosine wave *)
                    -xscale*waveamplitude,
                    yscale*(wavelocation-frombase+1), yscale*wavelength,
                    xscale*(wavebit+waveamplitude),
                    xmin,ymin,xmax,ymax);
         end;
      end;
      stoppic(outfile);
end;
(* end module themain *)
 
begin
      themain(rsdata, rsgrap, wave, marks, picture);
1: end.