MHfSH

I M m

By!IB

Eg© V-: ^H m

0**-*Ci.

Ban SCud rail

ft

..

*«li ^H

-rn

mKt

»s

M

^->:
NBS TECHNICAL NOTE 576

Computer Code for the Calculation

of Thermal Neutron Absorption

in Spherical and Cylindrical

Neutron Sources
NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards1 was established by an act of Congress March 3, 1901. The Bureau's overall goal is to strengthen and advance the Nation's science and

technology and facilitate their effective application for public benefit. To this end, the Bureau conducts research and provides: (1) a basis for the Nation's physical measure- ment system, (2) scientific and technological services for industry and government, (3) a technical basis for equity in trade, and (4) technical services to promote public safety.

The Bureau consists of the Institute for Basic Standards, the Institute for Materials

Research, the Institute for Applied Technology, the Center for Computer Sciences and

Technology, and the Office for Information Programs.

THE INSTITUTE FOR BASIC STANDARDS provides the central basis within the United States of a complete and consistent system of physical measurement; coordinates

that system with measurement systems of other nations; and furnishes essential services leading to accurate and uniform physical measurements throughout the Nation's scientific community, industry, and commerce. The Institute consists of a Center for Radiation Research, an Office of Measurement Services and the following divisions: Applied Mathematics—Electricity—Heat—Mechanics—Optical Physics—Linac

Radiation2—Nuclear Radiation2—Applied Radiation2—Quantum Electronics3—

Electromagnetics3—Time and Frequency3—Laboratory Astrophysics3 —Cryogenics3

. THE INSTITUTE FOR MATERIALS RESEARCH conducts materials research lead- ing to improved methods of measurement, standards, and data on the properties of

well-characterized materials needed by industry, commerce, educational institutions, and

Government; provides advisory and research services to other Government agencies;

and develops, produces, and distributes standard reference materials. The Institute consists of the Office of Standard Reference Materials and the following divisions: Analytical Chemistry—Polymers—Metallurgy—Inorganic Materials—Reactor Radiation—Physical Chemistry.

THE INSTITUTE FOR APPLIED TECHNOLOGY provides technical services to pro- mote the use of available technology and to facilitate technological innovation in industry and Government; cooperates with public and private organizations leading to the development of technological standards (including mandatory safety standards), codes

and methods of test; and provides technical advice and services to Government agencies

upon request. The Institute also monitors NBS engineering standards activities and

provides liaison between NBS and national and international engineering standards

bodies. The Institute consists of the following technical divisions and offices: Engineering Standards Services—Weights and Measures—Flammable Fabrics

Invention and Innovation—Vehicle Systems Research—Product Evaluation

Technology—Building Research—Electronic Technology—Technical Analysis

Measurement Engineering.

THE CENTER FOR COMPUTER SCIENCES AND TECHNOLOGY conducts re- search and provides technical services designed to aid Government agencies in improving cost effectiveness in the conduct of their programs through the selection, acquisition,

and effective utilization of automatic data processing equipment; and serves as the principal focus within the executive branch for the development of Federal standards for automatic data processing equipment, techniques, and computer languages. The Center

consists of the following offices and divisions: Information Processing Standards—Computer Information—Computer Services —Systems Development—Information Processing Technology.

THE OFFICE FOR INFORMATION PROGRAMS promotes optimum dissemination

and accessibility of scientific information generated within NBS and other agencies of

the Federal Government; promotes the development of the National Standard Reference

Data System and a system of information analysis centers dealing with the broader

aspects of the National Measurement System; provides appropriate services to ensure

that the NBS staff has optimum accessibility to the scientific information of the world,

and directs the public information activities of the Bureau. The Office consists of the following organizational units: Office of Standard Reference Data—Office of Technical Information and

Publications—Library—Office of Public Information—Office of International

Relations.

1 Headquarters and Laboratories at Gaithersburg, Maryland, unless otherwise noted; mailing address Washington, D.C. 20234.

2 Part of the Center for Radiation Research.

3 Located at Boulder, Colorado 80302.
JO u UNITED STATES DEPARTMENT OF COMMERCE

Maurice H. Stans, Secretary

NATIONAL BUREAU OF STANDARDS,* Lewis M. Branscornb, Director

» * o

NBS TECHNICAL NOTE 576

ISSUED MAY 1971

Nat. Bur. Stand. (U.S.), Tech. Note 576, 25 pages (May 1971) CODEN: NBTNA

Computer Code for the Calculation

of Thermal Neutron Absorption

in Spherical and Cylindrical Neutron Sources

V. Spiegel, Jr. and W. M. Murphey

Nuclear Radiation Division

Institute for Basic Standards

National Bureau of Standards

Washington, D.C. 20234

^°F %

\mj °<">eM of

'

NBS Technical Notes are designed to supplement the Bureau's regular publications program. They provide a means for making available scientific data that are of transient or limited interest. Technical Notes may be

listed or referred to in the open literature.

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C, 20402.

(Order by SD Catalog No. C 13.46:576). Price 35 cents. Stock No. 0303 0851
Abstract

A computer code has been written in FORTRAN IV for

the calculation of thermal neutron absorption in spherical

and cylindrical neutron sources. The formalism of the calculation, the structure of the computer code, a listing of

the code, and some sample results are presented. The comparison of the results of this calculation to experiment

appears elsewhere (l).

Key words: Manganous sulfate bath calibration of neutron

sources; neutron; neutron standards.

11
Computer Code for the Calculation of Thermal Neutron

Absorption in Spherical and Cylindrical Neutron Sources

V. Spiegel, Jr. and William M. Murphey

1. INTRODUCTION

This calculation has been carried out in connection with

a program to reduce the uncertainties in the corrections

applied to the manganous sulfate bath calibration of neutron

sources (2,3). The correction considered here is to account

for the reduction of the manganese activity due to the loss

of thermalized neutrons absorbed in the neutron source itself.

The source may be composed of up to three cladding and

one, possibly fissionable, core material. The calculation

is carried out in a single interaction approximation, i.e.,

the effects of elastic and inelastic scattering of thermal

neutrons are neglected. This approximation is adequate because

the neutrons are in thermal equilibrium and because the correction which is applied to the calibration is small (typically Y% or less). Two cases are available. The first is

for a spherically symmetric source and the second is for a

cylindrically symmetric source. Each consists of a core and

up to three cladding layers. The thicknesses of the ends and

side of a cladding cylinder may all be different. A measurement or knowledge of the thermal -neutron flux at the source

location is required. The thermal -neutron flux is assumed to

be isotropic, which enables one to carry out the computation

as a sum of mono-directional fluxes from different directions.

All integrations are performed with Weddle's formula (4).

2. DESCRIPTION OF THE CALCULATION

Part I. The probability of neutron loss for a given

neutron direction and position.

If the source materials are labeled A, B, C, and D from

the inside out, then the probability of a thermal neutron

interacting in passing through the source is

- Saa(t)-Ebb(t)- 2

cc(t)- £dd(t)

where a(t) is the thickness of "A" material for this particular direction and location of passing through the source,

b(t) is the total thickness of "B" material, etc., and "t"

denotes any particular path through the source. The E's are

the appropriate macroscopic cross sections, £

a being the sum of

Present address is: Technical Analysis Division, Institute for Applied

Technology, National Bureau of Standards, Washington, D.C 20234.
the absorption and fission cross sections for the "A" material,

2 being the absorption cross section for the "B" material,

etc.

The probability of fission in the fissionable material is

given by

, , Eaf , "E

a

a(t) w -(V^^V^^V^)} . \ P

f (t) »_i£ (l-e a

)(e d c ^ ) (2)

S

a

where d»(t) is the thickness of the "D" layer passed in this

direction and location going from the outside into the MA"

material, c'(t) the thickness of the "C" layer passed in

going into the source, etc., and zaf is the macroscopic

fission cross section of material "A".

Part II. Case 1. A spherically symmetric source.

In this case the probability of neutron disappearance

for a neutron striking the source is given by

R

%

2Tiy p(y) dy

(3)

TT R

2

where R is the outer radius of the source and y is the perpendicular distance from the center of the source to the

path through it. The probability of fission is similarly

,R

2TTy P (y) dy

'o f

P

f « (4)

TTR2

The probability of loss of neutrons from the bath per neutron

striking the source is therefore

PL

= P-kPf

(5)

where k is the number of neutrons per fission interaction.

Part II. Case 2. A cylindrically symmetric source.

The notation for the dimensions of the source is given in

Figure 1. Since an isotropic flux is assumed it is necessary

only to consider monodirectional fluxes from an appropriate

number of different directions. By symmetry it is most convenient to use cylindrical coordinates. For convenience in
normalization and avoidance of solid angle considerations an

imaginary sphere is placed about the source and all neutrons

striking this sphere are included in the integrals and probabilities. This procedure does not affect the results when

they are finally expressed in terms of neutrons lost as a

percent of the source sti'ength or in terras of a given thermal

flux in the vicinity of the source.

Given that a neutron enters this imaginary sphere, let

the probability of transmission without interaction be PT, probability of interaction be Pj, probability of absorption

be PA and the probability of fission be Pf

. The probability

of neutron loss in the source PL

If

PL

- Pj-kPf = l-PT

-kPf (6)

where k is the number of neutrons per fission. The program

computes PT and Pf by integration and p_ and PA by P

T = 1- PT

and PA = Pj -Pf

.

Part III. The integration and normalization.

Consider the cylindrical source surrounded with a sphere

centered about the core material with radius, R, just large

enough to enclose the outer cylinder as shown in Figure 1(a).

The thermal-neutron flux is taken to be isotropic. Neutron

flux is the number of neutrons per second entering a sphere

of 1 square centimeter cross section. The number of neutrons

per second entering the sphere surrounding the cylindrical

source is the cross sectional area of the surrounding sphere

times the thermal flux.

Since the geometry is symmetric about the x-axis, the probability for various results for an isotropic flux will equal

the probabilities evaluated assuming a flux uniformly distributed in z, w, and cp but all paths constrained to lie in

planes parallel to the x-y plane. The neutron paths are perpendicular to all points along the w-axis within the sphere

and cp is the angle between the y and w axes. The probability

of transmission of a neutron through the imaginary sphere

averaged over all directions is thus

R 2* r(z)

.2

T =

J

