NBSIR 87-3534

DCTDOS: Neutron and

Penetration in Composite Duct

Systems

L. V. Spencer

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Radiation Research

Gaithersburg, MD 20899

February 1987

epared for:

r\~ ;deral Emergency Management Agency

'0 ashington, DC 20472

T<e

H7

NBSIR 87-3534

DCTDOS: NEUTRON AND GAMMA

PENETRATION IN COMPOSITE DUCT

SYSTEMS

L. V. Spencer

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Radiation Research

Gaithersburg, MD 20899

February 1987

Prepared for: Federal Emergency Management Agency

Washington, DC 20472

U.S. DEPARTMENT OF COMMERCE, Malcolm Baidrige, Secretary

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

DCTOOS: NEUTRON AND GAMMA PENETRATION IN COMPOSITE DUCT SYSTEMS

page

List of Figures iii

I. Introduction 1

II. Structure of DCTDOS 5

III. Input Data and Options .............. . 7

A. Run-Description Card 7

B. Problem Description Cards . 8

C. Configuration Cards 10

IV. Output Selections and Elementary Examples 13

V. Comparison of Results with Some Classic Experiments 17

A. ORNL Studies 17

B. Armour Research Foundation Experiments Using Co 60 Gamma Rays. 23

VI. Procedures, Variables, and Basic Data 28

VII. Concluding Remarks and Acknowledgment 36

References and Bibliography . .............. 37

Appendix A. Some Notes on Subprograms 44

Appendix B. Sources of Data 52

Appendix C. Remarks on Inputs and Outputs with Examples 55

Appendix D. FORTRAN Listing of DCTDOS ....... 62
I
List of Figures

Page

Fig. 1. Diagram of DCTDOS ........ .......... 6

Fig. 2. Example of elementary input and output of DCTDOS ........ 14

Fig. 3. Example of elementary input and output of DCTDOS 15

Fig. 4. Detector pattern for ROOMX and DOSX, NW = 3, NL = 4 16

Fig. 5. Top view of source and duct geometry 18

Fig. 6. Thermal -neutron flux in various duct configurations due to

incident thermal neutrons 19

Fig. 7. Thermal -neutron fluxes in various duct configurations due to

incident epicadmium neutrons . 20

Fig. 8. Comparison of calculated and measured center-line singlecollision fast-neutron dose rates 21

Fig. 9. Cumulative epithermal neutron spectra ....... 24

Fig. 10. Gamma dose distribution measurements in the 6 by 6 foot

concrete entranceway ...................... 25

Fig. 11. Measured gamma dose rate in 6 by 6 foot concrete entranceway

with two right angle bends 27

Fig. Cl. Example of output from I PR I NT=3 , 56

Fig. C2. Input data for computations of the ORNL duct .... 60

Fig. C3. Input data for computations of the ARF duct 61

i i i

DCTDOS : NEUTRON AND GAMMA PENETRATION IN COMPOSITE DUCT SYSTEMS

L. V. Spencer

National Bureau of Standards

Washington, D.C. 20234

I . Introducti on

This paper describes computer methods for estimating neutron and

gamma ray fluence rate, dose, and even spectral features due to

penetration through a series of duct segments. The procedure links

together data for individual segments -- straight sections and bends --

in arbitrary combinations; and the resulting composite can include

computations for a room at the end, if there is one. This particular

method was developed for rapid estimates in the protection problems

against nuclear weapons; but the concepts which are employed are more broadly appl icabl e.

The study of radiation penetration in ducts and other cavities has

developed a large literature, which can be subdivided into a) studies

of neutron and gamma ray albedos, 5 ’ 13 ’ 18-24 ’ 29 ’ 38 ’ 43-45 ’ 58 ’ 61 ’ 78

, b)

experiments, l»3-4»6»8-9»12»15-16»27-28»30»33-34»37»39-40»48»54’59»

70-71

, c) Monte Carlo duct analysis, 17 ’ 25 ’ 33-37 ’ 43 ’ 46 " 50 ’ 52 " 59 ’ 75

, and

d) some form of analytical analysis. 7 ’ 12 " 13 ’ 26 ’ 31-32 ’ 41 ’ 53 ’ 60 ’ 62-70

Summary discussions of the various problems, data, and methods are to

be found in many places, among which are Refs. 15,41-42,55-58, and 74.

Many reactor shielding studies are not referenced above.
The program described here concludes a series of studies which

started in the middle and latter 1960's. In these studies various

aspects of radiation penetration through ducts were examined in regards

to the possibility of viewing complex duct configurations as composites

of elementary segments. 67 " 69 Many ideas appearing in this paper were described previously in Ref. 69. While this approach fits mainly with

category d) above, it is neither a quadrature nor a curve-fitting

method, nor does it emphasize center-line dose.

This report presents a usable computer program (DCTDOS) based on concepts described in that earlier study and gives applications to some

well known experimental duct studies. Several of our basic assumptions

are as follows:

k

1) That the relevant flux, which depends on both space and

direction variables, can be usefully described by introducing the

order-of-scatteri ng index n (or N=n+1). In DCTDOS we include up to

n = 20, which is barely large enough for complex systems with large

reduction, and may eventually be enlarged.

2) The albedo angular distribution for neutrons is well

approximated by the elementary distributions of Chandrasekhar 1

1

for

isotropic scatterers. 69 Since the sensitivity of most problems to

details of this angular distribution is rather low, we use it to obtain

solutions to elementary spatial configurations in terms of the variable

ru That is, we use one-velocity, i sotropic-scatterer albedos to

obtain data for simple segments as functions of ru

3) The neutron flux component which has been backscattered from

concrete walls (say) a total of n times has a spectral distribution

(for a given radiation source) which is assumed to depend on energy

only through _n. This takes advantage of a near-factori ng of the singlereflection albedo into a directional part and an energy part. 67-69 But

Here and throughout this manuscript the term flux stands for fluence

rate. (See ASTM Standard E170, 1986 Annual Book of ASTM Standards , Section 12.) Flux can also designate a time integral, i.e., a fluence

if the source is likewise integrated over time.

2
in turn this leads us to assume a factoring of the whole problem into

space-and-di rection studies, and mul ti -refl ection energy studies.

Tables of mul ti -refl ection energy albedos as functions of both initial

and eventual spectral energy are required. 69

4) In using the output of one segment as input into another, some information about the direction of particle currents is a requirement.

We have reduced this requirement to its least restrictive form by using

only a single parameter: For this we use the mean reciprocal cosine of

the obliquity, and which we refer to in the convenient form

(-1 + 2/<sece>) = s. Note that <sece> is also the ratio

fl ux/current . We record and use s as a matching parameter,

distinguishing between unreflected and reflected components in the

preceding of any two segments being matched, since these can be very

di fferent

.

(5) Radiation reflected in cylinders of circular cross section,

spheres, and between infinite plane surfaces has the property of

reentrance at the emergence obliquity when travel in the wall can be

ignored. Hence mul ti -refl ecti on albedos in these configurations are identical if the initial incidence distribution and spectrum are identical. We take advantage of the fact that cylinders with circular

and square cross sections have a rough equivalence, as do cubes and

spheres. Thus we believe that in a fairly usable first approximation

mul ti -refl ection distributions between infinite plane surfaces has a general applicability to far more complex configurations.

Combination of segments proceeds by simple convolutions of

successive functions of n. Spectra and different types of dose are extracted by additional convolutions with mul ti -refl ecti on energy data.

Hence DCTDOS marches to its solutions by a series of convolutions over data interpolated for elementary segments. As may be readily guessed,

such a process is fast and can be flexible.

Most of the data on duct sections which DCTDOS utilizes were developed in a series of unpublished studies by S. Woolf. 76 ” 77 The

3
mill ti -refl ection albedo energy data for neutrons is due to Hagan and

Simmons, 38 and correspondi ng data for gamma rays has been computed by

Spencer and Woolf. 78 Since all data are rather simple in concept and

readily obtained on file, we anticipate supplementation and replacement

if these methods are further developed.

Gamma ray albedos are more complex than neutron albedos, with strong

dependence on energy above 0.5 MeV [and for single scattering]. If the

single-scattering is given special attention, there is merit in treating

gamma rays as if they were neutrons, but using mul ti -refl ecti on albedo

energy data for gamma rays specifically. This simple approach turns out

to work rather well. This confirms that the n = 1 reflection component

for gamma rays dominates the difference from the correspondi ng neutron

angular distribution, and also that the grazing incidence case,

sometimes termed "skipping," is comparatively unimportant. Differences

of albedo magnitude in the basic segment data can be roughly corrected

by use of simple ratios raised to the n'th power. Recomputation of the

space-angular data for gamma rays specifically has not yet been done.

In Section V we present results for both neutrons and gamma rays in

comparison with experimental data and results from well known and

detailed computations: for neutron data in the Maerker and

Muckenthaler duct 46-49 and for gamma ray data in one of the AMF

ducts. 73 We discuss spectral and order-of- refl ecti on comparisons as well. A high degree of precision is neither attained nor expected; but

we feel that absolute values and trends are sati sfactori ly reproduced

for both types of radiation.

We have organized this report to put practical information necessary

to apply DCTDOS towards the front, with explanatory and

supplementary material located towards the rear. The comparison of

results for classic experiments constitutes a transition section

between these two types of information.

4
II. Structure of DCTDOS

In the diagram of Fig. 1, the capitalized words represent

subroutines, the function of each being suggested by choice of name. Data input are of three different types. (INPUT) represents problem

parameter data. The basic data on which the program depends are introduced by calling ENTDAT to set up elementary segment data and

SRSDET for spectra and mul ti -refl ection response weights, both making

use of types of mul ti -refl ection albedo energy data.

Parameter input begins with new Run data, and later is indicated by

four parallel lines designated Run, Problem, Length, and Weights. A new Length changes only a single (configuration) parameter, while a new

Run brings in all the data necessary for a new duct computation. A new Problem uses previous source and detector information, while "New

Weights" changes source and/or detector while retaining the previous

confi guration.

Branch points have been identified by lower-case letters a_ through i_, their significance being as follows:

a) Choice between DCTISO (isotropic source flux) and DCTSLT

(mono-obliquity source angular distribution).

b,c,d) The option here has to do with completion vs. an additional

segment of duct, according to the complexity of the

confi guration.

e) The configuration here may call for room computations or for

dose computations in the duct.

f) Choice of mul ti -refl ecti on albedo energy data and application

of these data to the dose weights: neutrons (SND) or gamma

rays (SGD).

g) This option permits recomputation with different dose weights.

The rerun through SRSDET gives wide latitude to this option.

h) This option permits modification of one of the straight segment

lengths, followed by recalculation with all other parameters

fixed

.

i) The choices here are for a new configuration with unchanged

source and detector data (16) or a fully new situation (New

Run, 13). Termination is also an option at this point.

5
b I AH I

WTLISTS » - DCTDOS

STOP

6
III. Input Data and Options

We use the term "card" here with reference to a line of input data,

which formerly was punched on an IBM card.

In preparing for a DCTDOS run, several different types of data card

are required. The first specifies mainly source and detector data to

be used, and we refer to this as a Run-Description card. Following

this is a card with index numbers mainly for specifying outputs, new or old dose weights, and a series of segment lengths or detector positions,

together with a title card. We refer to these as Probl em-Descri pti on cards. Finally, there is a set of cards giving dimensions for a sequence of duct segments, which we refer to here as Conf i gurati on cards. Sample run decks are indicated in Figs. 2 and 3; here we mainly

give descriptions, together with the options currently available.

A. Run-Description Card . The first run card specifies five

integers, followed by a decimal number in exponential form, using the

Fortran format ( 515, 1 PE10 .3) . The numbers, in their order, have the

following interpretations:

ISORS specifies the type of radiation spectrum selected,

IDOSE specifies the type of detector response selected,

ING specifies gamma rays or neutrons, with the inclusion or dismissal of wall -capture gamma rays,

LPIPE selects from several types of (straight) duct data available

for neutrons (only a single type of gamma ray data is used),

IEL first specifies the type of output by quick identification

with IPRINT. Later this parameter has a variety of other

uses through its COMMON storage in most subprograms.

DIMEN is a numerical factor which can change or rescale output

units.

7
Table 1 lists options currently available according to values of the

five index numbers referred to above. But special spectra or response

functions can be readily inserted, either by adding to the lists above,

or by insertion as a replacement of one of the cases listed in Table 1.

We have found this useful in substituting 137 Cs source data, say, for a 60 Co gamma spectrum.

Table 1. Options currently available, as specified by the different

input integers.

ISORS

1: Neutron INR spectrum (neutrons/group/incident dose), see Ref. 2,

Table 10)

2: Thermal neutron spectrum (normalized to unity)

3: Epithermal neutron spectrum (neutrons/sec/watt/cm 2

) incident. See

Ref. 48, Fig. 2, based on 34,200 n/sec/watt into (91.44 cm) 2

. 4. Gamma ray INR spectrum (gammas/group/incident dose). See Ref. 2,

Table 10).

5. 60 Cobalt gamma rays (normalized to unity). 50 Curies, 2 photons

per interaction; but only half moves forward into a duct of

6'x6'cross sectional area. This gives 3.319E9 gammas/cm 2 /mi n

IDOSE

