/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */
/* From input file "dnag.p" */


#include <p2c/p2c.h>


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

/* end of program */

/* begin module version */

#define 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 */
#define iascale         0.20   /* inches per angstrom across the plot */
#define fudgex          (-0.02)   /* fudge factor for x, in angstrom */
#define fudgey          (-0.04)   /* fudge factor for y, in angstrom */
/* end module dnag.const */

/* begin module interact.const */

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

/* begin module pic.const */

#define pi              3.14159265354
    /* circumference divided by diameter of circle */

#define picfield        12   /* width of numbers printed to the file */
#define picwidth        9   /* number of decimal places for numbers */

#define charwidth       0.05
/* the width of characters in the graphic space.
                    this allows centering of strings. */

#define 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 */

/* begin module interact.type */

typedef struct string {
  /* a string of characters */
  Char letters[maxstring];   /* the letters in the string */
  long length;   /* the number of characters in the string */
  long current;   /* the letter we are working on */
} string;

/* end module interact.type version = 1.28; (@ of domod, 1988 mar 2 */

/* begin module dnag.atom.type */

typedef struct atom {
  /* define an atom in space */
  Char group;   /* the chemical group the atom is a part of */
  Char id[2];
  /* identification string for the atom.
      first letter is atom type.  second is atom number.
      Me means methyl */
  double radius;   /* radial position, in angstrom */
  double angle;   /* radial position, in degrees */
  double z;   /* z axis coordinate, in angstrom */
} atom;

/* end module dnag.atom.type */

/* begin module dnag.group.type */

typedef struct groupatom {
  /* a group of atoms in space */
  atom element;   /* one of the atoms */
  struct groupatom *next;
} groupatom;

/* end module dnag.group.type */

/* begin module dnag.na.type */
/* definition of the structure of components of a nucleic acid */

typedef struct nastructure {
  groupatom *p;   /* phosphate */
  groupatom *d;   /* deoxyribose */
  groupatom *a;   /* adenine */
  groupatom *c;   /* cytosine */
  groupatom *g;   /* guanine */
  groupatom *t;   /* thymine */
} nastructure;

/* end module dnag.na.type */


/* begin module dnag.var */
Static double toradians;
    /* conversion factor to get to radians from degrees */

Static _TEXT bdna;   /* b form dna coordinates */
/* dnap: (@ parameters to control the program */
Static _TEXT dooin;   /* output of this program, input to pic program */
/* end module dnag.var */

/* begin module pic.var */
Static boolean inpicture;
/* true if we are drawing the picture,
                       ie, startpic has been called */
Static double picxglobal, picyglobal;   /* absolute location in the graph */
Static double pictolerance;
/* 10 raised to the picwidth,
   to detect values close to zero */
Static double scale;   /* scale factor.  graphic coordinate units per inch */


Static jmp_buf _JL1;


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

/* begin module halt */
Static Void 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. */
  printf(" program halt.\n");
  longjmp(_JL1, 1);
}


/* end module halt version = 1.28; (@ of domod, 1988 mar 2 */

/* begin module copyaline */
Static Void copyaline(fin, fout)
_TEXT *fin, *fout;
{
  /* copy a line from file fin to file fout */
  while (!P_eoln(fin->f)) {
    putc(P_peek(fin->f), fout->f);
    getc(fin->f);
  }
  fscanf(fin->f, "%*[^\n]");
  getc(fin->f);
  putc('\n', fout->f);
}  /* copyaline */


/* end module copyaline version = 1.28; (@ of domod, 1988 mar 2 */

/* begin module interact.clearstring */
Static Void clearstring(ribbon)
string *ribbon;
{
  /* empty the string */
  long index;   /* to the ribbon */

  for (index = 0; index < maxstring; index++)
    ribbon->letters[index] = ' ';
  ribbon->length = 0;
  ribbon->current = 0;
}  /* clearstring */


/* end module interact.clearstring version = 1.28; (@ of domod, 1988 mar 2 */

/* begin module interact.writestring */
Static Void writestring(tofile, s)
_TEXT *tofile;
string *s;
{
  /* write the string s to file tofile, no writeln */
  long i;   /* index to s */
  long FORLIM;

  FORLIM = s->length;
  for (i = 0; i < FORLIM; i++)
    putc(s->letters[i], tofile->f);
}  /* writestring */


/* end module interact.writestring version = 1.28; (@ of domod, 1988 mar 2 */

