LA-2120

PHYSICS AND MATHEMATICS

(TD-4500,13th Ed., Suppl.)

LOS ALAMOS SCIENTIFIC LABORATORY

OF THE UNIVERSITY OF CALIFORNIA LOS ALAMOS NEW MEXICO
REPORT WRITTEN: January 1957
REPORT DISTRIBUTED: December 18, 1957

A PRACTICAL MANUAL ON THE MONTE CARLO METHOD

FOR RANDOM WALK PROBLEMS

Work done bv: Report written by:

E . D. Cashwell E , D. Cashwell
C. J . Everett C. J . Everett
0. W. Rechard

Contract W-7405-ENG. 36 with the U. S. Atomic Energy Commission

r 1
 ABSTRACT

This r e p o r t is w r i t t e n t o serve as a guide t o those persons

who, having no previous experience with Monte Carlo methods, wish t o

apply these methods t o t h e i r own problems. P a r t i c u l a r emphasis i s

given t o techniques which a r e u s e f u l i n dealing with problems concerned

with the d i f f u s i o n of p a r t i c l e s (and gamma r a y s ) i n m a t e r i a l media of

some complexity, both from a geometrical and a nuclear standpoint.

Included as a n appendix a r e b r i e f summaries of a v a r i e t y of problems

of t h e above-mentioned type t o which t h e methods described herein have

been applied successfully.

-3-
 FOREWORD

The present r e p o r t i s a summary of the Monte Carlo method as it

a p p l i e s t o problems involving the i n t e r p l a y of neutrons and photons w i t h

bulk w t t e r i n g e o m t r i c systems of varying complexity, It is intended

t o serve as a n introduction and p r a c t i c a l guide f o r the fast-growing

group of people who are concerned with such systems,

A b r i e f res& of some of t h e u n c l a s s i f i e d problems s u c c e s s f u l l y

treated by the methods here o u t l i n e d i s included as a n appendix and may

serve t o emphasize the p r a c t i c a l i t y and f l e x i b i l i t y of these sampling

procedures. A similar res=’ of some of the classified problems which

have been handled i s issued s e p a r a t e l y as LAMS-2123,

The general method was o r i g i n a l l y developed by Fermi, U l a m and

von Neumann; many o t h e r s have contributed s p e c i a l techniques and devices.

No attempt has been made t o give all sources i n t h e t e x t , and n a t u r a l l y

no claims t o o r i g i n a l i t y a r e made f o r the procedures included.

- 5-
 CONTENTS
Page
ABSTRACT . * e 3
FOREWORD . . . . . . . . . . . . . . . . . . .
e

5
CHAPTER I . BASIC PRINCIPLES . . . . . . . . . . . . . 11

1 .. General nature of the problem

. . . . . .
. . . .. .. .. .. . ..
.
11
2
3 .. Outline of procedure
Production of random numbers . . . . . . . .
.
.. .. .. ..
.. .
..
12
15
4
5 . The fundamental p r i n c i p l e of Monte Carlo
Application of t h e p r i n c i p l e . . . .
16
21

CHAPTER I1 . TKE SOURCE ROUTIME . . . . . . . . . . . . 25

1.. Introduction . . .. .. .. .. .. .. .. .. .. .. .. .. .. 25

3.
27
2

4.
Remarks on u n i t s . . . . . . . . .. ..
P a r t i c l e parameters

Space coordinates f o r source particles .. .. .. .. 32

33
a
b .. Uniform source on a n annulus of radii R o < R1
Uniform source i n a s p h e r i c a l s h e l l . . . .. .. .. 34
34
e.
. . . . . . . . . .
Parallel-beam source incident on lateral c y l i n d r i c a l

d . surface. o r on a sphere
.
I s o t r o p i c p o i n t source e x t e r n a l t o cylinder
34
36
5 . Direction coordinates f o r source p a r t i c l e s . . . . . 36
a
b
.. . . . . .
I s o t r o p i c source; U.V. w d i r e c t i o n cosines
The cosine d i s t r i b u t i o n . .. .. .. .. 36
38
I s o t r o p i c and cosine sources i n s p h e r i c a l l y symmetric
. . . . . . . . . . . . . . .
d.
C.
systems 40
. . .
e.
I s o t r o p i c point source e x t e r n a l t o cylinder 40

f. . . . . . . .
General d i s t r i b u t i o n i n h a l f of direction-space
A prejudiced source . . .. . 42
43
. . . . . . . . . .
6.
7
8.
. Energy of particles
source
Other source pameters . . .
Source f o r u-type c a l c u l a t i o n s .. .. .. .. .. .. .. ..
45
46
46

-7-
CHAPTER I11 . THE MEAN FREE PATH AND TRANSMISSION v
Page
49

. . . .. .. . . . . . .
. . ... ...
2 ..
1. The c r o s s s e c t i o n concept
The =an free path
. . . . . . . . . . . . .
49
50

..
3 An example
4 rtSmU"systems and transmission . .. . .. .. .. .. 53
54
..
5 The "forced first c o l l i s i o n " r o u t i n e
. . .
6 Remark on the device i n s p h e r i c a l problems
. . .
.
.. ..
. 56
58
.
7 The transmission i n subsequent h i s t o r y
. . +

a Prejudiced f i r s t c o l l i s i o n i n "large" systems

61
61
CHAFTER N . THE COLLISION OR ESCAPE ROUTINE . . . . . . . 63
1. Introduction . . . . . . . . . . . . . . .
2 . A r o u t i n e f o r the s p h e r i c a l shell . . . . 63
3.
65
4. Flux problems i n s p h e r i c a l geometry . . ..
60
Reorientation formulas f o r t h e s p h e r i c a l s h e l l
.
5 . A r o u t i n e f o r the f i n i t e c y l i n d e r . . .. 72
6 . The f i n i t e c y l i n d r i c a l s h e l l w i t h c e n t r a l hole .. 72

. . . . . . . . .. .. . .. ..
74
7 . The s p h e r i c a l s h e l l i n a b s o l u t e space
8 . Slab geometry . 78

9 .
Problems run i n cycles of t i m e AT . . . . . . . 80
81
CHAPTER V. THE COLLISION ROUTINE FOR "I'RONS e 0 e
83 .
1..Introduction . . . . . . . . . . . . . . . 83
3.
2

4.
E l a s t i c c o l l i s i o n s i n general . 8 -
Capture and s e l e c t i o n of t h e type of c o l l i s i o n
9
84
89

5.
The d i f f e r e n t i a l e l a s t i c s c a t t e r i n g c r o s s s e c t i o n 99
.
.
A routine f o r e l a s t i c scattering
. .
6 D i f f e r e n t i a l e l a s t i c c r o s s s e c t i o n f o r the lab system
102

.. . . . . . . . . . . . . . . . ..
. . . . . 108
7 A weight device f o r e l a s t i c s c a t t e r i n g 109
8 Fission
. . . . . . 112

..
9 Inelastic (n-n) c o l l i s i o n s i n general
. . . . .
10 I n e l a s t i c (n-n) c o l l i s i o n s on heavy n u c l e i
.. .
114
118
.
11 I n e l a s t i c (n-n) c o l l i s i o n with Maxwell d i s t r i b u t i o n
.. .. .. .. .. .. .. ..
1 2 A combined t r a n s f e r matrix f o r f i s s i o n a b l e nuclei
119
121
.
1 3 C o l l i s i o n s s h a t t e r i n g a nucleus
14 The (n-2n) r e a c t i o n i n deuterium
124

. . . . .. . . . .. .. .. .. ..
129
.
1 5 An (n-2n) r e a c t i o n on heavy n u c l e i
16 Capture i n a small zone
132

. .. .. .. .. .. .. .. .. ..
133
Si
.
' Capture by a "point" d e t e c t o r
18 Remarks on thermal neutrons
19 Remark on determination of photon sources . . . . .
136
140
141

-8-
 . . . . . . . Page
CHAPTER V I . PHOTON COLLISIONS . . 143
. . . , . . ..
1.
2, Basic concepts and constants . . . .
Introduction

. . . . . . .
I

.
... .. .. . . 143
143
. . . . 149
3. Compton c o l l i s i o n s 144
4, The Klein-Nishina d i f f e r e n t i a l cross s e c t i o n
5. The photon energy d i s t r i b u t i o n and Compton cross s e c t i o n 150
6. Photoelectric e f f e c t and p a i r production . . . . 153
CHAPTER V I I . DIRECTION PARAMETERS AFTER COLLISION . 155
1. Introduction . . . . . . . . . . .. . .. 155
2. Formulas f o r t h e f i n a l d i r e c t i o n cosines . 155
.
e

3. Subroutine for t h e f i n a l d i r e c t i o n cosines 161

. 164
4. F i n a l d i r e c t i o n w i n slab o r s p h e r i c a l l y symmtric case 163
5. S c a t t e r i n g i s o t r o p i c i n the laboratory system .
CHAPTER VIII. T-NAL CLASSIFICATION . . . . . . . . . 165
1. Introduction . . . . . , . . . . . . 165
2. C l a s s i f i c a t i o n of escapes on nurnber of c o l l i s i o n s .
3. Energy and angle d i s t r i b u t i o n s of eacape . . e . 165
166
CHAPTER IX. REMARKS ON COMPUTATION . , . , . . . . . 175
Scaling . . . . , . . . . , , . . . . *175
1,
2.
3.
Debugging . . , . . , , , .
Special subroutines . '. . . , ... . . . . . . 177
175

. . . . .. .. .
b. Shifted random numbers . . . . . .. . . 177
a. A random nuniber r o u t i n e

A logarithm r o u t i n e . . . . . . . . . 178
C.
. . . .' . . . . . 179
..

. . . . . .
d. The exponential exp(- x/y) 182
e. A cosine r o u t i n e b . . *I83
4. A Monte Carlo device f o r r. . . . . . . . 186
5.
. ..
A Monte Carlo device f o r t h e cosine of a n e q u i d i s t r i b u t e d
angle..b,?**... .a87
CHAPTER X. STATISTICAL CONSIDERATIONS * 189
1. The limit theorem i n t h e Bernoulli case
. . . . . * 189
. . . . .. . . .
2. Application t o the terminal r a t i o s 190
3. The c e n t r a l l i m i t theorem + 193
195
. . .
4. Application t o problems using weights
5. I U u s t r a t i v e examples . . 197

-9-
APPENDIX . SUMMARY OF CERTAIN PROBUNS RUN ON MANIAC I . . . . . 204
Page

Problem1 . . . 205
Compton c o l l i s i o n s i n a spherical s h e l l

Problem 2 .
Compton c o l l i s i o n s i n a s o l i d sphere . . 205
Problem 3 . Neutrons i n a spherical s h e l l . . 206
Problem 4 . Energy independent s c a t t e r i n g i n a cylinder
0

. . 207
Psoblem 5 . Energy dependent s c a t t e r i n g i n a cylinder . . . 208

Problem 6 . Cylindrical shell w i t h p a r a l l e l source . . . . 209

Problem 7 . 1 4 Mev neutrons i n a c y l i n d r i c a l s h e l l . . . 210

Problem 8. Cylindrical sheU. with point source . . . 213

Problem 9 . A s c i n t i l l a t i o n counter . . . . . . . . . 215

Problem 1 0 . A "'long counterll . . . . . . . . . . . 216

Problem 11. problem connected neutrino detection . . 217
A with

Problem12 . determination . . . . . . . . . . 219

An u

Problem 13. Neutron flux i n air . . . . . . . . . . 220

Problem 14 . Escape d i s t r i b u t i o n from a carbon slab . . . . . 221

Problem 1 5. Escape d i s t r i b u t i o n from a spherical s h e l l . 221
Problem 16 . A rocket motor . . . . . . . . . . . . 222
Problem17 . Photon diffusion i n a reactor . . . . . . . 225
Problem 18. . . . .
A thick Z r t a r g e t problem . . . . 226
Problem19 . A thick C t a r g e t problem . . . . . . . . 226
Problem 20 . Heavy water experiment . . . . . . . . . 227

-10-
 CHAPTER I

BASIC PRINCIPLES

1. General nature of the problem. A l l problems treated in the

present manual involve estimation of what percentage of Wrticles emanate

ing from a given source, after undergoing specified processes in a

material medium of known geometry, can be expected to terminate in certain

stipulated categories.
If all relevant probabilities are known f o r the elementary events
in the "life history" of such a particle, the Monte Carlo method is
applicable, and indeed is usually the only method available.

Moreover, its technique is pre-eminenkly realistic, consisting in

actually following each of a large number of particles from the source
throughout its life history to its ''death" in some one of the terminal
categories, using the elementary probabilities at each stage of its
career in determining its fate.
The present state of development of high-speed digital, computers
permits the use of sasflples of a size sufficiently large to ensure satis-
factory accuracy in most practical problew.

-3.1-
 2. O u t l i n e of procedure. I n any p a r t i c u l a r problem, a p a r t i c l e

i s completely characterized by a s e t of parameters which a r e s u f f i c i e n t

t o determine i t s ( p r o b a b i l i t y ) behavior i n a l l s i t u a t i o n s it may encounter

during i t s h i s t o r y . These always include i t s p o s i t i o n and d i r e c t i o n

I

coordinates, and i n most cases i t s energy. Much more w i l l be s a i d about

these and other p a r t i c l e parameters i n the sequel ( t h e appendix t o t h i s

r e p o r t and LAMS-2121)
The Monte Carlo method of dealing w i t h problems of the kind we have

indicated breaks up n a t u r a l l y i n t o a well-defined s e t of subroutines,

which we s h a l l b r i e f l y describe here, postponing t h e i r d e t a i l e d treatment

t o subsequent chapters. They are schematized i n the following generalized

flow diagram (Fig. 1). This i s only intended as a general guide t o t h e

chapters of the t e x t , and i s s u b j e c t t o many r e v i s i o n s , depending on t h e

s p e c i a l circumstances of the problem.

(0). The proper assignment of a l l p a r t i c l e parameters t o a source

p a r t i c l e involves the spatial and angular d i s t r i b u t i o n of the source, and

i t s energy d i s t r i b u t i o n , i n case of a non-monoenergetic source.

(Po). A s p e c i a l r o u t i n e may be provided t o determine t h e p o s i t i o n

of f i r s t c o l l i s i o n i n cases of high transmission, where it i s d e s i r a b l e

'
t o d i s t i n g u i s h between t h i s and subsequent c o l l i s i o n s .

(B). A general r o u t i n e i s designed t o decide whether a p a r t i c l e ,

s t a r t i n g w i t h known parameters from a n a r b i t r a r y point of departure, i n

a p a r t i c u l a r zone of t h e system e i t h e r s u f f e r s a c o l l i s i o n within t h i s

zone o r reaches the boundary of t h i s zone on i t s l i n e of f l i g h t without

incident. The e s s e n t i a l physical concept involved i n t h i s decision is

-12-
@) . . - i+c Source parameters

I I

@I Parameters at forced 1st

collision (Ch. III)
Collision. Parameters at
I

point of colIision

I
I I
Passage to next
, zone

Angle of sgattering
New energy

Lost to energy

*I Direction after
scattering (ch. VII) t+~
Classification of
0-1 escapes (Ch. VIII)

Fig. 1

-13-
 the mean free path. The subsequent procedure depends on the nature of

this decision, as indicated i n Fig. 1. The particle parameters at the

point of collision ox escape are computed before proceeding. In case

of escape, one proceeds to ( 8 ) or to ( 6 ) according as to whether

the escape is to an adjacent zone, or from the system.

(r) . The collision routine naturally depends upon the physics of

the medium and is subject to the widest variation. Its objective is to

decide the exact nature of the collision and the immediate fate of the
particle after collision. This includes the eventuality that the particle
may terminate its career under the conditions of the problem; if this
occurs, a tally is made in the appropriate terminal category counter,
and one returns, as in all cases of termination, to an entry ( u ) which
leads to the source routine ( 0 ) for introducing a fresh particle. In
the event that the particle is to be followed further, the essential

. information required before proceeding to (6) is a knowledge of the

laboratory angle of deflection from the line of flight and the energy

and "weight" of the particle after deflection.

