X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=fluctuation_calc.c;fp=fluctuation_calc.c;h=210f3829ca7c6b495511043170969e62f06374b6;hb=25c22fe95e80365056c6a7fadc548119360ca8ce;hp=0000000000000000000000000000000000000000;hpb=0fbd532d8ddce9848df592ed586ffba439369284;p=physik%2Fposic.git

diff --git a/fluctuation_calc.c b/fluctuation_calc.c
new file mode 100644
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;
+}
+