dz [ dcp I* dw Pt

Pm (z,w,cp)/2tt(ttR ) (7) =

-R b -r(z)

The probability of transmission, Pt

(z,w,«p) is

-(2a

a(z,w,cp)+Sh

b(z,w,cp)+i: _c(z,w,cp)+Zd

d(z,w,cp))

PT

(z,w,cp)= e

b c a

(8)
where the £'s are the macroscopic cross sections. Parameters

are as in the spherical case and the a, b, c, and d are the

thicknesses of the A, B, C, and D materials as seen at the

particular z, w, and cp. For a transparent source the macroscopic cross sections are all zero and equation (7) integrates

to unity. The fission probability is similarly

R 2tt r(z)

Pp

* ^ dz ^

dcp ^ dw Pf (z,w,cp)/2rT(nR2

) (9)

-R o -r(z)

-P -Y

where Pf(z,w,cp) = (1-e )e (^ afAa

).

In P f

(z,w,cp).

p = Ea

a(z,w,cp) ,

-P

1-e is the probability for interaction in the "A"

material,

Y - 2b

b

1 (z,w,cp)+Sc c

i (z,w,cp)+E d

d

1 (z,w,cp)

where b^(z,w,<p) is the thickness of the "B" layer transversed passing into the "A" material and likewise for c^and d^;

-Y

e is the probability for passing through the cladding

material without interaction; a af is the fission cross section for the "A" material, and (<* af/a a ) is the probability

for fission, given that an interaction occurs in the "A"

material.

3. DESCRIPTION OF THE CODE

3.1 Spherical Case Calculation

3.1.1 Organization of the Code

The program first reads the title card and input data,

converts the source dimensions from inches to centimeters and

then prints out an echo check of the title, input data and

converted dimensions.

The approximate total number of subdivisions, NY, of

the radius is specified in the input. Then each layer is

subdivided in proportion to its thickness in the radial

direction. The boundary between the central core and
first encapsulation or the boundary between any of the

encapsulations are always made an endpoint of subdivision

so that any discontinuities occasioned by a change in

derivative at a boundary do not occur within an interval

of integration.

The integration of Equation (4) then begins with the

central core and continues to the outside radius, layer

by layer. All neutron paths at a given radius are perpendicular to a cross section through the center of the source

and are equal in length. Subroutine Grand calculates the

interaction probability in the A, B, C, and D materials and

the transmission probability along the path specified. The

probability of interaction in each material and the transmission probability with its error are summed and normalized

in the main program. The results are then printed below the

echo printout of the input data.

3.1.2 Approximations in the Numerical Integration

The numerical integration was performed using Weddle '

s

Rule (4). If the core material has been divided into N subdivisions, Weddle 's Rule replaces the integrand of the

integral by N fifth-degree parabolas and numerically integrates over that region. The accuracy of the integration

and estimate of the error depend upon the size of the integration intervals. Convergence of the integration procedure was checked for a plutonium-beryllium source with

encapsulations of tantalum and stainless steel by varying the

number of subdivisions, NY, of the radius between 10 and 1,000

In Table 1 the calculated transmissions and estimated errors

are shown together with the actual error based upon accepting

the finest subdivision as the correct answer. The estimated

and actual error appear to agree for 20 subdivisions of NY,

but this degree of accuracy is certainly not necessary. The

difference in net neutron loss calculated for NY equal 10

and 1000 is only 0.005$ out of 0.278$.

This program has been run in double precision, because

it was created, for the most part, from the deck of cards

to the cylindrical case, which had to be run in double precision. The number of additions involved in the transmission

integral is approximately equal to 6 times the number of

subdivisions, NY, of the radius specified in the input.

Therefore no rounding error should be introduced by performing the program in single precision. The program takes about

3 seconds to compile and about 2 seconds to execute in the

case of NY equal 100. It requires 1437 words of memory.
3.2 Cylindrical Case Calculation

3.2.1 Organization of the Code

Here too the title card and input data are read, the

source dimensions are converted to centimeters, and an echo

check of the input is printed.

According to the number, NZ, of subdivisions in the Zdirection, specified in the input, each layer in the Z-direction is subdivided in proportion to its thickness. The

subdivision of the angle of integration from zero to tt/2 is

performed according to the number, NPHI, specified in the

input. The subdivisions along the W-axis are set at the

interface of any two materials in a plane parallel to the

X-Y axis, distance Z above the origin. The D, C, and B

materials in the negative W-direction, the core A material,

and the B, C, and D materials in the positive W-direction

are all divided into six subintervals. The limits of integation for each of these subintervals along the W-axis are

passed to Subroutine IN. It checks the length of the interval and will further subdivide it by 6 or 36, if it is more

than 6 times the length of a subdivision in the Z -direction.

It sums and weights the contributions to the transmission,

the error for the transmission, and the interaction in the

core material for each path through the source perpendicular

to the W-axis between the specified limits of integration and

returns the result to the main program. These in turn are summed and weighted according to the Weddle Rule for the integration over angle, <p, and then over height, Z. The probabilities for transmission with its error, absorption, and

fission are then printed.

Subroutine C2 calculates the X and Y coordinates of the

points of incidence and the pathlength through any of the

cylindrical source encapsulations or core material and returns the result to Subroutine Grand. Subroutine Grand calculates the segments of the path in each of the source materials in order to compute the transmission and fission

probabilities along that path and returns the result to Subroutine PCE. Subroutine PCE subdivides by six the subinterval specified by Subroutine IN and performs the six point

Weddle integration of transmission and fission probabilities

with an estimate of the error for this segment of transmission. The seven values of the integrand are received from

Subroutine Grand and after proper normalization the results

are returned to Subroutine IN.

3.2.2 Approximations in the Numerical Integration

The comments about the use of the Weddle Rule in 3.1.2

also apply here. In this case we integrate over Z, <p, and W.

The size of the integration steps for Z and 9 are specified

in the input and, because the interval along the W-axis va-
ries with angle, <P , the thicknesses of the encapsulations,

and the inside or core length, the length of the W integration step is tested against the length of the Z step and

further subdivided, if necessary. The absolute values of

the estimated error for each segment of integration along

the W-axis are summed and added in quadrature with the estimated error in the <P integration. This error is then added

in quadrature with the estimated error in the Z integration.

It is also possible to test the relative error of any subinterval by testing against ET, which is specified in the

input, and to further subdivide, if necessary. It has been

found, however, that ET may safely be set as low as 1% or

0.5$ without causing many subdivisions, but setting ET as

low as o.l$ or 0.05$ will suddenly cause so many subdivisions that the program will be stopped automatically by

exceeding the maximum time specified on the run card and all

data are then lost. The transmission integral has also been

found to be insensitive to the value of ET. Each time an

interval in the W integration is subdivided once, twice or

three times a counter, KB, KC, or KD is incremented and

printed in the program output.

The convergence of the cylindrical program has been tested for two different plutonium-beryllium sources, one with

a single nickel encapsulation and the other with encapsulations of tantalum and stainless steel. In Figure 2 are shown

the calculated transmissions and estimated errors for the

nickel encapsulated source for various specified subdivisions

of the Z and <P integrals, NZ and NPHI, respectively. A solid

line joins those points for which NPHI is held constant at

10 and a dashed line joins those points for which NZ is held

constant at 10. Two points are also shown for both NPHI and

NZ equal to 5 and 15. Taking finer steps in the Z integral,

which forces finer steps in the W integral, increases the

transmission for this source, whereas finer steps in the <P integral decreases it. Taking equal steps in Z and <P appear

to result in a more rapid convergence, though this may be

fortuitous. The difference in net neutron loss in the source

for the two extreme points on the figure amounts to 0.01$,

which is of no consequence in such a calculation. The convergence of the program for the tantalum and stainless steel

encapsulated source was tested for NZ equal NPHI at values

of 3, 6, and 10. The transmissions were 0.1839, 0.1842 and

0.1843, respectively. The negligible difference in transmission corresponded to run times of 188, 516, and 1172 seconds, respectively.

This program was affected by rounding errors when run

in 'single precision by two different compilers. The results

were only independent of choice of compiler, when run in

double precision. The program takes about 5 seconds to compile and requires 2992 words of memory. Typical run times

are listed in the above paragraph.

4. LISTING OF THE PROGRAMS

4.1 SPHERICAL CASE

4.1.1 LIST OF THE PROGRAM

IMPLICIT DOUBLE PREC I S I ON ( A-H ,0-Z

)

1 DIMENSION V<6) .T< 14) .EH<6) , RZZ(5),KA(4) 2 DATA V/2. DO, 5. DO, 1. DO. 6. DO, 1. DO, 5. DO/, EH/ 1 . DO ,-6 . DO

,

3 ]15.D0,-2 0.D0.15.D0.-6.D0/,CC/l.D-8/ 4

1 FORMAT(13A6,A2/7G8.4,I4) 5

3 FORMATU2G6.3) 6 4 F0RMAT(1H1,13A6,A2/' SOURCE DlMENS IONS ' 25X ' / MACROSCOPIC CROSS SEC 7

1 '17X'/ INTEGRATION PARAMETERS'/' RA ='G14.6,' TB ='G14.6,' 8

1 / SIGAA='G14.6, ' AK ='G14.6.' / Ny ='16/' RB ='G14.6»« TC 9

1 ='G14.6,' / SIGAF='G14.6, « F ='G14.6» /' RC ='G1 10 14.6,' TD ='G14.6,' / SIGB ='G14.6» /' RD ='G 11 114.6.21X 1 / SIGC ='G14. 6.21X1 / NCY ='16/ 42X 12

1 ' / SIGD ='G14.6,21X' /' / 13

] 11X • IN INCHES, ' 2X • / COVER / WALL 14

1 THICKNESS '/13X'OUTSIDE DI AMETER= • G12 . 6

,

15

1' / B'4XG12.6» /43X'/ C4XG12.6. 16

1 /43X'/ D'4XG12.6) 17

5 FORMATI3X/ 46H FOR A SINGLE NEUTRON STRIKING THE SOURCE.../ 18 120X35HTHE TRANSMISSION PROBABILITY IS 22XG16.8/ 19 120X35HTHE INTERACTION PROBABILITY IS 22XG16.8/ 20

120X35HTHE ABSORPTION PROBABILITY IS 22XG16.8/ 21 120X44HTHE 'A' MATERIAL INTERACTION PROBABILITY IS 13XG16.8 / 22 120X44HTHE 'B' MATERIAL INTERACTION PROBABILITY IS 13XG16.8 / 23

120X44HTHE 'C MATERIAL INTERACTION PROBABILITY IS 13XG16.8 / 24

120X44HTHE 'D' MATERIAL INTERACTION PROBABILITY IS 13XG16.8 / 25

120X57HTHE PROBABILITY FOR A FISSION ABSORPTION IS G16 26

1.8/16X' NEUTRONS OUT PER NEUTRON STRIKE IS ' 22XG16.8 27

1/20X'NEUTRONS LOST PER NEUTRON STRIKE '22XG16.8,' .' 28

1//38H THE CROSS SECTION OF THE SOURCE IS G13.8.' CM**2.' 29

1/3X'THE ERROR IN THE TRANSMISSION PROBABILITY IS ESTIMATED TO BE 30

l'G12.3. 'PERCENT. ' ) 31

7 FORMAT) 3X/ 44H FOR A FLUX OF ONE NEUTRON PER CM**2-SEC * ( G 12'. 5 . 32

141H NEUTRONS PER SEC) STRIKING THE SOURCE.../ 33

120X26HTHE TRANSMISSION IS G16.8.20H NEUTRONS PER SECOND/ 34

120X26HTHE INTERACTION IS G16.8.20H NEUTRONS PER SECOND/ 35

120X26HTHE ABSORPTION IS G16.8.20H NEUTRONS PER SECOND/ 36

120X26HTHE FISSION EMISSION IS G16.8,20H NEUTRONS PER SECOND/ 37

116X26HAND THE RESULTANT LOSS IS 4XG16.8,' NEUTRONS PER SECOND.') 38

8 FORMATOX/ 35H IF THE MEASURED THERMAL FLUX IS G12.4, ' PER CE 39 INT OF Q PER CM**2 THEN...'/ 40

120X26HTHE TRANSMISSION IS G16.8.14H PER CENT OF Q/ 41 120X26HTHE INTERACTION IS G16.8.14H PER CENT OF Q/ 42

120X26HTHE ABSORPTION IS G16.8.14H PER CENT OF Q/ 43

120X26HTHE FISSION EMISSION IS G16.8.14H PER CENT OF Q/ 44

116X26HAND THE RESULTANT LOSS IS 4XG16.8,' PER CENT OF Q.') 45

9 FORMAT) lX'THE NUMBER OF GRAND CALLS WAS' 18.' THE LARGEST SUBDIVIS 46

HON OF THE RADIUS WAS 'G16.8) 47

100 READ(5.1 IT.SIGAA.SIGAF, SIGB, SIGC. SIGD. AK,F, NY 48

PEADI5.3) DIA.TBI ,TCI .TDI 'OUTSIDE D I AMETER+SHELL THICKNESSES 49

IGRAND=0 "COUNTER RESET FOR NUMBER OF SUB GRAND CALLS 50

NCY=4 'NCY EQUALS THE NUMBER OF LAYERS 51 IF(TDI.LT.CC) NCY=3 • 3 LAYERS 52 IF(TCI.LT.CC) NCY=2 • 2 LAYERS 53

IF(TBI.LT.CC) NCY=1 ' 1 LAYER 54

TB=TBI*2.54D0 'THE INPUT DIMENSIONS WERE IN INCHES 55 TC=TCI*2.54D0 'F=THE NEUTRON FLUX IN PERCENT OF Q 56

