#include <math.h>  
#include <stdio.h> 
#include <stdlib.h>
#define k 8.61738569226e-5 // eV/K
#define N 100

double M=1;
double Omega=1;
double Inter=0;
double A=0.025;
double T=294.15;
double INT_MAX1=0.31;
double INT_MAX2=0.31;

double tau(double x1[],double x2[],double x3[])
{
	int i,j;
	double rij=0,interaction=0,uij=0;
	for (i=1;i<=N;i++)
	{
		for (j=1;j<=N;j++)
		{
			rij=(pow((x1[i]-x1[j]),2)+pow((x2[i]-x2[j]),2)+pow((x3[i]-x3[j]),2));
			interaction+=Inter*0.5*rij*exp(-sqrt(rij)/A);
		}
		uij+=M*Omega*Omega*0.5*(pow(x1[i],2)+pow(x2[i],2)+pow(x3[i],2))-interaction;
	}
	return uij;
}

void main()
{
	FILE *stream;
	char datei[30];
	printf("\nEingabedatei: ");
	scanf("%s",datei);
	if ((stream=fopen(datei,"r"))==NULL)
	{
		printf("Kann Datei nicht öffnen!!!");
	}
	fscanf(stream,"%lf",&M);
	fscanf(stream,"%lf",&Omega);
	fscanf(stream,"%lf",&Inter);
	fscanf(stream,"%lf",&A);
	fscanf(stream,"%lf",&T);
	fscanf(stream,"%lf",&INT_MAX1);
	fscanf(stream,"%lf",&INT_MAX2);
	int nz=1,j;
	double sm1=0,sm2=0,sqm1=0,sqm2=0,stabw=0,stabw1=0,stabw2=0,prozabw=0;
	double pa;
	double t1=0,t2=0;
	double s1[100000]={0};
	double s2[100000]={0};
	double x1[N+1]={0};
	double x2[N+1]={0};
	double x3[N+1]={0};
	double energie=0;
	//char datei[50];
	printf("ProzAbw: "); scanf("%lf",&pa);
	//printf("Ausgabedateiname: "); scanf("%s",datei);*/
	static double beta=1/(k*T);
	nz=1;
	do {
		for (j=1;j<=N;j++) {
			x1[j]=(double(rand())/(RAND_MAX/2)-1)*INT_MAX1;
			x2[j]=(double(rand())/(RAND_MAX/2)-1)*INT_MAX1;
			x3[j]=(double(rand())/(RAND_MAX/2)-1)*INT_MAX1;}
		t1=tau(x1,x2,x3);
		for (j=1;j<=N;j++) {
			x1[j]=x1[j]*INT_MAX2/INT_MAX1;
			x2[j]=x2[j]*INT_MAX2/INT_MAX1;
			x3[j]=x3[j]*INT_MAX2/INT_MAX1;}
		t2=tau(x1,x2,x3);
		//printf("%le\N",t);
		s1[nz]=t1*exp(-beta*t1);
		s2[nz]=exp(-beta*t2);
		sm1+=s1[nz];
		sm2+=s2[nz];
		sqm1+=(s1[nz]*s1[nz]);
		sqm2+=(s2[nz]*s2[nz]);
		if (nz>10) {
			stabw1=sqrt(((sqm1/(nz))-((sm1/nz)*(sm1/nz)))/(nz-1));
			stabw2=sqrt(((sqm2/(nz))-((sm2/nz)*(sm2/nz)))/(nz-1));
			stabw=sm1/sm2*sqrt(pow(stabw1*nz/sm1,2)+pow(stabw2*nz/sm2,2));
			energie=(1.5*N/beta+sm1/sm2);
			prozabw=(stabw/energie)*100;
			if (!(nz % 100)) printf("Ergebnis: %1.5le +- %1.5le; Prozentuale Abw.: %2.2lf bei n=%d\n",energie,stabw,prozabw,nz);
			}
			else prozabw=100;
		if (stabw<0.00001) prozabw=100;
		nz++;
	} while (prozabw > pa);
	printf("%1.5le %1.5le\nErgebnis: %1.5le +- %1.5le; Prozentuale Abw.: %2.2lf bei n=%d\n",sm1,sm2,energie,stabw,prozabw,(nz-1));
}