estar.h File Reference

Core ESTAR density correction calculation interface. More...

#include <string>
#include "estar_formulaCalcs.h"
#include "estar_dataParser.h"

Go to the source code of this file.

Classes
struct	bspol
	Result struct returned by fbspol(). More...
struct	scof
	Spline coefficient struct returned by fscof(). More...

Functions
void	outputDensityFile (float mediaDensity, double *densityCorr, double *enGrid, float *meanIval, formula_calc fc, std::string outputFilename)
	Output a density correction file. More...
int	estarCalculation (int isCompound, int NEP, float mediaDensity, std::string *elementArray, double *massFraction, float *numOfAtoms, double *densityCorr, double *enGrid, float *meanIval, float *ipotval, int mediaNum, std::string outputFilename)
	Compute density correction factors for a medium. More...
bspol	fbspol (double s, const std::vector< double > &x, const std::vector< double > &a, const std::vector< double > &b, const std::vector< double > &c, const std::vector< double > &d)
	Evaluate the cubic spline interpolant at a given point. More...
scof	fscof (int nmax, const vector< double > &x, const vector< double > &f)
	Compute natural cubic spline coefficients. More...
double	objective_function (double tau, const vector< double > &f, const vector< double > &eps, int nmax, double x)
	Objective function for the bisection solver. More...
double	bisec (double lowerbound, double upperbound, double tolerance, double tau, const std::vector< double > &f, const std::vector< double > &eps, int nmax)
	Find the root of objective_function() by bisection. More...

Detailed Description

Core ESTAR density correction calculation interface.

This header exposes the public interface for computing electron stopping power density corrections using the ESTAR method described in:

Sternheimer, R.M., Berger, M.J., and Seltzer, S.M. (1984). "Density effect for the ionization loss of charged particles in various substances." Atomic Data and Nuclear Data Tables, 30, 261-271.

The main entry points for external callers are estar_() and compoundstoelements_(), which are called directly from Mortran/Fortran via the EGSnrc C/Fortran interface. All other functions in this header are internal implementation details exposed for unit testing purposes.

Definition in file estar.h.

Function Documentation

outputDensityFile()

void outputDensityFile	(	float	mediaDensity,
		double *	densityCorr,
		double *	enGrid,
		float *	meanIval,
		formula_calc	fc,
		std::string	outputFilename
	)

Output a density correction file.

Writes the density correction factors, energy grid, mean I-value, atomic numbers and mass fractions to a file in the standard EGSnrc density correction format (113-point energy grid).

Parameters

mediaDensity	Density of the medium (g/cm^3).
densityCorr	Array of 113 density correction factors delta, one per point on the standard energy grid.
enGrid	Array of 113 energies (MeV) defining the standard energy grid.
meanIval	Pointer to the mean excitation energy I (eV) of the medium.
fc	formula_calc struct containing the atomic numbers, weight fractions and element count for the medium.
outputFilename	Path to the output density correction file. Trailing whitespace is stripped before use.

estarCalculation()

int estarCalculation	(	int	isCompound,
		int	NEP,
		float	mediaDensity,
		std::string *	elementArray,
		double *	massFraction,
		float *	numOfAtoms,
		double *	densityCorr,
		double *	enGrid,
		float *	meanIval,
		float *	ipotval,
		int	mediaNum,
		std::string	outputFilename
	)

Compute density correction factors for a medium.

Implements the Sternheimer (1984) method for computing the density effect correction delta to the Bethe-Bloch stopping power formula. The correction factors are computed on the standard 113-point energy grid and stored in densityCorr. The corresponding energies are stored in enGrid.

The density-effect cutoff (the energy below which delta=0) was a feature of Sternheimer's parametric approximation (Sternheimer & Peierls, Phys. Rev. B3, 3681, 1971), introduced to reduce computational burden. Since EGSnrc uses the "exact" Sternheimer method (Sternheimer, Berger & Seltzer, Atom. Data Nucl. Data Tabl. 30, 261, 1984), in which delta emerges naturally from the oscillator- strength calculation and requires no imposed cutoff, this function is no longer needed.

Internally this function:

1. Calls getDataFromFormulae() to compute the mean excitation energy and mean Z/A for the medium.
2. Builds the dispersion oscillator model (Sternheimer eq. 5).
3. Solves for the adjustment factor using Newton's method (eq. 8).
4. Evaluates the density correction on the energy grid using a bisection solver and cubic spline interpolation (eq. 1).
5. Writes the results to a density correction file via outputDensityFile().

Parameters

