FPPAC88: A two-dimensional multispecies nonlinear Fokker-Planck package

Computer Physics Communications 51(1988) 373—380 North-Holland, Amsterdam

373

FPPAC88: A TWO-DIMENSIONAL MULTISPECIES NONLINEAR FOKKER—PLANCK PACKAGE * A.A. MIRIN, M.G. McCOY, G.P. TOMASCHKE and J. KILLEEN National Magnetic Fusion Energy Computer Center, Lawrence Livermore National Laboratory, Livermore, CA 94550, USA Received 27 May 1988

FPPAC88 solves the complete nonlinear multispecies Fokker—Planck collision operator for a plasma in two-dimensional velocity space, The operator is expressed in terms of spherical coordinates (speed and pitch angle) under the assumption of azimuthal symmetry. Provision is made for additional physics contributions (e.g. rf heating, electric field acceleration). The charged species, referred to as general species, are assumed to be in the presence of an arbitrary number of fixed Maxwellian species. The electrons may be treated either as one of these Maxwellian species or as a general species. Coulomb interactions among all charged species are considered. This program is a NEW VERSION of FPPAC, which was published in Computer Physics Communications in 1981. The most important new feature is the modification of the spatial differencing to accommodate situations in which advection dominates diffusion. Another important change is the discontinuance of the CDC-7600 version of FPPAC; thus the NEW VERSION, FPPAC88, is applicable to machines with a single fast memory, in particular the Cray-i with the CFT compiler. Any departures from standard FORTRAN are easily accommodated, so that FPPAC88 may be run on almost any supercomputer with a FORTRAN-77 compiler.

PROGRAM SUMMARY Title of new version: FPPAC88

Operating system or monitor under which the new version is executed: Cray Time-Sharing System (CTSS)

Catalogue number: ABFI Programming language used in the new version: FORTRAN-77 Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue) Reference to original program: Comput. Phys. Commun. 24 (1981) 37 Authors of original program: M.G. McCoy, A.A. Mirin and J. Killeen Computer for which the new version is designed and others on which iris operable: Designed for the Cray-i at the National MFE Computer Center, Lawrence Livermore National Laboratory, Livermore, California; applicable to any singlememory supercomputer with a FORTRAN-77 compiler. *

Work performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract W-7405-ENG-48.

High speed storage required: approximately 250000 words for Test Problem 1, which has a 71(v) by 31(8) mesh, 2 general species, 1 Maxwellian species and 9 Legendre polynomials No. of bits in a word: 64 No. of lines in combined program and testdeck: 7268

Keywords: two-dimensional, multispecies, implicit, nonlinear, Fokker—Planck Nature of the physical problem The complete nonlinear multispecies Fokker—Planck coffision operator for a plasma in two-dimensional velocity space is solved. The operator is expressed in terms of spherical coordinates (v = speed, 0= angle between velocity and magnetic field directions, 4~= azimuthal angle) under the assumption of azimuthal symmetry. Provision is made for additional physics contributions.

OO1O-4655/88/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

374

A.A. Mirin et al.

Method of solution

The Fokker—Planck ecluation is solved using finite differences. Spatial derivatives are approximated by central differences, with the exception of the advective terms which are approximated by combined central/upwind formulae designed to accommodate situations in which advection dominates diffusion. Time-advancement is accomplished through either implicit operator splitting, an alternating direction implicit (ADI) algorithm, or fully implicit differencing. (In the latter case the user must supply his own nine-banded linear systems solver.) The Fokker—Planck coefficients and their derivatives are computed by expanding the distribution functions and the Rosenbluth potentials in Legendre series, and equating the respective series coefficients. Reasons for the new version

The most important new feature is the modification of the spatial differencing to accommodate situations in which advection dominates diffusion (e.g. the slowing down of 14 MeV fusion protons in a deuteron background). Another important

/ FPPAC88 change is the discontinuance of the CDC-7600 version of FPPAC; thus the NEW VERSION, FPPAC88, is applicable to machines with a single fast memory, in particular the Cray-I with the CFT compiler. Restrictions on the complexity of the problem

The user must specify the number of meshpoints in the two coordinate directions as well as the number of Legendre polynomials used to calculate the Fokker-.Planck coefficients. Sufficient accuracy for most problems is attainable on the Cray-I. However, double precision should be used on a 32-bit-word machine. Ts’pical running time

