From facb6c3b0ab8d6a6dc1c9e15298fb2ec5e0317f5 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Tue, 6 Oct 2009 12:45:26 +0200
Subject: [PATCH] pc calc hopefully corrected now

---
 vasp_tools/pc_calc.c | 181 +++++++++++++++++++++++++++++++------------
 1 file changed, 132 insertions(+), 49 deletions(-)

diff --git a/vasp_tools/pc_calc.c b/vasp_tools/pc_calc.c
index f37506f..db12f50 100644
--- a/vasp_tools/pc_calc.c
+++ b/vasp_tools/pc_calc.c
@@ -43,30 +43,76 @@ int get_line(int fd,char *line,int max) {
         }
 }
 
+int get_max3_vals(char *buf,int amount,double *v1, double *v2,double *v3) {
+
+	int i,j,k;
+	char temp[128];
+
+	i=0;
+
+	for(k=0;k<amount;k++) {
+		while(buf[i]==' ')
+			i++;
+		j=0;
+		while((buf[i+j]!=' ')&&(buf[i+j]!=0))
+			j++;
+		memcpy(temp,buf+i,j);
+		temp[j]=0;
+		if(k==0) *v1=atof(temp);
+		if(k==1) *v2=atof(temp);
+		if(k==2) *v3=atof(temp);
+		i+=j;
+	}
+
+	return 0;
+}
+
+int get_2_ints(char *buf,int *i1,int *i2) {
+
+	int i,j,k;
+	char temp[128];
+
+	i=0;
+
+	for(k=0;k<2;k++) {
+		while(buf[i]==' ')
+			i++;
+		j=0;
+		while((buf[i+j]!=' ')&&(buf[i+j]!=0))
+			j++;
+		memcpy(temp,buf+i,j);
+		temp[j]=0;
+		if(k==0) *i1=atoi(temp);
+		if(k==1) *i2=atoi(temp);
+		i+=j;
+	}
+
+	return 0;
+}
+
 int main(int argc,char **argv) {
 
 	t_atom *atom;
 	double *aslot,*bslot,*cslot;
 	int cnt,count,fcnt,slots;
 	int fd;
-	char buf[256],*wptr;
+	char buf[256];
 	int i,j,k;
-	double dx2,dy2,dz2,dist,norm;
-	double sx,sy,sz;
+	int nsi,nc,ntot;
+	double dx,dy,dz,dist,norm;
+	double dxt,dyt,dzt;
+	double X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3;
+	double scale;
 
 	atom=NULL;
 
-	if(argc<5) {
-		printf("usage: %s sx sy sz file1 file2 ...\n",argv[0]);
+	if(argc<1) {
+		printf("usage: %s file1 file2 ...\n",
+		       argv[0]);
 		return -1;
 	}
 
-	sx=atof(argv[1]);
-	sy=atof(argv[2]);
-	sz=atof(argv[3]);
-
 	// prepare pc
-
 	slots=MAXR/DELTA;
 	aslot=malloc(slots*sizeof(double));
 	if(aslot==NULL) {
@@ -86,15 +132,14 @@ int main(int argc,char **argv) {
 		return -1;
 	}
 	memset(cslot,0,slots*sizeof(double));
-
 	printf("i: allocated 3 times %d slots ...\n",slots);
 
 	// use all given files ...
 	printf("using files:\n");
-	for(fcnt=4;fcnt<argc;fcnt++)
-		printf(" %d: %s\n",fcnt-4,argv[fcnt]);
+	for(fcnt=1;fcnt<argc;fcnt++)
+		printf(" %d: %s\n",fcnt,argv[fcnt]);
 
-	for(fcnt=4;fcnt<argc;fcnt++) {
+	for(fcnt=1;fcnt<argc;fcnt++) {
 
 	fd=open(argv[fcnt],O_RDONLY);
 	if(fd<0) {
@@ -105,25 +150,57 @@ int main(int argc,char **argv) {
 	// first line
 	cnt=get_line(fd,buf,256);
 
+	// second line -> scale
+	cnt=get_line(fd,buf,256);
+	get_max3_vals(buf,1,&scale,NULL,NULL);
+	printf("scale: %f\n",scale);
+
+	// third line -> X
+	cnt=get_line(fd,buf,256);
+	get_max3_vals(buf,3,&X1,&X2,&X3);
+	printf("X: %f %f %f\n",X1,X2,X3);
+	
+	// 4th line -> Y
+	cnt=get_line(fd,buf,256);
+	get_max3_vals(buf,3,&Y1,&Y2,&Y3);
+	printf("Y: %f %f %f\n",Y1,Y2,Y3);
+	
+	// 5th line -> Z
+	cnt=get_line(fd,buf,256);
+	get_max3_vals(buf,3,&Z1,&Z2,&Z3);
+	printf("Z: %f %f %f\n",Z1,Z2,Z3);
+	
+	// 6th line -> nsi + nc = ntot
+	cnt=get_line(fd,buf,256);
+	get_2_ints(buf,&nsi,&nc);
+	ntot=nsi+nc;
+	printf("Si: %d - C: %d - tot: %d\n",nsi,nc,ntot);
+	
+	// 7th line
+	cnt=get_line(fd,buf,256);
+
+	// 8th line
+	cnt=get_line(fd,buf,256);
+
 	// read in atoms
 	count=0;
-	while(1) {
+	while(count<ntot) {
 		cnt=get_line(fd,buf,256);
-		if(cnt<=0)
-			break;
-
-		// there is (another) atom	
+		if(cnt<=0) {
+			printf("Should never happen!\n");
+			return cnt;
+		}
 
 		atom=realloc(atom,(count+1)*sizeof(t_atom));
-		wptr=strtok(buf," ");
-		atom[count].type=wptr[0];
+		if(count<nsi)
+			atom[count].type='S';
+		else
+			atom[count].type='C';
 
-		wptr=strtok(NULL," ");
-		atom[count].x=atof(wptr);
-		wptr=strtok(NULL," ");
-		atom[count].y=atof(wptr);
-		wptr=strtok(NULL," ");
-		atom[count].z=atof(wptr);
+		get_max3_vals(buf,3,
+		              &(atom[count].x),
+		              &(atom[count].y),
+		              &(atom[count].z));
 
 		count+=1;
 	}
@@ -131,28 +208,33 @@ int main(int argc,char **argv) {
 	printf("i: read in %d atoms ...\n",count);
 
 	// calc pc
-
 	for(i=0;i<count;i++) {
 		for(j=0;j<i;j++) {
-			dx2=atom[i].x-atom[j].x;
-			if(dx2>sx/2.0)
-				dx2-=sx;
-			else if(dx2<-sx/2.0)
-				dx2+=sx;
-			dx2*=dx2;
-			dy2=atom[i].y-atom[j].y;
-			if(dy2>sy/2.0)
-				dy2-=sy;
-			else if(dy2<-sy/2.0)
-				dy2+=sy;
-			dy2*=dy2;
-			dz2=atom[i].z-atom[j].z;
-			if(dz2>sz/2.0)
-				dz2-=sz;
-			else if(dz2<-sz/2.0)
-				dz2+=sz;
-			dz2*=dz2;
-			dist=sqrt(dx2+dy2+dz2);
+			dx=atom[i].x-atom[j].x;
+			if(dx>0.5)
+				dx-=1;
+			else if(dx<-0.5)
+				dx+=1;
+
+			dy=atom[i].y-atom[j].y;
+			if(dy>0.5)
+				dy-=1;
+			else if(dy<-0.5)
+				dy+=1;
+
+			dz=atom[i].z-atom[j].z;
+			if(dz>0.5)
+				dz-=1;
+			else if(dz<-0.5)
+				dz+=1;
+
+			dxt=dx*X1+dy*Y1+dz*Z1;
+			dyt=dx*X2+dy*Y2+dz*Z2;
+			dzt=dx*X3+dy*Y3+dz*Z3;
+
+			dist=sqrt(dxt*dxt+dyt*dyt+dzt*dzt);
+			dist*=scale;
+
 			if(dist>=MAXR)
 				continue;
 			k=dist/DELTA;
@@ -168,13 +250,14 @@ int main(int argc,char **argv) {
 		}
 	}
 
-		close(fd);
+	close(fd);
+
 	}
 
 	// normalization and output
-
 	for(i=1;i<slots;i++) {
 		norm=4*M_PI*(i*DELTA)*(i*DELTA)*DELTA;
+		norm*=(argc-1);
 		printf("pc: %f %f %f %f\n",i*DELTA,
 		                           aslot[i]/norm,
 		                           bslot[i]/norm,
-- 
2.20.1