+/*
+ * pc_calc.c - calculate the pc
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+
+#include <math.h>
+
+
+#define MAXR 6.0
+#define DELTA 0.01
+
+typedef struct s_atom {
+ double x,y,z;
+ char type;
+} t_atom;
+
+int get_line(int fd,char *line,int max) {
+
+ int count,ret;
+
+ count=0;
+
+ while(1) {
+ if(count==max) return count;
+ ret=read(fd,line+count,1);
+ if(ret<=0) return ret;
+ if(line[count]=='\n') {
+ memset(line+count,0,max-count-1);
+ //line[count]='\0';
+ return count+1;
+ }
+ count+=1;
+ }
+}
+
+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];
+ int i,j,k;
+ 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<1) {
+ printf("usage: %s file1 file2 ...\n",
+ argv[0]);
+ return -1;
+ }
+
+ // prepare pc
+ slots=MAXR/DELTA;
+ aslot=malloc(slots*sizeof(double));
+ if(aslot==NULL) {
+ perror("slot a\n");
+ return -1;
+ }
+ memset(aslot,0,slots*sizeof(double));
+ bslot=malloc(slots*sizeof(double));
+ if(bslot==NULL) {
+ perror("slot a\n");
+ return -1;
+ }
+ memset(bslot,0,slots*sizeof(double));
+ cslot=malloc(slots*sizeof(double));
+ if(cslot==NULL) {
+ perror("slot a\n");
+ 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=1;fcnt<argc;fcnt++)
+ printf(" %d: %s\n",fcnt,argv[fcnt]);
+
+ for(fcnt=1;fcnt<argc;fcnt++) {
+
+ fd=open(argv[fcnt],O_RDONLY);
+ if(fd<0) {
+ perror("open file\n");
+ return fd;
+ }
+
+ // 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(count<ntot) {
+ cnt=get_line(fd,buf,256);
+ if(cnt<=0) {
+ printf("Should never happen!\n");
+ return cnt;
+ }
+
+ atom=realloc(atom,(count+1)*sizeof(t_atom));
+ if(count<nsi)
+ atom[count].type='S';
+ else
+ atom[count].type='C';
+
+ get_max3_vals(buf,3,
+ &(atom[count].x),
+ &(atom[count].y),
+ &(atom[count].z));
+
+ count+=1;
+ }
+
+ printf("i: read in %d atoms ...\n",count);
+
+ // calc pc
+ for(i=0;i<count;i++) {
+ for(j=0;j<i;j++) {
+ 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;
+ if((atom[i].type=='S')&&(atom[j].type=='S')) {
+ aslot[k]+=1;
+ }
+ else if((atom[i].type=='C')&&(atom[j].type=='C')) {
+ bslot[k]+=1;
+ }
+ else {
+ cslot[k]+=1;
+ }
+ }
+ }
+
+ 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,
+ cslot[i]/norm);
+ }
+
+ free(aslot);
+ free(bslot);
+ free(cslot);
+
+ free(atom);
+
+ return 0;
+
+}
+
projects / physik / posic.git / blobdiff