9.6 ~ss/meshpoint species on the Cray-I are required to compute the Fokker—Planck coefficients. 1.9 p. s/meshpoint are required to time-advance the distribution function for one species using implicit operator splitting. These times vary to some extent with the mesh configuration due to variations in vectorization efficiency.

LONG WRITE-UP 1. Introduction In the simulation of magnetically confined plasmas, there are a number of situations where the charged particle distribution functions are not Maxwellian. For cases where a precise knowledge of those distribution functions is important, kinetic equations must be solved. Examples include the heating of plasmas by energetic neutral beams or microwaves, the thermalization of a particles in D—T plasmas, runaway electrons in tokamaks, two-energy component fusion reactors, and mirror machines in which there is an expelling ambipolar potential. The kinetic equation of interest is the Boltzmann equation with Fokker—Planck collision terms [1]. The problem is to solve a nonlinear partial differential equation for the distribution function of each charged species in the plasma in terms of seven independent variables (three spatial coordinates, three velocity coordinates, and time). Since such an equation, even for a single species, exceeds the capability of any present computer, simplifications are required to make the problem tractable. The number of spatial dimensions can be reduced in the event that either the magnetic field is uniform (or square well), or that the magnetic

field inhomogeneity enters the problem principally through its variation along the direction of the field (i.e., the perpendicular drift velocity is neglected relative to the thermal velocity). Further reduction is realized by gyrophase averaging. This is warranted by the large difference between the time scale of the gyromotion and that of the collisional relaxation of the population of gyrocenter orbits. As a consequence, for many applications the number of independent variables may be reduced to four: one spatial variable, z; two velocity space coordinates, v and 0, the speed and pitch angle; and time, t. The increased tractability afforded by these basic assumptions, or more restrictive ones, has fostered an evolution of numerical Fokker—Planck calculations over the past thirty or so years [2]. One of the most useful numerical Fokker— Planck models applicable to a uniform or squarewell magnetic field over the past decade has been FPPAC [3], which time-advances the Fokker— Planck equations in two-dimensional velocity space (v, 9) for an arbitrary number of charged “general” species in the presence of an arbitrary number of fixed background species. Since FPPAC was first published, it has undergone a nurnber of revisions and improvements. The most im-

A.A. Mirin et al.

portant change has been the modification of the spatial differencing to accommodate situations in which advection dominates diffusion. This is accomplished by using a scheme similar to that of Chang and Cooper [4]. Another important change has been the discontinuance of the CDC-7600 version of FPPAC. Thus, the NEW VERSION, FPPAC88, is applicable to machines with a single fast memory, in particular the Cray-i with the CFT compiler. Note however that the departures from standard FORTRAN are easily accommodated, so that FPPAC may be run on almost any machine with a FORTRAN-77 compiler. FPPAC88 consists of a driver and a package. The driver/package may be run as a stand-alone code, thereby serving as the kernel of a new code, or the package itself may be stripped from the driver and incorporated into an existing code. Also, additional physics may be represented, within limits, by merely setting existing array elements to appropriate values. Because FPPAC88 is a NEW VERSION, this write-up will discuss only new issues. The reader is referred to the original publication [3] for further information.

/

FPPAC88

375

with 6,,J+1/2

=

