Extension of a numerical algorithm to proton dose calculations. I. Comparisons with Monte Carlo simulations

A numerical algorithm originally developed for electron dose calculations [Med. Phys. 21, 1591 (1994)] has been modified for use with proton beams. The algorithm recursively propagates the proton distribution in energy, angle, and space from one level in an absorbing medium to another at slightly greater depth until all protons stop. Vavilov's theory is used to predict, at any point in the absorber, the broadening of the primary proton energy spectrum. Moliere's theory is applied to describe the angular distribution, and it is shown that the Gaussian first term of Moliere's series expansion is of sufficient accuracy for dose calculations. These multiple scattering and energy loss distributions are sampled using equal probability spacing to optimize computational speed while maintaining calculational accuracy. Inelastic nuclear collisions along the proton trajectories are modeled by a simple exponential extinction. Predictions of the algorithm for absolute dose deposition by a 160 MeV initially monoenergetic proton beam are compared with the results of Monte Carlo simulations performed with the PTRAN code. The excellent level of agreement between the results of these two methods of dose calculation (<5% dose and <3 mm spatial deviations) demonstrate that dose deposition from proton beams may be computed to high accuracy using this algorithm without the need for extensive empirical measurement as input.