1: Snyder-Neufel d dose ( rads/( n/cm 2

)

2: Unit weights (for flux computations, e.g. neutrons/ ( n/cm 2

)

3: Henderson tissue dose data ( rads/ ( n/cm 2

)

4: Thermal neutron response only (neutrcns/( n/cm 2

)

5: Gamma-ray absorbed dose data (rads/(gamma/cm 2

)

ING LPIPE

-1: Wall -capture gammas only 1: Finite-duct (exit) neutron data

0: Neutrons only, or gammas only 2: Semi -i nfi ni te-duct data, forward

1: Neutrons + wall -capture gammas 3: Ditto, both forward and return

1-3: Any value for gamma rays

IEL

0: ( I PR I NT ) - - No print-out, first case;

1-4: ( I PR I NT ) — Print-out (Table 2, also see discussion in Section IV)

DIMEN

This constant factor can change dimensions for the weights and/or

renormalize the input current.

B. Problem Description Cards . Let us first dismiss the title card,

which is the last of these cards. Its format is presently (55H ... ).

8
The first input card has again the format ( 515 , 1PE10 .3) . The

integers specified have the options and i nterpretations listed in

Table 2 below, with entries in the same order as on the card.

Following the main input card is a card with a single integer

specifying the number (JDOS) of alternative sets of weights. And this

is followed by JDOS alternative-weights cards, each again using the

format (515, 1 PE10 . 3 ) . These cards list values for IPRBLM, ISORS,

IDOSE, ING, LP, and DIMEN in that order. It has turned out that of

these, IPRBLM is a dummy without application, and can be equated to its

value on the first input card.

Table 2. Options for the "Problem" Parameters

IPRBLM:

-1,-2,-

1: The after this after this one is to use the same run- weights card.

2: The computation after this one will have a new run-weights

card.

3: This is the last computation.

-3: The present computation will use weights already computed

and stored from the run-weights card which precedes.

NEL: This gives the number of duct elements (plus one, for the

source). This is also the number of case cards which follow.

IPRINT: Acceptable values are 0, 1, 2, 3, 4. Zero indicates

no printout; 1 means printout of input information

together with computed dose and dose ratio; 2 also gives

weights plus dose information at earlier junctions; 3 adds

group spectra and mean secant directional information; and 4

adds intermediate printout from the subroutines for use in

checking parameter sizes and other inquiry.

JVARY: This gives the index ( IC) for the straight section whose

length variation may be desired. Thus, JVARY = 2 would refer

to the first straight section of a duct.

JDEL: This index gives the number of values which the straight

section specified by JVARY is to have. JDEL = 1 means that

only one value is to be computed, the most common case: The

values of JVARY and DEL are then not relevant and DEL may be

omitted or set to zero.

DEL: This fixed point number gives the increment by which the

length of the JVARY segment is to be increased for each of

JDEL cases. It is important that DEL is length,

i.e., length divided by the square root of the cross sectional

area.

9
C. Configuration Cards . These cards in sequence describe the

configuration. We discuss them in terms of the names of the lists into

which they are put (in DATDCT). These lists are NEXT, WDTH, HGHT,

LNTH, OTHR, SELECT: The Format used is (15, 5F10.3)).

The first entry, NEXT, gives IEL values for a series of subprograms

to be called in the sequence given (DCTDOS, order 35). These are listed in Table 3 below for easy reference:

Table 3. NEXT indices used to call the main subprograms.

IEL SUBROUTINE Type of Section

0 Descriptive parameters, only used with IC = 1 1 IXTISO Straight segment, isotropic source

2 DCTSLT Straight segment, slant source

3 TURN Bend segment

4 BTODR Straight segment occurring after a bend

5 BTODR Straight segment after a bend

6 R00MX Room segment

7 DOSE Applies weights and computes total, detector in or at

the end of, a bend or straight section

8 OOSX Applies weights and computes total, detector pattern

in a room

The remaining entries on the case description cards include widths,

heights, and lengths. We summarize these in Table 4 using the

abbreviations, W, H, and L, respecti vely . We use dimensioned numbers,

usually in feet, in practice; the program does its own scaling. Those

parameters which are not W, H, or L require special explanations, which

fol 1 ow Tabl e 4.

Table 4. Input data assignments for the case description cards. (Note

Appendix B.) NEXT corresponds to IEL in Table 5.

NEXT WDTH HGHT XLTH OTHR ELECT

I ZD SO PEN ZD00R RMCASE

1 W H L ARAT S = 0

2 w H L ARAT S

= -1,0,1

3 w H RCSS F S = -1,0,1

4,5 w H L ARAT S = -1,0,1

6 W,WTH H L WD S = -1,0,1

7 RDLNTH S = -1,0,1

8 WDP NW-1 NL-1 AEX S = -1,0,1

10
The following interpretations apply to the various quantities in

Tables 3 and 4.

SELECT(IC) is normally zero or omitted. If it equals 1., the reflected

component is equated to zero for computations by the

subroutine selected by NEXT ( IC ) . If it equals -1., the

unreflected component is equated to zero. If 0., the normal

case, both contribute. Note that "unreflected" and

"reflected" here refer only to the preceding section and not

to N = 1 and N > 1 in general.

SO This index is zero for the cosine (current) source. In

general this quantity is = -1 + 2./<sece

0 > for slant

obliquity 9

0

.

PEN is an assumed effective scattering depth in side-walls prior

to reemergence. It has normally been set equal to zero; its

dimensions are those of W, H, and L.

ZDOOR is a door or wall thickness assumed (for IZD>0) at the

entrance to the shelter. Note: It is limited to dimensions

of feet of concrete presently, and assumes no modification

of source spectrum or directional distribution.

RMCASE has values +1., 0., -1.: The value +1. is used if the flux

entering the room (i.e., R00MX) and unreflected by the

preceding segment, is directed towards a nearby sidewall.

The value 0. is assigned if this flux component has no strong

bias towards or away from a sidewall. The value -1. is

assigned if this flux component (VELO) is directed towards

the array of detectors. Most commonly, RMCASE = 0.

RCSS gives the depth of the recess, in cases of a bend with a trap. Units are the same as for other XLTH values.

F is an azimuthal factor for the radiation entering the bend.

Normally it has the value zero. (1 - F) is applied as a factor to weight the component which traverses the bend

without wall reflections in the bend. Hence, for F=l, the source entering the bend cannot exit without reflections,

while for F=-l, the current exiting without reflections

(i.e., VELO) is doubled. (Note that the configuration can modify this component in ways not expressed by basic

bend data.)

RDLNTH Measures the distance from detector to the backscatteri ng

endwall of a straight duct segment. It is used in RETRN to

estimate and include backscattered current. It should be

omitted or set zero if this component is ignored, which is

the most common case.

WD is the width of the entrance door into a room section

(R00MX)

.

11
WDP

NW,NL

AEX

ARTO

RAT

ARAT

is the width of the pattern of detectors, i.e. the distance

from the center of the duct exit to the wall on the side

where the detectors are positioned, in the shelter room

( ROOMX)

(both fixed point numbers quickly converted to integers) are

lateral and longitudinal numbers of detectors in the array

(see ROOMX). NW*NL < 20, presently, by dimension statement.

is a ROOMX parameter. It is the area of the entrance to

the room (IEL=6). Note that IEL = 7,8 are alternatives;

and DOSX, rather than DOSE would be used with ROOMX.

is a ratio: (initial duct cross section area)/(final duct

cross section area). Here, initial and final refer to the

first and last duct segments. For DOSX, the final aperture

area is the same as AEX.

is the ratio of duct cross-sectional area taking account of

a mean scattering depth (PEN) to that with PEN = 0. It would normally be unity.

is the ratio of cross-sectional areas at a junction (the

value for the next segment divided by that for the

preceding) if there is a reduction.

12
IV. Output Selections and Elementary Examples

The input parameter I PR I NT has values 0-4 as al ready noted. The

value I PR I NT = 0 suppresses printout and has no interest to us here.

The value I PR I NT = 4 prints additional intermediate data and parameters

,

and is useful mainly to inquire into the causes of hard-to-expl ai

output; we will not attempt to discuss this option further. Desired

intermediate output may change with one or another particular

requirement or investigation.

I PR I NT = 1 causes only input data, together with minimum output, to

be printed. Input data have been described in the preceding section and are readily recognizable in the printout. Also these cards and their

sequence can be confirmed from the list INPUT lines added to the top of

Figs. 2 and 3. OUTPUT begins by recording the input data under suitable

headings. This is followed by values for the dose and dose reduction

factor as indicated by printed titles. The heading OUTPUT has been

added to clarify the figure and does not accompany the output. The

heading "NEXT CASE = .500," begins the printout of results for Fig. 2,

and this heading appears in general with output results. The number

.500 gives distance down the final straight segment specified by JVARY,

and may be zero if the dose is computed at the entering point for a new segment.

If JDEL = 1, the usual case, this distance will be the seal ed value

corresponding to XLTH for the segment identified by JVARY. Thus, the

duct of Maercker and Muckenthaler (see Fig. 5) has a mean lateral

extension of /(3‘ x 3

1 ) =3

' ; and for JVARY = 4 the value NEXT CASE = .500

identifies the position .5*3' = 1.5' down the segment identified by IC =

4. Since IC = 2 is the first straight segment and IC = 3 is the first

bend, IC = 4 is the second straight segment. The possibility for

incrementing the XLTH value specified by JVARY by one or more distances

DEL up to a total of JDEL-1 times exists. Then output of the form

described is repeated, occurring JDEL times in succession, with

appropriately changed headings, as shown in Fig. 2.

13
INPUT

215144 2031 1.973+000

1.333+000

0 *** M+M, 2-LEG, THERMAL FLUX (USES BTODR) , 8/29/86 ***

0 0.414 0.000 0.000 0.000 0.000

2 3.000 3.000 13.500 0.000 0.000

3 3.000 3.000 0.000 0.000 0.000

5 3.000 3.000 1.500 0.000 0.000

7 0.000 0.000 0.000 0.000 0.000

OUTPUT

•p ISORS I DOSE ING LPIPE IEL DIMEN

2 2 0 3 1 1. 973E+00

IPRBLM NEL IPRINT JVARY JDEL DEL

1 5 1 4 4 1. 333E+00

*** M, 2-LEG, THERMAL FLUX (USES BTODR), 8/29/86 ***

IEL WIDTH-FT HEIGHT-FT LENGTH-FT OTHER SELECT

0 .414 .000 .000 .000 .000

2 3.000 3.000 13.500 .000 .000

3 3.000 3.000 .000 .000 .000

5 3.000 3.000 1.500 .000 .000

7 .000 .000 .000 .000 .000

NEXT CASE .500

NEXT ( I )= 5 D0SE= 2. 361E-02 REDUX FACTR= 8.459E-03

NEXT CASE 1 .833

NEXT ( I

) = 5 D0SE= 4. 369E-03 REDUX FACTR= 1.565E-03

NEXT CASE 3.166

NEXT ( I

) = 5 DOSE= 1.215E-03 REDUX FACTR= 4.355E-04

NEXT CASE 4.499

NEXT ( I

) = 5 DOSE= 5.283E-04 REDUX FACTR= 1.893E-04

FIG. 2. Example of elementary input and output of DCTDOS, Data are for

a duct with concrete walls and 3'x3' passageway cross section.

Output is for the 2nd leg, at 1.5', 5.5', 9.5', and 13.5'

beyond the bend, which in this case is cubic in shape. The

(neutron) source is incident at 45° at all parts of the

entrance plane and all azimuths.

14
INPUT

2 2 0 3 1 1.973+000

3

0

**3X3X18

4 1 2 1 0.000

1-LEG DUCT INTO ROOM, THERMAL N FLUX, 8-15-83

0 0.000 0.000 0.000 0.000 0.000

1 3.000 3.000 18.000 0.000 0.000

6 10.000 8.000 15.000 3.000 0.000

8 5.000 1.000 2.000 9.000 0.000

OUTPUT

$ =$

ISORS IDOSE ING LPIPE IEL DIMEN

2 2 0 3 1 1.973E+00

IPRBLM NEL I PRINT JVARY JDEL DEL

3 5 1 2 1 .00CE+00

**3X3X18 1-LEG DUCT INTO ROOM, THERMAL N FLUX, 8-15-83

IEL WIDTH-FT HEIGHT-FT LENGTH-FT OTHER SELECT

0 .000 .000 .000 .000 .000

1 3.000 3.000 18.000 .000 .000

6 10.000 8.000 15.000 3.000 .000

8 5.000 1.000 2.500 9.000 .000

7 .000 .000 .000 .000 .000

NEXT CASE 6.000

DOSE/ENTRANCE DOSE W(FT)= 5.00 L(FT)= 15.00

1 . 0653E-02 9.0436E-03 5.3005E-03

1.4363E-03 1 . 5541E-03 1 .3522E-03

DOSE

4.2036E-02

5 .6678E-03

3.5686E-02

6.1 324E-03

2.0916E-02

5 . 3359E-03

Fig. 3 Example of elementary input and output of DCTDOS. A single

straight segment 3'x3'xl8' ends in a room of dimensions

lO'xS'xlS 1 . The source of thermal neutrons is isotropic at all

parts of the entrance plane. See Fig. 4 for the detector

pattern: For 6 detectors, 1, 3, 5, are opposite the entrance

and 2, 4, 6, are against the sidewall.

15
DOSX output for I PR I NT = 1 gives arrays rather than single numbers.

Thus in Fig. 3 output is in the form of two arrays -- for the dose

ratio as well as for the dose itself. The heading for the former is

DOSE/ENTRANCE DOSE, and for the latter, DOSE, as shown in Fig. 3. Both

arrays have, as their upper left value, the entrance value to the room. The top line gives a sequence along the entering axis in equidistant

positions crossing the room. Lower lines of detectors are displaced

laterally in the room by equal amounts from entrance to sidewall and

from front to back, as shown in Fig. 4 below.

1 4 7 10

2 5 8 11

3 6 9 12

Fig. 4. Detector pattern for ROOMX and DOSX, NW = 3, NL = 4.

This figure differs from that for the Fig. 3 data in its

proportions and in a non-centered entrance.

I PR I NT = 2 introduces several additional types of information.

Values are printed for the dose weights; also given are DOSE and

DOSE/ ( ENTRANCE DOSE) data for the successive segment junctions of the

system. Thus to obtain the dose weights used, one need run only the

simplest of programs using I PR I NT=2 . But dose weights are also to be

found in the file PREWTS. The types of data which are given with

IPRINT=2 may be required, or they may provide a better understanding as to why the dose or dose reduction turn out to have the values printed.

I PR I NT = 3 extends the output to include spectral values together

with minimum intermediate parameters, mostly mean secants at junctions.

I PR I NT - 4 prints additional intermediate values. Consult the

program listing in Appendix D to identify and possibly supplement this

seldom-used data selection.

16
V. Comparison of Results with Some Classic Experiments

The comparisons which follow in this section are absol ute , based on computations using known source strengths and absolute detector data.

No relative comparisons have been included.

A. ORNL Studies . The most extensive studies of neutrons in concrete

ducts are those carried out at ORNL by Maerker, Muckenthal er , and

others. 17 - 24 ’ 43 * 51 Their 3' x 3‘ duct eventually reached 45' in length,

with options of 1, 2, or 3 bends as shown in Fig. 5. In all their

studies the source was a beam, directed through the opening at an angle

of 45° relative to the duct axis. Measurements were made at points on the centerline.

The ORNL studies included computations of the albedo angular and

energy functions for concrete, with a follow-on computer Monte Carlo

study based on the resulting albedo data. These computations gave good

agreement with experiments; and they also permitted examination of

intermediate spectra and orders-of-refl ection distributions. Hence

there is a satisfying completeness about the information given in this

set of reports. (Note particularly Ref. 49.)

Data were developed for a thermal neutron source and for a fast

neutron source. In the latter case as well as in the former, thermal

neutron flux was measured and computed. For the latter case measurements and computations for a Hurst dosimeter were also

devel oped.

Fig. 6 shows the thermal neutron flux due to a thermal neutron

source. Fig. 7 gives the thermal neutron flux due to an epithermal

neutron source. Fig. 8 similarly gives the fast neutron dose due to

the "epithermal" neutron source. In each case we have added points

representing computations using DCTDOS, and the added points are readily identified. The following paragraphs discuss aspects of

these (DCTDOS) results.

17
Fig. 5. Top view Of source and duct geometry.

18
THERMALNEUTRON

FLUX

(neutrons

-cm /

-sec

ORNL-DWG 66-8671R

Fig. 6. Thermal -neutron flux in various duct configurations due to

incident thermal neutrons

19
SUBCADMIUM

FLUX

(neutrons

cm

*

sec'

0 5 tO 15 20 25 30 35 40 45 50

CENTER-LINE DISTANCE (ft)

Mg. 7. Thermal -neutron fluxes in various duct configurations due to

incident epicadmium neutrons.

20
FAST-NEUTRON

DOSE

RATE

(ergs

ORNL-DWG 66-1730R3

£

cn

0 6 12

,U 18 24 30 36 42

CENTER-LINE DISTANCE FROM MOUTH (ft)

Fig. 8. Comparison of calculated and measured center-line singlecollision fast-neutron dose rates.

21
A divergence near the source most likely reflects the difference

between the experimental beam source and our source, which covered the

aperture uniformly and at all azimuths, and which used interpolated data

for the 45° obliquity angle rather than expressing that obliquity

sharply.

An increasing divergence beyond 22' in the first straight segment

is, we believe, a failure of our model. But note that a somewhat

similar effect occurs in the ORNL computations. It may well indicate

that near-grazing trajectories are poorly described in both models. In our computations the discrepancy occurs at distances requiring

extrapolation beyond our data set; but a faulty extrapolation formula

does not seem to account fully for the effect.

A rise as a bend is approached is probably due to back-reflection

from surfaces marking the endwall of the bend. This can be estimated

if required.

High values computed by DCTOOS just following a bend can result in

part from differences between a centerline detector and the DCTDOS

detector, which averaged over the lateral plane. This will be evident

when lateral intensity variations are strong, as occurs when part

only of the duct cross section is exposed to intensities direct from

two segments preceding and possibly part of the source plane.

The general agreement with the experiments appears to be

comparable with that of the Monte Carlo computations, if account is

taken of the differences between source and detector configurations.

To pursue comparisons in greater detail, we computed both cumulative

group spectra and orders-of-refl ection components at several positions

in the duct. The ORNL studies published computed examples of both

types. Orders-of-refl ecti on distributions were compared for equal

centerline distances both down a duct continuation and around a bend.

Results show that the bend shifts the peak to values for the orders-of22
reflection which are larger by 4-5, for both first and second bends.

Our computations showed about the same shift. But our peaks were consistently located at smaller orders-of-scatteri ng by about the same amount. We think that this results partly from our use of a source configuration which permits neutron injection initially to greater

depths down the first leg. But in both ORNL computations and in ours,

peak location is not strongly sensitive to distance down a given leg.

We chose cumulative rather than differential spectra for a comparison because the ORNL energy group structure is quite different

from ours. 49 Fig. 9 gives the comparison between the two types of

computed spectra at 14' beyond the center of the first bend. Below

about .05 MeV the ORNL spectrum rises higher to a maximum excess factor

of about 4/3 at an energy just above the thermal component. This is

consistent with, although it may not fully result from, the larger

number of scatterings apparently undergone by the ORNL flux at the same position. Note that the thermal components are about equal at this

position, as shown in Fig. 7 at 29'.

B. Armour Research Foundation Experiments Using Co 60 Gamma Rays . These experiments were not supplemented with computations; hence

the types of comparison are much more limited than the above, although

these data have also be analyzed since by other authors.

In Fig. 10, we show the experimental duct with three legs used by

Terrell et al. 70-73 We have computed gamma attenuation in this concrete

duct and compared the data from DCTDOS with the experimental data of

Ref. 73. Our computations used a one-group approximation to the Co 60

spectrum. Similarly, we do not distinguish between three bends of the

Z shape as distinct from the U shape, a difference which seems to be

of minor importance here. Note that a poi nt , 50 Curie source was located at S 10 , 3' inside the first leg at a point which we take to

23
Cumulative

Fluence

Rate

Neutrons/cm

2

/sec/Watt

'III' I I r jTTTTT I I I JTTTTT I I I J

I I ? I I I | I JTTTT

3x10 - 4

o Maercker et al

• This Report

2x 10 -4

1 x 10 -4

o •

I ' ' ' i I I

| | | ! li i i i ' ' I I Inn i I| | 1 1 1 I 10 -1 10 -2 10-3

Lower Energy Limit

10

Fig. 9. Cumulative epithermal neutron fluence rates.

2.855

2.153 <-

o

10

" 5

‘

24
GkAPHiTC

Pace

Sio

4—3 —£44- !0

' -H

p\

Fig. 10. Top view of the 6 by 6 foot concrete entranceway.

25
be the duct entrance for the DCTDOS computations. Recall that these

are based on an isotropic input which is assumed uniform over the

lateral entrance plane. Again the centerline measurements are compared

with an average over the lateral plane used by DCTDOS. We have also

ignored source gammas back-reflected into the duct.

The comparison is given in Fig. 11. Our results show less

satisfactory agreement with the experiments than is the case for

neutrons, as might be expected. The most evident difference appears

wthin and just beyond the bend, where our computations do not

reproduce the feature that as long as the (point-like) source can be

viewed by the detector, the main fall off in the second leg is

postponed. The consequence of this is that our main fall off occurs 1-2

feet later than in the experiment. We believe here also that at least

pert of the explanation for this effect results from a centerline

approach vs an approach by lateral averaging. At the first bend, the

large experimental reduction occurs when the point detector is first

shielded from the point source. By contrast, a plane source is only

gradually shielded from a plane detector, in this case total shielding

occurring nearly 7' beyond the vertical lines which mark the turns to a new leg at the bend centers.

Regarding the factor of roughly 2 discrepancy in the third bend, we think that this is due to inadequacies of data expressing our model. In addition to differences of source and detector configuration, no attempt

has at this point been made to incorporate and apply gamma ray albedo

angular distributions to computations of basic segment data; neutron

albedo angular distributions have been used instead, as noted. The

substitution of more realistic data should improve accuracy somewhat, at

the expense of duplicating much of ENTDAT.

26
/min

1000 . 0

100. 0

10. 0

1.0 —

0 . 1

0 . 01

0 . 005

0 . 002 -

0 . 001

0

A

&

Normalized

here

O— Point Source C

• - DCTDOS

6 •A

O

&

O

Centerline Distance in Feet

_1 1 1 1

1 1 L_

12 16 20 24 28 32 36

Z

shape

40 44

Fig. 11. Measured gamma dose rate in 6 by 6 foot concrete entranceway

with two right angle bends (Z-Shape).

27
VI. Procedures, Variables, and Basic Data

1) Position, Energy, and Direction (N, SO, SP, IC, DLNGTH, SENVAL).

We identify position with the entrance or exit lateral planes in a straight duct segment, except in the special case of a detector pattern

in a room. _IC_ values sequence duct segments, and DLNGTH is straight

segment length or detector position.

The energy variable is expressed through the order of

backscatteri ng , _n or N=(n+1). While there is expected to be a dependence of the energy on the sequence of types of backscatteri ng

e.g., angles of incidence and emergence, in practice this dependence is

assumed to be weak and we keep track only of the number. In addition,

there is a dependence of the albedo on the neutron energy, which we neglect because it is not strong except in special cases which do not

concern us here. Our basic data (fcr concrete) is based on an assumed

albedo value of 0.696 throughout although more recent evaluations

suggest that 0.75 is a oetter figure. Rather than recompute the basic

data at this point, we neglect the differences of angular distribution

as well as spectrum and merely include a correction factor

(0.75/0. 696

n to take account of this, because these differences do not

appear to have a strong influence on segment outputs. Since the dose

weights have incorporated the correct albedo factor, only the

correction factor ( 1/0 . 696

n is applied in DOSE and D0SX.

The angular distribution is character!' zed , basically, by the ratio

f 1 ux/current . This is an average of the reciprocal cosine over the

radiation current. (It is also the reciprocal of the mean cosine averaged over the normalized flux.) We use this parameter

partly because it gives a simple way of transforming currents into

fluxes and vice versa. (Such a relationship can only be approximate

because dose or other energy weighting factors give variant flux-tocurrent ratios.

28
We have found it convenient to use a derived quantity

s 2 <sece>

1 . ( 1 )

Values of s for forward currents range from -1 to +1 ; and the largest

value (<sec 9 >=l) corresponds to neutrons travelling parallel to the

axis of penetration, which is nearly always the reference axis.

Dependence on only this single direction parameter is one of the basic

limitations of our approach; but we have found that it works rather

well. It is not difficult to identify configurations and directional

distributions which would be poorly described so as to lead to

unsatisfactory results. But such configurations are about as likely to

suggest a requirement for new data types as to suggest a change in

approach.

In DCTDOS, N_ or NREF always identifies (n + 1) rather than n_; and

S0(_IC_) , SP(JC_) correspond to ^values for unreflected and reflected

components, respectively, emerging from the segment J_C. (But note that

_IC = 2 is the first real segment of a duct system being evaluated.)

SENVAL is £ value for the current enteri ng a bend.

2) Dependent Variables (VELO, VELP, RMTX). Our computations are based on currents, rather than fluxes, and VELO, VELP identify currents

entering or leaving a duct segment. Specifically, VELP( I C , N

)

i s the

component of the current of neutrons reflected a total of (N-l) times

and emerging from the duct segment identified by IC, having had at

least one reflection in the IC'th segment . Correspondingly, VEL0( I C , N

)

is the component of the emerging current which passed through the IC'th

segment without any reflections. Its dependence on _N signifies that it

has no doubt been reflected in other segments, perhaps many times.

RMTX is intended to stand for Yefl ecti on matrix^' The current

entering any duct segment will be reflected 0, 1, 2, ... n times in

that segment. These are in addition to, and must be combined with,

previous reflections. The general procedure for doing this is the

29
evaluation of a simple folding computation. This is cast in the form

of a matrix operation in which the incoming currents, VELO and VELP,

are treated as vectors multiplied by an elementary matrix, RMTX, which

expresses the effect of the new segment. This product then yields new "vectors" VELO and VELP constituting the output current. If one has a vector with components V

Q

, V

x

, V

2

, ... together with a matrix which is zero to the right of the diagonal, has the value C

Q on the diagonal, C

2 in the positions just left of the diagonal, C

2 in all the positions

just left of the C

x elements, etc., then the product is

o

o 0 0 0 Vv o C V ^0 v 0

c

i o

o 0 0 Vw1 C V + C V

C

2 C1 oo 0 V2 CU V2

V 0 + Cu V

l v1 + ^0^2

C

3 C2 ci oo V3 C,V n3 0 + C V1 + C1 V2 + C0 V3

( 2 )

Let the subscripts now stand for numbers of reflections, the V be

n values for the current entering an element, and C

n for the probability

of emerging reflections in this element. Then the matrix to

the right has a form suitable for the corresponding current emerging

from the element. This assumes that the V are essentially

n interchangeable upon entrance. Sometimes we make this assumption; but

more often we perform different computations for VELO and VELP , since the

entering angular differences seem often to be important for these

components

RMTX always is given the general form of the matrix shown above, but

if VELO and VELP are treated separately, RMTX should more accurately be

computed separately for each, a refinement not in DCTDOS presently.

3) Data Forms for Structural Elements (BEND, DUCTO, DUCT1, DUCT2)

S. Woolf has performed many computations for duct segments of square

cross section and different lengths or different positions in a very

long, straight duct. 76 Other data from these computations explore more

general shapes: rectangular rather than square cross section shapes.

30
Still other computations are for bends of square cross section but with

lengths 1 x 1 (cubical), 1x2, etc.; and recessed bends were examined.

Study of Woolf's data leads to the empirical result that successive

orders of reflection beyond n = 1 nearly always follow rather well the

simple rule

C

n = A*a

n

, for n > 1 , (3)

where C

n is the n'th reflected current at the duct exit, say. We

assign to all of these components combined the angular parameter P. The n = 0 case is usually quite different, and hence is assigned the

different value SO, as already stated.

Thus, to render Woolf's data in a form suitable for our computations, we require five parameters:*

KX = 1: Co

2: A

3: a

4: SO

5: SP

*These types of data are used in the basic programs OCTISO, DCTSLT, and

TURN, where the sequence of orders-of-ref 1 ecti on computed is not that

of Eq. (3) above, but a somewhat more complicated form which agrees

better with the basic data. (Note the remarks about the subroutine

ORDST, in Appendix A.

Fixing the square cross section as basic, and the lxl bend shape

likewise as basic, we have, in the case of straight ducts, either the

duct length or the distance to the detector as variable: All the above

types of parameter are thus functions of this distance parameter. By

contrast, in the case of the duct bend the dimensions are fixed but

the angular input is not; and all of the above types of parameter are

31
functions of the value of the incoming current. Thus, straight duct

data are tabulated as a function of length or distance to the detector,

both identified by DLNGTH . Bend data are rendered as a function of the

s_ value ( SENVAL ) for the entering current, which covers the entrance

surface.

Straight duct data are based on an entering current proportional to cose relative to the duct axis in DUCTO , and on entering currents with

mono-obliquity angles cose = 0.25 and 0.75, in DUCT1 and DUCT2 . Other

data are available, but we have preferred to interpolate or combine

these, as required. In addition, data for finite ducts (LP = 1), for a. fixed distance in a semi -infinitely long duct (LP^ = 3), and for the

forward current only in the latter case (LP = 2) are included in ENTOAT.

LPIPE and _LP specify the input data type to be used.

4) Interpolation, Extrapolation of Segment Data . We rely mainly on a simple 2-point interpolation procedure applied to the logarithm of data

points. Thus, if we wish data for duct of length L which falls between

>

* tabulated values l

] and L

?

, and if the tabulated values for these two

points are y L and y 2 respectively, the value y or the desired value at L.

is given by

y = y i (y 2 /y i

)

^ L=” L

1

^

1-2-1

^

This is exact for exponential distributions.

In using such a simple interpolation procedure, various data sets

are modified to take account of the following problems:

a) reduce range of variation by known divisor functions for direct

use of DATA statements,

b) ensure that neither y L nor y 2 is zero in a tabulation,

c) counter a strongly non-exponential variation.

Such modifications are removed after an interpolation.

Also, an extrapolation procedure has been included to make possible

rough estimates for extreme cases such as very long and narrow ducts.

32
This is based on the assumption that the distribution being extrapolated

can be described by the general form

y

= Ax“ p + Bx

_

m

" £

(5)

for very large values of x. (See the discussion of the subroutine

XTRAP.)

5) Combination of Distributions . If one wishes to combine additive

distributions characteri zed by different values of s such as VELO and

VELP, a value for £ for the combination must be obtained by taking

proper account of the nature of s^ as a mean value. Let us designate

two currents with values s_

: and s_2 as ^(s_L ,cose) and C

2 (s_2 ,cose)

Total currents are designated c^, £2

, respectively. Then

<sece>

i

l

d(cose)sece[C

1 (s1 ,cose) + C2 (s2 ,cose)]

c**1 + c

2

c

1 *<sec9 1 > c

2 *<sec92 >

c

i + c2 Cx + C

2

( 6 )

Equation (1) can be used to express this equation in terms of s and s.

,

i=l ,2.

6) Correction for Eccentric Duct Cross Sections . Data due to

Woolf 76 show that for the reflected component (VELP) in a straight

duct, a rough correction for eccentricity (W not equal to H) can be

made which is relatively independent of duct length over the range of our data. This correction is not highly accurate, but it only becomes

large when the ratio of side widths becomes more extreme than 0.1. A

simple fit to this correction factor which is independent of n has been

i ncl uded.

For the unreflected component (VELO), the correction is more

complex. It turns out that duct length is more important for this case.

The empirical formula which we have incorporated attempts to modify the

simple form for the unreflected component to include the distance

33
parameter in about the simplest manner that can represent the data.

The results are not as dependable as in the case of the reflected

component; but the correction is smaller and again it does not become a major factor for duct eccentricities much less extreme than 0.1.

If e<l is the ratio of the smaller to the larger of the two lateral

dimensions, then our correction factors have the forms

K

d

= exp( - .109(1 ne) 2

) (7)

K

0 = exp( - ,06( 1 ne) 2 )L/ {V )

The factor Kp multiplies the output list VELP, while K

0 multiplies

the output list VELO, both in DCTISO and in DCTSLT.

No correction of this type is presently applied to TURN.

7) Sudden Changes of Duct Size . If the duct cross sectional area should suddenly increase, the total current would remain unchanged and

would rather quickly expand to fill the larger lateral space. By

contrast, a sudden reduction in the cross sectional area of a duct, say

by a factor DELR, would in first approximation decrease the total

current by DELR but would not otherwise be expected to change spectrum

or angular distribution.

8) Computations of Spectrum and Dose . At the conclusion of a sequence of duct elements the currents VELO and VELP are represented by

angular parameters SO, SP, respectively and by values for all orders of

reflection, and the problem then is to turn these data into spectra,

dose reduction factors, and possibly fluences. Let us consider the

problem of spectra first.

A basic requirement at this point is an evaluation of the neutron

spectrum as a function of initial energy and order of reflection.

To this end we have defined a special energy-loss problem in which

neutrons in a succession of 23 energy groups from 15 MeV to thermal are injected into the space between plane concrete walls with directions

34
proportional initially to cose. If the initial energy group is

identified with the index IE, and the final energy group by KE, the

resulting function A(IE,KE,N) is to be applied for most of our purposes.* The final spectrum of C(KE) of a current V(N) resulting

from an input source with spectrum S(IE) would then be simply

CO 00 C(KE) = l S( IE) I V(N)*A( IE.KE.N). (8)

I E= KE N=1

If the dose assigned to each spectral group is designated D(KE), it may

be obtained as usual with

Dose = l D( KE) C(KE) , (9)

KE=1

= l S( IE) l V(N) l

( D( KE ) A( IE ,KE ,N) . IE=1 N=1 KE=1

(9’)

In eq (9'), the brackets identify a two-index quantity which we designate WTSD. Rearrangement of (9') leads to an alternative form

with an interchange between positions of D and S. The corresponding

bracket term is then designated SWTS.

9) Flux and Flux Ratios . The entering current can readily be

weighted by dose factors to obtain a quantity which we refer to as "current dose". Note that N = 1 is the only component, and presumably

has the the input spectrum. Similarly one can compute a "current dose"

for other positions, as well as the ratios to that at the entrance.

The problem of the flux and flux reduction factors is more complex

in that currents with different angular distributions must be combined.

*Later considerations suggest that it may be desirable to supplement

with a function for mono-obliquity incidence cases.

35
Our approach is the simple procedure of computing current dose for

emerging SO and SP components which are then transformed by the simple

rul e

*current = <sec0>*current

( 10 )

Appropriate combinations and ratios then readily follow which have an accuracy which we think is consistent with many applications.

Additional information is given in the Appendices.

A listing of the computer program is presented in Appendix D.

Appendix A describes the various subroutines and their purpose.

Appendix B describes various types of data used and gives some information about data sources. Appendix C gives additional examples

of input and output of the program.

Albedo data has not been printed here, partly because it nearly

doubles the number of pages and partly because it is most practical in

magnetic disk form for direct use with computers.

Despite its apparent promise and potential usefulness, we do not

regard the DCTDOS program to be completed in the sense of some of

the best-known Monte Carlo programs currently available. But we hope

that the ideas incorporated will stimulate additional work of this type

leading to improvements and possibly applications to similar problems

in other technical fields such as reactor engineering.

In regard to ideas and data developed herein, the author appreciates

the many contributions from, and interest expressed by Dr. Stanley

Woolf, of ARCON Corp. in Wakefield, Massachusetts. The use of "we" and

"our" throughout this manuscript expresses the author's feeling of

partnership during the investigation. Dr. G. L. Simmons of SAIC in

San Diego assisted similarly with important aspects of the work.

VII. Concluding Remarks and Acknowledgment

36
REFERENCES AND BIBLIOGRAPHY

1. Albright, G.H., Beck, E.J., LeDoux, J.C., and Mitchell, R.A.,

Evaluation of Buried Conduits as Personnel Shelters , Report

No. WT-1421, Operation Plumbbob, Project 3.2, Defense Atomic

Support Agency (1960).

2. Auxier, J.A., Burson, Z.G., French, R.L., Haywood, F.F., Mooney,

L.G., and Straker, E.A., Nuclear Weapons Free-Field Environment

Recommended for Initial Radiation Shielding Calculations , Report

No. ORNL-TM-3396, Oak Ridge National Laboratory, Oak Ridge,

Tennessee (1972).

3. Beck, C. (Editor), Test of French Underground Personnel Shelters , Report No. WT-1453, Operation Plumbbob, Project 30.6, Defense Atomic

Support Agency (1962).

4. Beck, C. (Editor), Nuclear Weapons Effects Tests of Blast-Type

Shel ters , Report No. CEX-68.3, Civil Effects Test Operations,

U.S. Atomic Energy Commission (1969).

5. Berger, M.J., and Morris, E.E., Dose Albedo and Transmission

Coefficients for Cobalt-60 and Cesium-137 Gamma Rays Incident on Concrete Slabs, Report No. 9071, National Bureau of Standards

(1966).

6. Bramblett, R.M., Ewing, R.I., and Bonner, T.W., Nuclear Instruments

and Methods 9, 12 (1960).

7. Byoun, T.Y., Babel, P.J., and Danjani, A.T., Computer Program

SCAP-BR for Gamma-Ray Streaming Through Multi legged Ducts,

ANS Transactions 27, 777 (1977).

8. Cain, V.R., A Study of the Shielding Characteristics of Basic

Concrete Structures at the Tower Shielding Facility , ORNL-3464,

Oak Ridge National Laboratory, Oak Ridge, Tennessee (1964).

9. Cain, V.R., Clifford, C.E., and Hallard, L.B., Measurements of

Radiation Intensities in Vertical Concrete-Lined Holes and an Adjoining Tunnel at the Tower Shielding Facility , ORNL-3513,

Oak Ridge National Laboratory, Oak Ridge, Tennessee (1964).

10. Celnik, J., A Simple Albedo Method for Calculating Effects of

Single-Legged Ducts, Nuclear Reactor Shielding , p. 688,

Roussin, R.W., Abbott, L.S., Bartine, D.E., Editors (1977).

11. Chandrasekhar, S., Radiation Transfer, Dover Publications, Inc.,

New York, N.Y. (1960]“::

37
12. Chapman, J.M., Dose-Attenuation Variation with Incident Gamma-Ray

Energy in Two Legged Concrete~~and Steel Ducts, Nuclear Reactor

~~~

Shi el di ng , Tech. Report Y-F008-08-05-201 , Defense Atomic Support

Agency (1966)

13. Chilton, A.B. and Huddleston, C.M., A Semi empi rical Formula for

Differential Dose Albedo for Gamma Rays on Concrete, Nucl . Sci . Eng

27 , 419-424 (1962); earlier published as a Defense Atomic Support

Agency Report (1962).

14. Chilton, A.B., Shultis, J.K., and Faw, R.E., Principles of

Radiation Shielding, Chapter 9, Prentice-Hall, Inc., Englewood

Cliffs, NJ 07632 (1984)

15. Clifford, C.E., Gamrna Shielding Provided by Ducts , Report No. 370.

Defense Research Chemical Laboratories, Ottawa, Canada (1962).

16. Clifford, C.E.; Simons, G.G.; Faw, R.E.; and Shultis, J.K.,

Benchmark Gamma Ray Streaming Experiment , Report No. RRA-T7903

Radiation Research Associates, Fort Worth, Texas (1979).

17. Clifford, C.E.; Muckenthal er, F.J.; Maerker, R.E.; Straker, E.A.

and Mynatt, F.R.; Neutron Transport Studies in Ducts and Simple

Shields: Comparison of Monte Carlo and Discrete Ordinates

Calculations with Experimental Results from the Tower Shielding

Faci 1 i ty , Report No. ORNL-TM-1965, Oak Ridge National Laboratory,

Oak Ridge, Tennessee (1967).

18. Coleman, W.A., A Determination of Intermediate-Energy Neutron

Albedo Data for Concrete Using Monte Carlo , Report No. ORNL-TM1205, Oak Ridge National Laboratory, Oak Ridge, Tennessee (1965).

19. Coleman, W.A. and Maerker, R.E., Differential Neutron Current

Albedos for Concrete in the Incident Energy Range 0.5 eV to 200 keV . II. Results of Monte Carlo Calculations for an Incident Angle of~

~~

0 deg ., Report No. QRNL-3967, Oak Ridge National Laboratory, Oak

Ridge, Tennessee (1966).

20. Coleman, W.A. and Maerker, R.E., Differential Neutron Current

Albedos for Concrete in the Incident Energy Range 0.5 eV to 200 keV . III. Results of Monte Carlo Calculations for an Incident Angle of

45 deg . , Report No. ORNL-3967, Oak Ridge National Laboratory, Oak

Ridge, Tennessee (1966).

21. Coleman, W.A. and Maerker, R.E., Differential Neutron Current

Albedos for Concrete in the Incident Energy Range 0.5 eV to 200 keV . IV. Results of Monte Carlo Calculations for an Incident Angle of