~exp[

—

(RjJ+l/2/Ro)21.

R1~1/2 ~~j+ 1/2 I A11~1/2 I / I B1~1/2 I~ (4) The notation is the same as in section 8 of ref. [3], and R0 is a parameter. Using R0 equal to ~ recovers central differencing, whereas R0 equal to a very small number results in upwind differencing. We have found that a value of R0 equal to 3.5 yields quite satisfactory results. Discretization of the Df term is treated analogously: =

8(Df) I

[(Df)i+l/ 2,J

a

—

(Df)i_I/2.j]/~0,, (5)

.

where (Df).+l/2.

=

61÷l/21(Df)i,j

+ (i

—

6I+I/21)(Df)i+l.J,

(6)

with 6i+1/2,j

=

~exp[

—

(RE+l/ 21/Ro)2],

2. Modified spatial differencing

~

[(Af)1,j+l/2

—

(7)

and

The spatial differencing has been modified to more accurately treat situations in which advection dominates diffusion. Examples of such scenarios are the injection of a very high energy neutral beam or the slowing down of high energy fusion reaction products (e.g. 14 MeV protons). In this technique the advective terms are discretized by taking a weighted average of the previously used central differencing and of one-sided (upwind) differencing [5], according to the cell Reynolds number [6]. Eq. (56) of ref. [3] is replaced by ~(Af)

(3)

Here, R~1+ 1/2 is the cell Reynolds number:

(Af)i,J_l/2I/~, (1)

~°i+1/2 I DI+l/211/I ~+1/2,j I. (8) This technique, although developed independently, is very similar to that of ref. [4]. The main difference lies in the expression for 6(R) (see eqs. (3) and (7)). In both methods 6 is monotonically decreasing from 0.5 to 0 as R goes from 0 to oc. For the method here 6 is much smaller at very large R. This approach virtually eliminates the problem of the distribution function going substantially negative in situations where slowing down greatly exceeds scattering. Previously, the only alternative would have been to use an extremely large number of meshpoints, making the problem intractable.

R~~1/2J=

where

(Af )i,j+ 1/2 =6,J+1/2(Af)I,J+ (i —6,J±l/2)(Af)l,J÷I, (2)

3. Library subroutines FPPAC88 makes use of several optimized library and inline subprogram which reside on the

A,A. Mirin et al.

376

National MFE Computer Center system. For the benefit of users on other systems these routines (BCAST, SCOPY, SSCAL and CVMGT) are located toward the end of the source deck and are currently “commented out”. Users on other systems must remove the “ccc” from colunms 1—3.

/ FPPAC88 20, JXA 101, MXA 8 and KSYDMA 0. For non-NMFECC users, in order for the source to properly compile, the “.IF” and “.ENDIF” statements on lines 3309 and 3542 must be removed, and lines 3310—3341 should also be removed in the event PARAMETER NFCGA is not greater than 0; NMFECC users should precompile with PRECOMP. =

=

=

=

4. New COMMON variables FPPAC88 has additional COMMON variables, as follows: Variable R.Z in block NPARAM3 contains R0 (see eq. (3)); variable NFLG in block NXIN1 is described in user responsibility number 4 in the code listing; and arrays DLX and DLY in block NVALS2 stand for 6~~÷ 1/2 and 6,~ 1/2,1’ respectively (see eqs. (3) and (7)).

6. Test problems The same two test problems as in the original writeup are provided. The results are very slightly different because of the change in the difference algorithm. Because of this difference (although the densities and energies agree to at least three significant digits), updated output is provided.

5. Compile / load/ go deck References The compile/load/go deck makes use of the COSMOS job control language at the National MFE Computer Center. The actual FPPAC88 source comprises lines 3—7203, test problem number 1 input comprises lines 7208—7223, and test problem number 2 input comprises lines 7255—7266. The PARAMETER statements are set up for test problem number 1. In order to run test problem number 2, the source file must be changed to reflect the correct PARAMETER values, namely: MEQA 1, MEXA 0, NBOA 1, IYA =

=

=

[1] M.N. Rosenbluth, W.M. MacDonald and D.L. Judd, Phys. Rev. 107 (1957) 1. [21 J. Killeen, G.D. Kerbel, MG. McCoy and A.A. Mirin, Computational Methods for Kinetic Models of Magnetically Confined Plasma (Springer-Verlag, New York, 1986). [3] M.G. McCoy, A.A. Mirin and J. Killeen, Comput. Phys. Commun. 24 (1981) 37. [4] J.S. Chang and G. Cooper, J. Comput. Phys. 6 (1970) 1. [5] D. Potter, Comput. Phys. (John Wiley and Sons, London, 1973). [6] P.K. Roache, Computational Fluid Dynamics (Hermosa, Albuquerque, 1972).

A.A. Mirin et aL

/ FPPAC88

377

TEST RUN OUTPUT

Test 1

output to fokker-planck package with driver, for questions contact either m. g. mccoy or a. a. mirin of the national magnetic fusion energy computer center at the lawrence livermore national laboratory. see subroutine

initial for comments regarding input values.

number of velocity mesh points (jx> = 71 number of theta mesh points (iy) = 31 number of general,species inbo) = 2 number of maxwellian ions (meq) = 0 highest order legendre polynomial (mx) = 4 ksydm = 1, distributions functions not assumed symmetric about theta pi/2. mex = 1, electrons are maxweliian. itest= 1 nstop= 1000 nprint= 250 ncoef= 1 xmax= 0.l00000e+01 dtr= 0.800000e-02 secs vnorm= 0.600000e+09 cm/sec

general species follow species number mass= bnumb=

1

0.334330e-23 grams 0.l00000e+01

initial values for this species reden= 0.340000e+14 particles/cm**3 enit= 0.400000e102 key xnit= 0.1000008103 ynit= 0.l00000e+02 vnit 0.l00000e+O1 description of loss term for this species tau0 0.250000e-01 secs taul 0.200000e+03 secs description of source for this species asor= 0.650000e+14 particles/cm**3/sec esor 0.990000e+02 key xsor 0.940000e+03 ysor= 0.275000e+03 usor 0.l00000e+01 **** ** * * ***

** ** ** ******

species number

******

*

2

11 time step 1000 time= species no. 1 density species no. 2 density= species no. 3 density=

0.800000e+01 secs 0.729073e+14 particles/cm**3 0.527625e+14 particles/cm**3 0.700000e+14 particles/cm**3

energy= energy= energy=

0.364319.+02 key O.441727e+02 key O.800000e+01 key

378

A.A. Mirin et al.

/

FPPAC88

distribution function of species number I at normalized velocity xi 1)= 0. as a function of theta in radians from theta = 0. to theta = 0.314159e+0l 16359e+02 0.216359e-f02 0.216359e+02 O.216359e+02 0.216359e+02 0.216359e+02 0.216359e+02 0.2163598102 0.216359e402 0.216359e+02 0.216359e402 0.216359e402 0.2 0.216359e+02 0.216359e+02 0.216359e+02 0.216359e+02 0.216359e+02 0.216359e102 0.216359e+02 0.216359e+02 0.216359e+02 0.216359efOZ 0.216359e+02 0.216359e402 0.216359e+02 0.216359e102 0.216359e+02 0.216359ef02 0.216359e+02 0.216359e+02 0.21 6359e402 distribution function of species number 1 at normalized velocity xi 21= 0.142857e-01 as a function of theta in radians from theta = 0. to theta = 0.314159e+01 0.181144e+02 0.181193e+02 0.181273e-f02 0.181398e+02 0.181567e+02 0.181780ef02 0.182034e+02 0.182329e402 0.182660e+02 0.183025e+02 0.183419e+02 0.183839e-f02 0.184279e+02 0.184737e+02 0.185205e102 0.185681e102 0.186158e+02 0.186631e102 0.187095e+02 0.187545e+02 0.187977ei02 O.188385e102 0.188764e+02 0.189111e+02 0.189421e402 0.189691e402 0.189918e+02 0.190099e+02 O.190232e+02 0.190318e402 0.190371 e+02 distribution function of species number 1 at normalized velocity xi 31= 0.285714e-01 as a function of theta in radians from theta = 0. to theta 0.314159e+01 0.178707e+02 0.178801ef02 0.178951e+02 0.179185e+02 0.1795028102 0.179904e+02 0.180388e+02 O,180949e+02 0.181583e+02 0.182284e+02 0.183047e+02 0.183863e+02 0.184725e+02 0.185625e102 0.186552e+02 0.157496e+02 0.188449e+02 0.189399e+02 0.190337e+02 0.191250e-f02 0.192129e+02 0.192964e+02 O.193743e*02 O.194459e+02 0.195101e+02 0.195661e+02 0.196133e102 O.196509e+02 0.196788e102 0.196969e+02 0.1 97079e+02 distribution function of species number 1 at normalized velocity xi 41= O.428571e-01 as a function of theta in radians from theta = 0. to theta = 0.31 41 59e+01 0.178508e+02 0.178642e+02 0.178857e+02 0.179189e+02 0.179644e+02 0.180221e102 O.180918e+02 0.181731e+02 O.182654e+02 0.183681e+02 0.184803e+02 0.1860108102 0.1872918102 O.188634e+02 0~190026e+02 0.191451e+02 0.192895e+02 0.194342e+02 0.195776e+02 O.197179e+02 0.198536e+02 0.199828e+02 0.201040e+02 0.2021568102 0.203160e402 0.204038e+02 0.204780e+02 0.205373e+02 0.205813e+02 O.206098e+02 0. 206273e+02 distribution function of species number 1 at normalized velocity xi 51= 0.571429e-01 as a function of theta in radians from theta = 0. to theta 0.314159e+01 0.177422e+02 0.177594e+02 0,177865e+02 0.178286e+02 0.178862e+02 0.179597e+02 0.180489e+02 0.181533e+02 0.182726ef02 0.184059e-f02 0.185522e+02 0.187103e+02 0.188790e+02 0.190566e402 0.192415e+02 0.194318e-f02 0.196254e102 0.198202e+02 0.200140e+82 0.202044e+02 0.203892e+02 0.205659e-f02 O.207321e+02 0.208856e+02 0.210242e+02 0.211458e+02 0.212485e402 0.213310e102 0.213921e+02 0.214318e+02 0.21 4562e+02 distribution function of species number 1 at normalized velocity xi 61= 0.714286e-01 as a function of theta in radians from theta = 0. to theta = 8.31 41 59e+01 0.173867e+02 0.174071e+02 0.174389e+02 0.174881e+02 0.175558e+02 0.176424e+02 0.177480e402 0.178723e+02 0.180149e+02 0.181750e+02 0.183515e+02 0.185433e+02 0.187487e+02 0.189660e+02 0.191932e+02 0.194279e*02 0.196677e+02 0.199101e~02 0.201521e+02 0.203908e+02 0.206232e+02 0.208463e+02 0.210568ef02 0.212518e+02 0.214282e+02 0.215834e+02 0,217148e+02 0.218204e+02 0.218988e402 0.219499e102 0.21981 2e+02 distribution function of species number 1 at normalized velocity xi 7)= 0.857143e-01 as a function of theta in radians from theta = 0. to theta = 0.314159e+01 0.168052e+02 0.168281ef02 0.168633e102 0.169177e+02 0.169929e+02 O.170895e+02 0.172079e+02 0.173480ef02 0.175094ef02 0.176915e+02 0.178933e+02 0.181135e+02 0.183504e+02 0.186021e-f02 0.188664e+02 0.191406e*02 0.194218e+02 0.197071e102 0.199931e+02 0.202763e+02 0.205529e102 O.208193e+02 0.210715e+02 0.213051e+02 0.215182e+02 0.217055e102 0.218645e102 0.219925e+02 0.220876e+02 0.221496e+02 0.221877e+02 distribution function of species number 1 at normalized velocity xi 81= 0.l00000e+00 as a function of theta in radians from theta = 0. to theta = 0.314159e+01 O.160630e-f02 0.160878e+02 0.161250e+02 0.1618288102 0.162628e+02 0.163662e+02 0.164936e+02 0.166451e+02 0.168207e+02 0.170197e102 0.172414e+02 0.174843e+02 O.177469e+02 0.180272e+02 0.183226e+02 0.186303e+02 0.189474e+02 0.192702e+02 0.195950e+02 0.199177e+0Z 0.202341e+02 0.205397e+02 0.208300e+02 0.211003e+02 0.213462e+02 0.215634e+02 0.217481e+02 0.218970e+02 0.220079e+02 0.220802e+02 8.221 248e4 02 distribution function of species number 1 at normalized velocity xi 91= 0.114286e+00 as a function of theta in radians from theta = 0. to theta = 0.31 41 59e+01 0.152153e+02 0.152411e102 0.152793e+02 0.153384e+02 0.154207e+02 O.155277e+02 0.156603e+02 0.158191e+02 0.160041e+02 0.162150e102 0.164511e+02 0.167111e+02 0.169934e+02 0.172960e+02 0.176164e+02 0.179517e+02 0.182983e+02 0.186527e+02 0.190106e+02 0.193675e+02 0.197186e+02 0.200589e+02 0.203830e+02 0.206856e+02 0.209616e102 0.212060e+02 0.214141e102 0.215823e+02 0.217075e+02 0.217893e+02 0.21 0398e+02 distribution function of species number 1 at normalized velocity xi 10)= as a function of theta cn radcans from theta = 0. to theta =

