THE U N I V E R S I T Y OF MICHIGAN COLLEGE OF ENGINEERING Department of Nuclear Engineering Technical Report SOLUTIONS OF THE TWO-GROUP TRANSPORT EQUATION IN PLANE GEOMETRY Dale R. Metcalf Paul F. Zweifel, Project Director ORA Project 01046 supported by: NATIONAL SCIENCE FOUNDATION GRANT NO. GK- 1713 WASHINGTON, D.C. administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR February 1968

ACKNOWLEDGMENTS To my Doctoral Committee Chairman, Professor Paul F. Zweifel, I wish to express my sincere gratitude for his guidance and supervision. In particular, at the critical times of selecting a dissertation topic and of organizing the final report, Professor Zweifel provided outstanding assistance. Also I wish to express my appreciation to both Professors Paul F. Zweifel and K. M. Case for the opportunity to assist in the editing of their recent book on neutron transport theory. My real education in transport theory commenced with this undertaking. I express my appreciation to the other members of my doctoral committee, Professors K. M. Case, A. Z. Akcasu, G. C. Summerfield, and F. C. Shure. The financial support of the U.S. Atomic Energy commission through its fellowship program (1964-67) is gratefully acknowledged. The National Science Foundation provided research support during the summer of 1967 and during the final preparation and reproduction of this manuscript. For this assistance, I express my gratitude. I wish to express my appreciation to Dr. H. L. McMurry of Idaho Nuclear Corporation for his professional leadership and counsel in advising me to improve my education and for his continued interest in my progress. In conclusion, I wish to express my appreciation to my wife and family for their patience, understanding, and support. ii

TABLE OF CONTENTS Page LIST OF TABLES v LIST OF FIGURES vi ABSTRACT vii Chapter I. INTRODUCTION 1 II. EIGENFUNCTIONS AND EIGENVALUES 6 III. TWO-GROUP FULL-RANGE COMPLETENESS 15 IV. TWO-GROUP HALF-RANGE CASE 19 V. APPLICATIONS OF TWO-GROUP SINGULAR EIGENFUNCTION EXPANSIONS 29 A. Introduction 29 B. Milne Problem 29 C. Constant Isotropic Source Problem 36 VI. NUMERICAL METHODS AND COMPUTER SOLUTION TO HALF-SPACE PROBLEMS 41 A. Program for Finding Number, Type, and Value of the Discrete Eigenvalues 41 B. Discussion of Mesh Spacing 42 C. Numerical Evaluation of Singular Integrals 43 D. Iteration Procedure for the Expansion Coefficients 47 VII. SPHERICAL HARMONIC SOLUTIONS TO THE CONSTANT SOURCE AND MILNE PROBLEMS 48 VIII. NUMERICAL RESULTS AND COMPARISONS 67 IX. CONCLUSIONS 85 REFERENCES 87 iii

TABLE OF CONTENTS (Concluded) Appendix Page A. THE EVALUATION OF y(4) 89 B. NECESSARY MODIFICATIONS FOR TWO PAIRS OF DISCRETE ROOTS 93 C. AN ALTERNATE SOLUTION METHOD TO EQ. (4.29) 97 D. COMPUTER CODE LISTINGS 100 iv

LIST OF TABLES Table Page 2.1 The Zeros of the Dispersion Function 12 8.1 Two-Group Macroscopic Cross Sections 67 8.2 Eigenvalues of the Two-Group Transport Equation 69 8.3 Extrapolation Distances for Milne Problems 72 8.4 Total Fluxes for Group 1, Set I 82 8.5 Total Fluxes for Group 2, Set I 82 8,6 Total Fluxes for Group 1, Set IV 83

LIST OF FIGURES Figure Page 8.1. Expansion coefficients for Milne problem, Set I. 70 8.2. Expansion coefficients for Milne problem, Set IV. 71 8.3. Constant source angular distribution for Group 1, *1(0,G), Set I. 74 8.4. Constant source angular distribution for Group 2, *2(0,~), Set I. 75 8.5. Constant source angular distribution for Group 2, *2(1,1), Set I. 76 8.6. Modifications in exact angular distribution with distance into the medium, Set I. 77 8.7. Absorption effect on exact exit angular distribution for constant source problem, Group 2, Sets I and IV. 78 8.8. Exact total flux for the constant source vs. distance into the medium, Group 2, Set I. 80 8.9. Total flux for the constant source vs. distance into the medium, Group 2, Set IV. 81 8.10. Constant source current vs. distance for Group 2, Set IV. 84 vi

ABSTRACT The main objective of this thesis is to obtain exact solutions to the twogroup neutron transport equation utilizing the singular eigenfunction approach. By this method exact solutions to a limited set of problems can be obtained against which the accuracy of conventional approximations may be tested. A detailed study was made of the singular eigenfunctions of the two-group transport equation. A convenient linear combination of the continuum eigenfunctions was selected to simplify later analysis. Using these eigenfunctions, a full-range completeness proof for infinite medium applications is presented. For the half-range case, reliance is placed on a half-range theorem proved several years ago and (with this as authority) the general half-range problem is rederived in a form convenient for numerical analysis. Explicitly, two coupled equations are obtained where the unknowns are the discrete and one continuum mode expansion coefficients. Finally, the remaining continuum mode expansion coefficient is expressed directly in terms of the other two coefficients. This general result is then applied to the solution of two typical half-space problems, the Milne and constant source problems. For comparison purposes the solutions of the two-group Milne and constant source problems are also obtained in the P1, P3, and DP1 approximations. Special numerical methods developed for the computer solution to the exact problem are discussed in detail. Of particular importance is the development of a relatively fast and accurate method for calculating singular integrals. Computer codes are written in the Michigan Algorithm Decoder (MAD) language and a complete listing is given. With these codes anyone can (with different cross section sets) test the accuracy of conventional approximation methods. Calculations are performed for four light water systems with varying amounts of absorption. These calculations show that although the P1 approximation yields fairly good results for the total flux and current, the angular distribution is not so well represented especially at the interface. The P3 approximation continues to give negative inward angular fluxes but of a reduced magnitude. The DP1 method improves significantly on the P1 and P3 angular fluxes as well as on the total flux and current. The boundary condition at the vacuum-medium interface can be specified exactly in the DP1 approximation; therefore, the exit angular distribution at the boundary is well represented. An interesting result from this study was in the constant source case the exit angular distribution becomes less pronounced in the forward direction with increased absorption in the medium. This is in contrast to the Milne problem where the opposite effect is noted. Also of interest is the fact that the P1 approximation in two-group theory yields two eigenvalues in the convii

tinuum (as well as two outside the continuum) while in the one-group case both eigenvalues are outside the continuum. This means that in two-group theory more than just the asymptotic solution is obtained even in the lowest order approximation. In two-group neutron transport theory, the P3 and DP1 approximations present essentially the same computational difficulty. We conclude that in reactor physics design analysis (with reactor parameters similar to those studied here) the DP1 approximation is much to be preferred over the P3. viii

I. INTRODUCTION Essentially all of the (static) physical properties of a nuclear reactor are determined by a knowledge of the velocity and space dependent neutron flux. This quantity obeys to a very high approximation the neutron transport equation.l In its full generality, this equation is too difficult to solve even for the present generation of digital computers~ For this reason various approximation techniques have been developed through the years to determine the flux in nuclear reactor systems.2 The usual criterion for assessing the accuracy of the calculations has been to compare the results with critical experiments.3 The three main sources of error in these comparisons are the cross sections used, the geometric model of the reactor, and the simplification of the transport equation to a tractable approximate form~ The calculations have often given results that compare quite well with critical experimentsp3 But there remains a certain amount of justified skepticism in the calculations when one attempts to study a new and different system on which no experimental data exist~ Furthermore, the current requirements for economic performance of power reactors entail a great deal more than a simple knowledge of the static criticality factor. It has been well-documented that calculation of temperature coefficients, burnup, etc., require more sophisticated treatments of the transport equation>4 In the absence of detailed experimental information, it is possible to use exact solutions of the transport equation as test cases to assess the validity of various approximation schemes. Needless to say, the physical model 1

2 involved must be grossly simplified, but it is still possible to interpret the accuracy obtainable from the approximate methods when they are then applied to practical systems. The first detailed numerical efforts along these lines were made by Mendelson5 who treated the one-speed model in some detail. Our purpose is to extend these results to the "two-speed" (or as it is generally called the "two-group") approximation. This is the next logical extension to the grossly inadequate treatment of the neutron energy spectrum in a reactor yielded by the one-speed approximation. In fact, many reactors are treated reasonably well with two-group theory.3 For those systems which do require more groups (fast and intermediate systems particularity) the present work serves as a guide for the extension to a multigroup approximation. As a matter of fact, although there has been considerable effort expended in recent years to obtain exact analytical solutions to the transport equation (as is described in detail in Ref. 1), the reduction of these results to numerics has lagged. Obtaining numbers from the analytical formulae is far from a trivial problem, the major reason being the different and unusual nature of the functions involved. This point is clearly made in Mendelson's work.5 Thus, two contributions are made in this thesis. First, the (singular) eigenmodes of the two-group transport equation are studied in some detail and the two-group equations are reduced to a form convenient for numerical analysis4 Then, the numerical analysis is carried through and comparisons are made to standard approximation methods. Some attempt is made to draw

3 general conclusions as to the range of reactor parameters for which the common approximation schemes are valid. The two-group case has been studied previously by Zelazny and Kuszell-6 They succeeded in proving a "completeness theorem" which in effect tells us that the normal modes which we develop are adequate to solve all infinite and semi-infinite medium problems. We depend upon this theorem, but because the results of Ref. 6 are not presented in a form convenient for numerical or analytical computation much of their analysis is redone. Siewert and Zweife17 have obtained analytical solutions to a very special two-group problem, but their restrictions upon the parameters involved are invalid in the neutron transport case. (This work was applied to the problem of radiative transport.) Siewert and Shieh8 have carried out some preliminary work on the more general problem considered here, but their results can'be applied directly only to infinite medium problems. Recently some unpublished work along similiar lines to that described here has come to our attention and we are able to incorporate some of these results into our worko9 Aside from these very little work similiar to that presented here appears in the literature. Most of the effort expended upon the exact energydependent transport equation has been applied to neutron thermalization problems.10 Specifically, we treat the one-dimensional two-group transport equation with isotropic scattering. The method, based on Case's singular eigenfunction approach1'11 is a logical extension of that described in Ref 7.

4 In Chapter II we develop the appropriate eigenvectors and eigenvalues. In Chapter III we present a full-range completeness proof such as is necessary in infinite medium applications. As we have mentioned, similiar theorems have been proved in Refs. 6 and 7, so our treatment is brief, and is presented only because we use a different (and we believe more convenient) notation. In Chapter IV we re-derive the half-space problem in a form convenient for numerical analysis. Zelazny and Kuszell have proved completeness and we do not repeat their work. (Anyway, the completeness theorem is not constructivel2 so its utility is more or less academic.) In our final equation the presence of a nondegenerate Fredholm term prevents the usual Weiner-Hopf factorization. We develop a pair of coupled equations where the unknowns are the discrete and one continuum mode expansion coefficients. These equations can readily be solved with a computer program which is presented in Appendix D. In Chapter V we solve the two-group Milne and constant source problems. Some important simplifications are developed by using the X-function identities listed in Appendix A. Some of the special, unique numerical methods which were developed for the computer solution are discussed in Chapter VI. In particular, a relatively fast and accurate method for the evaluation of singular integrals is presented. In Chapter VII we develop the spherical harmonic solutions to the two-group Milne and constant source problems. These solutions are developed in the Pl, P3, and DP1 approximations.l3 To solve the Milne and constant source problems in these approximations most efficiently, a series of computer programs were written. A listing of these computer programs is given in Appendix Do Some special

problems, which are not of great importance to the overall goals of this thesis, are discussed in Appendices B and C.

II. EIGENFUNCTIONS AND EIGENVALUES The two-group transport equations with isotropic scattering in plane geometry can be written7 at (xl) + -(x4) = 11f rl*(x,)dt + t J2 2(X)d (2.1a) 6 I + a2t2(x,t) = 2 d f1 l(x,p)da + 22 f1 $2(x,4)dt (2.lb) 1(xt) and $2(x,p) are the angular fluxes in groups one and two respectively, al and 02 are the total macroscopic cross sections for each group and the aij are the macroscopic scattering transfer cross sections which describe the transfer of neutrons from group j to group i. All a's are assumed to'be spatially independent. We make the nonrestrictive assumption that y2 < Oa and we divide Eqs. (2,1) by a2. By defining a = al/a2 and using the optical length z = a2x, we write Eqs. (2ol) in the same convenient matrix form as in Ref. 6: L r(z,i) + *(z,) c Cf t(z,4)d4 (2o2) where (/ O\ A(Z\j\) EZ~ =, i (zp4) = C/a C 12= \ C CCS2L )2 2 and C ij *\C21 a22 ~ 2o2