/* begin module pic.functions */
/* ********************************************************************** */
/* begin module pic.await */
Static Void await()
{
  /* gutted procedure for now */
  /* Wait for user to type a carriage return.  the routine assumes that
there is a global file called input. */
  /* 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 module pic.await version = 1.28; (@ of domod, 1988 mar 2 */

/* begin module pic.startpic */
Static Void startpic(afile, setscale, x, y)
_TEXT *afile;
double setscale, x, y;
{
  /* 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 */
  fprintf(afile->f, ".PS %*.*f %*.*f %*.*f\n",
	  picfield, picwidth, setscale, picfield, picwidth, x, picfield,
	  picwidth, y);
  scale = setscale;   /* set the global scale */

  inpicture = true;
  picxglobal = 0.0;
  picyglobal = 0.0;
  pictolerance = (long)(exp(picwidth * log(10.0)) + 0.5);
  /*;writeln(output,'pictolerance = ',pictolerance:picfield:picwidth);*/
}


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

/* begin module pic.stoppic */
Static Void stoppic(afile)
_TEXT *afile;
{
  /* stop pic output to file afile */
  fprintf(afile->f, ".PE\n");
  inpicture = false;
}


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

/* begin module pic.drawr */
Static Void drawr(afile, dx, dy, visibility, spacing)
_TEXT *afile;
double dx, dy;
Char visibility;
double spacing;
{
  /* 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). */
  fprintf(afile->f, "drawr %*.*f %*.*f %c %*.*f\n",
	  picfield, picwidth, dx, picfield, picwidth, dy, visibility,
	  picfield, picwidth, spacing);
  picxglobal += dx;
  picyglobal += dy;
}


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

/* begin module pic.mover */
Static Void mover(afile, dx, dy)
_TEXT *afile;
double dx, dy;
{
  /* move relative the amount (dx, dy). */
  fprintf(afile->f, "mover %*.*f %*.*f\n",
	  picfield, picwidth, dx, picfield, picwidth, dy);
  picxglobal += dx;
  picyglobal += dy;
}


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

/* begin module pic.liner */
Static Void liner(afile, dx, dy)
_TEXT *afile;
double dx, dy;
{
  /* draw a line the relative amount (dx, dy). */
  fprintf(afile->f, "liner %*.*f %*.*f\n",
	  picfield, picwidth, dx, picfield, picwidth, dy);
  picxglobal += dx;
  picyglobal += dy;
}


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


/* begin module pic.drawa */
Static Void drawa(afile, x, y, visibility, spacing)
_TEXT *afile;
double x, y;
Char visibility;
double spacing;
{
  /* 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). */
  fprintf(afile->f, "drawa %*.*f %*.*f %c %*.*f\n",
	  picfield, picwidth, x, picfield, picwidth, y, visibility, picfield,
	  picwidth, spacing);
  picxglobal = x;
  picyglobal = y;
}


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

/* begin module pic.movea */
Static Void movea(afile, x, y)
_TEXT *afile;
double x, y;
{
  /* move to absolute x and y */
  fprintf(afile->f, "movea %*.*f %*.*f\n",
	  picfield, picwidth, x, picfield, picwidth, y);
  picxglobal = x;
  picyglobal = y;
}


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

/* begin module pic.linea */
Static Void linea(afile, x, y)
_TEXT *afile;
double x, y;
{
  /* draw a line from current position to absolute x and y */
  fprintf(afile->f, "linea %*.*f %*.*f\n",
	  picfield, picwidth, x, picfield, picwidth, y);
  picxglobal = x;
  picyglobal = y;
}


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

/* begin module pic.graphstring */
Static Void graphstring(tofile, s, centered)
_TEXT *tofile;
string *s;
boolean centered;
{
  /* 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 */
  long i;   /* index to s */
  boolean quoted;   /* true if the string is quoted */
  string sq;   /* s, with quotes around it to indicate centering */
  long FORLIM;

  if (s->length > 2) {
    if (s->letters[0] == '"' && s->letters[s->length - 1] == '"')
      quoted = true;
    else
      quoted = false;
  } else
    quoted = false;

  if (!quoted && centered) {
    clearstring(&sq);
    sq.length = s->length + 2;
    sq.letters[0] = '"';
    FORLIM = s->length;
    for (i = 1; i <= FORLIM; i++)
      sq.letters[i] = s->letters[i-1];
    sq.letters[sq.length - 1] = '"';
    writestring(tofile, &sq);
  } else
    writestring(tofile, s);
  /* just echo it */
  putc('\n', tofile->f);   /* complete the line */
}


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

/* begin module pic.stringinteger */
Static Void stringinteger(number, name, width, leadingzeros)
long number;
string *name;
long width;
boolean leadingzeros;
{
  /* 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. */
  long bigdigit;   /* the location of the biggest digit */
  long dig;   /* number of digits in the number */
  long place;   /* place to write the next digit of the number */
  long sign;   /* the sign of the number */

  if (number < 0) {
    sign = -1;
    name->length++;   /* provide room for the sign!! */
    number = -number;
    if (leadingzeros)
      printf(
	"WARNING: stringinteger: the sign of a negative number with leading zeros is lost\n");
  } 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)
    dig = (long)(log(number + 0.1) / log(10.0)) + 1;
  else
    dig = 1;

  if (dig > width) {
    printf("stringinteger: number width too small\n");
    printf("%ld digit number (%ld)\n", dig, number);
    printf("does not fit in %ld characters\n", width);
    halt();
  }
  if (leadingzeros)
    bigdigit = name->length + 1;   /* no sign if leading zeros */
  else {
    bigdigit = name->length + width - dig + 1;

    if (bigdigit <= name->length && sign < 0) {
      printf("stringinteger: no room for sign\n");
      halt();
    }
  }
  if (sign < 0)
    name->letters[bigdigit-2] = '-';

  for (place = name->length + width - 1; place >= bigdigit - 1; place--) {
/* p2c: dnag.p, line 436:
 * Note: Using % for possibly-negative arguments [317] */
    switch (number % 10) {

    case 0:
      name->letters[place] = '0';
      break;

    case 1:
      name->letters[place] = '1';
      break;

    case 2:
      name->letters[place] = '2';
      break;

    case 3:
      name->letters[place] = '3';
      break;

    case 4:
      name->letters[place] = '4';
      break;

    case 5:
      name->letters[place] = '5';
      break;

    case 6:
      name->letters[place] = '6';
      break;

    case 7:
      name->letters[place] = '7';
      break;

    case 8:
      name->letters[place] = '8';
      break;

    case 9:
      name->letters[place] = '9';
      break;
    }
    number /= 10;
  }
  name->length += width;
}


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

