/*
 * pdbio.c
 *
 * read positions from pdb file
 * write them out to another pdb file
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>


/* start with this size array allocated */
#define  INIT_MAX  10

typedef struct Vec_T {
  double x, y, z;
} Vec;

typedef char Line[96];


int main(int argc, char **argv)
{
  FILE *fin, *fout;
  const char *fname_in, *fname_out;
  Vec *p;
  int n = 0, max;
  Line buf;
  Line *line;
  int nlines = 0, maxlines = INIT_MAX;
  int i, j;
  double delta;

  if (argc != 4) {
    fprintf(stderr, "read atom coordinates from a PDB file\n"
        "write them to a new PDB file\n\n");
    fprintf(stderr, "syntax:  %s infile outfile delta\n", argv[0]);
  }
  fname_in = argv[1];
  fname_out = argv[2];
  delta = strtod(argv[3], NULL);
  delta *= 0.5773502692;

  srandom(time(NULL));

  /* init Line array */
  if ((line = (Line *) malloc(maxlines * sizeof(Line))) == NULL) {
    fprintf(stderr, "cannot malloc() memory for Line array\n");
    exit(1);
  }

  /* read in PDB file */
  if ((fin = fopen(fname_in, "r")) == NULL) {
    fprintf(stderr, "cannot open file \"%s\" for reading\n", fname_in);
    exit(1);
  }
  while (fgets(buf, sizeof(buf), fin) != NULL) {
    if (buf[strlen(buf)-1] != '\n') {
      fprintf(stderr, "line length exceeded buffer, "
          "badly formatted PDB file\n");
      exit(1);
    }
    if (nlines == maxlines) {
      void *tmp = realloc(line, 2 * maxlines * sizeof(Line));
      if (tmp == NULL) {
        fprintf(stderr, "cannot realloc() memory for Line array\n");
        exit(1);
      }
      line = (Line *) tmp;
      maxlines *= 2;
    }
    strcpy(line[nlines], buf);
    nlines++;
  }
  if (ferror(fin) || !feof(fin)) {
    fprintf(stderr, "error reading file \"%s\"\n", fname_in);
    exit(1);
  }
  else if (fclose(fin)) {
    fprintf(stderr, "error closing file \"%s\"\n", fname_in);
    exit(1);
  }

  /* parse coordinates from Line array (and blank coord substrings) */
  max = nlines;
  if ((p = (Vec *) calloc(max, sizeof(Vec))) == NULL) {
    fprintf(stderr, "cannot calloc() memory for Vec array\n");
    exit(1);
  }
  for (i = 0;  i < nlines;  i++) {
    if (strncmp(line[i], "ATOM ", 5) == 0) {
      if (sscanf(line[i] + 30, "%8lf%8lf%8lf",
                 &(p[n].x), &(p[n].y), &(p[n].z)) != 3) {
        fprintf(stderr, "incorrect formatting in PDB file\n");
        exit(1);
      }
      /* blank out the coordinate */
      memset(line[i] + 30, ' ', 24);
      n++;
    }
  }

/*
 * do something with coordinate array p[]
 */
  for (i = 0; i < nlines; i++) {
    p[i].x += (random() % 2 ? delta : -delta);
    p[i].y += (random() % 2 ? delta : -delta);
    p[i].z += (random() % 2 ? delta : -delta);
  } 

  /* write out coordinates */
  if ((fout = fopen(fname_out, "w")) == NULL) {
    fprintf(stderr, "cannot open file \"%s\" for writing\n", fname_out);
    exit(1);
  }
  j = 0;
  for (i = 0;  i < nlines;  i++) {
    if (strncmp(line[i], "ATOM ", 5) == 0) {
      sprintf(buf, "%8.3f%8.3f%8.3f", p[j].x, p[j].y, p[j].z);
      strncpy(line[i] + 30, buf, 24);
      j++;
    }
    if (fprintf(fout, "%s", line[i]) != strlen(line[i])) {
      fprintf(stderr, "error writing file \"%s\"\n", fname_out);
      exit(1);
    }
  }
  if (fclose(fout)) {
    fprintf(stderr, "error closing file \"%s\"\n", fname_out);
    exit(1);
  }
  free(p);
  free(line);
  return 0;
}