60 deg . , Report No. 0RNL-3967, Oak Ridge National Laboratory,

Oak Ridge Tennessee, (1966).

38
22. Coleman, W.A. and Maerker, R.E., Differential Neutron Current

Albedos for Concrete in the Incident Energy Range 0.5 eV to 200 keV . V. Results of Monte Carlo Calculations for an Incident Angle of

75 deg . , Report No. ORNL-3967, Oak Ridge National Laboratory, Oak

Ridge, Tennessee (1966).

23. Coleman, W.A. and Maerker, R.E., Differential Neutron Current

Albedos for Concrete in the Incident Energy Range 0,5 eV to 200 keV . VI. Results of Monte Carlo Calculations for an Incident Angle of

85 deg , Report No. ORNL-3967, Oak Ridge National Laboratory, Oak

Ridge, Tennessee (1966).

24. Coleman, W.A., Maerker, R.E., Muckenthal er, F.J., and Stevens,

P.N., Calculation of Doubly Differential Current Albedos for Monodirectional Beams of Epicadmium Neutrons Incident or. Concrete and

Comparison of the Reflected Subcadmium Component with Experimental

Results, Nucl . Sci . Eng . 27 , 411-422 (1967). Earlier published as Report No. ORNL-TM-1554, Oak Ridge National Laboratory, Oak Ridge,

Tennessee (1966).

25. Collins, D.G., A Monte Carlo Multibend Duct Procedure , NARF-62-13T

Nuclear Aerospace Research Facility, Fort Worth, Texas (1962).

26. Collins, D.G., French, R.L., Mooney, L.G., Tompkins, K.W., and

Wells, M.B., Simplified Methods for Calculating the Penetration of

Neutron and Gamma Radiation into Underground Concrete Structures , Report No. RRA-M45, Radiation Research Associated, Inc., Fort

Worth, Texas (1964).

27. Corps of Engineers, Huntsville Division (IJ.S. Army), 100 Man Blast

She! ter . Design Analysis Volume 1, March 1983.

28. Doty, D.R., Dose Measurements for Neutron Streaming in Ducts , Technical Report No. R282, U.S. Naval Civil Engineering Laboratory,

Port Hueneme, California (1964).

29. Doty, D.R., Differential Neutron Albedo for Finite Slabs, Nucl. Sci.

Eng. 27, 478 (1967).

30. Eisenhauer, C.M., Scattering of Cobalt-60 Gamma Radiation in Air

Ducts , Tech. Note 74, National Bureau of Standards (1960).

31. French, R.L., Wells, M.B., and Schaeffer, Penetration of Neutron and

Gamma Radiation Through the Openings of Underground Structures , Report No. RRA-T41, Radiation Research Associates, Inc., Fort Worth,

Texas (1964).

32. French, R.L., Wells, M.B., and Schaeffer, N.M., Neutron and Gamma

Ray Penetration into Concrete-Shielded Underground Structures , Report No. RRA-T51, Radiation Research Associates, Inc., Fort

Worth, Texas (1965).

39
33. Gardner, L.B., Neutron and Gamma-Ray Streaming through a TwoLegged, Thick-Wall Steel Duct , Technical Report No. R-558, Defense

Atomic Support Agency (1968).

34. Gardner, L.B., Neutron and Gamma-Ray Streaming through Thin-Steel

Shelter Entranceways , Technical Report No. R-569, Defense Atomic

Support Agency (1968)

35. Gardner, L.B. and Mettler, A.J., Monte Carlo Calculations of

Gamma-Ray Streaming through Coplanar and Noncoplanar Ducts , Report No. R412, Defense Atomic Support Agency (1966).

36. Gardner, L.B. and Mettler, A.J., Monte Carlo Calculation of Neutron

Streaming through Two-Legged Duct Entranceways , Report No. R379,

Defense Atomic Support Agency (1965).

37. Gardner, L.B.; Mettler, A..].; and Burdick, B., Experimental

Determination of Neutron Energy Spectra in Concrete Ducts , Report

No. R523, Defense Atomic Support Agency (1967).

38. Hagan, W.K. and Simmons, G.L., Neutron Spectra from Monoenerget i

Source Neutrons after Multiple-Reflection between Plane-Parallel

Concrete Interfaces , Report No. SA I-013-79-605-LJ , Science

Applications, Inc., La Jolla, California (1979).

39. Horton, C.C.; Halliday, D.B.; Lakey, J.R.A., and Harrison, J.R.,

The Transmission of Thermal Neutrons in a Straight Annular Duct , Report No. SWP/P.27, U.K. Atomic Energy Agency, Harwell ( 1956)

40. Huddleston, C.M. and Wilcoxson, W., Gamma Ray Streaming Through

Ducts , Report No. R289, Defense Atomic Support Agency (1964).

41. Huddleston, C.M., Gamma Ray Streaming through Ducts , Symposium on Protective Structures for Civilian Populations, National Academy of

Sciences, Washington, D.C. (1965).

42. Jaeger, R.G. Engineering Compendium on Radiation Shielding, Vol I . Shielding Fundamentals and Methods , Sections 4.4 and 5.4.1, and

"Chapter 8, Spri nger-Verl ag. New York Inc. (1968).

43. Maerker, R.E. and Muckenthal er, F.J., Cal cul ational and Measurement of the Fast-Neutron Differential Dose Albedo for Concrete,

Nucl . Sci . Eng . 22, 455-462 ( 1965).

44. Maerker, R.E. and Muckenthal er, F.J., Differential Fast-Neutron

Albedos for Concrete. I, Descriptions of Monte Carlo Calculation

and TSF Experiment and Comparison of Calculated and Experimental

Resul ts , Report No. ORNL-3822, Oak Ridge National Laboratory,

Oak Ridge, Tennessee (1965).

40
45. Maerker, R.E. and Muckenthal er , F.J., Measurements and SingleVelocity Calculations of Diff erential Angular Thermal -Neutron

Albedos for Concrete , Report No. 0RNL-4090, Oak Ridge National

Laboratory (1967). Also Report No. ORNL-TM-1421 (1966).

46. Maerker, R.E. and Muckenthal er, F.J., Monte Carlo Calculations,

Using the Albedo Concept, of the Fast-Neutron Dose Rates along the

Center Lines of One- and Two-Legged Square Concrete open Ducts and

Comparison with Experiment, Nucl . Sci . Eng . 27 , 423-432 (1967);

earlier published as Oak Ridge National Laboratory, ORNL-TM-1557

(1966).

47. Maerker, R.E. and Muckenthal er, F.J., Calculations, Using the

Albedo Concept, of Thermal -Neutron Fluxes along the Center Lines of

One-, Two- and Three-Legged Square Concrete open Ducts arising from

Thermal Neutrons Entering Duct Mouths: Comparison with Experiment , Report No. ORNL-TM-1758, Oak Ridge National Laboratory (1967).

48. Maerker, R.E. and Muckenthal er, F.J., Neutron Fluxes in Concrete

Ducts Due to Incident Epicadmium Neutrons: Calculations and

Experiments , Report No. ORNL-TM-1794, Oak Ridge National Laboratory

( 1967 ) . Also published in Nuclear Science and Engineering 30,

340-354 (1967) .

“ ~

49. Maerker, R.E. and Muckenthal er, F.J., Calculations, Using the

Albedo Concept, of Thermal -Neutron Fluxes, Epicadmi um-Neutron

Fluxes, and Fast-Neutron Dose Rates along the Center Lines of

Square Concrete Ducts: Comparison with Experiment , Report No.

ORNL-4147, Oak Ridge National Laboratory, (1967).

50. Maerker, R.E. and Cain, V.R., AMC: A Monte Carlo Code Utilizing

the Albedo Approach for Calculating Neutron and Capture Gamma-Ray

Distributions in Rectangular Concrete Ducts , Report No. ORNL-3964,

Oak Ridge National Laboratory, (1967).

51. Maerker, R.E., Claiborne, H.C., and Clifford, C.E., Neutron

Streaming in Rectangular Ducts , Oak Ridge National Laboratory,

( 1966) . Later published as Subsection 8.2.5, Engineering

Compendium on Radiation Shielding, Vol . 1 , Spri nger-Verl ag

,

New York Inc. ( 1968)

52. Moore, J.A., Eggen, J.B., and Leopard, F.O., Monte Carlo Procedure

for Analysis of Radiation in Ducts (RAID) , Report No.

AFWL-TR-67-9, Air Force Weapons Laboratory, Fort Worth, Texas

(1967).

53. Park, C.M., Agnihotri, C.B., and Silverman, J., Interim Report on Scattering of Gammas through Ducts , Report No. UMNE-2, Dept, of

Chemical Engineering, University of Maryland, College Park, MD

(1962).

41
54. Sayedahmed, F.M., Belov, S., Megahid, R., and Hamouda, I., Neutron

Transmission through Ducts Perpendicular to Neutron Beam Axis in ' Water Shields , Report No. U.A.R.A.E.E./Rep.-86, Atomic Energy Post

Office, United Arab Republic Atomic Energy Establishment, Cairo,

U.A.R. (1969).

55. Schaeffer, N.M. (Editor), Reactor Shielding for Nuclear Engineers , Chapter 7, Available in hardback as TID-25951 from the National

Technical Information Service, Springfield VA 22151 (1973).

56. Selph, W.E.; Newell, J.M.; and Bayer, C.E., Analysis of Neutron

Streaming through an Underground Hardened Installation , Report No.

RRA-T82, Radiation Research Associates, Fort Worth, Texas (1967).

57. Selph, W.E. and Claiborne, H.C., Methods for Calculating Effects of

Ducts, Access Ways, and Holes in Radiation Shields , Report No.

"0RNL-RSIC-20, Oak Ridge National Laboratory (1968).

58. Selph, W.E., Neutron and Gamma-Ray Albedos Report No. ORNL-RS IC -21

Oak Ridge National Laboratory (1968).

~~

59. Shin, K., Murakami, R., Taniuchi, H., Hyodo, T., and Oka, Y„,

Measurements of Neutron and Gamma-Ray Streaming in a Cavity-Duct

System and an Analysis by an Albedo Monte Carlo Nucl . Sci . Eng. 81,

161-171 (1982).

60. Shore, F.J. and Schamberge, R.D., The Transmission of Neutrons

through Ducts in Water , Report No. BNL 390 (T-74), Brookhaven

National Laboratory, New York (1956).

61. Simmons, G.L., Capture Gamma Ray Currents Emerging from a Type

Concrete Barrier, as Functions of Incident Neutron Energy and

Pi recti on . Private Communication 6/11/86 in fulfillment of a contract with NBS.

62. Song, Y.T., A Semiempi rical Formula for Differential Dose Albedo

for Neutrons on Concrete , Report No. TN-589, Naval Civil

Engineering Laboratory, Port Hueneme, California (1964).

63. Song, Y.T., Fast-Neutron Streaming through Concrete Ducts;

Theoretical Calculations Using the Albedo Concept Compared with

Experimental Results , Technical Report No. R-565, Naval Civil

Engineering Laboratory, Port Hueneme, California (1968).

64. Song, Y.T., Calculating the Dose Variation of Fast-Neutrons

Streaming through an Iron Duct , Technical Report No. R-570,

Naval Civil Engineering Laboratory (1968).

65. Song, Y.T.; Huddleston, C.M.; and Chilton, A.B., Neutron Streaming

through Ducts, ANS. Trans. 38, 182 (1969).

42
66. Song, Y.T., Huddleston, C.M., and Chilton, A.B., Neutron Streaming

through Ducts--Analysi s of Experimental Nucl . Sci . Eng. 44 252-275

(1971).

' ;

67. Spencer, L.V., Diaz, J.A., and Moses, E., Neutron Penetration in

Cylindrical Ducts, NBS Report 8542, National Bureau of Standards,

(1964).

68. Spencer, L.V., and Diaz, J.A., Neutron Penetration in Ducts and

Other Cavities, NBS Report 8541, National Bureau of Standards

(1964).

69. Spencer, L.V., and Woolf, S., On Neutron Penetration in Duct

Systems with Large Cross Section , ANS-SD-7, p. 199, Proceedings

of Special Sessions on Gamma-Ray Production and Transport and on

Civil Defense Shielding, American Nuclear Society, Chicago,

Illinois (1967).

70. Terrell, C.W., Radiation Streaming in Shelter Ent ranceways , p. 266,

Proceedings of the Shielding Symposium, Oct. 31-Nov. 1, U.S.

Radiological Defense Laboratory, San Francisco, California (1960).

71. Terrell, C.W., Jerri, A.J., Lyday, R.O., and Sperber, D. , Radi ati on Streaming in Shelter Ent ranceways , Report No. ARF 1158- 12, Armour

Research Foundation, (1960).

72. Terrell, C.W. and Jerri, A.J., Radiation Streaming in Shelter

Entranceways , Report No. ARF 1158A01-5, Armour Research Foundation

(1961).

73. Terrell, C.W. and Jerri , A.J., Radiation Streaming in Ducts and

Shelter Entranceways , Report No. ARF 1158A02-7, Armour Research

Foundation (1962).

74. Vanzant, B.W. and DeHart R.C., Penetration of Nuclear into Shafts,

Tunnels and Entranceways , Southwest Research Institute, San

Antoni o, Texas (1966).

75. Vogt, H.G., Intercompari son of Monte Carlo Calculations and

Measurements of the Neutron and Gamma-Ray-Streaming in Multi Bent

Nucl. Reactor Shielding, p. 857, Science Press, Princeton, N.J.

(1977).

76. Woolf, S., private communication of unpublished computations to

L.V. Spencer (May 20, 1968).

77. Woolf, S., private communication of unpublished computations to

L.V. Spencer (Nov. 22, 1967).

78. Woolf, S. and Spencer, L. V., unpublished computations (1983).

43
APPENDIX A. Some Notes on Subprograms

Figure 1 identifies the subprograms of DCTDOS, as well as indicating

the flow of control. The list below may also be useful for quick

reference. Following the list are brief descriptions.

a) Control Subprogram:

DCTDOS

b) Data Subprograms : ~ ENTDAT (storage)

DATDCT (input data preparation and scaling

SRSDET (source and detector data preparation)

SND (preparation of source and detector weight lists and

arrays from basic neutron data)

SGD (preparation of weight lists and arrays from basic gamma

ray data)

WTLISTS ( storage of computed weights in a permanent file called

PREWTS, to avoid recomputation with each new run)

c) Duct Segment Computations : DCTISO (isotropic source entering a straight duct segment)

DCTSLT (mono-obliquity source entering a straight duct segment)

TURN (bend segment for different source angular parameters)

BTODR (straight segment following a bend)

ROOMX (duct input, detector pattern in a room)

d) Dose Computations : DOSX (dose computations for ROOMX, i.e., in a room)

DOSE (dose computations in, or at the end of a duct)

e) Auxiliary Subprograms : RETRN (makes estimate of backrefl ection from end of a closed

duct)

XTRAP (a short program for extrapolation of data, mainly for

straight duct segments)

ORDST (prepares approximate orders-of-refl ecti on

di stributions.

a) Control Subprogram ( DCTDOS ) . The pattern shown by figure 1 gives

a fairly complete summary of the options presented by DCTDOS. For each

type of case or problem there is a loop indexed by IC which evaluates a sequence of duct elements. At the conclusion of each problem or case there are options to enter a full set of new data (Run Case), new configuration data only (Problem Case), a change of a straight section

(Length), or a new source and/or detector (Weights).

In the File #1, orders 82-85 set up a possible selection of reflected

or unreflected component, although usually S=0 . Order 94 applies LPIPE

segment data to the final segment only. Order 118 gives a new application of weights previously applied in the same run. Order 119

44
applies new alternative weights near the end of a problem. Orders 124

and 148 develop weights only once for cases in which a segment length is

varied. Order 121 develops gamma weights only once if there are at least

two bends. The orders from 127 to 132 relate to the computation of dose

for the segment preceding a room after a ROOMX dose pattern computation,

for comparisons.

b) Data Subprograms . ENTDAT stores output data for elementary duct segments (QUCTi and

BEND arrays). It also stores selected spectral source and dose

information in most of the remaining lists. The arrays AWTS and AMOD

represent obliquity incidence data (sidewall and endwal 1 , respecti vely

)

for gamma ray sources. These are required by SRSDET to obtain realistic

albedo weights for gamma ray problems.

DATDCT is mainly a buffer program in which input describing the

configuration is properly scaled and inserted in COMMON storage. Input

is summarized in Table 4, and output of DATDCT is summarized in Table

A5 below. The ratio PEN, usually equated to zero, can be given a small

value to take account of the mean wall depth from which reflection

occurs. (We have made almost no use of this option).

The segment eccentricity ECC replaces the duct width and height, and

COF is the length scaled by the square root of the cross sectional

area. DELR mostly gives a cross sectional area ratio (ARAT) between

successive segments. But this may not change, and then 0THR=0 so that

DELR=1 (order 230). DELR is also used to provide a required parameter

value for IEL = 6 and 8. DZRO also is given the value of a required

parameter for IEL = 3 and 6, otherwise parameters of occasional

application only.

Table Al . Scaled output produced by DATDCT and recorded with I PR I NT = 2.

For IEL = 6, ARAT repeats the value for the preceding segment.

IEL COF ECC DELR DZRO SCL

SO PEN ZD00R RMCASE( I EL==8)

1 L//A W/H ARAT,1 ARAT /RAT

2 L//A W/H ARAT ,

ARAT /RAT

3 L//A W/H 1 F /RAT

4,5 L//A W/H ARAT,1 ARAT /RAT

6 W H L WD AEX

7 RDLNTH//A W/H ARAT,1 RDLNTH-.01 ARTO

8 WDP NW-1 NL-1 ARAT ARTO

45
All the quantities named in Table A1 have the same interpretation as in Table 4, for the same value for IEL. Note that COF contains mainly

scaled lengths, ECC has eccentricity ratios, DELR has reduction factors

for constriction.

SCL is a list (in COMMON storage) of scale factors. Little use has

been made of /RAT, which is ordinarily unity. ARTO is the ratio of

entrance area to exit area; if there has been a change, either single

or progressive, this scale-factor can indicate the magnitude of the

total effect. But the effect is computed segment by segment using ARAT

val ues.

A constriction confined to a segment entrance only is not really

accounted for in DCTDOS, although it would not require much change.

Note that an increase of segment cross section area will spread the exit

flux over a larger area, while a reduction will reduce the total

entering flux without any other effect on flux/(exit area).

SRSDET computes source and detector weights. It also makes use of

three other programs: WTLISTS provides for more permanent storage of the

weights, and SND and SGD apply multiple albedo data to the computation

of orders-of-refl ecti on data for neutrons and gamma rays.

Orders 274-278 provide a selection parameter for neutrons, capture

gammas, or neutrons plus capture gammas.

Orders from 280 through 327 are directed to computation of orders-ofreflection weights for the gamma-ray case, which depends strongly on the

configuration of the first leg and first bend. The estimate here is

based on use of simple formulae for angular limits for a source at the

centerline point of the entrance to a cylindrical duct (unreflected in

the first segment, the CUT limit) and the complementary CUTP angular

range which applies mainly to scattered components from the end-surface

of first leg. Note again that the unreflected from the first segment

(TVO) may be reflected following near-perpendi cul ar incidence from an endwall while the reflected component (TVP) exits the first segment

after reflection(s) somewhere in the first leg, here assumed cylindrical

for convenience.

The main task for SRSDET is the computation which takes account of

spectral changes with multiple reflections, and the effect which this

has on the multiple-reflection weights. This requires SND and SGD and

their processing of the large tables of multiple-reflection (albedo)

46
energy data. Hence we give below a more detailed analysis of SND and

SGD, which are the heart of SRSDET.

SND Purpose: to generate the arrays WTSD, SWTS, and the list

ENGAM, for neutrons and capture gammas.

Arguments, Arrays, and Lists : SGRUPS: Input source list, neutrons.

SWTS: This begins as the source list which, multiplied by a multi

reflection albedo energy array gives dose weight per detector

response group-strength.

One should note that WTSD and SWTS are compl ementary

quantities. After constructing each with the orders-ofreflection data, one must still multiply the SWTS by the dose

weights in DGRUPS, or correspondingly the WTSD by the spectral

weights in SGRUPS, to obtain the dose weights required (see

Eqs.9,9').

DGRUPS: Dose weights for the different spectral groups.

ENGAM: This is a list of dose weights for capture gamma rays generated

in the duct walls. It is assumed that the ratio of gamma

number-current albedo to neutron number-current albedo is

.095/. 75; and it is assumed that all the gammas are about 1 MeV, with dose weight 4.5E-10. The indicated proportion of

gammas/neutrons is assumed to hold for each incident order of

reflection, with all orders of reflection contributing,

independently of spectral energy. (These assumptions probably

are overly simplistic, are untested, and should be reexamined.)

ALBEDO: This is the list of emergent neutron currents for each order

of reflection (2nd index) for each energy group (first index),

for each initial energy (MEG value). These basic data are included into the program by the READ orders.

Sequence of Computations : 405-416: Preparatory values assigned.

417-421: Begin main loop, first zeroing ALBEDO, which may not always

be necessary.

422-426: Read the data and assign the unref 1 ected components.

427-432: Cumulate the contributions to both main arrays, plus ENGAM,

for all source energies.

47
SGO Purpose: for gammas, to generate the arrays WTSD and SWTS.

Arguments, Arrays, and Lists : All are as described in SND except that ALB replaces ALBEDO of SND,

and that VEL1 differs from ENGAM. (One must interpret SGRUPS, DGRUPS,

and DWTS, as well as ALB, as applying to gamma rays in this subroutine).

VEL1: This is a weight list for different obliquity angles with which

the initiating gamma rays may strike the wall. It is used

because the first order of reflection for gamma rays will be

significantly different, and will affect all other orders of

reflection; hence some account must be taken whether most gamma

rays first strike normally, grazingly, or with a cosine

distribution.

Sequence of Computations : Closely similar to SND.

451-455: Zero the arrays.

456-462: Begin the source energy loop. KE is the energy index and

LOB is the obliquity index. Note that the energies

proceed from top down, each represented by the five

obliquity cases. 463-467: Read the (mass storage) data and take care of the

unreflected component, assumed to be unity and to be

incident in proportion to VEL1 obliquity weights.

468-476: Cumulate the two main arrays over all source energies and

obliquities.

WTLISTS is called by DCTDOS at the beginning of a run to set up two

files, IXFILE and WTFILE. The latter contains precomputed weights for

different combinations of IDOSE, ISORS, and LPIPE, currently a maximum

of 65 sets. The former contains an index which states whether the set

is, in fact, stored for retrieval or should be computed. A second index

confirms the location in the file, and a third index gives the number of

weights. At the end of DCTDOS a final call for WTFILE stores the file

again in PREWTS including any new weight lists that may have been added.

This subroutine is called twice in SRSDET. The first call either

accesses precomputed weights or exits the routine because they have

not as yet been stored. The second call either exits the routine

because they already are included or files a newly prepared set of

weights in PREWTS.

48
c) Duct Segment Computations . These constitute the main workiny

set for solutions. The three basic programs, DCTISQ , DCTSLT , and TURN

follow essentially the same pattern, namely

1) interpolate required data from tables in ENTDAT,

2) construct one or more suitable backscatteri ng matrices for the

segment (see Eq. 2),

3) fold data for the new segment over the input to obtain the

output as in Eq. 2, and

4) modify , total (to VEL) and possibly print the results.

The program ROOMX is much more complex, because it does not simply

focus on transmission of a current, but rather attempts to construct a dose pattern. It treats separately the component entering the room

which does not reflect before contributing to the dose. (Let us call

this the unreflected component even though it may consist of many

earlier reflections.) Note that we have a mean angular parameter for

each to assist our description of these two components.

Also given special treatment is the once-reflected component, which

is sensitive to the angular and energy properties of the entering

current. All other orders of reflection are lumped together and

treated with the assumption that these contributions follow a cosine

emission law from the room surfaces. The number albedo (.75) is

assumed to be modified (reduced) by an escape probability which is

taken to be proportional to the ratio of entrance area to half the

total wall area.

We believe that these principles are reasonably sound. Tests

comparing results with simple neutron Monte Carlo data have been

promising. Correspondi ng tests for gamma ray data have not been

carried out; and patterns for gamma ray cases should be less accurate.

Other ways to estimate these detector patterns can also be devised and

could prove better.

BTODR computes penetration down straight segments which follow

bends, and represents the key to successful evaluation of mul ti -segment

ducts. It is a second-generation program.

49
Consider a bend from the vantage point of a position well -down the

following straight section. Mainly visible is the endwal 1 . While

albedo from the endwal 1 does not have a strictly isotropic flux upon

emergence, the flux does not greatly differ from an isotropic distribution for penetration purposes down a straight section to follow. We

approximate this contribution with our DCTISO subprogram. All three

other walls of the bend, including the endwal 1 of the first segment,

decrease in visibility by a cosine factor with distance down the next

segment. We approximate their contribution with our DCTSLT program with

the aid of our estimate of the angular parameters SO, SP of exit from

the bend (entrance into the new segment) for components unscattered (SO) or scattered in the bend. It is mainly the SP from the bend which must

be apportioned to the walls which become the new endwal 1 (for DCTISO)

and the new sidewalls. The SO applies only to one (DCTSLT) component

which happens to bypass the bend but usually strikes a wall near the new segment entrance.

The remaining complexities relate to estimates of the number of

neutrons (or gammas) which strike the different walls, both reflected

in the bend and unreflected, to apportion suitably their relative

contributions and mean angular incident at the new segment entrance.

Rough symmetry between the exit flux and that striking the endwal 1 -to-be

aids in making this estimate.

d) Dose Computation Subprograms. DOSX is complicated by the requirement to evaluate the dose for the pattern used in RQOMX . Otherwise it

applies the weight lists from SRSDET to the orders of reflection

distributions at the different specified positions.

Both DOSX and DOSE have incorporated several options, but at present

the only option used computes an estimated flux dose rather than current

dose (using Eq. 10). Both also yield a reduction factor relative to the

initial position at the duct entrance.

e) Auxiliary Subprograms . RETRN is of limited application, since it

affects the final result only in the general vicinity of a straight

segment endwal 1. It is not a simple program, since it must evaluate

what strikes the endwall and subsequent back-attenuation. Possibility

50
of reduction due to a missing part of the duct sidewalls is ignored.

XTRAP is a short subprogram which has taken different forms with

little change in results once the general mathematical form for

reduction, i.e., a power law down straight duct segments, is fixed.

Present formula is something of a hodgepodge which works after a fashi on

,

Y = (AL1 + -^r ) * 4f . (Al)

x u.iX

xi

If the parameters of the fit give AL1<0, a simpler form is used instead,

namely

Y = ( AL 1 / X

X

) (Al 1 )

In the usual case, the constant XI is pre-assigned a value on the

basis of y a x“

XI ~* 05

for large X, and AL2, AL1 are then computed.

Various methods for specifying the constant X_I_ have been applied. This

present approach is based on the assumption that a single term is

sufficient for the large-X trend. The constants AL1 and L2^ are then

determined to fit the distribution well for the largest tabulated values

for X, and give reasonable extrapolated values at distances beyond the

