May 2014 

NASA/TP-2014-218271

A 3DHZETRN Code in a Spherical Uniform 

Sphere with Monte Carlo Verification 

John W. Wilson 

Old Dominion University, Norfolk, Virginia 

Tony C. Slaba 

Langley Research Center, Hampton, Virginia 

Francis F. Badavi 

Old Dominion University, Norfolk, Virginia 

Brandon D. Reddell and Amir A. Bahadori 

Johnson Space Center, Houston, Texas 
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NASA Aeronautics 

and Space Database and its public interface, the 

NASA Technical Report Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in 

both non-NASA channels and by NASA in the 

NASA STI Report Series, which includes the 

following report types: 

 TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

Programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but 

having less stringent limitations on manuscript 

length and extent of graphic presentations. 

 TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

 CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

 CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or cosponsored by NASA. 

 SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

 TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission. 

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI 

program, see the following: 

 Access the NASA STI program home page 

at http://www.sti.nasa.gov

 E-mail your question to help@sti.nasa.gov

 Fax your question to the NASA STI 

Information Desk at 443-757-5803 

 Phone the NASA STI Information Desk at 

443-757-5802 

 Write to: 

 STI Information Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320 
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

May 2014 

NASA/TP-2014-218271

A 3DHZETRN Code in a Spherical Uniform 

Sphere with Monte Carlo Verification 

John W. Wilson 

Old Dominion University, Norfolk, Virginia 

Tony C. Slaba 

Langley Research Center, Hampton, Virginia 

Francis F. Badavi 

Old Dominion University, Norfolk, Virginia 

Brandon D. Reddell and Amir A. Bahadori 

Johnson Space Center, Houston, Texas 
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 
iii

Contents

Abstract .......................................................................................................................................................... 1

Introduction .................................................................................................................................................... 1

Deterministic Code Development .................................................................................................................. 2

3D Marching Procedures .............................................................................................................................. 12

A Simple Geometry and Source ................................................................................................................... 14

Monte Carlo Benchmarks ............................................................................................................................. 19

Conclusions .................................................................................................................................................. 24

Acknowledgements ...................................................................................................................................... 24

Appendix A: Monte Carlo simulations ......................................................................................................... 24

Appendix B: Evaluation of i for N = 6, 10, 14, 18, 22 ............................................................................... 25

References .................................................................................................................................................... 27
iv

Figures 

1. Geometric relations of quantities useful in solving equation (1). The symbol 

n is a unit normal vector.

 ................................................................................................................................................................. 2

2. Range-energy relationship of ions in aluminum. Results for Z = 1 and Z = 2 are nearly overlapping. ... 4

3. Probability of nuclear reaction of ions in aluminum. Results for Z = 1 and Z = 6 are nearly overlapping.

 ................................................................................................................................................................. 5

4. Isotropic (iso) and forward (for) neutron spectra produced by 500 MeV protons in aluminum. ............. 5

5. Dose in water target protected by 20 g/cm2 aluminum shield (left pane) and iron shield (right pane) 

exposed to Webber SPE spectrum (from Wilson et al. [1991]). .............................................................. 7

6. Dose equivalent [ICRP 1977] in water target protected by 20 g/cm2 aluminum shield (left pane) and 

iron shield (right pane) exposed to Webber SPE spectrum (from Wilson et al. [1991]). ........................ 7

7. Observed dose rate versus calculated dose rate [Badhwar, Chapter 21 of Wilson et al. 1997]. .............. 8

8. Comparison of fluence spectra at 35 g/cm2 (left pane) and 40 g/cm2 (right pane) in a 40 g/cm2

aluminum slab (infinite in lateral dimensions) exposed to Webber SPE spectrum. .............................. 11

9. Geometry of 3D marching procedure. ................................................................................................... 12

10. Spherical geometry and external source orientation for benchmark comparisons (not all detectors 

shown). .................................................................................................................................................. 14

11. Isotropic neutron source (E > 1 MeV) distribution induced by the Webber SPE spectrum in the 

spherical benchmark geometry. ............................................................................................................. 15

12. Directional components for N = 1, 2, 6, 10, 14, 18, and 22. .................................................................. 16

13. Particle spectra at detector locations 35 g/cm2 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2) exposed to Webber SPE spectrum for different N in 3D transport solution. ........ 18

14. The neutron leakage ratio at detector locations 35 g/cm2 (left pane) and 40 g/cm2 (right pane) from top 

of aluminum sphere (diameter 40 g/cm2

) exposed to Webber SPE spectrum for different N in 3D 

transport solution. The largest relative difference between the N = 10 and N = 22 ratios is less than 3% 

at both depths. ....................................................................................................................................... 19

15. Secondary neutron fluence at the interface between a 20 g/cm2 aluminum shield and 30 g/cm2 water 

target exposed to the Webber SPE (from Slaba et al. [2010b]). ............................................................ 19

16. Particle spectra at detector locations 35 g/cm2 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2

) exposed to Webber SPE spectrum. ....................................................................... 20

17. The neutron leakage ratio factor at detector locations 35 g/cm2 (left pane) and 40 g/cm2 (right pane) 

from top of aluminum sphere (diameter 40 g/cm2) exposed to the Webber SPE spectrum. .................. 21

18. Particle spectra at detector locations 35 g/cm2 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2

) exposed to the 1977 solar minimum GCR hydrogen spectrum. ........................... 22

19. Particle spectra at detector locations 35 g/cm2 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2

) exposed to the 1977 solar minimum GCR helium spectrum. The proton and alpha 

results have been scaled by 0.1 and 0.01 to improve plot clarity. ......................................................... 22

20. Particle spectra at detector locations 35 g/cm2 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2

) exposed to the 1977 solar minimum GCR carbon spectrum. The proton and alpha 

results have been scaled by 0.1 and 0.01 to improve plot clarity. ......................................................... 23

21. Particle spectra at detector locations 35 g/cm2 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2

) exposed to the 1977 solar minimum GCR iron spectrum. The proton and alpha 

results have been scaled by 0.1 and 0.01 to improve plot clarity. ......................................................... 23
v

Tables 

1. Comparison of straight-ahead HZETRN transport model with Shuttle flight data [Wilson et al. 2003b].

 ................................................................................................................................................................. 8

2. Average thickness traversed at different detector locations along z-axis of sphere. The left-most 

column is the distance from the top of the sphere to the detector. ........................................................ 17

3. Dose (cGy/event) at different detector locations along z-axis of sphere for Webber SPE. The left-most 

column is the distance from the top of the sphere to the detector. ........................................................ 17

4. Dose equivalent (cSv/event) at different detector locations along z-axis of sphere for Webber SPE. The 

left-most column is the distance from the top of the sphere to the detector. ......................................... 18

5. Total CPU seconds required for benchmark calculations. Total number of histories ran in Monte Carlo 

simulations are shown in parentheses. ................................................................................................... 21
1

Abstract 

The computationally efficient HZETRN code has been used in recent trade studies for lunar and 

Martian exploration and is currently being used in the engineering development of the next 

generation of space vehicles, habitats, and extra vehicular activity equipment. A new version 

(3DHZETRN) capable of transporting High charge (Z) and Energy (HZE) and light ions 

(including neutrons) under space-like boundary conditions with enhanced neutron and light ion 

propagation is under development. In the present report, new algorithms for light ion and neutron 

propagation with well-defined convergence criteria in 3D objects is developed and tested against 

Monte Carlo simulations to verify the solution methodology. The code will be available through 

the software system, OLTARIS, for shield design and validation and provides a basis for personal 

computer software capable of space shield analysis and optimization. 

Introduction 

High charge (Z) and Energy (HZE) ion transport developed parallel to the experimental studies of 

Schimmerling et al. [1986] at the Lawrence Berkeley Laboratory and involved solving the Boltzmann transport 

equation for mono-energetic ion beams in the context of the continuous slowing down approximation (CSDA) 

[Wilson et al. 1984]. The main purpose was to identify the source of computational limitations. This was found in 

part to be the inadequacy of available nuclear data through comparing computational results with the experimental 

ionization data for a broad beam of 20Ne ions [Wilson et al. 1984]. Emphasis was soon overtaken by the need to 

establish the scope of the galactic cosmic ray (GCR) protection problem, and marching procedures were used to get 

first order estimates of shielding requirements [Wilson and Lamkin 1975, Lamkin 1974, Wilson and Badavi 1986, 

Wilson et al. 1991]. Testing the new computational marching model against atmospheric air shower data again 

pointed to the inadequacy of the available nuclear data [Wilson and Badavi 1986], and development of a semiempirical nuclear fragmentation model (NUCFRG) followed, leading to improved comparisons with experimental 

deep penetration data [Wilson et al. 1987a, 1987b]. This new NUCFRG model likewise improved comparisons with 

Schimmerling’s laboratory studies [Shavers et al. 1993, Wilson et al. 1991]. The next dozen years emphasized 

spaceflight validation of the marching solution and nuclear fragmentation model improvements [Wilson et al. 1995, 

1998; Badhwar et al. 1996; Shinn et al. 1998; Cucinotta et al. 1998; Hugger et al. 2003; Nealy et al. 2006; Wilson et 

al. 2005, 2006, 2007; Slaba et al. 2011a, 2013; Norman et al. 2012, 2013]. 

Advanced solution methods of the Boltzmann equation continued to develop [Wilson et al. 1994a, 1994b] 

but only slowly after NASA interest in deterministic transport code and related nuclear model development shifted 

in favor of Monte Carlo methods [Armstrong and Colburn 2001, Pinsky et al. 2001]. Focus of the Langley Research 

Center activity turned to trade studies, flight validation, and development of engineering design methods including 

shield optimization processes for which deterministic methods are very successful [Wilson et al. 2004a]. Renewed 

interest within NASA for deterministic code development is giving new emphasis for improved solution methods 

[Wilson et al. 2003a, 2006; Tweed et al. 2006a, 2006b; Slaba et al. 2010a, 2010b; Norman et al. 2013] with a focus 

on nuclear modeling activity as well [Cucinotta et al. 1994, 1998, 2007; Norbury et al. 2007; Adamczyk et al. 2012]. 

As a result, recent developments mainly utilized the older semi-empirical NUCFRG2 model with revisions [Wilson 

et al. 1994c, 1995, 2006] including energy downshifts and momentum dispersion [Tripathi et al. 1994]. Comparison 

of NUCFRG2 with other models and nuclear fragmentation experiments are given by Zeitlin et al. [1997, 2011], 

Golovchenko et al. [1999], Cucinotta et al. [2007], and Walker et al. [2006]. The present development will use the 

more recent NUCFRG3 [Adamcyzk et al. 2012] with an added light ion (Z < 2) coalescence model and additional 

electromagnetic dissociation channels. 

To support the need for keeping records of astronaut exposures in space operations, considerable time has 

been spent in validating and testing in Shuttle flight operations with an array of instruments [Badhwar et al. 1995, 

1996, 2001]. As a result of the need to design and test new systems for future exploration, there has been renewed 

interest in flight validation using the International Space Station (ISS) as a measurements platform in low Earth orbit 

(LEO) [Hugger et al. 2003; Wilson et al. 2005, 2006, 2007; Nealy et al. 2006; Slaba et al. 2011a, 2013]. 

In the present report, the current status of transport code development is reviewed with emphasis on 

extending these developments into a 3D version of the HZETRN code. These advances will use available Monte 

Carlo codes Geant4 [Agostinelli et al. 2003, Geant4 Collaboration 2012a, 2012b], FLUKA [Fasso et al. 2005, 

Battistoni et al. 2007], and PHITS [Sato et al. 2013] to judge the veracity of these developments especially with 

regard to their 3D aspects. Upon completion, the 3DHZETRN code will be integrated into the web based OLTARIS 
2

[Singleterry et al. 2011] software system for general testing, validation against spaceflight data, and use in spacecraft 

design. 

Deterministic Code Development 

The relevant transport equations are the time-independent linear Boltzmann equations derived on the basis 

of conservation principles [Wilson et al. 1991] for the flux density (, , ) fj x W E of type j particles as 

4 ( , , ) ( , ', , ') ( , ', ') ' ' ( ) ( , , ) j jk k j j E k

E E E E d dE E E p f s f sf ¥= - åò ò W W WW W W W  x xx , (1) 

and solved subject to a boundary condition over the enclosure of the solution domain as shown in Fig. 1. 

Fig. 1. Geometric relations of quantities useful in solving equation (1). The symbol 

n is a unit normal vector. 

The σj(E) and σjk(E,E ',Ω,Ω') are the media macroscopic cross sections. The double differential cross section,

