pc calc hopefully corrected now

diff --git a/vasp_tools/pc_calc.c b/vasp_tools/pc_calc.c
index f37506f6afe0cec749ff8859a3c64cb700281be0..db12f50c6e1a8f70cb758a6f17199a602e92d049 100644 (file)
--- 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,