#include "ReadTraj.h"

/*#####################################################################################*/
/*                          Stuff for get coord mdb                                    */
/*#####################################################################################*/

static int Recsize, Recpos;
static short Ixk[BiG_NuM_MaX];
static float Coords[BiG_NuM_MaX], Sclp, Sclw, Sclp1;

/*#####################################################################################*/



/* Gets the int size of a byte record from record header */
/* returns EOF if at End Of File                         */
static int rget_int(FILE *mdb, int *val)
{
  int i, t;
  char *c;

  c=(char *)val;
  for(i=0;i<INTSIZE;i++)  {
    t=getc(mdb);
    if(t==EOF) {
      return EOF;
    }
    c[i]=t;
  }
  return 0;
}


/* stores the size of the record in Recsize and sets Recpos to zero */
/* returns EOF if at End Of File                                    */
int newrecord(FILE *mdb)
{
  int stat;
  stat=rget_int(mdb, &Recsize);           /* read the record header */
  Recpos=0;
  return stat;
}


/* goes to the end of a record */
void endrecord(FILE *mdb)
{
  int t;
  while (Recpos<Recsize)  {
    t=getc(mdb);
    Recpos++;
  }
  rget_int(mdb, &t);                       /* what is this for??? */
}


/* reads a byte */
int ftget_byte(FILE *mdb, char *val)
{
  int t;
  t=getc(mdb);
  if (t==EOF)
    return EOF;
  Recpos++;
  val[0]=t;
  return 0;
}


/* reads a binary integer */
int ftget_int(FILE *mdb, int *val)
{
  int i, t;
  char *c;
  c=(char *)val;
  for(i=0;i<INTSIZE;i++)  {
    if(Recpos>=Recsize)
      return EOR;
    t=getc(mdb);
    if(t==EOF)
      return EOF;
    c[i]=t;
    Recpos++;
  }
  return 0;
}


/* read in a short integer */
int ftget_short(FILE *mdb, short *val)
{
  int i,t;
  char *c;
  c=(char *)val;
  for (i=0; i<2; i++)  {
    if (Recpos>=Recsize)
      return EOR;
    t=getc(mdb);
    if (t==EOF)
      return EOF;
    c[i]=t;
    Recpos++;
  }
  return 0;
}


/* read in a float */
int ftget_float(FILE *mdb, float *val)
{
  int i,t;
  char *c;
  c=(char *)val;
  for (i=0; i<4; i++)  {
    if (Recpos>=Recsize)
      return EOR;
    t=getc(mdb);
    if (t==EOF)
      return EOF;
    c[i]=t;
    Recpos++;
  }
  return 0;
}


/* do i need this or not? */
/* read a string */
int ftget_str(FILE *mdb, int n, char *p)
{
  int i,t;
  for (i=0; i<n; i++)  {
    if (Recpos>=Recsize)
      return EOR;
    t=getc(mdb);
    if (t==EOF)
      return EOF;
    p[i]=t;
    Recpos++;
  }
  p[i]==0;       /* Cap it off */
  return 0;
}



int get_mdb_coords (file_records *f, int * istep, float * min, float * max)
{
  int icode, nat3p, n, i, j;
  static int after_first_record = 0;
  signed char ixb;
  static FILE * mdb;


  if(after_first_record)  {

    if (newrecord(mdb)==EOF)
      return EOF;
    ftget_int(mdb,&icode);      /* icode specifies stored coords' format */
    endrecord(mdb);
    
    newrecord(mdb);
    ftget_int(mdb,istep);
    ftget_int(mdb,&n);
    if (icode==0) { /* icode = 0 indicates trajectory stored as a byte shift */
      for (i=0; i<n; i++)  {
	ftget_byte(mdb,(char *)&ixb);
	Ixk[i]+=ixb;
      }
    }
    else  {        /* otherwise coordinates written out completely in a short */
      for (i=0; i<n; i++)
	ftget_short(mdb,&(Ixk[i]));
    }
    endrecord(mdb);
    
    for (i=0; i<n; i+=3)  {       /* expands coordinates into floats */
      atom_record * a = f->atoms[i/3];

      a->x = (double)(Ixk[i])*Sclp1*0.01;
      a->y = (double)(Ixk[i+1])*Sclp1*0.01;
      a->z = (double)(Ixk[i+2])*Sclp1*0.01;
    }
  }
  /* this reads in first record in mdb */
  /* sets Sclp and Sclp1               */
  /* gets min and max arrays           */
  else {  

    if(!(mdb=fopen("fort.1", "r"))) {
      STDERR("Couldn't open fort.1");
      exit(1);
    }

    newrecord(mdb);

    ftget_int(mdb,istep);
    ftget_int(mdb,&n);

    for (i=0; i<n; i++)       /* getting raw coordinates */
      ftget_short(mdb,&(Ixk[i]));

    ftget_int(mdb,&nat3p);     /* old vs new style coords (see getix2) */
    ftget_float(mdb,&Sclp);    /* scale to multiply coordinates        */
    Sclp1=1.0/Sclp;           
    ftget_float(mdb,&Sclw);    /* something not used                   */

    /* getting min and max for x, y, & z */
    for (i=0; i<3; i++)
      ftget_float(mdb, &min[i]);
    for (i=0; i<3; i++) 
      ftget_float(mdb, &max[i]);

    /* expands coordinates into floats and writes into file record */
     for (i=0; i<n; i+=3)  {       /* expands coordinates into floats */
      atom_record * a = f->atoms[i/3];

      a->x = (double)(Ixk[i])*Sclp1*0.01;
      a->y = (double)(Ixk[i+1])*Sclp1*0.01;
      a->z = (double)(Ixk[i+2])*Sclp1*0.01;
    }
    endrecord(mdb);   

    after_first_record++;
  }  



  return 0;
}


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



bool BrittReadTrajWithSkip 
(long times, bool skip, file_records * f, long * istepP, double * min, double * max)
{
  float x[30000] ;
  long i;
  static short ix[30000];
  static char  ixb[30000];
  static float xmax[3], xmin[3];
  int istep,AtEOF ;

  fprintf(stderr,"ReadTrajWithSkip(): times=%d skip=%1d\n",times,skip);
  for (i=0; i<times; i++) {
    fprintf(stderr,"Before get_mdb_coords(): istep=%ld\n",istep);
    AtEOF = get_mdb_coords(f,&istep,xmin,xmax);
    fprintf(stderr,"After  get_mdb_coords(): istep=%ld\n",istep);
    
    *istepP = istep ;
    min[0] = xmin[0]; min[1] = xmin[1]; min[2] = xmin[2];
    max[0] = xmax[0]; max[1] = xmax[1]; max[2] = xmax[2];
    
    if ( AtEOF == EOF )  {
      STDERR("The trajectory is exhausted.\n");
      return false ;
    }
  }
  return true;
}

/*............................................................................*/

bool BrittSkip 
(long times, file_records * f, int * istepP)
{
  long i;
  static float xmax[3], xmin[3];
  int istep,AtEOF ;

  for (i=0; i<times; i++) {
    fprintf(stderr,"BrittSkip(): before get_mdb_coords() with istep=%ld\n",*istepP);
    if ( get_mdb_coords(f,istepP,xmin,xmax) == EOF )  {
      STDERR("The trajectory is exhausted.\n");
      return false ;
    }
  }
  return true;
}

/* .......................................................................... */