February 2013 

NASA/TP–2013-217803 

Faster Heavy Ion Transport for HZETRN 

Tony C. Slaba 

Langley Research Center, Hampton, Virginia 
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 

 

February 2013 

NASA/TP–2013-217803 

Faster Heavy Ion Transport for HZETRN 

Tony C. Slaba 

Langley Research Center, Hampton, Virginia 
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320

443-757-5802 
iii

Contents

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

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

2. Heavy Ion Transport Formalism................................................................................................................. 1

3. Heavy Ion Marching Equation ................................................................................................................... 3

4. New Numerical Method ............................................................................................................................. 6

5. Flux Comparison ........................................................................................................................................ 6

6. Convergence Test ....................................................................................................................................... 8

7. Run-time Comparisons..............................................................................................................................10

8. Conclusions...............................................................................................................................................10

9. References.................................................................................................................................................11
iv

Figures

1. The left pane shows the partial fragmentation cross section for 56Fe + Al as a function of proton range 

and kinetic energy. The right pane shows the scaled ion fluxes at 20 g/cm2 in aluminum exposed to 

solar minimum GCR conditions given by the Badhwar-O'Neill 2010 model [O'Neill 2010] with a solar 

modulation parameter of 475 MV. .......................................................................................................... 5

2. Relative difference between fluxes computed by the original and modified heavy ion transport 

algorithm as a function of kinetic energy in aluminum and water exposed to solar minimum GCR 

conditions. Light and dark error bars represent extreme and standard deviation of relative differences, 

respectively.............................................................................................................................................. 7

3. Relative difference between fluxes computed by the original and modified heavy ion transport 

algorithm as a function of shielding thickness in aluminum and water exposed to solar minimum GCR 

conditions. Light and dark error bars represent extreme and standard deviation of relative differences, 

respectively.............................................................................................................................................. 8

4. Relative difference between fluxes computed by the original and modified heavy ion transport 

algorithm as a function of charge in aluminum and water exposed to solar minimum GCR conditions. 

Light and dark error bars represent extreme and standard deviation of relative differences, respectively.

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

5. Dose equivalent versus depth in aluminum (left pane) and water (right pane) exposed to solar 

minimum GCR conditions computed with the original and modified heavy ion transport algorithms 

and all discretization parameters. ............................................................................................................ 9

6. Total discretization error for dose equivalent as a function of depth in aluminum (left pane) and water 

(right pane) exposed to a solar minimum GCR environment.................................................................. 9
v

Tables

1. Run-time (seconds) comparisons between original and new heavy ion transport algorithms. ...............10
1

Abstract

The deterministic particle transport code HZETRN was developed to enable fast and 

accurate space radiation transport through materials. As more complex transport 

solutions are implemented for neutrons, light ions (Z < 2), mesons, and leptons, it is 

important to maintain overall computational efficiency. In this work, the heavy ion (Z > 

2) transport algorithm in HZETRN is reviewed, and a simple modification is shown to 

provide an approximate 5x decrease in execution time for galactic cosmic ray transport. 

Convergence tests and other comparisons are carried out to verify that numerical 

accuracy is maintained in the new algorithm.

1. Introduction

The deterministic particle transport code, HZETRN (High charge (Z) and Energy TRaNsport)

[Wilson et al. 1991; Slaba et al. 2010a, 2010b], is a fast and accurate computational tool for transporting 

space radiation through shielding materials and human tissue. Recent updates to the code have provided 

significant reductions in run time [Slaba et al. 2010a], and further work reduced overall discretization error 

to less than 4% for shielding materials and thicknesses commonly found in space applications [Slaba et al. 

2012a]. This increased efficiency has enabled new calculations that would have otherwise been 

computationally prohibitive. First, a ray-by-ray transport methodology has been developed wherein 

transport is explicitly performed along each ray direction in a shielding thickness distribution, and no 

constraint is placed on the number or ordering of materials along a ray [Slaba et al. 2011a]. Second, the 

transport of pions, muons, electrons, positrons, and photons has been added into HZETRN at a moderate 

computational cost [Norman et al. 2012a, 2012b; Slaba et al. 2012b]. Third, calculations may be performed 

through detailed vehicle geometry on a minute-by-minute basis along a trajectory. This type of analysis has 

