Function theta_magic 

pub fn theta_magic(s: Screening, k: &MagicConstants, eps: f64, beta: f64) -> f64

Biersack-Haggmark “magic formula” (Nucl. Instrum. Methods 174 (1980) 257) for the centre-of-mass angle:

cos(theta/2) = (beta + rho + Delta) / (x0 + rho)
Delta = A (x0 - beta) / (1 + G)
A = 2 alpha eps beta^b,  alpha = 1 + C1 eps^(-1/2),
b = (C2 + eps^(1/2)) / (C3 + eps^(1/2)),
G = gamma / (sqrt(1 + A^2) - A),  gamma = (C4 + eps) / (C5 + eps)

with rho = 2 (1 - V(x0)/E) / (-V'(x0)/E) the radius of curvature of the trajectory at the distance of closest approach. Only the closest-approach root find is needed, no integral.

This form (no Delta in the denominator, gamma divided by the bracket) agrees with the MIT-licensed ir2-lab/screened_coulomb (commit f84c3c8). It is also the only arrangement that reproduces the quadrature to the accuracy expected of the fit. The tests check this to 1e-2 in cos(theta/2) for both constant sets.