#include "nr3.h"
#include "ran.h"

/* Monte-carlo test. integrate the mass and center-of-mass coordinates
of a hollow sphere, (inner radius rmin = 0.98 outer radius rmax= 1.0) of
Aluminum (rho 2.7), partly filled with water. The water sinks to the bottom,
the max water level in this example is z=-0.6.

The mass and C.o.m of the sphere is analytically given by
4/3pi (rmax^3-rmin^3) rho and (x,y,z)_com = 0,
and of the water level by 
mass = 2pi/3 rmin^3 (1-zmax)= 2pi/3 rmin^3 0.4, and 
x_com=y_com=0, z_com = rmin * (3/4 x (1+zmin)/2 x (1-zmin^4)

draw from spherical coordinates and from cartesian coordinates to show
the speed of the M.C. and to verify that the weighing is correct.

show advantage drawing with dk=dr^3 and dcostheta instead of dtheta -
in the latter case many variables are wasted.

this is an example only, you would not do an analytical integral in a MC
*/
int main() {
	const double PI = acos(-1);
	double rmin=0.98;
	double rmax=1.0;
	double z_level=-0.6;
	double theomassshell=(pow(rmax,3)-pow(rmin,3))*4*PI/3.0*2.7; // mass shell
	double theomassliq=pow(rmin,3)*2*PI/3.0+PI*z_level*pow(rmin,2)-PI/3.0*pow(z_level,3);
	double theozliq=PI/4*(2*pow(rmin*z_level,2)-pow(z_level,4)-pow(rmin,4));
	theozliq/=theomassliq;
	double theoz=theozliq*(theomassliq/(theomassliq+theomassshell));
	Ran ran(1234567); // random deviate
// integration over the volume, 1000000 randomly drawn results, for spherical coordinates
	double V1=4.0*PI*pow(rmax,3)/3.0; // total detection volume
	double mass=0;
        double xcm=0;
        double ycm=0;
	double zcm=0;
	for (int i=0; i<1000000; i++) {
		double costhet=(1-2*ran.doub()); // random number between 1 and -1
		double thet=acos(costhet);
		double y=ran.doub()*pow(rmax,3); // dy =dr^3
		double r=pow(y,1.0/3.0);
		double phi=2*PI*ran.doub(); // random phi
		if (r>rmin) {  // point inside the aluminum shell 
			mass+=2.7;
			xcm+=2.7*cos(phi)*sin(thet)*r;
			ycm+=2.7*sin(phi)*sin(thet)*r;
			zcm+=2.7*costhet*r;
		} else if (r*costhet<z_level) { // point inside the water
			mass+=1;
			xcm+=cos(phi)*sin(thet)*r;
			ycm+=sin(phi)*sin(thet)*r;
			zcm+=costhet*r;
		}
	}
	cout << "Monte-carlo results drawn spherically, from dcostheta and from dr^3 " << endl;
	cout << "mass : " << V1*mass/1000000 << " theoretical " << theomassliq+theomassshell << endl; 
	cout << " xcm : " << xcm/mass << endl;
	cout << " ycm : " << ycm/mass << endl;
	cout << " zcm : " << zcm/mass << " theoretical " << theoz << endl;
	cout << endl;
// try Cartesian coordinates, same number of counts
	V1=pow(2*rmax,3); // total detection volume, draw from -rmax to +rmax
	mass=0;
        xcm=0;
        ycm=0;
	zcm=0;
	for (int i=0; i<1000000; i++) {
		double x=(1-2*ran.doub())*rmax; // random number between -rmax and +rmax
		double y=(1-2*ran.doub())*rmax;
		double z=(1-2*ran.doub())*rmax;
		double r=sqrt(x*x+y*y+z*z);
		if (r>rmax) continue; // outside sphere
		if (r>rmin) {  // point inside the aluminum shell 
			mass+=2.7;
			xcm+=2.7*x;
			ycm+=2.7*y;
			zcm+=2.7*z;
		} else if (z<z_level) { // point inside the water
			mass+=1;
			xcm+=x;
			ycm+=y;
			zcm+=z;
		}
	}
	cout << "Monte-carlo results drawn Cartesian : " << endl;
	cout << "mass : " << V1*mass/1000000 << " theoretical " << theomassliq+theomassshell << endl; 
	cout << " xcm : " << xcm/mass << endl;
	cout << " ycm : " << ycm/mass << endl;
	cout << " zcm : " << zcm/mass << " theoretical " << theoz << endl;
	cout << endl;
	// now draw r, theta
	V1=2.0*PI*PI*rmax; // total detection volume
	mass=0;
        xcm=0;
        ycm=0;
	zcm=0;
	for (int i=0; i<1000000; i++) {
		double thet=PI*ran.doub(); // theta
		double r=ran.doub()*rmax; // r
		double phi=2*PI*ran.doub(); // random phi
		double w=r*r*sin(thet);
		if (r>rmin) {  // point inside the aluminum shell 
			mass+=2.7*w;
			xcm+=2.7*cos(phi)*sin(thet)*r*w;
			ycm+=2.7*sin(phi)*sin(thet)*r*w;
			zcm+=2.7*cos(thet)*r*w;
		} else if (r*cos(thet)<z_level) { // point inside the water
			mass+=w;
			xcm+=cos(phi)*sin(thet)*r*w;
			ycm+=sin(phi)*sin(thet)*r*w;
			zcm+=cos(thet)*r*w;
		}
	}
	cout << endl;
	cout << "Monte-carlo results drawn spherically, from r and theta " << endl;
	cout << "mass : " << V1*mass/1000000 << " theoretical " << theomassliq+theomassshell << endl; 
	cout << " xcm : " << xcm/mass << endl;
	cout << " ycm : " << ycm/mass << endl;
	cout << " zcm : " << zcm/mass << " theoretical " << theoz << endl;
	cout << endl;
}