tabulation which are not excessively large. The limiting case is most

commonly Y+0 for X->-«>. But a finite limit for Y is not ruled out and can

be implemented.

ORDST stands for order-of-refl ection distribution. Its purpose is

to change a simple power-law distribution fixed by two parameters into

a modified distribution which is more flexible and realistic, but which

can be readily computed from the same two parameters. More precisely,

a power law can be expressed by a single constant times

a, a

2

, a

3

, a , . .. . ,

( A 2

)

The more flexible and applicable form which is developed in this

subprogram is an overall factor times

a, a

2 +ab, a3 +a2 b+ab 2

, (A2 1 )

51
Assignment of a value for the parameter b suggested by Woolf data for

elementary segments 76 leads to formulae which give _a and the multiplicative factor in terms of a. and the original strength constant. Such a change is rendered useful because orders-of-refl ection distributions in

ducts are not simple power laws, and much better descriptions can be

obtained with little more added flexibility than this can give.

APPENDIX B. Sources of Data

a) S. Woolf Duct Data . The unpublished data of ref. 76 constituted a large number of configurations and sources, as listed in Table B1 below.

Two types of directional sources were computed, mono-obliquity sources

and continuous sources. The latter were sampled from a normalized

distribution of the form

f(s,cose) = 2( 1-s )cos9*d(cose)/[( 1+s) 2 -4scose] 3/2 . (Bl)

The value s=0 corresponds to the case of incident isotropic flux.

Output orders-of-refl ection data were given for the zero'th through

the 20 ' th order of reflection. Albedo data were given similarly. To

properly include higher orders the Woolf output also gave values of the

Fourier Transform,

oo 0(p) = l 0(n)e~np

. ( B2)

n=I

These data turned out to be extremely useful in obtaining simple

analytic approximations to the output distributions-in-n.

Only a small portion of the Woolf data were directly used to obtain

the ENTDAT data for DUCT0, DUCT1, DUCT2, and BEND, from which the first

three parameters of Table 4 interpolate. And as noted in Section VI,

corrections were later necessary for an albedo of 0.75 instead of

0.696.

52
It was early very clear that angular parameters would be necessary in

the joining of different segments. But it was not so clear which or how

many to use. To attempt perhaps the simplest option first, additional

computations were made to provide SO and SP, the (angular average)

parameter values listed in the DUCTi and BEND arrays, for many of the

most important cases.

The Table B1 list of cases were mostly recomputed 77 in order to

provide data for <sece>, and thus for SO, SP.

Computations were also made for eccentric ducts with rectangular

cross sections. The eccentricity ratio e is used which equals

width/length. Endpoints plus 5 internal planes of a duct of length 6

diameters were computed. Eccentricity ratios n were .01, .25, .50,

.75, and 1.0. Each computation was performed for at least 6 s values.

Fortunately the sensitivity to this parameter turned out to be rather

low, except for extreme cases. This confirms also that square cross section data has fairly general applicability.

Table Bl . Most of the cases listed below were computed with 5,000 to

10,000 histories.

FINITE, STRAIGHT, OPEN DUCTS, CONTINUOUS ANGULAR INPUT

8 lengths (0 to 8 diameters), each with 9 input s values

FINITE, STRAIGHT, OPEN DUCTS, M0N0DIRECTI0NAL ANGULAR INPUT

8 lengths (0 to 6 diameters), each with input obliquities

cose.j

n =.01, .25, .50, and .75. (These were redone to 8 diameters

with 30,000 histories.)

SEMI-INFINITE, STRAIGHT DUCTS, CONTINUOUS ANGULAR INPUT

13 depths (0 to 6 diameters), each with 6 s values, computation

of both flux and current.

OPEN, STRAIGHT, FINITE DUCTS INTERNALLY, CONTINUOUS ANGULAR INPUT

Midpoints, 1, 2, and 4 diameter lengths, each for 6 s values.

0.5 diameters depth, 4.5 diameters length, each for 6 s values.

^ q

«. ii 7 0

11 11 11 11 11 11 11 CLOSED, STRAIGHT, FINITE DUCTS, CONTINUOUS ANGULAR INPUT

4, 6 diameters length, each for 6 s values

OPEN BEND SEGMENTS, SQUARE CROSS SECTIONS

lxl, 2x2, 3x2, 2x1.5, 1.5x2, 7 s values (continuous angular input)

" " " " "

, input perpendicul ar to entrance.

53
b) SA I Neutron Multipl e -Reflection Data . 38 These tabulated data are currently in floppy-disk files. They correspond to various

monoenergeti c neutron fluxes which are backscattered many times between

two plane-parallel concrete barriers. At their start, the angular

distribution is isotropic. With each backscattering , the emergent

neutrons are classified into energy groups at and lower than the source

energy, in a fairly standard set of 23 energy groups. The concrete

composition is that referred to as Type 04, and has density 2.34 g/cm 3

and composition by weight as shown in Table D3 below:

Table B2. Type 04 concrete composition by fractional weight.

ELEMENT FRACTIONAL WEIGHT ELEMENT FRACTIONAL WEIGHT

H .0056 Si .3151

0 .4979 S .0013

Na .0171 Ca .0829

Mg .0026 Fe .0124

A1 .0457

mul ti pi e-•albedo neutron data are called and used by SND, as previously noted.

c) Gamma Albedo Data . 78 These data were the basis for the multiplealbedo data required for SGD . Original computations by Woolf provided

detailed albedo data for single incidence only. A program was then

written to apply these data to the space between plane, semi -i nf i ni te

concrete barriers. Source gamma rays of 17 different energy groups were incident at each of the following obliquity cosines (relative to the

perpendi cul ar)

cos 0

Q = .05, .25, .50, .75, 1 .0 .

Albedo energy groups and emergent obliquity cosines were then applied a total of 10 times to evaluate the mul ti -reflection data in ALB.

The emergent obliquity cosines were eventually not applied to the

gamma ray cases, with reliance instead on a similarity between low

energy gamma and neutron albedo angular distributions.

54
Our approach to computations for the different duct seyment data

types requires the application of approximate albedo anyular distributions in specialized Monte Carlo studies, plus evaluation of emeryent

current-dose anyular distributions characterized by the SO and SP

parameters. These permit ready estimates of fluxes from currents by the

simple rule previously given of

Flux

fl ux current

* current

a* 2

1 + S * current .

( 10 )

Thus the emergent flux and flux/current ratio are the main angular data

recorded, and are used to link the various duct components together.

APPENDIX C. Remarks on Inputs and Outputs with Examples

Figure Cl gives both input and output data for a run with I PR I NT =3

To save space let us note that if we had used I PR I NT=2 the first page

would have been identical, but would have been followed by the last five

lines of the 2nd page of these output data to conclude the run. The second page gives source listing and spectrum at the detector, as well as a breakdown of tne DOSE value into contributions from different

orders of reflection, termed DFLXN. VELO and VELP differ in that the

former is unreflected in the final segment while the latter is reflected

in the final segment: both usually have multiple orders of reflection