TD=TDI*2.54D0 'T=TITLE. AK=NEUTRONS PER FISSION 57 RD=DIA*2.54D0/2.D0 'OUTSIDE RADIUS 58 RC = RD-TD 'FOR ONE LAYER A, B.C. AND D DIMENS. ARE EQUAL 59

RB=RC-TC 'FOR TWO LAYERS B.C. AND D DIMENSIONS ARE EQUAL 60

RA=RB-TB 'FOR THREE LAYERS C AND D DIMENSIONS ARE EQUAL 61

DY=0.D0 'LARGEST DIV. OF RADIUS. Y. IS MAX VALUE OF DA 62

SIGAT=SIGAF+SI6AA 'SIGAT=TOTAL CROSS SECTION, 'A' MATERIAL 63
170

171

172

173

183

182

181

180

269

270

274

273

272

271

280

293

292

291

WRITE(6»4)T,RA»TB,SIGAA

1 NCY,SIGD,DIA,TBI ,TCI ,T RZZ< 1 )=0.D0

RZZ(2)=RA

RZZ(3)=RB

RZZ(4)=RC

RZZ(5 )=RD

EKI=O.DO

STI=O.DO

SYD=O.DO

SYC=O.DO

SYB=0.D0

SYA=0.D0

ANY=NY/RD

KA( 1)=RA*ANY+.5D0

KA(2)=(RB-RA)*ANY+.5D0

KA( 3)=(RC-RB)*ANY+.5D0

KA(4)=(RD-RC)*ANY+.5D0

DO 281 IAR=1,NCY

IF(KA( IAR).LT.l )KA( IAR)

DA=(RZZ( IAR+1)-RZZ( IAR)

IF(DA.GT.DY)DY=DA

NYP=6*KA( IAR)+1

IVY = 2 SPYT=0.D0

SPYD=0.D0

SPYC=0.D0

SPYB=0.D0

SPYA=0.D0

EK=0.D0

TEK=0.D0

DO 280 1=1, NYP Y=DA*DBLE ( I -1 ) +RZZ ( I AR PDY=0.D0

RCY=0.D0

RBY=0.D0

RAY=0.D0

GO TO (173,172,171,170)

IF(RD.GT,Y)RDY=DSQRT(RD

IF(RC.GT.Y) RCY=DSQRT(R

IF(RB.GT.Y) RBY=DSQRT(R

IF(RA.GT.Y) RAY=DSQRT(R

CALL GRAND( Y,SPIA ? SPIB»

IF( I.EQ.l.OR.I .EQ.NYP)

GO TO( 180,181,182,183)

,

SPYD=SPYD+V( IVY)*SPID

SPYC=SPYC+V( IVY)*SPIC

SPYB=SPYB+V( 1VY)*SPIB

SPYA=SPYA+V( IVY)*SPIA

SPYT=SPYT+V( IVY)*SPT

TEK=TEK+EH( IVY)*SPT SU

IVY=IVY+1

IF(IVY.GE.7) IVY=1

IF( IVY.EQ.2 )GO TO 269

GO TO 280

EK=EK+DABS(TEK)

TEK=SPT

GO TO 280

GO T01271, 272,273, 274)

i

SPYD=SPYD+SPID

SPYC=SPYC+SPIC

SPYB=SPYB+SPIB

SPYA=SPYA+SPIA

SPYT=SPYT+SPT

EK=EK+SPT

CONTINUE

EKI=EKI+EK*DA/140.D0

STI=.3D0*DA*SPYT+STI

GO TO(290, 291,292, 293)

,

SYD=SYD+0.3D0*DA*SPYD

SYC=SYC+0.3D0*DA*SPYC

SYB=SYB+0.3D0*DA*SPYB

,AK,NY,RB,TC,SIGAF,F,RC,TD,SIGB»RD,SIGC»

DI

•NY=APPROXIMATE NO.

'THIS IS TO SUBDIVI

•PROPORTION TO THE •FROM ZERO TO RA, F •TO RC, AND FROM RC

•ALWAYS OCCUR AT LA

=1 'DISCONTINUITIES

)/<6.D0*KA( IAR) ) •LARGEST DIVISION

'THE NUMBER OF PTS

•IVY= Y INTEGRATION

'Y PARTIAL SUM RESE

•INTERACTION IN 'D •INTERACTION IN 'C 'INTERACTION IN 'B •INTERACTION IN 'A 'ERROR RESET FOR Y •TEMPORARY EK •BEGIN Y INTEG

)

OF Y INTERVALS

DE THE Y-LAYERS IN THICKNESS OF EACH LAYER

ROM RA TO RB, FROM RB

TO RD. Y-INTEGRATION WILL

YER BOUNDARIES AND AVOID

F RADIUS, Y OF INTEGR. IN EACH LAYER

WEDDLE RULE COEF INDICATOR

T OF TRANSMISSION AND

• MATERIAL

• MATERIAL

• MATERIAL

• MATERIAL

INTEGRAL IN DO LOOP 280

R.FR O-RA-RB-RC-RD

IS

IS ONEONE

,NCY

**2-Y**2

)

C**2-Y**2) 'HIT C LAYER, IF THERE

B**2-Y**2) 'HIT B LAYER, IF THERE

A**2-Y**2) 'HIT A LAYER

SPICSPID.SPT)

GO TO 270

NCY

M ABSOLUTE VALUE OF ERROR EACH 6 SUBINTERVALS

•EK=ERROR IN Y INTEGRAL DUE TO WEDDLE RULE

NCYiCOEFF. OF FIRST AND LAST TERMS IS 1,

' Y NORMALIZATION, WEDDLE RULE

'TRANSMISSION INTEGRAL

NCY

64

65

66

67

68

69

70

71

7273

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

10
290 SYA=SYA+0.3D0*DA*SPYA

281 CONTINUE

VNORM=2.D0/RD**2

STIS=STI*VNORM

SAFIS=SYA*VNORM

SBFIS=SYB*VNORM

SCFIS=SYC*VNORM

SDFIS=SYD*VNORM

IF(SIGAT.LT.1.D-4)G0 TO

SFIS=SAFIS*SIGAF/SIGAT

GO TO 283

282 SFIS=0.D0

283 SNOS=SFIS*AK+STIS

ACROS=3.1416D0*RD**2

EM=EKI*1.D2/STI

FTSOU=STIS

FINSOU=l.D0-FTSOU

FMULSO=SFIS

FABSOU=FINSOU-FMULSO

FOUTSO=FTSOU+FMULSO*AK

FNOUT=l.D0-FOUTSO

FPT=STIS*ACROS

FPI=F.INSOU*ACROS

FPA=FABSOU*ACROS

FPF=SFIS*AK*ACROS

FNLS=FPI-FPF

FQT=F*FPT

FQI=F*FPI

FQA=F*FPA

FQF=F*FPF

FQL=F*FNLS

WRITE(6.9)IGRAND,DY

WRITE(6»5) FTS0U,FINSOU

1F0UTS0»FN0UT,ACR0S,EM

WRITE(6»7)ACROS,FPT.FPI

WRITE(6»8)F,FQT,FQI »FQA

GO TO 100

END OF Y INTEGRAL

VNORM=2/RD**2

THE INTEGRALS

INTERACTION IN "A" INTERACTION IN 'B' INTERACTION IN >C" INTERACTION IN «D» 282

FISSION PART OF «A" INTERACTION

AND UNIT NORMALIZES

MATERIAL

MATERIAL

MATERIAL

MATERIAL

NEUTRON

SOURCE

EM=ESTI

SINGLE

SINGLE

SINGLE

SINGLE

NEUTRON

RESULTA

UNIT FL UNIT FL UNIT FL UNIT FL UNIT FL PERCENT

PERCENT

PERCENT

PERCENT

PERCENT

S OUT PER NEUTRON HITTING SPHERE

CROSS SECTION

MATE OF PER CENT ERROR IN STI

STRIKE TRANSMISSION PROBABILITY

STRIKE INTERACTION PROBABILITY

STRIKE FISSION PROBABILITY

STRIKE ABSORPTION PROBABILITY

S OUT PER NEUTRON STRIKING SOURCE

NT LOSS

UX TRANSMISSION

UX INTERACTION

UX ABSORPTION

UX FISSION

UX. NET LOSS

Q OF TRANSMISSION

Q INTERACTION

Q ABSORPTION

Q FISSION GAIN

Q NET LOSS

,FABSOU,SAFIS,SBFIS,SCFIS,SDFIS,FMULSO,

,FPA,FPF,FNLS

,FQF,FQL

340

350

360

370

403

402

401

400

413

412

411

410

SUBROUTINE G IGRAND=IGRAN

DUMA=0.D0

DUMB=0.D0

DUMC=0.D0

DUMD=0.D0

A=0.D0

8=0. DO C=0.D0

D=0.D0

GO TO (370,3

D=(RDY-RCY)*

C=(RCY-RBY)*

B=(RBY-RAY)«

A=RAY*2.D0

DMD=-(SIGAT*

DMDEX=DEXP(D

DMDA=W*DMDEX

DEXD=DEXP(-S

DEXC=DEXP(-S

DEXB=DEXP(-S

DEXA=DEXP(-S

GO T0(400,40

DUMDA=W*(1-D

DUMCA=W*(DEX

DUMBA=W*(DEX

DUMA=W*(DEXD

GO T0(410,41

DUMD=DUMDA-D

DUMC=DUMCA-D

DUMB=DUMBA-D

RETURN

END

RAND (W, DUMA, DUMB, DUMC,DUMD,DMDA)

D+l 'COUNT OF SUBROUTINE CALL

•CONTAINS MULTIPLICATION PART OF INTEGRAND

60,350.340)

2. DO

2. DO

2. DO A+SIGB*B+SI

MD)

.NCY

•THE AMOUNT OF A.B.C, AND D MATERIAL

'TRAVERSED IS A,B,C, AND D, RESPECTIVELY.

GC*C+SIGD*D)

•TRANSMISSION PART OF INTEGRAND

IGD*D/

IGC*C/

IGB*B/

IGAT*A

1,402,

MDEX)

D-DMDE

D*DEXC

*DEXC*

1,412,

UMCA

UMBA

UMA

2. DO)

2. DO)

2. DO)

/2.D0

403) , X/DEX

-DMDE

DEXB413) ,

NCY

• IN

D)

X/( DMD

NCY

• IN

• IN

• IN

TERACTION IN A,B,C,AND D MATERIAL

•INTERACTION IN A,B, AND C MATERIAL

DEXD*DEXC) )• INTERACTION IN A AND B MATERIAL

EX/(DEXD*DEXC*DEXB) ) UNTERAC* IN A MAT.

TERACTION IN D MATERIAL

TERACTION IN C MATERIAL

TERACTION IN B MATERIAL

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

11
4.1.2 LIST OF THE INPUT DATA DECK

TITLE CARD. FORMAT t 13A6 »A2 ) EQUIVALENT SPHERICAL SOURCE FOR CYLINDRICAL SOURCE M-621 CARD 1 MACROSCOPIC CROSS SECTIONS AND INTEGRATION PARAMETERS. FORMAT ( 7G8 .4. 1 4

)

2 3456 789012 345678901234 567890123456789012 3456789012 3456789012345678901234567890

.865 4.305 1.16 .2811 .0 2.8

SIGAA SIGAF SIGB SIGC SIGD AK

I I I I

I I I I

I I I I

I I I I STRENGTH PER CM**2

I I I NUMBER OF NEUTRONS PER FISSION

I I ABSORPTION CROSS SECTION PER CM FOR D MATERIAL

I ABSORPTION CROSS SECTION PER CM FOR C MATERIAL

ABSORPTION CROSS SECTION PER CM FOR B MATERIAL

FISSION CROSS SECTION PER CM FOR A MATERIAL

ABSORPTION CROSS SECTION PER CM FOR A MATERIAL

.122 100 CARD 2 F NY

I 1/6 NUMBER OF INTEGRAI TION STEPS

FLUX IN PERCENT OF SOURCE

SOURCE DIMENSIONS IN INCHES. FORMAT t 12G6 . 3

)

