added mc_init

diff --git a/.cvsignore b/.cvsignore
index de08a3b8b637089616476cab181a784ab0050e00..d6485d3edb9b910eb770b4fe45b6aca2b821539f 100644 (file)
--- a/.cvsignore
+++ b/.cvsignore
@@ -10,3 +10,4 @@ bessel_1
 bessel_2
 check_rand
 kettenbruchentwicklung
+mc_init

diff --git a/Makefile b/Makefile
index 55038119c23cee666c2badd041e5ae895ba89b13..0af42366cdb042edf95631ff7ee6d5d6440de76c 100644 (file)
--- 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 (file)
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;
+}