added fluctuation calc code

[physik/posic.git] / fluctuation_calc.c

diff --git a/fluctuation_calc.c b/fluctuation_calc.c
new file mode 100644 (file)
index 0000000..210f382
--- /dev/null
+++ b/fluctuation_calc.c
@@ -0,0 +1,121 @@
+/*
+ * calcultae fluctuations and related values
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+
+#define AMU            1.6605388628e-27
+#define METER          1e10
+#define SECOND         1e15
+#define KILOGRAM       (1.0/AMU)
+#define NEWTON         (METER*KILOGRAM/(SECOND*SECOND))
+#define JOULE          (NEWTON*METER)
+#define K_BOLTZMANN    (1.380650524e-23*METER*NEWTON)
+#define K_B2           (K_BOLTZMANN*K_BOLTZMANN)
+#define EV             (1.6021765314e-19*METER*NEWTON)
+
+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') {
+                        line[count]='\0';
+                        return count+1;
+                }
+                count+=1;
+        }
+}
+
+int main(int argc,char **argv) {
+
+       double from,to;
+       int fd;
+       char file[256];
+       char buf[256];
+       double e,e2,de2,val,time;
+       char *ptr;
+       int i,count,lnr;
+       double np,m,t;
+
+       if(argc<5) {
+               printf("usage:\n");
+               printf("%s <file> <from> <to> <linenumber>\n",argv[0]);
+               printf("\n");
+               printf("optional add: <nr particles> <mass> <temperature>");
+               printf("\n");
+               return -1;
+       }
+       
+       strncpy(file,argv[1],255);
+       from=atof(argv[2]);
+       to=atof(argv[3]);
+       lnr=atoi(argv[4]);
+       np=0;
+       m=0.0;
+       t=0.0;
+
+       if(argc==8) {
+               np=atoi(argv[5]);
+               m=atof(argv[6])*AMU*np;
+               t=atof(argv[7])+273.0;
+       }
+
+       fd=open(file,O_RDONLY);
+       if(fd<2) {
+               perror("fd open");
+               return -1;
+       }
+
+       count=0;
+       e=0.0;
+       e2=0.0;
+
+       while(1) {
+               if(get_line(fd,buf,255)<=0)
+                       break;
+               ptr=strtok(buf," ");
+               time=atof(ptr);
+               if((time>=from)&(time<=to)) {
+                       count+=1;
+                       for(i=1;i<lnr;i++)
+                               ptr=strtok(NULL," ");
+                       val=atof(ptr);
+                       e+=val;
+                       e2+=(val*val);
+               }
+       }
+
+       de2=e2/count-e*e/(count*count);
+       printf("--> fluctuation [eV/atom]: %.20f\n",de2);
+       if(argc==8) {
+               de2*=(np*np)*(EV*EV);
+               printf("    specific heat capacity (T=%f K) [J/(kg K)]:\n",
+                      t);
+               printf("    nve: %f\n",
+                      (1.5*np*K_BOLTZMANN/(1.0-de2/(1.5*np*K_B2*t*t))/
+                      (m*JOULE)));
+               printf("    nvt: %f\n",
+                      (de2/(K_BOLTZMANN*t*t)+1.5*np*K_BOLTZMANN)/
+                      (m*JOULE));
+       }
+
+       close(fd);
+
+       return 0;
+}
+