12 3456 78 9012 345678901234 5678901234567890123456 7 8 9012 34567 89012 34567 8901234 567890

1.913 .11 .071 CARD 3 DIA TBI TCI TDI

I I I D SHELL THICKNESS IN INCHES

I I C SHELL THICKNESS IN INCHES

I B SHELL THICKNESS IN INCHES

OUTSIDE DIAMETER OF SOURCE IN INCHES

4.1.3 SIMPLIFIED SAMPLE OUTPUT

EQUIVALENT SPHERICAL S SOURCE DIMENSION INPUT

MACROSCOPIC CROSS SECT

SIGAA= 1.865 SIGAF=

NEUTRONS PER FISSION = FOR A SINGLE NEUTRON S THE TRANSMISSION PROBA

THE INTERACTION PROBAB

THE ABSORPTION PROBABI

THE 'A' MATERIAL INTER

THE <B' MATERIAL INTER

THE «C« MATERIAL INTER

THE PROBABILITY FOR A NEUTRONS OUT PER NEUTR

NEUTRONS LOST PER NEUT

THE CROSS SECTION OF T THE PERCENT ERROR IN T IF THE MEASURED THERMA

THE TRANSMISSION IN PE THE INTERACTION IN PER

THE ABSORPTION IN PERC

THE FISSION IN PERCENT

AND THE RESULTANT LOSS

OURCE FOR THE CYLINDRICAL PU-BE M-62

IN INCHES, OD=1.913» B SHELL=0.110»

IONS PER CM

4.305 » SIGB= 1.16

2.80 NY=100.

TRIKING THE SOURCE...

BILITY IS ILITY IS LITY IS ACTION PROBABILITY

ACTION PROBABILITY

ACTION PROBABILITY

FISSION ABSORPTION

ON STRIKE IS RON STRIKE IS HE SOURCE IN CM**2

HE TRANSMISSION PROBABILITY IS

L FLUX IN PERCENT OF Q PER CM**2 IS RCENT OF Q IS CENT OF Q IS ENT OF Q IS OF Q IS

IN PERCENT OF Q IS

1 SOURCE

C SHELL=0.071.

SIGC= 0.2811

IS

IS

IS

IS

IS

.13567724+0

.86432276+0

.59782176+0

.38195382+0

.36978466+0

.11258428+0

.26650101+0

.88188005+0

.11811995+0

.18543353+2

.0512

.122

.30694114+0

.19553480+1

.13524456+1

.16881265+1

.26722146+0

EQUIVALENT SPHERICAL SOURCE FOR THE CYLINDRICAL PU-BE A SOURCE

SOURCE DIMENSION INPUT IN INCHES, OD=1.180» B SHELL=0.044.

MACROSCOPIC CROSS SECTIONS PER CM

SIGAA= 1.214 , SIGAF= 2.802 . SIGB= 0.4106

NEUTRONS PER FISSION = 2.80 NY=100.

FOR A SINGLE NEUTRON STRIKING THE SOURCE...

THE TRANSMISSION PROBABILITY IS .26256084+0

THE INTERACTION PROBABILITY IS .73743916+0

THE ABSORPTION PROBABILITY IS .41973395+0

THE 'A' MATERIAL INTERACTION PROBABILITY IS .45535478+0

THE 'B' MATERIAL INTERACTION PROBABILITY IS .28208437+0

THE PROBABILITY FOR A FISSION ABSORPTION IS .31770520+0

NEUTRONS OUT PER NEUTRON STRIKE IS .11521354+1

NEUTRONS LOST PER NEUTRON STRIKE IS -.15213542+0

THE CROSS SECTION OF THE SOURCE IN CM**2 IS .70554114+1

12
THE PERCENT ERROR IN THE TRANSMISSION PROBABILITY IS .0185

IF THE MEASURED THERMAL FLUX IN PERCENT OF Q PER CM**2 IS .122 THE TRANSMISSION IN PERCENT OF Q IS .22600192+0

THE INTERACTION IN PERCENT OF Q IS .63475827+0

THE ABSORPTION IN PERCENT OF Q IS .36129028+0

THE r ISSION IN PERCENT OF Q IS .76571038+0

AND THE RESULTANT LOSS IN PERCENT OF Q IS -.13095211+0

EQUIVALENT SPHERICAL SOURCE FOR THE CYLINDRICAL AM-BE SOURCE

SOURCE DIMENSION INPUT IN INCHES. OD=1.551. B SHELL=0.080» C SHELL=0.042.

MACROSCOPIC CROSS SECTIONS PER CM

SIGAA= .1704 . SIGAF= 0.00082 . SIGB = 1.16 . SIGC= 0.281

NEUTRONS PER FISSION = 2.89 NY=100.

FOR A SINGLE NEUTRON STRIKING THE SOURCE...

THE TRANSMISSION PROBABILITY IS .34932104+0

THE INTERACTION PROBABILITY IS .65067896+0

THE ABSORPTION PROBABILITY IS .64994216+0

THE 'A' MATERIAL INTERACTION PROBABILITY IS ,15384751+0

THE 'B' MATERIAL INTERACTION PROBABILITY IS .41839021+0

THE «C« MATERIAL INTERACTION PROBABILITY IS .78441236^-1 THE PROBABILITY FOR A FISSION ABSORPTION IS .73680038-3

NEUTRONS OUT PER NEUTRON STRIKE IS .35145039+0

NEUTRONS LOST PER NEUTRON STRIKE IS .64854961+0

THE CROSS SECTION OF THE SOURCE IN CM#*2 IS .12189389+2

THE PERCENT ERROR IN THE TRANSMISSION PROBABILITY IS .0240

IF THE MEASURED THERMAL FLUX IN PERCENT OF Q PER CM**2 IS .121 THE TRANSMISSION IN PERCENT OF Q IS .51521920+0

THE INTERACTION IN PERCENT OF Q IS .95969683+0

THE ABSORPTION IN PERCENT OF Q IS .95861012+0

THE FISSION IN PERCENT OF Q IS .31406170-2

AND THE RESULTANT LOSS IN PERCENT OF Q IS .95655622+0

4.2 CYLINDRICAL CASE

4.2.1 LIST OF THE PROGRAM

IMPLICIT DOUBLE PREC I SI ON ( A-H ,0-Z ) 1 DIMENSION V(6) »T( 14) »EH<6> »AW(8) .RZZ(5) »KA(4) 2 ALIM(RZABCD.DLRDP) =RZABCD/CX+ < DLRDP-RZABCD*S/CX ) *S 3 DATA V/2.D0.5.D0.1.D0.6.D0.1.D0.5.D0/ » EH/ 1 .DO .-6 . DO

.

4 115.D0.-20.D0.15.D0.-6.D0/.CC/1.D-8/ 5

1 F0RMAT(13A6,A2/7G8.4,2I3,G8.4) 6

3 FORMATI12G6.3) 7 4 -0RMAT(1H1,13A6,A2/' SOURCE D IMENS I ONS • 25X • / MACROSCOPIC CROSS SEC 8

1 '17X'/ INTEGRATION PARAMETERS'/' RA ='G14.6,< TRB ='G14.6»' 9

1 / SIGAA='G14.6» • AK ='G14.6»' / NZ ='16/' RB ='G14.6»' TRC 10

1 ='G14.6»' / SIGAF='G14.6» • F ='G14.6»» / NPHI ='16/' RC ='G1 11 14. 6»' TRD ='G14.6»' / SIGB ='G14.6»21X< / ET ='G14.6/' RD = 'G 12 114.6.' TLB ='G14.6»' / SIGC =«G14.6,21X' / NCY .='16/' DLA ='G14 13 1.6.' TLC ='G14.6»' / SIGD ='G14.6»21X' /• /21X' TLD = 14 l'G14.6.' /•/• IN INCHES. OUTSIDE LENGTH ='G12.6.' / COVER / WALL 15

1 THICKNESS / LEFT END / RIGHT END ' / 12X ' OUTS I DE D I AMETER= ' G 12 . 6 . 16

1' / B'4XG12.6.4X'/'G12.6. ' / • Gl 2. 6/42X • / C • 4XG1 2 .6 . 4X ' / ' G12 .

6

17 1." /'G12.6/42X'/ D'4XG12.6,4X' /'G12.6. ' /'G12.6) 18

5 FORMAT13X/ 46H FOR A SINGLE NEUTRON STRIKING THE SOURCE.../ 19 120X35HTHE TRANSMISSION PROBABILITY IS 22XG16.8/ 20

120X35HTHE INTERACTION PROBABILITY IS 22XG16.8/ 21 120X35HTHE ABSORPTION PROBABILITY IS 22XG16.8/ 22 120X55HTHE 'A' MATERIAL ABSORPTION INTERACTION PROBABILITY IS 2XG16 23

1.8/20X52HTHE 'A' MATERIAL FISSION INTERACTION PROBABILITY IS 5XG16 24

1.8/20X44HTHE 'A' MATERIAL INTERACTION PROBABILITY IS 13XG16.8 / 25 120X49HTHE CLADDING MATERIAL INTERACTION PROBABILITY IS 8XG16.8 26

1 /16X' NEUTRONS OUT PER NEUTRON STRIKE IS ' 22XG16.8 27

1/20X'NEUTRONS LOST PER NEUTRON STRIKE '22XG16.8,' .' 28

1//46H THE AVERAGE CROSS SECTION OF THE SOURCE IS G13.8.' CM**2.' 29

1) 30

13
6 F0RMAT13X/ 62H FOR A SINGLE THERMAL NEUTRON STRIKING A SPHERE

10F RADIUS R=G13.8»' CM, WHICH JUST SURROUNDS THE SOURCE...'/

120X40HTHE PROBABILITY OF MISSING THE SOURCE IS17XG16.8/

120X'THE TRANSMISSION (INCLUDING MISS) PROBABILITY IS'9XG16.8/

120X35HTHE INTERACTION PROBABILITY IS 22XG16.8/

120X35HTHE ABSORPTION PROBABILITY IS 22XG16.8/

120X44HTHE 'A' MATERIAL INTERACTION PROBABILITY IS 13XG16.8 / 120X57HTHE PROBABILITY FOR A FISSION ABSORPTION IS G16

1.8/16X' NEUTRONS OUT PER NEUTRON INTO SPHERE ' 2 1XG1 6. 8 » 2H ./3XiT

1HE ERROR IN THE TRANSMISSION PROBABILITY IS ESTIMATED TO BE 'G12.3

l»i PER CENT. •

) 7 FORMAT13X/ 44H FOR A FLUX OF ONE NEUTRON PER CM**2-SEC . ( G12 . 5

151H NEUTRONS PER SEC) STRIKING THE THE ABOVE SPHERE.../

120X26HTHE TRANSMISSION IS G16.8.20H NEUTRONS PER SECOND/

120X26HTHE INTERACTION IS G16.8.20H NEUTRONS PER SECOND/

120X26HTHE ABSORPTION IS G16.8»20H NEUTRONS PER SECOND/

120X26HTHE FISSION EMISSION IS G16.8»20H NEUTRONS PER SECOND/

116X26HAND THE RESULTANT LOSS IS 4XG16.8,' NEUTRONS PER SECOND.')

8 FORMATI3X/ 35H IF THE MEASURED THERMAL FLUX IS G12.4, ' PER CE

INT OF PER CM**2 THEN...'/

120X26HTHE TRANSMISSION IS 120X26HTHE INTERACTION IS 120X26HTHE ABSORPTION IS 120X26HTHE FISSION EMISSION IS

G16.8,14H PER CENT OF Q/

G16.8»14H PER CENT OF Q/

G16.8»14H PER CENT OF Q/

G16.8»14H PER CENT OF Q/

116X26HAND THE RESULTANT LOSS IS 4XG16.8,' PER CENT OF Q.')

9 FORMAT! lX'THE NUMBER OF GRAND CALLS WAS'I8»'» 2ND 3RD, AND 4TH ORD

1ER COUNTS.KB, KC, AND KD WERE • I 8

,

