/*
   Project: CP
   Jan Hamaekers
   1999

   File: GetPdb.c
*/

#include<stdlib.h>
#include<stdio.h> 
#include<math.h>
#include"mergesort.h"
struct Coord {
  double coord[3];
  int Z;
};
double GetKey(void *C) { return(((struct Coord*)C)->Z); }




void GetPdb(FILE *source, FILE *target, int max) {
  struct Coord **sequence;
  int i;
  int dummyint,dummyint2;
  char atomchar;
  char dummystr[80];
  double center[3][2] = {{0.,0.},{0.,0.},{0.,0.}};
  double Bohr = 0.52919;
  sequence = (struct Coord **)
    malloc(sizeof(struct Coord*)*max);
  for (i=0; i < max; i++) {
    sequence[i] = (struct Coord *)
      malloc(sizeof(struct Coord));
    fscanf(source,"%s %i %c %i %lg %lg %lg",dummystr,&dummyint,&atomchar,&dummyint2, &sequence[i]->coord[0], &sequence[i]->coord[1], &sequence[i]->coord[2]);
    if (i==0) {
      center[0][0] = sequence[i]->coord[0];
      center[0][1] = sequence[i]->coord[0];
      center[1][0] = sequence[i]->coord[1];
      center[1][1] = sequence[i]->coord[1];
      center[2][0] = sequence[i]->coord[2];
      center[2][1] = sequence[i]->coord[2];
      switch (atomchar) {
      case 'H':
      	sequence[i]->Z = 1;
      	break;
      case 'B':
      	sequence[i]->Z = 5;
      	break;
      case 'N':
      	sequence[i]->Z = 7;
      	break;
      }
    } else {
      if (sequence[i]->coord[0] < center[0][0]) center[0][0] = sequence[i]->coord[0];
      if (sequence[i]->coord[0] > center[0][1]) center[0][1] = sequence[i]->coord[0];
      if (sequence[i]->coord[1] < center[1][0]) center[1][0] = sequence[i]->coord[1];
      if (sequence[i]->coord[1] > center[1][1]) center[1][1] = sequence[i]->coord[1];
      if (sequence[i]->coord[2] < center[2][0]) center[2][0] = sequence[i]->coord[2];
      if (sequence[i]->coord[2] > center[2][1]) center[2][1] = sequence[i]->coord[2];
      switch (atomchar) {
      case 'H':
      	sequence[i]->Z = 1;
      	break;
      case 'B':
      	sequence[i]->Z = 5;
      	break;
      case 'N':
      	sequence[i]->Z = 7;
      	break;
      }
    }
  }
  naturalmergesort((void **)sequence,0,max-1,&GetKey);
  fprintf(target,"Lattice in Bohr\n%g %g, %g %g, %g %g\nCenter in Angstroem %g %g %g\n",center[0][0]/Bohr,center[0][1]/Bohr,center[1][0]/Bohr,center[1][1]/Bohr,center[2][0]/Bohr,center[2][1]/Bohr, (center[0][0]+center[0][1])/2,(center[1][0]+center[1][1])/2 , (center[2][0]+center[2][1])/2);
  dummyint = 0;
  dummyint2 = 0;
  for (i =0;i < max;i++) {
    if (sequence[i]->Z != dummyint) {
      fprintf(target,"Atom: Following Z(%i) BeforeMax(%i)\n",sequence[i]->Z,dummyint2);
      dummyint2=0;
    }
    dummyint2++;
    dummyint = sequence[i]->Z;
    fprintf(target,"%g %g %g\n", sequence[i]->coord[0], sequence[i]->coord[1], sequence[i]->coord[2]);
    free(sequence[i]);
  }
  free(sequence);
  fprintf(target,"Atom: BeforeMax(%i)\n",dummyint2);
}

int main(int argc, char** argv) {
  int max;
  FILE *source;
  FILE *target;
  if(argc < 4) {
    fprintf(stderr,"Not enough arguments\nUsage: GetPdb source target nmax\n");
    exit(1);
  }
  source = fopen(argv[1], "r");
  if(source == 0) {
    fprintf(stderr,"Could not open source\n");
    exit(1);
  }
  target = fopen(argv[2], "w");
  if(target == 0) {
    fprintf(stderr,"Could not open target\n");
    exit(1);
  }
  max = atoi(argv[3]);
  fprintf(stderr,"Input: source: %s \ntarget: %s\nmax %i\n", argv[1], argv[2],max);
  GetPdb(source,target,max);
  fclose(source);
  fclose(target);
  return 0;
}