THE U N I V E R S I T Y OF M I C H I G A. N COLLEGE OF ENGINEERING Department of Electrical Engineering Space Physics Research Laboratory Final Report STUDIES TOWARD THE DEVELOPMENT OF A. D-REGION PROBE J. C. Pearl ORA. Project 05235 under contract with: DEPARTMENT OF THE ARMY DETROIT ORDNANCE DISTRICT CONTRA.CT NO. DA.-20-018-ORD-24992 DETROIT, MICHIGAN administered through: OFFICE OF RESEARCH A.DMINISTRATION ANN ARBOR Apri 1 1965

ACKNOWLEDGMENTS The author wishes to express his appreciation to Professor A. F. Nagy for his constant and willing advice and guidance throughout the period during which this work was performed. Further thanks go to Gary Poole and James Nichols for carrying out the extensive programming which the problem necessitated. For the design and construction of the electronics for the prototype probe assembly thanks are due to Tuck Bin Lee; the cooperation of the staff of the BRL Wind Tunnel during the actual testing of this probe is thankfully acknowledged. iii

TABLE OF CONTENTS Page LIST OF FIGURES vii NOMENC LATURE ix CHAPTER 1. INTRODUCTION 1 2. SELECTION OF THE PROBE CONFIGURATION 2 3. WIND TUNNEL TESTS 4 4. MATHEMATICAL INVESTIGATION OF THE CURRENT COLLECTING CHARACTERISTICS OF A MOVING PROBE 8 4.1 The Physical Problem 8 4.2 Mathematical Formulation 9 4.3 Discussion of Problem 10 5. THE "FLOW LINE" METHOD 12 5.1 Outline of Method 12 5.2 Determination of Flow Lines-Phase I 13 5.3 Determination of Charge Densities 16 5.4 Calculation of Potential 18 6. INITIAL APPROXIMATIONS USED IN PROGRAMMING THE FLOW LINE METHOD 26 6.1 Introduction 26 6.2 Initial Flow Lines 26 6.3 Charge Distribution Resulting from Starting Flow Lines 28 6.4 Termination of Flow Lines 28 7. COMPUTER PROGRAM AND NUMERICAL RESULTS 44 7.1 Computer Program 44 7.2 Presentation of Sample Run 45 8. EXTENSION AND CONCLUSIONS 51 8.1 Extension to Phase II 51 8.2 Further Extensions and Conclusions 53 BIBI;OGRAPHY 55 APPENDIX. FLOW DIAGRAM AND PROGRAM LISTING 57

LIST OF FIGURES Figure Page 1. Configuration of probe assembly. 3 2. Probe installed for wind tunnel test. 5 3. Measured voltage-current characteristic at M-3.02 for low range of applied voltages. 6 4. Measured voltage-current characteristic at M-3.02 for high range of applied voltages. 6 5. Schlieren photograph of probe at M=3,02. 7 6. Representation of the physical environment in front of a supersonic probe. 8 7. The difference scheme for potential gradient calculations. 14 8. Parametrization of flow lines. 17 9. Variables used in potential integration. 19 10. Variables relating charge with its image charge. 20 11. Section of charge toroid, 23 12. Construction for integration when field point lies within charge toroido 24 13. Flow past sphere, ( I1 < a). 27 14. Flow past sphere, ( ~I i1 > a). 27 15. Curve of t vs, (R=l). 36 16. Curve of t v's 7, ( =10). 37 17. Curve of t vs y, (r=100). 38 18. Curve of f(x) vs x. 39 vii

LIST OF FIGURES (Concluded) Figure Page 19. Curve of g(t) vs t, (y=O). 40 20. Curve of g(t) vs t, (y=0). 41 21. Curve of g(t) vs t, (y=O). 42 22. Curve of h(y,t) v' t, (y#O). 43 23. a flow lines (initial approximation). 46 24. C flow lines (initial approximation). 47 25. a flow lines (first iteration). 48 26. C flow lines (first iteration). 49 27. Flow lines crossing shock. 51 28. The three points of tabulated velocity field which are 52 nearest to point P.

NOMENCLATURE Alphabetic A surface area Ai see Eqn. (5.4-.10) a probe radius (= a ) Bi see Eqn. (5.4-10) D coefficient of diffusion e electronic charge f function representing ro (Eqn. (5.3-5)); function used in sheath cutoff (Eqn. (6.4-24)) g'function used in sheath cutoff (Eqn. (6.4-29)) H parameter for finite difference calculations (Eqn. (5.2-8)) h function used in sheath cutoff (Eqn. (6.4-30)) I current J current density (positive species) K mobility constant ( = e ) kT k index for solutions of cubic equation (Eqn. (6.4-14)); Boltzmann constant L current density (negative species) displacement n number of incremented steps in a calculation; number of particles in a Debye cube ix