in earl i er segments. (Note that the first column gives both an energy

group list and a list of orders of reflection (N=n+1).

The four VEL lists correspond to input and/or output for successive

duct segments. SO and SP values below give the angular distribution

parameters _S for the VELO and VELP components, respectively, which are,

however, recorded here only for the last VEL list. Note that for

forward currents, 1 < S < -1, with the value +1 correspond! ng to a current focused along the axis, and with the value -1 correspond! ng to

complete focus sideways , i.e., 90° from the forward direction. The

value .000 here in the IC=1 column has no meaning, being a default only.

55
INPUT35342 32033 1.000+000

3.000+000

0 *** M-M, 2-LEG, EPITH/TH ORDERS OF REFL., 9/2/86 ***

0 0.414 0.000 0.000 0.000 0.000

2 3.000 3.000 13.500 0.000 0.000

3 3.000 3.000 0.000 0.000 0.000

5 3.000 3.000 3.500 0.000 0.000

7 0.000 0.000 0.000 0.000 0.000

OUTPUT

$========$

ISORS32033 IDOSE ING LPIPE IEL DIMEN

l.OOOE+OO

IPRBLM 35342 NEL IPRINT JVARY JDEL DEL

3.000E+00

0 *** M-M, 2-LEG, EPITH/TH ORDERS OF REFL., 9/2/86 ***

IEL WIDTH-FT HEIGHT-FT LENGTH- FT OTHER SELECT

0 .414 .000 .000 .000 .000

2 3.000 3.000 13.500 .000 .000

3 3.000 3.000 .000 .000 .000

5 3.000 3.000 3.500 .000 .000

7 .000 .000 .000 .000 .000

IEL COF ECC DELR DZRO SCL

0 .41400 .00000 .00000 .00050 .00000

2 4.50000 1.00000 1.00000 .00000 1.00000

3 .00000 1.00000 1.00000 .00000 1.00000

5 1.16667 1.00000 1.00000 .00000 1.00000

7 .00000 1.00000 1.00000 -.01000 1.00000

NEXT CASE 1.167

DWTS//SWTD//ENGAM

4 . 124E+00 2.965E+00 2.252E+00 1.739E+00 1.320E+00 9.968E-01 7.457E-01

5.576E-01 4.178E-01 3 . 131E-01 2 . 334E-01 1.747E-01 1.302E-01 9 . 715E-02

7 . 214E-02 5.317E-02 3.968E-02 2.941E-02 2.207E-02 1 . 642E-02 1.214E-02

4 . 124E+00 2.965E+00 2.252E+00 1.739E+00 1.320E+00 9.968E-01 7.457E-01

5.576E-01 4. 178E-01 3. 131E-01 2.334E-01 1.747E-01 1 . 302E-01 9.715E-02

7.214E-02 5.317E-02 3.968E-02 2.941E-02 2.207E-02 1.642E-02 1.214E-02

.OOOE+OO 1.690E-10 1.283E-10 9.914E-11 7.523E-11 5.682E-11 4.250E-11

3.178E-11 2.381E-11 1.785E-11 1.331E-11 9.956E-12 7.421E-12 5.537E-12

4 . 112E-12 3.031E-12 2.262E-12 1.677E-12 1.258E-12 9.360E-13 6.919E-13

NEXT ( I ) = NEXT ( I

)

NEXT ( I

)

NEXT ( I ) =

0 DOSE=

2 DOSE=

3 DOSE=

5 DOSE=

5.833E+00

8.977E-02

4.889E-02

2.347E-02

REDUX FACTR= l.OOOE+OO

REDUX FACTR= 1.539E-02

REDUX FACTR= 8.381E-03

REDUX FACTR= 4.024E-03

Fig. Cl. Example of data from IPRINT=3.

56
Fig. Cl (Continued)

ARTO= l.OOOE+OO

MEG.N SOURCE SPCTRM DFLXN VELO VELP

1 €.OOOE+OO .OOOE+OO .OOOE+OO .OOOE+OO .OOOE+OO

2 1.OOOE+OO .OOOE+OO .OOOE+OO .OOOE+OO .OOOE+OO

3 2., 150E-02 2.010E-08 2.365E-04 5 .648E-05 5.605E-05

4 8,.000E-02 4.532E-07 4.912E-04 8.760E-05 1 .226E-04

5 1.405E-01 1.538E-06 6.543E-04 9.016E-05 1.640E-04

6 1., 341E-01 4.393E-06 7.069E-04 7.779E-05 1 .730E-04

7 2.. 67 3E-01 1.074E-05 6.775E-04 6 . 113E-05 1.594E-04

8 2.. 646E-01 1.986E-05 6.023E-04 4.550E-05 1.349E-04

9 9..440E-02 4.538E-06 5.095E-04 3.277E-05 1 .079E-04

10 3..038E-01 1.7 57E-05 4. 166E-04 2.312E-05 8.306E-05

11 4..430E-01 6.668E-05 3.327E-04 1 .610E-05 6.227E-05

12 3.. 922E-01 2.173E-04 2.613E-04 1 . 111E-05 4.581E-05

13 4.. 468E-01 2.341E-04 2.028E-04 7.630E-06 3.325E-05

14 2..464E-01 1.894E-04 1 . 561E-04 5.220E-06 2.389E-05

15 2.. 350E-01 4.362E-04 1.193E-04 3.563E-06 1 . 705E-05

16 2., 107E-01 4.200E-04 9.086E-05 2.429E-06 1.211E-05

17 2.. 099 E -01 4.961E-04 6.897E-05 1.654E-06 8.569E-06

18 1.. 492E-01 3.844E-04 5.223E-05 1 . 125E-06 6.048E-06

19 1.. 192E-01 3.246E-04 3.950E-05 7.655E-07 4.261E-06

20 1.. 497 E -0

4.463E-04 2.984E-05 5.206E-07 2.999E-06

21 1..202E-01 3.684E-04 1.609E-05 1.346E-07 1.586E-06

22 9., 570E-02 3.763E-04

23 «.OOOE+OO 1 . 946E-02

VEL(N)/IC=1 2 3 4 5

1 1., OOOE+OO .OOOE+OO .OOOE+OO .OOOE+OO

2 «.OOOE+OO 1 . 927 E -03 4.540E-04 .OOOE+OO

3 <.OOOE+OO 2. 163E-03 7.2 1 3E-04 1 . 125E-04

4 *.OOOE+OO 1.849E-03 7.729E-04 2.102E-04

5 t.OOOE+OO 1.424E-03 6.869E-04 2.541E-04

6 *.OOOE+OO 1.042E-03 5 .506E-04 2.508E-04

7 4.OOOE+OO 7.414E-04 4. 152E-04 2.205E-04

8 «.OOOE+OO 5.185E-04 3.016E-04 1.803E-04

.9 4.OOOE+OO 3.589E-04 2. 139E-04 1.406E-04

10 4.OOOE+OO 2.468E-04 1.495E-04 1.062E-04

11 4.OOOE+OO 1.690E-04 1.035E-04 7.837E-05

12 4.OOOE+OO 1 . 154E-04 7 . 114E-05 5.692E-05

13 4.OOOE+OO 7.870E-05 4.873E-05 4.088E-05

14 4.OOOE+OO 5.360E-05 3.328E-05 2.911E-05

15 4.OOOE+OO 3.648E-05 2.269E-05 2.062E-05

16 4.OOOE+OO 2.482E-05 1.546E-05 1.454E-05

17 .000E+00 1.688E-05 1 . 052E-05 1.022E-05

18 4.OOOE+OO 1.148E-05 7 .157E-06 7.173E-06

19 4.OOOE+OO 7.802E-06 4.868E-06 5.027E-06

20 4.OOOE+OO 5.304E-06 3.310E-06 3.519E-06

21 4.OOOE+OO 3.606E-06 2.250E-06 1 . 721E-06

S0= 4.. 140E-01 4 . 547 E -01 2.680E-01 3.602E-01

SP= 4.OOOE+OO 3. 532E-01 .OOOE+OO--7.022E-02

NEXT CASE 4.167

NEXT);i;1= 0 DOSE= 5.833E+00 REDUX FACTR= l.(100E+00

NEXTl i 1= 2 DOSE= 8.977E--02 REDUX FACTR= 1.!539E-02

NEXT! i 1= 3 DOSE= 4.889E--02 REDUX FACTR= 3.:381 E -03

NEXTl[i;1= 5 OOSE= 1.410E--03 REDUX FACTR= 2.418E-04

57
Fig. Cl (Continued)

ART0== 1.000E+00

MEG,N SOURCE SPCTRM DFLXN VELO VELP

1 .000E+00 .OOOE+OO .OOOE+OO .OOOE+OO .OOOE+OO

2 .OOOE+OO .OOOE+OO .OOOE+OO .OOOE+OO .OOOE+OO

3 2.150E-02 1.234E-09 1.472E-05 9 .337E-06 1 .242E-06

4 8.000E-02 2.641E-08 2. 598E-05 1.277E-05 3.863E-06

5 1.405E-01 9.033E-08 3.287E-05 1.236E-05 6 .445E-06

6 1.341E-01 2.550E-07 3.557E-05 1.031E-05 7.994E-06

7 2.673E-01 6.177E-07 3.51 IE-05 7.945E-06 8.410E-06

8 2.646E-01 1 . 131E-06 3.268E-05 5.842E-06 8.0Q3E-06

9 9.440E-02 2.576E-07 2.924E-05 4 . 176E-06 7.136E-06

10 3.038E-01 9.958E-07 2.547E-05 2.931E-06 6.088E-06

11 4.430E-01 3.728E-06 2.175E-05 2.035E-06 5.036E-06

12 3.922E-01 1 . 187E-05 1.832E-05 1.402E-06 4.073E-06

13 4.468E-01 1.267E-05 1.528E-05 9.613E-07 3.241E-06

14 2.464E-01 1.022E-05 1.264E-05 6.571E-07 2.546E-06

15 2.350E-01 2.334E-05 1.039E-05 4.483E-07 1.981E-06

16 2. 107E-01 2.24QE-05 8.498E-06 3.Q55E-07 1.529E-06

17 2.099E-01 2.645E-05 6.924E-06 2.079E-07 1.173E-06

18 1.492E-01 2.049E-05 5.624E-06 1.415E-07 8.950E-07

19 1.192E-01 1.731E-05 4.555E-06 9 .623E-08 6.801E-07

20 1.497E-01 2.378E-05 3.681E-06 6.544E-08 5. 149E-07

21 1.2Q2E-01 1.967E-05 2.116E-06 3 .207E-08 2.768E-07

22 9.570E-02 2.012E-05

23 •OOOE+OO 1.195E-03

VEL(N)/IC=1 2 3 4 5

1 1.000E+00 .OOOE+OO .OOOE+OO .OOOE+OO

2 .000E+00 1 . 927E-03 4 .540E-04 .OOOE+OO

3 .000E+00 2.163E-03 7.213E-04 1.058E-05

4 .000E+00 1.849E-03 7.729E-04 1.664E-05

5 .000E+00 1.424E-03 6.869E-04 1 . 88 IE -05

6 .000E+00 1.042E-03 5.506E-04 1.831E-05

7 .OOOE+OO 7.414E-04 4. 152E-04 1.636E-05

8 .000E+00 5 . 185E-04 3 .016E-04 1.385E-05

9 .000E+00 3.589E-04 2. 139E-04 1 . 131E-05

10 .OOOE+OO 2.468E-04 1.495E-04 9.019E-06

11 .OOOE+OO 1.690E-04 1.035E-04 7 .070E-06

12 .OOOE+OO 1 . 154E-04 7 . 1 14E-05 5.475E-06

13 .OOOE+OO 7.870E-05 4.873E-05 4. 202E-06

14 .OOOE+OO 5.360E-05 3.328E-05 3.203E-06

15 .OOOE+OO 3.648E-05 2.269E-05 2.429E-06

16 .OOOE+OO 2.482E-05 1.546E-05 1.834E-06

17 .OOOE+OO 1.688E-05 1 . 052E-05 1.381E-06

18 .OOOE+OO 1 . 148E-05 7.157E-06 1.037E-06

19 .OOOE+OO 7.802E-06 4.868E-06 7 . 763E-07

20 .OOOE+OO 5 .304E-06 3 . 31 OE -06 5.804E-07

21 .OOOE+OO 3.606E-06 2.250E-06 3.088E-07

S0= 4 . 140E-01 4.547E-01 2.680E-01 7.438E-01

SP= .OOOE+OO 3.532E-01 .OOOE+OO 2.636E-01

58
The second case (4.167) follows the same printout pattern.

Figure C2 reproduces the input for all the ORNL neutron duct

penetration sources of Figs. 6-9. Only a single run was requi red

because different sources and detectors were cycled at each detector

point. On an IBM XT, with the 8087 chip and storage of the ALBEDO,

ALB, PREWTS files in the RAM memory, this run required about 5 minutes.

It should be noted that the (unit) thermal neutron source required

additional normalizing of the data by a factor of very nearly 2 (neutrons/cm2

) to obtain the experimental source strength. With use of

IPRINT = 1 5 the many cases require about 9 pages of output.

Figure C3 shows the input for the ARF C060 gamma ray experiments.

Only a single detector response has been included. Note that in this case the first input data card suffices for 1-leg, 2-leg, and 3-leg cases. Time required was a little over 4 minutes. In every computation about 5/3 minutes is required at the start to set up the PREWTS

data file into the arrays consulted in SRSDET. Actual computation

time of problems can be better estimated by subtracting this setup

interval

The references include many other interesting cases, including

design data such as are in Refs. 3 and 4, additional experimental data

from many of the references listed, Monte Carlo studies, and the

various semi -analytic and semi -empi ri cal studies. Note that concretewalled cavities are presently required, or in the case of gamma rays

moderately low-Z walls. The program is also written only for near- right-angle bends. Because of the similarity of the 137Cs ARF results

to the 60 Co results we have not included Cesium source computations

in this report.

Preliminary applications have also been made to some duct configurations for possible Civil Defense structures. These have not been

included in this report because comparison data is lacking.

69
INPUT

2 2 0 3 1 1.000+000

2

•5 3 I 2 11 1.333+000

2 3 4 0 3 1.000+000

2 3 3 0 3 3.600+005

2 3 3 1 3 3.600+005

***MAERCKER-MUCKEN, 1-LEG, THERMAL 1MEUTS, 8/31/86***

0 0.414 0 e 000 0.000 0.000 0.000

2 3.000 3.000 0.000 0.000 0.000

7 0.000 0.000 0.000 0.000 0.000

2 2 0 3 1 1.000+000

2

o 5 1 4 8 1.333+000

•32 3 4 0 3 1.000+000

2 3 3 0 3 3.600+005

2 3 3 1 3 3.600+005

***M+M, 2-LEG, THERMAL FLUX, 8-15-133 ***

0 0.414 0.000 0.000 0.000 0.000

2 3.000 3.000 13.500 0.000 0.000

3 3.000 3.000 0.000 0.000 0.000

5 3.000 3.000 0.500 0.000 0.000

7 0.000 0.000 0.000 0.000 0.000

2 2 0 3 1 1.000+000

3

r> 7 1 6 6 1.333+000

>}3 3 4 0 3 1 . 000+000

3 3 3 0 3 3.600+005

3 3 3 1 3 3.600+005

***MAERCKER-MUCKEN, 3-LEGi, THERMAL NEUTS, 7/26/83[ ***

0 0.414 0.000 0.000 0.000 0.000

2 3.000 3.000 13.500 0.000 0.000

3 3.000 3.000 0.000 0.000 0.000

5 3.000 3.000 12.000 0.000 0.000

3 3.000 3.000 0.000 0.000 0.000

5 3.000 3.000 0.500 0.000 0.000

7 0.000 0.000 0.000 0.000 0.000

Fig* C2* Input Data for computations of the ORNL duct.

60
INPUT55031 3.319+009

1 3 1 2 12 3.000-001

0 *** ARF 1-LEG, 50 CURIE C0-60 , RAD DOSE, 9/1/86 ***

0 0.000 0.000 0.000 0.000 0.000

1 6.000 6,000 1.800 0.000 0.000

7 -1

0 *** ARF

0.000

5 1 0.000 0.000

4 10 3.000-001

0.000 0.000

2-LEG, 50 CURIE C0-60 , RAD DOSE, 9/1/86 ***

0 0.000 0.000 0.000 0.000 0.000

1 6.000 6.000 10.000 0.000 0.000

3 6.000 6.000 0.000 0.000 0.000

5 6.000 6.000 1.000 0.000 0.000

7 -3

0 *** ARF

0.000

7 1 0.000 0.0006 6 3.000-001

0.000 0.000

3-LEG, 50 CURIE C0-60 , RAD DOSE, 9/1/86 ***

0 0.000 0.000 0.000 0.000 0.000

1 6.000 6.000 10.000 0.000 0.000

3 6.000 6.000 0.000 0.000 0.000

5 6.000 6.000 8.000 0.000 0.000

3 6.000 6.000 0.000 0.000 0.000

5 6.000 6.000 1.000 0.000 0.000

7 0.000 0.000 0.000 0.000 0.000

Fig. C3. Input Data for computations of the ARF duct.

61
APPENDIX D. FORTRAN Listing of DCTDQS

The listing of DCTDOS subprograms which follows can be assembled

and run using Microsoft FORTRAN software with a microcomputer . It also requires three additional files, called NBS-TD-G, NBS-TD-N, and

PREWTS. The first of these file contains gamma ray multi-reflection

albedo data, the second contains neutron mul ti -refl ection data, and the

third must have a format acceptable to the subroutine WTLISTS, but is

initially full of zeros. All three files are between 45,000 and 56,000

bytes in size. Running times earlier quoted did not include preparatory

storage of the files in RAM memory.

For assembly purposes we have divided DCTDOS '’nto four parts roughly

equal in size because of a 64,000 byte memory limit. After assembly of

each part the object programs have been linked into a run program of

size roughly 180,000 bytes.

The subprograms have been little affected by the change from a mainframe to a microcomputer with the exception of ENTDAT, which

required a special format.

62
1

i 2:

3:4:

5:6:

7:

8:9: 10 : 11 : 12 : 13:

14:

15:

16:

17:

18:

19:

20: 21 : 22 : 23:

24:

25:

26:

27:

28: 29:

30:

31:

32:

33: 34:

35:

36:

37:

38:

39: 40:

41:

42:

43: 44:

45:

46:

47:

48:

49:

50: 51:

52:

53:

54:

55:

56: 57:

58:

59: 60:

61:

62:

63:

64:

65: 66 : 67: 68 : 69:

70:

71:

72:

73:

74:

75:

76:

77:

FILE 1

DCTDOS IS TO CALCULATE DOSE IN A DUCT-BEND-ROOM CONFIGURATION

COMMON DIAMS ( 18 ) , DUCTO ( 3 . 17 . 6 ) . DUCT1 ( 3 . 17 , 6 ) , DUCT2 (3,17,6) .IEL.LP. X BEND (10.6). SGRUPS ( 23 ) , DWTS ( 21 ) . CRRNTO (9.10). RCRRNT (9,10)

COMMON RVCT ( 21 ) , RMTX ( 21 , 21 ) . SO ( 9 ) . SP ( 9 ) . DZRO ( 9 ) , NEXT ( 9 ) . COF ( 9 )

X ECC ( 9 ) , DELR ( 9 ) ,SNYNEU(23) .SINVAL(IO) , COl , C02 , VEL (9,21) , IPRINT COMMON RCOSO . RCOSP , RCOS . DGRUPS ( 23 ) , VELO (9,21). VELP (9,21) X , VEL1 ( 5 ) .SCALE. SWTS (23. 21) ,WTSD(23.21)

COMMON IDOSE , HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) , SINR ( 23 ) . THMLN ( 23

)

X, ENGAM ( 21 ) , DIMEN, SCOBLT( 23) ,AWTS(6,5) ,SGAM(23) ,DGAM(23) ,AMOD(5,5)

DIMENSION SVELO (22.20), SVELP (22.20). SCL ( 9 ) , SELECT ( 9

)

DIMENSION IPR ( 6 ) . ISO ( 6 ) , IDO ( 6 ) . NGI ( 6 ) . ILP ( 6 ) . DIMI ( 6

)

1 FORMAT ( // ///20X, 10H$========$)

66 FORMAT (13,13,13,13,13, 1PE10 . 3 , 30H ( IPR, ISO . IDO , ING . ILP . DIMEN )

)

2 FORMAT (F10. 3)

3 FORMAT ( F10 . 3 , 7D13 . 3

)

4 FORMAT ( 7 ( 1PE10 . 3 )

)

5 FORMAT (1H ) 6 FORMAT (15, 5F10 . 5) 7 FORMAT (1H1)

8 FORMAT ( 2 5H NEXT CASE.F7.3)

9 FORMAT ( 515 , 1PE10 . 3

)

10 FORMAT ( 55H IPRBLM NEL IPRINT JVARY JDEL DEL ) 11 FORMAT (55H ISORS IDOSE ING LPIPE I EL DIMEN ) 12 FORMAT (55H I EL COF ECC DELR DZRO SCL ) OPEN ( 10 , FILE= ' B : NBS-DT-G '

)

OPEN ( 11 , FILE= ’ B : NBS-DT-N '

)

OPEN ( 05 , FILE= ' ') OPEN ( 06 , FILE= ' A : PAPER '

)

OPEN ( 07 , FILE= ' B : PREWTS '

)

IOX=-l

CALL WTLISTS ( JSTT . JXWT , IOX . NMAX . DWTS , SWTD . ENGAM , IDOSE . ISORS . LP

)

CALL ENTDAT

MODG=0

M0DN=0

13 READ (5, 9) ISORS , IDOSE . ING . LPIPE . IEL . DIMEN IPRINT=IEL

IF ( IPRINT. GT.O) WRITE (6. 7) 16 IDEL=0 READ (5.9) IPRBLM , NEL . IPRINT . JVARY , JDEL , DEL

IF ( IPRINT. EQ.O) GO TO 50 WRITE (6,1)

WRITE (6. 11)

WRITE (6.9) ISORS , IDOSE , ING . LPIPE , IEL . DIMEN WRITE (6. 10)

WRITE (6.9) IPRBLM . NEL . IPRINT . JVARY . JDEL , DEL

50 READ (5. 9) JDOS

WRITE ( 6 , 9 ) JDOS

IDOS=0

IF ( JDOS . LT . 1 ) GO TO 68

DO 67 1=1. JDOS

READ (5,9) IPR ( I ) , ISO ( I ) , IDO ( I ) . NGI ( I ) , ILP ( I ) , DIMI ( I

)

IF (IPRINT. EQ.O) GO TO 69

WRITE (6.9) IPR ( I ) , ISO ( I ) , IDO ( I ) , NGI ( I ) , ILP ( I ) . DIMI ( I

)

69 IPR ( JDOS+1 ) = IPRBLM

ISO( JDOS+1) =ISORS

IDO ( JDOS+1 )=IDOSE

NGI ( JDOS+1 )=ING ILP ( JDOS+1 )=LPIPE

DIMI( JDOS+1) =DIMEN

67 CONTINUE

68 CALL DATDCT(NEL. SCL. SELECT)

IF (IPRINT.LT. 2) GO TO 14 WRITE (6,12)

WRITE (6,6) (NEXT ( I ) ,COF(I) .ECC(I) , DELR(I) ,DZRO(I) .SCL(I) ,1=1, NEL)

14 IC=1

IF (IPRINT. EQ.O) GO TO 51 WRITE (6,8) COF ( JVARY

)

51 CONTINUE

VELO (IC, 1) =1

VEL ( IC , 1

) =1

VELP ( IC , 1

) =0

RVCT ( 1 ) =1

SO(IC) =COF(l)

DO 15 1=2.21

VELO ( IC , I

) =0

VELP(IC.I) =0. 63
7879

80

8182

83

8485

8687

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

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

VEL ( IC , I

) =0

.

15 RVCT ( I ) =0

.

30 IF (NEL-IC) 1000.1000.31

31 IC=IC+1

S= SELECT (1C)

IF(S.EQ.O.O) GO TO 33

V0=1.+0.5*S*(1.-S)

VP=1.-0.5*S*(1.+S)

DO 40 N=1 , 21 VELO ( IC-1 , N) =VO*VELO (IC-l.N)

40 VELP ( IC-1. N)=VP*VELP( IC-l.N)

33 IEL=NEXT ( IC

)

SCALE=SCL ( IC) LP = 2 RDLNTH = DZRO ( NEL

)

IF ( IC+l-NEL) 35,32,32

32 LP=LPIPE

C 35 GO TO (100.110,300,500,500,400,200,210) , IEL

C 100 DLNGTH=COF ( IC)

CALL DCTISO ( DLNGTH , IC)

GO TO 30

C 110 DLNGTH=COF ( IC

)

CALL DCTSLT ( DLNGTH, IC)

GO TO 30 C 300 CALL TURN ( SEN, IC)

GO TO 30

C 400 CALL ROOMX ( IC , SVELO , SVELP

)

GO TO 30

C 500 DLNGTH=COF ( IC

)

CALL BTODR ( DLNGTH. IC)

GO TO 30

C 200 IEL=NEL

IF (NEXT ( IC-1 ) .EQ.8) GO TO 204 IF ( IPRBLM . LT , 0 ) GO TO 205

IF( JDOS . GE. 1) GO TO 201 IF ( IDOSE . NE . 5 ) GO TO 202

IF (NEL. GE. 5) GO TO 202

201 CALL SRSDET ( ISORS , ING , MODG , MODN

)

GO TO 205

202 IF(IDEL.LT.l) CALL SRSDET ( ISORS , ING , MODG . MODN

)

GO TO 205

204 IF ( IDOS . GT . 0 ) GO TO 205

IC=IC-2

COF ( IC) =COF ( IC+2

)

SCL(IC) =SCL(IC+2)

DELR(IC) =DELR ( IC+2

)

ECC(IC) =ECC ( IC+2

)

DZRO ( IC ) =DZRO ( IC+2

)

205 CALL DOSE ( IC

)

IF( JDOS.LT.l) GO TO 220

IDOS=IDOS+l

IPRBLM=IPR ( IDOS

)

ISORS=ISO ( IDOS

)

IDOSE=IDO ( IDOS

)

ING=NGI ( IDOS

)

LP=ILP ( IDOS

)

DIMEN=DIMI ( IDOS

)

206 IF ( IDOS. GT. JDOS) GO TO 220

IF ( IPRINT . EQ . 0 ) GO TO 207

WRITE (6,66) IPRBLM . ISORS . IDOSE . ING , LP , DIMEN

207 IF (NEXT (IC-1) .EQ.8) GO TO 205

GO TO 200

210 IEL=NEL

IF(IDEL.LT.l) CALL SRSDET ( ISORS . ING , MODG . MODN)

CALL DOSX(IC, SVELO, SVELP)

IF (JDOS.LT.l) GO TO 212

IDOS=IDOS+l

IPRBLM=IPR ( IDOS

)

ISORS=ISO ( IDOS

)

IDOSE = IDO ( IDOS

)

ING=NGI ( IDOS

)

LP=ILP ( IDOS

)

DIMEN=DIMI (IDOS)

IF ( IPRINT. EQ.O) GO TO 211

IF ( IDOS . GT . JDOS ) GO TO 212

WRITE (6.66) IPRBLM . ISORS . IDOSE . ING . LP . DIMEN
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

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

FILE 1 Continued

211 IF ( IDOS . LE . JDOS ) GO TO 210

212 IF(IC.EQ.NEL) GO TO 240

IF(IC.LT.NEL) IC=IC+1

IF(IC.EQ.NEL) GO TO 233

220 IDEL=IDEL+1

IF (IDEL-JDEL) 230.240.240

230 COF ( JVARY ) = COF ( JVARY ) + DEL

IF ( JDOS . LT . 1 ) GO TO 14 233 IDOS=0

IF (NEXT ( IC-1 ) .EQ.8)GO TO 33

GO TO 14 240 IF ( IPRBLM**2-4 ) 16.13.1000

1000 IOX=l

CALL WTLISTS ( JSTT . JSWT , IOX . NMAX . DWTS . SWTD . ENGAM . I DOSE , ISORS . LP

)

STOP

END SUBROUTINE DATDCT ( NEL . SCL , SELECT

)

COMMON DIAMS ( 18 ) , DUCTO ( 3 . 17 , 6 ) , DUCT1 ( 3 , 17 . 6 ) , DUCT2 (3,17,6) .IEL.LP, X BEND (10.6). SGRUPS ( 23 ) . DWTS ( 21 ) . CRRNTO (9,10). RCRRNT (9.10)

COMMON RVCT ( 21 ) . RMTX ( 21 , 21 ) . SO ( 9 ) . SP ( 9 ) . DZRO ( 9 ) . NEXT ( 9 ) . COF ( 9 )

.

X ECC ( 9 ) , DELR ( 9 ) ,SNYNEU(23) .SINVAL(IO) , COl , C02 , VEL (9,21) . IPRINT COMMON RCOSO . RCOSP , RCOS . DGRUPS ( 23 ) , VELO (9,21), VELP (9.21) X , VEL1 ( 5 ) .SCALE, SWTS (23, 21) ,WTSD(23,21)

COMMON I DOSE , HNDRSN ( 23 ) . SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) . THMLN ( 23

)

X , ENGAM ( 21 ) , DIMEN , SCOBLT ( 23 ) . AWTS (6.5), SGAM ( 23 ) . DGAM ( 23 ) , AMOD (5.5)

DIMENSION SCL ( 9 ) . HGHT ( 9 ) , WDTH ( 9 ) . XLTH ( 9 ) . AREA ( 9 ) , OTHR ( 9 ) , SELECT ( 9

)

CHARACTER* 80 TITLE

2 FORMAT (F10. 3)

3 FORMAT ( F10 . 3 , 7D13 . 3

)

4 FORMAT ( 7 ( 1PE10 . 3 )

)

5 FORMAT (1H ) 6 FORMAT (15. 5F10.3)

7 FORMAT (1H1)

9

FORMAT ( 515, 2F5. 2, 15)

10 FORMAT (55H I EL WIDTH-FT HEIGHT-FT LENGTH-FT OTHER SELECT ) 1 FORMAT (A80) READ (5,1) TITLE READ (5. 6) (NEXT ( IC) .WDTH(IC) , HGHT(IC) .XLTH(IC) .OTHR(IC)

,

* SELECT ( IC ) ,IC=1. NEL)

IF ( IPRINT . EQ . 0 ) GO TO 500

WRITE (6, 5) WRITE (6,1) TITLE

WRITE (6. 5) WRITE (6. 10)

WRITE (6. 6) ( NEXT ( IC ) .WDTH(IC) , HGHT(IC) .XLTH(IC) .OTHR(IC)

.

* SELECT ( IC ) ,IC=1. NEL)

500 CONTINUE

COF ( 1

) =WDTH ( 1

)

PEN=HGHT ( 1

)

ZDOOR=XLTH ( 1

)

53 FORMAT ( 55H l=DOOR, PEN WALL (FT), DOOR THICKNESS (FT OF CONCRETE) ) IF ( IPRINT . LT . 4 ) GO TO 501

WRITE (6. 5) WRITE (6, 53)

WRITE (6, 6) NEXT(l) .PEN. ZDOOR

WRITE (6, 5) 501 CONTINUE

ECC ( 1 ) =PEN

DELR(l) = ZDOOR AREA ( 1 ) = 0

.

DO 25 IC=2 , NEL IEL=NEXT ( IC) IF ( IEL-7) 11,30,33

11 AREA ( IC) =WDTH ( IC) *HGHT ( IC

)

RAT= (2. *PEN+WDTH(IC) ) * ( 2

. *PEN+HGHT ( IC) ) /AREA(IC)

SCL ( IC ) =SQRT ( RAT)

ECC(IC) =WDTH(IC) /HGHT(IC)

DZRO ( IC) =OTHR ( IC

)

DELR( IC) =1

.

IF ( IEL . EQ . 3 ) GO TO 13 DELR(IC) =OTHR(IC)

IF ( OTHR ( IC ) . EQ . 0 . ) DELR(IC)=1.

IF ( IEL . EQ . 6 ) GO TO 22

13 COF ( IC )=XLTH(IC)/SQRT( AREA (IC)

)

GO TO 25

22 COF ( IC )= WDTH (IC)

ECC ( IC ) =HGHT ( IC

)

DELR ( IC) =XLTH ( IC

)

DZRO ( IC ) =OTHR ( IC

)

SCL(IC) =OTHR ( IC+1

)

GO TO 25

30 ICP=IC-1

65
244

245

246

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

288

289

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

319

320

321

322

323

324

325

FILE 1 IF (NEXT ( IC-l ) . EQ . 8 ) ICP=IC-3

COF ( IC) =XLTH ( IC ) /SQRT (AREA ( ICP )

)

SCL ( IC ) = AREA ( 2 ) /AREA ( ICP

)

IF (NEXT ( IC-l ) . EQ . 8 ) SCL ( IC ) =SCL ( IC-l

)

DELR ( IC) =DELR ( ICP

)

ECC ( IC) =ECC ( ICP)

DZRO ( IC) =XLTH ( IC) - . 01

GO TO 25

33 COF ( IC ) = WDTH ( IC

)

ECC(IC)=HGHT(IC)

DELR ( IC ) =XLTH ( IC

)

DZRO ( IC ) =OTHR ( IC ) /AREA ( IC- 2

)

SCL ( IC ) = AREA ( 2 ) /OTHR ( IC

)

DZRO ( 1

) =OTHR ( 1

)

25 CONTINUE

1000 RETURN

END SUBROUTINE SRSDET ( ISORS . ING . MODG , MODN

)

COMMON DIAMS ( 18 ) . DUCTO ( 3 . 17 . 6 ) .00011(3,17,6) , DUCT2 (3,17,6) .IEL.LP X BEND (10,6) , SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO (9,10) , RCRRNT (9,10)

COMMON RVCT ( 21 ) . RMTX ( 21 . 21 ) , SO ( 9 ) , SP ( 9 ) , DZRO ( 9 ) , NEXT ( 9 ) , COF ( 9 )

,

X ECC ( 9 ) , DELR ( 9 ) ,SNYNEU(23) ,SINVAL(10) , C01 , C02 , VEL (9,21) , IPRINT COMMON RCOSO , RCOSP , RCOS , DGRUPS ( 23 ) , VELO (9,21). VELP (9,21) X ,VEL1(5) .SCALE. SWTS(23. 21) ,WTSD(23,21)

COMMON IDOSE . HNDRSN ( 23 ) , SEPI ( 23 ) . STH ( 23 ) , SINR ( 23 ) , THMLN ( 23

)

X , ENGAM ( 21 ) , DIMEN , SCOBLT ( 23 ) . AWTS (6.5), SGAM ( 23 ) , DGAM ( 23 ) , AMOD (5,5)

DIMENSION SWTD ( 21 ) , TEMP ( 5

)

1 FORMAT (1H ) 2 FORMAT ( 5 5H DWTS //SWTD/ /ENGAM

3 FORMAT ( 7 ( 1PE10 . 3 )

)

CNG=ING

CG=CNG**2

CN=1 . + (CNG-CG) /2

.

CN=CN*DIMEN

CG=CG* DIMEN

JSTT=0

IOX=0

CALL WTLISTS ( JSTT . JXWT . IOX . NMAX . DWTS . SWTD , ENGAM , IDOSE . ISORS , LP

)

IF ( IPRINT, GT. 2) GO TO 100

IF ( JSTT . EQ . 0 ) GO TO 15

50 IF ( IDOSE.LT. 5) GO TO 100

IAG=IS0RS-3

IF(IAG.LT.l) GO TO 100

XL=COF ( 2) CUT=0 . 5/SQRT (XL**2+0 , 25

)

CUTP=SQRT ( 1 . -CUT**2

)

DO 97 J=1 ,

5

AMOD (3, J)=0.

97 AMOD ( 4 , J) =0

.

DO 98 J=1 ,

5

FJ=

J

IF ( CUT . GT . AMOD ( 1 , J ) ) GO TO 99 98 CONTINUE

99 NCUT=FJ+ .

1

DO 980 J=1 .

5

FJ= J IF ( CUTP GT . AMOD ( 1 , J ) ) GO TO 990

980 CONTINUE

990 NCUTP=FJ+ .

1

FR= ( AMOD ( 2 , NCUT ) -CUT ) / ( AMOD ( 2 . NCUT ) -AMOD ( 1 , NCUT )

)

FRP= ( AMOD ( 2 . NCUTP ) -CUTP ) / ( AMOD ( 2 , NCUTP ) -AMOD ( 1 . NCUTP )

)

DO 96 J=1 ,NCUT

96 AMOD ( 3 , J) =AWTS ( IAG. J) DO 960 J=l, NCUTP

960 AMOD ( 4 , J ) = AWTS ( IAG+ 2 . J

)

AMOD ( 3 , NCUT ) =FR*AMOD ( 3 . NCUT

)

AMOD ( 4 . NCUTP ) =FRP* AMOD ( 4 . NCUTP

)

SD = 0

.

ED=0

.

DO 102 J=1 ,

5

SD=SD+AMOD ( 3 . J)

102 ED=ED+AMOD (4 . J) DO 103 J=1 ,

5

AMOD ( 3 . J) = AMOD (3,J)/(SD+l.E-7)

103 AMOD ( 4 , J) =AMOD ( 4 , J) / ( ED+1 . E-7

)

TOT=0

.

TVP=VELP (2,2) + VELP (2.3) + VELP (2,4) +VELP (2.5) +VELP (2,6) +VELP (2.7)

TVO=VELO (2,1) IF ( IEL . LT . 4 ) TVO=0

.

DO 200 J=1 ,

5 Continued

66
326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

FILE 1 Continued

TEMP ( J

) =TVO*AMOD ( 4 . J) +TVP*AMOD ( 3 . J)

200 TOT = TOT+ TEMP ( J

)

DO 201 J=1 ,

5

201 VEL1( J) =TEMP(6-J) /TOT

IF (IPRINT.LT. 4) GO TO 100

WRITE (6,1)

WRITE (6. 2) WRITE (6,3) ( AMOD (3,J) ,J=1,5)

WRITE (6,3) ( AMOD (4,J) ,J=1,5)

WRITE (6. 3) TVO , TVP , TOT . TOT , (VEL1(J) , J=l,5)

WRITE (6, 2) 498 WRITE (6,1)

100 DO 12 MEG=1 , 23

GO TO (105,110,115,120,125) .ISORS

105 SGRUPS (MEG) =SINR (MEG)

GO TO 4 110 SGRUPS (MEG) =STH( MEG)

GO TO 4 115 SGRUPS (MEG) =SEPI (MEG)

GO TO 4 120 SGRUPS (MEG) =SGAM( MEG) GO TO 4 125 SGRUPS ( MEG )=SCOBLT( MEG)

4 GO TO (205.210,215,220,225) , IDOSE

205 DGRUPS ( MEG ) = SNYNEU ( MEG

)

GO TO 12 210 DGRUPS (MEG) =1. GO TO 12 215 DGRUPS (MEG) =HNDRSN( MEG)

GO TO 12 220 DGRUPS ( MEG ) =THMLN( MEG)

GO TO 12 225 DGRUPS (MEG) =DGAM( MEG)

12 CONTINUE

IF ( IDOSE . LT . 5 ) GO TO 13 NMAX=10

MMAX=17

DO 14 N=1 , 21 DWTS (N) =0

.

SWTD ( N) =0

14 ENGAM (N) =0

.

CALL SGD ( SGRUPS , SWTS , WTSD , DGRUPS , VEL1 , MODG

)

GO TO 11 13 CALL SND( SGRUPS, SWTS, WTSD. DGRUPS, ENGAM, MODN)

NMAX=21

MMAX=23

11 IF ( JSTT . EQ . 0 ) GO TO 15

DO 19 N=1,NMAX

DWTS (N) =0

.

SWTD ( N ) = 0

.

DO 19 M=1 , MMAX

SWTD(N) =SWTD (N) +SWTS(M,N) * DGRUPS (M) DWTS ( N ) =DWTS ( N ) + WTSD ( M . N

) * SGRUPS ( M

)

19 CONTINUE

15 CALL WTLISTS ( JSTT, JXWT, IOX. NMAX. DWTS, SWTD, ENGAM, IDOSE, ISORS, LP) DO 20 N=1 , NMAX

IF ( IDOSE. EQ. 5) DWTS (N) =SWTD (N) DWTS (N) = CN* DWTS (N) + CG* ENGAM ( N

)

20 CONTINUE

C 25 IF (IPRINT.LT. 2) GO TO 1000

WRITE (6.1)

WRITE (6. 2) WRITE (6,3) (DWTS(N) ,N=1.NMAX)

WRITE (6.1)

WRITE (6,3) (SWTD(N) ,N=1,NMAX)

WRITE (6.1)

WRITE (6,3) ( ENGAM (N) ,N=1,NMAX)

WRITE (6. 1) 1000 RETURN

END SUBROUTINE SND ( SGRUPS , SWTS . WTSD . DGRUPS . ENGAM . MODN

)

SLARGE ALBDO

DIMENSION SGRUPS (23) .SWTS (23. 21) .WTSD (23. 21) .DGRUPS (23)

.

* ALBDO (23. 21. 23) .ENGAM (21)

CHARACTER* 80 TITLE

1 FORMAT (A80)

7010 FORMAT ( 15

)

7020 FORMAT ( 8 ( 1PE10 . 3 )

)

IF ( MODN . EQ . 1 ) GO TO 12 IUT=11

REWIND IUT READ ( IUT . 1 ) TITLE

fi7
FILE 1 Continued

409

410

411

412

413

414

415

416

417

418

419

420

423 422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

IF (IPRINT.LT. 4) GO TO 12 WRITE (6.1) TITLE

12 FG=. 095/. 75*4. 5E-10

DO 25 N=1 , 21 ENGAM (N) =0

.

DO 25 M=1 » 23

SWTS(M.N) =0. 25 WTSD (M , N) =0

.

DO 350 MEG=1 , 23

IF(MODN.EQ.l) GO TO 35 DO 5 M=1 , 20 DO 5 J=1 , 23

5 ALBDO ( J , M , MEG ) = 0

.

3 READ (IUT. 7010) ISORS

DO 30 J=ISORS . 23

30 READ (IUT, 7020) (ALBDO (J.M. MEG) , M=l, 20)

35 SWTS ( ISORS, 1 ) =1 . *SGRUPS ( ISORS)

WTSD ( ISORS . 1

) =1 . *DGRUPS ( ISORS

)

DO 100 N=2 , 21 DO 100 M=ISORS , 23

SWTS (M . N) =SWTS (M . N) + ALBDO (M.N-l, MEG) *SGRUPS ( ISORS

)

ENGAM ( N ) = ENGAM ( N ) + ALBDO (M.N-l. MEG ) *FG* SGRUPS ( ISORS

)

100 WTSD ( ISORS. N)=WTSD( ISORS. N) + ALBDO (M.N-l .MEG) *DGRUPS (M) 350 CONTINUE

IF(MODN.EQ.O) CLOSE(ll)

MODN=l

RETURN

END SUBROUTINE SGD ( SGRUPS , SWTS . WTSD , DGRUPS , VEL1 , MODG

)

$ LARGE ALB DIMENSION SGRUPS ( 23 ) , SWTS (23,21), WTSD (23,21). DGRUPS ( 23 )

,

* ALB (17,9,85) ,VEL1(5)

CHARACTER* 80 TITLE

1 FORMAT (A80)

2 FORMAT (215)

3 FORMAT ( IX. 7 ( 1PE9 . 3 ) , 2 ( 1PE8 . 2 )) IF (MODG. EQ . 1 ) GO TO 12

IUT=10

REWIND IUT READ ( IUT. 1) TITLE

IF (IPRINT.LT. 4) GO TO 12 WRITE (6,1) TITLE

12 CONTINUE

DO 25 M=1 , 23 DO 25 N=1 , 21

SWTS(M.N) =0. 25 WTSD (M , N) =0

.

DO 350 K=1 , 85 IF (MODG . EQ . 1 ) GO TO 35

KE=(K-l)/5

LOB=K-5*KE

DO 5 M=1 .

9

DO 5 J=1 . 17

5 ALB ( J , M , K) =0

.

READ ( IUT , 2 ) IK. IE DO 30 J=IE. 17

30 READ (IUT, 3) ( ALB ( J. M , K) , M=1 . 9

)

35 SWTS (IE. 1) =SWTS ( IE , 1 ) +VEL1 ( LOB ) *SGRUPS ( IE

)

WTSD ( IE . 1 ) = WTSD ( IE . 1 ) +VEL1 ( LOB ) *DGRUPS ( IE

)

DO 100 N=2 , 10 DO 100 M=IE , 17

IF ( N . LE . 2 ) GO TO 40 IF ( M . GT . 16 ) GO TO 40

IF(N.GT.M-IO) GO TO 100

40 CONTINUE

SWTS(M.N) =SWTS (M , N) + ALB (M . N-l , K) * SGRUPS ( IE) *VEL1 (LOB

)

100 WTSD ( IE . N ) =WTSD ( IE , N) +ALB (M . N-l , K

) *DGRUPS (M ) *VEL1 ( LOB

)

350 CONTINUE

IF (MODG . EQ . 0 ) CLOSE (10)

MODG=l

RETURN

END SUBROUTINE ORDST ( RHO , A . AL , ORLIST

)

DIMENSION TEMP (21,21). BLST ( 21 ) , ORLIST ( 21

)

AB=A*RHO

B= ( A-AB ) / ( 1 . -AB* (1. 36788-. 36788*A)

)

ALB= ( 1 . -B ) * ( 1 . -AB ) *AL/RHO/ ( 1 . -A)

BLST ( 1

) =1

.

TEMP (1.1) =1. DO 10 1=2.21

BLST ( I

) =B*BLST ( 1-1

)

10 TEMP (1,1) =TEMP ( 1 , 1-1 ) *AB TEMP ( 1 , 1

) =0

.

68
492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

FILE 1 Continued

DO 20 J=2 , 21 DO 20 I=J, 21 20 TEMP ( J, I ) =TEMP ( J-l , 1-1

)

DO 30 1=1.21

ORLIST ( I ) =0

.

DO 30 J-l.I

30

ORLIST ( I ) =ORLIST ( I ) + TEMP ( J, I) *BLST ( J

)

DO 40 1=1.21

40 ORLIST ( I ) =ORLIST ( I

) *ALB

RETURN

END SUBROUTINE DCTISO ( DLNGTH , IC

)

COMMON DIAMS ( 18 ) , DUCTO ( 3 , 17 . 6 ) . DUCT1 ( 3 . 17 , 6 ) , DUCT2 (3.17,6) .IEL.LP, X BEND (10,6), SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO (9,10), RCRRNT (9,10)

COMMON RVCT ( 21 ) , RMTX ( 21 , 21 ) . SO ( 9 ) , SP ( 9 ) , DZRO ( 9 ) , NEXT ( 9 ) , COF ( 9 )

,

X ECC ( 9 ) , DELR ( 9 ) ,SNYNEU(23) .SINVAL(IO) , COl , C02 , VEL (9,21) , IPRINT COMMON RCOSO , RCOSP , RCOS , DGRUPS ( 23 ) , VELO (9,21), VELP (9.21) X , VEL1 ( 5 ) . SCALE , SWTS (23,21). WTSD (23,21)

COMMON I DOSE , HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) , SINR ( 23 ) , THMLN ( 23

)

X , ENGAM ( 21 ) , DIMEN , SCOBLT ( 23 ) , AWTS (6,5), SGAM ( 23 ) . DGAM ( 23 ) . AMOD (5,5)

DIMENSION XI ( 5 ) , ORLIST ( 21 ) , ESC ( 21

)

DO 4 1=1.5

4 XI (I) =0.

1 FORMAT (1H ) 2 FORMAT ( 13 , 11 ( 1PE10 . 2 )

)

3 FORMAT ( 13 , 7E13 . 3

)

IF ( DELR ( IC ) . LT . 0 . ) GO TO 131 DO 130 J=1 , 21

130 ESC( J) =DELR(IC)

GO TO 133

131 ESC(l) =1

.

DO 132 J=2 , 21 132 ESC ( J) =ESC ( J-l ) *

( -DELR ( IC )

)

133 CONTINUE

DLNGTH =DLNGTH /SCALE

EC=ECC(IC)

ECCO= - . 06 * LOG ( EC ) * * 2 * DLNGTH / SQRT (3.182 + DLNGTH * * 2

)

ECCP=- . 109*L0G(EC) **2

ECC0=EXP ( ECCO

)

ECCP=EXP ( ECCP

)

XI ( 1 ) =2

.

XI ( 4 ) =2

.

IF ( DLNGTH . GE . 8 . ) GO TO 10 LX= ( 2

. *DLNGTH+1 .

)

PWR= (DLNGTH-DIAMS(LX) )/( DIAMS ( LX+1 ) -DIAMS ( LX) ) DO 15 KX=1 .

5

Y2= DUCTO (LP , LX+1 , KX)

Yl= DUCTO ( LP , LX, KX)

RCRRNT (IC.KX) =0. 15 CRRNTO ( IC. KX) =Y1* (Y2/Y1) **PWR

GO TO 13

10 DO 12 KX=1 ,

5

AO=DUCTO (LP, 17 , KX)

A2= DUCTO ( LP , 13 , KX)

CALL XTRAP ( AO . A2 , XI ( KX ) , DLNGTH , AXTR

)

12 CRRNTO (IC,KX)=AXTR

13 CONTINUE

AL= ( - . 35+CRRNTO ( IC, 2 ) )/( . 5+DLNGTH**3

)

A=- . 30+CRRNTO ( IC , 3

)

RHO= . 95 RMTX (1,1) = CRRNTO ( IC , 1

)

CALL 0RDST(RH0. A, AL, ORLIST)

DO 19 1=2,21 RMTX (1,1) =ORLIST ( I

)

DO 19 J=2 ,

I

19 RMTX ( J , I

) =RMTX ( J-l , 1-1

)

20 SO ( IC ) = 1 . -CRRNTO ( IC , 4

)

SP (IC) =-l . +2 . *CRRNTO ( IC , 5

)

VELO ( IC , 1 ) =RMTX (1,1) *RVCT ( 1

) *ESC ( 1

) *ECCO

DO 31 J=2 , 21 VELO ( IC , J

) =RMTX ( J, J) *RVCT ( J

) *ESC ( 1) *ECCO

VELP ( IC , J) =0

.

IP=J-1

DO 30 1=1, IP

30 VELP ( IC , J) =VELP ( IC . J) +RMTX ( I , J) *RVCT ( I

) *ESC ( J-I+l

)

VELP ( IC , J

) =VELP ( IC , J) *ECCP

31 CONTINUE

DO 32 J=1 . 21 VEL ( IC , J ) = VELO ( IC , J ) +VELP ( IC, J)

32 RVCT ( J ) = VEL ( IC . J

)

DLNGTH = DLNGTH * SCALE

IF (IPRINT.LT. 4) GO TO 1000

69
r 1 Lt I Lonti nuea

574

575

576

577

578

579

580

581

582

583

584

40 WRITE (6.1)

WRITE (6,3) IC , DLNGTH , ( CRRNTO (IC.KX) ,KX=1,5)

WRITE (6,3) IC. CRRNTO ( IC , 1 ) , A , AL . SO ( IC ) , SP ( IC

)

WRITE (6, 2) I EL, (VELOUC. J) . J=l.ll)

