/*
	November 7,2011
	GET dR^2 vs dt plot from xyz values of a single particle
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

main(){
	FILE *fpr0,*fpr1,*fpw0;
	char buffer[256];
	float ts0,ts1;
	float r0[3],r1[3];
	float timestep = 0.05;

	fpr0=fopen("SingPart.xyz","r");
	fpr1=fopen("SingPart.xyz","r");
	fpw0=fopen("msd.out","w");

	int dn=0;
	int i,k;

	while(1){
		fgets(buffer,254,fpr0);
		fgets(buffer,254,fpr1);

		for(i=0;i<dn;i++)	fgets(buffer,254,fpr1);

		if(feof(fpr1)!=0){
			break;
		}
		else{
			double dr2=0;
			double ndr2=0;
			while(1){
				fgets(buffer,254,fpr0);
				sscanf(buffer,"%f %f %f %f",&ts0,&r0[0],&r0[1],&r0[2]);
				fgets(buffer,254,fpr1);
				sscanf(buffer,"%f %f %f %f",&ts1,&r1[0],&r1[1],&r1[2]);
				if(feof(fpr1)!=0){
					rewind(fpr0);
					rewind(fpr1);
					fprintf(fpw0,"%f %f\n",dn*timestep,dr2/ndr2);
					dn++;
					fflush(fpw0);
					break;
				}
				else{
					float dr2_t=0;
					for(k=0;k<3;k++) dr2_t+= (r1[k]-r0[k])*(r1[k]-r0[k]);
					dr2+=dr2_t;
					ndr2++;
				}

			}
		}

	}
	

	fclose(fpr0);
	fclose(fpr1);
	fclose(fpw0);
}