/* begin module pic.stringreal */
Static Void stringreal(number, name, width, decimal)
double number;
string *name;
long width, decimal;
{
  /* 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. */
  long abovezero;
  /* the number shifted above the decimal place, to
     'decimal' positions (and rounded) */
  long shift;
  /* power of ten used to shift a number around
     relative to the decimal point */
  long sign;   /* the sign of the number */
  long thedecimal;   /* integer version of the decimal part of the number */
  long theupper;   /* integer version of the upper part of the number */

  if (number < 0)
    sign = -1;
  else
    sign = 1;

  number = fabs(number);   /* make positive */

  /* the amount to shift the number above zero */
  shift = (long)floor(exp(decimal * log(10.0)) + 0.5);
      /* amount to move above zero */
  abovezero = (long)floor(number * shift + 0.5);
      /* move above zero, round off */
  theupper = (long)((double)abovezero / shift);
  thedecimal = abovezero - shift * theupper;

  /* create the actual real number */
  /* before decimal point */
  stringinteger(sign * theupper, name, width - decimal - 1, false);
  /* put in the decimal point */
  name->length++;
  name->letters[name->length - 1] = '.';
  stringinteger(thedecimal, name, decimal, true);   /* after decimal point */
}


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

/* begin module pic.picnumber */
Static Void picnumber(afile, dx, dy, number, width, decimal, centered)
_TEXT *afile;
double dx, dy, number;
long width, decimal;
boolean centered;
{
  /* 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 */
  fprintf(afile->f, "picnumber %*.*f %*.*f %*.*f %2ld %2ld",
	  picfield, picwidth, dx, picfield, picwidth, dy, picfield, picwidth,
	  number, width, decimal);
  if (centered)
    fprintf(afile->f, " true");
  else
    fprintf(afile->f, " false");
  putc('\n', afile->f);
}


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

/* begin module pic.xtic */
Static Void xtic(afile, length, dx, dy, number, width, decimal)
_TEXT *afile;
double length, dx, dy, number;
long width, decimal;
{
  /* 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. */
  fprintf(afile->f, "xtic %*.*f %*.*f %*.*f %*.*f %*ld %*ld\n",
	  picfield, picwidth, length, picfield, picwidth, dx, picfield,
	  picwidth, dy, picfield, picwidth, number, picfield, width, picfield,
	  decimal);
}


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