• , ' I 8 , • , AND ' I 3 , / IX ' THE LARGEST R 1ELATIVE ERRORS, ERL1, ERL2, AND ERL3, ENCOUNTERED AFTER THE 1ST, 2 1ND» AND 3RD SUBDIVISIONS, THE LARGEST • / 1 X ' SUBDI V I S I ON , DZ , ALONG

] THE Z-AXIS, AND THE LARGEST 3RD ORDER SUBDIVISION, DWD. OF THE W1AXIS, WERE'/ 5<4X,G16.8)>

100 KB=0 'SECOND ORDER COUNT, STATEMENT 60+2, SUB IN KC=0 'THIRD ORDER COUNT, STATEMENT 65+2 SUB IN KD=0 'FOURTH ORDER COUNT, STATEMENT 121 SUB IN READ (5,1 )T,SIGAA,SIGAF,SIGB,SIGC,SIGD,AK,F,NPHI ,NZ,ET

IF(ET.LT..1D-2)ET=.1D-1ET=SPEC I F IED RELATIVE ERROR IN SUBROUTINE IN READ (5 ,3) AA,AB,TRBI »TLBI »TSB»TRCI ,TLCI ,TSC,TRDI «TLDI ,TSD

IGRAND=0 'COUNTER RESET FOR NUMBER OF SUB GRAND CALLS

NCY=4 'AA=OUTSIDE LENGTH IN INCHES

IFtTSD. LT.CC.AND.TRDI.LT. CC.AND.TLDI.LT. CO NCY=3 ' 3 LAYERS

IFtTSC.LT.CC.AND.TRCI.LT.CC.AND.TLCI.LT.ee) NCY=2 ' 2 LAYERS

IF( TSB.LT.CC.AND.TRBI.LT.CC.AND.TLBI.LT.ee) NCY=1 ' 1 LAYER

TRB=2.54D0*TRBI

TLB=2.54D0*TLBI

TRC=2.54D0*TRCI

TLC=2.54D0*TLCI

TRD=2.54D0*TRDI

TLD=2.54D0*TLDI

RD=AB*2.54D0/2.D0

RC=RD-2.54D0*TSD

RB=RC-2.54D0*TSC

RA=RB-2.54D0*TSB

DLA=2.54D0*AA-TRB-TRC-TRD-TLB-TLC-TLD

BLA=DLA/2.D0

DRB=BLA+TRB

DRC=DRB+TRC

DRD=DRC+TRD

DLB=BLA+TLB

DLC=DLB+TLC

DLD=DLC+TLD

PTR=DSQRT(DRD**2+RD**2)

RTL=DSQRT(DLD**2+RD**2)

R = RTL

IF(RTR.GT.RTL) R=RTR

ELM=R«0.001D0

ERL1=O.DO

ERL2=0.D0

ERL3=0.D0

DWD=0.D0

DZ=0.D0

SIGAT=SIGAF+SIGAA

RK=DLD

THE INPUT DIMENSIONS WERE IN INCHES

F=THE NEUTRON FLUX IN PERCENT OF Q

AK=NEUTRONS PER FISSION

NCY=THE NUMBER OF DIFFERENT MATERIALS

T=TITLE

TRB= THICKNESS, RIGHT END, B CYLINDER

OUTSIDE RADIUS

FOR ONE LAYER A,B,C, AND D DIMENS. ARE EQUAL

FOR TWO LAYERS B,C, AND D DIMENSIONS ARE EQU

FOR THREE LAYERS C AND D DIMENSIONS ARE EQUA

ONE-HALF OF INSIDE LENGTH

AB=OUTSIDE DIAMETER OF SOURCE IN INCHES

SIGAA=MACROSCOPIC ABSORBTION CROSS SECTION OF MATERIAL A IN 'PER CM. < . SIGAF= FISSION IN 'A' IN PER CM.

SIGB=ABSORPTION IN 'B' TLD=THICKNESS, LEFT END, CYLINDER D RA=RADIUS OF CYLINDER A

R=SURROUNDING SPHERE RADIUS

SEE ST 210+4. SKIPS W INTEGR.IF LIM.TOO SMALl

LARGEST ERROR IN THE 1ST, 2ND, AND 3RD

ORDER WEDDLE DIVISIONS IN SUBROUTINE IN. LARGEST DIVISION OF W-AXIS, SUBROUTINE IN LARGEST DIV. OF Z-AXIS IS MAX VALUE OF DA SIGAT=TOTAL CROSS SECTION, 'A' MATERIAL

NZ=NO. OF DIVS. OF Z FOR WEDDLE INTEGRATION

31

3233

34

35

36

37

38

39

40

4142

43

44

45

46

47

48

49

50

5152

53

54

55

56

57

5859

60

61

62

63

64

65

66

67

68

69

70

7172

73

74

7576

77

78

79

80 AL 81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

14
IF(DRD.GT.DLD)RK=DRD 'RK=LARGER OF DRD AND DLD 103

U=DATAN2(RK,RD) • ARCTAN < RK/RD

)

104

AG=U*(R**2) 105

ARE=RD*RK 106

AD=AG-ARE 'AD = CROSS SECTIONAL AREA OF SPHERE ABOVE 107

NPHIP=6*NPHI+1 'CYLINDER AND NOT IN THE Z INTEGRAL RANGE. 108

VKB=6.D0*NPHI 'VKB = THE NUMBER OF PHI INTERVALS PER 90 DEGR 109

DB=1.5707963D0/VKB 'DB= PHI INTERVAL IN RADIANS 110

WRITE(6»4)T,RA,TRB,SIGAA,AK,NZ,RB,TRC,SIGAF»F,NPHI ,RC , TRD , S I GB , ET 111

1 »RD»TLB,SIGC,NCY,DLA,TLC,SIGD,TLD,AA,AB,TSB»TLBI »TRBI »TSC , TLCI 112

] TRCI ,TSD,TLDI »TRDI 113

RZZ(1)=0.D0 114

RZZ(2)=RA 115

RZZ(3)=RB 116

RZZ(4)=RC 117

RZZ(5)=RD 118

EJI=0.D0 119

EKI=0.D0 120

STI=0.D0 121

SFI=0.D0 122

SMI=0.D0 123

ANZ=NZ/RD 'NZ=APPROXIMATE NO. OF Z INTERVALS 124

KA< 1>=RA*ANZ+.5D0 'THIS IS TO SUBDIVIDE THE Z-LAYERS IN 125

KA(2)=(RB-RA)*ANZ+.5D0 'PROPORTION TO THE THICKNESS OF EACH LAYER 126

KA(3)=(RC-RB)*ANZ+.5DO 'FROM ZERO TO RA , FROM RA TO RB . FROM RB 127

KA(4) = (RD-RC)*ANZ+.5D0 'TO RC, AND FROM RC TO RD. Z-I NTEGRAT I ON WILL 128

DO 281 IAR=1.NCY 'ALWAYS OCCUR AT LAYER BOUNDARIES AND AVOID 129

IF(KA( IAR) ,LT.1)KA< IAR)=1 'DISCONTINUITIES. 130

IF( (RZZ( IAR+1)-RZZ( IAR) ) .LT.ELM)GO TO 281'IN CASE OF ZERO SIDEWALL 131

DA=(RZZ( IAR+1)-RZZ( IAR) )/(6.D0*KA( IAR) ) 132

IF(DA.GT.DZ)DZ=DA 'LARGEST DIVISION OF Z-AXIS 133

NZP=6*KA( IAR1+1 'THE NUMBER OF Z STEPS 134

IVZ=2 'IVZ= Z INTEGRATION WEDDLE RULE COEF INDICATOR 135

SPZT=0.D0 'Z PARTIAL SUM RESET OF TRANSMISSION, 136

SPZM=0.D0 'MULTIPLICATION, AND 137

SPZMIS=0.D0 'MISSES 138

EJ=0.D0 'ERROR RESET FOR Z INTEGRAND IN DO LOOP 280 139

EK=0.D0 'ERROR RESET FOR Z INTEGRAL IN DO LOOP 280 140

TEK=0.D0 'TEMPORARY EK 141

DO 280 1=1, NZP 'BEGIN. Z INTEGR.FROM O-RA-RB-RC-RD 142

IVP=2 'ANGLE INTEGRATION WEDDLE RULE COEFFICIENT INDICATOR. 143

SPPT = 0.D0 'ANGLE PARTIAL SUM RE-SET OF TRANSMISSION, 144

SPPM=0.D0 'MULTIPLICATION, AND 145

SPPMIS=0.D0 'MISSES 146

EF=0.D0 'ERROR RESET FOR PHI INTEGRAND IN DO LOOP 250 147

EG=0.D0 'ERROR RESET FOR PHI INTEGRAL IN DO LOOP 250 148

TEG=0.D0 'TEMPORARY EG 149

Z=DA*DBLE ( 1-1 )+RZZ( IAR) 150

RZ=DSQRT(R**2-Z**2) 151

RZD=0.D0 152

RZC=0.D0 153

RZB=0.D0 154

RZA=0.D0 155

GO TO ( 173,172,171 ,170) ,NCY 156

170 IF(RD.GT.Z)RZD=DSQRT(RD**2-Z**2) 157

171 IF(RC.GT.Z) RZC=DSQRT(RC**2-Z**2) 'HIT C LAYER, IF THERE IS ONE 158

172 IF(RB.GT.Z) RZB=DSQRT < RB**2~Z**2 ) 'HIT B LAYER, IF THERE IS ONE 159

173 IF(RA.GT.Z) RZA = DSQRT ( RA**2-Z**2 ) 'HIT A LAYER 160

DO 250 J=1,NPHIP 'BEGINNING OF PHI INTEGRATION FOR PHI FROM 161

PHI=DB*DBLE( J-l) '0 TO 90 DEGREES ONLY. PHI IS THE ANGLE 162

S=DSIN(PHI) 'BETWEEN THE Y-AXIS AND W-AXIS 163

CX=DCOS(PHI) 164

SPWT=0.D0 'W PARTIAL SUM RESET OF TRANSMISSION, 165

SPWM=0.D0 'AND MULTIPLICATION 166

IFLAG=3 'PHI NOT EQUAL TO OR 90 DEGREES 167

IF(J.EQ.l) GO TO 180 168

IF( J.EQ.NPHIP) GO TO 185 169

GO TO 200 170

180 IFLAG=1 'PHI=0 DEGREES 171

GO TO( 184,183,182,181) ,NCY 172

181 AW(1)=-RZD 'W(l) IS THE POSITION ALONG THE W-AXIS THAT A 173

AW(8)=RZD 'NEUTRON, INCIDENT PERPENDICULARLY FROM THE 174

15
182 AW(2)=-RZC 'RIGHT OR LEFT SIDE OF THE W-AXIS WILL JUST 175

AW<7)=RZC 'STRIKE THE LEFT OUTER CORNER OF THE D 176

183 AW(3)=-RZB 'MATERIAL. AW<2), AW(3) » AND AW<4) WILL JUST 177

AW(6)=RZB 'STRIKE THE C,B, AND A LAYERS* RESPECTIVELY. 178

184 AW(4)=-RZA 'AW(5)» AW(6)» AW ( 7 ) . AND AW(8) WILL JUST 179

AW(5)=RZA 'STRIKE THE UPPER RIGHT OUTER CORNER OF THE 180

GO TO 210 'A. B, C. AND D MATERIALS, RESPECTIVELY. THIS181

185 IFLAG=2'PHI=90 DEGREES. 'AVOIDS DISCONTINUITIES IN THE DERIVATIVES FOR182

GO TO( 189.188, 187*186) .NCY'THE WEDDLE RULE. 183

186 AW( 1)=-DLD 184

AW(8)=DRD 185

187 AW(2)=-DLC 186

AW(7)=DRC 187

188 AW(3)=-DLB 188

AW(6)=DRB 189

189 AW(4)=-BLA 190

AW(5)=BLA 191

GO TO 210 192

200 GO TO(204,203,202,201) ,NCY 193

201 AW( 1)=-ALIM(RZD»DLD) 194

