Merge branch 'leadoff'

[physik/posic.git] / pair_correlation_calc.c

diff --git a/pair_correlation_calc.c b/pair_correlation_calc.c
index efe81ad..26db76b 100644 (file)
--- a/pair_correlation_calc.c
+++ b/pair_correlation_calc.c
@@ -19,7 +19,7 @@
 int usage(char *prog) {
 
        printf("\nusage:\n");
-       printf("  %s <save file> <dr>\n\n",prog);
+       printf("  %s <dr> <save file 1> [<save file 2> ...]\n\n",prog);
 
        return -1;
 }
@@ -28,56 +28,81 @@ int main(int argc,char **argv) {
 
        t_moldyn moldyn;
        int ret;
-       double *stat;
+       double *stat,*total;
        int slots;
-       int i;
+       int i,j;
        double dr;
        int fd;
+       unsigned char first;
 
-       if(argc!=3) {
+       if(argc<3) {
                usage(argv[0]);
                return -1;
        }
 
-       memset(&moldyn,0,sizeof(t_moldyn));
+       dr=atof(argv[1]);
+
+       first=1;
+       stat=NULL;
+       total=NULL;
+
+       for(j=2;j<argc;j++) {
+
+               memset(&moldyn,0,sizeof(t_moldyn));
+
+               printf("[pair corr calc] reading save file ...\n");
+               ret=moldyn_read_save_file(&moldyn,argv[j]);
+               if(ret) {
+                       printf("[pair corr calc] exit!\n");
+                       return ret;
+               }
+
+               //moldyn.cutoff*=2;
+               //moldyn.cutoff_square*=4;
+               moldyn.cutoff=6.0;
+               moldyn.cutoff_square=36.0;
+
+               slots=moldyn.cutoff/dr;
+               printf("[pair corr calc]\n");
+               printf("  slots: %d\n",slots);
+               printf("  cutoff: %f\n",moldyn.cutoff);
+               printf("  dr: %f\n",dr);
+
+               if(first) {
+                       stat=(double *)malloc(3*slots*sizeof(double));
+                       total=(double *)malloc(3*slots*sizeof(double));
+                       first=0;
+                       if(stat==NULL) {
+                               perror("[pair corr calc] alloc mem (stat)");
+                               return -1;
+                       }
+                       if(total==NULL) {
+                               perror("[pair corr calc] alloc mem (toal)");
+                               return -1;
+                       }
+                       memset(total,0,3*slots*sizeof(double));
+               }
+
+               /* link cell init */
+               link_cell_init(&moldyn,VERBOSE);
+
+               calculate_pair_correlation(&moldyn,dr,stat);
+
+               for(i=0;i<3*slots;i++)
+                       total[i]+=stat[i];
 
-       printf("[pair corr calc] reading save file ...\n");
-       ret=moldyn_read_save_file(&moldyn,argv[1]);
-       if(ret) {
-               printf("[pair corr calc] exit!\n");
-               return ret;
        }
 
-       moldyn.cutoff*=2;
-       moldyn.cutoff_square*=4;
-
-       dr=atof(argv[2]);
-       slots=moldyn.cutoff/dr;
-       printf("[pair corr calc]\n");
-       printf("  slots: %d\n",slots);
-       printf("  cutoff: %f\n",moldyn.cutoff);
-       printf("  dr: %f\n",dr);
-
-       stat=(double *)malloc(3*slots*sizeof(double));
-       if(stat==NULL) {
-               perror("[pair corr calc] alloc mem");
-               return -1;
-       }
-
-       /* link cell init */
-       link_cell_init(&moldyn,VERBOSE);
-
-       calculate_pair_correlation(&moldyn,dr,stat);
-
        fd=open("pair_corr_func.txt",
                O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
        dprintf(fd,"#r #ab #aa #bb\n");
        for(i=0;i<slots;i++)
                dprintf(fd,"%f %f %f %f\n",
-                       i*dr,stat[i],stat[slots+i],stat[2*slots+i]);
+                       i*dr,total[i],total[slots+i],total[2*slots+i]);
        close(fd);
 
        free(stat);
+       free(total);
 
        moldyn_free_save_file(&moldyn);