/* begin module pic.ytic */
Static Void ytic(afile, length, dx, dy, number, width, decimal)
_TEXT *afile;
double length, dx, dy, number;
long width, decimal;
{
  /* 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. */
  fprintf(afile->f, "ytic %*.*f %*.*f %*.*f %*.*f %*ld %*ld\n",
	  picfield, picwidth, length, picfield, picwidth, dx, picfield,
	  picwidth, dy, picfield, picwidth, number, picfield, width, picfield,
	  decimal);
}


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

/* begin module pic.xaxis */
Static Void xaxis(afile, axlength, fromtic, interval, totic, length, dx, dy,
		  width, decimal)
_TEXT *afile;
double axlength, fromtic, interval, totic, length, dx, dy;
long width, decimal;
{
  /* 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. */
  fprintf(afile->f,
	  "xaxis %*.*f %*.*f %*.*f %*.*f %*.*f %*.*f %*.*f %*ld %*ld\n",
	  picfield, picwidth, axlength, picfield, picwidth, fromtic, picfield,
	  picwidth, interval, picfield, picwidth, totic, picfield, picwidth,
	  length, picfield, picwidth, dx, picfield, picwidth, dy, picfield,
	  width, picfield, decimal);
}


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

/* begin module pic.yaxis */
Static Void yaxis(afile, aylength, fromtic, interval, totic, length, dx, dy,
		  width, decimal)
_TEXT *afile;
double aylength, fromtic, interval, totic, length, dx, dy;
long width, decimal;
{
  /* 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. */
  fprintf(afile->f,
	  "yaxis %*.*f %*.*f %*.*f %*.*f %*.*f %*.*f %*.*f %*ld %*ld\n",
	  picfield, picwidth, aylength, picfield, picwidth, fromtic, picfield,
	  picwidth, interval, picfield, picwidth, totic, picfield, picwidth,
	  length, picfield, picwidth, dx, picfield, picwidth, dy, picfield,
	  width, picfield, decimal);
}


/* 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 */
Static Void circler(afile, radius)
_TEXT *afile;
double radius;
{
  /* make a circle at the current position of some radius. */
  fprintf(afile->f, "circler %*.*f\n", picfield, picwidth, radius);
}


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

/* begin module pic.polrec */
Static Void polrec(r, theta, x, y)
double r, theta, *x, *y;
{
  /* convert polar to rectangular coordinates,
theta is in radians */
  *x = r * cos(theta);
  *y = r * sin(theta);
}


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

/* ********************************************************************** */

/* begin module skipblanks */
Static Void skipblanks(thefile)
_TEXT *thefile;
{
  /* skip over blanks until a non-blank, or end of line, is found */
  while ((P_peek(thefile->f) == ' ') & (!P_eoln(thefile->f)))
    getc(thefile->f);
}


Static Void skipnonblanks(thefile)
_TEXT *thefile;
{
  /* skip over nonblanks until a blank, or end of line, is found */
  while ((P_peek(thefile->f) != ' ') & (!P_eoln(thefile->f)))
    getc(thefile->f);
}


/* 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 */
Static Void header(outfile)
_TEXT *outfile;
{
  /* do the header of the plot */
  if (*outfile->name != '\0') {
    if (outfile->f != NULL)
      outfile->f = freopen(outfile->name, "w", outfile->f);
    else
      outfile->f = fopen(outfile->name, "w");
  } else {
    if (outfile->f != NULL)
      rewind(outfile->f);
    else
      outfile->f = tmpfile();
  }
  if (outfile->f == NULL)
    _EscIO2(FileNotFound, outfile->name);
  SETUPBUF(outfile->f, Char);

  fprintf(outfile->f, ".nf\n");
  fprintf(outfile->f, "dnag %4.2f\n", version);
}


/* end module dnag.header */

/* begin module dnag.readatom */
Static Void readatom(afile, a, found)
_TEXT *afile;
atom *a;
boolean *found;
{
  /* 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.
*/
  boolean done = false;   /* done searching for a data line (or eof) */

  /*
       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
  */
  do {
    if (BUFEOF(afile->f)) {
      *found = false;
      done = true;
    } else if (P_peek(afile->f) == '*') {
      fscanf(afile->f, "%*[^\n]");
      getc(afile->f);
    } else {
      *found = true;
      done = true;
    }
  } while (!done);

  if (!*found)
    return;
  a->group = getc(afile->f);
  if (a->group == '\n')
    a->group = ' ';
  getc(afile->f);   /* skip the space */
  *a->id = getc(afile->f);
  a->id[1] = getc(afile->f);
  if (a->id[0] == '\n')
    a->id[0] = ' ';
  if (a->id[1] == '\n')
    a->id[1] = ' ';
  fscanf(afile->f, "%lg%lg%lg%*[^\n]", &a->radius, &a->angle, &a->z);
  getc(afile->f);
}