AW(8)=ALIM(RZD.DRD) 195

202 AW(2)=-ALIM(RZC.DLC) 196

AW(7)=ALIM(RZC»DRC) 197

203 AW(3)=-ALIM(RZB,DLB) 198

AW(6)=ALIM(RZB»DRB) 199

204 AW(4)=-ALIM(RZA,BLA) 200

AW(5)=-AW<4) 201

210 LL=5-NCY 202

LU=3+NCY 203

DELW=AW(LU+1 )-AW(LL) 'RQ=2* THE PORTION OF THE SPHERE CUT BEYOND 204

RQ=2.D0*(2.D0*RZ-DELW) 'CYLINDER BECAUSE OF INCIDENCE FROM TWO 205

IFIDELW.LT. ELM)GO TO 231 'DIRECTIONS 206

EZ=0.D0 207

DO 230 I2=LL»LU 208

CALL IN( AW( 12) ,AW( 12+1) ,E,SPT,SPM) 209

EZ=EZ+E 210

SPWT=SPWT+SPT 211

SPWM=SPWM+SPM 212

230 CONTINUE 'END OF W INTEGRAL DO LOOP, WHICH BECOMES 213

231 GO T0(240, 240, 235) ,IFLAG' INTEGRAND FOR PHI INTGR. WITH ABS. ERROR EZ 214

235 SPPT=SPPT+V( IVP)*SPWT 215

SPPM=SPPM+V( IVP)*SPWM 216

SPPMIS=SPPMIS+V( IVP)*RQ 'MISS CYLINDER CALCULATION 217

EF=EF+(V( IVP)* EZ)**2" 'ERROR**2 IN PHI IN. DUE TO INTEGRND. UNCRT. 218

TEG=TEG+EH( I VP ) *SPWT ' SUM ABSOLUTE VALUE OF ERROR EACH 6 SUBDIVISIONS 219

IVP=IVP+1 »EH( IVP)=WEDDLE RULE ERROR CONSTANTS 220

IFUVP.GE.7) IVP=1 221

IF( IVP.EQ.2JGO TO 239 222

GO TO 250 223

239 EG=EG+DABS(TEG) 'EG=ERROR IN PHI INTEGRAL DUE TO WEDDLE RULE 224

TEG=SPWT 225

GO TO 250 226

240 SPPT=SPPT+SPWT 'COEFF. FOR FIRST AND LAST TERMS IS ONE 227

SPPM=SPPM+SPWM 228

SPPMIS=SPPMIS+RQ 229

EG=EG+SPWT 230

EF=EF+EZ**2 231

250 CONTINUE 'END OF PHI INTEGRAL DO LOOP, WHICH BECOMES 232

SPPT=SPPT*0.3D0*DB 'INTEGRAND FOR THE Z INTEGRAL 233

SPPM=SPPM*0.3D0*DB 'THIS IS THE WEDDLE RULE NORMALIZATION 234

SPPMIS=SPPMIS*0.3D0*DB 'FOR THE PHI PART 235

EF=.3D0*DB*DSQRT(EF> 236

EG=EG*DB/140.D0 237

EI=DSQRT<EG**2+EF**2

)

'EI=TOTAL ABS. ERROR ESTIM. FOR PHI INTEGRAL 238

IF( I.EQ.l.OR.I.EQ.NZP) GO TO 270 239

SPZT=SPZT+V( IVZ)*SPPT 240

SPZM=SPZM+V< IVZ)*SPPM 241

SPZMIS=SPZMIS+SPPMIS*V< IVZ) 242

EJ=EJ+(V( IVZ)*EI )**2 'ERROR**2 DUE TO INTEGRAND UNCERT. IN Z INTEG. 243

TEK=TEK+EH( I VZ ) *SPPT ' SUM ABSOLUTE VALUE OF ERROR EACH 6 SUBDIVISIONS 244

IVZ=IVZ+1 245

IFUVZ.GE.7) IVZ=1 246

16
IF( IVZ.EQ.2)G0 TO 269

GO TO 280

269 EK=EK+DABS(TEK>

TEK=SPPT

GO TO 280

270 SPZT=SPZT+SPPT

SPZM=SPZM+SPPM

SPZMIS=SPZMIS+SPPMIS

EJ=EJ+EI**2

EK=EK+SPPT

280 CONTINUE

EJI=EJI+.3D0*DA*DSQRT(EJ)

EKI=EKI+EK*DA/140.D0

STI=.3D0*DA*SPZT+STI

SFI=.3D0*DA*SPZM+SFI

SMI=.3D0*DA#SPZMIS+SMI

281 CONTINUE

AR =3.1416D0*(R#*2)

VNORM=2.DO/(3.1416D0*AR)

EL=DSQRT(EJI**2+EKI**2)

STIS=STI*VNORM

SAFIS=SFI*VNORM 'INTERACTION

IF(SIGAT.LT.1.D-4)G0 TO 282

SFIS=SAFIS*SIGAF/SIGAT 'FISSION PART OF 'A' INTERACTION

GO TO 283

282 SFIS=0.D0

283 SMIS=SMI*VNORM

SLM=2.D0*AD/AR 'UNIT NORMALIZED. LUMP

SNOS=SFIS*AK+SLM+SMIS+STIS 'NEUTRONS OUT PER

'EK=ERROR IN Z INTEGRAL DUE TO WEDDLE RULE

•COEFF. OF FIRST AND LAST TERMS IS 1,

Z NORMALIZATION. WEDDLE RULE

•TRANSMISSION INTEGRAL

'FISSION INTEGRAL

'MISSES INTEGRAL

'END OF Z INTEGRAL

'SPHERE CROSS SECTION

'UNIT NORMALIZES THE INTEGRALS

IN 'A' MATERIAL

MISS PART

NEUTRON HITTING SPHERE

FMIS=SLM+SMIS

ACR0S=AR*(1.D0-FMIS)

FTRANS=SLM+STIS+SMIS

EM= EL* l.D2*VNORM/F TRANS

FINT=1.D0-FTRANS

FMUL=SFIS

FABS=FINT-FMUL

FNO=FMUL*AK+FTRANS

SINGLE NEUTRON

AVERAGE SOURCE

SINGLE NEUTRON

EM=ESTIMATE OF SINGLE NEUTRON

SINGLE NEUTRON

SINGLE NEUTRON

MISS PROBABILITY

CROSS SECTION

TRANSMISSION PROBABILITY

PER CENT ERROR IN FTRANS

INTERACTION PROBABILITY

FISSION PROBABILITY

ABSORPTION PROBABILITY

NEUTRONS OUT PER NEUTRON STRIKING SPHERE

SINGLE STRIKE INTERACTION PROBABILITY

'SINGLE STRIKE 'A' INTERACTION PROBABILITY

SINGLE STRIKE FISSION PROBABILITY

SINGLE STRIKE 'A' ABSORPTION INTERACT. PROBAB.288

SINGLE STRIKE CLADDING INTERACTION PROBABIL. 289

SINGLE STRIKE ABSORPTION PROBABILITY

NEUTRONS OUT PER RESULTANT LOSS

UNIT FLUX

FLUX

FLUX

FLUX

FLUX.

NEUTRON STRIKING SOURCE

FTSOU=(FTRANS-FMIS)*AR/ACROS 'SINGLE STRIKE TRANSMISSION PROBABILITY

FINSOU=l.D0-FTSOU

FSAFIS=SAFIS/(1.D0-FMIS

FMULSO=FMUL/ < 1.D0-FMIS)

AAIP=FSAFIS-FMULSO

BCDIP=FINSOU-FSAFIS

FABSOU=FINSOU-FMULSO

FOUTSO=FTSOU+FMULSO*AK

FNOUT=l.D0-FOUTSO

FPT=FTRANS*AR

FPI=FINT*AR

FPA=FABS*AR

FPF=FMUL*AR*AK

FNLS=FPI-FPF

FQT=F*FPT

FQI=F*FPI

FQA=F*FPA

FQF=F*FPF

FQL=F*FNLS

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

TRANSMISSION

INTERACTION

ABSORPTION

FISSION

NET LOSS

OF TRANSMISSION

INTERACTION

ABSORPTION

FISSION GAIN

NET LOSS

UNIT

'UNIT UNIT

'UNIT PERCENT

'PERCENT

'PERCENT

PERCENT

'PERCENT

WRITE (6 .9) IGRAND,KB.KC,KD.ERL1,ERL2.ERL3.DZ.DWD

WRITE (6. 5) FTSOU.FINSOU.FABSOU.AAIP.FMULSO.FSAFIS.BCDIP.FOUTSO.

1FNOUT.ACROS

WRITE(6.6)R,FMIS»FTRANS»FINT.FABS.SAFIS.FMUL.FNO.EM

WRITE(6.7)AR.FPT,FPI .FPA .FPF.FNLS

WRITE(6.8)F»FQT»FQI .FQA.FQF.FQL

GO TO 100

SUBROUTINE I N ( FLL , FUL »E . ANS.QE ) 'TO GIVE THE WEDDLE RULE- INTEGRAL TO

DW1=(FUL-FLL)/36.D0 'A SPECIFIED RELATIVE ERROR ET. DW2=DW1/6.D0 'IF THE SUBDIVISION ALONG THE W-AXIS IS 6 IF(DW2.GT.DA)GO TO 65 'TIMES LARGER THAN IN THE Z-DIRECTION, THEN

IF(DW1.GT.DA)G0 TO 60 'IT IS AUTOMATICALLY DIVIDED DOWN BY 6 OR 36

CALL PCE(FLL.FUL.ANS.AE.QE)

IF(ANS.LT.1.D-4)G0 TO 50

EE=DABS< AE/ANS)

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

17
IF(EE.GT.ERL1)ERL1=EE 'ERL1 IS LARGEST ERROR IN 1ST ORDER WEDDLE DIV.319

IF(EE.GT.ET)GO TO 60 'FLL=LOWER LIMIT 320

E=AE «FUL=UPPER LIMIT 321

RETURN • ANS=INTEGRAL 322

50 E=0.D0 323

RETURN 324

60 ADIV=6.D0 325

IA1=6 326

KB=KB+1 'RESET COUNTER IN STATEMENT 100 327

GO TO 70 328

65 ADIV=36.D0 329

IA1=36 330

KC=KC+1 '3RD ORDER COUNT 331

DW3=DW2/6.D0 332

IF(DW3.GT.DWD)DWD=DW3 'DWD IS LARGEST 3RD ORDER SUBDIV. OF W-AXIS 333

70 ED=(FUL-FLL) /ADIV ' PCE RETURNS ABSOLUTE ERROR. 334

P1I=0.D0 335

OE=0.D0 336

E=0.D0 337

FL2=FLL 338

FU2=FLL+ED 339

DO 120 IA=1»IA1 340

CALL PCE(FL2»FU2.PI tPEI »QF) "2ND ORDER OR 3RD ORDER 341

P1I=P1I+PI 342

QE=QE+QF 343

E=E+PEI 344

FL2=FU2 'INCREMENT LIMITS 345

120 FU2=FU2+ED 346

ANS=P1I 347

IF( ANS.LT.1.D-41G0 TO 122 348

EE=DABS( E/ANS) 349

IF( IA1.EQ.36JGO TO 121 350

IF(EE.GT.ERL2)ERL2=EE • ERL2 IS LARGEST ERROR IN 2ND ORDER WEDDLE DIV.351

IF(EE.GT.ET)GO TO 65 '4TH ORDER COUNT ,KD , USED HERE TO COUNT NO OF 352

RETURN 'TIMES 3RD ORDER WEDDLE ERROR FAILS TO 353

121 IF(EE.GT.ET)KD=KD+1 'SATISFY RELATIVE ERROR ET IN INPUT 354

IF(EE.GT.ERL3)ERL3=EE 'ERL3 IS LARGEST ERROR IN 3RD ORDER WEDDLE DIV.355

RETURN 356

122 E=0.DO 357

RETURN 358

359

SUBROUTINE PCE ( FL , FU , AI AE ,QC ) » GI VES THE WEDDLE CORRECT SUM BETWEEN FL360

DIMENSION VI I (7) ,VIE(7) 'AND FU USING 6 SEGMENTS. 361

DATA VII/l.D0»5.D0,l.D0»6.D0»l.D0»5.D0»l.D0/.VIE/l.D0»-6.D0»15.D0f 362

1-2G.D0»15.D0,-6.D0,1.D0/,CA/3.D-1/, CB/ . 71428571D-2/ 363

EC=(FU-FL)/6.D0 'AI IS THE FIRST INTEGRAL. TRANSMISSION PART 364

PI=0.D0 'QC IS THE SECOND INTEGRAL, MULTIPLICATION PART365

QB=0.D0 'AE IS THE ABSOLUTE ERROR IN THE FIRST INTEGRAL366

PE=0.D0 'EC IS THE INCREMENT SIZE 367

X=FL 'X IS THE INTEGRAND EVALUATION POINT 368

DO 100 ID=1,7 369

CALL GRAND(X,3A,QG) 'QA IS THE MULTIPLICATION PART. INTEGRAND 370

PI=PI+QG*VII (ID) 'QG IS THE TRANSMISSION PART, INTEGRAND 371

QB=QB+VI I ( ID)*QA 'PI IS THE TRANSMISSION PART PARTIAL SUM 372

PE=PE+QG*VIE< ID) 'QB IS THE MULTIPLICATION PARTIAL SUM 373

100 X=X+EC 374

AI=CA#EC*PI 'PE IS THE ERROR PARTIAL SUM (6TH DIFFERENCE) 375

AE=DABS(CB*EC*PE) 'QC IS THE MULTIPLICATION PART ANSWER 376

QC=CA*EC*QB 377

RETURN 378

379

SUBROUTINE GRANDt W.DUMC »DUMD

)

380

C FINAL EVALUATION OF THE INTEGRAND, WHICH IS PUT INTO PCE 381

C CONTRIBUTIONS ARE ADDED FOR NEUTRONS INCIDENT ON SOURCE FROM OPPOSITE 382

C DIRECTIONS, SO THAT THE PHI INTEGRAL IS CARRIED OUT FROM TO PI/2 383

C INSTEAD OF FROM TO PI. DUMD = TRANSM I SS I ON PART OF INTEGRAND. 384

IGRAND=IGRAND+1 385

DUMC=0.D0 'CONTAINS MULTIPLICATION PART OF INTEGRAND 386

A=0.D0 387

B=0.D0 388

C=0.D0 389

D=0.D0 390

18
GO TO(90, 60, 30,10) , NCY 'NCY=NUMBER OF LAYERS 391

10 CALL C2( IB,DRD»DLD,RZD,W,XID,YID,D> 392

GO TO (20,30) , IB 393

20 DUMD=2.D0 "TRANSMISSION. IT MISSES IN BOTH DIRECTIONS 394

RETURN 395

30 CALL C2< IB,DRC»DLC,RZC»W,XIC»YIC»C) 396

IFt IB.EQ.2) GO TO 60 397