been useful in providing more rigorous validation comparisons and uncertainty estimates [Slaba et al. 

2011b, 2012b]. 

Further developments that lead to more efficient numerical methods or algorithms are beneficial 

and may also enable new analyses to be completed. In this work, the heavy ion (Z > 2) transport algorithm 

in HZETRN is reviewed, and a simple modification is made that allows interpolation to be avoided in the 

ion collision source term. It is shown that the modification provides an approximate 5x decrease in 

execution time without sacrificing numerical accuracy. 

2. Heavy Ion Transport Formalism

In this section, the heavy ion transport formalism in HZETRN is reviewed. The Boltzmann 

transport equation with the continuous slowing down and straight ahead approximations is given as 

[Wilson et al. 1991] 

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

B xE   E E x E dE   , (1)

with the linear differential operator

1 [ ( , )] () () (, ) j j j j

j

B xE SE E x E x AE   

        

   

, (2)

and boundary condition

(0, ) ( ) j j E fE  . (3)
2

In equations (1) - (3), (, ) j  x E is the flux of type j particles at depth x with kinetic energy E, Aj is the 

atomic mass number of a type j particle, Sj(E) is the stopping power of a type j ion with kinetic energy E, 

σj(E) is the total macroscopic cross section for a type j particle with kinetic energy E, and σjk(E,E') is the 

differential macroscopic production cross section for interactions in which a type k particle with kinetic 

energy E' produces a type j particle with kinetic energy E. The summation limits in equation (1) are taken 

over all projectiles that might produce a given fragment. The boundary condition spectrum, ( ) j f E , is 

considered to be a known function over a broad energy spectrum. 

The ion stopping powers, Sj(E), in equation (2) may be approximated in terms of the proton 

stopping power, S(E), as [Bethe et al. 1930]

1 () () j j

j

S E SE

A   , (4)

where the scaling parameter is expressed in terms of the ion mass and charge (Aj, Zj) according to

2 / j jj   Z A . The transport equation may now be written as 

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

B xE   E E x E dE 

  , (5)

with the modified linear differential operator

[ ( , )] ( ) ( ) (, ) B xE j j jj SE E xE

x E   

   

   

    . (6)

Approximating the ion stopping powers with the scaled proton stopping power allows simplified numerical 

algorithms, as has been discussed elsewhere [Wilson et al. 2006].

For heavy ions, it is noted that projectile fragments have an energy and direction very near that of 

the projectile, while target fragments are produced nearly isotropically with low energy and travel only a 

short distance before being absorbed [Wilson et al. 2006]. The approximate decoupling of target and 

projectile fragments is discussed in detail by Wilson et al. [1995] and suggests that target fragments can be 

neglected in the heavy ion transport procedure (their contribution to dose is approximately accounted for 

after the transport procedure). The equal velocity assumption for heavy ions can be expressed in the 

production cross section as [Shinn et al. 1992]

( , ') ( ') ( ') jk EE E E E jk     , (7)

where σjk(E) is the partial fragmentation cross section for interactions in which a type k particle with kinetic 

energy E produces a type j particle. In this case, equation (5) becomes

[ ( , )] ( ) ( , ) j jk k

k

B xE E xE     . (8)

The absence of target fragments in the heavy ion transport procedure allows one to take the summation in 

equation (5) over all projectiles with mass greater than that of the fragment. If all the transported particles 

are ordered according to mass, then equation (5) can be succinctly written as 

[ ( , )] ( ) ( , ) j jk k

k j

B xE E xE  



  , (9)
3

which is the transport equation found in Wilson et al. [2006] and will be referred to as the heavy ion 

transport equation. The upper summation limit in equation (9) can vary, but it is common to use no fewer 

than 59 ions [Cucinotta et al. 2006]. 

For light particles (n, p, 2

H, 2H, 3He, 4He), both projectile and target fragments are included in the 

transport procedure. The broad energy distribution in collision events also indicates that the equal velocity 

assumption in equation (7) cannot be used. In this case, no simplifications to equation (5) are available, and 

the summation is taken over all light particles. It should be noted that for solar particle events with a 

negligible heavy ion component, only the solution to the light particle transport equation is required. For 

galactic cosmic ray (GCR) environments, there are both heavy ion and light ion components, and solutions 

to the light particle and heavy ion transport equations must be evaluated simultaneously. The light particle 

transport solutions and coupling to heavy ion transport has been discussed in detail elsewhere [Slaba et al. 

2010a] and will not be repeated here.

3. Heavy Ion Marching Equation

In this section, the mathematical formalism that allows equation (9) to be inverted and expressed 

as an explicit marching procedure is reviewed. The following derivation was taken from Wilson et al. 

[1977, 1986, 1991, 1995, 2006] and Shinn et al. [1992]. The scaled heavy ion flux is defined as

(, ) ( ) (, ) j jj   xr SE xE  , (10)

where r is the proton range defined in terms of the stopping power according to the equation

0

'

( ')

E dE r S E   . (11)

Equation (9) may be written in terms of the scaled flux and proton range as 

() (, ) () (, ) j

j j j jk k

k j k

r xr r xr

x r



   

 

       

     , (12)

where the quantities, ( ) j E and ( ) jk E have been replaced with ( ) j  r and ( ) jk  r since, for a given r, 

equation (11) may be inverted to find E. Equation (12) may be inverted using the method of characteristics

to obtain the Volterra integral equation 

     

(, )

, '

0

( , ) (0, )

' ', ' ' ,

j

j

r x

j j j

x j r x

jk j j k

k j k

xr e r x

e r x x xr x dx





 



   









 

   

(13)

with the exponential term

0 (, ) ( ) x

j jj   r x r t dt    . (14)

The complete details of how equation (9) is transformed into equation (12) as well as the inversion of 

equation (12) have been given by Slaba et al. [2010c]. Equation (13) may be written as a marching 

procedure in terms of a step-size, h, as 
4

     

(, )

, '

0

( , ) (, )

' ', ' ' .

j

j

r h

j j j

h j r x

jk j j k

k j k

x hr e xr h

e r x x h xr x dx





  



   









 

    

(15)

The integrand in equation (15) is approximated by considering the solution to the homogeneous form of 

equation (12). The homogeneous equation neglects secondary particle production through nuclear

interactions and accounts only for the slowing down of particles due to atomic interactions and the loss of 

particles due to nuclear absorption. If the step-size is taken to be sufficiently small such that 

1

( ) j

h  r , (16)

(i.e. much less than the nuclear mean free path), then the local truncation error will be negligible as the 

particles will not have travelled far enough to participate in a nuclear collision [Wilson et al. 1991]. The

homogeneous solution is given by the equation [Wilson et al. 2006]

(, ) ( , ) (, ) j r h

j j j x hr e xr h        , (17)

which may be equivalently written as 

( ', ') ( ', ') ( , ' ( ')) k j r xh x

k j k jk x h xr x e xr x h x               , (18)

and substituted into equation (15) to obtain

   

(, )

, ' ( ', ')

0

( , ) (, )

' ( , ' ( ')) ' .

j

j kj

r h

j j j

h j r x r xh x

jk j k j k

k j k

x hr e xr h

e r x x r x h x dx



 

  



   





  



 

    

(19)

The fragmentation cross section and particle flux appearing in the integrand of equation (19) are 

approximated as

  jk rx rh     j j '  jk  / 2, (20)

( , ' ( ')) ( , ( ) / 2) k j k k jk      xr x h x xr h      . (21)

Equation (20) is accurate since heavy ion fragmentation cross sections vary slowly as a function of energy 

(range), and most collisions occur at energies for which the ion range is large compared to typical stepsizes. Equation (21) is accurate since the scaled heavy ion fluxes are nearly constant at low range values. At 

higher range values, the step-size is comparatively small, and the step-size perturbation has a small affect. 

These statements are supported in Figure 1, showing representative fragmentation cross sections and scaled

ion fluxes as a function of proton range (bottom axis) and kinetic energy (top axis).
5

Figure 1. The left pane shows the partial fragmentation cross section for 56Fe + Al as a function of proton 

range and kinetic energy. The right pane shows the scaled ion fluxes at 20 g/cm2 in aluminum exposed to 

solar minimum GCR conditions given by the Badhwar-O'Neill 2010 model [O'Neill 2010] with a solar 

modulation parameter of 475 MV.

Upon substitution of equations (20) and (21), the heavy ion transport equation becomes

   

(, )

, ' ( ', ')

0

( , ) (, )

/ 2 ( , ( ) / 2) '.

j

j kj

r h

j j j

h j r x r xh x

jk j k j k

k j k

x hr e xr h

r h x r h e dx



 

  



   





  



 

     (22)

The remaining integral is evaluated by expanding the exponential terms in the integrand of equation (22). 

The complete details of the expansion and integral evaluation have been given by Slaba et al. [2010c]. The 

integral is approximately evaluated as 

  ( /2) ( ( ) /2)

, ' ( ', ')

0 ' ( ( ) / 2) ( / 2)

( ).

j j k jk

j kj

r hh r hh h r x r xh x

k jk j j

jk

e e e dx

r h rh

C r

     

   

   

       



 (23)

The final marching procedure is now given by the equation

 

(, ) ( , ) (, )

/ 2 ( , ( ) / 2) ( ),

j r h

j j j

j

jk j k j k jk

k j k

x hr e xr h

r h xr h C r

   



   







 

   

(24)

This is the heavy ion marching equation given by Wilson et al. [2006]. The heavy ion fluxes are evaluated 

on a pre-defined energy grid that is converted to proton range values. There are typically 100 grid points 

spaced logarithmically from 0.01 MeV/n to 50 GeV/n. In order to obtain the ion flux at a grid point, ri, 

interpolation is needed to evaluate the terms (, ) j j   xr h  and ( , ( ) / 2) k jk   xr h   . The 

interpolation method that has been typically used in HZETRN is log-log cubic Lagrange [Wilson et al. 

1995]. In this method, the logarithm of both the dependent and independent variables is evaluated. A cubic 

Lagrange interpolation scheme is applied, and the resulting value is exponentiated. 
6

4. New Numerical Method

In this section, the original numerical implementation of the heavy ion marching equation is 

modified to provide a more efficient algorithm. As stated in the previous section, interpolation is required

to evaluate the ion fluxes at values not on the prescribed energy/range grid. By examining the summation

term in equation (24), it can be determined that the number of interpolations per marching step scales 

linearly with the number of energy grid points and quadratically with the number of isotopes. For typical 

calculations, this number is on the order of 103 interpolations per step. For a given range value, ri, an 

alternative strategy to evaluate equation (24) is to use the nearest-neighbor approximation

* ( , ( ) / 2) ( , ) ki j k ki    x r   h xr , (25)

where *

i r is the nearest grid point value to ( ) /2 i jk r h     . For each ri, the heavy ion marching 

equation becomes 

 

(,)

*

( , ) (, )

/ 2 ( , ) ( ).

j i r h

j i ji j

j

jk i j k i jk i

k j k

x hr e xr h

r h xr C r

   



  







 

 

(26)

In this equation, no interpolation is needed to evaluate the ion flux within the production 

summation. Consequently, the only interpolation required is for (, ) ji j   xr h  , and the total number of 

interpolations per step is now directly proportional the number of energy grid points. In order to take 

advantage of the nearest neighbor approach, the mapping from each ri to the nearest neighbor of 

( ) /2 i jk r h     is evaluated once and used repeatedly throughout code execution. As in the previous 

section, equation (25) is accurate at low range values since the scaled ion fluxes are nearly constant in this 

region. At higher range values, the step-size perturbation has a small effect. Further, the number of grid 

points typically used for heavy ion transport provides sufficient fidelity so that the nearest-neighbor range 

value is not far from the actual perturbed value. 

It should be noted that at the first step from the front boundary, secondary ion fluxes are changing 

rapidly, and the previous heavy ion solution is more accurate than equation (26). Similar behavior is 

observed near material interfaces, but the impact is less noticeable than at the front boundary. In the 

updated heavy ion algorithm, a hybrid approach has been adopted. For the first step into any material, 

including material interfaces, the original heavy ion marching equation (24) is evaluated. All other steps 

use the simplified form in equation (26). Log-log cubic Lagrange interpolation is still used to evaluate the 

term (, ) ji j   xr h  . 

5. Flux Comparison

In this section, a flux-based comparison between the original and new heavy ion transport 

algorithms is given. In this comparison, both algorithms were used to transport a solar minimum GCR 

environment through a 100 g/cm2 aluminum slab and a 100 g/cm2 water slab, and fluxes were saved at 

various depths in both slabs. The solar minimum GCR environment was calculated with the BadhwarO'Neill 2010 GCR model [O'Neill 2010] using a solar modulation parameter of 475 MV. The relative 

