//calcul du champ gravitationnel genere par un ellipsoide dans le point X
#include <iostream>
#include <cmath>

using namespace std;

bool verifierVolume(double x[], double a[]);
void fieldComp(double xsub[], double xp[], double f[]);

//constante globale
const double GN = 6.673e-11;   //m^3 kg^-1 s^-2
  
int main() {
  //saisi des parametres
  double a[3];
  cout << "Champ gravitationnel genere par un ellipsoid dans le point X"
       << endl << endl;
  cout << "Entrez les dimensions de l'ellipsoide (axes en metres) :"
       << endl;
  cout << "a = ";   cin >> a[0];
  cout << "b = ";   cin >> a[1];
  cout << "c = ";   cin >> a[2];

  double M = 10e20;
  cout << "Entrez le masse de l'objet (~10e20 kg) :" << endl;
  cin >> M;

  double xp[3];
  cout << "Entrez les coordonnees du point X (en metres) :" << endl;
  cout << "(le point X peut etre situe a l'interieur de l'ellipsoide)"
       << endl;
  cout << "x = ";   cin >> xp[0];
  cout << "y = ";   cin >> xp[1];
  cout << "z = ";   cin >> xp[2];
  
  int div;
  cout << "Entrez le nombre de divisions de l'axe majeur (a) :" << endl;
  cout << "N = ";   cin >> div;
  double dx[3];
  for (int i=0; i<3; i++)
    dx[i] = a[i]/div;
    
  //boucle sur les "sous-volumes" de l'ellipsoide :
  //si le centre d'un sous-volume est dans l'ellipsoide
  //on additionne sa contribution au champ total
  double xsub[3];
  int sousVolumes = 0;
  double f[3];
  double field[3] = {0.};
  bool inside;
  //on commence par le cote negatif des axes
  for (int ix=0; ix<(2*div); ix++) {   //axe x
    for (int iy=0; iy<(2*div); iy++) {   //axe y
      for (int iz=0; iz<(2*div); iz++) {   //axe z
        //calcule de la position du suos-volume en 3 dim
        for (int i=0; i<3; i++)
          xsub[i] = -a[i] + dx[i]/2. + ix*dx[i];
          
        //on verifie que le point est dans le volume considere                        
        inside = verifierVolume(xsub,a);
        if (inside) {
          sousVolumes += 1;
          //contribution de ce sousvolume au champ gravitationnel :
          //on considere GN = 1 et m = 1: F = G*m/(r*r) = 1/(r*r)                         
          //il faut maintenant determiner les composantes du vecteur
          //pour pouvoir sommer toutes les contributions
          fieldComp(xsub,xp,f);
          for (int i=0; i<3; i++)
            field[i] += f[i];
        }
      }
    }
  }
  
  cout << sousVolumes << " sous-volumes sont definis." << endl;
  //on normalise le champ a la masse de l'objet ;
  //M / # sous-volumes = masse sous-volume
  double ff = 0.;
  for (int i=0; i<3; i++) {
    field[i] *= -GN*M/sousVolumes;
    ff += pow(field[i],2);
  }
  ff = sqrt(ff);
    
  cout << "Estimation du champ : |F| = " << ff << " m/s2 ." << endl;
  cout << " Fx = " << field[0] << endl;
  cout << " Fy = " << field[1] << endl;
  cout << " Fz = " << field[2] << endl;

  return 0;
}

bool verifierVolume(double x[], double a[]) {
  bool inside = false;

  //equation de l'ellipsoide  
  if ((pow(x[0]/a[0],2) + pow(x[1]/a[1],2) + pow(x[2]/a[2],2)) <= 1.)
  inside = true;
    
  return inside;
}

//calcul du champ gravitationnel en X
void fieldComp(double x[], double p[], double f[]) {
  //calcul de la distance entre le centre du sous-volume
  //et le point X
  double dist = 0.;
  for (int i=0; i<3; i++)
    dist += pow(p[i]-x[i],2);
  dist = sqrt(dist);
     
  //composantes de G non normalisees
  for (int i=0; i<3; i++)
    f[i] = (p[i]-x[i]) / pow(dist,3);

  return;
}