0.128571e100 0.314159e+O1

0.143016e+02 0.143279e+02 0.143657e+02 0.144243e+02 0.l45064e+02 0.146140e+02 0.147484e102 0.149104e+02 0.151003e+02 0.1531818102 0.155632e+02 0.158346e+02 0.161308e+02 0.164497e+02 0.167889e102 0.171453e+02 0.175153e+02 0.178951e+02 0.182801e+02 0.186655e+02 0.190458e+02 0.194156e+02 0.197689e102 0.200997e+02 0.204021e102 0.206704e+0Z 0.208994e+02 0.210846e+02 0.212228e+02 0.213131e+02 0.213689e402 distribution function of species number 1 at normalized velocity xi 11) 0.142857e100

A.A. Mirin et aL

/

FPPAC88

379

Test 2

output to fokker-planck package with driver, for questions contact either m. a. a. mirin of the national magnetic fusion energy computer center at the lawrence livermore national laboratory.

g. mccoy or

see subroutine initial for comments regarding input values. number of velocity mesh points ijx) = 101 number of theta mesh points iiy) = 20 number of general species inbo) = 1 number of maxwellian ions imeq) = 1 highest order le9endre polynomial (mx) = ksydm = 0. distribution functions assumed symmetric about pi/2. mex = 0, electrons are a general species. itest= 2 nstop= 30000 nprint= 7500 ncoef 100 xmax= 0.800000e+02 dtr= 0.200000e-03 secs vnorm= 0.600000e+09 cm/sec

