added pair correlation calc code

[physik/posic.git] / pair_correlation_calc.c

diff --git a/pair_correlation_calc.c b/pair_correlation_calc.c
new file mode 100644 (file)
index 0000000..2f93ab7
--- /dev/null
+++ b/pair_correlation_calc.c
@@ -0,0 +1,83 @@
+/*
+ * calcultae pair correlation function
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+
+#include "moldyn.h"
+
+int usage(char *prog) {
+
+       printf("\nusage:\n");
+       printf("  %s <save file> <dr>\n\n",prog);
+
+       return -1;
+}
+
+int main(int argc,char **argv) {
+
+       t_moldyn moldyn;
+       int ret;
+       double *stat;
+       int slots;
+       int i;
+       double dr;
+       int fd;
+
+       if(argc!=3) {
+               usage(argv[0]);
+               return -1;
+       }
+
+       ret=moldyn_read_save_file(&moldyn,argv[1]);
+       if(ret) {
+               printf("[pair corr calc] exit!\n");
+               return ret;
+       }
+
+       dr=atof(argv[2]);
+       slots=(int)(moldyn.cutoff/dr);
+
+       stat=(double *)malloc(3*slots*sizeof(double));
+       if(stat==NULL) {
+               perror("[pair corr calc] alloc mem");
+               return -1;
+       }
+
+       calculate_pair_correlation(&moldyn,dr,stat);
+
+       fd=open("pair_corr_func_ab.txt",
+               O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+       dprintf(fd,"# type a - type b bonds\n");
+       for(i=0;i<slots;i++)
+               dprintf(fd,"%f %f\n",i*dr,stat[i]);
+       close(fd);
+               
+       fd=open("pair_corr_func_aa.txt",
+               O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+       dprintf(fd,"# type a - type a bonds\n");
+       for(i=0;i<slots;i++)
+               dprintf(fd,"%f %f\n",i*dr,stat[slots+i]);
+       close(fd);
+               
+       fd=open("pair_corr_func_bb.txt",
+               O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+       dprintf(fd,"# type a - type b bonds\n");
+       for(i=0;i<slots;i++)
+               dprintf(fd,"%f %f\n",i*dr,stat[2*slots+i]);
+       close(fd);
+               
+       free(stat);
+
+       return 0;
+}