NOMENCLATURE (Continued) Alphabetic P arbitrary field point p parameter used in flow line computation (Eqn. (5.2-12)) Q charge; point on central axis of torus (Fig. 11) q charge R distance from coordinate origin R1,R2,R(P) distances identified in Fig. 9. r radial cylindrical coordinate; minor radius of torus (Eqn. (5.4-18)) S point in charge density field (Fig. 11) s distance normal to flow lines T point on torus (Fig. 11); temperature t time ( - tV' ); variable used in sheath cutoff (Eqn. (6.4-27)) V neutral gas flow velocity ( = X V ) D v composite velocity term (Eqn. (5.2-2)) x variable used in sheath cutoff (Eqno (6.4-23)) y see Eqn. (6.4-6) z axial cylindrical coordinate ~ti dimensionless ratio (Eqn, (5o4-13))

NOMENCLATURE (Continued) Alphabetic dimensionless parameter (Eqn. (6.2-5)) 7 parameter used in sheath cutoff (Eqn. (6.4-27)) A parameter in difference formula (Fig. 7) 6 indicator of finite increment dielectric constant; error due to approximations negative charge density ( = ~' ) parameter used in sheath cutoff (Eqn. (6.4-33)) complete elliptic integral of the first kind Debye length ( =,/ ) see Eqn. (6.4-.12) p charge density (Eqn. (5.4-9))' t,7" see Eqn. (5.4-26) a positive charge density ( = ) O00 T vo lume electrical potential ( - et kT azimuthal cylindrical coordinate Subscript s i counting or identifying index n evaluation after n increments; summation index xi

NOMENCLATURE (Concluded) Subscripts o identifies flow line (see Section 5.3, p. 16). evaluation for sphere in torus routine (Eqn. (5.4-24).) p condition at probe r radial component z axial component 0o condition at infinity Superscripts and Miscellaneous quantity before being nondimensionalized; image (Eqn. (5.4-5))X evaluation for flow line tangent to probe (Eqn. 6.4-3)) "1 evaluation for flow line target to probe (see p.30) ^ unit vector V gradient operator ( - v') l al absolute value of a [a] greatest integer in a. xii

1. INTRODUCTION Since 1946, when University of Michigan investigators first proposed such experiments, Langmuir probes have been widely used to study the E and F regions of the ionosphere. The theory of operation of Langmuir probes, for these regions, where the mean free path of the particles is large, is well understood and therefore their volt-ampere characteristics can be accurately interpreted in terms of ambient charged particle parameters. The growing interest in the physical properties of the lower parts of the ionosphere has provided impetus to study the feasibility of the use of such current collecting probes for the D region. Recently a number of Langmuir probe measurements have been made at these low altitudes. It is, however, difficult to interpret the data from such measurements, since to date no relevant theoretical work has been published dealing with the current collection characteristics of such probes in regions of high neutral gas density. This report deals with the development of a probe for D region charged particle density measurements and the work done to enable the interpretation of the current measured by such a probe in terms of the ambient density. A test model of the probe was constructed and placed in the wind tunnel of the Aberdeen Proving Grounds. The results of these tests are also discussed.

2. SELECTION OF THE PROBE CONFIGURATION A probe configuration enabling measurements to be made in the E as well as the D region was the aim of this study, since the sounding rockets which are likely to carry these experiments will pass through both of these regions. Because there is a great deal of confidence in E region Langmuir probe measurements, these results could also be used to provide some degree of calibration for the interpretation of the D region data. First the question of the probe position received a great deal of consideration. The velocity of the vehicle is expected to be supersonic, resulting in the formation of a bow shock wave. As even the neutral gas flow well behind the shock is not understood it is considered imperative that the probe be located on the tip of the nose cone or on a boom of sufficient length to place it outside the shock cone of the vehicle. In either case the probe itself will create its own shock disturbance. A nose cone probe was chosen since the technical problems involved in its construction are considerably less than those for a boom supported probe. Because the tip of the nose cone is preferred for a number of other experiments it is felt that if the results from the nose cone probe experiments are encouraging work should be initiated towards the design of boom mounted probes. Having established the nose tip as the preferred position of the probe, a configuration using a hemispherical tip was chosen. The main advantages of such a probe over a conical one, are: (a) the collector area projected in the direction perpendicular to the velocity vector is not very sensitive to small changes in angle of attack as would be the case for example with the usual conical nose tip; (b) the current collection theory for spherical collector has been worked out for the long mean free path operation to be encountered in the E region. The principal disadvantage of the spherical over the sharp conical arrangement is the increased aerodynamic drag, but this increase is not expected to be significant for most applications. For both D and E region application it is considered important to minimize the effects of the distortion in the sheath about the probe, caused by the rest of the nose cone. All these considerations led to the choice of probe configuration indicated in Fig. 1.

4f COLLECTOR INSULATOR -- A' - GUARD INSULATOR 4-I 4 BASE Fig. 1. Configuration of probe assembly.

3. WIND TUNNEL TESTS A decision was made at the beginning of the contract period to observe the behavior of the probe in the wind tunnel at the Aberdeen Proving Grounds. The relatively high operating pressure of this wind tunnel (total pressure 100 torr) clearly does not permit the simultion of conditions typical of the D region but still tests were carried out in the hope of obtaining certain information about the behavior of the probe (e.g., saturation, etc.). A radioactive source was employed to achieve partial ionization in the tunnel, but since there was no way of knowing the degree of ionization, the percentage of electrons and negative ions and what the recombination rates were along the flow, no real interpretation of the experimental results was possible. A photograph of the probe constructed for these tests shown in Fig. 2. The electrometer amplifier used in this test fixture was capable of measuring currents from 10-13 to 10-7 amperes in seven linear ranges. A remote controlled voltage stepping circuit allowed the application of eleven different voltages from -15 to +15 volts to the collector and guard. Current versus voltage readings were taken for a number of Mach numbers and source strength values. Figure 3 shows the current versus voltage characteristics for M=3.02. These data indicate a strong relation between the current collected and the radioactive source strength. These tests gave no indication of current saturation, so the test setup was changed to permit the application of potentials up to 200V. The results of such a test are shown in Fig. 4; even here there is no clear indication of saturation. The current detector and the voltage source for these runs had to be situated outside the tunnel and connected to the probe via shielded cable and therefore considerable noise voltage may have been present. This noise however is not believed to have been of sufficient magnitude to mask the real results. Numerous Schlieren photographs were taken during the test (Fig. 5). These show that the shock separation as expected, was very small. It was possible to vary the angle of attack in the wind tunnel from 0~ to 15~. Changes in the angle of attack over this range resulted in no detectable change in the current collected. However, since these tests were also run with the electronics outside the tunnel the results have to be treated with caution.

:2i(::::i::::::::: - -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~x Fig 2. Probe installed for wind tunnel test.:::;i:i:::::i::i:::i:.::i:: i::::~:~~ 8~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~::: —:::'::::::-~~:::i'i:i-:::::::::iii:::::::-:::::::i:i.i;~~~j::::i::::ii ~~ii~l:i:l:i::~::ilii:iliii~:i:.:::ii i:::i i.:: i~iiii:l::::::i:i,::::::::-~::::::::::: i I::::i:::i | - — ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~::.:_:::::jji:i: U ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~: -X-,- He ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~i~iiii:'l:lil~a~~~:-::::~::i.:j'~~~ji"" ~~I:-::::: $0QE;0{S 20'00?2j00g0000'B' ^>9>>'>XS: HE F:U: A:~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~::::::::~'':~':~::~~ ~TI BBfwBBg:L- ";~-'~i::-_iri, - L e)S~~~~~~~~~~~~~v~~~~ta >S'S~~~~~~~~~~~~~~j~;)'t;>;; i S E~~~~~~~~~~~~~~~~~~~~~i~iiES~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~iiiii:~:i:: Fig. 2. Probe installed for wind tunnel test.~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~i~iii'iini

* ~* * - Voltage (Volts) 0 -15 -10 -5 0 O0. 3 Curie * _-10 -20 E E M = 3.02 -30 - Po=80 cm c -40 u 1.0 Curie * -50 Fig. 3- Measured voltage-current characteristic at M=3.02 for low range of applied voltages. rl0 Source - I Currie M = 3.02 Po = 80 cm CURRENT (10-l AMP) 5 -200 -100 I 100 200 VOLTAGE (VOLTS) -5 0 Fig. 4. Measured voltage-current characteristic at M=3.02 for high range of applied voltages.

- -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:::::: 1~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-~::j: Fig.~~~~~:::::::: 5* Shlieen potorap fpoea =.2 7 iliiiii~ii~i

4. MATHEMATICAL INVESTIGATION OF THE CURRENT COLLECTING CHARACTERISTICS OF A MOVING PROBE 4.1 THE PHYSICAL PROBLEM To formulate a theoretical treatment for the-collection characteristics of the probe discussed in Section 3, we consider a rounded electrode moving with supersonic velocity through a dense neutral gas containing a low density plasma. The principal features of the physical environment near such a body are indicated in Fig. 6. / Neutral flow line I z Sh'eath Detached "Boundary edge shock layer edge Fig. 6. Representation of the physical environment in front of a supersonic probe. Here the flow of the nonionized component of the gas is uniform outside the shock, changes to a region of compressible flow behind this, and passes over into the boundary layer near the probe surface. Significant electrical effects on the ionized components of the flow are restricted to a region around the probe designated the sheath. The scale of this region of electrical influence is determined mainly by the temperatures and densities of the charged particle components and by the potential of the probe relative to the plasma; for sufficiently low densities and/or high probe potentials, this region may extend well beyond the shock. Consequently, within the sheath, strong ionic currents may be created, resulting in a net current to the probe. Our concern in the following work is to relate the charged particle densities in the undisturbed plasma to the current collected by the probe as a function of the probe voltage.

4.2 MA.THEMATICA.L FORMULATION Since we wish to apply our results to data from a probe in the D-region of the ionosphere, the following assumptions regarding the gas are considered justified: a. The degree of ionization of the gas is small so that the flow of the neutral particle component is unaffected by the ionic currents; b. The neutral particle density is high enough so that ion migration is governed by mobility and diffusion; c. The shock is nonionizing and the flow is frozen; d. The entire system is in a steady state; e. The Einstein relation between mobility and diffusion constants is valid; f. effects of magnetic fields are negligible; g. Each ionic species and the neutral gas may be considered to have a well defined temperature. The following two further assumptions are made to simplify calculations, though neither is essential in principle to the development: a. Only positive and negative ions of equal masses are present in the plasma; b. The plasma components are in thermal equilibrium with the neutral component. Under these assumptions we may consider the system governed by the equations of charge conservation for each species and by the Poisson equal tion, written respectively: V' (V'o'-Kc'V''-DV'a' ) = 0 V~(7t'V't'+KI'VtIDV1) = o (4'2-1) 9

The boundary conditions for the system are: on the probe surface2'5 a = =,' = (p; (4.2-2) at great distances from the probe: t = a oT = t T = O Written in terms of the appropriate dimensionless variables, the constants K, D and e can be absorbed (see nomenclature) and the system of differential equations may be written: V (Va- cV- V-c) = O V.(V + Vvo ):= o (4.2-3) V- VD = while the boundary conditions become: on the probe surface a =' = 0, O = (p, (4.2-4) at large distances from the probe C o = 1, = 0 Throughout the remainder of this report all quantities will be considered in dimensionless form. 4.3 DISCUSSION OF PROBLEM The problem as formulated is complicated in particular by the three flow regimes of the neutral gas which must be considered. Up to the present time, theoretical work has been done only on certain aspects of the problem as a whole. Talbot considers a small planar probe at the stagnation point of a blunt body in a low density supersonic flow. He assumes that the sheath is entirely within the boundary layer and that its thickness is much less than an ionic mean free path. 2 5 Chung and Pollin treat a similar problem, but consider ion motion in the sheath to be governed by mobility and diffusion. Again, the sheath is 10

taken as imbedded in the boundary layer. In addition, Pollin considers ionization caused by a strong shock. At the other extreme, Cohen,4 and Su and Lam5 have given the continuum treatment for a stationary sphere. A more general treatment for an arbitrary low velocity probe in a continuum plasma has been given by Lam. Considering only incompressible flow of the neutral gas, a complete solution to the problem is provided, making no assumptions about the range of the electrical effects of the probe. Unfortunately the various simplifying assumptions made in the above works render their results inextensible to the problem here considered. Indeed, an exact analytical solution to the entire problem seems out of the question. For this reason, this report describes work done toward its solution by a numerical method. 11

5. THE "FLOW LINE" METHOD 5.1 OUTLINE OF METHOD The object of this method is to determine the flow of the plasma components by a series of approximations to the charge and potential distributions about the probe. Starting with the potential field of the probe (taken as a spherical conductor) in free space, flow lines of small elements of the plasma are obtained. Considering the resulting space charge distribution, a new potential field is determined which is in turn used to find new flow lines. The process is to be continued until the variations in succeeding approximations are considered sufficiently small. The current to the probe can then be obtained by calculating the current in the flux tubes containing all flow lines which terminate on the probe surface. To determine the desired flow lines we apply the divergence Eqns. (4.2-3), to current flux tubes. These equations may be interpreted as expressing the constancy of current in any given tube. From these, we take the current densities for positive and negative species: J (V=-5)a-YV -L - (V+VO) V. (5.1-1) Here the velocity field V of the neutral gas is considered known; the flow ahead of the shock is uniform, while that between probe and shock can be determined numerically.718,9 and entered in the computer as a table. In the supersonic region of the flow the convective contribution in Eqns. (5.1-1) will dominate the others. Lam6 has shown, by his "outer" solution, that in such a region the charge densities are uniform, while the potential is harmonic. Thus, we can expect the diffusive terms in Eqns. (5ol-1) to be smaller than the mobility terms except near the probe. Consequently, as a first approximation, we can consider large diffusive effects to be confined to the region behind the shock where the neutral flow is most complex. This suggest that the program be constructed in two phases: I. The determination of flow lines in a region of uniform neutral flow, neglecting diffusive effects entirely; I:I.d Extension of the above to include the flow region behind the shock and density gradient terms over the entire flow field.

This breakdown has several advantages~ a. The phase I program alone should provide good estimates of the actual flow lines in. the region considered. Further because the shock standoff distance is small a reasonable estimate of expected probe current can then be obtained for sufficiently large probe potentials simply by extending the flow lines for the attracted ion species down to the probe surface and neglecting the retarded species altogether. b. By excluding diffusive contributions ahead of the shock analytical solutions for *the initial flow lines can be obtained. Comparing these with the corresponding results of the phase I program, the error in the computed flow lines can be estimated; c. To extend the phase I program to cover the entire flow field (phase II), two essentially similar interpolation routines must be introduced (see Section 8.1). Upon their insertion into the phase I program the general flow line program is completed. This report deals only with the phase I program. 5.2 DETERMINATION OF FLOW LINES-PHA.SE I In the following we treat flow lines of positive particles only; the development for negative species parallels this exactly with appropriate adjustments of signs. 1Neglecting the diffusive term in the expression for the positive ion current, density in Eqn. (5.1-1) we have: 7 (XTv)oD (5.2-L) From thisg the charges can be considered to move as the result of a composite velocity term: V_-V(o (5o2-2) Integrating this over a small but finite time increment bt we have: 5~ = (- 7v)bt + 0(6t2) (5.2-3) where bQ is the change in position of the charge element during time st. Taking a cylindrical coordinate system (r, B, z) with origin at the center of the probe and z-axis directed into the flow, the neutral flow vector and potential gradient become respectively: 13

= -v z (5.2-4) and VO r+- Z, (5.2-5) where r and ~z are unit vectors in the r and z directions respectively. Knowing the potential field (see Section 5.4) the gradients in Eqn. (5.2-5) are found in the program by using a difference method.* Desiring the gradients at some point P, and taking a difference parameter A (see Fig. 7) we have: =0~4-02 + 0(A ) Or p 2A 6( = 1- 3 + o(A2) (5.2-6) where pi represents the potential at a point i at a distance A along lines parallel to the coordinate axes. 40 P 3 52 z Fig. 7. The difference scheme for potential gradient calculations. The quantity A is chosen to be consistent with the point spacings in the trajectory field. Starting points for the trajectories (taken along a line of constant z far ahead of the probe) are spaced so that: br = V St (5.2-7) making the r and z spacings of the trajectory points of the same order. A is *The only variation in this procedure occurs in the determination of the initial trajectory field. In this case an analytical expression for ~ is assumed (see part 6). 14

then written as: A = H V5t (5.2-8) where H is an input constant, usually taken as 0.5. Thus, A is of order St. Consequently, using Eqns, (5.2-6), Eqn. (5.2-3) gives us: - F - r -+ t, (5.2-9) L 2a 2A which is accurate to terms in 5t2. Indeed, since the error incurred through the approximation (5.2-6) will contribute only a term in 5t3 to the above, the second order error term arises only from Eqno (5.2-3), and has the form: 2 =C (vV)v (5.2-10) To obtain the flow lines, Eqn. (5.2-9) is applied at a given point to determine the nextb Reapplying the equation at the new point another is found, and the process is continued until the flow line is terminated. If one considers the flow line to be generated over a fixed time interval tn = nbt, while letting n increase (by taking successively finer and finer steps) then the accumulated error behaves as 2 En =-V) t-, (5.2-11) n 2 so that the total error over the given time interval decreases inversely as the number of steps in the interval. Because of limited storage facilities in the computer, the region of study must be restricted to a finite volume of space about the probe. Consequently2 a region is considered which is bounded in the z direction by two planes normal to the z axis, one ahead of and one behind the probe, and radially by imposing a certain maximum on the initial radial coordinates of the flow lines. The positions of these planes and the value of the maximum radius are obtained below. The initial points for the flow lines can be chosen by considering the flow along the axis. Here r - o, and the initial z coordinate may be determined by requiring that the mobility term (VO) resulting from the starting potential be a certain fraction p of the free stream velocity V: V 15

The initial r coordinate for the flow line farthest from the axis may then be taken as, say, twice the value of z so determined, to allow for the greater range of influence of the potential field normal to the flow. Since the terminal value of z is fixed by the parameters of the problem, as explained below, the value of p is chosen to make the total number of points on the flow lines compatible with the available storage space in the computer. The terminal z coordinate of the flow line, Zminl is determined by other means. In using the first approximation to the potential, it is found that the flow lines of the repelled charge species are deflected so as never to close in again behind the probe. Since this does not occur with the flow lines of the attracted species, the result is an infinite sheath trailing behind the probe. This is clearly nonphysical, and its elimination provides the criterion used for terminating the flow lines. On physical grounds the probe must cause a disturbance in the plasma charge distribution which is commensurate with its own charge. Thus, the charge imbalance in the sheath must be equal to the charge residing on the sphere: Qsphere = Qsheath (5.2-13) Satisfying this condition requires that the charge "tail" behind the probe be cut off at a distance which can be calculated (see Section 6.4). This also restores charge neutrality in the plasma at infinity, as required by the boundary conditions. 553 DETERMINATION OF CHARGE DENSITIES Knowing the flow lines, the variations in charge density along them can be obtained. Since the flow far ahead of the probe is uniform each flow line in that region is parallel to and is characterized by its distance from the z axis. Let this characteristic distance be ro. Recalling the symmetry of the problem, each flow line is recognized as a section of a surface of revolution about the axis. Two such surfaces whose characteristic distances (radii) differ by an infinitesimal amount, dro, form a flux tube, By Gauss' Theorem and the first of Eqnso (5.1-1) we may write for such a tube: SVoj dT JdA = O. (5.3-1) T A

Here we consider a volume T obtained by twice slicing the flux tube transversely. The only portions of the surface of this volume (denoted A) which contribute to the above integral are these transverse surfaces. For convenience we take these surfaces as being everywhere orthogonal to the trajectories. Considering one such surface infinitely far upstream of the probe, and the other a finite distance from it, we may write Eqn. (5.3-1): J, dA, = J dA = I = constant, (5.3-2) where I is the current in the flux tube and the scalar product is neglected because J and dA are parallel. By symmetry and using Eqns. (5.2-1) and (5.2-2), Eqn. (5.3-2) may be written: I = aoV2trodro = a v 2Tr ds, (5.3-3) where ds is the distance between the two "flow line surfaces" at the point (r,z) (see Fig. 8). Since v is known, only ds needs to be determined in order to find the charge densities. r dk~~~~~~~~~~~~~~~ I zWt,r Fig. 8. Parametrization of flow lines. The equation of a given flow line may be written as: z = z(r,ro), (5 3-4) or, solving for ro in terms of r and z: ro = f(r,z). (.5 ~3-5) 17

Considering f as a function of its two independent variables and taking its gradient, we have: IVfI (5.3-6) ds where ds is the magnitude of an elemental displacement taken normal to the curves of constant f. From Eqn. (5.3-5) df- dro, and since we have taken our slicing surfaces normal to the flow lines, the ds of Eqn. (5.3-6) is seen to be the quantity we seek: ds = dro/ (5.3-7) Thus, to determine the charge densities from Eqn. (5-3-3) it is sufficient to know the equations of the flow lines and the local flow velocity v. In the phase I program the densities are not needed explicitly, but are used implicitly in calculating the potential. 5.4 CALCULATION OF POTENTIAL In order to make use of Eqn. (5.2-9), it is necessary to know the potential distribution in the space surrounding the probe. This field will be governed by the Poisson Equation. V24) = _ (5>4-1) and the boundary conditions. 0 = 0p at the probe + 0 at large distances from the probe (5.4-2) Assuming the charge densities known, the potential at any point P may be expressed in terms of the probe potential Op and integrals over the space charge distribution: D(p)a 1- -dq + - dq (.4-3) (P) = P R(P) 4,t Ri 4Tc R2 where dq = (o- )dT (5.4-4) dqt' = dq, 18

and a, R, R1, R2, R(P) are the nondimensionalized distances shown in Fig. 9. r R(P) R dq PROBE \ Fig. 9. Variables used in potential integration. The terms in Eqn. (5.4-3) may be identified as the potentials at P due to the probe, the space charge distribution, and the surface charge induced on the sphere respectively. In order to use this equation in the program for the calculation of potential, the integrals must be reduced to sums over the points composing the flow lines. We consider the integrals over the two charge densities in Eqn. (5.4-3) separately. For either species, the flow lines and the corresponding family of "slicing" surfaces we have chosen form an orthogonal net over the r,z plane. The differential volume element may then be conveniently expressed as: dT = r di ds do, (5.4-6) where 2 and s measure distances along the flow lines and the orthogonals, respectively, and r is the cylindrical radial coordinate of the volume element. Further, the elemental distance along a flow line may be expressed as: da = v dt. (5.4-7) 19

Combining Eqns. (5.3-3), (5.4-6), and (5.4-7) we obtain: a dT = Vrdr r v dt ds do v r ds = ao~V dt ro dro do, (5.4-8) so that the integrals of interest take the form: _ dT= fV P ro dt dro do (549) Ri Ri where integration extends along a flow line (dt), over all flow lines (dro) and around the symmetry axis (dy). Here p and Ri are to be taken respectively as the appropriate charge density and the distance between charge element and field point. Using the flow lines obtained previously (Section 5.2), Eqn. (5.4-9) is well suited for numerical evaluation. The ro of each flow line is known, as is the spacing between flow lines bro. The step size in the flow lines is determined by the input value of St. Thus, the integrations over ro and t in Eqn. (5.4-9) are easily transformed to sums over the points of the flow lines. Only the angular integration need be more thoroughly investigated. Consider a charge element dq and its image charge in the sphere, dq', as shown in Fig. 10. Let the potential be desired at the point P(r,z,o), let P(r, z, o) Fq(rczhi Fig. 10. Variables relating charge with its image charge. 20

dq be at (rl,zjl,), and let dq' be at (r2,z2,O3). It is easily established from the figure that: Ri (zzi)2 2 Ri = (-i) + r2 + ri - 2rri cos~ - Ai - Bi cos i = 1,2 (5.4-10) where Ai and Bi have the obvious definitions. From the position of the image charge (see Fig. 9) it is known that: 2 a r2 = rl 2 2 rl.z1 2 a Z2 = 2 zle (5.4-.11.) 2 2 rl+zl Using relations (5.4-10) and (5.4-11), all integrands in integrals of the form (5.4-9) can be expressed in terms of the coordinates of dq and the point P, i.e., in terms of the flow line point and the field point. Using the notation of Eqn. (5.4-10) the azimuthal integral in Eqn. (5.4-9) takes the form 2it do do (5.4-12) S R 4 o tl —coWs where Bi Ai The integral on the right side of Eqn. (5.4-12), may be expressed in terms of the complete elliptic integral of the first kind: 2T Tcl/2 O/1-Clicoso o C i_ li ir] o0 11 oi~ 0 l- 2i sin2~ 4 K'2 (5.4-14) Thus, Eqn. (5.4-12) may be written as: 21

do 4 / (K Z (54-15) Ri \A +B l+c In the program this elliptic integral is evaluated by the series: 4 4 = ) an 1 — - in (- n i -i (5.4-16) n=o n=o 10 where ao = 1.386 294 361 12 bo =.5 al =.096 663 442 59 bl =.124 985 935 97 a2 =.035 900 923 83 b2 =.068 802 485 76 a3 =.037 425 637 13 b3 =.033 283 553 46 a4 =.014 511 962 12 b4 =.004 417 870 12 Evaluating this series in the IBM 7090 computer using roughly 9 significant figures, the results are considered accurate to about one part in.10'. It is noted that the above integral diverges as the logarithm of (1-ci)-l when ci + 1, i.e., as the field point approaches a flow line point. Since this divergence is nonphysical, such situations must be treated by a special method. In evaluating the potential integrals we obtain Eqn. (5.4-12) by considering the charge associated with a given flow line point to be concentrated on a circle concentric with the symmetry axis, rather than distributed in a toroid of finite volume. Since the potential of such a ring diverges as one approaches it, the misbehavior of our integral is explained. Clearly then, our method breaks down when the field point P enters the torus associated with the flow line point (Fig. 11), This provides a criterion for changing the procedure used for finding the potential due to the toroid in question. Turning again to Fig. 11, the toroid associated with the point Q will generally have a roughly parallelogram-like cross-section with an area of order (vSt)2o Since the exact shape of this figure cannot be determined without greatly increased complexity, it is approximated in the program by a circle about the point Q having area (vbt)2. The radius of this circle is then given by: 22

r = vt (5.4.18) where the value of vSt is taken as the separation of the trajectory point Q and the preceding point S. If a point P falls within distance r of Q, the potential at that point is found in the following manner. o Flow line points vSt r I z Fig. 11. Section of charge toroid. The line PQ is constructed and extended through P to intersect the above circle at T. The potential at T is found in the usual manner, and the potential at Q is found as explained below. Then the potential at P is taken as that obtained by linear interpolation between the potentials at Q and T. To obtain the potential at a point Q on the central circle of the torus, the torus is first approximated by a series of spheres centered on this circle. The potential due to this new charge distribution is taken as an approximation to that of the original toroidal distribution. The torus is assumed to have major radius R, and minor radius r, given by Eqn. (5.4-18) (see Fig. 12). Measuring from the point Q the axis is cut into segments having length R=X 2iR (5.4-19) _r where [rtR/r] signifies the greatest integer in tR/r. At the end of each segment determined as above, a sphere of radius r is constructed. Each such sphere is considered to have charge given by: 23

q q (5.4-20o) LrR where q is the total charge of the complete torus. The flow lines have been generated in such a way that this charge is identically equal to the charge in a torus at the top of the flow line considered, i.e., q = 2t ro bro V St. (5.4-21) Fig. 12. Construction for integration when field point lies within charge toroid. Ca lling ~n = n, n = 0,1,..., [T (5.4-22) n b~ n =r the potential at Q due to this distribution is [E q(-1 = o+ X ( 5.4-23) n4-i l 2R2-2R2cos ~n where 3 ro 5ro V St = 5r R]Vt (5.4-24) 4r [7rc 24

represents the potential at Q due to the sphere of charge 3q centered at Q. Thus the total potential at Q due to the torus is: [" El-1 ro V ro t 3 + 1 (5.4-25) r n=l l-cosjn The above relations now allow us to calculate the potential at any point P due to the charge distribution given by the flow lines. The potential expression (5.4-4) becomes: )(p) a= p a V bt (7 ro t )l+l a ( I_+U P R(P) it L L R b4A2+ [I R]-l (5.4-26) v6t(ZIro bro 3 11 L _ 2 4( a )[2 r ] { 2 R 4/cosCni r n=l where all symbols are as previously indicated, and j are sums over all flow line points of the a and ~ species respectively except points which fall into the case treated by Eqn. (5.4-25), which are accounted for by the sums E and L The values of ro and 6ro which are used in the program are computed from the r coordinates at the tops of the flow line considered and of that with the next larger value of r; ro is taken as the average of these values, and 5ro as the difference. These calculations are performed for the farthest trajectory from the axis by considering another to be located 5ro farther out, where bro is taken as the last calculated value of that quano tity.

6. INITIAL APPROXIMATIONS USED IN PROGRAMMING THE FLOW LINE METHOD 6.1 INTRODUCTION In order to start the iterative procedure of the flow line method it is necessary to assume an initial potential field or charge distribution. We assume here an initial potential field: that due to a sphere of potential 0p in free space. This is not only reasonable physically but it also enables the flow lines and charge densities of the first iteration to be derived analytically, giving an insight into the operation of the procedure as well as an estimate of the errors involved due to finite step size. 6.2 INITIAL FLOW LINES In a cylindrical coordinate system the initial potential is given by c, =, a r2+z2 > a2 (6.2-1) P r2+Z 2 The paths followed by the charge elements will be described by a function z = z(rro) which, for given ro, will be everywhere tangent to the local velocity vector v, as given by Eqn. (5.2-2). This function is determined by integrating the differential equation: dr vr (6.2-2) dz where Vr and vz are the r and z components of the vector V. Using Eqns. (5.2-2), (5.2-4), (5.2-5), and (6.2-1), we find these velocity components to be: Op a r r (r 2+z2)3/2 (6.2-3) vp a z V = V+ -V.p. Z (r2+z2)3/2 Substituting these relations in Eqn. (6.2-2), the resulting differential equation is easily integrated. Requiring that r - ro as z -+ +oo the constant of integration is evaluated, and the flow line equation is found to be: o~2N 2 r2 6.2 (62-4) 26

where: Op a - a (6.2-5) V The equation for flow lines of the negative species is identical to the above except that the sign of the right side is changed. Taking P positive (positive probe) the flow lines of positive particles are deflected away from the probe, while those of negative particles are attracted, which results in a region where one species is excluded, i.e., a sheath. It can be shown from the equations of the flow lines that this region starts at a distance ahead of the origin equal to v fJI that it extends behind the probe an infinite distance, and that it is always of finite width, having a radius of 21 jIP at an infinite distance behind the probe. For the ionic species which is attracted to the probe, this radius also corresponds to the ro of the flow line farthest from the axis which is "collected." In fact, it is easily shown that the flow lines of the positive and negative species are merely the mirror images of one another across the plane z=o (see Figs. 13 and 14). ir Fig. 13. Flow past sphere, ( i? < a ). Fig. 14. Flow past sphere, (, 11 > a ). 27

6.3 CHARGE DISTRIBUTION RESULTING FROM STARTING FLOW LINES From Section 5.3 we know that the charge density distribution is given by Eqns. (5.3-3), (5.3-5), and (5.3-7). Combining these relations and using the fact that ox, = 1 we obtain V IV roI v. 2'(6.3-1) since by Eqn. (5.3-5), f - roSolving Eqn. (6.2-4) for ro2 and using Eqns. (6.2-3) we find that 2 Or V (6.2) o? v so that Eqn. (6.3-1) reduces to: ~~a 1~~~.1~ ~(6.3-3) A similar relation holds for ~, so that the charge distribution is everywhere uniform except in the region where one species of particles is excluded. The net density is therefore zero everywhere outside the sheath, and unity inside, where it has the sign opposite to that of the sphere potential. 6.4 TERMINATION OF FLOW LINES As mentioned in Section 5.2, the nonphysical nature of the infinite trailing sheath provides a criterion for terminating the flow lines. Since both charge densities are uniform everywhere outside the sheath, terminating the flow lines amounts to cutting off the sheath and restoring neutrality to the flow behind it, which also satisfies the boundary condition of charge neutrality at infinity. Our requirement is that the net charge in the sheath shall be equal in magnitude to the charge on the sphere due to its natural capacitance: Qsphere = Qsheath (6.4-1) Considering the spherical probe to be initially in free space, the none dimensionalized charge on it due to its own capacitance is easily found to be Qsphere = 4t an p (6. 4-2) 28

Here the charge is normalized with respect to the electronic charge and n, is the number of charged particles of either species in a cube at infinity having side equal to a Debye length, X. To determine the sheath charge for use in Eqn. (6.4-1) we must first find the boundary of the sheath. Figure 13 indicates some of the flow lines of the first approximation in the vicinity of the probe when 11J < a2. We note that the probe intercepts flow lines of both species. Assuming the probe to have a positive potential, the curves labeled ai and 5i represent the flow lines of positive and negative charges respectively which just graze the probe surface. All flow lines of both species having characteristic radii ro less than those of the corresponding flow line above are collected; all having values of ro greater than those of the above are not. Thus the sheath is the region (shaded) between the line Fi, and the line ~i or the probe, whichever is appropriate. In this region there is negative charge only. To determine the outer edge of the sheath we must find the flow line ail which has a point of tangency with the probe. Let this line have characteristic radius r'. At the required point of tangency, call it (r',z'), the following conditions are satisfied: r'2 + z'2 a a2 (6.4-3) and: r r2 z 2 rt - r i2- o (6.4-4) so that the quantity _ Z2 is equal for both also. Finding z2 from Eqns, (6.43) and (6.14-4), caculating these derivatives and equating the results we obtain: (4p2_y2)2 = y4.4y y2+ 8P2r12y (6425) where: y _r2 rT2 + 2P. (6.4-6) Solving the system of Eqns. (6.4-3) through (6.4-6) we obtain r' = Ba& L. a2 Z! a (6.4-7) r' = a -2 29

The outer boundary of the sheath is thus determined. The inner boundary of the sheath consists of two parts. Because of the symmetry of the flow lines of the different species about the plane z=o, the line Si of Fig. 13 is tangent to the sphere at the point (r', -z'). Thus, the inner boundary of the sheath is the probe surface between the points (r', z') and (r', z'), and is the flow line 5i for all values of z < -z' The value of ro for this flow line, call it ro, can be found rather simply in terms of r'o. We consider the flow line equation for negative ions: r2 - r2 - (6.4-8) Letting z go to minus infinity and denoting the value of r in this limit by r, we have: ro - r 4P. (6.4-9) Because of the previously mentioned symmetry of the flow lines we know that for the Ci trajectory r ro = ro, so that Eqn. (6.4-9) gives us: r112 4P + rt 2:. (64-10) Thus, the boundaries of the charge distribution in the sheath are known. With this information we can now turn to the calculation of the sheath charge as required in Eqn. (6.,4-1). Considering all flow lines terminated at some plane z = Zmin, the sheath charge may be expressed in terms of integrals, z' r2(z) 2 2 dz d r D< a; z > Zmin > a 2 Zmin a -z Qsheath = "n 2'6.4- 11 z' r2(z> -z' r2( z) 2 2 2 dz d r + dz d r d<a; zmin<-a 2 -z a -z Zmin r1(z) 2 where n represents the number of ions per Debye cube in the sheath, and rl (z) and r2 (z) are obtained from the equations for 5i and ai respectively, as shown below. Lettiing = r2 + 2 (6.4-12) 50

and expanding either of the equations for ai or Si we obtain a single cubic equation in r: (r2)3+(z2-2_) ( r2) 2+[( _-2z2) ( r2)+z2(24 _ 2) 0. (6.4-13) For values of the parameters in this equation which are of interest to us (IPI > 0; ro2, z2 > 0) there are three real solutions for r2: r2 k,z) = (rA+z2+2P) 1 +OS C - ((r+z2+2)3 l 3 3 k = 0,1,2 (6.4-14) Considering appropriate limits of these solutions we find that: 2 2 rL = r (0,z) (6.4-15) r2 - r2(2n\ r2 = r- 2,z while the solution r2 (l,z) is physically impossible. In the event that JIP > a2, Fig. 13 is altered (see Fig. 14). In this case all of the a flow lines miss the probe entirely. In addition, there is no region behind the probe which is free of charge, so that the range of integration for the sheath charge integral extends from the axis to ai except between a < z < -a where it is bounded on the inside by the sphere. In this case flow line ai has r' = 0, so that j in Eqn. (6.4-12) reduces to 2o. Equation (6.4-13) then becomes (r2) [(r2)2 (2+4z )(rr )+4p(5-z_)] = O (6,4-16) and the solutions of interest are: 2 rl(z) = 0 z2 z r2(z) = 2P - 2- - /8+z2 (6.4-1- Using these results, the sheath charge integral now takes the forms 31

r2(z) S dz S d r2 J1 >a2;! > Zmin > a Zmin 0 Qsheathd22 NJT(1' ) a -z Qsheath = S-dz = d r2 dz A r2 a; a > min > -a o...o Zmin Zmin (6.4-18) 222 1TB~ rpa(z) a a -z dz d r2 -fdz d r2 ]~ > a; Zmin < -a zmin 0 -a 0 2 where r2 (z) is given by Eqn. (6.4-17). 2 Thus, depending on the relative magnitudes of a and I, the appropriate form of Eqns. (6.4-11) or (6.4-18) is used with Eqn. (6.4-2) in Eqn. (6.4-1). This yields an equation for Zmin which involves only the parameters of the problem; when solved this provides the stopping point for the flow lines. In the case where i11 < a2 the evaluation of the charge integral is ex2 2 tremely difficult because of the complicated nature of the limits rl and r2. However, if Ij1 < a2 an approximation may be made which reduces the difficulty considerably and the integral may be evaluated. When J1J > a2 the sheath charge integral can be evaluated directly. We now consider these cases separately. If IP << a2 (a >> ) the two flow lines determining the sheath will be tangent to the probe very near (r,z) = (a,o). In the limit of p = 0 these flow lines will satisfy this condition exactly and will be everywhere parallel to the z-axis. We consider cases near this limit, taking!PJ to be small. Demanding that the two flow lines of interest pass through the point (a,O), the flow line equations give us the relations 2 2 ro - a + 2 (6.4-19) where the upper sign is valid for flow lines of positive species, the lower for negative species. Substituting these results in the flow line equations we obtain: 2 2 +2__z_ a -r = (6,4-20) In the limit we are conrsidering r will not deviate appreciably from its value at z = 0, so we set r = a in the radical of Eqn. (6,4-20) to obtain: 32

2 2- 2Pz r - a a. (6.4-21) ra2+ 2~ This gives a good approximation to r2 for all z, since near z = O, r is very nearly equal to a, and for large negative z the importance of the r2 in the radical is diminished so that the error introduced by the substitution r = a is very small. Using Eqn. (6.4-21), the sheath charge integrals (6.4-11) reduce to: 2 2Pz Qsheath min dz Zmin a + +2z la 2Liz - 4gn~a ( l1+x2 -1 ] (6.4-22) where: Zmin x a m (6.4-23) Substituting Eqns. (6.4-21), (6.4-22), and (6.4-23) into Eqn. (6.4-1) and rearranging terms we have: 2 = J f(x) x < 0 (6.4-24) which enables us to solve for Zmin. We note that the asympototic form of f(x) for large negative x, call it fO(x), is f,(x) = -x-l + - (6. 425) When IP1 > a2, we use Eqns. (6.4-18) for the sheath charge. These can be immediately integrated to give Qsheath' - - n [t3 + (8+t2) -12t-161 0 1> y l>t>/ + gnP B3/2~~- (6.4-26) + 53 / t3-3y2t+273 1> y2 y>t>-7y 4y 1> y2 t<-7 where 33

Zmin a (6.4-27) = (6.4-27) Substituting Eqns. (6.4-27) and (6.4-26) into Eqns. (6.4-1) and (6.4-2), and rearranging terms we have FD 2 247 = g(t) - 2h(Q,t) 1>7873 1> 72 t<-7 (6,4-28) where g(t) = t3+(8+t 2)3/ 2-12t-16 (6.4-29) 3 2 3 h(y,t) =t -3Y t+2y ~ (6.4-30) The asymptotic form of g(t) is ga(t) = -24t-16 t + -. (6.4-31) Having obtained solutions of Eqn. (6.4-1) for I|1 > a and 1n << a2 it is of interest to see whether the results can be adapted to the intermediate case of I<1 < a. Considering the relations for x, t and 7 (Eqns. (6.4-23) and (6.4-27)) we have t = 7x. (6.4-32) Defining~ _A = —, (6.4-33) Eqn. (6.4-24) can be solved for t: t = f= - ~(n+l) -l1. (6.4-34) Turning to Eqn. (6.4-28) we consider values of 7 beyond the stated range, i.e., 7/ > i, From the nature of h(y,t) we know that 24v1y < g(t) < 24TIZ+8y3, (6.4-55) and hence, by the monotone property of g(t), that: 34

Using Eqn. (6.4-34) and the inequalities (6.4-36), curves of t vs y are given in Figs. 15, 16, and 17 for three values of il. The curves have been extended beyond their ranges of definition in y for purposes of comparison. These indicate that for values of ~ > 10 the percentage variation between the values of t determined from f and g is quite small in the intermediate range of 1 < y < 3. Thus, extension of the use of g(t) to values of up to Y = 1.2 (see Fig. 15) in the form g(t) = 24ny+8y 1 < y < 1.2 (6.4-37) and the use of f(x) for values of y > 1.2 can be considered to give a good approximation for determining Zmin in the intermediate range of y. The restriction that ~ > 10 is not a serious one in practical cases of interest, since it will be satisfied by most rocket mounted probes, e.g., considering a probe at an altitude of 70 km with V' > 100 m/sec, V' = 10, a' = 1.3 cm, T' = 300~K, and n' = 1.3 x 104/cm3 of singly ionized E2, we find that ~ > 60. In order to obtain fast approximate solutions of the Zmin equations, graphs of the functions f(x), g(t), and g(t) - 2h(y,t) are given in Figs. 18 through 22. To obtain a graphical solution for Zmin the value of IJI is calculated from the given data of the probelm and compared with a. If Ii < a2, the parameter on the left side of Eqn. (6.4-24) is computed and the appropriate value of x determined from the graph of f(x) in Fig. 18, or from the asymptotic relation for f(x). Then Zmin is found from Eqn. (6.4-23). If JIP > a2, the parameter y is evaluated from Eqn. (6.4-27). This determines a unique "branch" of the g(t)-2h(y,t) curve of Fig. 22. Computing the parameter on the left side of Eqn. (6.4-28), this figure should be consulted to find the appropriate value of t, interpolating between the curves for different values of 7 if necessary. In the event that the parameter 247 yP is so large that the g(t)-2h(y,t) curves are not extended sufficiently far along the negative t axis to allow solution, 8y should be added to it, and the value of t read directly from the g(t) curve (Figs. 19, 20, and 21) or obtained from the asymptotic relation (6.4-31). When the values of a and P are such that y falls between the limiting cases discussed above,in the curves are used in accordance with the results discussed in the paragraph following Eqn. (6.4-36).

(T=U) 14 SA. JO GAins ~T * 2TJ 81lI= & - 91Zl( /+X + &tz)1 6 =.'018(&) A= \ 1 9~I- %I..- 2-~~~~~~~~~~~g

F'; F` ~ ~ ~ U0 r~ ~~ ~C~~~~~~~~~~~~~~~~C Omft I*Awwo~~~~~~~~.,U.~~~~~~~~~~~~~~~~~~~~~~U i'3 4-. ~ ~ ~ 4 00 0 11' S~~( 3 3 41 OF~~~~F OOj r CM (V PC).CCI OD~~~~~~~~~~~~ >*I',. >*.I'.~ P*~~~I (~ 1 (~ 1 r~~~~r) I~1') --- I I I ~ ~ ~ ~ ~ ~ ~~~~~~~~I I I I to 0 to 0 to 0 to~~~~~~~~~~~~~~I=

1 2 3 -50 -100 - -150 -200 -250 -I t=9 (241rqy) -300 Fig. 17.* Curve of t Vs r, (R=100).

4.5 4.0 - 3.5 3.0 f(x) f(x)=TVl+x2 -I'\ -~ 2.5 fc (x)=-x-I2. 2.0 1.5 1.0 0.5 -6 -5 -4 -3 -2 -I 0 Fig. 18. Curve of f(x) vs x. 39

60 50 g(t) = t3+ 8(2+ t2/4) 3/2-12t- 16 40 \g(t) -t=-~ 90D(t)= -24t- 1630 0 20 10 -3 -2 -I I Fig. 19. Curve of g(t) vs -t, (y=0).

220 200 180 160 g(t) 140 g(t) 120 g (t) 100 80 -9 -8 -7 -6 -5 -4 -3 Fig. 20. Curve of g(t) vs -t, (7=0).

340 320 300 280 g(t) 260 g(t) 240 a(t) 220 l I I I I I 200 -15 -14 -13 -12 -11 -10 -9 t Fig. 21. Curve of g(t) vs -t, (y=0).

16 h(Yt).85 y o1.09-95'0 BO75.70.60.50 12 25 -10 I I -1.0 -.8 -.6 -4 -.2 0.2.4.6.8 1.0 t Fig. 22. Curve of h(yt) Vs -t, (#o).

7. COMPUTER PROGRAM AND NUMERICAL RESULTS 7.1 COMPUTER PROGRAM The input parameters for the program are: 1. V = probe velocity (dimensionless) 2. PHI = 0p = probe potential (dimensionless) 3. DELTA.T = bt (dimensionless) 4. A = a = probe radius 5. ZMIN - Zmin (dimensionless) 6. H = spacing parameter for derivative calculations. 7. NUM = number of starting points to be used in generating flow lines. 8. NMAX = maximum number of points per flow line, set by computer storage restrictions. 9. NITER = number of iterations to be made after computing initial flow lines. 10. PSW = intermediate output switch. 11. (RI, ZI) = starting points for flow lines (NUM pairs will be read). From the starting points (RI, ZI) the program first computes an initial approximation to the flow lines from the neutral flow velocity and probe potential. These flow lines are stored in two pairs of arrays: RSA., ZSA and RZA, ZZA.. The first letter of each array name indicates SIGMA (positive charge) or ZETA (negative charge). Next a new approximation to the flow lines is made by using the same starting points and taking into account the charge distribution from the previous approximation. These new approximations are stored in arrays RSB, ZSB, RZB, and ZZB~ If the number of new approximations is less than NITER, the "B" arrays are shifted into the "A" arrays and the next iteration is calculated and stored in the "B" arrays.

During calculation of all approximations to the flow lines, a flow line is terminated at the last point before any of the following occur: 1. The calculated increment in Z is positive for the flow line along the vertical axis. 2. R is negative. 3o Z is less than ZMIN. 4. The flow line intersects the probe. 5. More than NMAX points are calculated. After each approximation to the flow lines is completed, it is written on the output tape. Because each iteration uses many minutes of computer time, it is sometimes desirable to print out each new flow line point as it is generated. Assigning a non-zero value to PSW accomplishes this, and allows the program to be stopped at any time without loss of output. For details of the method see the flow diagram and program listing in the Appendix. Subroutines which are not expanded in the flow diagram or the program listing are part of The University' of Michigan Computing Center Executive System. 7 2 PRESENTATION OF SAMPLE RUN At the time of writing this report, results of programming the procedure previously discussed are not satisfactory. The output of a run using tbhe program discussed in Section 4.1 is presented in Figs. 23 through 26; the dots indicate the computed points. Note also that the r and z axes are scaled differently. Herethe values of the parameters used are V=15, Op-+10 a=-35 t= - H= o It is found in this case that 1J <' a2, so the sheath cutoff distance is determined using Fig, 14; this gives Zmin = - 17,5 The starting coordinat;es for the flow Lines are arbitrarily chosen to be z=+309 r=0,2,4, o 120 Figures 23 and 24 show the flow lines as computed using only the free space potential due to the probe. Figures 25 and 26 show the first iteration, using the fields of the previous figures in addition to the potential of the probe. Total running time on an IBM 7090 computer for the two iterations was 15 minm The results shown in Figs. 25 and 26 have several salient features~

R v =15 15 )p =+ 10 a =3 St =1/6 30 20 10 0 -10 -20.Fig. 2.5. a flow lines (inital a1pproximation).

-0 o 0~~ -F 0.~ t 0~~~~~~~~~~~~~~~~~ /~~~~~ L. 0 N 1 I~~~~~~~~~~~~~~~~~~~~.~o~~~~~~a Ino to e I I I I II II 0 -40 ro~~~~~r

V =15 (pP= + 10 a=3 Bt = I/6 ~-I0 o i 30 20 10 0 -10 -20 z Fig. 25. a flow lines (first iteration).

V =15 -15 (p= + 10 a =3 St = 1/6 10 x * * ~* * x @.... v * 30 20 10 0 -10 -20 Fig. 26. ~ flow lines (first iteration).

1. The sigma flow line starting with r=6 has several irregularities near the lower end; 2. the magnitude of the overall interaction is markedly increased, as exemplified by the highly distorted flow lines; 3. near the probe the positive (negative) species are attracted (repelled) in the second iteration; and 4. the strength of the interaction seems to increase with distance from the probe, as indicated by the increasing distortion of the flow lines which begin with successively larger values of r. At present, these effects, all of which appear to be non-physical have not been satisfactorily explained. Until these are resolved, further iterations do not appear warranted.

8. EXTENSIONS AND CONCLUSIONS 8.1 EXTENSION TO PHASE II The extension of the phase I program to include nonuniform flow and diffusion can be achieved by considering these effects independently. In principle the complications necessary in the program because of the change in the neutral flow regime at and behind the shock are of a relatively simple nature. In passing through the shock both the flow velocity and the density change discontinuously. Thus, the first necessary alteration to the phase I program is the determination of the intercepts of the flow lines and the shock. This involves checking the points along a flow line to locate the first point falling behind the shock, and then by say, linear interpolation, finding the point of intersection P of the flow line and the shock front (Fig. 27). The point so located should then be used as an "initial" point for further determination of the flow line behind the shock. r o Flow line points o~o ~/SHOCK PROBE z Fig. 27. Flow lines crossing shock. The addition of an interpolation routine is also necessary in order to determine the neutral flow velocity at flow line points behind the shock. This will be required because such points cannot generally be expected to coincide with any of the points at which the neutral velocity is available, since the velocity field is entered in the computer as a table. Such an interpolation routine might be constructed as follows. 51

Consider a flow line point P at which the neutral stream velocity is desired (see Fig. 28). Let the three points nearest to P at which the velocity o Points in neutral flow field r * Points on flow line z Fig. 28. The three points of tabulated velocity field which are nearest to point P. field is tabulated be indicated by 1, a, and 3. Making a Taylor series expansion in the components of the velocity about the point P we have to first order: Vri =;+ + a6rzlp i = 1,2,3 (8.1-1) where Vri denotes the r component of the velocity vector at the point i and 5rin 5zi are the appropriate displacements of the points i from point P. This linear system can now be solved for the radial velocity components at the point P. A similar procedure can be used to find the required z-components of the vector V. It should be noted that the real difficulty in the above is the machine determination of the nearest neighbor points 1,2, and 3. This could perhaps be done by searching the velocity table for points within a specified distance of P, which is, say, twice the maximum distance between points in the velocity table. Decrementing this distance by small amounts, some of the points so determined could be excluded until the three necessary points were found. The contribution to the current density due to density gradients can be found by an interpolation technique similar to the above. Rearranging Eqns. (5.1-1) the current densities may be written:

J = (V - V ) _L = (V+o V _ ) 1, (8.1-2) We see then, that the contribution to the velocity of the charges due to diffusion is given by the ratio of the density gradient to the density itself (assuming the density is finite). Setting up a linear system analogous to Eqns. (8.1-1) for the charge densities, these contributions can be easily found. Since we assume the ionized flow to be frozen, the ion densities will change directly with the neutral density increase across the shock. For computer calculations we must therefore consider the neutral flow velocity and the densities to be double valued at the shock (point P of Fig. 23), using the shock as a bounding surface between the two regions of the flow. In calculating charge density gradients ahead of the shock, the "nearest" neighbor points considered must be restricted to the region ahead of the shock or on the shock itself as approached from the front. Similarly, behind the shock calculations must use only values in the region behind the shock or values along the shock which have been corrected for the changes caused by passing through the shock. Further corrections might also be included for such effects as the variation in diffusion (and thus mobility) constants due to the density variation behind the shock. 8.2 FURTHER EXTENSIONS AND CONCLUSIONS The major portion of this report has been concerned with a spherical probe, and with the free space potential of such a body as the potential approximation initiating the iterative procedure. In principle the extension of the program to utilize different configurations and potentials is straightforward. Changes would be necessary in sections of the program in which the initial flow lines are computed (by changing the initial potential function). Even if the spherical probe is retained, an extension of the process used for determining the flow line cutoff point zmin would be desirable. In general the charge distribution on and about the probe will be altered with each iteration. Consequently, new calculations of the charge imbalance in the plasma and of the net charge on the sphere (including for consistency the charge induced by the previous spatial distribution) should be used to obtain a new approximation to Zmin. This would require an additional section in the program for performing the necessary numerical integrations and for obtaining the dew sired cutoff value from them. The numerical integrations might easily utilize the density computations mentioned in connection with the phase II extension of the program. 53

Some serious points concerning the validity of the entire procedure must be mentioned. The whole method evolves from the observation that if the potential 0 is known, then the system (4.2-3) can be solved for a and A, or alternatively if the charge densities are known, the Poisson equation can be solved for the potential. This suggests solving the system by successive approximations, as has been done in this report. This process has the advantage of resolving the problem into the solution of a sequence of linear equations, rather than of a coupled system. However, this approach immediately raises the following questions: do there exist functions which can be used as initial potentials for which the process will ultimately converge to the solution of the system (4.2-3) with its boundary conditions, and if so, what are they? These questions remain unanswered; results of running the computer program evolved in this report (see Section 7.1) are not sufficient to indicate possible resolution of these important points at this time. 54

BIBLIOGRAPHY 1. Talbot, L., "Theory of the Stagnation-Point Langmuir Probe," Phys. Fluids, 3, 289, 1960. 2. Chung, P. M. Electrical Characteristics of Couette and Stagnation Boundary Layer Flows of Weakly Ionized Gases, Aerospace Corporation, Report No. TDR - 169 (3230-12)TN-2, October 1962. 3. Pollin, I., The Stagnation-Point Langmuir Probe in a Shock Tube-Theory and Measurements, Harry Diamond Laboratories, TR-1103, June 1963 (AD 416 346). 4, Cohen, I. M., "Asymptotic Theory of Spherical Electrostatic Probes in a Slightly Ionized, Collision-Dominated Gas," Phys. Fluids 6, 1492, 1963. 5. Su, C. H., and Lam, S. H., "Continuum Theory of Spherical Electrostatic Probes," Phys. Fluids 6, 1479, 1963. 6. Lam, S. H., "A General Theory for the Flow of Weakly Ionized Gases," A.IAA. Jour., 2, 256, 1964. 7. Lieberman, E., Description of IBM 709/90/94 Computer Programs and Analysis for Flow Fields About Bodies of Revolution in Hypersonic Flight, General Applied Science Laboratories, Technical Report Noo 34o, May 1963. 8. Garabedian, P. R,, Numerical Construction of Detached Shock Waves, Stanford- University, Applied Mathematics and Statistics Iaboratory, Technical Report No. 55, September 1956. 9. Nelson, PO, Jr., An Iterative Method for Approximating the Flow Field About a Supersonic Blunt Body, Brown Engineering Co., TN R-72, October 1963 (AD 435 078). 10. Hastings, C., Jr., "Approximations for Digital Computers," Princeton University Press, 1955.

APPENDIX FLOW DIAGRAM AND PROGRAM LISTING

START ) Tl —1 - ( RSAJ,K + ZSAJ, K-3/2 Read V, PHI, DELTAT, A, RSAJ, K+1'RSAJ, K T1 ZMIN, H, NUM, NITER, PSW ZSAJ,K+I-ZSAJ,K ~ T1 and( for J =1, 2,..., NUM) RSAj,1, ZSAj i Make program modifications and Is R SAJ,K+1 0 and initializations ZSAK <ZSAJK+1 or RSAJ,K+1 <0 or ZSAJ,K+1 < ZMIN or does the line segment J 1 I [(RSAJ,K,ZSAJ,K ), (RSAJ,K+IZSAJ,K+1 ) intersect the probe or is K = NMAX? Begin flow line J of RSA and ZSA Yes No K|- K1 |J I-J + 1 K- K + 1 Begin point K of flow lineJ Is J > NUM? No Ye s 1 4

Print initial SIGMA flow Begin flow line J of lines from RSA and ZSA RSB and ZSB K 1 Set J —-1 and repeat the steps from ( through 7 ( the preceding step for ZETA, using RZA, ZZA, line J ZETA; using RZA, 2ZZA, J t Begin point K of flow and T1-*-1 + ( RZA2 + ZZAJ,K2 )-3/2 JK No Is PSW n0- I.Print Yes Yes RSAJK Have NITER Zes SAJ,K iterations been Start DAYTIM made? No RTEMP —- RS BJK-DELTAI For J-1, 2,..., NUM ZTEM P —-Z S BJK VECJ-.(RSAJ+1,+RSA,l *(RSAJ+,1- RSAj,1) T2. -POTNTL (RSA, ZSA) -POTNTL (RZA, ZZA) VECNUM (2 RSANUMI+RSANuM,t ) (RSANUM,I -RSANUM-1,1) | 6.,-

RTEM P —RSBJ,K +DELTA TI —(RSBJK-ZSBJK )3/2 RSBJ,K+I-RSBJ,K + RGP ZSBJK+I,- ZSBJ,KT1 + ZGP RG P —R S B J,K - V ~ DELTAT POTNTL (RSA, ZSA) T.POTNTL (RZA, ZZA)-T 2. DELTA * 4ir Is RSBJ,K+1 = 0 and ZS BJ,K <ZSBJ,K+1 RTEMP —-RSBd,,K or RSBJ,K+1 < 0 |ZTEM P --- ZS BJ K + DELTA Ior does the line segment [URSBJ,K ZSBJ'K )] (R S BJ, K+I,ZS BJ, K+)] I T2 —POTNTL (RSA, ZSA) intersect the probe -POTNTL (RZA, ZZA) or is K = NMAX? Yes No ZTEMP —ZSBJ,K + DELTA J —J+ 1 K-' — K+I ZGP- ZS BJ,K rPOTNTL (RSA, ZSA) 1 1 S J>NUM POTNTL (RZA, ZZA)-T2 No Yes 2- DELTA ~ 4r1 ( (

Set J4 —1 and repeat the For J — l, 2,..., NUM steps from ( through and K —1, 2,..., NMAX the preceding step for RSAJK RSBJK RZB and ZZB using ZSAJ ZSBK RG P-RZ BJK K [POTNTL ( RSA, ZSA ) 1 RZAJK -RZ BJ,K + POTNTL ( RZA, ZZA )-T2J ZZAJ,K,-ZZBJ,K 2 ~ DELTA.4ir ZGP —- ZZBJ,K POTNTL ( RSA ZSA 1 Print the new flow line [Pf+l~( SA ZSAm~ iteration from the "A" 2 DELTA.4 r matrices 2 ~ DELTA.4ir RZ BJ,K+1 — RZBJ,K *TI+ RGP ZZBJ,K+li- -ZZBJK T1+ZGP -V * DELTAT

POTNTL TT20 —-A2 / (TEMPR2 + TEMPZ2) TEM PR. TEM PR TT20 Function of RXYA, ZXYA TEMPZ TEMPZ. TT20 12 TT4-a M - (TT20 ) 1/2,PONT TT 7 M w- M+1 r' ~~~ I ~IstheM-lpoint No the last point of 1 the L flow line? TT4'0 Yes I TT7 ~ —TT7 — + VECL TT4 TEM PR - RXYA L,M TEMPZ - ZXYAL,M J - J + 1 TEMPRR -- RXYA L,M -1 TEMPZZ *- ZXYA L, M-1 Is J >NUM? NO IYes TT44 —TT4 + PONT C (AC) -z1/22 V' DELTAT.TT7 TT20 - A/I (TEMPRR2 2 + TEMPZZ2) RETURN TEM PZZ - TEM PZZ * TT20 TEMPRR w-TEMPRR. TT20 (14)~~~~~6

PONT TT3' (RTEMP - TEMPR )2 + ( ZTEMP - TEMPZ )2 Yes Is TT3 > (V.DELTAT)? i) No TT8 —-[TEM PR r I TT5] DPHI-, — 2r / TT8 TT2 [(TEMPRR - TEMPR )2 TT110.11TT8 {3 I 2 (TT2) + (TEMPZZ- TEMPZ )2] 1r + 11V T EMPR Yes TT-[1 - COS (i. DPHI )]/2} Is TTr2 >TT3? ~ l No 15 | CAPB - 2'( TT2/TT3 )1/2 TEMPR ( RTEM P - TEMPR ) CAPB -2. RTEMP. TEMPR CAPA CAPB + TT2 CAPA-CA P B + TT3 INPONT - TT1 (2 C (AC) - INPONT C ( TT3 ) I/2 T I ( RETURN

INPONT ALPS- -ICAPA I CAPBI M-a( 1 - ALP) (1 + ALP) KM —- X AM' - i=l 4 ELOG (M) * i~ Bi M' C (AC)-4*KM I (CA PA + CAPB)1 RETURN

PAGE 1 UNIVERSITY OF MICHIA PROGRA? LE\CT1 (CCTAL) 522C7 PRCGRAM COPMCN LENGTF CCCCC LOWEST ERASABLE DEFINED OCCOO INITIAL TRANSFER LECTCR CCCCC TRANSFER VECTCR LENCTF CCCI1 FIRST LOCATION EXECUTED 00011 EXTERNAL RCUTINES CALLED EY TIIS PROGRAM CCCCC 6225(3254626 Cc SETECF CCCCI 335125212460 CC.REAC CCCC2 33CCC 1C3CCCC CC.C13C0 CCCO3 334751314563 CC.PRINT CCCC4 24217C633144 CC CAYTIM CCCC5 625C51636C6C CC SQRT CCCC6 234662606060 CC COS CCCC7 254346276060 CC ELCG CCC1C 25515146516C Cc ERRCR PR00RAfP LISTINC (NUMBERS IN CCTAL) BRIEF CN PMC CFF CCC11 O074 CC 4 CCCCC 010 CALL SETECFEND CCCL3 OC74 CC 4 CCCCI CiC START READ FIN REAC PARAMETERS AND STARTING POINTS CCC15 -1 CCCCC C C2670 CiC [OP V VELOCITY CCC16 -1 CCCCC C C2671 C1C ICP PHI PRCBE POTENTIAL CCC17 -1 CCCCC C 02672 C1C lOP DELTAT TIME INCREMENT CCC20 -1 CCCCC C C2673 C1C IDP A PROBE RADIUS CCC21 C56C CC C C2673 CIC LCQ A 0CC22 C26C CC C C2673 C1C FMP A CCC23 06C1 CC C C2674 C1C STO A2 CCC24 C560 CC C 02670 C01 LDQ V 000CCC25 C26C CC C C2672 C1C FMP DELTAT CCC26 0601 CC C 02676 C1C STO VCT CCC27 0131 CC C CCCCC CC XCA CCC30 C260 CC C 02616 CIC FMP VCT CC031 0601 CC C C2675 C1C STO V2 C0C32 0560 CC C 02671 010 LDQ PHI CCC33 C260 CC C 02672 010 FMP DELTAT 00034 0601 CC C C270C C01 STO K OCC35 -1 CCCCC C 02701 01C IOP ZMIvN MINIPUM Z COORDINATE CCC0 6 -1 0C(00 C 2~703 010C lP H FRACTION OF TOP OF TRAJECTORY SPACING FOR PCINTS FOR FINDING GRADIENTS CCC37 -1 CCCCC C C2717 010 lOP NLMv NLMBER OF TRAJECTORIES CCC4C -1 000C C C0272C C01C IP NIVAX MAXIPUM NUMBER OF POINTS PER TRAJECTORY CCC41 -1 00CCC C C2721 01C lOP NITER NUMBER OF INTERATIONS TC BE MADE CCC42 -1 CCC00 C 02641 C1C IP Psh CCC43 05CO CC C 272C C10 CLA NMAX CCC44 04CC CC C 52206 C1C ADD =1 00045 0601 CC C C272C C10 STO NMAX CCC46 0534 CC 1 0272C C01 LXA NMAX,1 SET CCC47 -0634 CC 1 CC222 010 SXD S5,1 UP CCC5C -0634 CC 1 CC241 C1C SXD S61 LOCCP CCC51 -0634 CC 1 CC373 010 SXD S7,1 TERMINATIONS CCC0 2 -0634 CC 1 CC6C6 010 SXD SB,1 AND CCC53 -0634 CC 1 CC727 010 SxD S1 INCREMENTS CCC54 -0634 CC 1 CC744 010 SXD L2CACtl CCC55 -0634 CC 1 01513 010 SXD S31,1 CCC56 -0634 CC I C1224 010 SXD S1C,1 CCC57 -0634 CC 1 C0252 010 SXD S11,1

PAfE 2 UNIVERSITY OF MICFIGAN CCC60 -0634 CC 1 CC072 C1C SXD S17,1 CCC61 -C634 CG 1 CCC75 C1C SXD sie,I CCC62 -0634 CC 1 CCC71 C1C SxC S16,1 CCC63 1 77776 1 CCC64 CIC TXI *+1,1,-2 CCC64 -0634 CC 1 CC364 C1C SXD S12,1 CCC65 -0634 CC 1 C0577 C1C SXC S13,1 OCC66 -0634 CC 1 C1215 C1C SXO S14,1 CCC67 -0634 CC 1 C1504 C1C 3XD SXD S15,1 CCC70 1 CCCC2 1 CCO71 C1C TXI *+1,1,2 CCC71 1 CCCCC 1 CC072 CIC S16 TXI *+1,1,** **=PhAX CCC72 I CCCCC 1 CCC73 C10 S17 TXI *+1,1,'* *'=NMAX CCC73 -0634 CC 1 C1242 C1C SXO S19,1 CCC74 -C0634 CC 1 C1531 C1C SXD S2C,1 CCC75 1 CCCCO 1 CC076 ClC S18 TXI *+1,1,* *'=hkMAX CCC76 -C634 CC 1 CC613 C1C SXD S21,1 CCC77 -0634 CC 1 CC4CC ClC SXO S22,1 CCICC 0535 CC 1 C272C C10 LAC NMAX,1 AC~RESS OCICI -0634 CC I CC125 C1C SXD S23,1'COIFICATION CC1C2 -0634 CC I CC142 C10 SXD S24,1 CC1C3 1 C7636 1 CC104 IClC TXI *+I,IRSA CC1C4 0634 CC 1 C0437 C1C SXA S25,1 0C1C5 0634 CC 1 C1625 C10C SXA S27,1 CC1C6 0634 CC 1 CC733 C1C SXA L2CA,1 CC1C7 0634 CC 1 CC736 C1C SXA S29,1 CCllC 1 11612 1 COlll C1C TXI *+1l,,RZA-RSA CCII0 C634 CC 1 CC652 C1G SXA S26,1 CC00112 C634 CC I C1704 ClC SXA S28,1 CC113 1 73C73 1 CC114 CO1C TXI *+1,1,ZSA-RZA OCN 0CC114 C634 CC 1 C044C C1C SXA S25+1,1 OC115 0634 CC 1 C1626 C1C SXA S27+1,1 CC116 1 11612 1 CC117 C1C TXI *+1,1,ZZA-ZSA CC117 0634 CC 1 CC653 ClC SXA S26+1,1 OC12C 0634 CC 1 017C5 C10 SXA S28+1,1 CC121 1 C47C5 1 CC122 C1C TXI *+I,I,RSB-ZZA OC122 0634 CC 1 CG755 C1C SXA L2CC,1 CC123 1 1161C 1 CC124 C1C TXI *+1,1,RZB-RSB OC124 0634 CC 1 C1246 C1C SXA L30C,1 OC125 1 CCCCC 1 CC126 ClC S23 TXI *+1,1,** **=-NPAX GC126 1 42544 1 CC127 C1C TXI *+1,1,RSA-RZB CC127 C634 CC 1 CC444 C1C SXA S25A,1 CC13C 0634 CC 1 C1632 C1C SXA S27A,1 CC131 1 11612 1 C0132 C1C TXI *+1,1,RZA-RSA CC132 0634 CC 1 CC657 C1C SXA S26A,1 CC133 0634 CC 1 C1711 C1C SXA S28A,1 C00134 1 73C73 1 CC135 C1C TXI *+1,1,ZSA-RZA CC135 C634 CC 1 CC445 CLC SXA S25A+1,1 CC136 C634 CC 1 C1633 C1C SXA S27A+1,1 0C137 1 11612 1 CC14C C1C TXI *+1,1,ZZA-ZSA CC14C 0634 CC 1 CC66C C1C SXA S26A+1,1 CC141 0634 CC 1 C1712 ClC SXA S28A+1,1 CC142 1 CCCCC 1 CC143 C1C S24 TXI *+1,1,.* **=-NiAX CC143 1 61261 1 CC144 C1C TXI *+1,1,RSA-ZZA 0C144 C634 CC 1 CC451 C1C SXA S258,1 CC145 0634 CC 1 C1637 C1C SXA S278,1 CC146 I 1Ili2 1 CC147 C0C TXI *+1,1,RZA-RSA CC147 0634 CC 1 CC664 C1C SXA S26B,1 OC15C 0634 CC 1 C1716 C1C SXA S288,1 CC151 1 73C03 1 CC152 C1C TXI *+1,1,ZSA-RZA

PAGE 3 UNIVERSITY OF MIC-IGAN OC152 C634 CC 1 CC452 C1C SXA S258+1,1 0C153 0634 CC 1 C164C C1C SXA S27B+1,1 CC154 1 11612 1 CC1SS C1C TXI *+1,1,ZZA-ZSA CC155 C634 CO 1 CC665 C1C SXA S268+1,1 CC156 CE34 CC 1 C1717 C1C SXA S2Ee+1,1 CC157 C5CC CC C C2717 C1C CLA NLt CC16C C621 CC C CC412 C1C STA L5A CC161 C621 CC C CC625 C1C STA L15A CC162 C621 CC C C1577 C C STA L25A CC163 C621 CC C C1656 CIC STA L35A CC164 04CC CC C 522C5 C1C ACD =3 CCMFLTE C0165 0131 CC C CCCCC CC XCA TCTAL CC166 -C754 CC C CCCCC CC ZAC PAGES PER CC167 0221 CC C 522C4 C1C ODP =4 TRAJECTCRY CC17C -C6CC CC C C2725 C1C STQ TP VATRIX CC171 -0634 CC C CC371 C10 SXO L3C,C CC172 -C634 CC C CC604 C1C SXD L13C,C 0C173 0774 CC 1 C47C4 CC AXT SIZE,1 ZERC CC174 06CC CC 1 C7636 C1C STZ RSA,1 CIT OC175 C6CC CC 1 14543 C1C STZ ZSA,1 TRAJECTCRY CC176 C6CC CC 1 2145C C1C STZ RZA,1 PATRICES CC177 C6CC CC 1 26355 C1C STZ ZZA,1 CC2CC 0600CC CC 1 33262 C1C STZ RSB,1 0C2C1 C6CC CC 1 40166 C1C STZ ZSB,1 CC2C2 C6CO CC 1 45C72 C1C STZ RZE,1 CC2C3 06CC CC 1 51776 C1C STZ ZZB,1 CC2C4 2 CCCC1 1 CC174 C1C TIX *-8,1,1 CC2C05 C774 CC 2 CCCC1 CO AXT 1,2 CC2C6 0534 CC 1 C2717 ClC LXA NLNt,1 REAE CC2C7 -1 CCCCC 2 C7636 C1C S4 ICP RSA,2 ANC CC21C -1 CCCCC 2 14543 C1C lOP ZSA,2 SET CC211 C05CC CC 2 07636 C1C CLA RSA,2 LP CC212 C6C1 CC 2 2145C C1C SitC RZA,2 INITIAL CC213 06C1 CC 2 C7637 C1C STO RSA+1,2 OC214 06C1 CC 2 21451 C1C STO RZA+1,2 CC215 C05CC CC 2 14543 C1C CLA ZSA,2 TRAJECTCRY CC216 C601 CC 2 26355 ClC STC ZZA,2 PCINTS CC217 C3CC CC C C267.6 C1C FAD VCT CC220 06C1 CC 2 14544 C1C STO ZSA+1,2 CC221 06C1 CC 2 26356 CIC STO ZZA+1,2 CC222 1 CCCCC 2 CC223 C1C S5 TXI *+1,2,'* *=hNPAX CC223 2 CCCC1 1 CC2C7 C1C TIX S4,1,1 CC224 -1 CCCCC C CCCCC CC ENDIO CC225 2 CCCC1 2 CC226 CLC TIX *+1,2,1 CC226 -0634 CC 2 CC374 C0C SXD L4,2 C(IR2)=NIPAX*NLIv CC227 -0634 CC 2 CC457 C1C SXD L8,2 MCRE CC23C -0634 CC 2 CC607 ClC SXD L14,2 LCCP CC23. 1 -C634 CC 2 CC672 CIC SXO L18,2 SET CC022 -0634 CC 2 C1226 C1C SXC L24A,2 CP CC233 -0634 CC 2 C1547 C1C SXO L24E,2 0C234 -C634 CC 2 C1645 C1C SXD L28,2 CC235 -0634 CC 2 C1515 C1C SXD L34A,2 CC236 -0634 CC 2 C1563 C1C SXD L34E,2 CC237 -0634 CC 2 C1724 C1C SXD L38,2 CC240 -0634 CC 2 C2C060C CIC SXC LL17,2 CC00241 2 CCC-CC 2 CC242 C1C S6 TIX *+1,2,** **=NPAX CC242 -0634 CC 2 CC745 C1C SXo L2CAE,2 C(IIR2)=M'AX*(NUF-1) CC243 -0634 CC 2 C1225 C1C SXD L24,2 LCCF

PAO F 4 UNIVERSITY OF MICFIGAN CC244 -C634 CC 2 C1514 CIC Sxi) L34,2 SET LP CC245 C77'. CC 2 CCCC1 CC AXT 1,2 CC246 C754 CC 2 CCCCC CC L1 PXA,2 CCUFLTE CC247 C734 CC I CCCCC CC PAX,1 INITIAL CC25C C77i CC 4 CCCCI CC AXT 1,4 AFPRCXIIATICNS CC251 C634 CC 4 Cc36C C1C SXA L2C,4 FCR 0C252 Ce34 CC 4 CC37C C1C SXA L3B,4 SIGtVA 0C253 C560 CC 1 C7636 (1C L2 LCQ RSA,i TRAJECTCRIES CC254 C260C CC 1 C7636 CIC FFP RSA,1 CC255 06C1 CC C C2642 C1C STC T1 CC256 056C CC 1 14543 C1C LCQ ZSA,1 CC257 C26C CC I 14543 C1C FMP ZSA,1 CC260C 03CC CC C C2642 C1C FAD TI CC261 CSLC CC C 522C3 C1C LCQ =-1.5 GC262 CC74 CC 4 CCCC2 ClC CALL.C13CC EXPCNIENTATICh SUBRCUTINE CC263 C131 CC C CCCCC CC XCA CC264 0260 (C C C 7CC CHC FMP K 0C265 C6C1 CC C C2642 CLC STI T1 OC266 C5CC CC C 522C2 CIC CLA =1. CC267 C3CC CC C C2642 C1C FAD TI CC27C GCC1 CC C C2642 C1C STO T1 CC271 C131 CC C CCOCC CC XCA CC272 C26u CC 1 C7636 C1C FMP RSA,1 CC273 C6Cl CC I C7635 ClC SIC RSA-1,1 CC274 -C12C CC C CC367 C1C TvI L3A CC275 C560 CC C C2642 C1C LCO T1 PARTICLLAR 00C276 C260C CC 1 14543 C1C FMP ZSA,1 TRAJECTCRY CC277 03C2 CC C C2676 C1C FSB VCT hILL BE CC3CC O6C1 CC 1 14542 ClC SO0 ZSA-I,1 TERPINATEC CC3C1 04C02 CC C C27C01 C1C SUB ZFIn IF R VALLE CC3C2 -C12C CC C CC366 CIC TMI L3 NEGATIVE CC3C3 CSCC CC 1 C7636 C1C CLA RSA,1 CC3C4 -CCO CC CCCCC C0 C TNZ *+4 CC3C5 C5CC CC 1 14542 CIC CLA ZSA-1,1 CC3C6 C3C2 CC 1 14543 C1C FSB ZSA,I CC3C7 012C CC C CC366 C1C TPL L3 CC31C 0560 CC 1 C7635 C1C LC( RSA-I,1 CC311 C260 CC 1 C7635 C1C FPP RSA-i,1 CC312 C6C1 C C C 2711 CLC STIC TEP CC313 C560 CC 1 14542 C1C LCQ ZSA-I,1 CC314 C260C CC 1 14542 C1C FVP ZSA-I,1 CC315 C3CC CC C C0711 C1 FAC TEMP CC316 C4Cc CC C C2674 C1C SUB A2 CC317 -0120 C C CC366 C1C TMI L3 CC32C CSCC CC 1 14543 C1C CLA ZSA,1 CC321 CC02 CO 1 14542 CIC FSB ZSA-1,1 CC322 C-C! CC C C,711 CIL STO TFIVP CC323 05CC CC 1 C7636 C1C CLA RSA,1 CC324 03C2 CC 1 C7635 CLC FSR RSA-1,1 CC325 C241 CC C C2711 C1C FCP TE P CC326-CCC CC C C2726 -C6CC CC C 02726 010 STQ B CC327 C260 CC C C2726 C C FMP B CC33C UCC1 CC C C2727 ClC ST RB2 CC331 C56C CC 1 14543 C1C LCQ ZSA,1 CC332 C2b0 CC C C2726 CiC FP'P P CC333 00C7 CC 1 C7636 CLC FSO RSA,1 CC334 06C1 CC C C273C C1C STO 1 CC335 0131 CC C CCCCC CC XCA

PAGE 5 UNIVERSITY OF MICFIGAN CC336 C26C CC C C2730 ClC FVP 81 CC337 06C1 CC C C2731 C1C STU B12 CC34C C3C2 CC C C2674 C1C FSB A2 CC341 0241 CC C C2731 ClC FCP 812 CC342 -06CC CC C C2711 C1C STC TEMP CC343 05CC CC C 522C2 C 1C CLA =1. CC344 0241 CC C C2127 C1C FCP B2 CC345 0131 CC C (CCCC CC XCA CC346 C3CC CC C 52C2 CIC FAD =1. CC347 0131 CC C CCCCC CC XCA CC35C C26C CC C C2711 C1C FVP TEMP CC351 04C2 CC C 522C2 C1C SLB =1. CC352 ClCC CC C CC366 CIC TZE L3 CC353 C12C CC C CC357 C1C TPL L28 CC354 C560 CC 1 14543 CIC LC0 ZSA,1 CC355 026C CC 1 14542 ClC FMP ZSA-1,1 CC356 -Clfe CC C CC366 CiC TMI L3 0C357 1 CCCC1 1 CC36C C1C L2e TXI *+41,1, CC36C C774 CC 4 CCCCC CC L2C AXT **,4 *:=INTERPECIATE IMAXSA CC361 1 CCCC1 4 CC362 C1C TXI *+1,4,1 CC362 C634 CC 4 C03.C C1C SXA L2C,4 CC363 C634 CO 4 CC37C C1C SXA L38,4 CC364 -3 CCCCC 4 CC253 C1C S12 TXL L2,4,** **:=fAX-1 CC365 OC2C CC C CC371 CLC TRA L3C CC366 C6CC CC 1 14542 C1C L3 STZ ZSA-1,1 CC367 C6C0 CC 1 C7635 C1C L3A STZ RSA-1,1 CC370C 774 CC 4 CCCCC CC L38 AXT **,4 *=INIERMECIATE IMAXSA 0 CC371 -3 CCCCC 4 CC373 CLC L3C TXL *+2,4,4* I*=IFAXSA CC372 -0634 CO 4 CC371 C1C SXD L3C,4 CC373 1 CCCCC 2 CC374 C1C S7 TXI *+1,2,** *+1=hPiAX CC374 -3 CCCCC 2 CC246 C1C L4 TXL L1,2,** **=hNAX*NUM CC375 -0534 CC 4 CC371 ClC LXD L3C,4 OC376 0634 CC 4 CC4C6 C1C SXA L5,4 CC377 -0535 CO 4 CC371 C1C LCC L3C,4 OC4CC 1 CCCCC 4 CC4C1 C1C 522 TXI *+1,4,,* *:=411*lIAX CC4C01 -0634 CC 4 CC456 C1C SX0 L7,4 CC4C2 C0500CC CC 0 522C4 C1C CLA =4 CC4C3 0601 CC C C2723 C0C STC VAR FCRPAT VARIAPLE CC4C4 0774 CC 2 CCCC1 CC AXT 1,2 CC4C5 UeCC CC C C2724 C1C STZ P C04C6 C774 CC 1 CCCCC CC L5 AXT **,1 *.=IIdAXSA CC4C7 05CO CC C C2724 CIC CLA P PAGE hNLBER CC41C 04CC CC C 522C6 C1C ACD =1 CC411 06C1 CC 0 C2724 C0lC STO P 00412 0774 CO 4 CCOCC CC L5A AXT **,4 *=S[IGFA TRAJECTORIES STILL TO BE PRINTED CC413 2 CCCC4 4 CC415 C'1C TIX *+2,4,4 CC414 C0634 CC 4 C2723 C1C SXA VAR,4 C0415 C634 CC 4 CC0412 C1C SXA L5A,4 00416 CC74 CO 4 CC0C3 ClC PRINT FVAR,...,FCLTI,P,TFV,PHI,CELTAT,A,H,ZMIN CC43C C5C0 CC C C2723 C1C CLA VAR 0C431 O1CC CC C CC453 C1C TZE L7-3 0C432 -1 CCCCC 2 C7636 C1C L6 ICP RSA,2 PRINT INITIAL CC433 -1 CCCCC 2 14543 10C ICP ZSA,2 SIGIA TRAJECTCRIES CC434 5C00 CC C C2723 C1C CLA VAR CC435 04C2 CC C 52206 C1C SUB =1 CC436, 01CC CC C C0453 C1C TZE L7-3 CC437 -1 C0CCC 2 0CCCC CC S25 ICP **,2 **=RSA-NMAX CC44C -1 CCCCC 2 CCCCC CC ICP **,2 **=ZSA-NMAX

PAGt 6 UNIVERSITY OF MICHIGAN OC441 05CC CC C C2723 C1C CLA VAR CC442 04C2 CC C 522C1 C1C SUB =2 CC443 O1CC CC C CC453 C1C TZE L7-3 CC444 -1 CCCCC 2 CCCCC CC S25A ICP **,2: *=RSA-2*N'AX CC445 -1 CCCCO 2 CCCCC CC ICP *4,2 **=ZSA-2 V'NFAX CC446 0500C CC C C2723 C C CLA VAR CC447 04C2 CC C 522C5 C1C SUB =3 CC45C ClCC CC C C0453 C1C TZE L7-3 00451 -i CCCCC 2 CCCCC CC S258 IGP **,2 *,=RSA-3'NMAX CC452 -1 CCCCC 2 CCCCC CC lOP **,2 **=ZSA-3*aNAX CC453 1 C~CC1 2 CC454 CIC TXI *+1,2,1 00C454 2 CCCC1 I C0432 C1C TIX L6,1,1l CC455 -1 CCCCC C CCCCC CC ENDIO CC456 1 CCCCC 2 CC457 C1C L7 TXI *+1,2,** *:=4hNMAX-IFAXSA CC457 -3 CCCCC 2 CC406 C1C L8 TXL L5,2,** **=hlAX*NUM CC46C 0774 CC 2 CCO01 CC AXT 1,2 CC00461 0754 CC 2 CCCCC CC L11 PXA,2 OC462 C734 CC 1 CCCCC CC PAX,1 CCMFLTE CC463 0774 CC 4 CCOC1 CC AXT 1,4 INITIAL 0C464 0634 CO 4 C0573 CLC SXA L12C,4 ZETA TRAJECTCRIES CC465 C634 CC 4 CC6C3 C1C SXA L138,4 CC466 056C CC 1 2145C C1C L12 LCQ RZA,1 CC46.7 C26C CC 1 2145C C1C FMP RZA,1 OC470 06C1 CC C C2642 C1C STO T1 CC471 C560 CC 1 26355 C1C LCQ ZZA,1 CC472 C260 CC 1 26355 C1C FMP ZZA,1 CC473 C3CC CC C C2642 C1C FAD T1 0C474 C560 CC 0 522C3 C1C LCQ =-1.5 o CC475 CC74 CC 4 CCCC2 C1C CALL.C13CC CC476 0131 CC C CCCCC CC XCA OC477 C26C CC C C27CC CLC FMP K OC5CC 06C1 CC C C2642 C C STO TI C0050C1 05CC CC C 522C2 C1C CLA =1. CC5C2 03C2 CC C C2642 ClC FSB T1 OC5C3 06C 1 CC C C2642 C1C STO T1 C0050C4 0131 CC C CCOCC CC XCA OC5C5 C26C CC 1 2145C C1C FMP RZA,1 CC5C6 06C1 CC 1 21447 C1C STO RZA-l,l OC5C7 -012C CC 0 CC6C2 C1C TVI L13A CC0051C 056C CC C C2642 C1C LCQ TI OC511 C26C CC 1 26355 C1C FMP ZZA,1 CC512 03C2 CC C C2676 C1C FSB VCT OC513 06C1 CC 1 26354 C1C STO ZZA-1,1 C00514 04C2 CC C C27C1 C1C SUB ZfINh CC515 -C12C CC C L'C6C1 C1C TMI L13 CC516 05CC CC 1 2145C CLC CLA RZA,1 OC517 -C1CC CC C CC523 CLC TNZ *+4 CC520 C5CC CC 1 26354 C C CLA ZZA-1,1 OC521 03C2 CC 1 26355 C1C FSB ZZA,1 CC522 0120 CC C CC6C1 CIC TPL L13 CC523 C560C CC 1 21447 C1C LOQ RZA-i,1 CC524 C26C CC 1 21447 C1C FMP RZA-1,1 CC525 06C1 CC C C2711 C1C STO TEMP CC526 C56C CC 1 26354 CIC LCQ ZZA-1,1 00C527 C26C CC 1 26354 C1C FMP ZZA-1,1 0CC53C 3CC C-C C C2711 C 1C FAD TEMP CC531 04C2 CC C C2674 C1C SLB A2 CC532 -C120C CC C CC6C1 C 1C TI L13

PAGE 7 UNIVERSITY OF MICFIGAN CC533 CSC CC 1 26355 C1C CLA ZZA,1 CC534 03C2 CC 1 26354 C1C FSB ZZA-1,l CC535 C601 CC C C2711 C1C STO TEMP CC536 CSCC CC 1 2145C C1C CLA RZA,1 CC537 03C2 CC i 21447 C1C FSB RZA-I,1 CC54C 0241 CC C C2711 C1C FOP TEMP 00541 -06CO CC C C2726 ClC STC B CC542 C26C CC C C2726 C1C FPMP B CC543 C6Cl CC C C2727 C1C STO B2 0C544 056C CC 1 26355 C1C LDQ ZZA,1 CC545 0260 CC C C2726 C1C FMP B CC546 03C2 CC 1 21450 C10 FSB RZP,1 CC547 060C1 CC C C2730 C1C STO BI CC55C 0131 CC C CCCCC CC XCA OC051 C260 CC C C2730 C1C FMP B1 CC552 C6CI CC C C2731 ClC STO B12 0C553 03C2 CC C C2674 t1C FSB A2 OC554 0241 CC C C2731 CIC FCP B12 00555 -0600 CC C C2711 C1C STQ TEMP CC556 C5CC CC 0 522C2 ClC CLA =1. CC557 0241 CC C C2727 C1C FCP B2 CC560 C131 CC C CCCCC CC XCA 0C561 C3CC CCC 522C2 C1C FAD =1. CC562 0131 CC C CCCCC CC XCA CC563 0260 CC C C2711 C1C FMP TEMP OC564 04C2 CC C 522C2 C1C SLB =1. 0C565 C1CO CC C CC6C1 ClC TZE L13 CC566 0120 CC C CC572 C1C TPL L12B CC567 C560 CC 1 26355 C1C LCO ZZA,1 OC57C C260 CC 1 26354 C1C FMP ZZA-1,1 CC571 -C12C CC C CC6C1 C1C TMI L13 OC572 1 CCCC1 1 CC573 C1C L12P TXI *+1,1,1 CC573 0774 CC 4 CCOCC CC L12C AXT **,4 ~*=INTEREDIIATE IMAXZA CC574 1 CCCC1 4 CC575 C1C TXI *+1,4,1 CC575 C0634 CC 4 CC573 C1C SXA L12C,4 CC576 C634 CC 4 CC6C3 C1C SXA L13B,4 CC517 -3 CCCCC 4 CC466 C1C S13 TXL L12,4,** **=hNA X-1 CC6CC CC2C CC C CC6C4 C1C TRA L13C CC6C1 06CC CC 1 21447 C1C L13 STZ RZA-i,l CC6C2 06CC CC 1 26354 C1C L13A STZ ZZA-1,l CC6C3 G774 LC 4 C0CCC CC L13E AXT **,4 **=INTERMECIATE IMAXZA CC6C4 -3 CCCCC 4 CC606 ClC L13C TXL *+2,4,**'*=IPAXZA CC065 -C634 CC 4 CC604 C0C SXD L13C,4 CC60C6 1 CCCCC 2 CC6C7 C1C Se TXI *+1,2,''* *a=NAX CC6C7 -3 CCCCC 2 CC461 ClC L14 TXL L1,2,** **=N:AX*NUM CC61C -0534 CC 4 CC604 C1C LXD L13C,4 C0611 0634 CO 4 CC621 CLC SXA L15,4 CC612 -0535 CC 4 CC604 C1C LCC L13C,4 CC613 1 CC(CC 4 CC614 C1C S21 TXI *+1,4,*4 **=4hMNAX CC614 -0634 CC 4 CC671 C1C SXD L17,4 CC15 C05CO CC C 522C4 C0C CLA =4 CC616 C06C CC C C2723 ClC STO VAR CC617 0774 CC 2 CCCC1 CC AXT 1,2 00C62C 06C CC C C2724 C1C STZ P CC621 0774 CC 1 CCCCC CC L15 AXT **,1 **=IIPXZA OC(22 05CC CC C C2724 ClC CLA P CC623 C4CC CC C 522C6 ClC AED =1 CC624 0CCl CC C C2724 C1C STO P

PAGE 8 UNIVERSITY OF MICI-IGAN 00625 C774 CC 4 CCOCC CC L15A AXT **,4 **=ZETA TRAJECTCRIES STILL TO BE PRINTED 0CC26 2 CCCC4 4 CC630 ClC TIX *+2,4,4 CC627 0634 CC 4 C2723 C1C SXA VAR,4 CC63C C634 CC 4 CC625 C1C SXA L15A,4 OC631 C0074 CC 4 CCC03 C1C PRINT FVAR,...,FCLT2,F,TP,V,PHI,CELTAT,A,I-,ZPIN CC0643 05CO CC O C2723 C C CLA VAR 0C644 01CC CC C CC.66 C1C TZE L17-3 CC645 -1 CCCCC 2 21450C C1C L1 IOP RZA,2 CC646 -1 CCCCC 2 26395 C1C IGP ZZA,2 OC647 C5CO CC C C2723 C C CLA VAR OC65C 04C2 CC C 522C6 C1C SUB =1 CC651 01CC CC C CC666 C1C TZE L17-3 CCt52 -1 CCCCC 2 CCOCG CC S26 ICP **,2 **=RZA-NMAX OC653 -1 CCCCC 2 CCCCC CC ICP **,2 *:=ZZA-hVIAX CC654 05CC CC C C2723 C1C CLA VAR CC655 04C2 CC 0 52201 C1C SUB =2 CC00656 O1CC CC C CC666 C1C TZE L17-3 00C657 -1 CCCCC 2 CC0CC CC S26A ICP **,2 **=RZA-2N1'AX OC6C -1 CCCCC 2 CCCCC CO IOP **,2 **=ZZA-2INMAX CC661 0500CC CC C C2723 CIC CLA VAR CC00662 04C2 CC C 522C5 C 1C SUB =3 C663 C010CC CC C CC666 C1C TZE L17-3 CC664 -1 CCCCC 2 CCCOC CC S268 IOP **,2 **=RZA-3lNMAX CC665 -1 CCCCC 2 CCCCC CC lIP **,2 **=ZZA-3lN\PAX CC666 1 CCCCI 2 CC667 C1C TXI *+1,2,1 CC667 2 CCCC1 1 CC645 C10C TIX L16,1,1 CC670C -1 CCCCC C CCCCC CC ENDIO CC671 1 CCCCC 2 CC672 ClC L17 TXI *+1,2,* *.*=4*PAX-IMAXZA CC672 -3 CCCCO 2 CCL21 C10 L18 TXL L15,2,9* *=:AX*NLM OC673 06CC CC C C2714 C1C STZ CCLNT CC674 -0534 CC 4 CC371 CIC LXD L3C,4 CC675 -0634 CC 4 C1231 C1C SXD L24AA,4 CC676 -0534 CO 4 CC604 C1C LXC L13C,4 CC677 -0634 CC 4 C152C C1C SXO L34AA,4 CC7CC 05CC CC C C C2714 C1C L2C CLA CCLNT CC7C1 04C2 CC C C2721 C1C SUB NITER 007C2 01CC CC C CC013 C1C TZE START CC7C3 -0634 CC C C1222 C1C SXD L23C,C CC7C4 -C634 CC C C1511 C1C SXC L33C,C CC7C5 C774 CC 2 CC001 CC AXT 1,2 CC7C6 C534 CC 1 C2717 C1C LXA NLM,1 CC7C7 C5C0 CC 2 C7636 C1C Se.S CLA RSA,2 CC71C C601 CC 2 23262 CIC STO RSB,2 PCINIS CC711 05CC CC 2 14543 C1C CLA ZSA,2 FCR NEXT OC712 060C1 CC 2 4C166 C1C STO ZSB,2 APPRCXIPATICN CC713 05CO CC 2 2145C C1C CLA RZA,2 TC TRAJECTCRIES CC714 C6Cl CC 2 45C72 ClC STO RZB,2 C0715 05CC CC 2 26355 C1C CLA ZZA,2 OC716 06C01 CC 2 51776 C1C STO ZZB,2 CC717 C5CO CC 2 C7637 C1C CLA RSA+1,2 CC72C 06C1 CC 2 33263 CC STO RSB+1,2 CC721 05CC CC 2 14544 C1C CLA ZSA+1,2 00C722 C6C1 CC 2 4C167 C1C STO ZSB+1,2 CC723 CSCC CC 2 21451 C1C CLA RZA+1,2 CC724 06C01 CC 2 45C73 C1C SC RZB+1,2 CC725 C050CG CC 2 26356 C1C CLA ZZA+1,2 CC726 06C1 CO 2 51777 C10 STO ZZB+1,2 CC727 1 CCCCC 2 CC73C C1C SS TXI *+1,2,* **I=hIAX

PAGE S UNIVERSITY CF MICFIGAN CC73C 2 CCCC1 1 CC7C7 C1C TIX SE.S,1,l CC731 0774 CC 1 CCCCl CC AXT 1, 1 CC732 0774 CC 2 CCCC1 CC AXT 1,2 CC733 C5CC CC 2 CCCCC CC L2CA CLA **,2 **=RSA-NMAX CC734 C C2 CC 2 C7636 CIC FSB RSA,2 CC735 C601 CC C C27C4 C 1C ST0 DA CC736 C5CC CC 2 CCCCC CC S2S CLA **,2 **=RSA-NVAX CC737 C3CC CC 2 C7636 CLC FAD RSA,2 CC74C C131 CC C CCCCC CC XCA CC741 C26C CC C C27C4 CIC FMP DA CC742 06C1 CC 1 52166 C1C STO VEC,1 CC743 1 CCCC1 1 C0744 CLC TXI *+1,1,1 CC744 1 CCCCC 2 CC745 ClC L2CAC TXI *+1,2,'* **=hPAX CC745 -3 CCCCC 2 CC733 C1C L2CAE TXL L2CA,2,** *=NIPAX*(NKUF- 1) CC746 05CO CC 2 C7636 C1C CLA RSA,2 CC747 04CC CC C 522CC C1C ACD =K1KS CC75C CSCC CC C C27C4 C1C FAD DA CC751 0131 CC C CCCCC CC XCA CC752 C26C CC C C2704 C1C FIP CA CC753 06C1 CC 1 52166 CLC STO VEC,1 CC754 0774 CC 2 CCOC1 CC AXT 1,2 CC755 05CC CC 2 CCCCC CC L2CC CLA **,2 e=PRSB-NMAX CC756 03C2 CC 2 33262 C1C FSB RSB,2 CCIFLTE CC757 0131 CC C CCCCC CC XCA NEXT CC76C C260 CC C C27C3 C1C FMP H APPRCXIMATI ChS CC761 06C1 CC C C2677 C1C STO CELTA TC TRAJECTCRIES CC762 0754 CC 2 CCOCC CC L21 PXA,2 0 CC763 C734 CC 1 CCCCC CC PAX,1 CC764 C774 CC 4 CCCC1 CC AXT 1,4 CC765 0634 CC 4 C1211 C1IC SXA L22C,4 CC766 0634 CC 4 C1221 C1C SXA L23B,4 CC767 -C520 CC C C2641 CLC L22 NZT PSh CC77C CC2C CC C CICC1 C1C TRA L22.1 OC771 CC74 CC 4 CCCC4 C1C CALL CAYTIM CC772 C6C1 CC C C264C CLC STO X2 CC773 CC74 CC 4 CCCC3 CLC PRINT FL CC775 -1 CCCCC 1 33262 C1C ICP RSe,1 CC776 -1 CCCCC 1 4C166 CLC ICP ZSB,1 CC777 -1 CCCCC C C264C CLC IOP X2 O1CCC -1 CCLCL C CCCCC CC ENCIO O1CCl CSCC CC 1 32262 CLC L22.1 CLA RSB,l C1CC2 O1CC CC C ClC40 C1C TZE L22.2 C1CC3 C3C2 CC C C2677 CLC FSB DELTA C1CC4 CeC2 CC C C2F6C CLC SLW RTEMP 01CC5 CSCC CC 1 4C166 C1C CLA ZSB,1 C1CC6 C6Cl CC C C2661 C1C STC ZTEVP C1CC7 CC74 CC 4 C1726 C1C TSX PCTNTL,4 OIClC 3 CCCCC C 2145C CIC PAR RZA ClCll 3 CCCCC C 26355 ClC PAR ZZA C01C12 OCLL CC C C2642 C1C STC T1 ClCl3 CC74 CC 4 C1726 C1C TSX PCTNTL,4 01C14 3 CCCCC C C7636 CIC PAR RSA C1C15 3 CCCCC C 14543 C1C PAR ZSA C1C16 C3C2 CC C C2642 C1C FSB T1 01C17 C6C1 CC C C2643 ClC SIG T2 CLC2C C5CC CC 1 33262 C1IC CLA RSB,1 CIC21 03CC (C C C2677 CIC FAD CELTA 01C22 06C1 CC C C266C CLC SIC RTEVP

PAGE 10 UNIVERSITY OF MICHIGAN C1C23 CC74 CO 4 C1726 C1C TSX PCTINITL,4 C1C024 3 CCCCC C 2145C C1C PAR RZA 01C25 3 CCCCC C 26355 CIC PAR ZZA C01C26 06C1 CC C C2642 C1C STO T1 C01C27 GC74 CC 4 C1726 C1C TSX PCThTL,4 ClC30 3 CCCCC C C7636 C1C PAR RSA 01C31 3 CCCCC C 14543 C1C PAR ZSA C1C32 03C2 CC C C2642 C1C FSB T1 C1C33 03C2 CC C C2643 C1C FSB T2 010C34 0241 CC C C2677 C1C FCP DELTA 01C35 026C CC C 52177 CIC FMP =.C3S7EE7 =.5/(4*PI) 01C36 C76C CC C CCCC2 CC CHS C1C037 03CC CC 1 33262 C1C FAC RSB,1 C1040 06C1 CC C C2666 CIC L22.2 STO RGP 01041 05CC CO 1 33262 ClC CLA RSB,1 01C42 06C1 CC C C260C C1C STO RTEIP 01C43 C5CC CC 1 46166 C1C CLA ZSB,1 01C44 C3C2 CC C C2677 C1C FSB DELTA 01C45 060C1 CC C C2661 CIC STO ZTEMP C1C46 0C74 CC 4 C1726 C1C TSX PCINTL,4 C1C47 3 0CCCC C 2145C C1C PAR RZA.1CC5C 3 CCCCC C 26355 C1C PAR ZZA 01C051 C6C1 CC C C2642 C1C STO Tl 01C052 CC74 CC 4 C1726 C1C TSX PCTNTL,4 01C53 3 CCCCC C C7636 C1C PAR RSA 01C54 3 CCCCCC 14543 C1C PAR ZSA 01C055 03C2 CC 0 C2642 C1C FSB T1 01C56 06C1 CC C C2643 C1C STO T2 C1C57 05C0 CC 1 4C166 C1C CLA ZSB,1 01C60 03CC CC C C2677 C1C FAD DELTA C1C61 06C1 CC C C2661 C1C STO ZTEFP 01C62 00C74 CC 4 C1726 C1C TSX PCTNTL,4 01063 3 CCCCC C 21450C CC PAR RZA 01C064 3 CCCCC C 26355 C1C PAR ZZA 01C65 C6C0 CC C C2642 C1C STO Ti 01C66 0C74 CC 4 C1726 C1C TSX PCTNTL,4 C1C67 3 CCCCC C C7636 C1C PAR RSA C107C 3 CCCCC C 14543 C1C PAR ZSA 01071 03C2 CC C C2642 C1C FSB T1 01C72 03C2 CC C C2643 C1C FSB T2 C1C73 0241 CC C C2677 ClC FCP CELIA C1C74 C26C CC C 52177 C1C FVP =.03S7EE7 01C75 076C CC C CCOC2 CC CHS C1C76 03C00 CC 1 4C0166 C1C FAD ZSB,1 C1C77 06C1 CC C C2667 C1C STO ZGP 01lCC 056C CC 1 ~3262 C1C LCQ RSB,1 01lC1 C260 CC 1 33262 C1C FMP RSe,1 011C2 0601 CC C C2642 C1C STO T1 011C3 0560 CC 1 4C166 C1C LCQ ZSB,1 01llC4 C260 CC 1 40166 C1C FMP ZSB,1 011C5 03CC CC C C2642 C1C FAD T1 CllC6 06C1 CC C C2642 C1C STO T1 011C7 0131 CC C CCCCC CC XCA C111C C26C CC C C2642 C1C FMP T1 Cl11l 06C1 CC C C2643 C1C STO T2 C1112 0C74 CC 4 CCCC5 C1C CALL SCRl,11 C1114 C241 CC C C2643 ClC FCP T2 01115 C26C CC C C27CC C1C FMP K

PAGE 11 UNIVERSITY OF MICFIGAN 01116 06C1 CC C C2642 C1C STO T1 01117 0131 CC C CCOCC CC XCA 0112C C26C CC 1 33262 C1C FMP RSB,1 C1121 03CC CC C C2666 C1C FAC RGP 01122 06C1 CC 1 33261 C1C STO RSB-1,1 01123 ClCC CC C C1125 C1C TZE *+2 01124 -0120 CC C C122C C1C TMI L23A 01125 0560 CC C C2642 C1C LCQ T1 01126 C26C CC 1 40166 C1C FMP ZSB,1 01127 C3CC CC C C2667 C1C FAC ZGP 0113C 03C2 CC C C2676 C1C FSB VCT 01131 06C1 CC 1 4C165 C1C STO ZSB-1,1l 01132 04C2 CC C C2701 C1C SUB ZMIN C01133 -C12C CC 0 C1217 C1C TMI L23 01134 05CC CC 1 33262 C1C CLA RSB,1 C1135 -C1CC CC C C1141 C1C TNZ *+4 01136 05CC CC 1 4C165 CI1C CLA ZSB- 1,1 01137 C3C2 CC 1 4C166 C1C FSB ZSB,1 C114C 012C CC C C1217 C1C TPL L23 C1141 0560 CC 1 33261 C1C LCQ RSB-I,1 01142 0260 CC 1 33261 C1C FMP RSB-1,1 01143 0601 CC C C2711 C1C STO TEMP C1144 0560 C C 40165 C 1C LCQ ZSB-1,1 C1145 C260 CC I 4C165 C1C FMP ZSB-I,1 01146 -03CC CC C C2711 C1C FAD TEMP C1.147 C4C2 CC C C2674 C1C SUB A2 01150 -0120 CC C C1217 C1C TMI L23 01151 0500 CC 1 4C166 C1C CLA ZSB,1 01152 03C2 CC 1 4C0165 ClC FSB ZSB-l,l 01153 06C01 CC C C2711 C1C STO TEFP 01154 05CC CC 1 33262 C1C CLA RSB,1 01155 03C2 CC 1 33261 C1C FSB RSB-,1l 01156 0241 CC 0 C2711 C1C FCP TEMP 01157 -C6C0 CC C C2726 C1C STQ B 0116C 0260 CC C C2726 C1C FMP B C1161 060C1 CC C C2727 C1C STO B2 C1162 0560 CC 1 4C166 C1C LCQ ZSB,1 01163 C260 CC C C2726 C1C FMP B C1164 C3C2 CC 1 33262 C1C FSB RSB,I C1165 06C1 CC C C2730 C 1C STO B1 01166 0131 CC C CCCCC CO XCA 01167 C260 CC C C273C C1C FMP B1 C1170 06C1 CC C C2731 C1C STC B12 C1171 03C0 CC C C2674 C1C FSB A2 01172 0241 CC C C2731 C1C FCP B12 01173 -C0600CO CC C C2711 C1C STQ TEMP C1174 G500 CC C 522C2 C1C CLA =1. 01175 C241 CC C C2727 ClC FCP B2 01176 0131 CC C CCOCC CC XCA C1177 C3CC CC C 522C2 C1i FAD =1. 012CC 0131 CC C CCOCC CC XCA 012C1 C260C CC C C2711 C1C FMP TEMP 012C2 04C2 CC 0 522C2 C10C SUB =1. C12C03 C1CC CC C C1217 C1C TZE L23 012C4 0120 CC C C121C C1C TPL L228 C12C05 0560C CC 1 4C166 C1C LCQ ZSB,1 C12C6 C26C CC 1 4C165 C1C FMP ZSe-I,l C12C7 -C12C CC C C1217 C 1C TMI L23

PAGE 12 UNIVERSITY CF MICFIGAN 0121C 1 CCCC1 1 C1211 CLC L22P TXI *+1,1,1 01211 0774 CC 4 CCCCC CC L22C AXT **,4 "*=INTERVECIATE ItMAXSE C1212 i CCCCl 4 C1213 C1C TXI *+1,4,1 C1213 0634 CC 4 C1211 C1C SXA L22C,4 01214 C634 CC 4 C1221 C1C SXA L238,4 C1215 -3 CCCCC 4 CC767 C1C S14 TXL L22,4,**.=NhPAX-1 C1216 CC2C CC C C1222 C1C TRA L23C 01217 C6CC CC 1 40165 C1C L23 STZ ZS8-1,1 C122C C6CC CC 1 33261 C1C L234 STZ RSE-1,1 C1221 C774 CC 4 CCCCC CC L23B AXT **,4 *.=IITERIECIATE IMAXSB C1222 -i CCCCC 4 C1224 C1C L230C TXL *+2,4,** *.=IIPAXSE 01223 -0634 CC 4 C1222 C1C SXC L23C,4 01224 1 CCCCC 2 C1225 C10 SiC TXI *+1,2,"* *=:hPAX 01225 -3 CCCCC 2 CC755 ClC L24 TXL L2CC,2,**.*=:pAX,(NLhP-1) C1226 -3 CCCCC 2 CC762 C1C L24A TXL L21,2,u* **=tAX*KUM 01227 -C0534 CC 2 C1222 C1C LXC L23C,2 C123C -0534 CC 4 C1222 C1C LXD L23C,4 C1231 3 CC(CC 4 C1233 C1C L24AA TXH *+2,4,.* *.=FRLVICLS IVAXSB C1232 -0534 CC 4 C1231 C1C LXD L24AA,4 01233 0634 CC 4 Gi535 C1C SXA L248,4 01234 C0634 CO 4 C1571 C1C SXA L25,4 M ICCLE C1235 0754 CC 4 CCOCC CC PXA,4 01236 -C760 CC C CCCC3 CC SSM C12:37 C4C CC C C2720C CIC ACD NtAX 0124C 0734 CC 4 CCCCC CC PAX,4 C1241 -0634 CO 4 Cl546 C1C SXC L24C,4 "",] 01242 1 CCCCO 4 C1243 C C S1S TXI *+1,4,.' 4=3*1\PAX 0c1 C01243 -0634 CC 4 C1644 C1C SXC L27,4 C1244 -0t34 (t 2 C!231 C1C SXC L24AA,2 C1245 C774 CC 2 CCCO1 CC AXT 1,2 01246 05CC CC 2 CCOCC CC L3CC CLA **,2 *"=RZB-INtvAX C1247 03C2 CC 2 450C72 C1C FSB RZB,2 C1250C C131 CC C CCCCC CC XCA 01251 C26C CC C C27C3 C 1C FMP H 01252 C06C1 CC C C2677 C1C STO DELTA 01253 0754 CC 2 CCCCC CO L31 PXA,2 01254 0734 CC 1 CCCCC CC PAX,1 01255 C774 CC 4 CCCC1 CC AXT 1,4 01256 C634 CO 4 C1500CC CC SXA L32C,4 01257 C634 CC 4 C151C C1C SXA L33E,4 C1260 -052C CC C C2641 C1C L32 NZT PSh C1261 CC74 CC 4 CCC04 C1C CALL DAYTIM C1262 C601 CC C C264C C1C STIC X2 01263 CC74 CC 4 CCCC3 CLC PRINT F1 C1265 -1 CCCCC 1 45C72 C1C ICP RZB,1 C1266 -1 CCCCC 1 51776 C1C ICP ZZB,1 C1267 -1 CCCCC C C264C C1C ICP X2 0127C -1 CCCCC C CCCCC CC ENOIO C1771 0C5CC CC 1 450C72 ClC L32.1 CLA RZ,1 C1272 ClCC CC C C1327 C1C TZE L32.2 C1273 C3C2 CC C C2677 C 1C FSB CELTA C1274 C06C2 CC C C266C C1C SLW RTEP C1275 05CC CC 1 51776 ClC CLA ZZB,1 C1276 060 CC C C2661 C1C ST] ZTEIVP 01277 CC74 CC 4 C1726 C1C TSX PCTNTL,4 C13CC 3 CCCCC C 2145C C 1C PAR RZA C13C1 3 CCCCC C 26355 C1C PAR ZZA

PAGE 13 UNIVERSITY OF MICFIGAN C13C2 06C 1 CC C C2642 C C STO T1 C13C3 0C74 CC 4 C 1726 C1C TSX PCTNTL,4 013C4 3 CCCCC C C7636 CLC PAR RSA C13C5 3 CCCCC C 14543 ClC PAR ZSA 013C6 03C2 CC C C2642 C1C FSB TI C13C7 C6C1 CC C C2643 ClC STU T2 0131C C5CC CC 1 45C 72 C 1C CLA RZB,1 C1311 C3CC CC C C2677 C1C FAD DELTA 01312 06C1 CC C C2660 ClC STO RTEVP C1313 0C74 CC 4 C1726 C1C TSX PCTNTL,4 01314 3 CCCCC 0 2145C C 1C PAR RZA 01315 3 CCCCC C 26355 C1C PAR ZZA 01316 06C1 CC C C2642 C1C STO T1 01317 CC74 CC 4 C1726 CIC TSX PCINTL,4 0132C 3 CCCCC C C7636 C C PAR RSA 01321 3 CCCCC C 14543 ClC PAR ZSA C1322 03C2 CC C C2642 C1C FSB TI 01323 03C2 CCC C2643 C 1C FSB T2 01324 C241 CC C C?677 C1C FCP DELTA C1325 0260 CC C 52177 C1C F'P =.C3S7EE7 C1326 C3CC CC 1 45C72 C1C FAD RZe,1 01327 06C1 CC C 02666 C1C L32.2 STO RGP C133C 05CC CC 1 45C72 C1C CLA RZB,1 C1331 C6C1 CC C C066C C10 STI RTE P C1332 05CC CC 1 51776 C1C CLA ZZ8,1 01333 C3C2 CC C C2677 C1C FSB CELTA C1334 06C1 CC C C2661 ClC STC ZTE'P 0 1335 0C74 CC 4 C1726 C1C TSX PCTNTL,4 C1336 3 CCCCC C 2145C C1C PAR RZA C13'7 3 CCCCC C 26355 ClC PAR ZZA C1340 C6C1 CC C C2642 C1C STO Tl C01341 CC74 CC 4 C1726 C1C TSX PCTNTL,4 C1342 3 CCCCC C C7636 C1C PAR RSA C1343 3 CCCCC C 14543 C 1C PAR ZSA 01344 03C2 CC C C2642 C1C FSB TI C1345 06C1 CC C C643 C 1C STO T2 C1346 05CO CC 1 51776 C1C CLA ZZR,1 C1347 03CC CC C C2677 C1C FAD CELTA 0135C 06C1 CC C C2661 C1C STi ZTE~P C1351 CC74 CC 4 C1726 C1C TSX PCTNTL,4 C1352 3 CCCCC C 2145C C1C PAR RZA 01353 3 CCCCC C 26355 C1C PAR ZZA C1354 C6CL CC C C2642 C1C STO TI 01355 CC74 CC 4 C1726 C1C TSX PCTNTL,4 C1356 3 CC(CC C C7636 C1C PAR RSA 01357 3 CCCCC C 14543 CLC PAR ZSA 0136C C3C2 CC C C2642 C1C FSB TI C1361 0362 CC C C2643 C1C FSB T2 C1362 0241 CC C C2677 C1C FCP CELTA 01363 C260 CC 0 52177 1C FMP =.C3S7667 C1354 C3CC CC 1 51776 C1C FAD ZZB,1 C1365 C6C1 CC C C2667 C1C SIC ZGP C1-26 0560 CC 1 45C72 C1C LCG RZB,1 C1367 C26C CC 1 45C72 C1C FMP RZB,1 C137C C6C1 CC C C2642 C0C STO T1 C1371 056C CC 1 51776 CIC L0Q ZZe,1 C1372 026C CC 1 51776 ClC FMP ZZe,I C1373 03CC CC C C 0 C2642 C1C FAC T1

PAGE 14 UNIVERSITY OF MICFIGAN 01374 0601 CC C C2642 ClC STO T1 C1375 C131 CC C CCCCC CC XCA 01376 026C CC C C2642 C1C FMP T1 C1377 06C1 C CC C2643 ClC STO 12 014CC OC74 CC 4 CCCC5 C1C CALL SCRIT1 014C2 0241 CC C C2643 C1C FCP T2 014C3 C26C CC C C27CC C1 FMP K C14C4 0601 CC C C2642 C1C STO Tl C14C5 0131 CC C CCCCO CC XCA C14C6 C260 CC 1 450C72 C1C FMP RZB,1 014C7 C760 CC C CCCC2 CC CHS C141C 03CC CCC C2666 C1C FAD RGP 01411 06C1 CC 1 45C71 C 1C STO RZB-1,1 01412 -C120C CC C C15C7 C1C TMI L33A C1413 0560 CC C C2642 C1C LCQ T1 C1414 026C CC 1 51776 C1C FtP ZZB,1 C1415 0760 CC C CCOC2 CC CHS 01416 C3CC CC C C2667 C1C FAC ZOP C1417 0C3C2 CC C C2676 C1C FSB VCT C1420 0601 CC 1 51775 C1C STO ZZB-1,1 C1421 04C2 CC C C2701 C1C SUB ZMIN 01422 -C12C CC C C15C6 C1C TMI L33 C1423 0500CO CC 1 45C72 C1C CLA RZe,1 C1424 -01CC CC C C1430 C1C TNZ *+4 01425 05CC 0C 1 51775 C1C CLA ZZB-1,1 01426 03C2 CC 1 51776 C1C FSB ZZB,1 C1427 012C CC C C15C06 C1C TPL L33 C143C 0560 CC 1 450C71 C1C LCQ RZB-1,1 01431 C260 CC 1 45C71 C1C FMP RZB-1,1 01432 06C01 CC C C2711 C10 STO TEP C1433 0560 CC 1 51775 C1C LCQ ZZB-1,1 C1434 026C CC 1 51775 C1C FMP ZZB-I,1 01435 C3CC CC C C2711 C1C FAD TEMP 01436 04C2 CC C C2674 C1C SUB A2 C1437 -C12C CC C C15C6 C10 TMI L33 0144C 05CC CC 1 51776 C1C CLA ZZo,1 C1441 03C02 CC 1 51775 C1C FSB ZZB-1,1 C1442 0601 CC C C2711 C1C STO TEMP C1443 0500CC CC 1 45C72 C1C CLA RZB,1 C1444 03C2 CC 1 45C71 C1C FSB RZB-1,1 01445 0241 CC C C2711 C1C FCP TEMP 01446 -C06CC CC C C2726 C1C STQ B C1447 C260 CC C C2726 C1C FvP B C1450 06C1 CC C C2727 C 1C STU 82 C1451 C0560 CC 1 51776 C1C LCO ZZ8,1 01452 C260 CC C C2726 C1C FMP B C1453 03C2 CC 1 45C72 ClC FS8 RZe,1 01454 C601 CC C C2730 ClC STO B1 01455 0131 CC C CCOCC CC XCA C1456 C260 CC C C273C ClC FMP B1 01457 06C1 CC C C2731 ClC STO 812 01460 03C2 CC C C2674 C1C FSB A2 C1461 C241 CC C C2731 C1C FOP B12 01462 -C06CO CC C C2711 C1C SC TEMP C1463 C0500CC CC C 522C2 CC CLA =1. 01464 C241 CC C C2727 C lC FCP 82 C1465 0131 CC C CCCCC CC XCA 01466 C3C0 CC C 522C2 C1C FAD =1.

PAGE 15 UNIVERSITY OF MICFIGAN 01467 C131 CC C CCCCC CC XCA C147C 0260 CC C C2711 C1C FMP TEMP C1471 04C2 CC C 522C2 C1C SUB =1. 01472 O1CC CC C C150C6 C1C TZE L33 01473 C120 CC C C1477 C1C TPL L32e 01474 C560 CC 1 51776 C1C LCQ ZZE,1 C1475 C26C CC 1 51775 C1C FMP ZZR-l,1 01476 -C12C CC 0 C15C6 C1C TMI L33 01477 1 CCCC1 1 C15CC C1C L328 TXI *+1,1,1 015CC 0774 CC 4 CCCCG CO L32C AXT **,4 **=INTERMEDIATE IMAXZB C150C1 1 CCCC1 4 C15C2 C1C TXI *+1,4,1 015C02 C634 CC 4 015C0 C1C SXA L32C,4 015C3 0634 CC 4 C151C C1C SXA L33e,4 C150C4 -3 CCCCC 4 C1260C CIC- S15 TXL L329,4,'* **={AX-1 C15C5 OC2C CC C C1511 C1C TRA L33C 01506 06C00 CC 1 51775 C10 L33 STZ ZZB-1,1 0150C7 C600CC CC 1 45C71 C1C L33A STZ RZe-,l1 0151C 0774 CC 4 CCCCC CO L33B AXT **,4 "*=INTERMECIATE IMAXZB 01511 -3 CCCCO 4 C1513 C1C L33C TXL *+2,4,**.*=I[AXZB 01512 -0634 CO 4 C1511 ClC SXD L33C,4 01513 1 CCCCC 2 C1514 C1C S31 TXI *+1,2,* **=hPtAX 01514 -3 cCCCO 2 C1246 CIC L34 TXL L3CC,2,**: *=NPhAX*(NuLM-1) 01515 -3 CCCCC 2 01253 C1C L34A TXL L31,2,** *'-=N'AX* UNu 01516 -0534 CO 2 C1511 C10 LXC L33C,2 01517 -0534 CC 4 C1511 C1C CXI LXD L330C,4 0152C 3 CCCCC 4 C1522 C1C L34AA TXH *+2,4,** **=PREVICUS IP'AXZB C1521 -0534 CC 4 C152C C1C LX0 L34AA,4 01522 0634 CO 4 01551 C1C SXA L348,4 01523 0634 CO 4 C1652 C1C SXA L35,4 01524 0754 CC 4 CCCCO CC PXA,4 C1525 -C760 CC C CCOC3 CC SSM 01526 C4CC CC C C2720 C1C ACD NMAX 01527 0734 CC 4 CCCCC CO PAX,4 0153C -0634 CC 4 C1562 C1C SXD L34C,4 01531 1 CCCCO 4 C1532 C1C S20C TXI *+1,4,** **=3*NAX 01632 -0634 CC 4 C1723 C1C SXC L37,4 C1533 -0634 CC 2 C152C CIC SXD L34AA,2 01534 C774 CC 2 CC001 CC AXT 1,2 01535 0774 CC 1 CCOCC CC L24B AXT **,1 **-PAXIMfM IMAXSB 01536 05CC CC 2 33262 C1C L24C CLA RSB,2 ShIFT NEW 01537 06C1 CC 2 C7636 C1C STO RSA,2 S IGMA C1540 0600 CC 2 33262 C1C STZ RSB,2 TRAJECTCRIES C1541 05CC00 CC 2 4C166 ClC CLA ZSB,2 TC A 01542 06C1 CC 2 14543 C1C STO ZSA,2 MATRICES 01543 06CC00 CC 2 4C166 C1C STZ ZS8,2 C1544 1 CCCC1 2 C1545 C1C TXI *+1,2,1 01545 2 CCCC1 1 C1536 C1C TIX L24C,1,1 01546 1 CCCCO 2 C1547 C1C L24C TXI *+1,2, * **=hKA-(AX-(AXUM IMAXSB) 01547 -3 CCCCC 2 C1535 C 1C L24E TXL L248,2,** *=hIpAX*NUM 0155C 0774 CC 2 CCCC001 CC AXT 1,2 01551 0774 CC 1 CCOCO CC L34B AXT **,1 *:=PAXItJUM IMAXZB 01552 0500CO CC 2 45C72 C1C L34C CLA RZB,2 01553 06C1 CC 2 21450 C1C STO RZA,2 01554 06CC CC 2 45C72 C1C STZ RZB,2 01555 05CO CC 2 51776 C1C CLA ZZB,2 01556 06C01 CC 2 26355 C1C STC ZZA,2 01557 C6C00 CC 2 51776 C1C STZ ZZB,2 0156C 1 CCCC1 2 C1561 C1C TXI *+1,2,1

PAGE 16 UNIVERSITY OF MICFIGAN 01561 2 CCCCI 1 C1552 01C TIx L34C,1,1 01562 1 CCCCC 2 015t3 CIC L34C TXI *+1,2,*. * =*=AX-(FAXIdUM IMAXZ8) 01563 -3 CCCO0 2 01551 CIC L34E TXL L348,2,** **=hPAX*NUM 01564 05Cn CC C C2714 CIO CLA CCLIT 01565 04CC CC C 522C6 C IC AED =1 01566 06C1 CC C 02714 CH STO COCNT 01567 C5CO CC 0 52204 CdO CLA =4 0157C 0601 CC C 02723 CIC STC VAR 01511 0774 CC 2 CCCC1 CC AXT 1,2 C1572 06CC CC C 02724 CLC STZ P 01573 0774 CC 1 CCCCC CO 125 AXT **,1 *=t=AXIVUt I"AXSB 01574 05CC CC C C2724 CLO CLA P 01575 04CC CC C 52206 CLC ACD =1 01576 0601 CC C C2724 CL1C STO P 01577 0774 CO 4 CCCCC CC L25A AX T **,4 **=SIGPA TRAJECTORIES STILL TO BE PRINTED C16CC 2 CCCC4 4 C1602 CLC TIX *+2,4,4 01601 0634 CC 4 02723 CHC SXA VAR,4 01602 0634 CC 4 C1577 C1C SXA L2SA,4 01603 CC74 CO 4 CCCC3 C10 PRINT FVAR,...,FCLT3,CCI\TPTPVPHIDELTATAtZMIN 016 16 05CO CC 0 C2723 CLC CLA VAR 016117 01C CC C C1641 10C TZE L27-3 01620 -1 CCCCC 2 07fa36 10C L26 lOP RSA,2 01621 -1 CCCCO 2 14543 C1C lHP ZSA,2 01622 05C0 CC C C2723 C1C CIA VAR 01623 0402 CC C 52206 C01 SUB =1 016 24 01CCC C C 01641 C 10 TZE L27-3 01625 -1 CCCCC 2 CCCO CO S27 ICP * *2 u**RSA-NPAX 01626 -1 CCCCC 2 (0CC0 CO lOP **42 **zZSA-NPAX C1627 05CC CC C 02723 010 CLA VAR 01630 C4C2 CC C 52201 010 SUB =2 016 3 1 0100 CC 0 01641 C 01 TZE L27-3 01632 -1 CCCCC 2 CC0CC CC S27A 7 OP + t 2 $*=RLA-2*N'MAX 01633 -1 CCCCC 2 CCCCC CC lOP **,2 **=ZSA-2*NMAX 01634 05CO CC C 02723 C1C CLA VAR 01635 0402 CC 0 52205 C01 SUB =3 01636 01CC CC C 01641 10C TiE L27-3 C1637 -1 CCCCC 2 CCCCC CC S278 ICP **,2 **=RSA-3*NVAX 01640 -1 CCCCC 2 CCCCO CC IP **,2.*=ZSA-3*NMAX 01641 1 CCCC1 2 01642 010 TXI *41,2,1 01642 2 CCCC0I 1 162C C1C TIX 126,1,1 C1643 -1 CCCCC C CCC00 CC ENOIO 01644 1 C0000 2 1645 C01C L27 TX I, 2 **=4*hPAX-IMAXIMUM IPAXSB) 0164.5 -3 CCCCC 2 01573 01C 12 TXL 125,2t** **=htAX*KUM 01646 0500 CC C 5204 C010 CIA =4 01647 0601 00 C 02723 010 STO VAR 01650 0774 CC 2 CC001 CC AXT 1,2 01651 06CC CC C 02724 010 STZ P 01652 0774 CC 1 CCCCC CC 135 AXT **,1 U**=AXIVLt IVAXZB 016 5.3 05CC 00 C C2724 C 10 CIA P 01654 04CC CC 0 52206 CIC AED =1 01655 0601 CC C C2724 010 STO P 01656 0774 CC 4 CCCCO CO L35A AXT **,4 **=ZETA TRAJECTORIES STILL TO BE PRINTED 01657 2 CCCC4 4 C1661 010 TIX *+2,4,4 0166C C634 CC 4 02723 (10 SXA VAR,4 01661 C634 00 4 01656 010 SXA 35A,4 01662 CC74 CC 4 C(003 010 PRIN T FVAR,...FCLT4,CCLNTPtTPVPHIOELTATAHZMIN 01675 05CC CC C 02723 010 CIA VAR 01676 0100 CC C 1172C C01 TZE 137-3

PAGE 17 UNIVERSITY OF MICFIGAN C16 7 -1 CCCCC 2 2145C C1C L36 ICP RZA,2 C17CC -1 CCCCC 2 26355 C C ICP ZZA,2 C17C1 C05CC CC C C2723 CIC CLA VAR C17C2 04C2 CC C 52206 C C SUB =1 C17C3 C1CC CC C C 1720 C1C TZE L37-3 C17C4 -1 CCCCC 2 CCCCC CC S2E ICP **,2 017C5 -1 CCCCC 2 CCCCC CC LOP **,2 C17C6 0500CC CC C C2723 C1C CLA VAR 017C7 04C2 CC C 522C1 C1C SUB =2 0171C ClCC CC 0 C]72C CIC TZE L37-3 C01711 -1 CCCCC 2 CCCCO CC S28A IOP **,2 C1712 -1 CCCCC 2 CCCOC CC ICP **,2 C1713 05CC CC C C2723 C1C CLA VAR C1714 04C2 CC C 522C5 C1C SUB =3 C1715 ClCC CC C 0172C C1C TZE L37-3 C1716 -1 CCCCC 2 CCCCC CC S2EB ICP **,2 01717 -1 CCCCC 2 CC0CO CO ICP **,2 0172C 1 CCCC1 2 C1721 C10 TXI *+1,2,1 C1721 2 CCCC1 1 C1677 C1G TIX L36,1,1 C1722 -1 CCCCC C CCOCC CC ENDIO C1723 1 CCCCC 2 Ci724 C1C L37 TXI *+1,2,** **=4hMNAX-(MAXIMUM IMAXZB) C1724 -3 CCCCC 2 C1652 C1C L3. TXL L35,2,** *=NphAX*KUM C1725 C0020 CC C CC7CC00 C1C TRA L2C C1726 0634 CC 4 C2064 C1C PCThTL SXA LLle,4 SAVE C1727 0634 CC 1 C2065 C10 SXA LLlS1 INDEX C173C 0634 CC 2 C2C66 C1C SXA LL2C,2 REGISTERS C1731 0500CC CC 4 CCOC CC1 CLA 1,4 TRANSFER CD C1732 0621 CC C C1754 C1C STA LL6 ARGLYENTS 01733 0621 CC C C2045 C1C STA LL14 01734 05C00 CC 4 CCCC2 CO CLA 2,4 C1735 0621 CC C C1756 ClC STA LL7 C1736 0621 CC C C2047 C1C STA LL15 C1737 0774 CC 2 CCCC1 C0 AXT 1,2 01740 0634 CC 2 C1747 C1 SXA LL2,2 C1741 0634 CC 2 C2C51 C1C SXA LL16,2 C1742 C6CC CC C C C2652 C1C STZ TT7 01743 0754 CC 2 CCOCC CO LL1 PXA,2 01744 0734 CC 1 CCCCC CC PAX,1 C174,5 06CC CC C C2647 C1C STZ TT4 C1746 C6CC CC C C2715 C1C STZ CCLNT2 01747 0774 CC 2 CCOCC CO LL2 AXT **,2 a*=1ST VEC SUBSCRIPT Ch CURRENT TRAJECTORY 0175C 0500CC CC 2 52166 C0C LL4 CLA VEC,2 C1751 C6C1 CC C C2644 C1C STO IT1 C1752 1 CCCC1 2 C1753 C1C TXI *+1,2,1 01753 0634 CC 2 C1747 C1C SXA LL2,2 01754 05CC CC 1 CCCCC CC LL6 CLA **,1 C1755 06C1 CC C C2662 CLC SI0 TEPPR C1756 0500CC CC 1 CCCCC CC LL7 CLA **,1 01757 C6C1 CC C C2663 C1C STO TEMPZ 01760C 1 77777 1 C1761 C1C TXI *+1,1,-1 01761 0522 CC C C1754 C1C XEC LL6 01762 C6Ci CC C C2664 ClC STO TEMPRR C1763 0522 CC C C1756 C1C XEC LL7 01764 0601 CC C C2665 ClC STC TEWPZZ C1765 CC74 CC 4 C2C7C C1C TSX PCNT,4 C1766 03CC CC C C2647 C1C FAD TT4 C1767 C6C1 CC C C2647 C1C STO TT4 0177C C56C CC C C2664 ClC LCC TEMPRP

PAGE 18 UNIVERSITY OF MICFIGAN C1771 C26C CC C C2664 C1C FMP TEIvPRR C1772 (:C1 CC C C2656 C1C STO TT2C C1'73 C560 CC C C2665 C1C LCO TEfPZZ C1774 C26C CC 0 C2665 CIC FMP TEFPZZ C1775 C3CC CC C C2656 CIC FAD TT2C C1776 C0C1 CC C C2656 C1C STO TT2C C1777 05C0 CC C C, C674 C1C CLA a2 C2CCC C241 CC C C2656 C1C FCP TT2C 02CC 1 -C06CC CC C C2656 C1C SIC TT2C C2CC2 C260 CC C C2665 C1C FVP TEtPZZ C2CC3 C6C1 CC C C2665 C1C STO TEIVPZZ C2CC4 C56C CC C C2656 C1C. LCC TT2C 02CC5 C260 CC C C2664 C1C FVP TEiPRR C020CC6 06C1 CC C C2664 C1C STIC TEIPRR C2CC7 C56C CC C C2662 C1C S32 LCQ TEIPR 0201C C26C CC C C2662 ClC FMP TtVPR C2C11 06C1 CC C C2656 C1C STI TT2C C2C012 0560 CC C C2663 C1C LDQ TEVFZ C2C13 C26C CC C C2663 C1C FVP TE!VPZ 02C14 03CC CC C C2656 CIC FAD TT2C C2C15 C6Ci CC C C2656 C1C STO TT2C 02016 C5CC CC C C?674 C1C CLA A2 C2C17 C241 CC C C2656 C1C FCP TT2C 02C2C -C6CC CC C C2656 C1C STC TT2C 02C21 C26C CC C C2662 C1C FMP TEPPR C2C22 06C1 CC C C2662 C1C STG TEPPR 02C23 0560C CC C C2656 C1C LCQ TT2C (]o 02C24 C26C CC 0 C2663 C1C FMP TEMPZ 02C25 C6C01 CC C C2663 C1C STIC TEfPZ C2C26 CC74 CC 4 CCCC5 CIC CALL SCRT,1T2C 02C30 06C CC C C02656 C1C STO TT2C 02C31 CC74 CC 4 C2C7C C1C TSX PCNT,4 02C32 C76C CC C CCCC2 CO CFS C2C33 0131 CC C CCCC0 CC XCA C2034 C26C CC C C2656 C1C FMP TT2C 02C35 03CC CC C 02647 C1C FAD TT4 02C36 C601 CC C C2647 C1C STO TT4 C2C37 05CC CC C C2715 C1C S'3 CLA CCLNT2 C204C C4CC CC CC 522C6 CIC ACD =1 02041 06C1 CC C 02715 C1C SIC COLKT2 C2C42 C4C2 CC C C272C C1C SL8 N/AX C2C43 O1CC CC C C2051 ClC TZE LL16 C2C44 1 CCCC2 1 C2C45 C1C TXI *+1,1,2 020C5 050C CC 1 CCCCC CC LL14 CLA **,1 **=RPTRX 02C46 -0100CC CC C C1754 C1C TNZ LL6 C2C47 05CO CC 1 CCCCC CC LL15 CLA **,1 **=ZPTRX 02C50 -01C0 CC C C1754 C1C TNZ LL6 02CS1 0774 CC 2 CCCCC CC LL16 AXT **,2 **=IST ITRX SUBSCRIPT Ch CURRENT TRAJECTORY 02C052 1 CCCCC 2 C2C53 CIC S11 TXI *+1,2,** **=N:AX C2C53 0634 CC 2 C2C051 C1C SXA LL16,2 02C054 0560 CC C C2644 C1C LCQ 1T1 02C055 026C CC C C2647 C1C FIP TT4 C2C56 C3CC CC C C2652 C1C FAD TT7 02C57 06C1 CC C C2652 C1C SfO TT7 02CEC -3 CCCCC 2 C1743 C1C LLI7 TXL LLI,2,** **=iAX*NLM 02C61 U5CC CC C C2652 C1C CLA TT7 C2C062 0241 CC C 52176 C1C FCP =2. C2C63 C26C CC C C2676 C1C FVP VCT

PAGE 19 UNIVERSITY OF MICHIGAN 02C64 0774 CC 4 CCCCC CC LL18 AXT **,4 **=ACCRESS CCPPLEMENT FOR TSX TO POTNTL C2C65 0774 CC 1 CC0OC CC LL19 AXT **,1 **=CCNTENTS GF XR1 BEFORE TSX TO POTNTL 02066 C774 CC 2 CCOOC CC LL20 AXT **,2 *=CCGhTENTS OF XR2 BEFOCRE TSX TO POTNTL 02C67 CC20 CO 4 CCCC3 CO TRA 3,4 02C7C 0634 CC 4 C2260 C1C PCNT SXA RET,4 02071 05CO CC C C2660 C1C CLA RTEMP 02072 0302 CC 0 C2662 C1C FSB TEMPR 02C73 06C.1 CC 0 C2646 C 1C STO TT3 C2C74 0131 CC C CCCOC CC XCA 02C75 0260 CC G C2646 CIC FMP TT3 02C076 0601 C0 C C2646 ClC STO TT3 02C77 05CO CC C C2661 C1C CLA ZTEIVP 02.1CC 03C02 CC C C2663 C1C FSB TEPPZ 021C1 06C1 CC C C2645 C1C STO TT2 021C2 0131 CC C CCOCO CO XCA C210C3 C260 CC 0 C2645 ClC FMP TT2 C2.104 03CC CC C C264f C1C FAD TT3 02105 0601 CC C C2646 C 1C STO TT3 02106 04C02 CC 0 C2675 C1C SUB V2 0210C7 C120 CC C C2131 C1C TPL LL2CA C211C 05CO CC C C2664 C1C CLA TEMPRR 02-111 03C2 CC C C2662 CIC FSB TEMPR 02112 0601 CC C C2645 C10 STO TT2 C2113 0'131 CC C CCCCC CO XCA C2114 0260 CC C C2645 C1C FMP 1T2 02115 0601 CC C C2645 CLC STO TT2 02.1 16 05C00 CC C 02665 C1C CLA TEMPZZ CO 02117 03C2 CO C C2663 C 1C FSB TEMPZ C02.120 0601 CC C C27C2 C1C STO ALP 02121 0131 CC C CC000C CC XCA 02122 0260 CC C C27,C2 CC10 FP ALP C212.3 03COC C C C2645 C1C FAD TT2 02124 0241 CC C 52175 C1C FOP =3.1415927 02125 0131 CC C CCOOC CO XCA 02126 0601 CC C C2645 C1C STO TT2 02127 04C2 CC C C2646 C10 SUB TT3 C2130 0120 CCC C2141 C1C TPL LL12 02131 05CC CC C C2660 C1C LL2CA CLA RTEIP NCRPAL C2132 C241 CC 0 52174 C1C FCP =.5 PCTENTIAL 02133 0260 CC C 2662 CIC FMP TEPPR CALCULATICN C2134 0601 CC C C2713 ClC STO CAPB 02135 03CC00 CC 0 C2646 C10 FAD TT3 02136 0601 CC C C2712 C10 STO CAPA 02.137 0074 CC 4 C2262 C1C TSX INPCNT,4 C2,140 CC20 CC C C2260 C1C TRA RET 02,141 OC74 CO 4 CC005 C10 LL12 CALL SQRT,TT2 CALCLLATICN 02143 06C1 CC 0 C265C C1C STO TT5 FCR PCINTS C2144 C600 CC C C2655 C10 STZ TT1C C2.145 05CC CC C C2662 C1C CLA TEMPR 02146 0241 CC C C2650C C1C FCP TT5 02.147 0260 CC C 52175 C10 FMP =3.1415927 SPHERE CF 02150 -C3C0 CC C 52173 C1C UFA =K233K9 EXCLLSICN C0151 0601 CC C C2653 C1C STO TT8 02152 C322 CC 0 52173 C C ERA =K233KS 02153 04C2 CC C 522C06 C10 SUB =1 02.154 C34C CC C 522C6 C1C CAS =1 C215-5 OC20 CC C C2162 C1C TRA *+5 C2156 CC20 CC 0 C2157 C 1C TRA *+1

PAGE 2C UNIVERSITY OF MICHIGAN 02157 05CO CC 0 52202 C LC CLA =1. 32160 06C1 CC 0 02653 CIO STO TT8 02161 OC20 CO C C2223 C1C TRA SS2+4 02162 06C1 CC C C2654 C1C STO TTS 02163 05CO CC 0 C2653 C IC CLA TT8 02.164 C3CO CC 0 52173 C10 FAD =K233KS 0216,5 0601 CC C C2653 C1C STO TT8 02166 05CO CC 0 52172 ClC CLA =6.2831853 0216,7 0241 CC C C2653 C C FOP TT8 C2170 -06CO CO C C27C6 C10 STQ DPHI 02171 0131 CC C CCOOC CO XCA 02.172 0601 CC C C2705 C1C SS3 STO PHIN 02173 OC74 CC 4 CCOC6 CIC CALL CGS,PFIh 02175 0760 CC C CCCC2 CO CHS 02176 03CC CC C 52202 C10 CAD =1. 02177 0601' CC C C2657 CLC STO JT4C C22CC 0074 CC 4 CCO05 C10 CALL SQRT,1T40 022C2 0601 CC C C2657 C1C STO TT4C 022C3 05CO CC C 522C2 C1C CLA =1. 022C4 0241 CC C C2657 CIC FDP TT40 02205 0'131 CO 0 CO0CO CO XCA 022C6 C3CC CC 0 C2655 C10 FAD TTIC 02207 060.1 CC 0 C2655 ClO STO TTC10 02210 05C CC Ca C2654 C10 CLA TT9 02211 04C2 CC 0 52206 C10 SUB =1 02212 01CO. CC C C2217 C10 TZE SS2 02213 0601 CC C C2654 C10 STO TTS Co C2214 05CO CC C C27CS05 C1C CLA PHIN 0221:5 C3CC CO 0 C2706 C10 FAD OPHI 02216 0-C20 CC C 02172 C10 TRA SS3 02117 05CO CC C C2655 C10 SS2 CLA TT1C C222C 0241 CC C C2662 C1C FOP TEHPR 02221 0260 CC 0 52171 C10 FMP =.7C711 02222 0601 CC C C2655 C10 STO TT1C 02223 05CO CC C 52170 C10 CLA =1.5 02224 0241 CC C C2650 ClC FOP TT5 0222-5 0131 CO 0 ccoC0 CO XCA 02226 03CO CG C C2655 CCO FAD TT1C 02227 0131 CC C CCCCC CC XCA 02230 0241 CC C C2653 ClC FOP TT8 02231 0260 CC 0 52167 C10 FMP =6.2831853072 02232 0601 CC 0 C2655 CO1 STO TT1C C2233 C074 CC 4 CCCC5 CLC CALL SQRTtT3 02235 06C..1 CC C C2646 ClC STO TT3 02236 C5CO CC 0 C265C C 1C CLA TT5 C2237 024.1 CC C C2646 C 1C FCP TT3 0224C -06CC CC C C2646 C C STQ TT3 02241 05CC CC C C2660 ClC CLA RTEIP 02242 03C2 CC C C2662 C10 FSB TEMPR 02243 0131 CC C CCCCO CO XCA 02244 0260 CC C C2646 C1C FMP TT3 02245 0300 CC C C2662 C 1C FAD TEMPR 02246 0241 CC 0 52174 C10 FCP =.5 02247 C260 CC C C2662 C1C FMP TEPPR 02250 0601 CC C C2713 CIC STO CAPB 02251 03CC CC C C2645 CIC FAD TT2 02252 0601 CC 0 C2712 C1C STO CAPA 02253 OC74 CC 4 C2262 C 1C TSX IKPCNT,4

PAGE 21 UNIVERSITY OF MIC-IGAN 02254 03C2 CC C C2655 CIC FSB TT1C 02255 0241 CC C C2646 C1C FOP TT3 02256 0131 CC C CCCCC CC XCA 02257 03CC CC 0 C2655 C10 FAD TT1C 0226C 0774 CC 4 CCCCC CO RET AXT **,4 C2261 CC20 CO 4 CCCC1 CO TRA 1,4 C2262 0634 CC 4 C2347 CIC INPCNT SXA LL29,4 02263 CSCC CC C C2713 C10 CLA CAPE 02264 024.1 CC C C2712 C 10 FOP CAPA 02265 0131 CC C CCCCC CO XCA 02266 06C2 CC C 027C2 C1C SLW ALP NEST MACRO X LCQ X+4 FMP M FAD X+3 XCA FMP M FAD X+2 XCA FMP M FAD X+1 XCA FMP M FAD X NEST END 02267 05CC CC C 52202 C 1C CLA =1. C2270 03CC CC C C27C2 C1C FAD ALP 0CD C2271 C6C1 CO 0 C271C 1 C STO TEMP 02272 05CO CC C 522C2 C1C CLA =1. C2273 03C2 CC C C27C2 C C FSB ALP C2274 0241 CC C C2711 C1C FDP TEMP 022a75 -C06CC CC C C27C7 C IC STQ M 02276 CC74 CC 4 CCC07 C1C CALL ELCGi, C23C0 0601 CC C C2711 C1C STO TEMP C23C1 NEST BB 02315 0131 CC 0 CCOCC CO XCA 02316 026C CC C C2711 C1C FMP TEMP C2317 OC1 CC C C2711 C1C STO TEuMP 02320 NEST AA 02334 03C2 CC C C2711 C1C FSB TEMP 0233,5 0601 CC C C271C ClC STO KM C2336 05CC CC C 2712 ClC CLA CAPA 02.337 03CC CC C C2713 C1C FAD CAPB C2340 06C1 CC C 02711 C1C STO TEMP 02341 0074 CC 4 CCCC5 C1C CALL SCRT,TEMP C2343 06C1 CC C C2711 C1C STO TEMP C2344 05CC CC C C2710 C1C CLA KM 02345 0241 CC C C2711 C1C FOP TEMP C2346 0260 CC C 52166 C1C FMP =4. C2347 0774 CC 4 CCCO0 CC LL2S AXT **,4 0235C C02C CC 4 CCC01 CC TRA 1,4 02351 0074 CC 4 CCC03 C1C ENC PRINT FMT,C C2354 0074 CC 4 CCC1C C1C CALL ERRCR C2355 6521516C6C6C CC BCI 1,VAR FCRFAT VARIABLE C2356 0 CC(CC 1 C2723 C1C PZE VAR,1 SYMeCL 02357 0 CCCCC C CCCC2 CC FVAR PZE 2 TAELE 02360 02C52CC1633C CO FMT BCI *,25H1TFE BEGINNING CF THE ENC* C2365 C626ClCC73C4 CO FIN BCI *,6FlC,4I5/(2FlC)

PAOE 22 UNIVERSITY OF MICHIGAN C237C 3C2CC1623127 CO FCLIl BCI *,H+1SIGtA FLCh LINES - INITIAL APPRCXIMATICN+S59,5HPAGE C24C1 6031C273C43C CC ETC I2,4H CF I2/1P-/S1C,3HV =F7.2,Sl0,5HPHI =F7.2,Sl0,9HDELTA 02412 254363216C63 CC ETC T =F7.4,S1C,3HA =F7.4,S1O,3HH =F7.4,SI0,6HZMIN =F7.2/lH02424 4C62C4731465 CO ETC S4,'VAR'(S9,11-RS12,1HZS8)//'VAR (S8,WF10.4,S3tWF10.4)* 02436 3C2CC16C7125 CC FCLT2 BCI *,H+1 ZETA FLCk LINES - INITIAL APPROXI'ATION+S59,5HPAGE C2447 6031C273C430 CC ETC I2,4H CF 12/IH-/S1C,3HV =F7.2,S10,5HPHI =F7.2,S1l0,9HDELTA 0246C 254363216C63 CO ETC T =F7.4,S10,3HA =F7.4,S0O,3HH =F7.4,S10,6HZMIN =F7.2/IH02472 4C62C4731465 CC ETC S4,'VAR'IS9,11RS12,1HZS8)//eVAR'(S8,hFIO.4,S3,WFIO.4). 025C4 302CC1623127 CC FCLT3 BCI *,H+1SIGPA FLCk LINES - ITERATION +12,S83,KHPAGE 12,4H OF 0251'5 266031C261C1 CO ETC 12/1I-/SlC,3HV =F7.2,SlO,5hPHI =F7.2,S10,9HDELTA T =F7.4 02-526 260733047362 CO ETC,SIC,3HA =F7.4,SIC,3H =F7.4,SIO1 6HZMIN =F7.2/1H-S4,'VAR'( 0254C 74621C736626 CC ETC S8,hFlC.4,S3,k Fl0.4)* 02:544 302CC16C7125 CO FCLT4 BC! *,H+1 ZETA FLCU LINES - ITERATION +I2,S83,KHPAGE I2,4H OF 02555 266C31C26101 CO ETC I2/1H-/SlC,3HV =F7.2,S1O,5HPHI =F7.2,S10,9HDELTA T =F7.4 02566 26C733047162 CO ETC,S1C,3HA =F7.4,SlC,3H =F7.4,S10,6HZMIN =F7.2/lH-S4,'VAR'( C26CC 74621C736626 Cc ETC Se,hF1C.4,S3,hFlC.4)* 026C4 02053C4C2225 CO F BC! *r,25H-BEGINNING WCRK CN R = F10.4,S3,4HZ = Fl0.4,S5,50 C2615 30633C256C63 CO BCI *,I-THE TIME IS (IN 1/60THS CF A SECOND PAST MIDNIGHT)I9* 02626 +2C154271C277 CO AA DEC 1.38629436112,.C9666344259,.03590092383,.03742563713 02632 +1727334155C05 CO CEC.C145119212 02633 +2004CCCCCCCC CO EB DEC.5,.124585S3557,.06880248576,.03328355346,.C0441787012 02640 ASSIGN X2,PSk 02642 ASSIGN TL,T2,TT1,TT2,TT3,TT4,TT5,TT6,TT57TT68TTTTTTIO,TT20,TT40 0266C ASSIGN RTEMP, ZTEMP, TEPRTEMPZTEPPRR, TEMPZZ RGP ZGP C2670 ASSIGN V, PHI,OELTAT,A A2,Y2,VCTtELTA,K,ZMIN, ALP,H, CAPPIN DPHI 02707 ASSIGN M,K,tTEPP,CAPA,CAP8 02714 ASSIGN COLNTCCUNT2,NA,tNL NAX tNITER, ISWVARPtTP D c02726 ASSIGN 8,82,81,812 C3N C47C4 SIZE SYN 25CC C7636 BES SIZE 07636 0 C0CCCC C CCCCC CO RSA PZE 1454 3 BES S I ZE 14543 0 CCCCC C CCCCO CC ZSA PZE 2145C BES SIZE 2145C 0 COCCC C CCOCC CC RZA PZE 20355 BES SIZE 26355 0 CCCCC C CCCCC CC ZZA PZE 33262 RSB BES SIZE 40166 ZSe BES SIZE 45C72 RZ8 BES SIZE 51776 ZZ8 BES SIZE 52166 VEC BES 120 BRIEF OFF END P-RGRAW LITERALS 52166 2C34CCCCCCCO CC 52167 203e22C77325 CC 5217C 2G1ECCCCCCCC CC 52171 2 CC 52024 1 CC 52,172 2C3E22077324 CC 52173 233CCCCCCCCC CC 52174 2CC4CCCCCCCC CC 52175 2C2622077326 CC 52176 2C24CCCCCCCC CC 52177 1745C5745716 CC 522CC OClCCCCCCCCC CC 52201 CCCCCCCCCCC2 CC

PAGE 23 UNIVERSITY OF MICHIGAN 52202 2C14CCCCCCCC CC 52203 6016CCCCCCCC CC 52204 CCCCCCCCCCC4 CC 52205 CCCCCCcCCCc3 Co 52206 CCCCCCCCCCC1 CO *.i* 52207 Is THE FIRST LCCATICN IOT LSED BY THIS PRCGRAI. OD

%5 o3f