/*
 * pdb.c
 *
 * Summary:  Read trajectories from PDB file.
 *
 * Author:  David Hardy
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mdtypes.h"
#include "pdb.h"
#include "pdbrec.h"
#include "error.h"


int pdb_startup(void)
{
  return pdbrec_startup();
}

void pdb_shutdown(void)
{
  pdbrec_shutdown();
}

int pdb_init(Pdb *p)
{
  memset(p, 0, sizeof(Pdb));
  return 0;
}

void pdb_destroy(Pdb *p)
{
  /* nothing to do! */
}

int pdb_open(Pdb *p, const char *filename, int filetype)
{
  p->f = fopen(filename, "r");
  MDIO_ASSERT(p->f != NULL);
  return pdbrec_init(&(p->pdbrec), filetype);
}

int pdb_close(Pdb *p)
{
  int ret = fclose(p->f);
  MDIO_ASSERT(ret == 0);
  memset(p, 0, sizeof(Pdb));
  return ret;
}

int pdb_read(Pdb *p, void *vec, void *a, void *b, int elemtype, int buflen,
             int *numread)
{
  int mask, i, k, ret;
  FILE *f = p->f;
  Pdbrec *pdbrec = &(p->pdbrec);
  Pdbrec_Data *data = &(p->data);
  MD_Fvec *fvec = NULL;
  float *fa = NULL, *fb = NULL;
  MD_Dvec *dvec = NULL;
  double *da = NULL, *db = NULL;
  char buf[128];

  mask = -1;           /* block everything */
  if (vec != NULL) {
    UNMASK(mask,  8);
    UNMASK(mask,  9);
    UNMASK(mask, 10);
  }
  if (a != NULL) {
    UNMASK(mask, 11);
  }
  if (b != NULL) {
    UNMASK(mask, 12);
  }
  p->data.d.atom.mask = mask;

  MDIO_ASSERT(elemtype == FLOATTYPE || elemtype == DOUBLETYPE);
  if (elemtype == FLOATTYPE) {
    fvec = (MD_Fvec *) vec;
    fa = (float *) a;
    fb = (float *) b;
  }
  else {
    dvec = (MD_Dvec *) vec;
    da = (double *) a;
    db = (double *) b;
  }

  k = 0;
  while (k < buflen && fgets(buf, sizeof(buf), f) != NULL) {
    /* pad out string with blanks */
    i = strlen(buf);
    if (buf[i-1] == '\n') i--;
    while (i < 80) buf[i++] = ' ';
    buf[80] = '\n';
    buf[81] = '\0';
    MDIO_ASSERT(strlen(buf) == 81);  /* 80 chars + newline */
    ret = pdbrec_extract(pdbrec, buf);
    MDIO_ASSERT(ret >= 0 && ret <= PDBREC_END);
    if (ret == PDBREC_ATOM || ret == PDBREC_HETATM) {
      ret = pdbrec_getdata(pdbrec, data);
      MDIO_ASSERT(ret == 0);

      if (dvec != NULL) {
        dvec[k].x = (double) data->d.atom.x;
        dvec[k].y = (double) data->d.atom.y;
        dvec[k].z = (double) data->d.atom.z;
      }
      else if (fvec != NULL) {
        fvec[k].x = data->d.atom.x;
        fvec[k].y = data->d.atom.y;
        fvec[k].z = data->d.atom.z;
      }

      if (da != NULL) {
        da[k] = (double) data->d.atom.occupancy;
      }
      else if (fa != NULL) {
        fa[k] = data->d.atom.occupancy;
      }

      if (db != NULL) {
        db[k] = (double) data->d.atom.tempFactor;
      }
      else if (fb != NULL) {
        fb[k] = data->d.atom.tempFactor;
      }

      k++;
    }
  }

  *numread = k;
  return 0;
}