8

general species follow species number mass bnumb=

1

0.910660e-27 grams 0.l00000e+01

initial values for this species reden 0.733000e+13 particles/cm**3 enit= 0. key xnit 0.511600e—02 ynit= 0. vnit= 0. description of loss term for this species tauo 0.l00000e+91 secs taul Q.l00000e+91 secs description of source for this species asor= 0. particles/cm**3/sec esor= 0.300000e+02 key xsor= 0.l00000e+02 ysor= 0. usor 0.

11 time step= species no. species no.

7500 1 2

time= density= denscty=

0.150000e+01 0.252686e413 0.733000e+13

secs particles/cm**3 particles/cm**3

energy= energy=

0.311152e+02 key 0.305000e+03 key

380

A.A. Mirin et aL

/

FPPAC88

distribution function of species number 1 at normalized velocity xi 11= 0. as a function of theta in radians from theta = 0. to theta = 0.157080e+01 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e-04 0.161794e—04 0.161794e—04 0.161794e—04 0.161794e-04 0.161794e—04 0. 161794e-04 0. 161794e-04 distribution function of species number 1 at normalized velocity xi 2)= as a function of theta in radians from theta = 0. to theta =

0.800000e400 0.157080e+01

0.161657e—04 Q.161657e—04 0.161657e—04 0.161657e—04 0.161657e-04 0.161657e—04 0.161657e—04 0.161656e—04 0.161656e—04 0.161656e—04 0.161656e—04 0.161656e—04 Q.161656e-04 0.161656e—04 0.161656e—04 0.161656e—04 0.161656e—04 0.161656e—04 0 .161656e-04 0.161656e-04 distribution function of species number 1 at normalized velocity xi 3) 0.160000e+01 as a function of theta in radians from theta = 0. to theta = 0.157080e+01 O.161388e-04 O.161388e-04 0.161388e-04 0.161388e-04 0.161388e-04 0.161388e-04 0.161388e-04 0.161388e—04 0.161387e-04 4 0.161387e-04 0.161387e-04 0.161387e-04 0.161387e-04 0.161387e—04 0.161387e-0 0.161387e-04 distribution function of specirs number 1 at normalized velocity xi 4) 0.2400008101 as a function of theta in radians from theta = 0. to theta = 0.157080e+01 0.160720e-04 0.160726e-04 0.160725e-04 0.160724e-04