difference between the results computed with both algorithms is computed as 

( ) ( )

( )

(, ) (, ) (, )

(, )

N O

j ik j ik

j ik O

j ik

Ex Ex

REx

E x

 



  , (27)
7

where ( )(, ) O j ik E x and ( )(, ) N j ik E x are the fluxes for each ion, j, at each energy, Ei, and depth, xk, 

computed with the original (O) and new (N) algorithms, respectively. It was found that the absolute value 

of the relative difference between the old and new algorithms never exceeded 5% over all energies, ions,

and depths. In order to summarize these results, the following quantities were computed as a function of 

energy, 

min ,

( ) min ( , ) i j ik j k

R E REx  (28)

max ,

( ) max ( , ) i j ik j k

R E REx  (29)

avg

1 1

1 ( ) (, )

Nx Nj

i j ik

j x k j

R E REx

N N  

  (30)

2

sd avg

1 1

1 ( ) (, ) ()

Nx Nj

i j ik i

j x k j

R E REx R E

N N  

      (31)

Equations (28) - (31) have also been computed as a function of particle type and depth. The results are 

shown in Figures 2 - 4. The symbol in each figure is the average of the relative differences, Ravg(Ei). The 

outer-most error bars represent the minimum and maximum relative differences defined in equations (28)

and (29). The inner-most error bars represent the standard deviation of relative differences defined in 

