#include <vector>
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <math.h>
#include <random>
#include <chrono>
#include <algorithm>
#include <cfloat>
#include <limits>
#include <fstream>
#include <cstdio>

#define LINES 1000000ul

using namespace std;

double eta07 = (2.639/26.)*6.022e23*1.; //Li or Li-7 per cm^3
double eta19 = eta07; //F-19 per cm^3
double eta06 = 0.;//0.0759 * eta07; //Li-6 per cm^3
//double sig07 = 1.; //barn
//double sig19 = 4.; //barns
unsigned long neutrons = (unsigned long)1e6;
double initKE = 17e6; //eV
double deadKE = 1e-2; //eV (thermal)
vector<double> sigmaT07, sigmaT19, sigmaT06, sigmaC07, sigmaC19, sigmaC06;

double crossSectionBarn ( signed short massNum, double kineticEnergy, bool load = false );

int main ( int argc, char **argv ) {
  
  if ( eta06 > 0. )
    eta07 -= eta06;
  if ( crossSectionBarn ( 7, 17e6, true ) != (1.40383-1.74425e-6) ) return 1;
  
  default_random_engine generator;
  uniform_real_distribution<double> distribution(0.0,1.0);
  
  double x0 = 0., y0 = 0., z0 = 0., rad3d = 0.;
  remove("LiF.txt");
  ofstream outFile;
  outFile.open("LiF.txt",ios::app);
  
  for ( unsigned long int i = 0; i < neutrons; ++i ) {
    
    double KE = initKE;
    double xi = x0, yi = y0, zi = z0;
    
    while ( KE > deadKE ) {
      
      double mfpS07 = 1. / ( eta07 * 1e-24 * crossSectionBarn ( 7, KE ) );
      double mfpS19 = 1. / ( eta19 * 1e-24 * crossSectionBarn ( 19, KE ) );
      double mfpA07 = 1. / ( eta07 * 1e-24 * crossSectionBarn (-7, KE ) );
      double mfpA19 = 1. / ( eta19 * 1e-24 * crossSectionBarn (-19, KE ) );
      double mfpS06 = 1. / ( eta06 * 1e-24 * crossSectionBarn ( 6, KE ) );
      double mfpA06 = 1. / ( eta06 * 1e-24 * crossSectionBarn (-6, KE ) );
      //cerr << mfpA07 << "\t" << mfpA19 << endl; return 1;
      
      double fpS07 = -mfpS07 * log(distribution(generator));
      double fpS19 = -mfpS19 * log(distribution(generator));
      double fpA07 = -mfpA07 * log(distribution(generator));
      double fpA19 = -mfpA19 * log(distribution(generator));
      double fpS06 = -mfpS06 * log(distribution(generator));
      double fpA06 = -mfpA06 * log(distribution(generator));
      vector<double> fp;
      fp.push_back(fpS07);
      fp.push_back(fpS19);
      fp.push_back(fpA07);
      fp.push_back(fpA19);
      fp.push_back(fpS06);
      fp.push_back(fpA06);
      std::sort(fp.begin(),fp.end());
      if ( fabs(fp[0]-fpA07) < 1e-6 || fabs(fp[0]-fpA19) < 1e-6 || fabs(fp[0]-fpA06) < 1e-6 )
	KE = 0.;
      
      double phi, theta;
      //if ( (xi == x0 && yi == y0 && zi == z0) || KE == initKE ) {
	phi = 2. * M_PI * distribution(generator);
	theta = acos ( 2. * distribution(generator) - 1. );
      //}
      
      double dx = fp[0] * sin(theta) * cos(phi);
      double dy = fp[0] * sin(theta) * sin(phi);
      double dz = fp[0] * cos(theta);
      double xf = xi + dx;
      double yf = yi + dy;
      double zf = zi + dz;
      xi = xf;
      yi = yf;
      zi = zf;
      if ( KE <= 0. )
	KE = 0.;
      else
	KE *= distribution(generator);
      rad3d = sqrt ( (xf-x0)*(xf-x0)+(yf-y0)*(yf-y0)+(zf-z0)*(zf-z0) );
      fp.clear();
      
    }
    
    //cout << rad3d << endl;
    outFile << rad3d << endl;
    
  }
  
  outFile.close();
  return 0;
  
}

double crossSectionBarn ( signed short int massNum, double kineticEnergy, bool load ) {

  if ( load ) {
    
    const char* xSectFileName = "Li7F19Li6.txt";
    FILE *xSectFile = fopen(xSectFileName,"r");
    double temp0 = 0., temp1 = 0., temp2 = 0., temp3 = 0., temp4 = 0., temp5 = 0.;
    
    for ( unsigned long int j = 0; j < LINES; ++j ) {
      fscanf(xSectFile, "%lf %lf %lf %lf %lf %lf", &temp0, &temp1, &temp2, &temp3, &temp4, &temp5 );
      sigmaT07.push_back(temp0);
      sigmaT19.push_back(temp1);
      sigmaT06.push_back(temp2);
      sigmaC07.push_back(temp3);
      sigmaC19.push_back(temp4);
      sigmaC06.push_back(temp5);
      //fprintf(stderr,"%f\t%f\t%f\t%f\t%f\t%f\n",temp0,temp1,temp2,temp3,temp4,temp5);
    }
    fclose(xSectFile);
    
  }
  
  long index = long(floor((log10(kineticEnergy)+2.)/9.3e-6+0.5));
  if ( index < 0 ) index = 0;
  if ( index >= LINES ) index = LINES - 1;
  
  double sigma = 0.0;
  switch ( massNum ) {
  case 7:
    sigma = sigmaT07[index] - sigmaC07[index];
    break;
  case 19:
    sigma = sigmaT19[index] - sigmaC19[index];
    break;
  case 6:
    sigma = sigmaT06[index] - sigmaC06[index];
    break;
  case -7:
    sigma = sigmaC07[index];
    break;
  case -19:
    sigma = sigmaC19[index];
    break;
  case -6:
    sigma = sigmaC06[index];
    break;
  default:
    sigma = 1e6;
  }
  
  if ( sigma <= 0. ) sigma = DBL_MIN;
  if ( sigma >= DBL_MAX ) sigma = 1.0;
  return sigma;
  
}