(6 1. A purely geometric routine determines the direction param-

eters of the deflected particle, and leads to (g). At this point, all
particle parameters should be evaluated as they obtain at the point of
collision, after deflection.

(e). In the event that escape *om the system occurs in the ( p )
routine, terminal classification is made and one returns to (a).
 This is the general scheme of the method as it operates in all

problems we propose to discuss. The chapters that follow are designed

to elaborate on each of these subroutines in detail, the necessary
physical concepts and their preparation for computation being developed
as they are needed.
Before this, however, must come a brief discussion of random
numbers and the "fundamental principle of Monte Carlo," upon which all
else rests.

3. Production of random numbers. It is necessary to have upon

call some source of random numbers equidistributed on the interval

Osr<l. Ideally, one might spin a wheel of uniform scale, and indeed

there exist of such numbers generated in truly random fashion.

There are computational algorithms adapted to digital computers which
appear to serve our purpose just as well. One can find descriptions of
such methods in the literature,
('j2) together with discussions of the

(l)Cf. articles in U. S, Department of Commerce, National Bureau o f t

Standards, Applied Mathematics Series #12, Monte Carlo Method,

Washington, D. C., 1951.
(2)The Rand Corporation, A Million .Random
, -
_ Digits with 100,000 Normal
Deviates, Free Press PGblishers, Glencoe Illinois, 1955
D; H, Lehmer, Mathematical Methods Large-Scale Computinq Units,
Proceedings -
of-a Second Symposium on Large-Scale DiKital Calculatin
Machinery, Harvard University Pres< Cambridge, Massachusetts, 1951:
Herbert A. Meyer, ed., S y y s i u m on Monte Carlo Methods, John Wiley &
Sons, Inc., New York, 195
"tests of randomness" which have been a p p l i e d t o them. We do not e n t e r

here upon such questions s i n c e t h e method adopted w i l l probably be

d i c t a t e d by t h e type of machine used and o t h e r p r a c t i c a l considerations.

A d e s c r i p t i o n of t h e method w e have employed w i l l be found i n Ch. I X ,

9 3a, b. We should l i k e t o c a l l a t t e n t i o n e s p e c i a l l y t o t h e s h o r t paper

by von Neumann, c i t e d i n footnote (l), which contains some words of

wisdom f o r those who may be troubled by t h e " s t a t e of s i n " accompanying

t h e use of d e t e r m i n i s t i c "random numbers .".

4. The fundamental p r i n c i p l e of Monte Carlo. Suppose t h a t a

homogeneous medium c o n s i s t s of n u c l e i of t h r e e d i f f e r e n t types A, B, and

C , and one knows t h a t , i n t h e event of a c o l l i s i o n of a neutron i n t h e

medium, t h e p r o b a b i l i t y of c o l l i s i o n with type A is .2, with B is

.3, and with C i s t h e remaining .5. It i s i n t u i t i v e l y c l e a r t h a t ,

i f a l a r g e number N of random numbers are produced, approximately

.2 N w i l l f a l l on t h e i n t e r v a l 0s r < . 2

.3 N w i l l fall. on the i n t e r v a l . 2 1 r < .5

.5 N w i l l f a l l on t h e i n t e r v a l . 5 6 r <1
and t h a t t h i s approximation w i l l improve with i n c r e a s i n g N; indeed,

t h i s i s t h e s a l i e n t f e a t u r e we demand of any scheme of random number

production. It i s c l e a r t h e r e f o r e how reference t o a random number can

be used t o decide which of t h e t h r e e types of n u c l e i i s h i t i n t h e event

of a colli'sion.

-16-
 One of the decisive features of digital computers now existing;
is their ability to make decisions at high speed, with no limitation on

the number of logical possibilities involved. The typicalMon$e Carlo

flow diagram is largely occupied with intricate decision features.

As a first example of how a flov diagram(3) schematizes the pro-
cedure in the above example we may note Fig. 2.
E r e generally, i f El, ..., En are n independent, mutually

exclusive events with probabilities pl, ..., pn, respectively, and

p1 + + pn = 1, we will agree that

P1 + ... + pi-ls r ( p l + ... + *i

detemnines Ei. This Is the "fundamental principle'' bsofw as it applies

to the discrete case of a finike number of eventualities.
1,
We may use it to furnish a heuristic approach t o the continuous
case, in which it is desired to determine from a random number one of a
continuum of values of a variable x. In this connection it is conven-

ient to bear in mind the fact that the latter case may always be regarded
as approximable by a large finite number of distinct Cases; indeed, in
computation with a fixed number of digits, we are always actually

concerned with the discrete approximation.

'3)1n all flow diagrams of this report we adopt the convention: x4-
means x30, x- means X C O . A box covtaining r always denotes
reference to the special subroutine generating the next random number
of the sequence.
Fig. 2

-18-
 Suppose that we arbitrarily assign a variable x on t h e m 1

O S x c n toth

i -1 x < i represents the event -Ei. Let us construct a grobability

density function p( x) by the definition

i = 1, ,.., n

Thus p(x) Will be a step f'un~tion'~)like that of Fig. 3. Now

suppose that we define the probability distribution f'unction

n
X

whose graph is a, broken line as indicated in F i g . 4. Note that P(0) = 0,

P(n) = 1. Since P(i) = p1 + ... t pi, we may interpret P(X) to mean

the probability of the inequality 4 sx, for x = 1, i = lc.*

* t - n.
Moreover, it is clear that the equation

determines x uniquely --
as a function of r, m sue,* a way that if ,, r
is uniformly distributed on the interval 0 d r e 1, x , falls,WlLh
.u--..----..""-.--A-

frequency p. in the interval

..- -~_yl
i ~ "
- 1-% x r i ,
F - L I - X - A
thereby determining the
event Ei under our agreement.

F)The figures are drawn for the simple A, B, C example.

.5
.4
.3
.2
.1
X
1 2 3

Fig. 3

Fig. 4

-20-
 We may state at once the fundamental principle as it applies to
the continuous case: If p(x) dx is the probability of x lying
between x and x -I-dx, with a f x < b, and

then
P
a
p(S) d0 = 1

determines x uniquely as a function of r; moreover, if r is

uniformly distributed on 0 r < 1, then x falls with frequency
p(x) dx in the interval (x, x I- d x ) .

p s a completely trivial first example, consider the case of

neutrons that are to be located uniformly on an interval a 8 x <b.

We have p(x) dx = dx/(b-a), [ , and
p(k) dS = 1 r = P(x) =
[ p( 6) dS = (x-a)/(b-a), whence x = a m (b-a) determines x as a
function of the random number r

5. Application of the principle, It may be noted that the

equation

-21-
can be expected to give rise to difficult implicit problems, since x

must be determined from r. At worst, some successive approximation

routine can be provided for the solution of the implicit equation

r = P(x) when P(x) is obtainable in closed analytic form. Such time-

consuming processes can be obviated in various ways, a few of which we

indicate here.
The simplest method, applicable in a l l cases, even when P(x)
is known only in experimental tabular form, consists in subdividing the
(a,b) interval, storing accurate values of P(xi) =, Pi for the end
points xo = a < xl < ... < xn = b of the subintervals, and using the

discrete method for determining the subinterval (xi - xi) on which

x falls, together with an interpolation for the actual value of x.
If i is the first value of the index for which r - Pi is negative,
r being the current random number, we may determine x from one of the
formulas

Pi - r
x = xi - Pi - Pi,l (Xi - xi-11

-22-
 The linear interpolation in (1) distributes x uniformly on the
interval (xi I, xi) and is strictly valid only when p(x) is a step
function. For sufficiently small subdivisions, (2) or ( 3 ) may give
better results at the cost of an additional square root and are appro-

priate when P(x) is concave up or concave doqrespectively,

One may contrast the latter formulas (2) and ( 3 ) with that
resulting from a linear assumption for p( x) on ( xiel, x~):

which is more complicated when solved for x, requires additional

storage of the pi, and uses a trapezoid rule for the Pi.
Very usefuL also is a device employed by von Neumann, especially
when p(x) is readily computable and storage space is at a premium.
This consists of "throwing" points ( 0 , V ) unifomily into the rectangle

bounded by the lines x = a, x = b, y = 0, y = 1 and rejecting the

points lying above the curve

x being assigned the value 6 whenever (6, TI) f a l l below the curve. In
many trials, the ratio of the number of points retained with 6 between

x and x I- dx to the number of points retained altogether will be

approximately the ratio of the areas

-23-
 The method is illustrated by the flow diagram of Fig. 5.
Obviously the device in this form is impractical if the area under the
curve y
*(x)
= p is a small fraction of the enclosing rectangle. How-
ever, modifications involving other enclosing areas can obviously be

devised.
*
One may also retain only points (S,.rl) above the curve y = p (x),

assigning the value 0 to x and adjusting the “weight“ of the particle

by a factor P(X>
a
1
b
(1 - p*(5))rl5/(1 - P”(x)).
Finally, a combination of the two methods may be used, a first
random number determining the interval (xi 1, xi) by reference to the
and the von Neumann device then used on p( x) on this interval.
Pi,
The method is then accurate, and the efficiency high.

Fig. 5
 cHAFTm I1

THE SOURCE ROUTINE

1. Introduction. It is the purpose of t h e present chapter t o

describe how a problem is i n i t i a t e d by t h e machine, how "print-outs"

are automatically e f f e c t e d , and how the p a r t i c l e s are drawn from t h e

source.

Suppose t h a t we l e t N denote t h e number of p a r t i c l e s already

processed, so t h a t N = 0 a t t h e s t a r t of t h e problem. A f t e r the

machine has processed any given number N of p a r t i c l e s it w i l l contain

i n various counters the numbers Ni of these N whose c a r e e r s have

terminated i n a set of d i s j o i n t , a l l - i n c l u s i v e c a t e g o r i e s Ci. Thus

the r a t i o s Ni/N c o n s t i t u t e t h e output of t h e problem and serve as

estimates of the p r o b a b i l i t i e s of a source particle terminating i n the

various c a t e g o r i e s Ci. It is o r d i n a r i l y d e s i r a b l e t o p r i n t t h e cumula-

,
tive totals Ni
p e r i o d i c a l l y during the course of the problem, say
~

every N* p a r t i c l e s , so t h a t convergence and "reasonableness" can be

observed

Thus t h e beginning of a flow diagram usually resembles t h a t i n

-25-
Fig. 6. It w i l l be noted t h a t N' denotes the number of p a r t i c l e s

processed s i n c e the l a s t p r i n t - o u t , being reset t o zero after each p r i n t ,

whereas N cumulates during the e n t i r e problem. The computation may be

stopped a t any t i m a t which t h e s t a b i l i t y of t h e Ni/N or other

s t a t i s t i c a l considerations i n d i c a t e t h a t s u f f i c i e n t accuracy has been

attained. The ( a ) e n t r y , as has been mentioned, i s t h a t p o i n t t o

which t h e machine r e t u r n s after t h e p a r t i c l e it was following has

terminated i t s c a r e e r i n some one of the c a t e g o r i e s

Ci.

pGV-1-

Fig. 6

-26-
 The (u) e x i t l e a d s t o that p a r t of t h e flow diagram devoted t o

assigning p a r t i c l e parameters t o a source p a r t i c l e .

There are two types of storage involved i n a l l machine problems.

"Permanent storage" places a r e reserved f o r constants l i k e 0, 1, N,

*
which do not change t h e i r values during t h e course of the problem,

while "dynamic storage" refers t o storage p o s i t i o n s i n the machine which

contain values of parameters l i k e N, N', Ni, which are s u b j e c t t o

change as the problem progresses.

It is d e s i r a b l e t o keep a record of a l l storage of both k i n d s

introduced i n t o t h e flow diagram as it i s constructed so t h a t none be

overlooked, and so t h a t some estimate of t h e "size" of t h e problem can

be gained as one proceeds. Occasionally it becomes c l e a r that tlp

memory of the machine is being exceeded, and various devices must be

introduced f o r reducing the permanent storage (e.g., by multiple

s t o r a g e ) or the length of the code i t s e l f

2, P a r t i c l e parameters, F'rom a consideration of t h e physical

and geometric f e a t u r e s of the problem, one fixes upon a s e t of p a r t i c l e

parameters whose values a t any time s u f f i c e t o completely c h a r a c t e r i z e

the p a r , t i c l e . We proceed t o discuss these i n detail.
We s h a l l limit our examples t o two types of coordinate systems
f Q r space and d i r e c t i o n ; namely, we s h a l l either employ space

coordinates x, y , z, together w i t h d i r e c t i o n cosines u 5: cos u ,

v = cos p , w = cos y f o r d i r e c t i o n of f l i g h t , where u, p, Y are

-27-
the angles made by the l i n e of f l i g h t w i t h the x, y, z axes, r e s p e c t i v e -

l y , o r ( i n cases of s p h e r i c a l symmetry only) we s h a l l use the r a d i a l

distance R of the p a r t i c l e from the center 0, together with t h e

cosine w = cos y of t h e angle y which t h e d i r e c t e d l i n e of f l i g h t

makes with t h e p o s i t i v e l y d i r e c t e d r a d i u s v e c t o r . These coordinates

are i n d i c a t e d i n Fig. 7. We may note t h a t 0 01, p ,y 6 IF and

on t h i s range the cosines assume a l l values on t h e range -16u, v,

w 6 + 1 once and only once. It i s a l s o h e l p f u l t o remember t h a t t h e

d
---i r e c t i o n coordinates (u,
_ _ ~-v, w) -_-may
_I - be regarded as defining a point on

t h e u n i t sphere
L
u2 + v2 + w2 = 1 i n d i r e c t i o n space U, V, W.
Considerable advantages a t t e n d t h e use of parameters R, w in

cases where s p h e r i c a l symmetry w a r r a n t s it. However, we have found

t h a t i n more complicated geometries, f o r example i n c y l i n d e r s , even

when symmetry o b t a i n s , the use of coordinates i n d i c a t e d by the symmetry

i s hardly worth while. The main reason f o r t h i s i s t h a t a p a r t i c l e

which proceeds from some point of departure t o a new place changes i t s

d i r e c t i o n a l coordinates i n the l a t t e r case. The x , y, z; u,v, w system

has t h e g r e a t advantage t h a t u, v , w remains unchanged under linear

displacements.

A l l problems r e q u i r e t h e use of s p a t i a l and d i r e c t i o n a l coordinates.

Other parameters are d i c t a t e d by the nature of the problem.

Most problems are energy dependent; t h a t i s t o say, the physical

processes involved have elementary p r o b a b i l i t i e s which are functions of

the p a r t i c l e energy. I n such cases we c a r r y a n energy parameter E

-28-
 Z

0 Y

Rk
Fig. 7

-29-
and usually an energy group index g. Cross s e c t i o n s are usually t o o

complicated as functions of energy t o be a c c u r a t e l y f i t t e d by simple

formulas, so that i n p r a c t i c e t h e whole range of energy involved i n

t h e problem i s subdivided i n t o s u i t a b l e groups g = 1, ..., G by lower

bounds E l > E2 > ... > EG, and a l l necessary functions of energy are

t a b u l a t e d f o r these i n t e r v a l s .

Moreover, systems encountered i n p r a c t i c e are frequently non-

homogeneous, being composed of zones 3= 1, ..., a- of varying d e n s i t i e s

o r of d i f f e r e n t materials. This usually involves s t o r i n g physic,al

q u a n t i t i e s as functions of 2 as w e l l as of g. Permanent storage

then contains numbers K with two independent i n d i c e s , which

,3+
n e c e s s i t a t e s a n a d d i t i o n a l p a r t i c l e parameter 3 t o i n d i c a t e the zone

p r e s e n t l y occupied by t h e p a r t i c l e .

I n problems involving high transmission, capture, f i s s i o n , (n-2x1)

r e a c t i o n s o r o t h e r such features, it i s d e s i r a b l e , although not

necessary, t o introduce a p a r t i c l e parameter W, c a l l e d i t s weight,

which i s i n i t i a l l y unity a t the source. To i l l u s t r a t e i t s use, consider

a problem i n which, upon c o l l i s i o n of a neutron with a nucleus, t h e r e

is a p r o b a b i l i t y p of capture. W e have the a l t e r n a t i v e s of ( a ) not

employing weights, using a random nuniber r i n case of c o l l i s i o n t o

determine whether capture occurs o r not, by reference t o t h e capture

probability p, and i f r 4 p, l o s i n g the neutron t o a capture cate-

gory, r e t u r n i n g t o (01); o r ( b ) using a weight parameter W, t a l l y i n g a

weight pW i n t h e capture category d e t e r m i n i s t i c a l l y , a n d continuing

-30-
with a neutron of weight (1-p)W, which now s c a t t e r s on t h e proper

nucleus. I n t h e l a t t e r case we l o s e no t r a j e c t o r i e s t o capture and

g e t a b e t t e r p i c t u r e of t h e capture d i s t r i b u t i o n i t s e l f .

-
Certain problems are concerned with the time a p a r t i c l e t a k e s t o

t r a v e l from the source t o i t s death, which c a l l s f o r a parameter T

giving t h e "age" of the p a r t i c l e a t a l l phases of i t s l i f e . A parameter

Y is employed f o r t h e number of c o l l i s i o n s suffered by a p a r t i c l e i n

some problems, e.g., i n those dealing with t h i c k t a r g e t corrections.

We include f o r easy reference a l i s t of these most f r e q u e n t l y

used p a r t i c l e parameters, and w i l l adhere throughout t o the notation

indicated here.

_I

P a r t i c l e Parameters

Space coordinates x,y,z or R

Direction coordinates u,v,w or w

E

Energy group index €3

Zone number

Number of c o l l i s i o n s Y
 3. Remarks on u n i t s . It i s perhaps d e s i r a b l e t o mention b r i e f l y

t h e matter of u n i t s a t t h i s point. We use t h e centimeter-gram-second

systems of u n i t s , with t h e following q u a l i f i c a t i o n s . ( a ) Neutron ener-

g i e s are expressed i n e l e c t r o n v o l t s (ev), Kev (103 e v ) , or Mev(106 ev).

For t h e u n i n i t i a t e d , voltage has dimensions of energy p e r u n i t charge;

one ev i s by d e f i n i t i o n t h e energy acquired by an e l e c t r o n which has

dropped through a p o t e n t i a l d i f f e r e n c e of one ( p r a c t i c a l u n i t ) v o l t .

Since t h e l a t t e r i s lo8/, of t h e e l e c t r o s t a t i c u n i t of voltage, ( 5 )

and t h e charge on t h e e l e c t r o n i s e = 4.8025 x 10-l' esu, it follows

t h a t one ev is eV = 4.8025 x lom1' x lo8/, = 1.60203 x erg.

(b) The energy of a photon i s measured i n (dimensionless) u n i t s of

moc2, where mo i s t h e rest-mass of t h e e l e c t r o n . This is explained

more f u l l y i n Ch. V I .

We include here t h e formula for computation of the time A T in

seconds f o r a neutron of energy E MeV t o t r a v e r s e a distance of d cm.

Letting k - 1.60203 x e r g s per Mev, we have i n t h e n o n - r e l a t i v i s t i c

range of neutron energies ( 6 )

Ek- -1 m v 2
2 1

(5) c= 2.99776 x lolo cm s e c - l , the v e l o c i t y of l i g h t i n vacuo.

(')The r e l a t i v i s t i c k i n e t i c energy i s (m-ml)c2, where m = m l / { z
and B = v/c. A s an e x e r c i s e one may determine the r e l a t i v e error
(1.1%) i n computing v from t h e t w o formulas when E = 14 Mev, the
highest neutron energy involved i n the present r e p o r t , which is thus
confined t o the n o n - r e l a t i v i s t i c range of neutron energies. W e do
give the r e l a t i v i s t i c treatment of the Compton e f f e c t i n Ch. V I ,
however.

-32 -
where ml(gm) i s t h e mass of t h e neutron and v i t s speed i n cm sec-'.

Thus

Now one ''gram atomic mass" of any physical p a r t i c l e contains

Avogadro's. number A = 6.0228 x 1023 p a r t i c l e s . Taking 1.00893 f o r t h e

atomic mass of t h e neutron, we have

v = k' v m cm sec -1

8
where k' = 13.83 x 10 numerically. Thus a 1 Mev neutron t r a v e l s about

14 meters per microsecond ( p sec). We have then

A 7 = k"d/ sec

f o r the t r a n s i t t i m e A T , where k" = 7.231 x 10-l'.

4. Space coordinates f o r source p a r t i c l e s . We have i n d i c a t e d i n

8 1 of the present chapter how t h e machine i s l e d from t h e "start" of

the problem t o the point ( u ) a t which it s e t s up a source p a r t i c l e ,

o r having f i n i s h e d processing a p a r t i c l e , it i s returned t o ( u ) via

(cy ), while i n 9 2 we have discussed the various kinds of parameters

as they a r e c a r r i e d throughout t h e c a r e e r of a p a r t i c l e . We are now

-33 -
ready t o consider how i n i t i a l values are assigned t o these parameters

as a p a r t i c l e issues from t h e source. W

e hope t o give enough examples

t o i n d i c a t e the nature of the procedure.

(a). Uniform source on an annulus o f r a d i i R o < R1. Determina-

t i o n of coordinates x,y ( c f . Fig. 8). The p r o b a b i l i t y density function

p(R) is 27r R/ 7r (R: - 2

Ro), so that t h e fundamental p r i n c i p l e sets

Having thus located the r a d i u s R , we note t h a t p ( 4 ) d4 = dq5 /2n

so t h a t the next random number determines 4 by

W e have therefore the r o u t i n e indicated i n Fig. 8.

v m
(b). Uniform source i n a s p h e r i c a l shell. The f i n a l result i s

R = i f s p h e r i c a l symmetry admits use of t h e s i n g l e

space coordinate R. If x,y,z must be s p e c i f i e d , we s h a l l have

x = Ru, y = Rv, z = Rw, where (u,v,w) i s a point uniformly d i s t r i b u t e d

on t h e u n i t sphere u2 + v2 + w2 = 1. How such points may be obtained

w i l l be discussed i n t h e next s e c t i o n (5a) on d i r e c t i o n parameters

u, v , w, where t h e same problem arises,

(c). Parallel-beam source incident on lateral surface of r i g h t

c i r c u l a r cylinder of r a d i u s R1, height H, o r on a sphere of r a d i u s R

1"

-34-
 Y

I t

Q, = n(2r - 1) x = R cos @

Fig, 8

-35-
If t h e beam i s t a k e n i n the p o s i t i v e y-direction (u = 0 , v = 1, w = O ) ,

and the cylinder i n the position shown i n Fig, 9, one may use the

r o u t i n e there indicated.

For a similar beam i n c i d e n t on a sphere of r a d i u s R we have

1’
x=RII/F; y = - 1-4 = - R~ i-), z = 0 , issuming symmetry

about the y a x i s .

(d), I s o t r o p i c point source a t distance d from r i g h t c i r c u l a r

cylinder. Here the determination of t h e x, y, z coordinates of the

point of e n t r y t o the cylinder i s c o r r e l a t e d with t h e d i r e c t i o n param-

eters u, v , W. We t h e r e f o r e postpone t h i s problem t o t h e next section.

5. Direction coordinates f o r source p a r t i c l e s . e next consider

W

t h e problem of assigning d i r e c t i o n cosines u, v , w, o r simply w in

the s p h e r i c a l l y symmetric case, t o source p a r t i c l e s . If t h e l a t t e r are

l i b e r a t e d i n a material medium, the d i r e c t i o n s may be expected t o be

drawn from a n i s o t r o p i c d i s t r i b u t i o n , whereas p a r t i c l e s emanating from

a surface are n a t u r a l l y l i m i t e d t o h a l f of d i r e c t i o n space and usually

are d i s t r i b u t e d i n some non-isotropic d i s t r i b u t i o n about t h e normal t o

the surface.

( a ) , I s o t r o p i c source; u,v,w d i r e c t i o n cosines. The problem

involved i s tantamount t o t h a t of choosing a point (u,v,w) uniformly

2
d i s t r i b u t e d on t h e u n i t sphere u2 + v2 + w = 1. The element of area
-
f o r t h i s sphere i n s p h e r i c a l coordinates Y , q5 is s i n y dr dq5 =

- dw a+ , where 4 i s the longitude, and we have used y instead of

-36-
 X

- 1)

-
x = R1(2r __c
Y =-

Fig. 9

-37-
 the traditional 0 f o r the remaining angle s i n c e Y has already been

introduced f o r t h i s angle by way of w = cos y ( c f . Fig. 7). The

p r o b a b i l i t y d e n s i t y function p(w) i s t h e r e f o r e given by p(w) dw =

- 2 3 ~s i n y dy/kfl =
I
dw. We may t h e r e f o r e first determine w by

w -
w: -1

and 4 subsequently by

Reference t o Fig. 10 t h e n completes the argumeat f o r t h e i s o t r o p i c

source; we need only remember t h a t P E 1-i =i-.

(b). The cosine d i s t r i b u t i o n . ( 7 ) %is r e f e r s t o a point source

emanating from a surface. If we agree t h a t the outer normal t o t h e

surface a t the point has the d i r e c t i o n u = 0, v = 0, w = 1, then t h e

cosine d i s t r i b u t i o n has by d e f i n i t i o n the p r o b a b i l i t y density function

P(W) = 2w, with w 3 0. Thus

r = 10
W
2w dw

‘7’For a discussion of the way i n which the cosine d i s t r i b u t i o n governs

p a r t i c l e s emanating from a surface one may r e f e r t o F. W. Sears, An
Introduction t o Thermodynamics, the Kinetic Theory of Gases, and
S t a t i s t i c a l Mechanics, Addison-Wesley Publishing Co. Inc Cambridge,, ,
Massachusetts, 1955. Cf. t h e chapter on k i n e t i c theory.

-38-
 t-
I r

u = p cos r#

v 5 p sin Q,

Fig. 10

-39-
and w = 1/F i n t h i s case. The flow diagram of Fig. 10 may be used

with t h e replacement of w = 2r - 1 by w = r/iF.

(c). -I s o t r o p i c and cosine sources i n s p h e r i c a l l y symmetric

systems. I n case s p h e r i c a l symmr?try i n d i c a t e s the use of coordinates

R, w defined i n Fig. 7, we may have t o assign the d i r e c t i o n parameter

w f o r point sources of t h e above k i n d s . Such sources w i l l be located

on s p h e r i c a l surfaces, and i n case of a cosine source, t h e la tte r w i l l

be r e l a t i v e t o the normal t o such a surface, i.e., t o the r a d i u s vector.

It i s therefore c l e a r t h a t the formula w = 2r - 1 is valid for the

i s o t r o p i c case, w h i l e w = a p p l i e s t o outwardly d i r e c t e d , and

w = -Ft o inwardly d i r e c t e d cosine sources.

(d). I s o t r o p i c p o i n t source a t d i s t a n c e ds f r o m c y l i n d e r of

radius R2, height H ( c f . Fig. 11). T h i s problem presents some

f e a t u r e s of interest. Clearly,if we proceed naively t o assign d i r e c t i o n

parameters u, v , w t o p a r t i c l e s as they leave the source s, using t h e

method of ( a ) above, most p a r t i c l e s w i l l f a i l t o h i t t h e cylinder, the

s i z e of which may be greatly exaggerated i n t h e f i g u r e . The physical

problem is n a t u r a l l y concerned w i t h questions r e l a t i v e t o the nuaiber of

incident p a r t i c l e s , not w i t h the problem (an i n t e r e s t i n g one , i n c i d e n t a l l y )

of what s o l i d angle i s subtended a t t h e source by the cylinder. How the

l a t t e r question may be attacked by Monte Carlo we leave as an exercise,

I n t h i s connection it i s worth mentioning t h a t many purely geometrical

problems which present forbidding a n a l y t i c a l d i f f i c u l t i e s are a p t t o

arise i n photon problems, and may be (many have been) successfully

Fig. 11
t r e a t e d by the simple sampling methods which we are discussing.

To r e t u r n t o our point source, we observe f i r s t t h a t t h e only

p a r t i c l e s which can h i t the cylinder are l i m i t e d t o t h e wedge defined

by t h e two tangent planes t o t h e cylinder through t h e source S. This

means t h a t d i r e c t i o n s a r e l i m i t e d , i n u, v, w space a t the source t o

those with longitude #J on t h e range - 4 2 & d, & 42 where 42 -

a r c s i n R2/d,, and + i s e q u i d i s t r i b u t e d on t h i s range, Moreover, for

a11 those p a r t i c l e s i s s u i n g from t h e source between 9 and 4 + d$

the d i r e c t i o n cosine w = cos y i s l i m i t e d by the end points of t h e

element of t h e cylinder determined by 6 , and on t h i s range w has a

constant p r o b a b i l i t y d e n s i t y function, since the element of surface a r e a

on the unit sphere is -dw d d .

Now from the f i g u r e , t h e half-chord c2 -
s o t h a t t h e distance d = d cos 9 - c 2 * Thus f o r t h i s 4 , the range

-1f.
2 s
of w is - <
w2 % w 4 w2, where w2 = ( H / 2 ) / (H/2) The

r o u t i n e for s e t t i n g up x, y, z, u, v , w a t t h e point of e n t r y there-

f o r e appears a s i n Fig. 12, where ,#2, c2>d2, w2 a r e given by the

above formulas. The l a s t box follows from t h e f a c t t h a t

z/d2 = tan(n/2 - Y ) = c o t y = cos y / s i n y = w/pe

(e ). General d i s t r i b u t i o n i n h a l f of direction-space. Occasion-

a l l y it i s necessary t o consider a point source w i t h w & 0 having

some experimentally determined d i s t r i b u t i o n i n w, which may be given

i n the form of a t a b l e . We may then t a b u l a t e Po = 0 < P < PI = 1,

1
w
0
= 1 > w1 > ... > wI = 0, where
'i
i s the p r o b a b i l i t y of a cosine

-42 -
 w = w (2r - 1)

u = p cos 9

Fig. 12

w & wi, and proceed as i n Fig. l 3 a (cf. a l s o formulas ( 2 ) and ( 3 ) ,

~
Ch. I, 85). When various d i s t r i b u t i o n s are t o be tried, it i s preferable
t o run a number of d i f f e r e n t problems, each f o r a s p e c i f i c w; t h e re-

I s u l t s may then be weighted t o give terminal percentages for any desired

source.

(f). A prejudiced source. It is sometimes desirable t o sample

c e r t a i n ranges of source d i r e c t i o n s more thoroughly t h a n o t h e r s . To

i l l u s t r a t e , suppose t h a t we have a n i s o t r o p i c source, with w uniformly
d i s t r i b u t e d on -14 w 4 1, b u t t h a t the position of a counter makes

-43-
 iL

-i 4-1-i

Fig. 13a

t
1-w-w

t f
w-24-w

- Fig. 13b

-44-
more important those p a r t i c l e s which leave the source on the range

-16 w 4 ? < 0. & m y then give source p a r t i c l e s equal likelihood

of s t a r t i n g on t h i s range or i t s complement, provided we assign weights

(? + 1) and (1 - w),hl respectivelJ I n t h i s way about half of a l l

source p a r t i c l e s , of r e l a t i v e l y smaller weight, originate i n the

important cone, the t o t a l weight processed for N p a r t i c l e s having

expectation N(? + 1) -t 2 N(1 - w)
rr)
= N. The procedure i s schematized

i n Fig. l3b.

6. Energy of source p a r t i c l e s , I n the case of energy dependent

I

problems, one usually decides upon a set of energy ranges with lower

bounds E l > E2 > ... > EG, with storage of physical q u a n t i t i e s f o r

each of these ranges, EG
being the lowest energy permitted t o p a r t i c l e s

i n the problem. P a r t i c l e s which by chance reach lower energies a r e

relegated t o a terminal category reserved for such losses. An index

g = 1, ..., G designates the group number. If the source is mono-

energetic, one simply s e t s Eo -+ E, g o 4Q , t h a t i s , E and g are

assigned thevalues of the i n i t i a l energy E, and the index go of the

group i n which t h i s energy f a l l s ; go may or may not be unity; one

may wish t o study the behavior of a s e r i e s of sources of various i n i t i a l

energies.

If the source p a r t i c l e s a r e not monoenergetic but are chosen

from some given energy d i s t r i b u t i o n , one tabulates Po = 0 < .. b

< PG
= 1 together w i t h El > ... > E G and proceeds exactly as indicated

-45-
in 5 5e above, reading g for i, E
g
for w
i'
and E for W.

7. Other source parameters. "he o t h e r parameters mentioned

i n Ch. 11, 52, i f c a l l e d f o r , have t h e following values a t the source:

(4 3 i s set equal t o t h e number of the zone i n t o which source

p a r t i c l e s enter (which may depend on t h e values assigned t o s p a t i a l a n d

d i r e c t i o n a l coordinates, and t h u s involve a decision routine ); (b ) the

weight W is usually unity a t the source; ( c ) age t = 0 initially;

and ( d ) Y = 0 f o r t h e number of c o l l i s i o n s undergone.

8. Source f o r a -type c a l c u l a t i o n s . I n a n important c l a s s of

problems, including those which a r e designed t o determine t h e rate of

growth of a neutron population i n a f i s s i o n a b l e m a t e r i a l , t h e source

c o n s i s t s of a s p e c i f i e d nuniber of neutrons N i n group g,

@;,pi
and d i r e c t i o n range e w i l l i l l u s t r a t e f o r the case of a
i. W
zone P
homogeneous sphere of r a d i u s RZ, where these ranges are determined by
bounds

Eo>E1> ..e
> EG

0 = R o < R1< ... < R Z

1 = wo > w17 ... > WI = -1

-46-
 The o b j e c t of such a problem is t h e determination of t h e d i s t r i -

bution N' i n which t h i s i n i t i a l d i s t r i b u t i o n r e s u l t s a t time

g ,3,i
A 7 , and then t o use the output N' as the input source d i s t r i -
63 $3,i
bution for the next cycle. A s u i t a b l e modification of the, (a) and

(0) r o u t i n e s is i n d i c a t e d i n Fig, 14. The determination of E, R,

and w on t h e i r ranges may be achieved by appropriate i n t e r p o l a t i o n ,

o r , i f preferred, by d e t e r m i n i s t i c assignment of s u i t a b l e mean values

The i n i t i a l values of are a r b i t r a r y i n a - c a l c u l a t i o n s , but

N
Q i ,a,
are guessed as well a s possible t o hasten convergence t o the limit

distribution, The
-g ,-p 3 refer t o t h e category being processed, and
mu@t be distinguished from t h e neutron parameters g 9 3 , i, which are

s u b j e c t t o change while t h e category from which the neutrons arose i s

being processed. The feedback of N' f o r N i s s u b j e c t t o s u i t a b l e

renormalization t o preserve a reasonable source s i z e

-47-
 I" /+ Det. E on ( E-, E- 1

I 1

(Cf. Ch.-IV, 09)

-
g+1-g 1 - F 1--i

Print all NI -
g,g,i
- I

Fig. 14

-48-
 C W T E R I11

THE MEAN FREE PATH AND !TRANSMISSION

1. The cross s e c t i o n concept. The cross s e c t i o n u of a

s t a t i o n a r y " t a r g e t " p a r t i c l e f o r p a r t i c l e s of energy E, r e l a t i v e t o

a given s i n g l e process may be thought of as t h e area presented by t h e

t a r g e t , assumed s t a t i o n a r y i n t h e l a b o r a t o r y system, t o a beam of '

( p o i n t ) p a r t i c l e s of t h i s energy, r e l a t i v e t o the l a b o r a t o r y system.

If we regard a t h i n slab of material of area a , thickness dl , numer-

i c a l d e n s i t y N ( t a r g e t p a r t i c l e s p e r cm3 ), t r a v e r s e d by a p a r a l l e l

beam of p a r t i c l e s (of energy E) normal t o a t h e t o t a l area presented

t o t h e beam by t a r g e t p a r t i c l e s i s u N a ( d j ) , assuming ai so s m a l l

t h a t no "shadowing" exists. The f r a c t i o n of beam-particles undergoing

t h e process i n t h i s volume should be 0 N a (d k ) / a . We, t h e r e f o r e ,

have the a t t e n u a t i o n l a w

dn = -n N o ( d b )

f o r t h e number n of p a r t i c l e s i n t h e beam, and

n = no exp ( -NuR)

-49,
represents the number of particles remaining in the beam after traversing

a distance k in such a uniform infinite medium, no being the number of

particles in the beam at 1= 0. Note that we are assuming no other

competing processes exist.

It is, therefore, suppoeed that

is the probability for a first collision between 1 and k + d k , and

P(a) = -
e
'
/ (-Nos) Nu(d1) = 1 - exp (-Not)
0

is the corresponding probability distribution function for a first

collision at distance 51.

2. The mean free path. The average distance to first collision

is defined as the first moment of the function p(I), i.e.,

---
and is called the mean free path for the process at this energy.
It follows that the Monte Carlo determination of distance 1 from
a n arbitrary point of departure to first collision, assuming the medium

-50-
 homogeneous and infinite must be

r = ~(9.) = 1 - exp (-RIA)

or

A= - A R ~(1 - r>

Since 1 - -
r is equidistributed on 0 5 r &l if r is, we may use

simply

on ,
The 6um of the cross sections 4
A (el.) + A
u (in.) + ... of all
types Tor a particular nucleus A is called its total cross section -
Q
A
(tot,). If the medium contains nuclei of types A, B, C, ... in
numerical densities NA, NB, Nc9 ..., respectively, the "total cross
section" for the medium is defined to be

-51-

i
 A
Z = NAu ( t o t . ) + B
NB u ( t o t . ) + ...
Reference t o t h e preceding discussion makes i t c l e a r t h a t t h e

t h e simple N u of t h a t argument should be replaced by C i n t h e general

case, and the mean free path f o r t h e medium i s , t h e r e f o r e ,

The Monte Carlo method i s always concerned with t h e d i s t a n c e 1

from point of departure t o c o l l i s i o n , and only i n case of c o l l i s i o n ,

t u r n s t o a consideration of t h e nature of t a r g e t h i t , and t h e type of

process involved. Thus it i s always the mean free path A = 1/ I; that

i s used and never t h e free path f o r any of t h e individual processes.

It must be remembered t h a t cross s e c t i o n s are, i n general, de-

pendent upon t h e energy of t h e p a r t i c l e (picturesquely, the s i z e . of t h e

t a r g e t depends on the speed of the arrow) ; thus we write 0

A
E(el.), ...
A
TE, and h E as functions of energy E, and, i n p r a c t i c e , t h e
0
E( t o t . ) ,
free path i s u s u a l l y tabulated as h g , g = 1, ..., G, where g is t h e

energy group index.

Moreover, i n systems c o n s i s t i n g of zones of d i f f e r i n g composition,

we w i l l have an a d d i t i o n a l zone index on t h e free path, thus A 3

@;*

-52-
 3. An example. Consider the problem of determining A for neu-
trons of energy E = 3 MeV in a medium of CH2 of density d = ,92 gm cm-3 .
C
At this energy, C has an elastic scattering cross section (el.) = 1.14
Q

barns (1 barn = lom2‘ ern2 ) while H has a similar cross section H (el.) Q

= 2,23 barns. No other processes are involved.

The atomic weight of C is 12 and of H is 1, so that the molecular

weight of CH2 is 12 + 2(1) = 14. In one gram molecular weight G of any

24 molecules. The mass of one CH2 molecule is,
compound are A = .6 x LO
therefore, G/A gm. Since 1 cm3 of CH2 has mass d gm, the number of
3
molecules of CH2 in 1 cm is N = d / ( G / A ) = dA/G = .0394 x 1024 cm-3
Hence, the numerical densities of C and H are NC = N, NB = 24 and
c = N~ uC (tot.) + N~ 0 H(tot.) = N(C C (el.) + 2 0H(el.)) = 221 cm-1 .
Thus finally A = 1/ 2. = 4.52 cm.
In problems involving many zones of the same material at different
densities, it may be necessary t o store only the basic nuclear constants,
together with zone densities, and provide for the machine to compute its

own h’when needed, by reference to the stored quantities. For instance,

g
in the preceding example, taking energy dependence into account, we should
have

where

-53-
and d are stored quantities. Such a procedure is time-consuming,
2
especially when it necessitates computation of the probabilities for

type of collision upon each collision.

The concept of free path as we have formulated it applies to
photons in their interaction with electrons and nuclei as well as to
neutrons interacting with nuclei. A discussion of photons will be

found in Chapter VI. Throughout the report, except in Chapters V and

VI, we speak in a general way of "particles" which may be photons or
neutrons.

4. "Small" systems and transmission. Consider a homogeneous

medium having mean free path h for particles of a monoenergetic
source. If the distance from source S to the boundary of the system

along a given direction is L, then, of a beam of N source particles

leaving the source in this direction, N exp ( - L / h ) will escape un-

deterred. It is clear that if the dimensions of such a system are

small compared to the free path A, most source particles will escape.
This is' especially undesirable if the assignment of source parameters

is complicated, and,in any case,requires needlessly large sources t o

produce an effective sample.
Now there is nothing to prevent us from regarding a single
mathematical particle leaving the source in a given direction as

representing a large set of W actual particles. Since the output

of a l l problems is a set of ratios Ni/N, where N is the total number

-54-
of source p a r t i c l e s , t h e i n i t i a l va,lue assigned t o W may be taken as

unity.

We may then argue t h a t , f o r a p a r t i c l e of weight W (=1)leaving

t h e source i n t h e s i t u a t i o n described above, a p a r t i a l weight W exp (-L/A)

i s transmitted, t h e l a t t e r weight being tallied i n a category T reserved

f o r t o t a l transmission (without c o l l i s i o n ) . Then we determine a p o s i t i o n

of f i r s t c o l l i s i o n on t h e i n t e r v a l 0 - d L- L within t h e medium f o r t h e

remaining p a r t i c l e of weight W (1 - exp ( - L / h ) ) according t o the formula

r = P((L)/P(L), where P(1) = 1 - exp (-!/A), as derived i n a preceding

section. Solving f o r 1, one o b t a i n s 1 = -A an(1 - r [l - exp (-L/A)]}.

If t h e medium is non-homogeneous, b u t may be regarded as c o n s i s t i n g

of a number of zones, each homogeneous i n i t s e l f , t h e portion of t h e l i n e

of f l i g h t l y i n g within t h e medium may be decomposed i n t o successive i n -

t e r v a l s of lengths LI, .*,, Lm, where L i s t h e segment l y i n g i n zone 7

3
w i t h free path A The transmission i s then c l e a r 1 y . t = exp (-Ll/Al)
7.
exp ( -L2/h2) . exp ( -Lm/Am) = exp - {Ll/AL + ... + L,/hm } . The

p r o b a b i l i t y of a f i r s t c o l l i s i o n at a d i s t a n c e 2k L from t h e p o i n t of

o r i g i n is, t h e r e f o r e , 1 - exp (-P), where P i s defined by

L1 + ..*e L5 .c -
LLl + .*. +L
3-1
+L
7
and

-55-
 Thus p is the number of free paths represented by 1. We
therefore may record the weight Wt as transmitted, and force a first

collision of the weight W ( l - t) on the line of flight at a distance

1<- L by means
&
of the formula

r = (1 - e
-P
)/(1 - t)
Thus P = - n
!
, [1 - r(1 - t)] determines 7 by means of the inequalities

- L
L1
+ ..a +=< p = - L1
1
+ ...
5 r-1

and 1 by the equation

5. The "forced first collision" routine. We illustrate the use

of the device in t w o t y p i c a l problems in this and the following section.

Consider a parallel-beam monoenergetic source incident on the lateral

surface of a cylindrical shell of radii R o < R1, and height H, composed
of homogeneous material having free path h at this energy. We suppose

that the source routine has already assigned to a source particle its
parameters at entry, say, u = 1, v = 0, w = 0, and x, y, z (in the manner
indicated in Chapter 11, 64c), together with E = Eo, g = go, II = 0, and

-56-
W = 1. The e x i t from ( 0 ) should then l e a d t o the ( B o ) r o u t i n e for

forced f i r s t c o l l i s i o n as indicated i n Fig. 16, based on t h e geometric

p r o p e r t i e s of Fig. 1.5.
J

c--)-ecL
I
I
S lo -x
___c 1
I
I
I
‘S-CI-l
c--$----c
I

- x

Fig. 16

-57 -
 The procedure naturally depends on whether or not the line of

flight crosses the hole. This is the reason for the y I I - R0 decision.
The distance L is the total path length of the line of flight lying
within the medium. The transmission T~ is based on the free path for
initial energy group g0 . The distance traversed within the medium from
point of entry to point of collision is denoted by 4. Note that the
exit leads to the collision routine (Y) with all parameters as they
exist at the point of collision, momentarily before impact.

6. Remark on the device in spherical problems. Although we

shall not include an example of the forced first collision device in
a system treated with spherical coordinates, we should mention that,
if it is used to determine a weight W and a radial distance R
a at the
point of first collision,(8) one should exit to an entry such as (7')

discussed in Chapter IV, 83, which sets up the direction w and

radius R as they exist at the point of collision, before entering
- - - A

(Y) itself.

A s an example of the forced first collision method in a non-

homogeneous medium, consider a parallel beam of particles directed

vertically upward (u = 0, v = 0, w = 1) and incident on a sphere
2 2
x2 + y2 + z = Ro subdivided into spherical homogeneous shells by

the radii Ro > R1 > ... > Rk = 0, zone 7 having total mean free path
I8)If R,w are the source parameters, and k the distance from source
to forced first collision, R2 = R2 + k2 + 2Rkw in the solid homo-
geneous case. B

-58-
 Fig, 16

A
5
f o r the i n c i d e n t energy, A set of s t o r a g e places PI, ..., pk are

reserved i n the machine f o r t h e numbers of free paths represented by

segments ( z ' , ~ " ) of the l i n e of f l i g h t i n t h e zones through which it

passes. The method i s i l l u s t r a t e d i n Fig. 16a. It i s assumed that

x = R o 6 , y = 0 have already been set f o r the p o i n t of e n t r y i n the

source r o u t i n e , symmetry obtaining about the z-axis.

-59-
0

d Ap = z"/A

1
1

T-
Wt+T-T t = exp (-2P)
AP -
P+Ap-P
pa-

w
w-wt-w

fi 9

Fig. 16a

-60-
 7. The transmission in subsequent history. The device of
forcing collisions of the non-transmitted weight may, of couxse, be
applied to collisions after the first. If applied consistently to
-
all collisions, then no trajectory ever terminates in escape, and one

would ordinarily rely on a weight cutoff for termination. Usually the

geometric complexity of paths after first collision renders use of the

device for further collisions impractical. We do not consider forced

collisions other than the first in the present manual.

8. Prejudiced first collision in "large" systems, In problems

concerned with "large" systems, such as those arising in shielding,

the transmission must be very small and yet one may have to obtain
energy-angle distribution of escape and space distribution of various
types of collision (e,g,, inelastic collision and radiative capture of
neutrons in determining Y -sources) throughout the system. Moreover,
the existence of energy cutoffs makes very unlikely the.arriva1 of
particles in the farther reaches of the system after many collisions.
In such cases, one may overcome the dwindling of first collisions due
to the exponential by prejudicing the distribution of first collisions
and weighting accordingly.
We illustrate with a simplified example. Consider a thick plane

slab of two layers (free paths Al, A2 at incident energy) bounded by

the planes z = 0, z = $ z = L2, with a source directed vertically

upward and incident on the surface z = 0. One may then determine the

-61-
position z of first collision, together with the weight at this point
by the scheme of Fig. 16b. Note that the formula z = rL2 distributes
first collisions uniformly throughout the slab, while the weighting

gives the correct expectation for first collision between z and z + dz,
namely,

-dz. -L2 e- P -
= e- P d z
L2 Am 4n

If N source particles are processed, the expected total weight

assigned to these at first collision is then N(l - where

T~),
To

=.,-{>+ 1
L2 - L1}
h2
is the transmission.

t
I I

b
Fig. 16b

-62-
 CHAPTER IV

THE COLLISION OR ESCAPE ROUTINE

1. Introduction. I n the last chapter, we discussed a s p e c i a l

r o u t i n e f o r "forced first c o l l i s i o n . " Conceptually, t h i s r o u t i n e con-

d u c t s a particle from a special p o i n t of departure, namely, t h e source,

t o a p o i n t of c o l l i s i o n within t h e system, with no escape a l t e r n a t i v e

f o r the t r a j e c t o r y i t s e l f . When t h i s ( B o ) r o u t i n e i s used, it i s

entered only once, namely, d i r e c t l y from the source, and has only one

e x i t , the collision routine ( 7 ) .

I n a l l problems, whether t h i s device is used o r not, a r o u t i n e

( e ) is required which is designed t o conduct a particle from a p e r f e c t l y

a r b i t r a r y p o i n t of departure i n an arbitrary zone, a t which the p a r t i c l e

parameters are known, t o i t s next p o i n t of c o l l i s i o n within t h e zone,

o r t o i t s p o i n t of departure from the zone, i n the event the boundary

i s reached without c o l l i s i o n . (Cf,, however, Bg of t h i s chapter.)

T h i s (p) r o u t i n e is e n t e r e d d i r e c t l y from the source, i n prob-

lems not using t h e (So) device, and i n various o t h e r s i t u a t i o n s , namely,

as a p q r t i c l e departs from a c o l l i s i o n , as i t e n t e r s a new zone, o r

-63-
upon re-entry into the same (central) zone, after crossing a central
hole, etc.
Aside from the initial determination of the distance k to

collision (calculated for an infinite homogeneous medium), this routine

is purely geometric, and consists essentially in comparing the distance

1 with the distance d to the boundary of the zone along the line of

flight. If 1\ <d, space and direction coordinates are set at the point
of collision, and one proceeds to ( ?’ ). If k - the particle is con-
sidered to reach the zone boundary. I n this case its parameters are

set at the boundary point, and one returns to ( B ).

It is clear that when a problem involves zones of different
geometric shapes, several such routines ( B A ) , ,&= 1,2,. .. may be
required, each designed for the geometric problem of its own type of
sector. In such cases the particle carries an additional parameter
which indicates the type of geometric zone it occupies at any given
time. Transfer is made from source points, points of collision, and
so on to the appropriate (fib).

One m y also note that, under our procedure, a particle reaching

the boundary of a zone is referred back to ( B ), unless escape from the
system is involved, and is then treated anew relative to the zone
entered. An alternative procedure using a single random number to
determine the eventual position of collision, or escape from the sptem,
by reference to all path segments defined on the line of flight by a l l

zone boundaries may be used but seem clumsier t o handle, and we do not

-64 -
d i s c u s s it, The method involved should be c l e a r from the discussion

a t the end of $4, Chapter 111.

2. A routine f o r the spherical shell. Consider t h e problem

of a p a r t i c l e with parameters R, w, E, g,$ a t a poSnt of departure i n

3-1< R8'
a s h e l l of r a d i i R The flow diagram of Fig. 17 w i l l be seen

t o keep computation a t a minimum. Here, Ra is t h e r a d i a l d i s t a n c e t o

t h e p o i n t of c o l l i s S a n i n a n i n f i n i t e medium of free path A' t i s the
g'
'tangential d i s t a n c e t o t h e inner boundary from t h e point of departure

a t R , and wt = -t/R is the cosine 00 t h e angle Yt from OR t o t h i s

tangent (cf. F i g . 18) Reference t o the cos Y curve of Fig. 19 shows
that, w being negative, wt2 - w2< 0 implies tha%the l i n e of f l i g h t

c u t s the inner bowdary. Moreover, if Fti2 - R2%-I is non-negative, %he

sign o f t2- t 2 d i s t i n g u i s h e s between a p o i n t of c o l l i s i o n on one o r

t h e o t h e r side of t h e i n n e r boundary.

It w i l l be seen that f o r " s o l i d spheres" with e o c e n t r a l hole

t h i s flow diagram w i l l work automatically i f a = 1 is the index of the
innermost zone, provided Ro = 0 be s t o r e d together with t h e o t h e r zone

radii R1 < .. < Rzo Thus no s p e c i a l ( B ) r o u t i n e i s required f o r t h e

c e n t r a l s p h e r i c a l zone. The exlts ( r ' ) indicate c o l l i s i o n within zone z

(cf. the next s e c t i o n ) .

-65-
 2 2
a t
w -w
0

Fig. 17

-66-
 0 R R
a -l
Fig. 18

w = cos y

1.

-1 --

Fig, 19

-67-
 3. Reorientation formulas for the spherical shell. ( 9 ) In case

of collision (Y'), one provides a routine for computing the parameters

R,w as they exist at the point of collision, before proceeding to (Y)

itself. This is indicated in Fig, 20. Moreover, in the event of escape

at R we have two alternatives. If 3 = 2 , the number of the outer-

5'
most zone, the particle leaves the system, and one proceeds to (e),

while if 7 ~ 3the
, particle enters a new zone 3-k 1. Classification
of an escaping particle may or may not involve its direction of escape.
In the latter case, one may by-pass the reorientation part of the escape
routine of Fig, 20 by putting the ?+ l++#and 2 - 7 boxes first. One

=i[l
obtains the w ' formula of the latter routine from the relation w'

= cos 71 - sin ?'I] and the law of sines: sin Y / R

7
= sin yl/R.

Note that cos Y ' takes the positive square root since an entry to a

zone from the inner boundary has acute angle of entry.

Similarly, one obtains w ' = -f[ 1 - (R/Rt-1)2(1 - w2)] for the
entry cosine at %-le For solid sphere problems, the escape contingency

at R occurs - , involves the substitutions w'+w,

only if ~ > 1 and
7-1
using the latter formula f o r w', R + R, 7 - l + ~and
, thence passage
8
to ($1.
If the central zone y = 1 is vacuum, escape at R1 from zone y = 2
necessitates re-entry of zone
a= 2 after crossing the hole. Thus for
the hollow sphere problem, we have the routine of Fig. 21.

'"All the geometric relations involved in the present section are

indicated in Fig. 22 and F i g , 23.

-68-
Fig. 20
 t

Fig. 21

-70-
Fig. 22

Fig. 23

-71-
 4. Flux problems i n s p h e r i c a l geometry. C e r t a i n problems, one

of which i s t h a t of a c e n t r a l neutron source i n homogeneous air, r e q u i r e

study of the number of neutrons of energy E which t r a v e r s e (imaginary)

s p h e r i c a l s u r f a c e s a t varying d i s t a n c e s from t h e source. Such s i t u a t i o n s

may be handled by subdivision of t h e medium i n t o spherical zones, and

incorporating i n t o t h e "Escape a t R ?and R7-1 I' r o u t i n e s a cumulative

t a l l y i n counters N r,g. -
These are not terminal c a t e g o r i e s , except f o r

?=?and do not e n t e r i n t o a "sum check," b u t count all neutrons of

a l l energy groups whenever they cross a spherical boundary

B.
5. A r o u t i n e f o r t h e f i n i t e cylinder. We consider now the case
af a p a r t i c l e w i t h parameters x, y, 2, u, v, w, E, g a t some p o i n t of

departure w i t h i n a f i n i t e homogeneous cyiinder of r a d i u s R1, height H

( c f . F i g . 24). The procedure i s i n d i c a t e d i n Fig. 25. Note t h a t before

e n t e r i n g ( Y ) i n the event of c o l l i s i o n , a l l parameters are s t o r e d i n t h e

machine as t h e y o b t a i n a t the p o i n t of c o l l i s i o n . The d i r e c t i o n coor-

d i n a t e s u, v, w are the same

. ,
a t the e n t r y t o ( Y ) as they were at the ( 0 )
e n t r y , and do not require new evaluation, as w a s the case f o r w i n

s p h e r i c a l geometry. The e x i t (e) denotes, as usual, escape from the

system and leads t o a c l a s s i f i c a t i o n r o u t i n e f o r escaping particles.

Such r o u t i n e s are discussed i n a later chapter.

-72-
Fig. 24

Fig. 25

-73-
 6. The f i n i t e c y l i n d r i c a l s h e l l with c e n t r a l hole. Let

x, y, z, u, v, w, E, g be t h e p a r t i c l e parameters as they obtain a t

a n a r b i t r a r y point of departure within a f i n i t e c y l i n d r i c a l s h e l l of

r a d i i Ro< Rl and height H. The procedure of 85 may be used unless

t h e l i n e of f l i g h t x' = x + ut, y t = y + vt, z' = z + wt, t20

c u t s t h e ( i n f i n i t e ) inner c y l i n d r i c a l surface x2 +'y 2 = Ro2 a t two

rea1,distinct distances

t = -I* V A O
2
1 - w

where 6 = ux + vy and A = I2 - (1 - w2 )(x2 + y2 - Ro 2 ). I n the

0

l a t t e r case one provides a n a d d i t i o n a l r o u t i n e as indicated i n Fig. 26.

Study of t h e flow diagram together with Fig. 27 should make

t h e method c l e a r . Note that t h e case w = 1 and t h e case of a l i n e

of flight tangent t o t h e i n n e r c y l i n d e r are handled automatically.

Nevertheless, some caution may be required a t e n t r y t o t h e "t"boxes

2
if 1 - w i s very small. Such a case may be handled e a s i l y by a pre-

liminary comparison of [ - I -f: &] - K(l - w2) with 0, where K is

a constant l a r g e r than t h e Largest d i s t a n c e within t h e s h e l l . I n case

t h e latter d i f f e r e n c e i s p o s i t i v e , one may consider 1

zt( as e s s e n t i a l l y
i n f i n i t e and may r o u t e t h e flow t o t h e proper box d i r e c t l y , by-passing

t h e e x p l i c i t computation of t and zt.

An alternative method which i s very convenient when s e v e r a l

c y l i n d r i c a l zones are involved i s that of Fig. 27a. This has t h e

-74-
 .

-6 +

1-w
m0
1
+ 2
-6 - mo
1 - w

= z + w t
5

-
(line of flight does not cut inner cylinder)

=z+wa

T XA =

y =y+vJ
x +u l

I x
y
i
.
i. v
ul- x
ZR -
l- y
z
I
(Escape through
hole bases)
-0

1 -I
""zi
@ (Re-entry across hole)

Fig. 26

-75-
 Z
I

Fig. 27 I
advantage of effecting the decision between escape or collieion with a
minimum of square roots, sets up the space parameters at escape position
with little repetition of code, and avoids the "infinity decision" re-

ferred to above, It i s understood,that the inner and outer radii

0 2 Ro< R1 and base plane z-coordinates H' and H" are properly set as
a particle enters one of the system of zones. Moreover, a parameter
2
R2 5 x + y2 is carried throughout with x, y, z, and an additional
Pi! f 1 - w2 together with u, v, W. An index X is also used, being
set equal to zero at the source. The latter is of a purely computa-

tional character.

-76-
 r--i
t = 1-6 + q rnlP A = 62 - p2W2- q) I

@
t
I xt=xtut
I
Y t = Y + d
2 = z + w t
t
2 2 2
R t = x + y
t t

Fig. 27a

-77-
 I t may be h e l p f u l t o note t h a t q = -1 or q = +1 according t o

whether the d i r e c t e d l i n e of f l i g h t does o r does not c u t t h e i n n e r ( i n -

f i n i t e ) surface, while d i s t t p l u s t to r "minus" according t o whether t h e

p o i n t of c o l l i s i o n i s outside o r inside t h e ( i n f i n i t e ) base planes.

2
After f i r s t t r a n s i t through t h e box i n which xt, yt, zt, R
t are corn-
puted, these q u a n t i t i e s refer t o t h e p o i n t of i n t e r s e c t i o n of t h e l i n e

of f l i g h t with t h e base plane z = H, when d 2 0, and t o t h e p o i n t of

c o l l i s i o n when d < 0. The e x i t s #? and /3 i n d i c a t e escape from t h e

R H
zone through the i n n e r o r o u t e r lateral surface, and from one of t h e

bases, r e s p e c t i v e l y . I n a n a c t u a l problem involving s e v e r a l s e c t o r s ,

such e x i t s must lead t o r a t h e r involved r o u t i n e s for deciding t h e s e c t o r

next entered and s e t t i n g up i t s geometric parameters Rf, HI, H", and

T o r c l a s s i f y i n g escape i f such i s the case.

7. The s p h e r i c a l s h e l l i n absolute space. Even when t h e medium

i s s p h e r i c a l l y symmetric, i t may be necessary t o keep t r a c k of t h e d i r e c -

t i o n of a p a r t i c l e i n absolute space, e.g., i n case emergent p a r t i c l e s

are t o be c l a s s i f i e d with r e s p e c t t o t h e i r d i r e c t i o n s relative t o a

given source d i r e c t i o n . I n such cases, i t is convenient t o use x, y,

z, u, v, w parameters. The method i s similar t o t h a t of t h e preceding

s e c t i o n and should r e q u i r e no f u r t h e r explanation. The procedure f o l -

lowing a "core h i t " depends on the problem and may involve an absorption,

passage t o an inner zone, or crossing of a c e n t r a l vacuum. W

e do not
e l a b o r a t e t h i s case f u r t h e r . The procedure i s indicated t o some extent

i n Fig. 28.

-78-
I A. = 62 - (x2 + y2 + f' + Rt I

Fig. 28

-79-
 8. Slab geometrx. I n problems on plane slabs, t h e medium can

be considered t o occupy t h e region of x, y, z space defined by Zo <=

&i zT9
t h e coordinate z and t h e cosine w of the angle Y which the l i n e of

f l i g h t makes w i t h t h e p o s i t i v e z-axis being t h e only r e l e v a n t coordinates.

The slab may c o n s i s t of zones of d i f f e r e n t kinds of media, with upper

boundaries defined by Z1 < ... < Z-.

P
The ( 8 ) r o u t i n e f o r such a problem

may be l i k e t h a t i n Fig. 29.

bt
z2-1- z - 8-1 - 2 -

-4
0

Fig. 29

-80-
 9. Problems run i n cycles of time A T . I n problems involving

t h e d i s t r i b u t i o n of p a r t i c l e s ( i n energy, p o s i t i o n , and d i r e c t i o n ) a t a

given t i m e A T from t h e i r o r i g i n i n some source d i s t r i b u t i o n , a p a r t i -

c l e i s followed u n t i l it i s l o s t t o some terminal category such as e s -

cape, capture, e t c . , ox u n t i l t h e a l l o t t e d time 4 7 has expired. Such

problems arise n a t u r a l l y i n a -determinations ( c f . Chapter 11, 5 8 ) , t h e

computation being performed i n successive cycles, t h e output d i s t r i b u -

t i o n of each cycle being used as t h e input of t h e next, The ( B ) routine

i n such cases must be modified s o t h a t t h e e s s e n t i a l decision between

c o l l i s i o n within t h e zone or a r r i v a l a t t h e zone boundary i s contingent

upon the q u a l i f y i n g condition: "if t i m e permits." Thus, t h e distance

d t o c o l l i s i o n o r boundary must be computed, and t h e t i m e ( c f , Chapter

11, § 3 )
TI = IT f k"d/ fi
of t h e s e events i s then compared with cycle time AT. If T'< AT one

substitutes T'+T and proceeds as usual, whereas, i f T' 2 A T , time

runs o u t before t h e event can occur, One then computes the distance

at = (A T - T )k1 fi
t h a t t h e p a r t i c l e can t r a v e l i n t h e time remaining, and t h e p o s i t i o n

and d i r e c t i o n a f t e r t h i s distance i s traversed. The p a r t i c l e i s then

c l a s s i f i e d i n an energy, position, and d i r e c t i o n category N g , p and

one r e t u r n s t o (a). We include i n Fig. 30 an example f o r a homogeneous

sphere with radius R- See §@,3 of the present chapter for (y') and
3..
Chapter 11, 58, for (a). (Lg r e f e r s t o l o s s from t h e boundary.)

-81-
 T' -AT

Fig. 30

-82-
 THE COLLISION ROUTINE FOR NEUTRONS

1, Introduction. We have seen i n t h e last chapter how t h e

geometry of t h e system determines the immediate fate of a particle, on

t h e basis of the equation g= - A l n r, as a c o l l i s i o n o r a n escape a t

t h e boundary of the zone. The present chapter is devoted t o the methods

involved i n d e a l i n g w i t h t h e former contingency i n the case of neutrons

c o l l i d i n g w i t h nuclei. I n t h e following chapter, we consider t h e c o l l i s i o n

r o u t i n e s f o r photons.

It may be t h a t t h e media occupying d i f f e r e n t zones T = 1, 2, ...,?

are so d i v e r s e i n t h e types of n u c l e i contained, and thus i n t h e types of

neutron processes involved, t h a t i t i s not worth w h i l e t o attempt a gen-

eral code f o r ( Y) covering a l l contingencies. I n such cases, d i f f e r e n t

c o l l i s i o n routines ( y.$ may be provided, each economically adapted t o

i t s own type of medium.

The basic o b j e c t of t h i s chapter and t h e next i s t o show how,

after c o l l i s i o n of a neutron o r photon i n a given medium, t h e new E, g,

u, and W are obtained, as w e l l as t h e cosine of t h e l a b o r a t o r y angle of

-03-
d e f l e c t i o n from the i n c i d e n t l i n e of f l i g h t . The e x i t (a) from t h e

(7) r o u t i n e refers t o a purely geometric procedure which determines

t h e new d i r e c t i o n parameters u, v, w ( o r simply w) from t h e incident

d i r e c t i o n parameters and t h e laboratory d e f l e c t i o n cosine, and thence

l e a d s back t o ( B ) , "he (6) r o u t i n e i s developed i n Ch. VI1

2. Capture and s e l e c t i o n of t h e type of c o l l i s i o n . It i s impos-

sible t o give a p e r f e c t l y general procedure f o r i h e (7) routine, so

d i v e r s e are t h e various types of processes. Each problem must be studied

i n i t s own i n d i v i d u a l i t y .

If t h e number v of c o l l i s i o n s i s among t h e neutron parameters, we

may begin with

Unless the m e di um i n question c o n s i s t s of only one type of nucleus,

and t h a t nucleus has only one type of cross s e c t i o n f o r neutrons, we

must next proceed t o decide on t h e type of nucleus h i t , and t h e type o f

collision. Recalling t h e discussion of cross s e c t i o n s (Ch. 111, f, 1,2,3), C

it i s c l e a r t h a t t h e t o t a l a r e a presented t o a beam of neutrons of energy

group g i n t h e t h i n slab t h e r e defined i s (N

A g
A (tot. ) +
u ... a,.
) 01

Hence, assuming a c o l l i s i o n i n the medium, t h e p r o b a b i l i t y t h a t the

c o l l i s i o n be with a nucleus of type A i s t h e r a t i o of areas:

-84-
where the a d k may be cancelled.
If one o r more of the nuclei A, '8, C, ... present admits
capture, or some other process which may be regarded as terminal for

purposes of the problem (e.g., inelastic collision in monoenergetic

problems), we may decide to use a weight parameter W, as discussed

previously. (Cf. Ch. 11, §2.) It is then economical to store deter-
'y

ministic fractions C of weight captured on collision, where

Q

-c
ET
= [N
A
0A (cap.) + ..*]/[N
g
o A (tot.)
A 6 3
+ ... J

the summations being over nuclear types A, B, ... . We should then

include at the outset the routine of Fig. 31. Here Lc is a "terminal"
category in the sense that capture is a final event in the life of a

physical neutron. However, the use of weights prevents the loss of the
geometric path being followed and greatly improves the statistics.

When weights are used, it is advisable to use a "weight cutoff" Wo,

below which weights are negligible, and an additional terminal category
to catch weights falling below the cutoff. This category is ter-
=w ,
minal in the mathematical sense that one returns to (a) f o r a new
source neutron in case of such a loss, the trajectory terminating at

this point. If the weight W of the uncaptured beam exceeds Wo, we

-85-
proceed to decide, if necessary, upon which type of nucleus A, B, C, ...
the scattering takes place, probabilities being notr dependent on the

assumption that a non-capture collision occws.

Suppose for simplicity only two nuclear types A and B are

present. For machine purposes, we assign to type A a "type parmeter"
value e = 1 and to B the value e = 2, We store in addition to the
ry

above C the probabilities

!3
A

Fig. 31

Fig. 32

-86-
 A = NA ( 0A (tot.) - uA (cap.))/[lTA (0: (tot.)
g Ei;

for non-capture collision on type A, and enter the decision routine of

Fig. 32 for determination of the parameter e.
We then proceed to routines designed to decide the type of colli-

sion on A (e = 1) or on B (e = 2). If the processes involved in the

two cases are sufficiently similar one m y join the exits of Fig. 32 as
shown and go to a comon routine which, by use of the variable e, can

handle both cases. Otherwise,one provides separate routines for the two
types. In any case there is at this point a real disjunction, the neutron
hitting one or the other type of nucleus.
It remains to decide which type of collision is undergone, assuming

it to be non-capture on nuclear type e, This involves reference of a

random number to an additional set of probabilities, the number of which
w i l l depend on the number of kinds of processes other than capture which

nuclear type e admits.

For example, if e = 1 and A is uranium, we might have three

possibilities: elastic scattering, inelastic scattering, i.e. , (n-n)

reaction with loss of energy, and fission. We should have to store in

this case the probabilities

-
-07
and proceed from the box =]I of Fig. 32 to the flow dLagram of Fig. 33.
It m y be preferred, even if capture is to be treated by weights,

to first decide on the type of nucleus hit, by reference to the total

probabilities

and then to store individual capture fractions

A A
0 (cap.) / u g (tot. )
Q

f o r each type of nucleus, allowtng the uncaptured weight to undergo one

of the remaining processes by use of probabilities such as

g
and ag
above.
If capture is not treated by weights, but is regarded as an event

terminating the trajectory, the procedure of this section is modified in

the obvious way, and we do not include it.
In any event, at the present stage of the flow diagram, we should
have effected the decision of the type of nucleus hit, and the type of
collision undergone by the incident neutron, with assignment, if necessaxy,
of a nuclear parameter e, and adjustment of the neutron weight.

-88 -
The p o s s i b i l i t i e s include, of course, r e t u r n t o (a) i f capture o r weight

cutoff terminates the c a r e e r of the neutron i n question.

We proceed t o discuss r o u t i n e s f o r various individual processes

commonly encountered i n p r a c t i c e .

c
0
c
el. coll.

f
0 0.)
in. coll.

Fig. 33

3. E l a s t i c c o l l i s i o n s i n general. We d i s c u s s i n t h i s s e c t i o n t h e

c o l l i s i o n between two p a r t i c l e s , assuming conservation of momentum and

energy. We reserve c a p i t a l s f o r vectors, and make the following conven-

tions: R = (x,y,z) denotes t h e p o s i t i o n vector of a p o i n t i n t h e lab-

o r a t o r y system, V = = (A,$,%), t h e v e l o c i t y vector of t h e p o i n t R.

If Ri = (xi,yi,zi) are the p o s i t i o n vectors of a set of p o i n t masses

-89-
mi, i = 1, ..*,n, and m = ';mi is their total mass, then R = m-1smiRi

defines the position of the center of mass of these particles. We also

recall that if A = (al, a2, a3), B = (bl, b2, b ) are any two vectors,
3
2
their inner product is AB = Caibi, A2 = p i , the norm of A is
defined as IAI = fl and the cosine of the angle 8 between A and B

is given by cos 0 = AB/ /AI J B ~ .

Now consider two point masses 1 and %, with laboratory velocity

vectors V1 and V2, respectively. The total momentum of the system is

the vector P =ZtmiVi' its total kinetic energy is the scalar

1
k = Z- m 3
2 i i' and its total mass is m = Zmi. Since the equation
mV =Emi Vi holds for the velocity V of the center of mass, we have
always P = mV.
It is customary to define the velocities, total momenta, and kinetic
energies, relative to the center of mass, by the equations

v; = vi - v

1 2
k' = C- m. V'
2 1 i

We find from these definitions the following relations between

absolute and relative quantities:

-90-
or, since P = mV,

Note that the first of these equations states that, relative to

the center of mass, the particles always travel on the same straight line

with opp'ositely directed velocities.

We assume now a collision between these two particles, adopting

the notations indicated below for the relevant quantities before and

after collision:
Before After

i' i'

P Q

V W

P' Q'
k' R'
and restate the fundamental relations

Now the defining relations for elastic collision are

P = Q

k = a
Since 3p = Q and P = mV, Q = mW, we have V = W, which states that
the center of mass of the system proceeds unperturbed by the collision.
This equality, together with the equations ( 3 ) , and the fact that k = 2,
implies that k' = jl, that is to say, the relative kinetic energy is
unchanged.
Now P t = 0 in general, and we have therefore the equations

yvi = - %v;p Svi = % 2 %I2, and +i* - %%&2= 0, with identical

results for the relative velocities Wi, W; after collision.

Hence the equations

-92-
w e satisfied by V i2, Vi2 and also (since k' = 1')by .
17; 2 and 1;r. 2

Since the determinant of the linear system is - S(ml

1 "22 t- y2 n$) # 0,

the solution is unique and Wi

2
= Vi2 , Wi2 = V2' 2 '

We have, therefore,the following very simple picture (cf. Fig. 34) ,

in the center of mass system, of elastic collision.
A - -I_

The upper part of Fig. 34 indicates the relations we have derived:

The center of mass velocity is unchanged, the relative incoming (and out-
going) velocities are collinear and oppositely directed, and the relative

speeds are unchanged. "he hshed lines are coplanar, as are the solid
lines; however, the two planes they determine are not necessarily the
same.

-
The conditions for elastic scattering do not determine the angle
of scattering, as defined (say) by the angle si f'rom W to Wi (nor
the angle between the two planes referred to above). This must be given
in the form of a distribution law which depends in our case upon the

neutron energy and the nucleus hit. Such distributions w i l l be discussed

in the next section. O m immediate object is the derivation of formulas
for the angles JI, of scattering in the laboratory system, that is, the
angles between 'W and the Wi, and for the ldnetic energies
1
1, = 5 mi WF of the two particles in the laboratory system, as functions
of an arbitrary angle 11.;
We first consider the $
1:

-93-
W;
Fig. 34

-94-
w* = (w -t W ' ) 2 = (v + w q 2 = v2 + mi + wi2 = v2 f vi2 -I- 2(vIIwi\ COB 111'
1 1 1 1
= v2 + v i 2 + 2Jvllv;J cos rl'
1

WE = (w + w$ = (v 4- w$ = v2 + 2vw' + w;2 = v2 + Wh2 +' 2\vI Iw; I cos $6

2

= v2 + v i 2 - elvl Ivg cos J , 1

'

Next, we compute t h e cos #i:

= (Ivl -. Iv;l COS rli)/lW21

where t h e lWil * are given above, These formulas provide the general

s o l u t i o n t o our problem. I n most cases t h e energies of neutrons are

suf'ficiently g r e a t e r than the (thermal) v e l o c i t i e s of nuclei t o admit

t h e assumption that t h e latter a r e e f f e c t i v e l y s t a t i o n a r y i n t h e

laboratory system.

-95-
 We proceed t h e r e f o r e t o s p e c i a l i z e t o the case of neutrons of mass
y,laboratory velocity V
1
s c a t t e r i n g on n u c l e i of mss m2, v e l o c i t y

V2 = 0. I n t h i s case we have

mV = yVl f m2V2 = y V l , Iv I= m-lml IV,~

Vi = V1 - V = V1 - m-1y V l = m
-1
m2Vl

Vi = V 2 - V = -V = -m -1m 1V 1

S u b s t i t u t i o n gives us

:w = v{:m -2 2
ml + -2 2 -2
m m2 + 2m y m 2 cos q1
f 1
1 2 1 2
and hence, i f = 2 y W 1 , kl = 2 y V l are t h e ( l a b o r a t o r y ) energies
of y after and before c o l l i s i o n , we obtain

l1/kl = m-2{< + Tm2 + m2

-T m2 + m cos
12
+.i

It is convenient t o introduce t h e q u a n t i t i e s A = m /m ,
2 1

-2 1 - -2
s I (1
= -2 + 'F) = 1 - 2m ym2, and T = 2 (1 - r) = 2m m 1m2' obtaining

finally

= kl ( S + T COS JIi)

-96-
 Moreover,

Although we make no use of them i n t h e sequel, we include for

completeness t h e r e s u l t s f o r t h e s c a t t e r e d nucleus
52

and COS *2 = f- - sing) , which the reader may v e r i f y

as an exercise.

The s a l i e n t r e s u l t s of t h i s s e c t i o n are t h e formulas

E ‘ = E (Se + T,P)

where w e agree on the following d e f i n i t i o n s :

me = mas of nucleus of type e

m = mss o f neutron

Ae = m /m
e
- +
re = (Ae - 1 ) 2 / ( A e 1)2
se = 2
1 (1 + Fe)

1 -
Te = (1 - re)
E = neutron energy before scattering
laboratory system
E ’ = neutron energy after scattering

-97-
 c c = cos Iq
II; = angle of deflection in center of mass system from original
line of flight (coincident with V = W since V2 = 0)

a =

q1 = angle of deflection in laboratory system from original line

of flight

The constants Ae, Se, Te are nuclear (energy independent) constants,

the subscript e, if required, being set in the routine of the preceding

section.
It may be noted that in the case of scattering on hydrogen, we
1
may assume A = 1, and the formulas become simply 5 = 0, S = T = 5'
E' = E (1 + P), and a = f-. Several remarks are of interest
2
here: (a) we have the relation E' = Ea ; (b) scattering on hydrogen is

v m
always forward in the laboratory system, viz., a = cos II12 0; (c)

a = cos = = cos ($i/2) implies G1 = *i/2; (d) the

formula a = f- should be used for hydrogen, since the general

formula is indeterminate at c = -1; (e) if scattering on H is isotropic

in the center of mass system (a good assumption, incidentally), we have

fl = 2 r - 1, and a = fl,which shows that the laboratory scattering i s
in the cosine distribution. (Cf. Chapter 11, 85b.)
F i n a l l y we observe that f o r "heavy1'elements (Ae large) we have
-
the following approximate results: re = 1, Se = 1, Te = 0, E' = E, a = cc.
If these relations a r e of acceptable accuracy, one may ignore the energy
change and determine a directly instead of N by the methods of the

-98-
 next section. I n t h i s connection one should r e f e r t o Chapter V I I , gpj,

f o r t h e case of s c a t t e r i n g i s o t r o p i c i n t h e l a b o r a t o r y system.

4. The d i f f e r e n t i a l e l a s t i c s c a t t e r i n g c r o s s s e c t i o n . We have

seen i n t h e last s e c t i o n that t h e preservation of t o t a l momentum and

k i n e t i c energy c h a r a c t e r i z i n g t h e e l a s t i c s c a t t e r i n g process does not

serve t o determine t h e angle of s c a t t e r i n g i n t h e c e n t e r of mass

1
' system, but does detemnine the new energy E ' of the s c a t t e r e d neutTon

and i t s d e f l e c t i o n cosine a i n t h e laboratory system i n terms of a

given LI = cos I)'

1 *
The determination of P is gcxverned by a physical d i s t r i b u t i o n
e
function u E ( P ) , whose u n i t s a r e barns per s t e r a d i a n , and which depends

upon t h e s c a t t e r i n g nucleus e and t h e energy E o f t h e i n c i d e n t neutron.

Specifically, 0; ( Q ) dS1 i s defined as the cross s e c t i o n ( i n t h e sense

of Chapter III,€j1)presented t o i n c i d e n t neutrons of energy E by a

nucleus of type e f o r the process of s c a t t e r i n g i n t h e c e n t e r o f a s s

system a t an angle i n t h e dSa neighborhood of t h e d i r e c t i o n Ja with t h e

i n c i d e n t l i n e of f l i g h t . Thus by d e f i n i t i o n ,

the i n t e g r a t i o n being over t h e e n t i r e surface of t h e u n i t sphere i n

e
d i r e c t i o n space U, V, W , and (J (el.) being t h e ( t o t a l ) e l a s t i c s c a t t e r -
E
ing cross s e c t i o n f o r t h e element e a t energy E. I n terms o f s p h e r i c a l

coordinates (azimuthal angle y , polar angle I):, 35),

Fig. we have

-99-
 W
tI
incident line of flight

C. M. system

Fig. 35

-100-
 In a11 problems of the text u e is independent of
E 9 and we
have

or, in the form usually given,

0 e (P) dlu e (el.)

= aE

It f o l l o w s that Pz (B) = 2n 0; (P)/'z (el.) is the probability

density function for elafatic scattering at the direction @.i,

f'rom the

line of flight in the center of mss system.

Hence, the Monte C a r l o procedure sets
nlu
r =J PE dp
-1

for the determination of from the random number r.

-101-
 5. A routine f o r e l a s t i c scattering. Suppose t h a t e l a s t i c

s c a t t e r i n g on t h e ( l i g h t ) element e i s i s o t r o p i c i n the c e n t e r of

mass system. This means t h a t t h e function oe(cl) i s a constant, and

t h e f i n a l formula of t h e preceding s e c t i o n gives simply cl= 2r - 1.

I n such a case, w e should e n t e r a r o u t i n e l i k e t h a t i n Fig. 36. Here

E r e p r e s e n t s t h e "energy c u t o f f , " namely, t h e l a s t of t h e lower bounds

G
E l > E2 > ... >E
G
of the energy groups adopted. Energies f a l l i n g below

t h i s n e c e s s i t a t e c l a s s i f i c a t i o n of the corresponding weight i n a t e r m i n a l

category reserved f o r l o s s t o energy c u t o f f . The loop on E

e: -E
begins by comparing t h e new energy E with t h e lower bound E of t h e
63
group which t h e neutron occupied before s c a t t e r i n g , s i n c e e l a s t i c

s c a t t e r i n g cannot raise t h e energy. Note a l s o t h a t EG - E < 0 at entry

t o t h i s loop ensures i t s termination a t some g 2 G. If EG > 0, t h e

previous E - E decision is mandatory.

G
The formulas f o r E' and a
are f u l l y discussed i n t h e preceding 83. A s remarked there, t h e pro-

cedure f o r t h e s p e c i a l case of hydrogen i s simpler and i s indicated i n

the lower h a l f of F i g , 36. For t h e case of s c a t t e r i n g i s o t r o p i c i n the

l a b o r a t o r y system, we again refer t o Chapter V I I , 85. This a p p l i e s

approximately t o the case of heavy elements with an i s o t r o p i c l a w i n

t h e c e n t e r of mass system.

The ( a ) of Fig. 36 refers t o a r o u t i n e f o r f i x i n g t h e f i n a l .

l a b o r a t o r y d i r e c t i o n parameters of t h e s c a t t e r e d neutron ( c f . Chapter

V I I ) before returning t o ( p ) .

-102-
Fig. 36

-103-
 Most elements have more complicated differential cross section

functions, usually characterized by preferential forward (and sometimes

backward) scattering, an effect which tends to increase with higher

incident energies. The problem of choosing CC from such a distribution

is essentially that (already discussed in Chapter 11, $5) of selecting

the direction parameter w from a given source distribution.

If the problem can be restricted to a single energy, or to

relatively few energy groups(lo)over each of which the differential
function bE (P) may be considered constant, one may well store tables
f o r the probabilities P of scattering at cosines p 2 Lc , where
g,J 3
fl
1> y-.., >gJ = -1 are suitably chosen bounds of subintervals of the

cosine range, and the

are pre-computed by numerical integration and stored. The routine for

determining u is then that of Fig. 13, where we read j for i, P
Q, j
for Pi, CC for wi, and CC for w.
3
Aside f'rom its demands on storage space, this is subject to the
errors of interpolation, which may be difficult to minimize in cases of
-. -.C
,I _- .-

(")These need not coincide withmr be as numerous as the energy groups

reserved for free path and scattering-type probabilities. It is not
unusual to carry two or even three different sets of energy classifications
with corresponding indices g , h,... as neutron parameters.

-104-
very forward scattering.
e
If the original curves (9
E (U) can be simply fitted, the simplest
scheme may be the von Neumann device (Cbpter 1 I 8 5) We should then

have a formula, for the function

and use the routine of Fig. 5 with a = -1,b

*
1, p =
*
= 0 (e, E, NU),
reading cc for e
Occasionally, the function 0
*(e, E, cc ) m y be fitted easily for
each energy E as a function of P , say,
*(e,
(I E, P ) =
e
AE -k
e
BEP + CE
e
CC
2 - 1 I P 6 1

but the coefficients m y prove intractable as functions of E, In such

a case we may be able to store the coefficients

for each of a reasonable number of energy groups, and use the routine
of Fig, 5 as indicated above, the machine computing o
*(e, g, P ) =

Ae + BeP + Ce CC from P and the stored coefficients.

@ ; & g
Finally, and in the worst cases, it may be necessary to store

u(vm, En) is the differential cross section for a particular element

(we drop the index e temporarily) evaluated at a suitable set of
...>EN, where
I

values
c1

p0 = 1> .. N

>P
M
= -1 and Eo.
cv N

Eo is at least

-105-
the highest energy encountered in the problem, and = EG the least
such energy. The energies and cosines used f o r this purpose may be
hand-picked, and need not coincide with those used f o r other purposes
elsewhere in the problem. We store a table of the form

t.m 0 1 ... N

0 (r*
00
1

.
0

U*
M M,N

N
En- 11 z0 Gl N

EN

Fig, 36a

and resort to the double interpolation routine of Fig. 37.

However we may determine C , the routine should lead to a

determination of E and g after collision and the deflection cosine

tt f o r the laboratory system as indicated in Fig. 36.

-106-
 I i

I gl = o'*
m-1, n
+ (.* m-1, n-1 - u*m-1, n

Fig. 37

-107-
 6. Differential elastic cross section for the laboratory system.
If a differential cross section
L(a) is given for elastic scattering
(J

relative to the laboratory cosine a = cos IL one may choose between

1'
the following two alternatives.

(a) One may employ (J,(a) directly to determine> a, just as

the center of mass differential cross section o(P) was used to determine
c1 . It is then necessary to use the formula

-1
P = A

to determine the corresponding center of m8s cosine P, which is needed

in computing the new energy E' = E(S + Tp).

The above formula is obtained by solving the equation

(5)

for CC in terms of a. The ( s g n a) choice of sign on t h e radical is

dictated by the following considerations: From equation ( 5 ) it is clear
that the sign of a is always that of 1 + A M . In solving equation (5)
for P we obtain

1+ = a2 f v [ ( l- a2)2 - (1 - a2) + .*A2]

so that the sign of the right side must always be that of a. It is

clear graphically that f o r the function f(a) = a2 f .rr;t?;lT to have the

sign of a, it is necessary and sufficient that one use the upper branch
for a > 0 and the lower for a C. 0.

-108-
 Formula (4) m y be r e w r i t t e n i n the form

where b2 = 1 - a2 .
( b ) A l t e r n a t i v e l y , we may use the defining r e l a t i o n

a ( P ) dP = uL (a) da

t o compute the c e n t e r of mass function a ( P ) i n advance, and then use

the usual procedures of $5 directly, This m y be done by means of the

formulas

where

and

7. A -
- weight device f o r e l---"
a s t i c scattering. I n problems involving

many c o l l i s i o n s i n a medium admitting e l a s t i c s c a t t e r i n g on l i g h t elements,

i t my prove worth while t o avoid l o s s of t r a j e c t o r i e s , which are u s u a l l y

followed through m n y c o l l i s i o n s , t o t h e energy cutoff Eo. These losses
m y be obviated by the following device, which makes further s t r a t e g i c

u8e of weights.

We have seen that t h e r e l a t i o n between emergent and i n c i d e n t energies

-109-
i n such c o l l i s i o n s i s E ' = E(S i- T P ) , where c1 i s t h e c e n t e r of mss

d e f l e c t i o n cosine, chosen randomly by t h e equation

Now,for a given incoming energy E E t h e --

least emergent energy
G'
is ETdn = E(S - T ) . If t h i s exceeds EG, the energy c u t o f f f o r the

problem, one proceeds as usual, with no p o s s i b i l i t y of loss t o energy.

However, i f E(S - T ) = E
<
that i s t o say, i f E d EG/(S - T) E E
*,
G'
t h e r e exists a c r i t i c a l cosine P C ' depending on E, such t h a t

E(S + Tuc) = EG, with -14 P

C
C 1, and

r e p r e s e n t s t h e p r o b a b i l i t y of s c a t t e r i n g with ,u on -15 ,u 5 c1 with

C

the r e s u l t i n g energy E ' 4 EG . Moreover, assuming P > pc

i s t h e proper Monte Carlo formula for determination of t h e cosine

on t h e range Pc 5 r(l 5 1, with a r e s u l t i n g energy E ' => EG. We may

t h e r e f o r e proceed as i n Fig. 38.

For t h e simple case of s c a t t e r i n g i s o t r o p i c i n t h e center o f

mass system, w e have t h e formulas

1
f = - (1 .t p C ) , P= %-1(r) = l-r(l-pc),
-1
P = pE (r) = 1 - 2 r .

F i n a l l y , w e note t h a t a neutron being followed by t h i s method

-110-
automatically has weight W > 0 and energy E ,Ec, if the random
number r = 1 does not a c t u a l l y occur i n our sequence.

(1 - f ) W -w
I
I I

E' = 5(S + Tp) - E

I
&

I a (Cf. 95)
I Fig, 38

-111-
 8. --
Fission. Suppose that a neutron causes f i s s i o n of some

(heavy) nucleus, r e s u l t i n g i n t h e production of a n average number u

0

of neutrons with energies of r e l e a s e independently d i s t r i b u t e d i n a

given f i s s i o n spectrum, and whose emergent d i r e c t i o n d i s t r i b u t i o n s a r e

i s o t r o p i c i n t h e laboratory system.

It i s possible t o invent a p r o b a b i l i t y d i s t r i b u t i o n p ( n ) with

V
0
= x n p ( n ) , decide on t h e number n o f progeny i n a given f i s s i o n

by chance, and f o l l o w t h e s e n neutrons individually. This may be

done even i n s u p e r c r i t i c a l systems i f recourse i s had t o time cycles

and census taking. However, t h i s i s unnecessary and undesirable. We

may i n s t e a d follow a single neutron, of weight vo times that of t h e

i n c i d e n t neutron, chosen from t h e f i s s i o n spectrum of energies, and

directed isotropically. This may reduce f l u c t u a t i o n s , enormously

s i m p l i f i e s t h e code, and decreases machine running t i m e .

Suppose that f ( E ) dE is t h e p r o b a b i l i t y of f i s s i o n energy E

between E and E + dE. We choose a set of energy i n t e r v a l s with

bounds Fo> Zl> . . ., %, where z0 i s t h e highest significant fission

energy and 2H t h e energy c u t o f f f o r t h e problem. These

P
Eh need

not agree with bounds of energy groups used for other purposes. Define

We m y then refer t o Fig. 39. Note t h a t i f t h e Xh are d i f f e r -

e n t from t h e used f o r ( s a y ) s t o r i n g free paths A
E one must d e t e r -
63 €3’
mine t h e index g, beginning t h e loop with g = 1 s i n c e f i s s i o n can

-112-
-.raise
-. energy above the incident energy E in group g . The exit (a)
refers to the routine of Chapter VII,$5,for scattering isotropic in the

laboratory system.

In some problems, designed to find the distribution, say, in

a number of cylindrical shell zones defined by radii R1 < R2< ...< Rm
and heights H1< H g < ...<Hn, of fissions resulting from a given spatial
distribution of fissions, fission is a terminal event, and counters
N are reserved for the number of neutrons terminating in fission in
i,J
radial zone i,height zone 3 . In the event of fission we should then

f o l l o w the simple classification routine of Fig. 40, rather than that

of Fig. 39.

- Eh
N N

E=: +(r-F) Eh-l

. h Fh-l - Fh

Fig. 39
 -i + 1- i 4 - j + 1 4 4-

-
7'

@ , 0
- R2 - - - l - - j -
R =x +y 1-i R2
io Z-H
_t

Fig. 40

9. I n e l a s t i c ( n - n ) c o l l i s i o n s i n general. Using t h e n o t a t i o n

and procedure of $ 3 of t h i s chapter, we now discuss a type of c o l l i s i o n

i n which a n i n c i d e n t neutron of k i n e t i c energy kl = ?VI2 = E strikes

a nucleus of mass m2, imparting t o t h e l a t t e r a known energy of e x c i t a -

t i o n e , and emerging with an energy l1 = E'. We obfain again

t h e fundamental r e l a t i o n s

(6)
P = mV & = M (7)
k = -' m V 2 + k ' t = - m1 W 2 +I' (8)
2 2

The momentum and energy c o n s e m t i o n l a w s now read

 P = Q

k=,!+@
however, and we f i n d V = W as before, but now
k'= 4' -+ E

We f i n d again that t h e equations

determine Wi2 and Wi2 uniquely, t h e d i f f e r e n c e being only that 4'

is not simply k' b u t k' - E. Thus t h e speeds of departure of ml
and m2 i n the c e n t e r of mass system are uniquely determined, although

they are not simply the incident relative speeds as before. The angle

$i must be determined according t o some physically determined d i s t r i -

bution l a w , and from t h i s the laboratory energy and d e f l e c t i o n cosine

may be computed.

We w i l l derive the r e l e v a n t formulas, F i r s t , solving the above

2
system f o r Wi , W; 2 , suppose tbt we first rewrite t h e system f o r

i n the form
where e
1
= 2m2e/mlm, 1 2m, and
E2 = 2me/m
1
%E1 + $ m2e2 = e ,
2
-$el + m2 c2 = 0. We see a t once t h a t

We d e a l only with t h e neutron %. Computing W2 as before

i n terms of cos +i, w e obtain

w2 = (w + wi)2 = (v + w;)2 = v2 + w;2 + 2IVlIWiI cos $i

1

and for cos JI,, the laboratory d e f l e c t i o n cosine,

cos *1 = vw,/[vl lWll = v(w + W i ) / l V l Iw,( = v(v + wi),

-1 -1
Again s p e c i a l i z i n g t o t h e case V2 = 0, V = m ?VI, V i = m m2V1,
we obtain

-116-
and theref o r e

( E , / v2 ~ )I- q m 2 m -2 cos #i-y

But B
1
= s o that rn2el/rn$f = mc/m2kl = m c / m 2 , and
hence, f i n a l l y , l! by E' and kl by E

Moreover,

or
m-'y + m-lm
2
cos $i
cos @l =
(.-'."
_I_

1 + m-2rng
.----

-
2
(me/%E) + 2m-2nipi,os#.$~)}'2

Now the r e a c t i o n can occur only i f !

I
, = k' - 0 => 0 ; since k' =
2; miVi2 = -2}V; = %m-h,t h e condition is simply

-117-
E 1 rn -1
m2 c . If we define t h i s c r i t i c a l energy as e c , we may write

1+ A COS $i\/L- (EC/E)

___ .-- _------
cos JI = --- ___I__,________ ....
1
~

(1 + A2 [l - (tc/E)] + 24 COS JIif-)} 2'1

About t h e s e formulas w e make t h e following remarks: ( a ) For

6 = 0, they reduce, as they should, to those f o r e l a s t i c s c a t t e r i n g ,

C

( b ) for e c 7 0 and A l a r g e (heavy n u c l e i ) , they become approximately

E' = E -e and cos I

)
1
= cos @i,as we should expect.

It must a l s o be pointed out t h a t a nucleus m y have various

excited states each with i t s own e x c i t a t i o n e n e r g y e , and a correspond-

ing i n e l a s t i c cross section 0

E
(c ) f o r neutrons of energy E 2 (m/m2) e =ec.

For l i g h t nuclei these must be d e a l t with individually, using t h e above

formulas with t h e appropriate e

C
.
For heavy n u c l e i , t h e states may be w e l l separated i n some cases,

and m y be d e a l t with individually,using t h e simple r e l a t i o n s of rerrark

( b ) above. If, however, t h e excited states are very c l o s e together,

one r e s o r t s t o such methods as those of t h e following two s e c t i o n s ,

10. I n e l a s t i c (n-n) c o l l i s i o n s on heavy

-..
nuclei. If t h e states
of t h e nucleus are c l o s e l y packed on the energy s c a l e , it i s c u s t o m r y

t o give a t a b l e of experimentally determined p r o b a b i l i t i e s for

'h,h'

-118-
 ... >-EH is a suitable
an incident neutron of energy group h to produce in (n-n) reaction a
2 Eht,where h' 3 h,
c c1

neutron of energy E' and Eo>

set of lower bounds for energy. One may then proceea to determine the
energy E ' in the usual way by uqe of a random number and linear inter-

polation, as in Fig. 41, which is drawn for the case of isotropic

scattering in the laboratory system, and therefore exits to (5) (Chap-

ter VII, 85). Note that the Ph,h, form a triangular matrix, with Ph,h'
non-zero only for h' 2 h, and Ph,H E: 1, h = 1, ..., H.
Some inaccuracy is unavoidable in such a method, since, strictly

speaking, 1 is a function not only of the group h, but of the energy

E in this group.

11. Inelastic (n-n) collision with Maxwell distribution. Con-

sider an (n-n) collision of a neutron of energy E < @ with a heavy
nucleus for which it may be assumed that the emergent energy E '

probability density function is proportional to

where T = @fm,
being a given constant.
a
 h' + 1 - ht

E'
hf
-- E
h

E - E
g 0-g + l - g
Fig. 41

P E ( E ' ) dE' = ( l / A ( E ) ) E ' exp(-E'/T) a'

and t h e Monte Carlo p r i n c i p l e involves determination of E' from E

and r by means of

-120-
 The s o l u t i o n f o r E' i s d i f f i c u l t and we use i n s t e a d t h e

device of Chapter 4 5 5 . The latter r e q u i r e s t h e function

p*(E') = n(E')/max N ( E ' )

t h e maximum being f o r t h e range 0 5- E' 6 E. Now t h e maximum of n ( E ' )

i s a t E ' = T, and we m y use t h i s provided T = a im C E, which

mans that the incident energy E shall exceed 01

2/E.
- L e t us suppose
2
f o r s i m p l i c i t y that c1 /Fc: EG, t h e energy cutoff (as i s usually the

case). We have then

*
p (E') =
e
(T)E'e
-E'/" <
-- 1

f o r a l l incident energies E on the range EG 4 E 2 E and we may pro-

ceed i n t h e usual way, with E ' = rE, e t c .

- f o r f i s s i o n a-.b- l e nuclei.
12. A- ----combined t r a n s f e r matrix I__---
If t h e

nature o f a problem involving a f i s s i o n a b l e nucleus does n o t demand

I
keeping t h e i n d i v i d u a l types of c o l l i s i o n separate ( f o r example, f o r

purposes ,of t a b u l a t i n g capture, f i s s i o n , etc.) and if t h e neutrons

emerging from each type of c o l l i s i o n are a11 emitted i a o t r o p i c a l l y i n

t h e laboratory system, a much simpler method can be used t o g r e a t

advantage, i n place of the several i n d i v i d u a l procedures i n d i c a t e d i n

previous sections. L e t us adopt t h e following notations, i n atidition

-121-
 p(h,h')= probability of scattering from group h -
to -group
- ~ - h'

(inelastic collision)
V = average number of emergent neutrons per fission
0

f(h') = probability of fission neutron energy ---

in group h'
= Kronecker b function, with value 1 f o r h=h', 0 otherwise
dh,h'
Now, assuming a collision of a neutron of group h w i t h the
-
nucleus, it is clear that the expected number v of neutrons emerging
h,h'
from the collision in group h' is

-
V h,h' ={oh(cap.) o+ oh(el.)dh,hl + oh(in.)p(h,h') +

ah( fiss ) f( h ' ) v O } / U h ( tot. )

Thus, the total. expected number of neutrons per collision is

Hence, the probability of an emergent neutron from such a colli-

sion being in group h' is

h'
We may therefore store a table of
%,h' = c q h , k
k=l
h = l,..., H, h' = l,.. ., H, and
Q
h,o
= 0, h = 1,. .., H, as shown
below,and proceed according to Fig. 42.

-122 -
 0 1 2 . . . h ' . . . H

1 0 1

2 0 1

h 0 1
'h, hf

Io 1

Fig. 42

-123-
 13. C o l l i s i o n.-s s h a t t e r i n g a nucleus. We consider i n t h i s

section a ty-pe o f c o l l i s i o n i n which a neutron of rest-mass m and

1
velocity V e n t e r s a nucleus of rest-mass m and v e l o c i t y V = 0,
1 2 2
nucleus system i n t o a s e t
r e s u l t i n g i n a s h a t t e r i n g of t h e neutron i-

o f fragments o f msses n
3
and v e l o c i t i e s W
s' j = 1,2,. .. , the

t o t a l r e s t mss n = En o f t h e fragments being greater than m = ?L+"2'

3
We treat t h e ~mechanics
- - . of t h e c o l l i s i o n n o n - r e l a t i v i s t i c a l l y , using t h e

formula -21 mV 2 f o r k i n e t i c energy of a p a r t i c l e of rest mss m, and

v e l o c i t y V, and w e suppose that m = n a t various points of t h e argu-

ment. The a c t u a l m s s - d i f f e r e n c e n - m corresponds t o an energy

e = (n - m)c2, where c i s t h e v e l o c i t y of light, and i n order f o r t h e

process t o occur, part of t h e k i n e t i c energy of t h e neutron must be

used t o supply t h i s energy e .

Proceeding as i n sf3 and 9 of t h e present chapter, we have t h e

general equations

> P = mV Q = W

k = k' + -1m V 2 + - r1r l w 2

2 2

(note t h e assumption m = n ) , and t h e conservation l a w s now read

P = Q

k = ,(+e

As usual, V = W, and so k ' =I'+ . J u s t as in Sg,

e we f i n d

that i n order f o r the r e a c t i o n t o take place, we must have k'= m

2
m-h 2E

-124-
or

Assuming such a n i n c i d e n t energy E, t h e conditions on t h e

relative velocities W' o f t h e fragments are

J

Cf=n.w(. = 0
J J

I n c o n t r a s t t o the s i t u a t i o n i n $ 9 , these r e l a t i o n s do not deter-

mine the W!2i f more than two fragments e x i s t (IncLuding the original
J
neutron). We propose i n s t e a d to s i n g l e out an a r b i t r a r y one of t h e

J and discover the maximum r e l a t i v e energy w i t h which it

fragments n

can emerge.

We adopt the convention that E' refers to summation on a l l

j # J. Thus we define n ' = x'ns t o be the mass of the remaining

fragments, and nfWR = c ' 11J.W j defines the v e l o c i t y W

R
of t h e i r

own c e n t e r of mss, Now WI; = WE - W denotes as usual t h e v e l o c i t y of

the residual c e n t e r of mass R r e l a t i v e t o that of t h e whole system.

Moreover, we define W" = w4 .- WB = ( ~+ jw)

j
- (vi c W) = W; ,- Wr; f o r
a l l j # J. Now the above system becomes
 But observe t h a t t h e f o l l o v i n g r e l a t i o n s hold f o r t h e r e s i d u a l

system:

z: 'n.Wt = ztnj( w i + W') = ntw; + zCln

jW"J (10)
J J J

Now we know by d e f i n i t i o n that

so that

and hence
n .W ! = ntWI;
J J

Reference t o (10 ) y i e l d s

and t h i s result with (11).gives us t h e familiar

Now t o return t o t h e system ( g ) , its second r e l a t i o n together

with ( 1 2 ) t e l l s us t h a t
nlW; = -nJW'J

-126-
 1 n'Wi2 = (nJ/n')(--1 n W ' ?-).
s o that - e m y therefore r e m i t e ( 1 3 ) a s
W
2 2 J J

The l a t t e r , s u b s t i t u t e d i n t o the f i r s t member of t h e system (g),

gives us f i n a l l y

I t follows that t h e mximum -..--..-

r e l a t i v e energy of nJ must be

Thus, i n sunwry, if t h e i n c i d e n t energy E exceeds t h e minim1

energy m m - I e: required f o r t h e process, the maximal relative speed

2
I of any fragment J i s given by

w -
J m n,
J

Suppose that we now construct t h e laboratory v e l o c i t y WJ and

the corresponding laboratory d e f l e c t i o n angle 3, by t h e usual a d d i t i o n

J
W
J
= W + W'
J
= V + W; . Reference t o Fig. 43 mkes it c l e a r that the
__
possible end p o i n t s of the vectors WJ occupy t h e sphere of r a d i u s wJ
about t h e end p o i n t of V. Moreover, it i s apparent that three e s s e n t i a l l y
d i f f e r e n t cases arise, depending on the r e l a t i v e m g n i t u d e s of qJ and
-
IVl, o r , equivalently (using t h e above formula for w ), the relative

magnitudes of E
--
and e = nte/(m2 - nJ).
J
If we are i n t e r e s t e d only i n

-127-
 Case 1. wJ( ] V I o r E ( 7. The angle @J i s l i m i t e d t o t h e
range 0 <=
J 5- arc s i n (wJ/ V I ), and t h e possible speeds IWJ I c o r r e s -
rl

ponding t o any given @J on t h i s range are bounded from zero.

-. ‘c

Case 2. w = IVl or E =e . The angle i s limited t o the

J
range 05
- @J <= 7r/2 b u t f o r a given such $ t h e possible lWJl range

f r o m 0 t o a maximum.

Case 3.
.-
wJ 7 I VI o r E 7 B
c
. The angle @ J has i t s f u l l range

0 <- rl,J <= T,, and IWJ I for a given rlJ ranges from 0 t o i t s upper bound.
I n any case, t h e maximum speed w a s s o c i a t e d with a given
J
possible
qJ i s given by t h e g r e a t e r r o o t of t h e quadratic equation

Fig. 43

-
--2
WJ = WJ
2
+ v2 - 2 WJI I cos qJ. Substituting f o r wJ and writing
J
f o r the maxim1 laboratory energy nYJ o f t h e Jth fragment a s s o c i a t e d
with t h e l a b o r a t c r y d e f l e c t i o n angle #J, we obtain

If we apply t h i s t o any neutron nJ, we have n = y and

nf = m2, so t h a t f i n a l l y

The discussion of t h i s s e c t i o n thus serves t o show why t h e i n -

elastic cross s e c t i o n U,(in. ) for such a process becomes zero f o r

u

i n c i d e n t energies E below m m-le

2
, and how t h e energy range of

neutrons emerging from t h e r e a c t i o n depends on the i n c i d e n t energy

and t h e cosine of t h e laboratory angle @J.

14. The (n-2n) r e a c t i o n i n deuterium. As a n example of t h e

preceding theory, w e consider the r e a c t i o n

n + ,H 2 --+ lH1 + 2n

i n which a neutron d i s r u p t s t h e deuterium nucleus i n t o i t s c o n s t i t u e n t

proton and neutron.

Adopting the atomic masses of t h e following table

neutron 1.00893

1 1.008123

LH2 2.014708

we f i n d a n excess of .002345. This gives a n actual mass of

.002345/A gm p e r deuterium atom. Multiplying by c2 and dividing

by t h e number 1.60203 x 10
-6 of e r g per Mev (Chapter 11, 53) gives

= 2.184 MeV. Hence m miL€ = 1.5 6 i s t h e minimal energy f o r t h e

process t o occur, and t h e three cases mentioned i n t h e preceding s e c t i o n

depend on the relative sizes of E and m26 /(% - ml) = 2 e ,

We w i l l suppose t h a t a "cross section" T ~ ( E ' ,J I ) i s defined i n

such a way that

i s t h e expected ---number of neutrons of energies between E' and E' + dE',

and d i r e c t i o n s between @ and $ + d $ r e s u l t i n g from such i n e l a s t i c

c o l l i s i o n s of a beam of B neutrons of energy E t r a v e r s i n g a medium

of numerical d e n s i t y N and thickness A t . (Contrast Chapter 111, $Lo)

Now

i s t h e --
t o t a l expected number of neutrons i n t h i s t r a v e r s a l , s i n c e each

such c o l l i s i o n produces two neutrons.

Hence w e have

where the (E, $ ) are the upper and lower energy bounds referred t o
E
2
i n the preceding s e c t i o n .
 Moreover,

-'(i.n.)T
nu E E ( E ' , $') dE' d ( c o s I))

i s t h e p r o b a b i l i t y of an emergent neutron being i n t h e range (E, E + dFI)

and (I),# + a # ) , assuming such a c o l l i s i o n occurs.

The Monte Carlo procedure is t h e r e f o r e clear. W

e decide i n the
usual way whether the c o l l i s i o n i s of t h i s kind by rererence of a random
nwnber t o QE(in.)/uE(tot.). If r i s less than t h i s r a t i o , i n e l a s t i c
c o l l i s i o n occurs. We then double the weight W of the i n c i d e n t neutron

and determine i t s c o s $ = a and i t s energy E' from the above p r o b a b i l i t y

,distribution. We may do this i n two steps, using first t h e p r o b a b i l i t y

d e n s i t y function

t o determine the a = cos rl, and then finding the energy E ' using t h e

d e n s i t y function f o r E ' :

f o r the a = c o s $ determined.
A more complicated example involving a n energy cutoff and

c r i t i c a l angles i n t h e same way as indicated i n 67 of t h i s chapter may

be found i n LA-1606 (not a v a i l a b l e ) .

-131-
 15. An (n-2n) reaction on heavy nuclei. A s a further illustra-
tion of the use of weights, consider a collision of a 14 Mev neutron with
a heavy nucleus which results in the emission of a pair of neutrons,
isotropically distributed in the laboratory system, one neutron being in
the energy distribution

E' exp(-E'/T1) dE' , 06E'(=e , T1<

and the other in a similar distribution

E' exp(-E'/T2) dE' , 0 5 E' 5 e , T2< e

where e, T1, T2 are constants.

for i = 1,2. We have then

-pi(E') d E 1 = AI1 E' exp(-E'/Ti) dE'

for the probability of the ith neutron emerging between E' and
E' + dE1. Since the expected number of neutrons in the latter range,

per collision, is pl(E1)dE' + p2(E')dF,', we may properly a8sume one

neutron of weight W suffering such a collision gives rise to a neutron

of weight 2W, chosen in the energy distribution

-132-
 Since t h e maximum mx p ( E ' ) on the range 0 <, E' (= e can be

computed i n advance, we may use the von Newnrznn device of Chapter I, s!jr

followed by a n induction loop for, determining t h e new energy group, and

the procedure of Chapter VI& @ , f o r the final d i r e c t i o n parameters,

16. Capture i n a smll zone. Certain problems involve determina-

t i o n of absorption i n a smll zone of material, surrounded by a r e l a t i v e l y

l a r g e system of moderator, which slows down neutrons by s c a t t e r i n g on l i g h t

nuclei. The s l i m chance of a neutron h i t t i n g t h i s amall capture zone

may be improved by various devices, of. which we i n d i c a t e only one.

Suppose f o r s i m p l i c i t y that a s p h e r i c a l l y symmetric system con-

t a i n s a small c e n t r a l core of capture mterial of r a d i u s R1. L e t R* 7 R1

be chosen comparatively c l o s e t o R.
I' If a neutron undergoes c o l l i s i o n
a t a point of r a d i u s R R
* we my proceed as usual. However, i f

R
*
5 R , we may process from this point t o termination 5 neutrons each

of weight W/z i n s t e a d of the customary s i n g l e neutron of w e i g h % W.

-x.
If the m u l t i p l i c i t y m i s s u f f i c i e n t l y large and R - R1 sufficiently

small, some of the h descendants are very l i k e l y t o h i t t h e absorbing

core.

Such a scheme r e q u i r e s modification of the o v e r - a l l flow diagram.

We c a l l a t t e n t i o n e s p e c i a l l y t o the following chawes: ( a ) the source

reutine ( 0 ) sets a new parameter m = 0 f o r each fresh neutron l a v i n g

t h e source j ( b ) t h e ( B ) o r ( P o ) r o u t i n e exits, i n event of c o l l i s i o n , t o

a new e n t r y ( 7 ) which i s preliminary t o the usual capture r o u t i n e (9');

-133-
(c) the ( 7 ) routine, in case the multiplication trick is indicated,
- .
- - - - - ...
stores all p m e t e r s R, w, E, g, W,.. of the "parent" neutron at

the point of collision in new positions designated by €3, w, E, g, W,

f o r reference in processing each of the progeny from this point to
termination; (d) the ( a ) entry, to which one returns on termination of

any particle, is modified to order complete processing of progeny

before starting out a new source neutron. The parameter m = l,2,...

-m
indicates the number of the descendant being processed. The following
flow diagram (Fig. 44) includes the essential modifications.
Note that the entry ( Y ) excludes using the multiplication trick
on a descendant. That is, we do not iterate the process. In practice
the device is usually of a more elaborate nature, using different
multiplicities m for different critical radii R
*. Moreover, the
capture cross section of the core usually becomes significant only at
low energies, so that the multiplication decision may also rest upon
a decision on E - E3c JC
where E is some stipulated low energy.
71

StoreNN N N N
R,w,E,g,W, ... at R,w,E,g,W ,...

Fig. 44
 17. Capture by a "point" d e t e c t o r . The method of t h e preceding

s e c t i o n may be i n e f f e c t i v e i f t h e capture zone i s extremely small. The

present s e c t i o n presents a much s i m p l i f i e d example of such a case.

Consider a s p h e r i c a l s h e l l ( r a d i i R 1 < R 2 ) of a s i n g l e l i g h t

element which can be considered (perhaps by use of t r a n s p o r t cross

s e c t i o n s ) t o s c a t t e r neutrons i s o t r o p i c a l l y i n the laboratory system.

Suppose f u r t h e r t h a t a d e t e c t o r of very small radius R o < < R1 i s located

a t t h e c e n t e r of t h e hole (vacuum), and i t i s desired t o f i n d t h e d i s -

t r i b u t i o n i n t o energy groups of neutrons impinging on the d e t e c t o r , con-

sidered as a p e r f e c t absorber.

Specifically, l e t Eo > El> ...> EG denote the bounds of the

energy groups adopted, EG being t h e "cutoff" energy, below which neu-

t r o n s are l o s t t o a category LE. Consider a c o l l i s i o n C a t x, y, z,

R, R1 4 R < Re of a neutron with incident d i r e c t i o n u, v , w. The

fraction f of t o t a l s o l i d angle subtended a t C by t h e d e t e c t o r i s

1 (1 -
- c o s u ) ~1- ( l - (1 - L o 2 ) ) = L o 2 z L s i n 2 a = L R 2 / R 2 where u ,
2 2 2 4 4 4 0
i s defined by s i n u =Ro/R. Observe t h a t we a r e only considering a case

where t h e s e approximations a r e very good.

Let t h e d i r e c t i o n from C t o the o r i g i n be denoted by u", v",

w", and l e t r~ = u u" + v VI' + w w" be t h e cosi,ne of t h e angle between

t h i s d i r e c t i o n and t h e i n c i d e n t d i r e c t i o n . Finally, l e t E" be t h a t

new energy which would result from a s c a t t e r i n g from the d i r e c t i o n

u, v, w t o the direction u", v", w".

-136-
 We proceed in two essentially different ways according as this

energy E" is below or above EG. If below, we a l l o w the colliding

neutron of weight W to scatter isotropically in the laboratory system
into a new direction ut, v ' , wl. If the new line of flight cuts the

detector, or if it f a l l s outside of the solid angle subtended at C by

the detector and corresponds to a scattered energy below cutoff, we de-

posit the weight W in LE. Otherwise, we follow the scattered neutron

further to escape or its next collision by the usual (#?) routine. It

will be seen that this is the orthodox way of treating a collision,

except that we assume that directions within the very s m a l l detector
solid angle have the same behavior as the direction u", v", w" insofar
as resultant energy is concerned.
However, in case the energy E" corresponding to u", VI', w" is
above cutoff, we deterministically add a weight Wft to the category D
J'
which records neutrons impinging on the detector with energies in the

group containing the energy E". Here f has its assigned meaning and

t = exp {- (R - Rl)/h (E")} is the transmission for neutrons of energy

E" in the direction toward the center of the detector. We then allow an
isotropic scattering into the direction u', VI, w'. If this direction
falls within the counter solid angle, we force a first collision of weight
W(l - t) at x + u"1, y + v'lk, z + w"k, where I = - A (E")kn [l - r(1 - t)] .
(Cf. Chapter 111, e5.) If the scattered direction u'v'w' is outside

this solid angle, we determine the corresponding new energy E' and

deposit weight W in LE if E t is below the cutoff EG or follow

the weight W further in the usual ( p ) routine if E' is above EG'
 In a large number N of collisions of the second kind (E" > Ec)
the expected number hitting the detector w i l l be NWft (with the proper
energy distribution), as it should be. Moreover (Nf) neutrons will

scatter into the detector direction, on the average, resulting in a

t o t a l weight of (Nf) W (1 - t) = (NW) f (1 - t) having first collisions
in this direction, and these are distributed spatially in the correct
exponential distribution. Finally,an expected number N( 1 - f) of
neutrons will scatter outside the detector solid angle, with total weight

U(l - f) W = (NW)(1 - f), and these are correctly processed in their

subsequent history,

The method has the great virtue of deterministically contributing

a correct positive weight to the detector on every collision of the second

kind. (Collisions of the first kind (E"< EG) cannot do ao physically.)

A flow diagram covering the method is given in Fig. k.For the
cc formula,'onemay refer to $6 of the present chapter. The computational
parameter J) is set to zero at the source (6).
I
1 I
\
\

U"
V"
--
v
u
w + LE - LE

I I
w1'- w
I t = exp -(R - Rl)/hj

I WfttD-Dj
1
I
8
P
J
(3 W(1 - t) -w
y +vl- y

L
2 + WJ- 2

R =x +y + &

Fig. 44a

-139-
 18. Remarks on thermal neutrons. I n a l l neutron problems w i t h

which we are concerned, t h e energy range of neutrons a c t u a l l y followed

i s bounded from zero by some minimum energy E

*> 0. This bound may be

some s t i p u l a t e d energy E
* above t h e mean thermal energy of t h e medium,

neutrons with energies f a l l i n g below t h i s upon any c o l l i s i o n being of no

i n t e r e s t i n t h e p a r t i c u l a r problem. The energy

*
E then s e r v e s as t h e

lower bound EG of t h e lowest energy group f o r which c r o s s s e c t i o n s are

stored.

However, when neutrons are t o be followed down t o t h e mean thermal

energy of the medium, t h i s energy serves as t h e minimum E

*, and several

remarks are i n order. A l l c r o s s s e c t i o n considerations and s c a t t e r i n g

formulas have been based on t h e assumption of t a r g e t n u c l e i which are a t

rest i n t h e laboratory system. I t i s c l e a r t h a t our procedure i s un-

j u s t i f i e d i n t h a t part of t h e energy range approaching t h e thermal energy

i n t h e case of l i g h t elements. Moreover, i f t h e mean energy of the medium

i s i n t h e range below t h e molecular binding energies involved, simple

nuclear c r o s s s e c t i o n s may no longer be applicable. (11)

If neutrons reaching "thermal energy" are dropped, w e again have

E* = EG as above, neutrons w i t h energies below E* being thrown i n t o

a counter f o r losses t o cutoff.

When t h e problem n e c e s s i t a t e s a c t u a l l y following neutrons which

have reached t h e thermal energy range, w e r e q u i r e a lowermost energy

__ - __

-
~~

(l')Cf. S. Glasstone, M. C . Edlund, The Elements of Nuclear Reactor

Theory, D. Van Nostrand Company, Inc., Prince%&, N. J., 1952,
on thermal neutron cross s e c t i o n s .

-140-
group f o r such "thermal neutrons," a s u i t a b l e "average free path," and

d i f f e r e n t i a l s c a t t e r i n g laws f o r t h i s group, as well as p r o b a b i l i t i e s f o r

d i f f e r e n t types of r e a c t i o n s . It i s impractical t o deal with t h e actual

case of a d i s t r i b u t i o n of neutron energies i n t h e thermal range, and a l l

methods assune neutrons within t h i s group have a fixed energy and upon

elastic c o l l i s i o n r e t a i n this energy. For heavy elements, isotropy i n

t h e laboratory system seems t o be a v a l i d assumption f o r e l a s t i c s c a t t e r -

ing. However, t h e choice of proper averages f o r t h e other nuclear con-

s t a n t s referred t o depends on t h e elements involved and on complicated

questions of t h e a c t u a l energy d i s t r i b u t i o n s obtaining, which are non-

Maxwellian i n media with strong thermal capture.

- souTces.
19. Remark on determination of photon
I - C o l l i s i o n s of

neutrons with p a r t i c u l a r types o f n u c l e i m y r e s u l t i n t h e production of

y -rays; radiative capture and c e r t a i n (n-n) i n e l a s t i c processes are

examples. One of t h e a p p l i c a t i o n s of Monte Carlo t o neutronics problems

which i s becoming of increasing importance i s the determination of

photon production by neutrons i n shields. The l a t t e r i n f o m a t i o n m y be

used i n turn f o r t h e input energy and spatial d i s t r i b u t i o n i n t h e Monte

Carlo treatment of photon d i f f u s i o n , discussed i n t h e following chapter.

Insofar as the neutronfcs of such shielding media i s concerned, one

need only record as they occur t h e number o f radiative captures, and t h e

number o f (n-n)radiative c o l l i s i o n s caused by incident neutrons of energy

group g, i n each of a set of spatial zones.

 This, together with a knowledge of t h e Y energies produced i n

capture and the production cross sections f o r Y 's by neutrons of energy

group g i n (n-n) reactions, yield8 the desired d i s t r i b u t i o n s f o r the

source i n the related photon problem.

-142-
 cmPTER VI

PHOTON COLLISIONS

1. Introduction. The procedures discusssd in the present

manual are adapted to problems involving the scatkering of photons on
electrons and nuclei insofar as photons may be considered as "particles"

subject to mean free path and differential scattering laws. It is for

this reason that we speak of particles in the text, It is only in t h e
collision routine ( Y ) tha.t the physical character of the particles

need be distinguished. We have discussed many of the collision proc-

esses for neutrons with nuclei in Chapter V. We now turn to problems
involving photons. The principal types of photon-electron collisions
are: (1) Compton scattering, whioh is the strict analogue of elastic
scattering of neutrons on nuclei; ( 2 ) pair-production; (3) photo-
electric effect. These we proceed to discuss from the point of view
of the Monte Carlo technique.

2. Basic concepts and constants, We shall use

c = 2.997'76 x lolo em seccl for the velocity of light in vacuo

 h = 6.624 x lo-= erg sec for Planck's constant
m0 = 9.10658 x gm for the rest-mass of the electron, and

e = 4.8025 x l
o
-
' esu for the electronic charge
A photon is characterized (for the purposes of this chapter) by

its energy E (ergs). To such a photon is ascribed a frequency Y (sec'l)

by means of the relation e = hu and a wave-length A (cm) by the equa-

tion h v = C. The "equivalent mss" m (gm) of such a photon is defined

by m 2 = e . The speed of a l l photons in vacuum being c, we m y

assign to a photon a momentum (vector) mV, where V is the velocity

vector of the photon in the laboratory system. Then

IVl = c, and ImV( = mc = h v/c

An electron is characterized for OUT purposes by its charge

e (esu), rest-mass mo (gm), and its velocity vector V. Its mas is

then m = mo/ fi- where B= IV I

/c. The momentum vector of
the electron is mV and its total energy mc
2
.
It is customary to express photon energies by means of a dimen-
sionlesa parameter E = hv/moc2, which gives the ratio of' photon energy
to the rest-mass energy of the electron, namely, m c2 = .81837 x 10-6 erg
0
or .51083 MeV, recalling that 1 MeV is 1.60203 x lom6erg (Chapter 11,§3).
This parameter E is the one we adopt for photon energy in the reminder
of the chapter.

3. Compton collisions. A Compton collision is by definition

a collision of a photon with an electron (the latter assumed free and
a t r e s t i n t h e laboratory system), with preservation of t o t a l momentum

( v e c t o r ) and t o t a l energy. We may, therefore, proceed i n much t h e same

way as w e d i d i n t h e case of elastic c o l l i s i o n between neutron and

nucleus (Chapter V, § 3 ) , except that we do not introduce a center of

mas. We agree on the n o t a t i o n s and r e l a t i o n s of t h e following table:

Before c o l l i s i o n A f t e r collision
Fhoton Electron Photon Electron

2
Mass "1 = hv/c In2 = mo n 1 = hvl/c 2 n2 = "0 /Vl7
B = w2/c a

Velocity
v2 = 0 w1 w2

Speed lVll = Iv2\ = O IWll = lwzl = w2

lvbmentum
mlVl m2v2 nlWl n2W2
2 2 e ' = hv'
2 2
!lbtal energy e = hv = y c
m2c = nlc n2C

We have t h e two conservation l a w s :

m1V 1 + m2V2 = n1W 1 -t n2W2

m c 2 + m c 2 = n c2 + n c 2
1 2 1 2

o r equivalently,

t hv = hv' + c 2 (n2 - m2)

2
Introducing the energy parameters E = hv/moc2 and E ' = hv'/moc ,
we read
 EVl = E'W1 + W2/ (3)

E - E' = l/fi -P2 - 1 (4)

remembering that, contrary t o our usual p r a c t :e, t h e capitals E, E' '

r e p r e s e n t scalars.

From ( 4 ) it follows that

[l + (E - E')] - 1 = 2 (E - E') + (E - E')2 = P2/(1 - B 2) (5)

while f r o m ( 3 ) we o b t a i n

E%; - 2 EE' IV1l(W cosUfl + Ef%f = We2 /(1 - B 2)

where 1L is t h e d e f l e c t i o n angle f r o m V1 t o W1. Hence,

1

(E2 + El2 - 2 EE' cos Ufl) = (E - E')2 + 2 EE' (1 - cos #1)

= 8*/(1 - f12) (6)

2
Eliminating p /(1 - B 2) from the two equations ( 5 ) and ( 6 ),

2(E - E') + (E - E')2 = (E - + 2 EE'(1 - cos al) or

E - E' = EE' (1 - COS 2)

Thus we have finally t h e r e l a t i o n

E' = E/[ 1 + ~ (- a)]

i
 where we use a = coswY as usual f o r the laboratory d e f l e c t i o n cosine

of t h e i n c i d e n t p a r t i c l e . This i s t h e analogue of t h e neutronics rela-

tion E! = E(S +T ,u) of Chapter V, $3.

I t i s noted that, as i n the case of neutrons s c a t t e r i n g on

n u c l e i heavier than hydrogen, there i s a p o s i t i v e minimal energy

E 'min I=: E/1 + 2E f o r the s c a t t e r e d photon, a t t a i n e d f o r a = cos*

1
= -1.

A s a s i d e r e m k , observe that the r e l a t i o n

leads a t once t o the familiar wave-length r e l a t i o n

i f we f i r s t write

Ei- - E - 2 sine (2)

and use t h e r e l a t i o n s E = hv/moc 2 and u h = c f o r both primed and

unprimed variables E,Y,h. The "Compton wave length" A. i s defined

t o be h/moc = ,024264 a (18 = angstrom u n i t = c m ) , and the

I r e l a t i o n i s sometimes w r i t t e n

We have seen how EL = cos% determines E ' , the r e s u l t i n g energy

of the photon. We have s t i l l t o show how the k i n e t i c energy and

defection cosine cos*2 of the s c a t t e r e d e l e c t r o n are found i n terms

of E and E'.

-3.47 -
 To be sure, equation ( 2 ) answers t h e first question, s i n c e

c 2(n2 - m2) = hu - hv' = moc2( E - E ' ) ; thus t h e k i n e t i c energy of the

scattered electron i s
1
II
Ee = .51083 (E - E ' ) Mev

We use equation ( 3 ) t o determine cos #2 a e follows:

= v l -B (EV: - E ' lVLl lWl[ cos Gl)/w,c

= \JL- (E - Efa)/B

where a/ 1/ 1 - B2 i s known from e i t h e r equation ( 5 ) o r ( 6 ).

SLurmrarizing t h e r e s u l t s of t h i s s e c t i o n , we have

photon d e f l e c t i o n cosine a = cos I P ~

s c a t t e r e d photon energy E ' = E/ [ 1 + E(l - a )1 (units of moc 2)

r e c o i l e l e c t r o n K.E.
Ee
= .51083 (E - E ' ) Mev

electron deflection cos G2 = (E - B'a,)/(i?/Vl -P2)

A s an example, suppose an 8 Mev photon s c a t t e r s a t an angle of

@ = 60'. Then
2

L.
 a = 1
v

E = 15.661 (8 MeV)
E' = 1.7735 ( 4 0 5 9 6 MeV)

Ee = 7.09404 Mev
8/ V l T = 14.854
COS
2
= .99464
'v2 = 5'56'

4. The Klein-Nishina d i f f e r e n t i a l c r o s s s e c t i o n . (u) This i s

t h e d i f f e r e n t i a l c r o s s s e c t i o n uE( Q ) i n cm2/steradian which governs

t h e d i s t r i b u t i o n of a = cos J /
1 f o r t h e s c a t t e r e d photon i n the Compton
e f f e c t , and i s t h e photon analogue of t h e neutron c r o s s s e c t i o n of

Chapter V, 954 and 6. It i s defined by t h e formula

bE( Q)dP = E 2(1 - a )2

2
(1 + a ) [ L + ~ ( -1

where (13)r = e2/moc2 = .28183 x cm.

0

Now w e might proceed j u s t as i n t h e case of e l a s t i c c o l l i s i o n

t o determine a by (say) yon Newnann ' s device (Chapter I, 55) and then
E ' from E and a using t h e r e l a t i o n

('*)For an account of t h i s and r e l a t e d functions see R . Latter and

H. Kahn, Gam-Ray Absorption Coefficients, P r o j e c t RAND, R-170 (1949)
_I_-

and t h e National Bureau of Strrndards Circular 542, Graphs of the Compton

--
Energy-Anale Relationshie and t h e Klein-Nishina Formula from 10 Kev t o
500 MeV.
----
(13)Thia i s t h e c l a s s i c a l Thomson r a d i u s of the2electron which n t e r s
i n t o t h e e l e c t r o n cross-section f o m u l a (8/3)m = ,6654 x 10-2E cm2
0
( = .6654 barns).
of the preceding section. In view of the complexity of (7)it is

fortunate that a simpler alternative exists. We describe this in the

following section.

5 . The photon energy distribution and Compton cross section.

In order to prepare for the alternative method, we first determine the
differential cross section ZE(Et) from the relation

UE(E')dE' = uE( Q ) d Q =-OE(Q ) 2 .(a)

Using the preceding equation (8)in the form

a = 1 + E1 - E1 t

we o b t a i n the result

We m y obtain easily from this the t o t a l cross section f o r

Compton scattering as follows:

= 2xro

-150-
which is graphed in the National Bureau of Standards Circular 542

( IC and m y be used to compute free paths A = l/Nebg (Compton)

t3
for a suitable set of energy groups (assuming Ne is the numerical

density of free electrons and no other processes exist).

We will actually use the Monte Carlo procedure

to determine E' from E and random number r, since a sufficiently

accurate fit for the inverse function is given by(14)

where A. = E/(1 .I- ,5625 E), and E 5 4 (-2 Mev). Addition of a term

1
-
2
(E - 4) r2 (1 - r)
2

yield's a reasonably good fit on the range 4< E 6 10.

I The equation ( 5 ) above then permits determination of a from
E and E ' .

Thus a Compton collision with incident energy E not exceeding

5 MeV may be satisfactorily handled by the flow dLa@;lramof Fig. 45.

(14)Bengt Cwlson, The Monte Carlo Method Applied to a Problem in

Y-ray Diffusion L o X b F S a t i f i c L a b o r a t o r y , ~ E ? % J ~ J ~ 3 .
Cf. H. Mayer,T: A . Burton, Tables -of -the Compton Effect Cross Sections
-
."
and Energies, LAMS-1199,Los Alamos Scientific Laboratory, 1953,for
goodness of fit.

-151-
s = E/(1 + .E625 E) bl r
I
E' = E/ 1 + sr + (2E - s) r3} I

Fig. 45

-152-
 6. Photoelectric effect and pair production. We have discussed

the Compton collision for photons with electrons and noted that it
corresponds to the elastic collision in neutronics. If a medium con-
sists of elements A, B , . . . , the contribution of the Cornpton effect

to the "total cross section" (cf. Chapter I I I , § 2 ) is NoE(Compton),

where N = NAZA + NBZB + ...
i B the total numerical electron density.

Here NA denotes the number of atoms of element A per cm3 and Z

A
i8 the atomic number (= number of electrons per atom) of element A ,
In determining the "total cross section''all collision processes
must be taken into account, and one must consider in this connection the

further contributions NAUE

A
(pe) .t ... (pe) of the photoelectric
E

A (pp) +
effect and NAoE ... z zE E
(pp) of pair production. Cross sections
2
f o r these processes for elements with 4 Z 5 92 and .10 $ E 20 (moc )
may be found in the RAND report R-170 referred to before. Thus the

free path determining collision positions is given by

= l/Z
where

In the present report we regard pair production and photoelectric

effect as absorptions, and have therefore, upon any collision, an

absorption probability [ - NoE(Compton)]/z , which m y be dealt

with by the alternative methods discussed in Chapter 11, 52.

-153-
 It must be remembered t h a t secondary photons (X-rays) may r e s u l t

from t h e p h o t o e l e c t r i c e f f e c t i n heavy elements, and t h a t t h e e l e c t r o n

and p o s i t r o n born i n pair production produce bremsstrahlung while being

slowed down by i o n i z a t i o n of t h e medium. Moreover, t h e inverse photo-

e l e c t r i c e f f e c t may come i n t o play, and t h e p o s i t r o n i s f i n a l l y annihi-

l a t e d with an e l e c t r o n t o produce s t i l l f u r t h e r y-rays. If t h e energy

cutoff of t h e problem i s not s u f f i c i e n t l y high t o exclude consideration

of these secondary photons, one must include i n t h e output of t h e o r i g -

i n a l problem s u f f i c i e n t information t o determine source d i s t r i b u t i o n s

f o r such secondary r a d i a t i o n . W e do not d e a l with t h i s case i n t h e

present manual.

-154-
 CHAPTER VI1

DIRECTION PARAMETERS AFTER COLLISION

1. Introduction. The purpose of the present chapter is to

develop formulas for the final direction parameters of a particle after

scattering through an angle *l of cosine a, in the laboratory system,
from an incident direction M, v, w (or w alone). While these formu-
las are somewhat long, they my be derived from the simple principle of
elementary complex variables which states that (x + iy)(cos €I+ i sin 0)
is a complex number whose vector is rotated through an angle 0 from that

of x + iy.

- - -
2. Formulas for the final direction cosines, Let u = cos 01, v =

cos $, w = cos -? be the direction cosines of the incident line of flight.

Consider (5,7 , w) as a point on the unit sphere u2 + v2 + w2 = 1 in
direction space U, V, W
-
(a?.Fig. 46). Letting P 7 be the polar
coordinates in the U, V plane of the point (u,v, 0) we see that

-155-
 sin = 7,’;
Our f i r s t o b j e c t i v e is t o derive formulas f o r a p a r t i c l e r o t a t i o n

of space i n t o i t s e l f which t a k e s the point (0,0,1) i n t o (u,v,w);

---
U,V,W
namely, t h a t r e s u l t i n g from i t e r a t i o n of two simple r o t a t i o n s : the
-
f i r s t about the V-axis through the angle Y , followed by a second

r o t a t i o n through t h e angle 7 about t h e W-axis.

U
Fig. 46

-156-
 We have for t h e first r o t a t i o n

w' + i u ' = (w e i u ) ( c o s 'i; + i sin ?)

v' = v
and f o r the second

ufl + ivll = (ul + iv')(cos .c i s i n 7)

w" = w '

Separating real and i m g i n a r y parts, we read

u ' = u cos 'i; +w sin T

v' = V

w' = -u s i n F +w c o s 'G

utl = u 1 cos - Y' s i n 7

VI' = ut s i n + VI cos (4)
w" = W'

S u b s t i t u t i o n o f ( 3 ) i n t o (4) y i e l d s

ult = u cos 7 cos g- v sin + w sin T COS

VI' = u cos F sin

w'' = -u s i n 'i; -t w cos r

Replacing t h e functions of 7, by t h e i r values i n terms of

-P = {v),
u, v, and G , we obtain from (5)

-157-
 1
u" =
--
(uwu-vv>/;+w~

V"

Under t h i s r o t a t i o n (6) it i s evident t h a t (u,v,w> = (0, 0, 1)goes over

Now we know that t h e d i r e c t i o n ( 5 , 7, G) of the d e f l e c t e d l i n e of

f l i g h t makes a n angle Q1 of cosine

- -
a w i t h t h e d i r e c t i o n (u, v, G)
of t h e i n c i d e n t l i n e of f l i g h t , It i s c l e a r t h a t t h i s r e s t r i c t s t h e

former t o a.cone of d i r e c t i o n s of opening about the l a t t e r d i r e c t i o n ,

?Pl
and t h e d i r e c t i o n s on t h i s cone are a l l equally likely. W e adopt t h e

following m b i t r a r y convention f o r f i x i n g the d e f l e c t e d d i r e c t i o n ;

namely, we always consider f i r s t a cone of opening w1 about the

W-axis OW, and a point P = ( s i n JI

1
cos d , s i n ?Pls i n 6 , cos wl)
located on t h i s cone and the u n i t sphere, determined by an azimuthal

angle of d uniformly d i s t r i b u t e d on - A 5 d 5 A (Fig. 47).

-158-
 1

/ Fig. 47

We then agree t h a t the final d i r e c t i o n (E, f, G ) is t h e image

-c_-

under r o t a t i o n ( 6 ) ---
of the point P -
so constructed. It should be

clear t h a t t h i s i s an appropriate convention f o r OUT purpose.

For s i m p l i c i t y we adopt the further notation

-159-
I n t h i s not.ation, t h e point P = (bc, bd, a ) = (u, v, w ) goes over i n t o

Summarizing the r e s u l t s of the present section, we have

w t = -bc v-*-
(1 -w ) + aw J
where

u, v, w a r e d i r e c t i o n cosines of t h e incident l i n e of f l i g h t
u’, v), w * a r e d i r e c t i o n cosines of t h e d e f i e c t e d l i n e of f l i g h t (lab.)

a = cos *,,where V, i s t h e angle between t h e two

C = COS 6 , where - 1~ I; d (= R uniformly

d = (sgn d )VI
7

- c2
Note t h a t t h e s e formulas should not be used if Iw I is too close t o

unity. They are then poor computationally and indeterminate a t l w l = 1.

I n such a case it i s b e t t e r t o by-pass t h e r o t a t i o n ( 6 ) and use f o r

-160 -
u t , v t , w t t h e coordinates

U' = bc

V' = bd

w' = aw

where a, b , c, d have t h e d e f i n i t i o n s above.

3. Subroutine f o r t h e f i n a l d i r e c t i o n cosines. We incorporate

t h e s e r e s u l t s i n t o t h e ( 8 ) routine of Fig. 48, The power 2-n used t o

determine whether t h e incident d i r e c t i o n should be considered v e r t i c a l

may w e l l depend on t h e accuracy required and on t h e number of " b i t s "

afforded by t h e particular machine being used,

-161-
 1-1
(1 - Iwl} - 2-"

Ul

vl

Wl =
=

= 1
Eqs (V)
U' =

v' = bd
bc

w f = aw

d
Fig. 48

-162-
 4. Final direction w -in s l a b or s p h e r i c a l l y s;mmetri.c case.
The present section i s concerned with determination of t h e f i n a l di.rectlon

cosine w' = cos Y 1 of t h e s c a t t e r e d l i n e of f l i g h t as it e x i s t s a t t h e

point of c o l l i s i o n , i n terms of t h e cosine a cf t h e laboratory angle Wl

determined i n ( Y ), as required i.n both slab geometry and s p h e r i c a l l y

symmetric geometry. (Cf. Chapter 11, 92, and Chapter I V , 5 8 . ) We r e c a l l

t h a t , i n t h e s p h e r i c a l case, t h e cosine w = cos y of t h e angle y which t h e

incident l h e o f f l i g h t makes with t h e radi.us vector - it e x i s t s -

as - at -
the

--
point of collj.sj.on has a l r e a d y been prepared (Chapter ITJ, $ 3 ) and s t o r e d a s

w before e n t r y a t (a).
It i s immediately evldent t h a t t h e w) formula o f equations (7), i n

42 of the present chapter$

w' = -bcvl - w2 + aw

I a p p l i e s t o t h e slab case, where a, b , and c have t h e d e f i n i t i o n s g-iven t h e r e .

I
Note t h a t t h e formula i s independent of t h e u t v used i.n t h e d e r i v a t i o n , as

it should be.
But it w i l l a l s o be c l e a r t h a t t h e i d e n t i c a l formula a p p l i e s t o t h e

s p h e r i c a l case, i f we imagine f o r t h e moment t h a t OW r e p r e s e n t s t h e outward1.y

d i r e c t e d radius vector d i r e c t i o n and u, v, w the d j r e c t i o n of t h e incident

l i n e of f l i g h t r e l a t i v e t o t h e radius vector, t h e u,v being immaterial,

I n both a p p l i c a t i o n s , it w i l l be observed t h a t t h e angle 6 involved

in c = cos d may be l i m i t e d t o t h e range 0 (= a 7r, since cos -6 = COS 6= c

i s t h e only functjon of d appearing i n t h e w' formula.

-163-
 0
1 w' = -bc (1 - w ) + aw

t
Fig. 49

I n both cases,therefore, we have t h e simple r o u t i n e of Fig. 49.

5, S c a t t e r i n g i s o t r o p i c i n t h e laboratory system. In any case

where p a r t i c l e s emerge from a c o l l i s i o n i n a d i s t r i b u t i o n i s o t r o p i c i n

t h e laboratory system, it i s p a t e n t l y f o o l i s h t o determine a = 2r -1

f o r the d e f l e c t i o n cosine r e l a t i v e t o t h e l i n e of f l i g h t and then determine

t h e . f i n a 1 d i r e c t i o n cosines by means of t h e preceding two s e c t i o n s . It

i s simpler t o assign t h e d i r e c t i o n of t h e d e f l e c t e d l i n e of f l i g h t by t h e

same method we used t o set up an i s o t r o p i c source i n Chapter 11, $5a, c,

ignoring t h e incident d i r e c t i o n completely. Thus t h e ( r ) routine referred

t o i s simply t h a t of Fig. 10, with e x i t t o ( fi 3 for t h e u, v, w case,and,

still more simply, boxes r and w = 2r - 1 f o r t h e slab and s p h e r i c a l case.

 CHAPTER VI11

TERMINAL CLASSIFICATION

1. -
Introduction. We have already indicated various terminal

events in the history of a particle, e.g., transmission, capture, loss

to enerey cutoff, and so on. We must still consider some of the many
ways in which it may be desired to classify a particle which escapes
the system. This is the general function of the ( 6 ) routines pre-

viously referred to. Such routines invariably exit to (a) for the
introduction of a new source particle. Of these many kinds of classi-

fication, we can hope to give only a brief indication, since the demands
of physicists on this score are frequently inv.olved and exacting, Indeed,
it is the ability of Monte Carlo methods to provide answers to the most

intricate questions of this kind which makes them an indispensable tool

in design,

2. Classification of escapes on number of collisions. A s the

simplest example, we may refer to a problem whose purpose it is to

determine the distribution of escaping neutrons with respect to the

number of collisions suffered within the system. Suppose we agree to

keep storage registers N,, v = 1, 2, ..., 10 for the total weight of

neutrons escaping with v 5 10 c o l l i s i o n s , and Nll f o r those escaping

with v 2 11. We suppose t h a t the forced f i r s t c o l l i s i o n device has been

used so t h a t t h e transmission counter T records t h e weight escaping

with v = 0 automatically ( c f , Ch. 111, e4). We may, t h e r e f o r e , follow

t h e r o u t i n e of Fig. 50.

W + N - N
V V

-+ Nll - N1l

Fig. 50
3. Energy and angle d i s t r i b u t i o n s of escape. Consider a problem

w i t h a source d i r e c t i o n u = 0, v = 0, w = 1 i n which one desires t h e

c o r r e l a t e d d i s t r i b u t i o n of escapes i n energy and angle with t h e source

direction. There must be provided i n permanent s t o r a g e a s u i t a b l e set

of cosines C1 > C2 > ... > CJ = -1 and energy bounds > ~2>*..>~
H
= EG, which need not be t h e same as those used f o r o t h e r purposes, while

i n dynamic s t o r a g e , we reserve JH p o s i t i o n s N f o r t h e t o t a l weights

j,h
escaping i n t h e corresponding categories. I n such a problem, an escaping

neutron e n t e r s (e) w i t h a known weight W, d i r e c t i o n cosine w, and

energy E, and i s c l a s s i f i e d as indicated i n Fig, 51.

-166-
 Fig. 5 1

Such a classification in angle is appropriate in a case like that

I
of Fig. 52,where a detector band is at essentially infinite distance

from a scattering medium symmetric about the Z-axis. Here infinite

distance means that all particles escaping from any point of the

medium surface in the same direction hit the band in the same angular

zone. Then the number N referred t o is the number of h-group

J ,h
escapes in the solid angle 2 R (C
J - CJ-1 ), so that
S,h
N
/ 2 f l (CJ CjWL) -
is the number in this category escaping per steradian, and the number

hitting a detector of given area on the circle indicated in the figure

can be predicted.
If the source direction is again u = 0, v = 0 , w = 1, and the
scattering medium is symmetric about the Z-axis ensuring symmetry of
escape about this axis, but the detector is not at infinity in the

-167-
 Z

I 180"

Fig. 52

-168-
sense r e f e r r e d t o , it i s necessary t o take i n t o account t h e d i s t a n c e
% from center 0 t o t h e d e t e c t o r band as i n d i c a t e d i n Fig. 53 and
Fig. 54. I n the l a t t e r , the e x i t leads t o t h e previous c l a s s i f i c a t i o n

r o u t i n e of Fig. 51 and thence t o (a). A t e n t r y t o t h e r o u t i n e of

Fig. 54, the x,y,z parameters are those of the last point of de-

p a r t u r e before escape. Again t h e r e s u l t s may be normed t o number p e r

s t e r a d i a n and i n t e r p r e t e d with reference t o t h e d e t e c t o r band indicated.

Z

I O0

I
I

I
I
I

:s

Fig, 53 180'

-169-
 2 2 2
@-+ 6=ux+vy+wz- A = 6" -(x + y + z

To Fig. 51

Fig. 54

The preceding discussion is limited to systems symmetrically

distributed about the source (+Z) direction. Problems of this kind

have the great advantage, from the Monte Carlo standpoint, that no
escaping particles are lost to classification.
Unfortunately, there are experimental considerations which, in
many cases, dictate a geometry in which such symmetry is lacking.
Consider, f o r example, the situation in Fig. 55, in which a source
with direction u = 0, v = 0, w = 1 impinges on the lateral surface

of a cylinder with axis on the X-axis. It is desired to count the

numbers N hitting a coaxial band zone j upon escape, the de-

J
tector band having radius dg and height h. Since the cylinder

-170-
Z

I
I
I
I
I
I
4
I
I S
I
I
180'

Fig. 55
lacks symmetry about the Z-axis, it is no longer possible t o classify
escaping particles according to the cosine w of the angle Y of the

line from origin to the point of intersection of the line of flight

with the sphere of radius $, since the escapes with this direction

are non-uniform about the Z-axis. Without rather complicated devices

for prejudicing scattering in the band direction, which we shall not

discuss, all one can do is to submit to running very much larger
samples in order that the rather small solid angle actually subtended
by the band shall receive a sufficient number of escapes to render the
N statistics reliable. This is a real difficulty, and it would be
J
of great value for experimentalists contemplating the use of Monte Carlo
methods f o r thick target corrections to consider the feasibility of the
symmetric geometry in the experimental setup.
A suitable routine for the case cited, without any special

importance sampling devices, is given in Fig. 56, the exit again re-
ferring to the routine of F i g . 51. The category % refers to all
escapes failing to hit the band. In computing the distance t to the
(infinite) cylindrical surface y
2
-f z2 = dg,
2
one avoids the 1 - u2 C 0

catastrophe,as indicated in Ch. IV, 86 (end of section).

Consider finally a case where it is desired to classify a
particle escaping w€th direction u,v,w from the Lateral surface of
2
a cylinder x + y2 = R2 according to the angle measured from the
normal to the surface. If (xt,yt,zt) is the point of escape, deter-
mined in the usual manner from u,v,w,R, and the position x,y,z of

-172-
 x =x+ut
t -

-
t
zddB w

+I
To Fig. 51

Fig. 56

-173-
l a s t departure, then t h e d i r e c t i o n of t h e normal a t t h i s point is given

by

xt - 0
R ’
Yt - 0’
R
Zt
R
- Zt

and t h e cosine w of the angle of escape w i t h t h e normal i s

x u + y v + z t . O
w = t t
R = (XtU -1- YtV)/R
 REMARKS ON COMPUTATION

1. Scaling. If the machine employed l a c k s the " f l o a t i n g decimal"

f e a t u r e o r i f considerations of machine speed exclude i t s use, a l l

computations indicated i n t h e flow diagram a t each step must be planned

i n advance t o remain on the i n t e r v a l - 1 C x < 1. The formulas and flow

diagrams of t h e present t e x t are unscaled and a r e expressed i n the usual

physical units, The s c a l i n g procedure i n kbnte Carlo problems usually

causes 190 d i f f i c u l t i e s , and may be l e f t t o the imagination, with the

reminder that one must be c a r e f u l t o avoid the l o s s of accuracy which

attends over-scaling, Remarks are made about s c a l i n g a t various points

i n the present chapter as it enters i n t o s p e c i a l subroutines, such as

exp( - x / y ) where x may exceed y, i n t h e J n x r o u t i n e where x

may be as small as 2'23, and i n the cosine r o u t i n e f o r - fl 5- x 5

- *.

2. Debugging. When a problem i s ready for actual running by

the machine, some method of d e t e c t i n g the i n e v i t a b l e human e r r o r s which

m y enter a t a l l phases of the coding must be employed t o assure a l l

concerned t h a t the machine w i l l perform exactly the routine intended.

It goes without saying that a l l permanent storage constants should be

printed out a t various i n t e r v a l s and checked against the o r i g i n a l list,

c e r t a i n l y a t the beginning and end of the problem. b r e o v e r , i t is

sometimes possible t o b u i l d i n a sum check of a l l code orders and

permanent constants which will automatically detect any changes due t o

electronic misbehavior once the machine is running. e refer here

W

rather t o the question of ensuring that the code itself i s correct.

This question presents r a t h e r greater d i f f i c u l t i e s i n Monte Carlo

problems than i n other types of computation because of the complexity

of. decisions involved.

We have found t h a t the most convincing guarantee of faithfulness

of code t o flow diagram consists i n the preparation by hand of a

deterministic sequence of nurnbers rl, r2, r3, ... designed t o process

automatically a hand-picked set of p a r t i c l e s which are so chosen t h a t

every l o g i c a l l y possible path through the f l o w diagram is traversed a t

least once with non-trivial parameters. A special debug routine m y be

used f o r t h i s purpose which i n s t r u c t s the machine t o follow the code

precisely except that, upon c a l l f o r the next random number, the

random number routine i s by-passed, and the prepared l i s t i s consulted

f o r the next number rn of the hand-picked sequence. A t every l i s t

consultation the machine p r i n t s out a l l s i g n i f i c a n t variables of the

problem, the printout then being compared w i t h the hand-computed

p a r t i c l e h i s t o r ie s
 This i s a laborious process i n t h e hand-computing phase, b u t i s

w e l l worth doing, even though, i n complicated problems, the processing

of 50 p a r t i c l e s or so may be necessary.

3. S p e c i a l subroutines. W
e include a number of s p e c i a l sub-

r o u t i n e s f o r some o f t h e functions commonly occurring i n Monte Carlo

calculations. No claim i s made that these are the b e s t a v a i l a b l e ; they

are only given f o r completeness i n case the reader knows of no b e t t e r

one.

a. A random number r o u t i n e . We have used i n a l l problems t h e

sequence rlj r2, r3, ... of 38 b i t diadic numbers, generated by the

algorithm defining rn+l t o be t h a t 38 b i t nuniber corresponding t o t h e

middle 38 b i t s of the square of rn' where r 1 is the 38 b i t number
defined by 10 BBB FA4DE i n the system with base 16, with t h e decimal

after the f i r s t b i t . Thus rl = 0.001, 0000, 1011, 1011, 1011, 1111,

1010, 0100, 1101, lllfi,t h e f i n a l zero being ignored. The sequence

automatically terminates i n zero a t about n = 750,000 and has been

thoroughly tested, not only f o r the usual s t a t i s t i c a l features b u t by

its a c t u a l use i n many problems. The Monte Carlo use of t h e sequence i s

of a very peculiar and unpredictable kind. If one fixes a t t e n t i o n on

any p a r t i c u l a r random number box of the f l o w diagram, it appears that t h e

nurnbers a c t u a l l y s e l e c t e d from the sequence f o r use i n t h i s box are a

subsequence s e l e c t e d from t h e main sequence by completely i n s c r u t a b l e

r u l e s depending on t h e sequence i t s e l f . Whenever a p r o b a b i l i t y of some

terminal event (e.g., transmission) was known, the corresponding Ni/N

-177-
 f o r t h a t category agreed even b e t t e r than one might expect i t t o on the

basis of a t r d y random sequencer

b. Shifted random numbeE. An average problem may involve

20,000 source particle6,each with perhaps f i v e collisions, and each

c o l l i s i o n may involvethree o r four random numbers. It i s c l e a r t h a t

one may expect t o exhaust the random nmiber sequence i n many problems

i f it i s used i n t h i s simple way. It i s t r u e t h a t the same sequence

may safely be used over a g a i n without r e p e t i t i o n of results provided

the i n i t i a l rl i s not used a t the i d e n t i c a l place it was first called

i n the flow diagram, and indeed, such re-use of the sequence m y be

resorted to.

We may IIlention however t h a t the usual problems do not require

anything l i k e 38 b i t resolution, and the number of random numbers

available may be increased by a f a c t o r 2 o r even 4 by a s h i f t i n g

routine, which i s incidentally f a s t e r than squaring and so helps t o

i reduce machine t i m e . This consists i n using i n sequence not only

rn,
2 rn, 2a2 (2alr ), etc.,as random
OL1
but a l s o the f r a c t i o n a l parts of
n
numbers, before squaring rn t o obtain the next rn+l. For example,
if three random numbers a r e t o be obtained from rn, one might use the
f r a c t i o n a l p a r t s of rn, 212rn, and 212(212rn). The sequences ob-

tained by such extensions of the basic routine have a l s o been t e s t e d

and proved satisfactory.

-178-
 c. A Logarithm routine. From the u s u a l series

&(I+ x ) = x - x2/2 + x3/3 - x4/4 + ... 1x1 < 1

one obtains

,tn(l - x > = - x - x2/2 - x3/3 - x4/4 - 1x1 L 1

and hence

which, under the transformation -

y = 1-x becomes
l+x

~
Convergence of the l a t t e r series i s rapid enough f o r the range

-12 < y < 1 t o permit use of a moderate number of terms; f i v e should

s u f f i c e f o r most purposes.
Now suppose 6 is a nurnber on the range 2-23 6 c 1 for

which j n t, i s t o be computed. e determine t h a t integer

W n for

which
 Then 2-1 5 2
n
6 TI < 1. and In C = - n In 2 + I n TI ,where
In 9 can be obtained from the s e r i e s .

Now - ,& 2'23 = 23 i n 2 < 23(.694) < 16, so that In 5 I < 1

for S on the range - 4 < 1. Hence a r o u t i n e suitable f o r a

scaled problem is given by Fig. 57.

O u r only use f o r the &I function has been i n c o l l i s i o n -

d i s t a n c e formulas such as

k' = - a tn r = - z4 A r)

Limikation of random numbers e n t e r i n g t h i s formula t o t h e range r &20z3

means t h a t we force a c o l l i s i o n w i t h i n 16 free paths. This is inaccurate
t o the e x t e n t that e-16 i s the chance of c o l l i s i o n a t a s t i l l f u r t h e r

distance, which is an acceptable e r r o r i n most cases. I n problems i n -

volving systems so l a r g e t h a t such e r r o r s are s i g n i f i c a n t , probably some

type of importance sampling w i l l be used i n place of t h e simple formula

above. I n any case it i s c l e a r how t o modify the routine f o r higher

powers, 2'23 being q u i t e a r b i t r a r y .

For the r o u t i n e as given, I t i s c l e a r that i f a l l 'distances (cm)

are s c a l e d down by a f a c t o r exceeding 16 times t h e greatest free path

A occurring i n t h e problem, one i s s a f e so far as t h i s p a r t i c u l a r

formula i s concerned.

-180-
1
-,Pnq=g+g
2

t
Fig. 57
 d. The exponential exp(- x / y ) . It i s frequently necessary t o

compute the value of exp(- x/y), where 0 6. x < 1, 0 < y < 1. Even
though x and y are scaled, the quotient x/y may exceed unity and

cannot be d e a l t w i t h d i r e c t l y by fixed decimal point machines. There

are excellent polynomial approximations(lg) t o exp(- p for o < p 5 c ~ n2,

so t h a t we may define the integer n by

and p on the prescribed range by P = (x/y) - n l n 2. We shall

therefore have

where e-' S a
. + P (al f pa2) with

a = l
0

al =
- 9664279798

a2 = 03535763634

is usually s u f f i c i e n t l y accurate for our purposes.

( 15)B. Carlson and M. Goldstein,"Rational &proximation of Iknctions " ,

Los Alamos S c i e n t i f i c Laboratory; LA-1943, 1955 Gives polynomials
of various degrees with bounds on the errors involved.

-182-
 e may t h e r e f o r e proceed a5 in Fig. 58.
W

e. A cosine routine. We i n d i c a t e a scaled r o u t i n e f o r computing

-12 cos & , where 6 i s t o be chosen uniformly on - R

- b g
2 R , as
required by Ch. V I I , 5 3. It i s based on t h e following s e r i e s of

formulas

Input 0 sx < 1, 0 < y < 1

Stored constants ao, al, a2,an 2

Fig. 58

-183-
 -1 cos = 22{2-3- s2 + s4}
2

the f i n a l formula a r i s i n g from

cos6 = COS 4af = 2 COS

2
2 a' - 1

= 1 - 0 s2 4.0 s4

It may perhaps be mentioned here that f o r scaled machines t h e

following procedure i s somewhat s h o r t e r than t h a t obtained by applying

the preceding method t o Fig. 49, the r o u t i n e f o r slab and s p h e r i c a l

f i n a l direction w'. I n place of choosing 6 uniformly on 0 5 6 6 r,
we may introduce the angle
N

d = d - "/2, choosing the l a t t e r uniformly

5
CI)

on -r/2 d r/2. Thus we may proceed as follows:

 d'

S s i n 6 ' as above

N
C cos d = s i n ,j = s i n 2 6 ' = 2 s i n d ' cos 6 '

aw - bc {-

We should t h e r e f o r e have i n place of Fig. 49 the r o u t i n e of

Fig. 59. I n t h i s way we avoid t h e square root f o r b = )!=and

make use of only one double angle formula instead of two.

Fig. 59

-185-
 4, A Monte Carlo device f o r 6,
There a r e many ingenious

devices f o r avoiding the computation of t h e value of various functions

with a r g w n t s involving a random number. A s a simple example, (16) the

process

involved i n throwing f o r the r a d i u s of a neutron i n a source d i s t r i b u t e d

uniformly i n area on a d i s k of radius may be replaced by the f o l -

R1
lowing :

That is t o say, one may use i n place of the square r o o t of one

random number the l a r g e r of any two successive ones, for t h e equation

4 = results from t h e standard Monte Carlo r e l a t i o n

and i s thus equivalent t o choosing 5 on t--e i n t e r v a l (x, x + d x ) with

frequency 2x dx. ,
But the a l t e r n a t i v e of throwing points (6 v ) into

the u n i t square 0 5 6 5 1, 0 5v5 1 uniformly i n area and choosing

the l a r g e r of the coordinates automatically determines a n x B max(6,v)

(16 'Cf. J. W. B u t l e r ' s paper i n Symposium on Monte Carlo Methods,

,
John Wiley & Sons, Inc. New York, 1956.
-_I_-

-186-
on ( 6 , 6 + d t ) with t h i s frequency, as i s apparent from Fig. 59a.

5. A Monte Carlo device for the cosine of an e q u i d i s t r i b u t e d

angle! 17) In various places we have needed t h e cosine and s i n e of a n

angle 0 e q u i d i s t r i b u t e d on the i n t e r v a l - T 50 7r, and have hereto-

fore given a r a t h e r cwnbersom b u t straightforward method involving

s e r i e s , square r o o t , and, i n case of fixed decimal r e q u i r e F n t s , a

multiple angle formula. We give i n Fig 5% an a l t e r n a t i v e method,

e a s i l y scaled i f necessary and avoiding series and square r o o t r o u t i n e s .

One observes t h a t t h e d e s i r e d determination of c = cos4 and

d = sin 4 i s equivalent t o choosing a point (c,d) uniformly on the u n i t

circle. This i n turn Is tantamount t o throwing p o i n t s (t,'~) uniformly

i n t h e square -1 6 x, y 2 1, r e j e c t i n g those l y i n g outside t h e u n i t

circle x
2
f y2 = 1, and taking f o r c and d the v a l u e s 5 /bq
and for the r e t a i n e d p o i n t s
q//m, (6, 'I) are uniforrt& d is tri-
buted i n area i n t h e u n i t c i r c l e and hence t h e i r p r o j e c t i o n s (c,d) on t h e

u n i t c i r c l e are uniformly d i s t r i b u t e d a l s o . But these square r o o t s may

be avoided by l i m i t i n g the s e l e c t i o n procedure t o the first quadrant and

using double angle formulas t o obtain ( c , d ) uniformly on the upper h a l f -

circle. F i n a l l y a n a d d i t i o n a l random nuniber can be used t o change t h e

s i g n of d with p r o b a b i l i t y 1/2. The e f f i c i e n c y of the r e t e n t i o n of

random number pairs is m/4.

( 17)The method i s described by von Neumann i n U, S, Department of

Commerce, National Bureau of Standards, Applied Mathematics S e r i e s #12,
--
Monte Carlo Method, Washington, Do C., 1951. Cf, footnote (l),

-187-
 Shaded area = 2x du

0 z

4
Fig. 59b

-188-
 CHAPTER X

STATISTICAL CONSIDERATIONS

1. The l i m i t theorem i n t h e Bernoulli case. Suppose t h a t a

c e r t a i n experiment can result i n c ways Of which a are favorable

t o an event
1
E while b = c -
a r e s u l t i n the event E
0
(not El).

-
Consider t h e set of all sequences of N trials of t h i s experiment.
N
I n t h i s set of c sequences, the number of sequences r e s u l t i n g i n

exactly M occurrences of E1 i s c l e a r l y

CM a
N M bN-M

N
where CM is the number of combinations of N things taken M at a

time. Hence, the p r o b a b i l i t y of e x a c t l y M occurrences of E1 i n a

sequence of N trials is

M
/C N = CNM pM qN-M
N N N-M
PM=CMab

where p = a/c and q = b/c are t h e p r o b a b i l i t i e s of El and Eo,

respectively.

-189-
 It i s meaningful, t h e r e f o r e , t o speak of t h e p r o b a b i l i t y t h a t ,

i n a sequence of N trials of t h i s experiment, t h e number M of

"successes" El s h a l l l i e between two given bounds; indeed,

There i s a fundamental p r o b a b i l i t y theorem ( 18) which states

that

where t = CF and
P9

2, Application t o t h e terminal. r a t i o s . The output of a non-

m u l t i p l i c a t i v e ( l9 ) Monte Carlo c a l c u l a t i o n without weights gives t h e

( "'J. V. Uspensky, Introduction t o Mathematical P r o b a b i l i t y , McGraw-

Hill Book Company, Inc., New Grk, 1937, p. 130.
( 19)I .e., not involving f i s s i o n , (n-2n) r e a c t i o n s , e t c . , i n which a
p a r t i c l e gives rise t o more than one p a r t i c l e .

-190-
number M out of N source particles which terminate in each of a set
of all-inclusive, mutually disjoint categories C. If we fix attention

upon some particular category C, we may regard the processing of a

source particle as an experiment which has as its outcome either the
event El of termination in this category o r Eo of non-termination

therein.

Now in any actual problem, there is a definite upper bound 1 on

the number of random numbers needed to process a particle due to the
existence of cutoffs based on energy, time, weight, number of collisions,

etc. We may then consider the class of all Y sequences of 1 random

numbers. A l l these sequences are equally likely, and a certain ~1 of

them determine a history terminating in category C. We may, therefore,

say that the probability of termination in C is p = a / y .

To be sure, this probability p of El is unknown; indeed, its
determination is precisely the object of the problem. We may, however,

gain some idea of the statistical reliability of the Monte Carlo result
by tentatively taking for p the value of M/N at some late stage of
the problem, when the latter ratio appears to have settled down, and
define l q= 1 - (M/N). Then the preceding theory states that the ratio

of the number of sequences of N trials resulting in a ratio of

I
~

satisfying the inequality

~

-191-
to the totality of all possible sequences of N trials is approximately

where t = 6 and erf(x) is the well-tabulated "error function"

Pq

for which we include a brief table.

X erf x

0 0
.2 .2227
.4 ,4284
-6 6039
.a .@21
1.0 . 8427
1.2 91°3
1.4 9523
1.6 9763
1.8 9891
2.0 09953
2.2 9981
2.4 9993
2 -6 09998
2.8 9999

Fig. 60

-192-
 By assigning to B the maximum error to be tolerated and to N
the number of source particles processed, one finds from the value of the
Integral the approximate chance of an error not exceeding e , which should

be close enough to unity for comfort and which approaches unity with in-

creasing N.
As an example, suppose a Monte Carlo run of 50,000 neutrons shows a

capture of 5000 neutrons in a certain zone, and the probability p of cap-

ture in this zone is desired with 5s accuracy. Let p = .l, e = .05p = .005.

Then t = 3.72678 and f(t) = erf(t/p) = erf(2,635) > ,9998. The

error I R I in using f(t) for P((M/N - pI < 6) does not exceed .0001.
This is, of course, far higher probability of safety than one can

usually hope for. Consider for contrast the extreme case of a capture p
Then
-4, and for N = 50,000,
= LO
of about .01 with a maximum error of 1$, B

t = ,225,and f(t) = erf(t/P) = erf(.175) E ,18, with IRI < .018. Here
it appears that a simple-minded Monte Carlo is quite ineffective, To be

sure, these requirements are far more stringent than are usually encoun-
tered, but such problems do arise, notably in counter design, and one sees
clearly the necessity for very large samples combined with ingenious de-
vices f o r improving statistics by use of weights in such cases.

3. The central limit theorem. The preceding remarks apply

strictly only to Monte Carlo procedures of the most straightforward
kind where no weights are employed, and a single, non-multiplying source

particle is followed until it ends its history in some terminal category.

-193-
It may then be considered as a physical particle undergoing an experiment
insofar as the random numbers employed may be regarded as truly random.

However, we have mentioned many devices employing weighted par-

ticles. In such cases, it is clear that the total weight M terminating
in a given category C after such procedures no longer represents the

number of successes out of N trials, that is, out of the processing of

N source particles, each initially of weight 1, and we can no longer

apply the theorem of §1. There is, however, a well-known generalization

called the central limit theorem, of which we shall give here a very

special case ,(20) which does apply t o the weight method,

Instead of the simple experiment which results in 1 with prob-
ability p and in 0 with probability q, let us consider an experiment
which can result in K different ways, to which we assign definite r e d
numbers W1, ..., WK, with probabilities pl, ..., pK, respectively,
pl+ ... + pK being unity. In a single trial of the experiment, therefore,
we may say that Wk has probability pk, k = 1, ..., K. In such a case,

we define the mean a = 2 pk Wk and t h e dispersion b = iT pk(Wk - 2

a)

= I2 p W2
k k
- a2 . Suppose now that N trials are made of this experiment,
and M is the sum of N weights so determined. Then our theorem states
that the probability

(
'
*
)
J
. V. Uspensky ( IC), page 294; cf. footnote (18).
where t = B and p N t 0 as N t 00, The estimation of PN is
given in terms of the third moment of the W-distribution.
Notice that the Bernoulli case is contained in this, since f o r

that case, we have p1 =: p, W1 = 1; p2 = q, We = 0 , a = P 1+ q 0 = P,

b = (p l2 + q
2
0 ) - a2 = p - p2 = p(1 - p) = pq, and the sum M of
the weights 1 and 0 recorded in the N trials is simply the number
of 1's (successes) in this sequence.

4. Application to problems using weights. Consider now a per-

fectly general Monte Carlo problem (21 ) in which weighted particles may

be used, each particle leaving the source with weight 1. For simplicity,

we congider the weights limited to a finite set of numbers, which is, in

fact, the case in machine computation. Each possible history assigns a

definite weight Wk to the particular terminal category C. Moreover,

it is clear that each such weight W, has a definite probability p,
of realization upon the trial of a random source particle, namelsthe
probability of a history which assigns W, to C, Thus the mean
a = pk Wk is the expected weight terminating in C pes source par-

ticle. If we process N particles in a Monte Carlo problem, tallying

the weight W(T) contribuwd to category C by each trial particle

T = 1, 2, ..., N, the final sum Z d T ) is the M of the theorem, and

M/N may be used as an approximation t o a. The difficulty in applying

' 21 'Multiplicative processes are not excluded, When they are involved,
we speak of the "expectation" a of termination in catebory C
rather than the probability whexher weights different from unity
are used or not.
t h e theorem l i e s i n t h e f a c t t h a t we do not o r d i n a r i l y have at hand an

estimate f o r t h e d i s p e r s i o n b. I n the simple Bernoulli case, we s a w

that b w a s expressible i n terms of p, but, i n general, a does not

N
determine b. Since b = C pk Wk
2
- a2, w e could use
7 =1
(WT)2/N -
as a guess f o r b i f we had provided i n t h e problem f o r t a l l y i n g t h e

(WT)2 as w e l l as t h e (W') f o r each trial p a r t i c l e 7 = 1, .,., N.

Note t h a t i n case of capture, one squares the t o t a l of capture

c o n t r i b u t i o n s from all c o l l i s i o n s of t h e p a r t i c l e ; i t i s not c o r r e c t t o

simply cumulate t h e squares of each weight captured.

I n a problem which does not involve m u l t i p l i c a t i o n ( f i s s i o n , (n-2n)

r e a c t i o n s , e t c . ) , t h e weight procedure i s b e t t e r than t h a t not employing

weights, a t a given s t a g e N, and f o r a given e r r o r c, t o the e x t e n t t h a t

t h e upper l i m i t of i n t e g r a t i o n exceeds LP,

P4
the probability i n

the l e f t member being identical, t h a t i s , the gain is based on the degree

t o which b is lesg than pq.

2
Now b = Cpk Wk - a2, pq = p - p2 , and
t h e expectation of category C i n such a problem i s i t s p r o b a b i l i t y ,

i.e., a = p. Hence, f i n a l l y t h e improvement rests on t h e f a c t that

zpk WES p = a = "

k k'

which i s c l e a r s i n c e Wk $ 1 f o r a11 k. However, the degree of

improvement can be c a l c u l a t e d only i n t h e most t r i v i a l examples. It

e<<
i s a t least c l e a r why t h e improvement i s so great i n cases of small

q u a n t i t i e s , since, then, t h e Wk are q u i t e small and Wk.

 It may a l s o be remarked t h a t t h e theorem of s e c t i o n ( 3 ) a p p l i e s

p e r f e c t l y w e l l i n finding t h e expected number a (per source p a r t i c l e )

of p a r t i c l e s terminating i n a p a r t i c u l a r category C, while t h e theorem

of 5 1 does not. The unfortunate f e a t u r e involved i n the a p p l i c a t i o n

of t h e former theorem i s t h a t one has no simple way of knowing the

d i s p e r s i o n b. Problems of t h e complexity i n p r a c t i c e usually f o r b i d

complete t a b u l a t i o n s of W"'s f o r purely s t a t i s t i c a l s t u d i e s Actually

one provides the b e s t weight t r i c k s one can t h i n k of and t r i e s t o judge

t h e r e l i a b i l i t y of an output r a t i o by i t s convergence and s t a b i l i t y as

N increases.

5. I l l u s t r a t i v e examples. We include some extremely t r i v i a l

examples which nevertheless may h e l p t o f i x t h e ideas involved i n the

preceding sections.

(1) Consider the non-multiplicative problem defined by t h e

-197-
 We may limit random numbers to be either 0 or 1. The upper

bound on lengths of sequences of random numbers required is 2, and

the total number Y of such sequences is 2* = 4. In fact, we have the

following results:

00 implies termination in A

01 implies termination in A

10 implies termination in B
11 implies termination in C

If we fix attention on C, we find that the number a of sequences

favorable to C is 1, and the probability of termination in C is
a/Y = 1/4. Hence, relative to this category, we have p = 1/4 for
W = 1 and q = 3/4 for W = 0, with dispersion pq = 3/16.
Suppose for comparison that we decide to use a weight method for

this problem as indicated below.

7iZ~~+-.*
- W + B - B
2
1 I
r--
2

-w+c-c

Q‘@
Now sequences of random numbers are limited to length Y = 1, and

we have the results

-198-
 0 implies A = 1/2 B = 1/2

1 implies A = 1/2 C = 1/2

Again f,ixingattention on C,

P1 = 1-12
for w1 = O

and

a = p W + p W
1 1 2 2
=1/4

b = p1 :W + p2 WE - a2 = 1/16 < pq = 3/16

which illustrates in this very trivial example the type of improvement

discussed in the preceding section.

(2) Now suppose we have the simple multiplicative system in-

dicated by the following scheme:

-199-
 Thus a source particle may terminate in A with probability 1/2

or,in the contrary case, it may either terminate in B with probability

-
1/2 or produce two particles with probability L/2. The two progeny

thereupon have equal chances of terminating in A or B.

Here the maximum length 1 of random number sequences of 0's and
4
1's is 4, Y = 2 = 16, and we have the following correspondence between
random number sequences and terminal weights:

0000 A = l B = O
0001 A = l B = O
0010 A = l B = O
----
0011 A = l B = O
0100 A = l B = O
0101 A = l B = O
0110 A = l B = O
----
0111 A = l B = O
1000 A = O B = l
1001 A = O B = l
1010 A = O B = l

----
1011 A = O B = l
1100 A = 2 B = O
1101 A = l B = l
1110 A = l B = l
1111 A = O B = 2

Let us consider only category B. We have by inspection the weights

and corresponding probabilities

-200-
 PO
= 9/16 wo = o

p1
= 6/16 w1 = i

p2 = 1/16 wg = 2.

The expectation and dispersion are

a = ~p w
k k
=1/2

2 2
b = X p W - 8 =3/8
k k

Let us compare this procedure with one which, in the event of

multiplication, processes one particle Of weight 2 rather thantwo

particles each of unit weight, as indicated below.

W t A - A

W t B - B

-201-
 We now obtain the results

000 A = l B = O
001 A = l B = O
010 A - 1 B = O
---
011 A - 1 B - 0
100 A - 0 B - 1
101 A = 0 B-1
110 A = 2 B = O
111 A - 0 B-2

and f o r category B,

Po = 5/8 wo = 0

= 210 w1 = 1

p2 = 110 w2 = 2

a = zpi wi = 112

b = C pk Wg - a* = 1/2

It will be noted that the revision has resulted In Increasing the

dispersion. That this is not an inevitable consequence of employing

weighted multiplication, the reader may verify by observing that the

dispersion for category B becomes 0 in the following treatment of

the same problem:

-202-
-203-
 APPENDIX

SUMMARY OF CERTAIN PROBLEMS RUN ON MANIAC I

Introduction. We conclude thi’s r e p o r t with a very b r i e f summary

of some of t h e Monte Carlo t y p e problems which w e have set up and run

on t h e MANIAC I a t Los Alamos during t h e p a s t t h r e e years. The coding

f o r t h e s e problems has been done by some extremely rapid and resourceful

coders a t Group T-7 of t h i s laboratory, whom we wish t o thank c o l l e c -

t i v e l y , and t o whom we refer i n d i v i d u a l l y i n t h i s appendix and LA-2121.

It i s hoped t h a t t h e reader unfamiliar with t h e method w i l l f i n d

in t h i s b r i e f l i s t of problems some i n d i c a t i o n of t h e a c t u a l s i t u a t i o n s

i n which t h e methods of t h e t e x t f i n d a p p l i c a t i o n . It may a l s o perhaps

h e l p those contemplating use of t h e method i n problems of t h e i r own t o

form some judgment of i t s p o s s i b i l i t i e s .

The c o l l e c t i o n included here n a t u r a l l y omits t h e many c l a s s i f i e d

problems which have been completed.

-204-
 Problem 1

COMPTON COLLISIONS IN A SPHERICAL SHELL

-
Coder: R. G. Schrandt

Requester: G. H. Best and H. C. Hoyt

Medium: Free-electron "gas" in spherical shell of radii Ro < R1.

Source: Monoenergetic photons directed radially outward from R
0
.
Various initial energies were studied.

I
Photon parameters: x, y, z, u, v, w, E, g.

Physical processes: Compton collisions of photons with free electrons.

Output: Correlated energy and angle distribution of photons escaping R1,

together with totals lost to core and lost to energy cutoff. Any photon
impinging on the inner boundary Ro was considered absorbed.
I
Remarks: "Forced first collision" device was used. The angular distsi-
bution desired was that of Ch. V I I I , Fig. 53.

Problem 2

COMPTON COLLISIONS IN A SOLID SPHERE

Coder: R. G. Schrandt
__li

Requester: G. H. Best and H. C. Hoyt

-205 -
Medium: Free e l e c t r o n gas i n s o l i d sphere of r a d i u s R1.

Source, photon parameters, -

and physical processes: As i n Problem 1.

Output: A s i n Problem 1 with exception of core losses.

Remarks: I t w a s d e s i r e d t o study the e f f e c t of t h e core on t h e escape

d i s t r i b u t i o n by comparison of t h e r e s u l t s of Problems 1 and 2. Actually,

by s u i t a b l e devices i n t h e f l o w diagram, it w a s po'ssible t o combine the

two problems i n t o a s i n g l e code, keeping two sets of terminal category

r e g i s t e r s and following a photon path f u r t h e r f a r purposes of Problem 2,

after i t had h i t the core and been l o s t f o r Problem 1. This i s a device

which c u t s machine time almost i n h a l f .

Problem 3

NEUTRONS I N A SPHERICAL SHELL

Coder: J. M. Kister

Requester: J. R. Beyster

Medium: Hollow s p h e r i c a l s h e l l of heavy nuclei.

Source: Monoenergetic i s o t r o p i c p o i n t source a t c e n t e r of sphere.

Neutron parameters: R , w, v .

Physical processes: I n e l a s t i c s c a t t e r i n g , regarded as terminal, e l a s t i c

s c a t t e r i n g according t o given l a b . system d i f f e r e n t i a l cross section,

attended by no energy loss.

-206-
Output: Number of neutrons lost to inelastic collision, number of
neutrons escaping after u collisions, V = 0, 1, ..., 10, and u 2 11,
Remarks: The distribution function P. was tabulated for scattering at
J
a cosine a => a3 f o r 32 strategically chosen intervals of the cosine

range. The problem was run as a check on an analytic method of'

H. A. Bethe based on transport cross sections, the results of the
latter being computed by J. R. Beyster. The Monte Carlo procedure
was coded in two different ways, one using the weight method f o r in-
elastic termination, the other not, The results were in excellent
agreement and checked Bethe's result. Transmission was obtained

without forced first collision and agreed closely with the analytic

value. The problem is written up in LA-1583, ----

A Monte Carlo Deter-

-- -- --
mination of the Escape Fraction for a Scattering Spherical Shell with

-
C e n t r a l Point Source,

Problem 4

ENERGY INDEPENDENT SCATTERING IN A CYLINDER

Coder: R. L. Bivins
7

Requester: M. Walt

Medium: Solid cylinder of homogeneous material. Problem run f o r many

heavy elements.

-207-
Source: Monoenergetic neutrons in parallel beam incident on lateral
surface of cylinder in direction normal to axis.

Neutron parameters: x, y, z, u, v, w, W, v .

Physical processes: Inelastic collision, treated as terminal, elastic

scattering at source energy governed by differential cross section. No

energy loss assumed on elastic scattering,

Output: Transmission, loss to inelastic collision, number of emergent

scattered neutrons not hitting detector band, classification of neutrons

hitting band zone i with u = 1, 2, 3, and v 2 4 collisions,' Detector
band was coaxial with cylinder and calibrated in 5' zones. Geometry
w a s khat of Fig. 55.

Remarks: Weights were used for transmission and inelastic losses,

Elastic scattering probabilities P were used as indicated in Problem 3.
J
It will be noted that many of the problems included are of this general
nature, the Monte Carlo results being used in correcting experimentally
determined differential cross sections and cross sections of other

processes for multiple collision effects in thick targets.

Problem 5

ENERGY DEPENDENT SCATTERING IN A CYLINDER

Coder: R. L. Bivins
Requester: M, Walt

-208-
Medium: Solid cylinder of homogeneous material. Many light elements
were run.

Source: Same as Problem 4.

Neutron parameters: x, y, z , u, V, w, E, g, W, V.

Physical processes: Inelastic collision regarded as terminal, elastic

scattering governed by differential cross sections ( P ) in C.M. system.

Output: Transmission, loss t o inelastic collision, Oloss

g I
to energy cutoff,

emergent neutrons failing to hit detector band (same geometry as in

Problem 4 and Fig, 55), and those hitting band in zone i with IJ =
> 1, 2,

and u 2 3 collisions.

Remarks: Weights used for transmission (forced first collision device)

and for inelastic losses. Differential elastic scattering handled by

the von Neumann device using
*( g , p )
(I = A + B Cr + C P 2 as indicated
g g g
in Ch, V, S5. Tables stored for free path and f o r coefficients
65
C f o r suitable energy groups g.
Bg’ 63

Problem 6

CYLINDRICAL SHELL WITH PARALLEL SOURCE

-
Coder: R , G. Schrandt

Requester: J. D, Seagrave
Medium: Cylindrical shell of homogeneous material. Problems run for
CD2, C, and CH2 shells.

-209-
Source: Monoenergetic parallel beam of neutrons incident on lateral
surface of cylinder in direction normal to axis. Various source energies

were studied.

Neutron parameters: x, y, z, u, v, w, E, g, W, V .

Physical processes: Elastic scattering on C, D, H according to C.M.

differential cross sections oe(p). Such scattering on H was assumed
g
*
isotropic; C scattering was treated by a fit u (g,fl) = A + C N 2 ;
@ ; &
while for D, the double interpolation method of Ch. V, 05 (Fig. 37) was
used on the basis of a stored table
*

Output: Transmission, loss of "degraded" neutrons to energy cutoff,

emergent undegraded scattered neutrons failing to hit band, and clas-
sification of scattered undegraded neutrons hitting the band (geometry
of Fig. 5 5 ) into the following categories:

(a) hits on band zone i with energy above given ,upperbound ri

(b) hits on band zone i with energy below given lower bound Li

(c) hits on band zone i with energy between ei and Ti with V = 1

(d) hits on band zone i with energy between Li and Ti with V > 1.

Problem 7

14 MEW NEUTRONS IN A CYLINDRICAL SHELL

-
Coder: R. G. Schrandt
Requester: L. Rosen and L. Stewart

-210-
-
Medium: Solid cylinder of homogeneous scattering material. Problems run

for two heavy elements (Bi, Ta) .

Source: Parallel beam of monoenergetic (14MeV) neutrons incident on

-
base of cylinder. Geometry of Fig. 52,

Neutron parameters: x, y, z, u, v, w, Y, E, f (group index for elastic

scattering differential cross section purposes), g (group index for free

,
path and probability of elastic scattering constants) h (group index

for escape classification) .

Physical processes: At most, two collisions were permitted, escape being
forced after a second collision, with the direction obtaining as a result
of that collision, The branching process of Fig, 61 was involved, The
assumptions made for the types of collisions were as follows:
(a) Elastic collision: no energy loss, lab. differential cross
sections of(a) f o r five energy ranges f.

(b) d-process. Energy of emergent neutron equidistributed between

5 and 12 MeV. Process occurs only at 14 MeV. One differential
cross section given f o r lab. angle a.

(c) (n-2n) process. Two neutrons emerge with energy distributions

nl -- E ' exp(- E'/T~)and n2 = E' ea(- E~/T*)on o I: E' 6 5 .

Isotropic emergence in lab. system. Process occurs only at

14 MeV. (Cf. Ch. V, 815#for method.)

 (d) Inelastic collision. One neutron emerges, isotropically in

lab. system, with energy distribution n = E' em(- E/T) where

T = k fw). Process occurs only below 12 MeV. (Cf. Ch. V,

§ll,for method ) .

E = 14 Mev
r

Fig, 61

Output: Energy and angle distribution of escape.

Remarks: Symmetry of system about source distribution obtains, and

counter was at infinite distance so that the optimum conditions of

Fig. 52 obtained. (Cf. discussion of Ch. V I I I , 83.)

-212-
 Problem 8

CYLINDRICAL SHELL W I T H POINT SOURCE

-
Coder: R. L. Bivins

Requester: L. Cranberg

Medium: C y l i n d r i c a l s h e l l of homogeneous material. Various heavy

elements run.

Source: Monoenergetic point source i n c i d e n t on l a t e r a l surface of cyl-

i n d e r from an e x t e r n a l p o i n t midway between base planes, ( C f . Ch. 11,

85d, and Fig, 11.)

Neutron parameters: x, y, z , u, v, w, E ( t h r e e d i s c r e t e values El, E*,

E$.

P h y s i c a l processes: A t y p i c a l element ( B i ) had two e x c i t e d state energy

levels at A = .9 and A 2 = 1.65 Mev above t h e normal state. A neutron

with energy E c o l l i d i n g with the nucleus can thus s c a t t e r e l a s t i c a l l y ,

with no energy loss, or i n e l a s t i c a l l y , emerging w i t h energy E - A l if

E => A 1 or with energy E - A if E 2 A2, Since t h e source energy was

2.5 MeV, t h e h i s t o r i e s of Fig. 62 were possible. Inelastic scattering

was assumed i s o t r o p i c i n t h e lab. system, while e l a s t i c s c a t t e r i n g was

assumed t o obey a given d i f f e r e n t i a l cross s e c t i o n a ( a ) applicable a t

a11 energies involved.

-213-
 E = 2.5 E' = 1.6

Fig, 62

Output: At most, two collisions were allowed. Neutrons having a third

collision were classified in a terminal category L At most, one in-
3'
elastic collision was permitted. Neutrons suffering a second inelastic
collision were classified in a terminal category 12. Thus a11 escaping

neutrons emerge with one of three possible energies El, E2, E3. Of
these, the neutrons hitting the detector band (geometry of Fig, 55) in
angular zone i and energy E were tabulated as Ni(E) for each E.
Also recorded were the transmission, escapes not hitting the band,
- -
total hitting band after one elastic collision, and the total hitting
-
t h e band after one collision of any kind.

-214-
 Problem 9

A SCINTILLATION COUNTER

Coder: R. L. Bivins
_.

Requester: R, F. Taschek and N. J, Terrell

Medium: Cylindrical s h e l l containing various compounds of H, C, 0, and Cd.

Source: Monoenergetic point source at center of cylindrical hole (vacuum)

with various deterministic directions.

Neutronjparameters: x, y, z, u, v, w, E, g .

Physical processes: H-capture, Cd-capture, elastic scattering isotropic

in C.M. system for H, 0, Cd, and according to a differential cross section
2
= A -t C P for C.
g g

Output: Loss from s h e l l surface, H-capture, Cd-capture, and losses to

energy cutoff classified according to 6 radiaL and 12 height zones.

Remarks: The purpose of the problem was the determination of the space
distribation of degraded neutrons referred to, for use as input in a
subsequent age theory calculation (not performed by us) connected with

sensitivity of a scintillation counter design.

 Problem 10

A "LONG COUNTER"

-
Coder: R . L. Bivins

Requester: A. W , Schardt

Medium: Cd cylinder (medium m = I, parts p = 1, 2, 4 of Fig. 63) con-

t a i n i n g B (gaseous boron compound under pressure) and vacuum space, and

surrounded by a hydrocarbon (medium m = 2, parts p = 3,5) i n cylindricaZ

outer container as shown.

m - 2 I p = 3

m - 2
p = 5

I
Im = 1 , p = 11

P P'

Fig, 63

-216-
Source: Monoenergetic neutrons in parallel beam impinging on various

radial zones ( p , p ' ) of base of outer container. Various source energies

studied

Neutron parameters: x, y, z, u, v, w, E, g, W, m, p.

Physical processes: E-capture, Cd-capture, B-capture, elastic scattering

on H, Cd isotropic in C.M. system, C as in Problem 9.

Output: Escape from Cd-base ( p = 1, p = 4), escape from hydrocarbon

surface (p = 3, 5), H-capture, Cd-capture, B-capture, losses to energy

and weight cutoff,

Remarks: The purpose was to estimate the efficiency of B-capture as a

function of source energy and distance from axis in connection with
design of a long counter. Weights were used throughout, neutrons
passing through B cylinder were attenuated exponentially, and the
multiplication device of Ch. V, 616,was used. T h i s problem probably
represents the most demanding one statistically that we have set up.
Enormous samples were involved, and the actual running of the machine

was done by Schardt and co-workers.

Problem 11

A PROBLEM CONNECTED WITH NEUTRINO DETECTION

-
Coder : R. L. Bivins

Requeste,r: F. Reines and C, L, Cowan

Medium: Plane s l a b of homogeneous CH
1.5
on t h e 2 - i n t e r v a l - H <= z 5 H,
with Cd-layer a t z = 0. Two d i f f e r e n t thicknesses of Cd were studied.

The s l a b thickness 2H was e s s e n t i a l l y i n f i n i t e .

Source: Monoenergetic neutrons a t various heights z and with various

d i r e c t i o n s from OZ. A series of source energies was studied.

Neutron pamameters: z, w, E, 7, W, and a parameter i n d i c a t i n g the status

of the neutron r e l a t i v e t o two problems run simultaneously. (See Remarks,)

Physical processes: H-capture a t thermal energy (only), Cd-capture,

e l a s t i c s c a t t e r i n g on H i s o t r o p i c i n C.M. system a t energies above

thermal, and i s o t r o p i c i n lab. system a t thermal with no f u r t h e r energy

loss. Elastic s c a t t e r i n g on C w a s assumed i s o t r o p i c i n t h e labe system

with no energy l o s s a t a11 energies.

Output: C l a s s i f i c a t i o n of Cd-captures i n .2 psec i n t e r v a l s from t i m e

( t = 0 ) of leaving source t o time cutoff of 30 psec, 150 intervals i n

a l l , hydrogen capture, l o s s t o t i m e c u t o f f .

Remarks: The following conventions on Cd w a l l capture were adopted.

The t h i c k e r wall captures all neutrons impinging on it a t energies

< .4 ev, while t h e t h i n n e r w a l l captures -

all neutrons impinging a t

energies <.34. The t h i c k and t h i n Cd-wall problems were run simulta-

neously by keeping separate counters N1, ...? N

150
and N i , ..., Ni5,-
f o r t h e two t i m e d i s t r i b u t i o n s a n d following the path of a neutron

f u r t h e r for purposes of the t h i n w a l l problem after it had been c l a s -

s i f i e d terminally as captured i n the t h i c k w a l l problem.

-218-
 Problem 1 2

AN Iy DETERMINATION

Coder: J. M. Kister
Requester: R . B. Lazarus

Medium: Bare homogeneous Oy sphere.

Source: I n i t i a l l y an a r b i t r a r y d i s t r i b u t i o n of neutrons i n energy

groups g = 1, 2, 3, radial zones z = 1, ..,, 10, and cosine i n t e r v a l s

j = 1, ...) 10, 300 c a t e g o r i e s i n a l l .
Neutron parameters: R, w, V = VE, g.

Physical processes: The usual processes f o r a fissionable element were

treated by the t r a n s f e r matrix method of Ch. V, g12.

Output: The problem w a s run i n cycles of f i x e d time length AT, the

output N of one cycle serving as t h e number of neutrons t o be

g,z,J
processed i n category ( g , ~ j, ) during t h e next. The problem run was

Close t o c r i t i c a l , so t h a t no re-normalization was involved. (cf* Ch. 11,

68,and ~ h ,SV, 59.)

-219-
 Problem 1 3

NEUTRON FLUX I N A I R

-
Coder: R. G. Schrandt
Requester: J. Hall, R . G. Wagner, and G. M. Wing
-
Medium: Sphere of homogeneous air zoned for computational purposes by

radii R1< R2 < R3.

Source: Monoenergetic isotropic point source at center R = 0. Various

initial energies were studied.

Neutron parameters: R, w, E, g, W, Z.

Physical processes: Elastic scattering on 0 and N assumed isotropic in

C.M. system.

Output: Loss to energy cutoff in each zone, loss from system (at R
3
>,
loss to weight cutoff. These were the only terminal categories. How-
ever, the purpose of the problem was to determine total flux at R1, R2,

R Counters were, therefore, provided t o tallr all neutrons crossing

3'
each boundary R
I in each energy group g in outward direction and in
inward direction. (Cf. Ch. IV, S4.)

Remarks: Weights were used in connection with the device discussed in

Ch. V, %7,for avoiding loss of trajectories to energy cutoff. Thus no

trajectory terminated until loss to weight cutoff or loss from R

3'

-220-
 Problem 14

ESCAPE DISTRIBUTION FROM A CARBON SLAB

Coders: M. B, Wells and R, G, Schranat

Requester: T. B. Taylor

Medium: Plane slab of C on Q <= z <= Z,

Source: Monoenergetic neutrons upwardly directed at z = 0 in cosine
___c

distribution.

Neutron parameters: z , w, E, g".

Physical processes: Elastic scattering on C governed by a differential.

cross section of form 0

96
(g,p) = A +C
@ ; G
c1
2
.
Output: 13 energy group classification escapes at z = 0 and z = Z.

Energy cutoff losses classified according to position in five sub-

intervals,

Problem 1 5

ESCAPE DISTRIBUTION FROM A SPHERICAL SHELL

Coders: M. B. Wells and R. G. Schrandt

.
_Lc

Requester: C L Longmire
Medium: Spherical shell of hydrocarbon on R1
_IC <= R 5 R2.

-221-
Source: Monoenergetic neutrons d i r e c t e d i n t o s h e l l a t R1 i n cosine

distribution.

Neutron parameters: R , w, E, g ” .

Physical processes: E l a s t i c s c a t t e r i n g on H i s o t r o p i c i n C.M. >ystem,

on carbon as i n Problem 1 4 ,

Output: 1 3 energy group c l a s s i f i c a t i o n of escapes a t R1, escapes a t

R2, l o s s t o energy c u t o f f .

Problem 16

A ROCKET MOTOR

Coder:
I__
R , L. Bivins

Requesters: G. I. B e l l and C. L. Longmire

Medium: Cylindrical container of C y U, H, with C of uniform d e n s i t y ,

H and U of prescribed d e n s i t i e s i n each of a set of s h e l l zones, de-

f i n e d by 7 r a d i a l and 15 height i n t e r v a l s . This cylinder was sur-

mounted by a coaxial cyclinder of H and Be a t s p e c i f i e d d e n s i t i e s ,

and t h e s e two c y l i n d e r s were surrounded by a c y l i n d r i c a l s h e l l of Be.

The U-density d i s t r i b u t i o n w a s s u b j e c t t o change from problem t o

problem, of which a whole series was run. ( C f . Fig. 64.)

-222 -
I
 Z

Be -'H
p = 12
I

I Be

' i
p = 3

C j
I
p = 1

i
I

Fig, 64

Source: I s o t r o p i c source i n f i s s i o n energy spectrum, s p a t i a l l y d i s t r i b -

uted as prescribed among the lo5 zones. The s p a t i a l d i s t r i b u t i o n also

v a r i e d with t h e problem and was so normalized t h a t i n t e g r a l numbers

of neutrons were involved i n each source.

-223-
Neutron parameters: x, y, z, u, v, w, E, g,y, W, and a sector parameter
(See Fig. 64.) The zone parameter
p = 1, 2, 3 .
7 is really a pair i,j
which define the associated radial and height zones.

Physical processes: Be-capture at "thermal" energy; U-capture, fission,

elastic and inelastic scattering; elastic scattering on C and Be assumed
isotropic in lab. system, using transport cross sections for free path
contributions but with correct energy losses. Elastic scattering on H
isotropic in C.M. system,

Output: Losses from system classified according to surface of escape;

Be-capture; U-capture; number N' of fissions in each of the lo5 spatial
7
zones. Fission was regarded as a terminal event.

Remarks: There was no loss to energy cutoff; neutrons degrading to

"thermal" energy were allowed to scatter isotropically in the lab.
system without further loss of energy. A neutron was followed until
capture, fission, or escape occurred. The purpose of the problem was

to determine a U-density distribution together with a fission source

distribution in space which would produce a steady state in the sense
that the output N' and input N should satisfy the equations
7 7
-1
'a' = u o
" Nr
where yo is the fission multiplication constant. Such problems arise

in the design of nuclear rocket motors, the hydrogen of the problem

being present in the form of a propellant gas.

-224-
 Problem 17

PHOTON DIFFUSION IN A REACTOR

-
Coder: R. G. Schrandt
Requesters: M, E. Battat and B. M. Carmichael

Medium: Cylindrical assembly consisting of inner fuel region, tamper,

and top and bottom reflectors of given electron densities.

Source: Prompt Y -rays in exponential energy distribution n(E)dE

= k exp(- 1.01 E)m, .2 5 E <= 5 MeV, isotropic in direction, and with

given radial and heightwise photon density distributions throughout

the fuel.

Photon parameters: x, y, z, u, v, w, E, Q, m (medium index indicating

region of occupancy: fuel, tamper, or reflectors), R " ~ , H ' , HI'
R I ~ ,

defining radial and height boundaries of region m.

Physical processes: Compton scattering of photons on electrons with

energy losses due to such scattering and to cutoffs (varying with m)
chosen at the point where the photoelectric effect becomes an

effective absorption,

Output: Distribution in energy and angle of photons escaping each of

outer surfaces, energy deposition within the system in 82 spatial
regions.

-225-
 Problem 18

A T H I C K Z r TARGET PROBLEM

Coder: R. G. Schrandt

Requesters: L. P. Stewart and L. P. Rossn

Medium: S o l i d Z r cylinder.

Source: P a r a l l e l beam of monoenergetic (14 MeV) neutrons i n c i d e n t on

base of cylinder.

Physical processes: Those of Problem 7 except t h a t t h e &2n) process

was replaced by a more general process involving production cross

s e c t i o n s f o r neutrons i n each of two energy d i s t r i b u t i o n s with a

t o t a l expectation of 2.32 neutrons.

Output: A s i n Problem 7.

Remarks: The problem was handled by a r e v i s i o n of the code of t h e

problem c i t e d .

Problem 19

A THICK C TARGET PROBLEM

-
Coder: R. G. Schrandt

Requesters: L. P. Stewart and L. P. Rosen

-226-
Medium: S o l i d carbon cylinder.

Source: Monoenergetic (14.1 MeV) neutrons i n parallel beam incident on

base of cylinder.

Neutron parameters: u, v, w, x, y, z , E, g, V .

Physical processes: E l a s t i c s c a t t e r i n g according t o given d i f f e r e n t i a l

c r o s s s e c t i o n s curves on f o u r ranges between 1 4 and 0 Mev; capture;

i n e l a s t i c s c a t t e r i n g f o r l e v e l s 4.4, 7.6, 9.6, and 11.2 MeV (the latter

an average f o r l e v e l s > 9.6), i n e l a s t i c s c a t t e r i n g a t 4.4 l e v e l f o r

i n c i d e n t energy >10 MeV governed by given u (E,€)), a l l o t h e r i n e l a s t i c

s c a t t e r i n g assumed i s o t r o p i c i n C,M. system.

Output: D i s t r i b u t i o n i n 28 energy groups and 1

0
' angular i n t e r v a l s

(with v e r t i c a l ) of neutrons escaping after one o r two c o l l i s i o n s ,

t o t a l number of second c o l l i s i o n s ,

Remarks: F i r s t c o l l i s i o n forced, neutrons forced out a f t e r a second

collision, Method of Ch. V, 69,used t o determine parameters after

inelastic collision.

Problem 20

HEAVY WATER EXPEXIMENT

Coder: R , G. Schrandt

Requester: J. N. Grundl

-227-
Medium: Heavy water in (A) spherical shell, (B) hemispherical shell.

Both problems run simultaneously. A series of runs involved various

inner and outer radii.

Source: Isotropic point source at center of system in fission-neutron

energy distribution, .1 5 E 5 10 MeV.

Neutron parameters: u, v, w, x, y, 2, E, g, I( (parameter indicating

traJectory lost relative to problem (B) or not).

Physical processes: Elastic scattering on D and 0 assumed isotropic in

C.M. system; transport cross section based on this assumption used f o r

free path.

Output: For problems (A) and (B), loss to outer surfaces, loss to energy

cutoff, classification in 25 energy groups of neutrons re-entering central

cavity with normal angles in each of four angular ranges.

-228-