IND=1

WRITE (6, 2) IND. (VELP(IC. J) , J=l.ll) IND = 2 WRITE (6,2) IND, (RVCT( J) , J=l,ll)

1000 RETURN

END

70
FILE 2

1 :* 2:

3:4:

5:6:

7:8:

9: 10: 11 : 12: 13:

14: 15:

16:

17:

18:

19:

20:

21: 22 : 23:

24:

25:

26: 27:

28:

29:

30: 31:

32: 33:

34:

35: 36:

37:

38:

39: 40:

41:

42:

43:

44:

45:

46:

47:

48:

49:

50: 51:

52:

53:

54:

55:

56:

57:

58: 59: 60:

61:

62:

63:

64: 65: 66 : 67: 68 : 69:

70: 71:

72:

73: 74:

75:

76:

77:

78:

79:

80:

81: 82:

SUBROUTINE DCTSLT ( DLNGTH . IC

)

COMMON DIAMS ( 18 ) . DUCTO ( 3 . 17 , 6 ) , DUCT1 ( 3 . 17 . 6 ) , DUCT2 (3,17,6) .IEL.LP. X BEND (10.6). SGRUPS ( 23 ) . DWTS ( 21 ) , CRRNTO (9.10). RCRRNT (9.10)

COMMON RVCT ( 21 ) . RMTX ( 21 . 21 ) . SO ( 9 ) . SP ( 9 ) . DZRO ( 9 ) . NEXT ( 9 ) . COF ( 9 )

,

X ECC ( 9 ) , DELR ( 9 ) . SNYNEU ( 23 ) , SINVAL ( 10 ) , COl . C02 . VEL (9.21). 1PRINT

COMMON RCOSO . RCOSP . RCOS . DGRUPS ( 23 ) , VELO (9.21). VELP (9,21) X , VEL1 ( 5 ) . SCALE , SWTS (23.21). WTSD (23.21)

COMMON IDOSE , HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) , THMLN ( 23

)

X , ENGAM ( 21 ) , DIMEN , SCOBLT ( 23 ) , AWTS(6.5) ,SGAM(23) ,DGAM(23) ,AMOD(5,5)

DIMENSION RMTX1 (21,21) , RMTX2 ( 21 . 21 ) ,V01(21) ,V02(21) ,VP1(21)

.

X VP2 ( 21 ) . XI ( 5 ) . ORLIST ( 21 ) , ESC ( 21

)

DO 4 1=1,5

4 XI (I) =0.

1 FORMAT (1H ) 2 FORMAT ( 13 , 11 ( 1PE10 . 2 )) 3 FORMAT ( 13 . 9E12 . 4

)

IF ( DELR ( IC ) . LT . 0 . ) GO TO 131 DO 130 J=1 , 21 130 ESC( J) =DELR(IC)

GO TO 133

131 ESC ( 1

) =1

.

DO 132 J=2 , 21 132 ESC ( J) =ESC ( J-l ) *

( -DELR ( IC)

)

133 CONTINUE

DLNGTH = DLNGTH /SCALE

C01= . 75/ (1 . +SO(IC-l) ) -0 .

5

C02=l . -COl

U= . 5* ( 1 . +SO ( IC-1 )

)

ZM=U/SQRT ( 1 . -U**2

)

VALO=l . - DLNGTH /ZM

IF ( VALO . LT . 0 . ) VALO=0

.

7 EC=ECC ( IC)

ECCO = -.06* LOG (EC) **2* DLNGTH / SQRT (3.182 + DLNGTH* * 2

)

ECCP=- . 109*LOG(EC) **2

ECCO=EXP(ECCO)

ECCP=EXP ( ECCP) IF ( DLNGTH. GE. 8. ) GO TO 10 LX= ( 2

. *DLNGTH+ 1 .

)

PWR= (DLNGTH-DIAMS (LX) )/( DIAMS ( LX+1 ) -DIAMS ( LX)

)

DO 15 KX=1 ,

5

Y2=DUCT1 (LP.LX+l.KX)

Y1=DUCT1 ( LP . LX , KX)

RCRRNT ( KX . 1 ) =Y1 *( Y2 /Y1 )** PWR Y2=DUCT2( LP.LX+l.KX)

Y1=DUCT2 ( LP , LX , KX) 15 RCRRNT ( KX, 2 ) = Y1* ( Y2/Y1 ) “PWR

GO TO 13

10 DO 12 KX=1 ,

5

A0=DUCT1 ( LP , 17 , KX

)

A2=DUCT1 ( LP , 13 , KX)

CALL XTRAP (AO.A2,XI(KX), DLNGTH . AXTR

)

RCRRNT ( KX. 1

) =AXTR AO=DUCT2(LP. 17, KX)

A2=DUCT2 ( LP . 13 , KX)

CALL XTRAP ( AO , A2 . XI ( KX ) . DLNGTH , AXTR

)

12 RCRRNT (KX. 2) =AXTR

13 CONTINUE

AL1= ( - . 35+RCRRNT ( 2 , 1 ) ) / ( . 5+DLNGTH**3

)

AL2= ( - . 35+RCRRNT (2,2) ) /( . 5+DLNGTH**3

)

Al=- . 3+RCRRNT (3.1)

A2=- . 3+ RCRRNT (3,2)

RHO= . 85

RMTX1 (1.1) =VALO

RMTX2 (1,1) =VALO

CALL 0RDST(RH0.A1.AL1. ORLIST)

DO 17 1=2.21

17 RMTX1 (1,1) =ORLIST ( I

)

CALL ORDST(RHO,A2.AL2. ORLIST)

DO 19 1=2.21

RMTX2 (1,1) =ORLIST ( I

)

DO 19 J=2 ,

I

RMTX1( J.I) =RMTX1( J-l.I-l)

19 RMTX2 (J.I) =RMTX2 ( J-l , 1-1

)

20 RC01= 2. /RCRRNT (4.1)

RC02=2 . / ( 1 . +RCRRNT (4,2)

)

RCP1=1. /RCRRNT (5,1)

RCP2=1 . /RCRRNT (5.2)

V01T=1 . 010 IE9

V02T=2 . 002E-9

VP1T=1 . 0101E-9

VP2T=2 . 0202E-9

DO 31 J=1 , 21

71
FILE 2 Continued

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:

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:

V01 ( J ) =RMTX1 ( J . J) *VELO ( IC-1 . J

) *ESC ( 1

)

V02 ( J) =RMTX2 ( J, J) *VELO ( IC-1 , J) *ESC ( 1) VP1 ( J) =0. VP2 ( J) =0

.

JP=J-1

IF(J-l) 28,28,29

29 DO 30 1=1, JP VP1(J)=VP1( J)+RMTX1(I, J)*VELO(IC-l.I)*ESC( J-Id)

30 VP2 ( J) =VP2 ( J) +RMTX2 ( I , J) *VELO ( IC-1 , I

) *ESC ( J-I+l

)

28 V01T=V01T+V01 ( J)

V02T=V02T+V02 ( J)

VP1T=VP1T+VP1 ( J)

VP2T=VP2T+VP2 ( J)

31 CONTINUE

DO 32 J=1 , 21 VELO ( IC, J) =C01*V01 ( J) +C02*V02 ( J)

32 VELP ( IC , J) =C01*VP1 ( J) +C02*VP2 ( J) DO 33 J=1 , 21 VELO(IC, J) =VELO(IC, J) *ECC0

VELP ( IC . J ) = VELP ( IC , J

) *ECCP

VEL ( IC , J) =VELO ( IC , J) +VELP ( IC , J

)

33 RVCT ( J ) = VEL ( IC , J

)

SO(IC) =-l. +2.* (C01*V01T+C02*V02T)/(C01*V01T*RC01+C02*V02T*RC02)

SP ( IC) =-l . +2 . * (C01*VP1T+C02*VP2T) / (C01*VP1T*RCP1+C02*VP2T*RCP2

)

DLNGTH = DLNGTH * SCALE

IF ( IPRINT . NE . 4 ) GO TO 1000

40 WRITE (6,1) IC .DLNGTH, (RCRRNT(KX.l) ,KX=1,5)

IC .DLNGTH. (RCRRNT(KX, 2) ,KX=1,5)

IC , AL1 , A1 , AL2 , A2 , COl , SO ( IC ) .SP(IC) ,C02 IC . V01T . V02T , VP1T , VP2T . RCOl . RC02 . RMTX1 (1.1), RMTX2 (1.1)

IEL, (VELO(IC. J) , J=l,ll)

IEL, (VELP ( IC. J) . J=l.ll)

IND. (VOl(J) , J=l,21.2)

IND. ( V02 ( J) , J=l,21.2)

WRITE (6.3)

WRITE (6. 3) WRITE (6. 3)

WRITE (6, 3) WRITE (6. 2) WRITE (6, 2)

IND=13

WRITE (6. 2)

WRITE (6, 2)

1000 RETURN

END SUBROUTINE XTRAP ( AO . A2 . XI , FA . AINT

)

SAVE=XI

FO = 8

.

F2 = 6 XI=L0G(A2/A0) /LOG ( FO/F2 ) - . 05

DL=0 .

1

FOL=FO**DL

F2L=F2**DL

FAL=FA**DL

DO=AO*FO**XI*FOL

D2=A2*F2**XI*F2L

AL1= ( DO-D2 ) / ( FOL-F2L

)

IF(ALl) 10.10.5

5 AL2= ( F0L*D2-F2L*D0 ) / (FOL-F2L) AINT= ( AL1+AL2/FAL) /FA**XI

GO TO 1000

10 XI=LOG ( A2/AO) /LOG ( FO/F2

)

ALX= ( AO*FO** (XI+1 . ) -A2*F2** (XI+1 ) )/(FO-F2)

AINT=AL1/FA**XI

XI = SAVE

1000 RETURN

END SUBROUTINE TURN ( SEN. IC)

COMMON DIAMS (18 ) . DUCTO (3.17.6), DUCT1 (3.17,6), DUCT2 (3.17,6), IEL . LP

.

X BEND (10.6), SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO (9.10). RCRRNT (9.10)

COMMON RVCT ( 21 ) , RMTX ( 21 . 21 ) , SO ( 9 ) , SP ( 9 ) . DZRO ( 9 ) . NEXT ( 9 ) , COF ( 9 )

.

X ECC ( 9 ) . DELR ( 9 ) . SNYNEU ( 23 ) . SINVAL ( 10 ) . COl , C02 . VEL (9.21), IPRINT COMMON RCOSO , RCOSP . RCOS , DGRUPS ( 23 ) . VELO (9,21), VELP (9,21)

X , VEL1 ( 5 ) , SCALE . SWTS (23.21), WTSD (23.21)

COMMON I DOSE , HNDRSN ( 23 ) . SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) , THMLN ( 23

)

X , ENGAM ( 21 ) . DIMEN , SCOBLT ( 23 ) , AWTS (6,5), SGAM ( 23 ) . DGAM ( 23 ) , AMOD (5.5)

DIMENSION RX1 ( 13 ) , RINF (13 ) , ORLIST ( 21 ) , ESC ( 21

)

DATA RX1/.88, .86, .84, .82, 79, . 74 , . 69 , . 62 , . 54 . .44, .26, .01, .001/ DATA RINF/ . 86 . . 83 . . 80 , . 76 . . 71 , . 66 , . 58 , . 48 , . 37 , . 24 , . 01

,

* .001.. 0001/

IF ( DELR (IC ) . LT . 0 . ) GO TO 131 DO 130 J=1 , 21 130 ESC ( J) =DELR (IC

)

GO TO 133

131 ESC (1 ) =1

.

DO 132 J=2 , 21 132 ESC ( J) =ESC ( J— 1

) *

( -DELR (IC)

)

133 CONTINUE

C0=0.

CP=0.

72
FILE 2 Continued

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

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

AZFAC=1 . -DZRO ( IC

)

CP1=0.

CP2=0

RC0S0=2 . / ( 1 . +S0 ( IC-1 )

)

RC0SP=2 . / ( 1 . +SP ( IC-1)

)

DO 305 1=1.21

CO=CO+VELO ( IC-1 , I

)

CP=CP+VELP ( IC-1 . I

)

VELO ( IC , I ) =0

VELP ( IC , I ) =0

305 RVCT ( I

) = VELO ( IC-1 . I

)

RCOSO= ( CO*RCOSO+CP*RCOSP ) / ( CO+CP

)

RCOSP=RCOSO

RCOS=RCOSO

IX=1

134 SEN=2 . / (RCOS+1 . E-6 ) -1

C RCSS=COF (IC)

KS= ( SEN* 10 . +1 .

)

IF(KS.LT.l) KS=1

RKS=KS

PWR= ( SEN*10 . +1 . 00001-RKS

)

C1=RX1 ( KS ) * (RX1 ( KS+1 ) /RX1 ( KS ) )**PWR

CINF=RINF ( KS ) * ( RINF (KS + 1) /RINF ( KS )

) * *PWR

BD= ( 1 . -Cl ) / (Cl-CINF

)

BN=CINF*BD

RX=RCSS**1 .

RDUX= ( 1 . +BN*RX) / ( 1 . +BD*RX)

C

1 FORMAT (1H ) 2 FORMAT ( 13 , 11 ( 1PE10 . 2 )

)

3 FORMAT ( 13 , 7E13 . 3

)

LS= ( 5

. *SEN+3 .

)

DO 15 KX=1 ,

Y2=BEND ( LS+1 , KX

)

Y1=BEND ( LS . KX) PWR= ( SEN-SINVAL ( LS ))/( SINVAL ( LS +1 ) -SINVAL ( LS )+. 00001

)

15 CRRNTO ( IC, KX) =Y1* ( Y2/Y1) **PWR

AL=CRRNTO ( IC , 2

)

A=CRRNTO ( IC , 3

)

SO ( IC) =CRRNTO ( IC , 4

) *SQRT ( 2 . -SEN) -1

SP ( IC) =2 . * CRRNTO ( IC . 5

) -1

IF(IX.GT.O) S01=S0(IC)

IF(IX.GT.O) SP1=SP(IC)

IF(IX.LT.O) S02=S0 ( IC)

IF(IX.LT.O) SP2=SP(IC)

RHO= . 75 RMTX (1.1) = CRRNTO ( IC . 1

)

CALL ORDST ( RHO , A , AL , ORLIST

)

DO 20 1=2.21 RMTX ( 1 . I ) =ORLIST ( I

)

DO 19 J=2 .

19 RMTX ( J, I ) =RMTX ( J-l . 1-1

)

20 CONTINUE

VELO ( IC . 1 ) =VELO ( IC , 1 ) +RMTX (1.1) *RVCT ( 1

) *ESC ( 1

)

IF(IX.GT.O) C01=VEL0 ( IC , 1

)

IF(IX.LT.O) C02=VEL0 ( IC, 1

)

DO 31 J=2 , 21 VELO(IC. J) =VELO(IC. J) +RMTX( J, J) *RVCT ( J

) *ESC ( 1

)

IP=J-1

DO 30 1=1, IP DVP=RMTX ( I . J) *RVCT ( I

) *ESC ( J-I+l ) *RDUX

30 VELP ( IC , J) =VELP ( IC , J) +DVP IF(IX.GT.O) C01=C01+VEL0 ( IC , J) IF(IX.GT.O) CPI =CP1+VELP ( IC , J) IF(IX.LT.O) C02=C02+VEL0 ( IC , J) IF(IX.LT.O) CP2=CP2+VELP(IC, J) 31 CONTINUE

C IF(IX.LT.O) GO TO 325

REN0RM=1 . + ( 1 . -AZFAC) *C01/ (C0+CP*1 . E-5

)

DO 315 J=1 , 21 VELO(IC, J) =VELO ( IC , J) * AZFAC

VELP ( IC . J ) = VELP ( IC , J

) *RENORM

315 RVCT ( J) =VELP (IC-1. J) C01=C01*AZFAC

CP1=CP1*REN0RM

RCOS=RCOSP

IX=IX-2

GO TO 134

325 CONTINUE

73
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

288

289

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

319

320

321

322

323

324

325

326

327

328

329

330

32

40

FILE 2

SO ( IC) =-l . + ( 1 . +S01 ) * ( 1 . +S02) *C02/ (C01* ( S02-S01 ) +C02* ( 1 . +S01 )

)

SP ( IC) =-l . + ( 1 . +SP1 ) * ( 1 . +SP2 ) *CP2/ (CPI* (SP2-SP1 ) +CP2* ( 1 . + SP1 )

)

DO 32 J=1 , 21 VEL(IC, J) =VELO ( IC , J) +VELP ( IC , J) RVCT ( J) =VELP ( IC , J) IF ( IPRINT . NE . 4 ) GO TO 1000

WRITE (6.1)

IC, SEN, ( CRRNTO ( IC , KX ) ,KX=1.5)

IC , RMTX (1,1) , AL , A , SO ( IC ) , SP ( IC

)

IEL, (VELO(IC. J) , J=l,ll)

WRITE (6, 3) WRITE (6, 3) WRITE ( 6 , 2) IND=1

WRITE (6,2) IND = 2 WRITE (6,2)

1000 RETURN

END SUBROUTINE BTODR ( DLNGTH

IND. (VELPdC. J) , J=l,ll)

IND, (RVCT ( J) . J=l.ll)

IC)

COMMON DIAMS ( 18 ) , DUCTO ( 3 , 17 , 6 ) , DUCT1 ( 3 , 17 . 6 ) , DUCT2 (3,17.6) .IEL.LP, X BEND ( 10 , 6 ) , SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO ( 9 , 10 ) .RCRRNT ( 9 . 10

)

COMMON RVCT (21) , RMTX (21. 21) .SO (9) ,SP(9) ,DZR0(9) ,NEXT(9) ,COF(9)

,

X ECC ( 9 ) , DELR ( 9 ) ,SNYNEU(23) ,SINVAL(10) ,C01.C02,VEL(9, 21) .IPRINT COMMON RCOSO , RCOSP . RCOS . DGRUPS ( 23 ) . VELO (9.21), VELP (9,21) X , VEL1 ( 5 ) , SCALE . SWTS (23.21). WTSD (23.21)

COMMON IDOSE , HNDRSN ( 23 ) . SEPI ( 23 ) . STH ( 23 ) , SINR ( 23 ) , THMLN ( 23

)

X. ENGAM ( 21 ) . DIMEN, SCOBLT( 23) , AWTS(6,5) ,SGAM(23) , DGAM ( 23 ) ,AMOD(5. 5) DIMENSION STRVEL(22) .TVLO (22) .TVLP (22) .SIDE (22, 22) .BACK (22. 22)

.

X VPB ( 22 ) . VPS ( 22 ) , VOS ( 22 ) , VOB ( 22

)

1 FORMAT (1H ) 2 FORMAT ( 13 , 11 ( 1PE10 . 2 )) 3 FORMAT ( 13 , 9E12 . 4

)

STRSO=SO ( IC-1

)

DLBND=1

.

DISTNS=DLNGTH+DLBND

DO 10 1=1.21

10 STRVEL ( I

) =VELO ( IC-1 , I

)

TVO=0

.

TVP=0

C

C DIRECT COMPONENT

CALL DCTSLT ( DLNGTH. IC) DO 20 1=1,21

TVLO (I ) =VELO ( IC , I

)

TVLP ( I ) =VELP ( IC , I

)

