/* matfract is a conversion program that takes a matrix file concatenated with a
   fractional coordinates file as input from stdin and converts it to a standard
   fractional file on stdout.
   
   format of the matrix file : unformatted input (only sequence important)
   number of symmetry operations : nsymm
   for every symmetry operation  : M11 M12 M13 T1 M21 M22 M23 T2 M31 M32 M33 T3
   
   The symmetry operations performed are :
   / M11 M12 M13 \ / x \   / T1 \
   | M21 M22 M23 | | y | + | T2 |
   \ M31 M32 M33 / \ z /   \ T3 /
   with (x,y,z) as the fractional coordinates listed in the appended fractional
   coordinates file.
   
   The symmetry operations are only performed on the input coordinates (not on the
   generated ones from previous symmetry operations).
   Atoms that end up with fractional coordinates within 0.00001 of existing one
   (input and previously generated ones) are deleted in the output file.
   
   format of the fractional coordinates file:
   first line : natoms alpha beta gamma a-axis b-axis c-axis (rest of line ignored)
   second line : title or comment line
   rest : atomname fract(x)  fract(y)  fract(z)              (rest of line ignored)

   atomname = elementname that may have a number concatenated to it
   i.g. C3 or H34b are valid : become C and H respectively
   fract(x) fract(y) fract(z) : floating point fractional coordinates

   to compile on SUN : cc matfract.c -o matfract -O -lm

   typical use :
     cat matrix.txt infile.fract | matfract >outfile.fract
   
   if you want to make it straight into an xyz file for xmol you can issue:
     cat matrix.txt infile.fract | matfract | fract2xyz >outfile.xyz

   bugs : report to KANTERS@URVAX.URICH.EDU

   Enjoy,
   Rene Kanters
   
*/

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


#define MINDIST 0.00001

typedef double vector[3];
typedef vector matrix[3];

typedef struct symmetry_operation {
  matrix m;
  vector t;
} SymmOp;

typedef struct atom {
  char name[4];
  vector xyz;
} Atom;

double vecxvec(v1,v2) vector v1,v2; { return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];}

void symmatom(o,s,d) SymmOp *o;Atom *s,*d;
{ /* note that the atom pointers should NOT be the same !!!! */
  d->xyz[0]=o->t[0]+vecxvec(o->m[0],s->xyz);
  d->xyz[1]=o->t[1]+vecxvec(o->m[1],s->xyz);
  d->xyz[2]=o->t[2]+vecxvec(o->m[2],s->xyz);
  strcpy(d->name,s->name);
}

void scanvec(v) vector v;  {scanf("%lf%lf%lf",&v[0],&v[1],&v[2]);}

main(argc,argv)
int argc;
char *argv[];
{ long natoms,nsymm,i,j,k;
  Atom *atoms,*a;
  SymmOp *symmops,*s;
  char buf[256],line1[256],line2[256],*ch;

  /* read the symmetry operations */
  scanf("%ld",&nsymm);
  if(!(symmops=(SymmOp *)malloc(nsymm*sizeof(struct symmetry_operation)))) {
    fprintf(stderr,"Not enough space to allocate memory for symmetry operations/n");
    exit(0);
  }
  for (i=0,s=symmops;i<nsymm;i++,s++) {
    scanvec(s->m[0]);scanf("%lf",&s->t[0]);
    scanvec(s->m[1]);scanf("%lf",&s->t[1]);
    scanvec(s->m[2]);scanf("%lf",&s->t[2]);
  }
  
  scanf("%ld",&natoms); gets(line1);
  gets(line2);
  if(!(atoms=(Atom *)malloc((nsymm+1)*natoms*sizeof(struct atom)))) {
    fprintf(stderr,"Not enough space to allocate memory for the atoms/n");
    exit(0);
  }
  /* read the atoms themselves */  
  for(i=0;i<natoms;i++) {
    scanf("%s",buf);
    ch=buf; while(isalpha(*ch)) ch++; *ch=0; /* remove dangling numbers */
    strncpy(atoms[i].name,buf,3);
    scanvec(atoms[i].xyz);
  }
  
  /* let's generate all the atoms */
  for(i=0;i<nsymm;i++) {
    s=&symmops[i];
    for(j=0,k=(i+1)*natoms;j<natoms;j++,k++) symmatom(s,&atoms[j],&atoms[k]);
  }
  
  /* delete doubly occupied spaces and count newnatoms */
  k=natoms*(nsymm+1);
  for(i=1,a=atoms+1;i<k;i++,a++) {
    for(j=0;j<i;j++)
      if(atoms[j].name[0] && /* check first if not deleted before */
         fabs(atoms[j].xyz[0]-a->xyz[0])<MINDIST &&
         fabs(atoms[j].xyz[1]-a->xyz[1])<MINDIST &&
         fabs(atoms[j].xyz[2]-a->xyz[2])<MINDIST) {
           atoms[i].name[0]=0;
           break; /* make sure that we don't test this one any further */
         }
  }
  
  natoms=0;
  for(i=0;i<k;i++) if(atoms[i].name[0]) natoms++;
  
  printf("%ld",natoms); puts(line1);
  puts(line2);
  for(i=0,a=atoms;i<k;i++,a++)
    if(a->name[0])
      printf("%3s %15.6lg %15.6lg %15.6lg\n",a->name,a->xyz[0],a->xyz[1],a->xyz[2]);
}