program dnag(bdna, dooin, output);
(* graphics of DNA

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

module libraries required: delman, domod *)

label 1; (* end of program *) 
  
const 
(* begin module version *)
version = 1.74; (* of dnag.p 1999 October 22
1.74 1999 October 22: add pointer to bdna
last previous change: 1993 January 26
origin 1986 feb 23 *)
(* end module version *)
  
(* begin module describe.dnag *)
(*
name
      dnag: graphics of DNA
  
synopsis
      dnag(bdna: in, dooin: out, output: out) 
  
files 
      bdna: b- form dna coordinates.  lines beginning with '*' are ignored.
         on each line following is the coordinate of one atom.
         the first character is the kind of group:
            * P = phosphate, D = deoxyribose
            * A = adenine, G = guanine, C = cytosine, T = thymine
         the next character is blank
         the next two characters are the atom and its number
         then the locations are given, separated by spaces:
            radius (angstrom) - angle (degree) - z axis (angstrom)
      dooin: graph of dna in doodle format
      output: messages to the user
  
description 
      dnag generates a graph of DNA.

documentation
      B-DNA Cylindrical Polar Coordinates from
      S. Arnott and D. W. L. Hukins
      Biochem. and Biophys. Res. Comm 47: 1504-1509 (1972)
      "Optimised Parameters for A-DNA and B-DNA"

      M. Karplus and R. N. Porter
      Atoms & Molecules
      Benjamin/Cummings Publishing Co., Menlo Park, Ca, 1970
      p. 204-7, crystal radii

see also
   dops.p, bdna

author
      Thomas D. Schneider 
  
bugs

      The location of the strings may not be centered exactly in the circles.
To make this easy to adjust, two fudge factors (fudgex and fudgey) are provided
as constants.

   The program's output is in the pic language for troff.  This can
   be converted to PostScript using the dops program.

*)
(* end module describe.dnag *)

(* begin module dnag.const *)
   iascale = 0.20;   (* inches per angstrom across the plot *)
   fudgex = -0.02; (* fudge factor for x, in angstrom *)
   fudgey = -0.04; (* fudge factor for y, in angstrom *)
(* end module dnag.const *)

(* begin module interact.const *)
      maxstring = 150; (* the maximum string *) 
(* end module interact.const version = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module pic.const *)
   pi = 3.14159265354; (* circumference divided by diameter of circle *)
   picfield = 12; (* width of numbers printed to the file *)
   picwidth = 9; (* number of decimal places for numbers *)
   charwidth = 0.05; (* the width of characters in the graphic space.
                       this allows centering of strings. *)
   defscale = 81; (* default scale factor.  coordinate units per inch *)

(*   scale = 1.252;*) (* scale factor.  81 pixles per inch.  the routines
                     use inches and this factor converts
                     for the sun pixle size (i think) *)

(* end module pic.const version = 1.28; (@ of domod, 1988 mar 2 *)

type
(* begin module interact.type *)
      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 *) 
      end;  
(* end module interact.type version = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module dnag.atom.type *)
     atom = record (* define an atom in space *)
        group: char; (* the chemical group the atom is a part of *)
        id: array[1..2] of char; (* identification string for the atom.
            first letter is atom type.  second is atom number.
            Me means methyl *)
        radius: real; (* radial position, in angstrom *)
        angle: real; (* radial position, in degrees *)
        z: real; (* z axis coordinate, in angstrom *)
     end;
(* end module dnag.atom.type *)

(* begin module dnag.group.type *)
     group = ^groupatom;
     groupatom = record (* a group of atoms in space *)
        element: atom; (* one of the atoms *)
        next: group
     end;
(* end module dnag.group.type *)

(* begin module dnag.na.type *)
   (* definition of the structure of components of a nucleic acid *)
   nastructure = record
      p: group; (* phosphate *)
      d: group; (* deoxyribose *)
      a: group; (* adenine *)
      c: group; (* cytosine *)
      g: group; (* guanine *)
      t: group; (* thymine *)
   end;
