Monte Carlo calculation of multiple scattering effects in thermal neutron scattering experiments

COMPUTER PHYSICS COMMUNICATIONS 7 (1974) 289-317. NORTH-HOLLAND PUBLISHING COMPANY

MONTE CARLO CALCULATiON OF MULTIPLE SCATTERING EFFECITS IN THERMAL NEUTRON SCATTERING EXPERIMENTS* J.R.D. COPLEY** Argonne National Laboratory, Argonne, Illinois 60439, USA Received 5 July 1973

PROGRAM SUMMARY Title of program: SLOW NEUTRON MULTIPLE SCATTERING Catalogue number: ACIC Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. heland (see application form in this issue) Computer: IBM 370/195; Installation: Argonne National Laboratory

a smoother function of scattering angle and energyand thus can be calculated using an approximate scattering function. Method of solution The program is adapted from that of Bischoff [1]. Given a scattering function the program tracks successive neutrons within the sample and/or sample container in a Monte Carlo

Operating system: OS/360, Release 21.7

fashion. At each scattering point, the response, within each time channel and for each detector, is calculated. A cut-off is imposed to avoid tracking neutrons indefinitely. Single and multiple scattering contributions for each time channel

Programming language used: FORTRAN IV

and for each angle, are separately accumulated.

High speed storage required: 75000 words. No. of bits in a word: 32 Overlay structure: None

Typical running time Once the scattering function has been set up, the Monte Carlo

No. of magnetic tapes required: None

loop takes 0.2 to 0.5 milliseconds per collision, per angle, per time channel. Typical running times for a complete problem are of order 1—10 mm on the IBM 370/195. If elastic coherent

Other peripherals used: Card reader, line printer, magnetic disk (up to three files)

scattering is to be calculated, the program may take considerably longer to obtain adequate statistics.

No. of cards in combined program and test deck: 5576 Card punching code: EBCDIC Keywords: Solid state, multiple scattering, thermal neutron scattering, elastic coherent scattering, elastic incoherent scattering, inelastic scattering, Monte Carlo method, alpha—beta sampling scheme. Nature of physical problem In a thermal neutron scattering experiment, the measured cross section includes both single and multiple scattering events. The former can be obtained from the measurements by subtracting an estimate for the latter, which is generally *

**

Based on work performed under the auspices of the U.S. Atomic Energy Commission. Present address: Institut Max von Laue-Paul Langevin, B.P. 156, 38042 Grenoble C~dex,France