/* end module dnag.readatom */

/* begin module dnag.readgroup */
Static Void readgroup(infile, c, g)
_TEXT *infile;
Char c;
groupatom **g;
{
  /* obtain the coordinates and data for all atoms of group c into structure g */
  groupatom *current;   /* the current atom being added to the group */
  boolean found;   /* true if an atom was found in infile */
  groupatom *previous;

  /* the previous atom just added to the group,
     kept track of to put nil when no more atoms exist */
  if (*infile->name != '\0') {
    if (infile->f != NULL)
      infile->f = freopen(infile->name, "r", infile->f);
    else
      infile->f = fopen(infile->name, "r");
  } else
    rewind(infile->f);
  if (infile->f == NULL)
    _EscIO2(FileNotFound, infile->name);
  RESETBUF(infile->f, Char);
  *g = (groupatom *)Malloc(sizeof(groupatom));
  current = *g;
  previous = current;

  do {
    readatom(infile, &current->element, &found);
    if (found) {
      if (current->element.group == c) {
	current->next = (groupatom *)Malloc(sizeof(groupatom));
	previous = current;
	current = current->next;
      }
    } else {
      previous->next = NULL;
      Free(current);
    }
  } while (found);
}


/* end module dnag.readgroup */

/* begin module dnag.readnastructure */
Static Void readnastructure(infile, b)
_TEXT *infile;
nastructure *b;
{
  /* read the atomic groups for b form dna from infile */
  readgroup(infile, 'P', &b->p);   /* phosphate */
  readgroup(infile, 'D', &b->d);   /* deoxyribose */
  readgroup(infile, 'A', &b->a);   /* adenine */
  readgroup(infile, 'C', &b->c);   /* cytosine */
  readgroup(infile, 'G', &b->g);   /* guanine */
  readgroup(infile, 'T', &b->t);   /* thymine */
}


/* Local variables for atomicradius: */
struct LOC_atomicradius {
  Char e1, e2;
} ;

Local Void unknown(LINK)
struct LOC_atomicradius *LINK;
{
  /* the particular element is not known to this procedure */
  printf("radius of element %c%c not known\n", LINK->e1, LINK->e2);
  halt();
}


/* end module dnag.readnastructure */

/* begin module dnag.atomicradius */
Static double atomicradius(e1_, e2_)
Char e1_, e2_;
{
  /* 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
  */
  struct LOC_atomicradius V;

  V.e1 = e1_;
  V.e2 = e2_;
  if (V.e2 == ' ') {
    if (V.e1 == 'H')
      return 0.25;
    if (V.e1 == 'C')
      return 0.70;
    if (V.e1 == 'N')
      return 0.65;
    if (V.e1 == 'O')
      return 0.50;
    if (V.e1 == 'P') {
      return 1.00;

      /* 'M ' means methyl group, ie, carbon */
    }
    if (V.e1 == 'M')
      return (atomicradius('C', ' '));
    unknown(&V);
  } else
    unknown(&V);
}


/* end module dnag.atomicradius */

/* begin module dnag.drawatom */
Static Void drawatom(outfile, a, zerox, zeroy, theta, scale, dyadic)
_TEXT *outfile;
atom a;
double zerox, zeroy, theta, scale;
boolean dyadic;
{
  /* 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) */
  double x, y;   /* coordinate in the plot */

  if (dyadic)
    polrec(a.radius * scale, (theta - a.angle) * toradians, &x, &y);
  else
    polrec(a.radius * scale, (theta + a.angle) * toradians, &x, &y);

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

  movea(outfile, x + zerox + fudgex, y + zeroy + fudgey);
  fprintf(outfile->f, "\"%c%c\"\n", a.id[0], a.id[1]);
}


/* end module dnag.drawatom */