(* end module dnag.na.type *)

var 
(* begin module dnag.var *)
   toradians: real; (* conversion factor to get to radians from degrees *)

   bdna, (* b form dna coordinates *)
   dooin: (* output of this program, input to pic program *)
   (* dnap: (@ parameters to control the program *)
      text;
(* end module dnag.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 *)
(* end module pic.var version = 1.28; (@ of domod, 1988 mar 2 *)

(* 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 = 1.28; (@ of domod, 1988 mar 2 *)
  
(* 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 = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module interact.clearstring *)
procedure clearstring(var ribbon: string);  
(* empty the string *)
var   index: integer; (* to the ribbon *) 
begin (* clearstring *) 
      with ribbon do begin  
         for index := 1 to maxstring do letters[index] := ' ';
         length := 0; 
         current := 0;
      end 
end; (* clearstring *)
(* end module interact.clearstring version = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module interact.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 interact.writestring version = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module pic.functions *)
(* ********************************************************************** *)
(* begin module pic.await *)
procedure await;
(* gutted procedure for now *)
(* Wait for user to type a carriage return.  the routine assumes that
there is a global file called input. *)
begin
(* the old way:
    writeln(output,'awaiting for a Return to continue');
    while not eoln(input) do begin get(input) end;
*)
    (* read past the input *)
(*
    readln(input)
*)
end;
(* end module pic.await version = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module pic.startpic *)
procedure startpic(var afile:text; setscale,x,y: real);
(* open the graphics field, with the given scale, and at (x,y)
in that scale.  scale is in device coordinates per inch. *)
(* open the graphics field *)
(* start pic output to file afile, set the globals *)
begin
   writeln(afile,'.PS',
                 ' ',setscale:picfield:picwidth,
                 ' ',       x:picfield:picwidth,
                 ' ',       y:picfield:picwidth);
   scale := setscale; (* set the global scale *)

   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 = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module pic.stoppic *)
procedure stoppic(var afile:text);
(* stop pic output to file afile *)
begin
   writeln(afile,'.PE');
   inpicture := false;
end;
(* end module pic.stoppic version = 1.28; (@ of domod, 1988 mar 2 *)

(* 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). *)
begin (* drawr *)
   writeln(afile,'drawr',
                 ' ',dx:picfield:picwidth,
                 ' ',dy:picfield:picwidth,
                 ' ',visibility,
                 ' ',spacing:picfield:picwidth);
   picxglobal := picxglobal + dx;
   picyglobal := picyglobal + dy;
end;
(* end module pic.drawr version = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module pic.mover *)
procedure mover(var afile: text; dx,dy: real);
(* move relative the amount (dx, dy). *)
begin
   writeln(afile,'mover',
                 ' ',dx:picfield:picwidth,
                 ' ',dy:picfield:picwidth);
   picxglobal := picxglobal + dx;
   picyglobal := picyglobal + dy;
end;
(* end module pic.mover version = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module pic.liner *)
procedure liner(var afile: text; dx,dy: real);
(* draw a line the relative amount (dx, dy). *)
begin
   writeln(afile,'liner',
                 ' ',dx:picfield:picwidth,
                 ' ',dy:picfield:picwidth);
   picxglobal := picxglobal + dx;
   picyglobal := picyglobal + dy;
end;
(* end module pic.liner version = 1.28; (@ of domod, 1988 mar 2 *)


(* 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). *)
begin
   writeln(afile,'drawa',
                 ' ',x:picfield:picwidth,
                 ' ',y:picfield:picwidth,
                 ' ',visibility,
                 ' ',spacing:picfield:picwidth);
   picxglobal := x;
   picyglobal := y;
end;
(* end module pic.drawa version = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module pic.movea *)
procedure movea(var afile: text; x,y: real);
(* move to absolute x and y *)
begin
   writeln(afile,'movea',
                 ' ',x:picfield:picwidth,
                 ' ',y:picfield:picwidth);
   picxglobal := x;
   picyglobal := y;
end;
(* end module pic.movea version = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module pic.linea *)
procedure linea(var afile: text; x,y: real);
(* draw a line from current position to absolute x and y *)
begin
   writeln(afile,'linea',
                 ' ',x:picfield:picwidth,
                 ' ',y:picfield:picwidth);
   picxglobal := x;
   picyglobal := y;
end;
(* end module pic.linea version = 1.28; (@ of domod, 1988 mar 2 *)

(* 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 *)
var i: integer; (* index to s *)
    quoted: boolean; (* true if the string is quoted *)
    sq: string; (* s, with quotes around it to indicate centering *)
begin
   with s do begin
      if length > 2 
      then if (letters[1]='"') and (letters[length]='"')
           then quoted := true
           else quoted := false
      else quoted := false;
 
      if (not quoted) and centered
      then begin
         clearstring(sq);
         sq.length := s.length + 2;
         sq.letters[1] := '"';
         for i := 1 to s.length do sq.letters[i+1] := s.letters[i];
         sq.letters[sq.length] := '"';
         writestring(tofile,sq)
      end
      else writestring(tofile,s); (* just echo it *)
      writeln(tofile) (* complete the line *)
   end
end;
(* end module pic.graphstring version = 1.28; (@ of domod, 1988 mar 2 *)

(* 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 = 1.28; (@ of domod, 1988 mar 2 *)

(* 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 = 1.28; (@ of domod, 1988 mar 2 *)

(* 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 *)
begin
   write(afile,'picnumber',
               ' ',dx:picfield:picwidth,
               ' ',dy:picfield:picwidth,
               ' ',number:picfield:picwidth,
               ' ',width:2,
               ' ',decimal:2);
   if centered then write(afile,' true')
               else write(afile,' false');
   writeln(afile);
end;
(* end module pic.picnumber version = 1.28; (@ of domod, 1988 mar 2 *)

(* 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
   writeln(afile,'xtic',
                 ' ',length:picfield:picwidth,
                 ' ',dx:picfield:picwidth,
                 ' ',dy:picfield:picwidth,
                 ' ',number:picfield:picwidth,
                 ' ',width:picfield,
                 ' ',decimal:picfield);
end;
(* end module pic.xtic version = 1.28; (@ of domod, 1988 mar 2 *)

(* 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
   writeln(afile,'ytic',
               ' ',length:picfield:picwidth,
               ' ',dx:picfield:picwidth,
               ' ',dy:picfield:picwidth,
               ' ',number:picfield:picwidth,
               ' ',width:picfield,
               ' ',decimal:picfield);
end;
(* end module pic.ytic version = 1.28; (@ of domod, 1988 mar 2 *)

(* 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. *)
begin
   writeln(afile,'xaxis',
                 ' ',axlength:picfield:picwidth,
                 ' ',fromtic:picfield:picwidth,
                 ' ',interval:picfield:picwidth,
                 ' ',totic:picfield:picwidth,
                 ' ',length:picfield:picwidth,
                 ' ',dx:picfield:picwidth,
                 ' ',dy:picfield:picwidth,
                 ' ',width:picfield,
                 ' ',decimal:picfield);
end;
(* end module pic.xaxis version = 1.28; (@ of domod, 1988 mar 2 *)

(* 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. *)
begin
   writeln(afile,'yaxis',
                 ' ',aylength:picfield:picwidth,
                 ' ',fromtic:picfield:picwidth,
                 ' ',interval:picfield:picwidth,
                 ' ',totic:picfield:picwidth,
                 ' ',length:picfield:picwidth,
                 ' ',dx:picfield:picwidth,
                 ' ',dy:picfield:picwidth,
                 ' ',width:picfield,
                 ' ',decimal:picfield);
end;
(* end module pic.yaxis version = 1.28; (@ of domod, 1988 mar 2 *)
(* ********************************************************************** *)
(* end module pic.functions version = 1.28; (@ of domod, 1988 mar 2 *)

(* begin module pic.circler *)
procedure circler(var afile: text; radius: real);
(* make a circle at the current position of some radius. *)
begin
   writeln(afile,'circler',
                 ' ',radius:picfield:picwidth);
end;
(* end module pic.circler version = 1.28; (@ of domod, 1988 mar 2 *)

(* 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 = 1.28; (@ of domod, 1988 mar 2 *)

(* ********************************************************************** *)
  
(* begin module skipblanks *) 
procedure skipblanks(var thefile: text);  
(* skip over blanks until a non-blank, or end of line, is found *)
begin 
      while (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 (thefile^ <> ' ') and not eoln(thefile) do get(thefile);  
end;  
(* end module skipblanks version = 1.28; (@ of domod, 1988 mar 2 *)
(* ********************************************************************** *)

(* begin module dnag.readparameters *)
(*procedure readparameters(var afile: text; var p: param);*)
(* read from afile the parameters p, as defined by type param *)
(*begin
      reset(afile);
      with p do begin
         if eof(afile)
         then wave := false
         else begin
            wave := true;
            readln(afile,extreme);
            readln(afile,wavelocation);
            readln(afile,wavebit);
            readln(afile,waveamplitude);
            readln(afile,wavelength);
         end
      end
end;*)
(* end module dnag.readparameters *)

(* begin module dnag.header *)
procedure header(var outfile: text);
(* do the header of the plot *)
begin
      rewrite(outfile);

      writeln(outfile,'.nf');
      writeln(outfile,'dnag ',version:4:2);
end;
(* end module dnag.header *)

(* begin module dnag.readatom *)
procedure readatom(var afile: text; var a: atom; var found: boolean);
(* read an atom's coordinates from afile.
if it was read, then found is true.  eof means found is false.
data is assumed to be in the form of 4 fixed characters left justified,
followed by the other atom data:
p o1 8.75 97.4 3.63
the data represent:
group - atom - radius (angstrom) - angle (degree) - z axis (angstrom)
lines begining with '*' are skipped.
*)
var
     done: boolean; (* done searching for a data line (or eof) *)
begin
     with a do begin
        done := false;
        repeat
           if eof(afile)
           then begin
              found := false;
              done := true
           end
           else if afile^='*'
                then readln(afile)
                else begin
                   found := true;
                   done := true
                end
        until done;

        if found then begin
           read(afile,group);
           get(afile); (* skip the space *)
           read(afile,id[1],id[2]);
           readln(afile,radius,angle,z);
        end
     end
(*
     if found then with a do begin
        write(output,group,'-',id[1],'-',id[2]);
        write(output,radius:6:2, angle:6:2,z:6:2);
     end
*)
end;
(* end module dnag.readatom *)

(* begin module dnag.readgroup *)
procedure readgroup(var infile: text; c: char; var g: group);
(* obtain the coordinates and data for all atoms of group c into structure g *)
var
     current: group; (* the current atom being added to the group *)
     found: boolean; (* true if an atom was found in infile *)
     previous: group; (* the previous atom just added to the group,
        kept track of to put nil when no more atoms exist *)
begin
      reset(infile);
      new(g);
      current := g;
      previous := current;

      repeat
         readatom(infile,current^.element,found);
         if found
         then begin
            if current^.element.group = c then begin
               new(current^.next);
               previous := current;
               current := current^.next
            end
         end
         else begin
            previous^.next := nil;
            dispose(current)
         end
      until not found
end;
(* end module dnag.readgroup *)

(* begin module dnag.readnastructure *)
procedure readnastructure(var infile: text; var b: nastructure);
(* read the atomic groups for b form dna from infile *)
begin
   with b do begin
      readgroup(infile,'P',p); (* phosphate *)
      readgroup(infile,'D',d); (* deoxyribose *)
      readgroup(infile,'A',a); (* adenine *)
      readgroup(infile,'C',c); (* cytosine *)
      readgroup(infile,'G',g); (* guanine *)
      readgroup(infile,'T',t); (* thymine *)
   end
end;
(* end module dnag.readnastructure *)

(* begin module dnag.atomicradius *)
function atomicradius(e1,e2: char): real;
(* give the radius corresponding to the element, given as two characters,
e1 and e2.  for example, helium is e1='H', e2='e'
source:
      M. Karplus and R. N. Porter
      Atoms & Molecules
      Benjamin/Cummings Publishing Co., Menlo Park, Ca, 1970
      p. 204-7, crystal radii
*) 
procedure unknown;
(* the particular element is not known to this procedure *)
begin
      writeln(output,'radius of element ',e1,e2,' not known');
      halt
end;
begin
      if e2=' '
      then    if e1 = 'H' then atomicradius := 0.25
         else if e1 = 'C' then atomicradius := 0.70
         else if e1 = 'N' then atomicradius := 0.65
         else if e1 = 'O' then atomicradius := 0.50
         else if e1 = 'P' then atomicradius := 1.00

          (* 'M ' means methyl group, ie, carbon *)
         else if e1 = 'M' then atomicradius := atomicradius('C',' ')
         else unknown
       else unknown
end;
(* end module dnag.atomicradius *)

(* begin module dnag.drawatom *)
procedure drawatom(var outfile: text; a: atom; zerox, zeroy, theta, scale: real;
                      dyadic: boolean);
(* draw the atom to file outfile, rotated by angle theta (degrees) given,
with the center of the polar coordinates at (zerox, zeroy).
if dyadic is true, draw at the dyadic position (negative angle) *)
var
      x,y: real; (* coordinate in the plot *)
begin
      with a do begin
         if dyadic
         then polrec(radius*scale,(theta-angle)*toradians,x,y)
         else polrec(radius*scale,(theta+angle)*toradians,x,y);

         movea(outfile,x+zerox,y+zeroy);
         circler(outfile,scale*atomicradius(id[1],' '));

         movea(outfile,x+zerox+fudgex,y+zeroy+fudgey);
         writeln(outfile,'"',id[1],id[2],'"');
      end
end;
(* end module dnag.drawatom *)

(* begin module dnag.drawgroup *)
procedure drawgroup(var outfile: text; g: group;
                    zerox, zeroy, theta: real;
                    scale: real; dyadic: boolean);
(* draw to outfile the atomic group, at some angle theta (degrees) and
on the given scale (inches per angstrom).  the center of the polar
coordinates is (zerox, zeroy).  if dyadic is true, draw
the group in the complementary position (negative angle).
   if the group is a base (a,c,g,t) then label it with its letter
at the average position of carbons C5 and C2, which is
the same on all bases. *)
var
   current: group; (* the current atom of the group *)
   labeling: boolean; (* true means attach a label to the group *)
   x2,y2: real; (* coordinates of the 2 carbon *)
   x5,y5: real; (* coordinates of the 5 carbon *)
begin
   current := g;
   (* label the 4 bases *)
   labeling := g^.element.group in ['A','C','G','T'];

   while current <> nil do begin
      drawatom(outfile, current^.element, zerox, zeroy, theta, scale, dyadic);
      if labeling
      then with current^.element do begin
         if id[1]='C' then begin
            if id[2]='2'
            then if dyadic
                 then polrec(scale*radius,(theta-angle)*toradians,x2,y2)
                 else polrec(scale*radius,(theta+angle)*toradians,x2,y2);
            if id[2]='5'
            then if dyadic
                 then polrec(scale*radius,(theta-angle)*toradians,x5,y5)
                 else polrec(scale*radius,(theta+angle)*toradians,x5,y5);
         end
      end;
      current := current^.next
   end;

   if labeling then begin
      movea(outfile,(x2+x5)/2+zerox+fudgex,(y2+y5)/2+zeroy+fudgey);
      (* put the letter in: *)
      writeln(outfile,'"',g^.element.group,'"');
   end;

end;
(* end module dnag.drawgroup *)

(* begin module dnag.drawbase *)
procedure drawbase(var outfile: text; b: char; dna: nastructure;
                   zerox, zeroy, angle: real; scale: real; dyadic: boolean);
(* draw to outfile the base defined by b, at some angle (degrees) and
on the given scale (inches per angstrom).  the center of the polar
coordinates is at (zerox, zeroy).  if dyadic is true, draw
the group in the complementary position (negative angle). *)
begin
   with dna do begin
      drawgroup(outfile,p,zerox,zeroy,angle,scale,dyadic);
      drawgroup(outfile,d,zerox,zeroy,angle,scale,dyadic);
      case b of
         'a': drawgroup(outfile,a,zerox,zeroy,angle,scale,dyadic);
         'c': drawgroup(outfile,c,zerox,zeroy,angle,scale,dyadic);
         'g': drawgroup(outfile,g,zerox,zeroy,angle,scale,dyadic);
         't': drawgroup(outfile,t,zerox,zeroy,angle,scale,dyadic);
      end
   end
end;
(* end module dnag.drawbase *)

(* begin module dnag.drawbp *)
procedure drawbp(var outfile: text; b: char; dna: nastructure;
                 zerox, zeroy, angle: real; scale: real);
(* draw a basepair to outfile.
place the group corresponding to b on one side of the helix,
and the other on the other side, in the dyad position.
use the coordinates in dna.  draw at some angle (degrees) and
on the given scale (inches per angstrom).  the center of the polar
polar coordinates is at (zerox, zeroy).  if dyadic is true, draw
the group in the complementary position (negative angle). *)
begin
   drawbase(outfile,b,dna,zerox,zeroy,angle,scale,false);
   case b of
      'a': drawbase(outfile,'t',dna,zerox,zeroy,angle,scale,true);
      'c': drawbase(outfile,'g',dna,zerox,zeroy,angle,scale,true);
      'g': drawbase(outfile,'c',dna,zerox,zeroy,angle,scale,true);
      't': drawbase(outfile,'a',dna,zerox,zeroy,angle,scale,true);
   end
end;
(* end module dnag.drawbp *)

(* begin module dnag.dnacircle *)
procedure dnacircle(var outfile: text; dnaradius: real;
                 zerox, zeroy, angle: real; scale: real);
(* draw a circle at the diameter of dna, with a line at angle (degrees)
to represent the zero angle (on the diagram produced).
the center of the polar polar coordinates is at (zerox, zeroy).  
put a circle of 1 angstrom diameter in the center.
use the given scale (inches per angstrom).
dnaradius is the radius of dna in angstroms (10.4 hits the o2 of the
phosphate group) *)
var
      x,y: real; (* coordinates for the end of the marker line *)
begin
      movea(outfile,zerox,zeroy);
      circler(outfile,scale/2); (* center of the axis, diameter 1 angstrom *)
      circler(outfile,dnaradius*scale); (* outer diameter *)
      polrec(dnaradius*scale,angle*toradians,x,y); (* zero axis *)
      linea(outfile,x+zerox,y+zeroy);
end;
(* end module dnag.dnacircle *)

(* begin module dnag.themain *)
procedure themain(var bdna, outfile: text); 
(* the main procedure of the program *) 
var
   dna: nastructure; (* the structure of dna *)
begin
      writeln(output,'dnag ',version:4:2);

(*      readparameters(dnag,params);*)

      header(outfile);
      readnastructure(bdna,dna);
      toradians := 2*pi/360;

      startpic(outfile,81.0,4.0,7.0);
      dnacircle(outfile,10.0,0.0,0.0,-90.0,iascale);
      drawbp(outfile,'a',dna,0.0,0.0,-90.0,iascale);
(*
      stoppic(outfile);

      writeln(outfile,'.bp');

as above      startpic(outfile,81.0,3.0,9.0);
*)
      dnacircle(outfile,10.0,0.0,-4.3,-90.0,iascale);
      drawbp(outfile,'g',dna,0.0,-4.3,-90.0,iascale);
      stoppic(outfile);
end;  
(* end module dnag.themain *)
 
begin
      themain(bdna, dooin);
1: end.