σjk(E,E ',Ω,Ω'), represents all those processes by which type k particles moving in direction Ω' with energy E ' 

produce a type j particle in direction Ω with energy E (including decay processes). Note that there may be several 

reactions that produce a particular product, and the appropriate cross sections for equation (1) are the inclusive ones. 

Exclusive processes are functions of the particle fields and may be included once the particle fields are known. The 

total cross section σj(E) for the medium may be written for each particle type as 

( ) () ( ) () () () () at el nuc

ss s s jj j j E EE E = ++ , (2) 

where the first term (at) refers to collision with atomic electrons, the second term (el) is for Coulomb elastic 

scattering on the atomic nuclei, and the third term (nuc) describes nuclear scattering and reactions where the minor 

nuclear inelastic processes (excitation) have been ignored. The corresponding differential cross sections are 

similarly ordered. Many atomic collisions (~106

) occur in a centimeter of ordinary matter, whereas ~103 nuclear 

Coulomb elastic collisions occur per centimeter, while nuclear processes are separated by a fraction to many 

centimeters depending on energy and particle type. Solution methods for equation (1) use perturbations based on 

the ordering of the physical cross sections with atomic interactions as the first physical perturbation; special methods 

are used for neutrons for which atomic cross sections are zero. Coulomb scattering is the second order physical 

perturbation process. The nuclear processes are treated in detail as the third order physical perturbation. 

The atomic/molecular interaction terms of equations (1) and (2) may be written as [Wilson et al. 1991] 
3

() ()

4

() () 1 1

,

() () () 1

, ,

( , ', , ') ( , ', ') ' ' ( ) ( , , )

( ) (, , ) ( ) (, , )

( ) (, , ) ( ) (, , ) (

at at

jk k j j E k

at at

jn j n j j n j j

n

at at at

jn j j n jn j j

n

E E E d dE E E

EA EA E E

E EA E E E

E

p

s f sf

s ef e s f

sf esf s

¥

- -

-

-

= + +-

=+ - ì ü ï ï ¶ é ù í ý ê ú ï ï ¶ ë û î þ

åò ò

å

å

WW W W W

W W

W W

x x

x x

x x 2

1 2

) (, , ) ( )

( ) ( , , ) ( ),

j n

j jj n

E O

A SE E O

E

f e

f e -

+

¶

= + é ù ë û ¶

W

W

x

x

 (3) 

where Aj is the mass of particle type j, and the stopping power Sj(E) is given by

( )

, () () atj n jn

n

SE E = åe s . (4) 

In equations (3) and (4), n denotes the electronic excitation levels, and εn refers to the corresponding excitation 

energies. Evaluation of the stopping power by equation (4) appears deceptively simple. However, all of the excited 

states including continuum states of the atomic/molecular system need to be evaluated. Furthermore, the projectile 

remains a bare ion except at low energies where the projectile ion atomic orbital states begin to resonate with the 

relative speed of the electrons of the media leading to electron capture and lowering of the ion charge (the stopping 

power is then related to RMS charge, called effective charge). These issues are further discussed in Wilson et al. 

[1991, 2005, 2006] and Tai et al. [1997]. The higher order terms of equation (3) are neglected in the CSDA and are 

further discussed elsewhere in connection to mono-energetic ion beams [Wilson et al. 2005, 2006]. The second 

physical perturbation, Coulomb scattering, results in angular diffusion about the ion direction leading to small 

angular corrections that are of little interest in the nearly omnidirectional space environment [Wilson et al. 1974, 

2005, 2006] and are therefore ignored. Henceforth, the total cross sections and the double differential cross sections 

appearing in equation (1) are the nuclear scattering and reaction values, and the superscript (nuc) will not be written. 

 Equation (1) can be written in the CSDA, ignoring multiple Coulomb scattering, as 

4 ( , , )] ( , ', , ') ( , ', ') ' ' j jk k E k

E E E E d dE p f sf ¥B[ W W W W W = å ò ò x x , (5) 

where the differential operator on the left hand side is defined as 

1 ( , , )] ( , , ) ( ) ( , , ) ( ) ( , , ) j j jj jj

j

E E SE E E E

A E f f f sf ¶ o - + é ù ë û ¶ B[ W W W W W xx xx  . (6) 

The right-hand side of equation (5) excludes the atomic/molecular processes now appearing on the left as the energy 

shifting operator in addition to the spatial drift and nuclear absorption terms. Neutral particles would have null 

atomic cross sections (except for negligible magnetic coupling) for which the stopping power term of B in equation 

(6) is set to zero. Application of CSDA in both laboratory and space shielding has been wide and the resulting errors 

discussed elsewhere [Tweed 2006a, 2006b]. Equation (5) can be rewritten as a Fredholm equation [Wilson 1977] 

given by 

' 4

( )( ) ( , , ) , ), , () ()

( ') ( ', '', , ') ( , ', '') ' '' ', () ()

j j

j j

j j

E j

j jk k E E

j j k

S E PE

E E

S EP E

A

P E E E E d dE dE

S EP E

g

g g

g

p

f f

s f ¥

= é ù ê ú ë û

+ åò òò W

