]> hackdaworld.org Git - physik/computational_physics.git/blob - mc_int.c
cdf9ad3cbbe6fc0c74e7fc751e97fa7c13e05278
[physik/computational_physics.git] / mc_int.c
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <time.h>
5
6 #define N 1000000
7 #define R1 5
8 #define R2 3
9 #define RHO 1
10 #define D 6
11
12 //#define DEBUG
13
14 int main(int argc,char **argv) {
15         double R12,R22,D2,D_2,YZ2;
16         double x[4],hlp[4],var[4];
17         double vol,vol_exakt,a1,a2;
18         int i,n;
19         double y_z_max,x_max;
20
21         n=N;
22         R12=R1*R1;
23         R22=R2*R2;
24         x_max=R1+R2-D;
25         a2=(R22+D*D-R12)/(2*D);
26         a1=D-a2;
27         y_z_max=2*sqrt(R12-a1*a1);
28         vol=y_z_max*y_z_max*x_max;
29         vol_exakt=M_PI*(R1-a1)*(R1-a1)*(3*R1-(R1-a1))/3;
30         vol_exakt+=M_PI*(R2-a2)*(R2-a2)*(3*R2-(R2-a2))/3;
31
32         srand(time(NULL));
33
34         /* init */
35         for(i=0;i<4;i++) x[i]=0;
36
37         while(n--) {
38                 hlp[0]=1.;
39                 hlp[1]=1.*rand()/RAND_MAX*x_max+D-R2;
40                 hlp[2]=1.*rand()/RAND_MAX*y_z_max-(y_z_max/2);
41                 hlp[3]=1.*rand()/RAND_MAX*y_z_max-(y_z_max/2);
42 #ifdef DEBUG
43                 printf("%f %f %f %f\n",hlp[0],hlp[1],hlp[2],hlp[3]);
44 #endif
45                 YZ2=hlp[2]*hlp[2]+hlp[3]*hlp[3];
46                 D2=(hlp[1]*hlp[1])+YZ2;
47                 D_2=(D-hlp[1])*(D-hlp[1])+YZ2;
48
49                 if((D2<=R12)&&(D_2<=R22)) {
50                         for(i=0;i<4;i++) {
51                                 hlp[i]*=RHO;
52                                 x[i]+=hlp[i];
53                                 var[i]+=hlp[i]*hlp[i];
54                         }
55                 }
56         }
57         for(i=0;i<4;i++) {
58                 x[i]/=N;
59                 var[i]=sqrt((var[i]/N-x[i]*x[i])/N)*vol;
60                 x[i]*=vol;
61                 if(i>0) x[i]/=vol_exakt;
62                 printf("x[%d] = %f +/- %f\n",i,x[i],var[i]);
63         }
64         printf("vol_exakt = %f\n",vol_exakt);
65
66         return 1;
67 }