//moment d'inerite d'une sphere
#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

//fonction densite
double densite(double r, double R, double rho0);
//distance du centre de la sphere
double r(const double point[]);
//distance a l'axe de rotation au carre
double perp2(const double point[]);

int main() {
  //saisi des parametres
  double rho0;
  cout << "Quelle est la densite de la sphere au centre (en kg/m^3) ? ";
  cin >> rho0;
  double R;    //rayon de la sphère
  cout << "et le rayon de la sphere (en m) ? ";
  cin >> R;
  int div; //divisions de l'axe x
  cout << "En combien de partie voulez vous diviser le rayon de la sphere ? ";
  cin >> div;

  double dx = R/div;
  double dV = pow(dx,3);   //sous-volume

  double mInertie = 0.;
  int iter = 0;   //# iterations
  double point[3];   //point a l'interieur de la sphere
  //on calcule le moment d'inertie d'un octant de la sphere 
  point[0] = dx/2.;
  for (int i=0; i<div; i++) {
    point[1] = dx/2.;
    for (int j=0; j<div; j++) {
      point[2] = dx/2.;
      for (int k=0; k<div; k++) {
      	//d'abord il faut verifier si le point est a l'interieur
        if (r(point)>R) break;
	    mInertie += densite(r(point),R,rho0) * dV * perp2(point) ;
	    point[2] += dx;
	    iter++;
      }
      point[1] += dx;
    }
    point[0] += dx;
  }
  
  cout << "Moment d'inertie de l'octant calcule : " << setw(10) << mInertie
       << " kg m^2" << endl
       << "et nombre d'iterations computees : " << iter << endl;
  
  cout << "Le moment d'inertie est egal a "
       << setw(15) << setprecision(10) << scientific << 8.*mInertie
       << " kg m^2" << endl;

  return 0;
}

//fonction densite
double densite(double r, double R, double rho0) {
  return rho0*exp(-r/R);
}

//distance du centre de la sphere
double r(const double point[]) {
  return sqrt(pow(point[0],2) + pow(point[1],2) + pow(point[2],2));
}

//distance a l'axe de rotation au carre
double perp2(const double point[]) {
  return pow(point[0],2) + pow(point[1],2);
}