TVO=TVO+VELO (IC , I 20 TVP=TVP+VELP ( IC , I

)

RCO=TVO*2 . / ( 1 . 001+S0 ( IC)

)

RCP=TVP*2 ./ (1.001+SP(IC)

C

C SIDE COMPONENT

DO 200 1=2,21

200 VELO (IC-1, I )=0.

VELO ( IC-1, 1)=1.

VELP ( IC, 1

) =0

.

SO(IC-l ) =0 .

5

CALL DCTSLT ( DISTNS, IC) DO 205 1=1,20

205 SIDE(1,I)=(VELP(IC,I+1)-.54*VELP(IC.I) )/.1162

DO 210 1=2,21

IP=I-1

DO 207 J=1.IP

207 SIDE ( I , J) =0

.

DO 210 J=I , 21 210 SIDE ( I , J) =SIDE (1-1 , J-l)

VPS(1)=0.

VOS(1)=0.

DO 220 J=2 , 20

VOS ( J) =SIDE ( J . J) * ( VELP ( IC-1 , J) - . 2055*STRVEL ( J-l )

)

VPS(J)=0.

JP=J-1

DO 220 1=1, JP VPS ( J) =VPS ( J) +SIDE (I , J) * VELP (IC-1. I) IF ( I . GT . 1 ) VPS( J)=VPS( J) -SIDE (I . J) * . 2055*STRVEL (1-1)

220 CONTINUE

SVO=0.

SVP=0.

DO 230 1=2.20

TVLO ( I

) =TVLO (I ) + . 829*VOS ( I

)

TVLP ( I

) =TVLP ( I ) + . 829*VPS ( I

)

SVO=SVO+ . 829*V0S (I

)

230 SVP=SVP+ . 829*VPS ( I

)

TVO=TVO+SVO

TVP=TVP+SVP

Continued

74
331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

FILE 2

RC0=RC0+SV0*2 . / ( 1 . 001+SP ( IC )

)

RCP=RCP+SVP*2 . / ( 1 . 001+SP ( IC)

)

C

C BACK COMPONENT

DO 300 1=2,21

VELO ( IC-1 , I ) =0

300 RVCT ( I ) =0

RVCT ( 1 ) = 1

VELO ( IC-1 , 1 ) =1

CALL DCTISO ( DISTNS , IC

)

BVP=0

BV0=0

DO 305 1=1,21

BVO=BVO+VELO ( IC . I

)

BVP=BVP+VELP ( IC , I

)

305 BACK (1,1) =VEL ( IC , I

)

RC0B=2 . / ( 1 . 001+SO( IC)

)

RCPB=2 ./( 1 . 001+SP ( IC)

)

C SO(IC-l) =-0.

CALL DCTSLT ( DISTNS, IC)

BVPM=0.

DO 310 1=2,21

BVPM=BVPM+VELP ( IC , I

)

310 BACK ( 1 , 1 ) = ( BACK ( 1 „ I

) - 1 . 3

7 *VELP ( IC , I ) )/.205

BACK (1,1)=BACK( 1,1)/. 205

RCSPB=BVP*RCPB-1 . 4*BVPM*2 . / ( 1 . 001+SP ( IC)

)

BVP=BVP-1.4*BVPM

RCSPB=RCSPB/BVP

C DO 315 1=2,21

IP=I-1

DO 312 J=1,IP

312 BACK ( I , J) =0

DO 315 J=I , 21 315 BACK(I, J) =BACK(I-1, J-l)

VOB ( 1

) =0

VPB(l) =0. DO 320 J=2 , 21 VOB( J)=BACK( J, J)* (VELP (IC-1, J) +STRVEL ( J-l )

)

VPB ( J ) =0. JP=J-1

DO 320 1=1, JP IF(I.GT.l) VPB ( J) =VPB ( J) +BACK ( I , J) *STRVEL ( 1-1

)

320 VPB ( J) =VPB ( J) +BACK ( I , J) *VELP ( IC-1 , I

)

SV0=0.

SVP=0

DO 330 1=1,21

SVO=SVO+ .171*VOB(I)

SVP=SVP+ . 171*VPB ( I

)

VELO ( IC , I

) =TVLO ( I ) + . 171*VOB ( I

)

VELP ( IC , I

) =TVLP ( I ) + , 171*VPB ( I

)

330 VEL ( IC , I

) =VELO ( IC , I ) +VELP ( IC , I

)

TVO=TVO+SVO

RCO=RCO+SVO*RCOB

TVP= TVP + SVP

RCP = RCP + SVP *RCSPB

SO(IC) =-l. +2, *TVO/RCO

SP ( IC) =-l . +2 . *TVP/RCP

VELP ( IC , 1

) =0

VEL ( IC , 1 ) =VELO ( IC , 1

)

DO 50 1=1,21

VELO ( IC-1 , I

) =STRVEL ( I

)

50 RVCT ( I

) =VEL ( IC , I

)

SO ( IC-1 ) =STRSO

IF ( IPRINT . NE . 4 ) GO TO 1000

41 WRITE (6,1)

IC,S0(IC) ,SP(IC) , TVO , RCO , TVP , RCP , VELO ( IC , 2 ) ,VELP(IC,2

I EL, (VELO(IC. J) , J=l,21,2)

IEL. (VELP ( IC , J) , J=l,21,2)

WRITE (6,3)

WRITE (6, 2) WRITE (6, 2) IND=17

WRITE (6, 2) WRITE (6, 2) 1000 RETURN

END

IND.

IND,

(TVLO(J) , J=1.21,2)

(TVLP(J) , J=1 , 21,2)

Continued

)

75
FILE 3

1 :*

2 : 3:4:

5:6:

7:

8 : 9: 10 : 11 : 12 : 13: 14:

15:

16:

17:

18: 19:

20: 21 : 22 : 23:

24:

25:

26:

27: 28:

29:

30:

31:

32:

33:

34:

35: C

36: 37:

38: C

39: 40:

41:

42:

43:

44:

45: C 46:

47:

48: 49: C

50:

51: 52:

53: C

54:

55:

56:

57:

58:

59: 60:

61:

62: 63:

64: 65: 66 : 67: 68 : 69:

70:

71 : 72:

73:

74: 75:

76: C

77:

78: C

79:

80:

81:

82:

SUBROUTINE ROOMX ( IC , SVELO , SVELP

)

COMMON DIAMS ( 18 ) . DUCTO ( 3 , 17 , 6 ) , DUCT1 ( 3 , 17 . 6 ) , DUCT2 (3,17,6) . IEL.LP, X BEND (10.6), SGRUPS ( 23 ) . DWTS ( 21 ) , CRRNTO (9.10). RCRRNT (9,10)

COMMON RVCT ( 21 ) . RMTX ( 21 , 21 ) , SO ( 9 ) , SP ( 9 ) , DZRO ( 9 ) , NEXT ( 9 ) . COF ( 9 )

,

X ECC ( 9 ) , DELR ( 9 ) , SNYNEU ( 23 ) , SINVAL ( 10 ) , C01 . C02 . VEL (9.21), IPRINT COMMON RCOSO . RCOSP . RCOS . DGRUPS ( 23 ) . VELO (9.21), VELP (9,21) X , VEL1 ( 5 ) . SCALE , SWTS (23.21), WTSD (23,21)

COMMON XDOSE . HNDRSN ( 23 ) . SEPI ( 23 ) . STH ( 23 ) , SINR ( 23 ) , THMLN ( 23

)

X, ENGAM ( 21 ) .DIMEN. SCOBLT( 23) . AWTS(6,5) ,SGAM(23) ,DGAM(23) ,AMOD(5,5)

DIMENSION D ( 20 ) . DL ( 20 ) , W ( 20 ) , OVEL (22.20), PVEL (22,20), SVELO (22,20).

* SVELP (22, 20)

1 FORMAT (1H ) 2 FORMAT ( 13 , 11 ( 1PE10 . 2 )

)

3 FORMAT ( 11 ( 1PE10 , 2 )) WTH=C0F ( IC)

HT=ECC(IC)

XL = DELR ( IC) AEX= SCALE

AEXR=DZRO ( IC+1

)

IF (AEXR . GT . 1 . ) AEXR=1

.

WD=DZRO ( IC)

RMCASE=DZRO ( 1

)

WDP=COF ( IC+1

)

NW= (ECC ( IC+1 ) +1 . 01)

NL=( DELR (IC+1) +1.01)

NT=NW*NL

DELW= WDP /ECC ( IC + 1

)

DELL=XL/DELR ( IC+1

)

TN=SQRT ( 1 . -SO ( IC-1 ) **2) /SO(IC-l)

A=SP ( IC-1 ) * (1 . +SO ( IC-1 ) ) / ( 1 . +SP ( IC-1 ) ) /SO ( IC-1

)

RR=AEX/(WTH*HT+HT*XL+XL*WTH)

AL= . 75 ALR=AL* ( 1 . -RR/2 .

)

TNR=SQRT ( HT*WTH/3 . 14159 ) /XL

POSITION 1,N=0 (UNREFLECTED, IN FRONT OF DUCT)

OVEL (1,1) =2 . / ( 1 . +SO ( IC-1 ) ) *AEXR

PVEL ( 1 , 1

) =2 . / ( 1 . +SP ( IC-1 ) ) *AEXR

POSITIONS 2 THRU NT. N IS 2 OR GREATER

DO 11 1=2, NT OVEL ( 3 . I ) = ( 1 . +SO ( IC-1 )

) *RR*ALR* *2 PVEL ( 3 , I

) = ( 1 . +SP ( IC-1 )

) *RR*ALR**2

DO 11 J=4 . 21 OVEL ( J. I

) =ALR*OVEL ( J-l , I

)

11 PVEL ( J , I

) =ALR*PVEL ( J-l , I

)

POSITION 1. N IS 2 OR GREATER

DO 12 J=3 , 21 OVEL ( J , 1 ) = . 5 *OVEL ( J , 2

)

12 PVEL ( J , 1 ) = . 5*PVEL ( J , 2

)

POSITIONS 2 THRU NW, N=0 (UNREFLECTED)

DO 8 1=2, NW

OVEL ( 1 , I

) =0

.

8 PVEL ( 1 , I ) =0

.

POSITIONS ( NW+1 ) THRU NT,N=0 (UNREFLECTED)

DO 13 1=1, NL

FI= ( 1-1

)

DO 13 J=1 , NW

F J= ( J-l

)

K= ( 1-1 ) *NW

DL ( K+ J) =DELL*FI

W ( K+ J) =DELW*FJ

13 D ( K+ J) =SQRT ( DL (K+J)**2+W(K+J)**2 +.0001)

K=NW+1

DO 25 I = K , NT

FR= ( DL ( I

) *TN-W ( I) + . 5*WD ) /WD

IF(FR) 16.16,17

16 FR=0.

17 IF ( FR-1 . ) 19.18,18

18 FR=1

.

19 SAF=SQRT ( 1 , +AEX/3 .14159/D(I)**2)

SAF=DL ( I ) /D ( I

) * ( 1 . -1 . /SAF) IF ( SAF- ( 1 . -SO ( IC-1 ) ) >20.20,21

21 SAF=1 . -SO ( IC-1

)

20 OVEL (1,1) =FR*SAF/ ( 1 . -SO ( IC-1 )

)

PVEL ( 1 , I

) = ( 1 . +SP ( IC-1 )

) * ( A*OVEL ( 1 , I ) / ( 1 . +SO ( IC-1 )

) + ( 1 . -A) *SAF)

25 CONTINUE

MODIFY FOR BEND INTO ROOM CASES

IF (RMCASE) 100.51,150

ENTERING UNREFLECTED MOVING TOWARDS DETECTORS

100 K=NW+1

DO 102 I=K , NT , NW OVEL ( 1 , I

) =0

.

DO 102 J=2 , NW 76
FILE 3 Continued

83

8485

8687

8889

90

9192

93

9495

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

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

102 OVEL ( 1 , 1+ J-l ) =2 . *OVEL ( 1 , 1+ J-l

)

GO TO 51 C ENTERING UNREFLECTED MOVING TOWARDS WALL

150 DO 152 1=1. NT TEMP=DL ( I

)

DL ( I ) =W ( I ) + ( WTH-WDP)

W ( I

) =TEMP

152 D ( I) =SQRT (DL( I )**2+W(I)**2+ .0001)

TEMP=NW

NW=NL

NL=TEMP

K=NW+1

DO 154 I=K , NT

154 OVEL ( 1 . 1

) = 0

GO TO 32

C ONCE-REFLECTED COMPONENT (N=l)

51 SPX=1.

RMX=0

IF(TNR.GT.TN) GO TO 30

SPX=0.

RMX=1

30 NLHF= (NL/2+1) DO 31 1=1 , NLHF

KI= (1-1) *NW KIR= (NL-I ) *NW

DO 31 J=1 , NW

TEMP=DL( J+KI

)

DL ( J+KI ) =DL ( J+KIR)

31 DL( J+KIR) =TEMP

DO 34 1=1. NT

34 D ( 1 ) =SQRT(DL (I)**2+W(I)**2 +.0001)

K=NT-NW+1

DO 33 J=1 , NT IF(J.LT.K) GO TO 35

OVEL (2. J) =OVEL ( 1 , J) *AL*2 . *DL ( J) /D( J) +RMX* ( 1 . +SO ( IC-1 )

) *AL*RR/2

PVEL ( 2, J) =PVEL ( 1 , J) *AL*2 . *DL ( J) /D ( J) +RMX* ( 1 . +SP ( IC-1 ) )*AL*RR/2.

GO TO 33

35 SAFP=SQRT ( 1 . +AEX/3 .1415926/D(J)**2)

3AFP=DL (J)/D(J)*(1.-1. /SAFP

)

OVEL ( 2 , J) = ( 1 . +SO ( IC-1 )

) *AL* ( SPX*SAFP+RMX*RR)

PVEL ( 2 , J) = ( 1 . +SP ( IC-1 )

) *AL* ( SPX* (A* SAFP/ (l.+SO(IC-l) )+() . -A) *RR)

* +RMX*RR)

33 CONTINUE

GO TO 57

C RMCASE=1 . (N=l

)

32 OVEL (2,1) =AL* ( 1 . +S0 ( IC-1 )

)

PVEL (2.1) =AL* ( 1 . +SP ( IC-1 )

)

DO 36 J=2 , NT SAFP=SQRT (1.+AEX/3. 1415926/D (J)**2)

SAFP=DL ( J) /D ( J) * ( 1 . -1 . /SAFP

)

OVEL ( 2 , J) =AL* (l.+SO(IC-l) )*SAFP

36 PVEL ( 2 , J) =AL* ( 1 . +SP ( IC-1 )

) *SAFP

TEMP=NW

NW=NL

NL=TEMP

C RENORMALIZATION

57 DO 55 1=1, NT

DO 55 J=1 . 21 IJ=I+J

IF ( I J . EQ . 2) GO TO 55 OVEL ( J, I

) =OVEL ( J. I

) *OVEL (1,1)

PVEL ( J, I

) =PVEL ( J , I

) *PVEL (1,1)

55 CONTINUE

C MATRIX GENERATION AND MULTIPLICATION ( S

)

DO 70 1=1. NT

DO 58 J=1 , 21

SVELP ( J , I

) =0

58 SVELO ( J, I) =OVEL (1,1) *VELO ( IC-1 . J) +PVEL (1,1) *VELP ( IC-1 . J) DO 65 IP=2 , 21 DO 65 J=2 . IP

65 SVELP ( IP, I

) = SVELP ( IP , I ) +OVEL ( J, I) *VELO ( IC-1 , IP+1- J) +PVEL ( J, I

)

* VELP (IC-l.IP+l-J)

70 CONTINUE

IF ( IPRINT . NE . 4 ) GO TO 1000

40 WRITE (6,1)

WRITE (6,2) IC . AEX , AEXR . WD , RMCASE . DELW , DELL , TN , TNR , A , RR , ALR WRITE (6, 1) WRITE (6,3) ( SVELO (J,l) . J=1 , 21 , 2) WRITE (6. 3) ( SVELO ( J, 2 ) , J=1 ,21,2)

WRITE (6,3) ( SVELO ( J , NT ) , J=1 , 21 , 2

)

WRITE (6,3) (OVEL(J.l) . J=l,10)

WRITE (6, 1) WRITE (6, 3) ( SVELP (J.l) . J=l,21,2)
r ILt i uontinuea

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

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

( SVELP ( J , 2 ) , J=l,21.2)

(SVELP(J.NT) . J=1.21.2)

( PVEL (J,l) , J=1 , 10

)

WRITE (6. 3)

WRITE (6, 3) WRITE (6. 3) WRITE (6.1)

1000 RETURN

END SUBROUTINE DOSX ( IC . SVELO , SVELP

)

COMMON DIAMS ( 18 ) . DUCTO ( 3 , 17 , 6 ) . DUCT1 ( 3 . 17 . 6 ) , DUCT2 (3.17,6) .IEL.LP. X BEND ( 10 , 6 ) , SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO (9,10) , RCRRNT (9,10)

COMMON RVCT ( 21 ) , RMTX ( 21 . 21 ) , SO ( 9 ) , SP ( 9 ) , DZRO ( 9 ) , NEXT ( 9 ) . COF ( 9 )

,

X ECC ( 9 ) , DELR ( 9 ) , SNYNEU ( 23 ) , SINVAL ( 10 ) . COl . C02 . VEL (9,21), IPRINT COMMON RCOSO . RCOSP . RCOS , DGRUPS ( 23 ) , VELO (9.21). VELP (9,21) X , VEL1 ( 5 ) , SCALE , SWTS (23,21) , WTSD (23,21)

COMMON IDOSE . HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) , THMLN ( 23

)

X , ENGAM ( 21 ) , DIMEN. SCOBLT( 23) . AWTS(6.5) ,SGAM(23) ,DGAM(23) ,AMOD(5.5)

DIMENSION SVELO (22. 20) , SVELP ( 22 . 20 ) .PRAT (23. 2) ,ORAT(20,2)

.

* DRAT (20.2). SPCTRM ( 23 ) . DOSE ( 20 ) , DFLXN ( 21

)

1 FORMAT (1H ) 2 FORMAT (5H IEL=.I2.5H IC=.I2)

3 FORMAT (8H DIR/ICD)

4 FORMAT (8H REF/ICD)

5 FORMAT (8H TOT/ICD)

6 FORMAT (8H DIR/IFD)

7 FORMAT (8H REF/IFD)

8 FORMAT (30H DOSE/ENTRANCE DOSE

10 FORMAT ( 9 ( 1PE12 . 4 )

)

11 FORMAT (8H RC0SIN= , 1PE10 . 3 , 7H X7H EXF= , 1PE10 . 3

)

12 FORMAT ( 55H) MEG,

N

FORMAT (5H DOSE)

FORMAT ( 55H VEL(N)/IC=1

FORMAT ( 15 . 7 ( 1PE10 . 3 )

)

FORMAT (5H SO= ,7(1PE10

FORMAT (5H SP = . 7 ( 1PE10

RCOSIN=2 . / ( 1 . + COF ( 1 )

)

CDIN=DWTS ( 1

)

IF (CDIN . EQ . 0 . ) CDIN=1

.

NW= ( ECC (IC) +1.01) NL= ( DELR (10+1.01)

NT=NW*NL

NMAX-- 21 IF ( IDOSE. NE. 5) GO TO 122

NMAX=10

122 DO 100 MEG=1 . 23

SPCTRM (MEG) =0

.

U=l.

DO 100 N=1 . NMAX

SPCTRM ( MEG ) = SPCTRM ( MEG ) + ( SVELO (N, 1) + SVELP (N.1))*SWTS(MEG.N)/U

100 U=U* . 696

W ( FT) = , F7 . 2 , 9H L(FT)=.F7.2)

CDIN= , 1PE10 . 3 , 7H DIN= , 1PE10 . 3

.

13

1415

1617

20

SOURCE

3) ) .3) )

SPCTRM DFLXN VELO

4

VELP )

DO 50 1=1. NT ORAT ( 1 , 1 ) = 0

.

PRAT ( I , 1

) =0

.

U=i.

DO 40 NREF=1 , NMAX

ORAT (1,1) =ORAT (1,1) + SVELO ( NREF . I

) *DWTS ( NREF ) /U PRAT ( 1 , 1

) = PRAT ( 1 , 1 ) + SVELP ( NREF . I

) * DWTS ( NREF ) /U DFLXN ( NREF ) = ( SVELP ( NREF . 1 ) + SVELO ( NREF , 1 )

) *DWTS ( NREF ) /U

40 U=U* . 696 DRAT (1,1) =PRAT (1,1) +ORAT ( I . 1

)

DOSE ( I

) =DRAT (1,1) 50 CONTINUE

IF (NEXT ( 1 ) ) 101.101.102

102 ZD=DELR ( 1

)

EXF= ( . 671 5+ ZD* ( . 924+ . 565*ZD) ) / ( . 0685+ZD* ( . 5445+ . 137*ZD)

)

CDIN=CDIN*EXP ( ZD*EXF

)

101 DIN=CDIN*RCOSIN

DO 120 1=1, NT ORAT (1,2) =0RAT (1,1) /DIN PRAT (1.2) =PRAT (1,1) /DIN DRAT ( I , 2

) = DRAT (1,1) /DIN

ORAT (1,1) =ORAT (1,1) /CDIN PRAT (1,1) =PRAT (1,1) /CDIN

120 DRAT (1,1) =DRAT (1,1) /CDIN

IF ( IPRINT. EQ.O) GO TO 1000

IF ( IPRINT . LE . 2 ) GO TO 128

WRITE (6. 2) I EL , IC

105 WRITE (6,1)

WRITE (6.12)

JC=IC-2

78
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

288

289

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

319

320

321

322

323

324

325

326

327

328

329

330

FILE 3 Continued

DO 109 1=1.21

WRITE (6,15)1, SGRUPS ( I ) .SPCTRM(I) .DFLXN(I) ,VELO( JC.I) ,VELP( JC,I)

109 CONTINUE

MEG=22

WRITE (6.15) MEG . SGRUPS ( MEG ) . SPCTRM ( MEG

)

MEG=23

WRITE (6,15) MEG . SGRUPS ( MEG ) . SPCTRM ( MEG

)

WRITE (6.1)

WRITE (6.1)

WRITE (6. 14) DO 113 N=1 , 21 WRITE (6,15) N. (VEL(ID.N) ,ID=1, JC)

113 CONTINUE

WRITE (6,16) (SO(ID) ,ID=1, JC)

WRITE (6,17) (SP(ID) ,ID=1. JC)

WRITE (6,1) IF (IPRINT.LT. 4) GO TO 128

WRITE (6. 11) RCOSIN, CDIN.DIN.EXF

WRITE (6, 3)

WRITE (6,10) ( ORAT (1,1) ,1=1, NT

)

WRITE (6,4)

WRITE (6. 10) ( PRAT (1,1) .1=1, NT

)

WRITE (6,1)

WRITE (6. 5) DO 125 J= 1 , NW

WRITE(6„ 10) ( DRAT ( 1+ J-l , 1) ,I=1,NT,NW)

125 CONTINUE

WRITE (6. 6)

WRITE (6. 10) ( ORAT (1.2) ,1=1. NT)

WRITE (6, 7) WRITE (6,10) ( PRAT (1,2) ,1=1, NT

)

WRITE (6,1)

128 WRITE (6. 8) COF ( IC ) , DELR ( IC-1

)

WRITE (6,1) DO 130 J=1 , NW WRITE (6. 10) ( DRAT (1+ J-l. 2) ,I=1,NT.NW)

130 CONTINUE

WRITE (6.1)

WRITE (6. 13) DO 135 J=1 , NW WRITE (6. 10) ( DOSE (1+ J-l) ,I=1.NT,NW)

135 CONTINUE

1000 RETURN

END

L'O

SUBROUTINE DOSE(IC)

COMMON DIAMS ( 18 ) . DUCTO ( 3 , 17 , 6 ) , DUCT1 ( 3 , 17 , 6 ) , DUCT2 (3,17,6) ,IEL,LP, X BEND (10,6), SGRUPS ( 23 ) , DWTS ( 21 ) , CRRNTO (9,10). RCRRNT (9,10)

COMMON RVCT ( 21 ) , RMTX ( 21 , 21 ) , SO ( 9 ) , SP ( 9 ) , DZRO ( 9 ) . NEXT ( 9 ) . COF ( 9 )

,

X ECC ( 9 ) , DELR ( 9 ) , SNYNEU ( 23 ) , SINVAL ( 10 ) , COl . C02 , VEL (9,21), IPRINT COMMON RCOSO . RCOSP . RCOS . DGRUPS ( 23 ) , VELO (9,21), VELP (9.21) X . VEL1 ( 5 ) , SCALE . SWTS (23,21), WTSD (23,21)

COMMON IDOSE , HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) , SINR ( 23 ) , THMLN ( 23

)

X , ENGAM ( 21 ) , DIMEN . SCOBLT ( 23 ) , AWTS (6.5), SGAM ( 23 ) , DGAM ( 23 ) , AMOD (5.5)

DIMENSION VELX ( 21 ) , SPCTRM ( 23 ) , DFLXN ( 21 ) , DELO ( 9 ) . DELP ( 9 )

,

X DOUT ( 9 ) , REDUXT ( 9 ) , REDUXD ( 9

)

1 FORMAT (1H ) 2 FORMAT (1H1)

3 FORMAT ( 15 , 7 ( 1PE11 . 3 )

)

4 FORMAT (55H IC DELP DELO DELX RCOSO RCOSP RCOSX ) 5 FORMAT ( 5 5H NEXT(JC) CDIN CDOUT DIN DOUT REDUXT (JC) REDUXC ) 7 FORMAT (9H NEXT ( I

) = , 12 . 7H DOSE= , 1PE10 . 3 . 14H REDUX FACTR= , 1PE10 . 3

)

8 FORMAT ( 15 , 5 ( 1PE10 . 3 )

)

9 FORMAT ( 5 5H MEG, N SOURCE SPCTRM DFLXN VELO VELP ) 10 FORMAT ( 15 , 4 ( 1PE10 . 3 ) . 110 . 4 ( 1PE10 . 3)

)

11 FORMAT ( 55H VEL(N)/IC=1 2 3 4 5 ) 12 FORMAT ( 15 , 10 ( 1PE10 . 3 )

)

13 FORMAT (5H SO= . 10 ( 1PE10 . 3)

)

14 FORMAT (5H SP= , 10 ( 1PE10 . 3 )

)

15 FORMAT (7H ARTO= .1PE10.3)

RCOSO=2 . / ( 1 . +SO ( IC-1 )

)

RCOSP=2 . / ( 1 . +SP ( IC-1 )

)

RCOSIN=2 . / ( 1 . +COF ( 1 )

)

CALL RETRN ( VELX , SX , IC)

RCOSX=2 . / (1 . +SX)

CDIN=0

JC=IC-1

DO 19 1=1, JC DELO ( I ) =0

19 DELP ( I

) =0

DELX=0

CDR1=0

NMAX=21

79
FILE 3 Continued

331 IF(IDOSE. LT. 5) GO TO 20

332 NMAX=10

333 20 CDIN=DWTS ( 1

)

334 IF (CDIN . EQ . 0 . ) CDIN=1

.

335 U=l.

336 CRD1=CRD1+VEL ( IC-1 , 2

) *DWTS ( 2 ) /U

337 DO 100 NREF= 1 . NMAX

338 DO 99 1 = 1. JC

339 DELO ( I

) =DELO ( I ) +VELO ( I . NREF ) *DWTS ( NREF ) /U

340 99 DELP ( I

) =DELP ( I ) +VELP ( I , NREF ) *DWTS ( NREF ) /U

341 DELX=DELX+VELX ( NREF ) * DWTS ( NREF ) /U

342 100 U=U*0.6y6

343 DELX=DELX* SCALE

344 DO 213 M=1 . 23

345 SPCTRM(M) =0

.

346 U=l.

347 V=l.

348 DO 212 N=1 , NMAX

349 DFLXN (N) = ( VELO ( JC . N

) *RCOSO+VELX ( N) *RC0SX ) *V/U

350 SPCTRM (M) =SPCTRM (M) + ( SWTS (M.N)*1.E10) *DFLXN (N ) /V

351 V=V* . 75 352 212 U=U* . 696

353 213 SPCTRM (M) =SPCTRM(M) /I .E10

3 54 IF (NEXT (1) ) 101.101,102

355 C

3 56 102 ZD=DELR ( 1

)

357 EXF= ( ,6715+ZD*( .924+.565*ZD) )/( ,0685+ZD*( . 5445+ . 137*ZD )

)

358 CDIN=CDIN*EXP ( ZD*EXF

)

359 C 360 101 DIN=CDIN*RCOSIN

361 CDOUT=DELO ( JC ) +DELX

362 DO 103 1=1, JC

363 DOUT ( I ) =DELO (I) *2 . / (1 . +SO (I) ) -s-DELP (I) *2 . / (1 . +SP(I) ) 3M REDUXT ( I

) =DOUT ( I ) /CDIN

3^5 103 REDUXD ( I

) =DOUT ( I ) /DIN

3t 6 REDUXC=CDOUT/CDIN

3 1 7 IF ( IPRINT . EQ . 0 ) GO TO 1000

36 8 IF ( IPRINT . EQ . 1 ) GO TO 106

36 9 DO 104 1=1, JC

370 WRITE (6,7) NEXT ( I ) , DOUT ( I ) . REDUXD ( I

)

371 104 CONTINUE

372 GO TO 107

373 106 WRITE (6, 7) NEXT( JC) ,DOUT( JC) , REDUXD (JC)

374 107 IF (NEXT ( JC+2 ) . EQ.8) GO TO 1000

375 WRITE (6, 1) 376 IF ( IPRINT. GT. 2) GO TO 105

3. 7 GO TO 1000

378 105 WRITE (6,1)

3/9 WRITE (6, 15) SCALE

360 WRITE (6. 9) 381 DO 109 1=1.21

382 WRITE (6,8)1, SGRUPS ( I ) , SPCTRM ( I ) . DFLXN ( I )

,

VELO ( JC , I ) , VELP ( JC , I

)

383 109 CONTINUE

384 MEG=22

385 WRITE (6.8) MEG , SGRUPS ( MEG ) . SPCTRM ( MEG

)

386 MEG=23

387 WRITE (6,8) MEG . SGRUPS ( MEG ) . SPCTRM ( MEG

)

388 WRITE (6.1)

389 WRITE (6. 11)

390 DO 113 N=1 , 21 391 WRITE (6,12) N. (VEL(ID.N) ,ID=1, JC)

392 113 CONTINUE

393 WRITE (6,13) ( SO ( ID) , ID=1 , JC)

394 WRITE (6,14) (SP(ID) ,ID=1. JC)

395 110 WRITE (6,1)

396 IF ( IPRINT. LE. 3) GO TO 1000

397 WRITE (6, 4)

398 WRITE (6,3) IC , DELP . DELO . DELX , RCOSO . RCOSP

,

RCOSX

399 WRITE (6. 5) 400 WRITE (6.3) NEXT ( JC ) , CDIN , CDOUT . DIN , DOUT ( JC ) , REDUXT ( JC ) . REDUXC

401 1000 RETURN

402 END

80
FILE 4

1 :* 2 : 3:4:

5:6 : 7:

8:

9: 10 : 11 : 12: C

13:

14:

15:

16:

17:

18:

19:

20:

21: 22 : 23:

24:

25:

26:

27: 28:

29: 30:

31:

32:

33:

34:

35: 36:

37: C

38:

39: 40:

41: 42:

43:

44:

45:

46:

47:

48:

49:

50:

51:

52:

53:

54: 55:

56:

57:

58:

59: 60:

61:

62: C

63:

64:

65: 66 : 67: 68 : 69:

70:

71: 72:

73:

74: 75:

76:

77:

78: 79: C

80:

81:

82:

SUBROUTINE RETRN (VELX, SX, IC)

COMMON DIAMS ( 18 ) . DUCTO ( 3 . 17 , 6 ) . DUCT1 ( 3 . 17 . 6 ) . DUCT2 (3.17,6) .IEL.LP. X BEND (10.6). SGRUPS ( 23 ) . DWTS ( 21 ) . CRRNTO (9.10). RCRRNT (9.10)

COMMON RVCT ( 21 ) . RMTX ( 21 . 21 ) . SO ( 9 ) . SP ( 9 ) . DZRO ( 9 ) . NEXT ( 9 ) . COF ( 9 )

,

X ECC ( 9 ) . DELR ( 9 ) . SNYNEU ( 23 ) , SINVAL ( 10 ) . COl . C02 . VEL (9.21). IPRINT COMMON RCOSO . RCOSP . RCOS . DGRUPS ( 23 ) . VELO (9,21), VELP (9,21) X , VEL1 ( 5 ) . SCALE . SWTS (23,21), WTSD (23.21)

COMMON IDOSE , HNDRSN ( 23 ) , SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) . THMLN ( 23

)