0.160728e-04 0.160726e—04 0.160725e-04 0.160724e-04

0.160727e-04 0.160726e-04 0.160724e-04

0.160727e-04 0.160726e-04 0.160724e-04

0.160727e—04 0.160725e-04 0.160724e-04

distribution function of species number 1 at normalized velocity xi 5i as a Function of theta in radians from theta = 0. to theta =

0.160727e-04 0.160725e-04 0.160724e-04

0.320000e+01 0.157080e+01

0.159310e—04 0.159310e—04 0.159310e-04 0.159309e-04 0.159309e—04 0.159309e-04 0.159308e—04 0.159308e—04 0.159307e—04 0.159307e-04 0.159307e-04 0.159306e-04 0.159306e—04 0.159305e—04 0.159305e—04 0.159305e-04 0.159305e—04 0.159305e-04 0.159304e—04 0.159304e-04 distribution function of species number 1 at normalized velocity xi 6) 0.400000e+01 as a function of theta in radians from theta 0. to theta = 0.157080e+01 0.156953e-04 0.156953e-04 0.156953e-Q4 0.156952e-04 0.156952e-04 0.156951e-04 0.156951e—04 Q.156950e—04 0.156950e-04 0.156949e—04 0.156949e-04 0.156949e—04 0.156948e-04 0.156948e-04 0.156947e-04 0.156947e-04 0.156947e-04 0,156947e—04 0.156947e-04 0.156947e-04 distribution function of species number 1 at normalized velocity xi 7i 0.480000e+01 as a function of theta in radians from theta = 0. to theta = 0.157080e+01 0.153651e-04 O.153651e-04 0.153651e-04 0.153650e-04 0.153649e-04 0.153649e-04 0.153648e-04 0.153648e-04 0.153647e-04 0.153647e-04 0.153647e-04 0.153647e-04 0.153646e-04 0.153646e—04 0.153646e-04 0.153645e-04 0.153645e-04 0.153645e-04 0.153645e-04 0.153645e-04 distribution function of species number 1 at normalized velocity xi 81= 0.560000e+01 as a function of theta in radians from theta = 0. to theta = 0.157080e401 0.149481e-04 0.149481e-04 0.149481e-04 0.149480e-04 0.149480e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149479e-04 0.149480e-04 0.149480e-04 distribution function of species number 1 at normalized velocity xi 9)= 0.640000e+01 as a function of theta in radians from theta 0. to theta = 0.157080e+01 0.144549e-04 0.144549e-04 0.144549e-04 0.144548e-04 0.144548e-04 0.144548e-04 0.144548e-04 0.144548e-04 0.144549e-04 0.144550e-04 0.144551e-04 0.144552e-04 0.144553e-04 0.144553e-04 0.144554e-04 0.144554e-04 0.144555e-04 0.144555e-04 0.144556e-04 0.144556e-04 distribution function of species number 1 at normalized velocity xi 10)= 0.720000e101 as a function of theta in radians from theta 0. to theta = 0.157080e+01 0.138963e-04 0.138964e-04 0.138976e-04 0.130984e-04

