/*
   Calculate inertia tensor of geometric object via random sampling.
   
   NB requires C++ V11 standard.
   
   http://www.cplusplus.com/reference/random/uniform_real_distribution/
*/
 
#include <iostream>
#include <random>

using namespace std;

// Function to decide whether point is inside object or not.
bool isInside(double x,double y,double z)
{
   // Cyllinder of radius 1 in the yz plane
   return (y*y+z*z)<1.0;

/*
   // Sphere with radius 1
   return (x*x+y*y+z*z)<1.0;
*/
}

int main()
{
// Number of points to sample
   const int N=10000009;
   mt19937 re;                     // using a Mersenne Twister

// define bounding box (BB) of object.
// Hence cyllindrical rod of length 20, and radius 1
   double xlo=-10,xhi=10;   
   double ylo=-1 ,yhi=1;   
   double zlo=-1 ,zhi=1;   
   double BBvol=(xhi-xlo)*(yhi-ylo)*(zhi-zlo);

// Pick random number uniform distributed in BB
   uniform_real_distribution<double> distx(xlo,xhi);
   uniform_real_distribution<double> disty(ylo,yhi);
   uniform_real_distribution<double> distz(zlo,zhi);

// Points inside object
   int count=0;

// Find centre-of-mass
   double Rx,Ry,Rz;
   Rx=Ry=Rz=0.0;

   for (int i=0 ; i<N ; i++)
     {
        double x=distx(re);
        double y=disty(re);
        double z=distz(re);
        
        if (isInside(x,y,z))
           {
               Rx+=x;
               Ry+=y;
               Rz+=z;
               count++;
           }
     }

   Rx/=count; Ry/=count; Rz/=count;
   cout << "Rcm=" << Rx << " " << Ry << " " << Rz << "\n";  
   cout << "Vol=" << BBvol*count/N << "\n";  

// Moments of inertia
   count=0;
   double Mxx,Myy,Mzz,Mxy,Mxz,Myz;
   Mxx=Myy=Mzz=Mxy=Mxz=Myz=0.0;

   for (int i=0 ; i<N ; i++)
     {
        double x=distx(re);
        double y=disty(re);
        double z=distz(re);
        
        if (isInside(x,y,z))
           {
               Mxx+=(x-Rx)*(x-Rx);
               Myy+=(y-Ry)*(y-Ry);
               Mzz+=(z-Rz)*(z-Rz);
               Mxy+=(x-Rx)*(y-Ry);
               Mxz+=(x-Rx)*(z-Rz);
               Myz+=(y-Ry)*(z-Rz);
               count++;
           }
     }
     
   Mxx/=count; Myy/=count; Mzz/=count; 
   Mxy/=count; Mxz/=count; Myz/=count; 
 
   cout << "Mxx=" << Mxx << " Myy=" << Myy << " Mzz=" << Mzz << "\n";  
   cout << "Mxy=" << Mxy << " Mxz=" << Mxz << " Myz=" << Myz << "\n";  
     
   
   return 0;
}