//
//
//
//
// Lab. Calcolo II - Esempio di codice
//
//
//
//
//
//
// <![CDATA[
// Example code: montecarlo integration. Let's compute the volume of
// a sphere slice.
// (20060911 francesco.prelz@mi.infn.it)
#include <iostream>
#include <random>
#include <cmath>
int
main (int argc, char *argv[])
{
double radius = 10,rsquare = radius*radius,rcube = rsquare*radius,
rcube8 = 8*rcube;
double zmin = 3.5;
double volume;
double est_sigma;
double req_sigma = 1;
double x,y,z;
unsigned int count_in, count_tot;
count_in = 0;
count_tot = 0;
std::random_device rndgen;
std::uniform_real_distribution<double> rndist(-radius,radius);
for (count_tot=1;;count_tot++)
{
// Pick the random coordinates inside a cube of side 2*radius
x = rndist(rndgen);
y = rndist(rndgen);
z = rndist(rndgen);
if ( z > zmin && ((x*x + y*y + z*z) <= rsquare) )
{
count_in++;
}
if (count_in > 0 && count_tot > count_in)
{
// One binomial distribution sigma
est_sigma = std::sqrt((double)count_in)/(double)count_tot*
std::sqrt((double)(count_tot - count_in)/(double)count_tot) *
rcube8;
if (est_sigma <= req_sigma) break;
}
volume = (static_cast<double>(count_in)/static_cast<double>(count_tot))
* rcube8;
if (!(count_tot%500000)) std::cout << "picks == " << count_tot
<< " sigma == " << est_sigma
<< " value == " << volume << std::endl;
}
std::cout << "Resulting volume == " << volume << " +/- " << est_sigma << std::endl;
std::cout << "Number of random picks == " << count_tot << std::endl;
return 0;
}
]]></plaintext></span>