Scripts for sampling fusion products (neutrons and alpha particles). The (R,Z) position is sampled 
according to the fusion reaction rate, the toroidal angle is sampled uniformly, the velocity is isotropic
and the energy is sampled in a manner consistent with the thermal Doppler broadening. 

Can either use prescribed Miller geometry flux surfaces or read a GEQDSK
numerical equilbrium. Uses prescribed 'Topeol style' kinetic profiles.

The method for sampling is described in L. Reali et al: <Cite our paper>.

The Miller geometry and profiles are described in Fausser et al, Fusion Engineering and Design
Volume 87, Issues 5–6, August 2012, Pages 787-792 https://doi.org/10.1016/j.fusengdes.2012.02.025

The Doppler Broadening calculation is derived in H Brysk 1973 Plasma Physics 15 611
DOI 10.1088/0032-1028/15/7/001

To sample particles edit the params dictionary in setup.py to define the profile parameters, geometry parameters and run settings.
Make sure the directory pointed to by the entry "output_dir" exists.

Then run:
> python sampling.py

This will generate nsamples of the chosen species and store them in a pickle file in the output directory. Reading the pickle file
will return a dictionary containing the following entries:

data = { 'R'   : R coordinates of sampled particles [m]
       	 'Phi' : Toroidal angles of sampled particles
	 'Z'   : Z coordinates of sampled particles [m]
	 'vx'  : x coordinate of velocity unit vector
	 'vy'  : y coordinate of velocity unit vector
	 'vz'  : z coordinate of velocity unit vector
	 'E'   : Energy of particle [MeV]
	 }

The output directory will also contain some diagnostic plots.

One use-case of these scripts is to generate neutrons with OpenMC. To do this set write_cdfs = True and run the code as normal.
Instead of generating samples, the code will calculate the cumulative distribution functions (CDFs) needed to perform sampling and
store them in a JSON file in the output directory named based on the cdf_filename setting described below. Using the json.hpp, 
plasma_source.hpp and plasma_source.cpp files in the cpp directory, OpenMC can read this JSON file and use the CDFs within to
sample neutrons. 

Possible future developments:
Parallelisation of CDF generation and sampling.
Expand OpenMC coupling to handle numerical equilibria and profiles (currently only works with Miller parameterisations and analytic profiles).

-------------------------------------------------------------------

Input settings required in params dictionary:

WARNING: RHO CAN MEAN DIFFERENT THINGS:
For Miller equilibria:    rho = normalised minor radius = a / a_sep
For numerical equilibria: rho = normalised poloidal flux 

# Read numerical profiles from a JSON file. 
# The JSON file is expected to contain arrays named 'rho', 'ni' 
# and 'Ti' for the flux coordinate, ion density [m**-3] and ion
# temperature [keV] respectively. These will then be interpolated.

read_profiles    : Set True to read profile information from JSON file.
profile_filename : The name of the JSON file to read.

# Profile parameters - see Fausser et al.
rho_ped : The pedestal rho 

# Density Profile Settings [m**-3]
ni_ped  : The pedestal density
ni_sep  : The separatrix density
ni_0    : The core density
alpha_n, beta_n : The density peaking factors

# Temperature Profile [keV]
Ti_ped  : The pedestal temperature
Ti_sep  : The separatrix temperature
Ti_0    : The core temperature
alpha_t, beta_t  : The temperature peaking factors

# Flag to set whether using a numerical equilibrium
use_geqdsk : True to read equilibrium from GEQDSK file
geqdsk_filename : Name of GEQDSK file to read

# Numerical Geometry parameters [ used if use_geqdsk = True ]
sample_sqrt_psi : Controls spacing of psi contours when constructing coordinate charts.
		  If True points are spaced evenly in sqrt(psi)
		  If False points are spaced evenly in psi.
radial_steps    : Number of radial steps to use when constructing R(psi) map for outer midplane.

rho_low         : Below rho_low a step size of step_low is used when tracing contours.
step_low          Above this value step_high is used. It is common to make step_high >> step_low
step_high         to accelerate contour tracking for larger flux surfaces.

nalpha_downsample : Downsample the number of points on the contours to be < 2 * nalpha_downsample.
		  This can be necessary due to the difficulty of interpolating s(theta) for contours
		  with very large numbers of points. 

npsi_map        : The number of flux surfaces to track when constructing the coordinate charts.
nalpha_map      : The number of alpha points to use when constructing the coordinate charts.

# Miller Geometry Parameters - Also described in Fausser et al. [ Used if use_geqdsk = False ]
R0      : The plasma geometric major radius [m]
delta   : The plasma triangularity
kappa   : The plasma elongation
shift   : The Shafranov Shift [m]
a_sep   : The plasma minor radius [m]

# Run Settings

# Grid Reolution
nrho    : The number of radial bins used to construct the radial sampling CDF
nalpha  : The number of poloidal angle bins used to construct the angular sampling CDFs (one per radial bin)

nsamples   : Number of samples to generate
seed       : The random number seed
bosch-hale : The form of the reactivity to use: 0 = Sadler-Van-Belle, 1 = Bosch-Hale
cdf_interp_kind : The interpolation type for CDFs - Cubic seems to throw bounds errors so best to leave as linear.
species    : Sets the species to sample : 0 = neutron, 1 = alpha

# Output settings

output_dir      : Output directory location to store files
write_cdfs      : Set True to write cdf information to a JSON file for use in OpenMC. Code exits after writing.  
cdf_filename    : This sets the filename for the JSON file containing the CDF information.
pickle_filename : The name of the pickle file to dump the sampled data to

show_plots      : Set True to view diagnostic plots during run. If False plots are just written to output directory.

# DEBUGGING ENTIES (development only) ----------------
DEBUG_JACOBIAN  
DEBUG_CDFS     
DEBUG_GEOMETRY 
DEBUG_REACTIVITY
DEBUG_PROFILES