/* concenxyz.c */

/* John H. Carpenter */
/* Created for Physics398B final project */

/* This program computes the concentration (number density) of the
 * two types of atoms present in the input lattice. For the calculation, 
 * each direction is cut up into bins having a width equal to the
 * lattice constant of the system
 */

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

int
main()
{
  FILE *infile,*outfile;
  char name[80];
  double xell[3],x[3],lconst;
  int type,btypeax[50],btypebx[50],btypeay[50],btypeby[50],
    btypeaz[50],btypebz[50],total=0,i,k[3];

  for (i=0;i<50;i++) { 
    btypeax[i] = 0; btypebx[i] = 0; 
    btypeay[i] = 0; btypeby[i] = 0; 
    btypeaz[i] = 0; btypebz[i] = 0;
  }

  printf("Enter input file:\n");
  scanf("%s",name);
  infile = fopen(name,"r");
  printf("Enter output file:\n");
  scanf("%s",name);
  outfile = fopen(name,"w");
  printf("Enter lattice constant:\n");
  scanf("%lf",&lconst);

  fscanf(infile,"%lf %lf %lf",&xell[0],&xell[1],&xell[2]);
  printf("a %g %g %g \n",xell[0],xell[1],xell[2]);
  
  while(fscanf(infile,"%d %lf %lf %lf",&type,&x[0],&x[1],&x[2]) != EOF){
    /* find layer number */
    for(i=0;i<3;i++){
      x[i] += xell[i]/2;
      k[i] = (int)floor(x[i]/lconst);
    }
    printf("b %d x-> %g %d y-> %g %d z-> %g %d\n",
	   type,x[0],k[0],x[1],k[1],x[2],k[2]);
    /* bin information */
    if (type==1) { 
      btypeax[k[0]]++; btypeay[k[1]]++; btypeaz[k[2]]++; total++;
    }
    else if (type==2) { 
      btypebx[k[0]]++; btypeby[k[1]]++; btypebz[k[2]]++; total++;
    }
  }

  
  for (i=0;i<50;i++) 
    fprintf(outfile,"%d %g %g %g %g %g %g\n",i,
	    btypeax[i]/(double)total,btypebx[i]/(double)total,
	    btypeay[i]/(double)total,btypeby[i]/(double)total,
	    btypeaz[i]/(double)total,btypebz[i]/(double)total);

  fclose(infile);
  fclose(outfile);

  return 0;
}