0.138963e-04 0.138966e-04 0.138978e-04 0.138984e-04

0.138963e-04 0.138968e-04 0.138979e-04

0.138963e-04 0.138970e-04 0.138981e-04

0.138963e-04 0.138972e-04 0.138982e-04

distribution function of species number 1 at normalized velocity xi 111= as a function of theta in radians from theta = 0. to theta =

0.138964e-04 0.138974e-04 0.138983e-04

0.800000e101 0.157080e+01

0.132831e-04 0.132832e-04 0.132832e-04 0.132833e-04 0.132833e-04 0.132835e-04 0.132837e—04 0.132840e-04 0.132843e-04 0.132847e-04 0.132851e-04 0.132854e-04 0.132858e-04 0.132861e-04 0.132864e-04 0.132866e-04 0.132868e-04 0.132870e-04 0.132871e-04 0.132872e-04 distribution function of species number 1 at normalized velocity xi 121= 0.880000e+01 as a function of theta in radians from theta = 0. to theta = 0.157080e+01 0.126259e-04 0.126261e-04 0.126261e-04 0.126262e-04 0.126264e-04 0.126266e-04 0.126270e-04 0.126274e-04 0.126280e-04 0.126285e-04 0.126291e-04 0.126297e-04 0.126302e-04 0.126307e-04 0.126312e-04 0.126315e-04 0.126319e-04 0.l26321e-04 0.126323e-04 0.126323e-04 distribution function of species number 1 at normalized velocity xi 13)= 0.960000e+01 as a function of theta in radians from theta = 0. to theta = 0.157080e+01 0.119350e-04 0.119352e-04 0.119353e-04 0.119354e-04 0.119357e-04 0.119361e-04 0.119366e-04 0.119372e-04 0.119379e-04 0.119387e-04 0.119395e-04 0.119403e-04 0.119411e-04 0.119417e-04 0.119424e-04 0.119429e-04 0.119433e-04 0.119437e-04