equation (31). The standard deviation error bars indicate where most of the relative differences occur, while 

the minimum/maximum error bars indicate the extreme values of relative differences. The absolute value of 

all relative differences is clearly within 5%, and most relative differences are within 2%. This indicates that 

the simpler nearest-neighbor interpolation approach gives results that are very similar to the more 

complicated and computationally expensive log-log cubic Lagrange method. It should also be noted that the 

relative differences shown in Figures 2 – 4 translate into relative differences for dose equivalent that are 

bounded by 0.3%. It will be shown later, that the increased speed provided by the nearest neighbor 

approach justifies this negligible loss in accuracy. 

Figure 2. Relative difference between fluxes computed by the original and modified heavy ion transport 

algorithm as a function of kinetic energy in aluminum and water exposed to solar minimum GCR 

conditions. Light and dark error bars represent extreme and standard deviation of relative differences, 

respectively.
8

Figure 3. Relative difference between fluxes computed by the original and modified heavy ion transport 

algorithm as a function of shielding thickness in aluminum and water exposed to solar minimum GCR 

conditions. Light and dark error bars represent extreme and standard deviation of relative differences, 

respectively.

Figure 4. Relative difference between fluxes computed by the original and modified heavy ion transport 

algorithm as a function of charge in aluminum and water exposed to solar minimum GCR conditions. Light 