X . ENGAM ( 21 ) , DIMEN , SCOBLT ( 23 ) . AWTS (6,5), SGAM ( 23 ) . DGAM ( 23 ) . AMOD (5,5)

DIMENSION VELX ( 21 ) , VELPX ( 21 ) . SVELP ( 21 ) , SVELO ( 21 ) . XVELO ( 21 )

X XVELP ( 21 ) . SVEL ( 21

)

1 FORMAT (1H ) 2 FORMAT ( 13 , 11 ( 1PE10 . 2 )

)

3 FORMAT ( 13 , 7E13 . 3

)

KSUP=DZRO ( IC

)

DO 8 IR=1 , 21 SVEL ( IR ) =VEL ( IC- 1 , IR

)

VELPX ( IR ) = VELP ( IC-1 , IR)

8 VELX ( IR )= VELP (IC-l.IR)

SX=SP(IC-1)

SXP=SP(IC-1)

IF(KSUP) 50.50.10

10 RDLNTH=COF ( IC) IF ( RDLNTH . GT . 5 . ) GO TO 50 DLNGTH=COF ( IC-1

)

IF ( RDLNTH. LT. DLNGTH + 1. ) GO TO 50 IF(RDLNTH.GT.DLNGTH+5. ) GO TO 50 TOTO=0

T0TP=0

DO 12 NREF=1 , 21 TOTO=TOTO+VELO ( IC-1 . NREF)

12 TOTP=TOTP+VELP( IC-1, NREF)

RC= (TOTO*RCOSO+TOTP*RCOSP) / (TOTO+TOTP)

IF (RC-1 . ) 13. 13. 14 13 RC=1 . 01

14 STOTO=0

ST0TP=0

XTOTO=0

XT0TP=0

SSP=0.

XSP=0.

XS0=0.

SO(IC+1)=0.

SP(IC+1)=0.

26 DO 27 IR=1 , 21 RVCT ( IR ) =VEL ( IC-2 , IR)

SVELO ( IR )= VELO ( IC- 1 . IR

)

27 SVELP (IR)=VELP( IC-l.IR)

SSO=SO(IC-l)

SSP=SP(IC-1)

IEL=NEXT ( IC-1

)

LP = 2 GO TO (100.110,300.500,500) ,IEL

100 CALL DCTISO (RDLNTH, IC-1) GO TO 30

110 CALL DCTSLT (RDLNTH, IC-1) GO TO 30

300 DLNGTH =0

500 IF ( IEL . GE . 4 ) CALL BT0DA( RDLNTH . IC-1

)

30 RISO=RDLNTH-DLNGTH

RSLT=RISO- . 75 LP=3

RVCT ( 1

) =0

DO 35 IR=1 , 20

35 RVCT (IR+1)=0.53*VEL( IC-l.IR)

CALL DCTISO (RISO.IC)

DO 40 IR=1 , 21 XVELO ( IR ) = VELO ( IC . IR

)

XVELP ( IR) = VELP (IC.IR) IF ( IR . GT . 20 ) GO TO 41 40 VELO (IC-1, IR+1) = .06*VEL( IC-1. IR) 41 XSO=SO(IC)

XSP=SP(IC)

SO(IC-1)=0.

CALL DCTSLT ( RSLT. IC)

DO 45 IR=1 . 21 VELO (IC-l.IR) = SVELO ( IR

)

VELP ( IC- 1 . IR ) = SVELP ( IR

)

81
FILE 4 Continued

83: VEL(IC-1.IR)=SVEL(IR)

84 : VELX ( IR ) =XVELO ( IR ) +XVELP ( IR ) + VELO ( IC , IR ) + VELP ( IC . IR ) +VELPX ( IR

)

85: XTOTO=XTOTO+XVELO ( IR) 86: XTOTP=XTOTP+XVELP(IR)

87:. STOTO=STOTO+VELO ( IC, IR) 88: 45 STOTP=STOTP+VELP ( IC . IR) 89: SO(IC-l)=SSO

90: SP(IC-1)=SSP

91 : TOTX= TOTP + STOTO + STOTP +XTOTO + XTOTP

92 : RCOS = 2

. * ( TOTP/ ( 1 . +SXP ) +XTOTO/ ( 1 . +XSO ) + XTOTP/ ( 1 . + XSP

)

93: X +STOTO/ ( 1 . +S0 ( IC) ) +STOTP/ ( 1 . +SP ( IC) ) ) /TOTX

94: SX=2 . /RCOS-1

.

95: C

96: IF ( IPRIMT . NE . 4 ) GO TO 50 97: 47 WRITE (6,1)

98: WRITE (6.2) IEL. (XVELO(IR) . IR=1« 21. 2) 99: WRITE (6,2) IC, ( VELPX(IR) . IR=1 - 21

.

2 ) 100: WRITE (6,2) IC, (XVELP (IR) , IR=1 , 21 . 2) 101: WRITE (6,2) IC, (VELX(IR) . IR=1, 21, 2) 102: WRITE (6,3) IC. SX, SSP, XSO, XSP , SO ( IC) , SP ( IC) , DLNGTH

103: WRITE (6, 3) IC, TOTP , XTOTO. XTOTP , STOTO , STOTP, TOTX .RDLNTH

104: WRITE (6.1)

105: 50 RETURN

106: END

107: SUBROUTINE WTLISTS ( JSTT, JXWT. IOX, NMAX, DWTS , SWTD , ENGAM, IDOSE,

108: X ISORS.LP)

109: DIMENSION DWTS (21) , ENGAM ( 21 ) , SWTD (21) . IXFILE ( 65 , 3 ) . WTFILE ( 65 . 21 , 3

)

110: 4 FORMAT (315. 3112)

111: 5 FORMAT (31 5) 112: 8 FORMAT ( 5 ( 1PE12 . 4 )

)

113: C - - - - -SETS UP AND READS PRECALCULATED WEIGHTS FROM A FILE (WTSFILE)

114: ITOT=65

115: 15 IF (IOX) 20.17.100

116: 17 IF ( IXALT) 80,1000.50

117: C (MAIN PROGT^AM: IOX=ql)

118: C - - - - - SET UP 3 INDEX FILES WITH TIE INTERNAL WTFILES

119: 20 DO 40 I=l.ITOT

120: READ (7.5) ( IXFILE ( I . J) . J=1 , 3) 121: DO 40 IX=1 , 21

122: READ (7,8) (WTFILE ( I. IX, J) . J=1 . 3

)

123: 40 CONTINUE

124: IOX=Q

125: IXALT=1

126 : REWIND 7 127: C - - WRITE (4, 4) ( IXFILE ( 12 , J) , J=1.3)

128: C - - WRITE (4,8) (WTFILE (12 , 2 , J) , J=1.3)

129: GO TO 1000

130: C - - - - - WRITE BACK THE INFORMATION INTO PREWTS

131: 100 DO 110 1=1 , ITOT

132: WRITE (7.5) ( IXFILE(I.J) ,J=1,3)

133: DO 110 NX=1 , 21

134: WRITE (7.8) ( WTFILE (I, NX. J) . J=l,3)

135: 110 CONTINUE

136: REWIND 7 137: GO TO 1000

138: C (IOX VALUE OUTPUT DOESN'T MATTER)

139: C - - - - - NEXT EITHER RIAD THE OLD WEIGHTS OR STORE NEW ONES

140: C ( JSTT=ql . CALC THEN WRITE: JSTT=0=ISTT, READ THEN SKIP)

141: 50 IXALT=- IXALT

142: JXWT=IDOSE»5*ISORS+20*LP-25

143: IF ( JSTT) 1000,51.1000

144: 51 JSTT= IXFILE ( JXWT, 2) 145: IF ( JSTT . NE . 0 ) GO TO 1000

146: 52 JMAX=IXFILE( JXWT. 3) 147: DO 55 NX=1.JMAX

148: DWTS ( NX ) =WTFILE ( JXWT . NX . 1

)

149: SWTD ( NX ) =WTFILE ( JXWT . NX , 2

)

150: 55 ENGAM ( NX )= WTFILE (JXWT. NX. 3) 151: NMAX= JMAX

152: GO TO 1000

153: C - - - - - WRITE WEIGHTS JUST COMPUTED

154: 80 IXALT = - IXALT

155: IF ( JSTT . EQ . 0 ) GO TO 1000

156: IXFILE ( JXWT. 2) =0 157: IXFILE ( JXWT. 3) =NMAX

158: DO 85 NX=1.NMAX

159: WTFILE ( JXWT . NX , 1 ) = DWTS ( NX

)

160: WTFILE (JXWT, NX, 2) = SWTD (NX) 161: WTFILE ( JXWT , NX , 3) = ENGAM (NX) 162: 85 CONTINUE

163: 1000 RETURN

164: END

165: SUBROUTINE ENTDAT 82
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

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

COMMON DIAMS ( 18 ) . DUCTO ( 3 . 17 . 6 ) , DUCT1 ( 3 . 17 , 6 ) . DUCT2 (3.17,6) .IEL.LP. X BEND (10.6). SGRUPS ( 23 ) , DWTS ( 21 ) . CRRNTO (9.10). RCRRNT (9.10)

COMMON RVCT ( 21 ) . RMTX ( 21 . 21 ) . SO ( 9 ) . SP ( 9 ) . DZRO ( 9 ) . NEXT ( 9 ) . COF ( 9 )

.

X ECC ( 9 ) . DELR ( 9 ) , SNYNEU( 23) , SINVAL ( 10 ) , COl , C02 , VEL (9,21) . IPRINT COMMON RCOSO . RCOSP , RCOS . DGRUPS ( 23 ) , VELO ( 9 . 21 ) . VELP ( 9 . 21

)

X , VEL1 ( 5 ) .SCALE. SWTS (23. 21) ,WTSD(23,21)

COMMON IDOSE , HNDRSN ( 23 ) . SEPI ( 23 ) , STH ( 23 ) . SINR ( 23 ) . THMLN ( 23

)

X,ENGAM(21) , DIMEN. SCOBLT ( 23) ,AWTS(6. 5) .SGAM(23) ,DGAM(23) . AMOD ( 5 . 5

)

DIMENSION TPX1 ( 17) . TPX2(17). TPX3(51), TPX4(51), TPX5(17).

X TPX6 ( 51 ) , TPX7 ( 17 ) , TPX8(51). TPX9(51), TP10(17), TP11(51).

X TP12(17) . TP13 ( 51 ) , TP14 ( 51 ) , TP15(17). TP16(51). TP17( 8). X TP18 ( 8). TP19 ( 8). TP20 ( 8), TP21 ( 8). TP22( 8). TP23(23).

X TP24(23) TP25 ( 23 ) , TP26(23). TP27(20). TP28(10). TP29(23).

X TP30 ( 23 ) , TP31 ( 23 ) . TP32(23). TP33(23)

C DATA TPX1 /0... 5,1. ,1.5. 2. .2. 5. 3. .3. 5. 4. .4. 5. 5. .5. 5. X 6. .6.5.7. .7.5.8./

DO 101 IL=1 , 17

101 DIAMS (IL)=TPX1(IL)

C DATA TPX2 /l. , .413, .205. .111. .0685. .0481. .0354. X .0292, .0226, .0191. .0155. .0127. .0108. .0091. .0077. .0065. .0055/

DATA TPX3 /. 35, .56. .66, .71, .74. .75. .76. X .753. .745. .73. .715, .698. .675. .655. . 635 . . 615 , . 60 . . 35 . .488. .595. X .675. .72. .747. . 76 . . 753 . . 745 . . 73 , . 715 , .698. . 675 . . 655 , . 635 . .615. X .60, .52, .555. .61, .725. .77. .79. .804. .80. .79. .775. .758. .733. .71. X .685. .66. .64. .62/

DATA TPX4 /. 30. .59. .75. .85. .93, .98.1.02.

X 1.04.1.05.1.05.1.05,1.05.1.05,1.05,1.05,1.05.1.05. .30. .705. .835, X .91, .975.1.02,1.05.1.06.1.07,1.08,1.08.1.08,1.08,1.08,1.08.1.08.

X 1.08. .74. .825. .90. .965.1.015.1.05,1.06,1.07.1.08,1.08.1.08.1.08.

X 1.08.1.08.1.08,1.08.1.08/

DATA TPX5 /l. . .517, .285. .154. .097. .063. .046. .037. X .028. .022. .019. .016. .013. .0108, .009. .0074. .0062/

DATA TPX3 /. 265. .38, .48, .55, .605, .65. .69. X .725, .75. .78, .80, .823, . 84 , . 85 . . 86 . . 87 . . 88 . .265. . 365 . . 46 . . 53 , . 58

,

X .625. .66 .693. .72. .745, . 765 . . 785 . . 8 . . 82 , . 83 , . 84 . . 85 . . 40 . . 40 . . 43

.

X .49, .54, 584, .62, .65. .675, .70. .72. .74. .75. .768. .78. .79. .80/

IDX=0

DO 106 J=1 ,

3

DO 106 IL=1 . 17

IDX=IDX+1

IF(J-l) 106.104.105

104 DUCTO ( 1 . IL , 1 ) =TPX2 ( IDX

)

DUCTO ( 1 , IL . 4

) =TPX5 ( IDX)

105 DUCTO (J, IL, 2) =TPX3 (IDX)

DUCTO ( J . IL , 3

) =TPX4 ( IDX

)

DUCTO ( J . IL . 5

) =TPX6 ( IDX

)

106 CONTINUE

C DATA TPX7 /2..1.,1.,1.,1.,1.,1.,1..1.,1.,1.,1.,

X 1..1..1..1..1./

DATA TPX3 /. 35. .62. .52, .51. .51, .512, .518, X .523. .533. .548, .562, .58, .6. .622, .645, . 665 , . 685 . . 35 . . 59 . . 515 , . 51

.

X .51, .512. .518, .523, .533. .548, .562. .58, .6. .622, . 645 . . 665 , . 685

,

X .615. .56. .52. .51. .51. .512, .518. .523, .533. .548. .562. .58. .6. .622. X .645, .665. .685/

DATA TPX9 / . 3 . . 61 . . 83 . . 945 , 1 . 00 , 1 . 04 . 1 . 06

X. 1.065, 1.067. 1.064. 1.06. 1.04. 1.02. 1.00. . 98 . . 96 . . 94 . . 3 . . 68 . . 865

.

X .965.1.025.1.06.1.08.1.084.1.08.1.07,1.05.1.04.1.02.1.00. .98. .96. X .94, .68, .79, .905,1. .1.07.1.11,1.125,1.13.1.12.1.11.1.1.1.07.1.05.

X 1.03.1.01. .99. .97/ DATA TP10 /. 5. .5. .5, .5. .5. .5. .5. .5. .5. .5, .5, .5, X .5, .5. .5, .5, .5/

DATA TP11 /. 35. .48. .58. .64. .69, .735. .766, X .79. .81, .83, .843. .85. .86, .866. .872, .878, .884. .33.-44, .533. .59. .64 X , .68, .71. .733. .75. .765. .775, .78. .786. .791. . 795 . . 8 . . 805 . . 35 . . 42

.

X .485. .545, . 587. .625. . 653 , . 68 , . 7 . . 716 , .728, . 737 . . 745 . . 753 . . 754

.

X .758, .76/

IDX=0

DO 111 J=1 ,

3

DO 111 IL-1,17

IDX=IDX+1

IF(J-l) 111,109.110

109 DUCT1 ( 1 , IL , 1

) =TPX7 ( IDX

)

DUCT1 ( 1 , IL . 4 ) =TP10 ( IDX)

110 DUCT1 ( J . IL . 2

) =TPX8 ( IDX)

DUCT1 ( J . IL . 3

) =TPX9 ( IDX

)

DUCT1 ( J . IL . 5 ) =TP11 ( IDX

)

111 CONTINUE

C DATA TP12 /2. .1.517,1.124.1. .1. .1. .1. .1. .1. ,1.

.
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

288

289

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

319

320

321

322

323

324

325

326

327

328

329

330

FILE 4 Continued

X 1. .1. .1. ,1. .1. ,i. .1-/

DATA TP13 /. 35. .565, .77. .81. .77. .71. .663. X .63. .606. .6, .61. .62. .635. .657. . 68 . . 7 , . 72 , . 35 , . 49 . . 685 . . 78 . . 76

X .712. .663, .63. .606. .6, .61. .62. .635. .657. . 68 . . 7 . . 72 . . 5 . . 61 , . 8

X .825. .804. .74, .685. .657. .64. .64, .643, . 655 , . 67 . . 69 , . 715 . . 745 . .764/

DATA TP14 /. 3. .56. .72, .845. .94,1.01,1.04,

X 1.07.1.085.1.10.1.10.1.09,1.077,1.06.1.04.1.02.1.00, .3. .62. .836, X .91, .97.1.01,1.04.1.07,1.085,1.10,1.10.1.09,1.077.1.06.1.04.1.02.

X 1.00, .767. .83, .907. .966.1.005.1.04.1.06,1.08.1.095,1.10.1,10,

X 1.09.1.077.1.06.1.04,1.02.1.00/

DATA TP15 /. 5, .5, .5. .5, .5, .5, .5. .5. .5. .5, .5, .5, X .5. .5. .5. .5, .5/

DATA TP16 /. 65. .57, .54, .568. .645, .73. X .78. .78, .78. .78, .78. .78. .78. .78, .78. .78. .78. .485. .456. .468. .513. X .59. .655. .675. .675. .675. .675. . 675 , . 675 . . 675 . .675, .675. .675. .675, X .41, .41. .43. .476. .536, .584, .607. .613. . 62 , . 62 . . 62 . . 62 , . 62 . . 62 . . 62

X .62. .62/

IDX=0

DO 116 J=l,3

DO 116 IL=1 . 17

IDX=IDX+1

IF(J-l) 116.114,115

114 DUCT2 ( 1 , IL, 1 ) = TP12 ( IDX)

DUCT2 ( 1 , IL , 4

) =TP15 ( IDX)

115 DUCT2 ( J . IL . 2 ) =TP13 ( IDX

)

DUCT2 ( J . IL . 3 ) =TP14 ( IDX

)

DUCT2 ( J, IL, 5

) =TP16 ( IDX

)

116 CONTINUE

C DO 21 K=1 ,4,3

DO 20 J=2 .

DO 20 IL»1 , 17

DUCTO ( J , IL . K ) = DUCTO ( 1 , IL , K

)

DUCT1 ( J, IL. K) =DUCT1 ( 1 , IL, K)

DUCT2 ( J , IL . K ) = DUCT2 ( 1 , IL . K

)

20 CONTINUE

21 CONTINUE

C DATA TP17 /-. 4, -.2,0. .0.2. 0.4. 0.6. 0.8.1./

DATA TP18 /. 152. .169. .189. .213, .243. .290. .361, .450/

DATA TP19 /. 241. .248, .257, .266, .275. .277, 279, .293/

DATA TP20 /. 41. .41. .41, .41, .41. .41, .41, .41/

DATA TP21 /. 988. .983. .980,1.105. .955. .9487,1.039,1./

DATA TP22 /. 5. .5. .5. .5. .5. .5. .5. .5/

DO 22 IL=1,

SINVAL ( IL) =TP17 ( IL)

BEND ( IL, 1

) =TP18 ( IL) BEND ( IL , 2 ) =TP19 ( IL

)

BEND ( IL . 3 ) =TP20 ( IL

)

BEND ( IL . 4

) =TP21 ( IL) BEND ( IL , 5 ) =TP22 ( IL) 22 CONTINUE

C DATA TP23 /0 ., 0 ... 0215 ,. 08 ,. 1405 .. 1341 ,. 2673 .. 2646 .. 0944

*, .3038, .443, .3922, .4468, .2464, .235. .2107, .2099, .1492. .1192. .1497.

* .1202. .0957,0./

DATA TP24 /0.,0.,0.,0.,0.,0.,0..0.,0.,0.,0..0..0..

0..0..0..0..0..0..0..0..0..1./

DATA TP25 /I . 898E+6 , 2 . 204E+6 . 3 . 777E+6 , 7 . 885E+6 , 1 . 05E+7

* 1.023E+7.1.027E+7.1.351E+7.5.54E+6.2.316E+7.4. 544E+7 , 1 . 054E+8

* 2.215E+8,2.E+8,1.74E+8.1.549E+8,1.538E+8,1.071E+8,8.454E+7,

* 1 . 009E+8 , 7 . 236E+7 , 6 . 108E+7 . 5 . 411E+7/

DATA TP26 /2 . 326E+8 . 7 . 325E+7 . 8 . 596E+7 . 4 . 877E+7 . 5 . 477E+7

* 4.258E+7,5.011E+7,6.085E+7,4.762E+7.6.542E+7.7.11E+7,7.11E+7,

* 9,493E+7,1.813E+8,5.777E+8.8.534E+8.1.398E+8,6*0./

DATA TP33 /0 . . 0 . , 0 . . 0 . , 0 . . 0 . . 0 . . 1 . 0 . 15*0 .

/

DO 26 M=1 , 23

SEPI (M) =TP23 (M)

STH (M) =TP24 (M)

SINR (M) =TP25 (M)

SGAM (M) =TP26 (M)

SCOBLT(M) =TP33 (M) 26 CONTINUE

C DATA TP27 /. 4359 ,. 324 .. 1566 .. 0722 ,. 0113 , 0 .. 1 .

* 3*0. , .1, .25. .25, .25. .15.0. ,1. ,3*0./ DATA TP28 / . 9 . . 65 , . 4 . . 15 . 0 . . 1 . , . 9 . . 65 . . 4 , . 15/ IDX = 0 DO 28 1=1.4

DO 28 J=1 ,

IDX=IDX+1

AWTS ( I . J) =TP27 ( IDX) IF ( 1-2 ) 27,27.28
FILE 4 Continued

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

27 AMOD ( I , J) =TP28 ( 1DX)

28 CONTINUE

DATA TP29 /7 . E-9 . 7 . E-9 . 7 . 08E-9 . 6 . 72E-9 . 6 . 03E-9

*5 . 43E-9 . 4 . 83E-9 . 4 . 48E-9 . 4 . 33E-9 . 4 . 23E-9 . 3 . 96E-9 . 3 . 3E-9 . 1 . 73E-9

* 7.E-10,6.5E-10,6.07E-10,6. 72E-10 , 5.35E-10,3.88E-10,3.42E-10,

* 3 . 27E-10 , 3 . 22E-10 , 3 . 2E-10/ DATA TP30 /5.46E-9.5.13E-9.4.84E-9.4.61E-9.4.44E-9.

* 4 . 13E-9 ,4.01E-9,3.39E-9,3.15E-9,3.09E-9,2.64E-9,1.97E-9,1.12E-9,

* 0 . 229E-9 , 9*0 .

/

DATA TP31 /22*0 . , 1 .

/

DATA TP32 / 1 . 66E-9 , 1 . 46E-9 . 1 . 27E-9 , 1 . 08E-9 , 8 . 75E-10

,

* 7 . 35E-10 , 6 . 44E-10 , 5 . 3E-10 , 4 . 45E-10 , 3 . 5E-10 , 2 . 56E-10 , 2 . 15E-10

,

* 1 . 77E-10 , 1 . 22E-10 , 6 . 6E-11 , 3 . 9E-11 , 8 . 37E-11 , 6*0 .

/

DO 32 M=1 , 23

SNYNEU(M) =TP29 (M) HNDRSN ( M ) =TP30 ( M

)

THMLN(M) =TP31 (M)

32 DGAM(M) =TP32(M)

RETURN

END

85
NBS-114A (rev. 2 »ec) UcS. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See ; nstruction s)

1„ PUBLICATION OR

REPORT NO.

NBSIR 87-3534

2. Performing Organ. Report No^ 3. Publication Date

MARCH 1987

4. TITLE AND SUBTITLE

OCTDOS : NEUTRON AND GAMMA PENETRATION IN COMPOSITE DUCT SYSTEMS

5. AUTHOR(S)

L. V. Spencer

6. PERFORMING ORGANIZATION (If joint or other than N BS, see in struction s)

national bureau of standards

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Contract/Grant No,

8 . Type of Report 8. Period Covered

9, SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS ( Street. City, State , ZIP)

10. SUPPLEMENTARY NOTES

j_T Document describes a compjter program; SF-185, FIPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most significant in formation. If document includes a significant

hi HI ioeraphy or literature survey, mention it here)

This paper describes computer methods for estimating neutron and

gamma ray fl uence rate, dose, and even spectral features due to

penetration through a series of duct segments. The procedure 1-nks

together data for individual segments -- straight sections and bends

in arbitrary combinations; and the resulting composite can include

computations for a room at the end, if there is one. This particular

method was developed for rapid estimates in the protection problems

against nuclear weapons; but the concepts which are employed are more

broadly applicable.

12. KEY WORDS (Six to twelve entries; alph a bet/ cal order; capi ta I i ze only proper names; and separate key words by semi colon s) computer software; gamma-ray albedos; gamma-ray penetration; multiple gamma-ray

albedos; multiple neutron albedos; neutron penetration; radiation in ducts;

weapons gamma-rav shielding: weapons neutron shielding

13. availability

'jcj Uni imited

S For Official Distribution. Do Not Release to NTIS

Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. “ 20402.

Order From National Technical Information Service (NTIS), Springfield, VA. 22161

14. NO. OF

PRINTED PAGES

90

15. Price

$13 . 95

USCOMM-OC 6043-P80

0f