+/*
+ * 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;
+}
+