X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fcomputational_physics.git;a=blobdiff_plain;f=mc_int.c;fp=mc_int.c;h=cdf9ad3cbbe6fc0c74e7fc751e97fa7c13e05278;hp=0000000000000000000000000000000000000000;hb=a649f7232da963470fd43e9a1a5384944ba5ed62;hpb=f5ca97914e79b4882001d2ae4517b203ee09fd6f

diff --git a/mc_int.c b/mc_int.c
new file mode 100644
index 0000000..cdf9ad3
--- /dev/null
+++ b/mc_int.c
@@ -0,0 +1,67 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+
+#define N 1000000
+#define R1 5
+#define R2 3
+#define RHO 1
+#define D 6
+
+//#define DEBUG
+
+int main(int argc,char **argv) {
+	double R12,R22,D2,D_2,YZ2;
+	double x[4],hlp[4],var[4];
+	double vol,vol_exakt,a1,a2;
+	int i,n;
+	double y_z_max,x_max;
+
+	n=N;
+	R12=R1*R1;
+	R22=R2*R2;
+	x_max=R1+R2-D;
+	a2=(R22+D*D-R12)/(2*D);
+	a1=D-a2;
+	y_z_max=2*sqrt(R12-a1*a1);
+	vol=y_z_max*y_z_max*x_max;
+	vol_exakt=M_PI*(R1-a1)*(R1-a1)*(3*R1-(R1-a1))/3;
+	vol_exakt+=M_PI*(R2-a2)*(R2-a2)*(3*R2-(R2-a2))/3;
+
+	srand(time(NULL));
+
+	/* init */
+	for(i=0;i<4;i++) x[i]=0;
+
+	while(n--) {
+		hlp[0]=1.;
+		hlp[1]=1.*rand()/RAND_MAX*x_max+D-R2;
+		hlp[2]=1.*rand()/RAND_MAX*y_z_max-(y_z_max/2);
+		hlp[3]=1.*rand()/RAND_MAX*y_z_max-(y_z_max/2);
+#ifdef DEBUG
+		printf("%f %f %f %f\n",hlp[0],hlp[1],hlp[2],hlp[3]);
+#endif
+		YZ2=hlp[2]*hlp[2]+hlp[3]*hlp[3];
+		D2=(hlp[1]*hlp[1])+YZ2;
+		D_2=(D-hlp[1])*(D-hlp[1])+YZ2;
+
+		if((D2<=R12)&&(D_2<=R22)) {
+			for(i=0;i<4;i++) {
+				hlp[i]*=RHO;
+				x[i]+=hlp[i];
+				var[i]+=hlp[i]*hlp[i];
+			}
+		}
+	}
+	for(i=0;i<4;i++) {
+		x[i]/=N;
+		var[i]=sqrt((var[i]/N-x[i]*x[i])/N)*vol;
+		x[i]*=vol;
+		if(i>0) x[i]/=vol_exakt;
+		printf("x[%d] = %f +/- %f\n",i,x[i],var[i]);
+	}
+	printf("vol_exakt = %f\n",vol_exakt);
+
+	return 1;
+}