IF(NCY.EQ.3) GO TO 20 398

GO TO 130 399

60 CALL C2

(

IB,DRB,DLB,RZB,W,XIB,YIB,B) 400

IF( IB.EQ.2) GO TO 90 401

IFtNCY.EQ.2) GO TO 20 402

GO TO 130 403

90 CALL C2( IB,BLA»BLA»RZA»W,XIA,YIA,A) 404

IFIIB.EQ.2) GO TO 130 405

IF(NCY.EQ.l) GO TO 20 406

130 GO TO ( 170,160.150,140) .NCY 407

140 D=D-C 'THE AMOUNT OF A,B,C» AND D MATERIAL 408

150 C=C-B 'TRAVERSED IS A,B,C, AND D, RESPECTIVELY. 409

160 B=B-A 410

170 DUM=-(SIGAT*A+SIGB*B+SIGC*C+SIGD*D) 411

DUMD=2.D0*DEXP(DUM) 'TRANSMISSION 412

IF(IB.EQ.l) RETURN 'DID NOT HIT A MATERIAL. NO FISSION 413

DIDI=0.D0 414

DIDO=0.D0 . 415

DICI=0.D0 416

DICO=0.D0 417

DBPBA=-SIGAT*A 418

DBPBE=(1.D0-DEXP(DBPBA)

)

419

GO TO (220,210,200,190) , NCY 420

190 DIDI=DSQRT< ( X I D-X I C ) **2+ ( Y I D-Y I C ) **2 ) 'ENTRY THICK. FOR FISSION CALC. 421

DIDO=D-DIDI 'SAME, FROM OTHER SIDE 422

200 DICI=DSQRT( ( X I C-X I B > **2+ ( YIC-YI3)**2) 423

DICO=C-DICI 424

210 DIBI=DSQRT( ( X I B-X I A ) **2+ ( Y I B- Y I A ) **2

>

425

DIBO=B-DIBI 426

DUMBI=-(SIGD*DIDI+SIGC*DICI+SIGB*DIBI ) 427

DUMBO=-(SIGD*DIDO+SIGC*DICO+SIGB*DIBO) 42 8 DUMCI=DEXP(DUMBI ) 'FISSION EFFECT 429

DUMCO=DEXP(DUMBO) 'FROM OPPOSITE SIDE 430

DUMC=( DUMCI+DUMCO)*DBPBE 431

RETURN 432

220 DUMC=2.D0*DBPBE 433

RETURN 434

435

SUBROUTINE C2 ( I B, DR , DL , C2R , C2W , X I , Y I , DTHRU

)

436

YPF(X)=C2W*CX-S*(X-C2W*S )/CX • Y-COORD I NATE , GIVEN X, PHI, AND C2W(=W)437

XPF( Y)=C2W*S + CX*(C2W*CX-Y) /S ' X-COORD I NATE , GIVEN Y, PHI, AND C2W(=W)438

IF( IFLAG.EQ.DGO TO 101 439

IF( IFLAG.EQ.21GO TO 102 440

TR=YPF(DR) 441

l"L = YPF(-DL) 442

DTHRU=0.D0 443

IF1C2R.GT.TR) GO TO 30 444

20 IB=1 'NEUTRON MISSES CYLINDER. 445

RETURN 446

30 IF(-C2R.GE.TL) GO TO 20 447

IB=2 'NEUTRON HITS CYLINDER. 448

IFI-C2R.GT.TR) GO TO 60 449

XI=DR 'NEUTRON HITS END 450

YI=TR 451

GO TO 70 452

60 XI=XPF(-C2R) 'NEUTRON HITS BOTTOM. 453

YI=-C2R 454

70 IFIC2R.GT.TL) GO TO 90 455

XO=XPF(C2R) 'NEUTRON EXITS FROM THE TOP. 456

YO=C2R 457

GO TO 100 458

90 XO=-DL 'NEUTRON EXITS FROM THE END 459

YO=TL 460

100 DTHRU=DSQRT( (XI-XO)**2+( YI-YO)**2) 461

RETURN 462

19
101 DTHRU=DR+DL

XI=DR

YI=C2W

IB = 2 RETURN

102 DTHRU=2.D0*C2R

XI=C2W

YI=-C2R

IB = 2 RETURN

END

463

464

465

466

467

468

469

470

471

472

473

4»2.2 LIST OF THE INPUT DATA DECK

TITLE CARD. FORMAT ( 13A6 , A2

)

SOURCE M-621 PU-BE 80 GRAM, TA AND STAINLESS STEEL ENCAPSUL

MACROSCOPIC CROSS SECTIONS AND INTEGRATION PARAMETERS. FORMA

234 567 89012345678901234567890123456789012 34 567890123 45678901

.865

SIGAA

4.305

SIGAF

1. 16 SIG3.2811

SIGC

.0SIGD

2.{ AK

I

I

I

I

I

I

I

I

I

I

.122

F 10 1 N- NZ

PHI I

I II

I1/

GR

1/6 N GRATI

FLUX IN PERCE

STRENGTH PER NUMBER OF NEUTRONS PE ABSORPTION CROSS SECTION PER

ABSORPTION CROSS SECTION PER CM FOR C ABSORPTION CROSS SECTION PER CM FOR B MATERIA

FISSION CROSS SECTION PER CM FOR A MATERIAL

ABSORPTION CROSS SECTION PER CM FOR A MATERIAL

ATED

T(7G8.4»2

234567890

0.01

ET

IMINIMUM

AL ERROR

IN SUBRO

6 NUMBER

ATION STE

UMBER OF

ON STEPS

NT OF SOU

CM**2

R FISSION

CM FOR D MATERIAL

L

CARD 1 I3.G8.4)

1234567890

CARD 2

FRACTIONALLOWED

UTINE PCE

OF Z INTEPS

PHI INTERCE

MATERIAL

SOURCE DIMENSIONS IN INCHES. FORMAT ( 12G6 . 3

)

12 34567890123456789012 34 56789012345678901234567890123 4567 89 012345678901234567890

IN.

2.72 1.310 .100 .250 .070 .250 .100 .03 .0 .0 .0 CA AA AB TRBI TLBI TSB TRCI TLCI TSC TRDI TLDI TSD

I I I I D SIDE WALL IN

I I I D LEFT END IN INCHES

I I D RIGHT END IN INCHES

I C SIDE WALL IN INCHES

C LEFT END IN INCHES

I I I I C RIGHT END IN INCHES

I I I B SIDE WALL IN INCHES

I I B LEFT END IN INCHES

I B RIGHT END IN INCHES

OUTSIDE DIAMETER IN INCHES

OUTSIDE LENGTH IN INCHES

4.2.3 SIMPLIFIED SAMPLE OUTPUT

CYLINDRICAL PU-BE M-621 SOURCE. ENCAPSULATED IN TANTALUM AND ST. STEEL,

SOURCE DIMENSION INPUT IN INCHES, L = 2.72, OD = 1.31, B WALL = 0.07,

B ENDS = 0.25 AND 0.1, C WALL = 0.03, C ENDS = 0.1 AND 0.25.

MACROSCOPIC CROSS SECTIONS PER CM

SIGAA= 1.865 , SIGAF= 4.305 , SIGB= 1.16 , SIGC=

NEUTRONS PER FISSION = 2.80. NZ=10. NPHI=10.

FOR A SINGLE NEUTRON STRIKING THE SOURCE...

THE TRANSMISSION PROBABILITY IS THE INTERACTION PROBABILITY IS THE ABSORPTION PROBABILITY IS THE 'A' MATERIAL ABSORPTION INTERACTION PROBABILITY IS THE «A' MATERIAL FISSION INTERACTION PROBABILITY IS THE 'A« MATERIAL INTERACTION PROBABILITY IS

0.2811

.18433302+0

.81566698+0

.54509177+0

.11721783+0

.27057522+0

.38779305+0

20
THE CLADDING MATERIAL INTERACTION PROBABILITY IS .42787393+0

NEUTRONS OUT PER NEUTRON STRIKE IS .94194362+0

NEUTRONS LOST PER NEUTRON STRIKE IS .58056380+0

THE AVERAGE CROSS SECTION OF THE SOURCE IN CM**2 IS .20170300+2

THE PERCENT ERROR IN THE TRANSMISSION PROBABILITY IS .0162

IF THE MEASURED THERMAL FLUX IN PERCENT OF Q PER CM**2 IS .122 THE INTERACTION IN PERCENT OF IS .20071742+1

THE XBSORPTION IN PERCENT OF Q IS .13413491+1

THE FISSION IN PERCENT OF Q IS .18643104+1

AND THE RESULTANT LOSS IN PERCENT OF Q IS .14286378+0

CYLINDRICAL SOURCE PU-BE A. ENCAPSULATED IN NICKEL.

SOURCE DIMENSION INPUT IN INCHES, L = 1.031. OD = 1.031, B WALL = 0.128,

B ENDS = 0.112 AND 0.129.

MACROSCOPIC CROSS SECTIONS PER CM

SIGAA= 1.214 , SIGAF= 2.802 , SIGB= 0.4106.

NEUTRONS PER FISSION = 2.80. NZ=15. NPHI=15.

FOR A SINGLE NEUTRON STRIKING THE SOURCE...

THE TRANSMISSION PROBABILITY IS .31458511+0

THE INTERACTION PROBABILITY IS .68541489+0

THE ABSORPTION PROBABILITY IS .38173937+0

THE «A» MATERIAL ABSORPTION INTERACTION PROBABILITY IS .13157105+0

THE »A' MATERIAL FISSION INTERACTION PROBABILITY IS .30367551+0

THE 'A' MATERIAL INTERACTION PROBABILITY IS .43524656+0

THE CLADDING MATERIAL INTERACTION PROBABILITY IS .25016833+0

NEUTRONS OUT PER NEUTRON STRIKE IS .11648765+1

NEUTRONS LOST PER NEUTRON STRIKE IS -.16487655+0

THE AVERAGE CROSS SECTION OF THE SOURCE IN CM**2 IS .77943026+1

THE PERCENT ERROR IN THE TRANSMISSION PROBABILITY IS .0270

IF THE MEASURED THERMAL FLUX IN PERCENT OF Q PER CM**2 IS .122 THE INTERACTION IN PERCENT OF Q IS .65176438+0

THE ABSORPTION IN PERCENT OF Q IS .36299785+0

THE FISSION IN PERCENT OF Q IS .80854631+0

AND THE RESULTANT LOSS IN PERCENT OF Q IS -.15678192+0

CYLINDRICAL AM-BE SOURCE. ENCAPSULATED IN TANTALUM AND STAINLESS STEEL.

SOURCE DIMENSION INPUT IN INCHES, L = 1.355, OD = 1.355, B WALL = 0.07

B ENDS = 0.07, C WALL = 0.03, C ENDS = 0.05.

MACROSCOPIC CROSS SECTIONS PER CM

SIGAA= .1704 , SIGAF= 0.00082 » SIGB= 1.16 , SIGC= 0.281

NEUTRONS PER FISSION = 2.89. NZ=10. NPHI=10.

FOR fi SINGLE NEUTRON STRIKING THE SOURCE...

THE TRANSMISSION PROBABILITY IS .40678453+0

THE INTERACTION PROBABILITY IS .59321547+0

THE ABSORPTION PROBABILITY IS .59253429+0

THE «A' MATERIAL ABSORPTION INTERACTION PROBABILITY IS .14155335+0

THE «A' MATERIAL FISSION INTERACTION PROBABILITY IS .68118396-3

THE 'A' MATERIAL INTERACTION PROBABILITY IS .14223453+0

THE CLADDING MATERIAL INTERACTION PROBABILITY IS .45098094+0

NEUTRONS OUT PER NEUTRON STRIKE IS .40875315+0

NEUTRONS LOST PER NEUTRON STRIKE IS .59124685+0

THE AVERAGE CROSS SECTION OF THE SOURCE IN CM**2 IS .13463194+2

THE PERCENT ERROR IN THE TRANSMISSION PROBABILITY IS .0356

IF THE MEASURED THERMAL FLUX IN PERCENT OF Q PER CM**2 IS .121

THE INTERACTION IN PERCENT OF Q IS .96637561+0

THE ABSORPTION IN PERCENT OF Q IS .96526593+0

THE FISSION IN PERCENT OF Q IS .32069762-2

AND THE RESULTANT LOSS IN PERCENT OF Q IS .96316863+0

21
Table 1.

Neutron Transmission Convergence Test for the Spherical

Source Program

NY Transmission Estimated

Error

Actual

Error*

10 0.14653888 0.819$ 1.205$

20 0.14542886 0.438$ 0.438$

50 0.14495702 0.113$ 0.113$

100 0.14485061 0.041$ 0.039$

500 0.14479752 0.008$ 0.002$

1000 0.14479397 0.004$ —

* Actual Error listed here assumes that transmission

for NY = 1000 is correct.

22
Y =

1(a) VIEW FROM -Y

x = o

1(b) VIEW FROM +X

z=o

1 (c) VIEW FROM +Z

z = z

1 (d) VIEW FROM +Z

Figure 1. Diagram for the cylindrical source program.

The cylindrical source is surrounded by an

imaginary sphere with center on the X axis

in the middle of the core material. Figure

1(d) shows the intersection of a plane parallel to the X-Y axis with the source at Z.

23
0.3170

0.3160

O 0.3150

if)

if)

