double Potenziale(double x, double y, double z) {
  // inizializzazione dell costanti
  int Npunti=10000;
  TRandom generatore(0);
  double G =    6.6742E-11 ;  // costante di gravitazione universale
  double massa= 5.9723E+24 ;  // massa della terra
  double raggio=6.5781E+06 ;  // raggio della terra
  double rho  =               // densità della terra
    massa/(4.*TMath::Pi()*pow(raggio,3)/3.);
  // volume in cui vengono distribuiti
  // i punti casuali:
  double volume=pow(2.*raggio,3);
  // somma delle valutazione della funzione 
  double potenziale=0; 
  
  for ( int i=0; i<Npunti; i++ ) {
    // scelgo un punto a caso all'interno di un cubo 
    // contenente al terra
    double x1 = generatore.Uniform(-raggio,raggio); 
    double y1 = generatore.Uniform(-raggio,raggio); 
    double z1 = generatore.Uniform(-raggio,raggio); 
    // all'interno del volume cubico, l'integranda
    // -         0     al di fuori della sfera
    // - -G*rho/(r-r') all'interno della sfera
    if ( (x1*x1+y1*y1+z1*z1) < (raggio*raggio) ) {
      double r=sqrt( pow(x-x1,2)+pow(y-y1,2)+pow(z-z1,2) );
      potenziale-= G*rho/r;
    }
  }
  // Resituisco l'integrale espresso come valor medio
  // della funzione moltiplicato per il volume della
  // regione di integrazione.
  return volume*(potenziale/Npunti);
}