From a649f7232da963470fd43e9a1a5384944ba5ed62 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Thu, 26 Feb 2004 13:07:37 +0000
Subject: [PATCH] added mc_init

---
 .cvsignore |  1 +
 Makefile   |  5 +++-
 mc_int.c   | 67 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 72 insertions(+), 1 deletion(-)
 create mode 100644 mc_int.c

diff --git a/.cvsignore b/.cvsignore
index de08a3b..d6485d3 100644
--- a/.cvsignore
+++ b/.cvsignore
@@ -10,3 +10,4 @@ bessel_1
 bessel_2
 check_rand
 kettenbruchentwicklung
+mc_init
diff --git a/Makefile b/Makefile
index 5503811..0af4236 100644
--- a/Makefile
+++ b/Makefile
@@ -5,7 +5,7 @@ CFLAGS = -O3 -Wall
 LIBS = -L/usr/lib -lm
 
 API = g_plot.o general.o
-OBJS = newton zentral homogen integral-1_2 integral-2_2 polynom_interpolation kettenbruchentwicklung bessel_1 bessel_2 check_rand
+OBJS = newton zentral homogen integral-1_2 integral-2_2 polynom_interpolation kettenbruchentwicklung bessel_1 bessel_2 check_rand mc_int
 
 all: $(OBJS)
 
@@ -39,6 +39,9 @@ bessel_2: $(API)
 check_rand: $(API)
 	$(CC) $(CFLAGS) -o $@ $(API) $(LIBS) check_rand.c
 
+mc_int: $(API)
+	$(CC) $(CFLAGS) -o $@ $(API) $(LIBS) mc_int.c
+
 clean:
 	rm $(API) $(OBJS)
 
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;
+}
-- 
2.39.2