|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
|
|
#define THETA 0.5
|
|
#define H 2.0
|
|
#define N_MC 100000
|
|
#define SQR(x) ((x) * (x))
|
|
|
|
double input();
|
|
double raggio(double y_p, double theta);
|
|
double V_teorico(double h, double R, double y_p, double r);
|
|
double V_MC(double h, double R, double y_p);
|
|
|
|
int main() {
|
|
double y_p;
|
|
double R = raggio(H, THETA);
|
|
double y_max = input();
|
|
double delta_y = (H + y_max) / 50;
|
|
|
|
FILE *output = fopen("volume.dat", "w");
|
|
for(y_p = -H; y_p < y_max; y_p += delta_y) {
|
|
double r = raggio(y_p, THETA);
|
|
double Vth = V_teorico(H, R, y_p, r);
|
|
double Vmc = V_MC(H, R, y_p);
|
|
|
|
fprintf(output, "%.3lf %.3lf %.3lf\n", y_p, Vmc, Vth);
|
|
}
|
|
fclose(output);
|
|
|
|
return 0;
|
|
}
|
|
|
|
double input() {
|
|
double y_max;
|
|
do {
|
|
printf("Inserire il valore massimo di y da considerare: ");
|
|
scanf("%lf", &y_max);
|
|
}
|
|
while(y_max > 0 || y_max < -H);
|
|
|
|
return y_max;
|
|
}
|
|
|
|
double raggio(double y_p, double theta) {
|
|
return fabs(y_p) * tan(theta);
|
|
}
|
|
|
|
double V_teorico(double h, double R, double y_p, double r) {
|
|
return M_PI / 3 * (SQR(R) * h - SQR(r) * fabs(y_p));
|
|
}
|
|
|
|
double V_MC(double h, double R, double y_p) {
|
|
int i;
|
|
int Np = 0;
|
|
double V_paral = SQR(2 * R) * H;
|
|
|
|
for(i = 0; i < N_MC; i++) {
|
|
double x = (drand48() - 0.5) * 2 * R;
|
|
double y = -drand48() * H;
|
|
double z = (drand48() - 0.5) * 2 * R;
|
|
|
|
if(y < y_p) {
|
|
double r = raggio(y, THETA);
|
|
if(SQR(x) + SQR(z) < SQR(r)) {
|
|
Np++;
|
|
}
|
|
}
|
|
}
|
|
|
|
return V_paral * Np / (double) N_MC;
|
|
}
|