if)

<

q:

0.3140

0.3130

0.3120

A-NZ=I0

?'

4 8

NZ OR NPHI

12 16

Figure 2. The calculated Transmission and Estimated

Error for a nickel encapsulated Pu-Be cylindrical source for various specified subdivisions of the Z and cp integrals, NZ and

NPHI. The solid line joins points for whicl'

NPHI equals 10 and the dashed line joins

points for which NZ equals 10. Two points

are shown for both NPHI and NZ equal to 5

and 15.

2k
REFERENCES

1. Spiegel, Jr., V. and Murphey, W. M. : Metrologia (1971)

(In press).

2. Murphey, W. M. : Nuclear Instr. Methods 37_, 13 (1965).

3. Mosburg, Jr., E. R. : J. Research Natl. Bur. Standards

62, 189 (1959).

4. Scarborough, J. B. : Numerical Mathematical Analysis , Baltimore: The Johns Hopkins Press: London: Humphrey

Milford Oxford University Press 1930; pp. 120 and 155.

25

Latest developments in the subject area of this publication, as well as

in other areas where the National Bureau of Standards is active, are

reported in the IMBS Technical News Bulletin. See following page.
HOW TO KEEP ABREAST OF IMBS ACTIVITIES

Your purchase of this publication indicates an interest in the research,

development, technology, or service activities of the National Bureau of

Standards.

The best source of current awareness in your specific area, as well as in other NBS programs of possible interest, is the TECHNICAL NEWS BULLETIN,

a monthly magazine designed for engineers, chemists, physicists, research

and product development managers, librarians, and company executives.

If you do not now receive the TECHNICAL NEWS BULLETIN and would like to subscribe, and/or to review some recent issues, please fill out and return the

form below.

Mail to: Office of Technical Information and Publications

National Bureau of Standards

Washington, D. C. 20234

Name

Affiliation

Address

City State Zip

| | Please send complimentary past issues of the Technical News Bulletin.

| | Please enter my 1— yr subscription. Enclosed is my check or money

order for $3.00 (additional $1.00 for foreign mailing).

Check is made payable to: SUPERINTENDENT OF DOCUMENTS.

TN 576
orm NBS-114A d-71)

U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OR REPORT NO. NBS-TN-576

2. Gov't Accession No.

3. Recipient's Accession No.

4. TITLE AND SUBTITLE

Computer Code for the Calculation of Thermal Neutron

5. Publication Date May 1971

Absorption in Spherical and Cylindrical Neutron

Sources

6. Performing Organization Code

7. AUTHOR(S)

V. Spiegel, Jr. and W. M. Murphey

8. Performing Organization

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS

10. Project/Task/ Work Unit No. 5030220

DEPARTMENT OF COMMERCE

WASHINGTON, D.C 20234

11. Contract/Grant No.

2. Sponsoring Organization Name and Address Nuclear Radiation Division

Institute for Basic Standards

13. Type of Report & Period teal

National Bureau of Standards

Washington, D.C. 20234

14. Sponsoring Agency Code

5. SUPPLEMENTARY NOTES

6. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here.) A computer code has been written in FORTRAN IV for

the calculation of thermal neutron absorption in spherical

and cylindrical neutron sources. The formalism of the calculation, the structure of the computer code, a listing of

the code, and some sample results are presented. The comparison of the results of this calculation to experiment

appears in Metrologia 7, No. 1, 34-38 (1971).

17. KEY WORDS (Alphabetical order, separated by semicolons)

Neutron; neutron standards; manganous sulfate bath calibration of neutron

sources.

18. AVAILABILITY STATEMENT

[XI UNLIMITED.

19. SECURITY CLASS

(THIS REPORT)

UNCLASSIFIED

21. NO. OF PAGES

Uo

1 I FOR OFFICIAL DISTRIBUTION. DO NOT RELEASE

TO NTIS.

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

22. Price 35 cents

USCOMM-DC 66244-P71

* U. S. GOVERNMENT PRINTING OFFICE : 1971 431-998/303

NBS TECHNICAL PUBLICATIONS

PERIODICALS

JOURNAL OF RESEARCH reports National

Bureau of Standards research and development in physics, mathematics, chemistry, and engineering.

Comprehensive scientific papers give complete details of the work, including laboratory data, experimental

procedures, and theoretical and mathematical analyses. Illustrated with photographs, drawings, and

charts. Published in three sections, available separately: • Physics and Chemistry

Papers of interest primarily to scientists working in these fields. This section covers a broad range of physical and chemical research, with major emphasis on standards of physical measurement, fundamental

constants, and properties of matter. Issued six times a year. Annual subscription: Domestic, $9.50; for- eign, $11.75*. • Mathematical Sciences

Studies and compilations designed mainly for the mathematician and theoretical physicist. Topics in mathematical statistics, theory of experiment design,

numerical analysis, theoretical physics and chemistry, logical design and programming of computers

and computer systems. Short numerical tables. Issued quarterly. Annual subscription: Domestic,

$5.00; foreign, $6.25*. • Engineering and Instrumentation

Reporting results of interest chiefly to the engineer

and the applied scientist. This section includes many

of the new developments in instrumentation resulting

from the Bureau's work in physical measurement,

data processing, and development of test methods.

It will also cover some of the work in acoustics,

applied mechanics, building research, and cryogenic

engineering. Issued quarterly. Annual subscription:

Domestic, $5.00 ; foreign, $6.25*.

TECHNICAL NEWS BULLETIN

The best single source of information concerning the Bureau's research, developmental, cooperative and

publication activities, this monthly publication is designed for the industry-oriented individual whose

daily work involves intimate contact with science and

technology

for engineers, chemists, physicists, re- search managers, product-development managers, and

company executives. Annual subscription: Domestic,

$3.00; foreign, $4.00*.

• Difference in price is due to extra cost of foreign mailing.

NONPERlODICALS

Applied Mathematics Series. Mathematical tables, manuals, and studies. Building Science Series. Research results, test methods, and performance criteria of building materials, components, systems, and structures. Handbooks. Recommended codes of engineering

and industrial practice (including safety codes) de- veloped in cooperation with interested industries,

professional organizations, and regulatory bodies. Special Publications. Proceedings of NBS conferences, bibliographies, annual reports, wall charts,

pamphlets, etc. Monographs. Major contributions to the technical

literature on various subjects related to the Bureau's

scientific and technical activities. National Standard Reference Data Series.

NSRDS provides quantitative data on the physical

and chemical properties of materials, compiled from

the world's literature and critically evaluated. Product Standards. Provide requirements for sizes, types, quality and methods for testing various industrial products. These standards are developed cooperatively with interested Government and industry groups and provide the basis for common understanding of product characteristics for both buyers and sellers. Their use is voluntary.

Technical Notes. This series consists of communications and reports (covering both other agency and

NBS-sponsored work) of limited or transitory interest. Federal Information Processing Standards Publications. This series is the official publication within

the Federal Government for information on standards

adopted and promulgated under the Public Law

89-306, and Bureau of the Budget Circular A-86

entitled, Standardization of Data Elements and Codes

in Data Systems.

Consumer Information Series. Practical information, based on NBS research and experience, covering

areas of interest to the consumer. Easily understandable language and illustrations provide useful background knowledge for shopping in today's technological marketplace.

NBS Special Publication 305, Supplement 1, Publications of the NBS, 1968-1969. When ordering, include Catalog No. C 13. 10: 305. Price $4.50;

foreign, $5.75.

Order NBS publications from: Superintendent of Documents

Government Printing Office

Washington, D.C. 20402
U.S. DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20230

OFFICIAL BUSINESS POSTAGE AND FEES PAID

U.S. DEPARTMENT OF COMMERCE

PENALTY FOR PRIVATE USE, $300