added pair correlation calc code
[physik/posic.git] / pair_correlation_calc.c
1 /*
2  * calcultae pair correlation function
3  *
4  * author: frank.zirkelbach@physik.uni-augsburg.de
5  *
6  */
7
8 #define _GNU_SOURCE
9 #include <stdio.h>
10 #include <stdlib.h>
11 #include <unistd.h>
12 #include <string.h>
13 #include <sys/types.h>
14 #include <sys/stat.h>
15 #include <fcntl.h>
16
17 #include "moldyn.h"
18
19 int usage(char *prog) {
20
21         printf("\nusage:\n");
22         printf("  %s <save file> <dr>\n\n",prog);
23
24         return -1;
25 }
26
27 int main(int argc,char **argv) {
28
29         t_moldyn moldyn;
30         int ret;
31         double *stat;
32         int slots;
33         int i;
34         double dr;
35         int fd;
36
37         if(argc!=3) {
38                 usage(argv[0]);
39                 return -1;
40         }
41
42         ret=moldyn_read_save_file(&moldyn,argv[1]);
43         if(ret) {
44                 printf("[pair corr calc] exit!\n");
45                 return ret;
46         }
47
48         dr=atof(argv[2]);
49         slots=(int)(moldyn.cutoff/dr);
50
51         stat=(double *)malloc(3*slots*sizeof(double));
52         if(stat==NULL) {
53                 perror("[pair corr calc] alloc mem");
54                 return -1;
55         }
56
57         calculate_pair_correlation(&moldyn,dr,stat);
58
59         fd=open("pair_corr_func_ab.txt",
60                 O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
61         dprintf(fd,"# type a - type b bonds\n");
62         for(i=0;i<slots;i++)
63                 dprintf(fd,"%f %f\n",i*dr,stat[i]);
64         close(fd);
65                 
66         fd=open("pair_corr_func_aa.txt",
67                 O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
68         dprintf(fd,"# type a - type a bonds\n");
69         for(i=0;i<slots;i++)
70                 dprintf(fd,"%f %f\n",i*dr,stat[slots+i]);
71         close(fd);
72                 
73         fd=open("pair_corr_func_bb.txt",
74                 O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
75         dprintf(fd,"# type a - type b bonds\n");
76         for(i=0;i<slots;i++)
77                 dprintf(fd,"%f %f\n",i*dr,stat[2*slots+i]);
78         close(fd);
79                 
80         free(stat);
81
82         return 0;
83 }