7 As usual we assume solutions of the form (z,) = eZ/ F(,T ). (2.3) This ansatz when substituted into Eq. (2.2), and after cancellation of the space dependence, yields the eigenvalue equation (B is the eigenvalue), A F(nq,) = nC f F(,ll)dL (204a) where crT -t O\ A (2,, b) 0 In contrast to the one-speed case where two eigenvalue spectrum consisting of the discrete eigenvalues of magnitude |1 > 1 and a continuous spectrum in the range re[-l,l], we must here consider three separate regions for the spectrum since in Eq. (2.4b) the term an-> does not vanish in the intervals cj[-1, 1] and Er[ 1,l]. (Recall that a > 1 ) Following Siewert and Zweifel, we consider Region 1 of the continuous spectra as le[- 1, 1]. The eigensolutions in this region are twofold degenerate and will be denoted by F (,) and F2 (,). To obtain linearly independent eigensolutions, we first normalize by fl F (1).4 )d = ((25a -l where Fl_ = 1,(.) (2.5b) F1l(YL

8 We substitute Eq. (2.5a) into Eq. (2.4a) and solve for F( )(n,p). Explicitly F(1)( ~) = - C1_ P + XL(1)+(a-_) (2.6a) and F(1n'O) C21P +X2(\ )F(f-) (2, 6b) n rl-2where ) \X2)l() is to be determined from Eq. (2.5a). By integrating Eqs. (2o6) over: d4 and using Eqo (2o5a),we find,(TI) = 1 2C1- C T(ao) (2.7a) }2(1 ) = -2C21 ('q) (2j7'b) where p1 -u = in (I+a\ = 2 tanh-l(a)l) 2r(T ) D (2o 2c) The symbol P in Eqs. (2.6) indicates that the integrals are to be interpreted in the sense of the Cauchy principal value, With Eqs. (2.7) inserted into Eqs. (2o6),we have C1rP + [12,TCL1.T_-(c)](arT 4) r_1) (T 1) = ~ (2.8) c-PI - 2 CaqT 2(]) T (q-n ) n _

9 A similar procedure with the normalization F2 (lbk)dp = (ON )(2.9) -1 F leads to the result C12T2P - 2C12TT (C)Y (C-kT) (210 F 2. (2.)0) \ 22JP + [l-2nc22T(n)l]5(1T-I) Thus the two linearly independent eigenvectors for pc[- 1, 1] are given -by Eq. (2.8) and Eq. (2.10). Next we consider [-l, - 1] and q,[,1], (Region 2). Defining the eigenvector as F(2)(n,kt) and noting that the normalization is arbitrary, we use J F( )(h, = = (2.11) -1 a2(r) -2C1lrT (l/r) where F2) -Ca)aj(1 + (2C 12a2() (12 (2)(T4 = C2jraj(i)P C22a2()P (212b We integrate each term in Eqsa (212) over d- and use the normalization given by Eq. (2.11) to derive

X(n) = 1 - 2rC227(1) - 2nC11T(l/an) + 4Cn2T (q )T (1/Cn) (2ol3) where C = det C Equations (2.11) and (2.13) are inserted into Eqs. (2.12) to derive the following eigenvector for Region 2: F(2)(~'1) = _ ) (_)k(2 - + b1rr-/2r) (2o14) Here t(n) = C22-2nCT(1//c) and X(n) is given by Eq. (2.13). Next we consider rj[-1,1] and define the associated eigenvector as \ F2(9,4)/ In this case, Eq. (2.4) can be solved directly to give Fj(T,4) = O CLll J F(r,y'T)det + C12 J F2(nq,')dj (2.15a) F2(f,4) = rl f1 Fl(,')d + C22 1 F2')d (215b) We integrate Eqs. (2.15) over d4 to obtain the following secular determinant whose zeros yield the discrete eigenvalues: [ -eC1 (1/an) e2C12ent(1/0 n) = 0,. (20L6) 2C21fT(1/r) l-2C22Tl/rj)

11 By defining the dispersion function O(z), we have ~(z) = 1-2cllzT(l/zz)-2C22zT(l/z)+4Cz2T(l/z)T(l/Z) = 0 (20X7) We have replaced q by the complex variable z to indicate that the roots can be real and/or complex. We can easily see from Eq. (2.17) that if qi is the positive root then -~i is also. And if ~i is a complex root of Eq. (2o17), is also. We have Q)= Q ) = 0 (2.18) O(Ti) = O(n*) = 0 (2.19) 9 ij Zelazny and Baran have made a detailed study of Eq. (2.17) to determine reproduce in Table 2.1 their results which take into account all possibilities For the two-group problems considered in this thesis we have only two discrete eigenvalueso By virtue of Eqso (2.18) and (2.19) these must be real or imaginary, In fact from Table 2.1 we note that the discrete roots are either real and/or purely imaginary. Therefore the eigenvalues occur in ~ pairs and we denote these by +rji with associated eigenvectors Fi~(p). The normalization in Eqs. (2153) is arbitrary. If we use the same normalization as in Region 2 (with Tri -' r),we derive the discrete eigenvectors for j[-l,l1; namely F,(.)= Cl2rt/csfi+) = i-1() (2.20)

TABLE 2.1 THE ZEROS OF THE DISPERSION FUNCTION Conditions Roots Cll+aC22 < a/2 2 Real C = 0 Cll+aC22 > a/2 2 Imaginary Cll+aC22 = c/2 2 Infinite Cll+aC22-2C < u/2 2 Real C < O C l1+C22-2C > a/2 2 Imaginary Cl1+GC22-2C = a/2 2 Infinite Cll+aC22-2C < a/2 2 Real C > O C11+CC22-2C > a/2 2 Imaginary C22 > 2C-T(l/a) Cl+oaC22-2C = </2 2 Infinite Cll < a/2 and C22 > 1/2 2 Real or and C11l > </2 and C22 < 1/2 2 Imaginary C -llaC22-2C < a/2 4 Real C > 0 Call < a/2 and C22 < 1/2 Cl1+oC22-2C > a/2 2 Real and 2 Imaginary C22 < 2CT(L/o) C1+laC22-2C = </2 2 Real and 2 Infinite Cil+aC22-2C < a/2 - 4 Imaginary C13, > o/2 and C22 > 1/2 Cl1+~C22-2C > a/2 2 Real and 2 Imaginary C11+C22-2C = a/2 2 Imaginary and 2 Infinite

13 where t(~i) = C22-2TiCT(l/Oai)) We note that ~(z) is analytic in the complex plane cut along the real axis from -1 to 1. Introducing the notation i (9) = 1, je~ Region i = 0, otherwise, we easily find +42 ()T ( 1/CA )G2 (p) -T242Col ( )+_i1 [ C22+C (4) -2C4 ( )G (4) -2C4T (4 )G) (4 ) -2C4T (l/94)92 (4 ) ] ~ (2,21) Here Q~ represent the boundary values of Q(z) as the'branch cut is approached from (abeow) That is, (4) = lim ~(ktii )o ~+0 It is convenient to use the following linear combination of Region 1 eigenvectors: -F( ) - -_ F( rIpJ. (2.22) ExplicCitly, Explicitly, C1 C 12 ('_-p ) Ci2CI1 The following linear combination is a convenient eigenvector valid over the full range, C[ 1,1]:

14 i,(t,4) = 2Cj2rjT(cj)Gj(r)F1 (ri) + [1l2CllTIT(a)]GL(rI>F2 (fl") + F(2)(TI,t)G2(n92). (2o24) Again the explicit form is C 12TP/YI -4\ (- Vtnghp rl + X(fl)6(TI- ) (2.25) Here we define g(? ) = C22'-2CT( a )1-( )-2Cn-T(1/an )82( ), (2,26) <(r1) = 1-2C22T-(T )-2C11,,T(( )-1,(n)-2C,-T-(1/cn )G2( )+4c2-T ()T,(:n )(T) +4C12T (1)T(l/a l )G2(] ). (2.27) In comparing Eq. (2.27) with Eq. (2.21), we note that we can write Eqo (2,27) in the compact form X(%) = Q +- (n) + r2 2CG (n) (2.28) Also (for later use) we have from Eq. (2O,2) f(1) 0 _ = C22+C11G1 (1)-2CT (CI)Gj(4 ) 2C+(c),-,)e ( )- 21C-,-( 2() e. (2,29) Thus the eigenvectors consist of the two continuum modes ~i(TI,1t) [Eq. (2~23)], q[- -, 1]; a2(T,ll) [Eq. (2.25)], TIE[-1,1]; and the discrete modes Fi(Li) [Eq. (2.20)].

III. TWO-GROUP FULL-RANGE COMPLETENESS8 The full-range theorem is stated as follows: 1 Any function *(t) of class G can be expanded using the eigenvectors _1(r4)2, _2(fl), and Fi+(4). We follow the standard procedure of attempting an expansion in terms of the continuum modes only; namely, 4(=) = cal(j()/1(r,4)dr + J o92(l)42(7jt )df (m1) We substitute the explicit form of 1(nI,.) and &2(%I,2) as given by Eqs. (2.23) and (2.25) into Eq. (3.1), and perform the integrations where possible to arrive at = -(1/') + C12P f1 rl2(Tl)dTr (3.2).(c() =-C 11(d - - 1 1( t2 (1)? /~_/af'qoa_ (nq)d _ oaj(~L)e:,. (0 [c1:1-2c7(~)] C11C12 -1/x } [-L C1 2Cl + p l I)g(d + -2()X(4) (303) We solve Eq. (3.2) for ajl(,/a) and then make the change of variable, 4/r = rt to obtain l(r) a= &C11v1(cr) - C11C112 1 2()d (4 The insertion of Eqo (3.4) for OL(r ) into Eq. (3o3) yields 15

24)+Ca /a P1f r (acr)dTI + ga V(ac)i_') [Cl-2C4T(4)] c ffl/a Pd ln'Ca2(1 ) - Gj(1)(2C4T([)-C21) () -1/cr -1 T + rf Prc2Qrj)g(r )dr + Q24 )0X(4) * (. 5) -1 -P The integration over ~ in the first term on the right-hand side of Eq. (3.5) can now be performed by using the Poincar4-Bertrand formuial4 and the partial fraction decomposition n 1 y u+'+ ) (3.6) The result is C fl/& Pid 1 pr2(r )d< = _-Cx21aa2()Gj_(4)-_2C p2(') -1/ cr- -l r<-rl -1 G- G2Qq)T( l/Crr) o(3y7) We define,) = +2( ) CCa p jl/I IfT1(cm)dq + a (,(l)a)Wm(tC)[CCl-2C 2C(T)] 2 C12 -1/acr C12 (3.8) By using Eq. (3.7) and Eqo (3.8) and recalling Eqs. (2o28) and (2o29), we write'( ) Q+(p)-2-(Il) pi /a C2()dri [2+(l)+2-(p)]cl2(~) 2iC1 -1 -, 2

17 As usual we introduce an auxiliary function N(z) = 2(ri)dL (3.10) with boundary values given by the Plemelj formulae14 N-() = 1 P a2()d X ) (.la) 2Ti -1 A-> 2 Therefore, Q2(t) = - (5 b) We write Eq. (3,9) in terms of these boundary values to obtain the factored form This has the immediate solution N(z) = _ IdL J1'()(33) 2Tri Q(z) l kZ From Eqo (3.10) we note that N(z) must be analytic in the complex plane cut from -1 to 1. In order that N(z) have this property for z = tMi we require f tY~g~d = 0 (3.14) ~ 1 1+ ~i where we recall the ~+i are the discrete eigenvalues. Our original expansion of 4r(i) was only in terms of the continuum modes. We have the discrete modes available to satisfy EqI (3 l4) Thus we make the replacement of i(l) by

18 - Ai+ Fi+(L) - E Ai- FLi- ~ (_ 15) We recall that _(t) contains components *1(4t) and ~22(~)' and *'(4) is a function of these components as given by Eq. (3.8). Thus V'(~) in Eq. (3.14) is adjusted according to (3.15) and with two pairs of roots we would have four equations with four unknowns from which the discrete coefficients, Ai~, could be found. With the discrete coefficient known we would then determine c2(p-) from Eq. (3.1lb) where the boundary value of N(z) are determined from Eqo (3.13). Also in Eq. (3.13) we must modify )'(i) according to (3o15). With e2(4) known, al(4) is calculated from Eq. (3.4). A better procedure for the full-range case is to prove an orthogonality theorem as Siewert and Shieh8 have done. This procedure makes possible the calculation of various normalization integrals from which the expansion coefficient can be derived directly for any particular full-range problem by taking scalar products.

IV. TWO-GROUP HALF-RANGE CASE From the completeness proof of Zelazny and Kuszell, we know that the normal modes of the transport equation provide us with sufficient eigenfunctions to solve half-space problems. By a procedure which is in many respects similar to the completeness theorem proof described in the previous chapter, we develop equations for the expansion coefficients in a form convenient for numerical calculations. The expansion for the half-range case, which we assume is valid from the work of Zelazny and Kuszell, is V(jL) = (-) fl// (n \ + 2()2(,)dT (4.1) *2 41) 0, We note that we are using the continuum modes only, the discrete mode will be introduced later. Using the explicit form of 1(q,,kI) and ~2(,n4) [Eqs. (2.3) and (2.25)] in Eq. (4.1) and performing the integration over the delta functions yields ( )( = / m + C12P Isl 2 (-)dn (4.2a) C11 o 0 r -~I and *2(4) =- C f (1/ r Ca ()dd -_ 1 [Cjj-2C4T(4)]aj([)() CllC12 o rl t CllC12 + P f n2(n)g(n)dn + x(4)Q2(4) (4.2b o rl-u. Making the change of variable p + pa in Eq. (4.2a) and solving for a1(p), we obtain 19

20 ( 4) = aC- 11a-() CllC12 + fJ a2(TI)dq (4e3) 10 0 - The insertion of Eq. (4.3) into Eq. (4.2b) gives + P 1 t1*(o()dCT + [C-2CiT(t)]( Cl2 o'q-l C12 C P f l/ TdTr P l 2(T)dO + [cj-2c?(C)]G1(4) P Jf n2(Tj)d o q'- o. ~ q o'o +p fi n2(n)g(~)dr + X(d)a2(c ). (4.4) o rl -[ The Poincare-Bertrand formula14 and the partial fraction decomposition n-' _ 1 f -4 -+ (4.5) (~'1 - ) ( -n') TI ) TI -[I l- 4 - -' Y enable us to perform an integration over dr' in the first term on the righthand side of Eq. (4.4)o The result is lg r dq1P 1 nla2(Ti)dq - p1 __a2(__)d__ Gl(,)C,,n( n o,-i ao + G2(Q)Cp4n 1 - )- I(r) el)Clqn - 1) (4.6) - G2(r1)Car n -1 -Ce2L2(4.)e1(4,) We can write Eqo (2.26) as g(r ) = C22-Cr1n ~n + a) +Cwe(r)Q.n - 1) +C2() un (1 - 1 (4,7) where we have used Eq. (2.7c)o The insertion of Eqs. (406) and (4~7) into Eq. (404) yields, after some cancellation and rearrangement, the compact form

21'(1+ C J 2(-)k~ r,4)di ()( p f1 nr2(f)dr + n+(2)+n-(G) o 0- 2Ji o 2 (4.8) Here we have defined 2(1) Ca P 1 l _(aq)d + l(t)Gl(Q)[Cll-2cT(44)] C2(1) + C12- + (4~9) and k(n,~) = nin + I- An 1 + 1 (4 10) We note that the second term on the left-hand side of Eq. (4.8) is a nondegenerate Fredholm term. If this term were absent, Eq. (4.8) could be solved directly by standard procedures.l If it were degenerate, the method of Shure and Natelson15 could be used. However, neither of these situations obtains, and no closed form solution of Eq. (4.8) has been found. Thus, we describe in Chapter VI an iterative procedure similiar to that used'by Mitsisl'l6 to solve the critical slab problem and by McCormich and Mendelson17 in treating the slab albedo problem. We define "(~) =-'() + C /'f nc(n)k(',~)n (4.11) Next we introduce a function N(z) defined by N(z) = 1 /L q2()da (4.12) 2Tti 0 Z Following Case and Zweifell if Q2(F) of class G exists then N(z) has the following properties: 1. N(z) is analytic in the complex plane cut along the real axis from 0 to 1;

22 2. N(z) - as zl + o; 3. N1() P f p 2(n)dn + 1 A2( ). (413) 2Tci 0 r1-1 2 We note from Property 3 that Q2( 4) = N+ (t)-N ([) (4.14) Inserting Eqso (4.11) and (4.13) into Eqo (4.8), we find (after some rearrangement) that At"(p) = -N+(p)QC+(p)) Ni(C)Q(YI). (4015) We recall that Q(z) is analytic in the complex plane cut from -1 to 1. We note from Property 1 that N(z) is analytic in the complex plane cut from 0 to 1. This requires (as in one-speed theoryl) the introduction of a function X(z) such that x+() = i+ -) > 0 (4.16) The X(z) function satisfying the necessary restrictions1 for the half-range and for one pair of discrete roots only is X(z) = exp 1 l ArgQ+(~)d~ (417) 1-z IT o [ Z The ratio condition, Eq. (4o16), is inserted into Eq. (4015) to yield Y(f),lr"() = N I(,)X+() - WN(,)X-( ) (4.1'8)

23 where 7+() (4.19) Assuming that the left-hand side of Eq. (4.18) is known, we can write the solution to Eq. (4.18) as N(z) = 1 (~)t(1)d. (4.20) 2tiX(z) 0o -z The Plemelj formulae give the boundary values of N(z) from Eq. (4o20). We insert these boundary values in Eq. (4.14) to obtain a2(i) = 1 1 1 > P p (' -)(IG)d4 + 1 + 1 2~i0 2 -+ () x-( oa yb ( )tr (11) It(4o21) (However, this is not a solution because we note from Eq. (4.11) that the unknown a2(p) is still contained in r"(()o) We next define Q1(1) = - 1 / 14, 1 ) = ()- ) (422a) 2~1ri x+ (4) x (k) 2~it ( )4+ (4)2- 40 and = - -- + x(%) = 2 +4 o) (4.222b) The last form of Eqs. (4.22) is derived from using Eqso (4.16) and (4.19), By defining the singular integral operator 0(%) m (%) P f'(r)d- +'p(2) o (4.2) 0 L4~J

24 We can write Eq. (4.21) in the compact form Q2 (k1) = 0() I"(Q(). (4,24) Property 2 requires that N(z) - as Izl + co, We note from Eq. (4.17) z that X(z) - 1 as z + oo. Thus from Eq. (4.20), we require z o f(l) L (l) + C f dif 2()k (i[)j d = O (4.25) where we have used Eq. (4,11) for r"([i). For the case of one pair of discrete roots we have the discrete eigenfunction available to satisfy Eqo (4.25)o We make the replacement of *(C1) in Eq. (4.25) by t(l) - A+Fl+([1) o (4.26) We recall that $'(U) is defined by Eq. (4.9) and is a functional of the components of ir(lL), i.e., $1([i) and 2()~. We define ~+(4) as the corresponding functional of the components of F1+()o. Thus the replacement given by the expression (4.26) is equivalent to replacing $!(~) in Eqo (4.25) by -At () 2A 4+([) o (4.27) With this replacement made in Eq. (4.25), we then solve for A+ to yield ( )(il) di + C ()d( )k( Y I)dn () ( )a4+ 1)d A+ = o o o n-u (4.28) Likewise with Eq. (4027) inserted into Eq. (4,11) and subsequently into Eqo (4.24), we obtain

~2QL) o(~)[t?(~>........ + a oQL) f' - ) 0 -.,, -- -o n A(4.29) Finally the replacement of jz(Cp) by $m(ta)-A+Fl+, 1(t) in Eq. (4.3) yields al(Qi) = aCll [t1 (ai)-A+Fl+,l( arp -a CllCa12 Pfl e2(n )d. 0 71-1 (4.30) An iteration procedure for solving Eqs. (4.28) and (4.29) for A+ and a2(1J) is discussed in Chapter VI. Providing this procedure converges, we can then insert A+ and 0z2(I) into Eq. (4.30) and solve for acl() competely determining all expansion coefficients. A number of important simplifications are now made for the terms containing 4+(Q) in Eqs. (4.28) and (4.29). First, from Eq. (4.9) we can write the explicit form +(4) = FL 2() + P f1'+C F+ (at) G(t)[C11-2C4T([)] C12 0 r1 l+,12 (4.31) Next we insert the discrete modes as given by Eq. (2.20) into Eqo (4.31) and perform the integration to obtain G+I(=) = r " LC22+ClleGl()+C4G(1)Qn - 1) +CG2 ()2n (1 - - 2CGQ(4)liT(p)-Crjjn + ). (4.32) It is easy to show that f(p) from Eq. (2.29) can be written as f() = C22+Cll.l(p) -2C()T(p)- Cl(p)~n (1 + "0 +C[Gl(p)n (1 -41'") - Ce2 (p)n (1 + 1-) +C2 (It.)Qn (1 - (4.3)

26 where we have used Eq. (2.7c). We insert Eq. (4.33) into Eq. (4.32) [recalling Eq. (4.10)] to write the compact form + ([) = ) [f(i)-Ck(l1,4)]. (4,34) The substitution of Eq. (4.34) into the integrand of the denominator of Eqo (4.28) yields f l= - riX(7i) - C f 7([)lk(TII)dI (435) o 0 rl -4 where we have used the X-function identity given by Eq. (A.1) to write X(l) = f1 y(I)f(4)dl (4.36) 0 t-Tl Next we consider the term 0 (p)~+(p) in Eqo (4.29), With o (p) given by Eq. (4.23) and Q+(t) by Eq. (4.34), we write =41) 1 y(t')r1 f(p)dP + f) c o ( _. J (4.37) The partial fraction decomposition = _. _ _ 1, 1 1,) inserted into the first term on right-hand side of Eq. (4~37) yields ~.(1)P 1 (')lf(' )d = 1(P) pl p l7(,(')f(dP - l')d 0 (.-)(..t) L1 0 o ~[ (4.38) The last term in brackets in Eq. (4.38) is -X(rl) Esee Eq. (4o362 and the

27 first term can be written in terms of the boundary values of X(z) see Eq. (A.1)]. Explicitly, Pf1 y(,')f(~,)dp = - xI+(-)+X[ -((,)) o 2 21 L439)1 The last form in Eq. (4.39) was derived by using Eqs. (4.16) and (4,19). We insert Eq. (4,39) into Eqo (4.38) to obtain 1(111 P) 1 P1 = Q () 2 LX(r)- () ( (. (p) (4.4o) We recall from Eq. (2.29) that f(,W) = 2*(i)- -(4.41) We insert Eqs. (4.22), (4.40) and (4.41) into Eqo (4.37) to obtain (after some cancellation) O (,{)+(~) = ~1()Tl)iX(1l) C O (4) k42) Now Eqso (4.35) and (4.42) can be substituted into Eqso (4.28) and (4,29) respectively to obtain a somewhat simpler form A+ lr(4).'(>), +Cl y(ll)djt rCi (4,43a) y- 7(It)nk(Tl,4)dj and o~(~) = o (1)i' (1) *+ C o (2) Jl.o2(r)k(.-,))d, -A+ i(I' )rlX( ) C0(-) L (4 z43-bl ) r"I \14 nS4

28 We recall that the kernel k(q,t) is defined by Eq. (4.10). The operator 0(p) is given by Eq. (4.23) and ~1(t) by (4.22a). Equations (4.43) are the set which will be treated numerically, as simultaneous equations for the expansion coefficients A+ and o2(i). After A+ and c2(4 ) are obtained, a1(t) is computed from a1 AC(ll (c)12 - CllC12P f.da2 (4.44) This completes the reduction of the general two-group expansion, In the next chapter, we shall obtain some specific forms for *'(i) in Eqs. (4~43) and for l1(ai) in Eq. (4.44). In fact with these forms for *'(0) and appropriate X-function identities, the integrals in Eq. (4.43) can be partially performed.

V. APPLICATIONS OF TWO-GROUP SINGULAR EIGENFUNCTION EXPANSIONS A. INTRODUCTION We have selected for detailed study the two-group Milne and constant source problems. These problems are representative of the type of two-group half-space problems that can be solved by our method. Two half-space problems appear also to be soluable. In fact one probably can solve any two group transport problem that is also soluable in the one-group case, although in the two-group case less analytical and more numerical analysis is required. Even in one-group theory the X-function must be generated by the computer or obtained from tables. Also in order to obtain the angular flux, total flux, and current in the medium, we must calculate principal value integrals on the computer. The point is that the computer is an essential tool to obtain complete results even in the one-group case. We shall find some helpful simplifications similar to those obtained in Chapter IV before beginning the numerical work. B. MILNE PROBLEM* We define the two-group Milne problem for one pair of discrete roots in a manner similiar to the one-group case.1 The solution must satisfy the conditions (o;C1) = O, i > 0 (C5la) *This definition is correct only for two discrete eigenvalues, 29

30 and z/,Al _(Z,) - FT ()e 0 (5 b) z+oo The solution which obeys Eq. (5.lb) is expanded in the two-group normal modes of the transport equation as A_ (z ) z/B1 + A+F+ (Q)e z/r + a _z/(,)e n 0 O - + f 1 2(2) (Tp leZ /lQ drj. (5.2) We use the boundary condition given by Eqo (5.1a) (normalize by setting AN = 1) to obtain - F1(>) = A+F1(+() + f l/uW (\,( )dn + fl a2(T )(P,4)d (553) O -- 0 Therefore the arbitrary expansion function, denoted by jn(k), becomes for the Milne problem [see Eq. (2o20)}, 4t () = - (5.4) }2 (i ~lt(1 )/A l+ The appropriate rn(p) for the Milne problem (see Eq. (4.9)) is closely related to the 4+(v) as defined by Eq. (4o31). In fact ~n(~) /= -+(B,-nl) (5 5) This means that r(~t) =- -I [f(2)-Ck(-lrp)] (5.6) T'2+p

Thus we have from Eqs. (4.35) and (4.42) the result that f1 y(it)>(n)d41 = -nX(-ri) + C f 7([ )yk(-T,_ )dT (7) 0 o 0 and 0(~)%*( =) - (I)nj x(-n) + co ((5.) k(] 8) We insert Eqs. (5.7) and (5.8) into Eqs. (4.43a) and (4.43b) to obtain -n1X(-ri)+C f1 y(4')jk(-nq )d4 + C f1y(4)dt( 11 nf2(n )k(n,+)dn A+ = (5,9a) 0~llX(~) ~(C | y(t),k(j,)-)d4 ~+ C0, + C 01(1).A OK2(Q>) = + C 0(v) L ( + " o(J)l + ( a2(\,)k(xnT4t TI +~1 Tt1+I 0 A+C1(4il)lx(l C O(,) Tlk( i (5i9'b) Sinc e We have from Eq. (4 30) [considering Eq. (2.20)] that ao(i) = -CLC12TI (1 + - CJJC12P J 2() di (59c) \jj1 +I 1- 0 1 Equations (59) are the final reduced equations for the Milne pro'blem. We note that Eqs, (5.9a) and (5.9b) are two coupled equations with unknowns ca2(kU) and A+, These equations are solved by an iteration process described

32 in Chapter VI. The operator 0([t) is a singular integral operator which requires special treatment otherwise;the solution is straightforward. Finally with A+ and cx2(4) known, we obtain jl(Jt) from Eq. (5.9c). The limiting case of C = 0 is of some interest. We note that the discrete mode expansion coefficient becomes (see Eq. (5.9a)) A+f = x(-_1) (5.10) X(n11) We recall that this is identical to the one-group Milne problem result.1 Thi~ is apparently different from Siewert and Zweifel's result but we recall that they considered the case for C = 0 and the discrete eigenvalues were infinite. When the discrete eigenvalues are written in the form of Eq. (2.20) we are implicitly assuming that the eigenvalues are finite. Also from Eq. (5.9b), we have a2 (41.) = _- 2 Ii)nX (-m) (5.11) In neutron transport theory, the physical problems for which C = 0 are quite specialized and we shall not consider this case further. We can prove from the form of the operator O(j) and Eq. (4.43b) that O2(Qi) + 0 for i = - and 1. (5.12) This result is important because otherwise the angular flux would be singular at these values of -. This follows from the fact that for any particular half-space problem the solution is written in terms of the expansion coefficients and the discrete and continuum mode eigenvectors with these eigen

vectors containing terms with T(ap), T(1/o4) and T(J). We note that T(go) and T(l/a4) are singular at z = 1/c and T(W) is singular for Ct = 1, To prove (5.12), we note from Eq. (4.22) that ~1(~) I 1 + 0 (5.13a) Xl(4i).1 + o (5 l3b) -+1/*J+ T(1/CT) 2X(4)' 1 2+0 (5.13c) 1- 2T()] Similar results follow for ~2(4). The operator O(4) as given'by Eq. (4,23) contains in each term either ~1(i) or ~2(4) and we note from Eqso (4~43) that these functions always appear multiplying principal value integrals or functions which are bounded at these values of 4. At least this statement is true for A = 1/a and we note from (5.13c) that things are still all right at ~ = 1 because of the squared function in the denominator. With the expansion coefficients known from the computer solution of Eqs. (5.9), the two-group angular fluxes are calculated from Eqo (5.2). The eigenvectors given'by Eqs. (2.20), (2~23) and (2o25) are substituted in Eqo (5.2) and the integration performed where possible to derive xlV(zll) =C12leZ/Th11 C12rnlA+e-z/ l+ 1 zl/4[G1(,)+G2(,)] ~'q z+4 c14.i C1o1 + C12P fl y (eZ/)efldy 0 Corn-_i (5 14a) where

34 and *2(z,) - = _lt___,_ez qlt(i)ez )A+e-Z/1 C fl/ qal(f)e-Z/dn TIl+t 1l-1 CellC12 0 - 4 l l() [ Cll_-2c T (u) ] ()e /[ + p f ng(T )q2(T )z/ Cl2Cll O + %(4)()G2~)(4)+ 1,2()] - 1 < 4 < 1 o (5.14b) For the angular distribution to the left we note that the integrals are not singular. A good criterion for the accuracy of the numerical methods is to see how closely condition (5.la) is satisfied. For convenience in computer calculation, the replacement in Eq. (5.14a) of p by 4a is made to yield (za~) - Cl2ne +, Ci21'A+e + 1 z/ 1 = + +.(~)~ () ~1 ao(nl+r ) 1(+- ) Clic + p f1 r(n)e-Z/T dn, - 1/ < < l/a (5.15) 0 nThe total flux and current for each group are derived by appropriate integrals over dp and 4d4., respectively. We shall use superscripts on the p's to indicate group number. Thus, O1 (Z) f.1 i(zp4)d = 2C12lT(-) [e Z/nl+A+e- Z/] 11 + 1 k-z/ p' 1 -z/n x J acl(4')e d' I - C12 f fqg2(T)e / T(g()-C22)dn (5o16a) o o

35 f~l)() = | V 2a (z,,u)djut = 2Ci2reeZ/rl1[l- (Z/+2 le (1) (z) = 111 0:+ 2C2 O1e 2Cl2fldnflo2 (q )e-Z/T1_nC12 Jl 2/a2 Z/)e (g (r ) -C22 )da (5.16b) o o ~p -(2)(z) = J-2(t) 21it(ql)T-) [e /l+A+e Z1] 2C f/Olcl(T)T ()ez//ld - Cl 2CT12 (l)]e/ (1 1z/1 -z/ 1 CllC12 ~ 11 12 1 + 2f gg(rq )2 ( )T ( )e d/di + I X( l)c2(l )eZ d/d (5 16c) o o p(2)(z),)dk = 2qlt(ql)eZ/l [ 1-nlT U )]+2qlt(ri)eZ/ ~ ~x A+~~[TIT 1 -i- C 11/ cl(:l)e/f [lT(r)-1] dq 1 _ 2 J l ( 1 1ca2 (q) [C:1-2CT () ]e fZd+2 1dg(')2T (rl )1[eT (d \)-1 ]+ CllC12 o x p-z, +1 qk(T )2 (l )e 7/ll d (5o16d) 0 The extrapolation distance (the distance to the left of the interface at which the asymptotic component of the density vanishes) is given by either Eq. (5.16a) or Eq. (5o16c). Thus we wish to determine zo such that ez / + A+e l= 0 (5o17) The solution for zo from this equation gives the same result as in the onespeed theory; namely, o~ I 2n (- A+) (5 18)

36 Again we emphasize that for one pair of discrete roots the extrapolation distances are equal for each group. For four discrete roots, we would calculate a different extrapolation distance for each group, Baran9 has studied this problem in detail, C. CONSTANT ISOTROPIC SOURCE PROBLEM Assume an isotropic source S02 in group two (the associated problem with a source in group one is virtually identical). We now seek solutions to the inhomogeneous transport equation which vanish at +oo subject to the condition 4(o,t) = 0 for t > 0. (5.19) The solution consists of a particular integral, 4p(z4,), plus a solution of the homogeneous equations, x (z,a). The latter consists only of those modes which vanish at +oo, i.eo. _ (z,=) A+e /lF (t) + f al/ ()L f l2(q)(l)e Z/ /dn C 0 0(2o) The particular integral can be found from the two-group transport equations (Eq. 2.2) in the limit as zroo. In this case, V(z,4) approaches a constant value denoted for each group by 41S and ~2So With t(z,i) constant, Eq. (2.2) reduces to a pair of simultaneous equations which are easily solved to yield lf = 2SC12 = S1 (5 2La ) /-,,-,,-, /_,-~-~ \ I,,,

37 and s (a-2C j) S(cx-2C = 1-S2 (5.21b) ~2S (l-2C22) (&-2C1 ) -4C 21C12 where S = S20/r2 We normalize all constant source problems to a unit source neutron in group two, ie., S20 = 1. The complete solution is written as Si1-1 =ZC _ ~(f2~d-z/fl +F ((p) + ) 1(zt) = ()+ A+e A/ i(41) + / -()3.(,-)ez/d+f2 x e -Z/ dq. (5522) By setting z = 0 in Eq. (5.22) and applying Eq. (5.19), we obtain pefrmn t inegat o-,)n th s e(o)n(dit 4)ddN A+ i-+ (4i) (5.23) We repeat here for convenience Eq. (4.8) of Chapter IV, ~,(~) 2(~)+C: f-/~ P~(:n)d + ~ 81(a~)ol(~)[CC.~-2C~T(~) ] o C12 o ~ _ C12 (4 8) And we note from Eq. (5.23) that l(~it) = = -. ( p 2 (T() s (5,24) 982 ( 41) S2 By inserting Eqs, (5024) into the appropriate terms of Eqo (4.8) and then performing the integration on the second term on the right-hand side, we find

38 r =S( -) = -S21 + C42(Q)in (1 - 1) S saC 12 LC cl / -2CtT (Q)1 ()+C11G1([I)j (5.25 The term in brackets in the second term in the right-hand side of Eq. (5.25) has terms similar to Eq. (4.33). In fact, we can write S- C12 -a + f(Q)+CI.n (1 + 1 ) - C22 (526) We recall from Chapter IV that we need 0()Qt(4) (5.27a) and I 7Q(2) ()dt (5.27b) O S in order to calculate the expansion coefficients A+ and a2(4). By the same method that was used in Chapter IV to simplify 0([a)+(pt), we can easily prove o(4)f( ) = o0 (5o28) Also from Eq. (A.1) in Appendix A, we have that 1 lim zX(z) = - I y(4)f(I)d4. (5.29) z-oo o But from Eq. (4.17), we see that X(z) - 1 (53o0) ZO z-+oo

39 Therefore from Eq. (5.29), we obtain o f1 y()f%,)d~ = 1. (5.51) Finally this gives the result that 0() c 1 1 Lw () ~) C iLQn(l + r 7i)j (5.32a) and J I~~tsza= - _ + J (~) sw_ St n 1 + (d (5.32b) o % C12 o C12 al where w = -S2 - I + aS1C22 (5.32c) C12 Cl2 S1 and S2 are defined by Eqs. (5o21). We insert Eqs. (5.32) into Eqs. (4.43) to obtain +cr1 1 F crSC 1 )l (1)k(,)d C- 12 o C12 l J +c ()d)( A+ + (5t 5 33a) -rllX(rll) - c 7(k(Tn1_,)dt o and = oQ() LW - c 1n 1 + ) + C 0(4) J 2(T)k(j,i)dD C rl- rllThe final equation for cz1(Q) is obtained from Eq. (4~44) where:L(aC ) = -s o

40 Explicitly, al() = -aCllS1 - CllClTIIA+ - CllC12 p j C()d (.3c) Equations (5533) are the final reduced equations for the constant source problem. We refer the reader to the paragraph following Eq, (5.9c) for comments concerning the solution to this set of equations. The angular fluxes for the constant source problem are given by Eqs. (5.14) where we replace the first term in the right-hand side by S1 and S2, respectively. The neutron current in each group is given by Eqs. (5.16b) and (5.16d) if we delete the first term on the right-hand side. For the total flux we replace the first term (term with positive exponential) on the righthand side of Eqs. (5.16a) and (5.16c) by 2S1 and 2S2, respectively.

VI. NUMERICAL METHODS AND COMPUTER SOLUTION TO HALF-SPACE PROBLEMS A. PROGRAM FOR FINDING TUMBER, TYPE, AND VALUE OF THE DISCRETE EIGENVALUES In solving Eqs. (5.9) and (5.33) for the expansion coefficients, the function y(t) must be calculated. But (as shown in Appendix A) the roots of the dispersion equation must first be found. A computer program was written which first finds the number and type of eigenvalues from Table 201 and then determines the magnitude of these eigenvalueso The method that was found quite effective was the false position method. For the problems studied in this thesis we have only one pair of real discrete roots. We concentrate on finding only the positive eigenvalue rj because we know from Eq. (2.18) that the other eigenvalue is -ru. The false position method require two starting values which we denote by j(1) and ). The subscript here means that these are two starting values while the superscript is the iteration index. To obtain good starting values we calculate Q(r) for T = 1.000001. We then increase T by.1 until Q(q) changes sign. We select T(1) and T(1) such that (1)) is opposite in sign to (W(1)). From these four values a value of,r3 is calculated according to (1) ()n () l)'13 is the intersection on the TI axis of a straight line which joins the points (~l)Q(l)))> and (T )( ))) For the next iteration the program selects

42 an(2) from and a 2 such that the corresponding functional values for Q( ) are smallest in absolute value. For the specific problems studied in this thesis only three iterations were required to find a root r. su.ch that o(rll) < LO7 10 < lo. This program is listed in Appendix D. B. DISCUSSION OF MESH SPACING In two-group, half-space problem there are always two basic integration intervals. We see from previous chapters that we have ordinary integrals and principal value intervals for integration ranges of O < <n and - < ~ < 1. We shall see in the next section that in the evaluati.on f pri.Lncipal valu.e integrals the singu:Larity is subtracted. Thi.s procedure gi.ves rise to a Logarithm term which is s:i.ngular at the end points of integration. Therefore we use the midpo:i.nt approximation in the evaluation of all integrals..In this approximation the interval from 0 to 1/o i.s di.v:ided into P equal1 intervals and all functions are evaluated at the midpoint of each of the P.intervals. A similar procedure is used for the mesh spacing in the interval from 1/n to 1.

43 C. NUMERICAL EVALUATION OF SINGULAR INTEGRALS In the solution of Eqs. (5.9) and (5.33) for A+ and c2(t) we must evaluate singular integrals of the type p f1 F([1')did' o?11' - The function F(l') in two-group, half-space problems studied in this thesis always contain factors of y(lt) as well as logarithmic functions. For example, from Eq. (4,43) all half-space problems have singular integrals where F(Ql) is of the form where we recall that k(-q,) = [n(ln 1 + ) -n 1+ -+ ) For the Milne problem we have in addition other integrals involving Wi the c(p)k(Tlh4a)the i g 6e While the constant source problem has another integral Where F) = w aS1C In (1 + (ce- (6o4) 1C2 The details for evaluating y(Qt) are given in Appendix A. The presence of logarithm factors and 7(i) in the function F( ) requires special care in order that an inordinate amount of computer time is not taken.

44 By usual proceduresl the singularity is subtracted to yield lp F(,C )dL' = l [F(g, )-F(4)]d~' + F(~)~n (6o9) o F~'4 O F n'-5 The functions F(Ol) and In(l-4)j/ are now evaluated at the midpoint of all the intervals for ~ between 0 and 1, and these values are stored in the fast access memory of the computero Let us assume that the singularity is at some midpoint value ti. Then we approximate the principal value integral in Eq. (6.5) by Di j1 F(I' )dR1' X [F( j)-F(4i)]hj + hiF, (i) + F(li)n(ipi 6) j=l 4j The prime on the summation indicates that in the sum the term j = i is not included. This particular term is given by the second term on the right-hand side of Eq. (6~6). Also hj is the width of interval j. We are using the trapezodial rule with the midpoint approximation in the evaluation of the integral term on the right-hand side of Eq. (6.5). We note that the second term on the right-hand side of Eq. (6.6) has a derivative factor. When the singular point pi is not in the interval adjacent to the end points or adjacent to the i.ntervals where an interval change is made, io.e, l = 1/0, this second term on the right-hand side of Eq. (6.6) is approximated by 11 hiF' (i) =2 [F( +1 )iF( When the singularity is at the midpoint of one of the other four remaining i.ntervals we use19

45 1 [F(42)+3F(i1)-4F(O)] i = 1 (6.7a) hiF'( i) [ 1 [4F(1/a)-3F(tp)-F(tp_i)]' i =P (6.7b) 3 1 [F(p+2 )+3F(ip+l)-F(l/a)], i = P+l (6.7c) I1 [4F(l)-3FQ(D)-F(Dl)], i = D. (6.7d) 3 We have P intervals from 0 to i/a and a total of D intervals from 0 to 1. Following to some extent Bareiss and Neumann,19 we shall derive in detail Eq. (6.7a), Let us first calculate the contribution of the integral term on the right-hand side of Eq. (6.5) to one general interval of width h about kt,' ie., h h h - < 4 < +-. 2 We have h + 2 f[F(p)-F(pdp fh/2 [F(y+4i)-F(4)]dy (6 8) h2 -h/2 Y 2 where we used the change of variable y = p'-4. We expand F(y+i) in a Taylor series about 4 to yield F(y+p.) = F() + yF'(4) + C y F"(t4) + Ho (6~9) 2 Inserting this expansion in Eq. (6.8) and integrating over dy yields h/2[F'() + 2 F"(1) + y2 F"'(p) +..o]dy =hF'(4) + h Ft'"() + O -h/2 3 (610o) Thus by keeping only the first term on the right-hand side of Eq. (6~10) we are omitting terms of 9(h3) or higher.

46 For Eq. (6.7a) we need to find a good approximation for hlF(i1) where l = 2L hi and h1 is the width of the first interval. To this end we write a Taylor series for F([+h); namely, F(Q+h) = F(p) + hF'(p) + - F"(Q) + S(h3). (6.11) 2' We insert 4 = - h1 in Eq. (6.11), and recalling from the midpoint values that = 3 h, we find F(42) = F(p1) + hlF'(l) + h- F"(t41) + Y(hh). (6.12) 2 Next we replace h by -h1/2 and 4 by hl/2 in Eq. (6.11) to arrive at F(O) = F(i1) - h1 F'(.1) + h_ F,"(1) + Y(h 3) (6[13) 2 8 2 Solving Eq. (6.13) for h- F"((41) yields 2 2 2h- F" (1) = 4F(O) - 4F(p) + 2hFI'(.l) 4'(h3) o (6.14) This result is substituted into Eq. (6.12) to eliminate the second derivative term. We then solve for hlFs (p1) to obtain Eq. (6.7a), Equations (67b), (6.7c) and (6.7d) are derived in a similar manner. The nonsingular integrals are also evaluated by the midpoint approximation. For example the integral in Eqo (5.9a) is approximated as follows: D 0 1-L.r y n4 -,r ) -,i, (p4j)q ( k(r ) p )hI (6 15 ) o I, Ji

47 D. ITERATION PROCEDURE FOR THE EXPANSION COEFFICIENTS To solve Eqs. (4.28) and (4.29) in the general case, we first approximate A+ by (o) o A+ = (6.16) f y(4))+( ()d 0 The first iterate solution for ac2(i) is given by 4 1)(1 ) = o()[,,()_(o)~+(~)].(617) (6._?) The superscript on A+ and a2(It) specifies the iteration index. Equation (6.17) is then used in Eqo (4.28) to yield A+ and then the second iterate on o2(4) is found; namely, (2) + () 2 (4) = o0()['(4)-A )+(1)] + C O)fl O (6.18) o Continuing in this manner, we have for the n'th iterative solution.n)(4) = O(i4)[ (i)-A +(! )]I + C 0(4).1 n2 ()k(,d (6.19) 0 We assume that for most physically interesting problems, e.g,, Milne and constant source, the above procedure converges in the sense that given an e > 0 we can find an integer M such that for all m > M and Be[O,] l( ) <m-l (6120) For the problems studied in this thesis, convergence was achieved in seven iterations where E = 107O An alternate method of solution by matrix inversion is discussed in Appendix C.

VII. SPHERICAL HARMONIC SOLUTIONS TO THE CONSTANT SOURCE AND MILNE PROBLEM In the present chapter some commonly used approximation methods are used to solve the MiSne and constant source problems which have been treated exactly in the earlier chapters of this thesis. We restrict our attention to the P1, P3, and DP1 cases,l and in the following chapter, compare results in order to derive a feeling for the validity of the various approximation schemes. In the P1 and P3 cases the angular flux is expanded in a series of Legs endre polynomials as 00 ~(z ~) = 2Q+1 P1(4)p_(z), (7.1).=0 with ()(Z (z) = (2 ) and r(z,.) = ( 1 ) (2)(z The superscript on the p's denotes the group numnber. We follow the standard procedurel of substituting Eq. (7.1) into Eq. (2~2), multiplying by Pk(4) and integrating over do. By including a source term in group two (consistent with Chapter V) and using the identity PP( = 2_+ [(1+l)p1 (4)+1P11() J 2~++1 +19 with the orthogonality relation 48

49 i P (4)Pkk([i)dF 2 k -1 ~ 2~+1 the two group P1 equations become' F dp (1) jz) ______ (1) k+k1 [+ (k+1) dk+zI + (pk (z) 2k+l dz dz = [2Cllp(1)(z)+2C12p(2)(z)1]ko (7.2a) (2) p(2) dpk- (z) +dp (+2)) 2k-l + (k+l) kd 1 | + P (z) 2k+l dz dz k = [22 p(1)(z )+2C22p(2 ) ( ))+2S] ko (7.2-b) where as before - S 20 02 The Pi approximation is made by truncating the expansion in EqO (7,1) i.e.; p (z) = 0 for m >. (7.3) In particular, the P1 equations are.dz ( + p)(z) = 2Crap(l)(z)+2C12p ( z), (7.4a) dz dpz-)(z) + 3,p (z) -= (7.4b) dz 1 dp (2) (z) = 2C21po() (z)+2C22p(2)(z )+2S 7 (7.4c) and

(2) dpo (Z) + 3p(2)(Z) = 0 (7.4d) dz These equations are solved by first finding the eigenvalues of the homogeneous set. Assume solutions of the form P(i)(z) = Aij e kj (7.5a) 1 and P(2)(z) = Bij e k (7.5b) The substitution of Eqs. (7.5) into Eqs. (7.4) and cancellation of the space dependence yields the following eigenvalue equations: Aljkj + (a-2C11)Aoj - 2C12Boj = 0 (7.6a) Aojkj + 3aAlj = 0 (7.6b) Bljkj + (1-2C22)Boj - 2C21Aoj = 0 (7.6c) Bojkj + 3Blj = 0. (7.6d) A nontrivial solution exists to Eqs. (7.6) if and only if the secular determinant vanishes: k c-2Cll 0 -2C12 3a k 0 0 = 0 (7?7) 0 -2C21 k 1-2C22 o o 3 k

51 The expansion of Eq. (7.7) yields k4 + Bk2 + D = 0 (7.8a) where B = -3[1-2C22+a(a-2C1)1] (7.8b) and D = 9a[a(1-2C22)-2C11+4C]. (7.8c) For the problems studied in this thesis, we obtain four real roots to Eqs. (7.8). We denote these roots by ~k1 and ~k2 where k1 < k2. The particular integral solution to Eqs. (7.4) is derived by realizing as z-oo the total flux in each group becomes constant. From Eqs. (74'b) and (7.4d), the current in each group is zero. Thus from Eqso (7.4a) and (7.4c) we have a pair of simultaneous equations whose solution is p(l)(z) - 4C2 = 2S1 (7.9a) ~ z-Doz (1-2C22) ( a-2Cll )-4C2C21 and p(2)(z) =2(gy2C)S = 2S2 ( (7-9b) z+oo (1-2C22)(a-2C1j.)-4C12C21 For the constant source problem the solutions as given by Eqso (7.5) must be finite as z-+oo. This requires that we consider only the negative real roots of Eq. (7.8a)o Also from Eqs. (7.6), Alj and Blj can be found in terms of Aoi and Bojo Thus we write solutions to Eqs. (7.4) as follows:

52 (z) = Ao2e z+ Aekz + 2S1 (7.10a) P Ao(z) + Ao4fl(k2)e-k2Z (7.10b) - -kez p(2)(z) = A02f2(kl)ekl + Ao4f2(k2)e k2 + 2S2 (7.10c) Pl2 (z) Ao2f3(kl)eklz + Ao4f3(k2)ekZZ (7.10d) The functions fn(k) are given by fl(k) = k (7.lla) 3a f2(k) -= o(i-2C11)-k2 (7.llb) 6 aC12 f3(k) = k[37(o-2C1L)-k2] (7.llc) 18yc12 The corresponding Milne problem solutions to Eqs. (7.4) are p(1)(z) = eklz + Ao2ek + Ao4e-k2Z (7.12a) 0 p(l)(z) = fl(-kl)eklZ + Aoefl(kl)e-klz + Ao4fl(k2)e-k2Z (7y12b) p(2)(z) = f2(kl)ekz + Ao2f2(kl)e-klZ + A04f2(k2)e-k2Z (7o12c) (2)(z) = f3(-kl)eklz + Ao2f3(kl)e-klz + Ao4f3(k2)e-k2Z. (7.l2d) Here we have used the same normalization as in Chapter V. We note that there are two arbitrary constants in Eqs. (7.10) and Eqs. (7.12) which are determined by the boundary conditions at z = 0. We approximate the rigorous boundary condition

55 r(o,ii) = 0 for Ct > 0 (7.13) by the lowest order Marshak boundary condition, 1 i1 4B(o,4)d = 0. (7.14) 0 - For the P? approximation the angular flux is given by Eq. (7.1) where p (z) = 0 for ~ > 1. Thus we have in this approximation that _=(z,) ) +(Z) 2 1(z) (7z15) The substitution of Eq. (7.15) into Eq. (7.14) and subsequent integration yields o(+ (o) = 0. (7.16) 2 For the constant source problem, the boundary condition given by Eq, (7.16) Lwhen related back to Eqs. (7o10)] yields Ao2g(kl) + Ao4g(k2) = -S (7o17a) and Ao2h(kl) + A04h(k2) = -S2 (7.17b) where g(k) =!+ fl(k) 2 and h(k) = 1 f2(k) + f3(k)

54 Likewise, from Eqs. (7.12) Ao2g(kl) + Ao4g(k2) = -g(-kl) (7.18a) and Ao2h(kL) + A04h(k2) = -h(-kl). (7.18b) With the coefficients determined, the complete solution in this approximation is given by Eq. (7.15) where the p's are given by Eqs. (7.10) and (7.12) for the constant source and Milne problem respectively. In the problems considered here the eigenvalues are real. One of these roots is less than unity (in absolute value) while the other is greater. Thus in two-group theory we obtain one of the continuum eigenvalues even in the P1 approximation. In the one-speed case1 the P1 approximation yields only the asymptotic solution, i.e., the eigenvalues are all greater than unity in absolute value. (At least the above statements were found to be true for the two-group problems considered in this thesis.) The two-group extrapolation distance is found from either Eq. (7.12a) or (7.12c). We set A04 = 0 and solve for the extrapolation distance Zo. This gives Zo = I n n( A1- (7.19) A2k AO We note that the extrapolation distance is the same for each group and is given by Eq. (7.19). This result can be contrasted with the usual3 procedure of using an extrapolation distance for each group given in the isotropic scattering case by.7104Xl and.7104X2 where Xl and \2 are the total mean free paths for Group 1 and Group 2, respectively.

55 In the P3 approximation, we truncate Eq. (7.1) by p,(z) = 0 for ~ > 3. (7~20) This condition along with Eqs. (7.2) and setting k = 0,1,2 and 3 yields the following set: dp (l) (z) (2) dp1( ) + (a-2Cll)p!)(z) - 2C12p (z) = 0 (7,21a) dz 2dp (z) + +(z 2dpl) (z) + dp l)(z) + 3Sp(l)(z) = 0 (7. 21b) dz dz 1 p()(z) 2dp(1) ( c) dz dz 23dp_ (z + 7(p) (z) = 21d) dz + (1-2C22)P (z) - 2C2lp l(z) = 2S (7021e) dz 2dp2)(z) dp+ (2)(z) = 0 (72Lf) 9 9 3P (Z) 0 (7o21f) dz dz 3dp (2)(z) 2dp )(z) 3dP2(z) 2+dp z) (2)(z ) = 0 (7.21g) dz dz 3dp(2)(z) (2) 2 p,. 7p3 (z) = 0 (7 2lh) dz By a procedure analogous to the previous Pa solution we obtain from Eqs. (y.2l) a secular determinant of order eight which, upon expansion, yields

Ck8 + Dk6 + Ek4 + Fk2 + G = 0. (7.22) The coefficients in this equation are given by C = 81 (7 23a) D = 990(C22+Clla) - 810(l+a2) (7.23b) E = 100la(llC11-9a)(11C22-9) + 945[1-2C22+C2(a-2C11)] - 1102 C21C12 (7.23c) F = 1050a(1-2C22)(11C11-9n) + 1050a3(a-2Cll)(llC22-9) + 210llOaC2lCl2(1+?2) (7.23d) G = (105)2a(3(a-2Cll)(1-2C22) - (210)2a3C21C12. (7.23e) In analogy with the P1 solutions, we can write immediately the solutions for both the constant source and Milne problem. The constant source problem solutions are 4 (1)(1z) = C Aojfn(kj)e iZ + 2S1no (7o24a) j=1 n = 0,1,2,3 P(2)(z) = j Aojfn+4(kj)ekz + 2S2no (7.24b) j=l And the Milne problem solutions are given by 4 p(l)(z) = X Aojfn(kj)e kz + fn(-kl)ek1z (7o25a) ji l n = 0,1,2,3 Pn= (X C ) Aojfn+4(kj)e 3 + fn+4 (-k)eklz ) (7 25b) ~on j=l

57 In Eqs. (7.24) and (7.25) we define fo(k) = 1 (7.26a) fl(k) = k(35-9k2) (7.26b) 5a(212-llk2) f2(k) = 14kofl(k) (7.26c) 35o2-9k2 f3(k) = k f2(k) ('726d) 7a f4(k) =1o)lOC2(21-llk2 (7 26e) k2(9k2( -35)+5(1-2C22) (21-11k2) f5(k) (35-9k) f4(k) (7.26f) 5(21-11k2) f6(k)4k (7.26g) fG5(k5-91k2 f- (k) 35-9k2 f7(k) = 5k f6(k) (7.26h) 7 The Milne problem is normalized as before and we identify the four possible positive roots of Eq. (7~22) as kl, k2, k3, and k4 where k. is the root of smallest magnitude. The unknown coefficients in Eqs. (7.24) or Eqs. (7.25) are now determined by two Marshak boundary conditions; namely, J Pa(4)r~~o>4l = 0,) (7.27a) O

58 f P3(i)*(o,j)d = 0. (7.27b) o - In this approximation the angular density at z = 0 is,(o,') = 2o(o~)+F2 Pi(" ) (o ) + 5 P2(~t)P2(o) + 7_ P3(2 )p3(o) * (7.28) We substitute Eq. (7.28) into Eqs. (7.27) and perform the integration over dt to derive 2p (o) + p (0 o) = o (7.29a) 2 — o ~ and PO(O) + 3 P (o) + P(o) +) = o. (729b) Thus for the constant source or the Milne problem, Eqs. (7~29) give us the four equations needed to determine the unknown coefficients. Specifically, we have for the constant source problem the following set of equations: 4 g(kj)Aoj = -Si (7-30a) j=l j=1 4 t (kj)Aoj = -2. (7.30c) j=l s2 And for the Milne problem,

59 4 g(k)Aoj = -g(-kl) (7.31a) j=l 4 X h(ks)Aoj = -h(-kl) (7l31b) j=s 4 t(kj)Aoj = -t(-kl) (7.31c) s(kj)Aoj = -s(-k1). (7.31d) j=jL Here we have g(k) = + fl(k) + 8 f2(k) (7o32a) 2 8 1 2 h(k) + fl(k) + 8 f2(k) + - f3(k) (7,32b) 45 8 5 t(k) = 1 f4(k) + f5(k) + 5 f6(k) (7.32c) 2 8 s(k.) = 1+ A f5(k) + 3 f6(k) + 2 ff(k) (7+32d) 4 5 8 5 With the coefficients, i.eo, A's, known from solving either Eqs. (7.31) or (7.32), we can obtain the total flux and current from Eqs. (7.24) or (7~25) or the angular flux from Eq. (7.1). Also the extrapolation distance is given by Eq. (7.19) with Ao2 derived, in this approximation, from Eqs. (7.31). An approximation method which permits the exact specification of the boundary condition at a free surface is the so-called double P1(DP1). This is in contrast to the ordinary PQ method in which the boundary condition can be fulfilled only approximately (for example, in the previous section we used the Marshak condition)~ The angular flux is expanded as

6o 00 = f (2i+l)P (2+-1)p~(z), 0 < i < 1 = (21+l)P(2k+1l)p-(z), -' < < 0. (7.33 ) 1=0 In the DP1 approximation, we set p(z) = p_(z) = 0 for I > 1 (7.34) The boundary condition at the vacuum-medium interface is p+(o) = P(o) = 0. (7.35) We shall see that this approximation is of the same computational difficulty as the P3 approximation, but apparently gives better accuracy. For the solution of Eq. (2.2) in this approximation we substitute the expansion given by Eq. (7.33) into Eqs. (202) and multiply successively by Po(24-l), P1(2t-l), Po(2[+l), Pl(2~1+l) and integrate over dl. The range of integration for the first two cases is from 0 to 1 while the latter two is from -1 to 0. We easily derive the following equations: d.pl)Z) + dp (z) + 2ap(l)+(z) 2Cllp(l)(z) + 2C12p(2)(z) (7.36a) dz dz 0 0 0 dp (l)+(z dp(i)+(z) 0 p+ 3 + 6ap(l)+(z) = 0 (7.36b) dz dz dl(z) dpl(z) = () 2C2p() (7+36c) dz dz 0 0

dp(l)-(,) dp~-)-(z) (1)dP~J-(z) 3 f1 + 6apl (z) = 0 (7.36d) dz dz dp)+() 2)+ () + 2po(2)+(z) 2C21po() (z) + 2C22p(2)(z) + 2S (7.36e) dz dz o p(2)+( dp + (Z) d + 3 5 p (z) + 6p(2)+(z) = 0 (7,36f) dz dz dp(2)-(z) dp(2 ) -(z) dz dz 1 p(1 Z) = 1 (4))( ) (7 (l37a)) Pf - o (7.3a) )2)(z) =) = ( +(2)(z) We could now obtain a secular determinant from Eqs., (7.36) in a manner quite similar to the P3 approximation. It has been found convenient 1 to derive an equivalent set of equations which are in form very similar to Eqs. (7.21). First we write the following identities pP (z ) =,, [P1 (2z,1-)+P0(2z4-l) 1 (7,38a) P:( 1 P ( =-l)+5Pp(21l)1 (z7 8c) 22l)=4 P(1+)3(2+)(78 pj(cl) = 1 [pl(2 +l)_ Po(24+!)1 (]8d)

62 P3(I) = 1 [P3(24-1 )+5P2(2_-i )+3Pl(2p-il)-Po(2C_-l)] (7.38e) P3(Q) = 8 [P3(24+1)-5P2(2l+1)+3Pl(24+1)+Po(24+1)]. (7 38f) We use identities given by Eqs. (7.38a) and (7.38b) to obtain the net current for each group, i.e., p (z) = f Pj(t)*j(z,4)d =- [P ( z)- (Z)P (Z)] + I [P (Z)+P (z)] 2 P o 2 F (7.39a) ( 2) (z) J1 P(4)2(z,)d = 2 [ (2)+(z)p(2)-(z)] + 1 [ (2)+(z)+p(2)-(z)] P! -1 o 22 (m ) (7.39b) Similiarly by using the appropriate identities above and considering the truncation given by Eq. (7.34), we obtain p(z ) =. [(z)- (z)l (7.40) and 3(Z) = 8 -P )+ 8 2e(z)+Pj(z)]. (7,41) 8 ( p( + 1 L p2(z) and p (z) are approximate full-range moments because we omit terms of p+ and p- for ~ > 1. From Eqs. (7.37), (7.39), (7.40), and (741), we can easily solve for the half-range moments to yield +4(z). P (z) - p (z) + I p (z) (7.42a) p-(z) P (-) (7' 42b) -o: Z P -o~z) - Pz~z) + P-3(z

63 P(z ) = ( p (z ) + p (z) (7042c) pi(z) ~ P(z) + (z) 3P (z) ) (7.42d) 1a 4 51 -- -32 By appropriate linear combinations of Eqs. (7.36),we derive an equivalent set of equations in terms of the full-range approximate moments; namely, dp(-) (z) dp (Z) + (a-2Cll)p(l)(z) = 2C12p(2)(z) (7.4y a) dp z + 2 + 3 oP (z) = 0 (7.43'b) dz dz 3 4........ (6) = O0 () 0 (7.43c) dz dz dp (z) + 6p(l)() = (743d) dz 3 dp ___ (2) d7o dp d2)( z) + (12C22)(2)(z ) 2C2p () 2S (7.43e) dz dzd dp, (2) () dp (z) d+ 2 + 3p(2)(z) = 0 (7.43g) dz dz 1 d3 2)(z) + 64 )8P (0) O (7.43h) dz dz 2 je ) (z ) + 8p (z) = 0 (7.43g) dz For example, to obtain Eqo (7.43a), Eq. (7.36a) is added to Eq. (7.36c) where the derivative term is simplified by Eq. (7.39a). The other equations are derived in a similar manner. We now solve Eqs. (7043) in a manner analogous to the P3 approximation. The characteristic equation obtained (on expansion of the secular determinant) is

64 Ck8 + Dk6 + Ek4 + Fk2 + G = 0 (7.44) where the coefficients are now given by C = -1 (7.45a) D = 24[a(a-Clj)+l-C22] (7.45b) E = 576a[C21C12+(l-C22)(C1l-a)] - 36[a3(a-2Cll) + 1-2C22] (7.45c) F = -864a[(Cll-a)(l-2C22) - o2(1-C22)(a-2Cll) ] - 1728C21C12Ca(C2+1) (7.45d) G = -1296u3(a-2Cll)(1-2C22) + 5184C21C123.~ (7.45e) The solutions to Eqs. (7.43) is exactly the same as the P3 solutions. For the constant source problem the solution is given by Eqs. (7.24) and for the Milne problem by Eqs. (7.25) where Eqs. (7.26) are replaced by fo(k) = 1 (7.46a) fl(k) k(12a2-k2) (7.46b) 12G(3oe-k2) f2(k) = 9kfl(k) (7.46c) 2(12ca2-k2) (7.46c) f3(k) = kf2(k) (7.46d) 6a f4(k) =24C2(k2-3) (7.46e) k2 (12-k2 )+12 (1-2C22) (k2-3) f5(k) = k(k2-12)f4(k) (7.46f) 12(k2-3)

65 f6(k) = 9kfs(k) (7.46g) 2 (12-k2) f7(k) = kf6,(k) (7.46g) 6 The boundary conditions in terms of the half-range moments are given'by Eq. (7.35). From Eqs. (7o42a) and (7.42c) we write the equivalent boundary conditions in terms of the full-range moments as p o(o) + 3 P (o) - p (o) = 0 (7.47a) 2 30 4-' -3 and p (o) + 2 P(o) + (0) = (7.47b) For the DP, approximation, we define g(k) = 1 + 2 fl(k) - f3(k) (7.48a) 2 4 h(k) -1 f (k) + 2 f2(k) + f3(k) (7.48b t(k) - f4(k) + (k)(k) - f7(k) (748c) s(k) = - f5(k) + - f6(k) + f7(k). (7.48d) 4 3 The equations for determination of the unknown coefficients is identical [except we use Eqs. (7048)] to either Eqs. (7030) or (7.31) depending on whether we are solving the constant source problem or the Milne problem. The total flux and current is given by Eqs. (7~24) or (7~25), respectively, the extrapolation distance by Eq. (7.19), and the angular distribution by Eq. (7?33)

with the truncation indicated by Eq. (7.34). In terms of the full-range moments the angular distribution from Eqs. (7.33) and (7.42) is (z,4) = _ pj(z)-2,(z) - 4 (Z) + 3J4P2(Z) + 3 p (z) + 6p (z)], 0 < i < 1 (7.49a) and _(zU):= 1 p (z)-2p (z)+4p (z)+4[-4P2(Z) + 2 pl(z)+6p3(z)1, -1 < 4 < O. (7.49b) We note the discontinuity in the angular distribution at p = O. In the next chapter, these solutions are evaluated numerically for some specific cases.

VIII. NUMERICAL RESULTS AND COMPARISONS We have selected four problems for detailed study. These are all light water systems. The first problem which we label Set I is for ordinary light water. Set II, III, and IV consist of boarted light water at concentrations of 1.025, 2.99, and 6.35 barns/hydrogen atom respectively. The ranges of the energy groups are chosen as Group 1: O. eV < E <.0253 eV and Group 2:.0253 eV < E <.532 eV The thermal spectra and cross section averaging procedures were performed'by using the INCITE22 code. The McMurry-Russell 3 H20 kernel was used at room temperature (293~K). The results of these cross section calculations are listed in Table 8.1. TABLE 8.1 TWO-GROUP MACROSCOPIC CROSS SECTIONS Set No. a a2 la a 2a 11 022 012 21 I 4.8822 3.2343.03166.01498 3.8180 2.8669.3524 1.0326 II 4.9270 351686.09725 o04424 3.7953 2,8005.3239 1.0345 III 5o0914 3.0707.28011.11738 3.7659 2.6828.2705 1.0454 IV 5. 3220 2.9738.58336.22326 3.6906 2.5341.2164 1.0481 Recall that these cross sections must be modified to 0 = 01/62 67

68 and Cij = ij/2G2 A computer program as described in Chapter VI is then used to find the eigenvalues of Eq. (2.17) i.e., the dispersion equation. For the four problem sets listed in Table 8.1, two real roots are obtained. In the process of solving the problems by approximate methods as detailed in Chapter VII,we also must obtain the roots of the characteristic equations i.e., Eqs. (7.8a), (7.22), and (7.44). The eigenvalues for the exact calculation and the approximate methods are given in Table 8.2 where from Eq. (2.3) and Eqs. (7.5),we see that kj = l/nj. The eigenvalues for ~ > 1 provide the asymptotic solution to the onegroup transport equation. This is also true in two-speed transport theory and we note in Table 8.2 that the P1 discrete eigenvalues are closer to the exact values for Set I and Set II than are the P3 or DP1 values. But we observe that as more absorption is added to the medium the P3 discrete eigenvalues'become almost identical to the exact eigenvalues (see Set III and Set IV in Table 8.2). With the discrete exact eigenvalues determined for each case, four Milne problem were solved by using the computer program listed in Appendix D. The behavior of the expansion coefficients for Set I and Set IV is shown in Fig. 8.1 and Fig. 8.2 respectively. For all problem sets, i.e., Milne and constant source, we always obtained the sharp resonance-like behavior in a2(Q). The "resonance" always appeared for. between 1/a and 1l In comparing Fig.

TABLE 8.2 EIGENVALUES OF THE TWO-GROUP TRANSPORT EQUATION Exact P1 P3 DPi Set No Discrete Continuum Discrete Continuum Discrete Continuum Discrete Continuum.8599.8359 I 7.1869 0 < < 1 7,2482.7435 7.2647.4835 7.2636.2786.3040.1760.8293.8054 II 4.1793 0 < < 1 4.1746.7158 4.2040.4823 4.2019.2780.2941.1703.7569.7333 III 2,5958 0 < < 1 2,5456.6494 2.5961.4788 2.5918.2765.2725.1581.6704.6466 IV 1,9361 0 < T < 1 1,8638.5683 1.9363.4730 1.9287.2740.2479.1444

70 -.004 -.003 -.002 (F) -.001 O -.001 A+= -.8309 -.3 -.2 0I.05 I t I I I 0.1.2.3.4.5.6.7.8.9 1.0 Fig. 8.1. Expansion coefficients for Milne problem, Set I.

71 -.008 -.006 al(F) o004 -.002 O\'l.002 A+ = -.4555 -24 -2.0 -1.6 a2(F) -1.2 -.8 -.4.4 I I I I I I I I I.1.2.3.4.5.6.7.8.9 1.0 Fig. 8.2. Expansion coefficients for Milne problem, Set IV.

72 8.1 with Fig. 8.2 we see that the added absorption causes the resonance to shift closer to 1/a and to become much narrower and higher (note the change of scale for (X2(i) in Fig. 8.1 and Fig. 8.2). In solving the Milne problem we obtain the extrapolation distances by using Eq. (5.18). The extrapolation distances were also caculated from the approximate methods using Eq. (7.19). The results are shown in Table 8.3. TABLE 8.3 EXTRAPOLATION DISTANCES FOR MILNE PROBLEMS Set No. Exact P1 P3 DP1 I.6658.6217.6736.6944 II.6783.6286.7111.7320 III.7121.6476 o7976.8177 IV.7613.6735.9279.9432 We see from Table 8.3 that the P1 results are consistently too low while the P3 results are two high and with increased absorption the Pj and P3 results become progressively less accurate. The DP1 extrapolation distances are about 3% higher than the P3. Angular distributions, total flux, and current are compared only for the constant source problems since the general conclusions should be valid for either problem. Also the constant source problem is of more interest for reactor applications. All distances into the medium are given in units of mean free path (m.f.p.) where we use the larger m.f.p. of the two energy groups, i.e~, l/a2.

73 In Fig. 8.3 the exact angular distribution of Group 1 at the interface is compared to that obtained by the three approximation methods. The P1 exit distribution fits the exact fairly well. As expected the P1 approximation gives both positive and negative inward components. The P3 calculation also gives a negative inward contribution which is difficult to show on the scale that we use in Fig. 8.3. We have indicated the DP1 results by points at +90~ and at ~101.5~. For G > 101.50 the points for the DP1 calculation were very close to the exact results. We note the discontinuity in the angular distribution at 9 = +90~ for the exact and DP1 calculations. In Fig. 8.4 we make a similar comparison of exit angular distribution for group 2. We observe the change of scale; otherwise the distributions are quite similar. The angular distribution at a distance of 1 m.f.p. in the medium is shown in Fig. 8.5. As can be seen from the figure the P3 and DP1 results are very close to the exact results. The DP1 results are slightly closer to the exact results than are the P3. We show in Fig. 8.6 the change in the exact angular distribution as a function of distance into the medium for group 2. We observe how the angular distribution becomes progressively more isotropic and increases in magnitude with distance into the medium. The effect of added adsorption on the exit angular distribution for group 2 is shown in Fig. 8.7. As a measure of the peaking of the exit angular distribution, we compute the ratio of the angular flux at k = -1 to that at P = O0 In the constant

74 VACUUM MEDIUM ~/. ~' *g~~ l.ed~~ A.6.8 1.0. LEGEND Exact. DFi Fig. 8.3. Constant source angular distribution for Group 1, f1(0,G), Set I.

VACUUM MEDIUM:'',"-.*- 2.0 4.0 5.0 \\ P3 *2(0,0), Set I. ~2(0, 0), Set I.

g ~, /e=coS't, 2.j0 6.0 10.0 LEGEND Exact P3 DPF Fig. 8.5. Constant source angular distribution for Group 2, *2(1,0), Set I.

77 (log,"() 5 I) \\t a Fig. 8.6. Modifications in exact angular distribution with Fig. 8.6. Modifications in exact angular distribution with distance into the medium, Set I.

VACUUM MEDIUM I IVn 1 --— 4 8=cs-l11 2.0 4.0 Fig. 8.7. Absorption effect on exact exit angular distribution for constant source problem, Group 2, Sets I and IV.

79 source problem, this ratio is 2.48 for Set I (pure water) while for Set IV (most poisoned case) it is 2.10. Thus with increased absorption the peaking in the forward direction becomes less pronounced. We observe the opposite effect for the Milne problem which is in agreement with the one-speed theory. Specifically, the ratios for the Milne problem are 2.91 and 3.77 for Set I and IV respectively. In reactor physics analysis the total flux and current are of most interest and utility. In Figs. 8.8 and 8.9 we show the exact total flux for Group 2 using Set I and Set IV, respectively. The approximate methods give a total flux quite close to the exact results, thus making a compariton of results difficult on the scale that we are using in these figures. Nevertheless we attempt to show in Fig. 8.9 the exact total flux and the approximate total flus. The DP1 calculation gives a total flux so close to the exact value that a separate line could not be drawn. We indicate in each figure with a horizontal line the asymptotic total flux. Numerical results for Set I and Set IV and Group 1 and Group 2 are given in Tables 8.4 to 8 6. Also the error is given in percent where we define [pe (approx)-p (exact)] ~100 error = p (exact)

80.0 7844 - ---- - 60.0 (2) 40.0 20.0 0 I I I I I 2 4 6 8 10 z Fig. 8.8. Exact total flux for the constant source vs. distance into the medium, Group 2, Set I.

16.107 12.0 8.0 LEGEND (2) Exact, DP,.'P 4.0 2.0 I I I I I 0 2 4 6 8 10 z Fig. 8.9. Total flux for the constant source vs. distance into the medium, Group 2, Set IV.

82 TABLE 8.4 TOTAL FLUXES FOR GROUP 1, SET I Di stanc e Exact P1 Error P3 Error DP1 Error (m.f.p. ) 0 1.4106 1,5966 +13.1 1.4630 +3.7 1.3975 -.9.5 3.4424. 3089 - 3.9 3.3768 -1.9 3.4364 -.2 1.0 5.0450 4.8649 - 3.6 4.9675 -1.5 5.0164 -.6 1.5 6.4926 6. 2946 - 3.0 6.4008 -1.4 6.4449 -.7 2.0 7.8245 7.6194 - 2.6 7.7207 -1.3 7.7613 -.8 3.0 10.2009 9.9946 - 2.0 10. 0829 -1.2 10. 1198 -.8 5.0 14.0417 13.8530 - 1.3 13.9153 -.8 13.9427 -.7 10.0 20.0287 19.8987 -.6 19.9202 -.5 19.9362 -. 5 TABLE 8.5 TOTAL FLUXES FOR GROUP 2, SET I Distance (m.f.p.) Exact P1 Error P3 Error DP1 Error 0 5.7604 6.5009 +12.9 6.2100 +7.8 5. 6877 -1.2 ~ 5 11.6354 11.2336 - 3.5 11.3838 -2.2 11.5746 -.5 1.0 16.2810 15. 6776 - 3.7 15.9721 -2.0 16.1772 -.6 1.5 20.5104 19.8403 - 3.3 20.1670 -1,7 20. 3304 -.9 2.0 24.4239 23.7335 - 2.8 24.0507 -1.5 24.1878 -1.0 3,0 31.4543 30.7663 - 2.2 31.0381 -1.4 31.1478 -1.0 5.0 42.8775 42.2395 - 1.5 42.4243 -1.1 42.5066 -.9 10,0 60,7205 60.2340 -.8 60,2984 -.7 60.3425 -.6

TABLE 8.6 TOTAL FLUXES FOR GROUP 1 SET IV Distance (m.f.p.) Exact P1 Error P3 Error DP1 Error 0 3.2683 3.6563 +11.8 3.4226 +4.6 3.2625.2 *5 8.3076 709137 - 4.7 8.2312 -.9 8.3723 +.8 1.0 11.5025 11.1189 - 3.3 11.4584 -.4 11.5392 +.3 1,5 13.8524 13.5487 - 2.2 13.8201 -.2 13.8656 +.1 2.0 15.6174 15. 3999 - 15.4 15.5924 -.2 15. 6174 0 3.0 17.9790 17.8879 -.5 17.9585 -.1 17.9702 0 5o0 20,1820 20.1908 + 0 20,1717 -.1 20.1776 0 10.0 21,2945 21.3033 + 0 21,2945 + 0 21.2959 0 We note the accuracy of the P3 calculation which is always within one or two percent of the exact value except at the interface. Even at the interface the error is always less than 8%. The error is always negative for the P1 and P3 approximations except at the boundary where the error is several percent positive. We observe that the DP1 calculation gives a significant improvement in accuracy over the P3 especially at the interface. Finally, we compare in Fig. 8.10 the current for the exact calculation with that obtained by the approximation methods, The P3 and DP1 currents were very close to the exact. Nevertheless, the DP1 calculation gave the'best results in every case studied.

16.0 (2) 1 0 Exact, P3,DP P (z)12.0 LEGEND -o-o- P 8.0 4.0 2 4 6 8 10 z Fig. 8.10. Constant source current vs. distance for Group 2, Set IV.

IX. CONCLUSIONS In the general two-group transport problem an exact analytical solution has not been found. Nevertheless by the Casell approach we have succeeded in deriving a pair of coupled equations (one is a singular integral equation) which are conveniently solved by numerical methods. Singular integral equations are often easier to solve than Fredholm equations. This is true in this thesis and may be generally valid,25 We develop in this work computer codes in the Michigan Algorithm Decoder Language (MAD) from which'anyone can (with different cross section sets) test the accuracy of conventional approximation methods. The numerical calculations show that although the P1 approximation yields fairly good results for the total flux and current, the angular distribution is not so well represented especially at the interface. The P3 approximation improves on the P1 but we continue (as in the PI) to obtain negative inward angular fluxes at the interface. The DP1 method improves significantly on the P3 angular fluxes as well as on the total flux and current, The DP, method has the important property that the boundary condition at a vacuummedium interface can be specified exactly even in the lowest order approximation. Thus the exit angular distribution at the boundary is well-represented. An interesting observation from this work is that in the constant source case the exit angular distribution becomes less pronounced in the forward direction with increased absorption in the medium~ This is in contrast to 85

86 the Milne problem where the opposite effect is noted. We recall that the P3 and DP1 approximation present the same computational difficulty. From this study we conclude that in two-group reactor physics design analysis (with reactor parameter similar to those studied here) the DP1 approximation is to be preferred over the P3 approximation. This conclusion has also been stated by Weinberg and Wigner24 in the one group case but never (to our knowledge) has it been tested against exact calculations in two-group problems. Of course, we realize that for complete core design calculations on complicated power or research reactors multigroup diffusion theory may still be the only feasible approach.

REFERENCES 1. K. M. Case and P. F. Zweifel, Linear Transport Theory, Addison-Wesley Publ. Co., Reading, Mass. (1967). 2. A review of some of these techniques is given in Ref. 1, Chapter 8; see also A. M. Weinberg and E. P. Wigner, The Physical Theory of Neutron Chain Reactors, University of Chicago Press, Chicago (1958) and R. V. Meghreblian and D. K. Holmes, Reactor Analysis, McGraw-Hill, New York (1960). 3. Reactor Physics Constants, ANL 5800, 2nd ed., Argonne National Laboratory (1965); see also J. L. Meem, Two Group Reactor Theory, Gordon and Breach (1964). 4. J. R. Beyster and J. M. Neill, "Status of Thermal Neutron Spectra," speech at I.A.E.A. Symposium on Neutron Thermalization and Reactor Spectra, Ann Arbor, Michigan (1967) (to be published). 5. M. R. Mendelson, Thesis, The University of Michigan, 1964. 6. R. Zelazny and A. Kuszell, Ann. Phys. 16, 81, 1961; also in Physics of Fast and Intermediate Reactors, I.A.E.A., Vienna, 1962. 7. C. E. Siewert and P. F. Zweifel, Ann. Phys. 36, 61, 1966; and J. Math. Phys. 7, 2092, 1966. 8. C. E. Siewert and P. S. Shieh, J. of Nucl. Energy, 21, 383 (1967). 9. R. Zelazny, Private communication, June 1967. 10. M.M.R. Williams, The Slowing Down and Thermalization of Neutrons, Interscience Publishers, New York (1966). 11. K. M. Case, Ann. Phys. 9, 1, 1960. 12. R. Zelazny and A. Kuszell, op. cit. (in the sense that the expansion coefficients are obtained). 13. See Chapter 8 of Ref. 1; also Joel H. Ferziger and P. F. Zweifel, The Theory of Neutron Slowing Down in Nuclear Reactors, The M.I.T. Press, Cambridge, Mass. (1966). 14. N. I. Muskelishvili, Singular Integral Equations, Noordhoff, Groningen, Holland (1953). 87

15. F. C. Shure and M, Natelson, Ann, Phys. 26, 274, 1964. 16. G. J. Mitsis, Nucl. Sci. Eng. 17, 55, 1963; see also G. J. Mitsis, Thesis, The University of Michigan, 1963. 17. N. J. McCormick and M. R. Mendelson, Nucl. Sci. Eng. 20, 462, 1964. 18. B. Carnahan, H. A. Luther and J. 0. Wilkes, Preliminary Edition of Applied Numerical Methods, John Wiley and Sons, Inc., New York, 1964. 19. E. H. Bareiss and C. P. Neumann, "Singular Integrals and Singular Integral Equations with a Cauchy Kernel and the Method of Symmetric Pairing," ANL 6988, Argonne National Laboratory (1965). 20. J. Yvon, J. Nucl. Energy 4, 305, 1957. 21. E. Gelbard, J. Davis, and J. Pearson, "Iterative Solutions to the P and the Double P1 Equations," WAPD-T-810, Westinghouse Atomic Power Development, 1958. 22. R. L. Curtis and R. A. Grimesey, "INCITE-A Fortran IV Program to Generate Neutron Spectra and Multigroup Constants Using Arbitrary Scattering Kernels," IN-1062, Idaho Nuclear Corporation (to be published). 235 H. L. McMurry, G. J. Russell, and R. M. Brugger, Nucl. Sci. Eng. 25, 248, 1966. 24. A. M. Weinberg and E. P. Wigner, The Physical Theory of Neutron Chain Reactors, University of Chicago Press, Chicago (1953). 25. P. F. Zweifel, Private communication, November 1967. 26. B. A. Galler, The Language of Computers, McGraw-Hill Book Co., Inc., New York (1962). 88

APPENDIX A THE EVALUATION OF y() Half-space identities quite similar to those in the one group case1 exist for the two group case. As the method of proving these identities parallels closely that for the one group case as given in Case and Zweifel, we shall omit the details in this work. The X-function identities that we use in this work are X(z) = Jf 7()f([L)d (A.1) o p-z X(Z)X(-Z) = Q(z) (A2) Q(oo) (,_-z2 X(z) = Lf(/.)d (A.3) o a(oo)( 2(_2)X(- 4)(4Z)-z We are considering only the case of one pair of discrete eigenvalues. Actually the only modification needed for the four eigenvalue cases is that we must make the replacement in identity (A.2) and (A3.) of;2z2 by (2-z2( 2-z2) where q, and BI2 are the two positive eigenvalues. Also in the above identities from Eq. (2o29), f(4) = ____(_)_ = C22 + [C11-2CCT(4)-C(a4)]01(4)-2C4T(1/G4)G2(4) 23-tiCi (A.4) and from the dispersion equation we have ~(oo) = 1 2C11_ 2c22 +. (A5) 89

90 The method developed by Shure and Natelson15 is now used to obtain a rapidly convergent integral equation which is then used to calculate the X-function and y(4). We write from Eq. (A.2), z2 1 - Q(z) Q(z)......... (A.6) X(z)X(-z) x2(o) We then define a new function K(z) such that K(z) = - X(z) (A.7) X (o) and note that K(O) = K(O) = le We solve Eq. (A.7) for X(z) and substitute in Eq. (A.6) to arrive at 1 - 2(z) I (A,8) K( Z ) (-Z) l-z2X2(o) Boundary values are now taken of Eq. (A.8) to S+(, - K-(11) = l-42X2(9 _____( _ _ (A.9) 1- 2 We observe that the function(K(z)-l)/z is analytic in the complex plane with a cut from 0 to 1 along the real axis. A direct application of Cauchy's theorem around this cut gives.(z )- 1 l1,+()_.(A.lO) z 2ti o ( (1-z) We substitute into Eq. (A.9) the boundary value of Q(z) as given in Eq. (A 4) and then insert this result. into Eq. (A.10) to obtain an integral equation for K(-z). This integral equation is

91 )= 1 _ Z [l-42X(o)]f(0 )da (A 1) K(-z) - )z fl (A. (, )) In the monoenergetic case this integral can be carried out explicitly in the first approximation,5 i.e., for K (-z) = 1, but in Eq. (A.11) above the function f(it) occurs and thus numerical procedures must be resorted to immediately. Actually Eq. (A.11) is a nonlinear singular integral equation'because from Eqo (A.4) and for C ~ 0 f() + ~oo as ~ + - a To simplify the notation, we define g(Iz) = [1-t2X2(o)] (A.12) (1 - i) (+Z)K(-4) By using Eq. (A.4) for f(Q) in Eq. (A.11) and by standard techniques of subtracting the singularity, we easily obtain (-z) = 1 - z f [gC22+ (C l-2CiT (O ))) G (4) ]di 0 4i + 2Cz f [g~l., Z) (aZ) I [T ( a4)Gj(4)+T (1/94)G2 40 ]d4 o + 2Cz g(,z): J1 [T(Ga4)e(4)+T(l/Gi4)G2(W)1d (A.L) This last integral in Eq. (A.13) can be performed to yield 1 1 2_1)1. (A.14) f [T(acm)GI(1)+T(l/a)2(,)]dt = [T(l/a) + a (Ao-4) o0 2 Because of the form of the integrands in Eq. (A.13) the integration is performed separately for O < y <_ l/a and 1/a < K < lo Eight point Gauss quadrature is now used to eva:luate the integrals and we ~iterate on s.(z) 1

92 for the first iteration. In a sample problem after five iterations the maximum difference in any of the 32 values of K was less than 106. From Eq. (Aj7) we see that X(-4) = K(-~) (A 15) 1 X(o) and if we recall that Y(~ = 4 X+ (k) (Ao16) ~+ then from Eq. (A.2), ((2 ( )+ ) X2(o) (A+ly) = o2 _4 () (A. 7) Once K(-z) is found then Eq. (A.17) gives y(p). The K-function for one-speed theory and isotropic scattering can be obtained from Eq. (A.11) by defining f( ) = 2 and X2(o) =,, 1,, where C has its usual one-speed definition. The computer code for calculating y(Ct) is an integral part of the computer code -which solves the exact Milne and constant source problems. This code is listed in Appendix Do

APPENDIX B NECESSARY MODIFICATIONS FOR TWO PAIRS OF DISCRETE ROOTS We have briefly noted in Appendix A the slight modification in the procedure to obtain y(p.) if there exists two pair of discrete, finite roots. If the modification is carried through we can easily derive for y(Q), 2 )2 + )X2(o) (B.) where (o) ) 1 (B,2) 2 2 Q00 1( rl l 2 and the integral equation for K(-z) has the form (-z)' = 1 z [l-t2X2(o)1f(ll)d4 (B53) ('\ - 1 (4+z)(-) Her 1,9 and 12 are the two positive discrete eigenvalues. The appropriate X(z) function for two pair of discrete roots is X(z): 1 exp 1 1 Ar-g 2+)di (B,4) (l-z)2 t o And in order that N(z) - - as z + co we require from Eq. (4.20) that () Lr( ) + C f d2( r)k(,: 0 (Bo ) (Bo6) o (f) Lr'(l) + C f d'2r (n)k( nI,) (B06) 93

94 We have two discrete modes available to satisfy Eqs. (B.5) and (B.6) and thus we make the replacement of'(,u) by'%(i)-Al+/1+(j)2-A2+2~(4) This gives two equations with two unknowns assuming a2(Qi) has been determined. The equations are J 7Q2YG'Q([i)d + (t - Al+ f 7(A2Y1f%2d1(1^)-2+ fl)dl + C + (o d o o 0 (x) fI.a-2(n0)k(n,) - o (Bo7) O T-, and o 0 + c l,(7)d4 j1 dfn2((n)k(r,L t) =. (B.8) o O T1-4 We solve these two equations for A1+ and A2+ to obtain Al+ f[ 7(l()X2+(4')d4 l f - I () )d2 (p d4 f d + fl )( 2+T)d4' f1 dkt,7(pi) f1 dn 2 (r)k(rj,p) o o - o o (X- ) ()d- l () (ll l ( (d (Bel) 2(+ )c f 0 - 11x) fl l. (BfLO) n n n -i i 0~~~~~~~~~~~~~~~~~

95 The only fundamental difference, other than more terms in this result, from the case treated in the text is the presence of terms of the type 0I 47(4)/n+(>)d~ ok(Bo1l) We recall that /1+(p) = 1 [f(t)-Ck(qiy4,1 (B.L2) Inserting this into the above integral we find { y()()d = J vx(k) "IL f(t)d - C f Ky(l[) -' k(,l4L)d4[(B.l1) The first integral can be written as &7(p) -( f(t)dli = -Tj1 f y( )f(t)dp + f y('Q) DL- f(-)d- d (B.l4) 0 1-1 0 0o o We recall the X-function defined by Eqo (Aol) is also valid for the two pair case, i.e, X(z) = l 7(I)f([)dA (A.L) o 0 - z We note that lim zX(z) = - I (Q)f(p)d o (B.ol) z-oo o0 But from Eq. (B.4) X(z) +- 1 (Bo6) z-+oo z thus I y(4)f()d = O0

96 And we find for the integral expression (B.1l) that 1 12 f Q(t) (l+()d = - x()-C f y(p) 1 k(rj4)d t (B. 1) o 1+ 0 T1 J7[ t1 Similar results can easily be derived for the other integrals of this type in Eqs. (Bo9) and (B.10). The basic equations for finding c2(4) would be the same as before except r'(Q) is replaced by *'(t)-A1+!+(~)-A2+~2+(t). Of course quite basic to the solution of a2(C) is the shape of the function y(C) which would certainly have a different functional form. Lastly the equation for al(t) is modified as follows: ca(Q) = aCll[tl(ca) - A1F+1+ (ca) - A2+F l)(Ca )]-CllC12 f1Pc a2( d)dn (Bol8) Here the functions F (cr1) and F(+)(a) are given by FI+ (1o) = (B.19a) 1F )(o) = C122 (B. 19b) F+(2-4) The simplification of 0(4)01+(P) and 0(e)w2+( ) would proceed in exactly the same manner as for 0(p)/+(4) in Chapter IV.

APPENDIX C AN ALTERNATE SOLUTION METHOD TO EQ. (4. 29) It is possible that there are values of C for problems of interest where the iteration technique described in Chapter VI will not converge at all or will converge very slowly. The method of solution outlined in this appendix will give a solution for these cases. For convenience we repeat Eqs. (4.28) and (4.29) of Chapter IV: J1 1(p)'(4)d+C f 7(4,d4f dTn12(T )k(, Hi) A+ 0 0 0 1 (4.28) and a2(4) = O(,)[t'(2)-A+/+(u)] + C O(2) f1 d lf2(n)k(,f)) (4.29) 0 We define Mo+ = 7(~)~+(~)d (C.1) A(o) o (C.2) + = (C.2) Mo+ We substitute Eq. (4.28) into Eq. (4.29) to arrive at + C O(1 ) fld2(2 ()k(.4) =0(l)*'(4) - A(~)O(,t)~+(4) - C ~ ~ Mo+ 97

98 Equation (C.3) is now reduced to a set of linear equations. We again use the midpoint approximation to evaluate the integrals. For example the last integral in Eq. (C.3) is D fo ana2(n)k(T,~) X q (a2(~i)k( i,.)hi (c.4) i=l i hi is the width of the ith interval. For the integral in the third term on the right-hand side of Eq. (Co3), we write D D D f X)i( i____ y(1 )ik(i4j)h When j=i in the sum over j, we use L'Hospital's rule to write [I j = n log(l + ) - i (c6) )k(sT ~) h i +i Thus we have for Eq. (C. 3) the set of equations D D D ( y(xyj)nik(qi',4j)hj + CQ) 2(i)jik(Ti4)hi (C) (x,) +..CO(1k) i (C67) li-~Ij =i =l For clarity, we write the explicit form for ) o(4k)~T()) = O(k-(k) A ~')O(+ [L ) +d +((ci).8)+

99 Equation (C.7) in matrix-vector form is A_ 2 = (C.9) where D = im C O(l h)+(l) 7(Yj)lmk(,4j)hj - CO( )H ( m k(fm,9)hm Mo~+ *1 4j =_ _,1 (C. 10) Al = 1 +C 0()d()) h 7([Sj)t k(j9jk j)hj _ co(l-) k(,9 )h MO+ j =1 - 0~(~)~' (I)-A( )0(~1)~+(~) x9,2 2~~~~~~~~~(C.12) / (o) o -.1 D D _C~2 ( 1D ) - _O ( ID) $'( )-A+ (1ID) +( 1) Equation (CO9) is now solved by standard matrix inversion. It was found that the iteration method (discussed in Chapter VI) was much faster with the same number of points. Nevertheless this matrix inversion method can always be applied for any value of C. One obvious requirement for a solution to Eqo (C.9) is det IAI p 0 o

APPENDIX D COMPUTER CODE LISTING Five computer programs are listed in this appendix. We denote these programs by Program 1, 2, 3, 4, and 5. For Program 1 we tabulate a list of principal variables and refer the reader to Table 2.l for information on the number and type of eigenvalueso The false position method for finding the roots of a transcendental equation is discussed in Ref. 18. We have inserted remark cards into the other programs to give details of the solution procedure. In many cases these cards refer to the appropriate equations in the text. We have listed the replacement cards at the end of Program 2, 3, and 4 which, when inserted into the previous program, will give a program for the solutions to the constant source problem. In Program 5 the changes for the DP1 are almost analogous to the P3, therefore, we have not listed these cards i.e,, we have listed only the solution to the constant source problem in the DP1 approximation. At the end of each program is a list of sample input data. All programming was done in the Michigan Algorithm Decoder (MAD) language 26 100

101 Program Number 1: List of Principal Variables Program Symbol Definition Dl det IC| Cl C1l (See Eq. (2,2)) C2 C12 C3 C21 C4 C22 SIG a OMEGAX. Subroutine for calculating 2(z) for z real (see Eq. (2.17)) OMEGAY. Subroutine for calculating Q(z) for z pure imaginary ITMAX Maximum number of iterations permitted for converging on the roots of Q(z).18 EPSl Convergence criterion for O(z)o EPS2 Criterion for terminating program when denominator is too small in false position methodo18 ITX Number of iterations, ko FX(l) O(xk) FX(2) n(xk) FX(3) 0(x3)

102 PROGRAM 1 R PROGRAM FOR FINDING THE REAL AND/OR IMAGINARY ROOTS OF THE R DISPERSION EQUATION FOR TWO GROUP TRANSPORT THEORY BY THE R FALSE POSITION METHOD DIMENSION X(3), FX(3), Y(3), FY(3) 1NTEGER ITMAX, ITX, ITY, I, Jp R START READ AND PRTNT DATA ITMAXP EPS1, EPS2, 1C1, C2, C3l C4, SIG D1=C1*C4'C2*C3 02=C1+SIG*C4 D3=CI+SIG*C4-2.*Dl D4=Dl*ELOG.((SIG+1.)/(SIG-1.)) WHENEVER D1.E.O. WHENEVER D2.L,SIG/2. R=1 PRINT COMMENT S TWO REAL ROOTS$ OR WHENEVER D2.G.SIG/2, R=2 PRINT COMMENT $ TWO IMAGINARY ROUTS$ OR WHENEVER D2.E.SIG/2. R=3 PRINT COMMENT S TWO INFINITE ROOTS$ END OF CONDITIONAL OR WHENEVER D1.L.O. WHENEVER D3.LSIG/2. R=1 PRINT COMMENT S TWO RPEAL ROOTS$ OR WHENEVER D3,..SIG/2. R=2 PRINT COMMENT S TWO IMAGINARY ROUTSS OR WHENEVER D3.FSIG/2. R=3 PRINT COMMENT $ TWO INFINITE ROOTS$ END OF CONDITIONAL OR WHENEVER Dl.G.O,.AND.C4.G.D4 WHENEVER 03.L.STG/2. R=1 PRINT COMMENT $ TWO REAL ROOTS$ OR WHENEVER D3,G.STG/2. R=2 PRINT COMMENT S TWO IMAGINARY ROOTS$ OR WHENEVER D3,E,SIG/2. R=3 PRINT COMMENT S TWO INFINITE ROOTS$ END OF CONDITIONAL OR WHENEVER D1.G.O.,ANDoC4.LED04 WHENEVER Cl.LE.SIG/2.AND*C4.GE.(1,/2,).OR. 1 Cl.GE,SIG/2..AND.C4.LE.(1,/2.) R=4 PRINT COMMENT S FOUR ROOTS TWO REAL AND TWO IMAGINARY$

103 PROGRAM 1 (Continued) OR WHENEVER Cl.L.SIG/2..AND.C4.L.(1./2.) WHENEVER D3.L.SIG/2. R=5 PRINT COMMENT S FOUR HEAL ROOTS$ OR WHENEVER D3.G.SIG/2. R=4 PRINT COMMENT S FOUR ROOTS TWO REAL AND TWO IMAGINARYS OR WHENEVER D3.E.SIG/2. R=l PRINT COMMENT S TWO REAL AND TWO INFINITE ROOTSS END OF CONDITIONAL OR WHENEVER Cl,GSIG/2,.ANDC4.G.(1,/2.) WHENEVER D3.L.SIG/2. R=6 PRINT COMMENT S FOUR IMAGINARY ROOTS$ OR WHENEVER D3,GSIG/2. R=4 PRINT COMMENT S TWO REAL AND TWO IMAGINARY ROOTSS OR WHENEVER D3.F.SIG/2. R=2 PRINT COMMENT S TWO IMAGINARY AND TWO INFINITE ROOTS$ END OF CUNDITIONAL END OF CONDITIONAL END OF CONDITIONAL WHENEVER R.E.1 X=1,000001 OMEGAR=OMEGAX.(X,C1,C4pDISIG) WHENEVER OMEGAR.L.O. THROUGH SERCH1, FOR X=1.100001O.1, 0MEGAR.G,0O.OR.X.G.209 WHENEVER X.G.19.75 PRINT COMMENT S NO SIGN CHANGE FROM NEG TO POS$ THROUGH PRINT1, FOR Z=1.000001po1, Z.G.20. PRINTI PRINT RESULTS Z. OMEGAX.(ZC1,C4pDlSIG) TRANSFER TO START END OF CONDITIONAL SERCHI OMEGAR=OMEGAX.C(XC1,C4,Dl1SIG) OR WHENEVER OMEGAR.G.O. THROUGH SERCH2P FOR X=1.200001,.2, OMEGARL.O..OR.XG,20, WHENEVER X.G.19.75 PRINT COMMENT $ NO SIGH CHANGE FROM POS TO NEGS THROUGH PRINT2, FOR Z=1.000001,.o2 Z. G.20 PRINT2 PRINT RESULTS Zp OMEGAX.(ZClC4D1l,SIG) TRANSFER TO START END OF CONDITIONAL SERCH2 OMEGAR=OMEGAX.(XC1lC4,D0ISIG) OTHERWISE PRINT COMMENT S X=1.000001 IS A ROOTS TRANSFER TO START END OF CONDITIONAL X(1)=XX.2

104 PROGRAM 1 (Continued) X(2)=X-. PRINT RESULTS X(1), OMEGAXo(X(1),C1,C4,DDISIG), 1X(2), OMEGAX.(X(2),C1,C4,DlSIG) FX(1)=OMEGAX.(X(1),CI,C4,DlISIG) FX(2)=OMEGAX.(X(2),C1,C4,D1,SIG) THROUGH LOOPXA, FOR ITX=I1,1ITX.G.ITMAX NUMER=X(2)*FX(l)-X(1)*FX(2) DENOM=FX(1)-FX(2) WHENEVER.ABS. DENOM.L..ABS.NUMER*EPS2 PRINT COMMENT S DENOMINATOR TOO SMALLS TRANSFER TO START END OF CONDITIONAL X(3)=NUMER/DENOM FX(3)=UMEGAX.(X(3),C1,C4,Dl,SIG) WHENEVER.ASS.FX(3).LE.EPS1 PRINT COMMENT S PROCEDURE CONVERGED$ PRINT RESULTS ITX,X(3),FX(3),R TRANSFER TO START END OF CONUITIONAL THROUGH LOOPXA, FOR I=lPlI.G,2 THROUGH LOOPXA, FOR J=I+1,lJ,(i.3 WHENEVER,ABS.FX(I).G..ABS.FX(J) TEMP=FX(I) FX(I)=FX(J) FX(J)=TEMP TEMP=X(I) X(I)=X(J) X(J)=TEMP LOOPXA END OF CONDITIONAL PRINT COMMENT $ NO CONVERGENCES TRANSFER TO START OR WHENEVER R.E.2 THROUGH SERCH3, FOR Y=.O0001.2# OMEGAI.L.O..OR.Y.G,40. WHENEVER Y.G.37.5 PRINT COMMENT S NO SIGN CHANGE FROM POS TO NEG$ THROUGH PRINT3, FOR M=.OOO01,2p MG.40. PRINT3 PRINT RESULTS MP OMEGAY.(MClC4,Dl,SIG) TRANSFER TO START END OF CONDITIUNAL SERCH3 OMEGAI=OMEGAY.(YPClC4,DISIG) WHENEVER OMEGAY.(.OOO1,ClC4,DlSIG).L.O. Y( 1 )=.0001 Y(2)=.2001 OTHERWISE Y(1)=Y-.4 Y(2)=Y.2 END OF CONDITIONAL PRINT RESULTS Y(1), OMEGAYC(Y(1),ClC4PD1,SIG), 1Y(2), OMEGAY.(Y(2),Cl1C4,DlSIG) FY 1)=OMEGAY (Y(l) 1)l,CIC41t)l1SIG)

105 PROGRAM 1 (Continued) FY(2)=OMEGAY (Y(2),CIC4,D1SIG) THROUGH LOOPYBFOR ITY=1,IITY.G ITMAX NUMER=Y(2)*FY(1)-Y(1)*FY(2) DENOM=FY(1)-FY(2) WHENEVER.ABS.DENOM,L,,..ABS.NUMER*EPS2 PRINT COMMENr $ DENOMINATOR TOO SMALLS TRANSFER TO START END OF CONDITIONAL Y(3)=NUMER/DENOM FY(3)=OMFGAY,(Y(3),,C,C4,D1,SIG) WHENEVER.ABS.FY(3).LE.EPS1 PRINT COMMENT $ PROCEDURE CONVERGEDS PRINT RESULTS ITY, Y(3),FY(3),* TRANSFER TU START ENO OF CONDITIONAL THROUGH LOOPYf, FOR 1zl1sI,G,2 THROUGH LOOPYB, FnR J=I+1,1,JG,3 WHENEVER,ABS. FY( I),G,ABSFY(J) TEMP=FY( I) FY(I)=FY(J) FY ( J)=TEMP TEMP=Y( I) Y( I )=Y(J) Y(J)=TEMP LOOPYB END OF CONDIIIONAL PRINT COMMENT $ NO CONVERGENCES OR WHENEVER R.E.4 THROUGH SEHCH4, FOR X=1.000001,.05,OMEGAR.L.O..OR,X,G.20. WHENEVER X.,G.1975 PRINT CUOMMENT NO SIGN CHANGE FOR REAL ROOTS$ fHROUGH PRINT4, FOP N=1.000001i.2, NG.20, PRINT4 PRINT RESULTS N, OMEGAX.(NpC1,C4,D1,SIG) TRANSFER TO IMAGR END OF CONDITIONAL SERCH4 OMEGAR=MEGAX X, C1C C4 DlSIG) WHENEVER OMEGAX.(1 000001,ClC4,D1,SIG) *LO 0 X(1)=1,.000001 X(2)=1 100001 OTHERWISE X( 1 )=X- 1 X ( 2)= X. 0 5 XC2)=X-.05 END OF CONDITIONAL PRINT RESULTS X(1), OMEGAX*(X(1),ClC4,D1I SIG), 1X(2), OMEGAX.(X(2),C1,C4,Dl,SIG) FX(1)=OMEGAX.(X( 1),C1,C4,DlSIG) FX(2)=OMEGAX.(X(2),C1IC4,D1SIG) THROUGH LOOPXC, FOR ITX=1lIITXG.ITMAX NUMER=X(2)*FX(1)-X(1)*FX(2) DENOM=FX(l)-FX(2) WHENEVER,ARS. DENOM,L,,ABS.NUMER*EPS2

106 PROGRAM 1 (Continued) PRINT COMMENT $ DENOMINATOR FOR REAL ROOTS TOO SMALL$ TRANSFER TO 1MAGR END OF CONDITIONAL X(3)=NUMER/DENOM FX(3)=OMFGAX.(X(3),ClC4,01,SIG) PRINT RESULTS X(3), FX(3) WHENEVER.ABS.FX(3).LE.EPS1 PRINT COMMENT $ PROCEDURE FOR REAL ROOTS CONVERGEO$ PRINT RESULTS ITX, Y(3), FX(3), R TRANSFER TU IMAGR END OF CnNOITIONAL rHROUGH LOOPXC, FOR I=1,II.G.2 THROUGH LOUPXCP FOR J=I+I,1,J,G.3 WHENEVER.ARS.FX(I).G..ABS.FX(J) TEMP=FX(I) FX(I)=FX(J) FX(J)=TEMP X(I)=X(J) X(J)=TEMP PRINT RESULTS X(1), FX(1), X(2), FX(2) LOOPXC END OF CONDITIONAL PRINT COMMENT S NO CONVERGENCE FOR REAL ROOTS$ IMAGR THROUGH SERCH5S FOR Y=.0001,.05,OMEGAI.L.O..OR.Y.G.40. WHENEVEi Y.G.40. PRINT COMMENT $ NO SIGN CHANGE FOR IMAG ROOTS$ THROUGH PRINTS, FOR T=.0001,.2, T.G,40, PRINT5 PRINT RESLULTS T, OMFGAY.(T,C1,C4,D1,SIG) TRANSFER TO STARS END OF CONDITIONAL SERCH5 OMEGAI=01MEGAY.(Y,C1,C4,D1,SIG) Y( 1 )=Y. 10 Y(2)=Y-.05 PRINT RESULTS Y(1), OMEGAY.(Y(1),C1,C4,D1,SIG), Y(2), OMEGAY. (Y(2),C 1,C4,01, SI G) FYC1)=OMEGAY.(Y(1),C1lC4,D1,SIG) FY(2)=OMEGAY.(Y(2),C1,C4,01,SIG) THROUIGH LOOPYD, FOR ITY=llITY*.GITMAX NUMER=Y(2)*FY(1)-Y(1)*FY(2) DENOM=FY(1)-FY(2) WHFNEVER,ABS. DENOM.L..ABSNUMER*EPS2 PRINT COMMENT $ DENOMINATOR FOR IMAG. ROOTS TOO SMALL$ iRANSFER TO START END OF CONDITIONAL Y(3)=NUMER/DENOM FY(3)=OMEGAY.(Y(3), C C4 01 SIG) PRINT RESULTS Y(3), FY(3) WHENEVER V ABS. FY(3).LE.EPS1 PRINT COMMENT $ PROCEDURE FOR IMAG, ROOTS CONVERGEDS PRINT RESULTS ITY,Y(3),FY(3),R

107 PROGRAM 1 (Concluded) TRANSFER TO START END OF CONDITIONAL THROUGH LOOPYD, FORJ=I+1,1,J.G.3 WHENEVER.ABS. FY(I).G..ABS.FY(J) THROUGH LOOPYDpFURI=,1 I.G.2 THROUGH LOOPYD, FOR J=I+I,1,J.G.3 TEMP=FY(I) FY( I )=FY(J) FY(J)=TEMP TEMP=Y( I) Y( I)=Y(J) Y ( J)=TEMP PRINT RESULTS Y(1), FY(1)J Y(2), FY(2) LOOPYD END OF CONDITIONAL PRINT COMMENT $ NO CONVERGENCE FUR IMAG ROOTS$ OTHERWISE PRINT COMMENT $ ROOTS ARE INFINITE OR WE HAVE 1 FOUR REAL OR FOUR IMAGINARY ROOTS$ END OF CONDITIONAL TRANSFER TO START END OF PROGRAM S COMPILE MAD, EXECUTE EXTERNAL FUNCTION (X,C 1C4,D1,SIG) INTERNAL FUNCTION T1,(X)=X*ELOG((SIG*X+1 * )/(SIG*X-1 )) INTERNAL FUNCTION T2.(X)=X*ELOG.((X+1,)/(X-1.)) ENTRY TO OMEGAX. OMEGAX=1,-C *T. (X)-C4*T2(X)+D1*TI. (X)*T2. X) FUNCTION RETURN OMEGAX END OF FUNCTION S COMPILE MAD, EXECUTE EXTERNAL FUNCTION(YC1i,C4,01DSIG) INTERNAL FUNCTION S,(Y)=2.*Y*ATAN.(, 1/(SIG*Y)) INTERNAL FUNCTION S2.(Y)=2.*Y*ATAN.(1./Y) ENTRY TO OMEGAY. OMEGAY=,-Cl*Sl. CY)C4*S2. (Y)+DI*Sl.(Y)*S2.(Y) FUNCTION RETURN OMEGAY END OF FUNCTION S DATA S DATA ITMAX=50, EPS1=I.E-7? EPS2=i.E-20, C1:.59023, C2=.05448, C3=.15963, C4=.44320 SIG=1.50951 *

108 PROGRAM 2 $ COMPILE MAD, EXECUTE R THIS IS THE SOLUTION FOR THE EXACT MILNE PROBLEM OM DIMENSION Y(16), S(16), W1(16), W2(16), T(16), T1(16), R(16), 1 R1(16), KY(16), KS(16), KY1(16), KS1(16), DIFF(32), 1 ZI(92)p KAPP(92), GU(96), ELG(92), ELGI(92), FK(8646,DIN), 1 FU(846,VIM), FUL(10), FUU(10), FUM(10), XFMI(92), OMAG1(92), 1 OMAG2(92), ALPHA2(828,VIN), SU(96), DIF(92)p FKK(92), 1 FKU(92), FTK(92), ALPHAl(92), ANGR1(164), FUMM(10), 1 FUMP(10), ANG(92), ANGO(92), ANGI(92), AGDIP2(92), 1 AGDIM2(92), AGDIP1(92), AGDIM1(92), Z2(92)',XFPI(92), 1 OPHI(92), OPSI(92), FBK(92), EXPO(92), Z3(20), ANGR(1B4) INTEGER P1, K, Pp J, I, L, ITMAX, ITMAXX, L1, P3, D, ZEE VECTOR VALUES VIM=2,0,O VECTOR VALUES VIN=2p0,0 VECTOR VALUES DlN=2,0,0 PROGRAM COMMON FU, FUL, FUM, FUU, N, ELG, ZiP H1, H2, P3, D, I SG R VECTOR VALUES FOR GAUSS QUADRATURE VECTOR VALUES X(1)=.0950125098,.2816035508,.4580167777, 1.6178762444,.7554044084S.8656312024,,9445750231, 2.9894009350 VECTOR VALUES W(1)=.1894506105,.1826034150,,1691565194, 1.1495959888,.1246289713,.0951585117,,0622535239, 2.0271524594 FTRAP. R INTERNAL FUNCTIONS FOR KAPPA FUNCTION CALCULATION INTERNAL FUNCTION F.(V1,Y1,KAP)=(l.-Vl*Vl*XO2)*Vl/((l1.Vl*V1/ 1 (ETA*ETA))*(Vl+Y1)*KAP) INTERNAL FUNCTION FG,(Vi,Yl,KAP)=1./((VI+Yl)*KAP) START READ AND PRINT DATA R NOTE THAT P, P3 AND D ARE DECLARED TO 8E INTEGERS P=N P3=N3 D=P+P3 VIM(2)=D+2 VIN(2)=D DIN(2)=0+2 R THIS PART OF PROGRAM OBTAINS A KAPPA FUNCTION WITHIN A R SPECIFIED ACCURACY EPS R PARAMETERS FOR GAUSS QUADRATURE INTERVAL FROM 1./SIG TO 1. A=(SIG-1,)/(2.*SIG) B=(SIG+1I)/(2.*SIG) RCALCULATION OF DETERMINANT OF C MATRIX SEE EQ.(2.2) D1=C *C4-C2*C3 R R,H,S. OF EQ.(A,14) MULTIPLIED BY C Q=D1/SIG*(SIG*(ELOG.((SIG+1.)/(SIG-I.)))+ELOG.(SIG*SIG-1.)) R SEE EQ.(A.2) AND EQ.(A.5),(X(O) SQUARED) X02=1,/((1.-2*C/1/ SIG-2.*C4+4.*D1/SIG)*(ETA*ETA)) THROUGH EVAL, FOR K=1,1' K.G.8

o109 PROGRAM 2 (Continued) R NEW VARIABLES CALCULATED FOR GAUSS QUADRATURE, REGION 1 R SEE REF, (18) Y(K)=XCK)/(2.*SIG)+1,/(2,*SIG) Y(K+8)=-X(K)/(2.*SIG)+o./(2,*SIG) W1(K)=W(K)/(2.*SIG) Wl(K+8)=W(K)/(2.*SIG) R NEW VARIABLES CALCULATED FOR GAUSS QUADRATURE, REGION 2 S(K)=X(K)*A+B S(K+8)=-X(K)*A+R W2(K)=W(K)*A EVAL W2(K+8)=W(K)*A THROUGH EVAL2, FOR J=l,1, J.G,16 R FUNCTIONS IN EQ.(A.13) THAT ARE DEPENDENT ON INTEGRATION INR TERVAL (MU) ONLY (NOTE GAUSS WEIGHTS ARE INCLUDED) T(J)=(C4+CI-Ul*Y(J)*ELOG.((l.+Y(J))/(1.-Y(J))))*( 1-Y(J)*Y(J) 1 *XO2)*W1(J)/(1.-Y(J)*Y(J)/(ETA*ETA)) Tl(J)=Wl(J)*01*ELOG.((1.+SIG*Y(J))/(1.-SIG*Y(J))) R(J)=W2(J)*C4*(1.-S(J)*S(J)*X02)/ (1.-S(J)*S(J)/(ETA*ETA)) EVAL2 R1(J)=W2(J)*Dl*ELOG.((SIG*S(J)+lo)/(SIG*S(J)-'1)) INTERNAL FUNCTION (7) ENTRY TO G. R THIS INTERNAL FUNCTION CALCULATES KAPPA FOR A SPECIFIED Z SUM=O THROUGH EVAL3, FOR L=1,1, L.G.16 EVAL3 SUM=SUM+T(L)*FG.(Y(L),ZKY(L))-Tl(L)*(F.(Y(L),ZKY(L))-F.(lo/ i SIG,ZKSG))+R(L)*FG,(S(L)ZKSC(L))-Rl(L)*(F.(S(L),ZKS(L))2 F.(1./SIG,Z,KSG)) FUNCTION RETURN 1.-Z*SUM+Q*Z*F.(l./SIGZ,KSG) END OF FUNCFION R AN ITERATION PROCESS COMMENCES AT THIS POINT WHERE FOR R INITIAL VALUS OF KAPPA WE SET KAPPA=1. THE ITERATION R CONTINUES UNTIL THF MAXIMUM DIFFERENCE BETWEEN ANY TWO R VALUES OF KAPPA IS LESS THAN EPS MAX=1. KSG=I. THROUGH EVAL4P FOR L=1,1, L.G,16 KY(L)=1. EVAL4 KS(L)=1. THROUGH EVAL5, FOR J=1,1. JGITMAXOR.MAX.L.EPS THROUGH EVAL6, FOR K=1,1, K,G,16 KY1(K)=G,(Y(K)) KS1(K)=G,(S(K)) DIFF(K)=,ABS.(KY(K)'KY1(K)) EVAL6 DIFF(K+16)=.ABS.(KS(K)-KS1(K)) KSG1=G.(1./SIG) MAX=DIFF(1) THROUGH ALPHA, FOR L.=2,l. L.G.32 ALPHA WHENEVER DIFF(L),G.MAX, MAX=DIFF(L) WHENEVER J.G.ITMAX PRINT COMMENT $ MAXIMUM NUMBER OF ITERATIONS FOR KAP EXCEEDS

110 PROGRAM 2 (Continued) PRINT RESULTS MAX, J TRANSFER TU START OR WHENEVER MAX.L.EPS PRINT COMMENT $ PROCEDURE FOR KAPPA FUNCTION CONVERGED$ PRINT RESULTS MAX, J, KY(1)...KY(16), KS(l).,.KS(16) TRANSFER TO GO OTHERWISE KSG=KSG 1 THROUGH SWITCH, FOR L1=1,1 L1.G.16 KY(L1)=KY1(Li) SWITCH KS(L1)=KS1 (CL) EVAL5 END OF CONDITIONAL R THIS CONCLUDES THE CALCULAIION UF KAPPA GO CONTINUE R NUMBER NEEDED FOR X-FUNCTION AND GAMMA FUNCTION SR=1./SQRT.(X02) R SET-UP MESH SPACING MID)POINT VALUES SG=1 /SIG R INTERVAL SPACING FOR REGION, 1 H1=SG/N R INTERVAL SPACING FOR RE(lION 2 H2=(1.-SG)/N3 K=1 R MIDPOINT VALUES OF MU VARIABLE DETERMINED THROUGH COMP1, FOR I1=1.,1. K.GD WHENEVER K.LE.P Zl(K)=(U-.5)*H1 OTHERW ISE ZI(K)=Zl(P)+.5*H1+(U-N-.5)*H2 END OF CONDITIONAL COMP K=K+I PRINT RESULTS SG, H1, H2, Z1(i)...Z1(D) R VARIOUS CONSTANTS NEEDED LATER R IF WE WANT TO SOLVF THE EXACT CONSTANT SOURCE PROBLEM WE R INSERT CARDS NO, 1S, 2S, 3S AND 4S PI=3,14159265 R HERE WE CALCULATE KAPPA FIRST AT 1,/SIG AND 1. AND THEN R GAMMA AT THESE SAMF VALUES KAPSG=G. (SG) GU(D+1)=SG*ETA*ETA*ETA(SR+SG)*XO2/((ETA*ETA-SG*SG)*KAPSG) KAONE=G,(1.) GU(D+2)=ETA*ETA*(SR+1.)*XO2/((ETA*ET A-1.)*KAONE) R KAPPA AT TWO DISCRETE ROOTS KAPM=G,(ETA) KAPP=G (-ETA) R FIRSf TERM IN EQ.,(4,10) FOR DISCRETE ROUTS ONLY ELGETP=ETA*ELOG*(1.+SG/ ETA) R MILNE CARD ONLY ELGETM=-ETA*ELOG,(1.-SG/ETA) R SEE EQUATION FOLLOWING EQ.(2.20)

111 PROGRAM 2 (Continued) ANGC=C4-Dl*(ELGETP+ELGETM) R THE LAST TERM IN EQ,(4.10) CALCULATED FOR 1. AND 1./SIG ELGONE=ELOG.(1, +SG) ELGSG=SG*ELOG.(2,) R EQS,(6.2) AND (6.3) EVALUATED AT 0, 1*/SIG AND 1, FUL(1)=O, FUL(2)=0, FUM(1)=GU(D+1)*ETA/(FTA+SG)*(ELGETM-ELGSG) 5M FUM(2)=GU(D+1)*ETA/(ETA-SG)*(ELGETP-ELGSG) FUU(1)=GU(0+2)*ETA/(ETA+le)*(ELGETM-ELGONE) 6M FUU(2)=GU(D+2)*ETA/(ETA'1.)*(ELGETP-ELGONE) R X-FUNCTION EVALUATED AT DISCRETE ROOTS (SEE EQ.(A,15) R MILNE CARD ONLY XFM=KAPM/(SR+ETA) XFP=KAPP/(SR-ETA) PRINT RESULTS KAPSGp GU(D+I), KAONE, GU(D+2), KAPM, ELGETM, 1 ELGONEP ELGSG, XFM, XFP R INTERNAL FUNCTION FOR THE BOUNDARY VALUES OF OMEGA EQ.(2,21) INTERNAL FUNCTION OM4,(Z)=Z INTERNAL FUNCTION 01,(Z)=Z*ELOG.((1,+Z)/(.1-Z)) INTERNAL FUNCTION 02,(Z)=Z*ELOG.((1,+SIG*Z)/(1.-SIG*Z)) INTERNAL FUNCTION tl3.(Z)=Z*ELOG,((SIG*Z+1,)/(SIG*Z-1.)) THROUGH COMP2, FOR K=:1,1 K,GD R COMPUTATION TO CHANGE MU VARIABLE TO ANGLE IN DEGREES FOR R EASIER PLOTIING ANRP=ARCCnS (Z1(K)) ANRM=ARCCOS (-Z1(K)) ANGR(K)=ANRP*360./(2.*PI) ANGR(K+D) =ANRM*360./(2.*PI) WHENEVER Zl(K).L,SG Zt 1=Z1(K)*SIG ANRP1=ARCC S,(Z11) ANRM1=ARCCOS.,( -Z11) ANGR1(K)=ANRP1*360,/(2.*Pi) ANGRI(K+P) =ANRM1*360,/(2.*PI) END OF CONDITIONAL R CALCULATION OF KAPPA THEN GAMMA FUNCTION FOR ALL MIDPOINT R VALUES (SEE EQ.(A,17)) KAPP(K)=G,(ZI(K)) GU(K)=ZL(K)*ETA*ETAETA*(SR+Z(K))*X2/((ETA*ETA-Zl(K)*Z(K))* 1 KAPP(K) ) R CALCULATION OF LOGARITHM TERM IN EQ.(6.6) ELG(K)=ELOG,((1.-Z1(K))/Z1(K)) R LAST TERM OF EQ,(4,10) EVALUATED AT ALL MIDPOINT VALUES ELG(K)=ZI(K)*ELOG (1 +SG/Z1 (K)) R CALCULATION OF EQ.(6.2) AND EQ.(6.3) AT ALL MIDPOINT VALUES FKU(K)=ETA/(ETA-Z(K) )*(ELGETPELGETPELG1 (K)) FU(C2K)=GU(K)*FKU(K) FKK(K)=ETA/(ETA+Z1(K))*(ELGETM-ELG1(K)) 7M FU(1pK)=GU(K)*FKK(K)

112 PROGRAM 2 (Continued) R EQS.(4.22) EVALUATED FOR EACH REGION WHENEVER Z1(K),L.SG S1=01.(Z1(K)) S2=02.(Z1(K)) S3=OM4.(Zl(K)) S4=S3'S3 DENOMI=(1.'C4*SI-I*S2+Dl*SI*S2-PI*PI*DI*S4).P.2+PI*PI*SA* 1 (C4+CI'Dl*(S2+S1)).P.2 NUMNI=S3*(C4+Cl-Dl*(Sl+S2)) NUMPI=1.-C4*S1-Ci*S2+D1*SI*S2-PI*PI*DI*S4 OMAG1(K) = NUMNI/(DENOHM*GlJ(K)) OMAG2(K) =NUP1l/DFNOM1 R FIRST TERM OF EQ9S(4.42) AND (5.8) EVALUATED R MILNE CARD ONLY XFMI(K)=ETA*XFM/(ETA+ZI(K))*OMAGI(K) XFPI(K)=OMAG1(K)*ETA*XFP/(ETA-ZI(K)) R FUNCTIONS NEEDED LATER FOR ANGULAR DISTRIBUTIONS ANG(K)=C4-Dl*S2 ANGI(K)=C1-Dl*Sl ANGO(K)=NUMP1+PI*PI*D1*S4 R SAME REMARKS AS ARnVE RUT FOR REGION 2 OTHERWISE T1=Ol.(Zl(K)) T2=O3.(ZL(K)) T3=OM4.(ZI(K)) T4=T3*T3 DENOM2=(1.-C4*T1-Cl*T2+DI*Tl*T2)P.2+PI*PI*T4*(C4-Dl*T2).P.2 NUMN2=T3*(C 4-D*'T2) NUMP2=1.-C4*T1-C1*T2+n1*Tl*T2 OMAG1(K) =NUMN2/(DENOM2*GU(K)) OMAG2(K) =NUMP2/DENOM2 R MILNE CARD ONLY XFMI(K)=ETAXFM/ETA*XFM/(T+Z1(K))*OMAG1(K) XFPI(K)=OMAG1(K)*ETA*XFP/(ETA-Z1(K)) ANG(K)=C4-D1*T2 ANGO(K)=NIJMP2 ANG1(K)=O. COMP2 END OF CONUITIONAL R INTEGRALS IN EQ.(5.9A) EVALUATED EXCEPT INTEGRAL WITH ALPHA2 R A SUBROUTINE IS USFD WHICH IS FOUND AT END OF PROGRAM INTEGM=INT.(1) INTEGP=INT.(2) R EQ.(4.35) EVALUATED APLDEN=-ETA*XFP-DI*INTEGP R EQ.(6.16) EVALUATED APLUS1=(-El'A*XFM+1D*INTEGM)/APLDEN 8M PRINT RESULTS KAPP(1)...KAPP(D), GU(1)...GU(D), 1 ELG(1)...ELG(D), AKIGOC1)...ANGO(D), FKK(1)*,.FKK(D), 1 FU(1I1)...FU(1JD), XFMI(1)..,XFMI(D), FUM(1), FUU(1), 1OMAGl(l)...OMAG1(D)0 OMAG2(1)...OMAG2(D),

113 PROGRAM 2 (Continued) 1 XFPI(1),..XFPI(D), APLUlS1, APLDEN, INTEGP, INTEGM THROUGH COMP3, FOR K=Il, K.GD R EQS.(4,42) AND (5.8) EVALUATED, PVI. IS A SUBROUTINE LISTED R AT THE END WHICH CALCULATES PRINCIPAL VALUE INTEGRALS OPHI(K)=D1*(OMA1GI(K)*PVI.(1,2K)'OMAG2(K)*FKU(K))+XFPI(K) OPSI(K)=-D1*(OMAGI(K)*PVI.(1,1,K)-OMAG2(K)*FKK(K))-XFMI(K) 9M R INITIAL VALUE OF ALPHA2 CALCULATED SEE EQ.(6.17) ALPHA2(C1K)=OPSI(K)-APLUS1*OPHI(K) COMP3 PRINT RESULTS Zt(K)p ALPHA2(1,K), OPHI(K), OPSI(K) R INTEGRAND IN SECOND TERM UF NUMERATOR OF EQ.(4.43A) R CALCULATED EXCLUDING ALPHA2 THROUGH COMP4, FOR J=lpl, JOG.D THROUGH COMP4, FOR K=Il, K.GU.+2 WHENEVER K.E.D+1 FK(JK)=Z1(J)/(Zl(J)-SG)*(ELGI(J)-ELGSG) OR WHENEVER K.E.D+2 FK(JK)=Zl(J)/(ZI(J)-I.)*(ELG1(J)-ELGONE) OR WHENEVER J.E.K FK(JK)=ELGI(K)-Z1(K)/(1.+SIG*Z1(K)) OTHERWISE FK(J,K)=ZI(J)/(Zi(J)-ZI(K))*(ELGI(J)-ELG1(K)) COMP4 END OF CONDITIONAL R HERE WE START THE ITERATION PROCEDURE FOR APLUS AND ALPHA2 MAXX=1. THROUGH COMP5S FOR J=2,1, J.G.ITMAXX.OR.MAXX.L.EPSI THROUGH COMP6, FOR K=I,1p K.GD+2 R FIRST INTEGhATION PERFORMED IN EQ.(4,43A) SUM10= SUM 2=0 THROUGH COMP7, FOR I=v1,1 I.G.D WHENEVER I.LE.P SUM1=SUMI+FK(I,K)*ALPHA2(J-l,I) OTHERWISE SUM2=SUM2+FK(I,K)*ALPHA2(J-l,I) COMP7 END OF CONDITIONAl. SU(K)=SUM1*Hi+SIJM2*N2 COMP6 FU(J,K)=SU(K)*GU(K) R SECOND INTEGRATION PERFORMED IN EQ.C(4.43A) INTG=INT.(J) R NEW VALUE OF A+ CALCULATED APLUS=APLUS1+INTG/APLDEN*D1 PRINT RESULTS FlU(J,l),..FU(J,D+2), SU(1)...SU(D+2), INTG FUL(J)=O, R FUNCTIONS NEEDED FOR NEXT PRINCIPAL VALUE INTEGRAL FUM(J)=FU(JPD+I) FUU(J)=FU(JO+2) PRINT RESULTS FUL(J), FUM(J), FUU(J), APLUS THROUGH COiAP8b FOR K=11, K.G.* R NEW VALUES OF ALPHA2 ALPHA2(JK)=UPSI(K)-APLUS*OPHI(K)+Dl*(-OMAGl(K)*PVI.(lJK)+

114 PROGRAM 2 (Continued) 1 OMAG2(K)*SU(K)) R DIFFERENCE 8ETWEEN PREVIOUS AND NEW VALUES OF ALPHA2 COMP8 DIF(K)=ALPHA2(JK)-ALPHA2(J-1 K) PRINT RESULTS DIF(1),,.DIF(D) R LOCATE MAXIMUM DIFFERENCE MAXX=.ABS,(DIF(1)) THROUGH COM1, FOR L=2,1J LwGD COM1 WHENEVER,ABS.(DIF(L)),G.MAXX, MAXX=.ABS.(DIF(L)) WHENEVER J.E.ITMAXX PRINT COMMENT $ MAXIMUM NUMBER OF ITERATIONS FOR ALPHA2 EXD$S PRINT RESULTS MAXX, J TRANSFER TO START R EPSI IS CRITERION FOR TERMINATION OF ITERATION PROCEDURE OR WHENEVER MAXX,L,EPSI PRINT COMMENT $ PROCEDURE FOR ALPHA2 CONVERGEDS TRANSFER TO GOO OTHERWISE PRINT RESULTS MAXX, J, ALPHA2(J,1).,.ALPHA2(JJD) COMP5 END OF CONDITIONAL GOO PRINT RESULTS ALPHA2(J,1)...ALPHA2(J,D)p Zl(1)...Z1(D), MAXX R FINAL CONVERGED CALCULATION OF A+ (SEE EQ,(5,9C)) THROUGH COM5, FOR K=1,p, K.G.D SUM 1=0 SUM2=O. THROUGH COM6P FUR I=1,1, I.G,D WHENEVER I.LE.P SUM1=SUM1+FK(IK)*ALPHA2(JI) OTHERWISE SUM2=SUM2+FK(I,K)*ALPHA2(JI) COM6 END OF CONDITIONAL SU(K)=SUM1*H1+SUM2*H2 COM5 FU(J,K)=SIJ(K)*GU(K) ALPAI=INT.(J) PRINT RESULTS SU(1)..,SU(D), FU(J,1)....FU(JD) THROUGH COM4, FOR K=1,1 K.G.D COM4 FU(J+1,K)=Z1(K)*ALPHA2(JK) APLUS=APLUS1+ALPAI/APLDEN*D1 FUL(J+1)=O. FUM(J+1)=O. FUU(J+1)=O, THROUGH COM3, FOR K=i,1. KJ G*D ALPHA1(K)=-ETA*C2*C1*(1./(ETA+Z1(K))+APLUS/(ETA-Z1(K)))-C2*C1 lOW 1*PVI.(1,J+1,K) COM3 PRINT RESULTS Z1(K), ALPHA1(K) R INTERNAL FUNCTION FOR QUADRATIC INTERPOLATION INTERNAL FUNCTION QUAD.( ZXOFOXlFl1,X2,F2)=FO+(Z-XO)*(Fl1 FO)/(XXO)+(ZXOO) XO)*Z-X1)/(X2XO)*((F2-F)/(X2Xl)(Fl-FO) 1 /(XI-XO)) R MILNE CARD ONLY EXTRAP=-ETA/2.*ELOG.(-1./APLUS)

115 PROGRAM 2 (Continued) PRINT RESULTS APLUSP EXTRAP, ALPAI R CALCULATION OF ANGULAR DISTRIBUTION FOLLOWS(SEE EQS.(5,14) R AND (5.15)) THROUGH CO4, FOR L=l,l1 L.G.ZEE THROUGH C06, FOR K=,ll, K.G.D C06 EXPO(K)=EXP.(-Z3(L)/ZI(K)) THROUGH COl FOR K=l,l. K.G.-D SUMI=O SUM2=0. SUM3=0. SUM4=0O SUM5=0O SUM6=0. FU(l,K)=ZI(K)*ALPHAI(K)*EXPO(K) FU(2pK)=Z1(K)*ANG(K)*ALPHA2(JK)*EXPO(K) FU(3,K)=Z1(K)*ALPHA2(J,K)*EXPO(K) THROUGH C03, FOR I=l,1, I.GD WHENEVER I.LE.P SUMI=SUMI+ZI(I)*ALPHAI(I)/(ZI(I)+Z1(K))*EXPO(I) SUM2=SUM2+Zl(I)*ANG(I)*ALPHA2(J,I)/(ZI(I)+ZI(K))*EXPO(I) SUM3=SUM3+Zl(I)*ALPHA2(J,I)/(Z1(I)+ZI(K))*EXPO(I) WHENEVER Z1(K).G*SG,SUM6 =SUM6 +Zl(I)*ALPHAl(I)/(Zl(I)-Zl(K)) 1 *EXPO(I) OTHERWISE SUM4=SUM4+Z1(I)*ANG(I)*ALPHA2(J,I)/(Z1(I)+ZI(K))*EXPO(I) SUM5 =SUM5 +Z1(I)*ALPHA2(J,I)/(ZI(I)+ZI(K))*EXPO(I) C03 END OF CONDITIONAL FKK(K)=SUM1*H1 FKU(K)=SUM2*HI+SUM4*H2 FTK(K)=SUM3*Hl+SUM5*H2 COI FBK(K)=H1*SUM6 FUL(l)=O, FUM(1)=SG*QUAD.(SGZI(P-1),ALPHAI(P-l),Z1(P),ALPHA1(P),ZI(P+I I),ALPHA1(P+l))*EXP,(-Z3(L)/SG) FUL(2)=0. SGl=ALPHA2(J,P'1)*ANG(P-1) SG2=ALPHA2(J,P)*ANG(P) SG3=ALPHA2(JP+I)*ANG(P+1) FUM(2)=SG*QUAD.(SG,ZI(P-I),SG1,Z1(P),SG2,Zl(P+l),SG3)*EXP.( 1 -Z3(L)/SG) PRINT RESULTS SG1, SG2, SG3, FUM(2) FUU(2)=O. FUL(3)=O. FUM(3)=O, FUU(3)=O. R MILNE CARD ONLY EXPETM=EXP.(Z3(L)/ETA) EXPETP=EXP.(-Z3(L)/ETA) THROUGH C02, FOR K=1l,l K.G.D AGDIM2(K)=ETA*ANGC/(ETA-Zl(K))*EXPETM+ETA*ANGC*APLUS/(ETA+Zl( l1M

116 PROGRAM 2 (Continued) 1 K))*EXPETP-L1/(CC*C2)*FKK(K)+FKU(K) WHENEVER Z1(K).L.SG AGDIP2(K)=ETA*ANGC/(ETA+Z1(K))*EXPETM+ETA*ANGC*APLUS/(ETA-Zl( 12M 1 K))*EXPETP-O1/(Cl*C2)*PVI.(2pl, K)-ALPHAl(K)/(C2*Cl)*ANGl(K)* 1 EXPO(K)+PVI.(1,2,K)+ANGO(K)*ALPHA2C(JK)*EXPO(K) AGDIPI(K)=ETM+C/(S*(TA+ZK *XPT+C2*ETA*APLUS/(SIG* 13M 1 (ETA-Zl(K)))*EXPETP+ALPHA1(K)/CI,*SIG)*EXPO(K)+C2/SIG*PVI.(( 1,3,K) AGDIMI(K)=C2*ETA/(SIG*(CTA-ZlCK)))*EXPETM+C2*ETA*APLUS/(SIG* 14H 1 (FTA+Z(K) ) )*EXPETP+C2/SIG*FTK(K) PRINT RESULTS ANGRI(K), AGDIP1(K)p ANGR1(K+P), AGDIM1(K) OTHERWISE AGDIP2(K)=ETA*ANGC/(ETA4Zl(K))*EXPETM+ETA*ANGC/(ETA-ZI(K))* 15M 1 APLUS*EXPETP-D1/(CI*C2)*FBK(K)+PVI.(C1,2K)+ANGO(K)*ALPHA2(J, 1 K)*EXPO(K) END OF CONDITIONAL C02 PRINT RESULTS ANGR(K), AGDIP2(K)p ANGR(K+D)p AGDIM2(K) AGDIP2(P)=QUAD.(Zl(P), Z1(P-2), AGDIP2(P-2),Zl(P-1)pAGDIP2(P 1 -1),Zl(P+2)vAGDIP?(P+2)) AGDIP2(P+1)=QUAD (Zl(P+1), ZI(P-1 )AGDlP2(P-1 )Zl(P+2),AGDIP2 1 (P+2)Z1 (P+3) AGDIP2(P+3)) PRINT RESULTS AGDIP2(P), AGDIP2(P+1) R APPROPRIATE INTEGRALS TO OBAIN THE TOTAL FLUX AND CURRENT R FOLLOWED BY A CALCULATION OF THE ASYMPTOTIC FLUX AND CURRENT R SEE EQS.(5.16) SUM 1=O SUM2=0. SUM3=0, SUM4=0. SUM5=0, SUM6=0. SUM7=O. SUM8=O0. SUM9=0. SUMIO 1 0 SUM11=O. SU M 12=O. THROUGH C05, FOR K=l,1, K.G.D WHENEVER K.L.P+I Z2(K)=Z1(K)*SIG SUM1=SUMI+AGDIMi(K) SUM2=SUM2+AGDIP1(K) SUM3=SUM3+AGDIM2(K) SUM4=SUM4+AG0IP2(K) SUM5=SUM5-Z2(K)*AGDIMI(K) SUM6=SUM6+Z2(K)*AGDIP1(K) SUM7=SUM7-Z1(K)*AGDIM2(K) SUM8=SUM6+Z(CK)*AGDIP2(K) OTHERWISE SUM9=SUM9+AGDIM2(K)

117 PROGRAM 2 (Continued) SUM10=SUM10O+AGDIP2(K) SUM1I=SUM11-Zl(K)*AGDIM2(K) SUM12=SUM12+Z1(K)*AGDIP2(K) C05 END OF CONDITIONAL FLNEG1=SUMI*HI*SIG FLPOSI=SUM2*H1*SIG 10TFL1=FLNEG1+FLPOS1 CUNEG1=SUM5*HI*SIG CUPOS1=SUM6*H1*SIG TOTCU1=CUNEG1+CUPOS1 FLNEG2=SUM3*HI+SUM9*H2 FLPOS2=SUM4*H1+SUIMO*H2 TOTFL2=FLNEG?+FLPOS2 CUNEG2=SUM7*H1+SUM1I*H2 CUPOS2=SUMS*Hi+SlM 12*H2 TOTCU2=CUNEG2+CUPOS? PRINT RESULTS FLNEG1, FLPOS1, TOTFL1I CUNEG1, CUPOSI, TOTCU1l 1 FLNEG?, FLPUS2P TOTFL2, CUNEG2, CUPOS2, TOTCU2 ARCSIG=ELGETP+ELGETM ASYFLI=C2*ARCSIG*(EXPETM+APLUS*EXPETP) 16M ASYCUI=C2*ETA*EXPETM*(2.-SIG*ARCSG*ACSIG)+C2*ETA*EXPETP*APLUS*( 17TM 1 SIG*ARCSIG_2.) ARCETA=ETA*ELOG.((ETA+1.)/(ETA-1,)) ASYFL2=ANGC*ARCETA*(EXPETM+APLUS*EXPETP) 18M ASYCU2=ETA*ANGC*EXPETM*(2,-ARCETA)+ETA*ANGC*APLUS*EXPETP*( 19M 1 ARCETA-2,) CO4 PRINT RESULTS ASYFLI, ASYCU1J ASYFL2, ASYCU2 TRANSFER TO START END OF PROGRAM S COMPILE MAD, EXECUTE R SUBROUTINE FOR THE CALCULATION OF INTEGRALS BY MIDPOINT R APPROXIMATION- TRAPEZOIDAL RULE EXTERNAL FUNCTION INT,(J) INTEGER K, I, Q4, P, P1, P3, 0D J DIMENSION FU(846,VIM). FUL(1O)p FUM(1O), FUU(1O)p ELG(92), 1 Z1(92) VECTOR VALUES VIM=2,0,O PROGRAM COMMON FUp FUL, FUM, FUUP N, ELG, ZI, H1, H2# P3, D, 1 SG VIM(2)=D+2 P=N SUMI=O0 SUM2=O0 THROUGH EVAL, FOR K=1,1, K.G.D WHENEVER K.LE.P SUM1=SUM1+FU(JPK) OTHERWISE SUM2=SUM2+FU(JK) EVAL END OF CONDITIONAL FUNCTION RETURN SUM1*H1+ SUM2*H2

118 PROGRAM 2 (Continued) END OF FUNCTION S COMPILE MAD, EXECUTE R SUBROUTINE FOR CALCULATING SINGULAR INTEGRALS SEE CHAPTER 6 EXTERNAL FUNCTION PVl.(J,I,Q4) DIMENSION FU(846,VIM), FUL(10O) FUM(1O), FUU(1O)p ELG(92), 1 Z1(92) VECTOR VALUES VIM=2PO,0 INTEGER K, IJ Q4, Pp P1J P3, DP J PROGRAM COMMON FUJ FUL, FUM, FJUr N, ELG, Z1P H1, H2, P3, D, 1 SG P=N VIM(2)=D+2 SUM=O, SUM1=O. SUM2=O. THROUGH COMPI, FOR K=1,1, K.G.D WHENEVER K.E.Q4 WHENEVER Q4,E.1 SUM=SUM+1e/3.*(FU(I,2)+3.*FIJ(I,1)-4.*FUL(I)) OR WHENEVER Q4.E.P SUM=SU+1 /3.*(4*FUM( I )-3.*FU(IP)-F U(IP-1)) OR WHENEVER Q4.E.P+1 SUM=SUM+1./3.*(FU(IP+2)+3.*FU( IP+I)-4*FUM(I)) OR WHENEVER Q4.E.D SUM=SUM+./3.*(4.*FUU( I) 3.*FU(I,D)'FU(ID-1)) OTHERWISE SUM=SUM+1./2.*(FU(I, Q4+1)'FU(I, Q4-1)) END OF CONDITIONAL OTHERWISE WHENEVER K.LE.P SUMI=SUMI+(FU(IK)-FU(IQ4))/(Zi(K)-ZI(Q4)) OTHERWISE SUM2=SUM2+(FU(IK)-FU(IQ4))/(Z1(K)-Zl(Q4)) END OF CONDITIONAL COMP1 END OF CONDITIONAL WHENEVER J.E.1 FUNCTION RETURN SUM+FU(IQ4)*ELG(Q4)+H1*SUM1+H2*SUM2 OTHERWISE R SINGULAR INTEGIAL FOR REGION 1 ONLY FUNCTION RETURN SUM+FIJ(IQ4)*ELOG.(SG/Z1(Q4)-I,)+HI*SUM1 END OF CONDITIONAL END OF FUNCTION S DATA C1=.61320P C2=.04404, C3=.17023, C4=,43685, SIG=1.65809r ETA=2.595825, EPS=*000001, N=22., ITMAXX=10. EPSI=.000001 N3=70. Z3(1)=0..5J,1.,1.5,2. 3.,5.,5,8. 1O,.20. ZEE=10, ITMAX=10 * R THE SOLUTION TO THF CONSTANT PROBLEM IS OBTAINED BY R REPLACING IN THE PREVIOUS PROGRAM THE CARDS WITH M IN R COLUMN 77 BY THE FOLLOWING CARDS WITH S IN COLUMN 77. THE R NUMBERS MUST ALSO MATCH. FOR EXAMPLE 5M IS REPLACED BY 5S5

119 PROGRAM 2 (Concluded) R ETC, ALSn WE MAY RFMOVE FROM THE PREVIOUS DECK ALL CARDS R PRECEDED BY THE REMARK $MILNE PROBLEM ONLY$. NOTE THAT R THE FIRST FOUR CARDS AFTER THE FIRST ARE PLACED IN THE R PREVIOUS DECK AT ONE DESIGNATED LOCATION, R THIS IS THE SOLUTION FOR THE EXACT CONSTANT SOURCE PROBLEM OS CS1=2,*C2/((1,-2.*C4)*(SIG-2.*C1)-4.*C2*C3) 1S CS2=(SIG-2.*Cl)/((1.-?.*C4)*(SIG'-2*Cl)-4,*C2*C3) 2S CSW=-CS2-CSI*D1/C2+SIG*CSI*C4/C2 3S CS3=SIG*CS1*D1/C2 I S FUM(1)=GU(I)+l)*(CSW-CS3*ELGSG) 5S FUU(I)=GU(D+2)*(CSW-CS3*ELGONE) 6S FKK(K)=CSW-CS3*ELG1(K) 7S APLUSI=(-SIG*CS1/C2+INTEGM)/APLDEN 8S OPSI(K)=-OMAG1(K)*PVI*(llK)+OMAG2(K)*FKK(K) 9S ALPHAI(K)='CSI*SIG*CI-ETA*C2*Cl*APLUS/(ETA-Zl(K))-C2*C1 10S AGDIM2(K)=CS2 +ETA*ANGC*APLUS/(ETA+Z1(K llS AGDIP2(K)=CS2 +ETA*ANGC*APLUS/(ETA-Zi( 12S AGDIP1(K)=CS1 +C2*ETA*APLUS/(SIG* 13S AGDIMI(K)=CS1 +C2*ETA*APLUS/(SIG* 14S AGDIP2(K)=CS2 +ETA*ANGC/(ETA-Z1(K))* 15S ASYFL1=2.*CSI+C2*ARCSIG*APLUS*EXPETP 16S ASYCUI= +C2*ETA*EXPETP*APLUS*( 17S ASYFL2=2,*CS2+ANGC*ARCETA*APLUS*EXPETP 18S ASYCU2= +ETA*ANGC*APLUS*EXPETP*( 19S

120 PROGRAM 3 S COMPILE MAD, EXECUTE, DUMP, PRINT OBJECT R PROGRAM FOR SOLUTION OF THE MILNE PROBLEM OM R SPHERICAL HARMONICS-Pi APPROXIMATION OM DIMENSION X(3), FX(3), RO(4), COP1(2), COP2(2), 1 F1(2). F2(2)J EXPO1(50), EXP02(50), RH01(50), RH11(50)p 1 RH02(50), RH12(50), Z3(50), PSI1(1050VIM), PSI2C1050,DIM), 1 COP3(2), ANGLE(21). EXP03(50), 1 ASYO1(50), ASY11(50), ASY02(50)p ASY12(50) INTEGER I,J, Q. ZEE, ITMAX, ITX, ITY VECTOR VALUES VIM=20O,21 VECTOR VALUES DIM=2,0,21 START READ AND PRINT DATA PI=3.14159265 R THE ROOTS OF EQ,(7.8A) ARE OBTAINED BY THE FALSE POSITION R METHOD (SEE REF. (28)) D1=C1*C4-C2*C3 B =-3.*(1.-2,*C4+SIG*SIG-2.*SIG*C1) D =9.*SIG*(SIG-2.*C1-2.*C4*SIG+4.*01) PRINT RESULTS B, D, D1 INTERNAL FUNC ION CHAR.(Z)=Z.P,4+8*Z.P.2+D T=O00000001 Q=1 ITER CHARO=CHAR.(T) WHENEVER CHARO.L.O. THROUGH SERCH1, FnR Y=T,O2,CHAR1.G,.O..ORY.G.40. WHENEVER Y.G.40. PRINT COMMENT S NO SIGN CHANGES TRANSFER TO START END OF CONDITIONAL SERCHi CHARI=CHAR.(Y) OTHERWISE THROUGH SFRCH2, FOR Y=T..O2,CHAR1.L.O..OR.Y.G.40. WHENEVER Y.G.40. PRINT COMMENT s NO SIGN CHANGES TRANSFER TO START END OF CONDITIONAL SERCH2 CHAR1=CHAR.(Y) END OF CONDITIONAL X(1)=Y-.04 XC2)=Y-.02 FX(1)=CHAR.(X(1)) FX(2)=CHAR.(X(2)) PRINT RESULTS X(1), FX(1)p X(2), FX(2) THROUGH LOOPXA, FOR ITX=1,1, ITX.G.ITMAX NUMER=X(2)*FX(1)-X(1)*FX(2) DENOM=FX(1)-FX(2) WHENEVER.ABS. DENOML..ABS.NUMER*EPS2 PRINT COMMENT $ DENOMINATOR TOO SMALLS TRANSFER TO START

121 PROGRAM 3 (Continued) ENO OF CONDITIONAL X ( 3 ) = N U ER / DEN NI1 FX(3)=CHAR. (X(3)) WHENEFVER.A8S,FX(3).LE.EPS1 PRINT COMMENT $ PROCEDURE CONVERGED$ TRANSFER TO ITER1 END OF CONDITIONAL THROUGH LOOPXA, FOR I=l,1, I.G.2 THROUGH LO[PXA, FOR J=I+1l1, J.G.3 WHENEVER.ABS.FX( I).G..AS.FXJ) TENP=FX(I) FX(I)=FX(J) FX(J)=TEMP TEMPX( I) X(I )=X(J) X (J)=TEMP LOOPXA END OF CONDITIONAL PRINT COMMENT $ NO CONIVERGENCES THROUGH EVALb, FOR Yl=l.E-8d, *2, Y1.G.30. CHAR3=CHA.(Y 1) EVAL6 PRINT RESULTS Yit CHAR3 TRANSFER TO START ITER1 RO(Q)=X(3) PRINT RESULTS ITX, X(3), FX(3), Q T=Y-, 02 Q=Q+ I WHENEVER Q.E.3, TRANSFER TO GO TRANSFER TOI ITER GO PRINT RESUILTS RO(1), RO(2) R END OF PROGRAM FOR O9T4INiNG THE ROOTS THROUJGH EVAL 1 FOR I = 1, 1 I1.5.2 R EQS.C(.11) EVALUATED COPl(I)=RO( I)/(3.*SIG) COP2(I)=(3.*SIG*(SIG2.k*Cl)-RO(I)*RO(I))/(6.*SIG*C2) COP3(I)=COP2(I) *RO()/3. R SEE EQS. FOLLOWING EQS.(7.17) F1(I)=.50+COP 1 (T) F2(I)=CUP2(I)/2.+COP3( I) EVAL1 PRINI RESIULTS COPi(T), COP2(I), COP3(I), F1(I)' F2(I) R NEXT THREE CARDS FfrR MILNE CASE ONLY (SEE ES.(17.18A) COP4=-HOL(1)/( 3.*SIG) COPS=COP2(1) COP6=-Rd(1) *ClOP5/3. F 3= 5+CUP4 1M F4=,5*C0P5+ COP6 2M R SEE EQS,(.1idA) BO2=(F4*F1(2)-F3*F2(2))/(Fl(1)*F2(2)-Fl(2)*F2(l)) 3]1 804=(F3*F2(1)-F4*F1l(1))/(F1(1)*F2(2)-F'1(2)*F2(1)) PRINT RESULTS F3, F43 Fa r302, B04 THROUGH EVAL2, FUR I=lt,1 I.G.ZEE

122 PROGRAM 3 (Continued) EXPOl(I)=FXP.(-RO(1)*Z3(I)) EXP02(I)=EXP.(-RO(2)*Z3(I)) R MILNE CARD ONLY EXPO3(I)=EXP.(R(l1)*Z3(I)) R SEE EQS,(7.12) RHOl(I)=EXPO3(I)+PO2*EXPOl(I)+B04*EXP02(I) 5M ASYOl(I)=EXPO3(I)+BO2*EXPOl(I) 6M RH(11I)=COP4*EXPO3(I)+BO2*C'OPl(1)*EXPOl(I)+BO4*COP1(2)*EXPO 7M R MILNE CARD ONLY 1 2(I) ASYl 1(I)=COP4*EXP3( I )+BO2*COP1 )*EXPO1 (I) 8M RH02(I)=CnP5*EXP03(I)+902*COP2( 1 )*EXPOI )+B4*COP2(2)* 9M R MILNE CARD ONLY 1 EXP02(I) ASYO2(I)=COP5*EXPO3(I)+B02*COP2(1)*EXPOl(I) 10M RH12(I)=COP6*EXPU3(I)+BO2*COP3(1)*EXPOl(I)+ B04*COP3(2)* 11M R MILNE CARD ONLY 1 EXP02(I) ASY12(I)=COP6*EXPO3(I)+302*COP3(1)*EXP0 O(I) 12M EVAL2 PRINT RESULTS Z3(I)J EXPO1(I), EXP02(I), RH01(I), RHII1(I), 1 RH02(I), RH12(I), ASYO1(I), ASYll(I), ASY02(I), ASY12(I) R ANGULAR DISTRIBuTION FOLLOWS Q=1 THROUGH EVAL4, FOR L=1.-.1, L.Ls-1. ANGLE(Q)=3bO./(2.,PI)*ARCCOS.(L) EVAL4 Q=Q+1 ANGP=360,/(2.*PI)*ARCCnS,(.975) ANGM=360,/(2.*PI)*ARCCOS.(-.975) THROUGH EVAL7, FOR J=1,lp J.G.ZEE PSP1=.5*RHO1(J)+1.5*,975*RH11(J) PSM1=.5*RHO1(J)l.5*,975*RH1I(J) PSP2=.5*RH02(J)+1,5*.975*NH12(J) PSM2=.5*RH02(J)-1,5*,975*RH12(J) Q=1 THROUGH EVAL3, FOR L=I.,-.1, L.L.-1. PSI1(JQ)=.5*RHO01(J)+15*L*RH11(J) PSI2(J,Q)=.5*RH02(J)+1.5*L*RH12(J) PRINT RESULTS ANGLE(Q), Z3(J), PSII(JQ), PSI2(JQ) EVAL3 Q=Q+1 EVAL7 PRINT RESULTS ANGP, PSP1, PSP2p ANGM, PSM1, PSM2 R NEXT TWO CARDS ARE MILNE CARDS UNLY R SEE EQ.(7.19) EXTRAP=-l./(2.*RO(1) )*ELOG.(-1//B02) PRINT RESULTS EXTRAP TRANSFER TO START END OF PROGRAM $ DATA C1=.59023, C2=.05448P C3=,15963, C4=.44320 Z3(1)=0,p*5,1. 1.52.p t.5,3,,3.5,4.p4.55.,*.8 10, 20., EPSI=1.E-6, EPS2=1.E-20D SIG=1,50951, ZEE=14, ITMAX=50 *

123 PROGRAM 3 (Concluded) R FOR THE CONSTANT SnURCE PROBLEM WE MAKE THE REPLACEMENT OF R THE FOLLOWING CARDS EXACTLY AS WE DID IN THE PREVIOUS CASE R PROGRAM FOR THE SOLUTION OF THE CONSTANT SOURCE PROBLEM OS R BY SPHERICAL HARMONICS-P1 APPROXIMATION OS S1=2.*C2/((1.-2.*C4)*(SIG-2.*Cl)-4,*C2*C3) 1S S2=(SIG-2.*C1)/((1.-2.*C4)*(SIG-2.*Cl)-4.*C2*C3) 25 A02=(S2*F1(2)-SI*F2(2))/(F 1()*F2(2)-Fl(2)*F2(1)) 3S AO4=(S1*F2(1)-S2*F(1))/))/(F1()*F2(2)Fl(2)*F2(l)) 4S RHOl(I)=AO2*EXPO1(I)+AO4*EXP02(I)+2.*S1 5S ASYO1(I)=AO2*EXPO1(I)+2.*S1 6S RH11(I)=COPl(1)*AO2*EXPO1(I)+COP1(2)*AO4*EXPO2(I) 7S ASYll(I)=COP1(1)*AO2*EXPO1(I) 8S RH02(I)=COP2(1)*AO2*EXPOl(I)+COP2(2)*AO4*EXP02(I)+2,*S2 9S ASY02(I)=COP2(1)*AO2*EXPOl(I)+2.*S2 10S RH12(I)=COP3(1)*AO2*EXPOl(I)+COP3(2)*AO4*EXP02(I) 1lS ASY12(I)=COP2(1)*A02*EXPOl(I) 12S

124 PROGRAM 4 S COMPILE MAD, EXECUTE, DUMP, PRINT OBJECT R PROGRAM FOR THE SOLUTION OF THE MILNE PROBLEM BY R SPHERICAL HARMONICS-P3 APPROXIMATION DIMENSION X(3), FX(3), RO(4), COP1(4), COP2(4), COP3(4), 1 COP4(4), COP5(4), COP6(4), COPT(4), F1(4), F2(4), F3(4), 1 F4(4), EXPO1(50), EXP02(50), EXP03(50), EXP04(50)p RH01(50), 1 RHI1(50), RH21(50), RH31C50), RH02C50), RH12(50), RH22(50)p 1 RH32(50), Z3(50), PSI1(1050,VIM), PSI2(1050,DIM), ANGLE(21), I Y(4), A(20,VIN), ASY01(50), ASY02(50), EXPO5(50) INTEGER IJ, Q, ZEE, ITMAX, ITXP ITY VECTOR VALUES VIM=2,0,21 VECTOR VALUES DIM=2O,,21 VECTOR VALUES VIN=2,0,5 START READ AND PRINT DATA Dl=CI*C4-C2*C3 PI=3.14159265 C=81, D=990,*(CA+Cl*SIG)-810,*(1.+SIG*SIG) E=(110.*SIG*Cl-90. *SIG*SIG)*(110 O*C4-90.)+945**(1,'2,*C4+SIG, 1 P.3*(SIG-2.*Cl))-10ll*110,*SIG*C2*C3 F=105*(1'.-2.*C4)*(110.*SIG*C1-90.*SIG*SIG)+105.*SIG.P.3*( 1 SIG-2.*C1)*(110.*C4-90.)+210.*110.*SIG*C2*C3*(1,+SIG*SIG) G=105,*105.*SIG.P.3*(SIG-2,*Cl)*(1,-2.*C4)-210,*210,*SIG,P,3 1 *C2*C3 PRINT RESULTS C, D, E, F, G INTERNAL FUNCTION CHAR,(Z)=C*Z.P.8+D*Z.P.6+E*Z.P.4+F*Z.P.2+G T=,00000001 Q=1 ITER CHARO=CHAR.(T) CHAR1=O WHENEVER CHARO.L.O. THROUGH SERCHl,FOR Yl=T,.OCHARl.GO..OR.Y1.G40. WHENEVER Y1.G.40, PRINT COMMENT S NO SIGN CHANGES TRANSFER TO START END OF CONDITIONAL SERCHi CHARI=CHAR.(Y1) OTHERWISE THROUGH SERCH2,FOR YI=T,,O1,CHAR1L.O..ORY1.G.40. WHENEVER Y1.G.40. PRINT COMMENT S NO SIGN CHANGES TRANSFER TO START END OF CONDITIONAL SERCH2 CHARI=CHAR.(Y1) END OF CONDITIONAL X(1)=Y1-,02 X(2)=Y1-,01 FX(1)=CHAR(X( 1)) FX(2)=CHAR (X(2))

125 PROGRAM 4 (Continued) THROUGH LOOPXA, FOR ITX=,11, ITXG.ITNMAX NUMER=X(?)*FX(l)-X(1)*FX(2) DENOM=FX1()-FX(2) WHENEVER.ABS. DENODML..ABS.NUMER*EPS2 PRINT COMMENT $ DENOMINATOR TOO SMALLS TRANSFER TO START END OF CONDITIONAL X(3)=NUMER/DENOM FX(3)=CHAR.(X(3)) WHENEVER,ABS.FX(3).LE.EPS1 PRINT COMMENT S PROCEODURE CONVERGEDS TRANSFER TO ITER1 OR WHENEVER Q.E.2.AND,ABS.FX(3).L..0003 TRANSFER TO ITER1 OR WHENEVER Q.E.3.AND..ABS.FX(3).L..003 TRANSFER TO ITER1 OR WHENEVER Q.E.4.AND..ABS.FX(3),L.,1 TRANSFER TO ITER1 END OF CONDITIONAL THROUGH LOOPXA, FOR I=1p,1 1.G.2 THROUGH LOOPXA, FOR J=I+l,1 J.G.3 WHENEVER.ABS.FX(I).G..ABS.FX(J) TEMP=FX(I) FX(I)=FX(J) FX(J)=TEMP TEMP=X(I) X(I)=X(J) X(J)=TEMP PRINT RESULTS X(1), FX(1), X(2), FX(2) LOOPXA END OF CONDITIONAL PRINT COMMENT S NO CONIVERGENCE$ THROUGH EVAL6, FOR Y1=O.,.l, Y1,G,30. CHAR3=CHAR.(Y1) EVAL6 PRINT RESULTS Y1, CHAR3 TRANSFER TO START ITERI RO(Q)=X(3) PRINT RESULTS ITX, X(3), FX(3), Q T=Yl-.Ol Q=Q+1 WHENEVER Q.E.5, TRANSFER TO GO TRANSFER TO ITER GO PRINT RESULTS RO(1)...,,RO(4) R END OF PROGRAM FOR CALCULATING THE ROOTS SEE CHAPTER 7 FOR R APPROPRIATE EQUATIONS FOR THE FOLLOWING PROGRAM THROUGH EVALl, FOR IT=,1, I.G,4 SQR=RO(I)*RO( I) SQSG=SIG*SIG COP1(I)=RO(I)*(35,*SQSG-9.*SQR)/(5.*SIG*(21.*SQSG-11.*SQR)) COP2(I)=14.*RO(I)*SIG*COP1(I)/(35**SQSG-9.*SQR) COP3(I)=3.*RO(I)*COP2(I)/(T.*SIG)

126 PROGRAM 4 (Continued) COP4(I)=10.*C3*(21.-11.*SQR)/(SQR*C9.*SQR-35*)+5.*(1.-2,*C4) 1 *(21.-ll11.*SQR)) COP5(I)=RO(I)*(35.-9.*SQR)*COP4(I)/(5.*(21.-11.*SQR)) COP6(I)=14.*RO(I)*COP5(I)/(35.-9,*SQR) COP7(I)=3.*RO(I)*COP6(I)/7. Fl(I)=.5+COP1(I)+5.*COP2(I)/8. F2(I)=.25+3,*COP1(I)/5,+5.*COP2(I)/8.+2.*COP3(I)/5. F3(I)=*5*COP4(I)+COP5(I)+5.*COP6(I)/8. F4(I)=,25*COP4(I)+3.*COP5(I)/5.+5.*COP6(I)/4.+2.*COP7(I)/5. EVAL1 PRINT RESULTS RO(I), COPI(I), COP2(I), COP3(I), COP4(I), 1 COP5(I), COP6(I), COP7(I)p F1(I), F2(I), F3(I), F4(I) R FOR THE CONSTANT SnURCE PROBLEM WE INSERT CARDS iS AND 2S A(1,5)='.5+COP1(1)+5./8.*COP2(1) 3M A(2,5)=-,25+3./5.*COP1(1 )+5./8.*COP2(1)+2./5.*COP3(1) 4M AC3,5)=-,5*COP4(1)+CIP5(1)+5o/8,*COP6(1) 5M A(4,5)=-.25*COP4(1)+3./5.*COP5(1)+5./8.*COP6C1)+2./5.*COP7(1) 6M THROUGH EVAL5, FOR I=il1, I.G.4 A(1I1)=Fl(I) A(2, I)=F2(I) A(3,I)=F3(I) EVAL5 A(4,I)=F4(I) PRINT RESULTS A(1,1)...A(4,5) EPSS=1.E-20 DETER=SIMUL.(4,A,Y,EPSS) PRINT RESULTS A(l,)...A(4, 5), Y(1)..Y(4), DETER THROUGH EVAL2, FOR I=l,l I.G.ZEE EXPOl(I)=EXP.(-RO(1)*Z3(I)) LXPO2(I)=EXP.(-R0(2)*Z3(I)) EXP03(I)=EXP.(-RO(3)*Z3(I)) EXP04(I)=EXP.(-RO(4)*Z3(I)) R MILNE CARD ONLY EXP05(1)=EXP.(R0(1)*Z3(I)) ASYO1(I)=Y(1)*FXPOI(I)+EXPO5(I) 7M RHO1(I)=Y(1)*EXPO1(I)+Y(2)*EXPO2(I)+Y(3)*EXPO3(I)+Y(4)* 8M 1 EXP4(I)+EXPOS5(I) 9M RHII(I)=Y(1)*COP1(1 )*EXPO1(I)+Y(2)*COP1(2)*EXPO2(I)+Y(3)* 10M i COPI(3)*EXPU3(I)+Y(4)*COP1(4)*EXP04(I)-EXP05(I)*COP1(1) 11M RH21(I)=Y(1)*COP2(1)*EXPOl(I)+Y(2)*COP2(2)*EXPO2(I)+Y(3)* 12M 1 COP2(3)*FXPO3(I)+Y(4)*COP2(4)*EXP04(I)-EXP05(I)*COP2(1) 13M RH31(1I)=Y(1)*COP3(1)*EXP1((I)+Y(2)*COP3(2)*EXP2(I)+Y(3)* 14M 1 COP3(3)*EXP03(I)+Y(4)*COP3(4)*EXP04(I)-EXP05(I)*COP3(1) 15M ASY02(1)=Y(1)*COP4(1)*EXP1C(I)+EXPO5( I)*COP4(1) 16M RH02(I)=Y(1)*COP4(1)*EXPO1(I)+Y(2)*COP4(2)*EXPO2(I)+Y(3)* 17M 1 COP4(3)*EXP03(I)+Y(4)*COP4(4)*EXP04(I)+EXPO5(I)*COP4(1) 18M RH12(I)=Y(1)*COP5(1)*EXP01(I)+Y(2)*COP5(2)*EXP02(I)+Y(3)* 19M 1 COP5(3)*EXP03(I)+Y(4)*COP5(4)*EXP04(I)-EXP05(I)*COP5(1) 20M RH22(I)=Y(1)*COP6(1)*EXPOl(I)+Y(2)*COP6( 2)*EXP02(I)+Y(3)* 21M 1 COP6(3)*EXP03(I)+Y(4)*COP6(4)*EXP04(I)-EXP05(I)*COP6(1) 22M RH32(I)=Y(1)*CDPT7(1)*EXPO(I)+Y(2)*COP7(2)*EXPO2(I)+Y(3)* 23M 1 COP7(3)*EXP3 )Y4*CP 4)*EXP4I)-EXP05(I)*COPT(1) 24M

127 PROGRAM 4 (Continued) EVAL2 PRINT RESULTS Z3(I), RHO1(I)p RHll(I), RH02(I)t RH12(I), 1 ASYO1(I), ASY02(I) R NEXT TvO CARDS ARE MILNE CARDS ONLY EXTRAP=-.5/RO(1)*ELOG.(-1./Y(l)) PRINT RESULTS EXTRAP R ANGULAR OISTRIBUTION FOLLOWS Q=1 THROUGH EVAL4P FOR L=l.,-.l, L.L.-1. ANGLE(Q)=360,/(2.*PI)*ARCCOS.(L) EVAL4 Q=Q+I THROUGH EVAL3P FOR J=l,lt J.G,ZEE Q=1 THROUGH EVAL3, FOR L=l1.-.l1 L.L.-1. PSIl(J,Q)=.5*RHOl(J)+1.5*L*RHll(J)+1.25*(3.*L*L-1.)*RH21(J) 1 +7,/4.*(5.*L.P.3-3.*L)*RH31(J) PST2(JQ)=.5*RH02(J)+1,5*L*RH12(J)+1.25*(3.*L*L-1.)*tH22(J) 1 +7,/4.*(5.*L.P.3-3.*L)*RH32(J) PRINT RESULTS ANGLE(Q), Z3(J), PSI1(JQ), PSI2(JAQ) EVAL3 Q=Q+1 TRANSFER TO START END OF PROGRAM S COMPILE MAD, EXECUTE, DUMP, PRINT OBJECT R THIS SUBROUTINE SOLVES THE SET OF SIMULTANEOUS EQUATIONS R SEE EQS,(7.30) AND (7.31) FOR DETAILS SEE REF. (18) EXTERNAL FUNCTION STMUL.(NA,XPEPS) DIMENSION IR(50), JC(50), CJ(50) NORMAL MODE IS INTEGER FLnATING POINT A, BIGA, XP DETERP EPS, AJCK NP = N+ 1 UETER=1. THROUGH L1, FOR K=1l,1 K.G.N BIGA=O. THROUGH L2J FOR I-1l1, I.G.N THROUGH L2P FOR J=llp, JG.N THROUGH L3P FOR I1=lp1, I1.E.K THROUGH L3, FOR J1=1,1, J1.E.K L3 WHENEVER I.E.IR(Il).OR.J,E.JC(J1). TRANSFER TO L2 WHENEVER.ABS.A(I,J).G.3IGA BIGA=.ABS.A(I,J) IR(K)=I JC(K)=J L2 END OF CONDITIONAL WHENEVER HIGA.L.EPSP FUNCTION RETURN 0 BIGA=A(IR(K),JC(K)) DETER = DETER*BIGA THROUGH L4o FOR J=,11, J.G.NP1 L4 A(IR(K),J)=A(IR(K),J)/BIGA THROUGH L1I FOR I=ll, I,G.N WHENEVER I.NE.IR(K) AJCK=A(I,JC(K))

128 PROGRAM 4 (Concluded) THROUGH L6, FOR J= l1, J.G.NP1 L6 A(I,J)=A(IJ)-AJCK*A(IR(K),J) L1 END OF CONDITIONAL THROUGH L7T FOR I=1,1, I.G.N CJ(IR(I) )JC(I) L7 X(JC(I))=A(IR(I),NP1) COUNT=O THROUGH L9, FOR I=1,1, I.E.N THROUGH LQ* FOR J=I+1,1, J.G.N WHENEVER CJ(J).L.CJ(I) TEMP=CJ(J) CJ(J)=CJ(I) CJ( I )=TEMP COUNT = cnUNT+l L9 END OF CONDITIONAL WHENEVER COUNT/2*2.NE.CUUNT, DETER=-DETER FUNCTION RETURN DETER END OF FUNCTION S DATA C1=.59023, C2=.05448, C3=.15963, C4=.44320, Z3(1)=0..5,1.,1.5,2., 2.5,3.,3.5S,4.4.5,5.,8.,10.,20., EPS1=I.E-6, EPS2=1.E-20, SIG=1.50951, RO(1)=.1376522, 1,162925, 2.068132, 3.289316, ZEE=14, ITMAX=50 * R FOR THE CONSTANT SnURCE PROBLEM WE REFER THE READER TO THE R PREVIOUS INSTRUCTIONS FOR THE P1 CASE R PROGRAM FOR THE SOLUTION OF THE CONSTANT SOURCE PROBLEM OS R BY SPHERICAL HARMONICS-P3 APPRUOXIMATION OS S1=2.*C2/((1.-2,*C4)*(SIG-2,*Cl)-4.*C2*C3) iS S2=(SIG-2.*Cl)/((1.-2.*C4)*(SIG-2.*Cl)-4.*C2*C3) 2S A(1,5)='S1 3S A(2,5)=-S1/2. 4 ASYO1(1)=Y(1)*EXPO1(I)+2)*S1 7S RHO1(I)=Y(l)*EXPO(II)+Y(2)*EXPO2(I)+Y(3)*EXPO3(I)+Y(4)* 8S 1 EXP04(I)+2.*S1 9S RH11(I)=Y(1)*COPl(1)*EXPO1(I)+Y(2)*C1OP(2)*EXPO2(I)+Y(3)* 10S 1 COPl(3)*EXP03(I)+Y(4)*COP1(4)*EXP04(I) 11S RH21(I)=Y(L)*COP2(1)*EXPO1(I)+Y(2)*COP2(2)*EXPO2(I)+YC3)* 12S 1 COP2(3)*EXPU3(I)+Y(4)*COP2(4)*EXP04(I) 13S RH31(I)=Y(1)*CUP3(1)*EXP01(I)+Y(2)*COP3(2)*EXPO2(I)+Y(3)* 14S i COP3(3)*FXP03(I)+Y(4)*COP3(4)*EXP04(I) 15S ASYO2(I)=Y(1)*COP4(1)*EXPO1(I)+2.*S2 16S RH02(I)=Y(1)*COP4(1)*EXPOl(I)+Y(2)*COP4(2)*EXP02(I)+Y(3)* 17S 1 COP4(3)*EXPO3(I)+Y(4)*COP4(4)*EXP04(I)+2.*S2 18S RH12(I)=Y(1)*COP5(1)*EXPO1(I)+Y(2)*cOP5(2)*EXPO2(I))Y(3)* 19S 1 COP5(3)*EXP03(I)+Y(4)*C0P5(4)*EXP04(I) 20S RH22(I)=Y(1)*COP6(1)*EXPOi(I)+Y(2)*COP6(2)*EXP02(I)+Y(3)* 21S 1 COP6(3)*EXP03(I)+Y(4)*COP6(4)*EXPO4(I) 22S RH32(I)=Y(1)*COP7(1)*EXPOi(I)+Y(2)*COP7C(2)*EXPO2()+Y(3)* 23S 1 COP7(3)*EXP03(I)+Y(4)*COP7(4)*EXPO4(I) 24S

129 PROGRAM 5 $ COMPILE MAD. EXECUTE, OUMP, PRINT OBJECT R PROGRAM FOR THE SOLUTION OF THE CONSTANT SOURCE PRObLEM BY R SPHERICAL HARMONICS-DP1 APPROXIMATION DIMENSION X(3), FX(3), RO(4), COPI(4), COP2(4), COP3(4), 1 COP4(4), COP5(4), COP6(4), COP7(4), F1(4), F2(4), F3(4), 1 F4(4), EXP01(50), FXP02(50), EXP03(50), EXP04(50), RHO1(50), 1 RH11(50), RH21(50), RH31(50), RH02(50), RH12(50), RH22(50), 1 RH32(50), Z3(50), PSI(1O50,VIM), PSI2(1050,OIM), ANGLE(21), 1 Y(4), A(20.,VIN), EXPO5C50), ASYU1(50), ASY02(50) INTEGER I,J, Q. ZEE, IT'AX, ITX, ITY VECTOR VALUES VIM=2O, 21 VECTOR VALUES DIM=2,0,21 VECTOR VALUES VIN=2,0,5 START READ AND PRINT DAIA PI=3,14159265 D1=C1 *C4-C2*C3 C=-1 D=24*(SIG*SIG-CI*SIG+i.-C4) E=576.*SIG*(C2*C3+(i.,-C4)*(C1-SIG))-36.*(SIG.P.3*(SIG-2,*C1)+ 1 1.-2.*C4) F=-864.*SIG*C-S*((C-SIG)*(I.-2.*C4)(1.-C4)*SIG*SIG*(SIG-2*C)) 1 -1728.*C2*C3*SIG*(IGSIG*SI+1) L=-1296.*SIG.P.3*(SIG-2.*Cl)*(l.-2.*C4)+5184.*C2*C3*SIGP.3 PRINT RESULTS Co D, E, F, G INTERNAL FUNCTION CHAR.(Z)=C*Z.P.8D+D*Z.P.6+E*ZP.4+F*ZP.2+G T=.00000001 Q=1 ITER CHARO=CHAR.(T) CH4RI=O0 WHENEVER CHARO.L.O, THROUGH SERCH1,FOR Y=T,.01,CHAR1 G,..R.Y1 G.40. WHENEVER Yl.G.40, PRINT COMMENT $ NO SIGN CHANGES TRANSFER TO START END OF CONDITIONAL SERCH1 CHAR1=CHAR.(Yl) OTHER I SE THROUGH SERCH2,FOR Yl=T,.OlCHAR.L.,..OR.Y1.G.40 WHENEVER Y1.G.40, PRINT COMMENT $ NO SIGN CHANGES TRANSFER TO START END OF CONDITIONAL SERCH2 CHARI.=CHAR,(Y1) END OF CONDITIONAL X(1)=Yl-.02 X(2)=Y-,01Ol FX(1)=CHAR.(X(1)) FX(2)=CHAR.(X(2)) THROUGH LOOPXA, FOR ITX=1,1' ITX.G.ITMAX

150 PROGRAM 5 (Continued) NUMER=X(2)*FX(1)-X C1)*FX(2) DENOM=FX( 1 )-F X(2) WHENEVER.ABS. DENOM.L...A8S.NUIER*EPS2 PRINT COMMENT S I)ENOMINATOR TOO SMALLS TRANSFER TO START END OF CONDITIONAL X(3)=NUMER/DENOM FX(3)=CHAR (X(3)) WHENEVER.ABS.FX(3),LE.EPS1 PRINT COMMENT $ PROCEDURE CONVERGEDS TRANSFER TO ITER1 OR WHENEVER.E.3.AND..ABs.FX(3),L..005 TRANSFER TO ITER1 OR WHENEVER Q.E.4.AND..AS.FX(3) L..01 TRANSFER TO ITER1 END OF CONDITIONAL THROUGH LOOPXA, FnR I=1,1, I.G.2 THROUGH LOOPXAP FOR J=I+l1lp JG.3 WHENEVER.ABS.FX(I),G..A9S.FX(J) TEMP=FX(I) FX(I)=FX(J) FX(J)=TEMP TEMP=X ( I ) X(I)=X(J) X(J)=TEMP PRINT RESULTS X(1), FX(1), X(2), FX(2) WHENEVER FX(1).E.FX(2) FX(3)=FX(1) X(3)=X(1) TRANSFER TO ITER1 END OF CONDITIONAL LOOPXA END OF CONDITIONAL PRINT COMMENT S NO CONVERGENCE$ THROUGH EVAL6, FOR Yl=O0,.1, Y1,G.30. CHAR3=CHAR.(Y1) EVAL6 PRINT RESULTS Y1, CHAR3 TRANSFER TO START ITER1 RO(Q)=X(3) PRINT RESULTS ITX, X(3), FX(3), Q T=Yl-,Ol WHENEVER Q.E.5, TRANSFER TO GO TRANSFER TO ITER GO PRINT RESULTS RO(1).,.RO(4) S1=2.*C2/((1.-2,*C4)*(SIG-2,*Cl)-4.*C2*C3) S2=(SIG-2.*C1)/((1.-2.*C4)*(SIG-2.*C1)-4.*C2*C3) THROUGH EVAL1, FOR I=1t1, I.G.4 SQR=RO(I)*RO(I) SQSG=SIb*SIG COP1(I)=RO(I)*(12.*SQSG-SQR)/(12*,SIG*(3.*SQSG-SQR))

131 PROGRAM 5 (Continued) COP2(I)=9:*SIG*RO(I)*COPI(I)/(2.*(12.*SQSG-SQR)) COP3(I)=RO(I)*COP2(I)/(6.*SIG) COP4(I)=24.*C3*(SQR-3.)/(SQR*(12,-SQR)+12,*(1.-2,*C4)*(SQR1 3.)) COP5(I)=RO(I)*(SQR-12.)*COP4(I)/(12.*CSQR-3.)) COP6(I)=9.*RU(I)*COP5( )/(2.*(12,-SQR)) COP7(I)=RO(I)*C0P6(I)/6. Fl(I)=.5+,75*COP1(I)-COP3(I) F2(I)=.25*COP1(I)+2./3.*CUP2(I)+COP3(I) F3(I)=.5*COP4(I)+.75*COP5(I)-CUP7(I) F4(I)=.25*COP5(I)+2./3.*COP6(I)+COP7(I) EVAL1 PRINT RESULTS RO(I), COPI(I), COP2(I), COP3(I)p COP4(I), 1 COP5S(I), COP6(I), COP7(I), FI(I), F2(I), F3(I), F4(I) At1 5)=-S1 A(2,5)=0. A(3,5)=-S2 A(4,5)=0, THROUGH EVALS5 FOR I=1,l, I.G,4 A( 1 I)=FI(I) A(2, I )=F2(I) A(3 I)=F3( I) EVAL5 A(4,I)=F4(I) PRINT RESULTS A(1,1)...A(4,5) EPSS=1.E-20 DETER=SIMUL.(4,A,Y,YPSS) PRINT RESULTS A(1,1)...A(4,5), Y(1)...Y(4), DETER THROUGH EVAL2, FOR I=l,1 I.G.ZEE EXPOl(I)=EXP.(-RO(1)*Z3(I)) EXPO2(I)=EXP.(-RO(2)*Z3(I)) EXP03(I)=EXP.(-RO(3)*Z3(I)) EXP04(I)=EXP.('RO(4)*Z3(I)) RHO I )=Y(1)*EXPO I )+Y(2)*EXP2 I )+Y3)*EXP03( I )+Y(4)* 1 EXPO4(I)+2.*S1 ASYOl(I)=Y(1)*EXP0(I)+2.*Sl RH11(I)=Y(1)*COPl(I)*EXPOl(I)+Y(2)*COP1(2)*EXPO2(I)+Y(3)* 1 COPI(3)*EXP03(I)+Y(4)*COP1((4)*EXP04(I) RH21(I)=Y(1)*COP2C1)*EXP01(I)+Y(2)*COP2(2)*EXP02(I)+Y(3)* i COP2(3)*EXPU3(I)+Y(4)*COP2(4)*EXP04(I) RH31(I)=Y(1)*COP3(1)*EXPO1(I)+Y(2)*COP3(2)*EXPO2(I)+Y(3)* 1 COP3(3)*EXP03(I)+Y(4)*COP3(4)*EXP04(I) RHO2(I)=Y(1)*COP4(1)*EXPOl(I)+Y(2)*COP4(2)*EXPO2(I)+Y(3)* 1 COP4(3)*EXPU3(I)+Y(4)*COP4(4)*EXP04(I)+2,*S2 ASY02(I)=Y(1)*COP4(1)*EXPO1(I)+2.*S2 RH12(I)=Y(1)*COP5(1)*EXP01(I)+Y(2)*COPS(2)*EXPO2(I)+Y(3)* 1 COP5(3)*EXP03(I)+Y(4)*COP5(4)*EXP04(I) RH22(I)=Y(1)*COP6(1)*EXPOl(I)+Y(2)*COP6(2)*EXP02(I)+Y(3)* 1 COP6(3)*EXP03(I)+Y(4)*COP6(4)*EXP04(I) RH32(I)=Y(1)*COP7(1)*EXPO1(I)+Y(2)*COP7( 2)*EXPO2(I)+Y(3)* 1 COP7(3)*EXPU3( )+Y(4)*COP7(4)*EXP04(1) EVAL2 PRINT RESULTS Z3(I), RH01(I), RH11(I), RH02(I)p RH12(I),

132 PROGRAM 5 (Continued) 1 ASY1I(I), ASY02(I) R ANGULAR DISTRIBUTION FOLLOWS Q=i THROUGH EVAL4P FOR L:=,,-.1o L.LI.-1 ANGLE(Q)=360./(2.*PI)*ARCCOS.(L) EVAL4 Q=Q+1 THROUGH EVAL3, FOR J=r1,l JG9ZEE Q=i THROUGH EVAL3P FOR L=1l.-.1, L.L.-I, WHENEVER L.G..001 PSI1(JQ)=,5*RHO1(J)-2.*RH21(J)-4**RH31(J)+L*(4,*RH21(J)+ 1 3./2.*RHll (J)+6.*RH31(J)) PSI2(JQ)=.5*RHO2(J)-2,*RH22(J)-4.*RH32(J)+L*(4.*RH22(J)+ 1 3./2.*RH12(J)+6.*RH32(J)) OR WHENEVER L.L.-.OO11 PSI1(JQ)=.5*RHO1(J)-2.*RH21(J)+4.*RH31(J)+L*(-4.*RH21(J)+ 1 3./2.*RH11(J)+6,*RH31(J)) PS12(JQ)=.5*RHO2(J)-2,*RH22(J)+4**RH32(J)+L*(-4**RH22(J)+ 1 3./2.*RH12(J)+6,*RH32(J)) R HERE WE CALCULATE THE DISCONTINUITY IN THE ANGULAR R DISTRIBUTION AT MU=O, OR WHENEVER Q.E.11 DISPI =.5*RHO1(J)'2.*RH21(J)'4.*RH31(J)+L*(4.*RH21(J)+ 1 3,/2,*RH11(J)+6.*RH31(J)) DISP2 =,5*RHO2(J)'2.*RH22(J)4,.*RH32(J)+L*(4,*RH22(J)+ 1 3./2.*RH12(J)+6.*RH32(J)) DISMi =.5*RHO1(J)-2.*RH21(J)+4.*RH31(J)+L*('4.*RH21(J)+ 1 3./2.*RH11(J)+6,*RH31(J)) DISM2 =,5*RHO2(J)-2.*RH22(J)+4.*RH32(J)+L*(-4.*RH22(J)+ 1 3./2.*RH12(J)+6,*RH32(J)) PRINT RESULTS DISP1i DISP2P DISMI1 DISM2, ANGLE(Q) END OF CONDITIONAL PRINT RESULTS ANGLE(Q), Z3(J), PSI1(J,Q)p PSI2(JQ) EVAL3 Q=Q+l TRANSFER TO START END OF PROGRAM S COMPILE MAD, EXECUTE, DUMP, PRINT OBJECT EXTERNAL FUNCTION SIMULI(NA,AXEPS) DIMENSION IR(50)p JC(50), CJ(50) NORMAL MODE IS INTEGER FLOATING POINT A, BIGA, X, DETER, EPS, AJCK NPI=N+I DETER=1. THROUGH LI1 FOR K=1)1, K.G.N BIGA=O. THROUGH L2, FOR I=1 1, I.G.N THROUGH L2J FOR J=1l,1 J.G.N THROUGH L3, FOR I1=1,1, I1.E.K THROUGH L3, FOR J1=1,1, J1.E.K L3 WHENEVER I.E.IR(I1).ORJ.E.*JC(J)(J TRANSFER TO L2

133 PROGRAM 5 (Concluded) WHENEVER.ABS.A(IJ).G.BIGA BIGA=,ABSA(IJ) IR(K)=I JC(K)=J L2 END OF CONDITIONAL WHENEVER BIGA*L.EPS, FUNCTION RETURN 0 BIGA=A(IR(K),JC(K)) DETER = DETER*BIGA THROUGH L4J FOR J=l,1, J,G.NP1'.4 A(IR(K),J)=A(IR(K),J)/BIGA THROUGH L1, FOR I=1'1, I.GN WHENEVER I.NE.IR(K) AJCK=A( IJC(K)) THROUGH L6P FOR J=1,l, JG.NP1 L6 A(IJ)=A(IJ)-AJCK*A(IR(K)PJ) L1 END OF CONDITIONAL THROUGH L7i FOR I=1,1, IoG.N CJ(IR(I))=JC(I) L7 X(JC(I))=A(IR(I),NP1) COUNT=O THROUGH L9' FOR I=1,1, IE.N THROUGH L9, FOR J=I+1,1, J.G.N WHENEVER CJ(J).L.CJ(I) TEMP=CJ(J) CJ(J)=CJ(I) CJ(I)=TEHP COUNT = COUNT+1 L9 END OF CONDITIONAL WHENEVER COUNT/2*2,NE,COUNT, DETER=-DETER FUNCTION RETURN DETER END OF FUNCTION S DATA C1=.59023, C2=,05448, C3=.15963, C4=.44320, Z3(1)=0.,.5,1.,15,2., 2.5,3,3.5,4,4.5,5.,8.,10.,~0., EPS1=1:E'6# EPS2=1:E'20P SIG=1,50951 ZEE=14, ITMAX=50 * R FOR THE SOLUTION TO THE MILNE PROBLEM IN THIS R APPROXIMATION WE REFER THE READER TO THE PREVIOUS P3 PROGRAM R FOR THE APPROPRIATE CHANGES, ALSO SEE CHAPTER 7

UNIVERSITY OF MICHIGAN III1 i'1 1ll IIIII11111111111111111111111111111 3 9015 03483 4047