/*
 * MONTE CARLO EVALUATION OF INTEGRALS
 *
 * Peter Signell, 3/19/2002
 */
#include"lib351.h"

void main(){

    // output initialization
	ofstream dataOutStr;
    outFileMaker("m355p1c.out", &dataOutStr);
    
    // program variable definitions
 	int intPoints;
	int intRuns;
	double dblMass;
	double dblR1;
	double dblR2;
	double dblDevAvMass;
	double dblRanNum;
	double dblAvMass;
	double dblAvMass2;
	double dblPoints;
	double dblRuns;
    double dblX[4];
	double dblCmX[4];
	double dblDevAvCmX[4];
	double dblAvCmX[4];
	double dblAvCmX2[4];
    int intRun;
	int k;
	int j;
    const double dblSeed = 13.;
    const double dblConst = 899.;
    const double dblMod = 32768.;
    const double dblMaxNum = dblMod - 1.;

	// cycle through the runs
  	bool goout = false;
  	while (!goout){

    	// data input
    	cout<<"Please, supply values :"<<endl;
    	cout<<"\nnumber of randomly located points to use, NUM ="<<endl;
    	cin>>intPoints;
    	cout<<"\nnumber of runs to make using this value of NUM, NRUNS ="<<endl;
    	cin>>intRuns;

		// cast some variables
		dblPoints = double(intPoints);
		dblRuns = double(intRuns);

		// print the header
		dataOutStr<<"\f"<<endl;
		dataOutStr<<"Monte Carlo Evaluation of Integrals"<<endl<<endl;
        dataOutStr<<"number of randomly located points to use, NUM = "<<intPoints<<endl;
		dataOutStr<<"number of runs to make using that number of points, NRUNS = "<<intRuns<<endl<<endl;

		// initialize the density for "inside"
    	const double dblRho = 1.;
    	
    	// initialize the random number generator
    	dblRanNum = dblSeed;
    	
    	// initialize the "run-averaged" variables
		dblAvMass = 0.;
		dblAvMass2 = 0.;
		for (k=1; k<=3;k++) {
			dblAvCmX[k] = 0.;
			dblAvCmX2[k] = 0.;
		}
		
		// loop over runs with different sets of random numbers
    	for (intRun=1; intRun <= intRuns; intRun++) {
    		//initialize the mass and cm coordinates
			dblMass = 0.;
    		for (k=1; k<=3; k++) {
    			dblCmX[k] = 0.;
    		}
    		// loop over the random points in this run
			for (j=1; j<=intPoints; j++) {
    			// generate 3 more random numbers in (-0.5,+0.5)
				// call the first random number the x-coordinate, the second y, etc.
    			for (k=1; k<=3; k++) {
    				dblRanNum = fmod(dblConst*dblRanNum,dblMod);
					dblX[k] = dblRanNum/dblMaxNum - 0.5;
				}

				// if the coordinates are inside the object, add to mass and cm coordinates
				dblR1 = pow(dblX[1],2) + pow(dblX[2],2) + pow(dblX[3],2);
				dblR2 = pow(dblX[1],2) + pow(dblX[2],2);
				if ((dblR1 <= 0.25) && (dblR2 >= 0.09)) {
					dblMass += dblRho/dblPoints;
					for (k=1; k<=3; k++) {
						dblCmX[k] += dblX[k]*dblRho/dblPoints;
					}
				}
			}
			dblAvMass += dblMass/dblRuns;
			dblAvMass2 += pow(dblMass,2)/dblRuns;
			for (k=1; k<=3; k++) {
				dblCmX[k] = dblCmX[k]/dblMass;
				dblAvCmX[k] += dblCmX[k]/dblRuns;
				dblAvCmX2[k] += pow(dblCmX[k],2)/dblRuns;
			}
			cout<<"Run "<<intRun<<": mass = "<<dblMass<<"; x-cm = "<<dblCmX[1]<<"; y-cm = "<<dblCmX[2]<<"; z-cm = "<<dblCmX[3]<<endl<<endl;
			dataOutStr<<"Run "<<intRun<<": mass = "<<dblMass<<"; x-cm = "<<dblCmX[1]<<"; y-cm = "<<dblCmX[2]<<"; z-cm = "<<dblCmX[3]<<endl<<endl;
		}
		dblDevAvMass = sqrt(dblAvMass2 - pow(dblAvMass,2));
		for (k=1; k<=3; k++) {
			dblDevAvCmX[k] = sqrt(dblAvCmX2[k] - pow(dblAvCmX[k],2));
		}
		cout<<"av mass = "<<dblAvMass<<"; av x-cm = "<<dblAvCmX[1]<<"; av y-cm = "<<dblAvCmX[2]<<"; av z-cm = "<<dblAvCmX[3]<<endl;
		dataOutStr<<"av mass = "<<dblAvMass<<"; av x-cm = "<<dblAvCmX[1]<<"; av y-cm = "<<dblAvCmX[2]<<"; av z-cm = "<<dblAvCmX[3]<<endl;
		cout<<"dev in av mass = "<<dblDevAvMass<<"; dev in av x-cm = "<<dblDevAvCmX[1]<<"; dev in av y-cm = "<<dblDevAvCmX[2]<<"; dev in av z-cm = "<<dblDevAvCmX[3]<<endl;
		dataOutStr<<"dev in av mass = "<<dblDevAvMass<<"; dev in av x-cm = "<<dblDevAvCmX[1]<<"; dev in av y-cm = "<<dblDevAvCmX[2]<<"; dev in av z-cm = "<<dblDevAvCmX[3]<<endl;

		// another run?
		goout = query();

	}
    // close the output file
	fileClose("m355p1c.out", &dataOutStr);
};