isCompound	1 if the medium is a compound, 0 otherwise.
NEP	Number of elements (or compounds) in elementArray.
mediaDensity	Density of the medium (g/cm^3). Must be > 0.
elementArray	Array of NEP element chemical symbols (e.g. "H", "O").
massFraction	Array of NEP mass fractions (weight fractions), one per element. Must sum to 1.
numOfAtoms	Array of NEP atom counts per formula unit (used for compounds).
densityCorr	Output array of 113 density correction factors. Must be pre-allocated by the caller.
enGrid	Output array of 113 energies (MeV). Must be pre-allocated by the caller.
meanIval	On input: if >= 0, overrides the computed I-value with this value (eV). On output: holds the I-value actually used.
ipotval	Pointer to the user-supplied I-value override (eV). Set to -1 if no override is desired.
mediaNum	Integer ID of the medium (used for diagnostic output).
outputFilename	Path to write the density correction file. Pass an empty string to suppress file output.

Returns

0 on success, 9 on any error (bad input, non-convergence, etc.). Error details are printed to stdout before returning.

fbspol()

bspol fbspol	(	double	s,
		const std::vector< double > &	x,
		const std::vector< double > &	a,
		const std::vector< double > &	b,
		const std::vector< double > &	c,
		const std::vector< double > &	d
	)

Evaluate the cubic spline interpolant at a given point.

Given the spline coefficients produced by fscof() and a query point s, performs a binary search to locate the bracketing interval and evaluates the cubic polynomial within that interval using Horner's method.

The result also carries the bracketing interval indices, which are used by the caller to set the bounds for the bisection solver (bisec()).

Parameters

s	The query point (log of y = tau*(tau+2), where tau = E/m_e c^2).
x	Array of n knot positions (log-energy values yql[]).
a	Spline coefficient array a[] from fscof().
b	Spline coefficient array b[] from fscof().
c	Spline coefficient array c[] from fscof().
d	Spline coefficient array d[] from fscof().
n	Number of knots.

Returns

A bspol struct containing the interpolated density correction and the bracketing interval indices.

Definition at line 980 of file estar.cpp.

References bspol::density_corr, egsFatal, egsWarning, bspol::lb_index, and bspol::ub_index.

fscof()

scof fscof	(	int	nmax,
		const vector< double > &	x,
		const vector< double > &	f
	)

Compute natural cubic spline coefficients.

Fits a natural cubic spline (zero second derivative at both endpoints) to the data points (x[i], f[i]) for i = 0, ..., nmax-1. The resulting coefficients are stored in the returned scof struct and are intended for use with fbspol().

Parameters

nmax	Number of data points. Must be >= 2.
x	Array of nmax knot positions (log-energy values yql[]). Must be strictly monotone.
f	Array of nmax function values (density corrections d[]).

Returns

A scof struct containing the four coefficient arrays a, b, c, d.

Definition at line 1062 of file estar.cpp.

References scof::a, scof::b, scof::c, scof::d, and egsFatal.

objective_function()

double objective_function	(	double	tau,
		const vector< double > &	f,
		const vector< double > &	eps,
		int	nmax,
		double	x
	)

Objective function for the bisection solver.

Evaluates the function whose root locates the exact density correction for a given kinetic energy. The root condition is:

log(y) - log(1 / sum_i( f[i] / (eps[i] + x) )) = 0

where y = tau*(tau+2) and tau = E / m_e c^2 (equation 1, Sternheimer 1984).

Parameters

tau	Reduced kinetic energy (E / m_e c^2) of the particle.
f	Array of nmax oscillator strengths.
eps	Array of nmax squared reduced oscillator energies (after adjustment factor has been applied).
nmax	Number of dispersion oscillators.
x	The trial value of l^2 (the variable being solved for).

Returns

The value of the objective function at x. A return value of zero indicates that x is the exact solution.

Definition at line 1119 of file estar.cpp.

References egsFatal.

bisec()

double bisec	(	double	lowerbound,
		double	upperbound,
		double	tolerance,
		double	tau,
		const std::vector< double > &	f,
		const std::vector< double > &	eps,
		int	nmax
	)

Find the root of objective_function() by bisection.

Locates l^2 such that objective_function(tau, f, eps, nmax, l^2) = 0, which gives the exact density correction via equation 1 of Sternheimer

1. The root must be bracketed by [lowerbound, upperbound] — i.e. the function must have opposite signs at the two bounds. The bracketing interval is provided by the lb_index and ub_index fields of the bspol struct returned by fbspol().

If the bounds do not bracket a root, an error is printed to stdout and lowerbound is returned as a sentinel value.

Parameters

lowerbound	Lower bound of the search interval (q[lb_index]).
upperbound	Upper bound of the search interval (q[ub_index]).
tolerance	Convergence tolerance on the absolute function value.
tau	Reduced kinetic energy (E / m_e c^2).
f	Array of nmax oscillator strengths.
eps	Array of nmax squared reduced oscillator energies.
nmax	Number of dispersion oscillators.

Returns

The value of l^2 at the root to within tolerance, or lowerbound if the bounds do not bracket a root.