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 }