W G(W W

WW W W

x x

x

 (7) 

with 
4

1

[ ( )] E R d RE g j j r - = -+ , (8) 

[ ( ) ( ')] =+ - RE RE j j x x W W . (9) 

In equations (7) and (8), Γ(Ω,x) is the point on the boundary connected to x along -Ω,  is the projection of x onto 

Ω (Fig. 1), d is the projection of Γ(Ω,x) onto Ω, Rj(E) is the distance a type j ion with energy E will travel before 

losing all of its energy to excitation/ionization of atomic electrons, and Pj(E) is the probability that a type j ion with 

energy E will not have a nuclear interaction in coming to rest in the media. The zero charge limit of equation (7) is 

discussed elsewhere [Wilson et al. 1991]. One obstacle to solving either equations (5) or (7) is the need to evaluate 

the integral dΩ' at arbitrary locations within the media and development of methods to handle this limitation 

expediently in ways amenable to efficient computational practice. 

The range-energy relationship is given (Fig. 2) by 

0

' ( ) ( ')

E

j j

j

dE RE A

S E = ò , (10) 

and the nuclear attenuation function is given (Fig. 3) by 

0

( ') ( ) exp ' ( ')

E j

j j

j

E

P E A dE

S E

é s ù

= -ê ú ê ú ê ú ë û ò , (11) 

where σj(E ') includes nuclear elastic and reactive processes. 

 The nuclear reactive and scattering differential cross sections can be written in the following form 

, ,

1 ( , ', , ') ( , ') ( , ', , ') 4 s ss jk jk iso jk for EE EE EE p WW WW = + , (12) 

where the first term (iso) is isotropic and associated with lower energy particles produced, including target 

fragments, and the second term (for) is highly peaked in the forward direction (Ω ≈ Ω') and is associated mainly 

with direct quasi-elastic events and projectile fragmentation products [Wilson 1977, Wilson et al. 1988]. 

Fig. 2. Range-energy relationship of ions in aluminum. Results for Z = 1 and Z = 2 are nearly overlapping. 
5

Fig. 3. Probability of nuclear reaction of ions in aluminum. Results for Z = 1 and Z = 6 are nearly overlapping. 

Fig. 4. Isotropic (iso) and forward (for) neutron spectra produced by 500 MeV protons in aluminum. 

Surprisingly, even nucleon-induced reactions follow this simple form, and the isotropic term extends to relatively 

high energies (see Fig. 4). For nucleon-induced reactions, the following form of the differential cross section has 

been used in versions of FLUKA [Ranft 1980] as follows 

( , ', , ') ( ') ( ') ( , ') ( , , ) s ns q jk jk k jk R T EE E E f EE g EA W W = , (13) 

where jk(E ') is average multiplicity, fjk(E,E ') is the energy spectrum of produced fragments, and the Ranft angular 

factor used in early versions of FLUKA is 

2

2

exp( / ) , 0 / 2 (, , ) exp( / 4 ) , /2

R R

R T

R R

N

g EA N

ql qp

q

p l p qp

ìï

ï - ££ = ï

í

ï

ï - < £ ïî

, (14) 

with cos θ = Ω·Ω' and 
6

(0.12 0.00036 ) / lR T = + A E . (15) 

In equations (14) and (15), AT is the mass of the target nucleus, and NR is a normalization constant. The Ranft 

formalism was modified by adding a quasi-elastic component [Wilson et al. 1988] and was instrumental in achieving 

the correct neutron KERMA [Wilson et al. 1991]. The isotropic component of equation (12) is determined by 

integrating equation (13) over the back hemisphere (i.e. cos(θ) < 0) according to 

, 2 ( , ') 2 ( , ', , ') ' jk iso jk B

EE EE d p s s = ò WW W , (16) 

where 2πB refers to the back hemisphere. The separation in phase space given by equations (12) - (16) is exploited 

in obtaining corrections to the straight-ahead approximation for light particles (Z < 2). The heavy ion (Z > 2) 

projectile fragment cross sections are further represented by 

2 2

, 2

( ') ( , ', , ') exp 2 '(1 ') / exp ( ' ) / 2

2

jk

jk t t jk jk jk

jk

E

EE N m E E E s s e le

pe

= - - --+ é ù é ù ê ú ê ú ë û ë û WW WW , (17) 

where m is the nucleon rest mass, jk is related to the momentum downshift, jk is related to the longitudinal 

momentum width, t,jk is related to the transverse momentum width, Nt is the transverse normalizing factor, and the 

fragmentation cross sections, σjk(E ') are obtained herein from NUCFRG3. Since the transverse width is small 

compared to the projectile-fragment energy, the transverse function is highly peaked about the forward direction

(Ω·Ω' ≈ 1). 

Atomic interactions limit the contributions of charged particles in the transport process. For example, the 

protons and alpha particles produced in aluminum below 100 A MeV contribute to the fluence only within a few 

centimeters of their collisional source and the heavier ions are even more restricted (see Fig. 2). This is an important 

factor in that the transported secondary charged particle flux tends to be small at low energies and the role of 

additional nuclear reactions by these low-energy ions are likewise limited (see Fig. 3), an important fact to be 

implemented later in the current development. 

Straight-ahead solution (N = 1) 

The approach to a practical solution of equation (5) is to develop a progression of solutions from the simple 

to increasingly complex, allowing early implementation of high-performance computational procedures and 

establishing a converging sequence of approximations with established accuracy criteria and means of 

verification/validation. The first step leading to the lowest order solution uses the straight-ahead approximation as 

guided by the nucleon transport studies of Alsmiller et al. [1965] using Monte Carlo methods in which the 

differential cross sections were approximated as 

( , ', , ') ( , ') ( ') s sd jk jk EE EE WW W-W = , (18) 

resulting in dose and dose equivalent per unit fluence to be within the statistical uncertainty of the Monte Carlo 

result obtained using the fully angle dependent cross sections in slab geometry. The relation of angular dependent 

cross sections to spacecraft geometry in space applications was further examined by Wilson and Kandelwal [1974] 

using asymptotic expansions about angular divergence parameters, demonstrating errors in the straight-ahead 

approximation to be on the order of the square of the ratio of distance of divergence (few centimeters) to radius of 

curvature of the shield (few to several m), resulting in a small error 10-4 in most human rated space systems 

[Wilson and Khandelwal 1974]. Introducing approximation (18) into equation (7) reduces the Fredholm equation to 

a Volterra equation that can be solved using marching procedures [Wilson and Lamkin 1975, Wilson 1977, Wilson 

and Badavi 1986]. Considering the success of the straight-ahead approximation for nucleons encourages its use in 

HZE transport where the mass of heavier ions reduces even further the angle of scatter. The success of equation (18) 

results from the fact that lateral diffusion from a given ray is fully compensated in a flat plate (i.e. slab with infinite 

lateral dimensions) by lateral diffusion along adjacent rays. The verification and validation processes are described 

elsewhere [Wilson et al. 2005, 2006]. 
7

Initial quantification of the usefulness of the code evaluated dose and dose equivalent (as given by the 

ICRP 26 [ICRP 1977]) in 30 cm slab of water behind a 20 g/cm2

 slab shield of aluminum (and alternately of iron) 

for the Webber [1966] approximation of the February 23, 1956 solar particle event (SPE) spectrum given as 

{ }

9 10 ( 938) ( ) exp 239.1 ( ) / 100 200 ( ) p

E E pE p E f + = - é ù ë û , (19) 

 

with pE EE ( ) ( 1876) = + . Comparisons were made against 3D Monte Carlo results (an early version of HETC) 

by Scott and Alsmiller [1967]. The results are shown in Fig. 5 for dose and Fig. 6 for dose equivalent. Statistical 

uncertainty of the Monte Carlo result is on the order of ten to fifteen percent, and the early version of HZETRN 

certainly lies within these bounds. 

Fig. 5. Dose in water target protected by 20 g/cm2

 aluminum shield (left pane) and iron shield (right pane) exposed 

to Webber SPE spectrum (from Wilson et al. [1991]). 

Fig. 6. Dose equivalent [ICRP 1977] in water target protected by 20 g/cm2

 aluminum shield (left pane) and iron 

shield (right pane) exposed to Webber SPE spectrum (from Wilson et al. [1991]). 
8

The newly developed dynamic LEO environmental models, shield models, and transport algorithm are 

shown to have an end-to-end uncertainty estimated to be 20 percent by Badhwar et al. [1995, 1996, 2001]. Aside 

from its accuracy, this version of HZETRN provides a computationally efficient procedure allowing rapid 

dosimetric evaluations and providing a basis for development of advanced engineering design methods capable of 

design optimization processes [Wilson et al. 2003a, 2004b]. A sample of results is given in Table 1 and Fig. 7. A 

sequence of corrections to the straight-ahead approximation is now taken up for light ions and neutrons. The heavy 

ion nuclear cross sections, represented by equation (17), indicate that the straight-ahead approximation is suitable for 

heavy ions for most applications, including space, as evidenced by laboratory validation with mono-energetic beams 

[Shavers et al. 1993, Shinn et al. 1998, Walker et al. 2006]. 

Fig. 7. Observed dose rate versus calculated dose rate [Badhwar, Chapter 21 of Wilson et al. 1997]. 

Table 1. Comparison of straight-ahead HZETRN transport model with Shuttle flight data [Wilson et al. 2003b]. 

Flight Date DRNM* DLOC TLD† (Gy/d) Calc. (Gy/d) 

STS-41A 11/83 6421 3 64.6 59.6 

STS-51D 4/85 6661 4 917.4 889.3 

STS-31 4/90 5701 1 2141 2290 

STS-43 8/91 5894 4 20.7 18.6 

STS-62 3/94 6771 1 94.3 89.2 

STS-65 7/94 6822 2 28.3 25.1 

STS-67 3/95 6925 3 250.8 238.1 

STS-80 11/96 6973 4 264.4 256.5 

STS-82 2/97 7074 1 2978 3080 

STS-91 6/98 6894 1 89.1 83.2 

STS-101 5/00 6460 2 140.8 131.1 

STS-92 10/00 6417 2 165.9 153.4 

* Deep River Neutron Monitor count rate, † GCR corrected TLD-100 data 

Bi-directional neutron transport (N = 2) 

It is clear that the straight-ahead approximation provides accurate dosimetric results and a good 

approximation to the particle fluence in many circumstances (see Badhwar et al. [1995]). One limitation of the 

straight-ahead approximation is near the front boundary on a flat plate where the calculated neutron fluence vanishes 

(unless they are part of the external radiation environment). This results from the straight-ahead approximation itself 

and can be improved by connecting to 3D corrections containing backward leakage to the front boundary [Marchuk 

and Lebedev 1986]. This limitation is resolved by treating the leakage that is an important process of neutron 

transport in space systems but of less importance to the charged particle transport that is limited by range/energy 
9

relations. This bi-directional approximation provides the first correction to the straight-ahead approximation and 

provides major improvements to the neutron fluence spectra. 

In the early transport studies of Wilson and Lamkin [1975] and Lamkin [1974], it was demonstrated that 

once neutrons are produced by a proton beam, the recoupling of neutrons back to the proton field is limited, but the 

neutron fields tend to build up with increasing penetration depth. This mainly occurs since many protons produced 

by neutron induced reactions are of lower energy than the forward components of equation (12) and are quickly 

removed by atomic interactions (especially these isotropically produced protons). The probability of even a 100 

MeV proton to undergo a nuclear reaction is only several percent (Fig. 3). The low energy neutrons on the other 

hand have no effective atomic interactions and propagate until a nuclear interaction occurs or they escape from the 

media through a boundary. Yet, these low energy protons and other light ions produced in tissue through neutron 

collisions contribute significantly to biological injury and must be included. This was the basis of studies starting 

with Clowdsley et al. [2000, 2001] and Heinbockel et al. [2003], with recent improvements provided by Slaba et al. 

[2010b]. Herein, we will concentrate on its extension to 3D transport methodologies. Future work will be a direct 

evaluation of equation (7) quantifying the propagated error of the 3D neutron code. 

 A bi-directional neutron transport code was first investigated by Clowdsley et al. [2000, 2001] using 

multigroup methods and recoupling to the light ions using an analytic solution as a quadrature over the light ion 

collision source term. These types of approximations were further studied and improved by Slaba et al. [2010b] and 

later compared to Monte Carlo simulations in slab geometry [Heinbockel et al. 2011a, 2011b] and on the lunar 

surface [Slaba et al. 2011b]. The model has been implemented into OLTARIS and is the basis of the ray-by-ray 

transport methodology [Slaba et al. 2011b, Singleterry et al. 2011]. In that transport model, the charged particle 

components were first solved in the straight-ahead approximation with the neutron coupling represented only by the 

forward component of the differential production cross section. The isotropically produced neutrons are then solved 

using a coupled bi-direction (forward/backward) transport to complete the neutron solution. These isotropically 

produced neutrons further collide with the media, giving rise to an isotropic source of ions that are transported 

assuming only atomic interactions as an approximation to the complete solution. Neglected in this solution is the 

minor coupling between low energy charged particles and any produced fragments. This assumption is generally 

accurate since the range of the low energy isotropic ions is much smaller than their nuclear mean free path lengths 

(see Figs. 2 and 3). 

The governing equation in this approach is first given by equation (5). The solution methodology requires 

particle fluxes to be separated into forward (φj,for) and isotropic components (φj,iso) according to 

, , (, , ) (, , ) (, , ) ff f j j for j iso x xx W WW E EE = + . (20) 

The differential cross section for neutron fragments (j = n) is also separated into forward and isotropic components 

according to equations (12) and (16). The differential cross section for charged fragments is treated in the straightahead approximation by equation (18). Upon substitution of the separated differential cross sections and particle 

fluxes into equation (5), the forward component of the charged particle (j ≠ n) fluxes is obtained by solving the 

straight-ahead equation 

, , ( , , )] ( , ') ( , , ') ' j for jk k for E k

f sf E E E E dE ¥B[ W W = åò x x , (21) 

and the forward neutron component is obtained by solving 

, ,, ( , , )] ( , ') ( , , ') ' n for nk for k for E k

f sf E E E E dE ¥B[ W W = åò x x , (22) 

for which the stopping power term in the operator B is zero for neutrons. Note that the full cross section of equation 

(18) appears in the ion transport of equation (21), but only the forward component of the cross sections appears in 

the transport of neutrons in equation (22). Equations (21) and (22) are solved simultaneously using methods 

previously developed [Wilson et al. 2006, Slaba et al. 2010a]. 

The forward fluxes give rise to an isotropic source of low energy neutrons which are treated by solving the 

equation 
10

, , 4

, , 4

( , , )] ( , ', , ') ( , ', ') ' '

( , ', , ') ( , ', ') ' ' .

n iso nn n iso E

nk iso k for E k

E E E E d dE

E E E d dE

p

p

f sf

s f

¥

¥

=

+

ò ò

åò ò

B[ W W W W W

WW W W

x x

x

 (23) 

Whereas equations (21) and (22) take on value only if an external source is present at the boundary, equation (23) 

has no external sources, and the resulting field develops from the internal isotropic neutron source induced by the 

forward solution. Development of an algorithm for solving equation (23) is expedited by separating the isotropic 

neutron fluence into forward and backward propagating components according to 

, , ( ) , () (, , ) (, , ) (, , ) ff f n iso n iso f n iso b x xx WW W EE E = +- , (24) 

and approximating the double differential cross sections as 

, , ( , ', , ') ( , ') ( ') ( , ') ( ') s sd sd nn nn f nn b EE EE EE WW W-W W+W = + , (25) 

, ,, 2 ( , ') ( , ', , ') ' ( , ') 2 ( , ') / 4 nn f nn nn for nn iso F

EE EE d EE EE p s s s ps p = =+ ò WW W , (26) 

, , 2 ( , ') ( , ', , ') ' 2 ( , ') / 4 nn b nn nn iso B

EE EE d EE p s s ps p = = ò WW W , (27) 

, () , , 2 ( , ') ( , ', , ') ' 2 ( , ') / 4 nk iso f nk iso nk iso F

EE EE d EE p s s ps p = = ò WW W , (28) 

, () , , 2 ( , ') ( , ', , ') ' 2 ( , ') / 4 nk iso b nk iso nk iso B

EE EE d EE p s s ps p = = ò WW W . (29) 

In equations (26) - (29), 2F and 2B represent integration over the forward and backward hemispheres with respect 

to Ω', respectively. The subscript (f) and (b) in equation (24) refers to the forward hemisphere and backward 

hemisphere components of the isotropic neutron flux, respectively. Using equations (25) - (29) in equation (23) 

results in 

, () , , ()

, , ()

, () ,

, () ,

[ ( )] ( , , ) ( , ') ( , , ') '

( , ') ( , , ') '

( , ') ( , , ') '

( , ') ( , , ') '

n n iso f nn f n iso f E

nn b n iso b E

nk iso f k for E k

nk iso b k for E k

E E E E E dE

E E E dE

E E E dE

E E E dE

sf s f

s f

s f

s f

¥

¥

¥

¥

+ =

+ -

+

+ -

ò

ò

åò

åò

WW W 

W

W

W

x x

x

x

x

 (30) 

for the forward component of the isotropic neutron flux. The backward moving isotropic neutron field satisfies a 

similar equation 

, () , , ()

, , ()

, () ,

, () ,

[ ( )] ( , , ) ( , ') ( , , ') '

( , ') ( , , ') '

( , ') ( , , ') '

( , ') ( , , ') '

n n iso b nn f n iso b E

nn b n iso f E

nk iso f k for E k

nk iso b k for E k

E E E E E dE

E E E dE

E E E dE

E E E dE

sf s f

s f

s f

s f

¥

¥

¥

¥

- + =

+ -

+ -

+

ò

ò

åò

åò

WW W 

W

W

W

x x

x

x

x

 (31) 
11

Note that , ( , , ') 0 fk for x - = W E in equations (30) and (31). Equations (30) and (31) represent a coupled set of 

differential equations that are solved efficiently using a Neumann series approach as shown by Slaba et al. [2010b]. 

This series is expected to converge quickly as σnn,f (E,E') ≫ σnn,b(E,E'). 

The isotropically produced neutrons initiate isotropic light ion production in collisions with the media 

resulting in 

, , ( , , )] ( , ') ( , , ') ' j iso jn n iso E f sf E E E E dE ¥

B[ W W = ò x x , (32) 

for which a semi-analytic solution can be efficiently implemented [Slaba et al. 2010b]. Equations (30) and (31) are 

solved simultaneously for no external fluence incident on the media directed along Ω (see Fig. 1). Note that the 

implicit assumption of equations (30) and (31) is that lateral diffusion from the ray Ω is fully compensated by 

diffusion from adjacent rays. This is strictly true for a flat plate slab with infinite lateral dimensions as will now be 

demonstrated and represents the use of periodic boundary conditions as an importance sampling method common to 

Monte Carlo simulation practice. 

 These two methods (straight-ahead (N = 1) and bi-directional (N = 2)) are now compared to values 

simulated by the Monte Carlo codes Geant4, FLUKA, and PHITS for the Webber SPE spectrum normally incident 

on a 40 g/cm2

 slab of aluminum with infinite lateral dimensions. All Monte Carlo fluence results presented herein 

are plotted as mean values with estimated statistical uncertainty shown as vertical error bars placed at the energy bin 

midpoint. The points are connected with lines to improve plot clarity. If not visible, vertical error bars are smaller 

than the symbols on the plot. Details of the Monte Carlo simulation setups are given in Appendix A. Results are 

shown in Fig. 8. The straight-ahead solution gives good agreement for charged components at most depths. Compare 

the left and right panes of Fig. 8. The neutron fluence reaches a broad maximum across 20 to 35 g/cm2

 wherein the 

neutron fluence is nearly isotropic at low energies and the straight-ahead (N = 1) fluence approaches the bidirectional (N = 2) fluence from below before the maximum and over estimates beyond the maximum. On the back 

face, leakage is an important process leaving the straight-ahead fluence as an over estimate at low neutron energies 

where the neutron fluence is reduced by about half of its interior value due to leakage through the back surface 

(compare left and right panes of Fig. 8). The resonance structure seen in the Monte Carlo neutron fluences below 1 

MeV is a result of elastic cross section databases that have not been included in 3DHZETRN. 

Fig. 8. Comparison of fluence spectra at 35 g/cm2

 (left pane) and 40 g/cm2 (right pane) in a 40 g/cm2 aluminum slab 

(infinite in lateral dimensions) exposed to Webber SPE spectrum. 
12

It is clear that the bi-directional solution for the semi-infinite slab where lateral leakage is negligible is a reasonable 

approximation. The success of the straight-ahead approximation in validation in human rated systems (Shuttle, ISS) 

clearly meets the requirements for a useful solution of the transport problem. Attention is now given to transport in 

finite objects in which lateral leakage plays an important role. 

3D Marching Procedures 

A sequence of approximations is now sought that will span the formalism from the simple bi-directional 

approximation to more complex propagation algorithms. With the demonstrated success of the bi-directional 

approximation, the solution is divided into forward and isotropic components, as before. The governing equation in 

the CSDA ignoring Coulomb scattering is given by equation (5). The particle flux is separated into forward and 

isotropic components, as in equation (20). Similarly, the double differential cross sections for neutron production are 

separated into forward and isotropic components given by equations (12) and (16). The charged particle production 

cross sections are assumed to take the form of equation (18). Upon substitution of the separated fluxes and 

production cross sections into equation (5), the forward components are obtained by solving 

, , ( , , )] ( , ') ( , , ') ' j for jk k for E k

f sf E E E E dE ¥B[ W W = åò x x , (33) 

for charged particles (j ≠ n), and for neutrons (j = n), 

, ,, [ ( )] ( , , ) ( , ') ( , , ') ' n n for nk for k for E k

sf s f E E E E E dE ¥+ = å ò WW W  x x . (34) 

 The space environment (external source) is nearly uniform in position and varies slowly in angle (if not 

isotropic), and for the forward flux, each directional component of the boundary condition may be propagated 

according to equations (33) and (34). Note that Ω enters equations (33) and (34) as a parameter and can be evaluated 

according to the boundary condition we have set as the incident flux from direction Ω0 (see Fig 9). Once the forward 

flux is evaluated for direction Ω0 over the domain of x and E, one must yet solve for the induced isotropic neutron 

flux within the shield configuration. Coupling to the isotropically produced neutron field is taken as 

, , ,0 4 [ ( )] ( , , ) ( , ', , ') ( , ', ') ' ' ( , , ) n n iso nn n iso n iso E

E E E E E d dE E p sf s f x ¥+ = + ò ò W W W W W W W, W  x xx , (35) 

where the isotropic neutron source is given by 

, 0 , 0, 0 ( , , ) ( , ', , ) ( , , ') ' n iso nk iso k for E k

x sf E E E E dE ¥W, W W W W = å ò x x . (36) 

Fig. 9. Geometry of 3D marching procedure. 
13

Note that the isotropic neutron source is generated by collisions between the forward flux and target media, for 

which the forward flux is a function of only the penetration depth z (Fig. 9) along the direction Ω0. This is denoted 

by writing the isotropic neutron source as 

, 0, 0 [ ( ), , ] ( , , ) n iso n iso c x zE E x x W, W W, W = , (37) 

where z(x) is the penetration depth along Ω0 to the point x. Equation (35) is now rewritten as 

, , ,0 4 [ ( )] ( , , ) ( , ', , ') ( , ', ') ' ' [ ( ), , ] n n iso nn n iso n iso E

E E E E E d dE z E p sf s f c ¥+ = + ò ò W W W W W W W, W  x xx (38) 

Equation (38) also provides a basis for introducing anisotropic source terms as well. Note that while the forward 

propagating solution has its source fixed by the boundary condition, the isotropic solution (having no inbound 

fluence) is driven by internal sources generated by the collisions of the inbound fluence (as modified by the forward 

propagator in penetrating to x along direction Ω0). 

 To solve equation (38), the neutron fluence , ( , ', ') fn iso x W E must be known at every x from all directions Ω'

approaching x from adjacent locations as was the case of equation (23). Although the use of approximation (25), 

based on the assumption that losses along the ray Ω are compensated by diffusion from adjacent rays, is a 

demonstrated useful assumption in a semi-infinite slab, it is an approximate overestimate in a sphere. Yet, using 

equation (25) in evaluation of equation (38) seems a reasonable approximation in most human rated vehicles leading 

to an extension of transport along an arbitrary ray Ω as a solution to the bi-directional equations. The forward 

component of the isotropic neutron field satisfies 

, () , , ()

, , ()

, 0

( )] ( , , ) ( , ') ( , , ') '

( , ') ( , , ') '

[ ( ), , ] ,

n n iso f nn f n iso f E

nn b n iso b E

n iso

E E E E E dE

E E E dE

z E

sf s f

s f

c

¥

¥

+ =

+ -

+

ò

ò

[W W W 

W

W, W

x x

x

x

 (39) 

and the backward moving field satisfies 

, () , , ()

, , ()

, 0

( )] ( , , ) ( , ') ( , , ') '

( , ') ( , , ') '

[ ( ), , ] .

n n iso b nn f n iso b E

nn b n iso f E

n iso

E E E E E dE

E E E dE

z E

sf s f

s f

c

¥

¥

+ =

+ -

+ -

ò

ò

[-W W W 

W

W, W

x x

x

x

 (40) 

It is clear that Ω enters equations (39) and (40) as a parameter and can be solved along an arbitrary number 

of rays used to construct a numerical representation of , (, , ) fn iso x W E whose Ω dependence comes from distance to 

the boundary t(Ω) and the penetration depth z along Ω0. The isotropic neutron source of equation (38) effectively 

decouples the incident direction Ω0 from the neutron direction of propagation Ω. 

These isotropically produced neutrons produce isotropic light ions in collisions with the media resulting in 

a light ion fluence given by 

, , ( , , )] ( , ') ( , , ') ' j iso jn n iso E f sf E E E E dE ¥

B[ W W = ò x x . (41) 

Equations (39) and (40) are solved simultaneously for no external fluence incident on the media directed along Ω

(see Fig. 9). Nonetheless, efficient numerical methods already developed for the bi-directional model may be 

utilized. Note also that the implicit assumption of equations (39) and (40) is that lateral diffusion from the ray Ω is 

fully compensated by diffusion from adjacent rays. This has been demonstrated for slab geometry with results in 
14

Fig. 8. The first order lateral leakage in this case is found by solving equations (39) and (40) along each Ωi over the 

ray limits {-t(-Ωi), t(Ωi)} for a given set of N directions Ωi. The next sections address convergence in the isotropic 

neutron and ion fields as a function of the number of directions (N) chosen to represent the fields. Numerical 

solutions for equations (39) - (41) are the same as in the bi-directional solution and described by Slaba et al. 

[2010b]. The only difference here is that the isotropic source term in each case is computed as a function of 

penetration depth with respect to the incident source direction, Ω0. 

 The present transport formalism is related, but quite different from SN transport theory [Marchuk and 

Lebedev 1986]. The main difference between the two formalisms is that solutions along a set of directions Ωi are 

fully coupled to all other discrete directions in the SN treatment. While in the present formalism, coupling is only 

included through the isotropic source term and for transport along anti-parallel rays (i.e. Ωi and -Ωi). It is 

demonstrated below that the present formalism compares well to Monte Carlo simulations where such coupling is 

implicitly included. 

A Simple Geometry and Source 

Simple geometry and source orientation are chosen to allow efficient Monte Carlo simulation for 

comparison with the calculation of quantities with the above formalism. In this way, the solution methodology is 

verified and experience is gained in the formalism in preparation for connecting to the general geometric 

descriptions generated by the engineering design process for which Monte Carlo simulation is impractical. The 

geometry herein is taken as a simple sphere of aluminum of 40 g/cm2

 diameter for comparison with the appropriate 

solutions of the current formalism. Target points (detectors) are placed at depths along the z–axis (x = 0, y = 0). The 

external radiation environment (external source boundary condition) is assumed to be anti-parallel with the z-axis, 

uniform in the x-y plane, positioned above the sphere, and directed down onto the top of the sphere as shown in Fig. 

10. 

Fig. 10. Spherical geometry and external source orientation for benchmark comparisons (not all detectors shown). 

The object geometry is specified at a point x where field values are to be evaluated, about which the 

thickness of material t is known over an array of directions Ωi (i.e. ti = t(Ωi)). Note that each Ωi for i odd is backed 

by i+1 = -i. The remainder of the solution methodology lies in selection of the arbitrary directions of Ωi of 

propagation within the target material and evaluation of the penetration depth z for any transport depth xT along each 

Ωi for the given Ω0 at the boundary. The N directions are taken as the central ray of an equal division of the unit 

sphere such that Ω = 4 / N. In this, the interest is in evaluating the omni-directional field function at an arbitrary 

location x about which the surrounding material is known as a thickness distribution of the various materials (in the 

present only homogenious aluminum). The evaluation will be made for a single direction of incidence Ω0 to 

highlight 3D transport effects, but this could include an array of incident directions for which the solutions would be 

integrated over all Ω0 values. 
15

Fig. 11. Isotropic neutron source (E > 1 MeV) distribution induced by the Webber SPE spectrum in the spherical 

benchmark geometry. 

In that low-energy neutron transport is dominated by diffusive processes, it is anticipated that the solution 

will rapidly converge for low values of N (i.e. less than 20). The source of neutrons induced within the sphere by the 

Webber SPE is shown in Fig. 11. The resulting neutron fluence depends not only on the distance t(Ω) to the 

boundary, but the intensity of the source along the ray in reaching the boundary. Clearly for SPEs, the most intense 

regions must be adequately represented in the domain of the solution. The next three sub-sections describe how the 

penetration depth, z, is determined for N = 1, 2 and the general case (N > 2). 

N = 1, Straight-ahead approximation

For each Ω0, only one value of Ωi is allowed and is taken to be that of the incident direction Ω0. The 

transport distance x along Ω1 and penetration depth z are the same, reducing complication in the formulation for N = 

1. In this case, leakage from the media is only through the distal face. The structure of this code is quite different 

than the original HZETRN code, although computational results are identical, and the original HZETRN code 

provides a verification test-case of the new 3DHZETRN code. For N = 1, the comparisons in the prior sections 

describing the verification and validation of HZETRN is fully applicable. If one restricts interest to N = 1, then the 

HZETRN code is preferred in that computational efficiency is higher but the added overhead of 3DHZETRN for N

= 1 is minimal. In the following sections, comparisons between 3DHZETRN and Monte Carlo simulations will also 

include results from the HZETRN code (N =1) to highlight shortcomings of the straight-ahead approximation and 

improvements of the 3D corrections presented herein. 

N = 2, Bi-directional approximation

The first step towards a full 3D code is the treatment of backward propagating particles. This was first 

accomplished by Clowdsley et al. [2000, 2001] and later improved by Slaba et al. [2010b] with demonstrated 

improved results [Slaba et al. 2010b, 2011b; Heinbockel et al. 2011a, 2011b]. In this case, propagation of the 

isotropic neutrons is represented in two directions by Ω1 = Ω0 and Ω2 = -Ω0. Therefore, the penetration depth z is 

equal to x along Ω1 and equal to t(Ω1)+ t(Ω2) - x along Ω2. Note that lateral leakage is not represented in this 

approximation although leakage at the forward and backward boundaries is so represented. The neutron transport in 

semi-infinite slab geometry has extensive Monte Carlo results for verification for which this two-stream 

approximation, ignoring lateral leakage, is most appropriate as demonstrated in Fig. 8. This approximation is 

appropriate for slab geometry and works well in applications in which the shield geometry is a collection of thin 

plates. The next step for treatment of lateral leakage is now discussed. 
16

N > 2, N-stream approximation 

For general N, the transport formalism requires the penetration depth, zp, to be known as a function of 

transport depth, xT, along the propagation direction Ω for each Ω0. In order to simplify algebraic manipulation and 

expressions, the center of the sphere is assumed to be positioned at the origin with the positive z-axis anti-parallel to 

the source direction, as shown in Fig. 10. Due to the simplified geometry and source orientation, an analytic 

expression for the penetration depth may be determined. The point x0 may be written as 

0 ) T x x =- - W(t x , (42) 

where t is the known thickness along direction -Ω from x to the boundary. The point xΓ may then be written in terms 

of x0 as 

0 0 p z x x G = - W . (43) 

The expression for zp is obtained by writing equation (43) in component form and substituting into the equation of a 

sphere with radius R = 20 g/cm2

. This yields 

2 2 2

00 0 00 ( ) p z R = +-+ x xx   W W . (44) 

 For preliminary testing and benchmark comparisons against Monte Carlo codes, the solutions for N = 6, 10, 

14, 18, and 22 are considered separately, in addition to the simple N = 1 and N = 2 cases discussed above. The 

distributions for all N are shown in Fig 12. Note that the N = 6, 14, and 22 cases include rays perpendicular to the zaxis (i.e. tangent to the sphere surface at the top and bottom detector locations). On the surface of the sphere, these 

rays correspond to zero shield thickness. This is corrected in the N = 10 and N = 18 distributions. The directional 

elements for each N are given in the Appendix B, along with a brief description of how the elements were 

determined. For each N, the omni-directional flux is obtained by integrating the fluxes along each ray direction after 

the transport procedure. 

Fig. 12. Directional components for N = 1, 2, 6, 10, 14, 18, and 22. 

Preliminary testing

A first quantity of interest is the average path length at each detector location where the solution is 

evaluated to the sphere boundary. For a given N, the average path length is computed over the material thicknesses 

traversed from the target point to the boundary along all rays. This quantity is expected to converge rapidly at most 

locations with increasing N. As shown in Table 2, convergence is reasonably rapid for all locations except on the 

sphere surface. It is clear that the results are reasonably converged at N = 18, especially in comparison to N = 10, as 

will be further demonstrated for transport quantities. 
17

Table 2. Average thickness traversed at different detector locations along z-axis of sphere. The left-most column is 

the distance from the top of the sphere to the detector. 

d (g/cm2

) Average thickness traversed (g/cm2

) 

N = 1 N = 2 N = 6 N = 10 N = 14 N = 18 N = 22

0 0.0 20.0 6.7 10.4 9.4 10.1 9.8 

0.3 0.3 20.0 9.0 10.9 10.6 10.7 10.7 

5 5.0 20.0 15.5 15.6 15.6 15.7 15.7 

10 10.0 20.0 18.2 18.2 18.2 18.2 18.2 

15 15.0 20.0 19.6 19.6 19.6 19.6 19.6 

20 20.0 20.0 20.0 20.0 20.0 20.0 20.0 

25 25.0 20.0 19.6 19.6 19.6 19.6 19.6 

30 30.0 20.0 18.2 18.2 18.2 18.2 18.2 

35 35.0 20.0 15.5 15.6 15.6 15.7 15.7 

40 40.0 20.0 6.7 10.4 9.4 10.1 9.8 

The sequence of solutions for N = 1 to 22 is used as a test on the convergence of the solution method. In 

addition, a Monte Carlo benchmark is defined for each of the indicated methods. Benchmarks for N = 1 and N = 2 

are given in slab geometry in Fig. 8 to be later compared with N = 18 benchmarks in spherical geometry. 

The first quantity of interest is the dose associated with the Webber SPE as given in Table 3 for all N. Dose 

is mainly driven by the proton fluence that still has a large contribution from the penetrating primary protons that 

improves the convergence rate (especially at shallower depths near the top of the sphere). The dose equivalent, 

receiving a large contribution from the other light ions for which the quality factor is large, will require convergence 

of all the light ions to ensure an accurate value for dose equivalent. The values of dose equivalent for the Webber 

SPE are given in Table 4 for various numbers of rays in the 3D evaluations. Again, rapid convergence is observed 

with increasing N. As expected, the dose equivalent converges most slowly near the bottom of the sphere. A second 

look at the convergence of the particle fluence spectra are shown in Fig. 13. The convergence of the fluence spectra 

at 35 g/cm2

 is typical for an interior location for which N = 18 is clearly well converged. Even on the back surface 

of the sphere at 40 g/cm2 the difference between N = 10 and N = 18 is very small in the graph even for alpha fluence; 

convergence is effectively obtained with N = 18. The rate of convergence of the alpha fluence mainly affected the 

convergence of the dose equivalent in Table 4. The N = 18 will be used in the Monte Carlo benchmarks discussed in 

the next section. 

 

Table 3. Dose (cGy/event) at different detector locations along z-axis of sphere for Webber SPE. The left-most 

column is the distance from the top of the sphere to the detector. 

d (g/cm2

) Dose (cGy/event) 

N = 1 N = 2 N = 6 N = 10 N = 14 N = 18 N = 22

0 43703.70 43703.77 43703.73 43703.84 43703.79 43703.82 43703.80 

0.3 1216.22 1216.28 1216.33 1216.35 1216.34 1216.34 1216.34 

5 69.95 69.98 70.03 70.02 70.02 70.02 70.02 

10 24.80 24.81 24.83 24.82 24.82 24.82 24.82 

15 12.23 12.23 12.23 12.23 12.23 12.23 12.23 

20 7.00 6.99 6.98 6.98 6.98 6.98 6.98 

25 4.37 4.36 4.34 4.34 4.34 4.34 4.34 

30 2.90 2.88 2.85 2.85 2.85 2.85 2.85 

35 2.01 1.98 1.95 1.95 1.95 1.95 1.95 

40 1.44 1.41 1.38 1.38 1.38 1.38 1.38 
18

Table 4. Dose equivalent (cSv/event) at different detector locations along z-axis of sphere for Webber SPE. The leftmost column is the distance from the top of the sphere to the detector. 

d (g/cm2

) Dose equivalent (cSv/event) 

N = 1 N = 2 N = 6 N = 10 N = 14 N = 18 N = 22

0 737745.84 737746.95 737746.21 737748.11 737747.20 737747.89 737747.45 

0.3 2871.17 2872.15 2873.07 2873.31 2873.10 2873.17 2873.12 

5 109.95 110.33 111.15 111.05 111.02 111.00 111.00 

10 37.94 38.10 38.45 38.39 38.36 38.38 38.38 

15 19.01 19.01 19.04 18.97 19.00 19.01 18.99 

20 11.33 11.20 10.97 10.97 10.98 10.97 10.98 

25 7.54 7.27 6.85 6.92 6.90 6.90 6.91 

30 5.42 5.03 4.58 4.55 4.59 4.58 4.58 

35 4.15 3.63 3.18 3.13 3.13 3.15 3.15 

40 3.34 2.69 2.26 2.21 2.19 2.19 2.21 

Fig. 13. Particle spectra at detector locations 35 g/cm2

 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2

) exposed to Webber SPE spectrum for different N in 3D transport solution. 

It is clear in the prior discussion that the main difference between N = 2 and N > 2 is the degree of lateral 

diffusion and the resulting leakage near boundaries. The ratio 

( , 2) () 1 ( , 2)

n

N

n

E N

r E

E N

f

f

> o - = (45) 

is directly related to the representation of lateral losses for neutrons. The quantity rN(E) is fixed by the shape of the 

object, the scattering angular distribution, and the fidelity of the transport solution represented by N. The quantity

rN(E) is shown in Fig. 14 for N = 6, 10, 14, 18, and 22. The value of zero at the highest energy is a result of the 

straight-ahead approximation for the forward component for which there is no lateral leakage. The largest relative 

difference between the N = 10 and N = 22 ratios is less than 3% at both depths, and the N = 18 solution appears 

converged based on the results in Tables 2-4 and Figs. 13 and 14. 
19

Fig. 14. The neutron leakage ratio at detector locations 35 g/cm2

 (left pane) and 40 g/cm2 (right pane) from top of 

aluminum sphere (diameter 40 g/cm2

) exposed to Webber SPE spectrum for different N in 3D transport solution. 

The largest relative difference between the N = 10 and N = 22 ratios is less than 3% at both depths. 

Monte Carlo Benchmarks 

In 2005, a systematic benchmark exercise was administered as a means of verifying various computational 

procedures resulting from the first phase of code development activity (1996-2004). Code verification uses 

established benchmarks as consisting of standard environments: the Webber SPE spectrum given in equation (19) 

and the 1977 solar minimum GCR environment. 

In earlier SPE benchmarking, the configuration was for a 30 cm semi-infinite slab of water shielded by 20 

g/cm2

 of either aluminum or iron utilizing only dose and dose equivalent (ICRP 26 quality factor [ICRP 1977]) as 

given in Figs. 5 and 6. For the moment, focus is mainly given to the results for the neutron spectra at the 

aluminum/water interface as shown in Fig. 15. It is seen from the figure that those codes with roots in the HETC 

code (HETC-HEDS and MCNPX) form a cluster of results including HZETRN at low energies (E < 40 MeV). In 

distinction, Geant4 and FLUKA form a cluster below HZETRN. All codes show improved agreement above 100 

MeV. 

Fig. 15. Secondary neutron fluence at the interface between a 20 g/cm2

 aluminum shield and 30 g/cm2 water target 

exposed to the Webber SPE (from Slaba et al. [2010b]). 
20

It was generally concluded that HZETRN agreeds with the Monte Carlo codes as well as the various Monte 

Carlo codes agreed among themselves. There are currently no compelling reasons to attempt to use Monte Carlo 

methods in design and operational applications. However, it is further demonstrated below, the Monte Carlo codes 

provide an indispensable guide to development of the 3DHEZTRN code (that does provide a useful basis for space 

radiation protection). Focus is now shifted towards benchmarking 3DHZETRN, in which the added complications of 

finite geometry are addressed. Details associated with Monte Carlo simulations presented in Figs. 16 - 21 are given 

in Appendix A. One needs to be mindful of the intercode comparisons of Fig. 15 when viewing succeeding 

benchmarks. 

Webber SPE Benchmark

The benchmark SPE spectrum is the Webber [1966] approximation of the February 23, 1956 SPE (not an 

accurate representation but a historically useful test spectrum, see Foelsche et al. [1974]) as given in equation (19). 

The induced neutron and proton fluence spectra were evaluated in a 40 g/cm2

 diameter aluminum sphere by Monte 

Carlo codes (Geant4, FLUKA, and PHITS) and 3DHZETRN with the results at 35 and 40 g/cm2 depths shown in 

Fig 16. The proton fluences among the three codes are in reasonable agreement. The larger descrepancy lies with the 

neutron spectra. One interesting note is that 3DHZETRN agrees with Geant4 below 10 MeV and agrees with 

FLUKA above 10 MeV. FLUKA and 3DHZETRN exhibit a similar inflection in the neutron spectrum near 10 to 60 

MeV that may be related to the use of an early FLUKA approximation to the secondary neutron spectrum used in 

HZETRN and carried over into 3DHZETRN. Some of the differences in proton spectra on the back surface may be 

associated with the simplified assumptions of the straight-ahead approximation in low energy isotropic light ion 

transport. A principle difference in the three codes is the computational efficiency, as discussed in more detail later 

in this section. It is also worth noting that the 3DHZETRN (N = 18) solution is a substantial improvement over the 

HZETRN (N = 1), indicating that lateral neutron leakage is more accurately described by the updated transport 

model. 

Fig. 16. Particle spectra at detector locations 35 g/cm2

 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2) exposed to Webber SPE spectrum. 

 The definition of the leakage factor given in equation (45) is specific to the 3DHZETRN formalism 

presented herein; however, an analogous quantity may also be computed for Monte Carlo results by comparing the 

ratio of the fluences in a slab (Fig. 8) to those obtained in a sphere (Fig. 16) provided the slab thickness and sphere 

diameters are identical. The leakage factor evaluated by three Monte Carlo codes (PHITS, FLUKA, and Geant4) are 

shown in comparison with 3DHZETRN (N = 18) at 35 g/cm2

 and 40 g/cm2 in Fig. 17. The leakage is related to the 

angular dependence of the cross section. For example, the near zero leakage in the 3DHZETRN results above 200 

MeV is a result of the straight-ahead approximation. It is clear that the Monte Carlo results show significant 
21

differences in leakage as determined mainly by differences in angular scattering cross sections. The 3DHZETRN 

results appear at a slight disadvantage with respect to Monte Carlo evaluations. It should also be clear that the 

improvement of the 3DHZETRN cross sections could result in more trustworthy results. This would not only 

improve leakage estimates but also give better representation of the resonance structure seen in the neutron fluences 

at lower energies. 

Fig. 17. The neutron leakage ratio factor at detector locations 35 g/cm2

 (left pane) and 40 g/cm2 (right pane) from 

top of aluminum sphere (diameter 40 g/cm2

) exposed to the Webber SPE spectrum. 

The time required by the three codes to run these benchmarks is given in CPU seconds used and listed in 

the first column of Table 5. The 5 CPU seconds of 3DHZETRN allows for the development of affordable means of 

design optimization software on a commercial laptop while Geant4, FLUKA, and PHITS software requires a 

multiprocessing cluster environment to do a single design evaluation and would still be impractical for design 

optimization. The main limitation of the Monte Carlo methods is the number of particle histories required for 

reasonable convergence. It should be noted that source biasing, physics biasing, or other optimization strategies 

could significantly reduce the Monte Carlo runtimes shown in Table 5. However, the fully optimized run times 

would still be orders of magnitude larger than the 3DHZETRN run times and still require a multi-processor cluster 

computing environment. 

Table 5. Total CPU seconds required for benchmark calculations. Total number of histories ran in Monte Carlo 

simulations are shown in parentheses. 

 Webber SPE GCR hydrogen GCR helium GCR carbon GCR iron 

3DHZETRN (N = 1) 2 2 2 2 2 

3DHZETRN (N = 18) 5 8 8 8 8 

Geant4 9×107

 (2×1011) 1×108 (2×1011) 3×108 (5×109) 2×108 (2×109) 4×108 (5×108) 

FLUKA 2×108

 (3×1011) 2×108 (3×1011) 1×108 (1×1010) 9×107 (7×109) 2×107 (5×108) 

PHITS 5×1010 (6×108

) 7×109 (1×108) 1×109 (3×108) 2×108 (1×108) 1×108 (3×108) 

 

Solar Minimum GCR benchmark. 

The present GCR benchmarks utilized the Badhwar-O’Neill 2010 [O’Neill 2010] model to represent the 1977 

solar minimum conditions (solar modulation parameter set to 475 MV). The hydrogen, helium, carbon, and iron 

spectra were considered separately as external source boundary conditions. The results for the neutron and proton 

fluence spectra within a 40 g/cm2

 sphere of aluminum for an incident solar minimum GCR hydrogen spectrum are 

shown in Fig. 18 at the two depths of 35 and 40 g/cm2

. The GCR hydrogen results are qualitatively similar to results 
22

obtained for the Webber SPE, only now the 3DHZETRN values favor Geant4 over the entire energy range except 

for the proton fluence below 100 MeV at the back surface (40 g/cm2

). The computation time used for transporting 

GCR hydrogen is slightly longer than for the Webber SPE transport in the same object. Still, the 3DHZETRN time 

is a very manageable 8 CPU seconds whereas the Monte Carlo codes take several orders of magnitude longer. 

Fig. 18. Particle spectra at detector locations 35 g/cm2

 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2

) exposed to the 1977 solar minimum GCR hydrogen spectrum. 

Fig. 19. Particle spectra at detector locations 35 g/cm2

 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2

) exposed to the 1977 solar minimum GCR helium spectrum. The proton and alpha results have 

been scaled by 0.1 and 0.01 to improve plot clarity. 

The neutron, proton, and alpha spectra generated by GCR helium ions are shown in Fig. 19. A large 

contribution to the alpha fluence spectrum is from the incident helium ions. There is generally better agreement with 

Geant4 results than those generated by FLUKA. The computer resources used in these calculations are given in 
23

Table 5. Similar results are shown for GCR carbon and iron ions in Figs. 20 and 21 with the corresponding resources 

used in Table 5. The light ion and neutron production nuclear database is more limited in HZETRN (and 

3DHZETRN) when the incident ion has Z > 2. First, the light ion and neutron fragments from the projectile are 

assumed to have the same velocity as the fragmenting ion itself. Second, the light ions and neutrons produced as 

target fragments have been ignored. Furthermore, the light ions are mostly treated within the straight-ahead 

approximation that mainly affects the light ion solution below 100 MeV/n. The spectral results will improve as these 

additional nuclear databases and corresponding transport procedures are added to the 3DHZETRN code. 

Fig. 20. Particle spectra at detector locations 35 g/cm2

 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2

) exposed to the 1977 solar minimum GCR carbon spectrum. The proton and alpha results have 

been scaled by 0.1 and 0.01 to improve plot clarity. 

Fig. 21. Particle spectra at detector locations 35 g/cm2

 (left pane) and 40 g/cm2 from top of aluminum sphere 

(diameter 40 g/cm2

) exposed to the 1977 solar minimum GCR iron spectrum. The proton and alpha results have 

been scaled by 0.1 and 0.01 to improve plot clarity. 
24

Although it is difficult to favor one code over another on the basis of the veracity of the solutions obtained, 

there are vast differences in computational time required to evaluate fluence spectra even in this simplest of 

geometries (homogeneous sphere of a single elemental material). Still, Monte Carlo methods provide an important 

test point for developing codes capable of meeting operational and especially design requirements. Computational 

speed has proven to be important in spaceflight validation in LEO where the time structure of the environment can 

be used to test various environmental components as was done using ISS and the Liulin instrument [Wilson et al. 

2007; Slaba et al. 2011a, 2013]. The main limitation on such studies remains the uncertainty in the environmental 

models (especially for the trapped environment) and uncertainty in nuclear cross sections. 

Conclusions 

The present version of 3DHZETRN shows improved agreement with Monte Carlo results in finite 

geometries where lateral leakage is non-negligible. Despite the improvements, the enhanced formalism and code has 

inherent limitations (aside from the simplified geometry of a homogeneous material) that carry over from the 

available bi-directional methods [Slaba et al. 2010b] in which the main share of light ions are still treated in the 

straight-ahead approximation as coupled to the forward propagating neutrons. Only the isotropically produced 

neutrons are first transported in 3D without the coupling to the isotropically produced light ion components. Once 

the isotropically produced neutrons are evaluated, the final coupling to the light ions is implemented (without further 

nuclear reactions) in a full 3D manner. Hence, only a fraction of the light ions have a 3D treatment, resulting in 

transport effects of unquantified magnitude. Future work will require extension to inhomogeneous media, complex 

geometries, and treatment of the cross-coupling of 3D neutron/light ion transport. Even so, the present expanded 

treatment of 3D effects shows great promise in the path to a fully 3D transport algorithm. 

The present benchmarks provide useful exercises for evaluation of computational procedures and 

atomic/nuclear database. It is now time for improving the physical description of the transport process that is still 

hampered by a lack of basic nuclear physics database developments. The codes for space design are still incomplete 

in both the physical description and corresponding implementation into computational procedures. Even so, the 

current development appears as state-of-the-art in computational procedures for spacecraft design. 

Acknowledgements 

This work was supported by the Human Research Program under the Human Exploration and Operations 

Mission Directorate of NASA and by NASA Grant number NNX09AR20A. 

Appendix A: Monte Carlo simulations 

 In this section, a brief description is given of the Monte Carlo codes Geant4, FLUKA, and PHITS used to 

simulate particle spectra shown in Figs. 8 and 16 - 21. A discussion of the source, geometry, and simulation setup 

used for the codes is also given. 

Geant4 version 9.4.6 [Agostinelli et al. 2003] was used in the present calculations with the physics list 

QGSP_BERT_HP. This physics list utilizes the quark gluon string model (QGS) for high energy interactions with 

nucleons, pions, and nuclei. Post-interaction nuclear de-excitation is handled by the precompound (P) model. The 

Bertini cascade (BERT) model is used for interactions below 10 GeV. A high-precision (HP) tracking model is used 

for neutrons below 20 MeV. Nucleus-nucleus collisions are represented by the native Quantum Molecular Dynamics 

(QMD) model [Koi 2008]. More information about the chosen physics list can be found at the Geant4 website 

[Geant4 2012a]. The QGSP_BERT_HP physics list has been recommended by the Geant4 collaboration for high 

energy physics applications [Geant4 2012b] and has been used, with some slight variation, in other space radiation 

related studies [Bernabeu and Casanova 2007, Hayatsu et al. 2008, Martinez and Kingston 2012, Slaba et al. 2013]. 

FLUKA is a Monte Carlo program that contains fully integrated physics and performs the transport of 

elementary particles and ions through materials. It has the ability to transport and calculate interactions with all 

elementary hadrons, light and heavy ions, and electrons and photons over an energy range which extends up to 104

TeV for all particles, and down to thermal energies for neutrons [Battistoni et al. 2007]. The nuclear models are 

integrated into the software for all particles and energies, with the exception of neutrons below 20 MeV, where 

standard international evaluated data files are used. This code has been extensively benchmarked against available 

accelerator and cosmic ray experimental data, at beam energies as low as a few MeV and as large as cosmic ray 

energies [Aarnio et al. 1993; Andersen et al. 2004]. 
25

PHITS has been used extensively in space radiation transport. It simulates ionization through transport 

processes and uses the mean free path to determine when collisions occur [Niita et al. 2006, Sihver et al. 2007]. 

Ionization transport can include angle and energy straggling, showing good agreement for Bragg peaks and 

fragmentation tails [Niita et al. 2006]. Tabulated nuclear cross section data are used for neutron interactions below 

20 MeV, and evaluated data are used for gamma and electron transport below 1 GeV [Niita et al. 2006]. High 

energy neutron and heavy charged particle interactions are simulated through the Jet AA Microscopic transport 

model (JAM), which addresses hadron-hadron interactions, and the Japan Atomic Energy Research Institute 

(JAERI) Quantum Molecular Dynamics model (JQMD), which addresses hadron-nucleus and nucleus-nucleus 

interactions, up to energies of 200 GeV/n [Niita et al. 2006, Sato et al. 2006, Sihver et al. 2007]. Both JAM and 

JQMD are used for the dynamic portion of the nuclear reaction, and once the dynamics are completed, the General 

Evaporation Model is used to address nuclear de-excitation through particle emission in a statistical manner [Niita et 

al. 2006, Sihver et al. 2007]. The total nucleus-nucleus cross sections used in PHITS are the same as those that have 

been used in versions of HZETRN [Niita et al. 2006]. 

In Fig. 8, the geometry setup was for a 40 g/cm2

 aluminum slab with infinite lateral dimensions and the 

external source was applied normally to the front face of the slab. Particle spectra were tallied at a thin (0.027 g/cm2) 

aluminum detector layer placed 35 g/cm2

 and 40 g/cm2 from the front of the slab. In Figs. 16 - 21, the geometry 

setup was for a sphere with diameter 40 g/cm2

, and the external source was applied uniformly onto the top of the 

sphere (covering the entire sphere) parallel with the z-axis as shown in Fig. 10. Particle spectra were tallied at 

spherical aluminum detectors (radius 0.675 g/cm2

) placed along the z-axis down through the sphere. Convergence 

tests on the detector radius were performed to arrive at the final choice. In general, the chosen detector radius is 

large enough to allow reasonable tally statistics, but small enough to effectively prevent fragmentation events from 

occurring within the detector. 

For all of the GCR simulations, the source energy spectrum was sampled directly from the input spectrum 

(i.e. no source energy biasing). For the SPE simulations, the source energy spectrum was biased so that source 

sampling was performed from a uniform distribution covering 0.01 MeV to 2500 MeV. This biasing technique 

places greater emphasis on the higher energy portion of the SPE spectrum improving tally statistics for secondary 

neutrons. 

Appendix B: Evaluation of i for N = 6, 10, 14, 18, 22 

The i are herein evaluated with 0 as the top hemispheric pole so that 1 = 0, while the bottom pole is at 

2 = -0 and are the central rays of the two polar regions. The remaining directions denoted by i for i = 3...N are 

symmetrically distributed about the mid-latitudes. Generally in application, the i’s must be rotated to the body 

coordinates of the transport media. This is accomplished by a rotation matrix R(0, 0) determined from the polar 

vector 0 relative to the body coordinate system (note,  denotes co-latitude with values on 0 to , and  denotes 

longitude running from 0 to 2). The regional boundaries are chosen such that 

4 / i

i d N p

Dò = W

W , (46) 

where Δi denotes the i

th region, N is the number of regions and is related to the number of N segments of cosi

and the number of N segments i plus the two polar regions such that 

N NN 2 = + f q , (47) 

and defined such that 

1 1

cos

N N

i ii

i i

f q

= =

å å D = DD W . (48) 

The polar caps are bound by 

1 cos 1 2 / q = - N , (49) 
26

2 cos 1 2 / q =- + N , (50) 

and the corresponding central rays 1 and 2 lie along the z axis. The middle latitudes are divided into NN

regions extending from the polar cap to the equator in each hemisphere. The corresponding longitudinal boundaries 

are set by the N such that the 

2 / i Nf D = f p . (51) 

Note that N is generally greater than or equal to 4. The remaining boundaries in latitude are given as 

1 2 cos (cos cos ) / i Nq D= - q qq . (52) 

For N = 6, there is only one midpoint in latitude taken as 

1 2 cos (cos cos ) / 2 0 i q qq =+ = , (53) 

for i  2. The four longitudinal points are given as 

f p pp = { } 0, / 2, , 3 / 2 . (54) 

The directional elements for each solid angle region are given by 

{ } cos sin , sin sin , cos i il il l W = fq fq q , (55) 

and are given for N = 6 in Table B.1. Similar results for other N used in this paper are given in Tables B.2 - B.5. 

Table B.1 Directional elements of the i for N = 6. 

i 1 2 3 4 5 6 

ex 0 0 1 -1 0 0 

ey 0 0 0 0 1 -1 

ez 1 -1 0 0 0 0 

Table B.2 Directional elements of the i for N = 10. 

i 1 2 3 4 5 6 7 8 9 10 

ex 0 0 0.9165 -0.9165 0 0 -0.9165 0.9165 0 0 

ey 0 0 0 0 0.9165 -0.9165 0 0 -0.9165 0.9165 

ez 1 -1 0.4 -0.4 0.4 -0.4 0.4 -0.4 0.4 -0.4 

Table B.3 Directional elements of the i for N = 14. 

i 1 3 5 7 9 11 13 

ex 0 0.8207 0 -0.8207 0 1 0 

ey 0 0 0.8207 0 -0.8207 0 1 

ez 1 0.5714 0.5714 0.5714 0.5714 0 0 

i 2 4 6 8 10 12 14 

ex 0 -0.8207 0 0.8207 0 -1 0 

ey 0 0 -0.8207 0 0.8207 0 -1 

ez -1 -0.5714 -0.5714 -0.5714 -0.5714 0 0 
27

Table B.4 Directional elements of the i for N = 18. 

i 1 3 5 7 9 11 13 15 17 

ex 0 0.7454 0 -0.7454 0 0.9750 0 -0.9750 0 

ey 0 0 0.7454 0 -0.7454 0 0.9750 0 -0.9750 

ez 1 0.6667 0.6667 0.6667 0.6667 0.2222 0.2222 0.2222 0.2222 

i 2 4 6 8 10 12 14 16 18 

ex 0 -0.7454 0 0.7454 0 -0.9750 0 0.9750 0 

ey 0 0 -0.7454 0 0.7454 0 -0.9750 0 0.9750 

ez -1 -0.6667 -0.6667 -0.6667 -0.6667 -0.2222 -0.2222 -0.2222 -0.2222 

Table B.5 Directional elements of the i for N = 22. 

i 1 3 5 7 9 11 13 15 17 19 21 

ex 0 0.6863 0 -0.6863 0 0.9315 0 -0.9315 0 1 0 

ey 0 0 0.6863 0 -0.6863 0 0.9315 0 -0.9315 0 1 

ez 1 0.7273 0.7273 0.7273 0.7273 0.3636 0.3636 0.3636 0.3636 0 0 

i 2 4 6 8 10 12 14 16 18 20 22 

ex 0 -0.6863 0 0.6863 0 -0.9315 0 0.9315 0 -1 0 

ey 0 0 -0.6863 0 0.6863 0 -0.9315 0 0.9315 0 -1 

ez -1 -0.7273 -0.7273 -0.7273 -0.7273 -0.3636 -0.3636 -0.3636 -0.3636 0 0 

References 

Aarnio, P.A., Fasso, A. Ferrari, A., Moehring, H.J., Ranft, J., Sala, P.R., Stevenson, G.R., Zazula, J.M., FLUKA: 

Hadronic benchmarks and applications. Proceedings of the MC93 International Conference on Monte Carlo 

Simulation in High Energy and Nuclear Physics, 1993. 

Adamczyk, A.M., Norman, R.B., Sriprisan, S.I., Townsend, L.W., Norbury, J.W., Blattnig, S.R., Slaba, T.C., 

NUCFRG3: Light ion improvements to the nuclear fragmentation model. Nucl. Instrum. & Methods A 678: 21-32; 

2012. 

Agostinelli, S. et al., Geant4--a simulation toolkit. Nucl. Instrum. & Methods A 506: 250-303; 2003. 

Alsmiller, R.G., Irving, D.C., Kinney, W.E., Moran, H.S., The validity of the straightahead approximation in space 

vehicle shielding shielding studies. In Second Symposium on Protection Against Radiations in Space, Arthur Reetz, 

ed. NASA SP-71, pp. 177-181, 1965. 

Andersen, V., Ballarini, F., Battistoni, G., Campanella, M., Carboni, M., Cerutti, F., Empl, A., Fasso, A., Ferrari, A., 

Gadioli, E., Garzelli, M.V., Lee, K., Ottolenghi, A., Pelliccioni, M., Pinsky, L.S., Ranft, J., Roesler, S., Sala, P.R., 

Wilson, T.L., The FLUKA code for space applications: recent developments. Adv. Space Res. 34: 1338–1346; 2004. 

Armstrong, T.W., Colburn, B.L., Predictions of secondary neutrons and their importance to radiation effects inside 

the International Space Station. Radiat. Meas. 33: 229-234; 2001. 

Badhwar, G.D., Patel, J.U., Cuccinotta, F.A., Wilson, J.W., Measurements of the secondary particle energy spectra 

in the space shuttle, Rad. Meas., 24: 129–138; 1995. 

Badhwar, G.D., O’Neill, P.M., Cucinotta, F.A., In-flight radiation measurements on STS-60. Radiat. Meas. 26:17-

34; 1996. 

Badhwar, G.D., Atwell, W., Badavi, F.F., Yang, T.C., Space radiation absorbed dose distribution in human 

phantom. Radiat. Meas. 33: 235-241; 2001. 
28

Battistoni, G., Muraro, S., Sala, P.R., Cerutti, F., Ferrari, A., Roesler, S.,Fasso, A., Ranft, J., The FLUKA code: 

Description and benchmarking. Proceedings of the Hadronic Shower Simulation Workshop 2006, 896: 31–49; 2007. 

Bernabeu, J. and Casanova, I., Geant4-based radiation hazard assessment for human exploration missions. Adv. 

Space Res. 40: 1368-1380; 2007. 

Clowdsley, M.S., Heinbockel, J.H., Kaneko, H., Wilson, J.W., Singleterry, R.C., Shinn, J.L., A comparison of the 

multigroup and collocation methods for solving the low-energy neutron Boltzmann equation. Can. J. Phys. 78: 45-

56; 2000. 

Clowdsley, M.S., Wilson, J.W., Shinn, J.L., Badavi, F.F., Heinbockel, J.H., Atwell, W., Neutron environment 

calculations for low Earth orbit. SAE 01ICES2327, 2001. 

Cucinotta, F.A. and Dubey, R.R., Alpha-cluster description of excitation energies in 12C(12C,3α)X at 2.1A GeV, 

Phys. Rev. C 50: 979–984; 1994. 

Cucinotta, F.A., Wilson, J.W., Shinn, J.L., Tripathi, R.K., Assessment and requirements of nuclear reaction data 

bases for GCR transport in the atmosphere and structures, Adv. Space Res., 21(12): 1753-1762; 1998. 

Cucinotta, F.A., Wilson, J.W., Saganti, P., Hu, X. Kim, M.Y., Cleghorn, T., Zeitlin, C.E., Tripathi, R.K., Isotopic 

dependence of GCR fluence behind shielding. Radiat. Meas. 41(9):1235-1249; 2007. 

Fasso, A., Ferrari, A., Ranft, J., Sala, P.R., FLUKA: A multi-particle transport code, CERN-2005-10, INFN/TC 

05/11, SLAC-R-773, 2005. 

Foelsche, T., Mendell, R.B., Wilson, J.W., Adams, R.R., Measured and calculated neutron spectra and dose 

equivalent rates at high altitudes: Relavance to SST operations and space research. NASA TN D-7715, 1974. 

Geant4 collaboration. Geant4 Reference Physics Lists. Website: 

http://www.geant4.org/geant4/support/physicsLists/referencePL/referencePL.shtml (2012a). 

Geant4 collaboration. Geant4 Recommended Physics Lists. Website: 

http://www.geant4.org/geant4/support/physicsLists/referencePL/useCases.shtml (2012b). 

Golovchenko, A.N., Skvarč, J., Ilić, R., Sihver, L., Bamblevski, V.P., Tretyakova, S.P., Schardt, D., Tripathi, R.K., 

Wilson, J.W., R. Bimbot, R., Fragmentation of 200 and 244 MeV/u carbon beams in thick tissue-like absorbers. 

Nucl. Instrum. & Methods B, 159: 233-240; 1999. 

Hayatsu, K., Hareyama, M., Kobayashi, S., Yamashita, N., Miyajima, M., Sakurai, K., Hasebe, N., Radiation doses 

for humans exposed to galactic cosmic rays and their secondary products on the lunar surface. Bio. Sci. in Space 22: 

59-66; 2008. 

Heinbockel, J.H., Feldman, G.A., Wilson, J.W., Singleterry, R.C., Leakeas, C.L., Clowdsley, M.S., Solutions to the 

low energy neutron Boltzmann equation for space applications. SAE 03ICES339, 2003. 

Heinbockel, J.H., Slaba, T.C., Blattnig, S.R., Tripathi, R.K., Townsend, L.W., Handler, T., Gabriel, T.A., Pinsky, 

L.S., Reddell, B., Clowdsley, M.S., Singleterry, R.C., Norbury, J.W., Badavi, F.F., Aghara, S.K., Comparison of the 

transport codes HZETRN, HETC and FLUKA for a solar particle event. Adv. Space Res. 47: 1079-1088; 2011a. 

Heinbockel, J.H., Slaba, T.C., Tripathi, R.K., Blattnig, S.R., Norbury, J.W., Badavi, F.F., Townsend, L.W., Handler, 

T., Gabriel, T.A., Pinsky, L.S., Reddell, B., Aumann, A.R., Comparison of the transport codes HZETRN, HETC and 

FLUKA for galactic cosmic rays. Adv. Space Res. 47:1089-1105; 2011b. 

Hugger, C.P., Qualls, G.D., Wilson, J.W., Cucinotta, F.A., Shavers, M.R., Zapp, N., ISS as a platform for 

environmental evaluation. AIP Conf. Proc. 654: 1011-1017; 2003. 
29

ICRP, Recommendations of the ICRP. ICRP Publication 26, 1977. 

Koi, T., New native QMD code in Geant4. IEEE Nuclear Science Symposium Conference Record, 2008. 

Lamkin, S.L., A theory for high-energy nucleon transport in one-dimension. Masters Thesis, Old Dominion 

University, June 1974. 

Nealy, J.E., Cucinotta, F.A., Wilson, J.W., Badavi, F.F., Dachev, Ts.P., Tomov, B.T., Walker, S.A., de Angelis, G., 

Blattnig, S.R., Atwell, W., Pre-engineering spaceflight validation of environmental models and the 2005 HZETRN 

simulation code. Adv. Space Res. 40(11): 1592-1610; 2006. 

Martinez, L.M. and Kingston, J., Space radiation analysis: Radiation effects and particle interaction outside the 

earth's magnetosphere using GRAS and GEANT4. Acta Astronautica 72: 156-164; 2012. 

Niita, K., Sato, T., Iwase, H., Nose, H., Nakashima, H., Sihver, L., 2006 PHITS - A particle and heavy ion transport 

code system, Rad. Meas., 41: 1080 – 1090; 2006. 

Marchuk, G.I. and Lebedev, V.I., Numerical Methods in the Theory of Neutron Transport, Harwood Academic 

Publishers, 1986. 

Norbury, J.W., Adamczyk, A., Differential cross sections for electromagnetic dissociation. Nucl. Instrum. & 

Methods B 254: 177-181; 2007. 

Norman, R.B., Blattnig, S.R., De Angelis, G., Badavi, F.F., Norbury, J.W., Deterministic pion and muon transport in 

Earth’s atmosphere. Adv. Space Res. 50: 146-155; 2012. 

Norman, R.B., Slaba, T.C., Blattnig, S.R., An extension of HZETRN for cosmic ray initiated electromagnetic 

cascades. Adv. Space Res. 51: 2251-2260; 2013. 

O'Neill, P.M., Badhwar-O'Neill galactic cosmic ray flux model – Revised. IEEE Trans. Nuc. Sci., 57: 3148–3153; 

2010. 

Pinsky, L., Carminati, F., Ferrari, A., Simulation of space shuttle neutron measurements with FLUKA. Radiat. 

Meas. 33: 335-339; 2001. 

Ranft, J., The FLUKA and KASPRO hadronic cascade codes. Computer Techniques in Radiation Transport and 

Dosimetry, W. R. Nelson and T.M. Jenkins, eds, Plenum Press, pp. 339-371, 1980 

Sato, T., Niita, K., Iwase, H., Nakashima, H., Yamaguchi, Y., Sihver, L., 2006 Applicability of particle and heavy 

ion transport code PHITS to the shielding design of spacecrafts, Rad. Meas., 41: 1142 – 1146; 2006. 

Sato, T., Niita, K., Matsuda, N., Hashimoto, S., Iwamoto, Y., Noda, S., Ogawa, T., Iwase, H.,Nakashima, H., 

Fukahori, T., Okumura, K., Kai, T., Chiba, S., Furuta T., Sihver, L., Particle and heavy ion transport code system 

PHITS, version 2.52, J. Nucl. Sci. Technol. 50: 913-923; 2013. 

Schimmerling, W., Rapkin, M., Wong, M., Howard, J., The propogation of relativistic heavy ions in multielement 

beam lines. Medical Phys. 13: 217-228; 1986. 

Scott, W.W., Alsmiller, R.G., Comparisons of results obtained with several proton penetration codes. ORNL-RSIC17. 1967. 

Shavers, M.R., Frankel, K., Miller, J., Schimmerling, W., Townsend, L.W., Wilson, J.W., The fragmentation of 

670 A MeV Neon-20 as a function of depth in water. III. Analytic Multigeneration Transport Theory. Radiat. Res.,

134(1): 1-14; 1993. 
30

Shinn, J.L., Cucinotta, F.A., Simonsen, L.C., Wilson, J.W., Badavi, F.F., Badhwar, G.D., Miller, J., Zeitlin, C., 

Heilbronn, L., Tripathi, R.K., Clowdsley, M.S., Heinbockel, J.H., Xapsos, M.A., Validation of a comprehensive 

space radiation transport code, IEEE Trans. Nuc. Sci., 45: 2711-2719; 1998. 

Sihver, L., Mancusi, D., Sato, T., Niita, K., Iwase, H., Iwamoto, Y., Matsuda, N., Nakashima, H., Sakamoto, Y., 

2007 Recent developments and benchmarking of the PHITS Code. Adv. Space Res., 40: 1320 – 1331; 2007. 

Singleterry, R.C., Blattnig, S.R., Clowdsley, M.S., Qualls, G.D., Sandridge, C.A., Simonsen, L.C., Slaba, T.C., 

Walker, S.A., Badavi, F.F., Spangler, J.L., Aumann, A.R., Zapp, E.N., Rutledge, R.D., Lee, K.T., Norman, R.B., 

Norbury, J.W., OLTARIS: On-line tool for the assessment of radiation in space. Acta Astronautica 68: 1086-1097; 

2011. 

Slaba, T.C., Blattnig, S.R., Badavi, F.F., Faster and more accurate transport procedures for HZETRN. J. Comp. 

Phys. 229: 9397-9417; 2010a. 

Slaba, T.C., Blattnig, S.R., Aghara, S.K., Townsend, L.W., Handler, T., Gabriel, T.A., Pinsky, L.S., Reddell, B., 

Coupled neutron transport for HZETRN. Radiat. Meas. 45: 173-182; 2010b. 

Slaba, T.C., Blattnig, S.R., Badavi, F.F., Stoffle, N.N., Rutledge, R.D., Lee, K.T., Zapp, E.N., Dachev, Ts. P., 

Tomov, B.T., Statistical validation of HZETRN as a function of vertical cutoff rigidity using ISS measurements. 

Adv. Space Res. 47: 600-610; 2011a. 

Slaba, T.C., Blattnig, S.R., Clowdsley, M.S., Variations in lunar neutron dose estimates. Rad. Res. 176: 827-841; 

2011b. 

Slaba, T.C., Blattnig, S.R., Reddell, B., Bahadori, A., Norman, R.B., Badavi, F.F., Pion and electromagnetic 

contribution to dose: Comparisons of HZETRN to Monte Carlo results and ISS data. Adv. Space Res. 52: 62-78; 

2013. 

Tai, H., Bichsel, H., Wilson, J.W., Shinn, J.L., Cucinotta, F.A., Badavi, F.F., Comparison of stopping power and 

range data bases for radiation transport study, NASA TP-3644, 1997. 

Tripathi, R.K., Townsend, L.W., Kahn, F., Role of intrinsic width in fragment momentum distributions in heavy ion 

collisions, Phys. Rev. C. 48(4): R1775-R1777; 1994. 

Tweed, J., Walker, S.A., Wilson, J.W., Cucinotta, F.A., Tripathi, R.K., Blattnig, S., Mertens, C.J., Computational 

methods for the HZETRN code. Adv. Space Res. 35: 194-201; 2006a. 

Tweed, J., Walker, S.A., Wilson, J.W., Tripathi, R.K., Cucinotta, F.A., Badavi, F.F., An improved Green’s function 

code for HZE transport. SAE/ICES paper 2006-01-2147, 2006b. 

Walker, S.A., Tweed, J., Wilson, J.W., Cucinotta, F.A., Tripathi, R.K., Blattnig, S., Zeitlin, C., Heilbronn, L., 

Miller, J., Validation of the HZETRN code for laboratory exposures with 1 A GeV iron ions in several targets. Adv. 

Space Res. 35: 202-207; 2006. 

Webber, W.R., An evaluation of solar-cosmic-ray events during solar minimum. D2-84274-1, Boeing Co. 1966. 

Wilson, J.W., Khandelwal, G.S., Proton dose approximation in convex geometry. Nucl. Tech. 23: 298-305; 1974. 

Wilson, J.W. and Lamkin, S.L, Perturbation theory for charged-particle transport in one dimension. Nucl. Sci. & 

Eng. 57(4): 292-299; 1975. 

Wilson, J.W., Analysis of the theory of high-energy ion transport. NASA TN D-8381, 1977. 

Wilson, J.W., Townsend, L.W., Bidasaria, H.B., Schimmerling, W., Wong, M. Howard, J., 20-Ne depth-dose 

relations in water. Health Phys. 46: 1101-1111; 1984. 
31

Wilson, J.W. and Badavi, F.F., Methods of galactic heavy ion transport. Radiat. Res. 108: 231-237; 1986. 

Wilson, J.W., Townsend, L.W., Badavi, F.F., A semiempirical nuclear fragmentation model. Nucl. Instrum. & 

Methods B 18: 225-231; 1987a. 

Wilson, J.W., Townsend, L.W., Badavi, F.F., Galactic cosmic ray propagation in earth's atmosphere. Radiat. Res.

109: 173-183; 1987b. 

Wilson, J.W., Chun, S. Y., Buck, W.W., Townsend, L.W., High energy nucleon data bases. Health Phys., 55: 817-

819; 1988. 

Wilson, J.W., Townsend, L.W., Schimmerling, W., Khandelwal, G.S., Khan, F., Nealy, J.E., Cucinotta, F.A., 

Simonsen, L.C., Shinn, J.L., Norbury, J.W., Transport methods and interactions for space radiations, NASA RP1257, 1991. 

Wilson, J. W., Shavers, M.R., Badavi, F.F., Miller, J., Shinn, J.L., Costen, R.C., Nonperturbative methods in hze 

propagazition. Radiat. Res. 140: 241-244; 1994a. 

Wilson, J.W., Townsend, L.W., Shinn, J.L., Badavi, F.F., Lamkin, S.L., Galactic cosmic ray transport methods: 

Past, present, and future. Adv. Space Res. 14: 841-852, 1994b. 

Wilson, J.W., Shinn, J.L., Townsend, L.W., Tripathi, R.K., Badavi, F.F., Chun, S.Y., NUCFRG2: A semiempirical 

nuclear fragmentation model. Nucl. Instrum. & Methods B 94: 95-102; 1994c. 

Wilson, J.W., Tripathi, R. K., Cucinotta , F. A., Shinn, J. L., Badavi, F. F., Chun, S. Y., Norbury, J. W., Zeitlin, 

C. J., Hielbronn, L. H., Miller, J., NUCFRG2, an evaluation of the semiempirical nuclear fragmentation database, 

NASA TP-3533, 1995. 

Wilson, J.W, Miller, J.M., Konradi, A., Cucinotta, F.A., Shielding strategies for human space exploration. NASA 

CP-3360, 1997. 

Wilson, J.W., Cucinotta, F.A., Tai, H., Shinn, J.L., Chun, S.Y., Tripathi, R.K., Sihver, L., Transport of light ions in 

matter. Adv. Space Res. 21(12): 1763-1771; 1998. 

Wilson, J.W., Tweed, J., Heinbockel, J.H., Tripathi, R.K., Cucinotta, F.A., Mertens, C.J., Advances in space 

radiation transport codes. 14th Space Radiation Health Investigators' Workshop, League City, TX; 2003a. 

Wilson, J.W., Nealy, J.E., De Angelis, G., Badavi, F.F., Hugger, C.P., Cucinotta, F.A., Kim, M.Y., 

Dynamic/anisotropic low earth orbit environmental models. AIAA Space Conference paper 2003-6221, 2003b. 

Wilson, J.W., Tripathi, R.K., Qualls, G.D., Cucinotta, F.A., Prael, R.E., Norbury, J.W., Heinbockel, J.H., Tweed, J., 

Space radiation transport methods development, Adv. Space Res. 34(6): 1319-27; 2004a. 

Wilson, J.W., Cucinotta, F.C., Schimmerling, W.S., Emerging radiation health-risk mitigation technologies. AIP 

Conf. Proc. 699, 913-924; 2004b. 

Wilson, J.W., Tripathi, R.K., Mertens, C.J., Blattnig, S.R., Clowdsley, M.S., Verification and validation: High 

charge and energy (HZE) transport codes and future development. NASA TP-2005-213784, 2005. 

Wilson, J.W., Tripathi, R.K., Badavi, F.F., Cucinotta, F.A., Standardized radiation shield design method: 2005 

HZETRN, SAE/ICES paper 2006-01-2109, 2006. 

Wilson, J.W., Nealy, J.E., Dachev, Ts.P., Tomov, B.T., Cucinotta, F.A., Badavi, F.F., De Angelis, G., Atwell, W., 

Leutke, N., Time serial analysis of the induced LEO environment within the ISS 6A. Adv. Space Res. 40(11): 1562-

1570; 2007. 
32

Zeitlin, C., Heilbronn, L., Miller, J., Rademacher, S.E., Borak, T., Carter, T.R., Frankel, K.A., Schimmerling, W., 

Stronach, C.E., Heavy fragment production cross sections from 1.05 GeV/nucleon 56Fe in C, Al, Cu, Pb, and CH2

targets. Phys. Rev. C 56: 388-397; 1997. 

Zeitlin, C., Miller, J., Guetersloh, S., Heilbronn, L., Fukumura, A., Iwata, Y., Murakami, T., Blattnig, S., Norman, 

R.. Mashnik, S., Fragmentation of 14N, 16O, 20Ne, and 24Mg nuclei at 290 to 1000 MeV/Nucleon. Phys. Rev. C 83: 

034909; 2011. 
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Publication

 4. TITLE AND SUBTITLE

A 3DHZETRN Code in a Spherical Uniform Sphere with Monte Carlo 

Verification

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Wilson, John W.; Slaba, Tony C.; Badavi, Francis F.; Reddell, Brandon; 

Bahadori, Amir A.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20415

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 93

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

The computationally efficient HZETRN code has been used in recent trade studies for lunar and Martian exploration and is 

currently being used in the engineering development of the next generation of space vehicles, habitats, and extra vehicular 

activity equipment. A new version (3DHZETRN) capable of transporting High charge (Z) and Energy (HZE) and light ions 

(including neutrons) under space-like boundary conditions with enhanced neutron and light ion propagation is under 

development. In the present report, new algorithms for light ion and neutron propagation with well-defined convergence 

criteria in 3D objects is developed and tested against Monte Carlo simulations to verify the solution methodology. The code 

will be available through the software system, OLTARIS, for shield design and validation and provides a basis for personal 

computer software capable of space shield analysis and optimization. 

15. SUBJECT TERMS

HZETRN; Particle transport; Radiation

18. NUMBER

 OF 

 PAGES

40

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

651549.02.07.10 

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TP-2014-218271

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 05 - 2014