From: hackbard Date: Thu, 26 Feb 2004 13:07:37 +0000 (+0000) Subject: added mc_init X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=a649f7232da963470fd43e9a1a5384944ba5ed62;p=physik%2Fcomputational_physics.git added mc_init --- 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 +#include +#include +#include + +#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; +}