and dark error bars represent extreme and standard deviation of relative differences, respectively.

6. Convergence Test

A computational algorithm is said to converge if the numerical solutions reach an asymptotic limit 

as the discretization parameters approach zero. HZETRN has two discretization parameters: the step size, h, 

and the number of energy grid points, N. In order to show convergence and quantify discretization error, 

discretization parameters are refined several times and the differences between the various solutions are 

compared. This process is referred to as a convergence test. 

In this section, convergence tests are applied to the original and modified heavy ion transport 

algorithms to show convergence and quantify overall discretization error. Five different energy grid sizes 

are considered ranging from N = 100 to N = 502, along with eleven different step-sizes ranging from h = 
9

0.5 g/cm2 to h = 2-11 g/cm2. These discretization parameters are used in both the original and modified 

heavy ion transport algorithms to transport a solar minimum GCR environment through 100 g/cm2 of 

aluminum and water. The solar minimum GCR environment was calculated with the Badhwar-O'Neill 2010 

GCR model [O'Neill 2010] using a solar modulation parameter of 475 MV. The same light particle 

transport algorithm was coupled to both the original and new heavy ion transport algorithm. Dose 

equivalent is computed at a number of depths in the slab and is used as the quantity of interest for the error 

analysis.

Dose equivalent as a function of depth in aluminum and water exposed to the solar minimum GCR 

environment is given in Figure 5 for both algorithms and all of the discretization parameters. The data in 

Figure 5 show that both algorithms with the various discretization parameters give very similar results over 

the thicknesses considered. The spread in dose equivalent values at 100 g/cm2 in aluminum and water is 

1.2% and 1.4%, respectively. This variation includes both the differences associated with the original and 

new heavy ion transport algorithm as well as the inherent discretization errors. 

Figure 5. Dose equivalent versus depth in aluminum (left pane) and water (right pane) exposed to solar 

minimum GCR conditions computed with the original and modified heavy ion transport algorithms and all 

discretization parameters.

Figure 6. Total discretization error for dose equivalent as a function of depth in aluminum (left pane) and 

water (right pane) exposed to a solar minimum GCR environment.
10

The total discretization error associated with the original and modified heavy ion transport 

algorithms for aluminum and water targets exposed to the solar minimum GCR environment is shown in 

Figure 6. For each algorithm, the discretization error is obtained by computing the relative percent 

difference between the converged results obtained with fine discretization parameters (h = 2-11 g/cm2 and N

= 502) and the solution obtained with typical discretization parameters (h = 0.5 g/cm2 and N = 100). A 

negative discretization error indicates that the results obtained with the typical discretization parameters 

underestimate the converged result. It is clear that the modified heavy ion transport algorithm results in a 

negligible increase in discretization error in both aluminum and water targets as thick as 100 g/cm2

. It will 

be shown in the next section that this minor loss of accuracy is justified by the significant reduction in code 

run time. 

7. Run-time Comparisons

The original and modified heavy ion transport algorithms are used to transport the same solar 

minimum GCR boundary condition through three different shielding configurations typically used with 

HZETRN. The solar minimum GCR environment was calculated with the Badhwar-O'Neill 2010 GCR 

model [O'Neill 2010] using a solar modulation parameter of 475 MV. These run-time tests were performed 

on a Linux computer with 8 quad-core AMD Opteron 8378 processors and 256 GB of RAM. Each test was 

run 25 times with the original and new algorithms, and the reported run-times are the average of the 25 

separate executions. Results are given in Table 1. The column labeled "1 layer" corresponds to a shielding 

configuration made of a single layer of aluminum. This type of calculation might be done to analyze the 

radiation protection properties of a slab of some material or to generate a simplified interpolation database 

for vehicle analysis. The column labeled "2 layer" corresponds to a shielding configuration made of an

aluminum shield followed by a tissue target. Fluxes are stored on a two-dimensional grid of points covering 

every combination of first and second layer thicknesses. This simulation is typically used to model a 

shielding structure (aluminum) surrounding a human phantom (tissue). The column labeled "3 layer" 

corresponds to an aluminum shield followed by a polyethylene shield and a tissue target. This simulation is 

typically used to model a shielding structure (aluminum) with a secondary shield (polyethylene) 

surrounding a human phantom (tissue). In all three cases, the maximum thickness of any layer is 100 g/cm2

, 

and each layer has 11 spatial grid points spaced nearly logarithmically at which fluxes are saved and 

printed to a file. Results are given in Table 1. 

Table 1. Run-time (seconds) comparisons between original and new heavy ion transport algorithms.

1 layer 2 layer 3 layer

Original 9.2 96.1 1042.3

New 2.9 20.1 212.4

Factor 3.2 4.8 4.9

The overhead computational cost of reading in cross section databases and writing out flux files is the same 

in both codes and is included in the total run-time. Consequently, the ratio between the original and new 

run times is not constant. Two additional 1-layer tests were run, but the slab thickness was increased to 

1000 g/cm2 and 5000 g/cm2

, and the fluxes were only printed out at the front and back boundaries of the 

slab to minimize the run-time dedicated to file writing. These tests help minimize the impact of reading and 

writing data on total run-time and more clearly show the increased speed associated with the new 

algorithm. It was found that the new code was 5.5x and 5.8x faster than the old code in the 1000 g/cm2 and 

5000 g/cm2 slabs, respectively. 

8. Conclusions

In this work, a modification to the heavy ion transport algorithm in HZETRN was shown to be 

approximately 5x faster than the original algorithm, and convergence tests showed that the loss in 

numerical accuracy was negligible (< 1%). This increase in overall code run time for GCR calculations will 
11

be beneficial for future analysis including risk assessment, ray-by-ray calculations, and geometries with 

large shielding thicknesses such as the International Space Station and martian atmosphere.

9. References

Bethe, Von H., Zur Theorie des Durchgangs schneller Korpuskularstraheln durch Materie. Annals of

Physics, 5 Folge, Bd. 5, pp. 325-400 (1930).

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

Isotopic Dependence of GCR Fluence behind Shielding. Radiation Measurements, Volume 41, pp. 1235-

1249 (2006).

Norman, R.B., Blattnig, S.R., De Angelis, G., Badavi, F.F., Norbury, J.W., Deterministic Pion and Muon 

Transport in Earth's Atmosphere. Advances in Space Research, Volume 50, pp. 146-155 (2012a).

Norman, R.B., Slaba, T.C., Blattnig, S.R., An Extension of HZETRN for Cosmic Ray Initiated 

Electromagnetic Cascades. Advances in Space Research, accepted for publication (2012b).

O'Neill, P.M., Badhwar-O'Neill Galactic Cosmic Ray Flux Model – Revised. IEEE Transactions on 

Nuclear Science, Volume 57, pp. 3148 – 3153 (2010).

Shinn, J.L., John, S., Tripathi, R.K., Wilson, J.W., Townsend, L.W., Norbury, J.W., A Fully EnergyDependent HZETRN (A Galactic Cosmic-Ray Transport Code). NASA Technical Paper 3243 (1992).

Slaba, T.C., Blattnig, S.R., Badavi, F.F., Faster and more Accurate Transport Procedures for HZETRN. 

Journal of Computational Physics, Volume 229, pp. 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. Radiation Measurements, Volume 45 pp. 173-182 

(2010b).

Slaba, T.C., Blattnig, S.R., Badavi, F.F., Faster and more Accurate Transport Procedures for HZETRN. 

NASA Technical Paper 2010 – 216213 (2010c).

Slaba, T.C., Blattnig, S.R., Clowdsley, M.S., Variations in Lunar Neutron Dose Estimates. Radiation 

Research, Volume 176, pp. 827-841 (2011a).

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. Advances in Space Research, Volume 47, pp. 600-610 (2011b).

Slaba, T.C., Blattnig, S.R., Tweed, J., Reduced Discretization Error in HZETRN. Journal of Computational 

Physics, Volume 234, pp. 217-229 (2012a).

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. 

Advances in Space Research, submitted (2012b).

Wilson, J.W., Analysis of the Theory of Heavy Ion Transport. NASA Technical Note D-8381 (1977).

Wilson, J.W. and Badavi, F.F., Methods of Galactic Heavy Ion Transport. Radiation Research, Volume 

108, pp. 231-237 (1986).
12

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 Reference Publication 1257 (1991).

Wilson, J.W., Badavi, F.F., Cucinotta, F.A., Shinn, J.L., Badhwar, G.D, Silberberg, R., Tsao, C.H., 

Townsend, L.W., Tripathi, R.K., HZETRN: Description of a Free-Space Ion and Nucleon Transport and 

Shielding Computer Program. NASA Technical Paper 3495 (1995).

Wilson, J.W., Tripathi, R.K., Badavi, F.F., Cucinotta, F.A., Standardized Radiation Shield Design Method: 

2005 HZETRN. SAE ICES 2006-18 (2006).
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Publication

 4. TITLE AND SUBTITLE

Faster Heavy Ion Transport for HZETRN

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Slaba, Tony C.

 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-20222

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 deterministic particle transport code HZETRN was developed to enable fast and accurate space radiation transport through 

materials. As more complex transport solutions are implemented for neutrons, light ions (Z < 2), mesons, and leptons, it is 

important to maintain overall computational efficiency. In this work, the heavy ion (Z > 2) transport algorithm in HZETRN is 

reviewed, and a simple modification is shown to provide an approximate 5x decrease in execution time for galactic cosmic ray 

transport. Convergence tests and other comparisons are carried out to verify that numerical accuracy is maintained in the new 

algorithm. 

15. SUBJECT TERMS

HZETRN; Heavy ions; Radiation transport; Space radiation

18. NUMBER

 OF 

 PAGES

20

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-2013-217803

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- 02 - 2013