/*
*                   SIMULATION OF A CHAIN REACTION
*/

#include<iostream.h>
#include<math.h>
#include"lib351.h"

void main(){

    // output initialization
	ofstream dataOutStr;
    outFileMaker("m356p1c.out", &dataOutStr);
    
    // program variable definitions
    int numNeutrons, insideReactions, nRun;
    double fissionNumbers[10];
    double mass, shape;
    double sizeA, sizeB, x0, y0, z0, x1, y1, z1, distance;
    double phi, costh, sinth, survivalFraction;
	double dblRanNum;
    const double randInit = 13.;
	const double m_pi = 3.141592654;
    const double dblSeed = 13.;
    const double dblConst = 899.;
    const double dblMod = 32768.;
    const double dblMaxNum = dblMod - 1.;

	// initialize the random number generator
    dblRanNum = dblSeed;

	// print the header
	dataOutStr<<"FISSION SURVIVAL FRACTION CALCULATION";
 
		// cycle through the runs
	nRun = 0;
	bool goout = false;
	while (!goout){

    	nRun++;
		
		// data input
  		cout <<"Please supply values for numNeutrons, shape and mass:"<< endl;
  		cout <<"\nnumber of randomly generated neutrons, numNeutrons ="<< endl;
        cin  >> numNeutrons;
        cout <<"\nshape of the slab, shape ="<< endl;
        cin  >> shape;
  		cout <<"\nmass of the rectangular slab, mass ="<< endl;
        cin  >> mass;

		// print the header
		dataOutStr<<endl<<endl;
		dataOutStr<<"run "<<nRun<<": "<<"numNeutrons= "<<numNeutrons<<"; shape="<<shape<<"; mass="<<mass;

        // do the calculation
        sizeA = pow((mass*shape),(1./3));
        sizeB = pow(mass/(shape*shape),(1./3));

        insideReactions = 0;
 
        for (int i = 1; i <= numNeutrons; i++){
        	for (int j = 1; j <= 9; j++){
    			dblRanNum = fmod(dblConst*dblRanNum,dblMod);
	    		fissionNumbers[j] = dblRanNum/dblMaxNum;
            }

            x0 = sizeA * (fissionNumbers[1] - 0.5);
            y0 = sizeA * (fissionNumbers[2] - 0.5);
  		    z0 = sizeB * (fissionNumbers[3] - 0.5);

            for (int k = 1; k <= 2; k++){
            	phi = 2 * m_pi *fissionNumbers[2*k+2];
            	costh=2*(fissionNumbers[2*k+3]-0.5);
            	sinth = sqrt(1 - pow(costh,2));

            	distance = fissionNumbers[k+7];
           		x1 = x0+distance * sinth * cos(phi);
            	y1 = y0+distance * sinth * sin(phi);
            	z1 = z0 + distance * costh;

        		if (fabs(x1)<=sizeA/2 && fabs(y1)<=sizeA/2&&fabs(z1)<= sizeB/2 ){
        			insideReactions += 1;
                }
        	}
        }
		survivalFraction = double(insideReactions)/double(numNeutrons);

		// data out
		dataOutStr <<"; survival fraction = "<<survivalFraction<<endl;

		// another run?
		goout = query();
	}
    // close the output file
	fileClose("m356p1c.out", &dataOutStr);
};