/* begin module dnag.drawgroup */
Static Void drawgroup(outfile, g, zerox, zeroy, theta, scale, dyadic)
_TEXT *outfile;
groupatom *g;
double zerox, zeroy, theta, scale;
boolean dyadic;
{
  /* 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. */
  groupatom *current = g;   /* the current atom of the group */
  boolean labeling;   /* true means attach a label to the group */
  double x2, y2;   /* coordinates of the 2 carbon */
  double x5, y5;   /* coordinates of the 5 carbon */
  atom *WITH;

  /* label the 4 bases */
  labeling = (g->element.group == 'T' || g->element.group == 'G' ||
	      g->element.group == 'C' || g->element.group == 'A');

  while (current != NULL) {
    drawatom(outfile, current->element, zerox, zeroy, theta, scale, dyadic);
    if (labeling) {
      WITH = &current->element;
      if (WITH->id[0] == 'C') {
	if (WITH->id[1] == '2') {
	  if (dyadic)
	    polrec(scale * WITH->radius, (theta - WITH->angle) * toradians,
		   &x2, &y2);
	  else
	    polrec(scale * WITH->radius, (theta + WITH->angle) * toradians,
		   &x2, &y2);
	}
	if (WITH->id[1] == '5') {
	  if (dyadic)
	    polrec(scale * WITH->radius, (theta - WITH->angle) * toradians,
		   &x5, &y5);
	  else
	    polrec(scale * WITH->radius, (theta + WITH->angle) * toradians,
		   &x5, &y5);
	}
      }
    }
    current = current->next;
  }

  if (labeling) {
    movea(outfile, (x2 + x5) / 2 + zerox + fudgex,
	  (y2 + y5) / 2 + zeroy + fudgey);
    /* put the letter in: */
    fprintf(outfile->f, "\"%c\"\n", g->element.group);
  }

}


/* end module dnag.drawgroup */

/* begin module dnag.drawbase */
Static Void drawbase(outfile, b, dna, zerox, zeroy, angle, scale, dyadic)
_TEXT *outfile;
Char b;
nastructure dna;
double zerox, zeroy, angle, scale;
boolean dyadic;
{
  /* 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). */
  drawgroup(outfile, dna.p, zerox, zeroy, angle, scale, dyadic);
  drawgroup(outfile, dna.d, zerox, zeroy, angle, scale, dyadic);
  switch (b) {

  case 'a':
    drawgroup(outfile, dna.a, zerox, zeroy, angle, scale, dyadic);
    break;

  case 'c':
    drawgroup(outfile, dna.c, zerox, zeroy, angle, scale, dyadic);
    break;

  case 'g':
    drawgroup(outfile, dna.g, zerox, zeroy, angle, scale, dyadic);
    break;

  case 't':
    drawgroup(outfile, dna.t, zerox, zeroy, angle, scale, dyadic);
    break;
  }
}


/* end module dnag.drawbase */

/* begin module dnag.drawbp */
Static Void drawbp(outfile, b, dna, zerox, zeroy, angle, scale)
_TEXT *outfile;
Char b;
nastructure dna;
double zerox, zeroy, angle, scale;
{
  /* 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). */
  drawbase(outfile, b, dna, zerox, zeroy, angle, scale, false);
  switch (b) {

  case 'a':
    drawbase(outfile, 't', dna, zerox, zeroy, angle, scale, true);
    break;

  case 'c':
    drawbase(outfile, 'g', dna, zerox, zeroy, angle, scale, true);
    break;

  case 'g':
    drawbase(outfile, 'c', dna, zerox, zeroy, angle, scale, true);
    break;

  case 't':
    drawbase(outfile, 'a', dna, zerox, zeroy, angle, scale, true);
    break;
  }
}


/* end module dnag.drawbp */

/* begin module dnag.dnacircle */
Static Void dnacircle(outfile, dnaradius, zerox, zeroy, angle, scale)
_TEXT *outfile;
double dnaradius, zerox, zeroy, angle, scale;
{
  /* 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) */
  double x, y;   /* coordinates for the end of the marker line */

  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 module dnag.dnacircle */

/* begin module dnag.themain */
Static Void themain(bdna, outfile)
_TEXT *bdna, *outfile;
{
  /* the main procedure of the program */
  nastructure dna;   /* the structure of dna */

  printf("dnag %4.2f\n", version);

  /*      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 module dnag.themain */

main(argc, argv)
int argc;
Char *argv[];
{
  PASCAL_MAIN(argc, argv);
  if (setjmp(_JL1))
    goto _L1;
  dooin.f = NULL;
  strcpy(dooin.name, "dooin");
  bdna.f = NULL;
  strcpy(bdna.name, "bdna");
  themain(&bdna, &dooin);
_L1:
  if (bdna.f != NULL)
    fclose(bdna.f);
  if (dooin.f != NULL)
    fclose(dooin.f);
  exit(EXIT_SUCCESS);
}



/* End. */