Unusual features of the program The minimum high speed storage requirement is about 49000 words. If inelastic scattering is included, further high speed storage is required, the amount depending on how accurately the inelastic scattering functions are defined. In the present case an extra 26000 words are used to store these functions. Restrictions on the complexity of the problem The present program is limited to one geometry and to no more than eight scattering angles. The sample (and sample container) must be isotropic materials. References [11 F.G. Bischoff, stitute (1970); Ph.D. Thesis, Rensselaer Polytechnic InF.G. Bischoff, M.L. Yeater and W.E. Moore, Nucl. Set. Eng. 48 (1972) 266.

290

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

LONG WRITE-UP I. Introduction The object of a thermal neutron scattering experiment is to measure the double differential scattering cross section, per unit solid angle &1, per unit scattered energy Ef, d2 U ‘k

-÷

k

\ 1.

d~dEf‘ 0 f This quantity is a measure of the probability for scattering from initial state k state kf (energy Ef) (fig. 1). 0 (energy E0) to final A typical experimental setup includes components which (i) monochromate the incident neutron beam and define its direction, and (ii) select scattered neutrons having a particular energy and final direction. In a triple-axis spectrometer, the monochromating and analyzing functions are performed using Bragg reflections from single crystals. On the other hand, some spectrometers utilize mechanical choppers and time-of-flight measurement techniques to achieve these functions. A reasonably comprehensive review of neutron spectrometry is contained in ref. [1] In many cases no analysis of scattered energy is performed, and one attempts instead to measure a single differential cross section .

da

d2cr

=

2cr 1k ~ k d d~ldEf” 0 f~S~ In an actual neutron scattering experiment, the measured cross section is reduced with respect to the single scattering cross section because of self shielding by the sample. On the other hand, it is augmented by multiple scattering within the sample. Furthermore, the sample is very often surrounded bya container which contributes single scattering as well as multiple scattering, and which further attenuates the beam. A system of programs which is routinely used to make corrections for these effects is described elsewhere [2] The present paper describes the program MSCAT, which is used to calculate multiple scattering contributions to the measured scattering cross section. The present program, MSCAT, is a modified version of the code MSC written by Bischoff [3,4] The .

f~~i~’ (k

0 -÷ kf) dEe,

P

where the path P of integration over final energy is performed in some well-defined manner. For example the path P, in the case of a conventional two-axis diffractometer is such that the scattering angle p (fig. 1) is fixed. In this work we restrict our attention to the case of an isotropic target, such as a liquid, an amorphous solid, or a polycrystalline material. In this case only the magnitudes of the wave vectors k 0 andWe kf, and the scattering angle p, need be considered. therefore write the cross section as

major improvement over Bischoff’s code is the provision for a container surrounding the sample. Container scattering contributes significantly to the meas-

TARGET ~.CONTAINER SAMPLE

E — —

4’ Ef kf Fig. 1. A typical scattering event. The neutron’s initial and final energies and wave-vectors are E0, It0 and Ef, kf respectively. The target consists of a sample enclosed in a container.

ured cross section in experiments where, for example, the container is made of a refractory metal (high-

temperature measurements), or where the container is necessarily thick (high-pressure measurements).

Other differences between MSCAT and MSC include (i) a more general treatment of the elastic coherent scattering into a detector, (ii) the ability to run a simulation with perfect resolution, and (iii) the provision for any combination of elastic coherent scattering, elastic incoherent scattering, and inelastic scattering, in sample and/or container. The user of MSCAT is strongly advised to consult ref. [3] ,both for a better understanding of various aspects of this program, and for a very useful introduction to the techniques of Monte Carlo simulation.

291

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

Bischoff’s code, MSC, allows the user to include in his simulation some (not all) of the uncertainties which contribute to instrumental resolution. These provisions remain in the present program, but in practice we normally make a separate correction for resolution as described in ref. [2] In some cases the dependence of the resolution function on energy transfer, which is neglected in ref. [2] ,must be included in the resolution correction. This subject is outside the scope of the present paper. Section 2 contains definitions of various quantities which are used in the program. The general structure of the program is described in Section 3, and the individual subprograms are introduced in Section 4. Input and output are explained in Section 5, and Section 6 contains a description of the COMMON blocks. The test run and other tests of the code are described in Section 7. Several aspects of the program are briefly discussed in the final section. .

2. Definitions

/

~‘i

1.5

~O5

/L~0

~L~-05

fL.-I

0.5 OC

1.25

~

-1.0

Fig. angle loci in (a,a) space, for e~= 1.25. The lines 2.~zConstant = 1 and ~i = —1 define the kinematically allowed region.

a~ 2~_i+ 1~ 2e~21(e1 +13k)’ 12pp or a = 2e 2(e 2p, 0 + 13 2e~/ 0+13)” which defines the allowed values of a and j3 for a particular initial energy ~0• This is illustrated in fig. 2, where constant angle loci are shown, for a particular choice of ~. The kinematically allowed region is de—

2.1. A scattering event The neutron’s energy and wave-vector are E 0 and k0 before it strikes the target, Ef and kf when it leaves the target, and E1 and k, after i collisions. The frequency transfer w,, wave-vector transfer Q,, scattering angle ~ and are given by scattering cosine .u~,at the ith collision,

fined by the lines ~.i = +1 and i~ = —1. The maximum and minimum values of a, for a particular choice of 13 and e~,are given ~y 2(eo+13)1/2. (1) a÷(f3,~ 2e~+ 13 ±2e~/

(..O~ =(E~—

~

12~=k~~_k~

DETHT

1, cosp~=p.=k.

DETWO

k./k.

Using the cosine rule, we obtain =

+

TARGET

k,~ 2k~_1k~p1.

Z.C..,1~~LECJ)

—

now define the following dimensionless quantities: 2Q,~/2mflkBT, 13~ =flw,/k~T, a~ il = Ei/kBT = 1l2k~/2mnkB 7’, We

where m~and kB are the neutron mass and Boltzmann’s constant respectively, and T is the temperature. Combining the above equations, we obtain the relation

WO

ZEUP ~-~--

“~

0MONITOR

XE

C3DIS

XZCPM Fig. 3. The geometry of the Monte Carlo “experiment”. Distances XZCPS and XZCPM are measured from the zero covariance point (Z.C.P.) in the direction of the incident beam, The detectors have height DETHT and width DETWD. The elastic coherent scattering coordinate system is defined by the axes XE, ~E and ZE. As drawn, the scattering angle ANGLE(J) is negative.

292

.J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

Note that w and Q are the natural choice of independent variables in a neutron scattering experiment. In the case of an isotropic target, the variables a and 13 are equally acceptable.

YZSEP

~ T

x

2.2. Geometry

TSEPWO

The geometry of the Monte Carlo “experiment” is illustrated in fig. 3. The incident beam of neutrons has direction cosines WO in the target coordinate systern (T.C.S.) and (1,0,0) in the elastic coherent scattering coordinate system (E.C.S.C.S.) (see below). The detectors are cylinders, with axes normal to the scattering plane, equidistant from the target. A detector angle is positive if neutrons bound for the detector were scattered to the right at the target position. Distances are measured from the zero covariance point (Z.C.P.), where the time and inverse velocity distributions of the incident beam are considered separable (see ref. [3] p. 154, and ref. [5]). A monitor detector is placed in the straight-through beam as shown in fig. 3; the time distribution of neutrons at the monitor is determined within MSCAT and may be compared with an experimental distribution, The (right-handed) elastic coherent scattering coordinate system is shown in fig. 3. This coordinate system is used within subroutine ELSCAT (Section 4.7.3). The target consist of a sample located within a container (fig. 1). Both the sample and the container are isotropic materials. The present program treats the multiple vertical cylinder geometry which is IIlustrated in fig. 4. This type of geometry is commonly used in situations where the ideal geometry would be a thin slab sample, but for engineering reasons such a sample cannot be made. The (right-handed) T.C.S. has its origin at the center of the left-hand tube (viewing along the incident beam direction) and the coordinate axes are shown in fig. 4. Note that WO > 0 always whereas —1
Y T

I T5CP

0\ \

_________________

N~~ON

SDIAN

THICK

IDIAM

Fig. 4. A plan view of the multiple vertical cylinder geometry used in the program. The target coordinate system is defined by the axes XT, ~T and ZT.

,

2.3. Scattering properties We divide the scattering by a material into three parts; elastic coherent scattering (E.C.S.) which is characterized by delta functions in both Q and w;

elastic incoherent scattering (E.I.S.) which contains a delta function in ~ and inelastic scattering (I.S.) which contains no delta functions. We shall consider each type of scattering in turn, and restrict our attention to isotropic systems. The double differential cross sections for E.I.S. and for I.S. are written as proportional to UB, the bound atom scattering cross section per scattering unit. Note that this scattering unit must be the same as in the definitionof p, the number density of scattering units. 2.3.1. Elastic coherent scattering The differential cross section for E.C.S. from an isotropic target is given by d a d~ZdEfE ~

~

=

p

0 Tj

~ 1 .

Xo(E —L )b[2irr.—2k sin~p]. f 0 i 0

The indexj runs over all reflections, r~is the magnitude 2 is the strucof a reciprocal lattice vector, and 1F ture factor for the jth reflection. The0(rj)1 second delta function expresses the Bragg condition. Note that the r’s of Bischoff et al. [4]are larger by a factor 2ir. We may instead write the cross section as a sum over sets of equivalent reflections i, as follows:

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

d2a

E.C.

=

IF(i~)I2exp[_2W(T~)} 2r~

AEC

293

The total cross section is then —

UEI(EO) —AEI

2~2GB 2

2

[1

2

exp(_BEIkO/27r )].

—

(4)

X S(Ef—EQ)5[2lrTf—2k 0 sin~p].

(2)

an average Debye—Waller factor exp[—2W(r,)], where W(r~)= ~BECT~, has been separated from the structure factor. Furthermore an “elastic coherent multiplicative factor” (E.C.M.F.)AEC has been introduced. This constant may used toAEC absorb un2.Thebeproduct 1F(r~)I2 wieldy factors in IF(r~)I is in barns. As an example, the reflections from a Bravais crystal may be described by writing AEC = Ucoh/41T and IF(r~)~2 = M~,whereUcoh is the coherent scattering cross section per scattering unit (in barns), and M is the multiplicity of the ith reflection. The total elastic coherent cross section UEC(Eo) is obtained by integrating eq. (2). We obtain (cf. ref. [61): Here

—

2ir2

UEC(EO —A~p

X ~

IF(r~)I2 exp[—2W(r~)],

2.3.3. Inelastic scattering In the present context the “inelastic scattering” (I.S.) is the total scattering, less the E.C.S. and the E.I.S. The differential cross section (for an isotropic system) is normally written 1/2 as d2u I a~(E~\ S(Q,w), =——~—~ d~dE~ 11 4ir \ E~i where S(Q, w) is the well-known scattering function (or scattering law, or dynamic structure factor), described by Van Hove [7] . Note that the scattering function is a function of the independent variables Q and w. It is related through Fourier transformation to the correlation function G(r, t) of Van Hove. We now introduce the dimensionless scattering function S(a,13), given by S(a,13) =1~~S(Q, w)k~T,

r~~ k

and the symmetrized scattering function (S.S.F.) 2S(a 13). S(a,13)e~ In terms of S (a,13) we define the alpha-integrated

0/ir, which may be written GEC(EO)=(AECIEO) ~P(T~),

BEIkO

T~~k

(3)

0fr,

scattering functionF(a,13)by F(a,13)=fS(a’,13)da’.

where

0

We introduce two intermediate scattering functions (I.S.F.’s),

~2

P(r1) = 2iii 2ir ~

exp[—2W(r1)].

2.3.2.elastic Elastic incoherent scatteringcross section is The incoherent differential simply given by d~2dEf~EJ = AEI

(5)

-~-

11(13, e0)

=

72 [F(a(13’

f

eo),13’)

0 e~~ —

F(a ~

dj3’,

(6)

exp(_BEI Q2 /8ir~)6(E~_ E 0),

where AEI is an “elastic incoherent multiplicative factor” (E.I.M.F.), the exponential term is the Debye— Wailer factor, and Q = 2k0 sin ~ ~.The E.I.M.F. is generally the ratio Ujflc/OB, where ~inc is the incoherent scattering cross section per scattering unit,

which is only defined for 13 D’(13, ~

=

J’ e

~

~

0, and

[F(a÷(j3’,e~)~ j3’)

13 —

~(a (13’,e0),13’)] d13’ —

which is only defined for 0 ~ j3 ~ —e0.

(7)

294

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

The marginal up- and down-scattering functions UQ3,

0) and

2.6. Linear interpolation

D(jl,e0) are given by We define the linear interpolation function LIN by

U(13,c0) and

=

U(13,c0)/U(°°,e0), with 13

?

0, LIN(X,x1,x2,y1,y2)

D(j3,e0) = D (13, e0)/D

(—~~

with 0 ~ 13 ~ —en.

,e0)

(x2—X)y1

+

(X—x1)y2

~ —x 2 1

Furthermore, the upscattering probability function is 3. Program structure

U’(°°, e~) —

[U’(oo,e0)

+D’(—e0,0)]

‘

and finally the total inelastic scattering cross section is simply a~e~=~a ~~u~°° e0,~+D(—e 1’ 0~’ “ B’/4e0’1 \ ‘~ 0’ e O~‘ That is to say,

(

=

a1(e0)

=

~f f GB

~*

—60

13 CO

) ,

,

,

,

S(a ,13 )da d13.

a(j3 ~

The functions S, F, U, D, U0 and a1 are collectively known as the scattering functions (S.F.’s).

2.4. Cross sections Atomic cross sections (in barns) are written as where the first subscript denotes the type of process, and the second subscript denotes the mate,rial. 1)

Similarly, absorption cross sections (in cm— are writtenlinear as ~xy’ Note that in general =

x,y

x,y’

where p~,is the number density of material y. The subscript x can be T, for total (i.e., scattering plus absorption)

S, for all scattering EC, for elastic coherent scattering El, for elastic incoherent scattering I, for inelastic scattering, or A, for absorption. The subscripty can be S, for the sample, or C, for the container, 2.5. Random number The symbol ~ denotes a random number chosen from a uniform distribution between 0 and 1.

In order to give a general feeling for the program, we present in fig. 5 a diagram which shows the relationships between the various subprograms. The subprograms are organized into groups, designated by rectangles of dashed lines. They are individually described in Section 4.. The main program, which we call MSCAT, is used to set the dimensions of various arrays (see Section 6) and to define the unit numbers of the various devices.

MSCAT calls MSC2, which maintains overall control

of the computations. This subroutine divides naturally

into three sections: Data Input, Monte Carlo Loop and Output. In the first section, Data Input, subroutines MSC2, TARGET, and RDTIME read in data for the run. Subroutine MSCIN2 is then called. This routine reads

several cards and then either generates, or reads from a storage device, the scattering functions. If the scatteringroutine, functions arePRELIM to be generated, SGENER is called. This with and SINPUT, generates the symmetrized scattering function S (a, 13); the remaining scattering functions are then computed within MSCIN2. If there is elastic scattering, the necessary information is read in; READEL reads elastic coherent scattering parameters. Control returns to Section II of MSC2.

The first step in the Monte Carlo loop is the mitialization of a neutron. Its direction and impact point are determined by entry IMPACT in TARGET. The removal probability, (l_exp[_~T cdc_~T ~ 1S determined by calling DIST (in TARGET) f~rds and dc (distances to travel through sample and container respectively), and by calling (as necessary) ESINC and/or SIGEL for the elastic incoherent and elastic coherent contributions to the removal cross sections ~T,S and ~T,C’ The inelastic and absorption contributions to ETS and ~T,C are computed within MSC2. A scattering point is now selected (by RANPOS in TARGET), and a type of scattering is chosen,

295

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

.~

.~

I

____________

U

i

Li

I I

i I

I

I

L

Ui

I ~

-

in

2 ~

i

0

0

_____

fl~fl

~-

~tI)

I

I

I

I-

—4

I I

L

-l

I I I

I-

C.)

I

U. 2

b

I I

I~~jII

j_4~j_ L___J

_---~ ____________

~

I~I I~I____________

I

~

_____i~ri~i r

-

- ~

‘~

L

a a

h.,L—’-~iuI I

1!i

i i L__...4

~ ____

______

________

________________

‘~ r.’I -

________ri

•E-~

.;.~i

‘Ui,

I

ICDI

~ ‘I

I

__________ Ui

=~~i~ri L.~J

2

.E

E.EB

a-I

lol 2 I I~I L~JLU a.

—.I~I I

2

I

a,

______________ 1~

L

II I I I

:

J

=

296

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

knowing the various contributions to the total scattering cross section ~ss or ~ c• Depending on the type of scattering, one’of the r~utinesELSEL (for E.C.S.), INCSEL (for E.I.S.) or SEL2 (for I.S.) is called, and a new energy and scattering direction for the neutron are selected in an appropriate manner. Subroutine SEL2 calls SELBET which selects a value for the energy transfer 13. At this stage, the program calculates the contribution of the neutron to the detector response, for each scattering angle, and for each time channel. These quantities are computed within ELSCOR (for E.C.S., E.I.S.) or SCORE (for I.S.). SCORE calls DFSCT to obtain the double differential scattering cross section: DFSCT calls SLAW2 to obtain the appropriate value of ELSCIN S(a,j3). ELSCOR calls ELSCAT (for E.C.S.) and (for E.I.S.). Subroutine BINFND determines the elastic scattering time channel. Both SCORE and ELSCOR call DIST (in TARGET) in order to determine d 5 and dc, the distances to travel through sample and container in the direction of a detector. These quantities are used to determine the attenuation of the beam leaving the target. The neutron now proceeds to further collisions in a similar fathion. At each collision its weight WGT (initially unity) is reduced. If WGT
.

.

Table 1 Alphabetical list of subprograms and additional entry points. The group number of the subprogram is also given. Subprogram Entries Group number BLNFND CYLGEO DIRECT ELSCAT ELSCIN ELSCOR ELSEL 1 ESINC 1 FLTRNF INCSEL 1 INDX KEEP

7 RANCGE

4 4 7 7 7 5

RANGET, RANSET

MSCAT2 MSCIN2

10 115 9 1 3

PRELIM RANDAN

3 4

RDTIME READEL SCORE SECOND SELBET SEL2 SETMS SGENI~R SIGEL SINPUT SLAW2 SREF TARGET

9 3 6 9 5 5 READMS,WRITMS 3 3 33 AINPUT, BINPUT 6 8 COORDS, DIST, IMPACT, 4 NEWDIR,RANPOS,TEFF

GET

TIMZRO Function subprogram. 2 Main program.

8

4. Subroutine and function subprogram descriptions The subroutine and function programs, with their entry points,are listed alphabetically in table 1. Unless otherwise stated, all computations are performed in single precision. 4.1. Group 1: the main program, MSCAT The main program MSCAT contains various COMMON statements which are used to set the dimensions of the scattering function arrays (Section 6). Each array is allocated a separate COMMON block,

and in the subroutines the arrays are given dimension unity. Furthermore, two-dimensional functions, such asS (a,13), are treated as one-dimensional arrays within the program. The remaining statements in MSCAT defme (i) the variables NAC and NBC (see Section 4.3), and (ii) the peripheral unit numbers NRD, NPRT, NST, NSTC and NST2. Subroutine SETMS is called, in order to define the dimensions of the arrays in the routine, and finally subroutine MSC2 is called. If and when control returns to MSCAT, the program terminates.

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

4.2. Group 2.• the control subroutine, MSC2 This routine divides naturally into three sections: Section I, Data Input and Initialization (Statement numbers between 1000 and 1999), Section II, The Monte Carlo Loop (2000—2999), and Section III, Final Data Manipulation and Output (3000—3999). A simplified flowchart is shown in fig. 6. In Section I, seven data cards are read. Subroutines TARGET and RDTIME are then called to read further data cards: the card input is described in Section 5.1. Control now passes to MSCIN2, which handles the ______

C_STA

~‘

READ

/

I

GENERATE __________

AND

I

~

NEUT

i

•

~

NO

1_~

—

can be a to troublesome correction it is clearly advantageous have the facility to runand additional histories.

EXIST - 0?

1NO NEUT YES

—

—

—

_____________________

~ INITIALIZE’

I

————n SAMPLE

~> P~?

ing (fig. 6). The time remaining, CPTIM, is returned by subroutine SECOND (Section 4.9.3), and is compared with the input parameter, TIMIN. If CPTIM < TIMIN, subroutine KEEP is called, and the necessary

2

program then terminates, but the may beunit. conarrays variables onrun a storage tinued and at some laterare timewritten by specifying IEXIST The 0

YES

______

I

MODIFY

Iii

_______ WGT

CHLSE I

I

SAMPLE

NO

_______ E.C.S.

NO > ~EC

~ SCORE c~WGTI

CcLLISiON~~_INELASTI

•

wGT

(Section 5.1). If on the other hand CPTIM ~ TIMIN, and if IEXIST * 0, entry GET (in KEEP) is called, and the necessary arrays and variables are read from

NO SCORE NNP?

L

~~~>l?

I

~C~Nv~

ELASTIC

~YES

________________

NDEVNDEV + i

ERRORS CALCULATE’

I

~NO

_~ —

The variable NNP (fig. 6) is NN/NSTD. Before entering the inner main loop, a test is performed to ensure that there is sufficient time remain-

SAMPLE

dAZIMUTHAL I ANGLE

NEUTRON -

.

are generally runThe in order optimizefrom the statistics of the simulation. elastictoscattering a container

TIMIN

•

scattering properties of both sample and container. It returns the scattering functions and the elastic scattering parameters for both materials. The remainder of Section 1 of MSC2 contains straightforward calculations of the scattered energy and the total cross section (scattering plus absorption) for both sample and container, for each time channel. The basic structure of the Monte Carlo section, Section II of MSC2, is due to Bischoff [3] The total number of neutrons to be tracked, NN + NEL (Section 5.1) is divided into groups of NSTD neutrons, where NSTD is the number of neutrons to be treated together in standard deviation calculations (ref. [3], p. 56).1 The from to NNPP outer=main (NN DO + NEL)/NSTD; loop variable, theNEUT, inner main runs DO variable, NDEV, additional runs from elastic ito NSTD. When thereloop is elastic scattering, histories

________

NDEV

~ NEUT+1 NEil-

PP EUT>

We now follow the history of a neutron. Entry the storage unit. IMPACT, in TARGET, initializes the neutron’s direction and impact point. Its initial energy and time-ofarrival are then chosen from truncated gaussian 2/F2], Ix x~I
NO

—

____J

01>SF,

L__1 ~

297

)\~J’ WRITE

______

Fig. 6. A simplified flow diagram for subroutine MSC2. The top and bottom lines correspond to Sections I and III respectively. The two rectangles, composed of dashed lines, represent the two main DO loops in Section II of MSC2.

where ing theP(x) valueis x, the x0(unnormalized) is the mean value probability of x, andofF choosis the standard deviation of the distribution. The neutron’s weight WGT (initially unity), is now multiplied by (I TRNS), where TRNS, the transmission probabili—

298

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

ty, is given by TRNS = exp[—~2TSdS ~T,CdCI, and d~and dC are the distances to travel through sample and container. A collision point is chosen by entry RANPOS, in TARGET, and the weight WGT is multiplied by (~S,S/~T,s) or (~S,C/~T,c)’ depending where the collision occurs. The neutron’s contributions to the single scattered detector response in each time channel, and for each scattering angle, are now evaluated. Subroutines ELSCOR and SCORE score the elastic and inelastic scattering responses respectively. The weight WGT is next compared with WCO, the cut-off weight (Section 5.1). If WGT > WCO, the neutron’s history is followed through another collision; if WGT
In Section III the dat~is normalized and converted to the employed S.S.F. form,using procedures analogous to those in the program ANALYS [2] ,which is used to reduce our measured data to that form. The data is finally printed as described in Section 5.3. Calculations of the standard deviations of the detector responses are performed in double precision. 4.3. Group 3: the scattering properties ofthe target. ESINC, MSCIN2, PRELIM, READEL, SETMS, SGENER, SIGEL, SINPUT This group of subroutines includes a control subroutine, MSCIN2; three routines, SGENER, PRELIM, and SINPUT, which handle the S.S.F. ~‘(s,13);READEL, which reads and processes the E.C.S. information; two function subprograms, SIGEL and ESINC, which evaluate elastic coherent and elastic incoherent cross sections respectively; and SETMS, which reads or writes the scattering functions on a peripheral device. 4.3.1. Subroutine MSCIN2 This subroutine first reads two cards (Section 5.1). The second card contains the variables IMAT, NSTOP, and NSTOR, which control the overall flow of con-

elastic coherent,

Read Cord I

<~
Read Card 2 _____

______

I_NOWALL[1—..~_MAT _______ IF

The neutron’s new energy, and the scattering cosine, are selected by calling ELSEL, INCSEL or SEL2, as

I

T

Error 9B81 STOP

FINISH 1T

appropriate. Finally, a random azimuthal angle is selected by NEWDIR (in TARGET). The flag, NCLF, (initially 1) is now set to 2, indicating that subsequent collisions constitute multiple scattering events. A new collision point is selected, and contributions to the elastic and inelastic scattering responses are computed. Again WGT is compared with WCO and the history is either terminated or followed through another collision. In the latter event a new energy and direction are chosen and the whole procedure is followed once again. Sooner or later the neutron’s history is terminated. (A test is performed, in Section I of MSC2, to ensure that WCO> l0~6). The above procedure is repeated until the supply of neutrons is exhausted, or until there is insufficient time remaining. If elastic scattering is included, extra neutrons are generally tracked, in order to optimize the statistics of the Monte Carlo run.

_______

o

-]

STOP

k-El

READ F

Read Card2C

TI

I

I

I I

IMAT IF

Error 994 STOP

L

lINE #0

I

Generate S.F. $

I

4

________

____________

lINE ~ ~

I

~t

________

~[

~r

I__I_Printout

$

IINEC

‘READINSTORE)

__________

IF

READ(NST)

I

NSTOR>O

l__i~IL_IMATC * 0

__________

IREADINSTC)

READF

4F

Error 9

]

o

________

WRITE(NSTORE)I

_____

RETURN

F’

Fj FINIS~j

Fig. 7. The overall logic of subroutine MSCIN2. The letters T and F stand for .TRUE. and .FALSE. respectively. The logical variables NOWALL, FINISH and READF are explained in the text. The error numbers refer to format statements in the subroutine.

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

299

Table 2 Types of run, defined by the logical variables READF, FINISH and NOWALL. The letters T and F stand for .TRUE. and .FALSE. respectively. A Generation Run includes the generation of the scattering functions. Designation of run

READF

FINISH

NOWALL

IMAT

Generation, sample (no container) F Generation, sample or container F Generation and production, sample (no container) F

T T F

T F T

+0 0, * 0 *0

Illegal combination Printout, sample (no container) Printout, sample or container Production, sample (no container) Production, container with or without sample

F T T F F

F T F T F

0+0 * 0 0, * 0 *0 *0

F T T T T

trol through the subroutine, as shown in fig. 7. If IMAT ~ 0 (= 0), the variable NSTORE is set equal to NST (NSTC). NST and NSTC are the peripheral unit numbers for the S.F. files for sample and container respectively (defined in MSCAT). Depending on ICOH and IINC, the E.C.S. and E.I.S. data are read from cards (the former by READEL). The logical variable NOWALL is now tested. This variable is .TRUE. if the wall thickness of the contamer is zero. If NOWALL is .TRUE., IMAT is then checked (IMAT must be nonzero for the sample). If NOWALL is .FALSE., the variable FINISH (= NSTOP .NE. 0) is tested. If FINISH is .FALSE., the run is a production run involving a container. In this case the variable READF ( NSTOR .LT.0) is tested. If READF is .FALSE., implying that the S.F.’s are to be generated, an error message is printed, since this combination is not permitted (table 2) and the program terminates. If on the other hand READF is .TRUE., a data card (MSCIN2, card 2C) is read for the container, and then IMAT and IMATC are checked (IMATC must be zero). The E.C.S. and/or E.I.S. data are then read (if this is specified by ICOHC and/or IINCC). If lINE ~ 0, entry READMS (in SETMS) is called to read in the S.F.’s for the sample from storage unit NST, and if IINEC ~ 0, READMS is called to read in the S.F.’s for the container from unit NSTC. Control is then returned to MSC2. We now consider the remaining situations: NOWALL = .TRUE., and/or FINISH = .TRUE. (see fig. 7 and table 2). The variable READF is checked, and if it is .TRUE. the S.F.’s are read from file NSTORE. If READF is .FALSE., and lINE 0, the S.F.’s are generated. This section of the program (Section ID is described below. Once the S.F.’s have been gener-

IMATC — —

— —

— —

0

ated, they are written on file NSTORE (if NSTOR> 0). The S.F.’s are then printed under control of the variables NOUT and NSONLY (MSCIN2, Card 2). If FINISH is .TRUE., the program terminates; if it is .FALSE., control returns to MSC2. Section II of subroutine MSCIN2 contains the calculation of the scattering functions. First the S.S.F. S(a,13) (designated SC in this subroutine) is generated by SGENER. The scattering function is stored in histogram form, such that SC(I,J)* is the value of S(a,J3) for AC(1— l)0, and BC(l)> 0. The arrays AC and BC, containing NAC and NBC values of a and j3 respectively, are also returned by SGENER. The function ~(a,13), designated F in the program, is calculated next. F(I,J) is simply the sum of terms SC(II,J)*(AC(II) AC(II—l)), for II = 1 to I. We now enter a ioop over the energy index IE. First the I.S.F.’s U’(13,e0) and D’(j3,c0), designated —

—

UPI and DNI, are computed. They are obtained [eqs. (6), (7)] by integrating the true scattering function S(a,13) over the region in (a’,13’) space bounded by the lines 13’ = 0,13’ = 13, and a’ = a~(/3,~ In the program, the lower alpha bound of the kinematically *

An element of SC, which is designated SC(l,J) in this writeup,is equivalent to the element SC(I+(J—1)~NAC)in the program, where NAC is the working dimension of the array AC. appliesand equally to the S.F.’s F, the UPMAR, SCC,This FC, comment and UPMARC, in modified form to S.F.’s DNMAR and DNMARC, which are triangular arrays.

300

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

~7 BC(4,

BC(2

L

—

The marginal upscattering function, for zero energy, is now evaluated. From eq. (1) we have a~(13,e0)~13±2(e013)~’2ase0~0,

IE 4

—

________

—

—

so that from eqs. (5) and (6),

—

U(13,0)~4~0 je’/2~(13’,13’)~d13’, as e~0.

—

(8) Since U(j3,0) is simply the ratio U’(13,0)/U’(°°,0), we drop the coefficient 4\/~in e9. (8), and evaluate UPI as the line integral of S(a, 13)131/2 along the line a = 13. UPMAR is then obtained. The final section of MSCIN2 contains statements which print the various scattering functions, as well as other quantities. 0

C

—

— —

—

—

—BC(2~—

—

—

—

—

\,

v.—

-BC(4.

—

—

—

,•/

~—

___.

_—.—-

—~

7/

~c~i

~c~3

4C(6)

AC(8)

Fig. 8. A schematic diagram showing how the kinematically allowed region, for a particular energy, is approximated. In this case C0 = BC(4).

allowed region (for energy C~)is approximated by a series of straight lines joining the points a (13,e0), for 13 = —BC(IE), —BC(IE—1), -- - —BC(l), 0, BC(l), • -. BC(NBC); this is illustrated in fig. 8. The upper alpha bound is similarly approximated. The IB’th element of the array UPI is then written as UPI(IB) = where TERM(J)

,E

=

TERM(J)*EXP(—0.25*(BC(J—l)+BC(J))),

(BC(P—BC(J—l)) \ /

NAC *

~ 1=1

FRAC(l,J)*(AC(I)—AC(I—l))*SC(I,J).

In this expression, FRAC(I,J) is the fraction of the rectangle bounded by the lines a = AC(I—l), a = AC(I), 13 = BC(J—l), 13 = BC(J), which lies within the (approximate) kmnematically allowed region (see fig. 8). A similar expression is used for DNI(IB). The arrays TERM and FLAG do not appear in MSCIN2. Having calculated U’ and D’ for energy 60, the remaining scattering functions are calculated. These functions are U(f3,eo),D(13,e0), and U0(e0), designated UPMAR, DNMAR, and UPROB. The inelastic scattering cross section, SIGMAC, is also calculated.

4.3.2. Subroutines SGENER, PRELIM and SINPUT Within MSCAT we represent the symmetrized scattering function in histogram form. The S.S.F. is designated SC (or SI), and the arrays defining the bounds of the histogram bins are designated AC (or Al) and BC (or BI). (The names in parentheses are the names used in SGENER and SINPUT.) On the other hand it is more convenient to specify the S.S.F. at a mesh ofpoints in (a,13) space; furthermore this meth need not coincide with the mesh defined by the arrays AC and BC. Within these three subroutmes, the S.S.F., designated SE, is specified on a mesh defined by the arrays AE and BE (such that AE(1) = BE(1) = 0). The purpose of these subroutines is to generate the arrays AE, BE, and SE, as well as Al and BI, and from them to calculate SI. The one-dimensional problem is illustrated in fig. 9. The “external” mesh, X for! = 1, -- -,NJ, (with X1 = 0) is the set of values of the abscissa for which the ordinate is defined. These points (X1, )~)are represented by open circles. The internal array, x1 for i = 1, - - - ,NI, (with x1 >0 and x0 = 0) defines the histogram channels, and we wish to know Y(X) dX

.1

y1. =

for i = ito NI.

—

i

i—i

We choose to represent the function Y(X) by its value in each interval (X1, Xj+1), as shown in fig. 9. The integral is then written as: average

.J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

Looking at fig. 9, and defining “internal region” i

__________________

I 2 k~ I I Ij I I m1 I 2

/ 0]

:~

3 I 2 3

4 2 2 5

NI’

5 6 3 4

to be the range x11
I

,

~

‘~

6

1. “External region”! comprises all or part of “internal regions” k1 through l~. 2. “External region”! contains “incremental regions”

4

NJ~~ NKL’ 6

_______

t

f

m1 through m1 + k1. 3. The length of the abscissa in”incremental region” ~ is The Kronecker delta in eq. (9) states that “external region”! only contributes to “internal regions” i fork.
.~

X

X1 X2 X3

X4

e2 e, e4

NT EXT

12

NCR I

2

3

I

3 4

3 I4

2

e7 ~

~

I

I

X5

5

5

I

6

I

4

1

6~

7

.

The two-dimensional internal function S~ is related to the external function ~ by

~j

Fig. 9. The conversion of a function, specified at a set of points, to a function represented by a histogram. In this example the input function is shown as open circles,with coordinates (X1, Y1), (X2, Y2), ---, (X7, Y7). The function Y(X), for other values of X, is represented by the dashed line. The histogram representation is shown by the solid lines. The various symbols are explained in the text.

—

~

—

r ~

‘—‘

~ “

1

2(~~+l)

/

(9)

h=k ~hi~m1+h —k,’

where

and
xk. 1 i

are defined by the inequalities and


/

x1. ~
j

i

m1 is defined by the recurrence relation m1

=

1

‘

m1~1=m1+l1—A~+l, for!= l,”-,NJ—l, and NKL is defined by t~tr,L

— —

—

1,

.~

ii

x

— ‘~‘

=~-—

1,

,

~

“l~

if xNI
where Jf is such that ~ < <~ Jf_ 1 XNI Jf Finally &~= ~,where the array ~ is the union of arrays x and X, with = 0. —

+~ +~ “j+l,j’ j,/’+l

j+1,j’+l

6hi~h—k+m

~II~

~lj~

h=k

~

~‘ h’=k~’

~‘

/

This calculation is performed within subroutine SGENER. The necessary arrays 5, A, AC,B and BC,

NKL

~

(x~—x~1)

j,j’

11

1 =

-

y~~_.1)

*

~

301

are obtained by calling SINPUT and its entries, AINPUT and BINPUT Subroutine PRELIM is called in order to calculate the arrays k, 1, m, and A new subroutine SINPUT must be written for each pçoblem using a different scattering function. In the present deck we include a subroutine which calculates S(a, 13) for a mass one ideal gas. In every case it is required that AE(l) = BE(l) = 0, and that ~.

AI(I) >0, BI(l)> Furthermore the final elements of these arrays must0.satisfy the condition M(NAI)~(3+2~~ñ)~BL(NBI). This requirement is explained in Section 5.5. In addition the elements of each array (AE, BE, Al and BI) must be arranged in increasing order. Two further restrictions apply to the element BI(NBI). It must be greater than (~o)m~, and it must be greater than (Ef)max - Here (Eo)max is the maximum value of the incident energy, denoted by EPMAX in subroutine MSC2, and (f)max is the maximum scattered energy, corresponding to the minimum time channel, and denoted by ENERGY(K) in subroutii3e MSC2

302

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

(between statements 1310 and 1320). The calculation of S(a,13), within subroutine SINPUT, is straightforward. Precautions are taken to avoid the singularity at a = 0. This is achieved by setting S(a1 ,I3~)= S(a2 ,i3~),for all i3~~ where a1 = 0, and a2 is the first nonzero value of a. 4.3.3. Subroutine READEL This subroutine is used to read in the E.C.S. parameters, which include the E.C.M.F.AEC, vectors 2,and thethe Debye— Ti, the structure factors F(r~)I Wailer coefficient BEC. These quantities are defined in Section 2.3.1. In addition the variable XNT is read. This variable is defined in Section 4.7.3. The quantities P(r~)and ~ P(r 1) are also calculated (Section 2.3.1). They are used in subroutines SIGEL, ELSEL and ELSCAT. 4.3.4. Function subprograms SIGEL and ESINC These function subprograms evaluate the total E.C.S. and E.I.S. cross sections respectively. Eqs. (2) and (3) are used, 4.3.5. Subroutine SETMS This subroutine handles the I/O of the scattering functions. Entry WRITMS handles the writing and entry READMS handles the reading. The dimensions of the various arrays are set by a call to SETMS from MSCAT, the main program. 4.4. Group 4: the geometry of the target. CYLGEO, DIRECT, RANDAN, TARGET The subroutine TARGET is written to treat vanous aspects of the sample geometry, which is ilustrated in fig. 4. CYLGEO is a more general routine which can be used with any multiple cylinder target geometry. Some of the calculations in TARGET, and all of the calculations in CYLGEO, are performed in double precision. Subroutines DIRECT and RANDAN complete this group. 4.4.1. Subroutine TARGET This subroutine contains seven entry points, TARGET, COORDS, IMPACT, DIST, RANPOS, NEWDIR and TEFF. All aspects of the target geometry are contained within this routine and within subroutine CYLGEO which it calls. To treat a different .

target geometry it is only necessary to write a new subroutine TARGET. If the different geometry involves a multiple cylinder arrangement, CYLGEO may be used as it stands. The present routine treats a system of up to 120 cylinders, with axes vertical, arranged at an angle to the incident beam (fig. 4). We may distinguish two situations. (i) TSEPWO (see fig. 4) TDIAM. In this case the tubes do not shadow each other. To make subsequent calculations WDINC width of the incident beam) is easier, adjusted so that(the an integral number of tubes is illuminated. (ii) TSEPWO < TDIAM. In this case the tubes shadow each other. The operation of the subroutine is described by corn~‘

ment cards in the deck, and no further explanation need be given here. 4.4.2. Subroutine CYLGEO This subroutine contains two entry points, CYLGEO and RANCGE. On calling CYLGEO (X, W, SSAM, SCAN, IT), the straight line trajectory of a neutron, starting at position X (in tube IT), traveling in direction LV, is traced. The distances from X to all interfaces (between air and container and between container and sample) are calculated and stored in the array S(4, 120). The total distance traveled through the sample (SSAM) and through the container (SCAN) are returned. The array S(i,j) is defined in the following manner: S(1 ,j) = distance from starting position to entry point into container, !th tube. S(2,j) = distance from starting position to entry point into sarnple,!th tube. S(3,j) = distance from starting position to exit point from sample,!th tube. S(4,/) = distance from starting position to exit point from container,!th tube. The tube within which the neutron starts (tube IT) has! = 1. The next tube (if there is one) through which the neutron passes (tube IT ±1)hasj = 2, and so on. In this notation, SSAM

=

~

[S(3,j~

—

S(2,!)],

—

S(l

and N

SCAN

=

I~[S(2,!) j=1

,/) +

S(4,j)

—

S(3,j)],

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

MI ~

CI

/ -

53

.~-

~ MG

(a)

(b)

(c)

Fig. 10. Various trajectories for a neutron travelling through a cylindrical sample surrounded by a container. In (a) the neutron starts within the sample; in (b) it starts within the container; in (c) it comes from outside. A heavy dot indicates that the neutron leaves through an end. Path Ml differs from path MS because a neutron following path MS would never enter the sample, no matter how long the cylinder was. This is not true of path Ml.

where N is the number of tubes intersected by the neutron. We now examine the various possible trajectories, and first consider the different ways the neutron can leave the initial tube (tube IT). (a) Initial tube.’ starting point in sample (fig. 1 Oa): S(1, 1) ‘S(2, 1) 0, Case Si - Neutron leaves through end: S(3, 1) = S(4, 1). Case S2. Neutron passes into container, leaves through end. Case S3. Neutron passes into container, leaves through side. (b) Initial tube: starting point in container (fig. lOb): S(1,1) —0. Case Cl - Neutron leaves through side: 5(2, 1) = 5(3, 1) = S(4, 1). Case C2. Neutron leaves through end: S(2, 1) = S(3, 1) =S(4, 1). Case C3. Neutron passes into sample, leaves through end: S(3,I)~’S(4,I). Case C4. Neutron passes into sample, back into container, and leaves through end, Case CS. Neutron passes into sample, back into contamer, and leaves through side. In all cases except S3, Cl and C5, the neutron leaves through the end of the tube and cannot pass through another tube. In cases S3, Cl and CS, tests are performed to determine whether the neutron will pass through another tube, (c)Subsequent tube (tube IT±(f-i)). (‘fig. lOc). Case Ml. Neutron passes into container, would pass .

303

into sample if tube were infinitely long, but leaves through end: S(2,/) = S(3,!) = S(4,j). Case M2. Neutron passes into container, into sample, and leaves through end: S(3,!) = 5(4,!). Case M3. Neutron passes into container, into sample, back into container, and leaves through end. Case M4. Neutron passes into container, into sax back into container, and leaves through side. Case MS. Neutron passes into container, misses sample, and leaves through end: S(2,!) = S(3,!) = S(4,!). Case M6. Neutron passes into container, misses sample, and leaves through side: S(2,j) = S(3,j) = S(4,i). The procedure for calculating the values ofS(i,!) is straightforward, once the particular case is identified. These calculations are of necessity performed in double precision. Having called CYLGEO to obtain the distances S(i,j), a subsequent call to RANCGE (X, SSS, SMU, CMU, IT) is used to choose a new collision point. Here SMU and CMU are linear absorption coefficients for the sample and the container respectively. The routine returns the new collision point (X), the new tube number (IT), and the distance traveled (SSS). Cumulative transmissions are first calculated: T(0) = 1, T(4!—3) = T(4!—4)*exp(.--CMU* [S(2J)—S(l T(4!—2) = T(41—3)*exp(--SMU* [S(3,j)—S(2,j)J), T(4j—l) = T(4j—2)*exp(—CMU*[S(4,/)—S(3,j)}), T(41) = for! = 1,2, - . , N. Cumulative removal probabilities are then obtained as ~!)J)~

CR’1~=

for i

—~———~~-

1 —T(41V)’

‘‘

=

0 1 2 ‘

-“

4N

‘

and the region in which the collision occurs is selected by solving the equation CR(k ~ < ~ < CR’k —

‘~

/

“

‘.

fork, where ~ is the usual random number. If k is even (odd) the collision has occurredin the sample (container). If ~k is even, an error message is printed. The position where the collision occurred is then given by ~ ~ ~k k d + *[ ~ 1’ 2~+

304

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

where k1

=

k

—

1, mod 4

+

1,

If W = (0,0, 1), no rotations are necessary, and W5

=

k2int(~k)+l, and

h’s.

d=~ln 1—~[1—T(4N)] T(k—l)

4.4.4. Subroutine RANDAN This subroutine is used to calculate the sine and cosine of a random angle. The method is due to von Neumann [9].

and ~(k)

=

SMU if k is even, CMU if k is odd. 4.5. Group 5: Selection of the result ofa scattering event. ELSEL, INCSEL, SELBET, SEL2

4.4.3. Subroutine DIRECT Given an incident direction LV, and a scattering cosine CO, this subroutine calculates a scattered direction WS, such that W- WS = CU, and the azimuthal

When a scattering event occurs, one of these routines is called in order to select the neutron’s new

angle U is a random angle between 0 and 2ir. The final expression for WS is obtained in the following way. A new coordinate system X’ is obtained from the original coordinate system X by two rotations: X’

=

energy and scattering cosine.

[1

0

0

4.5.1. Subroutine SEL2 This subroutine is called if an inelastic scattering event is to be selected. Given an incoming energy E, the scattering cosine CO and the outgoing energy ESC are chosen, (a) The first stage is to choose between up and

10

W3

—z

down scattering, The upscattering probability UP(E) is obtained as

Z

W3)

R1R2X,

where

R1

=

0

UP(E) = LIN [E,L~ , E~,UP(E~1),UP(E~)], where the index i is chosen such that E~i
and W2 —W1 0 W1 W2 0

B.2

0

0

Z)

(b) A value for 13 is now selected. for (i) energies Downscattering. E Generally the function DNMAR, 11 and E~,is interpolated to energy E, and a value for 13 is selected by choosing a random number ~and interpolating in beta. This sequence of

112.In this new coordinate

and Z denotes (W? + W~) system the incident direction W’ is (0,0,1). A scattered direction W~is then constructed: /

W~=

(

sin ~ ~ ~ sin ~ sin 0 cos~i

,

I

where cos p = CU, and the scattered direction in the original coordinate 1R 1 w~,system is then W~= R~ 1 which is -

I(W

~+ 1

wi

cos ~‘

\

operations is performed by subroutine SELBET. There is however one special case which must be considered separately. If i = 0, the maximum possible value for beta is E/kT, and a value for beta is selected by writing p = ~E/kT. (II) Upscattering. The function UPMAR, for energies E~i Et, is interpolated and a value for j3and is selected by choosingtoaenergy randomE,number ~and interpolating in beta. This sequence of operations is performed by subroutine SELBET. (c) The scattered energy ESC is now computed. (d) Finally a value for a is selected. The maximum

2 sin p cos 0 + W1 W3 sin ~ sin 0)(W~+W~) 2+W W~ ‘~(—W 1sin~cosO+ W2W3 sin~sinO)(W~+W~)~ 2cos~pI. and minimum values amax = a+(13,eo) and amin = —(W~-I-W~ sinp sino)+ W3 cos4 I a(13,e0), are first evaluated. Values of~(a,13),for

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

a = amax and a = amin, are obtained by interpolation, and a value of~(a,13) is then chosen at random by writing F(a,13) = F(a 13) + ~[F(a ,j3) —F(a 13)1. .

mm’

max

mm’

~(akl

1

,p),~(ak,p),ak 1

‘13) <~(a,0) ~ ~(ak ~

4.5.2. Subroutine SELBET This subroutine is called by SEL2. Given a section of the marginal up- or downscattering function for two adjacent energies, the distribution for an intermediate energy is obtained and a value for 13 is then selected. In fig. II we show the function for two adjacent energies, which we shall denote ~j1 and E1. For each value of 13 in the internal rne~hBC, we obtain D(E, f3,) = LIN [E,E11, E.,D(E~1, 13,~),D(E1,

i3~)1.

This function is shown as the dashed line. We now choose a random number ~, and find 1, such that D(E, ~/-l~< ~
=

LIN [~, D(E, 13/ —1) ,D(E, ~ P~_~ ~i3~].

E

10

~C/)L~

where ~=

J1

____

I~’ ,~

/

with P1~ = 0,

J1

and n (IP in the subroutine) is such that r~~ k0/ir < In other words n is the number of possible reflections for energy E0. The scattering angle is then calculated. The scattering energy is set equal to the incident energy. 4.5.4. Subroutine INCSEL This subroutine is called if an elastic incoherent scattering event is to be selected. In this case we solve the equation 2I 1 exp [_BElk~(l —p)/4ir for the scattering cosine p = cos ~, by writing 4~2

(~X\

~,

I /n ~P(r,)/~P(r1),

~1_exp[_BEJk~/2T2]

________________

4cr U z ~

4.5.3. Subroutine ELSEL This subroutine is called if an elastic- coherent scattering event is to be selected. The index i of the selected reflection is obtained by writing p; <~~
—

1

0

In practice D(E, 13,) is only computed in the region of interest, as illustrated in fig. 11. In the case of the downscattering function, the interpolation mD cannot be performed forj = i, since D(E~_1,/3~)does not exist, Instead the final value of the interpolated function, which13~_ we designate D(E, I3~), is defined such that the points ( 1,D(E, 13~_ (E/kT, 1.0), and (p~,D(E,13k)) are colinear.

Finally a is obtained by writing a = LIN[F(a,p),~(ak where k is such that

305

~E.

I-I

2]}). p

BEJk~ I +____~_In(1_-~{1_exp[BEIk~/27r

__________________

BETA — Fig. 11. A marginal scattering function is shown by crosses for two energies, E~ 1and E~.The solid circles and the dashed line represent the function for an intermediate energy E. The intersection between the dashed line, and the solid line labelled ~, determines beta. Given ~, it is only necessary to evaluate the marginal scattering function for energy E over the range of beta values indicated by the short vertical bars.

The scattered energy

is

set equal to the incident energy.

4.6. Group 6: the double differentialscattering cross section for inelastic scattering. DFSCT, SCORE, SLA W2 4.6.1, Subroutine SCORE This subroutine is used to calculate inelastic contributions to the intensity in each time channel, for

306

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

each detector angle. A typical contribution is f WGT* ~ cr‘d~2dE 1(E.) *A

where WGT is the current statistical weight of the neutron, and ATT is the probability that the neutron, having left the scattering point, will reach the detector. The expression between asterisks is the probability of scattering into the selected detector and time channel. The differential cross section is simply / d a t1B(~f’~ S(a,13) ex (_13/2~ d~ldE~ 4ir \E1/

kBT

~

“

4.6.2. Subroutine DFSCT This subroutine is used to evaluate the differential cross section. 4.6.3. Subroutine SLA W2 This subroutine is used to obtain S(a,13), given a and 13. 4. 7. Group 7: the differential cross sections for elastic scattering. BINFND, ELSCA T, ELSCIN, ELSCOR These routines are used to calculate elastic contrib. utions to the intensity in each time channel, for each detector angle. 4.7.1. Subroutine ELSCOR This subrotuine has overall control, and is used to compute quantities to bothofE.C,S. and E.I.S. These include ATT,common the attenuation the beam on leaving the scattering point; CO, the scattering cosine; arid TARV, the time-of-arrival at the detector, 4.7.2. Subroutine BINFND This subroutine is used to find the index of the elastic scattering time channel, given the time-ofarrival TARV. 4.7.3. Subroutine ELSCAT This subroutine is called to calculate the elastic coherent scattering cross section for a particular detector. The cross section is written DFEC

=

L~~ZAEC ~

~JP(r),

r-~~ k1, r,

where ~&2is the solid angle subtended by the detector. Thus DF is essentially ~c,/z~l (in barns per steraEC sum over reflections is necessary because dian). The more than one reflection can enter the detector. For the jth reflection,f1 is the fraction of the Debye— Scherrer cone which intersects the detector. This fraction is determined by Monte Carlo sampling. A direction hEC, lying on the Debye—Scherrer cone, is chosen at random. If this vector passes through the detector, the number of successes (ISCS) is incrernented. This procedure is repeated NT times: (NT is the same as the input parameter XNT, read in by subroutine READEL). The fractionf1 is then ISCS divided by NT. If NT ~ 0, DFEC is set to zero, and control returns from ELSCAT. Note that this routine assumes that the detectors are cylindrical, with axes normal to the scattering plane, and that they have 100% efficiency. The most probable Bragg reflection vector is tried first. Then the vectors and are tried, and so on, until two vectors have yielded zero fractions. Bischoff [31uses a different procedure to calculate the intensity DFEC, but we prefer the present method since it is more readily adapted to a different detector geometry, and it involves no approximations. On the other hand it is an inefficient method. 4.7.4. Subroutine ELSCIN This subroutine is called to calculate the elastic incoherent scattering cross section for a particular detector. The cross section is written 2(l—p)/4ir2], DFEI = AE I ~-exp[—B 4ir El k0 where p is the scattering cosine. 4.8. Group 8: the initial parameters of a neutron. SREF, TIMZRO 4.8.1. Subroutine TIMZRO This subroutine is called, in the case of an imperfect resolution run, to choose the neutron’s initial inverse velocity n 0 and its time at the zero covariance point, T0, from truncated gaussian distributions. 4.8.2. ThisSubroutine subroutineSREF is used to choose a value x from

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

a gaussian distribution centered at x0 of width i~.x. The method is due to Box and Muller [101. 4.9. Group 9.’ miscellany. KEEP, RDTIME, SECOND 4.9.1. Subroutine KEEP This subroutine handles the I/O of various quantities which are modified on each pass through the outer Monte Carlo loop in subroutine MSC2. When KEEP is called, sufficient data (including the current value of the random number generator) is dumped to storage device NST2 so that the run may be continued on a later occasion. When entry GET is called, the data is read from the dump and the current value of the random number generator is reset. Care must be exercised in redefining NN and NEL (MSC2, Card 2) if a run is to be extended (rather than continued on account of time limitations in a previous run). If NN’ and NEL’ further inelastic and elastic histories are desired, the new values of NN and NEL must be given by NN(new) = NN(old) NEL(new) = NEL’.

+

NEL (old) + NN’,

4.9.2. Subroutine RDTIME This subroutine is used to read in various quantities associated with the time-of-flight analysis scheme. The detector channel width is taken to be independent of channel number but this restriction is easily lifted simply by rewriting RDTIME. 4.9.3. Subroutine SECOND This subroutine returns the current time remaining for the job, in seconds. This routine calls an installation dependent routine, TLEFT, and must be rewritten for use at any other installation. If no suitable routine exists, a statement such as T = I .OE + 6 will make the program run, but the time remaining will not be correctly tested, 4.10. Group 10: the random number generator. FLTRNF The function FLTRNF returns a random number ~ such that 0 ( ~ 1, with probability F(~)= 1. The numbers are generated using the multiplicative con-

307

gruence 16+11),modulo231. rx~(2 This is a modified version of subroutine RANDU, described in the IBM System/360 Scientific Subroutine Package [111, page 77. If this routine is to be used on a computer having a word length of N bits (N> 32) it may be modified by including statements which mask out the first (N— 31) bits of the variable IX. This is described more fully by comment cards in the routine. Xn+l

4.11. Group 11.’ an interval-halving routine. INDX The function INDX returns the index i, such that ~ given X and the array x containing n elements which satisfy the conditions 0~ x 1, and x~1~ x~, for i = 2,”’, n. The interval halving technique is employed. See, for example, ref. [111, p. 250, where this technique (called “binary search”) is illustrated. 5. Input and output Table 3 lists the subprograms involved in the input of data from the card reader and from storage devices, and in the output of data to the line printer and to storage devices. We also indicate separately those subroutines which can print warning or fatal error messages. We shall consider each group in turn. 5.1. Input from cards The data read from cards is listed in table 4. Seven cards are read by MSC2. Card 1 is a comment card which typically identifies the run. Card 2 contains NN, the number of neutrons to be tracked in the run; IRANX, the initial value of the random number generator in subroutine FLTRNF (if IRANX is read as zero, it is set to 3125, and if the run is a continuation run, the current value of the generator is read in by entry GET in subroutine KEEP); NSTD, the number of neutrons to be treated as a unit for the purpose of calculating standard deviations (NN/NSTD must be an integer,> I); NEL, the number of additional elastic

308

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

Card 5 contains EPO, the incident energy in meV; GAMTOF, the standard deviation of the inverse yelocity distribution of the incident beam (in ps/cm);

Table 3 Subprograms which include input or output statements. The letter C opposite a subprogram name means that the subprogram includes statements to read from cards. Similarly the letters I, P, 0, and E respectively refer to other forms of input, output to the line printer (other than error messages), other forms of output, and output of error messages to the line printer,

.

Subroutme name

and GAMT, the standard deviation of the time of arrival at the Z.C.P. (in ps). Card 6 of MSC2 contains JANG, the number of detector angles (maximum 8), and card 7 contains the detector angles in degrees (the sign convention is described in Section 2.2). One card is read by TARGET. It contains TLGTH, the length of a tube; TRAD, its radius (including the container); TSEP, the separation between tube centers; THICK, the wall thickness; and NTB, the number of tubes. All dimensions are in cm. See fig. 4.

_____________

______________

BINFND CYLGEO ELSCAT ELSCIN

______

E E E E

S INCSEL INDX KEEP MSCIN2

E E I C

P P

0 E

RDTIME READEL SELBET SEL2 SETMS SGENER SIGEL TARGET

-

The first card read by RDTIME contains DELTM(1), the channel width (in ps) for all detectors; CHWD, the channel width (in ps) for the monitor; and KMX, the number of time channels (limit is 50). Card (or cards) 2 contain the energy transfers (Ef E0) in meV, arranged in decreasing order. Cart! 1 of subroutine MSCIN2 is a title card. Card 2 contains a number of parameters which normally apply to the sample. Under certain circumstances, --

C

p

I

C

P

0

E E E E E E E

histories, which must be a multiple of NSTD; TIMIN, a time in seconds such that if the time remaining is less than TIMIN, the necessary numbers are dumped by subroutine KEEP, and the run terminates; IEXIST, which is nonzero if the run is a continuation run, and is zero if it is a new run; IRES, which is zero if a perfect resolution run is desired, and nonzero if an imperfect resolution run is desired; and WCO, the neutron’s cutoff weight. WCO must be ~ 10_6 and should be made> 1 if only single scattering is desired (If WCO = 1, some multiple scattering is scored because of the “Russian Roulette” method of terminating a neutron’s history.) Card 3 contains WDINC and HTINC, the width and height of the incident beam (in cm) and WO, an array which contains the direction cosines of the incident beam (see Section 2.2). Card 4 contains the distances XZCPS (from Z.C.P. to target), XZCPM (Z.C.P. to monitor) and XSD (target to detectors), and the height and width (diameter) of a detector, DETHT and DETWD (all quantities in cm) (see fig. 3).

.

as described in Section 4.3.1 (see also table 2), they may apply to the container. SIGBA is the bound atom scattering cross section c~(in barns); TEMP is the temperature in Kelvin; SIGABS is the absorption cross section (inbarns) for neutron energy ENABS (in meV). The absorption cross section 112, RHONofis the material is assumed to vary as E~ the density in g/cc and WTMOL is the mass of a scattering unit in amu. If WTMOL = 0, RHON is instead the number density (scattering units/cc). The variable IMAT is nonzero (zero) if the numbers on this card refer to the sample (the container). The next three variables, lINE, IINC, and ICOH, are nonzero if and only if the material is an inelastic, an elastic incoherent, and/or an elastic coherent scatterer respectively. NSTOP = 0 if the run is to be a production run. NSONLY * 0 if only the S.S.F. is to be printed out. NSTOR = 0 if the S.F.’s are to be generated but not written on a storage file; NSTOR> 0 if they are to be generated and then written on a storage file; and NSTOR < 0 if the S.F.’s are to be read from the storage file. NOUT controls the degree to which the S.F.’s however,

are printed out (see Section 5.3).

If ICOH was nonzero, subroutine READEL reads in cards which describe the elastic coherent scattering.

J,R.D. Copley, Monte Carlo calculations for thermal neutron scattering

309

Table 4 Input data from cards. The variables are explained in Section 5.1. The variables RUNID, TITLE and WO are dimensioned 20, 20 and 3 respectively. Subroutine

Card(s)

Variables

Format

MSC2

1 2

RUNID NN, IRANX, NSTD, NEL, TIMIN, IEXIST, IRES, WCO

3 4 5 6 7

WDINC, HTINC, WO XZCPS, XZCPM, XSD, DETHT, DETWD EPO,GAMTOF,GAMT JANG (ANGLE (J),J=l,JANG)

(20A4) (18,112,215, FlO.5, 215, F10.5) (SF10.5) (SF10.5) (3F10.5) (IS) (8F10.5)

TARGET RDTIME MSCIN2

1

READEL

MSCIN2

1

TLGTH, TRAD, TSEP, THICK, NTB

(4F 10.5,15)

1 2ff.

DELTM(1),CHWD,KMX (EN(I),I=l,KMX)

(2F10.5,15) (8F10.5)

1 2

TITLE SIGBA, TEMP, SIGABS, ENABS, RHON, WTMOL, IMAT, lINE, IINC, ICOH, NSTOP, NSONLY, NSTOR, NOUT ITM (TA(I), 1=1, ITM) (SF(I), I’l, ITM) FELAST, DWINT,XNT DW1,FIN SIGBAC, TEMPC, SIGABC, ENABSC, RHONC, WTMOLC, IMATC, IINEC, IINCC, ICOHC

(20A4) (4F1 0.5, E10.4, F10.S,I2,5I1,12,l3) (12) (7F10.6) (7F10.6) (3F10.6) (2F10.S) (4F10.5, El0.4, Fl0.5, 12, 311)

1 2 ff. 3 ff. 4 32

2C3 READEL3’4

1C-4C

See READEL, Cards 1—4

MSCIN23’5

3C

DW1C,FINC

1 2

(2F10.5)

Card is only read if ICOH (MSCIN2, Card 2) ~ 0. Card is only read if IINC (MSCIN2, Card 2) 0.

~ Card is only read if THICK (TARGET, Card 1) > 0, NSTOP (MSCIN2, Card 2) = 0, and NSTOR (MSCIN2, Card 2) <0. See also table 2. Card is only read if ICOHC (MSCIN2, Card 2C) + 0. Card is only read if IINCC (MSCIN2, Card 2C) ~ 0.

Card 1 contains ITM, the number of reflections (limit 99). Card(s) 2 contains the tau-vectors r~(in A1), and card(s) 3 contains the structure factors IF(r 2 1)1 (see Section 2.3.1) for each reflection. Card 4 contains FELAST, the E.C.M.F.; DWINT, the coefficient BEC in the Debye—Waller factor in (A2); and XNT, which controls the scoring of E.C.S. by subroutine ELSCAT (Section 4.7.3). If IINC was nonzero, a third card is read by MSCIN2 (card 3). It contains DWI, which is the coefficient BEI (see Section 2.3.2), and FIN, which is the E.I.M~F,AEI. If the run is to be a

production run with a contamer, further cards are read in. Card 2C of MSCIN2

is similar to card 2, but the variables SIGBAC, TEMPC, SIGABC, ENABSC, RHONC and WTMOLC refer to the container, IMATC must be zero, and the three remaining variables IINEC, IINCC and ICOHC also refer to the container. If ICOHC is nonzero, cards 1C-4C are read in by READEL, describing the E.C.S. by the container. If IINCC is nonzero, card 3C of MSCIN2 is read, describing the E.I.S. by the container. 5.2. Other forms of input If the variable IEXIST (MSC2, Card 2) is nonzero, the run is a continuation run; entry GET in subroutine

310

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

KEEP is called, and sufficient data is read from storage unit NST2 so that the run can be continued, If NSTOR (MSCIN2, Card 2) <0, and lINE and/or IINEC is nonzero, entry READMS in subroutine SETMS is called, to read in the S.F.’s from storage unit NST and/or NSTC. 5.3. Output to the line printer (excluding error messages) In this section we give a general description of the printed output, excluding error messages. The bulk of the output is generated by MSC2 and by MSCIN2. First the contents of all the data cards are printed, The data cards are identified as in this write-up. The next section of printed output is from MSCIN2, under control of the parameters NOUT and NSONLY (MSCIN2, Card 2). If NOUT ( 2, there is no further printout from MSCIN2, Otherwise several input parameters are printed; then if there is E.C.S. the vectors Ti and the relative probabilities P(r 1) are printed, and —

if there is E.I.S. the quantity DW1 (i.e., BEI) is printed. Finally the S.F.’s are printed (if there is I.S.). First the internal a and 13 arrays, and the corresponding Q and t~ arrays, are written out. The functions S, F, U and D are now printed unless NOUT < 0. If NSONLY * 0 only Sis printed. Every kth value is printed, where k = NOUT + I. Thus if NOUT = 3, the 1st, 5th, 9th, etc. elements in each of the arrays is printed. The following functions are now printed: U0, a1, o~,tJA + 0EC + 0E1’ 0T~and the energy mesh. Appropriate messages are printed by subroutine KEEP as and when the routine is called to read or write the running arrays on storage unit NST2. Section II of MSC2 also contains several timing messages. Section III of MSC2 contains the bulk of the printout. The first statements print some of the input data. For each time channel, the total cross section, the time-of-flight from sample to detector, and the channel width are printed. The number of collisions NCOL, the number of collisions in the container, NCOLC, and the number of neutrons lost, NLOST, are then written out. (A neutron is considered lost if it is upscattered to the point where its energy is greater than the maximum value in the internal f3 mesh). The initial and final values of the random number generator are printed, followed by the computed transmission of the target, and the average computing

times (for the run) per collision and per neutron. The time distributions at the target position and at the monitor are then written out. The final results are now printed. A separate page is devoted to each angle, and the information for each time channel, printed on a separate line, includes: the scattered energy, the wave vector transfer Q, the frequency transfer w, the inelastic first scattering “response”, the elastic first scattering “response”, the total (elastic plusinelastic) multiple scattering “response” the total (single plus multiple) scattering “response”, the correction factor (i.e., the ratio of all multiple scattering to all single plus all multiple scattering), and finally the “ideal response”, which is simply S(Q, w). All the “responses” are presented in the symmetrized scattering function form, in units of l0~2 sec. The calculated “responses”, and the correction factor, are presented with their associated standard errors. 5.4. Other forms ofoutput If there is insufficient time remaining, or when the Monte Carlo 1oop is completed, subroutine KEEP is called and sufficient data is written on unit NST2 so that the run can be continued at some other time. If NSTOR (MSCIN2, Card 2) > 0, and lINE >0, entry WRITMS in subroutine SETMS is called, and the S.F’s are written on unit NST or NSTC. 5.5. Error messages The subroutines which can write error messages are listed in table 3. In every case the print statements are grouped at the end of the routine, and the error message itself includes the name of the routine. The message “EPMAX.GE.EN(KPX)”, in subroutine MSC2, indicates that EN(KPX), the maximum energy in the internal beta mesh, is smaller than the maximum possible incident neutron energy: the program cannot proceed because cross sections at energy EPMAX cannot be evaluated. Similarly the message “E .GE. EN(KPX)” in MSC2 indicates that there is a time channel for which the scattered energy E is> EN(KPX). The message “COMMENSURATE MESHES,..” requires some explanation. In a perfect resolution run, first scattering is scored for definite values of the energy transfer (the values read in by RDTIME); if one of

.J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

these values coincides with a value of 13 in the internal beta mesh, the scattering function is not well defined (the energy transfer lies at the edge of a histogram bin). This condition has caused problems in the past and for this reason it is not allowed. In subroutine MSCIN2, the messages “NOWALL .EQ.TRUE, AND LINE, IINC AND ICOH ALL .EQ. 0” and “NOWALL .EQ.FALSE AND IINEC, IINCC AND ICOHC ALL .EQ. 0” should be self-explanatory. The messages “NSTOP .EQ.0. . “IMAT .EQ. 0” and “IMAT .EQ. 0 AND/OR IMATC .NE. 0” correspond to errors 992,988 and 994 respectively in fig. 7. If there is inelastic scattering by both sample and container, the S.F ,‘s-for each material are read from a storage device. If the alpha meshes for the two materials differ, the message “ALPHA, . MESHES DIFFER” is printed. A similar message is printed if the beta meshes differ. The messages “CALCULATION OF UPI’ and “CALCULATION OF DNI’..,” appear if an error occurs in the calculation of the I.S.F.’s U’ orD’. Since these messages should never appear, we will not explain their significance. An error message is printed if there is no scattering from the sample, but SIGABS > 0. Such a situation would cause problems elsewhere in the program. Subroutine SGENER prints several self-explanatory error messages. The remaining message, “AI(NAI) LESS THAN appears when the maximum value in the internal a array is less than (3 + 2.J~)times the maximum value in the internal 13 array (I3max)’ Under these circumstances upscattering from e = 0max to = 20max cannot be properly treated because the a mesh does not extend to large enough values. Bischoff et al. [4] incorrectly state that the a mesh need only extend to 40max’ The subroutines READEL, SIGEL, ESINC, ELSEL, INCSEL, ELSCAT and ELSCIN use the variable K to denote sample scattering (K = 0) or container scattering (K = I). IfK takes any other value an error message is printed. READEL contains two other selfexplanatory error messages. When entry READMS (in SETMS) is called, the S.F,’s for a particular material are read from a storage unit. If the number of a values, the number of 13 values, the bound atom cross section, or the temperature, disagree with the expected value, error messages are printed. The expected values of the first two quantities are specified in MSCAT, and the expected values .“,

.

311

of the last two quantities are read in by MSCIN2 (card 2 or 2C). Subroutine TARGET contains three obvious error messages. Again CYLGEO contains several self-explanatory error messages, but these messages normally do not appear during a run. If they do, the problem may sometimes be alleviated by increasing the value of the variable SMALL (in a DATA statement). These messages generally occur in a rather specialized situation, such as the case of a neutron which grazes the edge of a cylinder. Subroutines SEL2, ELSEL, and INCSEL contain a message which appears if the computed scattering cosine, p, does not satisfy p~< 1. In addition SEL2 contains a message which appears if the computed scattered energy is negative, and SELBET writes error messages if the random number RAN2 lies outside the range of interpolated values of the marginal scattering function. An error message is printed by BINFND if the array TIME is not ordered properly. Two straightforward messages are printed by RDTIME, and INDX prints messages if the given value X0 lies outside the range of values of the array X. This program is routinely used by the neutron scattering group at Argonne for a variety of problems. We still find on occasion that unexpected error messages are printed. The cause is generally traced to a roundoff problem of some sort. Such problems would probably be avoided if the program were written in double precision. 6. Common blocks and variable dimensions There are 31 common blocks in the program. Those labelled BLOCK- are general purpose common blocks; within each of these blocks the variables are arranged in alphabetical order. BLOCK 1 ,BLOCK2, BLOCK3 and BLOCK4 contain real variables, whereas BLOCKA and BLOCKB contain logical and integer variables respectively. The variables within each of these blocks are described in comment statements in subroutine MSC2. Some of the variables in block BLOCK3 are in double precision.

Common blocks COMONI, COMON2, COMON3 and COMON4 are described in comment statements in subroutines SEL2, TARGET, TARGET and MSCAT respectively. Common block COMON3 contains double precision variables. Common block MULCOR only

312

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

appears in MSC2; it is included for compatibility the program MULCOR described in [2] . Common blocks VCB1 through VCB2O are variable common blocks. These blocks each contain one singly dimensioned array; VCB13 and VCB14 also contam an integer variable. In MSCAT the dimensions of the various arrays are defined (unless the arrays are not to be used), and in the subroutines the arrays are each dimensioned unity. Thus the common block statements in MSCAT set the dimensions of the arrays. The internal a and j3 histogram arrays contain NA! and NBI values respectively. Block VCBI contains the S,S.F. (for the sample) which should be dimen. sioned NAI*NBI. Similarly blocks VCB2 to VCB6

with

contain the arrays F, D, U~U0 and a~(for the sample) which should be dimensioned NAI*NBI, NBI*(NBI+l)/2 NBI*(NBI+1), NBI and NBI respectively. Blocks VCB7— 12 contain similar variables for the container, and should be dimensioned in the same way as blocks VCB1—6 respectively. If there is no I.S. by the sample (container) the blocks VCB1—6 (VCB7—12) need not be dimensioned since these arrays will not be used. VCB13, and respectively, 18 should beand dimensioned NBI, NB!14, and15NBI if there is NAL, a contamer, blocks VCB16 and 17 should be dimensioned NAL and NB! respectively. The dimensions of the arrays in VCB19 and VCB2O should be (NA! + NAE) and (NB! + NBE) respectively, where NAE and NBE are the numbers of values in the external alpha and beta arrays respectively. The arrays AE, KA, LA and MA in subroutine SGENER are dimensioned with 400 elements. The arrays BE, KB, LB, MB, SE1 and SE2 contain 500 elements. If NAE > 400 and/or NBE> 500, the appropriate dimension statements in this routine must be modified.

7. Tests of the program 7.1.

Test

run

It is impossible to test all sections of the program

therefore restrict our attention to a sample with no container. At the end of this subsection we outline a method for testing sections of the program which treat the container scattering. The program includes a subroutine named SINPUT (Section 4.3.2) which generates the mass one ideal ~s scattering function. This particular form for S(a,13) was chosen because it permits a test of those sections of subroutine MSCIN2 which perform the integrations outlined in Section 2.3.3. The computed cross section Ui(0) has been compared with the exact result [12] 0B u~(~)=

1

i—

j—

~ f(e

0+~)erf (ye0)

+

ve0/~exp (—e~)}.

In the worst case, where ~ 1.1, the difference is 0.6%. Bischoff [3] made a similar comparison, and obtained differences of the same order. The test run also calls for both elastic incoherent and elastic coherent scattering. The tau vectors and structure factors describe the Bragg scattering from a body centered cubic material having a lattice param3, eter of 5.7 A.with The this number density,The l.080X 1022 cm is consistent description. E.C.M.F., AEC = 0,79577, is simply OB/41r where °B = 10 barns. The structure factors F(T~)~2 are reflection multiplicities (see Section 2.3.1). The incident energy was chosen so that a set of Bragg reflections would occur at one of the scattering angles, 90°.The remaining input parameters were chosen somewhat arbitrarily. They require no further explanation. Several pages of the printed output are reproduced at the end of this paper. The first page shows the input data. The next page gives the upscattering probability, and various contributions to the total cross section, as a function of incident energy. The third page includes the number of collisions in the sample, the initial and final values of the random number generator, the computed transmission of the sample, and some timing information. There follow the results of the calculation for four angles and six scattered energies. In order to test out sections of the program which ‘~-

in a single run. For example, the program is written

treat the container scattering, the following procedure

way that one cannot, in a single run, compute the multiple scattering effects in a target con-

is suggested. (1) Run the program with NSTOP = and NSTOR = 1, in order to write the scattering functions on a data set. The necessary job control cards

in such a

sisting of a sample surrounded by a container. We

J,R.D, Copley, Monte Carlo calculations for thermal neutron scattering

must of course be included. (2) Run the program with NSTOP = 0 and NSTOR = —1, in order to check that the scattering functions are being read correctly from the data set. (3) If so, proceed to make the following changes: (i) Increase the core requirement by 100 kbytes (25 kwords). (ii) Insert COMMON statements into the main program: COMMON!VCB7/SCC(8000) COMMON!VCB8/FC(8000) COMMON!VCB9/DNMARC(3240) COMMON/VCB1O/UPMARC(6480) COMMON!VCB1 l/UPROBC(80) COMMON!VCB 1 2/SIGMCC(80) COMMON!VCB 1 6/ACC(-l00) COMMON!VCB 1 7/BCC(80) Also insert the statement NSTC = NST, following the definition of NST. (iii) Change the parameter THICK from 0.0 to 0.1 (TARGET, Card 1). (iv) After MSCIN2, Card 3, add the following cards: MSCIN2, Card 2C. Same as MSCIN2, Card 2, except column 62 contains 0 (variable IMATC). READEL, Cards 1C—4C. Same as READEL, Cards 1—4. MSCIN2, Card 3C. Same as MSCIN2, Card 3. (4) The results of this run should be identical to those of the test run, except that the “response functions” will be larger by a factor of 4, owing to the fact that the sample volume has been reduced by this factor (from 0,2 cm radius to 0.1 cm radius). 7.2. Tests against independent calculations In this section we describe two different tests of the program against independent calculations. The inelastic scattering has been tested by calculating the multiple scattering in a sample of 36A. The geometry of the sample, and the details of the “experiment” were chosen to match the experimental setup of Sk6ld et a!. [13] . The results of the calculation were compared with the results of an independent transport calculation, which is described in ref. [13] . The two sets of results show similar behaviour as a function of both angle and scattered energy. In general the transport calculation lies below the MSCAT calculation but there is no reason to doubt the correctness of either calculation. Different approximations

were used for each calculation, and in view of this fact the degree of agreement is quite satisfactory.

313

The elastic incoherent scattering has been tested by computing the multiple scattering in a plane slab of vanadium. Comparing with a transport calculation, using the expressions given in ref. [2] , Appendix 2, we find that the two calculations agree to better than 1% at all scattering angles. The tests outlined above were performed using target geometries which differ from the vertical cylinder geometry which is described in the present paper. Two other versions of TARGET, which treat (i) plane slab geometry and (ii) multiple horizontal cylinder geometry, will be described in separate papers in this journal.

8. Discussion In this final section we offer several comments on the usefulness of the program MSCAT. These remarks are based on more than two years’ experience with the program in various stages of evolution into its present form. The program realizes its full capability, and its maximum efficiency, when applied to inelastic scattering problems. Once the various scattering functions have been evaluated, the calculation of multiple scattering effects is essentially independent of the complexity of the scattering function. Since the calculation is statistical in nature, the accuracy of the results may be improved by continuing the calculation for a longer period of time. The calculation of elastic coherent scattering effects is by comparison very inefficient and time consuming. The most time-consuming section of the program is then the calculation of the intersection of a Debye— Scherrer cone with a detector (subroutine ELSCAT). For this reason we include an option which allows the user to follow Bragg reflections within the target, without calculating intensities of Bragg reflections which reach detectors. This saves considerably on computing time, and allows one to treat the scattering within the target in a proper fashion. It should be noted that if XNT ~ 0 (Section 4.7.3), no detector response is calculated if the final scattering process is elastic coherent, Thus an elastic coherent process, followed by an inelastic process, is evaluated but the reverse sequence is not evaluated. Since the two sequences are presumably equally probable, one might

314

J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

be able to include double scattering processes of the type (inelastic, elastic coherent) by separately storing detector responses to processes of the type (elastic coherent, inelastic) and then doubling the results. If a perfect resolution run is used (i.e., IRES = 0), self-shielding factors (ref. [2] , Appendix 1) may be evaluated by dividing the inelastic single scattering response by the “idea!” response. In general this is

an inefficient method of computing self-shielding factors. On the other hand we have found it useful to evaluate a total cross section UT(Eo) using MSCAT and then to calculate self-shielding factors, or sample attenuation factors, using a simple program which performs the integrals given in ref. [2] , Appendix 1. Note that the self-shielding factor properly includes the effects of shielding by both the sample and the container.

Acknowledgements It is a pleasure to thank F.G. Bischoff for sending a copy of his program, for his advice on several occasions, and for his interest in the modified version of his program. I am grateful to J.M, Rowe and D.L. Price for useful discussions, and to J.M. Rowe for communicating the results of his multiple scattering calculations in liquid argon.

References Ill 121

PA. Egelstaff (Ed.), Thermal neutron scattering (Acadeniic Press, London, 1965). J.R.D. Copley, D.L. Price and J,M. Rowe, NucI. Inst. Meth. 107 (1973) 501. An erratum, relating to the definition of self-shielding factors, is published in Nuci.

Inst. Meth. 114 (1974) 411, 13] F.G. Bischoff. Ph.D. Thesis, Rensselaer Polytechnic Institute (1970). Available from University Microfilms, Ann Arbor, Michigan, Catalogue No. 70-19931. 14] Eng. E.G. Bischoff, M.L. Yeater and W.E. Moore, Nucl. Sci. 48 (1972) 266. [5] R.J. Royston, Nuci. Inst. Meth. 30(1964)184. In this paper the zero covariance point is called the “effective chopping point”. [6] R. Weinstock, Phys. Rev. 65 (1944) 1. See eq. (32). 171 Van Hove, Phys. Rev. 95 (1954) 249. [81 L. H. Kahn,Applications of Monte Carlo, AECU-3259, ~,. 115. [9] J. von Neumann, NBS Appi. Math. Ser. 12 (1951) 36. [101 G.E.P. Box and ME. Muller, Ann. Math. Stat. 28 (1958) 610; See also M.E. Muller, J. Assoc. Comp. Mach. 6 (1959) 376. [11] IBM System/360 Scientific Subroutine Package, Version III. Edition No. GH2O-0205-4 (1970). [12] D.E. Parks, M.S. Nelkin, J.R. Beyster and N.F. Wikner, Slow neutron scattering and thermalization (W.A. Benjamin, Inc., NewGOYork, p. 67. a mass one gas, the 4parameters and A1970) in their eq. For (2.100) become and 1 respectively. oB/ [13] K. Sk~Id,J.M. Rowe, G.E. Ostrowski and P.D. Randolph, Phys. Rev. A6 (1972) 1107.

.J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

315

TEST RUN OUTPUT TIME REMAINING AT START OF PROGRAM ***** INPUT DATA READ IN BY MSCZ

424.79 SEC

CARD

1.

TEST CASE FOR PUBLICATION. S~E LONG WRITE—UP FOR DETAILS.

CARD

2.

NEUTRONS USED 100 INITIAL RANDOM NUMBER = NEUTRONS TREATED AS UNIT FOR STD. 0EV. CALCNS. ADDITIONAL ELASTIC HISTORIES 50 IEXIST = 0 IRES 1 CUTOFF WEIGHT =0.01000000

CARD

3.

WDINC

CARD

4.

01ST ZCP TO SAMPLE

2.50000

DETHT • 5.

EPO

6.

NUMBER OF ANGLES

CARD

7~

ANGLES IN DEGREES —120.000 —90.000

CARD

1.

*8*8*

CARD

1.

CARD

2 PP.

*****

CARD- 1. CARD 2.

*4*8*

CARD CARD

CARD

1. 2 FF.

3 FF.

32.73000

CARD

4. *89*8

CARD

3.

GAMTOF

01ST ZCP TO MONITOR •

300.000

01ST SAMPLE JU LIEIIrCTORS

250.000

5.00000 CM 0.01000

—

GAMY

•

10.00000

4

=

—60.000

—30.000

INPUT DATA READ IN BY TARGET TUBE LENGTH — WAIL THICKNESS

TORE RADIUS • 0.20000 NUMBER OF TUBES = 12

1.00000 0.0

TUBE SEPARATION

•

0.50000

INPUT DATA READ IN BY ROTIME DETECTOR CHANNEL WIDTH • 5.00000 MONITOR CHANNEL WIDTH NUMBER OF CHANNELS • 6 ENERGY TRANSFERS (IN MEV) 5.00000 2.00000 0.0 —2.00000 —5.00000 —15.00000 INPUT DATA READ IN BY MSCIN2 MASS ONE IDEAL GAS SCATTERING FUNCTION. SIGBA — 10.00000 TEMP • 30.00000 ENABS — 25.29999 RHON — 0.IOBOE 23 TMAT 1 lINE • 1 IINC • NSTOP • 0 NSONLY — 0 NSTOR •

—

4.00000

(MICKUSECSI

5490 SIGABS WTMOL 1 ICOH 0 NOUT

5.00000 0.0 • 1 • 9

INPUT DATA READ IN BY READEL ITM 72 ARRAY TA... 0.24811 0.85941 1.13691 1.33610 1.50918 1.64576

0.35088 0.89456 1.16373 1.38140 1.50918 1.66436

0.42973 0.89456 1.18988 1.38140 1.50918 1=66436

0.4962i 0.96092 1.21547 1.40351 1.52944 1.66436

0.55479 0.99243 1.24054 1.42527 5.54943 L.70094

ARRAY SF... 12.0000

6.0000

24.0000

12.0000

48.0000 24.0000 48.0000 48.0000 48.0000

24=0000 48.0000 48.0000 48.0000 48.0000

48.0000 8.0000 6=0000 24.0000 24.0000

24.0000 48.0000 24.0000 24.0000 24.0000

—

2.50000

DETW0

CARD

*4*8*

-

30.000

CARD

25

—

_________________________________________________________________________________

HTINC

50.00000

5480

_________________________

FELAST

=

0.79577

OWINT

•

0.60714 1.02297 1.24054 1.42527 1.56917 1.10094

0.65643 1.02291 1.24054 1.42527 1=58566 1.71894

24.0000

8.0000

12.0000 12.0000 24.0000 48.0000 48.0000

24.0000 48.0000 48.0000 24.0000 48.0000

0.80000

TNT

0.10175 1.05263 1.26510 i.44670 1.58566 1.73675

0.74432 1.05263 1.28920 1.44610 1.60192 1.73615

0.14432 1.08141 1.28920 1.46782 1.61695 1.73675

0.78458 1.08141 1.28920 1.46665 1.62595 1.15439

0.82288 1.10957 1.31286 1.48865 1.62695 1.75439

48.0000

6.0000

12.0000

24.0000

24.0000

24.0000

24.0000 24.0000 24.0000 24.0000 24.0000

24.0000 24.0000 24.0000 24.0000 12.0000

6.0000 24.0000 24.0000 48.0000 48.0000

46.0000 24.0000 48.0000 24.0000 48.0000

24.0000 48.0000 12.0000 48.0000 24.0000

24.0000 48.0000 24.0000 48.0000 6.0000

100.0

________

____________

INPUT DATA READ IN BY MSCIN2 OWl

—

0.52200

FIN

=

1.00000

_____________

FUNCTIONS S.F.DNMAR,UPMAR,UPROB, E SIUMAC CALCULATED THIS RUN UPSCATIERING PROBABILITY 9.9891—01 9.946E—01 9.8671—01 1.138E-01 7.45iF—01 7.1621-01 4.559F—01 4.301E—01 4.0941—01 2.5381-01 2.4261—01 2.320E—01 1.5421-01 1.4861—01 1.432E—01 5.023F—01 9.918E—02 9.620E—02 7.2501—02 7.O62E—02 6.881E—02

9.7561—01 6.873E—01 3.8981-01 2.221E—01 1.382E—01 9.3356—02 6.7011—02

9.616E—Oi 6.585E—01 3.TILE—01 2.1271—01 1.333E—01 9.0631-02 6.539E—O2

9.448E—O1 6.301E—01 3.5336—01 2.0381—01 1.2881—01 8.802E—02 6.3181—02

9.2566—01 6.022E—Oi 3.366E—01 1.955E—01 1.244E—01 8.552E—02 6.2221—02

-

9.0416—01 8.8086—01 5.7S0F—01~548AE—01 3.207E—01 3.0516—01 1.8761—01 1.8021—01 1.203E—01 1.1631—01 8.312E—02 8.0821—02 6.0726—02

8.5516—01 8.2946—01 8.0201—01 5.2306—01 4.9831—01 4.1461—01 2.9166—01 2.7821-01 2 .6561—01 1.7311-01 1.6656—01 1.6021—01 j.lZbE—0r1.0901—D1 1.usb6~or ~.661E—02 7.6506—02 7.4461—02

____________________________

________

-

316

.J.R.D. Copley, Monte Carlo calculations for thermal neutron scattering

INFIASTIC SCATTERING CROSS SECTION (BARNS) 5.650E 01 2.8301 01 1.8931 01 1.4271 01 1.1506 4.9011 00 4.633E 00 4.4051 00 4.2106 00 4.0426 3.2621 DO 3.2071 00 3.1511 00 3.1121 00 3.071E 2.851E 00 2.8321 00 2.815E 00 2.199E 00 2.785E 2.698E 00 2.6901 00 2.6821 00 2.67.RE 00 2.6681 2.626E 00 2.6221 00 2.61SF 00 2.6141 00 2.611E 2.588E 00 2.5851 002.5836 00 2.581E 00 2.519E

01 00 00 00 00 00 00

9.6591 3.8066 3.0346 2.1716 2.6621 2.6071 2.5711

00 00 00 00 ‘JO 00 00

8.3612 3.768!, 3.000E 2.1581 2.6562 2.6041 2.51SF

00 00 00 00 00 00 00

1.3982 3.656E 2.9706 2. 1461 2.6502 2.6012 2.5131

00 00 00 00 00 00 00

6.6582 3.5582 2.9411 2.13SF 2.645F 2.5981

ABSORPTION 3.12SF 02 2.4061 01 1.251E 01 8.4551 00 6.3841 00 5.128E 00 4.2851 00

ONLY (BARNS) 1.5641 02 1.043E 2.2351 (11 2.086F 1.203E 01 1.159E 8.2331 00 8.0211 6.251E 00 6.1341 5.046E 00 4.9661 4.228E 00 4.1111

7.8211 1.955E 1.1111 1.8211 6.016E 4.888F 4.116E

01 01 01 00 00 00 00

6.257E 1.840E 1.0191 1.6301 5.903E 4.BI3E 4.0631

01 01 01 00 00 00 00

5.2 141 1.1381 1.0431 7.4406 5.1932 4.140E 4.0111

01 01 01 00 00 00 00

4.4691 1.6471 1.0091 1.2152 5.688E 4.6691 3.9601

01 01 01 00 00 00 00

3.9101 1.564), 9.7161 7.1102 5.5861, 4.601E 3.9101

01 01 00 00 00 00 00

3.4762 1.4902 9.4801 6.9522 5=4882 4.5342

ABSORPTION 3.22SF 02 3.4001 01 3.230E 01 2.718E 01 2.3641 01 2.182E 01 1.9671 01

PLUS ELASTIC (BARNS) 1.664E 02 1.143E 02 8.820E 4.536E 01 4.2171 01 3.947E 3.1041 01 2.000E 01 3.03SF 2.7021 01 2.630E 01 2.587E 2.3151 01 2.3911 01 2.3421 2.144E 01 2.125E 01 2.1231 1.956E 01 1.9261 01 1.9061

01 01 01 01 01 01 01

7.2562 3.7171 2.933E 2.523E 2.2951 2.0811 1.8181

01 01 01 01 01 01 01

6.212E 3.5171 2.8361 2.5832 2.251E 2.0681 1.8801

01 01 01. 01 01 01 01

5.4672 3.343E 2.7471 2.520E 2.270E 2.0621 1.8541

01 01 Dl 01 01 01 DI

4.9081 3.4112 2.8696 2.4626 2.241E 2.028E 1.828E

01. 01 01 01 01 01 01

4.4131, 3.255E 2.7812 2.4111 2=2001 2.0001

TOTAL CROSS 3.193E 02 3.890E 01 3.556E 01 3.0631 01 2.6331 01 2.4451 01 2.226E 01

SECTION 1.947E 02 4.999F 01 3.4251 01 2.985E 01 2.584E 01 2.4061 01 2.2141 01

02 01 01 05 01 01 01

8.405E 4.12SF 3.240E 2.8021 2.562E 2.348E 2.1361

01 01 01 01 01 01 01

7.1782 3.906E 3.1392 2.8601 2.517E 2.3281 2.1311

01 01 01 01 01 01 01

6.303E 3.7201 3.0472 2.1961 2.536E 2.3226 2.1121

01 01 01 01 01 01 01

5.6486 3.1712 3. 1662 2. 1366 2.5071 2.288E 2.0851

01 01 01 01 01 01 01

5.1382 3.6101, 3.075), 2.7441, 2=4641 2.2601

(BARNS) 1.3321 4.65SF 3•305E 2.9121 2.659E 2.381E 2.18SF

02 01 01 00 00 00 00

02 01 01 01 01 01 01

1.025E 4.36SF 3.3491 2.8611 2.6091 2.384E 2.164F

ENERGY MESH (MFV) 6.4631—03 2.5851—02 5.B11E-,02 1.034E—01 1.6161—01 2.3272—01 3.161E—O1 4.1362—01 1.0921 00 1.267E 00 1.4541 00 1.6541 00 1.8686 00 2.0941 00 2.3331 00 2.58SF 00 4.039E 00 4.369E 00 4.711E 00 5.061E 00 5.4351 00 5.8172 00 6.2111 30 6.6181 00 8.548E 00 9.332E 00 9.8301 00 1.034E 01 1.0861 01 1.140E 01 1.1951 01 1.2511 01 1.5521 01 1.6161 Dl 1.681E 01 1.7481 01 1.8151 01 1.885E 01 1.9551 01 2.027E 01 2.405E 01 2.484E 01 2.565E 01 2.647F 01 2.1311 01 2.8151 01 2.9012 01 2.9881 01 3.4441 01_3.539E_01_3.635E_01_3.733E_01_3.532E_01_3.9321_01_4.O33F_01_4.136E_01

00 00 00 00 00 00

GAMT

=

‘JO 00 ‘70 ‘JO ‘JO (10

5.6062 3.3932 2.892!2.1152 2.6351 2.5932

00 00 00 00 00 00

01 01 00 00 ‘JO 00

3.1281 1.4222 9.201!, 6.8012 5.394!, 4=4692

______

01 01 00 JO 00 00

2.8441 1.3602 8.9381 6.6562 5.3021 4.4062

01. 01 00 (10 (10 JO

2.6011 1.3031 8.6901 6.5112 5.2141 4.3456

01 01 OU 00 00 00

4.1241 01 3.116F 01 2.1001 01 2.4152 01 2~2T01 01 2.0092 01

3.8392 2.9012 2.6151 2.4182 2.1702 1.9786

01 01 01 01 01 01

3.6012 2.8791 2.602), 2.3662 2.159), 1.9581

01 01 01 01 1W 01

—

01 01 01 01 01 01

________

01 01 01 01 01 01

4.1322 01 3.4632 0! 2.9912 01 2.6881- 01 2.414E0l 2.2692 01

___________

4.4001 3.331E 2.9642 2=6891 2.4332 2.2311

01 01 01 01 01 01

MICROSECS.

GAMTOF

0.01000

5.2352—01 2.8501 00 1.0381 00 1.3092 01 2.100F 01 3.0 liE 01

6.4632—01 3. 1281, 00 1.4711 00 1.3682 01 2.1141 01 3.1612 (11

1.8201—01 3.4102 00 1.911), 00 1.4282 01 2.250F 01 3.2582 01

___________________

5480 5490

MICROSECS/CM.

INCIDENT ENERGY~ 32.73 MEV INCIDENT DIRECTION = ( 0.7071068, 0.7071068, 0.0 100 INELASTIC HISTORIES 50 ADDITIONAL ELASTIC HISTORIES 25 NEUTRONS TREATED AS UNIT FOR STO 0EV CALCULATION TUBES SHADOW EACH OTHER TUBE LENGTH 1.000 CM. TUBE RADIUS 0.200 CM. TUBE SEPARATION 0.500 CM. WALL THICKNESS 0.0 CM. NUMBER OF TUBES 12 BEAM DIRECTION COSINES 0.70711 0.70711 0.0 DIMENSIONS OF BEAM 2.500 BY 1.000 CM. FLIGHT PATH= 250.000 CM NUMBIR OENSITY=1.O80000E 22CM**—3 SAMPLE TEMPERATURE=2.5851E 00MEV BOUND ATOM CROSS SECTION THERMAL

ABS0RBTION

CROSS

(SAMPLE)

ALFA

0.0

t.0000E 01 BARNS = 0.1080E 00 CM—i SECTION= 5.0000E 00 BARNS 0.54000—01 CM—I

—

—45.000

CM**—3 (CONTAINER) (SAMPLE) (SAMPLE) (SAMPLE) (SAMPLE)

4.1Z31 01 3.2122 01 2.8891 01 2.6372 01 2.42220T 2.2172 01

_____________

INELASTIC SCATTERING INCLUDED ELASTIC COHERENT SCATTERING INCLUDED ELASTIC INCOHERENT SCATTERING INCLUDED 10.00000

5.2216 00 3.3242 00 2.8116 00 2. 1061 00 2.6311 0O 2.5902 00

-

TEST CASE FOR PUBLICATION. SEE LONG WRITE—UP FOR DETAILS. MASS ONE IDEAL GAS SCATTERING FUNCTION. SAMPLE; SAMPLE; SAMPLE;

6.0152 3.471!, 2.9162 2. 1252 2.6406 2.5956

0.0 0.0 0.0 0.0

BARNS CM—S BARNS CM—I

(CONTAINER) (CONTAINER) (CONTAINER) (CONTAINER)

--

9.3061—01 3.7232 00 8.3162 00 1.4892 01 2.3271 01 3.3502 01

317

J.R.D. Copley, Monte Carlo calculationsfor thermal neutron scattering MACROSCOPIC TOTAL CROSS SECTION (SAMPLE), 2.3246E—0l 2.39421—01 2.41081—01 TIME CHANNELS 9.3051E 02

CM—I 2.44236—01

2.5334E—D1

2.19862—01

IMICROSECS) 9.6981E 02

9.99061 02

1.03111 03

1.08541 03

1.35742 03

CHANNEL WIDTHS (MICRDSECS) 5.0000E 00 5.00001 00

5.0000E 00

5.00001 00

5.00001 00

5.0000E 00

-

CUTOFF WEIGHT=1.0000E—02 TOTAL NUMBER OF COLLISIONS= 448 INCLUDING NUMBER DF NEUTRONS LOST OUT OF ENERGY RANGE= 0 INITIAL VALUE OF RANDOM NUMBER GENERATOR 3125 FINAL VALUE OF RANDOM NUMBER GENERATOR 1731165173 TRANSMISSION

=

0.91905

TIME, PER (COLLISION—ANGLE—TIME CHANNEL), TIME,

PER

(NEUTRON—ANGLE—TIME CHANNEL),

INCIDENT ENFRG1SCAT1 WAVE DELE E1i7 VEC. 11EV ANG—l PS—i 37.73 7.14 7.60 34.73 6.99 3.04 32.13 6.88 0.00 30.13 6.78 —3.04 21.73 6.61 —1.60 17.73 6.00 —22.19 SCATTERING SAMPLF; SAMPLE; SAMPLE;

MS

1.83614

MS

32.13

MEV

INELASTIC 510. SINGLE 0EV. SCATTERING 2.2351—07 7.3891—09 3.724E—W71T~TE—9 5.2321—07 1.040E—08 7.175E—07 1.O3SE—O8 1.1411—06 1.1381—OR 5.2011—06 8.3131—OR

ELASTIC ITO. SINGLE DIV. SCATTERING 0.0 0.0 0.0 0.0 1.0861—01 1.8462—02 1.4151—01 4.2906—02 0.0 0.0 0.0 0.0

ANGLE= —00.00 DEGREES INELASTIC SCATTERING INCLUDED ELASTIC COHERENT SCATTERING INCLUDED ELASTIC INCOHERENT SCATTERING INCLUDED

ENFRGY

SCA1T WAVE OELF EN. VET. MEV AND—i PS—I 31.13 5.83 7.60 34.73 5.71 3.04 32.13 5.62 0.00 30.73 5.53 —3.04 27.73 5.40 —7.60 17.13 4.93 —22.19 SCATTERING ANGLE

SAMPLE; SAMPLE; SAMPLE; INCIDENT

-0.61478

ANGLE= —120.00 DEGREES INELASTIC SCATTERING INCLUDED ELASTIC COHERENT SCATTERING INCLUDED ELASTIC INCOHERENT SCATTERING INCLUDED

SCATTERING SAMPLE; SAMPLE; SAMPLE;

INCITENT

0 WITHIN CONTAINER

32.13

ALL SID. MULTIPLE DIV. SCATTEFING 2.3441—04 2.057E—04 5.2441—05 i.644E—05 8.7502—03 3.1992—03 1.8541-02 1.3102—02 6.3402—04 5.811E—04 1.6472—05 1.9512—06

Sb. TOTAL ~DFV~ SCAITERING 2.3461-04 2.0511-04 5.2811—05 1.6441—05 l.113E—01 1.8136—02 1.6602—31 4.4042—02 6~3~Iff~E4 5.6112—04 2.1676—05 1.9522—06

________________________

_______

_______

ENERGY.

—60.00

32.73

I0EAC~~ 4.6562—01 6.5932—UU 1.4692 (10 1.2641—06 i.6911—0A

1.1812—06

_______________

11EV

INELASTIC STO. SINGLE 0EV. SCATTERING 9.647E—06 l.465E—O1 1.457E—O’59WI2E—07 l.BO4E—05 2.R1RE—O7 2.262E—O5 1.R73E—07 3~095E—O5 6.391E—05 4.BG2E—O7

ELASTIC SED. SINGLE DIV. SCATTERING 0.0 0.0 0.O~OO~ 1.OROE 00 3.O86E—0l 1.122E 00 2.336E—Ol

ALL Sb. MULTIPLE 0EV. SCATTERING 2.IAIE—O5 5.9376—06 3.441E—04 3.1421—04 4.0012—03 1.9461—03 8.352E—03 1.688E—O3

110. TOTAL ~ SCATTERING 3.1312—05 5.9301—06 3.5892—04 3.I4ZE—04 1.0842 00 3.QRbE—Oi l.I3OE 00 2.3362—01

RATIO 0

0.0

6.0081—05 4.3372—05

1.2402—04 4.337F—05

0.48454 0.38864

0.0

TO TOTAL 0.69138 U.9b05~ 0.00369 0.00730

SF0. 0.23060 1.21406 0.00208 0.00214

1.3532—05 1.105 1.3512 01 3.1612—05 4.(IU11-U5 1.3022—05

DEGREES

INELASTIC SCATTERING INCLUDED ELASTIC COHERENT SCATTERING INCLUDED ELASTIC INCOHERENT SCATTERING INCLUOED

SCATT WAVE DELI EN. VEC. 11EV AND—I PS—I 37.13 4.13 7.60 34.73 4.04 3.04 32.73 3.97 0.00 30.73 3.91 —3.04 27.13 3.83 —7.60 11.73 3.51 —22.79

RATIO UI’ 510. MULTIPLE 019 10 tOTAL 0.99905 1.23936 09~2950.43814 0.014S1 0.02013 0.1 1167 0.08451 0.99820 1.29411 0.75997 0.11311

-.________ -

—

11EV

INELASTIC STD. SINGLE 0EV. SCATTERING 4.400E—04 6.5281—06 5.E8RE—04 5.4ORE—06 6.R76E—O4 4.8ORE—06 1.3591—04 9.932E—O6 S.222E—04 5.4381—06 6.773E—04 6.2621—06

ELASTIC 510. SINGLF DIV. SCATTERING 0.0 0.0 0.0 0.0 2.815E—Ol 5.076E—02 5.314E—02 1.9051—02 0.0 0~0~ 0.0 0.0

ALL STO. MULTIPLE DIV. SCATTERING 5.093E—05 I.541E—05 1.163E—03 1.0872—03 l.109E—02 5.9432—03 1.3621—02 R.100E—03 6.4841~’O5 1.829E—05 1.392E—OS 2.734E—D6

510. RATIO Ii!- 510. TOTAL 0EV~~NW’EIPLEDEV~ IOEAL SCATTERING TO 1(7161 4.9181—04 1.6741—05 0.10355 0.03154 6.6011—04 1.732E—03 1.OR7E—D3 0.MI5601563F679W3F-O4 2.9332—01 5.1102—02 0.037(9 0.02130 1.8112 00 6.7492—32 2.1532—02 0.20111 0.13618 8.5866—04 S.8TOE~04T908FO307O73UW.02flE3 9.8711—04 6.912E—04 6.8322—06 0.02014 0.00396 1.9782—04

SCATTERING ANGLE=

—30.00 DECRIES INELASTIC SCATTERING INCLUDED ELASTIC COHERENT SCAT(EATNG INCLUDED ELASTIC INCOHERENT SCATTERING INCLUDED

SAMPLE; SAMPLE; SAMPLE; INCIDENT

ENERGY=

WAVE

EN. 11EV

VEC~~

SINGLE

ANG—1 PS—i 2.15 1.60 2.09 3.04 2.06 0.00 2.03 —3.04 2.00 —7.60 2.05 —22.79

SCATTERING 9.273E—03 1.3371—02 l.419E—02 1.3711—02 l.242E—02 I.2RRE—O3

37.73 34.73 32.13 30.73 21.73 17.73

DELE

32.73

SCATT

INELASTIC

________________

_______

________

11EV

STD. 0EV. 7.161E—05 1.0771—04 l.068E—O4 1.150E—04 l.171E—04 l.929F—05

ELASTIC S10 SINGLE 0EV. SCATTERING 0.0 0.0 0.0 ~ 3.857E—O1 6.2172—02 1.989E—O2 1.9891—02 0.0 0.0 0.0 0.0

ALL 510. ITO. MULTIPLE 0EV. IUTSL 0EV. SCATTERING SCATTERING l.223E—O4 5.5341—05 9.3941—03 9.5321—05 5.535E—O4 3.5631—04 1.3931—01 3.IIZE—O4 3.486E—02 1.036E—0Z 4.3481—01 6.3031—02 4.851E—03 3.7671—03 3.8516—02 2.0242—02 1.492E—05 2.0491—05 1.2491—02 1.1881-04 1.5901—05 l.953E—06 l.304E—03 l.939E—05

RATIO Of- 110. ~UJt,IPLE UEV. 10 IDEAL 0.01298 0.00589 1.1241—02 0.03914 O.OZ960~T75ST1-UT 0.08018 0.02651 1.9706 00 0.12598 0.11813 1.5812—02 0.00800 0.001R4 1.4266—02 0.01220 0.00151 1.3202—03